]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Improve performance in the typical use case of sparse matrix iterator by basing the...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 6 Feb 2013 10:44:40 +0000 (10:44 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 6 Feb 2013 10:44:40 +0000 (10:44 +0000)
git-svn-id: https://svn.dealii.org/trunk@28246 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/include/deal.II/lac/precondition.h
deal.II/include/deal.II/lac/sparse_matrix.h
deal.II/include/deal.II/lac/sparse_matrix.templates.h
deal.II/include/deal.II/lac/sparse_matrix_ez.h
deal.II/include/deal.II/lac/sparse_matrix_ez.templates.h
deal.II/source/lac/sparse_matrix.inst.in
deal.II/source/lac/sparse_matrix_ez.inst.in

index 497a0c153c1fcca993ee7ec3119068221cbe0726..a56e5554435a019cffc9ce93ea313f8887673734 100644 (file)
@@ -17,6 +17,7 @@
 #include <deal.II/base/config.h>
 #include <deal.II/base/smartpointer.h>
 #include <deal.II/base/utilities.h>
+#include <deal.II/base/parallel.h>
 #include <deal.II/base/template_constraints.h>
 #include <deal.II/lac/tridiagonal_matrix.h>
 #include <deal.II/lac/vector_memory.h>
@@ -25,6 +26,13 @@ DEAL_II_NAMESPACE_OPEN
 
 template <typename number> class Vector;
 template <typename number> class SparseMatrix;
+namespace parallel
+{
+  namespace distributed
+  {
+    template <typename number> class Vector;
+  }
+}
 
 /*! @addtogroup Preconditioners
  *@{
@@ -619,7 +627,7 @@ private:
    * row where the first position after
    * the diagonal is located.
    */
-  std::vector<unsigned int> pos_right_of_diagonal;
+  std::vector<std::size_t> pos_right_of_diagonal;
 };
 
 
@@ -1311,36 +1319,28 @@ PreconditionSSOR<MATRIX>::initialize (const MATRIX &rA,
 {
   this->PreconditionRelaxation<MATRIX>::initialize (rA, parameters);
 
-  // in case we have a SparseMatrix class,
-  // we can extract information about the
-  // diagonal.
+  // in case we have a SparseMatrix class, we can extract information about
+  // the diagonal.
   const SparseMatrix<typename MATRIX::value_type> *mat =
     dynamic_cast<const SparseMatrix<typename MATRIX::value_type> *>(&*this->A);
 
-  // calculate the positions first after
-  // the diagonal.
+  // calculate the positions first after the diagonal.
   if (mat != 0)
     {
-      const std::size_t   *rowstart_ptr =
-        mat->get_sparsity_pattern().get_rowstart_indices();
-      const unsigned int *const colnums =
-        mat->get_sparsity_pattern().get_column_numbers();
       const unsigned int n = this->A->n();
-      pos_right_of_diagonal.resize(n);
-      for (unsigned int row=0; row<n; ++row, ++rowstart_ptr)
+      pos_right_of_diagonal.resize(n, static_cast<std::size_t>(-1));
+      for (unsigned int row=0; row<n; ++row)
         {
-          // find the first element in this line
-          // which is on the right of the diagonal.
-          // we need to precondition with the
-          // elements on the left only.
-          // note: the first entry in each
-          // line denotes the diagonal element,
-          // which we need not check.
-          pos_right_of_diagonal[row] =
-            Utilities::lower_bound (&colnums[*rowstart_ptr+1],
-                                    &colnums[*(rowstart_ptr+1)],
-                                    row)
-            - colnums;
+          // find the first element in this line which is on the right of the
+          // diagonal.  we need to precondition with the elements on the left
+          // only. note: the first entry in each line denotes the diagonal
+          // element, which we need not check.
+          typename SparseMatrix<typename MATRIX::value_type>::const_iterator
+            it = mat->begin(row)+1;
+          for ( ; it < mat->end(row); ++it)
+            if (it->column() > row)
+              break;
+          pos_right_of_diagonal[row] = it - mat->begin();
         }
     }
 }
@@ -1648,7 +1648,8 @@ PreconditionChebyshev<MATRIX,VECTOR>::initialize (const MATRIX &matrix,
 
         // need at least two iterations to have
         // maximum and minimum eigenvalue
-        if (it > data.eig_cg_n_iterations || (it > 2 &&
+        if (res == 0. ||
+            it > data.eig_cg_n_iterations || (it > 2 &&
                                               res < data.eig_cg_residual))
           break;
 
@@ -1662,16 +1663,24 @@ PreconditionChebyshev<MATRIX,VECTOR>::initialize (const MATRIX &matrix,
         offdiagonal.push_back(std::sqrt(beta)/alpha);
       }
 
-    TridiagonalMatrix<double> T(diagonal.size(), true);
-    for (unsigned int i=0; i<diagonal.size(); ++i)
+    if (diagonal.size() == 0)
+      min_eigenvalue = max_eigenvalue = 1.;
+    else
       {
-        T(i,i) = diagonal[i];
-        if (i< diagonal.size()-1)
-          T(i,i+1) = offdiagonal[i];
+        TridiagonalMatrix<double> T(diagonal.size(), true);
+        for (unsigned int i=0; i<diagonal.size(); ++i)
+          {
+            T(i,i) = diagonal[i];
+            if (i< diagonal.size()-1)
+              T(i,i+1) = offdiagonal[i];
+          }
+        T.compute_eigenvalues();
+        min_eigenvalue = T.eigenvalue(0);
+        if (diagonal.size() > 1)
+          max_eigenvalue = T.eigenvalue(T.n()-1);
+        else
+          max_eigenvalue = min_eigenvalue;
       }
-    T.compute_eigenvalues();
-    min_eigenvalue = T.eigenvalue(0);
-    max_eigenvalue = T.eigenvalue(T.n()-1);
   }
 
   // include a safety factor since the CG
@@ -1687,6 +1696,161 @@ PreconditionChebyshev<MATRIX,VECTOR>::initialize (const MATRIX &matrix,
 
 
 
+namespace internal
+{
+  namespace PreconditionChebyshev
+  {
+    // for deal.II vectors, perform updates for Chebyshev preconditioner all
+    // at once to reduce memory transfer. Here, we select between general
+    // vectors and deal.II vectors where we expand the loop over the (local)
+    // size of the vector
+
+    // generic part for non-deal.II vectors
+    template <typename VECTOR>
+    inline
+    void
+    vector_updates (const VECTOR &src,
+                    const VECTOR &matrix_diagonal_inverse,
+                    const bool    start_zero,
+                    const double  factor1,
+                    const double  factor2,
+                    VECTOR &update1,
+                    VECTOR &update2,
+                    VECTOR &dst)
+    {
+      if (start_zero)
+        {
+          dst.equ (factor2, src);
+          dst.scale (matrix_diagonal_inverse);
+          update1.equ(-1.,dst);
+        }
+      else
+        {
+          update2 -= src;
+          update2.scale (matrix_diagonal_inverse);
+          if (factor1 == 0.)
+            update1.equ(factor2, update2);
+          else
+            update1.sadd(factor1, factor2, update2);
+          dst -= update1;
+        }
+    }
+
+    // worker loop for deal.II vectors
+    template <typename Number>
+    struct VectorUpdatesRange : public parallel::ParallelForInteger
+    {
+      VectorUpdatesRange (const size_t  size,
+                          const Number *src,
+                          const Number *matrix_diagonal_inverse,
+                          const bool    start_zero,
+                          const Number  factor1,
+                          const Number  factor2,
+                          Number       *update1,
+                          Number       *update2,
+                          Number       *dst)
+        :
+        src (src),
+        matrix_diagonal_inverse (matrix_diagonal_inverse),
+        start_zero (start_zero),
+        factor1 (factor1),
+        factor2 (factor2),
+        update1 (update1),
+        update2 (update2),
+        dst (dst)
+      {
+        if (size < internal::Vector::minimum_parallel_grain_size)
+          apply_to_subrange (0, size);
+        else
+          apply_parallel (0, size,
+                          internal::Vector::minimum_parallel_grain_size);
+      }
+
+      ~VectorUpdatesRange()
+      {}
+
+      virtual void
+      apply_to_subrange (const size_t begin,
+                         const size_t end) const
+      {
+        if (factor1 == Number())
+          {
+            if (start_zero)
+              for (unsigned int i=begin; i<end; ++i)
+                {
+                  dst[i] = factor2 * src[i] * matrix_diagonal_inverse[i];
+                  update1[i] = -dst[i];
+                }
+            else
+              for (unsigned int i=begin; i<end; ++i)
+                {
+                  update1[i] = ((update2[i]-src[i]) *
+                                factor2*matrix_diagonal_inverse[i]);
+                  dst[i] -= update1[i];
+                }
+          }
+        else
+          for (unsigned int i=begin; i<end; ++i)
+            {
+              const Number update2i = ((update2[i] - src[i]) *
+                                       matrix_diagonal_inverse[i]);
+              update1[i] = factor1 * update1[i] + factor2 * update2i;
+              dst[i] -= update1[i];
+            }
+      }
+
+      const Number *src;
+      const Number *matrix_diagonal_inverse;
+      const bool start_zero;
+      const Number factor1;
+      const Number factor2;
+      mutable Number *update1;
+      mutable Number *update2;
+      mutable Number *dst;
+    };
+
+    // selection for deal.II vector
+    template <typename Number>
+    inline
+    void
+    vector_updates (const ::dealii::Vector<Number> &src,
+                    const ::dealii::Vector<Number> &matrix_diagonal_inverse,
+                    const bool    start_zero,
+                    const double  factor1,
+                    const double  factor2,
+                    ::dealii::Vector<Number> &update1,
+                    ::dealii::Vector<Number> &update2,
+                    ::dealii::Vector<Number> &dst)
+    {
+      VectorUpdatesRange<Number>(src.size(), src.begin(),
+                                 matrix_diagonal_inverse.begin(),
+                                 start_zero, factor1, factor2,
+                                 update1.begin(), update2.begin(), dst.begin());
+    }
+
+    // selection for parallel deal.II vector
+    template <typename Number>
+    inline
+    void
+    vector_updates (const parallel::distributed::Vector<Number> &src,
+                    const parallel::distributed::Vector<Number> &matrix_diagonal_inverse,
+                    const bool    start_zero,
+                    const double  factor1,
+                    const double  factor2,
+                    parallel::distributed::Vector<Number> &update1,
+                    parallel::distributed::Vector<Number> &update2,
+                    parallel::distributed::Vector<Number> &dst)
+    {
+      VectorUpdatesRange<Number>(src.local_size(), src.begin(),
+                                 matrix_diagonal_inverse.begin(),
+                                 start_zero, factor1, factor2,
+                                 update1.begin(), update2.begin(), dst.begin());
+    }
+  }
+}
+
+
+
 template <class MATRIX, class VECTOR>
 inline
 void
@@ -1697,29 +1861,25 @@ PreconditionChebyshev<MATRIX,VECTOR>::vmult (VECTOR &dst,
   double rhok  = delta / theta,  sigma = theta / delta;
   if (data.nonzero_starting && !dst.all_zero())
     {
-      matrix_ptr->vmult (update1, dst);
-      update1 -= src;
-      update1 /= theta;
-      update1.scale (data.matrix_diagonal_inverse);
-      dst -= update1;
+      matrix_ptr->vmult (update2, dst);
+      internal::PreconditionChebyshev::vector_updates
+      (src, data.matrix_diagonal_inverse, false, 0., 1./theta, update1,
+       update2, dst);
     }
   else
-    {
-      dst.equ (1./theta, src);
-      dst.scale (data.matrix_diagonal_inverse);
-      update1.equ(-1.,dst);
-    }
+    internal::PreconditionChebyshev::vector_updates
+    (src, data.matrix_diagonal_inverse, true, 0., 1./theta, update1,
+     update2, dst);
 
   for (unsigned int k=0; k<data.degree; ++k)
     {
       matrix_ptr->vmult (update2, dst);
-      update2 -= src;
-      update2.scale (data.matrix_diagonal_inverse);
       const double rhokp = 1./(2.*sigma-rhok);
       const double factor1 = rhokp * rhok, factor2 = 2.*rhokp/delta;
       rhok = rhokp;
-      update1.sadd (factor1, factor2, update2);
-      dst -= update1;
+      internal::PreconditionChebyshev::vector_updates
+      (src, data.matrix_diagonal_inverse, false, factor1, factor2, update1,
+       update2, dst);
     }
 }
 
@@ -1735,29 +1895,25 @@ PreconditionChebyshev<MATRIX,VECTOR>::Tvmult (VECTOR &dst,
   double rhok  = delta / theta,  sigma = theta / delta;
   if (data.nonzero_starting && !dst.all_zero())
     {
-      matrix_ptr->Tvmult (update1, dst);
-      update1 -= src;
-      update1 /= theta;
-      update1.scale (data.matrix_diagonal_inverse);
-      dst -= update1;
+      matrix_ptr->Tvmult (update2, dst);
+      internal::PreconditionChebyshev::vector_updates
+      (src, data.matrix_diagonal_inverse, false, 0., 1./theta, update1,
+       update2, dst);
     }
   else
-    {
-      dst.equ (1./theta, src);
-      dst.scale (data.matrix_diagonal_inverse);
-      update1.equ(-1.,dst);
-    }
+    internal::PreconditionChebyshev::vector_updates
+    (src, data.matrix_diagonal_inverse, true, 0., 1./theta, update1,
+     update2, dst);
 
   for (unsigned int k=0; k<data.degree-1; ++k)
     {
       matrix_ptr->Tvmult (update2, dst);
-      update2 -= src;
-      update2.scale (data.matrix_diagonal_inverse);
       const double rhokp = 1./(2.*sigma-rhok);
       const double factor1 = rhokp * rhok, factor2 = 2.*rhokp/delta;
       rhok = rhokp;
-      update1.sadd (factor1, factor2, update2);
-      dst -= update1;
+      internal::PreconditionChebyshev::vector_updates
+      (src, data.matrix_diagonal_inverse, false, factor1, factor2, update1,
+       update2, dst);
     }
 }
 
index 0a2d7e08ed983cf3db0e2b48a5789ba17f299e2e..e061abc4f701e921738bbc293f2d5862784b8e61 100644 (file)
@@ -98,6 +98,12 @@ namespace SparseMatrixIterators
               const unsigned int  row,
               const unsigned int  index);
 
+    /**
+     * Constructor.
+     */
+    Accessor (MatrixType         *matrix,
+              const std::size_t   index_with_matrix);
+
     /**
      * Constructor. Construct the end accessor for the given matrix.
      */
@@ -234,6 +240,12 @@ namespace SparseMatrixIterators
               const unsigned int  row,
               const unsigned int  index);
 
+    /**
+     * Constructor.
+     */
+    Accessor (MatrixType         *matrix,
+              const std::size_t   index);
+
     /**
      * Constructor. Construct the end accessor for the given matrix.
      */
@@ -318,6 +330,13 @@ namespace SparseMatrixIterators
               const unsigned int row,
               const unsigned int index);
 
+    /**
+     * Constructor. Create an iterator into the matrix @p matrix for the given
+     * index in the complete matrix (counting from the zeroth entry).
+     */
+    Iterator (MatrixType        *matrix,
+              const std::size_t  index_within_matrix);
+
     /**
      * Constructor. Create the end iterator for the given matrix.
      */
@@ -1198,7 +1217,7 @@ public:
   void precondition_SSOR (Vector<somenumber>             &dst,
                           const Vector<somenumber>       &src,
                           const number                    omega = 1.,
-                          const std::vector<unsigned int> &pos_right_of_diagonal=std::vector<unsigned int>()) const;
+                          const std::vector<std::size_t> &pos_right_of_diagonal=std::vector<std::size_t>()) const;
 
   /**
    * Apply SOR preconditioning matrix to <tt>src</tt>.
@@ -1850,7 +1869,7 @@ inline
 number SparseMatrix<number>::diag_element (const unsigned int i) const
 {
   Assert (cols != 0, ExcNotInitialized());
-  Assert (cols->optimize_diagonal(),  ExcNotQuadratic());
+  Assert (m() == n(),  ExcNotQuadratic());
   Assert (i<m(), ExcInvalidIndex1(i));
 
   // Use that the first element in each row of a quadratic matrix is the main
@@ -1865,7 +1884,7 @@ inline
 number &SparseMatrix<number>::diag_element (const unsigned int i)
 {
   Assert (cols != 0, ExcNotInitialized());
-  Assert (cols->optimize_diagonal(),  ExcNotQuadratic());
+  Assert (m() == n(),  ExcNotQuadratic());
   Assert (i<m(), ExcInvalidIndex1(i));
 
   // Use that the first element in each row of a quadratic matrix is the main
@@ -1972,6 +1991,19 @@ namespace SparseMatrixIterators
 
 
 
+  template <typename number>
+  inline
+  Accessor<number,true>::
+  Accessor (const MatrixType   *matrix,
+            const std::size_t   index)
+    :
+    SparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern(),
+                                        index),
+    matrix (matrix)
+  {}
+
+
+
   template <typename number>
   inline
   Accessor<number,true>::
@@ -1999,7 +2031,8 @@ namespace SparseMatrixIterators
   number
   Accessor<number, true>::value () const
   {
-    return matrix->nth_entry_in_row(a_row, a_index);
+    AssertIndexRange(a_index, matrix->n_nonzero_elements());
+    return matrix->val[a_index];
   }
 
 
@@ -2028,8 +2061,8 @@ namespace SparseMatrixIterators
   inline
   Accessor<number, false>::Reference::operator number() const
   {
-    return accessor->matrix->nth_entry_in_row(accessor->a_row,
-                                              accessor->a_index);
+    AssertIndexRange(accessor->a_index, accessor->matrix->n_nonzero_elements());
+    return accessor->matrix->val[accessor->a_index];
   }
 
 
@@ -2039,8 +2072,8 @@ namespace SparseMatrixIterators
   const typename Accessor<number, false>::Reference &
   Accessor<number, false>::Reference::operator = (const number n) const
   {
-//TODO: one could optimize this by not going again through the mapping from row/col index to global index
-    accessor->matrix->set (accessor->row(), accessor->column(), n);
+    AssertIndexRange(accessor->a_index, accessor->matrix->n_nonzero_elements());
+    accessor->matrix->val[accessor->a_index] = n;
     return *this;
   }
 
@@ -2051,9 +2084,8 @@ namespace SparseMatrixIterators
   const typename Accessor<number, false>::Reference &
   Accessor<number, false>::Reference::operator += (const number n) const
   {
-//TODO: one could optimize this by not going again through the mapping from row/col index to global index
-    accessor->matrix->set (accessor->row(), accessor->column(),
-                           static_cast<number>(*this) + n);
+    AssertIndexRange(accessor->a_index, accessor->matrix->n_nonzero_elements());
+    accessor->matrix->val[accessor->a_index] += n;
     return *this;
   }
 
@@ -2064,9 +2096,8 @@ namespace SparseMatrixIterators
   const typename Accessor<number, false>::Reference &
   Accessor<number, false>::Reference::operator -= (const number n) const
   {
-//TODO: one could optimize this by not going again through the mapping from row/col index to global index
-    accessor->matrix->set (accessor->row(), accessor->column(),
-                           static_cast<number>(*this) - n);
+    AssertIndexRange(accessor->a_index, accessor->matrix->n_nonzero_elements());
+    accessor->matrix->val[accessor->a_index] -= n;
     return *this;
   }
 
@@ -2077,9 +2108,8 @@ namespace SparseMatrixIterators
   const typename Accessor<number, false>::Reference &
   Accessor<number, false>::Reference::operator *= (const number n) const
   {
-//TODO: one could optimize this by not going again through the mapping from row/col index to global index
-    accessor->matrix->set (accessor->row(), accessor->column(),
-                           static_cast<number>(*this)*n);
+    AssertIndexRange(accessor->a_index, accessor->matrix->n_nonzero_elements());
+    accessor->matrix->val[accessor->a_index] *= n;
     return *this;
   }
 
@@ -2090,9 +2120,8 @@ namespace SparseMatrixIterators
   const typename Accessor<number, false>::Reference &
   Accessor<number, false>::Reference::operator /= (const number n) const
   {
-//TODO: one could optimize this by not going again through the mapping from row/col index to global index
-    accessor->matrix->set (accessor->row(), accessor->column(),
-                           static_cast<number>(*this)/n);
+    AssertIndexRange(accessor->a_index, accessor->matrix->n_nonzero_elements());
+    accessor->matrix->val[accessor->a_index] /= n;
     return *this;
   }
 
@@ -2112,6 +2141,19 @@ namespace SparseMatrixIterators
 
 
 
+  template <typename number>
+  inline
+  Accessor<number,false>::
+  Accessor (MatrixType         *matrix,
+            const std::size_t   index)
+    :
+    SparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern(),
+                                        index),
+    matrix (matrix)
+  {}
+
+
+
   template <typename number>
   inline
   Accessor<number,false>::
@@ -2156,6 +2198,17 @@ namespace SparseMatrixIterators
 
 
 
+  template <typename number, bool Constness>
+  inline
+  Iterator<number, Constness>::
+  Iterator (MatrixType        *matrix,
+            const std::size_t  index)
+    :
+    accessor(matrix, index)
+  {}
+
+
+
   template <typename number, bool Constness>
   inline
   Iterator<number, Constness>::
@@ -2269,21 +2322,7 @@ namespace SparseMatrixIterators
 
     const SparsityPattern &sparsity = accessor.get_matrix().get_sparsity_pattern();
 
-    const unsigned int this_position
-      = (*this != (*this)->get_matrix().end()
-         ?
-         sparsity.get_rowstart_indices()[(*this)->row()] + (*this)->index()
-         :
-         sparsity.get_rowstart_indices()[sparsity.n_rows()]);
-
-    const unsigned int other_position
-      = (other != (*this)->get_matrix().end()
-         ?
-         sparsity.get_rowstart_indices()[other->row()] + other->index()
-         :
-         sparsity.get_rowstart_indices()[sparsity.n_rows()]);
-
-    return (this_position - other_position);
+    return (*this)->a_index - other->a_index;
   }
 
 
@@ -2310,14 +2349,7 @@ inline
 typename SparseMatrix<number>::const_iterator
 SparseMatrix<number>::begin () const
 {
-  // search for the first line with a nonzero number of entries
-  for (unsigned int r=0; r<m(); ++r)
-    if (cols->row_length(r) > 0)
-      return const_iterator(this, r, 0);
-
-  // alright, this matrix is completely empty. that's strange but ok. simply
-  // return the end() iterator
-  return end();
+  return const_iterator(this, 0);
 }
 
 
@@ -2335,14 +2367,7 @@ inline
 typename SparseMatrix<number>::iterator
 SparseMatrix<number>::begin ()
 {
-  // search for the first line with a nonzero number of entries
-  for (unsigned int r=0; r<m(); ++r)
-    if (cols->row_length(r) > 0)
-      return iterator(this, r, 0);
-
-  // alright, this matrix is completely empty. that's strange but ok. simply
-  // return the end() iterator
-  return end();
+  return iterator (this, 0);
 }
 
 
@@ -2351,7 +2376,7 @@ inline
 typename SparseMatrix<number>::iterator
 SparseMatrix<number>::end ()
 {
-  return iterator(this, m(), 0);
+  return iterator(this, cols->rowstart[cols->rows]);
 }
 
 
@@ -2362,10 +2387,7 @@ SparseMatrix<number>::begin (const unsigned int r) const
 {
   Assert (r<m(), ExcIndexRange(r,0,m()));
 
-  if (cols->row_length(r) > 0)
-    return const_iterator(this, r, 0);
-  else
-    return end (r);
+  return const_iterator(this, cols->rowstart[r]);
 }
 
 
@@ -2377,14 +2399,7 @@ SparseMatrix<number>::end (const unsigned int r) const
 {
   Assert (r<m(), ExcIndexRange(r,0,m()));
 
-  // place the iterator on the first entry past this line, or at the end of
-  // the matrix
-  for (unsigned int i=r+1; i<m(); ++i)
-    if (cols->row_length(i) > 0)
-      return const_iterator(this, i, 0);
-
-  // if there is no such line, then take the end iterator of the matrix
-  return end();
+  return const_iterator(this, cols->rowstart[r+1]);
 }
 
 
@@ -2396,10 +2411,7 @@ SparseMatrix<number>::begin (const unsigned int r)
 {
   Assert (r<m(), ExcIndexRange(r,0,m()));
 
-  if (cols->row_length(r) > 0)
-    return iterator(this, r, 0);
-  else
-    return end (r);
+  return iterator(this, cols->rowstart[r]);
 }
 
 
@@ -2411,14 +2423,7 @@ SparseMatrix<number>::end (const unsigned int r)
 {
   Assert (r<m(), ExcIndexRange(r,0,m()));
 
-  // place the iterator on the first entry past this line, or at the end of
-  // the matrix
-  for (unsigned int i=r+1; i<m(); ++i)
-    if (cols->row_length(i) > 0)
-      return iterator(this, i, 0);
-
-  // if there is no such line, then take the end iterator of the matrix
-  return end();
+  return iterator(this, cols->rowstart[r+1]);
 }
 
 
index 8b929f9e4519a27d71543e90608105b05b8b7a5d..30f7c18a67b7b8c4214d31294da3339ac4f7f8b5 100644 (file)
@@ -289,7 +289,7 @@ SparseMatrix<number>::symmetrize ()
     {
       // first skip diagonal entry
       number             *val_ptr = &val[cols->rowstart[row]];
-      if (cols->optimize_diagonal())
+      if (m() == n())
         ++val_ptr;
       const unsigned int *colnum_ptr = &cols->colnums[cols->rowstart[row]+1];
       const number       *const val_end_of_row = &val[cols->rowstart[row+1]];
@@ -450,11 +450,11 @@ SparseMatrix<number>::add (const unsigned int  row,
 #endif
 
       const unsigned int *this_cols =
-        &cols->get_column_numbers()[cols->get_rowstart_indices()[row]];
+        &cols->colnums[cols->rowstart[row]];
       const unsigned int row_length_1 = cols->row_length(row)-1;
-      number *val_ptr = &val[cols->get_rowstart_indices()[row]];
+      number *val_ptr = &val[cols->rowstart[row]];
 
-      if (cols->optimize_diagonal() == true)
+      if (m() == n())
         {
 
           // find diagonal and add it if found
@@ -521,9 +521,9 @@ SparseMatrix<number>::add (const unsigned int  row,
   // unsorted case: first, search all the
   // indices to find out which values we
   // actually need to add.
-  const unsigned int *const my_cols = cols->get_column_numbers();
-  unsigned int index = cols->get_rowstart_indices()[row];
-  const unsigned int next_row_index = cols->get_rowstart_indices()[row+1];
+  const unsigned int *const my_cols = cols->colnums;
+  unsigned int index = cols->rowstart[row];
+  const unsigned int next_row_index = cols->rowstart[row+1];
 
   for (unsigned int j=0; j<n_cols; ++j)
     {
@@ -580,9 +580,9 @@ SparseMatrix<number>::set (const unsigned int  row,
   // First, search all the indices to find
   // out which values we actually need to
   // set.
-  const unsigned int *my_cols = cols->get_column_numbers();
-  std::size_t index = cols->get_rowstart_indices()[row], next_index = index;
-  const std::size_t next_row_index = cols->get_rowstart_indices()[row+1];
+  const unsigned int *my_cols = cols->colnums;
+  std::size_t index = cols->rowstart[row], next_index = index;
+  const std::size_t next_row_index = cols->rowstart[row+1];
 
   if (elide_zero_values == true)
     {
@@ -902,10 +902,8 @@ SparseMatrix<number>::mmult (SparseMatrix<numberC>       &C,
   // clear previous content of C
   if  (rebuild_sparsity_C == true)
     {
-      // we are about to change the sparsity
-      // pattern of C. this can not work if
-      // either A or B use the same sparsity
-      // pattern
+      // we are about to change the sparsity pattern of C. this can not work
+      // if either A or B use the same sparsity pattern
       Assert (&C.get_sparsity_pattern() != &this->get_sparsity_pattern(),
               ExcMessage ("Can't use the same sparsity pattern for "
                           "different matrices if it is to be rebuilt."));
@@ -913,44 +911,35 @@ SparseMatrix<number>::mmult (SparseMatrix<numberC>       &C,
               ExcMessage ("Can't use the same sparsity pattern for "
                           "different matrices if it is to be rebuilt."));
 
-      // need to change the sparsity pattern of
-      // C, so cast away const-ness.
+      // need to change the sparsity pattern of C, so cast away const-ness.
       SparsityPattern &sp_C =
         *(const_cast<SparsityPattern *>(&C.get_sparsity_pattern()));
       C.clear();
       sp_C.reinit (0,0,0);
 
-      // create a sparsity pattern for the
-      // matrix. we will go through all the
-      // rows in the matrix A, and for each
-      // column in a row we add the whole row
-      // of matrix B with that row number. This
-      // means that we will insert a lot of
-      // entries to each row, which is best
-      // handled by the
+      // create a sparsity pattern for the matrix. we will go through all the
+      // rows in the matrix A, and for each column in a row we add the whole
+      // row of matrix B with that row number. This means that we will insert
+      // a lot of entries to each row, which is best handled by the
       // CompressedSimpleSparsityPattern class.
       {
         CompressedSimpleSparsityPattern csp (m(), B.n());
         for (unsigned int i = 0; i < csp.n_rows(); ++i)
           {
-            const unsigned int *rows =
-              &sp_A.get_column_numbers()[sp_A.get_rowstart_indices()[i]];
+            const unsigned int *rows = &sp_A.colnums[sp_A.rowstart[i]];
             const unsigned int *const end_rows =
-              &sp_A.get_column_numbers()[sp_A.get_rowstart_indices()[i+1]];
+              &sp_A.colnums[sp_A.rowstart[i+1]];
             for (; rows != end_rows; ++rows)
               {
                 const unsigned int col = *rows;
                 unsigned int *new_cols = const_cast<unsigned int *>
-                                         (&sp_B.get_column_numbers()
-                                          [sp_B.get_rowstart_indices()[col]]);
+                                         (&sp_B.colnums[sp_B.rowstart[col]]);
                 unsigned int *end_new_cols = const_cast<unsigned int *>
-                                             (&sp_B.get_column_numbers()
-                                              [sp_B.get_rowstart_indices()[col+1]]);
+                                             (&sp_B.colnums[sp_B.rowstart[col+1]]);
 
-                // if B has a diagonal, need to add that
-                // manually. this way, we maintain
-                // sortedness.
-                if (sp_B.optimize_diagonal() == true)
+                // if B has a diagonal, need to add that manually. this way,
+                // we maintain sortedness.
+                if (sp_B.n_rows() == sp_B.n_cols())
                   {
                     ++new_cols;
                     csp.add(i, col);
@@ -977,46 +966,37 @@ SparseMatrix<number>::mmult (SparseMatrix<numberC>       &C,
     max_n_cols_B = std::max (max_n_cols_B, sp_B.row_length(i));
   std::vector<numberC> new_entries(max_n_cols_B);
 
-  // now compute the actual entries: a
-  // matrix-matrix product involves three
-  // nested loops. One over the rows of A,
-  // for each row we then loop over all the
-  // columns, and then we need to multiply
-  // each element with all the elements in
-  // that row in B.
+  // now compute the actual entries: a matrix-matrix product involves three
+  // nested loops. One over the rows of A, for each row we then loop over all
+  // the columns, and then we need to multiply each element with all the
+  // elements in that row in B.
   for (unsigned int i=0; i<C.m(); ++i)
     {
-      const unsigned int *rows =
-        &sp_A.get_column_numbers()[sp_A.get_rowstart_indices()[i]];
-      const unsigned int *const end_rows =
-        &sp_A.get_column_numbers()[sp_A.get_rowstart_indices()[i+1]];
+      const unsigned int *rows = &sp_A.colnums[sp_A.rowstart[i]];
+      const unsigned int *const end_rows = &sp_A.colnums[sp_A.rowstart[i+1]];
       for (; rows != end_rows; ++rows)
         {
-          const double A_val = val[rows-&sp_A.get_column_numbers()[sp_A.get_rowstart_indices()[0]]];
+          const double A_val = val[rows-&sp_A.colnums[sp_A.rowstart[0]]];
           const unsigned int col = *rows;
           const unsigned int *new_cols =
-            (&sp_B.get_column_numbers()[sp_B.get_rowstart_indices()[col]]);
+            (&sp_B.colnums[sp_B.rowstart[col]]);
 
           // special treatment for diagonal
-          if (sp_B.optimize_diagonal())
+          if (sp_B.n_rows() == sp_B.n_cols())
             {
               C.add (i, *new_cols, A_val *
-                     B.val[new_cols-&sp_B.get_column_numbers()
-                           [sp_B.get_rowstart_indices()[0]]] *
+                     B.val[new_cols-&sp_B.colnums[sp_B.rowstart[0]]] *
                      (use_vector ? V(col) : 1));
               ++new_cols;
             }
 
-          // now the innermost loop that goes over
-          // all the elements in row 'col' of
-          // matrix B. Cache the elements, and then
-          // write them into C at once
+          // now the innermost loop that goes over all the elements in row
+          // 'col' of matrix B. Cache the elements, and then write them into C
+          // at once
           numberC *new_ptr = &new_entries[0];
           const numberB *B_val_ptr =
-            &B.val[new_cols-&sp_B.get_column_numbers()[sp_B.get_rowstart_indices()[0]]];
-          const numberB *const end_cols =
-            &B.val[&sp_B.get_column_numbers()[sp_B.get_rowstart_indices()[col+1]]-
-                   &sp_B.get_column_numbers()[sp_B.get_rowstart_indices()[0]]];
+            &B.val[new_cols-&sp_B.colnums[sp_B.rowstart[0]]];
+          const numberB *const end_cols = &B.val[sp_B.rowstart[col+1]];
           for (; B_val_ptr != end_cols; ++B_val_ptr)
             *new_ptr++ = A_val * *B_val_ptr * (use_vector ? V(col) : 1);
 
@@ -1049,10 +1029,8 @@ SparseMatrix<number>::Tmmult (SparseMatrix<numberC>       &C,
   // clear previous content of C
   if  (rebuild_sparsity_C == true)
     {
-      // we are about to change the sparsity
-      // pattern of C. this can not work if
-      // either A or B use the same sparsity
-      // pattern
+      // we are about to change the sparsity pattern of C. this can not work
+      // if either A or B use the same sparsity pattern
       Assert (&C.get_sparsity_pattern() != &this->get_sparsity_pattern(),
               ExcMessage ("Can't use the same sparsity pattern for "
                           "different matrices if it is to be rebuilt."));
@@ -1060,48 +1038,41 @@ SparseMatrix<number>::Tmmult (SparseMatrix<numberC>       &C,
               ExcMessage ("Can't use the same sparsity pattern for "
                           "different matrices if it is to be rebuilt."));
 
-      // need to change the sparsity pattern of
-      // C, so cast away const-ness.
+      // need to change the sparsity pattern of C, so cast away const-ness.
       SparsityPattern &sp_C =
         *(const_cast<SparsityPattern *>(&C.get_sparsity_pattern()));
       C.clear();
       sp_C.reinit (0,0,0);
 
-      // create a sparsity pattern for the
-      // matrix. we will go through all the
-      // rows in the matrix A, and for each
-      // column in a row we add the whole row
-      // of matrix B with that row number. This
-      // means that we will insert a lot of
-      // entries to each row, which is best
-      // handled by the
+      // create a sparsity pattern for the matrix. we will go through all the
+      // rows in the matrix A, and for each column in a row we add the whole
+      // row of matrix B with that row number. This means that we will insert
+      // a lot of entries to each row, which is best handled by the
       // CompressedSimpleSparsityPattern class.
       {
         CompressedSimpleSparsityPattern csp (n(), B.n());
         for (unsigned int i = 0; i < sp_A.n_rows(); ++i)
           {
             const unsigned int *rows =
-              &sp_A.get_column_numbers()[sp_A.get_rowstart_indices()[i]];
+              &sp_A.colnums[sp_A.rowstart[i]];
             const unsigned int *const end_rows =
-              &sp_A.get_column_numbers()[sp_A.get_rowstart_indices()[i+1]];
+              &sp_A.colnums[sp_A.rowstart[i+1]];
+            // cast away constness to conform with csp.add_entries interface
             unsigned int *new_cols = const_cast<unsigned int *>
-                                     (&sp_B.get_column_numbers()
-                                      [sp_B.get_rowstart_indices()[i]]);
+                                     (&sp_B.colnums[sp_B.rowstart[i]]);
             unsigned int *end_new_cols = const_cast<unsigned int *>
-                                         (&sp_B.get_column_numbers()
-                                          [sp_B.get_rowstart_indices()[i+1]]);
+                                         (&sp_B.colnums[sp_B.rowstart[i+1]]);
 
-            if (sp_B.optimize_diagonal() == true)
+            if (sp_B.n_rows() == sp_B.n_cols())
               ++new_cols;
 
             for (; rows != end_rows; ++rows)
               {
                 const unsigned int row = *rows;
 
-                // if B has a diagonal, need to add that
-                // manually. this way, we maintain
-                // sortedness.
-                if (sp_B.optimize_diagonal() == true)
+                // if B has a diagonal, need to add that manually. this way,
+                // we maintain sortedness.
+                if (sp_B.n_rows() == sp_B.n_cols())
                   csp.add(row, i);
 
                 csp.add_entries (row, new_cols, end_new_cols, true);
@@ -1125,47 +1096,37 @@ SparseMatrix<number>::Tmmult (SparseMatrix<numberC>       &C,
     max_n_cols_B = std::max (max_n_cols_B, sp_B.row_length(i));
   std::vector<numberC> new_entries(max_n_cols_B);
 
-  // now compute the actual entries: a
-  // matrix-matrix product involves three
-  // nested loops. One over the rows of A,
-  // for each row we then loop over all the
-  // columns, and then we need to multiply
-  // each element with all the elements in
-  // that row in B.
+  // now compute the actual entries: a matrix-matrix product involves three
+  // nested loops. One over the rows of A, for each row we then loop over all
+  // the columns, and then we need to multiply each element with all the
+  // elements in that row in B.
   for (unsigned int i=0; i<m(); ++i)
     {
-      const unsigned int *rows =
-        &sp_A.get_column_numbers()[sp_A.get_rowstart_indices()[i]];
-      const unsigned int *const end_rows =
-        &sp_A.get_column_numbers()[sp_A.get_rowstart_indices()[i+1]];
-      const unsigned int *new_cols =
-        (&sp_B.get_column_numbers()[sp_B.get_rowstart_indices()[i]]);
-      if (sp_B.optimize_diagonal())
+      const unsigned int *rows = &sp_A.colnums[sp_A.rowstart[i]];
+      const unsigned int *const end_rows = &sp_A.colnums[sp_A.rowstart[i+1]];
+      const unsigned int *new_cols = &sp_B.colnums[sp_B.rowstart[i]];
+      if (sp_B.n_rows() == sp_B.n_cols())
         ++new_cols;
 
-      const numberB *const end_cols =
-        &B.val[&sp_B.get_column_numbers()[sp_B.get_rowstart_indices()[i+1]]-
-               &sp_B.get_column_numbers()[sp_B.get_rowstart_indices()[0]]];
+      const numberB *const end_cols = &B.val[sp_B.rowstart[i+1]];
 
       for (; rows != end_rows; ++rows)
         {
           const unsigned int row = *rows;
-          const double A_val = val[rows-&sp_A.get_column_numbers()[sp_A.get_rowstart_indices()[0]]];
+          const double A_val = val[rows-&sp_A.colnums[sp_A.rowstart[0]]];
 
           // special treatment for diagonal
-          if (sp_B.optimize_diagonal())
+          if (sp_B.n_rows () == sp_B.n_cols())
             C.add (row, i, A_val *
-                   B.val[new_cols-1-&sp_B.get_column_numbers()
-                         [sp_B.get_rowstart_indices()[0]]] *
+                   B.val[new_cols-1-&sp_B.colnums[sp_B.rowstart[0]]] *
                    (use_vector ? V(i) : 1));
 
-          // now the innermost loop that goes over
-          // all the elements in row 'col' of
-          // matrix B. Cache the elements, and then
-          // write them into C at once
+          // now the innermost loop that goes over all the elements in row
+          // 'col' of matrix B. Cache the elements, and then write them into C
+          // at once
           numberC *new_ptr = &new_entries[0];
           const numberB *B_val_ptr =
-            &B.val[new_cols-&sp_B.get_column_numbers()[sp_B.get_rowstart_indices()[0]]];
+            &B.val[new_cols-&sp_B.colnums[sp_B.rowstart[0]]];
           for (; B_val_ptr != end_cols; ++B_val_ptr)
             *new_ptr++ = A_val * *B_val_ptr * (use_vector ? V(i) : 1);
 
@@ -1320,11 +1281,9 @@ SparseMatrix<number>::precondition_Jacobi (Vector<somenumber>       &dst,
 {
   Assert (cols != 0, ExcNotInitialized());
   Assert (val != 0, ExcNotInitialized());
-  Assert (cols->optimize_diagonal(),
-          typename SparsityPattern::ExcDiagonalNotOptimized());
-
-  Assert (dst.size() == n(), ExcDimensionMismatch (dst.size(), n()));
-  Assert (src.size() == n(), ExcDimensionMismatch (src.size(), n()));
+  AssertDimension (m(), n());
+  AssertDimension (dst.size(), n());
+  AssertDimension (src.size(), n());
 
   const unsigned int n = src.size();
   somenumber              *dst_ptr = dst.begin();
@@ -1358,7 +1317,7 @@ void
 SparseMatrix<number>::precondition_SSOR (Vector<somenumber>              &dst,
                                          const Vector<somenumber>        &src,
                                          const number                     om,
-                                         const std::vector<unsigned int> &pos_right_of_diagonal) const
+                                         const std::vector<std::size_t>  &pos_right_of_diagonal) const
 {
   // to understand how this function works
   // you may want to take a look at the CVS
@@ -1366,14 +1325,12 @@ SparseMatrix<number>::precondition_SSOR (Vector<somenumber>              &dst,
   // which is much clearer...
   Assert (cols != 0, ExcNotInitialized());
   Assert (val != 0, ExcNotInitialized());
-  Assert (cols->optimize_diagonal(),
-          typename SparsityPattern::ExcDiagonalNotOptimized());
-
-  Assert (dst.size() == n(), ExcDimensionMismatch (dst.size(), n()));
-  Assert (src.size() == n(), ExcDimensionMismatch (src.size(), n()));
+  AssertDimension (m(), n());
+  AssertDimension (dst.size(), n());
+  AssertDimension (src.size(), n());
 
   const unsigned int  n            = src.size();
-  const std::size_t *rowstart_ptr = &cols->rowstart[0];
+  const std::size_t  *rowstart_ptr = &cols->rowstart[0];
   somenumber         *dst_ptr      = &dst(0);
 
   // case when we have stored the position
@@ -1388,7 +1345,7 @@ SparseMatrix<number>::precondition_SSOR (Vector<somenumber>              &dst,
       for (unsigned int row=0; row<n; ++row, ++dst_ptr, ++rowstart_ptr)
         {
           *dst_ptr = src(row);
-          const unsigned int first_right_of_diagonal_index =
+          const std::size_t first_right_of_diagonal_index =
             pos_right_of_diagonal[row];
           Assert (first_right_of_diagonal_index <= *(rowstart_ptr+1),
                   ExcInternalError());
@@ -1492,9 +1449,6 @@ SparseMatrix<number>::precondition_SOR (Vector<somenumber> &dst,
 {
   Assert (cols != 0, ExcNotInitialized());
   Assert (val != 0, ExcNotInitialized());
-  Assert (cols->optimize_diagonal(),
-          typename SparsityPattern::ExcDiagonalNotOptimized());
-
 
   dst = src;
   SOR(dst,om);
@@ -1510,9 +1464,6 @@ SparseMatrix<number>::precondition_TSOR (Vector<somenumber> &dst,
 {
   Assert (cols != 0, ExcNotInitialized());
   Assert (val != 0, ExcNotInitialized());
-  Assert (cols->optimize_diagonal(),
-          typename SparsityPattern::ExcDiagonalNotOptimized());
-
 
   dst = src;
   TSOR(dst,om);
@@ -1527,10 +1478,8 @@ SparseMatrix<number>::SOR (Vector<somenumber> &dst,
 {
   Assert (cols != 0, ExcNotInitialized());
   Assert (val != 0, ExcNotInitialized());
-  Assert (cols->optimize_diagonal(),
-          typename SparsityPattern::ExcDiagonalNotOptimized());
-
-  Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
+  AssertDimension (m(), n());
+  AssertDimension (dst.size(), n());
 
   for (unsigned int row=0; row<m(); ++row)
     {
@@ -1556,10 +1505,8 @@ SparseMatrix<number>::TSOR (Vector<somenumber> &dst,
 {
   Assert (cols != 0, ExcNotInitialized());
   Assert (val != 0, ExcNotInitialized());
-  Assert (cols->optimize_diagonal(),
-          typename SparsityPattern::ExcDiagonalNotOptimized());
-
-  Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
+  AssertDimension (m(), n());
+  AssertDimension (dst.size(), n());
 
   unsigned int row=m()-1;
   while (true)
@@ -1590,8 +1537,7 @@ SparseMatrix<number>::PSOR (Vector<somenumber> &dst,
 {
   Assert (cols != 0, ExcNotInitialized());
   Assert (val != 0, ExcNotInitialized());
-  Assert (cols->optimize_diagonal(),
-          typename SparsityPattern::ExcDiagonalNotOptimized());
+  AssertDimension (m(), n());
 
   Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
   Assert (m() == permutation.size(),
@@ -1629,8 +1575,7 @@ SparseMatrix<number>::TPSOR (Vector<somenumber> &dst,
 {
   Assert (cols != 0, ExcNotInitialized());
   Assert (val != 0, ExcNotInitialized());
-  Assert (cols->optimize_diagonal(),
-          typename SparsityPattern::ExcDiagonalNotOptimized());
+  AssertDimension (m(), n());
 
   Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
   Assert (m() == permutation.size(),
@@ -1666,8 +1611,7 @@ SparseMatrix<number>::Jacobi_step (Vector<somenumber> &v,
 {
   Assert (cols != 0, ExcNotInitialized());
   Assert (val != 0, ExcNotInitialized());
-  Assert (cols->optimize_diagonal(),
-          typename SparsityPattern::ExcDiagonalNotOptimized());
+  AssertDimension (m(), n());
 
   Assert (m() == v.size(), ExcDimensionMismatch(m(),v.size()));
   Assert (m() == b.size(), ExcDimensionMismatch(m(),b.size()));
@@ -1698,9 +1642,7 @@ SparseMatrix<number>::SOR_step (Vector<somenumber> &v,
 {
   Assert (cols != 0, ExcNotInitialized());
   Assert (val != 0, ExcNotInitialized());
-  Assert (cols->optimize_diagonal(),
-          typename SparsityPattern::ExcDiagonalNotOptimized());
-
+  AssertDimension (m(), n());
   Assert (m() == v.size(), ExcDimensionMismatch(m(),v.size()));
   Assert (m() == b.size(), ExcDimensionMismatch(m(),b.size()));
 
@@ -1727,9 +1669,7 @@ SparseMatrix<number>::TSOR_step (Vector<somenumber> &v,
 {
   Assert (cols != 0, ExcNotInitialized());
   Assert (val != 0, ExcNotInitialized());
-  Assert (cols->optimize_diagonal(),
-          typename SparsityPattern::ExcDiagonalNotOptimized());
-
+  AssertDimension (m(), n());
   Assert (m() == v.size(), ExcDimensionMismatch(m(),v.size()));
   Assert (m() == b.size(), ExcDimensionMismatch(m(),b.size()));
 
@@ -1771,9 +1711,7 @@ SparseMatrix<number>::SSOR (Vector<somenumber> &dst,
 
   Assert (cols != 0, ExcNotInitialized());
   Assert (val != 0, ExcNotInitialized());
-  Assert (cols->optimize_diagonal(),
-          typename SparsityPattern::ExcDiagonalNotOptimized());
-
+  AssertDimension (m(), n());
   Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
 
   const unsigned int  n = dst.size();
index 08718a1c62bf52226b3c798ffc16e528426d41d5..0e0b63acb9d0707f6c96469e31ac5eec9bbc1292 100644 (file)
@@ -790,7 +790,7 @@ public:
   void precondition_SSOR (Vector<somenumber>       &dst,
                           const Vector<somenumber> &src,
                           const number              om = 1.,
-                          const std::vector<unsigned int> &pos_right_of_diagonal = std::vector<unsigned int>()) const;
+                          const std::vector<std::size_t> &pos_right_of_diagonal = std::vector<std::size_t>()) const;
 
   /**
    * Apply SOR preconditioning matrix to @p src.
index c3cb1ec036d57f9b2a2f13a1bc9e6b74b76a89dc..1c4ee5b8f644e3596e798d04ea9535aa1fc6c825 100644 (file)
@@ -328,7 +328,7 @@ void
 SparseMatrixEZ<number>::precondition_SSOR (Vector<somenumber>       &dst,
                                            const Vector<somenumber> &src,
                                            const number              om,
-                                           const std::vector<unsigned int> &) const
+                                           const std::vector<std::size_t> &) const
 {
   Assert (m() == n(), ExcNotQuadratic());
   Assert (dst.size() == n(), ExcDimensionMismatch (dst.size(), n()));
index 2ce3a691ab0f2439620c1745b394bdd44439960b..a958f2342405cedde3b40b0173bd6c55daf24de7 100644 (file)
@@ -67,7 +67,7 @@ for (S1, S2 : REAL_SCALARS)
       precondition_SSOR<S2> (Vector<S2> &,
                             const Vector<S2> &,
                             const S1,
-                            const std::vector<unsigned int>&) const;
+                            const std::vector<std::size_t>&) const;
 
     template void SparseMatrix<S1>::
       precondition_SOR<S2> (Vector<S2> &,
index f988b590d2c25d214fcad7de24685477ccfc0d3e..1747bc210b7063c8721c5ecd4a6f562fd86f6159 100644 (file)
@@ -37,7 +37,7 @@ for (S1, S2 : REAL_SCALARS)
       void SparseMatrixEZ<S1>::precondition_SSOR<S2> (Vector<S2> &,
                                                      const Vector<S2> &,
                                                      const S1,
-                                                     const std::vector<unsigned int>&) const;
+                                                     const std::vector<std::size_t>&) const;
     template
       void SparseMatrixEZ<S1>::precondition_SOR<S2> (Vector<S2> &,
                                                     const Vector<S2> &,

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.