Reference documentation for deal.II version 8.1.0
sparse_ilu.templates.h
1 // ---------------------------------------------------------------------
2 // @f$Id: sparse_ilu.templates.h 30036 2013-07-18 16:55:32Z maier @f$
3 //
4 // Copyright (C) 1999 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__sparse_ilu_templates_h
18 #define __deal2__sparse_ilu_templates_h
19 
20 
21 
22 #include <deal.II/base/config.h>
23 #include <deal.II/lac/vector.h>
24 #include <deal.II/lac/sparse_ilu.h>
25 
26 #include <algorithm>
27 #include <cmath>
28 
29 
31 
32 template <typename number>
34 {}
35 
36 
37 
38 template <typename number>
40 {
42 }
43 
44 
45 
46 
47 template <typename number>
48 template <typename somenumber>
50  const AdditionalData &data)
51 {
53 
54  decompose(matrix, data.strengthen_diagonal);
55 }
56 
57 
58 
59 template <typename number>
60 template <typename somenumber>
62  const double strengthen_diagonal)
63 {
64  Assert (matrix.m()==matrix.n(), ExcNotQuadratic ());
65  Assert (this->m()==this->n(), ExcNotQuadratic ());
66  Assert (matrix.m()==this->m(), ExcDimensionMismatch(matrix.m(), this->m()));
67 
68  Assert (strengthen_diagonal>=0,
69  ExcInvalidStrengthening (strengthen_diagonal));
70 
71  SparseLUDecomposition<number>::decompose (matrix, strengthen_diagonal);
72 
73  if (strengthen_diagonal>0)
74  this->strengthen_diagonal_impl();
75 
76  // in the following, we implement algorithm 10.4 in the book by Saad by
77  // translating in essence the algorithm given at the end of section 10.3.2,
78  // using the names of variables used there
79  const SparsityPattern &sparsity = this->get_sparsity_pattern();
80  const std::size_t *const ia = sparsity.rowstart;
81  const size_type *const ja = sparsity.colnums;
82 
83  number *luval = this->SparseMatrix<number>::val;
84 
85  const size_type N = this->m();
86  size_type jrow = 0;
87 
88  std::vector<size_type> iw (N, numbers::invalid_size_type);
89 
90  for (size_type k=0; k<N; ++k)
91  {
92  const size_type j1 = ia[k],
93  j2 = ia[k+1]-1;
94 
95  for (size_type j=j1; j<=j2; ++j)
96  iw[ja[j]] = j;
97 
98  // the algorithm in the book works on the elements of row k left of the
99  // diagonal. however, since we store the diagonal element at the first
100  // position, start at the element after the diagonal and run as long as
101  // we don't walk into the right half
102  size_type j = j1+1;
103 
104  // pathological case: the current row of the matrix has only the
105  // diagonal entry. then we have nothing to do.
106  if (j > j2)
107  goto label_200;
108 
109 label_150:
110 
111  jrow = ja[j];
112  if (jrow >= k)
113  goto label_200;
114 
115  // actual computations:
116  {
117  number t1 = luval[j] * luval[ia[jrow]];
118  luval[j] = t1;
119 
120  // jj runs from just right of the diagonal to the end of the row
121  size_type jj = ia[jrow]+1;
122  while (ja[jj] < jrow)
123  ++jj;
124  for (; jj<ia[jrow+1]; ++jj)
125  {
126  const size_type jw = iw[ja[jj]];
127  if (jw != numbers::invalid_size_type)
128  luval[jw] -= t1 * luval[jj];
129  }
130 
131  ++j;
132  if (j<=j2)
133  goto label_150;
134  }
135 
136 label_200:
137 
138  // in the book there is an assertion that we have hit the diagonal
139  // element, i.e. that jrow==k. however, we store the diagonal element at
140  // the front, so jrow must actually be larger than k or j is already in
141  // the next row
142  Assert ((jrow > k) || (j==ia[k+1]), ExcInternalError());
143 
144  // now we have to deal with the diagonal element. in the book it is
145  // located at position 'j', but here we use the convention of storing
146  // the diagonal element first, so instead of j we use uptr[k]=ia[k]
147  Assert (luval[ia[k]] != 0, ExcInternalError());
148 
149  luval[ia[k]] = 1./luval[ia[k]];
150 
151  for (size_type j=j1; j<=j2; ++j)
152  iw[ja[j]] = numbers::invalid_size_type;
153  }
154 }
155 
156 
157 
158 
159 template <typename number>
160 template <typename somenumber>
162  const Vector<somenumber> &src) const
163 {
164  Assert (dst.size() == src.size(), ExcDimensionMismatch(dst.size(), src.size()));
165  Assert (dst.size() == this->m(), ExcDimensionMismatch(dst.size(), this->m()));
166 
167  const size_type N=dst.size();
168  const std::size_t *const rowstart_indices
169  = this->get_sparsity_pattern().rowstart;
170  const size_type *const column_numbers
171  = this->get_sparsity_pattern().colnums;
172 
173  // solve LUx=b in two steps:
174  // first Ly = b, then
175  // Ux = y
176  //
177  // first a forward solve. since
178  // the diagonal values of L are
179  // one, there holds
180  // y_i = b_i
181  // - sum_{j=0}^{i-1} L_{ij}y_j
182  // we split the y_i = b_i off and
183  // perform it at the outset of the
184  // loop
185  dst = src;
186  for (size_type row=0; row<N; ++row)
187  {
188  // get start of this row. skip the
189  // diagonal element
190  const size_type *const rowstart = &column_numbers[rowstart_indices[row]+1];
191  // find the position where the part
192  // right of the diagonal starts
193  const size_type *const first_after_diagonal = this->prebuilt_lower_bound[row];
194 
195  somenumber dst_row = dst(row);
196  const number *luval = this->SparseMatrix<number>::val +
197  (rowstart - column_numbers);
198  for (const size_type *col=rowstart; col!=first_after_diagonal; ++col, ++luval)
199  dst_row -= *luval * dst(*col);
200  dst(row) = dst_row;
201  }
202 
203  // now the backward solve. same
204  // procedure, but we need not set
205  // dst before, since this is already
206  // done.
207  //
208  // note that we need to scale now,
209  // since the diagonal is not equal to
210  // one now
211  for (int row=N-1; row>=0; --row)
212  {
213  // get end of this row
214  const size_type *const rowend = &column_numbers[rowstart_indices[row+1]];
215  // find the position where the part
216  // right of the diagonal starts
217  const size_type *const first_after_diagonal = this->prebuilt_lower_bound[row];
218 
219  somenumber dst_row = dst(row);
220  const number *luval = this->SparseMatrix<number>::val +
221  (first_after_diagonal - column_numbers);
222  for (const size_type *col=first_after_diagonal; col!=rowend; ++col, ++luval)
223  dst_row -= *luval * dst(*col);
224 
225  // scale by the diagonal element.
226  // note that the diagonal element
227  // was stored inverted
228  dst(row) = dst_row * this->diag_element(row);
229  }
230 }
231 
232 
233 template <typename number>
234 template <typename somenumber>
236  const Vector<somenumber> &src) const
237 {
238  Assert (dst.size() == src.size(), ExcDimensionMismatch(dst.size(), src.size()));
239  Assert (dst.size() == this->m(), ExcDimensionMismatch(dst.size(), this->m()));
240 
241  const size_type N=dst.size();
242  const std::size_t *const rowstart_indices
243  = this->get_sparsity_pattern().rowstart;
244  const size_type *const column_numbers
245  = this->get_sparsity_pattern().colnums;
246 
247  // solve (LU)'x=b in two steps:
248  // first U'y = b, then
249  // L'x = y
250  //
251  // first a forward solve. Due to the
252  // fact that the transpose of U'
253  // is not easily accessible, a
254  // temporary vector is required.
255  Vector<somenumber> tmp (N);
256 
257  dst = src;
258  for (size_type row=0; row<N; ++row)
259  {
260  dst(row) -= tmp (row);
261  // scale by the diagonal element.
262  // note that the diagonal element
263  // was stored inverted
264  dst(row) *= this->diag_element(row);
265 
266  // get end of this row
267  const size_type *const rowend = &column_numbers[rowstart_indices[row+1]];
268  // find the position where the part
269  // right of the diagonal starts
270  const size_type *const first_after_diagonal = this->prebuilt_lower_bound[row];
271 
272  const somenumber dst_row = dst (row);
273  const number *luval = this->SparseMatrix<number>::val +
274  (first_after_diagonal - column_numbers);
275  for (const size_type *col=first_after_diagonal; col!=rowend; ++col, ++luval)
276  tmp(*col) += *luval * dst_row;
277  }
278 
279  // now the backward solve. same
280  // procedure, but we need not set
281  // dst before, since this is already
282  // done.
283  //
284  // note that we no scaling is required
285  // now, since the diagonal is one
286  // now
287  tmp = 0;
288  for (int row=N-1; row>=0; --row)
289  {
290  dst(row) -= tmp (row);
291 
292  // get start of this row. skip the
293  // diagonal element
294  const size_type *const rowstart = &column_numbers[rowstart_indices[row]+1];
295  // find the position where the part
296  // right of the diagonal starts
297  const size_type *const first_after_diagonal = this->prebuilt_lower_bound[row];
298 
299  const somenumber dst_row = dst (row);
300  const number *luval = this->SparseMatrix<number>::val +
301  (rowstart - column_numbers);
302  for (const size_type *col=rowstart; col!=first_after_diagonal; ++col, ++luval)
303  tmp(*col) += *luval * dst_row;
304  }
305 }
306 
307 
308 template <typename number>
309 std::size_t
311 {
313 }
314 
315 
316 
317 /*---------------------------- sparse_ilu.templates.h ---------------------------*/
318 
319 DEAL_II_NAMESPACE_CLOSE
320 
321 #endif
322 /*---------------------------- sparse_ilu.templates.h ---------------------------*/
const types::global_dof_index invalid_size_type
Definition: types.h:211
void initialize(const SparseMatrix< somenumber > &matrix, const AdditionalData &parameters=AdditionalData())
size_type m() const
size_type n() const
void initialize(const SparseMatrix< somenumber > &matrix, const AdditionalData parameters)
void decompose(const SparseMatrix< somenumber > &matrix, const double strengthen_diagonal=0.) DEAL_II_DEPRECATED
#define Assert(cond, exc)
Definition: exceptions.h:299
std::size_t size() const
virtual void reinit(const SparsityPattern &sparsity)
std::size_t * rowstart
virtual std::size_t memory_consumption() const
types::global_dof_index size_type
std::size_t memory_consumption() const
void Tvmult(Vector< somenumber > &dst, const Vector< somenumber > &src) const
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
::ExceptionBase & ExcInternalError()
void decompose(const SparseMatrix< somenumber > &matrix, const double strengthen_diagonal=0.) DEAL_II_DEPRECATED
void vmult(Vector< somenumber > &dst, const Vector< somenumber > &src) const