GetFEM  5.4.2
gmm_precond_ildltt.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2003-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_ildltt.h
33  @author Yves Renard <[email protected]>
34  @date June 30, 2003.
35  @brief incomplete LDL^t (cholesky) preconditioner with fill-in and threshold.
36 */
37 
38 #ifndef GMM_PRECOND_ILDLTT_H
39 #define GMM_PRECOND_ILDLTT_H
40 
41 // Store U = LT and D in indiag. On each line, the fill-in is the number
42 // of non-zero elements on the line of the original matrix plus K, except if
43 // the matrix is dense. In this case the fill-in is K on each line.
44 
45 #include "gmm_precond_ilut.h"
46 
47 namespace gmm {
48  /** incomplete LDL^t (cholesky) preconditioner with fill-in and
49  threshold. */
50  template <typename Matrix>
52  public :
53  typedef typename linalg_traits<Matrix>::value_type value_type;
54  typedef typename number_traits<value_type>::magnitude_type magnitude_type;
55 
57 
58  row_matrix<svector> U;
59  std::vector<magnitude_type> indiag;
60 
61  protected:
62  size_type K;
63  double eps;
64 
65  template<typename M> void do_ildltt(const M&, row_major);
66  void do_ildltt(const Matrix&, col_major);
67 
68  public:
69  void build_with(const Matrix& A, int k_ = -1, double eps_ = double(-1)) {
70  if (k_ >= 0) K = k_;
71  if (eps_ >= double(0)) eps = eps_;
72  gmm::resize(U, mat_nrows(A), mat_ncols(A));
73  indiag.resize(std::min(mat_nrows(A), mat_ncols(A)));
74  do_ildltt(A, typename principal_orientation_type<typename
75  linalg_traits<Matrix>::sub_orientation>::potype());
76  }
77  ildltt_precond(const Matrix& A, int k_, double eps_)
78  : U(mat_nrows(A),mat_ncols(A)), K(k_), eps(eps_) { build_with(A); }
79  ildltt_precond(void) { K=10; eps = 1E-7; }
80  ildltt_precond(size_type k_, double eps_) : K(k_), eps(eps_) {}
81  size_type memsize() const {
82  return sizeof(*this) + nnz(U)*sizeof(value_type) + indiag.size() * sizeof(magnitude_type);
83  }
84  };
85 
86  template<typename Matrix> template<typename M>
87  void ildltt_precond<Matrix>::do_ildltt(const M& A,row_major) {
88  typedef value_type T;
89  typedef typename number_traits<T>::magnitude_type R;
90 
91  size_type n = mat_nrows(A);
92  if (n == 0) return;
93  svector w(n);
94  T tmp;
95  R prec = default_tol(R()), max_pivot = gmm::abs(A(0,0)) * prec;
96 
97  gmm::clear(U);
98  for (size_type i = 0; i < n; ++i) {
99  gmm::copy(mat_const_row(A, i), w);
100  double norm_row = gmm::vect_norm2(w);
101 
102  for (size_type krow = 0, k; krow < w.nb_stored(); ++krow) {
103  typename svector::iterator wk = w.begin() + krow;
104  if ((k = wk->c) >= i) break;
105  if (gmm::is_complex(wk->e)) {
106  tmp = gmm::conj(U(k, i))/indiag[k]; // not completely satisfactory ..
107  gmm::add(scaled(mat_row(U, k), -tmp), w);
108  }
109  else {
110  tmp = wk->e;
111  if (gmm::abs(tmp) < eps * norm_row) { w.sup(k); --krow; }
112  else { wk->e += tmp; gmm::add(scaled(mat_row(U, k), -tmp), w); }
113  }
114  }
115  tmp = w[i];
116 
117  if (gmm::abs(gmm::real(tmp)) <= max_pivot)
118  { GMM_WARNING2("pivot " << i << " is too small"); tmp = T(1); }
119 
120  max_pivot = std::max(max_pivot, std::min(gmm::abs(tmp) * prec, R(1)));
121  indiag[i] = R(1) / gmm::real(tmp);
122  gmm::clean(w, eps * norm_row);
123  gmm::scale(w, T(indiag[i]));
124  std::sort(w.begin(), w.end(), elt_rsvector_value_less_<T>());
125  typename svector::const_iterator wit = w.begin(), wite = w.end();
126  for (size_type nnu = 0; wit != wite; ++wit) // copy to be optimized ...
127  if (wit->c > i) { if (nnu < K) { U(i, wit->c) = wit->e; ++nnu; } }
128  }
129  }
130 
131  template<typename Matrix>
132  void ildltt_precond<Matrix>::do_ildltt(const Matrix& A, col_major)
133  { do_ildltt(gmm::conjugated(A), row_major()); }
134 
135  template <typename Matrix, typename V1, typename V2> inline
136  void mult(const ildltt_precond<Matrix>& P, const V1 &v1, V2 &v2) {
137  gmm::copy(v1, v2);
138  gmm::lower_tri_solve(gmm::conjugated(P.U), v2, true);
139  for (size_type i = 0; i < P.indiag.size(); ++i) v2[i] *= P.indiag[i];
140  gmm::upper_tri_solve(P.U, v2, true);
141  }
142 
143  template <typename Matrix, typename V1, typename V2> inline
144  void transposed_mult(const ildltt_precond<Matrix>& P,const V1 &v1, V2 &v2)
145  { mult(P, v1, v2); }
146 
147  template <typename Matrix, typename V1, typename V2> inline
148  void left_mult(const ildltt_precond<Matrix>& P, const V1 &v1, V2 &v2) {
149  copy(v1, v2);
150  gmm::lower_tri_solve(gmm::conjugated(P.U), v2, true);
151  for (size_type i = 0; i < P.indiag.size(); ++i) v2[i] *= P.indiag[i];
152  }
153 
154  template <typename Matrix, typename V1, typename V2> inline
155  void right_mult(const ildltt_precond<Matrix>& P, const V1 &v1, V2 &v2)
156  { copy(v1, v2); gmm::upper_tri_solve(P.U, v2, true); }
157 
158  template <typename Matrix, typename V1, typename V2> inline
159  void transposed_left_mult(const ildltt_precond<Matrix>& P, const V1 &v1,
160  V2 &v2) {
161  copy(v1, v2);
162  gmm::upper_tri_solve(P.U, v2, true);
163  for (size_type i = 0; i < P.indiag.size(); ++i) v2[i] *= P.indiag[i];
164  }
165 
166  template <typename Matrix, typename V1, typename V2> inline
167  void transposed_right_mult(const ildltt_precond<Matrix>& P, const V1 &v1,
168  V2 &v2)
169  { copy(v1, v2); gmm::lower_tri_solve(gmm::conjugated(P.U), v2, true); }
170 
171 }
172 
173 #endif
174 
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::ildltt_precond
incomplete LDL^t (cholesky) preconditioner with fill-in and threshold.
Definition: gmm_precond_ildltt.h:51
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::copy
void copy(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:977