1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.apache.commons.math.linear;
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36 class BiDiagonalTransformer {
37
38
39 private final double householderVectors[][];
40
41
42 private final double[] main;
43
44
45 private final double[] secondary;
46
47
48 private RealMatrix cachedU;
49
50
51 private RealMatrix cachedB;
52
53
54 private RealMatrix cachedV;
55
56
57
58
59
60 public BiDiagonalTransformer(RealMatrix matrix) {
61
62 final int m = matrix.getRowDimension();
63 final int n = matrix.getColumnDimension();
64 final int p = Math.min(m, n);
65 householderVectors = matrix.getData();
66 main = new double[p];
67 secondary = new double[p - 1];
68 cachedU = null;
69 cachedB = null;
70 cachedV = null;
71
72
73 if (m >= n) {
74 transformToUpperBiDiagonal();
75 } else {
76 transformToLowerBiDiagonal();
77 }
78
79 }
80
81
82
83
84
85
86 public RealMatrix getU() {
87
88 if (cachedU == null) {
89
90 final int m = householderVectors.length;
91 final int n = householderVectors[0].length;
92 final int p = main.length;
93 final int diagOffset = (m >= n) ? 0 : 1;
94 final double[] diagonal = (m >= n) ? main : secondary;
95 cachedU = MatrixUtils.createRealMatrix(m, m);
96
97
98 for (int k = m - 1; k >= p; --k) {
99 cachedU.setEntry(k, k, 1);
100 }
101
102
103 for (int k = p - 1; k >= diagOffset; --k) {
104 final double[] hK = householderVectors[k];
105 cachedU.setEntry(k, k, 1);
106 if (hK[k - diagOffset] != 0.0) {
107 for (int j = k; j < m; ++j) {
108 double alpha = 0;
109 for (int i = k; i < m; ++i) {
110 alpha -= cachedU.getEntry(i, j) * householderVectors[i][k - diagOffset];
111 }
112 alpha /= diagonal[k - diagOffset] * hK[k - diagOffset];
113
114 for (int i = k; i < m; ++i) {
115 cachedU.addToEntry(i, j, -alpha * householderVectors[i][k - diagOffset]);
116 }
117 }
118 }
119 }
120 if (diagOffset > 0) {
121 cachedU.setEntry(0, 0, 1);
122 }
123
124 }
125
126
127 return cachedU;
128
129 }
130
131
132
133
134
135 public RealMatrix getB() {
136
137 if (cachedB == null) {
138
139 final int m = householderVectors.length;
140 final int n = householderVectors[0].length;
141 cachedB = MatrixUtils.createRealMatrix(m, n);
142 for (int i = 0; i < main.length; ++i) {
143 cachedB.setEntry(i, i, main[i]);
144 if (m < n) {
145 if (i > 0) {
146 cachedB.setEntry(i, i - 1, secondary[i - 1]);
147 }
148 } else {
149 if (i < main.length - 1) {
150 cachedB.setEntry(i, i + 1, secondary[i]);
151 }
152 }
153 }
154
155 }
156
157
158 return cachedB;
159
160 }
161
162
163
164
165
166
167 public RealMatrix getV() {
168
169 if (cachedV == null) {
170
171 final int m = householderVectors.length;
172 final int n = householderVectors[0].length;
173 final int p = main.length;
174 final int diagOffset = (m >= n) ? 1 : 0;
175 final double[] diagonal = (m >= n) ? secondary : main;
176 cachedV = MatrixUtils.createRealMatrix(n, n);
177
178
179 for (int k = n - 1; k >= p; --k) {
180 cachedV.setEntry(k, k, 1);
181 }
182
183
184 for (int k = p - 1; k >= diagOffset; --k) {
185 final double[] hK = householderVectors[k - diagOffset];
186 cachedV.setEntry(k, k, 1);
187 if (hK[k] != 0.0) {
188 for (int j = k; j < n; ++j) {
189 double beta = 0;
190 for (int i = k; i < n; ++i) {
191 beta -= cachedV.getEntry(i, j) * hK[i];
192 }
193 beta /= diagonal[k - diagOffset] * hK[k];
194
195 for (int i = k; i < n; ++i) {
196 cachedV.addToEntry(i, j, -beta * hK[i]);
197 }
198 }
199 }
200 }
201 if (diagOffset > 0) {
202 cachedV.setEntry(0, 0, 1);
203 }
204
205 }
206
207
208 return cachedV;
209
210 }
211
212
213
214
215
216
217
218 double[][] getHouseholderVectorsRef() {
219 return householderVectors;
220 }
221
222
223
224
225
226
227
228 double[] getMainDiagonalRef() {
229 return main;
230 }
231
232
233
234
235
236
237
238 double[] getSecondaryDiagonalRef() {
239 return secondary;
240 }
241
242
243
244
245
246 boolean isUpperBiDiagonal() {
247 return householderVectors.length >= householderVectors[0].length;
248 }
249
250
251
252
253
254
255 private void transformToUpperBiDiagonal() {
256
257 final int m = householderVectors.length;
258 final int n = householderVectors[0].length;
259 for (int k = 0; k < n; k++) {
260
261
262 double xNormSqr = 0;
263 for (int i = k; i < m; ++i) {
264 final double c = householderVectors[i][k];
265 xNormSqr += c * c;
266 }
267 final double[] hK = householderVectors[k];
268 final double a = (hK[k] > 0) ? -Math.sqrt(xNormSqr) : Math.sqrt(xNormSqr);
269 main[k] = a;
270 if (a != 0.0) {
271 hK[k] -= a;
272 for (int j = k + 1; j < n; ++j) {
273 double alpha = 0;
274 for (int i = k; i < m; ++i) {
275 final double[] hI = householderVectors[i];
276 alpha -= hI[j] * hI[k];
277 }
278 alpha /= a * householderVectors[k][k];
279 for (int i = k; i < m; ++i) {
280 final double[] hI = householderVectors[i];
281 hI[j] -= alpha * hI[k];
282 }
283 }
284 }
285
286 if (k < n - 1) {
287
288 xNormSqr = 0;
289 for (int j = k + 1; j < n; ++j) {
290 final double c = hK[j];
291 xNormSqr += c * c;
292 }
293 final double b = (hK[k + 1] > 0) ? -Math.sqrt(xNormSqr) : Math.sqrt(xNormSqr);
294 secondary[k] = b;
295 if (b != 0.0) {
296 hK[k + 1] -= b;
297 for (int i = k + 1; i < m; ++i) {
298 final double[] hI = householderVectors[i];
299 double beta = 0;
300 for (int j = k + 1; j < n; ++j) {
301 beta -= hI[j] * hK[j];
302 }
303 beta /= b * hK[k + 1];
304 for (int j = k + 1; j < n; ++j) {
305 hI[j] -= beta * hK[j];
306 }
307 }
308 }
309 }
310
311 }
312 }
313
314
315
316
317
318
319 private void transformToLowerBiDiagonal() {
320
321 final int m = householderVectors.length;
322 final int n = householderVectors[0].length;
323 for (int k = 0; k < m; k++) {
324
325
326 final double[] hK = householderVectors[k];
327 double xNormSqr = 0;
328 for (int j = k; j < n; ++j) {
329 final double c = hK[j];
330 xNormSqr += c * c;
331 }
332 final double a = (hK[k] > 0) ? -Math.sqrt(xNormSqr) : Math.sqrt(xNormSqr);
333 main[k] = a;
334 if (a != 0.0) {
335 hK[k] -= a;
336 for (int i = k + 1; i < m; ++i) {
337 final double[] hI = householderVectors[i];
338 double alpha = 0;
339 for (int j = k; j < n; ++j) {
340 alpha -= hI[j] * hK[j];
341 }
342 alpha /= a * householderVectors[k][k];
343 for (int j = k; j < n; ++j) {
344 hI[j] -= alpha * hK[j];
345 }
346 }
347 }
348
349 if (k < m - 1) {
350
351 final double[] hKp1 = householderVectors[k + 1];
352 xNormSqr = 0;
353 for (int i = k + 1; i < m; ++i) {
354 final double c = householderVectors[i][k];
355 xNormSqr += c * c;
356 }
357 final double b = (hKp1[k] > 0) ? -Math.sqrt(xNormSqr) : Math.sqrt(xNormSqr);
358 secondary[k] = b;
359 if (b != 0.0) {
360 hKp1[k] -= b;
361 for (int j = k + 1; j < n; ++j) {
362 double beta = 0;
363 for (int i = k + 1; i < m; ++i) {
364 final double[] hI = householderVectors[i];
365 beta -= hI[j] * hI[k];
366 }
367 beta /= b * hKp1[k];
368 for (int i = k + 1; i < m; ++i) {
369 final double[] hI = householderVectors[i];
370 hI[j] -= beta * hI[k];
371 }
372 }
373 }
374 }
375
376 }
377 }
378
379 }