Changeset 4068 for trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.5.0/ALGLIB-2.5.0/schur.cs
- Timestamp:
- 07/22/10 00:44:01 (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.5.0/ALGLIB-2.5.0/schur.cs
r3839 r4068 25 25 *************************************************************************/ 26 26 27 using System;28 27 29 namespace alglib 30 { 31 public class schur 32 { 33 /************************************************************************* 34 Subroutine performing the Schur decomposition of a general matrix by using 35 the QR algorithm with multiple shifts. 28 namespace alglib { 29 public class schur { 30 /************************************************************************* 31 Subroutine performing the Schur decomposition of a general matrix by using 32 the QR algorithm with multiple shifts. 36 33 37 38 39 34 The source matrix A is represented as S'*A*S = T, where S is an orthogonal 35 matrix (Schur vectors), T - upper quasi-triangular matrix (with blocks of 36 sizes 1x1 and 2x2 on the main diagonal). 40 37 41 42 43 44 38 Input parameters: 39 A - matrix to be decomposed. 40 Array whose indexes range within [0..N-1, 0..N-1]. 41 N - size of A, N>=0. 45 42 46 43 47 48 49 50 51 44 Output parameters: 45 A - contains matrix T. 46 Array whose indexes range within [0..N-1, 0..N-1]. 47 S - contains Schur vectors. 48 Array whose indexes range within [0..N-1, 0..N-1]. 52 49 53 54 55 56 50 Note 1: 51 The block structure of matrix T can be easily recognized: since all 52 the elements below the blocks are zeros, the elements a[i+1,i] which 53 are equal to 0 show the block border. 57 54 58 59 60 61 62 63 55 Note 2: 56 The algorithm performance depends on the value of the internal parameter 57 NS of the InternalSchurDecomposition subroutine which defines the number 58 of shifts in the QR algorithm (similarly to the block width in block-matrix 59 algorithms in linear algebra). If you require maximum performance on 60 your machine, it is recommended to adjust this parameter manually. 64 61 65 66 67 68 69 62 Result: 63 True, 64 if the algorithm has converged and parameters A and S contain the result. 65 False, 66 if the algorithm has not converged. 70 67 71 Algorithm implemented on the basis of the DHSEQR subroutine (LAPACK 3.0 library). 72 *************************************************************************/ 73 public static bool rmatrixschur(ref double[,] a, 74 int n, 75 ref double[,] s) 76 { 77 bool result = new bool(); 78 double[] tau = new double[0]; 79 double[] wi = new double[0]; 80 double[] wr = new double[0]; 81 double[,] a1 = new double[0,0]; 82 double[,] s1 = new double[0,0]; 83 int info = 0; 84 int i = 0; 85 int j = 0; 68 Algorithm implemented on the basis of the DHSEQR subroutine (LAPACK 3.0 library). 69 *************************************************************************/ 70 public static bool rmatrixschur(ref double[,] a, 71 int n, 72 ref double[,] s) { 73 bool result = new bool(); 74 double[] tau = new double[0]; 75 double[] wi = new double[0]; 76 double[] wr = new double[0]; 77 double[,] a1 = new double[0, 0]; 78 double[,] s1 = new double[0, 0]; 79 int info = 0; 80 int i = 0; 81 int j = 0; 86 82 87 88 // 89 // Upper Hessenberg form of the 0-based matrix 90 // 91 ortfac.rmatrixhessenberg(ref a, n, ref tau); 92 ortfac.rmatrixhessenbergunpackq(ref a, n, ref tau, ref s); 93 94 // 95 // Convert from 0-based arrays to 1-based, 96 // then call InternalSchurDecomposition 97 // Awkward, of course, but Schur decompisiton subroutine 98 // is too complex to fix it. 99 // 100 // 101 a1 = new double[n+1, n+1]; 102 s1 = new double[n+1, n+1]; 103 for(i=1; i<=n; i++) 104 { 105 for(j=1; j<=n; j++) 106 { 107 a1[i,j] = a[i-1,j-1]; 108 s1[i,j] = s[i-1,j-1]; 109 } 110 } 111 hsschur.internalschurdecomposition(ref a1, n, 1, 1, ref wr, ref wi, ref s1, ref info); 112 result = info==0; 113 114 // 115 // convert from 1-based arrays to -based 116 // 117 for(i=1; i<=n; i++) 118 { 119 for(j=1; j<=n; j++) 120 { 121 a[i-1,j-1] = a1[i,j]; 122 s[i-1,j-1] = s1[i,j]; 123 } 124 } 125 return result; 83 84 // 85 // Upper Hessenberg form of the 0-based matrix 86 // 87 ortfac.rmatrixhessenberg(ref a, n, ref tau); 88 ortfac.rmatrixhessenbergunpackq(ref a, n, ref tau, ref s); 89 90 // 91 // Convert from 0-based arrays to 1-based, 92 // then call InternalSchurDecomposition 93 // Awkward, of course, but Schur decompisiton subroutine 94 // is too complex to fix it. 95 // 96 // 97 a1 = new double[n + 1, n + 1]; 98 s1 = new double[n + 1, n + 1]; 99 for (i = 1; i <= n; i++) { 100 for (j = 1; j <= n; j++) { 101 a1[i, j] = a[i - 1, j - 1]; 102 s1[i, j] = s[i - 1, j - 1]; 126 103 } 104 } 105 hsschur.internalschurdecomposition(ref a1, n, 1, 1, ref wr, ref wi, ref s1, ref info); 106 result = info == 0; 107 108 // 109 // convert from 1-based arrays to -based 110 // 111 for (i = 1; i <= n; i++) { 112 for (j = 1; j <= n; j++) { 113 a[i - 1, j - 1] = a1[i, j]; 114 s[i - 1, j - 1] = s1[i, j]; 115 } 116 } 117 return result; 127 118 } 119 } 128 120 }
Note: See TracChangeset
for help on using the changeset viewer.
