GetFEM++  5.3
gmm_precond_ilu.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2002-2017 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 ilu.h from ITL.
33 // See http://osl.iu.edu/research/itl/
34 // Following the corresponding Copyright notice.
35 //===========================================================================
36 //
37 // Copyright (c) 1998-2001, 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 /**@file gmm_precond_ilu.h
65  @author Andrew Lumsdaine <lums@osl.iu.edu>
66  @author Lie-Quan Lee <llee@osl.iu.edu>
67  @author Yves Renard <yves.renard@insa-lyon.fr>
68  @date June 5, 2003.
69  @brief Incomplete LU without fill-in Preconditioner.
70 */
71 
72 #ifndef GMM_PRECOND_ILU_H
73 #define GMM_PRECOND_ILU_H
74 
75 //
76 // Notes: The idea under a concrete Preconditioner such
77 // as Incomplete LU is to create a Preconditioner
78 // object to use in iterative methods.
79 //
80 
81 #include "gmm_precond.h"
82 
83 namespace gmm {
84  /** Incomplete LU without fill-in Preconditioner. */
85  template <typename Matrix>
86  class ilu_precond {
87 
88  public :
89  typedef typename linalg_traits<Matrix>::value_type value_type;
90  typedef csr_matrix_ref<value_type *, size_type *, size_type *, 0> tm_type;
91 
92  tm_type U, L;
93  bool invert;
94  protected :
95  std::vector<value_type> L_val, U_val;
96  std::vector<size_type> L_ind, U_ind, L_ptr, U_ptr;
97 
98  template<typename M> void do_ilu(const M& A, row_major);
99  void do_ilu(const Matrix& A, col_major);
100 
101  public:
102 
103  size_type nrows(void) const { return mat_nrows(L); }
104  size_type ncols(void) const { return mat_ncols(U); }
105 
106  void build_with(const Matrix& A) {
107  invert = false;
108  L_ptr.resize(mat_nrows(A)+1);
109  U_ptr.resize(mat_nrows(A)+1);
110  do_ilu(A, typename principal_orientation_type<typename
111  linalg_traits<Matrix>::sub_orientation>::potype());
112  }
113  ilu_precond(const Matrix& A) { build_with(A); }
114  ilu_precond(void) {}
115  size_type memsize() const {
116  return sizeof(*this) +
117  (L_val.size()+U_val.size()) * sizeof(value_type) +
118  (L_ind.size()+L_ptr.size()) * sizeof(size_type) +
119  (U_ind.size()+U_ptr.size()) * sizeof(size_type);
120  }
121  };
122 
123  template <typename Matrix> template <typename M>
124  void ilu_precond<Matrix>::do_ilu(const M& A, row_major) {
125  typedef typename linalg_traits<Matrix>::storage_type store_type;
126  typedef value_type T;
127  typedef typename number_traits<T>::magnitude_type R;
128 
129  size_type L_loc = 0, U_loc = 0, n = mat_nrows(A), i, j, k;
130  if (n == 0) return;
131  L_ptr[0] = 0; U_ptr[0] = 0;
132  R prec = default_tol(R());
133  R max_pivot = gmm::abs(A(0,0)) * prec;
134 
135 
136  for (int count = 0; count < 2; ++count) {
137  if (count) {
138  L_val.resize(L_loc); L_ind.resize(L_loc);
139  U_val.resize(U_loc); U_ind.resize(U_loc);
140  }
141  L_loc = U_loc = 0;
142  for (i = 0; i < n; ++i) {
143  typedef typename linalg_traits<M>::const_sub_row_type row_type;
144  row_type row = mat_const_row(A, i);
145  typename linalg_traits<typename org_type<row_type>::t>::const_iterator
146  it = vect_const_begin(row), ite = vect_const_end(row);
147 
148  if (count) { U_val[U_loc] = T(0); U_ind[U_loc] = i; }
149  ++U_loc; // diagonal element
150 
151  for (k = 0; it != ite && k < 1000; ++it, ++k) {
152  // if a plain row is present, retains only the 1000 firsts
153  // nonzero elements. ---> a sort should be done.
154  j = index_of_it(it, k, store_type());
155  if (j < i) {
156  if (count) { L_val[L_loc] = *it; L_ind[L_loc] = j; }
157  L_loc++;
158  }
159  else if (i == j) {
160  if (count) U_val[U_loc-1] = *it;
161  }
162  else {
163  if (count) { U_val[U_loc] = *it; U_ind[U_loc] = j; }
164  U_loc++;
165  }
166  }
167  L_ptr[i+1] = L_loc; U_ptr[i+1] = U_loc;
168  }
169  }
170 
171  if (A(0,0) == T(0)) {
172  U_val[U_ptr[0]] = T(1);
173  GMM_WARNING2("pivot 0 is too small");
174  }
175 
176  size_type qn, pn, rn;
177  for (i = 1; i < n; i++) {
178 
179  pn = U_ptr[i];
180  if (gmm::abs(U_val[pn]) <= max_pivot) {
181  U_val[pn] = T(1);
182  GMM_WARNING2("pivot " << i << " is too small");
183  }
184  max_pivot = std::max(max_pivot,
185  std::min(gmm::abs(U_val[pn]) * prec, R(1)));
186 
187  for (j = L_ptr[i]; j < L_ptr[i+1]; j++) {
188  pn = U_ptr[L_ind[j]];
189 
190  T multiplier = (L_val[j] /= U_val[pn]);
191 
192  qn = j + 1;
193  rn = U_ptr[i];
194 
195  for (pn++; pn < U_ptr[L_ind[j]+1] && U_ind[pn] < i; pn++) {
196  while (qn < L_ptr[i+1] && L_ind[qn] < U_ind[pn])
197  qn++;
198  if (qn < L_ptr[i+1] && U_ind[pn] == L_ind[qn])
199  L_val[qn] -= multiplier * U_val[pn];
200  }
201  for (; pn < U_ptr[L_ind[j]+1]; pn++) {
202  while (rn < U_ptr[i+1] && U_ind[rn] < U_ind[pn])
203  rn++;
204  if (rn < U_ptr[i+1] && U_ind[pn] == U_ind[rn])
205  U_val[rn] -= multiplier * U_val[pn];
206  }
207  }
208  }
209 
210  L = tm_type(&(L_val[0]), &(L_ind[0]), &(L_ptr[0]), n, mat_ncols(A));
211  U = tm_type(&(U_val[0]), &(U_ind[0]), &(U_ptr[0]), n, mat_ncols(A));
212  }
213 
214  template <typename Matrix>
215  void ilu_precond<Matrix>::do_ilu(const Matrix& A, col_major) {
216  do_ilu(gmm::transposed(A), row_major());
217  invert = true;
218  }
219 
220  template <typename Matrix, typename V1, typename V2> inline
221  void mult(const ilu_precond<Matrix>& P, const V1 &v1, V2 &v2) {
222  gmm::copy(v1, v2);
223  if (P.invert) {
224  gmm::lower_tri_solve(gmm::transposed(P.U), v2, false);
225  gmm::upper_tri_solve(gmm::transposed(P.L), v2, true);
226  }
227  else {
228  gmm::lower_tri_solve(P.L, v2, true);
229  gmm::upper_tri_solve(P.U, v2, false);
230  }
231  }
232 
233  template <typename Matrix, typename V1, typename V2> inline
234  void transposed_mult(const ilu_precond<Matrix>& P,const V1 &v1,V2 &v2) {
235  gmm::copy(v1, v2);
236  if (P.invert) {
237  gmm::lower_tri_solve(P.L, v2, true);
238  gmm::upper_tri_solve(P.U, v2, false);
239  }
240  else {
241  gmm::lower_tri_solve(gmm::transposed(P.U), v2, false);
242  gmm::upper_tri_solve(gmm::transposed(P.L), v2, true);
243  }
244  }
245 
246  template <typename Matrix, typename V1, typename V2> inline
247  void left_mult(const ilu_precond<Matrix>& P, const V1 &v1, V2 &v2) {
248  copy(v1, v2);
249  if (P.invert) gmm::lower_tri_solve(gmm::transposed(P.U), v2, false);
250  else gmm::lower_tri_solve(P.L, v2, true);
251  }
252 
253  template <typename Matrix, typename V1, typename V2> inline
254  void right_mult(const ilu_precond<Matrix>& P, const V1 &v1, V2 &v2) {
255  copy(v1, v2);
256  if (P.invert) gmm::upper_tri_solve(gmm::transposed(P.L), v2, true);
257  else gmm::upper_tri_solve(P.U, v2, false);
258  }
259 
260  template <typename Matrix, typename V1, typename V2> inline
261  void transposed_left_mult(const ilu_precond<Matrix>& P, const V1 &v1,
262  V2 &v2) {
263  copy(v1, v2);
264  if (P.invert) gmm::upper_tri_solve(P.U, v2, false);
265  else gmm::upper_tri_solve(gmm::transposed(P.L), v2, true);
266  }
267 
268  template <typename Matrix, typename V1, typename V2> inline
269  void transposed_right_mult(const ilu_precond<Matrix>& P, const V1 &v1,
270  V2 &v2) {
271  copy(v1, v2);
272  if (P.invert) gmm::lower_tri_solve(P.L, v2, true);
273  else gmm::lower_tri_solve(gmm::transposed(P.U), v2, false);
274  }
275 
276 
277 }
278 
279 #endif
280 
Incomplete LU without fill-in Preconditioner.
void copy(const L1 &l1, L2 &l2)
*/
Definition: gmm_blas.h:977
gmm preconditioners.