[2806] | 1 | /*************************************************************************
|
---|
| 2 | Copyright (c) 2005-2007, Sergey Bochkanov (ALGLIB project).
|
---|
| 3 |
|
---|
| 4 | >>> SOURCE LICENSE >>>
|
---|
| 5 | This program is free software; you can redistribute it and/or modify
|
---|
| 6 | it under the terms of the GNU General Public License as published by
|
---|
| 7 | the Free Software Foundation (www.fsf.org); either version 2 of the
|
---|
| 8 | License, or (at your option) any later version.
|
---|
| 9 |
|
---|
| 10 | This program is distributed in the hope that it will be useful,
|
---|
| 11 | but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
| 12 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
| 13 | GNU General Public License for more details.
|
---|
| 14 |
|
---|
| 15 | A copy of the GNU General Public License is available at
|
---|
| 16 | http://www.fsf.org/licensing/licenses
|
---|
| 17 |
|
---|
| 18 | >>> END OF LICENSE >>>
|
---|
| 19 | *************************************************************************/
|
---|
| 20 |
|
---|
| 21 | using System;
|
---|
| 22 |
|
---|
| 23 | namespace 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 hasnt 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 | }
|
---|