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