1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23 package org.hipparchus.linear;
24
25 import org.hipparchus.exception.LocalizedCoreFormats;
26 import org.hipparchus.exception.MathIllegalArgumentException;
27 import org.hipparchus.util.FastMath;
28 import org.hipparchus.util.Precision;
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47 public class HessenbergTransformer {
48
49 private final double[][] householderVectors;
50
51 private final double[] ort;
52
53 private RealMatrix cachedP;
54
55 private RealMatrix cachedPt;
56
57 private RealMatrix cachedH;
58
59
60
61
62
63
64
65 public HessenbergTransformer(final RealMatrix matrix) {
66 if (!matrix.isSquare()) {
67 throw new MathIllegalArgumentException(LocalizedCoreFormats.NON_SQUARE_MATRIX,
68 matrix.getRowDimension(), matrix.getColumnDimension());
69 }
70
71 final int m = matrix.getRowDimension();
72 householderVectors = matrix.getData();
73 ort = new double[m];
74 cachedP = null;
75 cachedPt = null;
76 cachedH = null;
77
78
79 transform();
80 }
81
82
83
84
85
86
87
88 public RealMatrix getP() {
89 if (cachedP == null) {
90 final int n = householderVectors.length;
91 final int high = n - 1;
92 final double[][] pa = new double[n][n];
93
94 for (int i = 0; i < n; i++) {
95 for (int j = 0; j < n; j++) {
96 pa[i][j] = (i == j) ? 1 : 0;
97 }
98 }
99
100 for (int m = high - 1; m >= 1; m--) {
101 if (householderVectors[m][m - 1] != 0.0) {
102 for (int i = m + 1; i <= high; i++) {
103 ort[i] = householderVectors[i][m - 1];
104 }
105
106 for (int j = m; j <= high; j++) {
107 double g = 0.0;
108
109 for (int i = m; i <= high; i++) {
110 g += ort[i] * pa[i][j];
111 }
112
113
114 g = (g / ort[m]) / householderVectors[m][m - 1];
115
116 for (int i = m; i <= high; i++) {
117 pa[i][j] += g * ort[i];
118 }
119 }
120 }
121 }
122
123 cachedP = MatrixUtils.createRealMatrix(pa);
124 }
125 return cachedP;
126 }
127
128
129
130
131
132
133
134 public RealMatrix getPT() {
135 if (cachedPt == null) {
136 cachedPt = getP().transpose();
137 }
138
139
140 return cachedPt;
141 }
142
143
144
145
146
147
148 public RealMatrix getH() {
149 if (cachedH == null) {
150 final int m = householderVectors.length;
151 final double[][] h = new double[m][m];
152 for (int i = 0; i < m; ++i) {
153 if (i > 0) {
154
155 h[i][i - 1] = householderVectors[i][i - 1];
156 }
157
158
159 System.arraycopy(householderVectors[i], i, h[i], i, m - i);
160 }
161 cachedH = MatrixUtils.createRealMatrix(h);
162 }
163
164
165 return cachedH;
166 }
167
168
169
170
171
172
173
174
175 double[][] getHouseholderVectorsRef() {
176 return householderVectors;
177 }
178
179
180
181
182
183 private void transform() {
184 final int n = householderVectors.length;
185 final int high = n - 1;
186
187 for (int m = 1; m <= high - 1; m++) {
188
189 double scale = 0;
190 for (int i = m; i <= high; i++) {
191 scale += FastMath.abs(householderVectors[i][m - 1]);
192 }
193
194 if (!Precision.equals(scale, 0)) {
195
196 double h = 0;
197 for (int i = high; i >= m; i--) {
198 ort[i] = householderVectors[i][m - 1] / scale;
199 h += ort[i] * ort[i];
200 }
201 final double g = (ort[m] > 0) ? -FastMath.sqrt(h) : FastMath.sqrt(h);
202
203 h -= ort[m] * g;
204 ort[m] -= g;
205
206
207
208
209 for (int j = m; j < n; j++) {
210 double f = 0;
211 for (int i = high; i >= m; i--) {
212 f += ort[i] * householderVectors[i][j];
213 }
214 f /= h;
215 for (int i = m; i <= high; i++) {
216 householderVectors[i][j] -= f * ort[i];
217 }
218 }
219
220 for (int i = 0; i <= high; i++) {
221 double f = 0;
222 for (int j = high; j >= m; j--) {
223 f += ort[j] * householderVectors[i][j];
224 }
225 f /= h;
226 for (int j = m; j <= high; j++) {
227 householderVectors[i][j] -= f * ort[j];
228 }
229 }
230
231 ort[m] = scale * ort[m];
232 householderVectors[m][m - 1] = scale * g;
233 }
234 }
235 }
236 }