MyraMath
minres_complex


Source: tests/iterative/minres.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/utility/eprintf.h>
14 #include <myramath/expression/Expression.h>
15 #include <myramath/dense/Vector.h>
16 #include <myramath/dense/VectorRange.h>
17 #include <myramath/dense/Matrix.h>
18 #include <myramath/dense/MatrixRange.h>
19 #include <myramath/dense/DiagonalMatrix.h>
20 #include <myramath/dense/DiagonalMatrixRange.h>
21 #include <myramath/sparse/SparseMatrix.h>
22 #include <myramath/sparse/SparseMatrixRange.h>
23 #include <myramath/sparse/SparseMatrixBuilder.h>
24 
25 // Algorithms.
26 #include <myramath/expression/generators.h>
27 #include <myramath/dense/gemv.h>
28 #include <myramath/dense/gemm.h>
29 #include <myramath/dense/dimm.h>
30 #include <myramath/dense/euclidean.h>
31 #include <myramath/dense/frobenius.h>
32 #include <myramath/sparse/diag.h>
33 #include <myramath/sparse/phasor.h>
34 
35 // Iterative solve/action stuff.
36 #include <myramath/iterative/Action.h>
37 #include <myramath/iterative/GemmAction.h>
38 #include <myramath/iterative/DimmAction.h>
39 #include <myramath/iterative/minres.h>
40 
41 // Reporting.
42 #include <myramath/utility/stlprint.h>
43 #include <tests/myratest.h>
44 
45 using namespace myra;
46 
47 namespace {
48 
49 // Implementation detail, makes tridiagonal indefinite input A for minres() tests below.
50 template<class Number> SparseMatrix<Number> make_minres_A(int N)
51  {
52  // Form A.
53  if ( (N % 2) != 0)
54  throw eprintf("make_minres_A(N), must have N even [%d]",N);
55  int K = N/2;
56  Number one(1);
57  Number dmin(5);
58  Number dmax(10);
59  auto D = DiagonalMatrix<Number>::evaluate( make_LinspaceExpression(dmax,dmin,K) );
60  SparseMatrixBuilder<Number> A_builder(N,N);
61  for (int k = 0; k < K; ++k)
62  {
63  A_builder(k,k) = -D(k);
64  A_builder(N-k-1,N-k-1) = D(k);
65  }
66  for (int n = 0; n < N-1; ++n)
67  {
68  A_builder(n,n+1) = -one;
69  A_builder(n+1,n) = -one;
70  }
71  return A_builder.make_SparseMatrix();
72  }
73 
74 // Tests minres() for real symmetric indefinite A.
75 template<class Precision> void test1(int N, Precision tol)
76  {
77  using myra_stlprint::operator<<;
78  // Form A.
79  myra::out() << typestring<Precision>() << std::endl;
80  auto A = make_minres_A<Precision>(N);
81  // Precondition using diagonal scaling. Note M must be positive.
82  auto M = diag(A);
83  Precision one(1);
84  for (int n = 0; n < N; ++n)
85  M(n) = std::abs(one/M(n));
86  // Make forcing data, zero valued initial guess.
87  auto b = Vector<Precision>::ones(N);
88  auto x = Vector<Precision>::zeros(N);
89  // Solve A*x=b using minres()
90  auto M_action = make_DimmAction(M);
91  auto A_action = make_GemmAction(A);
92  auto output1 = minres(M_action,A_action,b,x,tol);
93  // Examine residual.
94  Precision error1 = euclidean( M_action*(A_action*x-b) ) / euclidean(b);
95  myra::out() << "history = " << output1.history << std::endl;
96  myra::out() << "|M*(A*x-b)|/|b| = " << error1 << std::endl;
97  REQUIRE( error1 < tol );
98  // Solve A*X=B using minres()
99  auto B = Matrix<Precision>::random(N,10);
100  auto X = Matrix<Precision>::zeros(N,10);
101  auto output2 = minres(M_action,A_action,B,X,tol);
102  // Examine residual.
103  Precision error2 = frobenius( M_action*(A_action*X-B) ) / frobenius(B);
104  myra::out() << "history = " << output2.history << std::endl;
105  myra::out() << "|M*(A*X-B)|/|B| = " << error2 << std::endl;
106  REQUIRE( error2 < tol );
107  }
108 
109 // Tests minres() for complex hermitian indefinite A.
110 template<class Precision> void test2(int N, Precision tol)
111  {
112  // Form A, make it nontrivially hermitian.
113  typedef typename ReflectComplex<Precision>::type Number;
114  myra::out() << typestring<Number>() << std::endl;
115  auto A = make_minres_A<Number>(N);
116  phasor_hermitian(A);
117  // Precondition using diagonal scaling. Note M must be positive.
118  auto M = diag(A);
119  Number one(1);
120  for (int n = 0; n < N; ++n)
121  M(n) = std::abs(one/M(n));
122  // Make forcing data, zero valued initial guess.
123  auto b = Vector<Number>::ones(N);
124  auto x = Vector<Number>::zeros(N);
125  // Solve A*x=b using minres()
126  auto M_action = make_DimmAction(M);
127  auto A_action = make_GemmAction(A);
128  auto output1 = minres(M_action,A_action,b,x,tol);
129  // Examine residual.
130  Precision error1 = euclidean( M_action*(A_action*x-b) ) / euclidean(b);
131  using myra_stlprint::operator<<;
132  myra::out() << "history = " << output1.history << std::endl;
133  myra::out() << "|M*(A*x-b)|/|b| = " << error1 << std::endl;
134  REQUIRE( error1 < tol );
135  // Solve A*X=B using minres()
136  auto B = Matrix<Number>::random(N,10);
137  auto X = Matrix<Number>::zeros(N,10);
138  auto output2 = minres(M_action,A_action,B,X,tol);
139  // Examine residual.
140  Precision error2 = frobenius( M_action*(A_action*X-B) ) / frobenius(B);
141  myra::out() << "history = " << output2.history << std::endl;
142  myra::out() << "|M*(A*X-B)|/|B| = " << error2 << std::endl;
143  REQUIRE( error2 < tol );
144  }
145 
146 } // namespace
147 
148 ADD_TEST("minres_real","[iterative]")
149  {
150  test1<float>(100,1.0e-3f);
151  test1<double>(100,1.0e-6);
152  }
153 
154 ADD_TEST("minres_complex","[iterative]")
155  {
156  test2<float>(100,1.0e-3f);
157  test2<double>(100,1.0e-6);
158  }
159 
euclidean.h
Routines for computing euclidean norm of a Vector/VectorRange, or normalizing a Vector/VectorRange to...
DiagonalMatrixRange.h
Interface class for representing subranges of DiagonalMatrix&#39;s.
eprintf.h
Returns a std::runtime_error() whose message has been populated using printf()-style formatting...
MatrixRange.h
Interface class for representing subranges of dense Matrix&#39;s.
generators.h
Generators for basic Expression&#39;s (constant, random, linspace, etc).
myra::Vector::ones
static Vector< Number > ones(int N)
Generates a ones Vector of specified size.
Definition: Vector.cpp:284
VectorRange.h
Interface class for representing subranges of dense Vector&#39;s.
Action.h
Applies the "Action" of a linear operator, b := A*x, used in iterative solution algorithms.
myra::Matrix::zeros
static Matrix< Number > zeros(int I, int J)
Generates a zeros Matrix of specified size.
Definition: Matrix.cpp:357
frobenius.h
Routines for computing Frobenius norms of various algebraic containers.
myra::Matrix::random
static Matrix< Number > random(int I, int J)
Generates a random Matrix of specified size.
Definition: Matrix.cpp:353
dimm.h
Routines for multiplying by a DiagonalMatrix.
Expression.h
An interface used to fill containers from Expression&#39;s (see Matrix::evaluate(), for example)...
SparseMatrix.h
General purpose compressed-sparse-column (CSC) container.
gemv.h
Variety of routines all for dense Matrix*Vector multiplies. All just delegate to gemm() ...
myra
Definition: syntax.dox:1
myra::ReflectComplex
Reflects a Precision type into a complex type.
Definition: Number.h:46
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.
myra::Vector::zeros
static Vector< Number > zeros(int N)
Generates a zeros Vector of specified size.
Definition: Vector.cpp:280
Matrix.h
General purpose dense matrix container, O(i*j) storage.
Vector.h
Container for either a column vector or row vector (depends upon the usage context) ...
myra::SparseMatrixBuilder
Convenience type for building SparseMatrix&#39;s, uses coordinate/triplet format.
Definition: SparseMatrix.h:32
GemmAction.h
An Action for multiplying by a dense Matrix or SparseMatrix using gemm().
SparseMatrixBuilder.h
Convenience type for building SparseMatrix&#39;s, uses coordinate/triplet format. Note that SparseMatrixB...
diag.h
Returns the diagonal of a SparseMatrix.
DimmAction.h
An Action for multiplying by a DiagonalMatrix using dimm()
DiagonalMatrix.h
Container for a diagonal matrix, O(n) storage. Used by SVD, row/column scaling, etc.
myra::SparseMatrix
Stores an IxJ matrix A in compressed sparse column format.
Definition: bothcat.h:23
myra::DiagonalMatrix::evaluate
static DiagonalMatrix< Number > evaluate(const Expression< 1, Number > &e)
Generates a DiagonalMatrix by evaluating an arity-1 Expression of Number.
Definition: DiagonalMatrix.cpp:258
gemm.h
Variety of routines all for dense Matrix*Matrix multiplies. Delegates to the BLAS.
phasor.h
Applies random phase shifts to a complex square SparseMatrix.
minres.h
Linear system solution via minimum residual method (symmetric action A, can be indefinite) ...
SparseMatrixRange.h
Range/Iterator types associated with SparseMatrix.


Results: [PASS]

std::complex<float>
history = [ 0.370527 0.089801 0.0881976 0.0274258 0.0268242 0.0100926 0.00990909 0.00374722 0.00369557 0.0014238 0.00140456 0.000520099 ] (12)
|M*(A*x-b)|/|b| = 0.00037029
history = [ 0.355903 0.0608088 0.0605837 0.0202432 0.0201083 0.00733785 0.0072879 0.00280633 0.00279597 0.00104344 0.00103963 0.000370742 ] (12)
|M*(A*X-B)|/|B| = 0.000378027
std::complex<double>
history = [ 0.371143 0.090826 0.0897114 0.0311435 0.0303255 0.0118155 0.011417 0.00443448 0.00438151 0.00159859 0.00157597 0.000520615 0.000513226 0.000175812 0.000173959 5.63975e-05 5.54209e-05 1.78631e-05 1.76299e-05 6.36678e-06 6.31986e-06 2.16545e-06 2.11853e-06 6.80687e-07 ] (24)
|M*(A*x-b)|/|b| = 5.20129e-07
history = [ 0.353645 0.0633671 0.0631533 0.0210318 0.0209579 0.00767153 0.00764574 0.00290698 0.00289398 0.00109641 0.00109227 0.000392763 0.000390975 0.000137503 0.000137032 4.70734e-05 4.69283e-05 1.64386e-05 1.63875e-05 5.47276e-06 5.45639e-06 1.7953e-06 1.78959e-06 5.85076e-07 ] (24)
|M*(A*X-B)|/|B| = 6.0409e-07


Go back to Summary of /test programs.

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