Reference documentation for deal.II version 8.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
sparse_decomposition.templates.h
1 // ---------------------------------------------------------------------
2 // @f$Id: sparse_decomposition.templates.h 30040 2013-07-18 17:06:48Z 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 
18 #ifndef __deal2__sparse_decomposition_templates_h
19 #define __deal2__sparse_decomposition_templates_h
20 
21 
22 #include <deal.II/base/memory_consumption.h>
23 #include <deal.II/base/template_constraints.h>
24 #include <deal.II/base/utilities.h>
25 #include <deal.II/lac/sparse_decomposition.h>
26 #include <algorithm>
27 #include <cstring>
28 
29 DEAL_II_NAMESPACE_OPEN
30 
31 template<typename number>
33  :
34  SparseMatrix<number>(),
35  decomposed(false),
36  own_sparsity(0)
37 {}
38 
39 
40 
41 template<typename number>
44  SparseMatrix<number>(sparsity),
45  decomposed(false),
46  own_sparsity(0)
47 {}
48 
49 
50 
51 template<typename number>
53 {
54  clear();
55 }
56 
57 
58 template<typename number>
60 {
61  decomposed = false;
62 
63  std::vector<const size_type *> tmp;
64  tmp.swap (prebuilt_lower_bound);
65 
67 
68  if (own_sparsity)
69  {
70  delete own_sparsity;
71  own_sparsity=0;
72  }
73 }
74 
75 
76 
77 template<typename number>
78 template <typename somenumber>
80  const SparseMatrix<somenumber> &matrix,
81  const AdditionalData data)
82 {
83  const SparsityPattern &matrix_sparsity=matrix.get_sparsity_pattern();
84 
85  const SparsityPattern *sparsity_pattern_to_use = 0;
86 
87  if (data.use_this_sparsity)
88  sparsity_pattern_to_use = data.use_this_sparsity;
89  else if (data.use_previous_sparsity &&
90  !this->empty() &&
91  (this->m()==matrix.m()))
92  {
93  // Use the sparsity that was
94  // previously used. This is
95  // particularly useful when
96  // matrix entries change but
97  // not the sparsity, as for the
98  // case of several Newton
99  // iteration steps on an
100  // unchanged grid.
101  sparsity_pattern_to_use = &this->get_sparsity_pattern();
102  }
103  else if (data.extra_off_diagonals==0)
104  {
105  // Use same sparsity as matrix
106  sparsity_pattern_to_use = &matrix_sparsity;
107  }
108  else
109  {
110  // Create new sparsity
111 
112  // for the case that
113  // own_sparsity wasn't deleted
114  // before (e.g. by clear()), do
115  // it here
116  if (own_sparsity)
117  {
118  // release the sparsity
120  // delete it
121  delete own_sparsity;
122  }
123 
124  // and recreate
125  own_sparsity = new SparsityPattern(matrix_sparsity,
126  matrix_sparsity.max_entries_per_row()
127  +2*data.extra_off_diagonals,
128  data.extra_off_diagonals);
129  own_sparsity->compress();
130  sparsity_pattern_to_use = own_sparsity;
131  }
132 
133  // now use this sparsity pattern
134  Assert (sparsity_pattern_to_use->n_rows()==sparsity_pattern_to_use->n_cols(),
135  typename SparsityPattern::ExcDiagonalNotOptimized());
136  decomposed = false;
137  {
138  std::vector<const size_type *> tmp;
139  tmp.swap (prebuilt_lower_bound);
140  }
141  SparseMatrix<number>::reinit (*sparsity_pattern_to_use);
142 }
143 
144 
145 template<typename number>
146 template<typename somenumber>
147 void
150  const double strengthen_diagonal)
151 {
152  decomposed = false;
153 
154  this->strengthen_diagonal = strengthen_diagonal;
155  prebuild_lower_bound ();
156  copy_from (matrix);
157  decomposed = true;
158 }
159 
160 
161 
162 template <typename number>
164 {
165  Assert (sparsity.n_rows() == sparsity.n_cols(),
166  typename SparsityPattern::ExcDiagonalNotOptimized());
167  decomposed = false;
168  {
169  std::vector<const size_type *> tmp;
170  tmp.swap (prebuilt_lower_bound);
171  }
172  SparseMatrix<number>::reinit (sparsity);
173 }
174 
175 
176 
177 template<typename number>
178 void
180 {
181  const size_type *const
182  column_numbers = this->get_sparsity_pattern().colnums;
183  const std::size_t *const
184  rowstart_indices = this->get_sparsity_pattern().rowstart;
185  const size_type N = this->m();
186 
187  prebuilt_lower_bound.resize (N);
188 
189  for (size_type row=0; row<N; row++)
190  {
191  prebuilt_lower_bound[row]
192  = Utilities::lower_bound (&column_numbers[rowstart_indices[row]+1],
193  &column_numbers[rowstart_indices[row+1]],
194  row);
195  }
196 }
197 
198 template <typename number>
199 template <typename somenumber>
200 void
202 {
203  // check whether we use the same sparsity
204  // pattern as the input matrix
205  if (&this->get_sparsity_pattern() == &matrix.get_sparsity_pattern())
206  {
207  const somenumber *input_ptr = matrix.val;
208  number *this_ptr = this->val;
209  const number *const end_ptr = this_ptr + this->n_nonzero_elements();
211  std::memcpy (this_ptr, input_ptr, this->n_nonzero_elements()*sizeof(number));
212  else
213  for ( ; this_ptr != end_ptr; ++input_ptr, ++this_ptr)
214  *this_ptr = *input_ptr;
215  return;
216  }
217 
218  // preset the elements by zero. this needs to be written in a slightly
219  // awkward way so that we find the corresponding function in the base class.
221 
222  // both allow more and less entries in the new matrix
223  for (size_type row=0; row<this->m(); ++row)
224  {
225  typename SparseMatrix<number>::iterator index = this->begin(row);
227  in_index = matrix.begin(row);
228  index->value() = in_index->value();
229  ++index, ++in_index;
230  while (index < this->end(row) && in_index < matrix.end(row))
231  {
232  while (index->column() < in_index->column() && index < this->end(row))
233  ++index;
234  while (in_index->column() < index->column() && in_index < matrix.end(row))
235  ++in_index;
236 
237  index->value() = in_index->value();
238  ++index, ++in_index;
239  }
240  }
241 }
242 
243 
244 
245 template <typename number>
246 void
248 {
249  for (size_type row=0; row<this->m(); ++row)
250  {
251  // get the global index of the first
252  // non-diagonal element in this row
253  Assert (this->m() == this->n(), ExcNotImplemented());
255  diagonal_element = this->begin(row);
256 
257  number rowsum = 0;
258  for (typename SparseMatrix<number>::iterator
259  p = diagonal_element + 1;
260  p != this->end(row); ++p)
261  rowsum += std::fabs(p->value());
262 
263  diagonal_element->value() += this->get_strengthen_diagonal (rowsum, row) *
264  rowsum;
265  }
266 }
267 
268 
269 
270 template <typename number>
271 std::size_t
273 {
275  MemoryConsumption::memory_consumption(prebuilt_lower_bound));
276 }
277 
278 
279 DEAL_II_NAMESPACE_CLOSE
280 
281 #endif // __deal2__sparse_decomposition_templates_h
virtual void clear()
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:731
types::global_dof_index size_type
void copy_from(const SparseMatrix< somenumber > &matrix)
size_type m() const
const_iterator begin() const
void initialize(const SparseMatrix< somenumber > &matrix, const AdditionalData parameters)
void decompose(const SparseMatrix< somenumber > &matrix, const double strengthen_diagonal=0.) DEAL_II_DEPRECATED
SparseMatrix< number > & operator=(const SparseMatrix< number > &)
#define Assert(cond, exc)
Definition: exceptions.h:299
std::size_t memory_consumption(const T &t)
const SparsityPattern & get_sparsity_pattern() const
size_type n_cols() const
void reinit(const SparsityPattern &sparsity) DEAL_II_DEPRECATED
const_iterator end() const
virtual void reinit(const SparsityPattern &sparsity)
size_type max_entries_per_row() const
virtual std::size_t memory_consumption() const
size_type n_rows() const
::ExceptionBase & ExcNotImplemented()