1 /** @file gsMinimalResidual.hpp
2
3 @brief Preconditioned iterative solver using the minimal residual method.
4
5 This file is part of the G+Smo library.
6
7 This Source Code Form is subject to the terms of the Mozilla Public
8 License, v. 2.0. If a copy of the MPL was not distributed with this
9 file, You can obtain one at http://mozilla.org/MPL/2.0/.
10
11 Author(s): J. Sogn, S. Takacs
12 */
13
14 namespace gismo
15 {
16
17 template<class T>
initIteration(const typename gsMinimalResidual<T>::VectorType & rhs,typename gsMinimalResidual<T>::VectorType & x)18 bool gsMinimalResidual<T>::initIteration( const typename gsMinimalResidual<T>::VectorType& rhs, typename gsMinimalResidual<T>::VectorType& x )
19 {
20 Base::initIteration(rhs,x);
21 //if (Base::initIteration(rhs,x)) return true; // z will not be initialized!
22
23 int n = m_mat->cols();
24 int 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; gamma = math::sqrt(z.col(0).dot(v.col(0))); gammaNew = 1;
42 eta = gamma;
43 sPrev = 0; s = 0; sNew = 0;
44 cPrev = 1; c = 1; cNew = 1;
45
46 return false;
47 }
48
49 template<class T>
step(typename gsMinimalResidual<T>::VectorType & x)50 bool gsMinimalResidual<T>::step( typename gsMinimalResidual<T>::VectorType& x )
51 {
52 z /= gamma;
53 m_mat->apply(z,Az);
54
55 T delta = z.col(0).dot(Az.col(0));
56 vNew = Az - (delta/gamma)*v - (gamma/gammaPrev)*vPrev;
57 m_precond->apply(vNew, zNew);
58 gammaNew = math::sqrt(zNew.col(0).dot(vNew.col(0)));
59 const T a0 = c*delta - cPrev*s*gamma;
60 const T a1 = math::sqrt(a0*a0 + gammaNew*gammaNew);
61 const T a2 = s*delta + cPrev*c*gamma;
62 const T a3 = sPrev*gamma;
63 cNew = a0/a1;
64 sNew = gammaNew/a1;
65 wNew = (z - a3*wPrev - a2*w)/a1;
66 if (!m_inexact_residual)
67 AwNew = (Az - a3*AwPrev - a2*Aw)/a1;
68 x += cNew*eta*wNew;
69 if (!m_inexact_residual)
70 negResidual += cNew*eta*AwNew;
71
72 if (m_inexact_residual)
73 m_error *= math::abs(sNew); // see https://eigen.tuxfamily.org/dox-devel/unsupported/MINRES_8h_source.html
74 else
75 m_error = negResidual.norm() / m_rhs_norm;
76
77 eta = -sNew*eta;
78
79 // Test for convergence
80 if (m_error < m_tol)
81 return true;
82
83 //Update variables
84 vPrev.swap(v); v.swap(vNew); // for us the same as: vPrev = v; v = vNew;
85 wPrev.swap(w); w.swap(wNew); // for us the same as: wPrev = w; w = wNew;
86 if (!m_inexact_residual)
87 { AwPrev.swap(Aw); Aw.swap(AwNew); } // for us the same as: AwPrev = Aw; Aw = AwNew;
88 z.swap(zNew); // for us the same as: z = zNew;
89 gammaPrev = gamma; gamma = gammaNew;
90 sPrev = s; s = sNew;
91 cPrev = c; c = cNew;
92 return false;
93 }
94
95 template<class T>
finalizeIteration(typename gsMinimalResidual<T>::VectorType &)96 void gsMinimalResidual<T>::finalizeIteration( typename gsMinimalResidual<T>::VectorType& )
97 {
98 // cleanup temporaries
99 negResidual.clear();
100 vPrev.clear(); v.clear(); vNew.clear();
101 wPrev.clear(); w.clear(); wNew.clear();
102 AwPrev.clear(); Aw.clear(); AwNew.clear();
103 zNew.clear(); z.clear(); Az.clear();
104 }
105
106
107
108 }
109