Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.5.0/ALGLIB-2.5.0/spdgevd.cs @ 3839

Last change on this file since 3839 was 3839, checked in by mkommend, 14 years ago

implemented first version of LR (ticket #1012)

File size: 19.7 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 = evd.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            matinv.matinvreport rep = new matinv.matinvreport();
248            int info = 0;
249            int i_ = 0;
250            int i1_ = 0;
251
252            System.Diagnostics.Debug.Assert(n>0, "SMatrixGEVDReduce: N<=0!");
253            System.Diagnostics.Debug.Assert(problemtype==1 | problemtype==2 | problemtype==3, "SMatrixGEVDReduce: incorrect ProblemType!");
254            result = true;
255           
256            //
257            // Problem 1:  A*x = lambda*B*x
258            //
259            // Reducing to:
260            //     C*y = lambda*y
261            //     C = L^(-1) * A * L^(-T)
262            //     x = L^(-T) * y
263            //
264            if( problemtype==1 )
265            {
266               
267                //
268                // Factorize B in T: B = LL'
269                //
270                t = new double[n-1+1, n-1+1];
271                if( isupperb )
272                {
273                    for(i=0; i<=n-1; i++)
274                    {
275                        for(i_=i; i_<=n-1;i_++)
276                        {
277                            t[i_,i] = b[i,i_];
278                        }
279                    }
280                }
281                else
282                {
283                    for(i=0; i<=n-1; i++)
284                    {
285                        for(i_=0; i_<=i;i_++)
286                        {
287                            t[i,i_] = b[i,i_];
288                        }
289                    }
290                }
291                if( !trfac.spdmatrixcholesky(ref t, n, false) )
292                {
293                    result = false;
294                    return result;
295                }
296               
297                //
298                // Invert L in T
299                //
300                matinv.rmatrixtrinverse(ref t, n, false, false, ref info, ref rep);
301                if( info<=0 )
302                {
303                    result = false;
304                    return result;
305                }
306               
307                //
308                // Build L^(-1) * A * L^(-T) in R
309                //
310                w1 = new double[n+1];
311                w2 = new double[n+1];
312                r = new double[n-1+1, n-1+1];
313                for(j=1; j<=n; j++)
314                {
315                   
316                    //
317                    // Form w2 = A * l'(j) (here l'(j) is j-th column of L^(-T))
318                    //
319                    i1_ = (0) - (1);
320                    for(i_=1; i_<=j;i_++)
321                    {
322                        w1[i_] = t[j-1,i_+i1_];
323                    }
324                    sblas.symmetricmatrixvectormultiply(ref a, isuppera, 0, j-1, ref w1, 1.0, ref w2);
325                    if( isuppera )
326                    {
327                        blas.matrixvectormultiply(ref a, 0, j-1, j, n-1, true, ref w1, 1, j, 1.0, ref w2, j+1, n, 0.0);
328                    }
329                    else
330                    {
331                        blas.matrixvectormultiply(ref a, j, n-1, 0, j-1, false, ref w1, 1, j, 1.0, ref w2, j+1, n, 0.0);
332                    }
333                   
334                    //
335                    // Form l(i)*w2 (here l(i) is i-th row of L^(-1))
336                    //
337                    for(i=1; i<=n; i++)
338                    {
339                        i1_ = (1)-(0);
340                        v = 0.0;
341                        for(i_=0; i_<=i-1;i_++)
342                        {
343                            v += t[i-1,i_]*w2[i_+i1_];
344                        }
345                        r[i-1,j-1] = v;
346                    }
347                }
348               
349                //
350                // Copy R to A
351                //
352                for(i=0; i<=n-1; i++)
353                {
354                    for(i_=0; i_<=n-1;i_++)
355                    {
356                        a[i,i_] = r[i,i_];
357                    }
358                }
359               
360                //
361                // Copy L^(-1) from T to R and transpose
362                //
363                isupperr = true;
364                for(i=0; i<=n-1; i++)
365                {
366                    for(j=0; j<=i-1; j++)
367                    {
368                        r[i,j] = 0;
369                    }
370                }
371                for(i=0; i<=n-1; i++)
372                {
373                    for(i_=i; i_<=n-1;i_++)
374                    {
375                        r[i,i_] = t[i_,i];
376                    }
377                }
378                return result;
379            }
380           
381            //
382            // Problem 2:  A*B*x = lambda*x
383            // or
384            // problem 3:  B*A*x = lambda*x
385            //
386            // Reducing to:
387            //     C*y = lambda*y
388            //     C = U * A * U'
389            //     B = U'* U
390            //
391            if( problemtype==2 | problemtype==3 )
392            {
393               
394                //
395                // Factorize B in T: B = U'*U
396                //
397                t = new double[n-1+1, n-1+1];
398                if( isupperb )
399                {
400                    for(i=0; i<=n-1; i++)
401                    {
402                        for(i_=i; i_<=n-1;i_++)
403                        {
404                            t[i,i_] = b[i,i_];
405                        }
406                    }
407                }
408                else
409                {
410                    for(i=0; i<=n-1; i++)
411                    {
412                        for(i_=i; i_<=n-1;i_++)
413                        {
414                            t[i,i_] = b[i_,i];
415                        }
416                    }
417                }
418                if( !trfac.spdmatrixcholesky(ref t, n, true) )
419                {
420                    result = false;
421                    return result;
422                }
423               
424                //
425                // Build U * A * U' in R
426                //
427                w1 = new double[n+1];
428                w2 = new double[n+1];
429                w3 = new double[n+1];
430                r = new double[n-1+1, n-1+1];
431                for(j=1; j<=n; j++)
432                {
433                   
434                    //
435                    // Form w2 = A * u'(j) (here u'(j) is j-th column of U')
436                    //
437                    i1_ = (j-1) - (1);
438                    for(i_=1; i_<=n-j+1;i_++)
439                    {
440                        w1[i_] = t[j-1,i_+i1_];
441                    }
442                    sblas.symmetricmatrixvectormultiply(ref a, isuppera, j-1, n-1, ref w1, 1.0, ref w3);
443                    i1_ = (1) - (j);
444                    for(i_=j; i_<=n;i_++)
445                    {
446                        w2[i_] = w3[i_+i1_];
447                    }
448                    i1_ = (j-1) - (j);
449                    for(i_=j; i_<=n;i_++)
450                    {
451                        w1[i_] = t[j-1,i_+i1_];
452                    }
453                    if( isuppera )
454                    {
455                        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);
456                    }
457                    else
458                    {
459                        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);
460                    }
461                   
462                    //
463                    // Form u(i)*w2 (here u(i) is i-th row of U)
464                    //
465                    for(i=1; i<=n; i++)
466                    {
467                        i1_ = (i)-(i-1);
468                        v = 0.0;
469                        for(i_=i-1; i_<=n-1;i_++)
470                        {
471                            v += t[i-1,i_]*w2[i_+i1_];
472                        }
473                        r[i-1,j-1] = v;
474                    }
475                }
476               
477                //
478                // Copy R to A
479                //
480                for(i=0; i<=n-1; i++)
481                {
482                    for(i_=0; i_<=n-1;i_++)
483                    {
484                        a[i,i_] = r[i,i_];
485                    }
486                }
487                if( problemtype==2 )
488                {
489                   
490                    //
491                    // Invert U in T
492                    //
493                    matinv.rmatrixtrinverse(ref t, n, true, false, ref info, ref rep);
494                    if( info<=0 )
495                    {
496                        result = false;
497                        return result;
498                    }
499                   
500                    //
501                    // Copy U^-1 from T to R
502                    //
503                    isupperr = true;
504                    for(i=0; i<=n-1; i++)
505                    {
506                        for(j=0; j<=i-1; j++)
507                        {
508                            r[i,j] = 0;
509                        }
510                    }
511                    for(i=0; i<=n-1; i++)
512                    {
513                        for(i_=i; i_<=n-1;i_++)
514                        {
515                            r[i,i_] = t[i,i_];
516                        }
517                    }
518                }
519                else
520                {
521                   
522                    //
523                    // Copy U from T to R and transpose
524                    //
525                    isupperr = false;
526                    for(i=0; i<=n-1; i++)
527                    {
528                        for(j=i+1; j<=n-1; j++)
529                        {
530                            r[i,j] = 0;
531                        }
532                    }
533                    for(i=0; i<=n-1; i++)
534                    {
535                        for(i_=i; i_<=n-1;i_++)
536                        {
537                            r[i_,i] = t[i,i_];
538                        }
539                    }
540                }
541            }
542            return result;
543        }
544    }
545}
Note: See TracBrowser for help on using the repository browser.