Reference documentation for deal.II version 8.1.0
sparse_mic.templates.h
1 // ---------------------------------------------------------------------
2 // @f$Id: sparse_mic.templates.h 30036 2013-07-18 16:55:32Z maier @f$
3 //
4 // Copyright (C) 2002 - 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_mic_templates_h
18 #define __deal2__sparse_mic_templates_h
19 
20 
21 #include <deal.II/base/memory_consumption.h>
22 #include <deal.II/lac/sparse_mic.h>
23 #include <deal.II/lac/vector.h>
24 
26 
27 template <typename number>
29  :
30  diag(0),
31  inv_diag(0),
32  inner_sums(0)
33 {}
34 
35 
36 
37 template <typename number>
39  :
40  diag(0),
41  inv_diag(0),
42  inner_sums(0)
43 {
45 }
46 
47 
48 template <typename number>
50 {
51  clear();
52 }
53 
54 
55 template <typename number>
57 {
58  {
59  std::vector<number> tmp;
60  tmp.swap (diag);
61  }
62  {
63  std::vector<number> tmp;
64  tmp.swap (inv_diag);
65  }
66  {
67  std::vector<number> tmp;
68  tmp.swap (inner_sums);
69  }
70 
72 }
73 
74 
75 template <typename number>
76 template <typename somenumber>
77 inline
79  const AdditionalData &data)
80 {
82 
83  decompose(matrix, data.strengthen_diagonal);
84 }
85 
86 
87 
88 template <typename number>
90 {
91  {
92  std::vector<number> tmp;
93  tmp.swap (diag);
94  }
95  {
96  std::vector<number> tmp;
97  tmp.swap (inv_diag);
98  }
99  {
100  std::vector<number> tmp;
101  tmp.swap (inner_sums);
102  }
103 
105  this->decomposed = false;
106 }
107 
108 
109 
110 template <typename number>
111 template <typename somenumber>
113  const double strengthen_diagonal)
114 {
115  SparseLUDecomposition<number>::decompose(matrix, strengthen_diagonal);
116 
117  Assert (matrix.m()==matrix.n(), ExcNotQuadratic ());
118  Assert (this->m()==this->n(), ExcNotQuadratic ());
119  Assert (matrix.m()==this->m(), ExcDimensionMismatch(matrix.m(), this->m()));
120 
121  Assert (strengthen_diagonal>=0, ExcInvalidStrengthening (strengthen_diagonal));
122 
123  if (strengthen_diagonal > 0)
124  this->strengthen_diagonal_impl ();
125 
126  // MIC implementation: (S. Margenov lectures)
127  // x[i] = a[i][i] - sum(k=1, i-1,
128  // a[i][k]/x[k]*sum(j=k+1, N, a[k][j]))
129 
130  // TODO: for sake of simplicity,
131  // those are placed here. A better
132  // implementation would store this
133  // values in the underlying sparse
134  // matrix itself.
135  diag.resize (this->m());
136  inv_diag.resize (this->m());
137  inner_sums.resize (this->m());
138 
139  // precalc sum(j=k+1, N, a[k][j]))
140  for (size_type row=0; row<this->m(); row++)
141  inner_sums[row] = get_rowsum(row);
142 
143  for (size_type row=0; row<this->m(); row++)
144  {
145  const number temp = this->begin(row)->value();
146  number temp1 = 0;
147 
148  // work on the lower left part of the matrix. we know
149  // it's symmetric, so we can work with this alone
151  p = matrix.begin(row)+1;
152  (p != matrix.end(row)) && (p->column() < row);
153  ++p)
154  temp1 += p->value() / diag[p->column()] * inner_sums[p->column()];
155 
156  Assert(temp-temp1 > 0, ExcStrengthenDiagonalTooSmall());
157  diag[row] = temp - temp1;
158 
159  inv_diag[row] = 1.0/diag[row];
160  }
161 }
162 
163 
164 
165 template <typename number>
166 inline number
168 {
169  Assert(this->m()==this->n(), ExcNotQuadratic());
170 
171  number rowsum = 0;
173  p = this->begin(row)+1;
174  p != this->end(row); ++p)
175  if (p->column() > row)
176  rowsum += p->value();
177 
178  return rowsum;
179 }
180 
181 
182 
183 template <typename number>
184 template <typename somenumber>
185 void
187  const Vector<somenumber> &src) const
188 {
190  Assert (dst.size() == src.size(), ExcDimensionMismatch(dst.size(), src.size()));
191  Assert (dst.size() == this->m(), ExcDimensionMismatch(dst.size(), this->m()));
192 
193  const size_type N=dst.size();
194  // We assume the underlying matrix A is: A = X - L - U, where -L and -U are
195  // strictly lower- and upper- diagonal parts of the system.
196  //
197  // Solve (X-L)X{-1}(X-U) x = b in 3 steps:
198  dst = src;
199  for (size_type row=0; row<N; ++row)
200  {
201  // Now: (X-L)u = b
202 
203  // get start of this row. skip
204  // the diagonal element
206  p = this->begin(row)+1;
207  (p != this->end(row)) && (p->column() < row);
208  ++p)
209  dst(row) -= p->value() * dst(p->column());
210 
211  dst(row) *= inv_diag[row];
212  }
213 
214  // Now: v = Xu
215  for (size_type row=0; row<N; row++)
216  dst(row) *= diag[row];
217 
218  // x = (X-U)v
219  for (int row=N-1; row>=0; --row)
220  {
221  // get end of this row
223  p = this->begin(row)+1;
224  p != this->end(row);
225  ++p)
226  if (p->column() > static_cast<size_type>(row))
227  dst(row) -= p->value() * dst(p->column());
228 
229  dst(row) *= inv_diag[row];
230  }
231 }
232 
233 
234 
235 template <typename number>
236 std::size_t
238 {
243 }
244 
245 
246 
247 DEAL_II_NAMESPACE_CLOSE
248 
249 #endif // __deal2__sparse_mic_templates_h
void initialize(const SparseMatrix< somenumber > &matrix, const AdditionalData &parameters=AdditionalData())
void vmult(OutVector &dst, const InVector &src) const
size_type m() const
size_type n() const
const_iterator begin() const
void reinit(const SparsityPattern &sparsity) DEAL_II_DEPRECATED
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 memory_consumption(const T &t)
std::size_t memory_consumption() const
std::size_t size() const
void vmult(Vector< somenumber > &dst, const Vector< somenumber > &src) const
const_iterator end() const
types::global_dof_index size_type
Definition: sparse_mic.h:53
virtual void reinit(const SparsityPattern &sparsity)
virtual void clear()
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)
virtual ~SparseMIC()