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  package org.hipparchus.linear;
23  
24  import java.math.BigDecimal;
25  import java.util.Arrays;
26  import java.util.List;
27  
28  import org.hipparchus.CalculusFieldElement;
29  import org.hipparchus.Field;
30  import org.hipparchus.UnitTestUtils;
31  import org.hipparchus.exception.LocalizedCoreFormats;
32  import org.hipparchus.exception.MathIllegalArgumentException;
33  import org.hipparchus.exception.MathIllegalStateException;
34  import org.hipparchus.exception.NullArgumentException;
35  import org.hipparchus.fraction.BigFraction;
36  import org.hipparchus.fraction.Fraction;
37  import org.hipparchus.fraction.FractionField;
38  import org.hipparchus.util.Binary64;
39  import org.hipparchus.util.Binary64Field;
40  import org.hipparchus.util.FastMath;
41  import org.hipparchus.util.Precision;
42  import org.junit.Assert;
43  import org.junit.Test;
44  
45  /**
46   * Test cases for the {@link MatrixUtils} class.
47   *
48   */
49  
50  public final class MatrixUtilsTest {
51  
52      protected double[][] testData = { {1d,2d,3d}, {2d,5d,3d}, {1d,0d,8d} };
53      protected double[][] testData3x3Singular = { { 1, 4, 7, }, { 2, 5, 8, }, { 3, 6, 9, } };
54      protected double[][] testData3x4 = { { 12, -51, 4, 1 }, { 6, 167, -68, 2 }, { -4, 24, -41, 3 } };
55      protected double[][] nullMatrix = null;
56      protected double[] row = {1,2,3};
57      protected BigDecimal[] bigRow =
58          {new BigDecimal(1),new BigDecimal(2),new BigDecimal(3)};
59      protected String[] stringRow = {"1", "2", "3"};
60      protected Fraction[] fractionRow =
61          {new Fraction(1),new Fraction(2),new Fraction(3)};
62      protected double[][] rowMatrix = {{1,2,3}};
63      protected BigDecimal[][] bigRowMatrix =
64          {{new BigDecimal(1), new BigDecimal(2), new BigDecimal(3)}};
65      protected String[][] stringRowMatrix = {{"1", "2", "3"}};
66      protected Fraction[][] fractionRowMatrix =
67          {{new Fraction(1), new Fraction(2), new Fraction(3)}};
68      protected double[] col = {0,4,6};
69      protected BigDecimal[] bigCol =
70          {new BigDecimal(0),new BigDecimal(4),new BigDecimal(6)};
71      protected String[] stringCol = {"0","4","6"};
72      protected Fraction[] fractionCol =
73          {new Fraction(0),new Fraction(4),new Fraction(6)};
74      protected double[] nullDoubleArray = null;
75      protected double[][] colMatrix = {{0},{4},{6}};
76      protected BigDecimal[][] bigColMatrix =
77          {{new BigDecimal(0)},{new BigDecimal(4)},{new BigDecimal(6)}};
78      protected String[][] stringColMatrix = {{"0"}, {"4"}, {"6"}};
79      protected Fraction[][] fractionColMatrix =
80          {{new Fraction(0)},{new Fraction(4)},{new Fraction(6)}};
81  
82      @Test
83      public void testCreateRealMatrix() {
84          Assert.assertEquals(new BlockRealMatrix(testData),
85                  MatrixUtils.createRealMatrix(testData));
86          try {
87              MatrixUtils.createRealMatrix(new double[][] {{1}, {1,2}});  // ragged
88              Assert.fail("Expecting MathIllegalArgumentException");
89          } catch (MathIllegalArgumentException ex) {
90              // expected
91          }
92          try {
93              MatrixUtils.createRealMatrix(new double[][] {{}, {}});  // no columns
94              Assert.fail("Expecting MathIllegalArgumentException");
95          } catch (MathIllegalArgumentException ex) {
96              // expected
97          }
98          try {
99              MatrixUtils.createRealMatrix(null);  // null
100             Assert.fail("Expecting NullArgumentException");
101         } catch (NullArgumentException ex) {
102             // expected
103         }
104     }
105 
106     @Test
107     public void testcreateFieldMatrix() {
108         Assert.assertEquals(new Array2DRowFieldMatrix<Fraction>(asFraction(testData)),
109                      MatrixUtils.createFieldMatrix(asFraction(testData)));
110         Assert.assertEquals(new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), fractionColMatrix),
111                      MatrixUtils.createFieldMatrix(fractionColMatrix));
112         try {
113             MatrixUtils.createFieldMatrix(asFraction(new double[][] {{1}, {1,2}}));  // ragged
114             Assert.fail("Expecting MathIllegalArgumentException");
115         } catch (MathIllegalArgumentException ex) {
116             // expected
117         }
118         try {
119             MatrixUtils.createFieldMatrix(asFraction(new double[][] {{}, {}}));  // no columns
120             Assert.fail("Expecting MathIllegalArgumentException");
121         } catch (MathIllegalArgumentException ex) {
122             // expected
123         }
124         try {
125             MatrixUtils.createFieldMatrix((Fraction[][])null);  // null
126             Assert.fail("Expecting NullArgumentException");
127         } catch (NullArgumentException ex) {
128             // expected
129         }
130     }
131 
132     @Test
133     public void testCreateRealVector() {
134         RealVector v1 = MatrixUtils.createRealVector(new double[] { 0.0, 1.0, 2.0, 3.0 });
135         Assert.assertEquals(4, v1.getDimension());
136         for (int i = 0; i < v1.getDimension(); ++i) {
137             Assert.assertEquals(i, v1.getEntry(i), 1.0e-15);
138         }
139         RealVector v2 = MatrixUtils.createRealVector(7);
140         Assert.assertEquals(7, v2.getDimension());
141         for (int i = 0; i < v2.getDimension(); ++i) {
142             Assert.assertEquals(0.0, v2.getEntry(i), 1.0e-15);
143         }
144     }
145 
146     @Test
147     public void testCreateFieldVector() {
148         FieldVector<Binary64> v1 = MatrixUtils.createFieldVector(new Binary64[] {
149             new Binary64(0.0), new Binary64(1.0), new Binary64(2.0), new Binary64(3.0)
150         });
151         Assert.assertEquals(4, v1.getDimension());
152         for (int i = 0; i < v1.getDimension(); ++i) {
153             Assert.assertEquals(i, v1.getEntry(i).getReal(), 1.0e-15);
154         }
155         FieldVector<Binary64> v2 = MatrixUtils.createFieldVector(Binary64Field.getInstance(), 7);
156         Assert.assertEquals(7, v2.getDimension());
157         for (int i = 0; i < v2.getDimension(); ++i) {
158             Assert.assertEquals(0.0, v2.getEntry(i).getReal(), 1.0e-15);
159         }
160     }
161 
162     @Test
163     public void testCreateRowRealMatrix() {
164         Assert.assertEquals(MatrixUtils.createRowRealMatrix(row),
165                      new BlockRealMatrix(rowMatrix));
166         try {
167             MatrixUtils.createRowRealMatrix(new double[] {});  // empty
168             Assert.fail("Expecting MathIllegalArgumentException");
169         } catch (MathIllegalArgumentException ex) {
170             // expected
171         }
172         try {
173             MatrixUtils.createRowRealMatrix(null);  // null
174             Assert.fail("Expecting NullArgumentException");
175         } catch (NullArgumentException ex) {
176             // expected
177         }
178     }
179 
180     @Test
181     public void testCreateRowFieldMatrix() {
182         Assert.assertEquals(MatrixUtils.createRowFieldMatrix(asFraction(row)),
183                      new Array2DRowFieldMatrix<Fraction>(asFraction(rowMatrix)));
184         Assert.assertEquals(MatrixUtils.createRowFieldMatrix(fractionRow),
185                      new Array2DRowFieldMatrix<Fraction>(fractionRowMatrix));
186         try {
187             MatrixUtils.createRowFieldMatrix(new Fraction[] {});  // empty
188             Assert.fail("Expecting MathIllegalArgumentException");
189         } catch (MathIllegalArgumentException ex) {
190             // expected
191         }
192         try {
193             MatrixUtils.createRowFieldMatrix((Fraction[]) null);  // null
194             Assert.fail("Expecting NullArgumentException");
195         } catch (NullArgumentException ex) {
196             // expected
197         }
198     }
199 
200     @Test
201     public void testCreateColumnRealMatrix() {
202         Assert.assertEquals(MatrixUtils.createColumnRealMatrix(col),
203                      new BlockRealMatrix(colMatrix));
204         try {
205             MatrixUtils.createColumnRealMatrix(new double[] {});  // empty
206             Assert.fail("Expecting MathIllegalArgumentException");
207         } catch (MathIllegalArgumentException ex) {
208             // expected
209         }
210         try {
211             MatrixUtils.createColumnRealMatrix(null);  // null
212             Assert.fail("Expecting NullArgumentException");
213         } catch (NullArgumentException ex) {
214             // expected
215         }
216     }
217 
218     @Test
219     public void testCreateColumnFieldMatrix() {
220         Assert.assertEquals(MatrixUtils.createColumnFieldMatrix(asFraction(col)),
221                      new Array2DRowFieldMatrix<Fraction>(asFraction(colMatrix)));
222         Assert.assertEquals(MatrixUtils.createColumnFieldMatrix(fractionCol),
223                      new Array2DRowFieldMatrix<Fraction>(fractionColMatrix));
224 
225         try {
226             MatrixUtils.createColumnFieldMatrix(new Fraction[] {});  // empty
227             Assert.fail("Expecting MathIllegalArgumentException");
228         } catch (MathIllegalArgumentException ex) {
229             // expected
230         }
231         try {
232             MatrixUtils.createColumnFieldMatrix((Fraction[]) null);  // null
233             Assert.fail("Expecting NullArgumentException");
234         } catch (NullArgumentException ex) {
235             // expected
236         }
237     }
238 
239     /**
240      * Verifies that the matrix is an identity matrix
241      */
242     protected void checkIdentityMatrix(RealMatrix m) {
243         for (int i = 0; i < m.getRowDimension(); i++) {
244             for (int j =0; j < m.getColumnDimension(); j++) {
245                 if (i == j) {
246                     Assert.assertEquals(m.getEntry(i, j), 1d, 0);
247                 } else {
248                     Assert.assertEquals(m.getEntry(i, j), 0d, 0);
249                 }
250             }
251         }
252     }
253 
254     @Test
255     public void testCreateIdentityMatrix() {
256         checkIdentityMatrix(MatrixUtils.createRealIdentityMatrix(3));
257         checkIdentityMatrix(MatrixUtils.createRealIdentityMatrix(2));
258         checkIdentityMatrix(MatrixUtils.createRealIdentityMatrix(1));
259         try {
260             MatrixUtils.createRealIdentityMatrix(0);
261             Assert.fail("Expecting MathIllegalArgumentException");
262         } catch (MathIllegalArgumentException ex) {
263             // expected
264         }
265     }
266 
267     /**
268      * Verifies that the matrix is an identity matrix
269      */
270     protected void checkIdentityFieldMatrix(FieldMatrix<Fraction> m) {
271         for (int i = 0; i < m.getRowDimension(); i++) {
272             for (int j =0; j < m.getColumnDimension(); j++) {
273                 if (i == j) {
274                     Assert.assertEquals(m.getEntry(i, j), Fraction.ONE);
275                 } else {
276                     Assert.assertEquals(m.getEntry(i, j), Fraction.ZERO);
277                 }
278             }
279         }
280     }
281 
282     @Test
283     public void testcreateFieldIdentityMatrix() {
284         checkIdentityFieldMatrix(MatrixUtils.createFieldIdentityMatrix(FractionField.getInstance(), 3));
285         checkIdentityFieldMatrix(MatrixUtils.createFieldIdentityMatrix(FractionField.getInstance(), 2));
286         checkIdentityFieldMatrix(MatrixUtils.createFieldIdentityMatrix(FractionField.getInstance(), 1));
287         try {
288             MatrixUtils.createRealIdentityMatrix(0);
289             Assert.fail("Expecting MathIllegalArgumentException");
290         } catch (MathIllegalArgumentException ex) {
291             // expected
292         }
293     }
294 
295     @Test
296     public void testBigFractionConverter() {
297         BigFraction[][] bfData = {
298                 { new BigFraction(1), new BigFraction(2), new BigFraction(3) },
299                 { new BigFraction(2), new BigFraction(5), new BigFraction(3) },
300                 { new BigFraction(1), new BigFraction(0), new BigFraction(8) }
301         };
302         FieldMatrix<BigFraction> m = new Array2DRowFieldMatrix<BigFraction>(bfData, false);
303         RealMatrix converted = MatrixUtils.bigFractionMatrixToRealMatrix(m);
304         RealMatrix reference = new Array2DRowRealMatrix(testData, false);
305         Assert.assertEquals(0.0, converted.subtract(reference).getNorm1(), 0.0);
306     }
307 
308     @Test
309     public void testFractionConverter() {
310         Fraction[][] fData = {
311                 { new Fraction(1), new Fraction(2), new Fraction(3) },
312                 { new Fraction(2), new Fraction(5), new Fraction(3) },
313                 { new Fraction(1), new Fraction(0), new Fraction(8) }
314         };
315         FieldMatrix<Fraction> m = new Array2DRowFieldMatrix<Fraction>(fData, false);
316         RealMatrix converted = MatrixUtils.fractionMatrixToRealMatrix(m);
317         RealMatrix reference = new Array2DRowRealMatrix(testData, false);
318         Assert.assertEquals(0.0, converted.subtract(reference).getNorm1(), 0.0);
319     }
320 
321     public static final Fraction[][] asFraction(double[][] data) {
322         Fraction d[][] = new Fraction[data.length][];
323         try {
324             for (int i = 0; i < data.length; ++i) {
325                 double[] dataI = data[i];
326                 Fraction[] dI  = new Fraction[dataI.length];
327                 for (int j = 0; j < dataI.length; ++j) {
328                     dI[j] = new Fraction(dataI[j]);
329                 }
330                 d[i] = dI;
331             }
332         } catch (MathIllegalStateException fce) {
333             Assert.fail(fce.getMessage());
334         }
335         return d;
336     }
337 
338     public static final Fraction[] asFraction(double[] data) {
339         Fraction d[] = new Fraction[data.length];
340         try {
341             for (int i = 0; i < data.length; ++i) {
342                 d[i] = new Fraction(data[i]);
343             }
344         } catch (MathIllegalStateException fce) {
345             Assert.fail(fce.getMessage());
346         }
347         return d;
348     }
349 
350     @Test
351     public void testSolveLowerTriangularSystem(){
352         RealMatrix rm = new Array2DRowRealMatrix(
353                 new double[][] { {2,0,0,0 }, { 1,1,0,0 }, { 3,3,3,0 }, { 3,3,3,4 } },
354                        false);
355         RealVector b = new ArrayRealVector(new double[] { 2,3,4,8 }, false);
356         MatrixUtils.solveLowerTriangularSystem(rm, b);
357         UnitTestUtils.assertEquals( new double[]{1,2,-1.66666666666667, 1.0}  , b.toArray() , 1.0e-12);
358     }
359 
360 
361     /*
362      * Taken from R manual http://stat.ethz.ch/R-manual/R-patched/library/base/html/backsolve.html
363      */
364     @Test
365     public void testSolveUpperTriangularSystem(){
366         RealMatrix rm = new Array2DRowRealMatrix(
367                 new double[][] { {1,2,3 }, { 0,1,1 }, { 0,0,2 } },
368                        false);
369         RealVector b = new ArrayRealVector(new double[] { 8,4,2 }, false);
370         MatrixUtils.solveUpperTriangularSystem(rm, b);
371         UnitTestUtils.assertEquals( new double[]{-1,3,1}  , b.toArray() , 1.0e-12);
372     }
373 
374     /**
375      * This test should probably be replaced by one that could show
376      * whether this algorithm can sometimes perform better (precision- or
377      * performance-wise) than the direct inversion of the whole matrix.
378      */
379     @Test
380     public void testBlockInverse() {
381         final double[][] data = {
382             { -1, 0, 123, 4 },
383             { -56, 78.9, -0.1, -23.4 },
384             { 5.67, 8, -9, 1011 },
385             { 12, 345, -67.8, 9 },
386         };
387 
388         final RealMatrix m = new Array2DRowRealMatrix(data);
389         final int len = data.length;
390         final double tol = 1e-14;
391 
392         for (int splitIndex = 0; splitIndex < 3; splitIndex++) {
393             final RealMatrix mInv = MatrixUtils.blockInverse(m, splitIndex);
394             final RealMatrix id = m.multiply(mInv);
395 
396             // Check that we recovered the identity matrix.
397             for (int i = 0; i < len; i++) {
398                 for (int j = 0; j < len; j++) {
399                     final double entry = id.getEntry(i, j);
400                     if (i == j) {
401                         Assert.assertEquals("[" + i + "][" + j + "]",
402                                             1, entry, tol);
403                     } else {
404                         Assert.assertEquals("[" + i + "][" + j + "]",
405                                             0, entry, tol);
406                     }
407                 }
408             }
409         }
410     }
411 
412     @Test(expected=MathIllegalArgumentException.class)
413     public void testBlockInverseNonInvertible() {
414         final double[][] data = {
415             { -1, 0, 123, 4 },
416             { -56, 78.9, -0.1, -23.4 },
417             { 5.67, 8, -9, 1011 },
418             { 5.67, 8, -9, 1011 },
419         };
420 
421         MatrixUtils.blockInverse(new Array2DRowRealMatrix(data), 2);
422     }
423 
424     @Test
425     public void testIsSymmetric() {
426         final double eps = Math.ulp(1d);
427 
428         final double[][] dataSym = {
429             { 1, 2, 3 },
430             { 2, 2, 5 },
431             { 3, 5, 6 },
432         };
433         Assert.assertTrue(MatrixUtils.isSymmetric(MatrixUtils.createRealMatrix(dataSym), eps));
434 
435         final double[][] dataNonSym = {
436             { 1, 2, -3 },
437             { 2, 2, 5 },
438             { 3, 5, 6 },
439         };
440         Assert.assertFalse(MatrixUtils.isSymmetric(MatrixUtils.createRealMatrix(dataNonSym), eps));
441     }
442 
443     @Test
444     public void testIsSymmetricTolerance() {
445         final double eps = 1e-4;
446 
447         final double[][] dataSym1 = {
448             { 1,   1, 1.00009 },
449             { 1,   1, 1       },
450             { 1.0, 1, 1       },
451         };
452         Assert.assertTrue(MatrixUtils.isSymmetric(MatrixUtils.createRealMatrix(dataSym1), eps));
453         final double[][] dataSym2 = {
454             { 1,   1, 0.99990 },
455             { 1,   1, 1       },
456             { 1.0, 1, 1       },
457         };
458         Assert.assertTrue(MatrixUtils.isSymmetric(MatrixUtils.createRealMatrix(dataSym2), eps));
459 
460         final double[][] dataNonSym1 = {
461             { 1,   1, 1.00011 },
462             { 1,   1, 1       },
463             { 1.0, 1, 1       },
464         };
465         Assert.assertFalse(MatrixUtils.isSymmetric(MatrixUtils.createRealMatrix(dataNonSym1), eps));
466         final double[][] dataNonSym2 = {
467             { 1,   1, 0.99989 },
468             { 1,   1, 1       },
469             { 1.0, 1, 1       },
470         };
471         Assert.assertFalse(MatrixUtils.isSymmetric(MatrixUtils.createRealMatrix(dataNonSym2), eps));
472     }
473 
474     @Test
475     public void testCheckSymmetric1() {
476         final double[][] dataSym = {
477             { 1, 2, 3 },
478             { 2, 2, 5 },
479             { 3, 5, 6 },
480         };
481         MatrixUtils.checkSymmetric(MatrixUtils.createRealMatrix(dataSym), Math.ulp(1d));
482     }
483 
484     @Test(expected=MathIllegalArgumentException.class)
485     public void testCheckSymmetric2() {
486         final double[][] dataNonSym = {
487             { 1, 2, -3 },
488             { 2, 2, 5 },
489             { 3, 5, 6 },
490         };
491         MatrixUtils.checkSymmetric(MatrixUtils.createRealMatrix(dataNonSym), Math.ulp(1d));
492     }
493 
494     @Test(expected=MathIllegalArgumentException.class)
495     public void testInverseSingular() {
496         RealMatrix m = MatrixUtils.createRealMatrix(testData3x3Singular);
497         MatrixUtils.inverse(m);
498     }
499 
500     @Test(expected=MathIllegalArgumentException.class)
501     public void testInverseNonSquare() {
502         RealMatrix m = MatrixUtils.createRealMatrix(testData3x4);
503         MatrixUtils.inverse(m);
504     }
505 
506     @Test
507     public void testInverseDiagonalMatrix() {
508         final double[] data = { 1, 2, 3 };
509         final RealMatrix m = new DiagonalMatrix(data);
510         final RealMatrix inverse = MatrixUtils.inverse(m);
511 
512         final RealMatrix result = m.multiply(inverse);
513         UnitTestUtils.assertEquals("MatrixUtils.inverse() returns wrong result",
514                 MatrixUtils.createRealIdentityMatrix(data.length), result, Math.ulp(1d));
515     }
516 
517     @Test
518     public void testInverseRealMatrix() {
519         RealMatrix m = MatrixUtils.createRealMatrix(testData);
520         final RealMatrix inverse = MatrixUtils.inverse(m);
521 
522         final RealMatrix result = m.multiply(inverse);
523         UnitTestUtils.assertEquals("MatrixUtils.inverse() returns wrong result",
524                 MatrixUtils.createRealIdentityMatrix(testData.length), result, 1e-12);
525     }
526 
527     @Test
528     public void testMatrixExponentialNonSquare() {
529         double[][] exponentArr = {
530                 {0.0001, 0.001},
531                 {0.001, -0.0001},
532                 {0.001, -0.0001}
533         };
534         RealMatrix exponent = MatrixUtils.createRealMatrix(exponentArr);
535 
536         try {
537             MatrixUtils.matrixExponential(exponent);  // ragged
538             Assert.fail("Expecting MathIllegalArgumentException");
539         } catch (MathIllegalArgumentException ex) {
540             // expected
541         }
542     }
543 
544     @Test
545     public void testMatrixExponential3() {
546         double[][] exponentArr = {
547                 {0.0001, 0.001},
548                 {0.001, -0.0001}
549         };
550         RealMatrix exponent = MatrixUtils.createRealMatrix(exponentArr);
551 
552         double[][] expectedResultArr = {
553                 {1.00010050501688, 0.00100000016833332},
554                 {0.00100000016833332, 0.999900504983209}
555         };
556         RealMatrix expectedResult = MatrixUtils.createRealMatrix(expectedResultArr);
557 
558         UnitTestUtils.assertEquals("matrixExponential pade3 incorrect result",
559                 expectedResult, MatrixUtils.matrixExponential(exponent), 32.0 * Math.ulp(1.0));
560     }
561 
562 
563     @Test
564     public void testMatrixExponential5() {
565         double[][] exponentArr = {
566                 {0.1, 0.1},
567                 {0.001, -0.1}
568         };
569         RealMatrix exponent = MatrixUtils.createRealMatrix(exponentArr);
570 
571         double[][] expectedResultArr = {
572                 {1.10522267021001, 0.100168418362112},
573                 {0.00100168418362112, 0.904885833485786}
574         };
575         RealMatrix expectedResult = MatrixUtils.createRealMatrix(expectedResultArr);
576 
577         UnitTestUtils.assertEquals("matrixExponential pade5 incorrect result",
578                 expectedResult, MatrixUtils.matrixExponential(exponent), 2.0 * Math.ulp(1.0));
579     }
580 
581     @Test
582     public void testMatrixExponential7() {
583         double[][] exponentArr = {
584                 {0.5, 0.1},
585                 {0.001, -0.5}
586         };
587         RealMatrix exponent = MatrixUtils.createRealMatrix(exponentArr);
588 
589         double[][] expectedResultArr = {
590                 {1.64878192423569, 0.104220769814317},
591                 {0.00104220769814317, 0.606574226092523}
592         };
593         RealMatrix expectedResult = MatrixUtils.createRealMatrix(expectedResultArr);
594 
595         UnitTestUtils.assertEquals("matrixExponential pade7 incorrect result",
596                 expectedResult, MatrixUtils.matrixExponential(exponent), 32.0 * Math.ulp(1.0));
597     }
598 
599     @Test
600     public void testMatrixExponential9() {
601         double[][] exponentArr = {
602                 {1.8, 0.3},
603                 {0.001, -0.9}
604         };
605         RealMatrix exponent = MatrixUtils.createRealMatrix(exponentArr);
606 
607         double[][] expectedResultArr = {
608                 {6.05008743087114, 0.627036746099251},
609                 {0.00209012248699751, 0.406756715977872}
610         };
611         RealMatrix expectedResult = MatrixUtils.createRealMatrix(expectedResultArr);
612 
613         UnitTestUtils.assertEquals("matrixExponential pade9 incorrect result",
614                 expectedResult, MatrixUtils.matrixExponential(exponent), 16.0 * Math.ulp(1.0));
615     }
616 
617     @Test
618     public void testMatrixExponential13() {
619         double[][] exponentArr1 = {
620                 {3.4, 1.2},
621                 {0.001, -0.9}
622         };
623         RealMatrix exponent1 = MatrixUtils.createRealMatrix(exponentArr1);
624 
625         double[][] expectedResultArr1 = {
626                 {29.9705442872504, 8.2499077972773},
627                 {0.00687492316439775, 0.408374680340048}
628         };
629         RealMatrix expectedResult1 = MatrixUtils.createRealMatrix(expectedResultArr1);
630 
631         UnitTestUtils.assertEquals("matrixExponential pade13-1 incorrect result",
632                 expectedResult1, MatrixUtils.matrixExponential(exponent1), 16.0 * Math.ulp(30.0));
633 
634 
635         double[][] exponentArr2 = {
636                 {1.0, 1e5},
637                 {0.001, -1.0}
638         };
639         RealMatrix exponent2 = MatrixUtils.createRealMatrix(exponentArr2);
640 
641         double[][] expectedResultArr2 = {
642                 {12728.3536593144, 115190017.08756},
643                 {1.1519001708756, 10424.5533175632}
644         };
645         RealMatrix expectedResult2 = MatrixUtils.createRealMatrix(expectedResultArr2);
646 
647         UnitTestUtils.assertEquals("matrixExponential pade13-2 incorrect result",
648                 expectedResult2, MatrixUtils.matrixExponential(exponent2), 65536.0 * Math.ulp(1e8));
649 
650 
651         double[][] exponentArr3 = {
652                 {-1e4, 1e4},
653                 {1.0, -1.0}
654         };
655         RealMatrix exponent3 = MatrixUtils.createRealMatrix(exponentArr3);
656 
657         double[][] expectedResultArr3 = {
658                 {9.99900009999e-05, 0.999900009999},
659                 {9.99900009999e-05, 0.999900009999}
660         };
661         RealMatrix expectedResult3 = MatrixUtils.createRealMatrix(expectedResultArr3);
662 
663         UnitTestUtils.assertEquals("matrixExponential pade13-3 incorrect result",
664                 expectedResult3, MatrixUtils.matrixExponential(exponent3), 4096.0 * Math.ulp(1.0));
665     }
666 
667     @Test
668     public void testOrthonormalize1() {
669 
670         final List<RealVector> basis =
671                         MatrixUtils.orthonormalize(Arrays.asList(new ArrayRealVector(new double[] {  1, 2, 2 }),
672                                                                  new ArrayRealVector(new double[] { -1, 0, 2 }),
673                                                                  new ArrayRealVector(new double[] {  0, 0, 1 })),
674                                                    Precision.EPSILON, DependentVectorsHandler.GENERATE_EXCEPTION);
675         Assert.assertEquals(3, basis.size());
676         checkBasis(basis);
677         checkVector(basis.get(0),  1.0 / 3.0,  2.0 / 3.0, 2.0 / 3.0);
678         checkVector(basis.get(1), -2.0 / 3.0, -1.0 / 3.0, 2.0 / 3.0);
679         checkVector(basis.get(2),  2.0 / 3.0, -2.0 / 3.0, 1.0 / 3.0);
680 
681     }
682 
683     @Test
684     public void testOrthonormalize2() {
685 
686         final List<RealVector> basis =
687                         MatrixUtils.orthonormalize(Arrays.asList(new ArrayRealVector(new double[] { 3, 1 }),
688                                                                  new ArrayRealVector(new double[] { 2, 2 })),
689                                                    Precision.EPSILON, DependentVectorsHandler.GENERATE_EXCEPTION);
690         final double s10 = FastMath.sqrt(10);
691         Assert.assertEquals(2, basis.size());
692         checkBasis(basis);
693         checkVector(basis.get(0),  3 / s10,  1 / s10);
694         checkVector(basis.get(1), -1 / s10,  3 / s10);
695 
696     }
697 
698     @Test
699     public void testOrthonormalize3() {
700 
701         final double small = 1.0e-12;
702         final List<RealVector> basis =
703                         MatrixUtils.orthonormalize(Arrays.asList(new ArrayRealVector(new double[] { 1, small, small }),
704                                                                  new ArrayRealVector(new double[] { 1, small, 0     }),
705                                                                  new ArrayRealVector(new double[] { 1, 0,     small })),
706                                                    Precision.EPSILON, DependentVectorsHandler.GENERATE_EXCEPTION);
707         Assert.assertEquals(3, basis.size());
708         checkBasis(basis);
709         checkVector(basis.get(0), 1,  small, small);
710         checkVector(basis.get(1), 0,  0,     -1   );
711         checkVector(basis.get(2), 0, -1,      0   );
712 
713     }
714 
715     @Test
716     public void testOrthonormalizeIncompleteBasis() {
717 
718         final double small = 1.0e-12;
719         final List<RealVector> basis =
720                         MatrixUtils.orthonormalize(Arrays.asList(new ArrayRealVector(new double[] { 1, small, small }),
721                                                                  new ArrayRealVector(new double[] { 1, small, 0     })),
722                                                    Precision.EPSILON, DependentVectorsHandler.GENERATE_EXCEPTION);
723         Assert.assertEquals(2, basis.size());
724         checkBasis(basis);
725         checkVector(basis.get(0), 1,  small, small);
726         checkVector(basis.get(1), 0,  0,     -1   );
727 
728     }
729 
730     @Test
731     public void testOrthonormalizeDependent() {
732         final double small = 1.0e-12;
733         try {
734             MatrixUtils.orthonormalize(Arrays.asList(new ArrayRealVector(new double[] { 1, small, small }),
735                                                      new ArrayRealVector(new double[] { 1, small, small })),
736                                        Precision.EPSILON, DependentVectorsHandler.GENERATE_EXCEPTION);
737             Assert.fail("an exception should have been thrown");
738         } catch (MathIllegalArgumentException miae) {
739             Assert.assertEquals(LocalizedCoreFormats.ZERO_NORM, miae.getSpecifier());
740         }
741     }
742 
743     @Test
744     public void testOrthonormalizeDependentGenerateException() {
745         try {
746             MatrixUtils.orthonormalize(Arrays.asList(new ArrayRealVector(new double[] { 2, 3, 0 }),
747                                                      new ArrayRealVector(new double[] { 2, 7, 0 }),
748                                                      new ArrayRealVector(new double[] { 4, 5, 0 }),
749                                                      new ArrayRealVector(new double[] { 0, 0, 1 })),
750                                        7 * Precision.EPSILON, DependentVectorsHandler.GENERATE_EXCEPTION);
751             Assert.fail("an exception should have been thrown");
752         } catch (MathIllegalArgumentException miae) {
753             Assert.assertEquals(LocalizedCoreFormats.ZERO_NORM, miae.getSpecifier());
754         }
755 
756     }
757 
758     @Test
759     public void testOrthonormalizeDependentAddZeroNorm() {
760         List<RealVector> basis = MatrixUtils.orthonormalize(Arrays.asList(new ArrayRealVector(new double[] { 2, 3, 0 }),
761                                                                           new ArrayRealVector(new double[] { 2, 7, 0 }),
762                                                                           new ArrayRealVector(new double[] { 4, 5, 0 }),
763                                                                           new ArrayRealVector(new double[] { 0, 0, 1 })),
764                                                             7 * Precision.EPSILON, DependentVectorsHandler.ADD_ZERO_VECTOR);
765         Assert.assertEquals(4, basis.size());
766         Assert.assertEquals(0, basis.get(2).getEntry(0), 1.0e-15);
767         Assert.assertEquals(0, basis.get(2).getEntry(1), 1.0e-15);
768         Assert.assertEquals(0, basis.get(2).getEntry(2), 1.0e-15);
769     }
770 
771     @Test
772     public void testOrthonormalizeDependentReduceToSpan() {
773         List<RealVector> basis = MatrixUtils.orthonormalize(Arrays.asList(new ArrayRealVector(new double[] { 2, 3, 0 }),
774                                                                           new ArrayRealVector(new double[] { 2, 7, 0 }),
775                                                                           new ArrayRealVector(new double[] { 4, 5, 0 }),
776                                                                           new ArrayRealVector(new double[] { 0, 0, 1 })),
777                                                             7 * Precision.EPSILON, DependentVectorsHandler.REDUCE_BASE_TO_SPAN);
778         Assert.assertEquals(3, basis.size());
779         Assert.assertEquals(0, basis.get(2).getEntry(0), 1.0e-15);
780         Assert.assertEquals(0, basis.get(2).getEntry(1), 1.0e-15);
781         Assert.assertEquals(1, basis.get(2).getEntry(2), 1.0e-15);
782     }
783 
784     @Test
785     public void testFieldOrthonormalize1() {
786         doTestOrthonormalize1(Binary64Field.getInstance());
787     }
788 
789     @Test
790     public void testFieldOrthonormalize2() {
791         doTestOrthonormalize2(Binary64Field.getInstance());
792     }
793 
794     @Test
795     public void testFieldOrthonormalize3() {
796         doTestOrthonormalize3(Binary64Field.getInstance());
797     }
798 
799     @Test
800     public void testFieldOrthonormalizeIncompleteBasis() {
801         doTestOrthonormalizeIncompleteBasis(Binary64Field.getInstance());
802     }
803 
804     @Test
805     public void testFieldOrthonormalizeDependent() {
806         doTestOrthonormalizeDependent(Binary64Field.getInstance());
807     }
808 
809     @Test
810     public void testFieldOrthonormalizeDependentGenerateException() {
811         doTestOrthonormalizeDependentGenerateException(Binary64Field.getInstance());
812     }
813 
814     @Test
815     public void testFieldOrthonormalizeDependentAddZeroNorm() {
816         doTestOrthonormalizeDependentAddZeroNorm(Binary64Field.getInstance());
817     }
818 
819     @Test
820     public void testFieldOrthonormalizeDependentReduceToSpan() {
821         doTestOrthonormalizeDependentReduceToSpan(Binary64Field.getInstance());
822     }
823 
824     private <T extends CalculusFieldElement<T>> void doTestOrthonormalize1(final Field<T> field) {
825 
826         final List<FieldVector<T>> basis =
827                         MatrixUtils.orthonormalize(field,
828                                                    Arrays.asList(convert(field, new ArrayRealVector(new double[] {  1, 2, 2 })),
829                                                                  convert(field, new ArrayRealVector(new double[] { -1, 0, 2 })),
830                                                                  convert(field, new ArrayRealVector(new double[] {  0, 0, 1 }))),
831                                                    field.getZero().newInstance(Precision.EPSILON),
832                                                    DependentVectorsHandler.GENERATE_EXCEPTION);
833         Assert.assertEquals(3, basis.size());
834         checkBasis(field, basis);
835         checkVector(basis.get(0),  1.0 / 3.0,  2.0 / 3.0, 2.0 / 3.0);
836         checkVector(basis.get(1), -2.0 / 3.0, -1.0 / 3.0, 2.0 / 3.0);
837         checkVector(basis.get(2),  2.0 / 3.0, -2.0 / 3.0, 1.0 / 3.0);
838 
839     }
840 
841     private <T extends CalculusFieldElement<T>> void doTestOrthonormalize2(final Field<T> field) {
842 
843         final List<FieldVector<T>> basis =
844                         MatrixUtils.orthonormalize(field,
845                                                    Arrays.asList(convert(field, new ArrayRealVector(new double[] { 3, 1 })),
846                                                                  convert(field, new ArrayRealVector(new double[] { 2, 2 }))),
847                                                    field.getZero().newInstance(Precision.EPSILON),
848                                                    DependentVectorsHandler.GENERATE_EXCEPTION);
849         final double s10 = FastMath.sqrt(10);
850         Assert.assertEquals(2, basis.size());
851         checkBasis(field, basis);
852         checkVector(basis.get(0),  3 / s10,  1 / s10);
853         checkVector(basis.get(1), -1 / s10,  3 / s10);
854 
855     }
856 
857     private <T extends CalculusFieldElement<T>> void doTestOrthonormalize3(final Field<T> field) {
858 
859         final double small = 1.0e-12;
860         final List<FieldVector<T>> basis =
861                         MatrixUtils.orthonormalize(field,
862                                                    Arrays.asList(convert(field, new ArrayRealVector(new double[] { 1, small, small })),
863                                                                  convert(field, new ArrayRealVector(new double[] { 1, small, 0     })),
864                                                                  convert(field, new ArrayRealVector(new double[] { 1, 0,     small }))),
865                                                    field.getZero().newInstance(Precision.EPSILON),
866                                                    DependentVectorsHandler.GENERATE_EXCEPTION);
867         Assert.assertEquals(3, basis.size());
868         checkBasis(field, basis);
869         checkVector(basis.get(0), 1,  small, small);
870         checkVector(basis.get(1), 0,  0,     -1   );
871         checkVector(basis.get(2), 0, -1,      0   );
872 
873     }
874 
875     private <T extends CalculusFieldElement<T>> void doTestOrthonormalizeIncompleteBasis(final Field<T> field) {
876 
877         final double small = 1.0e-12;
878         final List<FieldVector<T>> basis =
879                         MatrixUtils.orthonormalize(field,
880                                                    Arrays.asList(convert(field, new ArrayRealVector(new double[] { 1, small, small })),
881                                                                  convert(field, new ArrayRealVector(new double[] { 1, small, 0     }))),
882                                                    field.getZero().newInstance(Precision.EPSILON),
883                                                    DependentVectorsHandler.GENERATE_EXCEPTION);
884         Assert.assertEquals(2, basis.size());
885         checkBasis(field, basis);
886         checkVector(basis.get(0), 1,  small, small);
887         checkVector(basis.get(1), 0,  0,     -1   );
888 
889     }
890 
891     private <T extends CalculusFieldElement<T>> void doTestOrthonormalizeDependent(final Field<T> field) {
892         final double small = 1.0e-12;
893         try {
894             MatrixUtils.orthonormalize(field,
895                                        Arrays.asList(convert(field, new ArrayRealVector(new double[] { 1, small, small })),
896                                                      convert(field, new ArrayRealVector(new double[] { 1, small, small }))),
897                                        field.getZero().newInstance(Precision.EPSILON),
898                                        DependentVectorsHandler.GENERATE_EXCEPTION);
899             Assert.fail("an exception should have been thrown");
900         } catch (MathIllegalArgumentException miae) {
901             Assert.assertEquals(LocalizedCoreFormats.ZERO_NORM, miae.getSpecifier());
902         }
903     }
904 
905     private <T extends CalculusFieldElement<T>> void doTestOrthonormalizeDependentGenerateException(final Field<T> field) {
906         try {
907             MatrixUtils.orthonormalize(field,
908                                        Arrays.asList(convert(field, new ArrayRealVector(new double[] { 2, 3, 0 })),
909                                                      convert(field, new ArrayRealVector(new double[] { 2, 7, 0 })),
910                                                      convert(field, new ArrayRealVector(new double[] { 4, 5, 0 })),
911                                                      convert(field, new ArrayRealVector(new double[] { 0, 0, 1 }))),
912                                        field.getZero().newInstance(7 * Precision.EPSILON),
913                                        DependentVectorsHandler.GENERATE_EXCEPTION);
914             Assert.fail("an exception should have been thrown");
915         } catch (MathIllegalArgumentException miae) {
916             Assert.assertEquals(LocalizedCoreFormats.ZERO_NORM, miae.getSpecifier());
917         }
918     }
919 
920     private <T extends CalculusFieldElement<T>> void doTestOrthonormalizeDependentAddZeroNorm(final Field<T> field) {
921         List<FieldVector<T>> basis = MatrixUtils.orthonormalize(field,
922                                                                 Arrays.asList(convert(field, new ArrayRealVector(new double[] { 2, 3, 0 })),
923                                                                               convert(field, new ArrayRealVector(new double[] { 2, 7, 0 })),
924                                                                               convert(field, new ArrayRealVector(new double[] { 4, 5, 0 })),
925                                                                               convert(field, new ArrayRealVector(new double[] { 0, 0, 1 }))),
926                                                                 field.getZero().newInstance(7 * Precision.EPSILON),
927                                                                 DependentVectorsHandler.ADD_ZERO_VECTOR);
928        Assert.assertEquals(4, basis.size());
929        Assert.assertEquals(0, basis.get(2).getEntry(0).getReal(), 1.0e-15);
930        Assert.assertEquals(0, basis.get(2).getEntry(1).getReal(), 1.0e-15);
931        Assert.assertEquals(0, basis.get(2).getEntry(2).getReal(), 1.0e-15);
932     }
933 
934     private <T extends CalculusFieldElement<T>> void doTestOrthonormalizeDependentReduceToSpan(final Field<T> field) {
935         List<FieldVector<T>> basis = MatrixUtils.orthonormalize(field,
936                                                                 Arrays.asList(convert(field, new ArrayRealVector(new double[] { 2, 3, 0 })),
937                                                                               convert(field, new ArrayRealVector(new double[] { 2, 7, 0 })),
938                                                                               convert(field, new ArrayRealVector(new double[] { 4, 5, 0 })),
939                                                                               convert(field, new ArrayRealVector(new double[] { 0, 0, 1 }))),
940                                                                 field.getZero().newInstance(7 * Precision.EPSILON),
941                                                                 DependentVectorsHandler.REDUCE_BASE_TO_SPAN);
942         Assert.assertEquals(3, basis.size());
943         Assert.assertEquals(0, basis.get(2).getEntry(0).getReal(), 1.0e-15);
944         Assert.assertEquals(0, basis.get(2).getEntry(1).getReal(), 1.0e-15);
945         Assert.assertEquals(1, basis.get(2).getEntry(2).getReal(), 1.0e-15);
946     }
947 
948     private void checkVector(final RealVector v, double... p) {
949         Assert.assertEquals(p.length, v.getDimension());
950         for (int i = 0; i < p.length; ++i) {
951             Assert.assertEquals(p[i], v.getEntry(i), 1.0e-15);
952         }
953     }
954 
955     private void checkBasis(final List<RealVector> basis) {
956         for (int i = 0; i < basis.size(); ++i) {
957             for (int j = i; j < basis.size(); ++j) {
958                 Assert.assertEquals(i == j ? 1.0 : 0.0, basis.get(i).dotProduct(basis.get(j)), 1.0e-12);
959             }
960         }        
961     }
962 
963     private <T extends CalculusFieldElement<T>> void checkVector(final FieldVector<T> v, double... p) {
964         Assert.assertEquals(p.length, v.getDimension());
965         for (int i = 0; i < p.length; ++i) {
966             Assert.assertEquals(p[i], v.getEntry(i).getReal(), 1.0e-15);
967         }
968     }
969 
970     private <T extends CalculusFieldElement<T>> void checkBasis(final Field<T> field, final List<FieldVector<T>> basis) {
971         for (int i = 0; i < basis.size(); ++i) {
972             for (int j = i; j < basis.size(); ++j) {
973                 Assert.assertEquals(i == j ? 1.0 : 0.0, basis.get(i).dotProduct(basis.get(j)).getReal(), 1.0e-12);
974             }
975         }        
976     }
977 
978     private <T extends CalculusFieldElement<T>> FieldVector<T> convert(final Field<T> field, final RealVector v) {
979         ArrayFieldVector<T> c = new ArrayFieldVector<T>(v.getDimension(), field.getZero());
980         for (int k = 0; k < v.getDimension(); ++k) {
981             c.setEntry(k, field.getZero().newInstance(v.getEntry(k)));
982         }
983         return c;
984     }
985 
986 }