1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
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
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}});
88 Assert.fail("Expecting MathIllegalArgumentException");
89 } catch (MathIllegalArgumentException ex) {
90
91 }
92 try {
93 MatrixUtils.createRealMatrix(new double[][] {{}, {}});
94 Assert.fail("Expecting MathIllegalArgumentException");
95 } catch (MathIllegalArgumentException ex) {
96
97 }
98 try {
99 MatrixUtils.createRealMatrix(null);
100 Assert.fail("Expecting NullArgumentException");
101 } catch (NullArgumentException ex) {
102
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}}));
114 Assert.fail("Expecting MathIllegalArgumentException");
115 } catch (MathIllegalArgumentException ex) {
116
117 }
118 try {
119 MatrixUtils.createFieldMatrix(asFraction(new double[][] {{}, {}}));
120 Assert.fail("Expecting MathIllegalArgumentException");
121 } catch (MathIllegalArgumentException ex) {
122
123 }
124 try {
125 MatrixUtils.createFieldMatrix((Fraction[][])null);
126 Assert.fail("Expecting NullArgumentException");
127 } catch (NullArgumentException ex) {
128
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[] {});
168 Assert.fail("Expecting MathIllegalArgumentException");
169 } catch (MathIllegalArgumentException ex) {
170
171 }
172 try {
173 MatrixUtils.createRowRealMatrix(null);
174 Assert.fail("Expecting NullArgumentException");
175 } catch (NullArgumentException ex) {
176
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[] {});
188 Assert.fail("Expecting MathIllegalArgumentException");
189 } catch (MathIllegalArgumentException ex) {
190
191 }
192 try {
193 MatrixUtils.createRowFieldMatrix((Fraction[]) null);
194 Assert.fail("Expecting NullArgumentException");
195 } catch (NullArgumentException ex) {
196
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[] {});
206 Assert.fail("Expecting MathIllegalArgumentException");
207 } catch (MathIllegalArgumentException ex) {
208
209 }
210 try {
211 MatrixUtils.createColumnRealMatrix(null);
212 Assert.fail("Expecting NullArgumentException");
213 } catch (NullArgumentException ex) {
214
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[] {});
227 Assert.fail("Expecting MathIllegalArgumentException");
228 } catch (MathIllegalArgumentException ex) {
229
230 }
231 try {
232 MatrixUtils.createColumnFieldMatrix((Fraction[]) null);
233 Assert.fail("Expecting NullArgumentException");
234 } catch (NullArgumentException ex) {
235
236 }
237 }
238
239
240
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
264 }
265 }
266
267
268
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
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 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 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
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
376
377
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
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);
538 Assert.fail("Expecting MathIllegalArgumentException");
539 } catch (MathIllegalArgumentException ex) {
540
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 }