1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16 package org.apache.commons.math.linear;
17
18 import junit.framework.Test;
19 import junit.framework.TestCase;
20 import junit.framework.TestSuite;
21
22 /***
23 * Test cases for the {@link RealMatrixImpl} class.
24 *
25 * @version $Revision: 1.20 $ $Date: 2004/10/12 06:19:50 $
26 */
27
28 public final class RealMatrixImplTest extends TestCase {
29
30
31 protected double[][] id = { {1d,0d,0d}, {0d,1d,0d}, {0d,0d,1d} };
32
33
34 protected double[][] testData = { {1d,2d,3d}, {2d,5d,3d}, {1d,0d,8d} };
35 protected double[][] testDataLU = {{2d, 5d, 3d}, {.5d, -2.5d, 6.5d}, {0.5d, 0.2d, .2d}};
36 protected double[][] testDataPlus2 = { {3d,4d,5d}, {4d,7d,5d}, {3d,2d,10d} };
37 protected double[][] testDataMinus = { {-1d,-2d,-3d}, {-2d,-5d,-3d},
38 {-1d,0d,-8d} };
39 protected double[] testDataRow1 = {1d,2d,3d};
40 protected double[] testDataCol3 = {3d,3d,8d};
41 protected double[][] testDataInv =
42 { {-40d,16d,9d}, {13d,-5d,-3d}, {5d,-2d,-1d} };
43 protected double[] preMultTest = {8,12,33};
44 protected double[][] testData2 ={ {1d,2d,3d}, {2d,5d,3d}};
45 protected double[][] testData2T = { {1d,2d}, {2d,5d}, {3d,3d}};
46 protected double[][] testDataPlusInv =
47 { {-39d,18d,12d}, {15d,0d,0d}, {6d,-2d,7d} };
48
49
50 protected double[][] luData = { {2d,3d,3d}, {0d,5d,7d}, {6d,9d,8d} };
51 protected double[][] luDataLUDecomposition = { {6d,9d,8d}, {0d,5d,7d},
52 {0.33333333333333,0d,0.33333333333333} };
53
54
55 protected double[][] singular = { {2d,3d}, {2d,3d} };
56 protected double[][] bigSingular = {{1d,2d,3d,4d}, {2d,5d,3d,4d},
57 {7d,3d,256d,1930d}, {3d,7d,6d,8d}};
58 protected double[][] detData = { {1d,2d,3d}, {4d,5d,6d}, {7d,8d,10d} };
59 protected double[][] detData2 = { {1d, 3d}, {2d, 4d}};
60
61
62 protected double[] testVector = {1,2,3};
63 protected double[] testVector2 = {1,2,3,4};
64
65
66 protected double[][] subTestData = {{1, 2, 3, 4}, {1.5, 2.5, 3.5, 4.5},
67 {2, 4, 6, 8}, {4, 5, 6, 7}};
68
69 protected double[][] subRows02Cols13 = { {2, 4}, {4, 8}};
70 protected double[][] subRows03Cols12 = { {2, 3}, {5, 6}};
71 protected double[][] subRows03Cols123 = { {2, 3, 4} , {5, 6, 7}};
72
73 protected double[][] subRows20Cols123 = { {4, 6, 8} , {2, 3, 4}};
74 protected double[][] subRows31Cols31 = {{7, 5}, {4.5, 2.5}};
75
76 protected double[][] subRows01Cols23 = {{3,4} , {3.5, 4.5}};
77 protected double[][] subRows23Cols00 = {{2} , {4}};
78 protected double[][] subRows00Cols33 = {{4}};
79
80 protected double[][] subRow0 = {{1,2,3,4}};
81 protected double[][] subRow3 = {{4,5,6,7}};
82
83 protected double[][] subColumn1 = {{2}, {2.5}, {4}, {5}};
84 protected double[][] subColumn3 = {{4}, {4.5}, {8}, {7}};
85
86
87 protected double entryTolerance = 10E-16;
88 protected double normTolerance = 10E-14;
89
90 public RealMatrixImplTest(String name) {
91 super(name);
92 }
93
94 public void setUp() {
95
96 }
97
98 public static Test suite() {
99 TestSuite suite = new TestSuite(RealMatrixImplTest.class);
100 suite.setName("RealMatrixImpl Tests");
101 return suite;
102 }
103
104 /*** test dimensions */
105 public void testDimensions() {
106 RealMatrixImpl m = new RealMatrixImpl(testData);
107 RealMatrixImpl m2 = new RealMatrixImpl(testData2);
108 assertEquals("testData row dimension",3,m.getRowDimension());
109 assertEquals("testData column dimension",3,m.getColumnDimension());
110 assertTrue("testData is square",m.isSquare());
111 assertEquals("testData2 row dimension",m2.getRowDimension(),2);
112 assertEquals("testData2 column dimension",m2.getColumnDimension(),3);
113 assertTrue("testData2 is not square",!m2.isSquare());
114 }
115
116 /*** test copy functions */
117 public void testCopyFunctions() {
118 RealMatrixImpl m = new RealMatrixImpl(testData);
119 RealMatrixImpl m2 = new RealMatrixImpl(m.getData());
120 assertEquals(m2,m);
121 }
122
123 /*** test add */
124 public void testAdd() {
125 RealMatrixImpl m = new RealMatrixImpl(testData);
126 RealMatrixImpl mInv = new RealMatrixImpl(testDataInv);
127 RealMatrixImpl mPlusMInv = (RealMatrixImpl)m.add(mInv);
128 double[][] sumEntries = mPlusMInv.getData();
129 for (int row = 0; row < m.getRowDimension(); row++) {
130 for (int col = 0; col < m.getColumnDimension(); col++) {
131 assertEquals("sum entry entry",
132 testDataPlusInv[row][col],sumEntries[row][col],
133 entryTolerance);
134 }
135 }
136 }
137
138 /*** test add failure */
139 public void testAddFail() {
140 RealMatrixImpl m = new RealMatrixImpl(testData);
141 RealMatrixImpl m2 = new RealMatrixImpl(testData2);
142 try {
143 RealMatrixImpl mPlusMInv = (RealMatrixImpl)m.add(m2);
144 fail("IllegalArgumentException expected");
145 } catch (IllegalArgumentException ex) {
146 ;
147 }
148 }
149
150 /*** test norm */
151 public void testNorm() {
152 RealMatrixImpl m = new RealMatrixImpl(testData);
153 RealMatrixImpl m2 = new RealMatrixImpl(testData2);
154 assertEquals("testData norm",14d,m.getNorm(),entryTolerance);
155 assertEquals("testData2 norm",7d,m2.getNorm(),entryTolerance);
156 }
157
158 /*** test m-n = m + -n */
159 public void testPlusMinus() {
160 RealMatrixImpl m = new RealMatrixImpl(testData);
161 RealMatrixImpl m2 = new RealMatrixImpl(testDataInv);
162 assertClose("m-n = m + -n",m.subtract(m2),
163 m2.scalarMultiply(-1d).add(m),entryTolerance);
164 try {
165 RealMatrix a = m.subtract(new RealMatrixImpl(testData2));
166 fail("Expecting illegalArgumentException");
167 } catch (IllegalArgumentException ex) {
168 ;
169 }
170 }
171
172 /*** test multiply */
173 public void testMultiply() {
174 RealMatrixImpl m = new RealMatrixImpl(testData);
175 RealMatrixImpl mInv = new RealMatrixImpl(testDataInv);
176 RealMatrixImpl identity = new RealMatrixImpl(id);
177 RealMatrixImpl m2 = new RealMatrixImpl(testData2);
178 assertClose("inverse multiply",m.multiply(mInv),
179 identity,entryTolerance);
180 assertClose("inverse multiply",mInv.multiply(m),
181 identity,entryTolerance);
182 assertClose("identity multiply",m.multiply(identity),
183 m,entryTolerance);
184 assertClose("identity multiply",identity.multiply(mInv),
185 mInv,entryTolerance);
186 assertClose("identity multiply",m2.multiply(identity),
187 m2,entryTolerance);
188 try {
189 RealMatrix a = m.multiply(new RealMatrixImpl(bigSingular));
190 fail("Expecting illegalArgumentException");
191 } catch (IllegalArgumentException ex) {
192 ;
193 }
194 }
195
196
197
198 private double[][] d3 = new double[][] {{1,2,3,4},{5,6,7,8}};
199 private double[][] d4 = new double[][] {{1},{2},{3},{4}};
200 private double[][] d5 = new double[][] {{30},{70}};
201
202 public void testMultiply2() {
203 RealMatrix m3 = new RealMatrixImpl(d3);
204 RealMatrix m4 = new RealMatrixImpl(d4);
205 RealMatrix m5 = new RealMatrixImpl(d5);
206 assertClose("m3*m4=m5", m3.multiply(m4), m5, entryTolerance);
207 }
208
209 /*** test isSingular */
210 public void testIsSingular() {
211 RealMatrixImpl m = new RealMatrixImpl(singular);
212 assertTrue("singular",m.isSingular());
213 m = new RealMatrixImpl(bigSingular);
214 assertTrue("big singular",m.isSingular());
215 m = new RealMatrixImpl(id);
216 assertTrue("identity nonsingular",!m.isSingular());
217 m = new RealMatrixImpl(testData);
218 assertTrue("testData nonsingular",!m.isSingular());
219 }
220
221 /*** test inverse */
222 public void testInverse() {
223 RealMatrixImpl m = new RealMatrixImpl(testData);
224 RealMatrix mInv = new RealMatrixImpl(testDataInv);
225 assertClose("inverse",mInv,m.inverse(),normTolerance);
226 assertClose("inverse^2",m,m.inverse().inverse(),10E-12);
227
228
229 m = new RealMatrixImpl(testData2);
230 try {
231 m.inverse();
232 fail("Expecting InvalidMatrixException");
233 } catch (InvalidMatrixException ex) {
234
235 }
236
237
238 m = new RealMatrixImpl(singular);
239 try {
240 m.inverse();
241 fail("Expecting InvalidMatrixException");
242 } catch (InvalidMatrixException ex) {
243
244 }
245 }
246
247 /*** test solve */
248 public void testSolve() {
249 RealMatrixImpl m = new RealMatrixImpl(testData);
250 RealMatrix mInv = new RealMatrixImpl(testDataInv);
251
252 assertClose("inverse-operate",mInv.operate(testVector),
253 m.solve(testVector),normTolerance);
254 try {
255 double[] x = m.solve(testVector2);
256 fail("expecting IllegalArgumentException");
257 } catch (IllegalArgumentException ex) {
258 ;
259 }
260 RealMatrix bs = new RealMatrixImpl(bigSingular);
261 try {
262 RealMatrix a = bs.solve(bs);
263 fail("Expecting InvalidMatrixException");
264 } catch (InvalidMatrixException ex) {
265 ;
266 }
267 try {
268 RealMatrix a = m.solve(bs);
269 fail("Expecting IllegalArgumentException");
270 } catch (IllegalArgumentException ex) {
271 ;
272 }
273 try {
274 RealMatrix a = (new RealMatrixImpl(testData2)).solve(bs);
275 fail("Expecting illegalArgumentException");
276 } catch (IllegalArgumentException ex) {
277 ;
278 }
279 try {
280 (new RealMatrixImpl(testData2)).luDecompose();
281 fail("Expecting InvalidMatrixException");
282 } catch (InvalidMatrixException ex) {
283 ;
284 }
285 }
286
287 /*** test determinant */
288 public void testDeterminant() {
289 RealMatrix m = new RealMatrixImpl(bigSingular);
290 assertEquals("singular determinant",0,m.getDeterminant(),0);
291 m = new RealMatrixImpl(detData);
292 assertEquals("nonsingular test",-3d,m.getDeterminant(),normTolerance);
293
294
295 m = new RealMatrixImpl(detData2);
296 assertEquals("nonsingular R test 1",-2d,m.getDeterminant(),normTolerance);
297 m = new RealMatrixImpl(testData);
298 assertEquals("nonsingular R test 2",-1d,m.getDeterminant(),normTolerance);
299
300 try {
301 double a = new RealMatrixImpl(testData2).getDeterminant();
302 fail("Expecting InvalidMatrixException");
303 } catch (InvalidMatrixException ex) {
304 ;
305 }
306 }
307
308 /*** test trace */
309 public void testTrace() {
310 RealMatrix m = new RealMatrixImpl(id);
311 assertEquals("identity trace",3d,m.getTrace(),entryTolerance);
312 m = new RealMatrixImpl(testData2);
313 try {
314 double x = m.getTrace();
315 fail("Expecting illegalArgumentException");
316 } catch (IllegalArgumentException ex) {
317 ;
318 }
319 }
320
321 /*** test sclarAdd */
322 public void testScalarAdd() {
323 RealMatrix m = new RealMatrixImpl(testData);
324 assertClose("scalar add",new RealMatrixImpl(testDataPlus2),
325 m.scalarAdd(2d),entryTolerance);
326 }
327
328 /*** test operate */
329 public void testOperate() {
330 RealMatrix m = new RealMatrixImpl(id);
331 double[] x = m.operate(testVector);
332 assertClose("identity operate",testVector,x,entryTolerance);
333 m = new RealMatrixImpl(bigSingular);
334 try {
335 x = m.operate(testVector);
336 fail("Expecting illegalArgumentException");
337 } catch (IllegalArgumentException ex) {
338 ;
339 }
340 }
341
342 /*** test transpose */
343 public void testTranspose() {
344 RealMatrix m = new RealMatrixImpl(testData);
345 assertClose("inverse-transpose",m.inverse().transpose(),
346 m.transpose().inverse(),normTolerance);
347 m = new RealMatrixImpl(testData2);
348 RealMatrix mt = new RealMatrixImpl(testData2T);
349 assertClose("transpose",mt,m.transpose(),normTolerance);
350 }
351
352 /*** test preMultiply by vector */
353 public void testPremultiplyVector() {
354 RealMatrix m = new RealMatrixImpl(testData);
355 assertClose("premultiply",m.preMultiply(testVector),preMultTest,normTolerance);
356 m = new RealMatrixImpl(bigSingular);
357 try {
358 m.preMultiply(testVector);
359 fail("expecting IllegalArgumentException");
360 } catch (IllegalArgumentException ex) {
361 ;
362 }
363 }
364
365 public void testPremultiply() {
366 RealMatrix m3 = new RealMatrixImpl(d3);
367 RealMatrix m4 = new RealMatrixImpl(d4);
368 RealMatrix m5 = new RealMatrixImpl(d5);
369 assertClose("m3*m4=m5", m4.preMultiply(m3), m5, entryTolerance);
370
371 RealMatrixImpl m = new RealMatrixImpl(testData);
372 RealMatrixImpl mInv = new RealMatrixImpl(testDataInv);
373 RealMatrixImpl identity = new RealMatrixImpl(id);
374 RealMatrixImpl m2 = new RealMatrixImpl(testData2);
375 assertClose("inverse multiply",m.preMultiply(mInv),
376 identity,entryTolerance);
377 assertClose("inverse multiply",mInv.preMultiply(m),
378 identity,entryTolerance);
379 assertClose("identity multiply",m.preMultiply(identity),
380 m,entryTolerance);
381 assertClose("identity multiply",identity.preMultiply(mInv),
382 mInv,entryTolerance);
383 try {
384 RealMatrix a = m.preMultiply(new RealMatrixImpl(bigSingular));
385 fail("Expecting illegalArgumentException");
386 } catch (IllegalArgumentException ex) {
387 ;
388 }
389 }
390
391 public void testGetVectors() {
392 RealMatrix m = new RealMatrixImpl(testData);
393 assertClose("get row",m.getRow(0),testDataRow1,entryTolerance);
394 assertClose("get col",m.getColumn(2),testDataCol3,entryTolerance);
395 try {
396 double[] x = m.getRow(10);
397 fail("expecting MatrixIndexException");
398 } catch (MatrixIndexException ex) {
399 ;
400 }
401 try {
402 double[] x = m.getColumn(-1);
403 fail("expecting MatrixIndexException");
404 } catch (MatrixIndexException ex) {
405 ;
406 }
407 }
408
409 public void testGetEntry() {
410 RealMatrix m = new RealMatrixImpl(testData);
411 assertEquals("get entry",m.getEntry(0,1),2d,entryTolerance);
412 try {
413 m.getEntry(10, 4);
414 fail ("Expecting MatrixIndexException");
415 } catch (MatrixIndexException ex) {
416
417 }
418 }
419
420 public void testLUDecomposition() throws Exception {
421 RealMatrixImpl m = new RealMatrixImpl(testData);
422 RealMatrix lu = m.getLUMatrix();
423 assertClose("LU decomposition", lu, (RealMatrix) new RealMatrixImpl(testDataLU), normTolerance);
424 verifyDecomposition(m, lu);
425 m = new RealMatrixImpl(luData);
426 lu = m.getLUMatrix();
427 assertClose("LU decomposition", lu, (RealMatrix) new RealMatrixImpl(luDataLUDecomposition), normTolerance);
428 verifyDecomposition(m, lu);
429 m = new RealMatrixImpl(testDataMinus);
430 lu = m.getLUMatrix();
431 verifyDecomposition(m, lu);
432 m = new RealMatrixImpl(id);
433 lu = m.getLUMatrix();
434 verifyDecomposition(m, lu);
435 try {
436 m = new RealMatrixImpl(bigSingular);
437 lu = m.getLUMatrix();
438 fail("Expecting InvalidMatrixException");
439 } catch (InvalidMatrixException ex) {
440
441 }
442 try {
443 m = new RealMatrixImpl(testData2);
444 lu = m.getLUMatrix();
445 fail("Expecting InvalidMatrixException");
446 } catch (InvalidMatrixException ex) {
447
448 }
449 }
450
451 /*** test examples in user guide */
452 public void testExamples() {
453
454 double[][] matrixData = { {1d,2d,3d}, {2d,5d,3d}};
455 RealMatrix m = new RealMatrixImpl(matrixData);
456
457 double[][] matrixData2 = { {1d,2d}, {2d,5d}, {1d, 7d}};
458 RealMatrix n = new RealMatrixImpl(matrixData2);
459
460 RealMatrix p = m.multiply(n);
461 assertEquals(2, p.getRowDimension());
462 assertEquals(2, p.getColumnDimension());
463
464 RealMatrix pInverse = p.inverse();
465 assertEquals(2, pInverse.getRowDimension());
466 assertEquals(2, pInverse.getColumnDimension());
467
468
469 double[][] coefficientsData = {{2, 3, -2}, {-1, 7, 6}, {4, -3, -5}};
470 RealMatrix coefficients = new RealMatrixImpl(coefficientsData);
471 double[] constants = {1, -2, 1};
472 double[] solution = coefficients.solve(constants);
473 assertEquals(2 * solution[0] + 3 * solution[1] -2 * solution[2], constants[0], 1E-12);
474 assertEquals(-1 * solution[0] + 7 * solution[1] + 6 * solution[2], constants[1], 1E-12);
475 assertEquals(4 * solution[0] - 3 * solution[1] -5 * solution[2], constants[2], 1E-12);
476
477 }
478
479
480 public void testSubMatrix() {
481 RealMatrix m = new RealMatrixImpl(subTestData);
482 RealMatrix mRows23Cols00 = new RealMatrixImpl(subRows23Cols00);
483 RealMatrix mRows00Cols33 = new RealMatrixImpl(subRows00Cols33);
484 RealMatrix mRows01Cols23 = new RealMatrixImpl(subRows01Cols23);
485 RealMatrix mRows02Cols13 = new RealMatrixImpl(subRows02Cols13);
486 RealMatrix mRows03Cols12 = new RealMatrixImpl(subRows03Cols12);
487 RealMatrix mRows03Cols123 = new RealMatrixImpl(subRows03Cols123);
488 RealMatrix mRows20Cols123 = new RealMatrixImpl(subRows20Cols123);
489 RealMatrix mRows31Cols31 = new RealMatrixImpl(subRows31Cols31);
490 assertEquals("Rows23Cols00", mRows23Cols00,
491 m.getSubMatrix(2 , 3 , 0, 0));
492 assertEquals("Rows00Cols33", mRows00Cols33,
493 m.getSubMatrix(0 , 0 , 3, 3));
494 assertEquals("Rows01Cols23", mRows01Cols23,
495 m.getSubMatrix(0 , 1 , 2, 3));
496 assertEquals("Rows02Cols13", mRows02Cols13,
497 m.getSubMatrix(new int[] {0,2}, new int[] {1,3}));
498 assertEquals("Rows03Cols12", mRows03Cols12,
499 m.getSubMatrix(new int[] {0,3}, new int[] {1,2}));
500 assertEquals("Rows03Cols123", mRows03Cols123,
501 m.getSubMatrix(new int[] {0,3}, new int[] {1,2,3}));
502 assertEquals("Rows20Cols123", mRows20Cols123,
503 m.getSubMatrix(new int[] {2,0}, new int[] {1,2,3}));
504 assertEquals("Rows31Cols31", mRows31Cols31,
505 m.getSubMatrix(new int[] {3,1}, new int[] {3,1}));
506 assertEquals("Rows31Cols31", mRows31Cols31,
507 m.getSubMatrix(new int[] {3,1}, new int[] {3,1}));
508
509 try {
510 m.getSubMatrix(1,0,2,4);
511 fail("Expecting MatrixIndexException");
512 } catch (MatrixIndexException ex) {
513
514 }
515 try {
516 m.getSubMatrix(-1,1,2,2);
517 fail("Expecting MatrixIndexException");
518 } catch (MatrixIndexException ex) {
519
520 }
521 try {
522 m.getSubMatrix(1,0,2,2);
523 fail("Expecting MatrixIndexException");
524 } catch (MatrixIndexException ex) {
525
526 }
527 try {
528 m.getSubMatrix(1,0,2,4);
529 fail("Expecting MatrixIndexException");
530 } catch (MatrixIndexException ex) {
531
532 }
533 try {
534 m.getSubMatrix(new int[] {}, new int[] {0});
535 fail("Expecting MatrixIndexException");
536 } catch (MatrixIndexException ex) {
537
538 }
539 try {
540 m.getSubMatrix(new int[] {0}, new int[] {4});
541 fail("Expecting MatrixIndexException");
542 } catch (MatrixIndexException ex) {
543
544 }
545 }
546
547 public void testGetRowMatrix() {
548 RealMatrix m = new RealMatrixImpl(subTestData);
549 RealMatrix mRow0 = new RealMatrixImpl(subRow0);
550 RealMatrix mRow3 = new RealMatrixImpl(subRow3);
551 assertEquals("Row0", mRow0,
552 m.getRowMatrix(0));
553 assertEquals("Row3", mRow3,
554 m.getRowMatrix(3));
555 try {
556 m.getRowMatrix(-1);
557 fail("Expecting MatrixIndexException");
558 } catch (MatrixIndexException ex) {
559
560 }
561 try {
562 m.getRowMatrix(4);
563 fail("Expecting MatrixIndexException");
564 } catch (MatrixIndexException ex) {
565
566 }
567 }
568
569 public void testGetColumnMatrix() {
570 RealMatrix m = new RealMatrixImpl(subTestData);
571 RealMatrix mColumn1 = new RealMatrixImpl(subColumn1);
572 RealMatrix mColumn3 = new RealMatrixImpl(subColumn3);
573 assertEquals("Column1", mColumn1,
574 m.getColumnMatrix(1));
575 assertEquals("Column3", mColumn3,
576 m.getColumnMatrix(3));
577 try {
578 m.getColumnMatrix(-1);
579 fail("Expecting MatrixIndexException");
580 } catch (MatrixIndexException ex) {
581
582 }
583 try {
584 m.getColumnMatrix(4);
585 fail("Expecting MatrixIndexException");
586 } catch (MatrixIndexException ex) {
587
588 }
589 }
590
591 public void testEqualsAndHashCode() {
592 RealMatrixImpl m = new RealMatrixImpl(testData);
593 RealMatrixImpl m1 = (RealMatrixImpl) m.copy();
594 RealMatrixImpl mt = (RealMatrixImpl) m.transpose();
595 assertTrue(m.hashCode() != mt.hashCode());
596 assertEquals(m.hashCode(), m1.hashCode());
597 assertEquals(m, m);
598 assertEquals(m, m1);
599 assertFalse(m.equals(null));
600 assertFalse(m.equals(mt));
601 assertFalse(m.equals(new RealMatrixImpl(bigSingular)));
602 }
603
604 public void testToString() {
605 RealMatrixImpl m = new RealMatrixImpl(testData);
606 assertEquals("RealMatrixImpl{{1.0,2.0,3.0},{2.0,5.0,3.0},{1.0,0.0,8.0}}",
607 m.toString());
608 m = new RealMatrixImpl();
609 assertEquals("RealMatrixImpl{}",
610 m.toString());
611 }
612
613
614
615 /*** verifies that two matrices are close (1-norm) */
616 protected void assertClose(String msg, RealMatrix m, RealMatrix n,
617 double tolerance) {
618 assertTrue(msg,m.subtract(n).getNorm() < tolerance);
619 }
620
621 /*** verifies that two vectors are close (sup norm) */
622 protected void assertClose(String msg, double[] m, double[] n,
623 double tolerance) {
624 if (m.length != n.length) {
625 fail("vectors not same length");
626 }
627 for (int i = 0; i < m.length; i++) {
628 assertEquals(msg + " " + i + " elements differ",
629 m[i],n[i],tolerance);
630 }
631 }
632
633 /*** extracts the l and u matrices from compact lu representation */
634 protected void splitLU(RealMatrix lu, double[][] lowerData, double[][] upperData) throws InvalidMatrixException {
635 if (!lu.isSquare() || lowerData.length != lowerData[0].length || upperData.length != upperData[0].length ||
636 lowerData.length != upperData.length
637 || lowerData.length != lu.getRowDimension()) {
638 throw new InvalidMatrixException("incorrect dimensions");
639 }
640 int n = lu.getRowDimension();
641 for (int i = 0; i < n; i++) {
642 for (int j = 0; j < n; j++) {
643 if (j < i) {
644 lowerData[i][j] = lu.getEntry(i, j);
645 upperData[i][j] = 0d;
646 } else if (i == j) {
647 lowerData[i][j] = 1d;
648 upperData[i][j] = lu.getEntry(i, j);
649 } else {
650 lowerData[i][j] = 0d;
651 upperData[i][j] = lu.getEntry(i, j);
652 }
653 }
654 }
655 }
656
657 /*** Returns the result of applying the given row permutation to the matrix */
658 protected RealMatrix permuteRows(RealMatrix matrix, int[] permutation) {
659 if (!matrix.isSquare() || matrix.getRowDimension() != permutation.length) {
660 throw new IllegalArgumentException("dimension mismatch");
661 }
662 int n = matrix.getRowDimension();
663 int m = matrix.getColumnDimension();
664 double out[][] = new double[m][n];
665 for (int i = 0; i < n; i++) {
666 for (int j = 0; j < m; j++) {
667 out[i][j] = matrix.getEntry(permutation[i], j);
668 }
669 }
670 return new RealMatrixImpl(out);
671 }
672
673 /*** Extracts l and u matrices from lu and verifies that matrix = l times u modulo permutation */
674 protected void verifyDecomposition(RealMatrix matrix, RealMatrix lu) throws Exception{
675 int n = matrix.getRowDimension();
676 double[][] lowerData = new double[n][n];
677 double[][] upperData = new double[n][n];
678 splitLU(lu, lowerData, upperData);
679 RealMatrix lower =new RealMatrixImpl(lowerData);
680 RealMatrix upper = new RealMatrixImpl(upperData);
681 int[] permutation = ((RealMatrixImpl) matrix).getPermutation();
682 RealMatrix permuted = permuteRows(matrix, permutation);
683 assertClose("lu decomposition does not work", permuted, lower.multiply(upper), normTolerance);
684 }
685
686
687 /*** Useful for debugging */
688 private void dumpMatrix(RealMatrix m) {
689 for (int i = 0; i < m.getRowDimension(); i++) {
690 String os = "";
691 for (int j = 0; j < m.getColumnDimension(); j++) {
692 os += m.getEntry(i, j) + " ";
693 }
694 System.out.println(os);
695 }
696 }
697
698 }
699