]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Write a mmult function also for Trilinos matrices (it works in parallel, too).
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Tue, 5 May 2009 14:41:21 +0000 (14:41 +0000)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Tue, 5 May 2009 14:41:21 +0000 (14:41 +0000)
git-svn-id: https://svn.dealii.org/trunk@18814 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/trilinos_sparse_matrix.h
deal.II/lac/include/lac/trilinos_vector.h
deal.II/lac/source/trilinos_sparse_matrix.cc

index b370ad2d3e6d3d1e22f2b28726bbc71543759156..15e1f611068aeffdc9a75b6e7225da5762b40a6c 100755 (executable)
@@ -1534,6 +1534,34 @@ namespace TrilinosWrappers
                               const VectorBase &x,
                               const VectorBase &b) const;
 
+                                    /**
+                                     * Perform the matrix-matrix
+                                     * multiplication <tt>C = A * B</tt>,
+                                     * or, if an optional vector argument
+                                     * is given, <tt>C = A * diag(V) *
+                                     * B</tt>, where <tt>diag(V)</tt>
+                                     * defines a diagonal matrix with the
+                                     * vector entries.
+                                     *
+                                     * This function assumes that the
+                                     * calling matrix <tt>A</tt> and
+                                     * <tt>B</tt> have compatible
+                                     * sizes. The size of <tt>C</tt> will
+                                     * be set within this function.
+                                     *
+                                     * The content as well as the sparsity
+                                     * pattern of the matrix C will be
+                                     * changed by this function, so make
+                                     * sure that the sparsity pattern is
+                                     * not used somewhere else in your
+                                     * program. This is an expensive
+                                     * operation, so think twice before you
+                                     * use this function.
+                                     */
+    void mmult (SparseMatrix       &C,
+               const SparseMatrix &B,
+               const VectorBase   &V = VectorBase()) const;
+
 //@}
 /**
  * @name Matrix norms
@@ -1545,8 +1573,8 @@ namespace TrilinosWrappers
                                         * <i>l</i><sub>1</sub>-norm of
                                         * the matrix, that is
                                         * $|M|_1=
-                                       * \max_{\mathrm{all columns } j}
-                                       * \sum_{\mathrm{all rows } i} 
+                                       * \max_{\mathrm{all\ columns\ } j}
+                                       * \sum_{\mathrm{all\ rows\ } i} 
                                        * |M_{ij}|$, (max. sum
                                         * of columns).  This is the
                                         * natural matrix norm that is
@@ -1561,9 +1589,9 @@ namespace TrilinosWrappers
                                        /**
                                         * Return the linfty-norm of the
                                         * matrix, that is
-                                        * $|M|_\infty=\max_{\mathrm{all
-                                        * rows} i}\sum_{\mathrm{all
-                                        * columns} j} |M_{ij}|$,
+                                        * $|M|_\infty=\max_{\mathrm{all
+                                        * rows\ } i}\sum_{\mathrm{all\ 
+                                        * columns} j} |M_{ij}|$,
                                         * (max. sum of rows).  This is
                                         * the natural matrix norm that
                                         * is compatible to the
index b87cc53849e55250638ae08a6619d95824f9caaa..55831373409c41a919ea5e114709741fb8ecf69f 100755 (executable)
@@ -305,7 +305,7 @@ namespace TrilinosWrappers
                                        * exception will be thrown.
                                        */
        Vector &
-         operator = (const ::dealii::TrilinosWrappers::Vector &V);
+       operator = (const ::dealii::TrilinosWrappers::Vector &V);
 
                                        /**
                                        * Another copy function. This
index 876ca49fb9017d86fa6e09590e52839fd656de3c..0e942b33d13fb8fd2734046c9dff8bbfd289eebf 100755 (executable)
@@ -22,6 +22,9 @@
 
 #ifdef DEAL_II_USE_TRILINOS
 
+#include <ml_epetra_utils.h>
+#include <ml_struct.h>
+
 DEAL_II_NAMESPACE_OPEN
 
 namespace TrilinosWrappers
@@ -710,16 +713,11 @@ namespace TrilinosWrappers
 
     matrix->FillComplete (col_map, row_map, true);
 
-    int length;
-    int *row_indices;
-    TrilinosScalar *values;
-    for (unsigned int row=0; row<m(); ++row)
-      if (row_map.MyGID(row))
-       {
-         const int local_row = row_map.LID(row);
-         input_matrix.ExtractMyRowView(local_row, length, values, row_indices);
-         matrix->ReplaceMyValues(local_row, length, values, row_indices);
-       }
+    const TrilinosScalar * in_values = input_matrix[0];
+    TrilinosScalar * values = (*matrix)[0];
+    const unsigned int my_nonzeros = input_matrix.NumMyNonzeros();
+    for (unsigned int i=0; i<my_nonzeros; ++i)
+      values[i] = in_values[0];
 
     compress();
   }
@@ -1265,6 +1263,69 @@ namespace TrilinosWrappers
 
 
 
+  void
+  SparseMatrix::mmult (SparseMatrix       &C,
+                      const SparseMatrix &B,
+                      const VectorBase   &V) const
+  {
+    const bool use_vector = V.size() == n() ? true : false;
+    Assert (n() == B.m(), ExcDimensionMismatch(n(), B.m()));
+    Assert (this->matrix->DomainMap().SameAs(B.matrix->RangeMap()),
+           ExcMessage ("Parallel partitioning of A and B does not fit."));
+
+    C.clear();
+                                  // use ML built-in method for performing
+                                  // the matrix-matrix product.
+
+                                  // create ML operators on top of the
+                                  // Epetra matrix.
+    ML_Comm* comm;
+    ML_Comm_Create(&comm);
+    ML_Operator *A_ = ML_Operator_Create(comm);
+    ML_Operator *B_ = ML_Operator_Create(comm);
+    ML_Operator *C_ = ML_Operator_Create(comm);
+
+    ML_Operator_WrapEpetraCrsMatrix(&*matrix,A_,false);
+    ML_Operator_WrapEpetraCrsMatrix(&*B.matrix,B_,false);
+
+                                  // create an Epetra_CrsMatrix with the
+                                  // vector content on the diagonal.
+    Epetra_CrsMatrix Vmat (Copy, this->matrix->DomainMap(), 
+                          B.matrix->RangeMap(), 1, true);
+    const int min_my_gid = this->matrix->DomainMap().MinMyGID();
+    const int max_my_gid = this->matrix->DomainMap().MaxMyGID() + 1;
+    for (int i=min_my_gid; i<max_my_gid; ++i)
+      {
+       double V_val = use_vector ? V(i) : 1;
+       Vmat.InsertGlobalValues (i, 1, &V_val, &i);
+      }
+    Vmat.FillComplete();
+    Vmat.OptimizeStorage();
+
+                                  // wrap even this matrix into ML format
+    ML_Operator *V_ = ML_Operator_Create (comm);
+    ML_Operator_WrapEpetraCrsMatrix (&Vmat, V_, false);
+
+                                  // perform triple matrix-vector product
+    ML_rap (A_, V_, B_, C_, ML_CSR_MATRIX);
+
+    ML_Operator_Destroy (&V_);
+
+                                  // create an Epetra matrix from the ML
+                                  // matrix that we got as a result.
+    Epetra_CrsMatrix * C_mat;
+    ML_Operator2EpetraCrsMatrix(C_, C_mat);
+    C_mat->FillComplete();
+    C_mat->OptimizeStorage();
+    C.reinit (*C_mat);
+
+    ML_Operator_Destroy (&A_);
+    ML_Operator_Destroy (&B_);
+    ML_Operator_Destroy (&C_);
+  }
+
+
+
   void
   SparseMatrix::add (const TrilinosScalar  factor,
                     const SparseMatrix   &rhs)

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.