17 #ifndef __deal2__sparse_ilu_templates_h
18 #define __deal2__sparse_ilu_templates_h
22 #include <deal.II/base/config.h>
23 #include <deal.II/lac/vector.h>
24 #include <deal.II/lac/sparse_ilu.h>
32 template <
typename number>
38 template <
typename number>
47 template <
typename number>
48 template <
typename somenumber>
59 template <
typename number>
60 template <
typename somenumber>
62 const double strengthen_diagonal)
64 Assert (matrix.
m()==matrix.
n(), ExcNotQuadratic ());
65 Assert (this->m()==this->n(), ExcNotQuadratic ());
68 Assert (strengthen_diagonal>=0,
69 ExcInvalidStrengthening (strengthen_diagonal));
73 if (strengthen_diagonal>0)
74 this->strengthen_diagonal_impl();
80 const std::size_t *
const ia = sparsity.
rowstart;
117 number t1 = luval[j] * luval[ia[jrow]];
122 while (ja[jj] < jrow)
124 for (; jj<ia[jrow+1]; ++jj)
128 luval[jw] -= t1 * luval[jj];
149 luval[ia[k]] = 1./luval[ia[k]];
159 template <
typename number>
160 template <
typename somenumber>
168 const std::size_t *
const rowstart_indices
169 = this->get_sparsity_pattern().rowstart;
171 = this->get_sparsity_pattern().colnums;
190 const size_type *
const rowstart = &column_numbers[rowstart_indices[row]+1];
193 const size_type *
const first_after_diagonal = this->prebuilt_lower_bound[row];
195 somenumber dst_row = dst(row);
197 (rowstart - column_numbers);
198 for (
const size_type *col=rowstart; col!=first_after_diagonal; ++col, ++luval)
199 dst_row -= *luval * dst(*col);
211 for (
int row=N-1; row>=0; --row)
214 const size_type *
const rowend = &column_numbers[rowstart_indices[row+1]];
217 const size_type *
const first_after_diagonal = this->prebuilt_lower_bound[row];
219 somenumber dst_row = dst(row);
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);
228 dst(row) = dst_row * this->diag_element(row);
233 template <
typename number>
234 template <
typename somenumber>
242 const std::size_t *
const rowstart_indices
243 = this->get_sparsity_pattern().rowstart;
245 = this->get_sparsity_pattern().colnums;
260 dst(row) -= tmp (row);
264 dst(row) *= this->diag_element(row);
267 const size_type *
const rowend = &column_numbers[rowstart_indices[row+1]];
270 const size_type *
const first_after_diagonal = this->prebuilt_lower_bound[row];
272 const somenumber dst_row = dst (row);
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;
288 for (
int row=N-1; row>=0; --row)
290 dst(row) -= tmp (row);
294 const size_type *
const rowstart = &column_numbers[rowstart_indices[row]+1];
297 const size_type *
const first_after_diagonal = this->prebuilt_lower_bound[row];
299 const somenumber dst_row = dst (row);
301 (rowstart - column_numbers);
302 for (
const size_type *col=rowstart; col!=first_after_diagonal; ++col, ++luval)
303 tmp(*col) += *luval * dst_row;
308 template <
typename number>
319 DEAL_II_NAMESPACE_CLOSE
const types::global_dof_index invalid_size_type
void initialize(const SparseMatrix< somenumber > &matrix, const AdditionalData ¶meters=AdditionalData())
types::global_dof_index size_type
double strengthen_diagonal
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)
virtual void reinit(const SparsityPattern &sparsity)
virtual std::size_t memory_consumption() const
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