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