From 43e1584b93fafe984d6052d2ced8dea4f97e29ca Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Tue, 5 May 2009 14:41:21 +0000 Subject: [PATCH] Write a mmult function also for Trilinos matrices (it works in parallel, too). git-svn-id: https://svn.dealii.org/trunk@18814 0785d39b-7218-0410-832d-ea1e28bc413d --- .../lac/include/lac/trilinos_sparse_matrix.h | 38 +++++++-- deal.II/lac/include/lac/trilinos_vector.h | 2 +- deal.II/lac/source/trilinos_sparse_matrix.cc | 81 ++++++++++++++++--- 3 files changed, 105 insertions(+), 16 deletions(-) diff --git a/deal.II/lac/include/lac/trilinos_sparse_matrix.h b/deal.II/lac/include/lac/trilinos_sparse_matrix.h index b370ad2d3e..15e1f61106 100755 --- a/deal.II/lac/include/lac/trilinos_sparse_matrix.h +++ b/deal.II/lac/include/lac/trilinos_sparse_matrix.h @@ -1534,6 +1534,34 @@ namespace TrilinosWrappers const VectorBase &x, const VectorBase &b) const; + /** + * Perform the matrix-matrix + * multiplication C = A * B, + * or, if an optional vector argument + * is given, C = A * diag(V) * + * B, where diag(V) + * defines a diagonal matrix with the + * vector entries. + * + * This function assumes that the + * calling matrix A and + * B have compatible + * sizes. The size of C 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 * l1-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 diff --git a/deal.II/lac/include/lac/trilinos_vector.h b/deal.II/lac/include/lac/trilinos_vector.h index b87cc53849..5583137340 100755 --- a/deal.II/lac/include/lac/trilinos_vector.h +++ b/deal.II/lac/include/lac/trilinos_vector.h @@ -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 diff --git a/deal.II/lac/source/trilinos_sparse_matrix.cc b/deal.II/lac/source/trilinos_sparse_matrix.cc index 876ca49fb9..0e942b33d1 100755 --- a/deal.II/lac/source/trilinos_sparse_matrix.cc +++ b/deal.II/lac/source/trilinos_sparse_matrix.cc @@ -22,6 +22,9 @@ #ifdef DEAL_II_USE_TRILINOS +#include +#include + 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; rowReplaceMyValues(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; imatrix->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; iFillComplete(); + 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) -- 2.39.5