GetFEM++  5.3
gmm_precond_ildlt.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2003-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 cholesky.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 #ifndef GMM_PRECOND_ILDLT_H
65 #define GMM_PRECOND_ILDLT_H
66 
67 /**@file gmm_precond_ildlt.h
68  @author Andrew Lumsdaine <lums@osl.iu.edu>
69  @author Lie-Quan Lee <llee@osl.iu.edu>
70  @author Yves Renard <yves.renard@insa-lyon.fr>
71  @date June 5, 2003.
72  @brief Incomplete Level 0 ILDLT Preconditioner.
73 */
74 
75 #include "gmm_precond.h"
76 
77 namespace gmm {
78 
79  /** Incomplete Level 0 LDLT Preconditioner.
80 
81  For use with symmetric real or hermitian complex sparse matrices.
82 
83  Notes: The idea under a concrete Preconditioner such as Incomplete
84  Cholesky is to create a Preconditioner object to use in iterative
85  methods.
86 
87 
88  Y. Renard : Transformed in LDLT for stability reason.
89 
90  U=LT is stored in csr format. D is stored on the diagonal of U.
91  */
92  template <typename Matrix>
93  class ildlt_precond {
94 
95  public :
96  typedef typename linalg_traits<Matrix>::value_type value_type;
97  typedef typename number_traits<value_type>::magnitude_type magnitude_type;
98  typedef csr_matrix_ref<value_type *, size_type *, size_type *, 0> tm_type;
99 
100  tm_type U;
101 
102  protected :
103  std::vector<value_type> Tri_val;
104  std::vector<size_type> Tri_ind, Tri_ptr;
105 
106  template<typename M> void do_ildlt(const M& A, row_major);
107  void do_ildlt(const Matrix& A, col_major);
108 
109  public:
110 
111  size_type nrows(void) const { return mat_nrows(U); }
112  size_type ncols(void) const { return mat_ncols(U); }
113  value_type &D(size_type i) { return Tri_val[Tri_ptr[i]]; }
114  const value_type &D(size_type i) const { return Tri_val[Tri_ptr[i]]; }
115  ildlt_precond(void) {}
116  void build_with(const Matrix& A) {
117  Tri_ptr.resize(mat_nrows(A)+1);
118  do_ildlt(A, typename principal_orientation_type<typename
119  linalg_traits<Matrix>::sub_orientation>::potype());
120  }
121  ildlt_precond(const Matrix& A) { build_with(A); }
122  size_type memsize() const {
123  return sizeof(*this) +
124  Tri_val.size() * sizeof(value_type) +
125  (Tri_ind.size()+Tri_ptr.size()) * sizeof(size_type);
126  }
127  };
128 
129  template <typename Matrix> template<typename M>
130  void ildlt_precond<Matrix>::do_ildlt(const M& A, row_major) {
131  typedef typename linalg_traits<Matrix>::storage_type store_type;
132  typedef value_type T;
133  typedef typename number_traits<T>::magnitude_type R;
134 
135  size_type Tri_loc = 0, n = mat_nrows(A), d, g, h, i, j, k;
136  if (n == 0) return;
137  T z, zz;
138  Tri_ptr[0] = 0;
139  R prec = default_tol(R());
140  R max_pivot = gmm::abs(A(0,0)) * prec;
141 
142  for (int count = 0; count < 2; ++count) {
143  if (count) { Tri_val.resize(Tri_loc); Tri_ind.resize(Tri_loc); }
144  for (Tri_loc = 0, i = 0; i < n; ++i) {
145  typedef typename linalg_traits<M>::const_sub_row_type row_type;
146  row_type row = mat_const_row(A, i);
147  typename linalg_traits<typename org_type<row_type>::t>::const_iterator
148  it = vect_const_begin(row), ite = vect_const_end(row);
149 
150  if (count) { Tri_val[Tri_loc] = T(0); Tri_ind[Tri_loc] = i; }
151  ++Tri_loc; // diagonal element
152 
153  for (k = 0; it != ite; ++it, ++k) {
154  j = index_of_it(it, k, store_type());
155  if (i == j) {
156  if (count) Tri_val[Tri_loc-1] = *it;
157  }
158  else if (j > i) {
159  if (count) { Tri_val[Tri_loc] = *it; Tri_ind[Tri_loc]=j; }
160  ++Tri_loc;
161  }
162  }
163  Tri_ptr[i+1] = Tri_loc;
164  }
165  }
166 
167  if (A(0,0) == T(0)) {
168  Tri_val[Tri_ptr[0]] = T(1);
169  GMM_WARNING2("pivot 0 is too small");
170  }
171 
172  for (k = 0; k < n; k++) {
173  d = Tri_ptr[k];
174  z = T(gmm::real(Tri_val[d])); Tri_val[d] = z;
175  if (gmm::abs(z) <= max_pivot) {
176  Tri_val[d] = z = T(1);
177  GMM_WARNING2("pivot " << k << " is too small [" << gmm::abs(z) << "]");
178  }
179  max_pivot = std::max(max_pivot, std::min(gmm::abs(z) * prec, R(1)));
180 
181  for (i = d + 1; i < Tri_ptr[k+1]; ++i) Tri_val[i] /= z;
182  for (i = d + 1; i < Tri_ptr[k+1]; ++i) {
183  zz = gmm::conj(Tri_val[i] * z);
184  h = Tri_ind[i];
185  g = i;
186 
187  for (j = Tri_ptr[h] ; j < Tri_ptr[h+1]; ++j)
188  for ( ; g < Tri_ptr[k+1] && Tri_ind[g] <= Tri_ind[j]; ++g)
189  if (Tri_ind[g] == Tri_ind[j])
190  Tri_val[j] -= zz * Tri_val[g];
191  }
192  }
193  U = tm_type(&(Tri_val[0]), &(Tri_ind[0]), &(Tri_ptr[0]),
194  n, mat_ncols(A));
195  }
196 
197  template <typename Matrix>
198  void ildlt_precond<Matrix>::do_ildlt(const Matrix& A, col_major)
199  { do_ildlt(gmm::conjugated(A), row_major()); }
200 
201  template <typename Matrix, typename V1, typename V2> inline
202  void mult(const ildlt_precond<Matrix>& P, const V1 &v1, V2 &v2) {
203  gmm::copy(v1, v2);
204  gmm::lower_tri_solve(gmm::conjugated(P.U), v2, true);
205  for (size_type i = 0; i < mat_nrows(P.U); ++i) v2[i] /= P.D(i);
206  gmm::upper_tri_solve(P.U, v2, true);
207  }
208 
209  template <typename Matrix, typename V1, typename V2> inline
210  void transposed_mult(const ildlt_precond<Matrix>& P,const V1 &v1,V2 &v2)
211  { mult(P, v1, v2); }
212 
213  template <typename Matrix, typename V1, typename V2> inline
214  void left_mult(const ildlt_precond<Matrix>& P, const V1 &v1, V2 &v2) {
215  copy(v1, v2);
216  gmm::lower_tri_solve(gmm::conjugated(P.U), v2, true);
217  for (size_type i = 0; i < mat_nrows(P.U); ++i) v2[i] /= P.D(i);
218  }
219 
220  template <typename Matrix, typename V1, typename V2> inline
221  void right_mult(const ildlt_precond<Matrix>& P, const V1 &v1, V2 &v2)
222  { copy(v1, v2); gmm::upper_tri_solve(P.U, v2, true); }
223 
224  template <typename Matrix, typename V1, typename V2> inline
225  void transposed_left_mult(const ildlt_precond<Matrix>& P, const V1 &v1,
226  V2 &v2) {
227  copy(v1, v2);
228  gmm::upper_tri_solve(P.U, v2, true);
229  for (size_type i = 0; i < mat_nrows(P.U); ++i) v2[i] /= P.D(i);
230  }
231 
232  template <typename Matrix, typename V1, typename V2> inline
233  void transposed_right_mult(const ildlt_precond<Matrix>& P, const V1 &v1,
234  V2 &v2)
235  { copy(v1, v2); gmm::lower_tri_solve(gmm::conjugated(P.U), v2, true); }
236 
237 
238 }
239 
240 #endif
241 
conjugated_return< L >::return_type conjugated(const L &v)
return a conjugated view of the input matrix or vector.
void copy(const L1 &l1, L2 &l2)
*/
Definition: gmm_blas.h:977
gmm preconditioners.
Incomplete Level 0 LDLT Preconditioner.