Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/ALGLIB/spdgevd.cs @ 2575

Last change on this file since 2575 was 2563, checked in by gkronber, 15 years ago

Updated ALGLIB to latest version. #751 (Plugin for for data-modeling with ANN (integrated into CEDMA))

File size: 32.2 KB
Line 
1/*************************************************************************
2Copyright (c) 2005-2007, Sergey Bochkanov (ALGLIB project).
3
4>>> SOURCE LICENSE >>>
5This program is free software; you can redistribute it and/or modify
6it under the terms of the GNU General Public License as published by
7the Free Software Foundation (www.fsf.org); either version 2 of the
8License, or (at your option) any later version.
9
10This program is distributed in the hope that it will be useful,
11but WITHOUT ANY WARRANTY; without even the implied warranty of
12MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13GNU General Public License for more details.
14
15A copy of the GNU General Public License is available at
16http://www.fsf.org/licensing/licenses
17
18>>> END OF LICENSE >>>
19*************************************************************************/
20
21using System;
22
23namespace alglib
24{
25    public class spdgevd
26    {
27        /*************************************************************************
28        Algorithm for solving the following generalized symmetric positive-definite
29        eigenproblem:
30            A*x = lambda*B*x (1) or
31            A*B*x = lambda*x (2) or
32            B*A*x = lambda*x (3).
33        where A is a symmetric matrix, B - symmetric positive-definite matrix.
34        The problem is solved by reducing it to an ordinary  symmetric  eigenvalue
35        problem.
36
37        Input parameters:
38            A           -   symmetric matrix which is given by its upper or lower
39                            triangular part.
40                            Array whose indexes range within [0..N-1, 0..N-1].
41            N           -   size of matrices A and B.
42            IsUpperA    -   storage format of matrix A.
43            B           -   symmetric positive-definite matrix which is given by
44                            its upper or lower triangular part.
45                            Array whose indexes range within [0..N-1, 0..N-1].
46            IsUpperB    -   storage format of matrix B.
47            ZNeeded     -   if ZNeeded is equal to:
48                             * 0, the eigenvectors are not returned;
49                             * 1, the eigenvectors are returned.
50            ProblemType -   if ProblemType is equal to:
51                             * 1, the following problem is solved: A*x = lambda*B*x;
52                             * 2, the following problem is solved: A*B*x = lambda*x;
53                             * 3, the following problem is solved: B*A*x = lambda*x.
54
55        Output parameters:
56            D           -   eigenvalues in ascending order.
57                            Array whose index ranges within [0..N-1].
58            Z           -   if ZNeeded is equal to:
59                             * 0, Z hasn’t changed;
60                             * 1, Z contains eigenvectors.
61                            Array whose indexes range within [0..N-1, 0..N-1].
62                            The eigenvectors are stored in matrix columns. It should
63                            be noted that the eigenvectors in such problems do not
64                            form an orthogonal system.
65
66        Result:
67            True, if the problem was solved successfully.
68            False, if the error occurred during the Cholesky decomposition of matrix
69            B (the matrix isn’t positive-definite) or during the work of the iterative
70            algorithm for solving the symmetric eigenproblem.
71
72        See also the GeneralizedSymmetricDefiniteEVDReduce subroutine.
73
74          -- ALGLIB --
75             Copyright 1.28.2006 by Bochkanov Sergey
76        *************************************************************************/
77        public static bool smatrixgevd(double[,] a,
78            int n,
79            bool isuppera,
80            ref double[,] b,
81            bool isupperb,
82            int zneeded,
83            int problemtype,
84            ref double[] d,
85            ref double[,] z)
86        {
87            bool result = new bool();
88            double[,] r = new double[0,0];
89            double[,] t = new double[0,0];
90            bool isupperr = new bool();
91            int j1 = 0;
92            int j2 = 0;
93            int j1inc = 0;
94            int j2inc = 0;
95            int i = 0;
96            int j = 0;
97            double v = 0;
98            int i_ = 0;
99
100            a = (double[,])a.Clone();
101
102           
103            //
104            // Reduce and solve
105            //
106            result = smatrixgevdreduce(ref a, n, isuppera, ref b, isupperb, problemtype, ref r, ref isupperr);
107            if( !result )
108            {
109                return result;
110            }
111            result = sevd.smatrixevd(a, n, zneeded, isuppera, ref d, ref t);
112            if( !result )
113            {
114                return result;
115            }
116           
117            //
118            // Transform eigenvectors if needed
119            //
120            if( zneeded!=0 )
121            {
122               
123                //
124                // fill Z with zeros
125                //
126                z = new double[n-1+1, n-1+1];
127                for(j=0; j<=n-1; j++)
128                {
129                    z[0,j] = 0.0;
130                }
131                for(i=1; i<=n-1; i++)
132                {
133                    for(i_=0; i_<=n-1;i_++)
134                    {
135                        z[i,i_] = z[0,i_];
136                    }
137                }
138               
139                //
140                // Setup R properties
141                //
142                if( isupperr )
143                {
144                    j1 = 0;
145                    j2 = n-1;
146                    j1inc = +1;
147                    j2inc = 0;
148                }
149                else
150                {
151                    j1 = 0;
152                    j2 = 0;
153                    j1inc = 0;
154                    j2inc = +1;
155                }
156               
157                //
158                // Calculate R*Z
159                //
160                for(i=0; i<=n-1; i++)
161                {
162                    for(j=j1; j<=j2; j++)
163                    {
164                        v = r[i,j];
165                        for(i_=0; i_<=n-1;i_++)
166                        {
167                            z[i,i_] = z[i,i_] + v*t[j,i_];
168                        }
169                    }
170                    j1 = j1+j1inc;
171                    j2 = j2+j2inc;
172                }
173            }
174            return result;
175        }
176
177
178        /*************************************************************************
179        Algorithm for reduction of the following generalized symmetric positive-
180        definite eigenvalue problem:
181            A*x = lambda*B*x (1) or
182            A*B*x = lambda*x (2) or
183            B*A*x = lambda*x (3)
184        to the symmetric eigenvalues problem C*y = lambda*y (eigenvalues of this and
185        the given problems are the same, and the eigenvectors of the given problem
186        could be obtained by multiplying the obtained eigenvectors by the
187        transformation matrix x = R*y).
188
189        Here A is a symmetric matrix, B - symmetric positive-definite matrix.
190
191        Input parameters:
192            A           -   symmetric matrix which is given by its upper or lower
193                            triangular part.
194                            Array whose indexes range within [0..N-1, 0..N-1].
195            N           -   size of matrices A and B.
196            IsUpperA    -   storage format of matrix A.
197            B           -   symmetric positive-definite matrix which is given by
198                            its upper or lower triangular part.
199                            Array whose indexes range within [0..N-1, 0..N-1].
200            IsUpperB    -   storage format of matrix B.
201            ProblemType -   if ProblemType is equal to:
202                             * 1, the following problem is solved: A*x = lambda*B*x;
203                             * 2, the following problem is solved: A*B*x = lambda*x;
204                             * 3, the following problem is solved: B*A*x = lambda*x.
205
206        Output parameters:
207            A           -   symmetric matrix which is given by its upper or lower
208                            triangle depending on IsUpperA. Contains matrix C.
209                            Array whose indexes range within [0..N-1, 0..N-1].
210            R           -   upper triangular or low triangular transformation matrix
211                            which is used to obtain the eigenvectors of a given problem
212                            as the product of eigenvectors of C (from the right) and
213                            matrix R (from the left). If the matrix is upper
214                            triangular, the elements below the main diagonal
215                            are equal to 0 (and vice versa). Thus, we can perform
216                            the multiplication without taking into account the
217                            internal structure (which is an easier though less
218                            effective way).
219                            Array whose indexes range within [0..N-1, 0..N-1].
220            IsUpperR    -   type of matrix R (upper or lower triangular).
221
222        Result:
223            True, if the problem was reduced successfully.
224            False, if the error occurred during the Cholesky decomposition of
225                matrix B (the matrix is not positive-definite).
226
227          -- ALGLIB --
228             Copyright 1.28.2006 by Bochkanov Sergey
229        *************************************************************************/
230        public static bool smatrixgevdreduce(ref double[,] a,
231            int n,
232            bool isuppera,
233            ref double[,] b,
234            bool isupperb,
235            int problemtype,
236            ref double[,] r,
237            ref bool isupperr)
238        {
239            bool result = new bool();
240            double[,] t = new double[0,0];
241            double[] w1 = new double[0];
242            double[] w2 = new double[0];
243            double[] w3 = new double[0];
244            int i = 0;
245            int j = 0;
246            double v = 0;
247            int i_ = 0;
248            int i1_ = 0;
249
250            System.Diagnostics.Debug.Assert(n>0, "SMatrixGEVDReduce: N<=0!");
251            System.Diagnostics.Debug.Assert(problemtype==1 | problemtype==2 | problemtype==3, "SMatrixGEVDReduce: incorrect ProblemType!");
252            result = true;
253           
254            //
255            // Problem 1:  A*x = lambda*B*x
256            //
257            // Reducing to:
258            //     C*y = lambda*y
259            //     C = L^(-1) * A * L^(-T)
260            //     x = L^(-T) * y
261            //
262            if( problemtype==1 )
263            {
264               
265                //
266                // Factorize B in T: B = LL'
267                //
268                t = new double[n-1+1, n-1+1];
269                if( isupperb )
270                {
271                    for(i=0; i<=n-1; i++)
272                    {
273                        for(i_=i; i_<=n-1;i_++)
274                        {
275                            t[i_,i] = b[i,i_];
276                        }
277                    }
278                }
279                else
280                {
281                    for(i=0; i<=n-1; i++)
282                    {
283                        for(i_=0; i_<=i;i_++)
284                        {
285                            t[i,i_] = b[i,i_];
286                        }
287                    }
288                }
289                if( !cholesky.spdmatrixcholesky(ref t, n, false) )
290                {
291                    result = false;
292                    return result;
293                }
294               
295                //
296                // Invert L in T
297                //
298                if( !trinverse.rmatrixtrinverse(ref t, n, false, false) )
299                {
300                    result = false;
301                    return result;
302                }
303               
304                //
305                // Build L^(-1) * A * L^(-T) in R
306                //
307                w1 = new double[n+1];
308                w2 = new double[n+1];
309                r = new double[n-1+1, n-1+1];
310                for(j=1; j<=n; j++)
311                {
312                   
313                    //
314                    // Form w2 = A * l'(j) (here l'(j) is j-th column of L^(-T))
315                    //
316                    i1_ = (0) - (1);
317                    for(i_=1; i_<=j;i_++)
318                    {
319                        w1[i_] = t[j-1,i_+i1_];
320                    }
321                    sblas.symmetricmatrixvectormultiply(ref a, isuppera, 0, j-1, ref w1, 1.0, ref w2);
322                    if( isuppera )
323                    {
324                        blas.matrixvectormultiply(ref a, 0, j-1, j, n-1, true, ref w1, 1, j, 1.0, ref w2, j+1, n, 0.0);
325                    }
326                    else
327                    {
328                        blas.matrixvectormultiply(ref a, j, n-1, 0, j-1, false, ref w1, 1, j, 1.0, ref w2, j+1, n, 0.0);
329                    }
330                   
331                    //
332                    // Form l(i)*w2 (here l(i) is i-th row of L^(-1))
333                    //
334                    for(i=1; i<=n; i++)
335                    {
336                        i1_ = (1)-(0);
337                        v = 0.0;
338                        for(i_=0; i_<=i-1;i_++)
339                        {
340                            v += t[i-1,i_]*w2[i_+i1_];
341                        }
342                        r[i-1,j-1] = v;
343                    }
344                }
345               
346                //
347                // Copy R to A
348                //
349                for(i=0; i<=n-1; i++)
350                {
351                    for(i_=0; i_<=n-1;i_++)
352                    {
353                        a[i,i_] = r[i,i_];
354                    }
355                }
356               
357                //
358                // Copy L^(-1) from T to R and transpose
359                //
360                isupperr = true;
361                for(i=0; i<=n-1; i++)
362                {
363                    for(j=0; j<=i-1; j++)
364                    {
365                        r[i,j] = 0;
366                    }
367                }
368                for(i=0; i<=n-1; i++)
369                {
370                    for(i_=i; i_<=n-1;i_++)
371                    {
372                        r[i,i_] = t[i_,i];
373                    }
374                }
375                return result;
376            }
377           
378            //
379            // Problem 2:  A*B*x = lambda*x
380            // or
381            // problem 3:  B*A*x = lambda*x
382            //
383            // Reducing to:
384            //     C*y = lambda*y
385            //     C = U * A * U'
386            //     B = U'* U
387            //
388            if( problemtype==2 | problemtype==3 )
389            {
390               
391                //
392                // Factorize B in T: B = U'*U
393                //
394                t = new double[n-1+1, n-1+1];
395                if( isupperb )
396                {
397                    for(i=0; i<=n-1; i++)
398                    {
399                        for(i_=i; i_<=n-1;i_++)
400                        {
401                            t[i,i_] = b[i,i_];
402                        }
403                    }
404                }
405                else
406                {
407                    for(i=0; i<=n-1; i++)
408                    {
409                        for(i_=i; i_<=n-1;i_++)
410                        {
411                            t[i,i_] = b[i_,i];
412                        }
413                    }
414                }
415                if( !cholesky.spdmatrixcholesky(ref t, n, true) )
416                {
417                    result = false;
418                    return result;
419                }
420               
421                //
422                // Build U * A * U' in R
423                //
424                w1 = new double[n+1];
425                w2 = new double[n+1];
426                w3 = new double[n+1];
427                r = new double[n-1+1, n-1+1];
428                for(j=1; j<=n; j++)
429                {
430                   
431                    //
432                    // Form w2 = A * u'(j) (here u'(j) is j-th column of U')
433                    //
434                    i1_ = (j-1) - (1);
435                    for(i_=1; i_<=n-j+1;i_++)
436                    {
437                        w1[i_] = t[j-1,i_+i1_];
438                    }
439                    sblas.symmetricmatrixvectormultiply(ref a, isuppera, j-1, n-1, ref w1, 1.0, ref w3);
440                    i1_ = (1) - (j);
441                    for(i_=j; i_<=n;i_++)
442                    {
443                        w2[i_] = w3[i_+i1_];
444                    }
445                    i1_ = (j-1) - (j);
446                    for(i_=j; i_<=n;i_++)
447                    {
448                        w1[i_] = t[j-1,i_+i1_];
449                    }
450                    if( isuppera )
451                    {
452                        blas.matrixvectormultiply(ref a, 0, j-2, j-1, n-1, false, ref w1, j, n, 1.0, ref w2, 1, j-1, 0.0);
453                    }
454                    else
455                    {
456                        blas.matrixvectormultiply(ref a, j-1, n-1, 0, j-2, true, ref w1, j, n, 1.0, ref w2, 1, j-1, 0.0);
457                    }
458                   
459                    //
460                    // Form u(i)*w2 (here u(i) is i-th row of U)
461                    //
462                    for(i=1; i<=n; i++)
463                    {
464                        i1_ = (i)-(i-1);
465                        v = 0.0;
466                        for(i_=i-1; i_<=n-1;i_++)
467                        {
468                            v += t[i-1,i_]*w2[i_+i1_];
469                        }
470                        r[i-1,j-1] = v;
471                    }
472                }
473               
474                //
475                // Copy R to A
476                //
477                for(i=0; i<=n-1; i++)
478                {
479                    for(i_=0; i_<=n-1;i_++)
480                    {
481                        a[i,i_] = r[i,i_];
482                    }
483                }
484                if( problemtype==2 )
485                {
486                   
487                    //
488                    // Invert U in T
489                    //
490                    if( !trinverse.rmatrixtrinverse(ref t, n, true, false) )
491                    {
492                        result = false;
493                        return result;
494                    }
495                   
496                    //
497                    // Copy U^-1 from T to R
498                    //
499                    isupperr = true;
500                    for(i=0; i<=n-1; i++)
501                    {
502                        for(j=0; j<=i-1; j++)
503                        {
504                            r[i,j] = 0;
505                        }
506                    }
507                    for(i=0; i<=n-1; i++)
508                    {
509                        for(i_=i; i_<=n-1;i_++)
510                        {
511                            r[i,i_] = t[i,i_];
512                        }
513                    }
514                }
515                else
516                {
517                   
518                    //
519                    // Copy U from T to R and transpose
520                    //
521                    isupperr = false;
522                    for(i=0; i<=n-1; i++)
523                    {
524                        for(j=i+1; j<=n-1; j++)
525                        {
526                            r[i,j] = 0;
527                        }
528                    }
529                    for(i=0; i<=n-1; i++)
530                    {
531                        for(i_=i; i_<=n-1;i_++)
532                        {
533                            r[i_,i] = t[i,i_];
534                        }
535                    }
536                }
537            }
538            return result;
539        }
540
541
542        public static bool generalizedsymmetricdefiniteevd(double[,] a,
543            int n,
544            bool isuppera,
545            ref double[,] b,
546            bool isupperb,
547            int zneeded,
548            int problemtype,
549            ref double[] d,
550            ref double[,] z)
551        {
552            bool result = new bool();
553            double[,] r = new double[0,0];
554            double[,] t = new double[0,0];
555            bool isupperr = new bool();
556            int j1 = 0;
557            int j2 = 0;
558            int j1inc = 0;
559            int j2inc = 0;
560            int i = 0;
561            int j = 0;
562            double v = 0;
563            int i_ = 0;
564
565            a = (double[,])a.Clone();
566
567           
568            //
569            // Reduce and solve
570            //
571            result = generalizedsymmetricdefiniteevdreduce(ref a, n, isuppera, ref b, isupperb, problemtype, ref r, ref isupperr);
572            if( !result )
573            {
574                return result;
575            }
576            result = sevd.symmetricevd(a, n, zneeded, isuppera, ref d, ref t);
577            if( !result )
578            {
579                return result;
580            }
581           
582            //
583            // Transform eigenvectors if needed
584            //
585            if( zneeded!=0 )
586            {
587               
588                //
589                // fill Z with zeros
590                //
591                z = new double[n+1, n+1];
592                for(j=1; j<=n; j++)
593                {
594                    z[1,j] = 0.0;
595                }
596                for(i=2; i<=n; i++)
597                {
598                    for(i_=1; i_<=n;i_++)
599                    {
600                        z[i,i_] = z[1,i_];
601                    }
602                }
603               
604                //
605                // Setup R properties
606                //
607                if( isupperr )
608                {
609                    j1 = 1;
610                    j2 = n;
611                    j1inc = +1;
612                    j2inc = 0;
613                }
614                else
615                {
616                    j1 = 1;
617                    j2 = 1;
618                    j1inc = 0;
619                    j2inc = +1;
620                }
621               
622                //
623                // Calculate R*Z
624                //
625                for(i=1; i<=n; i++)
626                {
627                    for(j=j1; j<=j2; j++)
628                    {
629                        v = r[i,j];
630                        for(i_=1; i_<=n;i_++)
631                        {
632                            z[i,i_] = z[i,i_] + v*t[j,i_];
633                        }
634                    }
635                    j1 = j1+j1inc;
636                    j2 = j2+j2inc;
637                }
638            }
639            return result;
640        }
641
642
643        public static bool generalizedsymmetricdefiniteevdreduce(ref double[,] a,
644            int n,
645            bool isuppera,
646            ref double[,] b,
647            bool isupperb,
648            int problemtype,
649            ref double[,] r,
650            ref bool isupperr)
651        {
652            bool result = new bool();
653            double[,] t = new double[0,0];
654            double[] w1 = new double[0];
655            double[] w2 = new double[0];
656            double[] w3 = new double[0];
657            int i = 0;
658            int j = 0;
659            double v = 0;
660            int i_ = 0;
661            int i1_ = 0;
662
663            System.Diagnostics.Debug.Assert(n>0, "GeneralizedSymmetricDefiniteEVDReduce: N<=0!");
664            System.Diagnostics.Debug.Assert(problemtype==1 | problemtype==2 | problemtype==3, "GeneralizedSymmetricDefiniteEVDReduce: incorrect ProblemType!");
665            result = true;
666           
667            //
668            // Problem 1:  A*x = lambda*B*x
669            //
670            // Reducing to:
671            //     C*y = lambda*y
672            //     C = L^(-1) * A * L^(-T)
673            //     x = L^(-T) * y
674            //
675            if( problemtype==1 )
676            {
677               
678                //
679                // Factorize B in T: B = LL'
680                //
681                t = new double[n+1, n+1];
682                if( isupperb )
683                {
684                    for(i=1; i<=n; i++)
685                    {
686                        for(i_=i; i_<=n;i_++)
687                        {
688                            t[i_,i] = b[i,i_];
689                        }
690                    }
691                }
692                else
693                {
694                    for(i=1; i<=n; i++)
695                    {
696                        for(i_=1; i_<=i;i_++)
697                        {
698                            t[i,i_] = b[i,i_];
699                        }
700                    }
701                }
702                if( !cholesky.choleskydecomposition(ref t, n, false) )
703                {
704                    result = false;
705                    return result;
706                }
707               
708                //
709                // Invert L in T
710                //
711                if( !trinverse.invtriangular(ref t, n, false, false) )
712                {
713                    result = false;
714                    return result;
715                }
716               
717                //
718                // Build L^(-1) * A * L^(-T) in R
719                //
720                w1 = new double[n+1];
721                w2 = new double[n+1];
722                r = new double[n+1, n+1];
723                for(j=1; j<=n; j++)
724                {
725                   
726                    //
727                    // Form w2 = A * l'(j) (here l'(j) is j-th column of L^(-T))
728                    //
729                    for(i_=1; i_<=j;i_++)
730                    {
731                        w1[i_] = t[j,i_];
732                    }
733                    sblas.symmetricmatrixvectormultiply(ref a, isuppera, 1, j, ref w1, 1.0, ref w2);
734                    if( isuppera )
735                    {
736                        blas.matrixvectormultiply(ref a, 1, j, j+1, n, true, ref w1, 1, j, 1.0, ref w2, j+1, n, 0.0);
737                    }
738                    else
739                    {
740                        blas.matrixvectormultiply(ref a, j+1, n, 1, j, false, ref w1, 1, j, 1.0, ref w2, j+1, n, 0.0);
741                    }
742                   
743                    //
744                    // Form l(i)*w2 (here l(i) is i-th row of L^(-1))
745                    //
746                    for(i=1; i<=n; i++)
747                    {
748                        v = 0.0;
749                        for(i_=1; i_<=i;i_++)
750                        {
751                            v += t[i,i_]*w2[i_];
752                        }
753                        r[i,j] = v;
754                    }
755                }
756               
757                //
758                // Copy R to A
759                //
760                for(i=1; i<=n; i++)
761                {
762                    for(i_=1; i_<=n;i_++)
763                    {
764                        a[i,i_] = r[i,i_];
765                    }
766                }
767               
768                //
769                // Copy L^(-1) from T to R and transpose
770                //
771                isupperr = true;
772                for(i=1; i<=n; i++)
773                {
774                    for(j=1; j<=i-1; j++)
775                    {
776                        r[i,j] = 0;
777                    }
778                }
779                for(i=1; i<=n; i++)
780                {
781                    for(i_=i; i_<=n;i_++)
782                    {
783                        r[i,i_] = t[i_,i];
784                    }
785                }
786                return result;
787            }
788           
789            //
790            // Problem 2:  A*B*x = lambda*x
791            // or
792            // problem 3:  B*A*x = lambda*x
793            //
794            // Reducing to:
795            //     C*y = lambda*y
796            //     C = U * A * U'
797            //     B = U'* U
798            //
799            if( problemtype==2 | problemtype==3 )
800            {
801               
802                //
803                // Factorize B in T: B = U'*U
804                //
805                t = new double[n+1, n+1];
806                if( isupperb )
807                {
808                    for(i=1; i<=n; i++)
809                    {
810                        for(i_=i; i_<=n;i_++)
811                        {
812                            t[i,i_] = b[i,i_];
813                        }
814                    }
815                }
816                else
817                {
818                    for(i=1; i<=n; i++)
819                    {
820                        for(i_=i; i_<=n;i_++)
821                        {
822                            t[i,i_] = b[i_,i];
823                        }
824                    }
825                }
826                if( !cholesky.choleskydecomposition(ref t, n, true) )
827                {
828                    result = false;
829                    return result;
830                }
831               
832                //
833                // Build U * A * U' in R
834                //
835                w1 = new double[n+1];
836                w2 = new double[n+1];
837                w3 = new double[n+1];
838                r = new double[n+1, n+1];
839                for(j=1; j<=n; j++)
840                {
841                   
842                    //
843                    // Form w2 = A * u'(j) (here u'(j) is j-th column of U')
844                    //
845                    i1_ = (j) - (1);
846                    for(i_=1; i_<=n-j+1;i_++)
847                    {
848                        w1[i_] = t[j,i_+i1_];
849                    }
850                    sblas.symmetricmatrixvectormultiply(ref a, isuppera, j, n, ref w1, 1.0, ref w3);
851                    i1_ = (1) - (j);
852                    for(i_=j; i_<=n;i_++)
853                    {
854                        w2[i_] = w3[i_+i1_];
855                    }
856                    for(i_=j; i_<=n;i_++)
857                    {
858                        w1[i_] = t[j,i_];
859                    }
860                    if( isuppera )
861                    {
862                        blas.matrixvectormultiply(ref a, 1, j-1, j, n, false, ref w1, j, n, 1.0, ref w2, 1, j-1, 0.0);
863                    }
864                    else
865                    {
866                        blas.matrixvectormultiply(ref a, j, n, 1, j-1, true, ref w1, j, n, 1.0, ref w2, 1, j-1, 0.0);
867                    }
868                   
869                    //
870                    // Form u(i)*w2 (here u(i) is i-th row of U)
871                    //
872                    for(i=1; i<=n; i++)
873                    {
874                        v = 0.0;
875                        for(i_=i; i_<=n;i_++)
876                        {
877                            v += t[i,i_]*w2[i_];
878                        }
879                        r[i,j] = v;
880                    }
881                }
882               
883                //
884                // Copy R to A
885                //
886                for(i=1; i<=n; i++)
887                {
888                    for(i_=1; i_<=n;i_++)
889                    {
890                        a[i,i_] = r[i,i_];
891                    }
892                }
893                if( problemtype==2 )
894                {
895                   
896                    //
897                    // Invert U in T
898                    //
899                    if( !trinverse.invtriangular(ref t, n, true, false) )
900                    {
901                        result = false;
902                        return result;
903                    }
904                   
905                    //
906                    // Copy U^-1 from T to R
907                    //
908                    isupperr = true;
909                    for(i=1; i<=n; i++)
910                    {
911                        for(j=1; j<=i-1; j++)
912                        {
913                            r[i,j] = 0;
914                        }
915                    }
916                    for(i=1; i<=n; i++)
917                    {
918                        for(i_=i; i_<=n;i_++)
919                        {
920                            r[i,i_] = t[i,i_];
921                        }
922                    }
923                }
924                else
925                {
926                   
927                    //
928                    // Copy U from T to R and transpose
929                    //
930                    isupperr = false;
931                    for(i=1; i<=n; i++)
932                    {
933                        for(j=i+1; j<=n; j++)
934                        {
935                            r[i,j] = 0;
936                        }
937                    }
938                    for(i=1; i<=n; i++)
939                    {
940                        for(i_=i; i_<=n;i_++)
941                        {
942                            r[i_,i] = t[i,i_];
943                        }
944                    }
945                }
946            }
947            return result;
948        }
949    }
950}
Note: See TracBrowser for help on using the repository browser.