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