17 #ifndef __deal2__sparse_mic_templates_h
18 #define __deal2__sparse_mic_templates_h
21 #include <deal.II/base/memory_consumption.h>
22 #include <deal.II/lac/sparse_mic.h>
23 #include <deal.II/lac/vector.h>
27 template <
typename number>
37 template <
typename number>
48 template <
typename number>
55 template <
typename number>
59 std::vector<number> tmp;
63 std::vector<number> tmp;
67 std::vector<number> tmp;
68 tmp.swap (inner_sums);
75 template <
typename number>
76 template <
typename somenumber>
88 template <
typename number>
92 std::vector<number> tmp;
96 std::vector<number> tmp;
100 std::vector<number> tmp;
101 tmp.swap (inner_sums);
105 this->decomposed =
false;
110 template <
typename number>
111 template <
typename somenumber>
113 const double strengthen_diagonal)
117 Assert (matrix.
m()==matrix.
n(), ExcNotQuadratic ());
118 Assert (this->m()==this->n(), ExcNotQuadratic ());
121 Assert (strengthen_diagonal>=0, ExcInvalidStrengthening (strengthen_diagonal));
123 if (strengthen_diagonal > 0)
124 this->strengthen_diagonal_impl ();
135 diag.resize (this->m());
136 inv_diag.resize (this->m());
137 inner_sums.resize (this->m());
140 for (
size_type row=0; row<this->m(); row++)
141 inner_sums[row] = get_rowsum(row);
143 for (
size_type row=0; row<this->m(); row++)
145 const number temp = this->begin(row)->value();
151 p = matrix.
begin(row)+1;
152 (p != matrix.
end(row)) && (p->column() < row);
154 temp1 += p->value() / diag[p->column()] * inner_sums[p->column()];
156 Assert(temp-temp1 > 0, ExcStrengthenDiagonalTooSmall());
157 diag[row] = temp - temp1;
159 inv_diag[row] = 1.0/diag[row];
165 template <
typename number>
169 Assert(this->m()==this->n(), ExcNotQuadratic());
173 p = this->begin(row)+1;
174 p != this->end(row); ++p)
175 if (p->column() > row)
176 rowsum += p->value();
183 template <
typename number>
184 template <
typename somenumber>
206 p = this->begin(row)+1;
207 (p != this->end(row)) && (p->column() < row);
209 dst(row) -= p->value() * dst(p->column());
211 dst(row) *= inv_diag[row];
216 dst(row) *= diag[row];
219 for (
int row=N-1; row>=0; --row)
223 p = this->begin(row)+1;
226 if (p->column() >
static_cast<size_type>(row))
227 dst(row) -= p->value() * dst(p->column());
229 dst(row) *= inv_diag[row];
235 template <
typename number>
247 DEAL_II_NAMESPACE_CLOSE
249 #endif // __deal2__sparse_mic_templates_h
void initialize(const SparseMatrix< somenumber > &matrix, const AdditionalData ¶meters=AdditionalData())
void reinit(const SparsityPattern &sparsity) DEAL_II_DEPRECATED
double strengthen_diagonal
virtual void reinit(const SparsityPattern &sparsity)
#define Assert(cond, exc)
std::size_t memory_consumption(const T &t)
std::size_t memory_consumption() const
void vmult(Vector< somenumber > &dst, const Vector< somenumber > &src) const
const_iterator end() const
void decompose(const SparseMatrix< somenumber > &matrix, const double strengthen_diagonal=0.) DEAL_II_DEPRECATED
number get_rowsum(const size_type row) const
void decompose(const SparseMatrix< somenumber > &matrix, const double strengthen_diagonal=0.) DEAL_II_DEPRECATED
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
types::global_dof_index size_type
void vmult(OutVector &dst, const InVector &src) const
void initialize(const SparseMatrix< somenumber > &matrix, const AdditionalData parameters)
const_iterator begin() const