Reference documentation for deal.II version 8.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
chunk_sparse_matrix.templates.h
1 // ---------------------------------------------------------------------
2 // @f$Id: chunk_sparse_matrix.templates.h 31932 2013-12-08 02:15:54Z heister @f$
3 //
4 // Copyright (C) 2008 - 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__chunk_sparse_matrix_templates_h
18 #define __deal2__chunk_sparse_matrix_templates_h
19 
20 
21 #include <deal.II/base/template_constraints.h>
22 #include <deal.II/base/parallel.h>
23 #include <deal.II/lac/chunk_sparse_matrix.h>
24 #include <deal.II/lac/vector.h>
25 #include <deal.II/lac/full_matrix.h>
26 
27 
28 #include <ostream>
29 #include <iomanip>
30 #include <algorithm>
31 #include <functional>
32 #include <cmath>
33 
34 #include <vector>
35 #include <numeric>
36 
37 #include <deal.II/base/thread_management.h>
38 #include <deal.II/base/multithread_info.h>
39 
40 DEAL_II_NAMESPACE_OPEN
41 
42 
43 namespace internal
44 {
45 //TODO: the goal of the ChunkSparseMatrix class is to stream data and use
46 // the vectorization features of modern processors. to make this happen,
47 // we will have to vectorize the functions in the following namespace, either
48 // by hand or by using, for example, optimized BLAS versions for them.
49  namespace ChunkSparseMatrix
50  {
55 
63  template <typename MatrixIterator,
64  typename SrcIterator,
65  typename DstIterator>
66  inline
67  void
68  chunk_vmult_add (const size_type chunk_size,
69  const MatrixIterator matrix,
70  const SrcIterator src,
71  DstIterator dst)
72  {
73  MatrixIterator matrix_row = matrix;
74 
75  for (size_type i=0; i<chunk_size;
76  ++i, matrix_row += chunk_size)
77  {
78  typename std::iterator_traits<DstIterator>::value_type
79  sum = 0;
80 
81  for (size_type j=0; j<chunk_size; ++j)
82  sum += matrix_row[j] * src[j];
83 
84  dst[i] += sum;
85  }
86  }
87 
88 
89 
94  template <typename MatrixIterator,
95  typename SrcIterator,
96  typename DstIterator>
97  inline
98  void
99  chunk_vmult_subtract (const size_type chunk_size,
100  const MatrixIterator matrix,
101  const SrcIterator src,
102  DstIterator dst)
103  {
104  MatrixIterator matrix_row = matrix;
105 
106  for (size_type i=0; i<chunk_size;
107  ++i, matrix_row += chunk_size)
108  {
109  typename std::iterator_traits<DstIterator>::value_type
110  sum = 0;
111 
112  for (size_type j=0; j<chunk_size; ++j)
113  sum += matrix_row[j] * src[j];
114 
115  dst[i] -= sum;
116  }
117  }
118 
119 
125  template <typename MatrixIterator,
126  typename SrcIterator,
127  typename DstIterator>
128  inline
129  void
130  chunk_Tvmult_add (const size_type chunk_size,
131  const MatrixIterator matrix,
132  const SrcIterator src,
133  DstIterator dst)
134  {
135  for (size_type i=0; i<chunk_size; ++i)
136  {
137  typename std::iterator_traits<DstIterator>::value_type
138  sum = 0;
139 
140  for (size_type j=0; j<chunk_size; ++j)
141  sum += matrix[j*chunk_size+i] * src[j];
142 
143  dst[i] += sum;
144  }
145  }
146 
147 
152  template <typename result_type,
153  typename MatrixIterator,
154  typename SrcIterator1,
155  typename SrcIterator2>
156  inline
157  result_type
158  chunk_matrix_scalar_product (const size_type chunk_size,
159  const MatrixIterator matrix,
160  const SrcIterator1 u,
161  const SrcIterator2 v)
162  {
163  result_type result = 0;
164 
165  MatrixIterator matrix_row = matrix;
166 
167  for (size_type i=0; i<chunk_size;
168  ++i, matrix_row += chunk_size)
169  {
170  typename std::iterator_traits<SrcIterator2>::value_type
171  sum = 0;
172 
173  for (size_type j=0; j<chunk_size; ++j)
174  sum += matrix_row[j] * v[j];
175 
176  result += u[i] * sum;
177  }
178 
179  return result;
180  }
181 
182 
183 
192  template <typename number,
193  typename InVector,
194  typename OutVector>
195  void vmult_add_on_subrange (const ChunkSparsityPattern &cols,
196  const unsigned int begin_row,
197  const unsigned int end_row,
198  const number *values,
199  const std::size_t *rowstart,
200  const size_type *colnums,
201  const InVector &src,
202  OutVector &dst)
203  {
204  const size_type m = cols.n_rows();
205  const size_type n = cols.n_cols();
206  const size_type chunk_size = cols.get_chunk_size();
207 
208  // loop over all chunks. note that we need to treat the last chunk row
209  // and column differently if they have padding elements
210  const size_type n_filled_last_rows = m % chunk_size;
211  const size_type n_filled_last_cols = n % chunk_size;
212 
213  const size_type last_regular_row = n_filled_last_rows > 0 ?
214  std::min(m/chunk_size,
215  static_cast<size_type>(end_row)) :
216  end_row;
217  const size_type irregular_col = n/chunk_size;
218 
219  typename OutVector::iterator dst_ptr = dst.begin()+chunk_size*begin_row;
220  const number *val_ptr= &values[rowstart[begin_row]*chunk_size*chunk_size];
221  const size_type *colnum_ptr = &colnums[rowstart[begin_row]];
222  for (unsigned int chunk_row=begin_row; chunk_row<last_regular_row;
223  ++chunk_row)
224  {
225  const number *const val_end_of_row = &values[rowstart[chunk_row+1] *
226  chunk_size * chunk_size];
227  while (val_ptr != val_end_of_row)
228  {
229  if (*colnum_ptr != irregular_col)
230  chunk_vmult_add (chunk_size,
231  val_ptr,
232  src.begin() + *colnum_ptr * chunk_size,
233  dst_ptr);
234  else
235  // we're at a chunk column that has padding
236  for (size_type r=0; r<chunk_size; ++r)
237  for (size_type c=0; c<n_filled_last_cols; ++c)
238  dst_ptr[r] += (val_ptr[r*chunk_size + c] *
239  src(*colnum_ptr * chunk_size + c));
240 
241  ++colnum_ptr;
242  val_ptr += chunk_size * chunk_size;
243  }
244 
245  dst_ptr += chunk_size;
246  }
247 
248  // now deal with last chunk row if necessary
249  if (n_filled_last_rows > 0 && end_row == (m/chunk_size+1))
250  {
251  const size_type chunk_row = last_regular_row;
252 
253  const number *const val_end_of_row = &values[rowstart[chunk_row+1] *
254  chunk_size * chunk_size];
255  while (val_ptr != val_end_of_row)
256  {
257  if (*colnum_ptr != irregular_col)
258  {
259  // we're at a chunk row but not column that has padding
260  for (size_type r=0; r<n_filled_last_rows; ++r)
261  for (size_type c=0; c<chunk_size; ++c)
262  dst_ptr[r]
263  += (val_ptr[r*chunk_size + c] *
264  src(*colnum_ptr * chunk_size + c));
265  }
266  else
267  // we're at a chunk row and column that has padding
268  for (size_type r=0; r<n_filled_last_rows; ++r)
269  for (size_type c=0; c<n_filled_last_cols; ++c)
270  dst_ptr[r]
271  += (val_ptr[r*chunk_size + c] *
272  src(*colnum_ptr * chunk_size + c));
273 
274  ++colnum_ptr;
275  val_ptr += chunk_size * chunk_size;
276  }
277  }
278  Assert(std::size_t(colnum_ptr-&colnums[0]) == rowstart[end_row],
279  ExcInternalError());
280  Assert(std::size_t(val_ptr-&values[0]) ==
281  rowstart[end_row] * chunk_size * chunk_size,
282  ExcInternalError());
283  }
284  }
285 }
286 
287 
288 
289 template <typename number>
291  :
292  cols(0, "ChunkSparseMatrix"),
293  val(0),
294  max_len(0)
295 {}
296 
297 
298 
299 template <typename number>
301  :
302  Subscriptor (m),
303  cols(0, "ChunkSparseMatrix"),
304  val(0),
305  max_len(0)
306 {
310 }
311 
312 
313 
314 template <typename number>
317 {
321 
322  return *this;
323 }
324 
325 
326 
327 template <typename number>
329  :
330  cols(0, "ChunkSparseMatrix"),
331  val(0),
332  max_len(0)
333 {
334  reinit (c);
335 }
336 
337 
338 
339 template <typename number>
341  const IdentityMatrix &id)
342  :
343  cols(0, "ChunkSparseMatrix"),
344  val(0),
345  max_len(0)
346 {
347  Assert (c.n_rows() == id.m(), ExcDimensionMismatch (c.n_rows(), id.m()));
348  Assert (c.n_cols() == id.n(), ExcDimensionMismatch (c.n_cols(), id.n()));
349 
350  reinit (c);
351  for (size_type i=0; i<n(); ++i)
352  this->set(i,i,1.);
353 }
354 
355 
356 
357 template <typename number>
359 {
360  cols = 0;
361 
362  if (val != 0)
363  delete[] val;
364 }
365 
366 
367 
368 namespace internal
369 {
370  namespace ChunkSparseMatrix
371  {
372  template<typename T>
373  void zero_subrange (const unsigned int begin,
374  const unsigned int end,
375  T *dst)
376  {
377  std::memset (dst+begin,0,(end-begin)*sizeof(T));
378  }
379  }
380 }
381 
382 
383 
384 template <typename number>
387 {
389 
390  Assert (cols != 0, ExcNotInitialized());
391  Assert (cols->sparsity_pattern.compressed || cols->empty(),
392  ChunkSparsityPattern::ExcNotCompressed());
393 
394  // do initial zeroing of elements in parallel. Try to achieve a similar
395  // layout as when doing matrix-vector products, as on some NUMA systems, a
396  // memory block is assigned to memory banks where the first access is
397  // generated. For sparse matrices, the first operations is usually the
398  // operator=. The grain size is chosen to reflect the number of rows in
399  // minimum_parallel_grain_size, weighted by the number of nonzero entries
400  // per row on average.
401  const unsigned int matrix_size = cols->sparsity_pattern.n_nonzero_elements()
402  * cols->chunk_size * cols->chunk_size;
403  const unsigned int grain_size =
404  internal::SparseMatrix::minimum_parallel_grain_size *
405  (matrix_size+m()) / m();
406  if (matrix_size>grain_size)
407  parallel::apply_to_subranges (0U, matrix_size,
408  std_cxx1x::bind(&internal::ChunkSparseMatrix::template
409  zero_subrange<number>,
410  std_cxx1x::_1, std_cxx1x::_2,
411  val),
412  grain_size);
413  else if (matrix_size > 0)
414  std::memset (&val[0], 0, matrix_size*sizeof(number));
415 
416  return *this;
417 }
418 
419 
420 
421 template <typename number>
424 {
425  Assert (cols->n_rows() == id.m(),
426  ExcDimensionMismatch (cols->n_rows(), id.m()));
427  Assert (cols->n_cols() == id.n(),
428  ExcDimensionMismatch (cols->n_cols(), id.n()));
429 
430  *this = 0;
431  for (size_type i=0; i<n(); ++i)
432  this->set(i,i,1.);
433 
434  return *this;
435 }
436 
437 
438 
439 template <typename number>
440 void
442 {
443  cols = &sparsity;
444 
445  if (cols->empty())
446  {
447  if (val != 0)
448  delete[] val;
449  val = 0;
450  max_len = 0;
451  return;
452  }
453 
454  // allocate not just m() * n() elements but enough so that we can store full
455  // chunks. this entails some padding elements
456  const size_type chunk_size = cols->get_chunk_size();
457  const size_type N = cols->sparsity_pattern.n_nonzero_elements() *
458  chunk_size * chunk_size;
459  if (N > max_len || max_len == 0)
460  {
461  if (val != 0)
462  delete[] val;
463  val = new number[N];
464  max_len = N;
465  }
466 
467  // fill with zeros. do not just fill N elements but all that we allocated to
468  // ensure that also the padding elements are zero and not left at previous
469  // values
470  this->operator=(0.);
471 }
472 
473 
474 
475 template <typename number>
476 void
478 {
479  cols = 0;
480  if (val) delete[] val;
481  val = 0;
482  max_len = 0;
483 }
484 
485 
486 
487 template <typename number>
488 bool
490 {
491  if (cols == 0)
492  return true;
493  else
494  return cols->empty();
495 }
496 
497 
498 
499 template <typename number>
502 {
503  Assert (cols != 0, ExcNotInitialized());
504  return cols->n_nonzero_elements ();
505 }
506 
507 
508 
509 template <typename number>
512 {
513  Assert (cols != 0, ExcNotInitialized());
514 
515  // count those elements that are nonzero, even if they lie in the padding
516  // around the matrix. since we have the invariant that padding elements are
517  // zero, nothing bad can happen here
518  const size_type chunk_size = cols->get_chunk_size();
519  return std::count_if(&val[0],
520  &val[cols->sparsity_pattern.n_nonzero_elements () *
521  chunk_size * chunk_size],
522  std::bind2nd(std::not_equal_to<double>(), 0));
523 }
524 
525 
526 
527 template <typename number>
528 void
530 {
531  Assert (cols != 0, ExcNotInitialized());
532  Assert (cols->rows == cols->cols, ExcNotQuadratic());
533 
534  Assert (false, ExcNotImplemented());
535 }
536 
537 
538 
539 template <typename number>
540 template <typename somenumber>
543 {
544  Assert (cols != 0, ExcNotInitialized());
545  Assert (val != 0, ExcNotInitialized());
546  Assert (cols == matrix.cols, ExcDifferentChunkSparsityPatterns());
547 
548  // copy everything, including padding elements
549  const size_type chunk_size = cols->get_chunk_size();
550  std::copy (&matrix.val[0],
551  &matrix.val[cols->sparsity_pattern.n_nonzero_elements()
552  * chunk_size * chunk_size],
553  &val[0]);
554 
555  return *this;
556 }
557 
558 
559 
560 template <typename number>
561 template <typename somenumber>
562 void
564 {
565  // first delete previous content
566  *this = 0;
567 
568  // then copy old matrix
569  for (size_type row=0; row<matrix.m(); ++row)
570  for (size_type col=0; col<matrix.n(); ++col)
571  if (matrix(row,col) != 0)
572  set (row, col, matrix(row,col));
573 }
574 
575 
576 
577 template <typename number>
578 template <typename somenumber>
579 void
580 ChunkSparseMatrix<number>::add (const number factor,
581  const ChunkSparseMatrix<somenumber> &matrix)
582 {
583  Assert (cols != 0, ExcNotInitialized());
584  Assert (val != 0, ExcNotInitialized());
585  Assert (cols == matrix.cols, ExcDifferentChunkSparsityPatterns());
586 
587  // add everything, including padding elements
588  const size_type chunk_size = cols->get_chunk_size();
589  number *val_ptr = &val[0];
590  const somenumber *matrix_ptr = &matrix.val[0];
591  const number *const end_ptr = &val[cols->sparsity_pattern.n_nonzero_elements()
592  * chunk_size * chunk_size];
593 
594  while (val_ptr != end_ptr)
595  *val_ptr++ += factor **matrix_ptr++;
596 }
597 
598 
599 template <typename number>
600 template <class OutVector, class InVector>
601 void
603  const InVector &src) const
604 {
605  Assert (cols != 0, ExcNotInitialized());
606  Assert (val != 0, ExcNotInitialized());
607  Assert(m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
608  Assert(n() == src.size(), ExcDimensionMismatch(n(),src.size()));
609 
610  Assert (!PointerComparison::equal(&src, &dst), ExcSourceEqualsDestination());
611 
612  // set the output vector to zero and then add to it the contributions of
613  // vmults from individual chunks. this is what vmult_add does
614  dst = 0;
615  vmult_add (dst, src);
616 }
617 
618 
619 
620 template <typename number>
621 template <class OutVector, class InVector>
622 void
624  const InVector &src) const
625 {
626  Assert (val != 0, ExcNotInitialized());
627  Assert (cols != 0, ExcNotInitialized());
628  Assert(n() == dst.size(), ExcDimensionMismatch(n(),dst.size()));
629  Assert(m() == src.size(), ExcDimensionMismatch(m(),src.size()));
630 
631  Assert (!PointerComparison::equal(&src, &dst), ExcSourceEqualsDestination());
632 
633  Assert (cols != 0, ExcNotInitialized());
634  Assert (val != 0, ExcNotInitialized());
635  Assert(m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
636  Assert(n() == src.size(), ExcDimensionMismatch(n(),src.size()));
637 
638  Assert (!PointerComparison::equal(&src, &dst), ExcSourceEqualsDestination());
639 
640  // set the output vector to zero and then add to it the contributions of
641  // vmults from individual chunks. this is what vmult_add does
642  dst = 0;
643  Tvmult_add (dst, src);
644 }
645 
646 
647 
648 template <typename number>
649 template <class OutVector, class InVector>
650 void
652  const InVector &src) const
653 {
654  Assert (cols != 0, ExcNotInitialized());
655  Assert (val != 0, ExcNotInitialized());
656  Assert(m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
657  Assert(n() == src.size(), ExcDimensionMismatch(n(),src.size()));
658 
659  Assert (!PointerComparison::equal(&src, &dst), ExcSourceEqualsDestination());
661  std_cxx1x::bind (&internal::ChunkSparseMatrix::vmult_add_on_subrange
662  <number,InVector,OutVector>,
663  std_cxx1x::cref(*cols),
664  std_cxx1x::_1, std_cxx1x::_2,
665  val,
668  std_cxx1x::cref(src),
669  std_cxx1x::ref(dst)),
670  internal::SparseMatrix::minimum_parallel_grain_size/cols->chunk_size+1);
671 
672 }
673 
674 
675 template <typename number>
676 template <class OutVector, class InVector>
677 void
679  const InVector &src) const
680 {
681  Assert (cols != 0, ExcNotInitialized());
682  Assert (val != 0, ExcNotInitialized());
683  Assert(m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
684  Assert(n() == src.size(), ExcDimensionMismatch(n(),src.size()));
685 
686  Assert (!PointerComparison::equal(&src, &dst), ExcSourceEqualsDestination());
687 
688  const size_type n_chunk_rows = cols->sparsity_pattern.n_rows();
689 
690  // loop over all chunks. note that we need to treat the last chunk row and
691  // column differently if they have padding elements
692  const bool rows_have_padding = (m() % cols->chunk_size != 0),
693  cols_have_padding = (n() % cols->chunk_size != 0);
694 
695  const size_type n_regular_chunk_rows
696  = (rows_have_padding ?
697  n_chunk_rows-1 :
698  n_chunk_rows);
699 
700  // like in vmult_add, but don't keep an iterator into dst around since we're
701  // not traversing it sequentially this time
702  const number *val_ptr = val;
703  const size_type *colnum_ptr = cols->sparsity_pattern.colnums;
704 
705  for (size_type chunk_row=0; chunk_row<n_regular_chunk_rows; ++chunk_row)
706  {
707  const number *const val_end_of_row = &val[cols->sparsity_pattern.rowstart[chunk_row+1]
708  * cols->chunk_size
709  * cols->chunk_size];
710  while (val_ptr != val_end_of_row)
711  {
712  if ((cols_have_padding == false)
713  ||
714  (*colnum_ptr != cols->sparsity_pattern.n_cols()-1))
715  internal::ChunkSparseMatrix::chunk_Tvmult_add
716  (cols->chunk_size,
717  val_ptr,
718  src.begin() + chunk_row * cols->chunk_size,
719  dst.begin() + *colnum_ptr * cols->chunk_size);
720  else
721  // we're at a chunk column that has padding
722  for (size_type r=0; r<cols->chunk_size; ++r)
723  for (size_type c=0; c<n() % cols->chunk_size; ++c)
724  dst(*colnum_ptr * cols->chunk_size + c)
725  += (val_ptr[r*cols->chunk_size + c] *
726  src(chunk_row * cols->chunk_size + r));
727 
728  ++colnum_ptr;
729  val_ptr += cols->chunk_size * cols->chunk_size;
730  }
731  }
732 
733  // now deal with last chunk row if necessary
734  if (rows_have_padding)
735  {
736  const size_type chunk_row = n_chunk_rows - 1;
737 
738  const number *const val_end_of_row = &val[cols->sparsity_pattern.rowstart[chunk_row+1]
739  * cols->chunk_size
740  * cols->chunk_size];
741  while (val_ptr != val_end_of_row)
742  {
743  if ((cols_have_padding == false)
744  ||
745  (*colnum_ptr != cols->sparsity_pattern.n_cols()-1))
746  {
747  // we're at a chunk row but not column that has padding
748  for (size_type r=0; r<m() % cols->chunk_size; ++r)
749  for (size_type c=0; c<cols->chunk_size; ++c)
750  dst(*colnum_ptr * cols->chunk_size + c)
751  += (val_ptr[r*cols->chunk_size + c] *
752  src(chunk_row * cols->chunk_size + r));
753  }
754  else
755  // we're at a chunk row and column that has padding
756  for (size_type r=0; r<m() % cols->chunk_size; ++r)
757  for (size_type c=0; c<n() % cols->chunk_size; ++c)
758  dst(*colnum_ptr * cols->chunk_size + c)
759  += (val_ptr[r*cols->chunk_size + c] *
760  src(chunk_row * cols->chunk_size + r));
761 
762  ++colnum_ptr;
763  val_ptr += cols->chunk_size * cols->chunk_size;
764  }
765  }
766 }
767 
768 
769 template <typename number>
770 template <typename somenumber>
771 somenumber
773 {
774  Assert (cols != 0, ExcNotInitialized());
775  Assert (val != 0, ExcNotInitialized());
776  Assert(m() == v.size(), ExcDimensionMismatch(m(),v.size()));
777  Assert(n() == v.size(), ExcDimensionMismatch(n(),v.size()));
778 
779  somenumber result = 0;
780 
782  // like matrix_scalar_product, except that the two vectors are now the same
783 
784  const size_type n_chunk_rows = cols->sparsity_pattern.n_rows();
785 
786  // loop over all chunks. note that we need to treat the last chunk row and
787  // column differently if they have padding elements
788  const bool rows_have_padding = (m() % cols->chunk_size != 0),
789  cols_have_padding = (n() % cols->chunk_size != 0);
790 
791  const size_type n_regular_chunk_rows
792  = (rows_have_padding ?
793  n_chunk_rows-1 :
794  n_chunk_rows);
795 
796  const number *val_ptr = val;
797  const size_type *colnum_ptr = cols->sparsity_pattern.colnums;
798  typename Vector<somenumber>::const_iterator v_ptr = v.begin();
799 
800  for (size_type chunk_row=0; chunk_row<n_regular_chunk_rows; ++chunk_row)
801  {
802  const number *const val_end_of_row = &val[cols->sparsity_pattern.rowstart[chunk_row+1]
803  * cols->chunk_size
804  * cols->chunk_size];
805  while (val_ptr != val_end_of_row)
806  {
807  if ((cols_have_padding == false)
808  ||
809  (*colnum_ptr != cols->sparsity_pattern.n_cols()-1))
810  result +=
811  internal::ChunkSparseMatrix::
812  chunk_matrix_scalar_product<somenumber>
813  (cols->chunk_size,
814  val_ptr,
815  v_ptr,
816  v.begin() + *colnum_ptr * cols->chunk_size);
817  else
818  // we're at a chunk column that has padding
819  for (size_type r=0; r<cols->chunk_size; ++r)
820  for (size_type c=0; c<n() % cols->chunk_size; ++c)
821  result
822  +=
823  v(chunk_row * cols->chunk_size + r)
824  * (val_ptr[r*cols->chunk_size + c] *
825  v(*colnum_ptr * cols->chunk_size + c));
826 
827  ++colnum_ptr;
828  val_ptr += cols->chunk_size * cols->chunk_size;
829  }
830 
831 
832  v_ptr += cols->chunk_size;
833  }
834 
835  // now deal with last chunk row if necessary
836  if (rows_have_padding)
837  {
838  const size_type chunk_row = n_chunk_rows - 1;
839 
840  const number *const val_end_of_row = &val[cols->sparsity_pattern.rowstart[chunk_row+1]
841  * cols->chunk_size
842  * cols->chunk_size];
843  while (val_ptr != val_end_of_row)
844  {
845  if ((cols_have_padding == false)
846  ||
847  (*colnum_ptr != cols->sparsity_pattern.n_cols()-1))
848  {
849  // we're at a chunk row but not column that has padding
850  for (size_type r=0; r<m() % cols->chunk_size; ++r)
851  for (size_type c=0; c<cols->chunk_size; ++c)
852  result
853  +=
854  v(chunk_row * cols->chunk_size + r)
855  * (val_ptr[r*cols->chunk_size + c] *
856  v(*colnum_ptr * cols->chunk_size + c));
857  }
858  else
859  // we're at a chunk row and column that has padding
860  for (size_type r=0; r<m() % cols->chunk_size; ++r)
861  for (size_type c=0; c<n() % cols->chunk_size; ++c)
862  result
863  +=
864  v(chunk_row * cols->chunk_size + r)
865  * (val_ptr[r*cols->chunk_size + c] *
866  v(*colnum_ptr * cols->chunk_size + c));
867 
868  ++colnum_ptr;
869  val_ptr += cols->chunk_size * cols->chunk_size;
870  }
871  }
872 
873  return result;
874 }
875 
876 
877 
878 template <typename number>
879 template <typename somenumber>
880 somenumber
882  const Vector<somenumber> &v) const
883 {
884  Assert (cols != 0, ExcNotInitialized());
885  Assert (val != 0, ExcNotInitialized());
886  Assert(m() == u.size(), ExcDimensionMismatch(m(),u.size()));
887  Assert(n() == v.size(), ExcDimensionMismatch(n(),v.size()));
888 
889  // the following works like the vmult_add function
890  somenumber result = 0;
891 
892  const size_type n_chunk_rows = cols->sparsity_pattern.n_rows();
893 
894  // loop over all chunks. note that we need to treat the last chunk row and
895  // column differently if they have padding elements
896  const bool rows_have_padding = (m() % cols->chunk_size != 0),
897  cols_have_padding = (n() % cols->chunk_size != 0);
898 
899  const size_type n_regular_chunk_rows
900  = (rows_have_padding ?
901  n_chunk_rows-1 :
902  n_chunk_rows);
903 
904  const number *val_ptr = val;
905  const size_type *colnum_ptr = cols->sparsity_pattern.colnums;
906  typename Vector<somenumber>::const_iterator u_ptr = u.begin();
907 
908  for (size_type chunk_row=0; chunk_row<n_regular_chunk_rows; ++chunk_row)
909  {
910  const number *const val_end_of_row = &val[cols->sparsity_pattern.rowstart[chunk_row+1]
911  * cols->chunk_size
912  * cols->chunk_size];
913  while (val_ptr != val_end_of_row)
914  {
915  if ((cols_have_padding == false)
916  ||
917  (*colnum_ptr != cols->sparsity_pattern.n_cols()-1))
918  result +=
919  internal::ChunkSparseMatrix::
920  chunk_matrix_scalar_product<somenumber>
921  (cols->chunk_size,
922  val_ptr,
923  u_ptr,
924  v.begin() + *colnum_ptr * cols->chunk_size);
925  else
926  // we're at a chunk column that has padding
927  for (size_type r=0; r<cols->chunk_size; ++r)
928  for (size_type c=0; c<n() % cols->chunk_size; ++c)
929  result
930  +=
931  u(chunk_row * cols->chunk_size + r)
932  * (val_ptr[r*cols->chunk_size + c] *
933  v(*colnum_ptr * cols->chunk_size + c));
934 
935  ++colnum_ptr;
936  val_ptr += cols->chunk_size * cols->chunk_size;
937  }
938 
939 
940  u_ptr += cols->chunk_size;
941  }
942 
943  // now deal with last chunk row if necessary
944  if (rows_have_padding)
945  {
946  const size_type chunk_row = n_chunk_rows - 1;
947 
948  const number *const val_end_of_row = &val[cols->sparsity_pattern.rowstart[chunk_row+1]
949  * cols->chunk_size
950  * cols->chunk_size];
951  while (val_ptr != val_end_of_row)
952  {
953  if ((cols_have_padding == false)
954  ||
955  (*colnum_ptr != cols->sparsity_pattern.n_cols()-1))
956  {
957  // we're at a chunk row but not column that has padding
958  for (size_type r=0; r<m() % cols->chunk_size; ++r)
959  for (size_type c=0; c<cols->chunk_size; ++c)
960  result
961  +=
962  u(chunk_row * cols->chunk_size + r)
963  * (val_ptr[r*cols->chunk_size + c] *
964  v(*colnum_ptr * cols->chunk_size + c));
965  }
966  else
967  // we're at a chunk row and column that has padding
968  for (size_type r=0; r<m() % cols->chunk_size; ++r)
969  for (size_type c=0; c<n() % cols->chunk_size; ++c)
970  result
971  +=
972  u(chunk_row * cols->chunk_size + r)
973  * (val_ptr[r*cols->chunk_size + c] *
974  v(*colnum_ptr * cols->chunk_size + c));
975 
976  ++colnum_ptr;
977  val_ptr += cols->chunk_size * cols->chunk_size;
978  }
979  }
980 
981  return result;
982 }
983 
984 
985 
986 template <typename number>
989 {
990  Assert (cols != 0, ExcNotInitialized());
991  Assert (val != 0, ExcNotInitialized());
992 
993  const size_type n_chunk_rows = cols->sparsity_pattern.n_rows();
994 
995  // loop over all rows and columns; it is safe to also loop over the padding
996  // elements (they are zero) if we make sure that the vector into which we
997  // sum column sums is large enough
998  Vector<real_type> column_sums(cols->sparsity_pattern.n_cols() *
999  cols->chunk_size);
1000 
1001  for (size_type chunk_row=0; chunk_row<n_chunk_rows; ++chunk_row)
1002  for (size_type j=cols->sparsity_pattern.rowstart[chunk_row];
1003  j<cols->sparsity_pattern.rowstart[chunk_row+1] ; ++j)
1004  for (size_type r=0; r<cols->chunk_size; ++r)
1005  for (size_type s=0; s<cols->chunk_size; ++s)
1006  column_sums(cols->sparsity_pattern.colnums[j] *
1007  cols->chunk_size + s) +=
1009  cols->chunk_size +
1010  r * cols->chunk_size +
1011  s]);
1012 
1013  return column_sums.linfty_norm();
1014 }
1015 
1016 
1017 
1018 template <typename number>
1021 {
1022  Assert (cols != 0, ExcNotInitialized());
1023  Assert (val != 0, ExcNotInitialized());
1024 
1025  // this function works like l1_norm(). it can be made more efficient
1026  // (without allocating a temporary vector) as is done in the SparseMatrix
1027  // class but since it is rarely called in time critical places it is
1028  // probably not worth it
1029  const size_type n_chunk_rows = cols->sparsity_pattern.n_rows();
1030 
1031  // loop over all rows and columns; it is safe to also loop over the padding
1032  // elements (they are zero) if we make sure that the vector into which we
1033  // sum column sums is large enough
1034  Vector<real_type> row_sums(cols->sparsity_pattern.n_rows() *
1035  cols->chunk_size);
1036 
1037  for (size_type chunk_row=0; chunk_row<n_chunk_rows; ++chunk_row)
1038  for (size_type j=cols->sparsity_pattern.rowstart[chunk_row];
1039  j<cols->sparsity_pattern.rowstart[chunk_row+1] ; ++j)
1040  for (size_type r=0; r<cols->chunk_size; ++r)
1041  for (size_type s=0; s<cols->chunk_size; ++s)
1042  row_sums(chunk_row * cols->chunk_size + r) +=
1044  cols->chunk_size +
1045  r * cols->chunk_size +
1046  s]);
1047 
1048  return row_sums.linfty_norm();
1049 }
1050 
1051 
1052 
1053 template <typename number>
1056 {
1057  // simply add up all entries in the sparsity pattern, without taking any
1058  // reference to rows or columns
1059  //
1060  // padding elements are zero, so we can add them up as well
1061  real_type norm_sqr = 0;
1062  for (const number *ptr = &val[0]; ptr != &val[max_len]; ++ptr)
1064 
1065  return std::sqrt (norm_sqr);
1066 }
1067 
1068 
1069 
1070 template <typename number>
1071 template <typename somenumber>
1072 somenumber
1074  const Vector<somenumber> &u,
1075  const Vector<somenumber> &b) const
1076 {
1077  Assert (cols != 0, ExcNotInitialized());
1078  Assert (val != 0, ExcNotInitialized());
1079  Assert(m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
1080  Assert(m() == b.size(), ExcDimensionMismatch(m(),b.size()));
1081  Assert(n() == u.size(), ExcDimensionMismatch(n(),u.size()));
1082 
1083  Assert (&u != &dst, ExcSourceEqualsDestination());
1084 
1085  // set dst=b, then subtract the result of A*u from it. since the purpose of
1086  // the current class is to promote streaming of data rather than more random
1087  // access patterns, breaking things up into two loops may be reasonable
1088  dst = b;
1089 
1091  // the rest of this function is like vmult_add, except that we subtract
1092  // rather than add A*u
1094  const size_type n_chunk_rows = cols->sparsity_pattern.n_rows();
1095 
1096  // loop over all chunks. note that we need to treat the last chunk row and
1097  // column differently if they have padding elements
1098  const bool rows_have_padding = (m() % cols->chunk_size != 0),
1099  cols_have_padding = (n() % cols->chunk_size != 0);
1100 
1101  const size_type n_regular_chunk_rows
1102  = (rows_have_padding ?
1103  n_chunk_rows-1 :
1104  n_chunk_rows);
1105 
1106  const number *val_ptr = val;
1107  const size_type *colnum_ptr = cols->sparsity_pattern.colnums;
1108  typename Vector<somenumber>::iterator dst_ptr = dst.begin();
1109 
1110  for (size_type chunk_row=0; chunk_row<n_regular_chunk_rows; ++chunk_row)
1111  {
1112  const number *const val_end_of_row = &val[cols->sparsity_pattern.rowstart[chunk_row+1]
1113  * cols->chunk_size
1114  * cols->chunk_size];
1115  while (val_ptr != val_end_of_row)
1116  {
1117  if ((cols_have_padding == false)
1118  ||
1119  (*colnum_ptr != cols->sparsity_pattern.n_cols()-1))
1120  internal::ChunkSparseMatrix::chunk_vmult_subtract
1121  (cols->chunk_size,
1122  val_ptr,
1123  u.begin() + *colnum_ptr * cols->chunk_size,
1124  dst_ptr);
1125  else
1126  // we're at a chunk column that has padding
1127  for (size_type r=0; r<cols->chunk_size; ++r)
1128  for (size_type c=0; c<n() % cols->chunk_size; ++c)
1129  dst(chunk_row * cols->chunk_size + r)
1130  -= (val_ptr[r*cols->chunk_size + c] *
1131  u(*colnum_ptr * cols->chunk_size + c));
1132 
1133  ++colnum_ptr;
1134  val_ptr += cols->chunk_size * cols->chunk_size;
1135  }
1136 
1137 
1138  dst_ptr += cols->chunk_size;
1139  }
1140 
1141  // now deal with last chunk row if necessary
1142  if (rows_have_padding)
1143  {
1144  const size_type chunk_row = n_chunk_rows - 1;
1145 
1146  const number *const val_end_of_row = &val[cols->sparsity_pattern.rowstart[chunk_row+1]
1147  * cols->chunk_size
1148  * cols->chunk_size];
1149  while (val_ptr != val_end_of_row)
1150  {
1151  if ((cols_have_padding == false)
1152  ||
1153  (*colnum_ptr != cols->sparsity_pattern.n_cols()-1))
1154  {
1155  // we're at a chunk row but not column that has padding
1156  for (size_type r=0; r<m() % cols->chunk_size; ++r)
1157  for (size_type c=0; c<cols->chunk_size; ++c)
1158  dst(chunk_row * cols->chunk_size + r)
1159  -= (val_ptr[r*cols->chunk_size + c] *
1160  u(*colnum_ptr * cols->chunk_size + c));
1161  }
1162  else
1163  // we're at a chunk row and column that has padding
1164  for (size_type r=0; r<m() % cols->chunk_size; ++r)
1165  for (size_type c=0; c<n() % cols->chunk_size; ++c)
1166  dst(chunk_row * cols->chunk_size + r)
1167  -= (val_ptr[r*cols->chunk_size + c] *
1168  u(*colnum_ptr * cols->chunk_size + c));
1169 
1170  ++colnum_ptr;
1171  val_ptr += cols->chunk_size * cols->chunk_size;
1172  }
1173 
1174 
1175  dst_ptr += cols->chunk_size;
1176  }
1177 
1178  // finally compute the norm
1179  return dst.l2_norm();
1180 }
1181 
1182 
1183 
1184 template <typename number>
1185 template <typename somenumber>
1186 void
1188  const Vector<somenumber> &src,
1189  const number /*om*/) const
1190 {
1191  Assert (cols != 0, ExcNotInitialized());
1192  Assert (val != 0, ExcNotInitialized());
1193  Assert (m() == n(), ExcMessage("This operation is only valid on square matrices."));
1194 
1195  Assert (dst.size() == n(), ExcDimensionMismatch (dst.size(), n()));
1196  Assert (src.size() == n(), ExcDimensionMismatch (src.size(), n()));
1197 
1198  Assert (false, ExcNotImplemented());
1199 }
1200 
1201 
1202 
1203 template <typename number>
1204 template <typename somenumber>
1205 void
1207  const Vector<somenumber> &src,
1208  const number /*om*/) const
1209 {
1210  // to understand how this function works you may want to take a look at the
1211  // CVS archives to see the original version which is much clearer...
1212  Assert (cols != 0, ExcNotInitialized());
1213  Assert (val != 0, ExcNotInitialized());
1214  Assert (m() == n(), ExcMessage("This operation is only valid on square matrices."));
1215 
1216  Assert (dst.size() == n(), ExcDimensionMismatch (dst.size(), n()));
1217  Assert (src.size() == n(), ExcDimensionMismatch (src.size(), n()));
1218 
1219  Assert (false, ExcNotImplemented());
1220 }
1221 
1222 
1223 template <typename number>
1224 template <typename somenumber>
1225 void
1227  const Vector<somenumber> &src,
1228  const number om) const
1229 {
1230  Assert (cols != 0, ExcNotInitialized());
1231  Assert (val != 0, ExcNotInitialized());
1232  Assert (m() == n(), ExcMessage("This operation is only valid on square matrices."));
1233 
1234 
1235  dst = src;
1236  SOR(dst,om);
1237 }
1238 
1239 
1240 template <typename number>
1241 template <typename somenumber>
1242 void
1244  const Vector<somenumber> &src,
1245  const number om) const
1246 {
1247  Assert (cols != 0, ExcNotInitialized());
1248  Assert (val != 0, ExcNotInitialized());
1249  Assert (m() == n(), ExcMessage("This operation is only valid on square matrices."));
1250 
1251 
1252  dst = src;
1253  TSOR(dst,om);
1254 }
1255 
1256 
1257 template <typename number>
1258 template <typename somenumber>
1259 void
1261  const number /*om*/) const
1262 {
1263  Assert (cols != 0, ExcNotInitialized());
1264  Assert (val != 0, ExcNotInitialized());
1265  Assert (m() == n(), ExcMessage("This operation is only valid on square matrices."));
1266  Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
1267 
1268  Assert (false, ExcNotImplemented());
1269 }
1270 
1271 
1272 template <typename number>
1273 template <typename somenumber>
1274 void
1276  const number /*om*/) const
1277 {
1278  Assert (cols != 0, ExcNotInitialized());
1279  Assert (val != 0, ExcNotInitialized());
1280  Assert (m() == n(), ExcMessage("This operation is only valid on square matrices."));
1281  Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
1282 
1283  Assert (false, ExcNotImplemented());
1284 }
1285 
1286 
1287 template <typename number>
1288 template <typename somenumber>
1289 void
1291  const std::vector<size_type> &permutation,
1292  const std::vector<size_type> &inverse_permutation,
1293  const number /*om*/) const
1294 {
1295  Assert (cols != 0, ExcNotInitialized());
1296  Assert (val != 0, ExcNotInitialized());
1297  Assert (m() == n(), ExcMessage("This operation is only valid on square matrices."));
1298 
1299  Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
1300  Assert (m() == permutation.size(),
1301  ExcDimensionMismatch(m(), permutation.size()));
1302  Assert (m() == inverse_permutation.size(),
1303  ExcDimensionMismatch(m(), inverse_permutation.size()));
1304 
1305  Assert (false, ExcNotImplemented());
1306 }
1307 
1308 
1309 template <typename number>
1310 template <typename somenumber>
1311 void
1313  const std::vector<size_type> &permutation,
1314  const std::vector<size_type> &inverse_permutation,
1315  const number /*om*/) const
1316 {
1317  Assert (cols != 0, ExcNotInitialized());
1318  Assert (val != 0, ExcNotInitialized());
1319  Assert (m() == n(), ExcMessage("This operation is only valid on square matrices."));
1320 
1321  Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
1322  Assert (m() == permutation.size(),
1323  ExcDimensionMismatch(m(), permutation.size()));
1324  Assert (m() == inverse_permutation.size(),
1325  ExcDimensionMismatch(m(), inverse_permutation.size()));
1326 
1327  Assert (false, ExcNotImplemented());
1328 }
1329 
1330 
1331 
1332 template <typename number>
1333 template <typename somenumber>
1334 void
1336  const Vector<somenumber> &b,
1337  const number /*om*/) const
1338 {
1339  Assert (cols != 0, ExcNotInitialized());
1340  Assert (val != 0, ExcNotInitialized());
1341  Assert (m() == n(), ExcMessage("This operation is only valid on square matrices."));
1342 
1343  Assert (m() == v.size(), ExcDimensionMismatch(m(),v.size()));
1344  Assert (m() == b.size(), ExcDimensionMismatch(m(),b.size()));
1345 
1346  Assert (false, ExcNotImplemented());
1347 }
1348 
1349 
1350 
1351 template <typename number>
1352 template <typename somenumber>
1353 void
1355  const Vector<somenumber> &b,
1356  const number /*om*/) const
1357 {
1358  Assert (cols != 0, ExcNotInitialized());
1359  Assert (val != 0, ExcNotInitialized());
1360  Assert (m() == n(), ExcMessage("This operation is only valid on square matrices."));
1361 
1362  Assert (m() == v.size(), ExcDimensionMismatch(m(),v.size()));
1363  Assert (m() == b.size(), ExcDimensionMismatch(m(),b.size()));
1364 
1365  Assert (false, ExcNotImplemented());
1366 }
1367 
1368 
1369 
1370 template <typename number>
1371 template <typename somenumber>
1372 void
1374  const Vector<somenumber> &b,
1375  const number om) const
1376 {
1377  SOR_step(v,b,om);
1378  TSOR_step(v,b,om);
1379 }
1380 
1381 
1382 
1383 template <typename number>
1384 template <typename somenumber>
1385 void
1387  const number /*om*/) const
1388 {
1389  Assert (cols != 0, ExcNotInitialized());
1390  Assert (val != 0, ExcNotInitialized());
1391  Assert (m() == n(), ExcMessage("This operation is only valid on square matrices."));
1392 
1393  Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
1394 
1395  Assert (false, ExcNotImplemented());
1396 }
1397 
1398 
1399 
1400 template <typename number>
1401 void ChunkSparseMatrix<number>::print (std::ostream &out) const
1402 {
1403  AssertThrow (out, ExcIO());
1404 
1405  Assert (cols != 0, ExcNotInitialized());
1406  Assert (val != 0, ExcNotInitialized());
1407 
1408  Assert (false, ExcNotImplemented());
1409 
1410  AssertThrow (out, ExcIO());
1411 }
1412 
1413 
1414 template <typename number>
1416  const unsigned int precision,
1417  const bool scientific,
1418  const unsigned int width_,
1419  const char *zero_string,
1420  const double denominator) const
1421 {
1422  AssertThrow (out, ExcIO());
1423 
1424  Assert (cols != 0, ExcNotInitialized());
1425  Assert (val != 0, ExcNotInitialized());
1426 
1427  unsigned int width = width_;
1428 
1429  Assert (false, ExcNotImplemented());
1430 
1431  std::ios::fmtflags old_flags = out.flags();
1432  unsigned int old_precision = out.precision (precision);
1433 
1434  if (scientific)
1435  {
1436  out.setf (std::ios::scientific, std::ios::floatfield);
1437  if (!width)
1438  width = precision+7;
1439  }
1440  else
1441  {
1442  out.setf (std::ios::fixed, std::ios::floatfield);
1443  if (!width)
1444  width = precision+2;
1445  }
1446 
1447  for (size_type i=0; i<m(); ++i)
1448  {
1449  for (size_type j=0; j<n(); ++j)
1451  out << std::setw(width)
1452  << val[cols->sparsity_pattern(i,j)] * denominator << ' ';
1453  else
1454  out << std::setw(width) << zero_string << ' ';
1455  out << std::endl;
1456  };
1457  AssertThrow (out, ExcIO());
1458 
1459  // reset output format
1460  out.precision(old_precision);
1461  out.flags (old_flags);
1462 }
1463 
1464 
1465 
1466 template <typename number>
1468  const double threshold) const
1469 {
1470  AssertThrow (out, ExcIO());
1471 
1472  Assert (cols != 0, ExcNotInitialized());
1473  Assert (val != 0, ExcNotInitialized());
1474 
1475  const size_type chunk_size = cols->get_chunk_size();
1476 
1477  // loop over all chunk rows and columns, and each time we find something
1478  // repeat it chunk_size times in both directions
1479  for (size_type i=0; i<cols->sparsity_pattern.n_rows(); ++i)
1480  {
1481  for (size_type d=0; d<chunk_size; ++d)
1482  for (size_type j=0; j<cols->sparsity_pattern.n_cols(); ++j)
1484  {
1485  for (size_type e=0; e<chunk_size; ++e)
1486  out << '.';
1487  }
1488  else if (std::fabs(val[cols->sparsity_pattern(i,j)]) > threshold)
1489  {
1490  for (size_type e=0; e<chunk_size; ++e)
1491  out << '*';
1492  }
1493  else
1494  {
1495  for (size_type e=0; e<chunk_size; ++e)
1496  out << ':';
1497  }
1498  out << std::endl;
1499  }
1500 
1501  AssertThrow (out, ExcIO());
1502 }
1503 
1504 
1505 
1506 template <typename number>
1507 void
1509 {
1510  AssertThrow (out, ExcIO());
1511 
1512  // first the simple objects, bracketed in [...]
1513  out << '[' << max_len << "][";
1514  // then write out real data
1515  out.write (reinterpret_cast<const char *>(&val[0]),
1516  reinterpret_cast<const char *>(&val[max_len])
1517  - reinterpret_cast<const char *>(&val[0]));
1518  out << ']';
1519 
1520  AssertThrow (out, ExcIO());
1521 }
1522 
1523 
1524 
1525 template <typename number>
1526 void
1528 {
1529  AssertThrow (in, ExcIO());
1530 
1531  char c;
1532 
1533  // first read in simple data
1534  in >> c;
1535  AssertThrow (c == '[', ExcIO());
1536  in >> max_len;
1537 
1538  in >> c;
1539  AssertThrow (c == ']', ExcIO());
1540  in >> c;
1541  AssertThrow (c == '[', ExcIO());
1542 
1543  // reallocate space
1544  delete[] val;
1545  val = new number[max_len];
1546 
1547  // then read data
1548  in.read (reinterpret_cast<char *>(&val[0]),
1549  reinterpret_cast<char *>(&val[max_len])
1550  - reinterpret_cast<char *>(&val[0]));
1551  in >> c;
1552  AssertThrow (c == ']', ExcIO());
1553 }
1554 
1555 
1556 
1557 template <typename number>
1558 std::size_t
1560 {
1561  return sizeof(*this) + max_len*sizeof(number);
1562 }
1563 
1564 
1565 DEAL_II_NAMESPACE_CLOSE
1566 
1567 #endif
size_type n_rows() const
void TSOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
void TSOR(Vector< somenumber > &v, const number om=1.) const
void SSOR(Vector< somenumber > &v, const number omega=1.) const
void print(std::ostream &out) const
::ExceptionBase & ExcMessage(std::string arg1)
void print_pattern(std::ostream &out, const double threshold=0.) const
types::global_dof_index size_type
size_type n_nonzero_elements() const
#define AssertThrow(cond, exc)
Definition: exceptions.h:362
static real_type abs(const number &x)
Definition: numbers.h:316
void SOR(Vector< somenumber > &v, const number om=1.) const
size_type n_actually_nonzero_elements() const
size_type n_nonzero_elements() const
void apply_to_subranges(const RangeType &begin, const typename identity< RangeType >::type &end, const Function &f, const unsigned int grainsize)
Definition: parallel.h:464
void vmult_add(OutVector &dst, const InVector &src) const
SparsityPattern sparsity_pattern
size_type n_nonzero_elements() const
void SSOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
void vmult(OutVector &dst, const InVector &src) const
::ExceptionBase & ExcIO()
size_type get_chunk_size() const
void SOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
size_type n() const
numbers::NumberTraits< number >::real_type real_type
T sum(const T &t, const MPI_Comm &mpi_communicator)
Definition: mpi.h:447
unsigned int global_dof_index
Definition: types.h:100
#define Assert(cond, exc)
Definition: exceptions.h:299
void print_formatted(std::ostream &out, const unsigned int precision=3, const bool scientific=true, const unsigned int width=0, const char *zero_string=" ", const double denominator=1.) const
void block_write(std::ostream &out) const
somenumber matrix_scalar_product(const Vector< somenumber > &u, const Vector< somenumber > &v) const
size_type n() const
static bool equal(const T *p1, const T *p2)
size_type n_cols() const
void PSOR(Vector< somenumber > &v, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const number om=1.) const
void precondition_Jacobi(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1.) const
iterator begin()
ChunkSparseMatrix< number > & operator=(const ChunkSparseMatrix< number > &)
::ExceptionBase & ExcInvalidConstructorCall()
ChunkSparseMatrix< number > & copy_from(const ChunkSparseMatrix< somenumber > &source)
real_type l2_norm() const
void precondition_TSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
size_type m() const
std::size_t * rowstart
size_type n_cols() const
somenumber residual(Vector< somenumber > &dst, const Vector< somenumber > &x, const Vector< somenumber > &b) const
void precondition_SOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
size_type n_rows() const
void precondition_SSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
void block_read(std::istream &in)
std::size_t size() const
::ExceptionBase & ExcNotImplemented()
::ExceptionBase & ExcNotInitialized()
void set(const size_type i, const size_type j, const number value)
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
std::size_t memory_consumption() const
virtual void reinit(const ChunkSparsityPattern &sparsity)
::ExceptionBase & ExcScalarAssignmentOnlyForZeroValue()
void add(const size_type i, const size_type j, const number value)
SmartPointer< const ChunkSparsityPattern, ChunkSparseMatrix< number > > cols
somenumber matrix_norm_square(const Vector< somenumber > &v) const
::ExceptionBase & ExcInternalError()
void Tvmult(OutVector &dst, const InVector &src) const
bool empty() const
void TPSOR(Vector< somenumber > &v, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const number om=1.) const
void Tvmult_add(OutVector &dst, const InVector &src) const
static const size_type invalid_entry