1 | /*************************************************************************
|
---|
2 | Copyright (c) 1992-2007 The University of Tennessee. All rights reserved.
|
---|
3 |
|
---|
4 | Contributors:
|
---|
5 | * Sergey Bochkanov (ALGLIB project). Translation from FORTRAN to
|
---|
6 | pseudocode.
|
---|
7 |
|
---|
8 | See subroutines comments for additional copyrights.
|
---|
9 |
|
---|
10 | >>> SOURCE LICENSE >>>
|
---|
11 | This program is free software; you can redistribute it and/or modify
|
---|
12 | it under the terms of the GNU General Public License as published by
|
---|
13 | the Free Software Foundation (www.fsf.org); either version 2 of the
|
---|
14 | License, or (at your option) any later version.
|
---|
15 |
|
---|
16 | This program is distributed in the hope that it will be useful,
|
---|
17 | but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
18 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
19 | GNU General Public License for more details.
|
---|
20 |
|
---|
21 | A copy of the GNU General Public License is available at
|
---|
22 | http://www.fsf.org/licensing/licenses
|
---|
23 |
|
---|
24 | >>> END OF LICENSE >>>
|
---|
25 | *************************************************************************/
|
---|
26 |
|
---|
27 | using System;
|
---|
28 |
|
---|
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.
|
---|
36 |
|
---|
37 | The source matrix A is represented as S'*A*S = T, where S is an orthogonal
|
---|
38 | matrix (Schur vectors), T - upper quasi-triangular matrix (with blocks of
|
---|
39 | sizes 1x1 and 2x2 on the main diagonal).
|
---|
40 |
|
---|
41 | Input parameters:
|
---|
42 | A - matrix to be decomposed.
|
---|
43 | Array whose indexes range within [0..N-1, 0..N-1].
|
---|
44 | N - size of A, N>=0.
|
---|
45 |
|
---|
46 |
|
---|
47 | Output parameters:
|
---|
48 | A - contains matrix T.
|
---|
49 | Array whose indexes range within [0..N-1, 0..N-1].
|
---|
50 | S - contains Schur vectors.
|
---|
51 | Array whose indexes range within [0..N-1, 0..N-1].
|
---|
52 |
|
---|
53 | Note 1:
|
---|
54 | The block structure of matrix T can be easily recognized: since all
|
---|
55 | the elements below the blocks are zeros, the elements a[i+1,i] which
|
---|
56 | are equal to 0 show the block border.
|
---|
57 |
|
---|
58 | Note 2:
|
---|
59 | The algorithm performance depends on the value of the internal parameter
|
---|
60 | NS of the InternalSchurDecomposition subroutine which defines the number
|
---|
61 | of shifts in the QR algorithm (similarly to the block width in block-matrix
|
---|
62 | algorithms in linear algebra). If you require maximum performance on
|
---|
63 | your machine, it is recommended to adjust this parameter manually.
|
---|
64 |
|
---|
65 | Result:
|
---|
66 | True,
|
---|
67 | if the algorithm has converged and parameters A and S contain the result.
|
---|
68 | False,
|
---|
69 | if the algorithm has not converged.
|
---|
70 |
|
---|
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;
|
---|
86 |
|
---|
87 |
|
---|
88 | //
|
---|
89 | // Upper Hessenberg form of the 0-based matrix
|
---|
90 | //
|
---|
91 | hessenberg.rmatrixhessenberg(ref a, n, ref tau);
|
---|
92 | hessenberg.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;
|
---|
126 | }
|
---|
127 |
|
---|
128 |
|
---|
129 | public static bool schurdecomposition(ref double[,] a,
|
---|
130 | int n,
|
---|
131 | ref double[,] s)
|
---|
132 | {
|
---|
133 | bool result = new bool();
|
---|
134 | double[] tau = new double[0];
|
---|
135 | double[] wi = new double[0];
|
---|
136 | double[] wr = new double[0];
|
---|
137 | int info = 0;
|
---|
138 |
|
---|
139 | hessenberg.toupperhessenberg(ref a, n, ref tau);
|
---|
140 | hessenberg.unpackqfromupperhessenberg(ref a, n, ref tau, ref s);
|
---|
141 | hsschur.internalschurdecomposition(ref a, n, 1, 1, ref wr, ref wi, ref s, ref info);
|
---|
142 | result = info==0;
|
---|
143 | return result;
|
---|
144 | }
|
---|
145 | }
|
---|
146 | }
|
---|