GetFEM  5.4.2
gmm_precond_mr_approx_inverse.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2002-2020 Yves Renard
5 
6  This file is a part of GetFEM
7 
8  GetFEM is free software; you can redistribute it and/or modify it
9  under the terms of the GNU Lesser General Public License as published
10  by the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version along with the GCC Runtime Library
12  Exception either version 3.1 or (at your option) any later version.
13  This program is distributed in the hope that it will be useful, but
14  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  License and GCC Runtime Library Exception for more details.
17  You should have received a copy of the GNU Lesser General Public License
18  along with this program; if not, write to the Free Software Foundation,
19  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
20 
21  As a special exception, you may use this file as it is a part of a free
22  software library without restriction. Specifically, if other files
23  instantiate templates or use macros or inline functions from this file,
24  or you compile this file and link it with other files to produce an
25  executable, this file does not by itself cause the resulting executable
26  to be covered by the GNU Lesser General Public License. This exception
27  does not however invalidate any other reasons why the executable file
28  might be covered by the GNU Lesser General Public License.
29 
30 ===========================================================================*/
31 
32 
33 // This file is a modified version of approximate_inverse.h from ITL.
34 // See http://osl.iu.edu/research/itl/
35 // Following the corresponding Copyright notice.
36 //===========================================================================
37 //
38 // Copyright (c) 1998-2020, University of Notre Dame. All rights reserved.
39 // Redistribution and use in source and binary forms, with or without
40 // modification, are permitted provided that the following conditions are met:
41 //
42 // * Redistributions of source code must retain the above copyright
43 // notice, this list of conditions and the following disclaimer.
44 // * Redistributions in binary form must reproduce the above copyright
45 // notice, this list of conditions and the following disclaimer in the
46 // documentation and/or other materials provided with the distribution.
47 // * Neither the name of the University of Notre Dame nor the
48 // names of its contributors may be used to endorse or promote products
49 // derived from this software without specific prior written permission.
50 //
51 // THIS SOFTWARE IS PROVIDED BY THE TRUSTEES OF INDIANA UNIVERSITY AND
52 // CONTRIBUTORS ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
53 // BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
54 // FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE TRUSTEES
55 // OF INDIANA UNIVERSITY AND CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
56 // INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
57 // NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
58 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
59 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
60 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
61 // THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
62 //
63 //===========================================================================
64 
65 /**@file gmm_precond_mr_approx_inverse.h
66  @author Andrew Lumsdaine <[email protected]>
67  @author Lie-Quan Lee <[email protected]>
68  @author Yves Renard <[email protected]>
69  @date June 5, 2003.
70  @brief Approximate inverse via MR iteration.
71 */
72 
73 #ifndef GMM_PRECOND_MR_APPROX_INVERSE_H
74 #define GMM_PRECOND_MR_APPROX_INVERSE_H
75 
76 
77 #include "gmm_precond.h"
78 
79 namespace gmm {
80 
81  /** Approximate inverse via MR iteration (see P301 of Saad book).
82  */
83  template <typename Matrix>
84  struct mr_approx_inverse_precond {
85 
86  typedef typename linalg_traits<Matrix>::value_type value_type;
87  typedef typename number_traits<value_type>::magnitude_type magnitude_type;
88  typedef typename principal_orientation_type<typename
89  linalg_traits<Matrix>::sub_orientation>::potype sub_orientation;
90  typedef wsvector<value_type> VVector;
91  typedef col_matrix<VVector> MMatrix;
92 
93  MMatrix M;
94  size_type nb_it;
95  magnitude_type threshold;
96 
97  void build_with(const Matrix& A);
98  mr_approx_inverse_precond(const Matrix& A, size_type nb_it_,
99  magnitude_type threshold_)
100  : M(mat_nrows(A), mat_ncols(A))
101  { threshold = threshold_; nb_it = nb_it_; build_with(A); }
102  mr_approx_inverse_precond(void)
103  { threshold = magnitude_type(1E-7); nb_it = 5; }
104  mr_approx_inverse_precond(size_type nb_it_, magnitude_type threshold_)
105  { threshold = threshold_; nb_it = nb_it_; }
106  const MMatrix &approx_inverse(void) const { return M; }
107  };
108 
109  template <typename Matrix, typename V1, typename V2> inline
110  void mult(const mr_approx_inverse_precond<Matrix>& P, const V1 &v1, V2 &v2)
111  { mult(P.M, v1, v2); }
112 
113  template <typename Matrix, typename V1, typename V2> inline
114  void transposed_mult(const mr_approx_inverse_precond<Matrix>& P,
115  const V1 &v1,V2 &v2)
116  { mult(gmm::conjugated(P.M), v1, v2); }
117 
118  template <typename Matrix>
119  void mr_approx_inverse_precond<Matrix>::build_with(const Matrix& A) {
120  gmm::resize(M, mat_nrows(A), mat_ncols(A));
121  typedef value_type T;
122  typedef magnitude_type R;
123  VVector m(mat_ncols(A)),r(mat_ncols(A)),ei(mat_ncols(A)),Ar(mat_ncols(A));
124  T alpha = mat_trace(A)/ mat_euclidean_norm_sqr(A);
125  if (alpha == T(0)) alpha = T(1);
126 
127  for (size_type i = 0; i < mat_nrows(A); ++i) {
128  gmm::clear(m); gmm::clear(ei);
129  m[i] = alpha;
130  ei[i] = T(1);
131 
132  for (size_type j = 0; j < nb_it; ++j) {
133  gmm::mult(A, gmm::scaled(m, T(-1)), r);
134  gmm::add(ei, r);
135  gmm::mult(A, r, Ar);
136  T nAr = vect_sp(Ar,Ar);
137  if (gmm::abs(nAr) > R(0)) {
138  gmm::add(gmm::scaled(r, gmm::safe_divide(vect_sp(r, Ar), vect_sp(Ar, Ar))), m);
139  gmm::clean(m, threshold * gmm::vect_norm2(m));
140  } else gmm::clear(m);
141  }
142  if (gmm::vect_norm2(m) == R(0)) m[i] = alpha;
143  gmm::copy(m, M.col(i));
144  }
145  }
146 }
147 
148 #endif
149 
gmm::resize
void resize(M &v, size_type m, size_type n)
*‍/
Definition: gmm_blas.h:231
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:59
gmm::vect_sp
strongest_value_type< V1, V2 >::value_type vect_sp(const V1 &v1, const V2 &v2)
*‍/
Definition: gmm_blas.h:263
gmm::mat_euclidean_norm_sqr
number_traits< typename linalg_traits< M >::value_type >::magnitude_type mat_euclidean_norm_sqr(const M &m)
*‍/
Definition: gmm_blas.h:626
bgeot::alpha
size_type alpha(short_type n, short_type d)
Return the value of which is the number of monomials of a polynomial of variables and degree .
Definition: bgeot_poly.cc:47
gmm_precond.h
gmm preconditioners.
gmm::mat_trace
linalg_traits< M >::value_type mat_trace(const M &m)
Trace of a matrix.
Definition: gmm_blas.h:528
gmm::mr_approx_inverse_precond
Approximate inverse via MR iteration (see P301 of Saad book).
Definition: gmm_precond_mr_approx_inverse.h:84
gmm::wsvector
sparse vector built upon std::map.
Definition: gmm_def.h:487