GetFEM  5.4.2
gmm_precond_ilutp.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2004-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 /**@file gmm_precond_ilutp.h
33  @author Yves Renard <Yves.Renard@insa-lyon.fr>
34  @date October 14, 2004.
35  @brief ILUTP: Incomplete LU with threshold and K fill-in Preconditioner and
36  column pivoting.
37 
38 
39 */
40 #ifndef GMM_PRECOND_ILUTP_H
41 #define GMM_PRECOND_ILUTP_H
42 
43 #include "gmm_precond_ilut.h"
44 
45 namespace gmm {
46 
47  /**
48  ILUTP: Incomplete LU with threshold and K fill-in Preconditioner and
49  column pivoting.
50 
51  See Yousef Saad, Iterative Methods for
52  sparse linear systems, PWS Publishing Company, section 10.4.4
53 
54  TODO : store the permutation by cycles to avoid the temporary vector
55  */
56  template <typename Matrix>
57  class ilutp_precond {
58  public :
59  typedef typename linalg_traits<Matrix>::value_type value_type;
60  typedef wsvector<value_type> _wsvector;
61  typedef rsvector<value_type> _rsvector;
62  typedef row_matrix<_rsvector> LU_Matrix;
63  typedef col_matrix<_wsvector> CLU_Matrix;
64 
65  bool invert;
66  LU_Matrix L, U;
67  gmm::unsorted_sub_index indperm;
68  gmm::unsorted_sub_index indperminv;
69  mutable std::vector<value_type> temporary;
70 
71  protected:
72  size_type K;
73  double eps;
74 
75  template<typename M> void do_ilutp(const M&, row_major);
76  void do_ilutp(const Matrix&, col_major);
77 
78  public:
79  void build_with(const Matrix& A, int k_ = -1, double eps_ = double(-1)) {
80  if (k_ >= 0) K = k_;
81  if (eps_ >= double(0)) eps = eps_;
82  invert = false;
83  gmm::resize(L, mat_nrows(A), mat_ncols(A));
84  gmm::resize(U, mat_nrows(A), mat_ncols(A));
85  do_ilutp(A, typename principal_orientation_type<typename
86  linalg_traits<Matrix>::sub_orientation>::potype());
87  }
88  ilutp_precond(const Matrix& A, size_type k_, double eps_)
89  : L(mat_nrows(A), mat_ncols(A)), U(mat_nrows(A), mat_ncols(A)),
90  K(k_), eps(eps_) { build_with(A); }
91  ilutp_precond(int k_, double eps_) : K(k_), eps(eps_) {}
92  ilutp_precond(void) { K = 10; eps = 1E-7; }
93  size_type memsize() const {
94  return sizeof(*this) + (nnz(U)+nnz(L))*sizeof(value_type);
95  }
96  };
97 
98 
99  template<typename Matrix> template<typename M>
100  void ilutp_precond<Matrix>::do_ilutp(const M& A, row_major) {
101  typedef value_type T;
102  typedef typename number_traits<T>::magnitude_type R;
103 
104  size_type n = mat_nrows(A);
105  CLU_Matrix CU(n,n);
106  if (n == 0) return;
107  std::vector<T> indiag(n);
108  temporary.resize(n);
109  std::vector<size_type> ipvt(n), ipvtinv(n);
110  for (size_type i = 0; i < n; ++i) ipvt[i] = ipvtinv[i] = i;
111  indperm = unsorted_sub_index(ipvt);
112  indperminv = unsorted_sub_index(ipvtinv);
113  _wsvector w(mat_ncols(A));
114  _rsvector ww(mat_ncols(A));
115 
116  T tmp = T(0);
117  gmm::clear(L); gmm::clear(U);
118  R prec = default_tol(R());
119  R max_pivot = gmm::abs(A(0,0)) * prec;
120 
121  for (size_type i = 0; i < n; ++i) {
122 
123  copy(sub_vector(mat_const_row(A, i), indperm), w);
124  double norm_row = gmm::vect_norm2(mat_const_row(A, i));
125 
126  typename _wsvector::iterator wkold = w.end();
127  for (typename _wsvector::iterator wk = w.begin();
128  wk != w.end() && wk->first < i; ) {
129  size_type k = wk->first;
130  tmp = (wk->second) * indiag[k];
131  if (gmm::abs(tmp) < eps * norm_row) w.erase(k);
132  else { wk->second += tmp; gmm::add(scaled(mat_row(U, k), -tmp), w); }
133  if (wkold == w.end()) wk = w.begin(); else { wk = wkold; ++wk; }
134  if (wk != w.end() && wk->first == k)
135  { if (wkold == w.end()) wkold = w.begin(); else ++wkold; ++wk; }
136  }
137 
138  gmm::clean(w, eps * norm_row);
139  gmm::copy(w, ww);
140 
141  std::sort(ww.begin(), ww.end(), elt_rsvector_value_less_<T>());
142  typename _rsvector::const_iterator wit = ww.begin(), wite = ww.end();
143  size_type ip = size_type(-1);
144 
145  for (; wit != wite; ++wit)
146  if (wit->c >= i) { ip = wit->c; tmp = wit->e; break; }
147  if (ip == size_type(-1) || gmm::abs(tmp) <= max_pivot)
148  { GMM_WARNING2("pivot " << i << " too small"); ip=i; ww[i]=tmp=T(1); }
149  max_pivot = std::max(max_pivot, std::min(gmm::abs(tmp) * prec, R(1)));
150  indiag[i] = T(1) / tmp;
151  wit = ww.begin();
152 
153  size_type nnl = 0, nnu = 0;
154  L[i].base_resize(K); U[i].base_resize(K+1);
155  typename _rsvector::iterator witL = L[i].begin(), witU = U[i].begin();
156  for (; wit != wite; ++wit) {
157  if (wit->c < i) { if (nnl < K) { *witL++ = *wit; ++nnl; } }
158  else if (nnu < K || wit->c == i)
159  { CU(i, wit->c) = wit->e; *witU++ = *wit; ++nnu; }
160  }
161  L[i].base_resize(nnl); U[i].base_resize(nnu);
162  std::sort(L[i].begin(), L[i].end());
163  std::sort(U[i].begin(), U[i].end());
164 
165  if (ip != i) {
166  typename _wsvector::const_iterator iti = CU.col(i).begin();
167  typename _wsvector::const_iterator itie = CU.col(i).end();
168  typename _wsvector::const_iterator itp = CU.col(ip).begin();
169  typename _wsvector::const_iterator itpe = CU.col(ip).end();
170 
171  while (iti != itie && itp != itpe) {
172  if (iti->first < itp->first)
173  { U.row(iti->first).swap_indices(i, ip); ++iti; }
174  else if (iti->first > itp->first)
175  { U.row(itp->first).swap_indices(i,ip);++itp; }
176  else
177  { U.row(iti->first).swap_indices(i, ip); ++iti; ++itp; }
178  }
179 
180  for( ; iti != itie; ++iti) U.row(iti->first).swap_indices(i, ip);
181  for( ; itp != itpe; ++itp) U.row(itp->first).swap_indices(i, ip);
182 
183  CU.swap_col(i, ip);
184 
185  indperm.swap(i, ip);
186  indperminv.swap(ipvt[i], ipvt[ip]);
187  std::swap(ipvtinv[ipvt[i]], ipvtinv[ipvt[ip]]);
188  std::swap(ipvt[i], ipvt[ip]);
189  }
190  }
191  }
192 
193  template<typename Matrix>
194  void ilutp_precond<Matrix>::do_ilutp(const Matrix& A, col_major) {
195  do_ilutp(gmm::transposed(A), row_major());
196  invert = true;
197  }
198 
199  template <typename Matrix, typename V1, typename V2> inline
200  void mult(const ilutp_precond<Matrix>& P, const V1 &v1, V2 &v2) {
201  if (P.invert) {
202  gmm::copy(gmm::sub_vector(v1, P.indperm), v2);
203  gmm::lower_tri_solve(gmm::transposed(P.U), v2, false);
204  gmm::upper_tri_solve(gmm::transposed(P.L), v2, true);
205  }
206  else {
207  gmm::copy(v1, P.temporary);
208  gmm::lower_tri_solve(P.L, P.temporary, true);
209  gmm::upper_tri_solve(P.U, P.temporary, false);
210  gmm::copy(gmm::sub_vector(P.temporary, P.indperminv), v2);
211  }
212  }
213 
214  template <typename Matrix, typename V1, typename V2> inline
215  void transposed_mult(const ilutp_precond<Matrix>& P,const V1 &v1,V2 &v2) {
216  if (P.invert) {
217  gmm::copy(v1, P.temporary);
218  gmm::lower_tri_solve(P.L, P.temporary, true);
219  gmm::upper_tri_solve(P.U, P.temporary, false);
220  gmm::copy(gmm::sub_vector(P.temporary, P.indperminv), v2);
221  }
222  else {
223  gmm::copy(gmm::sub_vector(v1, P.indperm), v2);
224  gmm::lower_tri_solve(gmm::transposed(P.U), v2, false);
225  gmm::upper_tri_solve(gmm::transposed(P.L), v2, true);
226  }
227  }
228 
229  template <typename Matrix, typename V1, typename V2> inline
230  void left_mult(const ilutp_precond<Matrix>& P, const V1 &v1, V2 &v2) {
231  if (P.invert) {
232  gmm::copy(gmm::sub_vector(v1, P.indperm), v2);
233  gmm::lower_tri_solve(gmm::transposed(P.U), v2, false);
234  }
235  else {
236  copy(v1, v2);
237  gmm::lower_tri_solve(P.L, v2, true);
238  }
239  }
240 
241  template <typename Matrix, typename V1, typename V2> inline
242  void right_mult(const ilutp_precond<Matrix>& P, const V1 &v1, V2 &v2) {
243  if (P.invert) {
244  copy(v1, v2);
245  gmm::upper_tri_solve(gmm::transposed(P.L), v2, true);
246  }
247  else {
248  copy(v1, P.temporary);
249  gmm::upper_tri_solve(P.U, P.temporary, false);
250  gmm::copy(gmm::sub_vector(P.temporary, P.indperminv), v2);
251  }
252  }
253 
254  template <typename Matrix, typename V1, typename V2> inline
255  void transposed_left_mult(const ilutp_precond<Matrix>& P, const V1 &v1,
256  V2 &v2) {
257  if (P.invert) {
258  copy(v1, P.temporary);
259  gmm::upper_tri_solve(P.U, P.temporary, false);
260  gmm::copy(gmm::sub_vector(P.temporary, P.indperminv), v2);
261  }
262  else {
263  copy(v1, v2);
264  gmm::upper_tri_solve(gmm::transposed(P.L), v2, true);
265  }
266  }
267 
268  template <typename Matrix, typename V1, typename V2> inline
269  void transposed_right_mult(const ilutp_precond<Matrix>& P, const V1 &v1,
270  V2 &v2) {
271  if (P.invert) {
272  copy(v1, v2);
273  gmm::lower_tri_solve(P.L, v2, true);
274  }
275  else {
276  gmm::copy(gmm::sub_vector(v1, P.indperm), v2);
277  gmm::lower_tri_solve(gmm::transposed(P.U), v2, false);
278  }
279  }
280 
281 }
282 
283 #endif
284 
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_precond_ilut.h
ILUT: Incomplete LU with threshold and K fill-in Preconditioner.
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::ilutp_precond
ILUTP: Incomplete LU with threshold and K fill-in Preconditioner and column pivoting.
Definition: gmm_precond_ilutp.h:57
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::copy
void copy(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:977
gmm::wsvector
sparse vector built upon std::map.
Definition: gmm_def.h:487