deal.II/include/deal.II/lac/sparse_matrix.templates.h

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