MyraMath
normal2_solver


Source: tests/multifrontal/normal/normal2_solver.cpp

1 // ========================================================================= //
2 // This file is part of MyraMath, copyright (c) 2014-2019 by Ryan A Chilton //
3 // and distributed by MyraCore, LLC. See LICENSE.txt for license terms. //
4 // ========================================================================= //
5 
11 // Containers.
12 #include <myramath/utility/Number.h>
13 #include <myramath/dense/Vector.h>
14 #include <myramath/dense/VectorRange.h>
15 #include <myramath/dense/Matrix.h>
16 #include <myramath/dense/MatrixRange.h>
17 #include <myramath/sparse/Pattern.h>
18 #include <myramath/sparse/PatternRange.h>
19 #include <myramath/sparse/SparseMatrix.h>
20 #include <myramath/sparse/SparseMatrixRange.h>
21 
22 // Algorithms.
23 #include <myramath/dense/frobenius.h>
24 #include <myramath/sparse/laplacian2.h>
25 #include <myramath/sparse/transpose.h>
26 #include <myramath/sparse/hermitian.h>
27 #include <myramath/sparse/conjugate.h>
28 #include <myramath/sparse/gemm.h>
29 #include <myramath/sparse/gemv.h>
30 
31 // Solver under test.
32 #include <myramath/multifrontal/normal/SparseNormalSolver.h>
33 
34 // For testing serialization roundtrip.
35 #include <myramath/io/VectorStream.h>
36 
37 // Reporting.
38 #include <myramath/utility/stlprint.h>
39 #include <tests/myratest.h>
40 
41 using namespace myra;
42 using namespace myra_stlprint;
43 
44 namespace {
45 
46 template<class Number> void test_solve(int I, int J, typename ReflectPrecision<Number>::type tolerance)
47  {
48  myra::out() << typestring<Number>() << std::endl;
49  typedef typename ReflectPrecision<Number>::type Precision;
50  // Make random A with unsymmetric pattern.
51  auto P = stencil2_unsymmetric(I,J);
52  auto A = SparseMatrix<Number>::random(P);
53  int N = I*J;
54  // Make SparseNormalSolver.
55  typedef SparseNormalSolver<Number> Solver;
56  Solver solver(A);
57  // Test A*x = b.
58  {
59  auto b = Vector<Number>::random(N);
60  auto x = b;
61  solver.solve(x.column(),'L','N');
62  Precision residual = frobenius(A*x-b) / frobenius(b);
63  myra::out() << " |A*x-b| (solve) = " << residual << std::endl;
64  }
65  // Test A^T*x = b.
66  {
67  auto b = Vector<Number>::random(N);
68  auto x = b;
69  solver.solve(x.column(),'L','T');
70  Precision residual = frobenius(transpose(A)*x-b) / frobenius(b);
71  myra::out() << " |A^T*x-b| (solve) = " << residual << std::endl;
72  }
73  // Test A^H*x = b.
74  {
75  auto b = Vector<Number>::random(N);
76  auto x = b;
77  solver.solve(x.column(),'L','H');
78  Precision residual = frobenius(hermitian(A)*x-b) / frobenius(b);
79  myra::out() << " |A^H*x-b| (solve) = " << residual << std::endl;
80  }
81  // Test A^C*x = b.
82  {
83  auto b = Vector<Number>::random(N);
84  auto x = b;
85  solver.solve(x.column(),'L','C');
86  Precision residual = frobenius(conjugate(A)*x-b) / frobenius(b);
87  myra::out() << " |A^C*x-b| (solve) = " << residual << std::endl;
88  }
89  // Test x*A = b.
90  {
91  auto b = Vector<Number>::random(N);
92  auto x = b;
93  solver.solve(x.row(),'R','N');
94  Precision residual = frobenius(x*A-b) / frobenius(b);
95  myra::out() << " |x*A-b| (solve) = " << residual << std::endl;
96  }
97  // Test x*A^T = b.
98  {
99  auto b = Vector<Number>::random(N);
100  auto x = b;
101  solver.solve(x.row(),'R','T');
102  Precision residual = frobenius(x*transpose(A)-b) / frobenius(b);
103  myra::out() << " |x*A^T-b| (solve) = " << residual << std::endl;
104  }
105  // Test x*A^H = b.
106  {
107  auto b = Vector<Number>::random(N);
108  auto x = b;
109  solver.solve(x.row(),'R','H');
110  Precision residual = frobenius(x*hermitian(A)-b) / frobenius(b);
111  myra::out() << " |x*A^H-b| (solve) = " << residual << std::endl;
112  }
113  // Test x*A^C = b.
114  {
115  auto b = Vector<Number>::random(N);
116  auto x = b;
117  solver.solve(x.row(),'R','C');
118  Precision residual = frobenius(x*conjugate(A)-b) / frobenius(b);
119  myra::out() << " |x*A^C-b| (solve) = " << residual << std::endl;
120  }
121  }
122 
123 template<class Number> void test_refine(int I, int J, typename ReflectPrecision<Number>::type tolerance)
124  {
125  myra::out() << typestring<Number>() << std::endl;
126  typedef typename ReflectPrecision<Number>::type Precision;
127  // Make random A with unsymmetric pattern.
128  auto P = stencil2_unsymmetric(I,J);
129  auto A = SparseMatrix<Number>::random(P);
130  int N = I*J;
131  // Make SparseNormalSolver.
132  typedef SparseNormalSolver<Number> Solver;
133  Solver solver(A);
134  // Test A*x = b, using refine.
135  {
136  auto b = Vector<Number>::random(N);
137  auto x = b;
138  auto history = solver.refine(x.column(),'L','N');
139  Precision residual = frobenius(A*x-b) / frobenius(b);
140  myra::out() << " |A*x-b| (refine) = " << residual << std::endl;
141  myra::out() << " history = " << history << std::endl;
142  // Note, we do check the residual on this one (and all the rest)
143  REQUIRE(residual < tolerance);
144  }
145  // Test A^T*x = b, using refine.
146  {
147  auto b = Vector<Number>::random(N);
148  auto x = b;
149  auto history = solver.refine(x.column(),'L','T');
150  Precision residual = frobenius(transpose(A)*x-b) / frobenius(b);
151  myra::out() << " |A^T*x-b| (refine) = " << residual << std::endl;
152  myra::out() << " history = " << history << std::endl;
153  // Note, we do check the residual on this one (and all the rest)
154  REQUIRE(residual < tolerance);
155  }
156  // Test A^H*x = b, using refine.
157  {
158  auto b = Vector<Number>::random(N);
159  auto x = b;
160  auto history = solver.refine(x.column(),'L','H');
161  Precision residual = frobenius(hermitian(A)*x-b) / frobenius(b);
162  myra::out() << " |A^H*x-b| (refine) = " << residual << std::endl;
163  myra::out() << " history = " << history << std::endl;
164  // Note, we do check the residual on this one (and all the rest)
165  REQUIRE(residual < tolerance);
166  }
167  // Test A^C*x = b, using refine.
168  {
169  auto b = Vector<Number>::random(N);
170  auto x = b;
171  auto history = solver.refine(x.column(),'L','C');
172  Precision residual = frobenius(conjugate(A)*x-b) / frobenius(b);
173  myra::out() << " |A^C*x-b| (refine) = " << residual << std::endl;
174  myra::out() << " history = " << history << std::endl;
175  // Note, we do check the residual on this one (and all the rest)
176  REQUIRE(residual < tolerance);
177  }
178  // Test x*A = b, using refine.
179  {
180  auto b = Vector<Number>::random(N);
181  auto x = b;
182  auto history = solver.refine(x.row(),'R','N');
183  Precision residual = frobenius(x*A-b) / frobenius(b);
184  myra::out() << " |x*A-b| (refine) = " << residual << std::endl;
185  myra::out() << " history = " << history << std::endl;
186  // Note, we do check the residual on this one (and all the rest)
187  REQUIRE(residual < tolerance);
188  }
189  // Test x*A^T = b, using refine.
190  {
191  auto b = Vector<Number>::random(N);
192  auto x = b;
193  auto history = solver.refine(x.row(),'R','T');
194  Precision residual = frobenius(x*transpose(A)-b) / frobenius(b);
195  myra::out() << " |x*A^T-b| (refine) = " << residual << std::endl;
196  myra::out() << " history = " << history << std::endl;
197  // Note, we do check the residual on this one (and all the rest)
198  REQUIRE(residual < tolerance);
199  }
200  // Test x*A^H = b, using refine.
201  {
202  auto b = Vector<Number>::random(N);
203  auto x = b;
204  auto history = solver.refine(x.row(),'R','H');
205  Precision residual = frobenius(x*hermitian(A)-b) / frobenius(b);
206  myra::out() << " |x*A^H-b| (refine) = " << residual << std::endl;
207  myra::out() << " history = " << history << std::endl;
208  // Note, we do check the residual on this one (and all the rest)
209  REQUIRE(residual < tolerance);
210  }
211  // Test x*A^C = b, using refine.
212  {
213  auto b = Vector<Number>::random(N);
214  auto x = b;
215  auto history = solver.refine(x.row(),'R','C');
216  Precision residual = frobenius(x*conjugate(A)-b) / frobenius(b);
217  myra::out() << " |x*A^C-b| (refine) = " << residual << std::endl;
218  myra::out() << " history = " << history << std::endl;
219  // Note, we do check the residual on this one (and all the rest)
220  REQUIRE(residual < tolerance);
221  }
222  // Make a deep copy of solver, test A*x = b again.
223  {
224  Solver copy = solver;
225  auto b = Vector<Number>::random(N);
226  auto x = b;
227  auto history = solver.refine(x.column(),'L','N');
228  Precision residual = frobenius(A*x-b) / frobenius(b);
229  myra::out() << " |inv(A)*b-x| (copy) = " << residual << std::endl;
230  REQUIRE(residual < tolerance);
231  }
232  // Roundtrip solver to file, test A*x = b again.
233  {
234  VectorStream inout;
235  solver.write(inout);
236  Solver saved(inout);
237  auto b = Vector<Number>::random(N);
238  auto x = b;
239  auto history = solver.refine(x.column(),'L','N');
240  Precision residual = frobenius(A*x-b) / frobenius(b);
241  myra::out() << " |A*x-b| (saved) = " << residual << std::endl;
242  REQUIRE(residual < tolerance);
243  }
244  }
245 
246 } // namespace
247 
248 ADD_TEST("normal2_solver","[multifrontal][parallel]")
249  {
250  int I = 64;
251  int J = 64;
252  // Note single precision SparseNormalSolver is not very reliable.
253  test_solve<NumberD>(I,J,1.0e-8);
254  test_solve<NumberZ>(I,J,1.0e-8);
255  test_refine<NumberD>(I,J,1.0e-8);
256  test_refine<NumberZ>(I,J,1.0e-8);
257  }
transpose.h
Returns a transposed copy of a SparseMatrix.
MatrixRange.h
Interface class for representing subranges of dense Matrix&#39;s.
VectorRange.h
Interface class for representing subranges of dense Vector&#39;s.
gemm.h
Variety of routines for mixed dense*sparse or dense*sparse matrix multiplies. The dense*dense case is...
frobenius.h
Routines for computing Frobenius norms of various algebraic containers.
myra::SparseMatrix::random
static SparseMatrix< Number > random(int I, int J, int N)
Generates a random SparseMatrix with size IxJ and (approximately) N nonzeros.
Definition: SparseMatrix.cpp:493
myra::VectorStream
Wraps a std::vector<char>, presents it as both an InputStream and OutputStream. Useful for hygienic u...
Definition: VectorStream.h:22
SparseMatrix.h
General purpose compressed-sparse-column (CSC) container.
myra::Vector::random
static Vector< Number > random(int N)
Generates a random Vector of specified size.
Definition: Vector.cpp:276
myra
Definition: syntax.dox:1
myra::SparseNormalSolver
Definition: SparseNormalSolver.h:65
myra_stlprint
Definition: stlprint.h:32
stlprint.h
Routines for printing the contents of various std::container&#39;s to a std::ostream using operator <<...
Number.h
Various utility functions/classes related to scalar Number types.
gemv.h
Signatures for sparse matrix * dense vector multiplies. All delegate to gemm() under the hood...
VectorStream.h
A stream that serialize/deserializes to std::vector<char> buffer.
Matrix.h
General purpose dense matrix container, O(i*j) storage.
PatternRange.h
Range/Iterator types associated with Pattern.
Vector.h
Container for either a column vector or row vector (depends upon the usage context) ...
SparseNormalSolver.h
Sparse direct solver for general systems, solves the normal equations.
myra::ReflectPrecision
Reflects Precision trait for a Number, scalar Number types should specialize it.
Definition: Number.h:33
Pattern.h
Container class for a sparse nonzero pattern, used in reordering/symbolic analysis.
conjugate.h
Returns a conjugated copy of a SparseMatrix.
hermitian.h
Returns a hermitian copy of a SparseMatrix.
laplacian2.h
Helper routines for reordering/filling 2D structured grids. Used by many unit tests.
SparseMatrixRange.h
Range/Iterator types associated with SparseMatrix.


Results: [PASS]

double
|A*x-b| (solve) = 1.25194e-08
|A^T*x-b| (solve) = 4.09645e-08
|A^H*x-b| (solve) = 1.41502e-07
|A^C*x-b| (solve) = 8.50439e-09
|x*A-b| (solve) = 1.0355e-07
|x*A^T-b| (solve) = 3.16521e-09
|x*A^H-b| (solve) = 7.73489e-09
|x*A^C-b| (solve) = 2.46258e-08
std::complex<double>
|A*x-b| (solve) = 1.43095e-12
|A^T*x-b| (solve) = 9.9526e-11
|A^H*x-b| (solve) = 1.2041e-10
|A^C*x-b| (solve) = 1.45929e-12
|x*A-b| (solve) = 2.96004e-11
|x*A^T-b| (solve) = 1.53725e-12
|x*A^H-b| (solve) = 1.19461e-12
|x*A^C-b| (solve) = 1.23877e-10
double
|A*x-b| (refine) = 1.82258e-12
history = [ 1.19777e-07 1.82258e-12 ] (2)
|A^T*x-b| (refine) = 9.20735e-12
history = [ 1.55939e-06 9.20735e-12 ] (2)
|A^H*x-b| (refine) = 2.86293e-12
history = [ 4.47109e-07 2.86293e-12 ] (2)
|A^C*x-b| (refine) = 4.75954e-13
history = [ 3.07473e-08 4.75954e-13 ] (2)
|x*A-b| (refine) = 1.4148e-11
history = [ 2.05893e-06 1.4148e-11 ] (2)
|x*A^T-b| (refine) = 2.36662e-13
history = [ 1.62772e-08 2.36662e-13 ] (2)
|x*A^H-b| (refine) = 8.45199e-13
history = [ 6.86263e-08 8.45199e-13 ] (2)
|x*A^C-b| (refine) = 1.4028e-12
history = [ 2.46126e-07 1.4028e-12 ] (2)
|inv(A)*b-x| (copy) = 2.10464e-12
|A*x-b| (saved) = 1.03203e-13
std::complex<double>
|A*x-b| (refine) = 1.20593e-12
history = [ 1.20593e-12 ] (1)
|A^T*x-b| (refine) = 8.96022e-15
history = [ 7.30743e-11 8.96022e-15 ] (2)
|A^H*x-b| (refine) = 1.31296e-14
history = [ 9.92012e-11 1.31296e-14 ] (2)
|A^C*x-b| (refine) = 8.36052e-15
history = [ 2.64522e-12 8.36052e-15 ] (2)
|x*A-b| (refine) = 9.33219e-15
history = [ 6.61999e-11 9.33219e-15 ] (2)
|x*A^T-b| (refine) = 1.81592e-12
history = [ 1.81592e-12 ] (1)
|x*A^H-b| (refine) = 1.31926e-12
history = [ 1.31926e-12 ] (1)
|x*A^C-b| (refine) = 1.52026e-14
history = [ 1.27944e-10 1.52026e-14 ] (2)
|inv(A)*b-x| (copy) = 1.40238e-14
|A*x-b| (saved) = 1.08513e-14


Go back to Summary of /test programs.

Generated on Sun Aug 4 2024 12:26:14 for MyraMath by   doxygen 1.8.13