]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Added a block sparse matrix class based on Trilinos (parallel) matrices.
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 22 Aug 2008 16:17:11 +0000 (16:17 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 22 Aug 2008 16:17:11 +0000 (16:17 +0000)
git-svn-id: https://svn.dealii.org/trunk@16653 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/trilinos_block_sparse_matrix.h [new file with mode: 0644]
deal.II/lac/source/trilinos_block_sparse_matrix.cc [new file with mode: 0644]
deal.II/lac/source/trilinos_sparse_matrix.cc

diff --git a/deal.II/lac/include/lac/trilinos_block_sparse_matrix.h b/deal.II/lac/include/lac/trilinos_block_sparse_matrix.h
new file mode 100644 (file)
index 0000000..c9a5a81
--- /dev/null
@@ -0,0 +1,407 @@
+//---------------------------------------------------------------------------
+//    $Id: trilinos_block_sparse_matrix.h 14783 2007-06-18 14:52:01Z bangerth $
+//    Version: $Name$
+//
+//    Copyright (C) 2004, 2005, 2006, 2007 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+#ifndef __deal2__trilinos_block_sparse_matrix_h
+#define __deal2__trilinos_block_sparse_matrix_h
+
+
+#include <base/config.h>
+#include <base/table.h>
+#include <lac/block_matrix_base.h>
+#include <lac/trilinos_sparse_matrix.h>
+#include <lac/trilinos_block_vector.h>
+#include <lac/exceptions.h>
+
+#include <cmath>
+
+
+#ifdef DEAL_II_USE_TRILINOS
+
+DEAL_II_NAMESPACE_OPEN
+
+
+
+namespace TrilinosWrappers
+{
+  
+/*! @addtogroup TrilinosWrappers
+ *@{
+ */
+
+/**
+ * Blocked sparse matrix based on the TrilinosWrappers::SparseMatrix class. This
+ * class implements the functions that are specific to the Trilinos SparseMatrix
+ * base objects for a blocked sparse matrix, and leaves the actual work
+ * relaying most of the calls to the individual blocks to the functions
+ * implemented in the base class. See there also for a description of when
+ * this class is useful.
+ *
+ * In contrast to the deal.II-type SparseMatrix class, the Trilinos matrices do
+ * not have external objects for the sparsity patterns. Thus, one does not
+ * determine the size of the individual blocks of a block matrix of this type
+ * by attaching a block sparsity pattern, but by calling reinit() to set the
+ * number of blocks and then by setting the size of each block separately. In
+ * order to fix the data structures of the block matrix, it is then necessary
+ * to let it know that we have changed the sizes of the underlying
+ * matrices. For this, one has to call the collect_sizes() function, for much
+ * the same reason as is documented with the BlockSparsityPattern class.
+ *
+ * @ingroup Matrix1
+ * @author Wolfgang Bangerth, 2004
+ */
+  class BlockSparseMatrix : public BlockMatrixBase<TrilinosWrappers::SparseMatrix>
+  {
+    public:
+                                       /**
+                                        * Typedef the base class for simpler
+                                        * access to its own typedefs.
+                                        */
+      typedef BlockMatrixBase<SparseMatrix> BaseClass;
+    
+                                       /**
+                                        * Typedef the type of the underlying
+                                        * matrix.
+                                        */
+      typedef BaseClass::BlockType  BlockType;
+
+                                       /**
+                                        * Import the typedefs from the base
+                                        * class.
+                                        */
+      typedef BaseClass::value_type      value_type;
+      typedef BaseClass::pointer         pointer;
+      typedef BaseClass::const_pointer   const_pointer;
+      typedef BaseClass::reference       reference;
+      typedef BaseClass::const_reference const_reference;
+      typedef BaseClass::size_type       size_type;
+      typedef BaseClass::iterator        iterator;
+      typedef BaseClass::const_iterator  const_iterator;
+
+                                       /**
+                                        * Constructor; initializes the
+                                        * matrix to be empty, without
+                                        * any structure, i.e.  the
+                                        * matrix is not usable at
+                                        * all. This constructor is
+                                        * therefore only useful for
+                                        * matrices which are members of
+                                        * a class. All other matrices
+                                        * should be created at a point
+                                        * in the data flow where all
+                                        * necessary information is
+                                        * available.
+                                        *
+                                        * You have to initialize the
+                                        * matrix before usage with
+                                        * reinit(BlockSparsityPattern). The
+                                        * number of blocks per row and
+                                        * column are then determined by
+                                        * that function.
+                                        */
+      BlockSparseMatrix ();
+
+                                       /**
+                                        * Destructor.
+                                        */
+      ~BlockSparseMatrix ();
+
+                                       /** 
+                                        * Pseudo copy operator only copying
+                                        * empty objects. The sizes of the block
+                                        * matrices need to be the same.
+                                        */
+      BlockSparseMatrix &
+      operator = (const BlockSparseMatrix &);
+
+                                       /**
+                                        * This operator assigns a scalar to a
+                                        * matrix. Since this does usually not
+                                        * make much sense (should we set all
+                                        * matrix entries to this value? Only
+                                        * the nonzero entries of the sparsity
+                                        * pattern?), this operation is only
+                                        * allowed if the actual value to be
+                                        * assigned is zero. This operator only
+                                        * exists to allow for the obvious
+                                        * notation <tt>matrix=0</tt>, which
+                                        * sets all elements of the matrix to
+                                        * zero, but keep the sparsity pattern
+                                        * previously used.
+                                        */
+      BlockSparseMatrix &
+      operator = (const double d);
+
+                                       /**
+                                        * Resize the matrix, by setting
+                                        * the number of block rows and
+                                        * columns. This deletes all
+                                        * blocks and replaces them by
+                                        * unitialized ones, i.e. ones
+                                        * for which also the sizes are
+                                        * not yet set. You have to do
+                                        * that by calling the @p reinit
+                                        * functions of the blocks
+                                        * themselves. Do not forget to
+                                        * call collect_sizes() after
+                                        * that on this object.
+                                        *
+                                        * The reason that you have to
+                                        * set sizes of the blocks
+                                        * yourself is that the sizes may
+                                        * be varying, the maximum number
+                                        * of elements per row may be
+                                        * varying, etc. It is simpler
+                                        * not to reproduce the interface
+                                        * of the @p SparsityPattern
+                                        * class here but rather let the
+                                        * user call whatever function
+                                        * she desires.
+                                        */
+      void reinit (const unsigned int n_block_rows,
+                   const unsigned int n_block_columns);
+
+                                       /**
+                                        * Resize the matrix, by using an
+                                       * array of Epetra maps to determine
+                                       * the distribution of the 
+                                       * individual matrices. This function
+                                       * assumes that a quadratic block
+                                       * matrix is generated.
+                                        */
+      //void reinit (const std::vector<Epetra_Map> &input_maps);
+
+                                      /** 
+                                       * This function calls the compress()
+                                       * command of all matrices after 
+                                       * the sparsity pattern has been 
+                                       * generated or the assembly is 
+                                       * completed. Note that all MPI
+                                       * processes need to call this 
+                                       * command (in contrast to sparsity
+                                       * pattern generation, which is 
+                                       * probably only called on one 
+                                       * single process) before any 
+                                       * can complete it.
+                                       */
+      void compress ();
+
+                                       /**
+                                        * This function collects the
+                                        * sizes of the sub-objects and
+                                        * stores them in internal
+                                        * arrays, in order to be able to
+                                        * relay global indices into the
+                                        * matrix to indices into the
+                                        * subobjects. You *must* call
+                                        * this function each time after
+                                        * you have changed the size of
+                                        * the sub-objects.
+                                        */
+      void collect_sizes ();
+
+                                       /**
+                                        * Matrix-vector multiplication:
+                                        * let $dst = M*src$ with $M$
+                                        * being this matrix.
+                                        */
+      void vmult (BlockVector       &dst,
+                  const BlockVector &src) const;
+
+                                       /**
+                                        * Matrix-vector
+                                        * multiplication. Just like the
+                                        * previous function, but only
+                                        * applicable if the matrix has
+                                        * only one block column.
+                                        */
+      void vmult (BlockVector          &dst,
+                  const Vector &src) const;
+
+                                       /**
+                                        * Matrix-vector
+                                        * multiplication. Just like the
+                                        * previous function, but only
+                                        * applicable if the matrix has
+                                        * only one block row.
+                                        */
+      void vmult (Vector    &dst,
+                  const BlockVector &src) const;
+
+                                       /**
+                                        * Matrix-vector
+                                        * multiplication. Just like the
+                                        * previous function, but only
+                                        * applicable if the matrix has
+                                        * only one block.
+                                        */
+      void vmult (Vector       &dst,
+                  const Vector &src) const;
+    
+                                       /**
+                                        * Matrix-vector multiplication:
+                                        * let $dst = M^T*src$ with $M$
+                                        * being this matrix. This
+                                        * function does the same as
+                                        * vmult() but takes the
+                                        * transposed matrix.
+                                        */
+      void Tvmult (BlockVector       &dst,
+                   const BlockVector &src) const;
+  
+                                       /**
+                                        * Matrix-vector
+                                        * multiplication. Just like the
+                                        * previous function, but only
+                                        * applicable if the matrix has
+                                        * only one block row.
+                                        */
+      void Tvmult (BlockVector  &dst,
+                   const Vector &src) const;
+
+                                       /**
+                                        * Matrix-vector
+                                        * multiplication. Just like the
+                                        * previous function, but only
+                                        * applicable if the matrix has
+                                        * only one block column.
+                                        */
+      void Tvmult (Vector    &dst,
+                   const BlockVector &src) const;
+
+                                       /**
+                                        * Matrix-vector
+                                        * multiplication. Just like the
+                                        * previous function, but only
+                                        * applicable if the matrix has
+                                        * only one block.
+                                        */
+      void Tvmult (Vector       &dst,
+                   const Vector &src) const;    
+
+                                       /**
+                                        * Make the clear() function in the
+                                        * base class visible, though it is
+                                        * protected.
+                                        */
+      using BlockMatrixBase<SparseMatrix>::clear;
+      
+                                      /** @addtogroup Exceptions
+                                       * @{
+                                       */
+      
+                                       /**
+                                        * Exception
+                                        */
+      DeclException4 (ExcIncompatibleRowNumbers,
+                      int, int, int, int,
+                      << "The blocks [" << arg1 << ',' << arg2 << "] and ["
+                      << arg3 << ',' << arg4 << "] have differing row numbers.");
+                                       /**
+                                        * Exception
+                                        */
+      DeclException4 (ExcIncompatibleColNumbers,
+                      int, int, int, int,
+                      << "The blocks [" << arg1 << ',' << arg2 << "] and ["
+                      << arg3 << ',' << arg4 << "] have differing column numbers.");
+                                      ///@}
+  };
+
+
+
+/*@}*/
+
+// ------------- inline and template functions -----------------
+
+  inline
+  void
+  BlockSparseMatrix::vmult (BlockVector       &dst,
+                            const BlockVector &src) const
+  {
+    BaseClass::vmult_block_block (dst, src);
+  }
+  
+
+
+  inline
+  void
+  BlockSparseMatrix::vmult (BlockVector  &dst,
+                            const Vector &src) const
+  {
+    BaseClass::vmult_block_nonblock (dst, src);
+  }
+  
+
+
+  inline
+  void
+  BlockSparseMatrix::vmult (Vector            &dst,
+                            const BlockVector &src) const
+  {
+    BaseClass::vmult_nonblock_block (dst, src);
+  }
+  
+
+
+  inline
+  void
+  BlockSparseMatrix::vmult (Vector       &dst,
+                            const Vector &src) const
+  {
+    BaseClass::vmult_nonblock_nonblock (dst, src);
+  }
+
+
+  inline
+  void
+  BlockSparseMatrix::Tvmult (BlockVector       &dst,
+                            const BlockVector &src) const
+  {
+    BaseClass::Tvmult_block_block (dst, src);
+  }
+  
+
+
+  inline
+  void
+  BlockSparseMatrix::Tvmult (BlockVector  &dst,
+                            const Vector &src) const
+  {
+    BaseClass::Tvmult_block_nonblock (dst, src);
+  }
+  
+
+
+  inline
+  void
+  BlockSparseMatrix::Tvmult (Vector            &dst,
+                            const BlockVector &src) const
+  {
+    BaseClass::Tvmult_nonblock_block (dst, src);
+  }
+  
+
+
+  inline
+  void
+  BlockSparseMatrix::Tvmult (Vector       &dst,
+                            const Vector &src) const
+  {
+    BaseClass::Tvmult_nonblock_nonblock (dst, src);
+  }
+  
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif    // DEAL_II_USE_TRILINOS
+
+#endif    // __deal2__trilinos_block_sparse_matrix_h
diff --git a/deal.II/lac/source/trilinos_block_sparse_matrix.cc b/deal.II/lac/source/trilinos_block_sparse_matrix.cc
new file mode 100644 (file)
index 0000000..314c8a1
--- /dev/null
@@ -0,0 +1,108 @@
+//---------------------------------------------------------------------------
+//    $Id: trilinos_block_sparse_matrix.cc 15631 2008-01-17 23:47:31Z bangerth $
+//    Version: $Name$
+//
+//    Copyright (C) 2004, 2005, 2006, 2008 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+#include <lac/trilinos_block_sparse_matrix.h>
+
+
+#ifdef DEAL_II_USE_TRILINOS
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace TrilinosWrappers
+{
+  BlockSparseMatrix::BlockSparseMatrix ()
+  {}
+  
+  
+
+  BlockSparseMatrix::~BlockSparseMatrix ()
+  {}
+
+
+
+  BlockSparseMatrix &
+  BlockSparseMatrix::operator = (const BlockSparseMatrix &m)
+  {
+    BaseClass::operator = (m);
+
+    return *this;
+  }
+
+
+
+  BlockSparseMatrix &
+  BlockSparseMatrix::operator = (const double d)
+  {
+    for (unsigned int r=0; r<this->n_block_rows(); ++r)
+      for (unsigned int c=0; c<this->n_block_cols(); ++c)
+        this->block(r,c) = d;
+
+    return *this;
+  }
+
+
+  void
+  BlockSparseMatrix::
+  reinit (const unsigned int n_block_rows,
+          const unsigned int n_block_columns)
+  {
+                                     // first delete previous content of
+                                     // the subobjects array
+    for (unsigned int r=0; r<this->n_block_rows(); ++r)
+      for (unsigned int c=0; c<this->n_block_cols(); ++c)
+        {
+          BlockType *p = this->sub_objects[r][c];
+          this->sub_objects[r][c] = 0;
+          delete p;
+        }
+  
+    this->sub_objects.reinit (0,0);
+
+                                     // then resize. set sizes of blocks to
+                                     // zero. user will later have to call
+                                     // collect_sizes for this
+    this->sub_objects.reinit (n_block_rows,
+                              n_block_columns);
+    this->row_block_indices.reinit (n_block_rows, 0);
+    this->column_block_indices.reinit (n_block_columns, 0);
+
+                                     // and reinitialize the blocks
+    for (unsigned int r=0; r<this->n_block_rows(); ++r)
+      for (unsigned int c=0; c<this->n_block_cols(); ++c)
+        {
+          BlockType *p = new BlockType();
+          this->sub_objects[r][c] = p;
+        }
+  }
+
+  
+  void
+  BlockSparseMatrix::compress()
+  {
+    for (unsigned int r=0; r<this->n_block_rows(); ++r)
+      for (unsigned int c=0; c<this->n_block_cols(); ++c)
+       this->block(r,c).compress();
+  }
+
+  void
+  BlockSparseMatrix::collect_sizes ()
+  {
+    BaseClass::collect_sizes ();
+  }
+  
+}
+
+  
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index 6e3d0afa0385912a0954c32fe5ee5c803738592e..f8f8872550db8364da9084b316769eee54d7b8d1 100755 (executable)
@@ -194,6 +194,14 @@ namespace TrilinosWrappers
     Assert (matrix->NumGlobalRows() == (int)sparsity_pattern.n_rows(),
            ExcDimensionMismatch (matrix->NumGlobalRows(),
                                  sparsity_pattern.n_rows()));
+    
+                                // Trilinos seems to have a bug for
+                                // rectangular matrices at this point,
+                                // so do not check for consistent 
+                                // column numbers here.
+//    Assert (matrix->NumGlobalCols() == (int)sparsity_pattern.n_cols(),
+//         ExcDimensionMismatch (matrix->NumGlobalCols(),
+//                               sparsity_pattern.n_cols()));
 
     std::vector<int> n_entries_per_row(n_rows);
     
@@ -216,6 +224,7 @@ namespace TrilinosWrappers
     unsigned int n_rows = sparsity_pattern.n_rows();
     matrix.reset();
     row_map = input_map;
+    col_map = row_map;
     
     Assert (input_map.NumGlobalElements() == (int)sparsity_pattern.n_rows(),
            ExcDimensionMismatch (input_map.NumGlobalElements(),
@@ -266,6 +275,15 @@ namespace TrilinosWrappers
     matrix = std::auto_ptr<Epetra_FECrsMatrix>
              (new Epetra_FECrsMatrix(Copy, row_map, col_map,
                                      &n_entries_per_row[0], true));
+      
+//      Assert (matrix->NumGlobalCols() == col_map.NumGlobalElements(),
+//           ExcInternalError());
+//      Assert (matrix->NumGlobalRows() == row_map.NumGlobalElements(),
+//           ExcInternalError());
+//      Assert (matrix->NumGlobalCols() == sparsity_pattern.n_cols(),
+//           ExcInternalError());
+//      Assert (matrix->NumGlobalRows() == sparsity_pattern.n_rows(),
+//           ExcInternalError());
     
     const unsigned int n_max_entries_per_row = *std::max_element (
                    &n_entries_per_row[0], &n_entries_per_row[n_rows-1]);

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.