]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Fix a few deprecation warnings.
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 6 Feb 2013 10:45:52 +0000 (10:45 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 6 Feb 2013 10:45:52 +0000 (10:45 +0000)
git-svn-id: https://svn.dealii.org/trunk@28248 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/include/deal.II/lac/constraint_matrix.templates.h
deal.II/include/deal.II/lac/sparse_decomposition.templates.h
deal.II/include/deal.II/lac/sparse_ilu.templates.h
deal.II/include/deal.II/lac/sparse_mic.templates.h
deal.II/source/lac/constraint_matrix.cc

index c95e72b0a046f0e53312f713f9b42568706d5b9e..56154db94943aabda16a2ea335e1382e502d0539 100644 (file)
@@ -1729,28 +1729,18 @@ namespace internals
   {
     template <typename SparseMatrixIterator>
     static inline
-    void add_value (const double value,
-                    const unsigned int row,
-                    const unsigned int column,
-                    const unsigned int *col_ptr,
-                    const bool   are_on_diagonal,
-                    unsigned int &counter,
-                    SparseMatrixIterator      val_ptr)
+    void add_value (const double          value,
+                    const unsigned int    row,
+                    const unsigned int    column,
+                    SparseMatrixIterator &matrix_values)
     {
       if (value != 0.)
         {
-          Assert (col_ptr != 0,
-                  typename SparseMatrix<typename SparseMatrixIterator::MatrixType::value_type>::ExcInvalidIndex (row, column));
-          if (are_on_diagonal)
-            {
-              val_ptr->value() += value;
-              return;
-            }
-          while (col_ptr[counter] < column)
-            ++counter;
-          Assert (col_ptr[counter] == column,
+          while (matrix_values->column() < column)
+            ++matrix_values;
+          Assert (matrix_values->column() == column,
                   typename SparseMatrix<typename SparseMatrixIterator::MatrixType::value_type>::ExcInvalidIndex(row, column));
-          (val_ptr+counter)->value() += value;
+          matrix_values->value() += value;
         }
     }
   }
@@ -1780,20 +1770,12 @@ namespace internals
     if (sparsity.n_nonzero_elements() == 0)
       return;
 
-    const std::size_t *row_start = sparsity.get_rowstart_indices();
-    const unsigned int *sparsity_struct = sparsity.get_column_numbers();
-
     const unsigned int row = global_rows.global_row(i);
     const unsigned int loc_row = global_rows.local_row(i);
 
-    const unsigned int *col_ptr = sparsity.row_length(row) == 0 ? 0 :
-                                  &sparsity_struct[row_start[row]];
     typename SparseMatrix<number>::iterator
-    val_ptr = (sparsity.row_length(row) == 0 ?
-               sparse_matrix->end() :
-               sparse_matrix->begin(row));
-    const bool optimize_diagonal = sparsity.optimize_diagonal();
-    unsigned int counter = optimize_diagonal;
+      matrix_values = sparse_matrix->begin(row);
+    const bool optimize_diagonal = sparsity.n_rows() == sparsity.n_cols();
 
     // distinguish three cases about what can
     // happen for checking whether the diagonal is
@@ -1814,8 +1796,7 @@ namespace internals
                 const double col_val = matrix_ptr[loc_col];
                 dealiiSparseMatrix::add_value (col_val, row,
                                                global_rows.global_row(j),
-                                               col_ptr, false, counter,
-                                               val_ptr);
+                                               matrix_values);
               }
           }
         else
@@ -1825,55 +1806,57 @@ namespace internals
                 double col_val = resolve_matrix_entry (global_rows, global_rows, i, j,
                                                        loc_row, local_matrix);
                 dealiiSparseMatrix::add_value (col_val, row,
-                                               global_rows.global_row(j), col_ptr,
-                                               false, counter, val_ptr);
+                                               global_rows.global_row(j),
+                                               matrix_values);
               }
           }
       }
     else if (i>=column_start && i<column_end) // case 2: can split loop
       {
+        ++matrix_values; // jump over diagonal element
         if (global_rows.have_indirect_rows() == false)
           {
             AssertIndexRange (loc_row, local_matrix.m());
             const double *matrix_ptr = &local_matrix(loc_row, 0);
 
+            sparse_matrix->begin(row)->value() += matrix_ptr[loc_row];
             for (unsigned int j=column_start; j<i; ++j)
               {
                 const unsigned int loc_col = global_rows.local_row(j);
                 const double col_val = matrix_ptr[loc_col];
                 dealiiSparseMatrix::add_value(col_val, row,
-                                              global_rows.global_row(j), col_ptr,
-                                              false, counter, val_ptr);
+                                              global_rows.global_row(j),
+                                              matrix_values);
               }
-            val_ptr->value() += matrix_ptr[loc_row];
             for (unsigned int j=i+1; j<column_end; ++j)
               {
                 const unsigned int loc_col = global_rows.local_row(j);
                 const double col_val = matrix_ptr[loc_col];
                 dealiiSparseMatrix::add_value(col_val, row,
-                                              global_rows.global_row(j), col_ptr,
-                                              false, counter, val_ptr);
+                                              global_rows.global_row(j),
+                                              matrix_values);
               }
           }
         else
           {
+            sparse_matrix->begin(row)->value() +=
+              resolve_matrix_entry (global_rows, global_rows, i, i,
+                                    loc_row, local_matrix);
             for (unsigned int j=column_start; j<i; ++j)
               {
                 double col_val = resolve_matrix_entry (global_rows, global_rows, i, j,
                                                        loc_row, local_matrix);
                 dealiiSparseMatrix::add_value (col_val, row,
-                                               global_rows.global_row(j), col_ptr,
-                                               false, counter, val_ptr);
+                                               global_rows.global_row(j),
+                                               matrix_values);
               }
-            val_ptr->value() += resolve_matrix_entry (global_rows, global_rows, i, i, loc_row,
-                                                      local_matrix);
             for (unsigned int j=i+1; j<column_end; ++j)
               {
                 double col_val = resolve_matrix_entry (global_rows, global_rows, i, j,
                                                        loc_row, local_matrix);
                 dealiiSparseMatrix::add_value (col_val, row,
-                                               global_rows.global_row(j), col_ptr,
-                                               false, counter, val_ptr);
+                                               global_rows.global_row(j),
+                                               matrix_values);
               }
           }
       }
@@ -1881,6 +1864,7 @@ namespace internals
     // the loop
     else if (global_rows.have_indirect_rows() == false)
       {
+        ++matrix_values; // jump over diagonal element
         AssertIndexRange (loc_row, local_matrix.m());
         const double *matrix_ptr = &local_matrix(loc_row, 0);
 
@@ -1888,22 +1872,27 @@ namespace internals
           {
             const unsigned int loc_col = global_rows.local_row(j);
             const double col_val = matrix_ptr[loc_col];
-            dealiiSparseMatrix::add_value(col_val, row,
-                                          global_rows.global_row(j), col_ptr,
-                                          row==global_rows.global_row(j),
-                                          counter, val_ptr);
+            if (row==global_rows.global_row(j))
+              sparse_matrix->begin(row)->value() += col_val;
+            else
+              dealiiSparseMatrix::add_value(col_val, row,
+                                            global_rows.global_row(j),
+                                            matrix_values);
           }
       }
     else
       {
+        ++matrix_values; // jump over diagonal element
         for (unsigned int j=column_start; j<column_end; ++j)
           {
-            double col_val = resolve_matrix_entry (global_rows, global_rows, i, j,
-                                                   loc_row, local_matrix);
-            dealiiSparseMatrix::add_value (col_val, row,
-                                           global_rows.global_row(j), col_ptr,
-                                           row==global_rows.global_row(j),
-                                           counter, val_ptr);
+            double col_val = resolve_matrix_entry (global_rows, global_rows, i,
+                                                   j, loc_row, local_matrix);
+            if (row==global_rows.global_row(j))
+              sparse_matrix->begin(row)->value() += col_val;
+            else
+              dealiiSparseMatrix::add_value (col_val, row,
+                                             global_rows.global_row(j),
+                                             matrix_values);
           }
       }
   }
index 0f4ca7017fe358809f06af4daa67db60244c1d4b..33a47c6e2226238ad36593b3216617e1c00f6c56 100644 (file)
@@ -125,7 +125,7 @@ void SparseLUDecomposition<number>::initialize (
     }
 
   // now use this sparsity pattern
-  Assert (sparsity_pattern_to_use->optimize_diagonal(),
+  Assert (sparsity_pattern_to_use->n_rows()==sparsity_pattern_to_use->n_cols(),
           typename SparsityPattern::ExcDiagonalNotOptimized());
   decomposed = false;
   {
@@ -156,7 +156,7 @@ decompose (const SparseMatrix<somenumber> &matrix,
 template <typename number>
 void SparseLUDecomposition<number>::reinit (const SparsityPattern &sparsity)
 {
-  Assert (sparsity.optimize_diagonal(),
+  Assert (sparsity.n_rows() == sparsity.n_cols(),
           typename SparsityPattern::ExcDiagonalNotOptimized());
   decomposed = false;
   {
@@ -173,19 +173,19 @@ void
 SparseLUDecomposition<number>::prebuild_lower_bound()
 {
   const unsigned int *const
-  column_numbers = this->get_sparsity_pattern().get_column_numbers();
+    column_numbers = this->get_sparsity_pattern().colnums;
   const std::size_t *const
-  rowstart_indices = this->get_sparsity_pattern().get_rowstart_indices();
+    rowstart_indices = this->get_sparsity_pattern().rowstart;
   const unsigned int N = this->m();
-
-  prebuilt_lower_bound.resize (N);
-
-  for (unsigned int row=0; row<N; row++)
-    {
-      prebuilt_lower_bound[row]
-        = Utilities::lower_bound (&column_numbers[rowstart_indices[row]+1],
-                                  &column_numbers[rowstart_indices[row+1]],
-                                  row);
+   prebuilt_lower_bound.resize (N);
+   for (unsigned int row=0; row<N; row++)
+     {
+       prebuilt_lower_bound[row]
+         = Utilities::lower_bound (&column_numbers[rowstart_indices[row]+1],
+                                   &column_numbers[rowstart_indices[row+1]],
+                                   row);
     }
 }
 
@@ -213,33 +213,24 @@ SparseLUDecomposition<number>::copy_from (const SparseMatrix<somenumber> &matrix
   // awkward way so that we find the corresponding function in the base class.
   SparseMatrix<number>::operator= (number(0));
 
-  // note: pointers to the sparsity
-  // pattern of the old matrix!
-  const std::size_t *const in_rowstart_indices
-    = matrix.get_sparsity_pattern().get_rowstart_indices();
-  const unsigned int *const in_cols
-    = matrix.get_sparsity_pattern().get_column_numbers();
-  const unsigned int *cols = this->get_sparsity_pattern().get_column_numbers();
-  const std::size_t *rowstart_indices =
-    this->get_sparsity_pattern().get_rowstart_indices();
-
-  // both allow more and less entries
-  // in the new matrix
+  // both allow more and less entries in the new matrix
   std::size_t in_index, index;
   for (unsigned int row=0; row<this->m(); ++row)
     {
-      index = rowstart_indices[row];
-      in_index = in_rowstart_indices[row];
-      this->val[index++] = matrix.val[in_index++];
-      while (in_index < in_rowstart_indices[row+1] &&
-             index < rowstart_indices[row+1])
+      typename SparseMatrix<number>::iterator index = this->begin(row);
+      typename SparseMatrix<somenumber>::const_iterator
+        in_index = matrix.begin(row);
+      index->value() = in_index->value();
+      ++index, ++in_index;
+      while (index < this->end(row) && in_index < matrix.end(row))
         {
-          while (cols[index] < in_cols[in_index] && index < rowstart_indices[row+1])
+          while (index->column() < in_index->column() && index < this->end(row))
             ++index;
-          while (in_cols[in_index] < cols[index] && in_index < in_rowstart_indices[row+1])
+          while (in_index->column() < index->column() && in_index < matrix.end(row))
             ++in_index;
 
-          this->val[index++] = matrix.val[in_index++];
+          index->value() = in_index->value();
+          ++index, ++in_index;
         }
     }
 }
@@ -254,7 +245,7 @@ SparseLUDecomposition<number>::strengthen_diagonal_impl ()
     {
       // get the global index of the first
       // non-diagonal element in this row
-      Assert (this->cols->optimize_diagonal(),  ExcNotImplemented());
+      Assert (this->m() == this->n(),  ExcNotImplemented());
       typename SparseMatrix<number>::iterator
       diagonal_element = this->begin(row);
 
index 6bb71fa990bc110c5ea93e9840cc2e55caef269e..f80470bc3c58bc10892cc095947932928456f49c 100644 (file)
@@ -66,15 +66,12 @@ void SparseILU<number>::decompose (const SparseMatrix<somenumber> &matrix,
   if (strengthen_diagonal>0)
     this->strengthen_diagonal_impl();
 
-  // in the following, we implement
-  // algorithm 10.4 in the book by
-  // Saad by translating in essence
-  // the algorithm given at the end
-  // of section 10.3.2, using the
-  // names of variables used there
+  // in the following, we implement algorithm 10.4 in the book by Saad by
+  // translating in essence the algorithm given at the end of section 10.3.2,
+  // using the names of variables used there
   const SparsityPattern     &sparsity = this->get_sparsity_pattern();
-  const std::size_t   *const ia       = sparsity.get_rowstart_indices();
-  const unsigned int *const ja       = sparsity.get_column_numbers();
+  const std::size_t   *const ia      = sparsity.rowstart;
+  const unsigned int *const ja       = sparsity.colnums;
 
   number *luval = this->SparseMatrix<number>::val;
 
@@ -91,22 +88,14 @@ void SparseILU<number>::decompose (const SparseMatrix<somenumber> &matrix,
       for (unsigned int j=j1; j<=j2; ++j)
         iw[ja[j]] = j;
 
-      // the algorithm in the book
-      // works on the elements of row
-      // k left of the
-      // diagonal. however, since we
-      // store the diagonal element
-      // at the first position, start
-      // at the element after the
-      // diagonal and run as long as
-      // we don't walk into the right
-      // half
+      // the algorithm in the book works on the elements of row k left of the
+      // diagonal. however, since we store the diagonal element at the first
+      // position, start at the element after the diagonal and run as long as
+      // we don't walk into the right half
       unsigned int j = j1+1;
 
-      // pathological case: the current row
-      // of the matrix has only the
-      // diagonal entry. then we have
-      // nothing to do.
+      // pathological case: the current row of the matrix has only the
+      // diagonal entry. then we have nothing to do.
       if (j > j2)
         goto label_200;
 
@@ -121,9 +110,7 @@ label_150:
         number t1 = luval[j] * luval[ia[jrow]];
         luval[j] = t1;
 
-        // jj runs from just right of
-        // the diagonal to the end of
-        // the row
+        // jj runs from just right of the diagonal to the end of the row
         unsigned int jj = ia[jrow]+1;
         while (ja[jj] < jrow)
           ++jj;
@@ -141,25 +128,15 @@ label_150:
 
 label_200:
 
-      // in the book there is an
-      // assertion that we have hit
-      // the diagonal element,
-      // i.e. that jrow==k. however,
-      // we store the diagonal
-      // element at the front, so
-      // jrow must actually be larger
-      // than k or j is already in
+      // in the book there is an assertion that we have hit the diagonal
+      // element, i.e. that jrow==k. however, we store the diagonal element at
+      // the front, so jrow must actually be larger than k or j is already in
       // the next row
       Assert ((jrow > k) || (j==ia[k+1]), ExcInternalError());
 
-      // now we have to deal with the
-      // diagonal element. in the
-      // book it is located at
-      // position 'j', but here we
-      // use the convention of
-      // storing the diagonal element
-      // first, so instead of j we
-      // use uptr[k]=ia[k]
+      // now we have to deal with the diagonal element. in the book it is
+      // located at position 'j', but here we use the convention of storing
+      // the diagonal element first, so instead of j we use uptr[k]=ia[k]
       Assert (luval[ia[k]] != 0, ExcInternalError());
 
       luval[ia[k]] = 1./luval[ia[k]];
@@ -181,10 +158,10 @@ void SparseILU<number>::vmult (Vector<somenumber>       &dst,
   Assert (dst.size() == this->m(), ExcDimensionMismatch(dst.size(), this->m()));
 
   const unsigned int N=dst.size();
-  const std::size_t   *const rowstart_indices
-    = this->get_sparsity_pattern().get_rowstart_indices();
+  const std::size_t  *const rowstart_indices
+    = this->get_sparsity_pattern().rowstart;
   const unsigned int *const column_numbers
-    = this->get_sparsity_pattern().get_column_numbers();
+    = this->get_sparsity_pattern().colnums;
 
   // solve LUx=b in two steps:
   // first Ly = b, then
@@ -256,9 +233,9 @@ void SparseILU<number>::Tvmult (Vector<somenumber>       &dst,
 
   const unsigned int N=dst.size();
   const std::size_t   *const rowstart_indices
-    = this->get_sparsity_pattern().get_rowstart_indices();
+    = this->get_sparsity_pattern().rowstart;
   const unsigned int *const column_numbers
-    = this->get_sparsity_pattern().get_column_numbers();
+    = this->get_sparsity_pattern().colnums;
 
   // solve (LU)'x=b in two steps:
   // first U'y = b, then
index fc322ef10271558b1d23a300aae0a1f532bfe6c4..f760b384900a8fff4c1c722d9961f9576c9ea4dc 100644 (file)
@@ -132,12 +132,9 @@ void SparseMIC<number>::decompose (const SparseMatrix<somenumber> &matrix,
   for (unsigned int row=0; row<this->m(); row++)
     inner_sums[row] = get_rowsum(row);
 
-  const unsigned int *const col_nums = this->get_sparsity_pattern().get_column_numbers();
-  const std::size_t *const rowstarts = this->get_sparsity_pattern().get_rowstart_indices();
-
   for (unsigned int row=0; row<this->m(); row++)
     {
-      const number temp = this->diag_element(row);
+      const number temp = this->begin(row)->value();
       number temp1 = 0;
 
       // work on the lower left part of the matrix. we know
@@ -187,15 +184,10 @@ SparseMIC<number>::vmult (Vector<somenumber>       &dst,
   Assert (dst.size() == this->m(), ExcDimensionMismatch(dst.size(), this->m()));
 
   const unsigned int N=dst.size();
-  const std::size_t   *const rowstart_indices = this->get_sparsity_pattern().get_rowstart_indices();
-  const unsigned int *const column_numbers   = this->get_sparsity_pattern().get_column_numbers();
-  // We assume the underlying matrix A is:
-  // A = X - L - U, where -L and -U are
-  // strictly lower- and upper- diagonal
-  // parts of the system.
+  // We assume the underlying matrix A is: A = X - L - U, where -L and -U are
+  // strictly lower- and upper- diagonal parts of the system.
   //
-  // Solve (X-L)X{-1}(X-U) x = b
-  // in 3 steps:
+  // Solve (X-L)X{-1}(X-U) x = b in 3 steps:
   dst = src;
   for (unsigned int row=0; row<N; ++row)
     {
index 550a0a393c23460986d8628aa7b333cbe0c964a1..a03a4791499e09f623bf77fc4039a21798185173 100644 (file)
@@ -768,8 +768,7 @@ void ConstraintMatrix::condense (const SparsityPattern &uncondensed,
 {
   Assert (sorted == true, ExcMatrixNotClosed());
   Assert (uncondensed.is_compressed() == true, ExcMatrixNotClosed());
-  Assert (uncondensed.n_rows() == uncondensed.n_cols(),
-          ExcNotQuadratic());
+  Assert (uncondensed.n_rows() == uncondensed.n_cols(), ExcNotQuadratic());
 
 
   // store for each line of the matrix
@@ -811,28 +810,22 @@ void ConstraintMatrix::condense (const SparsityPattern &uncondensed,
 
 
   next_constraint = lines.begin();
-  // note: in this loop we need not check
-  // whether @p{next_constraint} is a valid
-  // iterator, since @p{next_constraint} is
-  // only evaluated so often as there are
-  // entries in new_line[*] which tells us
-  // which constraints exist
+  // note: in this loop we need not check whether @p{next_constraint} is a
+  // valid iterator, since @p{next_constraint} is only evaluated so often as
+  // there are entries in new_line[*] which tells us which constraints exist
   for (unsigned int row=0; row<uncondensed.n_rows(); ++row)
     if (new_line[row] != -1)
-      // line not constrained
-      // copy entries if column will not
-      // be condensed away, distribute
-      // otherwise
-      for (unsigned int j=uncondensed.get_rowstart_indices()[row];
-           j<uncondensed.get_rowstart_indices()[row+1]; ++j)
-        if (new_line[uncondensed.get_column_numbers()[j]] != -1)
-          condensed.add (new_line[row], new_line[uncondensed.get_column_numbers()[j]]);
+      // line not constrained copy entries if column will not be condensed
+      // away, distribute otherwise
+      for (SparsityPattern::iterator j=uncondensed.begin(row);
+           j<uncondensed.end(row); ++j)
+        if (new_line[j->column()])
+          condensed.add (new_line[row], new_line[j->column()]);
         else
           {
-            // let c point to the constraint
-            // of this column
+            // let c point to the constraint of this column
             std::vector<ConstraintLine>::const_iterator c = lines.begin();
-            while (c->line != uncondensed.get_column_numbers()[j])
+            while (c->line != j->column())
               ++c;
 
             for (unsigned int q=0; q!=c->entries.size(); ++q)
@@ -841,14 +834,14 @@ void ConstraintMatrix::condense (const SparsityPattern &uncondensed,
     else
       // line must be distributed
       {
-        for (unsigned int j=uncondensed.get_rowstart_indices()[row];
-             j<uncondensed.get_rowstart_indices()[row+1]; ++j)
+        for (SparsityPattern::iterator j=uncondensed.begin(row);
+             j<uncondensed.end(row); ++j)
           // for each entry: distribute
-          if (new_line[uncondensed.get_column_numbers()[j]] != -1)
+          if (new_line[j->column()] != -1)
             // column is not constrained
             for (unsigned int q=0; q!=next_constraint->entries.size(); ++q)
               condensed.add (new_line[next_constraint->entries[q].first],
-                             new_line[uncondensed.get_column_numbers()[j]]);
+                             new_line[j->column()]);
 
           else
             // not only this line but
@@ -857,7 +850,7 @@ void ConstraintMatrix::condense (const SparsityPattern &uncondensed,
               // let c point to the constraint
               // of this column
               std::vector<ConstraintLine>::const_iterator c = lines.begin();
-              while (c->line != uncondensed.get_column_numbers()[j]) ++c;
+              while (c->line != j->column()) ++c;
 
               for (unsigned int p=0; p!=c->entries.size(); ++p)
                 for (unsigned int q=0; q!=next_constraint->entries.size(); ++q)
@@ -877,8 +870,7 @@ void ConstraintMatrix::condense (SparsityPattern &sparsity) const
 {
   Assert (sorted == true, ExcMatrixNotClosed());
   Assert (sparsity.is_compressed() == false, ExcMatrixIsClosed());
-  Assert (sparsity.n_rows() == sparsity.n_cols(),
-          ExcNotQuadratic());
+  Assert (sparsity.n_rows() == sparsity.n_cols(), ExcNotQuadratic());
 
   // store for each index whether it must be
   // distributed or not. If entry is
@@ -898,28 +890,17 @@ void ConstraintMatrix::condense (SparsityPattern &sparsity) const
     {
       if (distribute[row] == numbers::invalid_unsigned_int)
         {
-          // regular line. loop over cols all
-          // valid cols. note that this
-          // changes the line we are
-          // presently working on: we add
-          // additional entries. these are
-          // put to the end of the
-          // row. however, as constrained
-          // nodes cannot be constrained to
-          // other constrained nodes, nothing
-          // will happen if we run into these
-          // added nodes, as they can't be
-          // distributed further. we might
-          // store the position of the last
-          // old entry and stop work there,
-          // but since operating on the newly
-          // added ones only takes two
-          // comparisons (column index valid,
-          // distribute[column] necessarily
-          // ==numbers::invalid_unsigned_int),
-          // it is cheaper to not do so and
-          // run right until the end of the
-          // line
+          // regular line. loop over cols all valid cols. note that this
+          // changes the line we are presently working on: we add additional
+          // entries. these are put to the end of the row. however, as
+          // constrained nodes cannot be constrained to other constrained
+          // nodes, nothing will happen if we run into these added nodes, as
+          // they can't be distributed further. we might store the position of
+          // the last old entry and stop work there, but since operating on
+          // the newly added ones only takes two comparisons (column index
+          // valid, distribute[column] necessarily
+          // ==numbers::invalid_unsigned_int), it is cheaper to not do so and
+          // run right until the end of the line
           for (SparsityPattern::iterator entry = sparsity.begin(row);
                ((entry != sparsity.end(row)) &&
                 entry->is_valid_entry());

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.