17 #ifndef __deal2__sparse_vanka_templates_h
18 #define __deal2__sparse_vanka_templates_h
21 #include <deal.II/base/memory_consumption.h>
22 #include <deal.II/base/thread_management.h>
23 #include <deal.II/lac/sparse_vanka.h>
24 #include <deal.II/lac/full_matrix.h>
25 #include <deal.II/lac/sparse_matrix.h>
26 #include <deal.II/lac/vector.h>
33 template<
typename number>
35 const std::vector<bool> &selected,
36 const bool conserve_mem,
37 const unsigned int n_threads)
39 matrix (&M, typeid(*this).name()),
40 conserve_mem (conserve_mem),
42 n_threads (n_threads),
45 Assert (M.
m() == M.
n(), ExcNotQuadratic ());
48 if (conserve_mem ==
false)
53 template<
typename number>
57 for (i=inverses.begin(); i!=inverses.end(); ++i)
66 template <
typename number>
70 #ifndef DEAL_II_WITH_THREADS
71 compute_inverses (0, matrix->m());
73 const size_type n_inverses = std::count (selected.begin(),
77 const size_type n_inverses_per_thread = std::max(n_inverses / n_threads,
78 static_cast<size_type> (1U));
96 std::vector<std::pair<size_type, unsigned int> > blocking (n_threads);
99 unsigned int thread = 0;
100 blocking[0].first = 0;
102 for (
size_type i=0; (i<matrix->m()) && (thread+1<n_threads); ++i)
104 if (selected[i] ==
true)
106 if (c == n_inverses_per_thread)
108 blocking[thread].second = i;
109 blocking[thread+1].first = i;
115 blocking[n_threads-1].second = matrix->m();
123 for (
unsigned int i=0; i<n_threads; ++i)
132 template <
typename number>
140 std::vector<size_type> local_indices;
144 for (
size_type row=begin; row<end; ++row)
145 if (selected[row] ==
true)
146 compute_inverse (row, local_indices);
151 template <
typename number>
154 std::vector<size_type> &local_indices)
160 = matrix->get_sparsity_pattern();
168 local_indices.resize (row_length);
173 inverses[row]->extract_submatrix_from (*matrix,
178 inverses[row]->gauss_jordan();
182 template<
typename number>
183 template<
typename number2>
192 apply_preconditioner (dst, src);
197 template<
typename number>
198 template<
typename number2>
202 const std::vector<bool> *
const dof_mask)
const
213 = matrix->get_sparsity_pattern();
221 const bool range_is_restricted = (dof_mask != 0);
234 std::map<size_type, size_type> local_index;
240 if ((selected[row] ==
true) &&
241 ((range_is_restricted ==
false) || ((*dof_mask)[row] ==
true)))
250 if (conserve_mem ==
true)
252 inverses[row] = &local_matrix;
253 inverses[row]->reinit (row_length, row_length);
256 b.reinit (row_length);
257 x.reinit (row_length);
271 local_index.clear ();
273 local_index.insert(std::pair<size_type, size_type>
277 for (std::map<size_type, size_type>::const_iterator is=local_index.
begin();
278 is!=local_index.
end(); ++is)
296 p != matrix->end(row); ++p)
303 const std::map<size_type, size_type>::const_iterator js
304 = local_index.find(p->column());
310 if (js == local_index.
end())
312 if (!range_is_restricted ||
313 ((*dof_mask)[p->column()] ==
true))
314 b(i) -= p->value() * dst(p->column());
319 if (conserve_mem ==
true)
320 (*inverses[row])(i,js->second) = p->value();
325 if (conserve_mem ==
true)
326 inverses[row]->gauss_jordan();
329 inverses[row]->vmult(x,b);
332 for (std::map<size_type, size_type>::const_iterator is=local_index.
begin();
333 is!=local_index.
end(); ++is)
338 if (!range_is_restricted ||
339 ((*dof_mask)[irow] ==
true))
348 if (conserve_mem ==
true)
355 template <
typename number>
359 std::size_t mem = (
sizeof(*this) +
361 for (
size_type i=0; i<inverses.size(); ++i)
370 template <
typename number>
372 const std::vector<bool> &selected,
373 const unsigned int n_blocks,
375 const bool conserve_memory,
376 const unsigned int n_threads)
378 SparseVanka<number> (M, selected, conserve_memory, n_threads),
381 std::vector<
bool>(M.m(), false))
387 template <
typename number>
390 const std::vector<bool> &selected,
395 const size_type n_inverses = std::count (selected.begin(),
399 const size_type n_inverses_per_block = std::max(n_inverses / n_blocks,
400 static_cast<size_type> (1U));
403 std::vector<std::pair<size_type, size_type> > intervals (n_blocks);
422 unsigned int block = 0;
423 intervals[0].first = 0;
425 for (
size_type i=0; (i<M.
m()) && (block+1<n_blocks); ++i)
427 if (selected[i] ==
true)
429 if (c == n_inverses_per_block)
431 intervals[block].second = i;
432 intervals[block+1].first = i;
438 intervals[n_blocks-1].second = M.
m();
447 switch (blocking_strategy)
449 case index_intervals:
451 for (
unsigned int block=0; block<n_blocks; ++block)
452 std::fill_n (dof_masks[block].begin()+intervals[block].first,
453 intervals[block].second - intervals[block].first,
479 std::map<size_type, size_type> local_index;
483 if (selected[row] ==
true)
488 unsigned int block_number = 0;
489 while (row>=intervals[block_number].second)
501 ++access_count[block_number][structure.
column_number(row, i)];
517 if (selected[row] ==
true)
519 unsigned int block_number = 0;
520 while (row>=intervals[block_number].second)
523 dof_masks[block_number][row] =
true;
531 unsigned int max_access_block = 0;
532 for (
unsigned int block=0; block<n_blocks; ++block)
533 if (access_count[block][row] > max_accesses)
535 max_accesses = access_count[block][row];
536 max_access_block = block;
538 dof_masks[max_access_block][row] =
true;
551 template <
typename number>
552 template <
typename number2>
561 this->apply_preconditioner (dst, src);
565 #ifdef DEAL_II_WITH_THREADS
581 const Vector<number2> &,
582 const std::vector<bool> *
const)
const;
586 for (
unsigned int block=0; block<n_blocks; ++block)
589 dst, src,&dof_masks[block]);
592 for (
unsigned int block=0; block<n_blocks; ++block)
593 this->apply_preconditioner (dst, src,
601 template <
typename number>
606 for (
size_type i=0; i<dof_masks.size(); ++i)
613 DEAL_II_NAMESPACE_CLOSE
types::global_dof_index size_type
SparseBlockVanka(const SparseMatrix< number > &M, const std::vector< bool > &selected, const unsigned int n_blocks, const BlockingStrategy blocking_strategy, const bool conserve_memory=false, const unsigned int n_threads=multithread_info.n_threads())
void vmult(Vector< number2 > &dst, const Vector< number2 > &src) const
void apply_preconditioner(Vector< number2 > &dst, const Vector< number2 > &src, const std::vector< bool > *const dof_mask=0) const
Thread< RT > new_thread(const std_cxx1x::function< RT()> &function)
void compute_dof_masks(const SparseMatrix< number > &M, const std::vector< bool > &selected, const BlockingStrategy blocking_strategy)
#define Assert(cond, exc)
std::size_t memory_consumption(const T &t)
const SparsityPattern & get_sparsity_pattern() const
unsigned int row_length(const size_type row) const
size_type column_number(const size_type row, const unsigned int index) const
std::size_t memory_consumption() const
const_iterator begin() const
size_type max_entries_per_row() const
SparseVanka(const SparseMatrix< number > &M, const std::vector< bool > &selected, const bool conserve_memory=false, const unsigned int n_threads=multithread_info.n_threads())
void compute_inverse(const size_type row, std::vector< size_type > &local_indices)
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
const_iterator end() const
::ExceptionBase & ExcInternalError()
void vmult(Vector< number2 > &dst, const Vector< number2 > &src) const
std::size_t memory_consumption() const