Tvmult for PreconditionBlockSOR

git-svn-id: https://svn.dealii.org/trunk@3152 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/lac/include/lac/precondition_block.h b/deal.II/lac/include/lac/precondition_block.h
index ea913a9f611d07fe36e73df96ef1c9da31daf027..7bbbcf9e101393e68d01494f7c59c73c61ce5987 100644 (file)
--- a/deal.II/lac/include/lac/precondition_block.h
+++ b/deal.II/lac/include/lac/precondition_block.h
@@ -365,6 +365,12 @@ class PreconditionBlockSOR : public  Subscriptor,
     template <typename number2>
     void vmult (Vector<number2>&, const Vector<number2>&) const;
 
+                                    /**
+                                     * Transpose of @ref{vmult}
+                                     */
+    template <typename number2>
+    void Tvmult (Vector<number2>&, const Vector<number2>&) const;
+
   private:
                                     /**
                                      * Relaxation parameter.

diff --git a/deal.II/lac/include/lac/precondition_block.templates.h b/deal.II/lac/include/lac/precondition_block.templates.h
index 22a34040841151681285531c3c1af3d0a420b557..6d53a57a2a031e483909ec8692154d63a92cf7b7 100644 (file)
--- a/deal.II/lac/include/lac/precondition_block.templates.h
+++ b/deal.II/lac/include/lac/precondition_block.templates.h
@@ -274,7 +274,6 @@ PreconditionBlockSOR<number,inverse_type>::set_omega(number om)
 }
 
 
-//TODO: implement Tvmult
 
 template <typename number, typename inverse_type>
 template <typename number2>
@@ -368,4 +367,105 @@ void PreconditionBlockSOR<number,inverse_type>::vmult (Vector<number2>       &ds
 }
 
 
+template <typename number, typename inverse_type>
+template <typename number2>
+void PreconditionBlockSOR<number,inverse_type>::Tvmult (Vector<number2>       &dst,
+                                                      const Vector<number2> &src) const
+{
+                                  // introduce the following typedef
+                                  // since in the use of exceptions,
+                                  // strict C++ requires us to
+                                  // specify them fully as they are
+                                  // from a template dependent base
+                                  // class. thus, we'd have to write
+                                  // PreconditionBlock<number,inverse_type>::ExcNoMatrixGivenToUse,
+                                  // which is lengthy, but also poses
+                                  // some problems to the
+                                  // preprocessor due to the comma in
+                                  // the template arg list. we could
+                                  // then wrap the whole thing into
+                                  // parentheses, but that creates a
+                                  // parse error for gcc for the
+                                  // exceptions that do not take
+                                  // args...
+  typedef PreconditionBlock<number,inverse_type> BaseClass;
+
+  Assert (A!=0, ExcNotInitialized());
+  
+  const SparseMatrix<number> &M=*A;
+  const unsigned int n_cells=M.m()/blocksize;
+
+  const SparsityPattern &spars    = M.get_sparsity_pattern();
+  const unsigned int    *rowstart = spars.get_rowstart_indices();
+  const unsigned int    *columns  = spars.get_column_numbers();
+
+  Vector<number2> b_cell(blocksize), x_cell(blocksize);
+
+                                      // cell_row, cell_column are the
+                                      // numbering of the blocks (cells).
+                                      // row_cell, column_cell are the local
+                                      // numbering of the unknowns in the
+                                      // blocks.
+                                      // row, column are the global numbering
+                                      // of the unkowns.
+  unsigned int row, column, row_cell, end_diag_block=blocksize *n_cells;
+  number2 b_cell_row;
+
+  if (!inverses_ready())
+    {
+      FullMatrix<number> M_cell(blocksize);
+      for (int icell=n_cells-1; icell>=0 ; --icell)
+       {
+         unsigned int cell = (unsigned int) icell;
+                                          // Collect upper triangle
+         for (row=cell*blocksize, row_cell=0; row_cell<blocksize; ++row_cell, ++row)
+           {
+             b_cell_row=src(row);
+             for (unsigned int j=rowstart[row]; j<rowstart[row+1]; ++j)
+               {
+                 column=columns[j];
+                 if (column >= end_diag_block)
+                   b_cell_row -= M.global_entry(j) * dst(column);
+               }
+             
+             b_cell(row_cell)=b_cell_row;
+                                              // Collect diagonal block
+             for (unsigned int column_cell=0, column=cell*blocksize;
+                  column_cell<blocksize; ++column_cell, ++column)
+                 M_cell(row_cell,column_cell)=M(row,column);
+           }
+         M_cell.householder(b_cell);
+         M_cell.backward(x_cell,b_cell);
+                                          // distribute x_cell to dst
+         for (row=cell*blocksize, row_cell=0; row_cell<blocksize; ++row_cell, ++row)
+           dst(row)=omega*x_cell(row_cell);
+         
+         end_diag_block-=blocksize;
+       }
+    }
+  else
+    for (int icell=n_cells-1; icell >=0 ; --icell)
+      {
+       unsigned int cell = (unsigned int) icell;
+       for (row=cell*blocksize, row_cell=0; row_cell<blocksize; ++row_cell, ++row)
+         {
+           b_cell_row=src(row);
+           for (unsigned int j=rowstart[row]; j<rowstart[row+1]; ++j)
+             {
+               column=columns[j];
+               if (column >= end_diag_block)
+                 b_cell_row -= M.global_entry(j) * dst(column);
+             }
+           b_cell(row_cell)=b_cell_row;
+         }
+       inverse(cell).vmult(x_cell, b_cell);
+                                        // distribute x_cell to dst
+       for (row=cell*blocksize, row_cell=0; row_cell<blocksize; ++row_cell, ++row)
+         dst(row)=omega*x_cell(row_cell);
+       
+       end_diag_block-=blocksize;
+      }
+}
+
+
 #endif

diff --git a/deal.II/lac/source/precondition_block.cc b/deal.II/lac/source/precondition_block.cc
index 426c93ce768235700abc3b836b4b3b7e8625bd65..eed558300dd1707009f1ec94309d187025bad751 100644 (file)
--- a/deal.II/lac/source/precondition_block.cc
+++ b/deal.II/lac/source/precondition_block.cc
@@ -76,6 +76,10 @@ template void PreconditionBlockSOR<float, float>::vmult (
   Vector<float> &, const Vector<float> &) const;
 template void PreconditionBlockSOR<float, float>::vmult (
   Vector<double> &, const Vector<double> &) const;
+template void PreconditionBlockSOR<float, float>::Tvmult (
+  Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockSOR<float, float>::Tvmult (
+  Vector<double> &, const Vector<double> &) const;
 
 
 // the instantiation for class PreconditionBlockSOR<float, double> is skipped
@@ -91,6 +95,10 @@ template void PreconditionBlockSOR<double, float>::vmult (
   Vector<float> &, const Vector<float> &) const;
 template void PreconditionBlockSOR<double, float>::vmult (
   Vector<double> &, const Vector<double> &) const;
+template void PreconditionBlockSOR<double, float>::Tvmult (
+  Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockSOR<double, float>::Tvmult (
+  Vector<double> &, const Vector<double> &) const;
 
 template class PreconditionBlockSOR<double, double>;
 
@@ -98,5 +106,9 @@ template void PreconditionBlockSOR<double, double>::vmult (
   Vector<float> &, const Vector<float> &) const;
 template void PreconditionBlockSOR<double, double>::vmult (
   Vector<double> &, const Vector<double> &) const;
+template void PreconditionBlockSOR<double, double>::Tvmult (
+  Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockSOR<double, double>::Tvmult (
+  Vector<double> &, const Vector<double> &) const;