Reference documentation for deal.II version 8.1.0
sparse_matrix.templates.h
1 // ---------------------------------------------------------------------
2 // @f$Id: sparse_matrix.templates.h 31932 2013-12-08 02:15:54Z heister @f$
3 //
4 // Copyright (C) 1999 - 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_matrix_templates_h
19 #define __deal2__sparse_matrix_templates_h
20 
21 
22 #include <deal.II/base/config.h>
23 #include <deal.II/base/template_constraints.h>
24 #include <deal.II/base/parallel.h>
25 #include <deal.II/base/thread_management.h>
26 #include <deal.II/base/multithread_info.h>
27 #include <deal.II/base/utilities.h>
28 #include <deal.II/lac/sparse_matrix.h>
29 #include <deal.II/lac/trilinos_sparse_matrix.h>
30 #include <deal.II/lac/vector.h>
31 #include <deal.II/lac/full_matrix.h>
32 #include <deal.II/lac/compressed_simple_sparsity_pattern.h>
33 #include <deal.II/lac/vector_memory.h>
34 
35 #include <ostream>
36 #include <iomanip>
37 #include <algorithm>
38 #include <functional>
39 #include <cmath>
40 #include <vector>
41 #include <numeric>
42 #include <deal.II/base/std_cxx1x/bind.h>
43 
44 
45 
47 
48 
49 template <typename number>
51  :
52  cols(0, "SparseMatrix"),
53  val(0),
54  max_len(0)
55 {}
56 
57 
58 
59 template <typename number>
61  :
62  Subscriptor (m),
63  cols(0, "SparseMatrix"),
64  val(0),
65  max_len(0)
66 {
70 }
71 
72 
73 
74 template <typename number>
77 {
81 
82  return *this;
83 }
84 
85 
86 
87 template <typename number>
89  :
90  cols(0, "SparseMatrix"),
91  val(0),
92  max_len(0)
93 {
94  reinit (c);
95 }
96 
97 
98 
99 template <typename number>
101  const IdentityMatrix &id)
102  :
103  cols(0, "SparseMatrix"),
104  val(0),
105  max_len(0)
106 {
107  Assert (c.n_rows() == id.m(), ExcDimensionMismatch (c.n_rows(), id.m()));
108  Assert (c.n_cols() == id.n(), ExcDimensionMismatch (c.n_cols(), id.n()));
109 
110  reinit (c);
111  for (size_type i=0; i<n(); ++i)
112  this->set(i,i,1.);
113 }
114 
115 
116 
117 template <typename number>
119 {
120  cols = 0;
121 
122  if (val != 0)
123  delete[] val;
124 }
125 
126 
127 
128 namespace internal
129 {
130  namespace SparseMatrix
131  {
133 
134  template<typename T>
135  void zero_subrange (const size_type begin,
136  const size_type end,
137  T *dst)
138  {
139  std::memset (dst+begin,0,(end-begin)*sizeof(T));
140  }
141  }
142 }
143 
144 
145 
146 template <typename number>
149 {
151 
152  Assert (cols != 0, ExcNotInitialized());
153  Assert (cols->compressed || cols->empty(), SparsityPattern::ExcNotCompressed());
154 
155  // do initial zeroing of elements in parallel. Try to achieve a similar
156  // layout as when doing matrix-vector products, as on some NUMA systems, a
157  // memory block is assigned to memory banks where the first access is
158  // generated. For sparse matrices, the first operations is usually the
159  // operator=. The grain size is chosen to reflect the number of rows in
160  // minimum_parallel_grain_size, weighted by the number of nonzero entries
161  // per row on average.
162  const size_type matrix_size = cols->n_nonzero_elements();
163  const size_type grain_size =
164  internal::SparseMatrix::minimum_parallel_grain_size *
165  (cols->n_nonzero_elements()+m()) / m();
166  if (matrix_size>grain_size)
167  parallel::apply_to_subranges (0U, matrix_size,
168  std_cxx1x::bind(&internal::SparseMatrix::template
169  zero_subrange<number>,
170  std_cxx1x::_1, std_cxx1x::_2,
171  val),
172  grain_size);
173  else if (matrix_size > 0)
174  std::memset (&val[0], 0, matrix_size*sizeof(number));
175 
176  return *this;
177 }
178 
179 
180 
181 template <typename number>
184 {
185  Assert (cols->n_rows() == id.m(),
186  ExcDimensionMismatch (cols->n_rows(), id.m()));
187  Assert (cols->n_cols() == id.n(),
188  ExcDimensionMismatch (cols->n_cols(), id.n()));
189 
190  *this = 0;
191  for (size_type i=0; i<n(); ++i)
192  this->set(i,i,1.);
193 
194  return *this;
195 }
196 
197 
198 
199 template <typename number>
200 void
202 {
203  cols = &sparsity;
204 
205  if (cols->empty())
206  {
207  if (val != 0)
208  delete[] val;
209  val = 0;
210  max_len = 0;
211  return;
212  }
213 
214  const std::size_t N = cols->n_nonzero_elements();
215  if (N > max_len || max_len == 0)
216  {
217  if (val != 0)
218  delete[] val;
219  val = new number[N];
220  max_len = N;
221  }
222 
223  *this = 0.;
224 }
225 
226 
227 
228 template <typename number>
229 void
231 {
232  cols = 0;
233  if (val) delete[] val;
234  val = 0;
235  max_len = 0;
236 }
237 
238 
239 
240 template <typename number>
241 bool
243 {
244  if (cols == 0)
245  return true;
246  else
247  return cols->empty();
248 }
249 
250 
251 
252 template <typename number>
254 SparseMatrix<number>::get_row_length (const size_type row) const
255 {
256  Assert (cols != 0, ExcNotInitialized());
257  return cols->row_length(row);
258 }
259 
260 
261 
262 template <typename number>
265 {
266  Assert (cols != 0, ExcNotInitialized());
267  return cols->n_nonzero_elements ();
268 }
269 
270 
271 
272 template <typename number>
275 {
276  Assert (cols != 0, ExcNotInitialized());
277  Assert (threshold >= 0, ExcMessage ("Negative threshold!"));
278  size_type nnz = 0;
279  const size_type nnz_alloc = n_nonzero_elements();
280  for (size_type i=0; i<nnz_alloc; ++i)
281  if (std::fabs(val[i]) > threshold)
282  ++nnz;
283  return nnz;
284 }
285 
286 
287 
288 template <typename number>
289 void
291 {
292  Assert (cols != 0, ExcNotInitialized());
293  Assert (cols->rows == cols->cols, ExcNotQuadratic());
294 
295  const size_type n_rows = m();
296  for (size_type row=0; row<n_rows; ++row)
297  {
298  // first skip diagonal entry
299  number *val_ptr = &val[cols->rowstart[row]];
300  if (m() == n())
301  ++val_ptr;
302  const size_type *colnum_ptr = &cols->colnums[cols->rowstart[row]+1];
303  const number *const val_end_of_row = &val[cols->rowstart[row+1]];
304 
305  // treat lower left triangle
306  while ((val_ptr != val_end_of_row) && (*colnum_ptr<row))
307  {
308  // compute the mean of this
309  // and the transpose value
310  const number mean_value = (*val_ptr +
311  val[(*cols)(*colnum_ptr,row)]) / 2.0;
312  // set this value and the
313  // transpose one to the
314  // mean
315  *val_ptr = mean_value;
316  set (*colnum_ptr, row, mean_value);
317 
318  // advance pointers
319  ++val_ptr;
320  ++colnum_ptr;
321  };
322  };
323 }
324 
325 
326 
327 template <typename number>
328 template <typename somenumber>
331 {
332  Assert (cols != 0, ExcNotInitialized());
333  Assert (val != 0, ExcNotInitialized());
334  Assert (cols == matrix.cols, ExcDifferentSparsityPatterns());
335 
336  std::copy (&matrix.val[0], &matrix.val[cols->n_nonzero_elements()],
337  &val[0]);
338 
339  return *this;
340 }
341 
342 
343 
344 template <typename number>
345 template <typename somenumber>
346 void
348 {
349  // first delete previous content
350  *this = 0;
351 
352  // then copy old matrix
353  for (size_type row=0; row<matrix.m(); ++row)
354  for (size_type col=0; col<matrix.n(); ++col)
355  if (matrix(row,col) != 0)
356  set (row, col, matrix(row,col));
357 }
358 
359 
360 
361 #ifdef DEAL_II_WITH_TRILINOS
362 
363 template <typename number>
366 {
367  Assert (m() == matrix.m(), ExcDimensionMismatch(m(), matrix.m()));
368  Assert (n() == matrix.n(), ExcDimensionMismatch(n(), matrix.n()));
369 
370  // first delete previous content
371  *this = 0;
372 
373  std::vector < TrilinosScalar > value_cache;
374  std::vector<size_type> colnum_cache;
375 
376  for (size_type row = 0; row < matrix.m(); ++row)
377  {
378  value_cache.resize(matrix.n());
379  colnum_cache.resize(matrix.n());
380 
381  // copy column indices and values and at the same time enquire about the
382  // length of the row
383  int ncols;
384  int ierr
385  = matrix.trilinos_matrix().ExtractGlobalRowCopy
386  (row, matrix.row_length(row), ncols,
387  &(value_cache[0]),
388  reinterpret_cast<TrilinosWrappers::types::int_type *>(&(colnum_cache[0])));
389  Assert (ierr==0, ExcTrilinosError(ierr));
390 
391  // resize arrays to the size actually used
392  value_cache.resize(ncols);
393  colnum_cache.resize(ncols);
394 
395  // then copy everything in one swoop
396  this->set(row,
397  colnum_cache,
398  value_cache);
399  }
400 
401  return *this;
402 }
403 
404 #endif
405 
406 
407 template <typename number>
408 template <typename somenumber>
409 void
410 SparseMatrix<number>::add (const number factor,
411  const SparseMatrix<somenumber> &matrix)
412 {
413  Assert (cols != 0, ExcNotInitialized());
414  Assert (val != 0, ExcNotInitialized());
415  Assert (cols == matrix.cols, ExcDifferentSparsityPatterns());
416 
417  number *val_ptr = &val[0];
418  const somenumber *matrix_ptr = &matrix.val[0];
419  const number *const end_ptr = &val[cols->n_nonzero_elements()];
420 
421  while (val_ptr != end_ptr)
422  *val_ptr++ += factor **matrix_ptr++;
423 }
424 
425 
426 
427 namespace internal
428 {
429  namespace SparseMatrix
430  {
439  template <typename number,
440  typename InVector,
441  typename OutVector>
442  void vmult_on_subrange (const size_type begin_row,
443  const size_type end_row,
444  const number *values,
445  const std::size_t *rowstart,
446  const size_type *colnums,
447  const InVector &src,
448  OutVector &dst,
449  const bool add)
450  {
451  const number *val_ptr = &values[rowstart[begin_row]];
452  const size_type *colnum_ptr = &colnums[rowstart[begin_row]];
453  typename OutVector::iterator dst_ptr = dst.begin() + begin_row;
454 
455  if (add == false)
456  for (size_type row=begin_row; row<end_row; ++row)
457  {
458  typename OutVector::value_type s = 0.;
459  const number *const val_end_of_row = &values[rowstart[row+1]];
460  while (val_ptr != val_end_of_row)
461  s += *val_ptr++ * src(*colnum_ptr++);
462  *dst_ptr++ = s;
463  }
464  else
465  for (size_type row=begin_row; row<end_row; ++row)
466  {
467  typename OutVector::value_type s = *dst_ptr;
468  const number *const val_end_of_row = &values[rowstart[row+1]];
469  while (val_ptr != val_end_of_row)
470  s += *val_ptr++ * src(*colnum_ptr++);
471  *dst_ptr++ = s;
472  }
473  }
474  }
475 }
476 
477 
478 template <typename number>
479 template <typename number2>
480 void
481 SparseMatrix<number>::add (const size_type row,
482  const size_type n_cols,
483  const size_type *col_indices,
484  const number2 *values,
485  const bool elide_zero_values,
486  const bool col_indices_are_sorted)
487 {
488  Assert (cols != 0, ExcNotInitialized());
489 
490  // if we have sufficiently many columns
491  // and sorted indices it is faster to
492  // just go through the column indices and
493  // look whether we found one, rather than
494  // doing many binary searches
495  if (elide_zero_values == false && col_indices_are_sorted == true &&
496  n_cols > 3)
497  {
498  // check whether the given indices are
499  // really sorted
500 #ifdef DEBUG
501  for (size_type i=1; i<n_cols; ++i)
502  Assert (col_indices[i] > col_indices[i-1],
503  ExcMessage("List of indices is unsorted or contains duplicates."));
504 #endif
505 
506  const size_type *this_cols =
507  &cols->colnums[cols->rowstart[row]];
508  const size_type row_length_1 = cols->row_length(row)-1;
509  number *val_ptr = &val[cols->rowstart[row]];
510 
511  if (m() == n())
512  {
513 
514  // find diagonal and add it if found
515  Assert (this_cols[0] == row, ExcInternalError());
516  const size_type *diag_pos =
517  Utilities::lower_bound (col_indices,
518  col_indices+n_cols,
519  row);
520  const size_type diag = diag_pos - col_indices;
521  size_type post_diag = diag;
522  if (diag != n_cols && *diag_pos == row)
523  {
524  val_ptr[0] += *(values+(diag_pos-col_indices));
525  ++post_diag;
526  }
527 
528  // add indices before diagonal
529  size_type counter = 1;
530  for (size_type i=0; i<diag; ++i)
531  {
532  while (this_cols[counter]<col_indices[i] && counter<row_length_1)
533  ++counter;
534 
535  Assert (this_cols[counter] == col_indices[i] || values[i] == 0,
536  ExcInvalidIndex(row,col_indices[i]));
537 
538  val_ptr[counter] += values[i];
539  }
540 
541  // add indices after diagonal
542  for (size_type i=post_diag; i<n_cols; ++i)
543  {
544  while (this_cols[counter]<col_indices[i] && counter<row_length_1)
545  ++counter;
546 
547  Assert (this_cols[counter] == col_indices[i] || values[i] == 0,
548  ExcInvalidIndex(row,col_indices[i]));
549 
550  val_ptr[counter] += values[i];
551  }
552 
553  Assert (counter < cols->row_length(row),
554  ExcMessage("Specified invalid column indices in add function."));
555  }
556  else
557  {
558  size_type counter = 0;
559  for (size_type i=0; i<n_cols; ++i)
560  {
561  while (this_cols[counter]<col_indices[i] && counter<row_length_1)
562  ++counter;
563 
564  Assert (this_cols[counter] == col_indices[i] || values[i] == 0,
565  ExcInvalidIndex(row,col_indices[i]));
566 
567  val_ptr[counter] += values[i];
568  }
569  Assert (counter < cols->row_length(row),
570  ExcMessage("Specified invalid column indices in add function."));
571  }
572  return;
573  }
574 
575  // unsorted case: first, search all the
576  // indices to find out which values we
577  // actually need to add.
578  const size_type *const my_cols = cols->colnums;
579  size_type index = cols->rowstart[row];
580  const size_type next_row_index = cols->rowstart[row+1];
581 
582  for (size_type j=0; j<n_cols; ++j)
583  {
584  const number value = values[j];
586 
587 #ifdef DEBUG
588  if (elide_zero_values==true && value == 0)
589  continue;
590 #else
591  if (value == 0)
592  continue;
593 #endif
594 
595  // check whether the next index to add is
596  // the next present index in the sparsity
597  // pattern (otherwise, do a binary
598  // search)
599  if (index < next_row_index && my_cols[index] == col_indices[j])
600  goto add_value;
601 
602  index = cols->operator()(row, col_indices[j]);
603 
604  // it is allowed to add elements to
605  // the matrix that are not part of
606  // the sparsity pattern, if the
607  // value we add is zero
608  if (index == SparsityPattern::invalid_entry)
609  {
610  Assert (value == 0., ExcInvalidIndex(row,col_indices[j]));
611  continue;
612  }
613 
614 add_value:
615  val[index] += value;
616  ++index;
617  }
618 }
619 
620 
621 
622 template <typename number>
623 template <typename number2>
624 void
625 SparseMatrix<number>::set (const size_type row,
626  const size_type n_cols,
627  const size_type *col_indices,
628  const number2 *values,
629  const bool elide_zero_values)
630 {
631  Assert (cols != 0, ExcNotInitialized());
632  Assert (row < m(), ExcInvalidIndex1(row));
633 
634  // First, search all the indices to find
635  // out which values we actually need to
636  // set.
637  const size_type *my_cols = cols->colnums;
638  std::size_t index = cols->rowstart[row], next_index = index;
639  const std::size_t next_row_index = cols->rowstart[row+1];
640 
641  if (elide_zero_values == true)
642  {
643  for (size_type j=0; j<n_cols; ++j)
644  {
645  const number value = values[j];
647 
648  if (value == 0)
649  continue;
650 
651  // check whether the next index to set is
652  // the next present index in the sparsity
653  // pattern (otherwise, do a binary
654  // search)
655  if (index != next_row_index && my_cols[index] == col_indices[j])
656  goto set_value;
657 
658  next_index = cols->operator()(row, col_indices[j]);
659 
660  // it is allowed to set elements in
661  // the matrix that are not part of
662  // the sparsity pattern, if the
663  // value to which we set it is zero
664  if (next_index == SparsityPattern::invalid_entry)
665  {
666  Assert (false, ExcInvalidIndex(row,col_indices[j]));
667  continue;
668  }
669  index = next_index;
670 
671 set_value:
672  val[index] = value;
673  ++index;
674  }
675  }
676  else
677  {
678  // same code as above, but now check for zeros
679  for (size_type j=0; j<n_cols; ++j)
680  {
681  const number value = values[j];
683 
684  if (index != next_row_index && my_cols[index] == col_indices[j])
685  goto set_value_checked;
686 
687  next_index = cols->operator()(row, col_indices[j]);
688 
689  if (next_index == SparsityPattern::invalid_entry)
690  {
691  Assert (value == 0., ExcInvalidIndex(row,col_indices[j]));
692  continue;
693  }
694  index = next_index;
695 
696 set_value_checked:
697  val[index] = value;
698  ++index;
699  }
700  }
701 }
702 
703 
704 
705 template <typename number>
706 template <class OutVector, class InVector>
707 void
709  const InVector &src) const
710 {
711  Assert (cols != 0, ExcNotInitialized());
712  Assert (val != 0, ExcNotInitialized());
713  Assert(m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
714  Assert(n() == src.size(), ExcDimensionMismatch(n(),src.size()));
715 
716  Assert (!PointerComparison::equal(&src, &dst), ExcSourceEqualsDestination());
717 
719  std_cxx1x::bind (&internal::SparseMatrix::vmult_on_subrange
720  <number,InVector,OutVector>,
721  std_cxx1x::_1, std_cxx1x::_2,
722  val,
723  cols->rowstart,
724  cols->colnums,
725  std_cxx1x::cref(src),
726  std_cxx1x::ref(dst),
727  false),
728  internal::SparseMatrix::minimum_parallel_grain_size);
729 }
730 
731 
732 
733 template <typename number>
734 template <class OutVector, class InVector>
735 void
737  const InVector &src) const
738 {
739  Assert (val != 0, ExcNotInitialized());
740  Assert (cols != 0, ExcNotInitialized());
741  Assert(n() == dst.size(), ExcDimensionMismatch(n(),dst.size()));
742  Assert(m() == src.size(), ExcDimensionMismatch(m(),src.size()));
743 
744  Assert (!PointerComparison::equal(&src, &dst), ExcSourceEqualsDestination());
745 
746  dst = 0;
747 
748  for (size_type i=0; i<m(); i++)
749  {
750  for (size_type j=cols->rowstart[i]; j<cols->rowstart[i+1] ; j++)
751  {
752  const size_type p = cols->colnums[j];
753  dst(p) += val[j] * src(i);
754  }
755  }
756 }
757 
758 
759 
760 template <typename number>
761 template <class OutVector, class InVector>
762 void
764  const InVector &src) const
765 {
766  Assert (cols != 0, ExcNotInitialized());
767  Assert (val != 0, ExcNotInitialized());
768  Assert(m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
769  Assert(n() == src.size(), ExcDimensionMismatch(n(),src.size()));
770 
771  Assert (!PointerComparison::equal(&src, &dst), ExcSourceEqualsDestination());
772 
774  std_cxx1x::bind (&internal::SparseMatrix::vmult_on_subrange
775  <number,InVector,OutVector>,
776  std_cxx1x::_1, std_cxx1x::_2,
777  val,
778  cols->rowstart,
779  cols->colnums,
780  std_cxx1x::cref(src),
781  std_cxx1x::ref(dst),
782  true),
783  internal::SparseMatrix::minimum_parallel_grain_size);
784 }
785 
786 
787 
788 template <typename number>
789 template <class OutVector, class InVector>
790 void
792  const InVector &src) const
793 {
794  Assert (val != 0, ExcNotInitialized());
795  Assert (cols != 0, ExcNotInitialized());
796  Assert(n() == dst.size(), ExcDimensionMismatch(n(),dst.size()));
797  Assert(m() == src.size(), ExcDimensionMismatch(m(),src.size()));
798 
799  Assert (!PointerComparison::equal(&src, &dst), ExcSourceEqualsDestination());
800 
801  for (size_type i=0; i<m(); i++)
802  for (size_type j=cols->rowstart[i]; j<cols->rowstart[i+1] ; j++)
803  {
804  const size_type p = cols->colnums[j];
805  dst(p) += val[j] * src(i);
806  }
807 }
808 
809 
810 namespace internal
811 {
812  namespace SparseMatrix
813  {
825  template <typename number,
826  typename InVector>
827  number matrix_norm_sqr_on_subrange (const size_type begin_row,
828  const size_type end_row,
829  const number *values,
830  const std::size_t *rowstart,
831  const size_type *colnums,
832  const InVector &v)
833  {
834  number norm_sqr=0.;
835 
836  for (size_type i=begin_row; i<end_row; ++i)
837  {
838  number s = 0;
839  for (size_type j=rowstart[i]; j<rowstart[i+1] ; j++)
840  s += values[j] * v(colnums[j]);
841  norm_sqr += v(i)*numbers::NumberTraits<number>::conjugate(s);
842  }
843  return norm_sqr;
844  }
845  }
846 }
847 
848 
849 
850 template <typename number>
851 template <typename somenumber>
852 somenumber
854 {
855  Assert (cols != 0, ExcNotInitialized());
856  Assert (val != 0, ExcNotInitialized());
857  Assert(m() == v.size(), ExcDimensionMismatch(m(),v.size()));
858  Assert(n() == v.size(), ExcDimensionMismatch(n(),v.size()));
859 
860  return
861  parallel::accumulate_from_subranges<number>
862  (std_cxx1x::bind (&internal::SparseMatrix::matrix_norm_sqr_on_subrange
863  <number,Vector<somenumber> >,
864  std_cxx1x::_1, std_cxx1x::_2,
865  val, cols->rowstart, cols->colnums,
866  std_cxx1x::cref(v)),
867  0, m(),
868  internal::SparseMatrix::minimum_parallel_grain_size);
869 }
870 
871 
872 
873 namespace internal
874 {
875  namespace SparseMatrix
876  {
888  template <typename number,
889  typename InVector>
890  number matrix_scalar_product_on_subrange (const size_type begin_row,
891  const size_type end_row,
892  const number *values,
893  const std::size_t *rowstart,
894  const size_type *colnums,
895  const InVector &u,
896  const InVector &v)
897  {
898  number norm_sqr=0.;
899 
900  for (size_type i=begin_row; i<end_row; ++i)
901  {
902  number s = 0;
903  for (size_type j=rowstart[i]; j<rowstart[i+1] ; j++)
904  s += values[j] * v(colnums[j]);
905  norm_sqr += u(i)*numbers::NumberTraits<number>::conjugate(s);
906  }
907  return norm_sqr;
908  }
909  }
910 }
911 
912 
913 
914 template <typename number>
915 template <typename somenumber>
916 somenumber
918  const Vector<somenumber> &v) const
919 {
920  Assert (cols != 0, ExcNotInitialized());
921  Assert (val != 0, ExcNotInitialized());
922  Assert(m() == u.size(), ExcDimensionMismatch(m(),u.size()));
923  Assert(n() == v.size(), ExcDimensionMismatch(n(),v.size()));
924 
925  return
926  parallel::accumulate_from_subranges<number>
927  (std_cxx1x::bind (&internal::SparseMatrix::matrix_scalar_product_on_subrange
928  <number,Vector<somenumber> >,
929  std_cxx1x::_1, std_cxx1x::_2,
930  val, cols->rowstart, cols->colnums,
931  std_cxx1x::cref(u),
932  std_cxx1x::cref(v)),
933  0, m(),
934  internal::SparseMatrix::minimum_parallel_grain_size);
935 }
936 
937 
938 
939 template <typename number>
940 template <typename numberB, typename numberC>
941 void
943  const SparseMatrix<numberB> &B,
944  const Vector<number> &V,
945  const bool rebuild_sparsity_C) const
946 {
947  const bool use_vector = V.size() == n() ? true : false;
948  Assert (n() == B.m(), ExcDimensionMismatch(n(), B.m()));
949  Assert (cols != 0, ExcNotInitialized());
950  Assert (B.cols != 0, ExcNotInitialized());
951  Assert (C.cols != 0, ExcNotInitialized());
952 
953  const SparsityPattern &sp_A = *cols;
954  const SparsityPattern &sp_B = *B.cols;
955 
956  // clear previous content of C
957  if (rebuild_sparsity_C == true)
958  {
959  // we are about to change the sparsity pattern of C. this can not work
960  // if either A or B use the same sparsity pattern
961  Assert (&C.get_sparsity_pattern() != &this->get_sparsity_pattern(),
962  ExcMessage ("Can't use the same sparsity pattern for "
963  "different matrices if it is to be rebuilt."));
965  ExcMessage ("Can't use the same sparsity pattern for "
966  "different matrices if it is to be rebuilt."));
967 
968  // need to change the sparsity pattern of C, so cast away const-ness.
969  SparsityPattern &sp_C =
970  *(const_cast<SparsityPattern *>(&C.get_sparsity_pattern()));
971  C.clear();
972  sp_C.reinit (0,0,0);
973 
974  // create a sparsity pattern for the matrix. we will go through all the
975  // rows in the matrix A, and for each column in a row we add the whole
976  // row of matrix B with that row number. This means that we will insert
977  // a lot of entries to each row, which is best handled by the
978  // CompressedSimpleSparsityPattern class.
979  {
980  CompressedSimpleSparsityPattern csp (m(), B.n());
981  for (size_type i = 0; i < csp.n_rows(); ++i)
982  {
983  const size_type *rows = &sp_A.colnums[sp_A.rowstart[i]];
984  const size_type *const end_rows =
985  &sp_A.colnums[sp_A.rowstart[i+1]];
986  for (; rows != end_rows; ++rows)
987  {
988  const size_type col = *rows;
989  size_type *new_cols = const_cast<size_type *>
990  (&sp_B.colnums[sp_B.rowstart[col]]);
991  size_type *end_new_cols = const_cast<size_type *>
992  (&sp_B.colnums[sp_B.rowstart[col+1]]);
993 
994  // if B has a diagonal, need to add that manually. this way,
995  // we maintain sortedness.
996  if (sp_B.n_rows() == sp_B.n_cols())
997  {
998  ++new_cols;
999  csp.add(i, col);
1000  }
1001 
1002  csp.add_entries (i, new_cols, end_new_cols, true);
1003  }
1004  }
1005  sp_C.copy_from (csp);
1006  }
1007 
1008  // reinit matrix C from that information
1009  C.reinit (sp_C);
1010  }
1011 
1012  Assert (C.m() == m(), ExcDimensionMismatch(C.m(), m()));
1013  Assert (C.n() == B.n(), ExcDimensionMismatch(C.n(), B.n()));
1014 
1015  // create an array that caches some
1016  // elements that are going to be written
1017  // into the new matrix.
1018  unsigned int max_n_cols_B = 0;
1019  for (size_type i=0; i<B.m(); ++i)
1020  max_n_cols_B = std::max (max_n_cols_B, sp_B.row_length(i));
1021  std::vector<numberC> new_entries(max_n_cols_B);
1022 
1023  // now compute the actual entries: a matrix-matrix product involves three
1024  // nested loops. One over the rows of A, for each row we then loop over all
1025  // the columns, and then we need to multiply each element with all the
1026  // elements in that row in B.
1027  for (size_type i=0; i<C.m(); ++i)
1028  {
1029  const size_type *rows = &sp_A.colnums[sp_A.rowstart[i]];
1030  const size_type *const end_rows = &sp_A.colnums[sp_A.rowstart[i+1]];
1031  for (; rows != end_rows; ++rows)
1032  {
1033  const double A_val = val[rows-&sp_A.colnums[sp_A.rowstart[0]]];
1034  const size_type col = *rows;
1035  const size_type *new_cols =
1036  (&sp_B.colnums[sp_B.rowstart[col]]);
1037 
1038  // special treatment for diagonal
1039  if (sp_B.n_rows() == sp_B.n_cols())
1040  {
1041  C.add (i, *new_cols, A_val *
1042  B.val[new_cols-&sp_B.colnums[sp_B.rowstart[0]]] *
1043  (use_vector ? V(col) : 1));
1044  ++new_cols;
1045  }
1046 
1047  // now the innermost loop that goes over all the elements in row
1048  // 'col' of matrix B. Cache the elements, and then write them into C
1049  // at once
1050  numberC *new_ptr = &new_entries[0];
1051  const numberB *B_val_ptr =
1052  &B.val[new_cols-&sp_B.colnums[sp_B.rowstart[0]]];
1053  const numberB *const end_cols = &B.val[sp_B.rowstart[col+1]];
1054  for (; B_val_ptr != end_cols; ++B_val_ptr)
1055  *new_ptr++ = A_val **B_val_ptr * (use_vector ? V(col) : 1);
1056 
1057  C.add (i, new_ptr-&new_entries[0], new_cols, &new_entries[0],
1058  false, true);
1059  }
1060  }
1061 }
1062 
1063 
1064 
1065 
1066 template <typename number>
1067 template <typename numberB, typename numberC>
1068 void
1070  const SparseMatrix<numberB> &B,
1071  const Vector<number> &V,
1072  const bool rebuild_sparsity_C) const
1073 {
1074  const bool use_vector = V.size() == m() ? true : false;
1075  Assert (m() == B.m(), ExcDimensionMismatch(m(), B.m()));
1076  Assert (cols != 0, ExcNotInitialized());
1077  Assert (B.cols != 0, ExcNotInitialized());
1078  Assert (C.cols != 0, ExcNotInitialized());
1079 
1080  const SparsityPattern &sp_A = *cols;
1081  const SparsityPattern &sp_B = *B.cols;
1082 
1083  // clear previous content of C
1084  if (rebuild_sparsity_C == true)
1085  {
1086  // we are about to change the sparsity pattern of C. this can not work
1087  // if either A or B use the same sparsity pattern
1088  Assert (&C.get_sparsity_pattern() != &this->get_sparsity_pattern(),
1089  ExcMessage ("Can't use the same sparsity pattern for "
1090  "different matrices if it is to be rebuilt."));
1092  ExcMessage ("Can't use the same sparsity pattern for "
1093  "different matrices if it is to be rebuilt."));
1094 
1095  // need to change the sparsity pattern of C, so cast away const-ness.
1096  SparsityPattern &sp_C =
1097  *(const_cast<SparsityPattern *>(&C.get_sparsity_pattern()));
1098  C.clear();
1099  sp_C.reinit (0,0,0);
1100 
1101  // create a sparsity pattern for the matrix. we will go through all the
1102  // rows in the matrix A, and for each column in a row we add the whole
1103  // row of matrix B with that row number. This means that we will insert
1104  // a lot of entries to each row, which is best handled by the
1105  // CompressedSimpleSparsityPattern class.
1106  {
1107  CompressedSimpleSparsityPattern csp (n(), B.n());
1108  for (size_type i = 0; i < sp_A.n_rows(); ++i)
1109  {
1110  const size_type *rows =
1111  &sp_A.colnums[sp_A.rowstart[i]];
1112  const size_type *const end_rows =
1113  &sp_A.colnums[sp_A.rowstart[i+1]];
1114  // cast away constness to conform with csp.add_entries interface
1115  size_type *new_cols = const_cast<size_type *>
1116  (&sp_B.colnums[sp_B.rowstart[i]]);
1117  size_type *end_new_cols = const_cast<size_type *>
1118  (&sp_B.colnums[sp_B.rowstart[i+1]]);
1119 
1120  if (sp_B.n_rows() == sp_B.n_cols())
1121  ++new_cols;
1122 
1123  for (; rows != end_rows; ++rows)
1124  {
1125  const size_type row = *rows;
1126 
1127  // if B has a diagonal, need to add that manually. this way,
1128  // we maintain sortedness.
1129  if (sp_B.n_rows() == sp_B.n_cols())
1130  csp.add(row, i);
1131 
1132  csp.add_entries (row, new_cols, end_new_cols, true);
1133  }
1134  }
1135  sp_C.copy_from (csp);
1136  }
1137 
1138  // reinit matrix C from that information
1139  C.reinit (sp_C);
1140  }
1141 
1142  Assert (C.m() == n(), ExcDimensionMismatch(C.m(), n()));
1143  Assert (C.n() == B.n(), ExcDimensionMismatch(C.n(), B.n()));
1144 
1145  // create an array that caches some
1146  // elements that are going to be written
1147  // into the new matrix.
1148  unsigned int max_n_cols_B = 0;
1149  for (size_type i=0; i<B.m(); ++i)
1150  max_n_cols_B = std::max (max_n_cols_B, sp_B.row_length(i));
1151  std::vector<numberC> new_entries(max_n_cols_B);
1152 
1153  // now compute the actual entries: a matrix-matrix product involves three
1154  // nested loops. One over the rows of A, for each row we then loop over all
1155  // the columns, and then we need to multiply each element with all the
1156  // elements in that row in B.
1157  for (size_type i=0; i<m(); ++i)
1158  {
1159  const size_type *rows = &sp_A.colnums[sp_A.rowstart[i]];
1160  const size_type *const end_rows = &sp_A.colnums[sp_A.rowstart[i+1]];
1161  const size_type *new_cols = &sp_B.colnums[sp_B.rowstart[i]];
1162  if (sp_B.n_rows() == sp_B.n_cols())
1163  ++new_cols;
1164 
1165  const numberB *const end_cols = &B.val[sp_B.rowstart[i+1]];
1166 
1167  for (; rows != end_rows; ++rows)
1168  {
1169  const size_type row = *rows;
1170  const double A_val = val[rows-&sp_A.colnums[sp_A.rowstart[0]]];
1171 
1172  // special treatment for diagonal
1173  if (sp_B.n_rows () == sp_B.n_cols())
1174  C.add (row, i, A_val *
1175  B.val[new_cols-1-&sp_B.colnums[sp_B.rowstart[0]]] *
1176  (use_vector ? V(i) : 1));
1177 
1178  // now the innermost loop that goes over all the elements in row
1179  // 'col' of matrix B. Cache the elements, and then write them into C
1180  // at once
1181  numberC *new_ptr = &new_entries[0];
1182  const numberB *B_val_ptr =
1183  &B.val[new_cols-&sp_B.colnums[sp_B.rowstart[0]]];
1184  for (; B_val_ptr != end_cols; ++B_val_ptr)
1185  *new_ptr++ = A_val **B_val_ptr * (use_vector ? V(i) : 1);
1186 
1187  C.add (row, new_ptr-&new_entries[0], new_cols, &new_entries[0],
1188  false, true);
1189  }
1190  }
1191 }
1192 
1193 
1194 
1195 template <typename number>
1198 {
1199  Assert (cols != 0, ExcNotInitialized());
1200  Assert (val != 0, ExcNotInitialized());
1201 
1202  Vector<real_type> column_sums(n());
1203  const size_type n_rows = m();
1204  for (size_type row=0; row<n_rows; ++row)
1205  for (size_type j=cols->rowstart[row]; j<cols->rowstart[row+1] ; ++j)
1206  column_sums(cols->colnums[j]) += numbers::NumberTraits<number>::abs(val[j]);
1207 
1208  return column_sums.linfty_norm();
1209 }
1210 
1211 
1212 
1213 template <typename number>
1216 {
1217  Assert (cols != 0, ExcNotInitialized());
1218  Assert (val != 0, ExcNotInitialized());
1219 
1220  const number *val_ptr = &val[cols->rowstart[0]];
1221 
1222  real_type max=0;
1223  const size_type n_rows = m();
1224  for (size_type row=0; row<n_rows; ++row)
1225  {
1226  real_type sum = 0;
1227  const number *const val_end_of_row = &val[cols->rowstart[row+1]];
1228  while (val_ptr != val_end_of_row)
1229  sum += numbers::NumberTraits<number>::abs(*val_ptr++);
1230  if (sum > max)
1231  max = sum;
1232  }
1233  return max;
1234 }
1235 
1236 
1237 
1238 template <typename number>
1241 {
1242  // simply add up all entries in the
1243  // sparsity pattern, without taking any
1244  // reference to rows or columns
1245  real_type norm_sqr = 0;
1246  const size_type n_rows = m();
1247  for (const number *ptr = &val[0];
1248  ptr != &val[cols->rowstart[n_rows]]; ++ptr)
1250 
1251  return std::sqrt (norm_sqr);
1252 }
1253 
1254 
1255 
1256 namespace internal
1257 {
1258  namespace SparseMatrix
1259  {
1271  template <typename number,
1272  typename InVector,
1273  typename OutVector>
1274  number residual_sqr_on_subrange (const size_type begin_row,
1275  const size_type end_row,
1276  const number *values,
1277  const std::size_t *rowstart,
1278  const size_type *colnums,
1279  const InVector &u,
1280  const InVector &b,
1281  OutVector &dst)
1282  {
1283  number norm_sqr=0.;
1284 
1285  for (size_type i=begin_row; i<end_row; ++i)
1286  {
1287  number s = b(i);
1288  for (size_type j=rowstart[i]; j<rowstart[i+1] ; j++)
1289  s -= values[j] * u(colnums[j]);
1290  dst(i) = s;
1292  }
1293  return norm_sqr;
1294  }
1295  }
1296 }
1297 
1298 
1299 template <typename number>
1300 template <typename somenumber>
1301 somenumber
1303  const Vector<somenumber> &u,
1304  const Vector<somenumber> &b) const
1305 {
1306  Assert (cols != 0, ExcNotInitialized());
1307  Assert (val != 0, ExcNotInitialized());
1308  Assert(m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
1309  Assert(m() == b.size(), ExcDimensionMismatch(m(),b.size()));
1310  Assert(n() == u.size(), ExcDimensionMismatch(n(),u.size()));
1311 
1312  Assert (&u != &dst, ExcSourceEqualsDestination());
1313 
1314  return
1315  std::sqrt (parallel::accumulate_from_subranges<number>
1316  (std_cxx1x::bind (&internal::SparseMatrix::residual_sqr_on_subrange
1318  std_cxx1x::_1, std_cxx1x::_2,
1319  val, cols->rowstart, cols->colnums,
1320  std_cxx1x::cref(u),
1321  std_cxx1x::cref(b),
1322  std_cxx1x::ref(dst)),
1323  0, m(),
1324  internal::SparseMatrix::minimum_parallel_grain_size));
1325 }
1326 
1327 
1328 
1329 template <typename number>
1330 template <typename somenumber>
1331 void
1333  const Vector<somenumber> &src,
1334  const number om) const
1335 {
1336  Assert (cols != 0, ExcNotInitialized());
1337  Assert (val != 0, ExcNotInitialized());
1338  AssertDimension (m(), n());
1339  AssertDimension (dst.size(), n());
1340  AssertDimension (src.size(), n());
1341 
1342  const size_type n = src.size();
1343  somenumber *dst_ptr = dst.begin();
1344  const somenumber *src_ptr = src.begin();
1345  const std::size_t *rowstart_ptr = &cols->rowstart[0];
1346 
1347  // optimize the following loop for
1348  // the case that the relaxation
1349  // factor is one. In that case, we
1350  // can save one FP multiplication
1351  // per row
1352  //
1353  // note that for square matrices,
1354  // the diagonal entry is the first
1355  // in each row, i.e. at index
1356  // rowstart[i]. and we do have a
1357  // square matrix by above assertion
1358  if (om != 1.)
1359  for (size_type i=0; i<n; ++i, ++dst_ptr, ++src_ptr, ++rowstart_ptr)
1360  *dst_ptr = om **src_ptr / val[*rowstart_ptr];
1361  else
1362  for (size_type i=0; i<n; ++i, ++dst_ptr, ++src_ptr, ++rowstart_ptr)
1363  *dst_ptr = *src_ptr / val[*rowstart_ptr];
1364 }
1365 
1366 
1367 
1368 template <typename number>
1369 template <typename somenumber>
1370 void
1372  const Vector<somenumber> &src,
1373  const number om,
1374  const std::vector<std::size_t> &pos_right_of_diagonal) const
1375 {
1376  // to understand how this function works
1377  // you may want to take a look at the CVS
1378  // archives to see the original version
1379  // which is much clearer...
1380  Assert (cols != 0, ExcNotInitialized());
1381  Assert (val != 0, ExcNotInitialized());
1382  AssertDimension (m(), n());
1383  AssertDimension (dst.size(), n());
1384  AssertDimension (src.size(), n());
1385 
1386  const size_type n = src.size();
1387  const std::size_t *rowstart_ptr = &cols->rowstart[0];
1388  somenumber *dst_ptr = &dst(0);
1389 
1390  // case when we have stored the position
1391  // just right of the diagonal (then we
1392  // don't have to search for it).
1393  if (pos_right_of_diagonal.size() != 0)
1394  {
1395  Assert (pos_right_of_diagonal.size() == dst.size(),
1396  ExcDimensionMismatch (pos_right_of_diagonal.size(), dst.size()));
1397 
1398  // forward sweep
1399  for (size_type row=0; row<n; ++row, ++dst_ptr, ++rowstart_ptr)
1400  {
1401  *dst_ptr = src(row);
1402  const std::size_t first_right_of_diagonal_index =
1403  pos_right_of_diagonal[row];
1404  Assert (first_right_of_diagonal_index <= *(rowstart_ptr+1),
1405  ExcInternalError());
1406  number s = 0;
1407  for (size_type j=(*rowstart_ptr)+1; j<first_right_of_diagonal_index; ++j)
1408  s += val[j] * dst(cols->colnums[j]);
1409 
1410  // divide by diagonal element
1411  *dst_ptr -= s * om;
1412  Assert(val[*rowstart_ptr]!= 0., ExcDivideByZero());
1413  *dst_ptr /= val[*rowstart_ptr];
1414  };
1415 
1416  rowstart_ptr = &cols->rowstart[0];
1417  dst_ptr = &dst(0);
1418  for ( ; rowstart_ptr!=&cols->rowstart[n]; ++rowstart_ptr, ++dst_ptr)
1419  *dst_ptr *= om*(2.-om)*val[*rowstart_ptr];
1420 
1421  // backward sweep
1422  rowstart_ptr = &cols->rowstart[n-1];
1423  dst_ptr = &dst(n-1);
1424  for (int row=n-1; row>=0; --row, --rowstart_ptr, --dst_ptr)
1425  {
1426  const size_type end_row = *(rowstart_ptr+1);
1427  const size_type first_right_of_diagonal_index
1428  = pos_right_of_diagonal[row];
1429  number s = 0;
1430  for (size_type j=first_right_of_diagonal_index; j<end_row; ++j)
1431  s += val[j] * dst(cols->colnums[j]);
1432 
1433  *dst_ptr -= s * om;
1434  Assert(val[*rowstart_ptr]!= 0., ExcDivideByZero());
1435  *dst_ptr /= val[*rowstart_ptr];
1436  };
1437  return;
1438  }
1439 
1440  // case when we need to get the position
1441  // of the first element right of the
1442  // diagonal manually for each sweep.
1443  // forward sweep
1444  for (size_type row=0; row<n; ++row, ++dst_ptr, ++rowstart_ptr)
1445  {
1446  *dst_ptr = src(row);
1447  // find the first element in this line
1448  // which is on the right of the diagonal.
1449  // we need to precondition with the
1450  // elements on the left only.
1451  // note: the first entry in each
1452  // line denotes the diagonal element,
1453  // which we need not check.
1454  const size_type first_right_of_diagonal_index
1455  = (Utilities::lower_bound (&cols->colnums[*rowstart_ptr+1],
1456  &cols->colnums[*(rowstart_ptr+1)],
1457  row)
1458  -
1459  &cols->colnums[0]);
1460 
1461  number s = 0;
1462  for (size_type j=(*rowstart_ptr)+1; j<first_right_of_diagonal_index; ++j)
1463  s += val[j] * dst(cols->colnums[j]);
1464 
1465  // divide by diagonal element
1466  *dst_ptr -= s * om;
1467  Assert(val[*rowstart_ptr]!= 0., ExcDivideByZero());
1468  *dst_ptr /= val[*rowstart_ptr];
1469  };
1470 
1471  rowstart_ptr = &cols->rowstart[0];
1472  dst_ptr = &dst(0);
1473  for (size_type row=0; row<n; ++row, ++rowstart_ptr, ++dst_ptr)
1474  *dst_ptr *= (2.-om)*val[*rowstart_ptr];
1475 
1476  // backward sweep
1477  rowstart_ptr = &cols->rowstart[n-1];
1478  dst_ptr = &dst(n-1);
1479  for (int row=n-1; row>=0; --row, --rowstart_ptr, --dst_ptr)
1480  {
1481  const size_type end_row = *(rowstart_ptr+1);
1482  const size_type first_right_of_diagonal_index
1483  = (Utilities::lower_bound (&cols->colnums[*rowstart_ptr+1],
1484  &cols->colnums[end_row],
1485  static_cast<size_type>(row)) -
1486  &cols->colnums[0]);
1487  number s = 0;
1488  for (size_type j=first_right_of_diagonal_index; j<end_row; ++j)
1489  s += val[j] * dst(cols->colnums[j]);
1490  *dst_ptr -= s * om;
1491  Assert(val[*rowstart_ptr]!= 0., ExcDivideByZero());
1492  *dst_ptr /= val[*rowstart_ptr];
1493  };
1494 }
1495 
1496 
1497 template <typename number>
1498 template <typename somenumber>
1499 void
1501  const Vector<somenumber> &src,
1502  const number om) const
1503 {
1504  Assert (cols != 0, ExcNotInitialized());
1505  Assert (val != 0, ExcNotInitialized());
1506 
1507  dst = src;
1508  SOR(dst,om);
1509 }
1510 
1511 
1512 template <typename number>
1513 template <typename somenumber>
1514 void
1516  const Vector<somenumber> &src,
1517  const number om) const
1518 {
1519  Assert (cols != 0, ExcNotInitialized());
1520  Assert (val != 0, ExcNotInitialized());
1521 
1522  dst = src;
1523  TSOR(dst,om);
1524 }
1525 
1526 
1527 template <typename number>
1528 template <typename somenumber>
1529 void
1531  const number om) const
1532 {
1533  Assert (cols != 0, ExcNotInitialized());
1534  Assert (val != 0, ExcNotInitialized());
1535  AssertDimension (m(), n());
1536  AssertDimension (dst.size(), n());
1537 
1538  for (size_type row=0; row<m(); ++row)
1539  {
1540  somenumber s = dst(row);
1541  for (size_type j=cols->rowstart[row]; j<cols->rowstart[row+1]; ++j)
1542  {
1543  const size_type col = cols->colnums[j];
1544  if (col < row)
1545  s -= val[j] * dst(col);
1546  }
1547 
1548  Assert(val[cols->rowstart[row]]!= 0., ExcDivideByZero());
1549  dst(row) = s * om / val[cols->rowstart[row]];
1550  }
1551 }
1552 
1553 
1554 template <typename number>
1555 template <typename somenumber>
1556 void
1558  const number om) const
1559 {
1560  Assert (cols != 0, ExcNotInitialized());
1561  Assert (val != 0, ExcNotInitialized());
1562  AssertDimension (m(), n());
1563  AssertDimension (dst.size(), n());
1564 
1565  size_type row=m()-1;
1566  while (true)
1567  {
1568  somenumber s = dst(row);
1569  for (size_type j=cols->rowstart[row]; j<cols->rowstart[row+1]; ++j)
1570  if (cols->colnums[j] > row)
1571  s -= val[j] * dst(cols->colnums[j]);
1572 
1573  Assert(val[cols->rowstart[row]]!= 0., ExcDivideByZero());
1574  dst(row) = s * om / val[cols->rowstart[row]];
1575 
1576  if (row == 0)
1577  break;
1578 
1579  --row;
1580  }
1581 }
1582 
1583 
1584 template <typename number>
1585 template <typename somenumber>
1586 void
1588  const std::vector<size_type> &permutation,
1589  const std::vector<size_type> &inverse_permutation,
1590  const number om) const
1591 {
1592  Assert (cols != 0, ExcNotInitialized());
1593  Assert (val != 0, ExcNotInitialized());
1594  AssertDimension (m(), n());
1595 
1596  Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
1597  Assert (m() == permutation.size(),
1598  ExcDimensionMismatch(m(), permutation.size()));
1599  Assert (m() == inverse_permutation.size(),
1600  ExcDimensionMismatch(m(), inverse_permutation.size()));
1601 
1602  for (size_type urow=0; urow<m(); ++urow)
1603  {
1604  const size_type row = permutation[urow];
1605  somenumber s = dst(row);
1606 
1607  for (size_type j=cols->rowstart[row]; j<cols->rowstart[row+1]; ++j)
1608  {
1609  const size_type col = cols->colnums[j];
1610  if (inverse_permutation[col] < urow)
1611  {
1612  s -= val[j] * dst(col);
1613  }
1614  }
1615 
1616  Assert(val[cols->rowstart[row]]!= 0., ExcDivideByZero());
1617  dst(row) = s * om / val[cols->rowstart[row]];
1618  }
1619 }
1620 
1621 
1622 template <typename number>
1623 template <typename somenumber>
1624 void
1626  const std::vector<size_type> &permutation,
1627  const std::vector<size_type> &inverse_permutation,
1628  const number om) const
1629 {
1630  Assert (cols != 0, ExcNotInitialized());
1631  Assert (val != 0, ExcNotInitialized());
1632  AssertDimension (m(), n());
1633 
1634  Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
1635  Assert (m() == permutation.size(),
1636  ExcDimensionMismatch(m(), permutation.size()));
1637  Assert (m() == inverse_permutation.size(),
1638  ExcDimensionMismatch(m(), inverse_permutation.size()));
1639 
1640  for (size_type urow=m(); urow != 0;)
1641  {
1642  --urow;
1643  const size_type row = permutation[urow];
1644  somenumber s = dst(row);
1645  for (size_type j=cols->rowstart[row]; j<cols->rowstart[row+1]; ++j)
1646  {
1647  const size_type col = cols->colnums[j];
1648  if (inverse_permutation[col] > urow)
1649  s -= val[j] * dst(col);
1650  }
1651 
1652  Assert(val[cols->rowstart[row]]!= 0., ExcDivideByZero());
1653  dst(row) = s * om / val[cols->rowstart[row]];
1654  }
1655 }
1656 
1657 
1658 
1659 template <typename number>
1660 template <typename somenumber>
1661 void
1663  const Vector<somenumber> &b,
1664  const number om) const
1665 {
1666  Assert (cols != 0, ExcNotInitialized());
1667  Assert (val != 0, ExcNotInitialized());
1668  AssertDimension (m(), n());
1669 
1670  Assert (m() == v.size(), ExcDimensionMismatch(m(),v.size()));
1671  Assert (m() == b.size(), ExcDimensionMismatch(m(),b.size()));
1672 
1674  typename VectorMemory<Vector<somenumber> >::Pointer w(mem);
1675  w->reinit(v);
1676 
1677  if (!v.all_zero())
1678  {
1679  vmult (*w, v);
1680  *w -= b;
1681  }
1682  else
1683  w->equ (-1.,b);
1684  precondition_Jacobi (*w, *w, om);
1685  v -= *w;
1686 }
1687 
1688 
1689 
1690 template <typename number>
1691 template <typename somenumber>
1692 void
1694  const Vector<somenumber> &b,
1695  const number om) const
1696 {
1697  Assert (cols != 0, ExcNotInitialized());
1698  Assert (val != 0, ExcNotInitialized());
1699  AssertDimension (m(), n());
1700  Assert (m() == v.size(), ExcDimensionMismatch(m(),v.size()));
1701  Assert (m() == b.size(), ExcDimensionMismatch(m(),b.size()));
1702 
1703  for (size_type row=0; row<m(); ++row)
1704  {
1705  somenumber s = b(row);
1706  for (size_type j=cols->rowstart[row]; j<cols->rowstart[row+1]; ++j)
1707  {
1708  s -= val[j] * v(cols->colnums[j]);
1709  }
1710  Assert(val[cols->rowstart[row]]!= 0., ExcDivideByZero());
1711  v(row) += s * om / val[cols->rowstart[row]];
1712  }
1713 }
1714 
1715 
1716 
1717 template <typename number>
1718 template <typename somenumber>
1719 void
1721  const Vector<somenumber> &b,
1722  const number om) const
1723 {
1724  Assert (cols != 0, ExcNotInitialized());
1725  Assert (val != 0, ExcNotInitialized());
1726  AssertDimension (m(), n());
1727  Assert (m() == v.size(), ExcDimensionMismatch(m(),v.size()));
1728  Assert (m() == b.size(), ExcDimensionMismatch(m(),b.size()));
1729 
1730  for (int row=m()-1; row>=0; --row)
1731  {
1732  somenumber s = b(row);
1733  for (size_type j=cols->rowstart[row]; j<cols->rowstart[row+1]; ++j)
1734  {
1735  s -= val[j] * v(cols->colnums[j]);
1736  }
1737  Assert(val[cols->rowstart[row]]!= 0., ExcDivideByZero());
1738  v(row) += s * om / val[cols->rowstart[row]];
1739  }
1740 }
1741 
1742 
1743 
1744 template <typename number>
1745 template <typename somenumber>
1746 void
1748  const Vector<somenumber> &b,
1749  const number om) const
1750 {
1751  SOR_step(v,b,om);
1752  TSOR_step(v,b,om);
1753 }
1754 
1755 
1756 
1757 template <typename number>
1758 template <typename somenumber>
1759 void
1761  const number om) const
1762 {
1763 //TODO: Is this called anywhere? If so, multiplication with om(2-om)D is missing
1764  Assert(false, ExcNotImplemented());
1765 
1766  Assert (cols != 0, ExcNotInitialized());
1767  Assert (val != 0, ExcNotInitialized());
1768  AssertDimension (m(), n());
1769  Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
1770 
1771  const size_type n = dst.size();
1772  size_type j;
1773  somenumber s;
1774 
1775  for (size_type i=0; i<n; i++)
1776  {
1777  s = 0.;
1778  for (j=cols->rowstart[i]; j<cols->rowstart[i+1] ; j++)
1779  {
1780  const size_type p = cols->colnums[j];
1782  {
1783  if (i>j) s += val[j] * dst(p);
1784  }
1785  }
1786  dst(i) -= s * om;
1787  Assert(val[cols->rowstart[i]]!= 0., ExcDivideByZero());
1788  dst(i) /= val[cols->rowstart[i]];
1789  }
1790 
1791  for (int i=n-1; i>=0; i--) // this time, i is signed, but always positive!
1792  {
1793  s = 0.;
1794  for (j=cols->rowstart[i]; j<cols->rowstart[i+1] ; j++)
1795  {
1796  const size_type p = cols->colnums[j];
1798  {
1799  if (static_cast<size_type>(i)<j) s += val[j] * dst(p);
1800  }
1801  }
1802  Assert(val[cols->rowstart[i]]!= 0., ExcDivideByZero());
1803  dst(i) -= s * om / val[cols->rowstart[i]];
1804  }
1805 }
1806 
1807 
1808 
1809 template <typename number>
1810 const SparsityPattern &
1812 {
1813  Assert (cols != 0, ExcNotInitialized());
1814  return *cols;
1815 }
1816 
1817 
1818 
1819 template <typename number>
1821  const unsigned int precision,
1822  const bool scientific,
1823  const unsigned int width_,
1824  const char *zero_string,
1825  const double denominator) const
1826 {
1827  Assert (cols != 0, ExcNotInitialized());
1828  Assert (val != 0, ExcNotInitialized());
1829 
1830  unsigned int width = width_;
1831 
1832  std::ios::fmtflags old_flags = out.flags();
1833  unsigned int old_precision = out.precision (precision);
1834 
1835  if (scientific)
1836  {
1837  out.setf (std::ios::scientific, std::ios::floatfield);
1838  if (!width)
1839  width = precision+7;
1840  }
1841  else
1842  {
1843  out.setf (std::ios::fixed, std::ios::floatfield);
1844  if (!width)
1845  width = precision+2;
1846  }
1847 
1848  for (size_type i=0; i<m(); ++i)
1849  {
1850  for (size_type j=0; j<n(); ++j)
1851  if ((*cols)(i,j) != SparsityPattern::invalid_entry)
1852  out << std::setw(width)
1853  << val[cols->operator()(i,j)] * denominator << ' ';
1854  else
1855  out << std::setw(width) << zero_string << ' ';
1856  out << std::endl;
1857  };
1858  AssertThrow (out, ExcIO());
1859 
1860  // reset output format
1861  out.precision(old_precision);
1862  out.flags (old_flags);
1863 }
1864 
1865 
1866 
1867 template <typename number>
1868 void SparseMatrix<number>::print_pattern (std::ostream &out,
1869  const double threshold) const
1870 {
1871  Assert (cols != 0, ExcNotInitialized());
1872  Assert (val != 0, ExcNotInitialized());
1873 
1874  for (size_type i=0; i<m(); ++i)
1875  {
1876  for (size_type j=0; j<n(); ++j)
1877  if ((*cols)(i,j) == SparsityPattern::invalid_entry)
1878  out << '.';
1879  else if (std::fabs(val[cols->operator()(i,j)]) > threshold)
1880  out << '*';
1881  else
1882  out << ':';
1883  out << std::endl;
1884  };
1885  AssertThrow (out, ExcIO());
1886 }
1887 
1888 
1889 
1890 template <typename number>
1891 void
1892 SparseMatrix<number>::block_write (std::ostream &out) const
1893 {
1894  AssertThrow (out, ExcIO());
1895 
1896  // first the simple objects,
1897  // bracketed in [...]
1898  out << '[' << max_len << "][";
1899  // then write out real data
1900  out.write (reinterpret_cast<const char *>(&val[0]),
1901  reinterpret_cast<const char *>(&val[max_len])
1902  - reinterpret_cast<const char *>(&val[0]));
1903  out << ']';
1904 
1905  AssertThrow (out, ExcIO());
1906 }
1907 
1908 
1909 
1910 template <typename number>
1911 void
1913 {
1914  AssertThrow (in, ExcIO());
1915 
1916  char c;
1917 
1918  // first read in simple data
1919  in >> c;
1920  AssertThrow (c == '[', ExcIO());
1921  in >> max_len;
1922 
1923  in >> c;
1924  AssertThrow (c == ']', ExcIO());
1925  in >> c;
1926  AssertThrow (c == '[', ExcIO());
1927 
1928  // reallocate space
1929  delete[] val;
1930  val = new number[max_len];
1931 
1932  // then read data
1933  in.read (reinterpret_cast<char *>(&val[0]),
1934  reinterpret_cast<char *>(&val[max_len])
1935  - reinterpret_cast<char *>(&val[0]));
1936  in >> c;
1937  AssertThrow (c == ']', ExcIO());
1938 }
1939 
1940 
1941 template <typename number>
1942 void
1943 SparseMatrix<number>::compress (::VectorOperation::values)
1944 {}
1945 
1946 
1947 template <typename number>
1948 std::size_t
1950 {
1951  return max_len*static_cast<std::size_t>(sizeof(number)) + sizeof(*this);
1952 }
1953 
1954 
1955 DEAL_II_NAMESPACE_CLOSE
1956 
1957 #endif
virtual void clear()
somenumber matrix_norm_square(const Vector< somenumber > &v) const
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:731
void precondition_Jacobi(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1.) const
static const number & conjugate(const number &x)
Definition: numbers.h:298
void Tvmult(OutVector &dst, const InVector &src) const
void vmult(OutVector &dst, const InVector &src) const
void PSOR(Vector< somenumber > &v, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const number om=1.) const
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:858
void Tvmult_add(OutVector &dst, const InVector &src) const
numbers::NumberTraits< number >::real_type real_type
void mmult(SparseMatrix< numberC > &C, const SparseMatrix< numberB > &B, const Vector< number > &V=Vector< number >(), const bool rebuild_sparsity_pattern=true) const
const Epetra_CrsMatrix & trilinos_matrix() const
void SSOR(Vector< somenumber > &v, const number omega=1.) const
size_type m() const
size_type get_row_length(const size_type row) const
::ExceptionBase & ExcMessage(std::string arg1)
void set(const size_type i, const size_type j, const number value)
size_type n() const
void block_read(std::istream &in)
size_type n_actually_nonzero_elements(const double threshold=0.) const
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 Tmmult(SparseMatrix< numberC > &C, const SparseMatrix< numberB > &B, const Vector< number > &V=Vector< number >(), const bool rebuild_sparsity_pattern=true) const
void SOR(Vector< somenumber > &v, const number om=1.) 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
bool is_finite(const double x)
bool all_zero() const
void print_pattern(std::ostream &out, const double threshold=0.) const
::ExceptionBase & ExcIO()
real_type l1_norm() const
void precondition_TSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
size_type n() const
SparseMatrix< number > & operator=(const SparseMatrix< number > &)
unsigned int global_dof_index
Definition: types.h:100
#define Assert(cond, exc)
Definition: exceptions.h:299
std::size_t memory_consumption() const
real_type linfty_norm() const
const SparsityPattern & get_sparsity_pattern() const
static bool equal(const T *p1, const T *p2)
void precondition_SSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1., const std::vector< std::size_t > &pos_right_of_diagonal=std::vector< std::size_t >()) const
size_type n_cols() const
unsigned int row_length(const size_type row) const
std::size_t size() const
void precondition_SOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
::ExceptionBase & ExcInvalidConstructorCall()
size_type n_nonzero_elements() const
size_type m() const
somenumber residual(Vector< somenumber > &dst, const Vector< somenumber > &x, const Vector< somenumber > &b) const
somenumber matrix_scalar_product(const Vector< somenumber > &u, const Vector< somenumber > &v) const
iterator begin()
void add(const size_type i, const size_type j, const number value)
void TPSOR(Vector< somenumber > &v, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const number om=1.) const
SparseMatrix< number > & copy_from(const SparseMatrix< somenumber > &source)
virtual void reinit(const SparsityPattern &sparsity)
void Jacobi_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
std::size_t * rowstart
void vmult_add(OutVector &dst, const InVector &src) const
types::global_dof_index size_type
::ExceptionBase & ExcNumberNotFinite()
void TSOR(Vector< somenumber > &v, const number om=1.) const
void SSOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
size_type n_rows() const
SmartPointer< const SparsityPattern, SparseMatrix< number > > cols
void block_write(std::ostream &out) const
::ExceptionBase & ExcNotImplemented()
::ExceptionBase & ExcNotInitialized()
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
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
::ExceptionBase & ExcScalarAssignmentOnlyForZeroValue()
void SOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
::ExceptionBase & ExcInternalError()
::ExceptionBase & ExcDivideByZero()
real_type frobenius_norm() const
VectorizedArray< Number > max(const ::VectorizedArray< Number > &x, const ::VectorizedArray< Number > &y)
unsigned int row_length(const size_type row) const
std::size_t max_len
void compress(::VectorOperation::values)
real_type linfty_norm() const
static const size_type invalid_entry
void TSOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const