Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.1.2.2591/ALGLIB-2.1.2.2591/htridiagonal.cs @ 2645

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

extracted external libraries and adapted dependent plugins (ticket #837)

File size: 25.4 KB
Line 
1/*************************************************************************
2Copyright (c) 1992-2007 The University of Tennessee.  All rights reserved.
3
4Contributors:
5    * Sergey Bochkanov (ALGLIB project). Translation from FORTRAN to
6      pseudocode.
7
8See subroutines comments for additional copyrights.
9
10>>> SOURCE LICENSE >>>
11This program is free software; you can redistribute it and/or modify
12it under the terms of the GNU General Public License as published by
13the Free Software Foundation (www.fsf.org); either version 2 of the
14License, or (at your option) any later version.
15
16This program is distributed in the hope that it will be useful,
17but WITHOUT ANY WARRANTY; without even the implied warranty of
18MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
19GNU General Public License for more details.
20
21A copy of the GNU General Public License is available at
22http://www.fsf.org/licensing/licenses
23
24>>> END OF LICENSE >>>
25*************************************************************************/
26
27using System;
28
29namespace alglib
30{
31    public class htridiagonal
32    {
33        /*************************************************************************
34        Reduction of a Hermitian matrix which is given  by  its  higher  or  lower
35        triangular part to a real  tridiagonal  matrix  using  unitary  similarity
36        transformation: Q'*A*Q = T.
37
38        Input parameters:
39            A       -   matrix to be transformed
40                        array with elements [0..N-1, 0..N-1].
41            N       -   size of matrix A.
42            IsUpper -   storage format. If IsUpper = True, then matrix A is  given
43                        by its upper triangle, and the lower triangle is not  used
44                        and not modified by the algorithm, and vice versa
45                        if IsUpper = False.
46
47        Output parameters:
48            A       -   matrices T and Q in  compact form (see lower)
49            Tau     -   array of factors which are forming matrices H(i)
50                        array with elements [0..N-2].
51            D       -   main diagonal of real symmetric matrix T.
52                        array with elements [0..N-1].
53            E       -   secondary diagonal of real symmetric matrix T.
54                        array with elements [0..N-2].
55
56
57          If IsUpper=True, the matrix Q is represented as a product of elementary
58          reflectors
59
60             Q = H(n-2) . . . H(2) H(0).
61
62          Each H(i) has the form
63
64             H(i) = I - tau * v * v'
65
66          where tau is a complex scalar, and v is a complex vector with
67          v(i+1:n-1) = 0, v(i) = 1, v(0:i-1) is stored on exit in
68          A(0:i-1,i+1), and tau in TAU(i).
69
70          If IsUpper=False, the matrix Q is represented as a product of elementary
71          reflectors
72
73             Q = H(0) H(2) . . . H(n-2).
74
75          Each H(i) has the form
76
77             H(i) = I - tau * v * v'
78
79          where tau is a complex scalar, and v is a complex vector with
80          v(0:i) = 0, v(i+1) = 1, v(i+2:n-1) is stored on exit in A(i+2:n-1,i),
81          and tau in TAU(i).
82
83          The contents of A on exit are illustrated by the following examples
84          with n = 5:
85
86          if UPLO = 'U':                       if UPLO = 'L':
87
88            (  d   e   v1  v2  v3 )              (  d                  )
89            (      d   e   v2  v3 )              (  e   d              )
90            (          d   e   v3 )              (  v0  e   d          )
91            (              d   e  )              (  v0  v1  e   d      )
92            (                  d  )              (  v0  v1  v2  e   d  )
93
94        where d and e denote diagonal and off-diagonal elements of T, and vi
95        denotes an element of the vector defining H(i).
96
97          -- LAPACK routine (version 3.0) --
98             Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
99             Courant Institute, Argonne National Lab, and Rice University
100             October 31, 1992
101        *************************************************************************/
102        public static void hmatrixtd(ref AP.Complex[,] a,
103            int n,
104            bool isupper,
105            ref AP.Complex[] tau,
106            ref double[] d,
107            ref double[] e)
108        {
109            int i = 0;
110            AP.Complex alpha = 0;
111            AP.Complex taui = 0;
112            AP.Complex v = 0;
113            AP.Complex[] t = new AP.Complex[0];
114            AP.Complex[] t2 = new AP.Complex[0];
115            AP.Complex[] t3 = new AP.Complex[0];
116            int i_ = 0;
117            int i1_ = 0;
118
119            if( n<=0 )
120            {
121                return;
122            }
123            for(i=0; i<=n-1; i++)
124            {
125                System.Diagnostics.Debug.Assert((double)(a[i,i].y)==(double)(0));
126            }
127            if( n>1 )
128            {
129                tau = new AP.Complex[n-2+1];
130                e = new double[n-2+1];
131            }
132            d = new double[n-1+1];
133            t = new AP.Complex[n-1+1];
134            t2 = new AP.Complex[n-1+1];
135            t3 = new AP.Complex[n-1+1];
136            if( isupper )
137            {
138               
139                //
140                // Reduce the upper triangle of A
141                //
142                a[n-1,n-1] = a[n-1,n-1].x;
143                for(i=n-2; i>=0; i--)
144                {
145                   
146                    //
147                    // Generate elementary reflector H = I+1 - tau * v * v'
148                    //
149                    alpha = a[i,i+1];
150                    t[1] = alpha;
151                    if( i>=1 )
152                    {
153                        i1_ = (0) - (2);
154                        for(i_=2; i_<=i+1;i_++)
155                        {
156                            t[i_] = a[i_+i1_,i+1];
157                        }
158                    }
159                    creflections.complexgeneratereflection(ref t, i+1, ref taui);
160                    if( i>=1 )
161                    {
162                        i1_ = (2) - (0);
163                        for(i_=0; i_<=i-1;i_++)
164                        {
165                            a[i_,i+1] = t[i_+i1_];
166                        }
167                    }
168                    alpha = t[1];
169                    e[i] = alpha.x;
170                    if( taui!=0 )
171                    {
172                       
173                        //
174                        // Apply H(I+1) from both sides to A
175                        //
176                        a[i,i+1] = 1;
177                       
178                        //
179                        // Compute  x := tau * A * v  storing x in TAU
180                        //
181                        i1_ = (0) - (1);
182                        for(i_=1; i_<=i+1;i_++)
183                        {
184                            t[i_] = a[i_+i1_,i+1];
185                        }
186                        hblas.hermitianmatrixvectormultiply(ref a, isupper, 0, i, ref t, taui, ref t2);
187                        i1_ = (1) - (0);
188                        for(i_=0; i_<=i;i_++)
189                        {
190                            tau[i_] = t2[i_+i1_];
191                        }
192                       
193                        //
194                        // Compute  w := x - 1/2 * tau * (x'*v) * v
195                        //
196                        v = 0.0;
197                        for(i_=0; i_<=i;i_++)
198                        {
199                            v += AP.Math.Conj(tau[i_])*a[i_,i+1];
200                        }
201                        alpha = -(0.5*taui*v);
202                        for(i_=0; i_<=i;i_++)
203                        {
204                            tau[i_] = tau[i_] + alpha*a[i_,i+1];
205                        }
206                       
207                        //
208                        // Apply the transformation as a rank-2 update:
209                        //    A := A - v * w' - w * v'
210                        //
211                        i1_ = (0) - (1);
212                        for(i_=1; i_<=i+1;i_++)
213                        {
214                            t[i_] = a[i_+i1_,i+1];
215                        }
216                        i1_ = (0) - (1);
217                        for(i_=1; i_<=i+1;i_++)
218                        {
219                            t3[i_] = tau[i_+i1_];
220                        }
221                        hblas.hermitianrank2update(ref a, isupper, 0, i, ref t, ref t3, ref t2, -1);
222                    }
223                    else
224                    {
225                        a[i,i] = a[i,i].x;
226                    }
227                    a[i,i+1] = e[i];
228                    d[i+1] = a[i+1,i+1].x;
229                    tau[i] = taui;
230                }
231                d[0] = a[0,0].x;
232            }
233            else
234            {
235               
236                //
237                // Reduce the lower triangle of A
238                //
239                a[0,0] = a[0,0].x;
240                for(i=0; i<=n-2; i++)
241                {
242                   
243                    //
244                    // Generate elementary reflector H = I - tau * v * v'
245                    //
246                    i1_ = (i+1) - (1);
247                    for(i_=1; i_<=n-i-1;i_++)
248                    {
249                        t[i_] = a[i_+i1_,i];
250                    }
251                    creflections.complexgeneratereflection(ref t, n-i-1, ref taui);
252                    i1_ = (1) - (i+1);
253                    for(i_=i+1; i_<=n-1;i_++)
254                    {
255                        a[i_,i] = t[i_+i1_];
256                    }
257                    e[i] = a[i+1,i].x;
258                    if( taui!=0 )
259                    {
260                       
261                        //
262                        // Apply H(i) from both sides to A(i+1:n,i+1:n)
263                        //
264                        a[i+1,i] = 1;
265                       
266                        //
267                        // Compute  x := tau * A * v  storing y in TAU
268                        //
269                        i1_ = (i+1) - (1);
270                        for(i_=1; i_<=n-i-1;i_++)
271                        {
272                            t[i_] = a[i_+i1_,i];
273                        }
274                        hblas.hermitianmatrixvectormultiply(ref a, isupper, i+1, n-1, ref t, taui, ref t2);
275                        i1_ = (1) - (i);
276                        for(i_=i; i_<=n-2;i_++)
277                        {
278                            tau[i_] = t2[i_+i1_];
279                        }
280                       
281                        //
282                        // Compute  w := x - 1/2 * tau * (x'*v) * v
283                        //
284                        i1_ = (i+1)-(i);
285                        v = 0.0;
286                        for(i_=i; i_<=n-2;i_++)
287                        {
288                            v += AP.Math.Conj(tau[i_])*a[i_+i1_,i];
289                        }
290                        alpha = -(0.5*taui*v);
291                        i1_ = (i+1) - (i);
292                        for(i_=i; i_<=n-2;i_++)
293                        {
294                            tau[i_] = tau[i_] + alpha*a[i_+i1_,i];
295                        }
296                       
297                        //
298                        // Apply the transformation as a rank-2 update:
299                        // A := A - v * w' - w * v'
300                        //
301                        i1_ = (i+1) - (1);
302                        for(i_=1; i_<=n-i-1;i_++)
303                        {
304                            t[i_] = a[i_+i1_,i];
305                        }
306                        i1_ = (i) - (1);
307                        for(i_=1; i_<=n-i-1;i_++)
308                        {
309                            t2[i_] = tau[i_+i1_];
310                        }
311                        hblas.hermitianrank2update(ref a, isupper, i+1, n-1, ref t, ref t2, ref t3, -1);
312                    }
313                    else
314                    {
315                        a[i+1,i+1] = a[i+1,i+1].x;
316                    }
317                    a[i+1,i] = e[i];
318                    d[i] = a[i,i].x;
319                    tau[i] = taui;
320                }
321                d[n-1] = a[n-1,n-1].x;
322            }
323        }
324
325
326        /*************************************************************************
327        Unpacking matrix Q which reduces a Hermitian matrix to a real  tridiagonal
328        form.
329
330        Input parameters:
331            A       -   the result of a HMatrixTD subroutine
332            N       -   size of matrix A.
333            IsUpper -   storage format (a parameter of HMatrixTD subroutine)
334            Tau     -   the result of a HMatrixTD subroutine
335
336        Output parameters:
337            Q       -   transformation matrix.
338                        array with elements [0..N-1, 0..N-1].
339
340          -- ALGLIB --
341             Copyright 2005, 2007, 2008 by Bochkanov Sergey
342        *************************************************************************/
343        public static void hmatrixtdunpackq(ref AP.Complex[,] a,
344            int n,
345            bool isupper,
346            ref AP.Complex[] tau,
347            ref AP.Complex[,] q)
348        {
349            int i = 0;
350            int j = 0;
351            AP.Complex[] v = new AP.Complex[0];
352            AP.Complex[] work = new AP.Complex[0];
353            int i_ = 0;
354            int i1_ = 0;
355
356            if( n==0 )
357            {
358                return;
359            }
360           
361            //
362            // init
363            //
364            q = new AP.Complex[n-1+1, n-1+1];
365            v = new AP.Complex[n+1];
366            work = new AP.Complex[n-1+1];
367            for(i=0; i<=n-1; i++)
368            {
369                for(j=0; j<=n-1; j++)
370                {
371                    if( i==j )
372                    {
373                        q[i,j] = 1;
374                    }
375                    else
376                    {
377                        q[i,j] = 0;
378                    }
379                }
380            }
381           
382            //
383            // unpack Q
384            //
385            if( isupper )
386            {
387                for(i=0; i<=n-2; i++)
388                {
389                   
390                    //
391                    // Apply H(i)
392                    //
393                    i1_ = (0) - (1);
394                    for(i_=1; i_<=i+1;i_++)
395                    {
396                        v[i_] = a[i_+i1_,i+1];
397                    }
398                    v[i+1] = 1;
399                    creflections.complexapplyreflectionfromtheleft(ref q, tau[i], ref v, 0, i, 0, n-1, ref work);
400                }
401            }
402            else
403            {
404                for(i=n-2; i>=0; i--)
405                {
406                   
407                    //
408                    // Apply H(i)
409                    //
410                    i1_ = (i+1) - (1);
411                    for(i_=1; i_<=n-i-1;i_++)
412                    {
413                        v[i_] = a[i_+i1_,i];
414                    }
415                    v[1] = 1;
416                    creflections.complexapplyreflectionfromtheleft(ref q, tau[i], ref v, i+1, n-1, 0, n-1, ref work);
417                }
418            }
419        }
420
421
422        /*************************************************************************
423
424          -- LAPACK routine (version 3.0) --
425             Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
426             Courant Institute, Argonne National Lab, and Rice University
427             October 31, 1992
428        *************************************************************************/
429        public static void hermitiantotridiagonal(ref AP.Complex[,] a,
430            int n,
431            bool isupper,
432            ref AP.Complex[] tau,
433            ref double[] d,
434            ref double[] e)
435        {
436            int i = 0;
437            AP.Complex alpha = 0;
438            AP.Complex taui = 0;
439            AP.Complex v = 0;
440            AP.Complex[] t = new AP.Complex[0];
441            AP.Complex[] t2 = new AP.Complex[0];
442            AP.Complex[] t3 = new AP.Complex[0];
443            int i_ = 0;
444            int i1_ = 0;
445
446            if( n<=0 )
447            {
448                return;
449            }
450            for(i=1; i<=n; i++)
451            {
452                System.Diagnostics.Debug.Assert((double)(a[i,i].y)==(double)(0));
453            }
454            tau = new AP.Complex[Math.Max(1, n-1)+1];
455            d = new double[n+1];
456            e = new double[Math.Max(1, n-1)+1];
457            t = new AP.Complex[n+1];
458            t2 = new AP.Complex[n+1];
459            t3 = new AP.Complex[n+1];
460            if( isupper )
461            {
462               
463                //
464                // Reduce the upper triangle of A
465                //
466                a[n,n] = a[n,n].x;
467                for(i=n-1; i>=1; i--)
468                {
469                   
470                    //
471                    // Generate elementary reflector H(i) = I - tau * v * v'
472                    // to annihilate A(1:i-1,i+1)
473                    //
474                    alpha = a[i,i+1];
475                    t[1] = alpha;
476                    if( i>=2 )
477                    {
478                        i1_ = (1) - (2);
479                        for(i_=2; i_<=i;i_++)
480                        {
481                            t[i_] = a[i_+i1_,i+1];
482                        }
483                    }
484                    creflections.complexgeneratereflection(ref t, i, ref taui);
485                    if( i>=2 )
486                    {
487                        i1_ = (2) - (1);
488                        for(i_=1; i_<=i-1;i_++)
489                        {
490                            a[i_,i+1] = t[i_+i1_];
491                        }
492                    }
493                    alpha = t[1];
494                    e[i] = alpha.x;
495                    if( taui!=0 )
496                    {
497                       
498                        //
499                        // Apply H(i) from both sides to A(1:i,1:i)
500                        //
501                        a[i,i+1] = 1;
502                       
503                        //
504                        // Compute  x := tau * A * v  storing x in TAU(1:i)
505                        //
506                        for(i_=1; i_<=i;i_++)
507                        {
508                            t[i_] = a[i_,i+1];
509                        }
510                        hblas.hermitianmatrixvectormultiply(ref a, isupper, 1, i, ref t, taui, ref tau);
511                       
512                        //
513                        // Compute  w := x - 1/2 * tau * (x'*v) * v
514                        //
515                        v = 0.0;
516                        for(i_=1; i_<=i;i_++)
517                        {
518                            v += AP.Math.Conj(tau[i_])*a[i_,i+1];
519                        }
520                        alpha = -(0.5*taui*v);
521                        for(i_=1; i_<=i;i_++)
522                        {
523                            tau[i_] = tau[i_] + alpha*a[i_,i+1];
524                        }
525                       
526                        //
527                        // Apply the transformation as a rank-2 update:
528                        //    A := A - v * w' - w * v'
529                        //
530                        for(i_=1; i_<=i;i_++)
531                        {
532                            t[i_] = a[i_,i+1];
533                        }
534                        hblas.hermitianrank2update(ref a, isupper, 1, i, ref t, ref tau, ref t2, -1);
535                    }
536                    else
537                    {
538                        a[i,i] = a[i,i].x;
539                    }
540                    a[i,i+1] = e[i];
541                    d[i+1] = a[i+1,i+1].x;
542                    tau[i] = taui;
543                }
544                d[1] = a[1,1].x;
545            }
546            else
547            {
548               
549                //
550                // Reduce the lower triangle of A
551                //
552                a[1,1] = a[1,1].x;
553                for(i=1; i<=n-1; i++)
554                {
555                   
556                    //
557                    // Generate elementary reflector H(i) = I - tau * v * v'
558                    // to annihilate A(i+2:n,i)
559                    //
560                    i1_ = (i+1) - (1);
561                    for(i_=1; i_<=n-i;i_++)
562                    {
563                        t[i_] = a[i_+i1_,i];
564                    }
565                    creflections.complexgeneratereflection(ref t, n-i, ref taui);
566                    i1_ = (1) - (i+1);
567                    for(i_=i+1; i_<=n;i_++)
568                    {
569                        a[i_,i] = t[i_+i1_];
570                    }
571                    e[i] = a[i+1,i].x;
572                    if( taui!=0 )
573                    {
574                       
575                        //
576                        // Apply H(i) from both sides to A(i+1:n,i+1:n)
577                        //
578                        a[i+1,i] = 1;
579                       
580                        //
581                        // Compute  x := tau * A * v  storing y in TAU(i:n-1)
582                        //
583                        i1_ = (i+1) - (1);
584                        for(i_=1; i_<=n-i;i_++)
585                        {
586                            t[i_] = a[i_+i1_,i];
587                        }
588                        hblas.hermitianmatrixvectormultiply(ref a, isupper, i+1, n, ref t, taui, ref t2);
589                        i1_ = (1) - (i);
590                        for(i_=i; i_<=n-1;i_++)
591                        {
592                            tau[i_] = t2[i_+i1_];
593                        }
594                       
595                        //
596                        // Compute  w := x - 1/2 * tau * (x'*v) * v
597                        //
598                        i1_ = (i+1)-(i);
599                        v = 0.0;
600                        for(i_=i; i_<=n-1;i_++)
601                        {
602                            v += AP.Math.Conj(tau[i_])*a[i_+i1_,i];
603                        }
604                        alpha = -(0.5*taui*v);
605                        i1_ = (i+1) - (i);
606                        for(i_=i; i_<=n-1;i_++)
607                        {
608                            tau[i_] = tau[i_] + alpha*a[i_+i1_,i];
609                        }
610                       
611                        //
612                        // Apply the transformation as a rank-2 update:
613                        // A := A - v * w' - w * v'
614                        //
615                        i1_ = (i+1) - (1);
616                        for(i_=1; i_<=n-i;i_++)
617                        {
618                            t[i_] = a[i_+i1_,i];
619                        }
620                        i1_ = (i) - (1);
621                        for(i_=1; i_<=n-i;i_++)
622                        {
623                            t2[i_] = tau[i_+i1_];
624                        }
625                        hblas.hermitianrank2update(ref a, isupper, i+1, n, ref t, ref t2, ref t3, -1);
626                    }
627                    else
628                    {
629                        a[i+1,i+1] = a[i+1,i+1].x;
630                    }
631                    a[i+1,i] = e[i];
632                    d[i] = a[i,i].x;
633                    tau[i] = taui;
634                }
635                d[n] = a[n,n].x;
636            }
637        }
638
639
640        /*************************************************************************
641
642          -- ALGLIB --
643             Copyright 2005, 2007 by Bochkanov Sergey
644        *************************************************************************/
645        public static void unpackqfromhermitiantridiagonal(ref AP.Complex[,] a,
646            int n,
647            bool isupper,
648            ref AP.Complex[] tau,
649            ref AP.Complex[,] q)
650        {
651            int i = 0;
652            int j = 0;
653            AP.Complex[] v = new AP.Complex[0];
654            AP.Complex[] work = new AP.Complex[0];
655            int i_ = 0;
656            int i1_ = 0;
657
658            if( n==0 )
659            {
660                return;
661            }
662           
663            //
664            // init
665            //
666            q = new AP.Complex[n+1, n+1];
667            v = new AP.Complex[n+1];
668            work = new AP.Complex[n+1];
669            for(i=1; i<=n; i++)
670            {
671                for(j=1; j<=n; j++)
672                {
673                    if( i==j )
674                    {
675                        q[i,j] = 1;
676                    }
677                    else
678                    {
679                        q[i,j] = 0;
680                    }
681                }
682            }
683           
684            //
685            // unpack Q
686            //
687            if( isupper )
688            {
689                for(i=1; i<=n-1; i++)
690                {
691                   
692                    //
693                    // Apply H(i)
694                    //
695                    for(i_=1; i_<=i;i_++)
696                    {
697                        v[i_] = a[i_,i+1];
698                    }
699                    v[i] = 1;
700                    creflections.complexapplyreflectionfromtheleft(ref q, tau[i], ref v, 1, i, 1, n, ref work);
701                }
702            }
703            else
704            {
705                for(i=n-1; i>=1; i--)
706                {
707                   
708                    //
709                    // Apply H(i)
710                    //
711                    i1_ = (i+1) - (1);
712                    for(i_=1; i_<=n-i;i_++)
713                    {
714                        v[i_] = a[i_+i1_,i];
715                    }
716                    v[1] = 1;
717                    creflections.complexapplyreflectionfromtheleft(ref q, tau[i], ref v, i+1, n, 1, n, ref work);
718                }
719            }
720        }
721    }
722}
Note: See TracBrowser for help on using the repository browser.