1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.apache.commons.math4.legacy.linear;
18
19 import java.util.Arrays;
20
21 import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
22 import org.apache.commons.math4.core.jdkmath.JdkMath;
23 import org.junit.Assert;
24 import org.junit.Test;
25
26 public class SymmLQTest {
27
28 public void saundersTest(final int n, final boolean goodb,
29 final boolean precon, final double shift,
30 final double pertbn) {
31 final RealLinearOperator a = new RealLinearOperator() {
32
33 @Override
34 public RealVector operate(final RealVector x) {
35 if (x.getDimension() != n) {
36 throw new DimensionMismatchException(x.getDimension(), n);
37 }
38 final double[] y = new double[n];
39 for (int i = 0; i < n; i++) {
40 y[i] = (i + 1) * 1.1 / n * x.getEntry(i);
41 }
42 return new ArrayRealVector(y, false);
43 }
44
45 @Override
46 public int getRowDimension() {
47 return n;
48 }
49
50 @Override
51 public int getColumnDimension() {
52 return n;
53 }
54 };
55 final double shiftm = shift;
56 final double pertm = JdkMath.abs(pertbn);
57 final RealLinearOperator minv;
58 if (precon) {
59 minv = new RealLinearOperator() {
60 @Override
61 public int getRowDimension() {
62 return n;
63 }
64
65 @Override
66 public int getColumnDimension() {
67 return n;
68 }
69
70 @Override
71 public RealVector operate(final RealVector x) {
72 if (x.getDimension() != n) {
73 throw new DimensionMismatchException(x.getDimension(),
74 n);
75 }
76 final double[] y = new double[n];
77 for (int i = 0; i < n; i++) {
78 double d = (i + 1) * 1.1 / n;
79 d = JdkMath.abs(d - shiftm);
80 if (i % 10 == 0) {
81 d += pertm;
82 }
83 y[i] = x.getEntry(i) / d;
84 }
85 return new ArrayRealVector(y, false);
86 }
87 };
88 } else {
89 minv = null;
90 }
91 final RealVector xtrue = new ArrayRealVector(n);
92 for (int i = 0; i < n; i++) {
93 xtrue.setEntry(i, n - i);
94 }
95 final RealVector b = a.operate(xtrue);
96 b.combineToSelf(1.0, -shift, xtrue);
97 final SymmLQ solver = new SymmLQ(2 * n, 1E-12, true);
98 final RealVector x = solver.solve(a, minv, b, goodb, shift);
99 final RealVector y = a.operate(x);
100 final RealVector r1 = new ArrayRealVector(n);
101 for (int i = 0; i < n; i++) {
102 final double bi = b.getEntry(i);
103 final double yi = y.getEntry(i);
104 final double xi = x.getEntry(i);
105 r1.setEntry(i, bi - yi + shift * xi);
106 }
107 final double enorm = x.subtract(xtrue).getNorm() / xtrue.getNorm();
108 final double etol = 1E-5;
109 Assert.assertTrue("enorm=" + enorm + ", " +
110 solver.getIterationManager().getIterations(), enorm <= etol);
111 }
112
113 @Test
114 public void testSolveSaunders1() {
115 saundersTest(1, false, false, 0., 0.);
116 }
117
118 @Test
119 public void testSolveSaunders2() {
120 saundersTest(2, false, false, 0., 0.);
121 }
122
123 @Test
124 public void testSolveSaunders3() {
125 saundersTest(1, false, true, 0., 0.);
126 }
127
128 @Test
129 public void testSolveSaunders4() {
130 saundersTest(2, false, true, 0., 0.);
131 }
132
133 @Test
134 public void testSolveSaunders5() {
135 saundersTest(5, false, true, 0., 0.);
136 }
137
138 @Test
139 public void testSolveSaunders6() {
140 saundersTest(5, false, true, 0.25, 0.);
141 }
142
143 @Test
144 public void testSolveSaunders7() {
145 saundersTest(50, false, false, 0., 0.);
146 }
147
148 @Test
149 public void testSolveSaunders8() {
150 saundersTest(50, false, false, 0.25, 0.);
151 }
152
153 @Test
154 public void testSolveSaunders9() {
155 saundersTest(50, false, true, 0., 0.10);
156 }
157
158 @Test
159 public void testSolveSaunders10() {
160 saundersTest(50, false, true, 0.25, 0.10);
161 }
162
163 @Test
164 public void testSolveSaunders11() {
165 saundersTest(1, true, false, 0., 0.);
166 }
167
168 @Test
169 public void testSolveSaunders12() {
170 saundersTest(2, true, false, 0., 0.);
171 }
172
173 @Test
174 public void testSolveSaunders13() {
175 saundersTest(1, true, true, 0., 0.);
176 }
177
178 @Test
179 public void testSolveSaunders14() {
180 saundersTest(2, true, true, 0., 0.);
181 }
182
183 @Test
184 public void testSolveSaunders15() {
185 saundersTest(5, true, true, 0., 0.);
186 }
187
188 @Test
189 public void testSolveSaunders16() {
190 saundersTest(5, true, true, 0.25, 0.);
191 }
192
193 @Test
194 public void testSolveSaunders17() {
195 saundersTest(50, true, false, 0., 0.);
196 }
197
198 @Test
199 public void testSolveSaunders18() {
200 saundersTest(50, true, false, 0.25, 0.);
201 }
202
203 @Test
204 public void testSolveSaunders19() {
205 saundersTest(50, true, true, 0., 0.10);
206 }
207
208 @Test
209 public void testSolveSaunders20() {
210 saundersTest(50, true, true, 0.25, 0.10);
211 }
212
213 @Test(expected = NonSquareOperatorException.class)
214 public void testNonSquareOperator() {
215 final Array2DRowRealMatrix a = new Array2DRowRealMatrix(2, 3);
216 final IterativeLinearSolver solver;
217 solver = new SymmLQ(10, 0., false);
218 final ArrayRealVector b = new ArrayRealVector(a.getRowDimension());
219 final ArrayRealVector x = new ArrayRealVector(a.getColumnDimension());
220 solver.solve(a, b, x);
221 }
222
223 @Test(expected = DimensionMismatchException.class)
224 public void testDimensionMismatchRightHandSide() {
225 final Array2DRowRealMatrix a = new Array2DRowRealMatrix(3, 3);
226 final IterativeLinearSolver solver;
227 solver = new SymmLQ(10, 0., false);
228 final ArrayRealVector b = new ArrayRealVector(2);
229 solver.solve(a, b);
230 }
231
232 @Test(expected = DimensionMismatchException.class)
233 public void testDimensionMismatchSolution() {
234 final Array2DRowRealMatrix a = new Array2DRowRealMatrix(3, 3);
235 final IterativeLinearSolver solver;
236 solver = new SymmLQ(10, 0., false);
237 final ArrayRealVector b = new ArrayRealVector(3);
238 final ArrayRealVector x = new ArrayRealVector(2);
239 solver.solve(a, b, x);
240 }
241
242 @Test
243 public void testUnpreconditionedSolution() {
244 final int n = 5;
245 final int maxIterations = 100;
246 final RealLinearOperator a = new HilbertMatrix(n);
247 final InverseHilbertMatrix ainv = new InverseHilbertMatrix(n);
248 final IterativeLinearSolver solver;
249 solver = new SymmLQ(maxIterations, 1E-10, true);
250 final RealVector b = new ArrayRealVector(n);
251 for (int j = 0; j < n; j++) {
252 b.set(0.);
253 b.setEntry(j, 1.);
254 final RealVector x = solver.solve(a, b);
255 for (int i = 0; i < n; i++) {
256 final double actual = x.getEntry(i);
257 final double expected = ainv.getEntry(i, j);
258 final double delta = 1E-6 * JdkMath.abs(expected);
259 final String msg = String.format("entry[%d][%d]", i, j);
260 Assert.assertEquals(msg, expected, actual, delta);
261 }
262 }
263 }
264
265 @Test
266 public void testUnpreconditionedInPlaceSolutionWithInitialGuess() {
267 final int n = 5;
268 final int maxIterations = 100;
269 final RealLinearOperator a = new HilbertMatrix(n);
270 final InverseHilbertMatrix ainv = new InverseHilbertMatrix(n);
271 final IterativeLinearSolver solver;
272 solver = new SymmLQ(maxIterations, 1E-10, true);
273 final RealVector b = new ArrayRealVector(n);
274 for (int j = 0; j < n; j++) {
275 b.set(0.);
276 b.setEntry(j, 1.);
277 final RealVector x0 = new ArrayRealVector(n);
278 x0.set(1.);
279 final RealVector x = solver.solveInPlace(a, b, x0);
280 Assert.assertSame("x should be a reference to x0", x0, x);
281 for (int i = 0; i < n; i++) {
282 final double actual = x.getEntry(i);
283 final double expected = ainv.getEntry(i, j);
284 final double delta = 1E-6 * JdkMath.abs(expected);
285 final String msg = String.format("entry[%d][%d)", i, j);
286 Assert.assertEquals(msg, expected, actual, delta);
287 }
288 }
289 }
290
291 @Test
292 public void testUnpreconditionedSolutionWithInitialGuess() {
293 final int n = 5;
294 final int maxIterations = 100;
295 final RealLinearOperator a = new HilbertMatrix(n);
296 final InverseHilbertMatrix ainv = new InverseHilbertMatrix(n);
297 final IterativeLinearSolver solver;
298 solver = new SymmLQ(maxIterations, 1E-10, true);
299 final RealVector b = new ArrayRealVector(n);
300 for (int j = 0; j < n; j++) {
301 b.set(0.);
302 b.setEntry(j, 1.);
303 final RealVector x0 = new ArrayRealVector(n);
304 x0.set(1.);
305 final RealVector x = solver.solve(a, b, x0);
306 Assert.assertNotSame("x should not be a reference to x0", x0, x);
307 for (int i = 0; i < n; i++) {
308 final double actual = x.getEntry(i);
309 final double expected = ainv.getEntry(i, j);
310 final double delta = 1E-6 * JdkMath.abs(expected);
311 final String msg = String.format("entry[%d][%d]", i, j);
312 Assert.assertEquals(msg, expected, actual, delta);
313 Assert.assertEquals(msg, x0.getEntry(i), 1., Math.ulp(1.));
314 }
315 }
316 }
317
318 @Test(expected = NonSquareOperatorException.class)
319 public void testNonSquarePreconditioner() {
320 final Array2DRowRealMatrix a = new Array2DRowRealMatrix(2, 2);
321 final RealLinearOperator m = new RealLinearOperator() {
322
323 @Override
324 public RealVector operate(final RealVector x) {
325 throw new UnsupportedOperationException();
326 }
327
328 @Override
329 public int getRowDimension() {
330 return 2;
331 }
332
333 @Override
334 public int getColumnDimension() {
335 return 3;
336 }
337 };
338 final PreconditionedIterativeLinearSolver solver;
339 solver = new SymmLQ(10, 0., false);
340 final ArrayRealVector b = new ArrayRealVector(a.getRowDimension());
341 solver.solve(a, m, b);
342 }
343
344 @Test(expected = DimensionMismatchException.class)
345 public void testMismatchedOperatorDimensions() {
346 final Array2DRowRealMatrix a = new Array2DRowRealMatrix(2, 2);
347 final RealLinearOperator m = new RealLinearOperator() {
348
349 @Override
350 public RealVector operate(final RealVector x) {
351 throw new UnsupportedOperationException();
352 }
353
354 @Override
355 public int getRowDimension() {
356 return 3;
357 }
358
359 @Override
360 public int getColumnDimension() {
361 return 3;
362 }
363 };
364 final PreconditionedIterativeLinearSolver solver;
365 solver = new SymmLQ(10, 0d, false);
366 final ArrayRealVector b = new ArrayRealVector(a.getRowDimension());
367 solver.solve(a, m, b);
368 }
369
370 @Test(expected = NonPositiveDefiniteOperatorException.class)
371 public void testNonPositiveDefinitePreconditioner() {
372 final Array2DRowRealMatrix a = new Array2DRowRealMatrix(2, 2);
373 a.setEntry(0, 0, 1d);
374 a.setEntry(0, 1, 2d);
375 a.setEntry(1, 0, 3d);
376 a.setEntry(1, 1, 4d);
377 final RealLinearOperator m = new RealLinearOperator() {
378
379 @Override
380 public RealVector operate(final RealVector x) {
381 final ArrayRealVector y = new ArrayRealVector(2);
382 y.setEntry(0, -x.getEntry(0));
383 y.setEntry(1, -x.getEntry(1));
384 return y;
385 }
386
387 @Override
388 public int getRowDimension() {
389 return 2;
390 }
391
392 @Override
393 public int getColumnDimension() {
394 return 2;
395 }
396 };
397 final PreconditionedIterativeLinearSolver solver;
398 solver = new SymmLQ(10, 0d, true);
399 final ArrayRealVector b = new ArrayRealVector(2);
400 b.setEntry(0, -1d);
401 b.setEntry(1, -1d);
402 solver.solve(a, m, b);
403 }
404
405 @Test
406 public void testPreconditionedSolution() {
407 final int n = 8;
408 final int maxIterations = 100;
409 final RealLinearOperator a = new HilbertMatrix(n);
410 final InverseHilbertMatrix ainv = new InverseHilbertMatrix(n);
411 final RealLinearOperator m = JacobiPreconditioner.create(a);
412 final PreconditionedIterativeLinearSolver solver;
413 solver = new SymmLQ(maxIterations, 1E-15, true);
414 final RealVector b = new ArrayRealVector(n);
415 for (int j = 0; j < n; j++) {
416 b.set(0.);
417 b.setEntry(j, 1.);
418 final RealVector x = solver.solve(a, m, b);
419 for (int i = 0; i < n; i++) {
420 final double actual = x.getEntry(i);
421 final double expected = ainv.getEntry(i, j);
422 final double delta = 1E-6 * JdkMath.abs(expected);
423 final String msg = String.format("coefficient (%d, %d)", i, j);
424 Assert.assertEquals(msg, expected, actual, delta);
425 }
426 }
427 }
428
429 @Test
430 public void testPreconditionedSolution2() {
431 final int n = 100;
432 final int maxIterations = 100000;
433 final Array2DRowRealMatrix a = new Array2DRowRealMatrix(n, n);
434 double daux = 1.;
435 for (int i = 0; i < n; i++) {
436 a.setEntry(i, i, daux);
437 daux *= 1.2;
438 for (int j = i + 1; j < n; j++) {
439 if (i != j) {
440 final double value = 1.0;
441 a.setEntry(i, j, value);
442 a.setEntry(j, i, value);
443 }
444 }
445 }
446 final RealLinearOperator m = JacobiPreconditioner.create(a);
447 final PreconditionedIterativeLinearSolver prec;
448 final IterativeLinearSolver unprec;
449 prec = new SymmLQ(maxIterations, 1E-15, true);
450 unprec = new SymmLQ(maxIterations, 1E-15, true);
451 final RealVector b = new ArrayRealVector(n);
452 final String pattern = "preconditioned SymmLQ (%d iterations) should"
453 + " have been faster than unpreconditioned (%d iterations)";
454 String msg;
455 for (int j = 0; j < 1; j++) {
456 b.set(0.);
457 b.setEntry(j, 1.);
458 final RealVector px = prec.solve(a, m, b);
459 final RealVector x = unprec.solve(a, b);
460 final int np = prec.getIterationManager().getIterations();
461 final int nup = unprec.getIterationManager().getIterations();
462 msg = String.format(pattern, np, nup);
463 for (int i = 0; i < n; i++) {
464 msg = String.format("row %d, column %d", i, j);
465 final double expected = x.getEntry(i);
466 final double actual = px.getEntry(i);
467 final double delta = 5E-5 * JdkMath.abs(expected);
468 Assert.assertEquals(msg, expected, actual, delta);
469 }
470 }
471 }
472
473 @Test
474 public void testEventManagement() {
475 final int n = 5;
476 final int maxIterations = 100;
477 final RealLinearOperator a = new HilbertMatrix(n);
478 final IterativeLinearSolver solver;
479
480
481
482
483
484
485 final int[] count = new int[] {0, 0, 0, 0};
486 final RealVector xFromListener = new ArrayRealVector(n);
487 final IterationListener listener = new IterationListener() {
488
489 @Override
490 public void initializationPerformed(final IterationEvent e) {
491 ++count[0];
492 }
493
494 @Override
495 public void iterationPerformed(final IterationEvent e) {
496 ++count[2];
497 Assert.assertEquals("iteration performed",
498 count[2],
499 e.getIterations() - 1);
500 }
501
502 @Override
503 public void iterationStarted(final IterationEvent e) {
504 ++count[1];
505 Assert.assertEquals("iteration started",
506 count[1],
507 e.getIterations() - 1);
508 }
509
510 @Override
511 public void terminationPerformed(final IterationEvent e) {
512 ++count[3];
513 final IterativeLinearSolverEvent ilse;
514 ilse = (IterativeLinearSolverEvent) e;
515 xFromListener.setSubVector(0, ilse.getSolution());
516 }
517 };
518 solver = new SymmLQ(maxIterations, 1E-10, true);
519 solver.getIterationManager().addIterationListener(listener);
520 final RealVector b = new ArrayRealVector(n);
521 for (int j = 0; j < n; j++) {
522 Arrays.fill(count, 0);
523 b.set(0.);
524 b.setEntry(j, 1.);
525 final RealVector xFromSolver = solver.solve(a, b);
526 String msg = String.format("column %d (initialization)", j);
527 Assert.assertEquals(msg, 1, count[0]);
528 msg = String.format("column %d (finalization)", j);
529 Assert.assertEquals(msg, 1, count[3]);
530
531
532
533
534
535 for (int i = 0; i < n; i++){
536 msg = String.format("row %d, column %d", i, j);
537 final double expected = xFromSolver.getEntry(i);
538 final double actual = xFromListener.getEntry(i);
539 Assert.assertEquals(msg, expected, actual, 0.0);
540 }
541 }
542 }
543
544 @Test(expected = NonSelfAdjointOperatorException.class)
545 public void testNonSelfAdjointOperator() {
546 final RealLinearOperator a;
547 a = new Array2DRowRealMatrix(new double[][] {
548 {1., 2., 3.},
549 {2., 4., 5.},
550 {2.999, 5., 6.}
551 });
552 final RealVector b;
553 b = new ArrayRealVector(new double[] {1., 1., 1.});
554 new SymmLQ(100, 1., true).solve(a, b);
555 }
556
557 @Test(expected = NonSelfAdjointOperatorException.class)
558 public void testNonSelfAdjointPreconditioner() {
559 final RealLinearOperator a = new Array2DRowRealMatrix(new double[][] {
560 {1., 2., 3.},
561 {2., 4., 5.},
562 {3., 5., 6.}
563 });
564 final Array2DRowRealMatrix mMat;
565 mMat = new Array2DRowRealMatrix(new double[][] {
566 {1., 0., 1.},
567 {0., 1., 0.},
568 {0., 0., 1.}
569 });
570 final DecompositionSolver mSolver;
571 mSolver = new LUDecomposition(mMat).getSolver();
572 final RealLinearOperator minv = new RealLinearOperator() {
573
574 @Override
575 public RealVector operate(final RealVector x) {
576 return mSolver.solve(x);
577 }
578
579 @Override
580 public int getRowDimension() {
581 return mMat.getRowDimension();
582 }
583
584 @Override
585 public int getColumnDimension() {
586 return mMat.getColumnDimension();
587 }
588 };
589 final RealVector b = new ArrayRealVector(new double[] {
590 1., 1., 1.
591 });
592 new SymmLQ(100, 1., true).solve(a, minv, b);
593 }
594
595 @Test
596 public void testUnpreconditionedNormOfResidual() {
597 final int n = 5;
598 final int maxIterations = 100;
599 final RealLinearOperator a = new HilbertMatrix(n);
600 final IterativeLinearSolver solver;
601 final IterationListener listener = new IterationListener() {
602
603 private void doTestNormOfResidual(final IterationEvent e) {
604 final IterativeLinearSolverEvent evt;
605 evt = (IterativeLinearSolverEvent) e;
606 final RealVector x = evt.getSolution();
607 final RealVector b = evt.getRightHandSideVector();
608 final RealVector r = b.subtract(a.operate(x));
609 final double rnorm = r.getNorm();
610 Assert.assertEquals("iteration performed (residual)",
611 rnorm, evt.getNormOfResidual(),
612 JdkMath.max(1E-5 * rnorm, 1E-10));
613 }
614
615 @Override
616 public void initializationPerformed(final IterationEvent e) {
617 doTestNormOfResidual(e);
618 }
619
620 @Override
621 public void iterationPerformed(final IterationEvent e) {
622 doTestNormOfResidual(e);
623 }
624
625 @Override
626 public void iterationStarted(final IterationEvent e) {
627 doTestNormOfResidual(e);
628 }
629
630 @Override
631 public void terminationPerformed(final IterationEvent e) {
632 doTestNormOfResidual(e);
633 }
634 };
635 solver = new SymmLQ(maxIterations, 1E-10, true);
636 solver.getIterationManager().addIterationListener(listener);
637 final RealVector b = new ArrayRealVector(n);
638 for (int j = 0; j < n; j++) {
639 b.set(0.);
640 b.setEntry(j, 1.);
641 solver.solve(a, b);
642 }
643 }
644
645 @Test
646 public void testPreconditionedNormOfResidual() {
647 final int n = 5;
648 final int maxIterations = 100;
649 final RealLinearOperator a = new HilbertMatrix(n);
650 final JacobiPreconditioner m = JacobiPreconditioner.create(a);
651 final RealLinearOperator p = m.sqrt();
652 final PreconditionedIterativeLinearSolver solver;
653 final IterationListener listener = new IterationListener() {
654
655 private void doTestNormOfResidual(final IterationEvent e) {
656 final IterativeLinearSolverEvent evt;
657 evt = (IterativeLinearSolverEvent) e;
658 final RealVector x = evt.getSolution();
659 final RealVector b = evt.getRightHandSideVector();
660 final RealVector r = b.subtract(a.operate(x));
661 final double rnorm = p.operate(r).getNorm();
662 Assert.assertEquals("iteration performed (residual)",
663 rnorm, evt.getNormOfResidual(),
664 JdkMath.max(1E-5 * rnorm, 1E-10));
665 }
666
667 @Override
668 public void initializationPerformed(final IterationEvent e) {
669 doTestNormOfResidual(e);
670 }
671
672 @Override
673 public void iterationPerformed(final IterationEvent e) {
674 doTestNormOfResidual(e);
675 }
676
677 @Override
678 public void iterationStarted(final IterationEvent e) {
679 doTestNormOfResidual(e);
680 }
681
682 @Override
683 public void terminationPerformed(final IterationEvent e) {
684 doTestNormOfResidual(e);
685 }
686 };
687 solver = new SymmLQ(maxIterations, 1E-10, true);
688 solver.getIterationManager().addIterationListener(listener);
689 final RealVector b = new ArrayRealVector(n);
690 for (int j = 0; j < n; j++) {
691 b.set(0.);
692 b.setEntry(j, 1.);
693 solver.solve(a, m, b);
694 }
695 }
696 }
697