Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/ALGLIB/hevd.cs @ 2636

Last change on this file since 2636 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: 8.8 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 hevd
26    {
27        /*************************************************************************
28        Finding the eigenvalues and eigenvectors of a Hermitian matrix
29
30        The algorithm finds eigen pairs of a Hermitian matrix by  reducing  it  to
31        real tridiagonal form and using the QL/QR algorithm.
32
33        Input parameters:
34            A       -   Hermitian matrix which is given  by  its  upper  or  lower
35                        triangular part.
36                        Array whose indexes range within [0..N-1, 0..N-1].
37            N       -   size of matrix A.
38            IsUpper -   storage format.
39            ZNeeded -   flag controlling whether the eigenvectors  are  needed  or
40                        not. If ZNeeded is equal to:
41                         * 0, the eigenvectors are not returned;
42                         * 1, the eigenvectors are returned.
43
44        Output parameters:
45            D       -   eigenvalues in ascending order.
46                        Array whose index ranges within [0..N-1].
47            Z       -   if ZNeeded is equal to:
48                         * 0, Z hasn’t changed;
49                         * 1, Z contains the eigenvectors.
50                        Array whose indexes range within [0..N-1, 0..N-1].
51                        The eigenvectors are stored in the matrix columns.
52
53        Result:
54            True, if the algorithm has converged.
55            False, if the algorithm hasn't converged (rare case).
56
57        Note:
58            eigen vectors of Hermitian matrix are defined up to multiplication  by
59            a complex number L, such as |L|=1.
60
61          -- ALGLIB --
62             Copyright 2005, 23 March 2007 by Bochkanov Sergey
63        *************************************************************************/
64        public static bool hmatrixevd(AP.Complex[,] a,
65            int n,
66            int zneeded,
67            bool isupper,
68            ref double[] d,
69            ref AP.Complex[,] z)
70        {
71            bool result = new bool();
72            AP.Complex[] tau = new AP.Complex[0];
73            double[] e = new double[0];
74            double[] work = new double[0];
75            double[,] t = new double[0,0];
76            AP.Complex[,] q = new AP.Complex[0,0];
77            int i = 0;
78            int k = 0;
79            double v = 0;
80            int i_ = 0;
81
82            a = (AP.Complex[,])a.Clone();
83
84            System.Diagnostics.Debug.Assert(zneeded==0 | zneeded==1, "HermitianEVD: incorrect ZNeeded");
85           
86            //
87            // Reduce to tridiagonal form
88            //
89            htridiagonal.hmatrixtd(ref a, n, isupper, ref tau, ref d, ref e);
90            if( zneeded==1 )
91            {
92                htridiagonal.hmatrixtdunpackq(ref a, n, isupper, ref tau, ref q);
93                zneeded = 2;
94            }
95           
96            //
97            // TDEVD
98            //
99            result = tdevd.smatrixtdevd(ref d, e, n, zneeded, ref t);
100           
101            //
102            // Eigenvectors are needed
103            // Calculate Z = Q*T = Re(Q)*T + i*Im(Q)*T
104            //
105            if( result & zneeded!=0 )
106            {
107                work = new double[n-1+1];
108                z = new AP.Complex[n-1+1, n-1+1];
109                for(i=0; i<=n-1; i++)
110                {
111                   
112                    //
113                    // Calculate real part
114                    //
115                    for(k=0; k<=n-1; k++)
116                    {
117                        work[k] = 0;
118                    }
119                    for(k=0; k<=n-1; k++)
120                    {
121                        v = q[i,k].x;
122                        for(i_=0; i_<=n-1;i_++)
123                        {
124                            work[i_] = work[i_] + v*t[k,i_];
125                        }
126                    }
127                    for(k=0; k<=n-1; k++)
128                    {
129                        z[i,k].x = work[k];
130                    }
131                   
132                    //
133                    // Calculate imaginary part
134                    //
135                    for(k=0; k<=n-1; k++)
136                    {
137                        work[k] = 0;
138                    }
139                    for(k=0; k<=n-1; k++)
140                    {
141                        v = q[i,k].y;
142                        for(i_=0; i_<=n-1;i_++)
143                        {
144                            work[i_] = work[i_] + v*t[k,i_];
145                        }
146                    }
147                    for(k=0; k<=n-1; k++)
148                    {
149                        z[i,k].y = work[k];
150                    }
151                }
152            }
153            return result;
154        }
155
156
157        /*************************************************************************
158
159          -- ALGLIB --
160             Copyright 2005, 23 March 2007 by Bochkanov Sergey
161        *************************************************************************/
162        public static bool hermitianevd(AP.Complex[,] a,
163            int n,
164            int zneeded,
165            bool isupper,
166            ref double[] d,
167            ref AP.Complex[,] z)
168        {
169            bool result = new bool();
170            AP.Complex[] tau = new AP.Complex[0];
171            double[] e = new double[0];
172            double[] work = new double[0];
173            double[,] t = new double[0,0];
174            AP.Complex[,] q = new AP.Complex[0,0];
175            int i = 0;
176            int k = 0;
177            double v = 0;
178            int i_ = 0;
179
180            a = (AP.Complex[,])a.Clone();
181
182            System.Diagnostics.Debug.Assert(zneeded==0 | zneeded==1, "HermitianEVD: incorrect ZNeeded");
183           
184            //
185            // Reduce to tridiagonal form
186            //
187            htridiagonal.hermitiantotridiagonal(ref a, n, isupper, ref tau, ref d, ref e);
188            if( zneeded==1 )
189            {
190                htridiagonal.unpackqfromhermitiantridiagonal(ref a, n, isupper, ref tau, ref q);
191                zneeded = 2;
192            }
193           
194            //
195            // TDEVD
196            //
197            result = tdevd.tridiagonalevd(ref d, e, n, zneeded, ref t);
198           
199            //
200            // Eigenvectors are needed
201            // Calculate Z = Q*T = Re(Q)*T + i*Im(Q)*T
202            //
203            if( result & zneeded!=0 )
204            {
205                work = new double[n+1];
206                z = new AP.Complex[n+1, n+1];
207                for(i=1; i<=n; i++)
208                {
209                   
210                    //
211                    // Calculate real part
212                    //
213                    for(k=1; k<=n; k++)
214                    {
215                        work[k] = 0;
216                    }
217                    for(k=1; k<=n; k++)
218                    {
219                        v = q[i,k].x;
220                        for(i_=1; i_<=n;i_++)
221                        {
222                            work[i_] = work[i_] + v*t[k,i_];
223                        }
224                    }
225                    for(k=1; k<=n; k++)
226                    {
227                        z[i,k].x = work[k];
228                    }
229                   
230                    //
231                    // Calculate imaginary part
232                    //
233                    for(k=1; k<=n; k++)
234                    {
235                        work[k] = 0;
236                    }
237                    for(k=1; k<=n; k++)
238                    {
239                        v = q[i,k].y;
240                        for(i_=1; i_<=n;i_++)
241                        {
242                            work[i_] = work[i_] + v*t[k,i_];
243                        }
244                    }
245                    for(k=1; k<=n; k++)
246                    {
247                        z[i,k].y = work[k];
248                    }
249                }
250            }
251            return result;
252        }
253    }
254}
Note: See TracBrowser for help on using the repository browser.