1 /** @file gsMinimalResidual.h
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 #pragma once
14 
15 #include <gsSolver/gsIterativeSolver.h>
16 
17 namespace gismo
18 {
19 
20 /** @brief The minimal residual (MinRes) method.
21   *
22   * \ingroup Solver
23   */
24 template<class T=real_t>
25 class gsMinimalResidual : public gsIterativeSolver<T>
26 {
27 
28 public:
29     typedef gsIterativeSolver<T> Base;
30 
31     typedef gsMatrix<T>  VectorType;
32 
33     typedef typename Base::LinOpPtr LinOpPtr;
34 
35     typedef memory::shared_ptr<gsMinimalResidual> Ptr;
36     typedef memory::unique_ptr<gsMinimalResidual> uPtr;
37 
38     /// @brief Constructor using a matrix (operator) and optionally a preconditionner
39     ///
40     /// @param mat     The operator to be solved for, see gsIterativeSolver for details
41     /// @param precond The preconditioner, defaulted to the identity
42     template< typename OperatorType >
43     explicit gsMinimalResidual( const OperatorType& mat,
44                                 const LinOpPtr& precond = LinOpPtr())
Base(mat,precond)45         : Base(mat, precond), m_inexact_residual(false) {}
46 
47     /// @brief Make function using a matrix (operator) and optionally a preconditionner
48     ///
49     /// @param mat     The operator to be solved for, see gsIterativeSolver for details
50     /// @param precond The preconditioner, defaulted to the identity
51     template< typename OperatorType >
52     static uPtr make( const OperatorType& mat, const LinOpPtr& precond = LinOpPtr() )
53     { return uPtr( new gsMinimalResidual(mat, precond) ); }
54 
55     bool initIteration( const VectorType& rhs, VectorType& x );
56     void finalizeIteration( VectorType& x );
57 
58     bool step( VectorType& x );
59 
60     /// @brief Returns a list of default options
defaultOptions()61     static gsOptionList defaultOptions()
62     {
63         gsOptionList opt = Base::defaultOptions();
64         opt.addSwitch( "InexactResidual",
65                        "If true, the residual is estimated, not accurately computed.",
66                        false );
67         return opt;
68     }
69 
70     /// @brief Set the options based on a gsOptionList
setOptions(const gsOptionList & opt)71     gsMinimalResidual& setOptions(const gsOptionList & opt)
72     {
73         Base::setOptions(opt);
74         m_inexact_residual = opt.askSwitch("InexactResidual", m_inexact_residual);
75         return *this;
76     }
77 
78     /// @brief If true, the residual is estimated, not accurately computed.
setInexactResidual(bool flag)79     void setInexactResidual( bool flag )     { m_inexact_residual = flag; }
80 
81     /// Prints the object as a string.
print(std::ostream & os)82     std::ostream &print(std::ostream &os) const
83     {
84         os << "gsMinimalResidual\n";
85         return os;
86     }
87 
88 private:
89     using Base::m_mat;
90     using Base::m_precond;
91     using Base::m_max_iters;
92     using Base::m_tol;
93     using Base::m_num_iter;
94     using Base::m_rhs_norm;
95     using Base::m_error;
96 
97     gsMatrix<T> negResidual,
98                      vPrev, v, vNew,
99                      wPrev, w, wNew, AwPrev, Aw, AwNew,
100                      zNew, z, Az;
101 
102     T eta,
103            gammaPrev, gamma, gammaNew,
104            sPrev, s, sNew,
105            cPrev, c, cNew;
106 
107     bool m_inexact_residual;
108 };
109 
110 } // namespace gismo
111 
112 #ifndef GISMO_BUILD_LIB
113 #include GISMO_HPP_HEADER(gsMinimalResidual.hpp)
114 #endif
115