]> https://gitweb.dealii.org/ - dealii.git/commitdiff
new const_iterator and conjugate_add improved
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Sat, 19 Jul 2003 21:03:00 +0000 (21:03 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Sat, 19 Jul 2003 21:03:00 +0000 (21:03 +0000)
git-svn-id: https://svn.dealii.org/trunk@7884 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/precondition_block.h
deal.II/lac/include/lac/precondition_block.templates.h
deal.II/lac/include/lac/sparse_matrix_ez.h

index 86fd63a54027d47265aa41f8f17791b2071fe38f..18eb95603e765247f2f8e563fc3f0654e689c713 100644 (file)
@@ -24,6 +24,9 @@
 template <typename number> class FullMatrix;
 template <typename number> class Vector;
 
+template<class MATRIX, typename inverse_type>
+class PreconditionBlockJacobi;
+
 /*! @addtogroup Preconditioners
  *@{
  */
@@ -197,6 +200,12 @@ class PreconditionBlock : public virtual Subscriptor
                                      * yields a preconditioner.
                                      */
     void set_same_diagonal ();
+
+                                    /**
+                                     * Does the matrix use only one
+                                     * diagonal block?
+                                     */
+    bool same_diagonal () const;
     
                                     /**
                                      * Stores the inverse of the
@@ -232,6 +241,12 @@ class PreconditionBlock : public virtual Subscriptor
                                      */
     unsigned int block_size () const;
 
+                                    /**
+                                     * The number of blocks of the
+                                     * matrix.
+                                     */
+    unsigned int n_blocks() const;
+    
                                     /**
                                      * Determine an estimate for the
                                      * memory consumption (in bytes)
@@ -324,7 +339,12 @@ class PreconditionBlock : public virtual Subscriptor
                                      */
     bool store_diagonals;
 
+                                    /**
+                                     * Number of blocks.
+                                     */
+    unsigned int nblocks;
   private:
+    
                                     /**
                                      * Storage of the inverse
                                      * matrices of the diagonal
@@ -350,7 +370,7 @@ class PreconditionBlock : public virtual Subscriptor
                                      * Flag for diagonal compression.
                                      * @ref set_same_diagonal()
                                      */
-    bool same_diagonal;
+    bool var_same_diagonal;
 };
 
 
@@ -358,7 +378,7 @@ class PreconditionBlock : public virtual Subscriptor
 /**
  * Block Jacobi preconditioning.
  * See @ref{PreconditionBlock} for requirements on the matrix.
- * @author Ralf Hartmann, Guido Kanschat, 1999, 2000
+ * @author Ralf Hartmann, Guido Kanschat, 1999, 2000, 2003
  */
 template<class MATRIX, typename inverse_type = typename MATRIX::value_type>
 class PreconditionBlockJacobi : public virtual Subscriptor,
@@ -371,6 +391,148 @@ class PreconditionBlockJacobi : public virtual Subscriptor,
     typedef typename MATRIX::value_type number;
     
   public:
+                                    /**
+                                     * STL conforming iterator.
+                                     */
+    class const_iterator
+    {
+      private:
+                                         /**
+                                          * Accessor class for iterators
+                                          */
+        class Accessor
+        {
+          public:
+                                             /**
+                                              * Constructor. Since we use
+                                              * accessors only for read
+                                              * access, a const matrix
+                                              * pointer is sufficient.
+                                              */
+            Accessor (const PreconditionBlockJacobi<MATRIX, inverse_type> *matrix,
+                      const unsigned int row);
+
+                                             /**
+                                              * Row number of the element
+                                              * represented by this
+                                              * object.
+                                              */
+            unsigned int row() const;
+
+                                             /**
+                                              * Column number of the
+                                              * element represented by
+                                              * this object.
+                                              */
+            unsigned int column() const;
+
+                                             /**
+                                              * Value of this matrix entry.
+                                              */
+            inverse_type value() const;
+       
+          protected:
+                                             /**
+                                              * The matrix accessed.
+                                              */
+            const PreconditionBlockJacobi<MATRIX, inverse_type>* matrix;
+
+                                            /**
+                                             * Save block size here
+                                             * for further reference.
+                                             */
+           unsigned int bs;
+           
+                                            /**
+                                             * Current block number.
+                                             */
+           unsigned int a_block;
+
+                                            /**
+                                             * Current block for extracting data.
+                                             *
+                                             * This is either
+                                             * @p{a_block} or zero,
+                                             * depending on
+                                             * @p{same_diagonal}.
+                                             */
+           unsigned int u_block;
+
+                                            /**
+                                             * Iterator inside block.
+                                             */
+           typename FullMatrix<inverse_type>::const_iterator b_iterator;
+
+                                            /**
+                                             * End of current block.
+                                             */
+           typename FullMatrix<inverse_type>::const_iterator b_end;
+           
+                                             /**
+                                              * Make enclosing class a
+                                              * friend.
+                                              */
+            friend class const_iterator;
+        };
+        
+      public:
+                                         /**
+                                          * Constructor.
+                                          */ 
+       const_iterator(const PreconditionBlockJacobi<MATRIX, inverse_type>* matrix,
+                      const unsigned int row);
+         
+                                         /**
+                                          * Prefix increment.
+                                          */
+       const_iterator& operator++ ();
+
+                                         /**
+                                          * Postfix increment.
+                                          */
+       const_iterator& operator++ (int);
+
+                                         /**
+                                          * Dereferencing operator.
+                                          */
+       const Accessor& operator* () const;
+
+                                         /**
+                                          * Dereferencing operator.
+                                          */
+       const Accessor* operator-> () const;
+
+                                         /**
+                                          * Comparison. True, if
+                                          * both iterators point to
+                                          * the same matrix
+                                          * position.
+                                          */
+       bool operator == (const const_iterator&) const;
+                                         /**
+                                          * Inverse of @p{==}.
+                                          */
+       bool operator != (const const_iterator&) const;
+
+                                         /**
+                                          * Comparison
+                                          * operator. Result is true
+                                          * if either the first row
+                                          * number is smaller or if
+                                          * the row numbers are
+                                          * equal and the first
+                                          * index is smaller.
+                                          */
+       bool operator < (const const_iterator&) const;
+
+      private:
+                                         /**
+                                          * Store an object of the
+                                          * accessor class.
+                                          */
+        Accessor accessor;
+    };
+
                                     /**
                                      * Make initialization function
                                      * publicly available.
@@ -443,6 +605,29 @@ class PreconditionBlockJacobi : public virtual Subscriptor,
     template <typename number2>
     void Tvmult_add (Vector<number2>&, const Vector<number2>&) const;
 
+                                    /**
+                                     * STL-like iterator with the
+                                     * first entry.
+                                     */
+    const_iterator begin () const;
+
+                                    /**
+                                     * Final iterator.
+                                     */
+    const_iterator end () const;
+    
+                                    /**
+                                     * STL-like iterator with the
+                                     * first entry of row @p{r}.
+                                     */
+    const_iterator begin (const unsigned int r) const;
+
+                                    /**
+                                     * Final iterator of row @p{r}.
+                                     */
+    const_iterator end (const unsigned int r) const;
+    
+
   private:
                                   /**
                                    * Actual implementation of the
@@ -723,6 +908,14 @@ PreconditionBlock<MATRIX, inverse_type>::empty () const
 }
 
 
+template<class MATRIX, typename inverse_type>
+inline unsigned int
+PreconditionBlock<MATRIX, inverse_type>::n_blocks () const
+{
+  return nblocks;
+}
+
+
 template<class MATRIX, typename inverse_type>
 inline inverse_type
 PreconditionBlock<MATRIX, inverse_type>::el (
@@ -745,4 +938,204 @@ PreconditionBlock<MATRIX, inverse_type>::el (
   return B(ib, jb);
 }
 
+//----------------------------------------------------------------------//
+
+template<class MATRIX, typename inverse_type>
+inline
+PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::Accessor::
+Accessor (const PreconditionBlockJacobi<MATRIX, inverse_type>* matrix,
+          const unsigned int row)
+               :
+               matrix(matrix),
+               b_iterator(&matrix->inverse(0), 0, 0),
+               b_end(&matrix->inverse(0), 0, 0)
+{
+  bs = matrix->block_size();
+  a_block = row / bs;
+
+                                  // This is the end accessor, which
+                                  // does not hava a valid block.
+  if (a_block == matrix->n_blocks())
+    return;
+
+  const unsigned int r = row % bs;
+
+  u_block = (matrix->same_diagonal()) ? 0 : a_block;
+
+  b_iterator = matrix->inverse(u_block).begin(r);
+  b_end = matrix->inverse(u_block).end();
+
+  Assert (a_block < matrix->n_blocks(),
+         ExcIndexRange(a_block, 0, matrix->n_blocks()));
+}
+
+
+template<class MATRIX, typename inverse_type>
+inline
+unsigned int
+PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::Accessor::row() const
+{
+  Assert (a_block < matrix->n_blocks(),
+         ExcIteratorPastEnd());
+  
+  return bs * a_block + b_iterator->row();
+}
+
+
+template<class MATRIX, typename inverse_type>
+inline
+unsigned int
+PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::Accessor::column() const
+{
+  Assert (a_block < matrix->n_blocks(),
+         ExcIteratorPastEnd());
+  
+  return bs * a_block + b_iterator->column();
+}
+
+
+template<class MATRIX, typename inverse_type>
+inline
+inverse_type
+PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::Accessor::value() const
+{
+  Assert (a_block < matrix->n_blocks(),
+         ExcIteratorPastEnd());
+  
+  return b_iterator->value();
+}
+
+
+template<class MATRIX, typename inverse_type>
+inline
+PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::
+const_iterator(const PreconditionBlockJacobi<MATRIX, inverse_type> *matrix,
+               const unsigned int row)
+               :
+               accessor(matrix, row)
+{}
+
+
+template<class MATRIX, typename inverse_type>
+inline
+typename PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator &
+PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::operator++ ()
+{
+  Assert (*this != accessor.matrix->end(), ExcIteratorPastEnd());
+  
+  ++accessor.b_iterator;
+  if (accessor.b_iterator == accessor.b_end)
+    {
+      ++accessor.a_block;
+      if (!accessor.matrix->same_diagonal())
+       ++accessor.u_block;
+      
+      if (accessor.a_block < accessor.matrix->n_blocks())
+       {
+         accessor.b_iterator = accessor.matrix->inverse(accessor.a_block).begin();
+         accessor.b_end = accessor.matrix->inverse(accessor.a_block).end();
+       }
+    }
+  return *this;
+}
+
+
+template<class MATRIX, typename inverse_type>
+inline
+const typename PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::Accessor &
+PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::operator* () const
+{
+  return accessor;
+}
+
+
+template<class MATRIX, typename inverse_type>
+inline
+const typename PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::Accessor *
+PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::operator-> () const
+{
+  return &accessor;
+}
+
+
+template<class MATRIX, typename inverse_type>
+inline
+bool
+PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::
+operator == (const const_iterator& other) const
+{
+  if (accessor.a_block == accessor.matrix->n_blocks() &&
+      accessor.a_block == other.accessor.a_block)
+    return true;
+
+  if (accessor.a_block != other.accessor.a_block)
+    return false;
+  
+  return (accessor.row() == other.accessor.row() &&
+          accessor.column() == other.accessor.column());
+}
+
+
+template<class MATRIX, typename inverse_type>
+inline
+bool
+PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::
+operator != (const const_iterator& other) const
+{
+  return ! (*this == other);
+}
+
+
+template<class MATRIX, typename inverse_type>
+inline
+bool
+PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::
+operator < (const const_iterator& other) const
+{
+  return (accessor.row() < other.accessor.row() ||
+         (accessor.row() == other.accessor.row() &&
+           accessor.column() < other.accessor.column()));
+}
+
+
+template<class MATRIX, typename inverse_type>
+inline
+typename PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator
+PreconditionBlockJacobi<MATRIX, inverse_type>::begin () const
+{
+  return const_iterator(this, 0);
+}
+
+
+template<class MATRIX, typename inverse_type>
+inline
+typename PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator
+PreconditionBlockJacobi<MATRIX, inverse_type>::end () const
+{
+  return const_iterator(this, n_blocks() * block_size());
+}
+
+
+template<class MATRIX, typename inverse_type>
+inline
+typename PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator
+PreconditionBlockJacobi<MATRIX, inverse_type>::begin (
+  const unsigned int r) const
+{
+  Assert (r<A->m(), ExcIndexRange(r,0,A->m()));
+  return const_iterator(this, r);
+}
+
+
+
+template<class MATRIX, typename inverse_type>
+inline
+typename PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator
+PreconditionBlockJacobi<MATRIX, inverse_type>::end (
+  const unsigned int r) const
+{
+  Assert (r<A->m(), ExcIndexRange(r,0,A->m()));
+  return const_iterator(this, r+1);
+}
+
 #endif
index f339aca3cbf75698daa8cc89acf6a2c6bce9b5ab..634bfb86ab1f44297dc3095acf80a57fdb9aba72 100644 (file)
@@ -28,7 +28,7 @@ PreconditionBlock<MATRIX,inverse_type>::PreconditionBlock ():
                blocksize(0),
                A(0),
                store_diagonals(false),
-               same_diagonal(false)
+               var_same_diagonal(false)
 {}
 
 
@@ -50,7 +50,7 @@ void PreconditionBlock<MATRIX,inverse_type>::clear ()
   if (var_diagonal.size()!=0)
     var_diagonal.erase(var_diagonal.begin(), var_diagonal.end());
   blocksize     = 0;
-  same_diagonal = false;
+  var_same_diagonal = false;
   A = 0;
 }
 
@@ -69,7 +69,9 @@ void PreconditionBlock<MATRIX,inverse_type>::initialize (
   Assert (A->m()%bsize==0, ExcWrongBlockSize(bsize, A->m()));
   blocksize=bsize;
   relaxation = parameters.relaxation;
-  same_diagonal = parameters.same_diagonal;
+  var_same_diagonal = parameters.same_diagonal;
+  nblocks = A->m()/bsize;
+  
   if (parameters.invert_diagonal)
     invert_diagblocks();
 }
@@ -79,7 +81,7 @@ template <class MATRIX, typename inverse_type>
 const FullMatrix<inverse_type>&
 PreconditionBlock<MATRIX,inverse_type>::inverse(unsigned int i) const
 {
-  if (same_diagonal)
+  if (same_diagonal())
     return var_inverse[0];
   
   Assert (i < var_inverse.size(), ExcIndexRange(i,0,var_inverse.size()));
@@ -93,7 +95,7 @@ PreconditionBlock<MATRIX,inverse_type>::diagonal(unsigned int i) const
 {
   Assert(store_diagonals, ExcDiagonalsNotStored());
   
-  if (same_diagonal)
+  if (same_diagonal())
     return var_diagonal[0];
   
   Assert (i < var_diagonal.size(), ExcIndexRange(i,0,var_diagonal.size()));
@@ -113,7 +115,15 @@ void
 PreconditionBlock<MATRIX,inverse_type>::set_same_diagonal()
 {
   Assert(var_inverse.size()==0, ExcInverseMatricesAlreadyExist());
-  same_diagonal = true;
+  var_same_diagonal = true;
+}
+
+
+template <class MATRIX, typename inverse_type>
+bool
+PreconditionBlock<MATRIX,inverse_type>::same_diagonal() const
+{
+  return var_same_diagonal;
 }
 
 
@@ -134,11 +144,9 @@ void PreconditionBlock<MATRIX,inverse_type>::invert_diagblocks()
   const MATRIX &M=*A;
   Assert (var_inverse.size()==0, ExcInverseMatricesAlreadyExist());
 
-  const unsigned int n_cells = M.m()/blocksize;
-
   FullMatrix<inverse_type> M_cell(blocksize);
 
-  if (same_diagonal)
+  if (same_diagonal())
     {
       deallog << "PreconditionBlock uses only one diagonal block" << std::endl;
                                       // Invert only the first block
@@ -176,7 +184,7 @@ void PreconditionBlock<MATRIX,inverse_type>::invert_diagblocks()
 
                                       // set the @p{var_inverse} array to the right
                                       // size. we could do it like this:
-                                      // var_inverse = vector<>(n_cells,FullMatrix<>())
+                                      // var_inverse = vector<>(nblocks,FullMatrix<>())
                                       // but this would involve copying many
                                       // FullMatrix objects.
                                       //
@@ -185,14 +193,14 @@ void PreconditionBlock<MATRIX,inverse_type>::invert_diagblocks()
       if (store_diagonals)
        {
          std::vector<FullMatrix<inverse_type> >
-           tmp(n_cells, FullMatrix<inverse_type>(blocksize));
+           tmp(nblocks, FullMatrix<inverse_type>(blocksize));
          var_diagonal.swap (tmp);
        }
 
       if (true)
        {
          std::vector<FullMatrix<inverse_type> >
-           tmp(n_cells, FullMatrix<inverse_type>(blocksize));
+           tmp(nblocks, FullMatrix<inverse_type>(blocksize));
          var_inverse.swap (tmp);
                                           // make sure the tmp object
                                           // goes out of scope as
@@ -201,7 +209,7 @@ void PreconditionBlock<MATRIX,inverse_type>::invert_diagblocks()
 
       M_cell.clear ();
       
-      for (unsigned int cell=0; cell<n_cells; ++cell)
+      for (unsigned int cell=0; cell<nblocks; ++cell)
        {
          const unsigned int cell_start = cell*blocksize;
          for (unsigned int row_cell=0; row_cell<blocksize; ++row_cell)
@@ -276,7 +284,6 @@ void PreconditionBlockJacobi<MATRIX,inverse_type>
   Assert(this->A!=0, ExcNotInitialized());
   
   const MATRIX &M=*this->A;
-  const unsigned int n_cells=M.m()/this->blocksize;
 
   Vector<number2> b_cell(this->blocksize), x_cell(this->blocksize);
 
@@ -292,7 +299,7 @@ void PreconditionBlockJacobi<MATRIX,inverse_type>
   if (!this->inverses_ready())
     {
       FullMatrix<number> M_cell(this->blocksize);
-      for (unsigned int cell=0; cell<n_cells; ++cell)
+      for (unsigned int cell=0; cell<nblocks; ++cell)
        {
          for (row=cell*this->blocksize, row_cell=0;
               row_cell<this->blocksize;
@@ -318,7 +325,7 @@ void PreconditionBlockJacobi<MATRIX,inverse_type>
        }
     }
   else
-    for (unsigned int cell=0; cell<n_cells; ++cell)
+    for (unsigned int cell=0; cell<nblocks; ++cell)
       {
        for (row=cell*this->blocksize, row_cell=0;
             row_cell<this->blocksize;
@@ -412,7 +419,6 @@ void PreconditionBlockSOR<MATRIX,inverse_type>::do_vmult (
   Assert (this->A!=0, ExcNotInitialized());
   
   const MATRIX &M=*this->A;
-  const unsigned int n_cells=M.m()/this->blocksize;
 
   Vector<number2> b_cell(this->blocksize), x_cell(this->blocksize);
 
@@ -429,7 +435,7 @@ void PreconditionBlockSOR<MATRIX,inverse_type>::do_vmult (
   if (!this->inverses_ready())
     {
       FullMatrix<number> M_cell(this->blocksize);
-      for (unsigned int cell=0; cell<n_cells; ++cell)
+      for (unsigned int cell=0; cell<nblocks; ++cell)
        {
          for (row = block_start, row_cell = 0;
               row_cell < this->blocksize;
@@ -464,7 +470,7 @@ void PreconditionBlockSOR<MATRIX,inverse_type>::do_vmult (
        }
     }
   else
-    for (unsigned int cell=0; cell<n_cells; ++cell)
+    for (unsigned int cell=0; cell<nblocks; ++cell)
       {
        for (row = block_start, row_cell = 0;
             row_cell<this->blocksize;
@@ -522,7 +528,6 @@ void PreconditionBlockSOR<MATRIX,inverse_type>::do_Tvmult (
   Assert (this->A!=0, ExcNotInitialized());
   
   const MATRIX &M=*this->A;
-  const unsigned int n_cells=M.m()/this->blocksize;
 
   Vector<number2> b_cell(this->blocksize), x_cell(this->blocksize);
 
@@ -534,13 +539,13 @@ void PreconditionBlockSOR<MATRIX,inverse_type>::do_Tvmult (
                                       // row, column are the global numbering
                                       // of the unkowns.
   unsigned int row, row_cell, block_start = 0;
-  unsigned int block_end=this->blocksize *n_cells;
+  unsigned int block_end=this->blocksize *nblocks;
   number2 b_cell_row;
 
   if (!this->inverses_ready())
     {
       FullMatrix<number> M_cell(this->blocksize);
-      for (int icell=n_cells-1; icell>=0 ; --icell)
+      for (int icell=nblocks-1; icell>=0 ; --icell)
        {
          unsigned int cell = (unsigned int) icell;
                                           // Collect upper triangle
@@ -578,7 +583,7 @@ void PreconditionBlockSOR<MATRIX,inverse_type>::do_Tvmult (
        }
     }
   else
-    for (int icell=n_cells-1; icell >=0 ; --icell)
+    for (int icell=nblocks-1; icell >=0 ; --icell)
       {
        unsigned int cell = (unsigned int) icell;
        for (row=cell*this->blocksize, row_cell=0;
@@ -672,10 +677,8 @@ void PreconditionBlockSSOR<MATRIX,inverse_type>::vmult (Vector<number2>       &d
   Vector<inverse_type> cell_src(this->blocksize);
   Vector<inverse_type> cell_dst(this->blocksize);
   
-  const unsigned int n_cells = this->A->m()/this->blocksize;
-
                                   // Multiply with diagonal blocks
-  for (unsigned int cell=0; cell<n_cells; ++cell)
+  for (unsigned int cell=0; cell<nblocks; ++cell)
     {
       unsigned int row = cell*this->blocksize;
       
@@ -704,10 +707,8 @@ void PreconditionBlockSSOR<MATRIX,inverse_type>::Tvmult (Vector<number2>       &
   Vector<inverse_type> cell_src(this->blocksize);
   Vector<inverse_type> cell_dst(this->blocksize);
   
-  const unsigned int n_cells = this->A->m()/this->blocksize;
-
                                   // Multiply with diagonal blocks
-  for (unsigned int cell=0; cell<n_cells; ++cell)
+  for (unsigned int cell=0; cell<nblocks; ++cell)
     {
       unsigned int row = cell*this->blocksize;
       
index a5a456116504e0e98eeeaa2fbdebc872122d2301..ebd4faa72b8554432a45227a2f8386cfa68d0240 100644 (file)
@@ -1635,9 +1635,9 @@ SparseMatrixEZ<number>::conjugate_add (const MATRIXA& A,
                                       // Determine minimal and
                                       // maximal row for a column in
                                       // advance.
-/*
-      std::vector<unsigned int>::minrow(B.n(), B.m());
-      std::vector<unsigned int>::maxrow(B.n(), 0);
+
+      std::vector<unsigned int> minrow(B.n(), B.m());
+      std::vector<unsigned int> maxrow(B.n(), 0);
       while (b1 != b_final)
        {
          const unsigned int r = b1->row();
@@ -1645,29 +1645,55 @@ SparseMatrixEZ<number>::conjugate_add (const MATRIXA& A,
            minrow[b1->column()] = r;
          if (r > maxrow[b1->column()])
            maxrow[b1->column()] = r;
+         ++b1;
        }
-*/
-      
-      b1 = B.begin();
-      while (b1 != b_final)
+
+      typename MATRIXA::const_iterator ai = A.begin();
+      const typename MATRIXA::const_iterator ae = A.end();
+
+      while (ai != ae)
        {
-         const unsigned int i = b1->row();
-         const unsigned int k = b1->column();
-         typename MATRIXB::const_iterator b2 = B.begin();
-         while (b2 != b_final)
+         const typename MATRIXA::value_type a = ai->value();
+                                          // Don't do anything if
+                                          // this entry is zero.
+         if (a == 0.) continue;
+
+                                          // Now, loop over all rows
+                                          // having possibly a
+                                          // nonzero entry in column
+                                          // ai->row()
+         b1 = B.begin(minrow[ai->row()]);
+         const typename MATRIXB::const_iterator
+           be1 = B.end(maxrow[ai->row()]);
+         const typename MATRIXB::const_iterator
+           be2 = B.end(maxrow[ai->column()]);
+
+         while (b1 != be1)
            {
-             const unsigned int j = b2->row();
-             const unsigned int l = b2->column();
-             
-             const typename MATRIXA::value_type a = A.el(k,l);
-             
-             if (a != 0.)
+             const double b1v = b1->value();
+                                              // We need the product
+                                              // of both. If it is
+                                              // zero, we can save
+                                              // the work
+             if (b1->column() == ai->row() && (b1v != 0.))
                {
-                 add (i, j, a * b1->value() * b2->value());
+                 const unsigned int i = b1->row();
+                 
+                 typename MATRIXB::const_iterator
+                   b2 = B.begin(minrow[ai->column()]);
+                 while (b2 != be2)
+                   {
+                     if (b2->column() == ai->column())
+                       {
+                         const unsigned int j = b2->row();
+                         add (i, j, a * b1v * b2->value());
+                       }
+                     ++b2;
+                   }
                }
-             ++b2;
+             ++b1;
            }
-         ++b1;
+         ++ai;
        }
     }
 }

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.