G+Smo  25.01.0
Geometry + Simulation Modules
 
Loading...
Searching...
No Matches
gsMinimalResidual.hpp
Go to the documentation of this file.
1
14namespace gismo
15{
16
17template<class T>
19{
20 Base::initIteration(rhs,x);
21 //if (Base::initIteration(rhs,x)) return true; // z will not be initialized!
22
23 index_t n = m_mat->cols();
24 index_t m = 1; // = rhs.cols();
25
26 vPrev.setZero(n,m); vNew.setZero(n,m);
27 wPrev.setZero(n,m); w.setZero(n,m); wNew.setZero(n,m);
28 if (!m_inexact_residual)
29 { AwPrev.setZero(n,m); Aw.setZero(n,m); AwNew.setZero(n,m); }
30
31 m_mat->apply(x,negResidual);
32 negResidual -= rhs;
33
34 m_error = negResidual.norm() / m_rhs_norm;
35 if (m_error < m_tol)
36 return true;
37
38 v = -negResidual;
39 m_precond->apply(v, z);
40
41 gammaPrev = 1;
42 T ip = z.col(0).dot(v.col(0));
43 GISMO_ASSERT(ip >= T(0), "gsMinimalResidual::initIteration(...), preconditioner not positive semi-definite");
44 gamma = math::sqrt(ip);
45 gammaNew = 1;
46 eta = gamma;
47 sPrev = 0; s = 0; sNew = 0;
48 cPrev = 1; c = 1; cNew = 1;
49
50 return false;
51}
52
53template<class T>
55{
56 z /= gamma;
57 m_mat->apply(z,Az);
58
59 T delta = z.col(0).dot(Az.col(0));
60 vNew = Az - (delta/gamma)*v - (gamma/gammaPrev)*vPrev;
61 m_precond->apply(vNew, zNew);
62 T ip = zNew.col(0).dot(vNew.col(0));
63 GISMO_ASSERT(ip >= T(0), "gsMinimalResidual::step(...), preconditioner not positive semi-definite");
64 gammaNew = math::sqrt(ip);
65 const T a0 = c*delta - cPrev*s*gamma;
66 const T a1 = math::sqrt(a0*a0 + gammaNew*gammaNew);
67 const T a2 = s*delta + cPrev*c*gamma;
68 const T a3 = sPrev*gamma;
69 cNew = a0/a1;
70 sNew = gammaNew/a1;
71 wNew = (z - a3*wPrev - a2*w)/a1;
72 if (!m_inexact_residual)
73 AwNew = (Az - a3*AwPrev - a2*Aw)/a1;
74 x += cNew*eta*wNew;
75 if (!m_inexact_residual)
76 negResidual += cNew*eta*AwNew;
77
78 if (m_inexact_residual)
79 m_error *= math::abs(sNew); // see https://eigen.tuxfamily.org/dox-devel/unsupported/MINRES_8h_source.html
80 else
81 m_error = negResidual.norm() / m_rhs_norm;
82
83 eta = -sNew*eta;
84
85 // Test for convergence
86 if (m_error < m_tol)
87 return true;
88
89 //Update variables
90 vPrev.swap(v); v.swap(vNew); // for us the same as: vPrev = v; v = vNew;
91 wPrev.swap(w); w.swap(wNew); // for us the same as: wPrev = w; w = wNew;
92 if (!m_inexact_residual)
93 { AwPrev.swap(Aw); Aw.swap(AwNew); } // for us the same as: AwPrev = Aw; Aw = AwNew;
94 z.swap(zNew); // for us the same as: z = zNew;
95 gammaPrev = gamma; gamma = gammaNew;
96 sPrev = s; s = sNew;
97 cPrev = c; c = cNew;
98 return false;
99}
100
101template<class T>
103{
104 // cleanup temporaries
105 negResidual.clear();
106 vPrev.clear(); v.clear(); vNew.clear();
107 wPrev.clear(); w.clear(); wNew.clear();
108 AwPrev.clear(); Aw.clear(); AwNew.clear();
109 zNew.clear(); z.clear(); Az.clear();
110}
111
112
113
114}
void finalizeIteration(VectorType &x)
Some post-processing might be required.
Definition gsMinimalResidual.hpp:102
bool step(VectorType &x)
Perform one step, requires initIteration.
Definition gsMinimalResidual.hpp:54
bool initIteration(const VectorType &rhs, VectorType &x)
Init the iteration.
Definition gsMinimalResidual.hpp:18
#define index_t
Definition gsConfig.h:32
#define GISMO_ASSERT(cond, message)
Definition gsDebug.h:89
The G+Smo namespace, containing all definitions for the library.