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 import java.util.ArrayList;
21 import java.util.Arrays;
22 import java.util.List;
23
24 import org.apache.commons.math.ConvergenceException;
25 import org.apache.commons.math.MathRuntimeException;
26 import org.apache.commons.math.MaxIterationsExceededException;
27 import org.apache.commons.math.util.MathUtils;
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58 public class EigenDecompositionImpl implements EigenDecomposition {
59
60
61 private static final double TOLERANCE = 100 * MathUtils.EPSILON;
62
63
64 private static final double TOLERANCE_2 = TOLERANCE * TOLERANCE;
65
66
67 private double splitTolerance;
68
69
70 private double[] main;
71
72
73 private double[] secondary;
74
75
76 private double[] squaredSecondary;
77
78
79 private TriDiagonalTransformer transformer;
80
81
82 private double lowerSpectra;
83
84
85 private double upperSpectra;
86
87
88 private double minPivot;
89
90
91 private double sigma;
92
93
94 private double sigmaLow;
95
96
97 private double tau;
98
99
100 private double[] work;
101
102
103 private int pingPong;
104
105
106 private double qMax;
107
108
109 private double eMin;
110
111
112 private int tType;
113
114
115 private double dMin;
116
117
118 private double dMin1;
119
120
121 private double dMin2;
122
123
124 private double dN;
125
126
127 private double dN1;
128
129
130 private double dN2;
131
132
133 private double g;
134
135
136 private double[] realEigenvalues;
137
138
139 private double[] imagEigenvalues;
140
141
142 private ArrayRealVector[] eigenvectors;
143
144
145 private RealMatrix cachedV;
146
147
148 private RealMatrix cachedD;
149
150
151 private RealMatrix cachedVt;
152
153
154
155
156
157
158
159
160
161
162 public EigenDecompositionImpl(final RealMatrix matrix,
163 final double splitTolerance)
164 throws InvalidMatrixException {
165 if (isSymmetric(matrix)) {
166 this.splitTolerance = splitTolerance;
167 transformToTridiagonal(matrix);
168 decompose();
169 } else {
170
171
172 throw new InvalidMatrixException("eigen decomposition of assymetric matrices not supported yet");
173 }
174 }
175
176
177
178
179
180
181
182
183
184
185
186 public EigenDecompositionImpl(final double[] main, double[] secondary,
187 final double splitTolerance)
188 throws InvalidMatrixException {
189
190 this.main = main.clone();
191 this.secondary = secondary.clone();
192 transformer = null;
193
194
195 squaredSecondary = new double[secondary.length];
196 for (int i = 0; i < squaredSecondary.length; ++i) {
197 final double s = secondary[i];
198 squaredSecondary[i] = s * s;
199 }
200
201 this.splitTolerance = splitTolerance;
202 decompose();
203
204 }
205
206
207
208
209
210
211 private boolean isSymmetric(final RealMatrix matrix) {
212 final int rows = matrix.getRowDimension();
213 final int columns = matrix.getColumnDimension();
214 final double eps = 10 * rows * columns * MathUtils.EPSILON;
215 for (int i = 0; i < rows; ++i) {
216 for (int j = i + 1; j < columns; ++j) {
217 final double mij = matrix.getEntry(i, j);
218 final double mji = matrix.getEntry(j, i);
219 if (Math.abs(mij - mji) > (Math.max(Math.abs(mij), Math.abs(mji)) * eps)) {
220 return false;
221 }
222 }
223 }
224 return true;
225 }
226
227
228
229
230
231
232 private void decompose() {
233
234 cachedV = null;
235 cachedD = null;
236 cachedVt = null;
237 work = new double[6 * main.length];
238
239
240 computeGershgorinCircles();
241
242
243 findEigenvalues();
244
245
246 eigenvectors = null;
247
248 }
249
250
251 public RealMatrix getV()
252 throws InvalidMatrixException {
253
254 if (cachedV == null) {
255
256 if (eigenvectors == null) {
257 findEigenVectors();
258 }
259
260 final int m = eigenvectors.length;
261 cachedV = MatrixUtils.createRealMatrix(m, m);
262 for (int k = 0; k < m; ++k) {
263 cachedV.setColumnVector(k, eigenvectors[k]);
264 }
265
266 }
267
268
269 return cachedV;
270
271 }
272
273
274 public RealMatrix getD()
275 throws InvalidMatrixException {
276 if (cachedD == null) {
277
278 cachedD = MatrixUtils.createRealDiagonalMatrix(realEigenvalues);
279 }
280 return cachedD;
281 }
282
283
284 public RealMatrix getVT()
285 throws InvalidMatrixException {
286
287 if (cachedVt == null) {
288
289 if (eigenvectors == null) {
290 findEigenVectors();
291 }
292
293 final int m = eigenvectors.length;
294 cachedVt = MatrixUtils.createRealMatrix(m, m);
295 for (int k = 0; k < m; ++k) {
296 cachedVt.setRowVector(k, eigenvectors[k]);
297 }
298
299 }
300
301
302 return cachedVt;
303
304 }
305
306
307 public double[] getRealEigenvalues()
308 throws InvalidMatrixException {
309 return realEigenvalues.clone();
310 }
311
312
313 public double getRealEigenvalue(final int i)
314 throws InvalidMatrixException, ArrayIndexOutOfBoundsException {
315 return realEigenvalues[i];
316 }
317
318
319 public double[] getImagEigenvalues()
320 throws InvalidMatrixException {
321 return imagEigenvalues.clone();
322 }
323
324
325 public double getImagEigenvalue(final int i)
326 throws InvalidMatrixException, ArrayIndexOutOfBoundsException {
327 return imagEigenvalues[i];
328 }
329
330
331 public RealVector getEigenvector(final int i)
332 throws InvalidMatrixException, ArrayIndexOutOfBoundsException {
333 if (eigenvectors == null) {
334 findEigenVectors();
335 }
336 return eigenvectors[i].copy();
337 }
338
339
340
341
342
343 public double getDeterminant() {
344 double determinant = 1;
345 for (double lambda : realEigenvalues) {
346 determinant *= lambda;
347 }
348 return determinant;
349 }
350
351
352 public DecompositionSolver getSolver() {
353 if (eigenvectors == null) {
354 findEigenVectors();
355 }
356 return new Solver(realEigenvalues, imagEigenvalues, eigenvectors);
357 }
358
359
360 private static class Solver implements DecompositionSolver {
361
362
363 private double[] realEigenvalues;
364
365
366 private double[] imagEigenvalues;
367
368
369 private final ArrayRealVector[] eigenvectors;
370
371
372
373
374
375
376
377 private Solver(final double[] realEigenvalues, final double[] imagEigenvalues,
378 final ArrayRealVector[] eigenvectors) {
379 this.realEigenvalues = realEigenvalues;
380 this.imagEigenvalues = imagEigenvalues;
381 this.eigenvectors = eigenvectors;
382 }
383
384
385
386
387
388
389
390
391
392 public double[] solve(final double[] b)
393 throws IllegalArgumentException, InvalidMatrixException {
394
395 if (!isNonSingular()) {
396 throw new SingularMatrixException();
397 }
398
399 final int m = realEigenvalues.length;
400 if (b.length != m) {
401 throw MathRuntimeException.createIllegalArgumentException(
402 "vector length mismatch: got {0} but expected {1}",
403 b.length, m);
404 }
405
406 final double[] bp = new double[m];
407 for (int i = 0; i < m; ++i) {
408 final ArrayRealVector v = eigenvectors[i];
409 final double[] vData = v.getDataRef();
410 final double s = v.dotProduct(b) / realEigenvalues[i];
411 for (int j = 0; j < m; ++j) {
412 bp[j] += s * vData[j];
413 }
414 }
415
416 return bp;
417
418 }
419
420
421
422
423
424
425
426
427
428 public RealVector solve(final RealVector b)
429 throws IllegalArgumentException, InvalidMatrixException {
430
431 if (!isNonSingular()) {
432 throw new SingularMatrixException();
433 }
434
435 final int m = realEigenvalues.length;
436 if (b.getDimension() != m) {
437 throw MathRuntimeException.createIllegalArgumentException(
438 "vector length mismatch: got {0} but expected {1}",
439 b.getDimension(), m);
440 }
441
442 final double[] bp = new double[m];
443 for (int i = 0; i < m; ++i) {
444 final ArrayRealVector v = eigenvectors[i];
445 final double[] vData = v.getDataRef();
446 final double s = v.dotProduct(b) / realEigenvalues[i];
447 for (int j = 0; j < m; ++j) {
448 bp[j] += s * vData[j];
449 }
450 }
451
452 return new ArrayRealVector(bp, false);
453
454 }
455
456
457
458
459
460
461
462
463
464 public RealMatrix solve(final RealMatrix b)
465 throws IllegalArgumentException, InvalidMatrixException {
466
467 if (!isNonSingular()) {
468 throw new SingularMatrixException();
469 }
470
471 final int m = realEigenvalues.length;
472 if (b.getRowDimension() != m) {
473 throw MathRuntimeException.createIllegalArgumentException(
474 "dimensions mismatch: got {0}x{1} but expected {2}x{3}",
475 b.getRowDimension(), b.getColumnDimension(), m, "n");
476 }
477
478 final int nColB = b.getColumnDimension();
479 final double[][] bp = new double[m][nColB];
480 for (int k = 0; k < nColB; ++k) {
481 for (int i = 0; i < m; ++i) {
482 final ArrayRealVector v = eigenvectors[i];
483 final double[] vData = v.getDataRef();
484 double s = 0;
485 for (int j = 0; j < m; ++j) {
486 s += v.getEntry(j) * b.getEntry(j, k);
487 }
488 s /= realEigenvalues[i];
489 for (int j = 0; j < m; ++j) {
490 bp[j][k] += s * vData[j];
491 }
492 }
493 }
494
495 return MatrixUtils.createRealMatrix(bp);
496
497 }
498
499
500
501
502
503 public boolean isNonSingular() {
504 for (int i = 0; i < realEigenvalues.length; ++i) {
505 if ((realEigenvalues[i] == 0) && (imagEigenvalues[i] == 0)) {
506 return false;
507 }
508 }
509 return true;
510 }
511
512
513
514
515
516 public RealMatrix getInverse()
517 throws InvalidMatrixException {
518
519 if (!isNonSingular()) {
520 throw new SingularMatrixException();
521 }
522
523 final int m = realEigenvalues.length;
524 final double[][] invData = new double[m][m];
525
526 for (int i = 0; i < m; ++i) {
527 final double[] invI = invData[i];
528 for (int j = 0; j < m; ++j) {
529 double invIJ = 0;
530 for (int k = 0; k < m; ++k) {
531 final double[] vK = eigenvectors[k].getDataRef();
532 invIJ += vK[i] * vK[j] / realEigenvalues[k];
533 }
534 invI[j] = invIJ;
535 }
536 }
537 return MatrixUtils.createRealMatrix(invData);
538
539 }
540
541 }
542
543
544
545
546
547 private void transformToTridiagonal(final RealMatrix matrix) {
548
549
550 transformer = new TriDiagonalTransformer(matrix);
551 main = transformer.getMainDiagonalRef();
552 secondary = transformer.getSecondaryDiagonalRef();
553
554
555 squaredSecondary = new double[secondary.length];
556 for (int i = 0; i < squaredSecondary.length; ++i) {
557 final double s = secondary[i];
558 squaredSecondary[i] = s * s;
559 }
560
561 }
562
563
564
565
566 private void computeGershgorinCircles() {
567
568 final int m = main.length;
569 final int lowerStart = 4 * m;
570 final int upperStart = 5 * m;
571 lowerSpectra = Double.POSITIVE_INFINITY;
572 upperSpectra = Double.NEGATIVE_INFINITY;
573 double eMax = 0;
574
575 double eCurrent = 0;
576 for (int i = 0; i < m - 1; ++i) {
577
578 final double dCurrent = main[i];
579 final double ePrevious = eCurrent;
580 eCurrent = Math.abs(secondary[i]);
581 eMax = Math.max(eMax, eCurrent);
582 final double radius = ePrevious + eCurrent;
583
584 final double lower = dCurrent - radius;
585 work[lowerStart + i] = lower;
586 lowerSpectra = Math.min(lowerSpectra, lower);
587
588 final double upper = dCurrent + radius;
589 work[upperStart + i] = upper;
590 upperSpectra = Math.max(upperSpectra, upper);
591
592 }
593
594 final double dCurrent = main[m - 1];
595 work[lowerStart + m - 1] = dCurrent - eCurrent;
596 work[upperStart + m - 1] = dCurrent + eCurrent;
597 minPivot = MathUtils.SAFE_MIN * Math.max(1.0, eMax * eMax);
598
599 }
600
601
602
603
604
605 private void findEigenvalues()
606 throws InvalidMatrixException {
607
608
609 List<Integer> splitIndices = computeSplits();
610
611
612 realEigenvalues = new double[main.length];
613 imagEigenvalues = new double[main.length];
614 int begin = 0;
615 for (final int end : splitIndices) {
616 final int n = end - begin;
617 switch (n) {
618
619 case 1:
620
621 process1RowBlock(begin);
622 break;
623
624 case 2:
625
626 process2RowsBlock(begin);
627 break;
628
629 case 3:
630
631 process3RowsBlock(begin);
632 break;
633
634 default:
635
636
637 final double[] range = eigenvaluesRange(begin, n);
638 final double oneFourth = 0.25 * (3 * range[0] + range[1]);
639 final int oneFourthCount = countEigenValues(oneFourth, begin, n);
640 final double threeFourth = 0.25 * (range[0] + 3 * range[1]);
641 final int threeFourthCount = countEigenValues(threeFourth, begin, n);
642 final boolean chooseLeft = (oneFourthCount - 1) >= (n - threeFourthCount);
643 final double lambda = chooseLeft ? range[0] : range[1];
644
645 tau = (range[1] - range[0]) * MathUtils.EPSILON * n + 2 * minPivot;
646
647
648 ldlTDecomposition(lambda, begin, n);
649
650
651 processGeneralBlock(n);
652
653
654 if (chooseLeft) {
655 for (int i = 0; i < n; ++i) {
656 realEigenvalues[begin + i] = lambda + work[4 * i];
657 }
658 } else {
659 for (int i = 0; i < n; ++i) {
660 realEigenvalues[begin + i] = lambda - work[4 * i];
661 }
662 }
663
664 }
665 begin = end;
666 }
667
668
669 Arrays.sort(realEigenvalues);
670 for (int i = 0, j = realEigenvalues.length - 1; i < j; ++i, --j) {
671 final double tmp = realEigenvalues[i];
672 realEigenvalues[i] = realEigenvalues[j];
673 realEigenvalues[j] = tmp;
674 }
675
676 }
677
678
679
680
681
682 private List<Integer> computeSplits() {
683
684 final List<Integer> list = new ArrayList<Integer>();
685
686
687 double absDCurrent = Math.abs(main[0]);
688 for (int i = 0; i < secondary.length; ++i) {
689 final double absDPrevious = absDCurrent;
690 absDCurrent = Math.abs(main[i + 1]);
691 final double max = splitTolerance * Math.sqrt(absDPrevious * absDCurrent);
692 if (Math.abs(secondary[i]) <= max) {
693 list.add(i + 1);
694 secondary[i] = 0;
695 squaredSecondary[i] = 0;
696 }
697 }
698
699 list.add(secondary.length + 1);
700 return list;
701
702 }
703
704
705
706
707
708
709 private void process1RowBlock(final int index) {
710 realEigenvalues[index] = main[index];
711 }
712
713
714
715
716
717
718
719 private void process2RowsBlock(final int index)
720 throws InvalidMatrixException {
721
722
723
724 final double q0 = main[index];
725 final double q1 = main[index + 1];
726 final double e12 = squaredSecondary[index];
727
728 final double s = q0 + q1;
729 final double p = q0 * q1 - e12;
730 final double delta = s * s - 4 * p;
731 if (delta < 0) {
732 throw new InvalidMatrixException("cannot solve degree {0} equation", 2);
733 }
734
735 final double largestRoot = 0.5 * (s + Math.sqrt(delta));
736 realEigenvalues[index] = largestRoot;
737 realEigenvalues[index + 1] = p / largestRoot;
738
739 }
740
741
742
743
744
745
746
747 private void process3RowsBlock(final int index)
748 throws InvalidMatrixException {
749
750
751
752 final double q0 = main[index];
753 final double q1 = main[index + 1];
754 final double q2 = main[index + 2];
755 final double e12 = squaredSecondary[index];
756 final double q1q2Me22 = q1 * q2 - squaredSecondary[index + 1];
757
758
759 final double b = -(q0 + q1 + q2);
760 final double c = q0 * q1 + q0 * q2 + q1q2Me22 - e12;
761 final double d = q2 * e12 - q0 * q1q2Me22;
762
763
764 final double b2 = b * b;
765 final double q = (3 * c - b2) / 9;
766 final double r = ((9 * c - 2 * b2) * b - 27 * d) / 54;
767 final double delta = q * q * q + r * r;
768 if (delta >= 0) {
769
770
771
772 throw new InvalidMatrixException("cannot solve degree {0} equation", 3);
773 }
774 final double sqrtMq = Math.sqrt(-q);
775 final double theta = Math.acos(r / (-q * sqrtMq));
776 final double alpha = 2 * sqrtMq;
777 final double beta = b / 3;
778
779 double z0 = alpha * Math.cos(theta / 3) - beta;
780 double z1 = alpha * Math.cos((theta + 2 * Math.PI) / 3) - beta;
781 double z2 = alpha * Math.cos((theta + 4 * Math.PI) / 3) - beta;
782 if (z0 < z1) {
783 final double t = z0;
784 z0 = z1;
785 z1 = t;
786 }
787 if (z1 < z2) {
788 final double t = z1;
789 z1 = z2;
790 z2 = t;
791 }
792 if (z0 < z1) {
793 final double t = z0;
794 z0 = z1;
795 z1 = t;
796 }
797 realEigenvalues[index] = z0;
798 realEigenvalues[index + 1] = z1;
799 realEigenvalues[index + 2] = z2;
800
801 }
802
803
804
805
806
807
808
809
810
811
812
813
814 private void processGeneralBlock(final int n)
815 throws InvalidMatrixException {
816
817
818 double sumOffDiag = 0;
819 for (int i = 0; i < n - 1; ++i) {
820 final int fourI = 4 * i;
821 final double ei = work[fourI + 2];
822 sumOffDiag += ei;
823 }
824
825 if (sumOffDiag == 0) {
826
827 return;
828 }
829
830
831 flipIfWarranted(n, 2);
832
833
834 initialSplits(n);
835
836
837 tType = 0;
838 dMin1 = 0;
839 dMin2 = 0;
840 dN = 0;
841 dN1 = 0;
842 dN2 = 0;
843 tau = 0;
844
845
846 int i0 = 0;
847 int n0 = n;
848 while (n0 > 0) {
849
850
851 sigma = (n0 == n) ? 0 : -work[4 * n0 - 2];
852 sigmaLow = 0;
853
854
855 double eMin = (i0 == n0) ? 0 : work[4 * n0 - 6];
856 double eMax = 0;
857 double qMax = work[4 * n0 - 4];
858 double qMin = qMax;
859 i0 = 0;
860 for (int i = 4 * (n0 - 2); i >= 0; i -= 4) {
861 if (work[i + 2] <= 0) {
862 i0 = 1 + i / 4;
863 break;
864 }
865 if (qMin >= 4 * eMax) {
866 qMin = Math.min(qMin, work[i + 4]);
867 eMax = Math.max(eMax, work[i + 2]);
868 }
869 qMax = Math.max(qMax, work[i] + work[i + 2]);
870 eMin = Math.min(eMin, work[i + 2]);
871 }
872 work[4 * n0 - 2] = eMin;
873
874
875 dMin = -Math.max(0, qMin - 2 * Math.sqrt(qMin * eMax));
876
877 pingPong = 0;
878 int maxIter = 30 * (n0 - i0);
879 for (int k = 0; i0 < n0; ++k) {
880 if (k >= maxIter) {
881 throw new InvalidMatrixException(new MaxIterationsExceededException(maxIter));
882 }
883
884
885 n0 = goodStep(i0, n0);
886 pingPong = 1 - pingPong;
887
888
889
890 if ((pingPong == 0) && (n0 - i0 > 3) &&
891 (work[4 * n0 - 1] <= TOLERANCE_2 * qMax) &&
892 (work[4 * n0 - 2] <= TOLERANCE_2 * sigma)) {
893 int split = i0 - 1;
894 qMax = work[4 * i0];
895 eMin = work[4 * i0 + 2];
896 double previousEMin = work[4 * i0 + 3];
897 for (int i = 4 * i0; i < 4 * n0 - 11; i += 4) {
898 if ((work[i + 3] <= TOLERANCE_2 * work[i]) &&
899 (work[i + 2] <= TOLERANCE_2 * sigma)) {
900
901 work[i + 2] = -sigma;
902 split = i / 4;
903 qMax = 0;
904 eMin = work[i + 6];
905 previousEMin = work[i + 7];
906 } else {
907 qMax = Math.max(qMax, work[i + 4]);
908 eMin = Math.min(eMin, work[i + 2]);
909 previousEMin = Math.min(previousEMin, work[i + 3]);
910 }
911 }
912 work[4 * n0 - 2] = eMin;
913 work[4 * n0 - 1] = previousEMin;
914 i0 = split + 1;
915 }
916 }
917
918 }
919
920 }
921
922
923
924
925
926 private void initialSplits(final int n) {
927
928 pingPong = 0;
929 for (int k = 0; k < 2; ++k) {
930
931
932 double d = work[4 * (n - 1) + pingPong];
933 for (int i = 4 * (n - 2) + pingPong; i >= 0; i -= 4) {
934 if (work[i + 2] <= TOLERANCE_2 * d) {
935 work[i + 2] = -0.0;
936 d = work[i];
937 } else {
938 d *= work[i] / (d + work[i + 2]);
939 }
940 }
941
942
943 d = work[pingPong];
944 for (int i = 2 + pingPong; i < 4 * n - 2; i += 4) {
945 final int j = i - 2 * pingPong - 1;
946 work[j] = d + work[i];
947 if (work[i] <= TOLERANCE_2 * d) {
948 work[i] = -0.0;
949 work[j] = d;
950 work[j + 2] = 0.0;
951 d = work[i + 2];
952 } else if ((MathUtils.SAFE_MIN * work[i + 2] < work[j]) &&
953 (MathUtils.SAFE_MIN * work[j] < work[i + 2])) {
954 final double tmp = work[i + 2] / work[j];
955 work[j + 2] = work[i] * tmp;
956 d *= tmp;
957 } else {
958 work[j + 2] = work[i + 2] * (work[i] / work[j]);
959 d *= work[i + 2] / work[j];
960 }
961 }
962 work[4 * n - 3 - pingPong] = d;
963
964
965 pingPong = 1 - pingPong;
966
967 }
968
969 }
970
971
972
973
974
975
976
977
978
979
980
981
982 private int goodStep(final int start, final int end) {
983
984 g = 0.0;
985
986
987 int deflatedEnd = end;
988 for (boolean deflating = true; deflating;) {
989
990 if (start >= deflatedEnd) {
991
992 return deflatedEnd;
993 }
994
995 final int k = 4 * deflatedEnd + pingPong - 1;
996
997 if ((start == deflatedEnd - 1) ||
998 ((start != deflatedEnd - 2) &&
999 ((work[k - 5] <= TOLERANCE_2 * (sigma + work[k - 3])) ||
1000 (work[k - 2 * pingPong - 4] <= TOLERANCE_2 * work[k - 7])))) {
1001
1002
1003 work[4 * deflatedEnd - 4] = sigma + work[4 * deflatedEnd - 4 + pingPong];
1004 deflatedEnd -= 1;
1005
1006 } else if ((start == deflatedEnd - 2) ||
1007 (work[k - 9] <= TOLERANCE_2 * sigma) ||
1008 (work[k - 2 * pingPong - 8] <= TOLERANCE_2 * work[k - 11])) {
1009
1010
1011 if (work[k - 3] > work[k - 7]) {
1012 final double tmp = work[k - 3];
1013 work[k - 3] = work[k - 7];
1014 work[k - 7] = tmp;
1015 }
1016
1017 if (work[k - 5] > TOLERANCE_2 * work[k - 3]) {
1018 double t = 0.5 * ((work[k - 7] - work[k - 3]) + work[k - 5]);
1019 double s = work[k - 3] * (work[k - 5] / t);
1020 if (s <= t) {
1021 s = work[k - 3] * work[k - 5] / (t * (1 + Math.sqrt(1 + s / t)));
1022 } else {
1023 s = work[k - 3] * work[k - 5] / (t + Math.sqrt(t * (t + s)));
1024 }
1025 t = work[k - 7] + (s + work[k - 5]);
1026 work[k - 3] *= work[k - 7] / t;
1027 work[k - 7] = t;
1028 }
1029 work[4 * deflatedEnd - 8] = sigma + work[k - 7];
1030 work[4 * deflatedEnd - 4] = sigma + work[k - 3];
1031 deflatedEnd -= 2;
1032 } else {
1033
1034
1035 deflating = false;
1036
1037 }
1038
1039 }
1040
1041 final int l = 4 * deflatedEnd + pingPong - 1;
1042
1043
1044 if ((dMin <= 0) || (deflatedEnd < end)) {
1045 if (flipIfWarranted(deflatedEnd, 1)) {
1046 dMin2 = Math.min(dMin2, work[l - 1]);
1047 work[l - 1] =
1048 Math.min(work[l - 1],
1049 Math.min(work[3 + pingPong], work[7 + pingPong]));
1050 work[l - 2 * pingPong] =
1051 Math.min(work[l - 2 * pingPong],
1052 Math.min(work[6 + pingPong], work[6 + pingPong]));
1053 qMax = Math.max(qMax, Math.max(work[3 + pingPong], work[7 + pingPong]));
1054 dMin = -0.0;
1055 }
1056 }
1057
1058 if ((dMin < 0) ||
1059 (MathUtils.SAFE_MIN * qMax < Math.min(work[l - 1],
1060 Math.min(work[l - 9],
1061 dMin2 + work[l - 2 * pingPong])))) {
1062
1063 computeShiftIncrement(start, deflatedEnd, end - deflatedEnd);
1064
1065
1066 for (boolean loop = true; loop;) {
1067
1068
1069 dqds(start, deflatedEnd);
1070
1071
1072 if ((dMin >= 0) && (dMin1 > 0)) {
1073
1074 updateSigma(tau);
1075 return deflatedEnd;
1076 } else if ((dMin < 0.0) &&
1077 (dMin1 > 0.0) &&
1078 (work[4 * deflatedEnd - 5 - pingPong] < TOLERANCE * (sigma + dN1)) &&
1079 (Math.abs(dN) < TOLERANCE * sigma)) {
1080
1081 work[4 * deflatedEnd - 3 - pingPong] = 0.0;
1082 dMin = 0.0;
1083 updateSigma(tau);
1084 return deflatedEnd;
1085 } else if (dMin < 0.0) {
1086
1087 if (tType < -22) {
1088
1089 tau = 0.0;
1090 } else if (dMin1 > 0.0) {
1091
1092 tau = (tau + dMin) * (1.0 - 2.0 * MathUtils.EPSILON);
1093 tType -= 11;
1094 } else {
1095
1096 tau *= 0.25;
1097 tType -= 12;
1098 }
1099 } else if (Double.isNaN(dMin)) {
1100 tau = 0.0;
1101 } else {
1102
1103 loop = false;
1104 }
1105 }
1106
1107 }
1108
1109
1110 dqd(start, deflatedEnd);
1111
1112 return deflatedEnd;
1113
1114 }
1115
1116
1117
1118
1119
1120
1121
1122
1123 private boolean flipIfWarranted(final int n, final int step) {
1124 if (1.5 * work[pingPong] < work[4 * (n - 1) + pingPong]) {
1125
1126 for (int i = 0, j = 4 * n - 1; i < j; i += 4, j -= 4) {
1127 for (int k = 0; k < 4; k += step) {
1128 final double tmp = work[i + k];
1129 work[i + k] = work[j - k];
1130 work[j - k] = tmp;
1131 }
1132 }
1133 return true;
1134 }
1135 return false;
1136 }
1137
1138
1139
1140
1141
1142
1143
1144 private double[] eigenvaluesRange(final int index, final int n) {
1145
1146
1147 final int lowerStart = 4 * main.length;
1148 final int upperStart = 5 * main.length;
1149 double lower = Double.POSITIVE_INFINITY;
1150 double upper = Double.NEGATIVE_INFINITY;
1151 for (int i = 0; i < n; ++i) {
1152 lower = Math.min(lower, work[lowerStart + index +i]);
1153 upper = Math.max(upper, work[upperStart + index +i]);
1154 }
1155
1156
1157 final double tNorm = Math.max(Math.abs(lower), Math.abs(upper));
1158 final double relativeTolerance = Math.sqrt(MathUtils.EPSILON);
1159 final double absoluteTolerance = 4 * minPivot;
1160 final int maxIter =
1161 2 + (int) ((Math.log(tNorm + minPivot) - Math.log(minPivot)) / Math.log(2.0));
1162 final double margin = 2 * (tNorm * MathUtils.EPSILON * n + 2 * minPivot);
1163
1164
1165 double left = lower - margin;
1166 double right = upper + margin;
1167 for (int i = 0; i < maxIter; ++i) {
1168
1169 final double range = right - left;
1170 if ((range < absoluteTolerance) ||
1171 (range < relativeTolerance * Math.max(Math.abs(left), Math.abs(right)))) {
1172
1173 break;
1174 }
1175
1176 final double middle = 0.5 * (left + right);
1177 if (countEigenValues(middle, index, n) >= 1) {
1178 right = middle;
1179 } else {
1180 left = middle;
1181 }
1182
1183 }
1184 lower = Math.max(lower, left - 100 * MathUtils.EPSILON * Math.abs(left));
1185
1186
1187 left = lower - margin;
1188 right = upper + margin;
1189 for (int i = 0; i < maxIter; ++i) {
1190
1191 final double range = right - left;
1192 if ((range < absoluteTolerance) ||
1193 (range < relativeTolerance * Math.max(Math.abs(left), Math.abs(right)))) {
1194
1195 break;
1196 }
1197
1198 final double middle = 0.5 * (left + right);
1199 if (countEigenValues(middle, index, n) >= n) {
1200 right = middle;
1201 } else {
1202 left = middle;
1203 }
1204
1205 }
1206 upper = Math.min(upper, right + 100 * MathUtils.EPSILON * Math.abs(right));
1207
1208 return new double[] { lower, upper };
1209
1210 }
1211
1212
1213
1214
1215
1216
1217
1218
1219 private int countEigenValues(final double t, final int index, final int n) {
1220 double ratio = main[index] - t;
1221 int count = (ratio > 0) ? 0 : 1;
1222 for (int i = 1; i < n; ++i) {
1223 ratio = main[index + i] - squaredSecondary[index + i - 1] / ratio - t;
1224 if (ratio <= 0) {
1225 ++count;
1226 }
1227 }
1228 return count;
1229 }
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242 private void ldlTDecomposition(final double lambda, final int index, final int n) {
1243 double di = main[index] - lambda;
1244 work[0] = Math.abs(di);
1245 for (int i = 1; i < n; ++i) {
1246 final int fourI = 4 * i;
1247 final double eiM1 = secondary[index + i - 1];
1248 final double ratio = eiM1 / di;
1249 work[fourI - 2] = ratio * ratio * Math.abs(di);
1250 di = (main[index + i] - lambda) - eiM1 * ratio;
1251 work[fourI] = Math.abs(di);
1252 }
1253 }
1254
1255
1256
1257
1258
1259
1260
1261 private void dqds(final int start, final int end) {
1262
1263 eMin = work[4 * start + pingPong + 4];
1264 double d = work[4 * start + pingPong] - tau;
1265 dMin = d;
1266 dMin1 = -work[4 * start + pingPong];
1267
1268 if (pingPong == 0) {
1269 for (int j4 = 4 * start + 3; j4 <= 4 * (end - 3); j4 += 4) {
1270 work[j4 - 2] = d + work[j4 - 1];
1271 final double tmp = work[j4 + 1] / work[j4 - 2];
1272 d = d * tmp - tau;
1273 dMin = Math.min(dMin, d);
1274 work[j4] = work[j4 - 1] * tmp;
1275 eMin = Math.min(work[j4], eMin);
1276 }
1277 } else {
1278 for (int j4 = 4 * start + 3; j4 <= 4 * (end - 3); j4 += 4) {
1279 work[j4 - 3] = d + work[j4];
1280 final double tmp = work[j4 + 2] / work[j4 - 3];
1281 d = d * tmp - tau;
1282 dMin = Math.min(dMin, d);
1283 work[j4 - 1] = work[j4] * tmp;
1284 eMin = Math.min(work[j4 - 1], eMin);
1285 }
1286 }
1287
1288
1289 dN2 = d;
1290 dMin2 = dMin;
1291 int j4 = 4 * (end - 2) - pingPong - 1;
1292 int j4p2 = j4 + 2 * pingPong - 1;
1293 work[j4 - 2] = dN2 + work[j4p2];
1294 work[j4] = work[j4p2 + 2] * (work[j4p2] / work[j4 - 2]);
1295 dN1 = work[j4p2 + 2] * (dN2 / work[j4 - 2]) - tau;
1296 dMin = Math.min(dMin, dN1);
1297
1298 dMin1 = dMin;
1299 j4 = j4 + 4;
1300 j4p2 = j4 + 2 * pingPong - 1;
1301 work[j4 - 2] = dN1 + work[j4p2];
1302 work[j4] = work[j4p2 + 2] * (work[j4p2] / work[j4 - 2]);
1303 dN = work[j4p2 + 2] * (dN1 / work[j4 - 2]) - tau;
1304 dMin = Math.min(dMin, dN);
1305
1306 work[j4 + 2] = dN;
1307 work[4 * end - pingPong - 1] = eMin;
1308
1309 }
1310
1311
1312
1313
1314
1315
1316
1317
1318 private void dqd(final int start, final int end) {
1319
1320 eMin = work[4 * start + pingPong + 4];
1321 double d = work[4 * start + pingPong];
1322 dMin = d;
1323
1324 if (pingPong == 0) {
1325 for (int j4 = 4 * start + 3; j4 < 4 * (end - 3); j4 += 4) {
1326 work[j4 - 2] = d + work[j4 - 1];
1327 if (work[j4 - 2] == 0.0) {
1328 work[j4] = 0.0;
1329 d = work[j4 + 1];
1330 dMin = d;
1331 eMin = 0.0;
1332 } else if ((MathUtils.SAFE_MIN * work[j4 + 1] < work[j4 - 2]) &&
1333 (MathUtils.SAFE_MIN * work[j4 - 2] < work[j4 + 1])) {
1334 final double tmp = work[j4 + 1] / work[j4 - 2];
1335 work[j4] = work[j4 - 1] * tmp;
1336 d *= tmp;
1337 } else {
1338 work[j4] = work[j4 + 1] * (work[j4 - 1] / work[j4 - 2]);
1339 d *= work[j4 + 1] / work[j4 - 2];
1340 }
1341 dMin = Math.min(dMin, d);
1342 eMin = Math.min(eMin, work[j4]);
1343 }
1344 } else {
1345 for (int j4 = 4 * start + 3; j4 < 4 * (end - 3); j4 += 4) {
1346 work[j4 - 3] = d + work[j4];
1347 if (work[j4 - 3] == 0.0) {
1348 work[j4 - 1] = 0.0;
1349 d = work[j4 + 2];
1350 dMin = d;
1351 eMin = 0.0;
1352 } else if ((MathUtils.SAFE_MIN * work[j4 + 2] < work[j4 - 3]) &&
1353 (MathUtils.SAFE_MIN * work[j4 - 3] < work[j4 + 2])) {
1354 final double tmp = work[j4 + 2] / work[j4 - 3];
1355 work[j4 - 1] = work[j4] * tmp;
1356 d *= tmp;
1357 } else {
1358 work[j4 - 1] = work[j4 + 2] * (work[j4] / work[j4 - 3]);
1359 d *= work[j4 + 2] / work[j4 - 3];
1360 }
1361 dMin = Math.min(dMin, d);
1362 eMin = Math.min(eMin, work[j4 - 1]);
1363 }
1364 }
1365
1366
1367 dN2 = d;
1368 dMin2 = dMin;
1369 int j4 = 4 * (end - 2) - pingPong - 1;
1370 int j4p2 = j4 + 2 * pingPong - 1;
1371 work[j4 - 2] = dN2 + work[j4p2];
1372 if (work[j4 - 2] == 0.0) {
1373 work[j4] = 0.0;
1374 dN1 = work[j4p2 + 2];
1375 dMin = dN1;
1376 eMin = 0.0;
1377 } else if ((MathUtils.SAFE_MIN * work[j4p2 + 2] < work[j4 - 2]) &&
1378 (MathUtils.SAFE_MIN * work[j4 - 2] < work[j4p2 + 2])) {
1379 final double tmp = work[j4p2 + 2] / work[j4 - 2];
1380 work[j4] = work[j4p2] * tmp;
1381 dN1 = dN2 * tmp;
1382 } else {
1383 work[j4] = work[j4p2 + 2] * (work[j4p2] / work[j4 - 2]);
1384 dN1 = work[j4p2 + 2] * (dN2 / work[j4 - 2]);
1385 }
1386 dMin = Math.min(dMin, dN1);
1387
1388 dMin1 = dMin;
1389 j4 = j4 + 4;
1390 j4p2 = j4 + 2 * pingPong - 1;
1391 work[j4 - 2] = dN1 + work[j4p2];
1392 if (work[j4 - 2] == 0.0) {
1393 work[j4] = 0.0;
1394 dN = work[j4p2 + 2];
1395 dMin = dN;
1396 eMin = 0.0;
1397 } else if ((MathUtils.SAFE_MIN * work[j4p2 + 2] < work[j4 - 2]) &&
1398 (MathUtils.SAFE_MIN * work[j4 - 2] < work[j4p2 + 2])) {
1399 final double tmp = work[j4p2 + 2] / work[j4 - 2];
1400 work[j4] = work[j4p2] * tmp;
1401 dN = dN1 * tmp;
1402 } else {
1403 work[j4] = work[j4p2 + 2] * (work[j4p2] / work[j4 - 2]);
1404 dN = work[j4p2 + 2] * (dN1 / work[j4 - 2]);
1405 }
1406 dMin = Math.min(dMin, dN);
1407
1408 work[j4 + 2] = dN;
1409 work[4 * end - pingPong - 1] = eMin;
1410
1411 }
1412
1413
1414
1415
1416
1417
1418
1419
1420 private void computeShiftIncrement(final int start, final int end, final int deflated) {
1421
1422 final double cnst1 = 0.563;
1423 final double cnst2 = 1.010;
1424 final double cnst3 = 1.05;
1425
1426
1427
1428 if (dMin <= 0.0) {
1429 tau = -dMin;
1430 tType = -1;
1431 return;
1432 }
1433
1434 int nn = 4 * end + pingPong - 1;
1435 switch (deflated) {
1436
1437 case 0 :
1438 if (dMin == dN || dMin == dN1) {
1439
1440 double b1 = Math.sqrt(work[nn - 3]) * Math.sqrt(work[nn - 5]);
1441 double b2 = Math.sqrt(work[nn - 7]) * Math.sqrt(work[nn - 9]);
1442 double a2 = work[nn - 7] + work[nn - 5];
1443
1444 if (dMin == dN && dMin1 == dN1) {
1445
1446 final double gap2 = dMin2 - a2 - dMin2 * 0.25;
1447 final double gap1 = a2 - dN - ((gap2 > 0.0 && gap2 > b2) ? (b2 / gap2) * b2 : (b1 + b2));
1448 if (gap1 > 0.0 && gap1 > b1) {
1449 tau = Math.max(dN - (b1 / gap1) * b1, 0.5 * dMin);
1450 tType = -2;
1451 } else {
1452 double s = 0.0;
1453 if (dN > b1) {
1454 s = dN - b1;
1455 }
1456 if (a2 > (b1 + b2)) {
1457 s = Math.min(s, a2 - (b1 + b2));
1458 }
1459 tau = Math.max(s, 0.333 * dMin);
1460 tType = -3;
1461 }
1462 } else {
1463
1464 tType = -4;
1465 double s = 0.25 * dMin;
1466 double gam;
1467 int np;
1468 if (dMin == dN) {
1469 gam = dN;
1470 a2 = 0.0;
1471 if (work[nn - 5] > work[nn - 7]) {
1472 return;
1473 }
1474 b2 = work[nn - 5] / work[nn - 7];
1475 np = nn - 9;
1476 } else {
1477 np = nn - 2 * pingPong;
1478 b2 = work[np - 2];
1479 gam = dN1;
1480 if (work[np - 4] > work[np - 2]) {
1481 return;
1482 }
1483 a2 = work[np - 4] / work[np - 2];
1484 if (work[nn - 9] > work[nn - 11]) {
1485 return;
1486 }
1487 b2 = work[nn - 9] / work[nn - 11];
1488 np = nn - 13;
1489 }
1490
1491
1492 a2 = a2 + b2;
1493 for (int i4 = np; i4 >= 4 * start + 2 + pingPong; i4 -= 4) {
1494 if(b2 == 0.0) {
1495 break;
1496 }
1497 b1 = b2;
1498 if (work[i4] > work[i4 - 2]) {
1499 return;
1500 }
1501 b2 = b2 * (work[i4] / work[i4 - 2]);
1502 a2 = a2 + b2;
1503 if (100 * Math.max(b2, b1) < a2 || cnst1 < a2) {
1504 break;
1505 }
1506 }
1507 a2 = cnst3 * a2;
1508
1509
1510 if (a2 < cnst1) {
1511 s = gam * (1 - Math.sqrt(a2)) / (1 + a2);
1512 }
1513 tau = s;
1514
1515 }
1516 } else if (dMin == dN2) {
1517
1518
1519 tType = -5;
1520 double s = 0.25 * dMin;
1521
1522
1523 final int np = nn - 2 * pingPong;
1524 double b1 = work[np - 2];
1525 double b2 = work[np - 6];
1526 final double gam = dN2;
1527 if (work[np - 8] > b2 || work[np - 4] > b1) {
1528 return;
1529 }
1530 double a2 = (work[np - 8] / b2) * (1 + work[np - 4] / b1);
1531
1532
1533 if (end - start > 2) {
1534 b2 = work[nn - 13] / work[nn - 15];
1535 a2 = a2 + b2;
1536 for (int i4 = nn - 17; i4 >= 4 * start + 2 + pingPong; i4 -= 4) {
1537 if (b2 == 0.0) {
1538 break;
1539 }
1540 b1 = b2;
1541 if (work[i4] > work[i4 - 2]) {
1542 return;
1543 }
1544 b2 = b2 * (work[i4] / work[i4 - 2]);
1545 a2 = a2 + b2;
1546 if (100 * Math.max(b2, b1) < a2 || cnst1 < a2) {
1547 break;
1548 }
1549 }
1550 a2 = cnst3 * a2;
1551 }
1552
1553 if (a2 < cnst1) {
1554 tau = gam * (1 - Math.sqrt(a2)) / (1 + a2);
1555 } else {
1556 tau = s;
1557 }
1558
1559 } else {
1560
1561
1562 if (tType == -6) {
1563 g += 0.333 * (1 - g);
1564 } else if (tType == -18) {
1565 g = 0.25 * 0.333;
1566 } else {
1567 g = 0.25;
1568 }
1569 tau = g * dMin;
1570 tType = -6;
1571
1572 }
1573 break;
1574
1575 case 1 :
1576 if (dMin1 == dN1 && dMin2 == dN2) {
1577
1578
1579 tType = -7;
1580 double s = 0.333 * dMin1;
1581 if (work[nn - 5] > work[nn - 7]) {
1582 return;
1583 }
1584 double b1 = work[nn - 5] / work[nn - 7];
1585 double b2 = b1;
1586 if (b2 != 0.0) {
1587 for (int i4 = 4 * end - 10 + pingPong; i4 >= 4 * start + 2 + pingPong; i4 -= 4) {
1588 final double oldB1 = b1;
1589 if (work[i4] > work[i4 - 2]) {
1590 return;
1591 }
1592 b1 = b1 * (work[i4] / work[i4 - 2]);
1593 b2 = b2 + b1;
1594 if (100 * Math.max(b1, oldB1) < b2) {
1595 break;
1596 }
1597 }
1598 }
1599 b2 = Math.sqrt(cnst3 * b2);
1600 final double a2 = dMin1 / (1 + b2 * b2);
1601 final double gap2 = 0.5 * dMin2 - a2;
1602 if (gap2 > 0.0 && gap2 > b2 * a2) {
1603 tau = Math.max(s, a2 * (1 - cnst2 * a2 * (b2 / gap2) * b2));
1604 } else {
1605 tau = Math.max(s, a2 * (1 - cnst2 * b2));
1606 tType = -8;
1607 }
1608 } else {
1609
1610
1611 tau = 0.25 * dMin1;
1612 if (dMin1 == dN1) {
1613 tau = 0.5 * dMin1;
1614 }
1615 tType = -9;
1616 }
1617 break;
1618
1619 case 2 :
1620
1621
1622 if (dMin2 == dN2 && 2 * work[nn - 5] < work[nn - 7]) {
1623 tType = -10;
1624 final double s = 0.333 * dMin2;
1625 if (work[nn - 5] > work[nn - 7]) {
1626 return;
1627 }
1628 double b1 = work[nn - 5] / work[nn - 7];
1629 double b2 = b1;
1630 if (b2 != 0.0){
1631 for (int i4 = 4 * end - 9 + pingPong; i4 >= 4 * start + 2 + pingPong; i4 -= 4) {
1632 if (work[i4] > work[i4 - 2]) {
1633 return;
1634 }
1635 b1 *= work[i4] / work[i4 - 2];
1636 b2 += b1;
1637 if (100 * b1 < b2) {
1638 break;
1639 }
1640 }
1641 }
1642 b2 = Math.sqrt(cnst3 * b2);
1643 final double a2 = dMin2 / (1 + b2 * b2);
1644 final double gap2 = work[nn - 7] + work[nn - 9] -
1645 Math.sqrt(work[nn - 11]) * Math.sqrt(work[nn - 9]) - a2;
1646 if (gap2 > 0.0 && gap2 > b2 * a2) {
1647 tau = Math.max(s, a2 * (1 - cnst2 * a2 * (b2 / gap2) * b2));
1648 } else {
1649 tau = Math.max(s, a2 * (1 - cnst2 * b2));
1650 }
1651 } else {
1652 tau = 0.25 * dMin2;
1653 tType = -11;
1654 }
1655 break;
1656
1657 default :
1658 tau = 0.0;
1659 tType = -12;
1660 }
1661
1662 }
1663
1664
1665
1666
1667
1668 private void updateSigma(final double tau) {
1669
1670
1671
1672 if (tau < sigma) {
1673 sigmaLow += tau;
1674 final double t = sigma + sigmaLow;
1675 sigmaLow -= t - sigma;
1676 sigma = t;
1677 } else {
1678 final double t = sigma + tau;
1679 sigmaLow += sigma - (t - tau);
1680 sigma = t;
1681 }
1682 }
1683
1684
1685
1686
1687 private void findEigenVectors() {
1688
1689 final int m = main.length;
1690 eigenvectors = new ArrayRealVector[m];
1691
1692
1693 final double[] d = new double[m];
1694 final double[] l = new double[m - 1];
1695 double di = main[0];
1696 d[0] = di;
1697 for (int i = 1; i < m; ++i) {
1698 final double eiM1 = secondary[i - 1];
1699 final double ratio = eiM1 / di;
1700 di = main[i] - eiM1 * ratio;
1701 l[i - 1] = ratio;
1702 d[i] = di;
1703 }
1704
1705
1706 for (int i = 0; i < m; ++i) {
1707 eigenvectors[i] = findEigenvector(realEigenvalues[i], d, l);
1708 }
1709
1710 }
1711
1712
1713
1714
1715
1716
1717
1718
1719
1720
1721 private ArrayRealVector findEigenvector(final double eigenvalue,
1722 final double[] d, final double[] l) {
1723
1724
1725
1726 final int m = main.length;
1727 stationaryQuotientDifferenceWithShift(d, l, eigenvalue);
1728 progressiveQuotientDifferenceWithShift(d, l, eigenvalue);
1729
1730
1731
1732 int r = m - 1;
1733 double minG = Math.abs(work[6 * r] + work[6 * r + 3] + eigenvalue);
1734 for (int i = 0, sixI = 0; i < m - 1; ++i, sixI += 6) {
1735 final double g = work[sixI] + d[i] * work[sixI + 9] / work[sixI + 10];
1736 final double absG = Math.abs(g);
1737 if (absG < minG) {
1738 r = i;
1739 minG = absG;
1740 }
1741 }
1742
1743
1744
1745 double[] eigenvector = new double[m];
1746 double n2 = 1;
1747 eigenvector[r] = 1;
1748 double z = 1;
1749 for (int i = r - 1; i >= 0; --i) {
1750 z *= -work[6 * i + 2];
1751 eigenvector[i] = z;
1752 n2 += z * z;
1753 }
1754 z = 1;
1755 for (int i = r + 1; i < m; ++i) {
1756 z *= -work[6 * i - 1];
1757 eigenvector[i] = z;
1758 n2 += z * z;
1759 }
1760
1761
1762 final double inv = 1.0 / Math.sqrt(n2);
1763 for (int i = 0; i < m; ++i) {
1764 eigenvector[i] *= inv;
1765 }
1766
1767 return (transformer == null) ?
1768 new ArrayRealVector(eigenvector, false) :
1769 new ArrayRealVector(transformer.getQ().operate(eigenvector), false);
1770
1771 }
1772
1773
1774
1775
1776
1777
1778
1779
1780
1781 private void stationaryQuotientDifferenceWithShift(final double[] d, final double[] l,
1782 final double lambda) {
1783 final int nM1 = d.length - 1;
1784 double si = -lambda;
1785 for (int i = 0, sixI = 0; i < nM1; ++i, sixI += 6) {
1786 final double di = d[i];
1787 final double li = l[i];
1788 final double diP1 = di + si;
1789 final double liP1 = li * di / diP1;
1790 work[sixI] = si;
1791 work[sixI + 1] = diP1;
1792 work[sixI + 2] = liP1;
1793 si = li * liP1 * si - lambda;
1794 }
1795 work[6 * nM1 + 1] = d[nM1] + si;
1796 work[6 * nM1] = si;
1797 }
1798
1799
1800
1801
1802
1803
1804
1805
1806
1807 private void progressiveQuotientDifferenceWithShift(final double[] d, final double[] l,
1808 final double lambda) {
1809 final int nM1 = d.length - 1;
1810 double pi = d[nM1] - lambda;
1811 for (int i = nM1 - 1, sixI = 6 * i; i >= 0; --i, sixI -= 6) {
1812 final double di = d[i];
1813 final double li = l[i];
1814 final double diP1 = di * li * li + pi;
1815 final double t = di / diP1;
1816 work[sixI + 9] = pi;
1817 work[sixI + 10] = diP1;
1818 work[sixI + 5] = li * t;
1819 pi = pi * t - lambda;
1820 }
1821 work[3] = pi;
1822 work[4] = pi;
1823 }
1824
1825 }