GetFEM  5.4.2
gmm_precond_ilut.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 // This file is a modified version of ilut.h from ITL.
33 // See http://osl.iu.edu/research/itl/
34 // Following the corresponding Copyright notice.
35 //===========================================================================
36 //
37 // Copyright (c) 1998-2020, University of Notre Dame. All rights reserved.
38 // Redistribution and use in source and binary forms, with or without
39 // modification, are permitted provided that the following conditions are met:
40 //
41 // * Redistributions of source code must retain the above copyright
42 // notice, this list of conditions and the following disclaimer.
43 // * Redistributions in binary form must reproduce the above copyright
44 // notice, this list of conditions and the following disclaimer in the
45 // documentation and/or other materials provided with the distribution.
46 // * Neither the name of the University of Notre Dame nor the
47 // names of its contributors may be used to endorse or promote products
48 // derived from this software without specific prior written permission.
49 //
50 // THIS SOFTWARE IS PROVIDED BY THE TRUSTEES OF INDIANA UNIVERSITY AND
51 // CONTRIBUTORS ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
52 // BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
53 // FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE TRUSTEES
54 // OF INDIANA UNIVERSITY AND CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
55 // INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
56 // NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
57 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
58 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
59 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
60 // THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
61 //
62 //===========================================================================
63 
64 #ifndef GMM_PRECOND_ILUT_H
65 #define GMM_PRECOND_ILUT_H
66 
67 /**@file gmm_precond_ilut.h
68  @author Andrew Lumsdaine <[email protected]>, Lie-Quan Lee <[email protected]>
69  @date June 5, 2003.
70  @brief ILUT: Incomplete LU with threshold and K fill-in Preconditioner.
71 */
72 
73 /*
74  Performane comparing for SSOR, ILU and ILUT based on sherman 5 matrix
75  in Harwell-Boeing collection on Sun Ultra 30 UPA/PCI (UltraSPARC-II 296MHz)
76  Preconditioner & Factorization time & Number of Iteration \\ \hline
77  SSOR & 0.010577 & 41 \\
78  ILU & 0.019336 & 32 \\
79  ILUT with 0 fill-in and threshold of 1.0e-6 & 0.343612 & 23 \\
80  ILUT with 5 fill-in and threshold of 1.0e-6 & 0.343612 & 18 \\ \hline
81 */
82 
83 #include "gmm_precond.h"
84 
85 namespace gmm {
86 
87  template<typename T> struct elt_rsvector_value_less_ {
88  inline bool operator()(const elt_rsvector_<T>& a,
89  const elt_rsvector_<T>& b) const
90  { return (gmm::abs(a.e) > gmm::abs(b.e)); }
91  };
92 
93  /** Incomplete LU with threshold and K fill-in Preconditioner.
94 
95  The algorithm of ILUT(A, 0, 1.0e-6) is slower than ILU(A). If No
96  fill-in is arrowed, you can use ILU instead of ILUT.
97 
98  Notes: The idea under a concrete Preconditioner such as ilut is to
99  create a Preconditioner object to use in iterative methods.
100  */
101  template <typename Matrix>
102  class ilut_precond {
103  public :
104  typedef typename linalg_traits<Matrix>::value_type value_type;
105  typedef wsvector<value_type> _wsvector;
106  typedef rsvector<value_type> _rsvector;
107  typedef row_matrix<_rsvector> LU_Matrix;
108 
109  bool invert;
110  LU_Matrix L, U;
111 
112  protected:
113  size_type K;
114  double eps;
115 
116  template<typename M> void do_ilut(const M&, row_major);
117  void do_ilut(const Matrix&, col_major);
118 
119  public:
120  void build_with(const Matrix& A, int k_ = -1, double eps_ = double(-1)) {
121  if (k_ >= 0) K = k_;
122  if (eps_ >= double(0)) eps = eps_;
123  invert = false;
124  gmm::resize(L, mat_nrows(A), mat_ncols(A));
125  gmm::resize(U, mat_nrows(A), mat_ncols(A));
126  do_ilut(A, typename principal_orientation_type<typename
127  linalg_traits<Matrix>::sub_orientation>::potype());
128  }
129  ilut_precond(const Matrix& A, int k_, double eps_)
130  : L(mat_nrows(A), mat_ncols(A)), U(mat_nrows(A), mat_ncols(A)),
131  K(k_), eps(eps_) { build_with(A); }
132  ilut_precond(size_type k_, double eps_) : K(k_), eps(eps_) {}
133  ilut_precond(void) { K = 10; eps = 1E-7; }
134  size_type memsize() const {
135  return sizeof(*this) + (nnz(U)+nnz(L))*sizeof(value_type);
136  }
137  };
138 
139  template<typename Matrix> template<typename M>
140  void ilut_precond<Matrix>::do_ilut(const M& A, row_major) {
141  typedef value_type T;
142  typedef typename number_traits<T>::magnitude_type R;
143 
144  size_type n = mat_nrows(A);
145  if (n == 0) return;
146  std::vector<T> indiag(n);
147  _wsvector w(mat_ncols(A));
148  _rsvector ww(mat_ncols(A)), wL(mat_ncols(A)), wU(mat_ncols(A));
149  T tmp;
150  gmm::clear(U); gmm::clear(L);
151  R prec = default_tol(R());
152  R max_pivot = gmm::abs(A(0,0)) * prec;
153 
154  for (size_type i = 0; i < n; ++i) {
155  gmm::copy(mat_const_row(A, i), w);
156  double norm_row = gmm::vect_norm2(w);
157 
158  typename _wsvector::iterator wkold = w.end();
159  for (typename _wsvector::iterator wk = w.begin();
160  wk != w.end() && wk->first < i; ) {
161  size_type k = wk->first;
162  tmp = (wk->second) * indiag[k];
163  if (gmm::abs(tmp) < eps * norm_row) w.erase(k);
164  else { wk->second += tmp; gmm::add(scaled(mat_row(U, k), -tmp), w); }
165  if (wkold == w.end()) wk = w.begin(); else { wk = wkold; ++wk; }
166  if (wk != w.end() && wk->first == k)
167  { if (wkold == w.end()) wkold = w.begin(); else ++wkold; ++wk; }
168  }
169  tmp = w[i];
170 
171  if (gmm::abs(tmp) <= max_pivot) {
172  GMM_WARNING2("pivot " << i << " too small. try with ilutp ?");
173  w[i] = tmp = T(1);
174  }
175 
176  max_pivot = std::max(max_pivot, std::min(gmm::abs(tmp) * prec, R(1)));
177  indiag[i] = T(1) / tmp;
178  gmm::clean(w, eps * norm_row);
179  gmm::copy(w, ww);
180  std::sort(ww.begin(), ww.end(), elt_rsvector_value_less_<T>());
181  typename _rsvector::const_iterator wit = ww.begin(), wite = ww.end();
182 
183  size_type nnl = 0, nnu = 0;
184  wL.base_resize(K); wU.base_resize(K+1);
185  typename _rsvector::iterator witL = wL.begin(), witU = wU.begin();
186  for (; wit != wite; ++wit)
187  if (wit->c < i) { if (nnl < K) { *witL++ = *wit; ++nnl; } }
188  else { if (nnu < K || wit->c == i) { *witU++ = *wit; ++nnu; } }
189  wL.base_resize(nnl); wU.base_resize(nnu);
190  std::sort(wL.begin(), wL.end());
191  std::sort(wU.begin(), wU.end());
192  gmm::copy(wL, L.row(i));
193  gmm::copy(wU, U.row(i));
194  }
195 
196  }
197 
198  template<typename Matrix>
199  void ilut_precond<Matrix>::do_ilut(const Matrix& A, col_major) {
200  do_ilut(gmm::transposed(A), row_major());
201  invert = true;
202  }
203 
204  template <typename Matrix, typename V1, typename V2> inline
205  void mult(const ilut_precond<Matrix>& P, const V1 &v1, V2 &v2) {
206  gmm::copy(v1, v2);
207  if (P.invert) {
208  gmm::lower_tri_solve(gmm::transposed(P.U), v2, false);
209  gmm::upper_tri_solve(gmm::transposed(P.L), v2, true);
210  }
211  else {
212  gmm::lower_tri_solve(P.L, v2, true);
213  gmm::upper_tri_solve(P.U, v2, false);
214  }
215  }
216 
217  template <typename Matrix, typename V1, typename V2> inline
218  void transposed_mult(const ilut_precond<Matrix>& P,const V1 &v1,V2 &v2) {
219  gmm::copy(v1, v2);
220  if (P.invert) {
221  gmm::lower_tri_solve(P.L, v2, true);
222  gmm::upper_tri_solve(P.U, v2, false);
223  }
224  else {
225  gmm::lower_tri_solve(gmm::transposed(P.U), v2, false);
226  gmm::upper_tri_solve(gmm::transposed(P.L), v2, true);
227  }
228  }
229 
230  template <typename Matrix, typename V1, typename V2> inline
231  void left_mult(const ilut_precond<Matrix>& P, const V1 &v1, V2 &v2) {
232  copy(v1, v2);
233  if (P.invert) gmm::lower_tri_solve(gmm::transposed(P.U), v2, false);
234  else gmm::lower_tri_solve(P.L, v2, true);
235  }
236 
237  template <typename Matrix, typename V1, typename V2> inline
238  void right_mult(const ilut_precond<Matrix>& P, const V1 &v1, V2 &v2) {
239  copy(v1, v2);
240  if (P.invert) gmm::upper_tri_solve(gmm::transposed(P.L), v2, true);
241  else gmm::upper_tri_solve(P.U, v2, false);
242  }
243 
244  template <typename Matrix, typename V1, typename V2> inline
245  void transposed_left_mult(const ilut_precond<Matrix>& P, const V1 &v1,
246  V2 &v2) {
247  copy(v1, v2);
248  if (P.invert) gmm::upper_tri_solve(P.U, v2, false);
249  else gmm::upper_tri_solve(gmm::transposed(P.L), v2, true);
250  }
251 
252  template <typename Matrix, typename V1, typename V2> inline
253  void transposed_right_mult(const ilut_precond<Matrix>& P, const V1 &v1,
254  V2 &v2) {
255  copy(v1, v2);
256  if (P.invert) gmm::lower_tri_solve(P.L, v2, true);
257  else gmm::lower_tri_solve(gmm::transposed(P.U), v2, false);
258  }
259 
260 }
261 
262 #endif
263 
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_norm2
number_traits< typename linalg_traits< V >::value_type >::magnitude_type vect_norm2(const V &v)
Euclidean norm of a vector.
Definition: gmm_blas.h:557
gmm::rsvector
sparse vector built upon std::vector.
Definition: gmm_def.h:488
gmm::nnz
size_type nnz(const L &l)
count the number of non-zero entries of a vector or matrix.
Definition: gmm_blas.h:68
gmm_precond.h
gmm preconditioners.
gmm::copy
void copy(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:977
gmm::ilut_precond
Incomplete LU with threshold and K fill-in Preconditioner.
Definition: gmm_precond_ilut.h:102
gmm::wsvector
sparse vector built upon std::map.
Definition: gmm_def.h:487