View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      https://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  
23  package org.hipparchus.linear;
24  
25  import static org.junit.Assert.assertNotNull;
26  
27  import java.util.Arrays;
28  import java.util.Random;
29  
30  import org.hipparchus.UnitTestUtils;
31  import org.hipparchus.complex.Complex;
32  import org.hipparchus.complex.ComplexComparator;
33  import org.hipparchus.complex.ComplexField;
34  import org.hipparchus.random.RandomDataGenerator;
35  import org.hipparchus.random.RandomGenerator;
36  import org.hipparchus.random.Well1024a;
37  import org.hipparchus.util.FastMath;
38  import org.hipparchus.util.Precision;
39  import org.junit.After;
40  import org.junit.Assert;
41  import org.junit.Before;
42  import org.junit.Test;
43  
44  public class EigenDecompositionNonSymmetricTest {
45  
46      private Complex[] refValues;
47      private RealMatrix matrix;
48  
49      /** test dimensions */
50      @Test
51      public void testDimensions() {
52          final int m = matrix.getRowDimension();
53          EigenDecompositionNonSymmetric ed = new EigenDecompositionNonSymmetric(matrix);
54          Assert.assertEquals(m, ed.getV().getRowDimension());
55          Assert.assertEquals(m, ed.getV().getColumnDimension());
56          Assert.assertEquals(m, ed.getD().getColumnDimension());
57          Assert.assertEquals(m, ed.getD().getColumnDimension());
58          Assert.assertEquals(m, ed.getVInv().getRowDimension());
59          Assert.assertEquals(m, ed.getVInv().getColumnDimension());
60      }
61  
62      /** test eigenvalues */
63      @Test
64      public void testEigenvalues() {
65          EigenDecompositionNonSymmetric ed = new EigenDecompositionNonSymmetric(matrix);
66          Complex[] eigenValues = ed.getEigenvalues();
67          Assert.assertEquals(refValues.length, eigenValues.length);
68          for (int i = 0; i < refValues.length; ++i) {
69              Assert.assertEquals(refValues[i].getRealPart(),      eigenValues[i].getRealPart(), 3.0e-15);
70              Assert.assertEquals(refValues[i].getImaginaryPart(), eigenValues[i].getImaginaryPart(), 3.0e-15);
71          }
72      }
73  
74      @Test
75      public void testNonSymmetric() {
76          // Vandermonde matrix V(x;i,j) = x_i^{n - j} with x = (-1,-2,3,4)
77          double[][] vData = { { -1.0, 1.0, -1.0, 1.0 },
78                               { -8.0, 4.0, -2.0, 1.0 },
79                               { 27.0, 9.0,  3.0, 1.0 },
80                               { 64.0, 16.0, 4.0, 1.0 } };
81          checkNonSymmetricMatrix(MatrixUtils.createRealMatrix(vData));
82  
83          RealMatrix randMatrix = MatrixUtils.createRealMatrix(new double[][] {
84                  {0,  1,     0,     0},
85                  {1,  0,     2.e-7, 0},
86                  {0, -2.e-7, 0,     1},
87                  {0,  0,     1,     0}
88          });
89          checkNonSymmetricMatrix(randMatrix);
90  
91          // from http://eigen.tuxfamily.org/dox/classEigen_1_1RealSchur.html
92          double[][] randData2 = {
93                  {  0.680, -0.3300, -0.2700, -0.717, -0.687,  0.0259 },
94                  { -0.211,  0.5360,  0.0268,  0.214, -0.198,  0.6780 },
95                  {  0.566, -0.4440,  0.9040, -0.967, -0.740,  0.2250 },
96                  {  0.597,  0.1080,  0.8320, -0.514, -0.782, -0.4080 },
97                  {  0.823, -0.0452,  0.2710, -0.726,  0.998,  0.2750 },
98                  { -0.605,  0.2580,  0.4350,  0.608, -0.563,  0.0486 }
99          };
100         checkNonSymmetricMatrix(MatrixUtils.createRealMatrix(randData2));
101     }
102 
103     @Test
104     public void testRandomNonSymmetricMatrix() {
105         for (int run = 0; run < 100; run++) {
106             RandomGenerator r = new Well1024a(0x171956baefeac83el);
107 
108             // matrix size
109             int size = r.nextInt(20) + 4;
110 
111             double[][] data = new double[size][size];
112             for (int i = 0; i < size; i++) {
113                 for (int j = 0; j < size; j++) {
114                     data[i][j] = r.nextInt(100);
115                 }
116             }
117 
118             RealMatrix m = MatrixUtils.createRealMatrix(data);
119             checkNonSymmetricMatrix(m);
120         }
121     }
122 
123     /**
124      * Tests the porting of a bugfix in Jama-1.0.3 (from changelog):
125      *
126      *  Patched hqr2 method in Jama.EigenvalueDecomposition to avoid infinite loop;
127      *  Thanks Frederic Devernay &lt;frederic.devernay@m4x.org&gt;
128      */
129     @Test
130     public void testMath1051() {
131         double[][] data = {
132                 {0,0,0,0,0},
133                 {0,0,0,0,1},
134                 {0,0,0,1,0},
135                 {1,1,0,0,1},
136                 {1,0,1,0,1}
137         };
138 
139         RealMatrix m = MatrixUtils.createRealMatrix(data);
140         checkNonSymmetricMatrix(m);
141     }
142 
143     @Test
144     public void testNormalDistributionNonSymmetricMatrix() {
145         for (int run = 0; run < 100; run++) {
146             final RandomGenerator r = new Well1024a(0x511d8551a5641ea2l);
147             final RandomDataGenerator gen = new RandomDataGenerator(100);
148 
149             // matrix size
150             int size = r.nextInt(20) + 4;
151 
152             double[][] data = new double[size][size];
153             for (int i = 0; i < size; i++) {
154                 for (int j = 0; j < size; j++) {
155                     data[i][j] = gen.nextNormal(0.0, r.nextDouble() * 5);
156                 }
157             }
158 
159             RealMatrix m = MatrixUtils.createRealMatrix(data);
160             checkNonSymmetricMatrix(m);
161         }
162     }
163 
164     @Test
165     public void testMath848() {
166         double[][] data = {
167                 { 0.1849449280, -0.0646971046,  0.0774755812, -0.0969651755, -0.0692648806,  0.3282344352, -0.0177423074,  0.2063136340},
168                 {-0.0742700134, -0.0289063030, -0.0017269460, -0.0375550146, -0.0487737922, -0.2616837868, -0.0821201295, -0.2530000167},
169                 { 0.2549910127,  0.0995733692, -0.0009718388,  0.0149282808,  0.1791878897, -0.0823182816,  0.0582629256,  0.3219545182},
170                 {-0.0694747557, -0.1880649148, -0.2740630911,  0.0720096468, -0.1800836914, -0.3518996425,  0.2486747833,  0.6257938167},
171                 { 0.0536360918, -0.1339297778,  0.2241579764, -0.0195327484, -0.0054103808,  0.0347564518,  0.5120802482, -0.0329902864},
172                 {-0.5933332356, -0.2488721082,  0.2357173629,  0.0177285473,  0.0856630593, -0.3567126300, -0.1600668126, -0.1010899621},
173                 {-0.0514349819, -0.0854319435,  0.1125050061,  0.0063453560, -0.2250000688, -0.2209343090,  0.1964623477, -0.1512329924},
174                 { 0.0197395947, -0.1997170581, -0.1425959019, -0.2749477910, -0.0969467073,  0.0603688520, -0.2826905192,  0.1794315473}};
175         RealMatrix m = MatrixUtils.createRealMatrix(data);
176         checkNonSymmetricMatrix(m);
177     }
178 
179     /**
180      * Checks that the eigen decomposition of a general (non-symmetric) matrix is valid by
181      * checking: A*V = V*D
182      */
183     private void checkNonSymmetricMatrix(final RealMatrix m) {
184         try {
185             EigenDecompositionNonSymmetric ed = new EigenDecompositionNonSymmetric(m, 1.0e-15);
186 
187             RealMatrix d = ed.getD();
188             RealMatrix v = ed.getV();
189 
190             RealMatrix x = m.multiply(v);
191             RealMatrix y = v.multiply(d);
192 
193             double diffNorm = x.subtract(y).getNorm1();
194             Assert.assertTrue("The norm of (X-Y) is too large: " + diffNorm + ", matrix=" + m.toString(),
195                     x.subtract(y).getNorm1() < 1000 * Precision.EPSILON * FastMath.max(x.getNorm1(), y.getNorm1()));
196 
197             RealMatrix invV = new LUDecomposition(v).getSolver().getInverse();
198             double norm = v.multiply(d).multiply(invV).subtract(m).getNorm1();
199             Assert.assertEquals(0.0, norm, 1.0e-10);
200         } catch (Exception e) {
201             Assert.fail("Failed to create EigenDecomposition for matrix " + m.toString() + ", ex=" + e.toString());
202         }
203     }
204 
205     /** test eigenvectors */
206     @Test
207     public void testEigenvectors() {
208         EigenDecompositionNonSymmetric ed = new EigenDecompositionNonSymmetric(matrix);
209         FieldMatrix<Complex> cMatrix = MatrixUtils.createFieldMatrix(ComplexField.getInstance(),
210                                                                      matrix.getRowDimension(),
211                                                                      matrix.getColumnDimension());
212         for (int i = 0; i < matrix.getRowDimension(); ++i) {
213             for (int j = 0; j < matrix.getColumnDimension(); ++j) {
214                 cMatrix.setEntry(i, j, new Complex(matrix.getEntry(j, i)));
215             }
216         }
217 
218         for (int i = 0; i < matrix.getRowDimension(); ++i) {
219             final Complex              lambda = ed.getEigenvalue(i);
220             final FieldVector<Complex> v      = ed.getEigenvector(i);
221             final FieldVector<Complex> mV     = cMatrix.operate(v);
222             for (int k = 0; k < v.getDimension(); ++k) {
223                 Assert.assertEquals(0, mV.getEntry(k).subtract(v.getEntry(k).multiply(lambda)).norm(), 1.0e-13);
224             }
225         }
226     }
227 
228     /** test A = VDV⁻¹ */
229     @Test
230     public void testAEqualVDVInv() {
231         EigenDecompositionNonSymmetric ed = new EigenDecompositionNonSymmetric(matrix);
232         RealMatrix v  = ed.getV();
233         RealMatrix d  = ed.getD();
234         RealMatrix vI = MatrixUtils.inverse(v);
235         double norm = v.multiply(d).multiply(vI).subtract(matrix).getNorm1();
236         Assert.assertEquals(0, norm, 6.0e-13);
237     }
238 
239     /**
240      * Matrix with eigenvalues {2, 0, 12}
241      */
242     @Test
243     public void testDistinctEigenvalues() {
244         RealMatrix distinct = MatrixUtils.createRealMatrix(new double[][] {
245                 {3, 1, -4},
246                 {1, 3, -4},
247                 {-4, -4, 8}
248         });
249         EigenDecompositionNonSymmetric ed = new EigenDecompositionNonSymmetric(distinct);
250         checkEigenValues((new Complex[] { new Complex( 2), new Complex( 0), new Complex(12) }), ed, 1E-12);
251         checkEigenVector((new Complex[] { new Complex( 1), new Complex(-1), new Complex( 0) }), ed, 1E-12);
252         checkEigenVector((new Complex[] { new Complex( 1), new Complex( 1), new Complex( 1) }), ed, 1E-12);
253         checkEigenVector((new Complex[] { new Complex(-1), new Complex(-1), new Complex( 2) }), ed, 1E-12);
254     }
255 
256     /**
257      * Verifies operation on indefinite matrix
258      */
259     @Test
260     public void testZeroDivide() {
261         RealMatrix indefinite = MatrixUtils.createRealMatrix(new double [][] {
262                 { 0.0, 1.0, -1.0 },
263                 { 1.0, 1.0, 0.0 },
264                 { -1.0,0.0, 1.0 }
265         });
266         EigenDecompositionNonSymmetric ed = new EigenDecompositionNonSymmetric(indefinite);
267         checkEigenValues((new Complex[] { new Complex(2), new Complex(1), new Complex(-1) }), ed, 1E-12);
268         Complex isqrt3 = new Complex(1 / FastMath.sqrt(3.0));
269         checkEigenVector((new Complex[] {isqrt3, isqrt3, isqrt3.negate() }), ed, 1E-12);
270         Complex isqrt2 = new Complex(1 / FastMath.sqrt(2.0));
271         checkEigenVector((new Complex[] { new Complex(0.0), isqrt2.negate(), isqrt2.negate() }), ed, 1E-12);
272         Complex isqrt6 = new Complex(1 / FastMath.sqrt(6.0));
273         checkEigenVector((new Complex[] { isqrt6.multiply(2), isqrt6.negate(), isqrt6 }), ed, 1E-12);
274     }
275 
276     /** test eigenvalues for a big matrix. */
277     @Test
278     public void testBigMatrix() {
279         Random r = new Random(17748333525117l);
280         double[] bigValues = new double[200];
281         for (int i = 0; i < bigValues.length; ++i) {
282             bigValues[i] = 2 * r.nextDouble() - 1;
283         }
284         Arrays.sort(bigValues);
285         EigenDecompositionNonSymmetric ed = new EigenDecompositionNonSymmetric(createTestMatrix(r, bigValues));
286         Complex[] eigenValues = ed.getEigenvalues();
287         Arrays.sort(eigenValues, new ComplexComparator());
288         Assert.assertEquals(bigValues.length, eigenValues.length);
289         for (int i = 0; i < bigValues.length; ++i) {
290             Assert.assertEquals(bigValues[i], eigenValues[i].getRealPart(), 2.0e-14);
291         }
292     }
293 
294     /**
295      * Verifies operation on very small values.
296      * Matrix with eigenvalues {2e-100, 0, 12e-100}
297      */
298     @Test
299     public void testTinyValues() {
300         final double tiny = 1.0e-100;
301         RealMatrix distinct = MatrixUtils.createRealMatrix(new double[][] {
302                 {3, 1, -4},
303                 {1, 3, -4},
304                 {-4, -4, 8}
305         });
306         distinct = distinct.scalarMultiply(tiny);
307 
308         final EigenDecompositionNonSymmetric ed = new EigenDecompositionNonSymmetric(distinct);
309         checkEigenValues(new Complex[] { new Complex(2).multiply(tiny), new Complex(0).multiply(tiny), new Complex(12).multiply(tiny) },
310                          ed, 1e-12 * tiny);
311         checkEigenVector(new Complex[] { new Complex( 1), new Complex(-1), new Complex(0) }, ed, 1e-12);
312         checkEigenVector(new Complex[] { new Complex( 1), new Complex( 1), new Complex(1) }, ed, 1e-12);
313         checkEigenVector(new Complex[] { new Complex(-1), new Complex(-1), new Complex(2) }, ed, 1e-12);
314     }
315 
316     @Test
317     public void testDeterminantWithCompleEigenValues() {
318         final RealMatrix m = MatrixUtils.createRealMatrix(new double[][] { { -3, -1.5, -3 }, { 0, -1, 0 }, { 1, 0, 0 } });
319         EigenDecompositionNonSymmetric decomposition = new EigenDecompositionNonSymmetric(m);
320         Assert.assertEquals(-3.0, decomposition.getDeterminant().getRealPart(),      1.0e-15);
321         Assert.assertEquals( 0.0, decomposition.getDeterminant().getImaginaryPart(), 1.0e-15);
322     }
323 
324     @Test
325     public void testReal() {
326         // AA = [1 2;1 -3];
327 
328         RealMatrix A = MatrixUtils.createRealMatrix(new double[][] {
329             { 1, 2 }, { 1, -3 }
330         });
331 
332         EigenDecompositionNonSymmetric ed = new EigenDecompositionNonSymmetric(A);
333 
334         assertNotNull(ed.getD());
335         assertNotNull(ed.getV());
336 
337         final double s2 = FastMath.sqrt(2);
338         final double s3 = FastMath.sqrt(3);
339         RealMatrix D_expected = MatrixUtils.createRealMatrix(new double[][] {
340             { s2 * s3 - 1.0, 0 }, { 0,  -1 - s2 * s3 }
341         });
342 
343         RealMatrix V_expected = MatrixUtils.createRealMatrix(new double[][] {
344             { (s2 + 2 * s3) / 5, (s3 - 3 * s2)     /  6 },
345             { (2 * s2 - s3) / 5, (4 * s3 + 3 * s2) / 12 }
346         });
347 
348         UnitTestUtils.assertEquals("D", D_expected, ed.getD(), 1.0e-15);
349         UnitTestUtils.assertEquals("V", V_expected, ed.getV(), 1.0e-15);
350 
351         // checking definition of the decomposition A = V*D*inv(V)
352         UnitTestUtils.assertEquals("A", A,
353                                    ed.getV().multiply(ed.getD()).multiply(MatrixUtils.inverse(ed.getV())),
354                                    2.0e-15);
355 
356     }
357 
358     @Test
359     public void testImaginary() {
360         // AA = [3 -2;4 -1];
361 
362         RealMatrix A = MatrixUtils.createRealMatrix(new double[][] {
363             { 3, -2 }, { 4, -1 }
364         });
365 
366         EigenDecompositionNonSymmetric ordEig = new EigenDecompositionNonSymmetric(A);
367 
368         assertNotNull(ordEig.getD());
369         assertNotNull(ordEig.getV());
370 
371         RealMatrix D_expected = MatrixUtils.createRealMatrix(new double[][] {
372             { 1, 2 }, { -2, 1 }
373         });
374 
375         RealMatrix V_expected = MatrixUtils.createRealMatrix(new double[][] {
376             { -0.5, 0.5 }, { 0, 1 }
377         });
378 
379         UnitTestUtils.assertEquals("D", D_expected, ordEig.getD(), 1.0e-15);
380         UnitTestUtils.assertEquals("V", V_expected, ordEig.getV(), 1.0e-15);
381 
382         // checking definition of the decomposition A = V*D*inv(V)
383         UnitTestUtils.assertEquals("A", A, ordEig.getV().multiply(ordEig.getD()).multiply(MatrixUtils.inverse(ordEig.getV())), 1.0e-15);
384 
385         // checking A*V = D*V
386         FieldMatrix<Complex> ac = new Array2DRowFieldMatrix<>(ComplexField.getInstance(), 2, 2);
387         for (int i = 0; i < ac.getRowDimension(); ++i) {
388             for (int j = 0; j < ac.getColumnDimension(); ++j) {
389                 ac.setEntry(i, j, new Complex(A.getEntry(i, j)));
390             }
391         }
392         for (int i = 0; i < ordEig.getEigenvalues().length; ++i) {
393             final Complex              li  = ordEig.getEigenvalue(i);
394             final FieldVector<Complex> ei  = ordEig.getEigenvector(i);
395             final FieldVector<Complex> aei = ac.operate(ei);
396             for (int j = 0; j < ei.getDimension(); ++j) {
397                 Assert.assertEquals(aei.getEntry(j).getRealPart(),      li.multiply(ei.getEntry(j)).getRealPart(),      1.0e-10);
398                 Assert.assertEquals(aei.getEntry(j).getImaginaryPart(), li.multiply(ei.getEntry(j)).getImaginaryPart(), 1.0e-10);
399             }
400         }
401 
402     }
403 
404     @Test
405     public void testImaginary33() {
406         // AA = [3 -2 0;4 -1 0;1 1 1];
407 
408         RealMatrix A = MatrixUtils.createRealMatrix(new double[][] {
409             { 3, -2, 0 }, { 4, -1, 0 }, { 1, 1, 1 }
410         });
411 
412         EigenDecompositionNonSymmetric ordEig = new EigenDecompositionNonSymmetric(A);
413 
414         assertNotNull(ordEig.getD());
415         assertNotNull(ordEig.getV());
416 
417         RealMatrix D_expected = MatrixUtils.createRealMatrix(new double[][] {
418             { 1, 2, 0 }, { -2, 1, 0}, {0, 0, 1 }
419         });
420 
421         final double a = FastMath.sqrt(8.0 / 17.0);
422         final double b = FastMath.sqrt(17.0) / 2.0;
423         RealMatrix V_expected = MatrixUtils.createRealMatrix(new double[][] {
424             { 0,        a, 0 },
425             { a,        a, 0 },
426             { a, -0.5 * a, b }
427         });
428         UnitTestUtils.assertEquals("D", D_expected, ordEig.getD(), 1.0e-15);
429         UnitTestUtils.assertEquals("V", V_expected, ordEig.getV(), 1.0e-15);
430 
431         // checking definition of the decomposition A = V*D*inv(V)
432         UnitTestUtils.assertEquals("A", A,
433                                    ordEig.getV().multiply(ordEig.getD()).multiply(MatrixUtils.inverse(ordEig.getV())),
434                                    8.0e-15);
435 
436         // checking A*V = D*V
437         FieldMatrix<Complex> ac = new Array2DRowFieldMatrix<>(ComplexField.getInstance(), 3, 3);
438         for (int i = 0; i < ac.getRowDimension(); ++i) {
439             for (int j = 0; j < ac.getColumnDimension(); ++j) {
440                 ac.setEntry(i, j, new Complex(A.getEntry(i, j)));
441             }
442         }
443         for (int i = 0; i < ordEig.getEigenvalues().length; ++i) {
444             final Complex              li  = ordEig.getEigenvalue(i);
445             final FieldVector<Complex> ei  = ordEig.getEigenvector(i);
446             final FieldVector<Complex> aei = ac.operate(ei);
447             for (int j = 0; j < ei.getDimension(); ++j) {
448                 Assert.assertEquals(aei.getEntry(j).getRealPart(),      li.multiply(ei.getEntry(j)).getRealPart(),      1.0e-10);
449                 Assert.assertEquals(aei.getEntry(j).getImaginaryPart(), li.multiply(ei.getEntry(j)).getImaginaryPart(), 1.0e-10);
450             }
451         }
452 
453     }
454 
455     @Test
456     public void testImaginaryNullEigenvalue() {
457         // AA = [3 -2 0;4 -1 0;3 -2 0];
458 
459         RealMatrix A = MatrixUtils.createRealMatrix(new double[][] {
460             { 3, -2, 0 }, { 4, -1, 0 }, { 3, -2, 0 }
461         });
462 
463         EigenDecompositionNonSymmetric ordEig = new EigenDecompositionNonSymmetric(A);
464 
465         assertNotNull(ordEig.getD());
466         assertNotNull(ordEig.getV());
467 
468         RealMatrix D_expected = MatrixUtils.createRealMatrix(new double[][] {
469             { 1, 2, 0 }, { -2, 1, 0 }, { 0, 0, 0 }
470         });
471 
472         final double a  = FastMath.sqrt(11.0 / 50.0);
473         final double s2 = FastMath.sqrt(2.0);
474         RealMatrix V_expected = MatrixUtils.createRealMatrix(new double[][] {
475             {      -a, -2 * a,  0.0 },
476             {  -3 * a,     -a,  0.0 },
477             {      -a, -2 * a,   s2 }
478         });
479 
480         UnitTestUtils.assertEquals("D", D_expected, ordEig.getD(), 1.0e-15);
481         UnitTestUtils.assertEquals("V", V_expected, ordEig.getV(), 2.0e-15);
482 
483         // checking definition of the decomposition A = V*D*inv(V)
484         UnitTestUtils.assertEquals("A", A,
485                                    ordEig.getV().multiply(ordEig.getD()).multiply(MatrixUtils.inverse(ordEig.getV())),
486                                    1.0e-14);
487 
488         // checking A*V = D*V
489         FieldMatrix<Complex> ac = new Array2DRowFieldMatrix<>(ComplexField.getInstance(), 3, 3);
490         for (int i = 0; i < ac.getRowDimension(); ++i) {
491             for (int j = 0; j < ac.getColumnDimension(); ++j) {
492                 ac.setEntry(i, j, new Complex(A.getEntry(i, j)));
493             }
494         }
495 
496         for (int i = 0; i < ordEig.getEigenvalues().length; ++i) {
497             final Complex              li  = ordEig.getEigenvalue(i);
498             final FieldVector<Complex> ei  = ordEig.getEigenvector(i);
499             final FieldVector<Complex> aei = ac.operate(ei);
500             for (int j = 0; j < ei.getDimension(); ++j) {
501                 Assert.assertEquals(aei.getEntry(j).getRealPart(),      li.multiply(ei.getEntry(j)).getRealPart(),      1.0e-10);
502                 Assert.assertEquals(aei.getEntry(j).getImaginaryPart(), li.multiply(ei.getEntry(j)).getImaginaryPart(), 1.0e-10);
503             }
504         }
505 
506     }
507 
508     /**
509      * Verifies that the given EigenDecomposition has eigenvalues equivalent to
510      * the targetValues, ignoring the order of the values and allowing
511      * values to differ by tolerance.
512      */
513     protected void checkEigenValues(Complex[] targetValues,
514                                     EigenDecompositionNonSymmetric ed, double tolerance) {
515         Complex[] observed = ed.getEigenvalues();
516         for (int i = 0; i < observed.length; i++) {
517             Assert.assertTrue(isIncludedValue(observed[i], targetValues, tolerance));
518             Assert.assertTrue(isIncludedValue(targetValues[i], observed, tolerance));
519         }
520     }
521 
522 
523     /**
524      * Returns true iff there is an entry within tolerance of value in
525      * searchArray.
526      */
527     private boolean isIncludedValue(Complex value, Complex[] searchArray, double tolerance) {
528        boolean found = false;
529        int i = 0;
530        while (!found && i < searchArray.length) {
531            if (value.subtract(searchArray[i]).norm() < tolerance) {
532                found = true;
533            }
534            i++;
535        }
536        return found;
537     }
538 
539     /**
540      * Returns true iff eigenVector is a scalar multiple of one of the columns
541      * of ed.getV().  Does not try linear combinations - i.e., should only be
542      * used to find vectors in one-dimensional eigenspaces.
543      */
544     protected void checkEigenVector(Complex[] eigenVector, EigenDecompositionNonSymmetric ed, double tolerance) {
545         Assert.assertTrue(isIncludedColumn(eigenVector, ed.getV(), tolerance));
546     }
547 
548     /**
549      * Returns true iff there is a column that is a scalar multiple of column
550      * in searchMatrix (modulo tolerance)
551      */
552     private boolean isIncludedColumn(Complex[] column, RealMatrix searchMatrix, double tolerance) {
553         boolean found = false;
554         int i = 0;
555         while (!found && i < searchMatrix.getColumnDimension()) {
556             Complex multiplier = Complex.ONE;
557             boolean matching = true;
558             int j = 0;
559             while (matching && j < searchMatrix.getRowDimension()) {
560                 double colEntry = searchMatrix.getEntry(j, i);
561                 // Use the first entry where both are non-zero as scalar
562                 if (multiplier.subtract(1).norm() <= FastMath.ulp(1.0) && FastMath.abs(colEntry) > 1E-14
563                         && column[j].norm() > 1e-14) {
564                     multiplier = column[j].reciprocal().multiply(colEntry);
565                 }
566                 if (column[j].multiply(multiplier).subtract(colEntry).norm() > tolerance) {
567                     matching = false;
568                 }
569                 j++;
570             }
571             found = matching;
572             i++;
573         }
574         return found;
575     }
576 
577     @Before
578     public void setUp() {
579         double[] real = {
580         		0.001, 1.000, 1.001, 2.001, 2.002, 2.003
581         };
582         refValues = new Complex[real.length];
583         for (int i = 0; i < refValues.length; ++i) {
584             refValues[i] = new Complex(real[i]);
585         }
586         matrix = createTestMatrix(new Random(35992629946426l), real);
587     }
588 
589     @After
590     public void tearDown() {
591         refValues = null;
592         matrix    = null;
593     }
594 
595     static RealMatrix createTestMatrix(final Random r, final double[] eigenValues) {
596         final int n = eigenValues.length;
597         final RealMatrix v = createOrthogonalMatrix(r, n);
598         final RealMatrix d = MatrixUtils.createRealDiagonalMatrix(eigenValues);
599         return v.multiply(d).multiplyTransposed(v);
600     }
601 
602     public static RealMatrix createOrthogonalMatrix(final Random r, final int size) {
603 
604         final double[][] data = new double[size][size];
605 
606         for (int i = 0; i < size; ++i) {
607             final double[] dataI = data[i];
608             double norm2 = 0;
609             do {
610 
611                 // generate randomly row I
612                 for (int j = 0; j < size; ++j) {
613                     dataI[j] = 2 * r.nextDouble() - 1;
614                 }
615 
616                 // project the row in the subspace orthogonal to previous rows
617                 for (int k = 0; k < i; ++k) {
618                     final double[] dataK = data[k];
619                     double dotProduct = 0;
620                     for (int j = 0; j < size; ++j) {
621                         dotProduct += dataI[j] * dataK[j];
622                     }
623                     for (int j = 0; j < size; ++j) {
624                         dataI[j] -= dotProduct * dataK[j];
625                     }
626                 }
627 
628                 // normalize the row
629                 norm2 = 0;
630                 for (final double dataIJ : dataI) {
631                     norm2 += dataIJ * dataIJ;
632                 }
633                 final double inv = 1.0 / FastMath.sqrt(norm2);
634                 for (int j = 0; j < size; ++j) {
635                     dataI[j] *= inv;
636                 }
637 
638             } while (norm2 * size < 0.01);
639         }
640 
641         return MatrixUtils.createRealMatrix(data);
642 
643     }
644 }