[Trilinos-Users] EpetraExt::MatrixMatrix::Multiply()

[Roberts, Nathan V.]

Hello all,

I'm working on a geometric multigrid preconditioner.  I have the fine mesh's stiffness matrix in an Epetra_CrsMatrix* A, and I compute a prolongation operator to the coarse mesh, and store this in an Epetra_CrsMatrix* _P.  I would like to compute my coarse grid operator P^T A P and store this in another Epetra_CrsMatrix.

Two calls to EpetraExt::MatrixMatrix::Multiply() would seem to fit the bill:
1. Compute AP := A * P
2. Compute PT_A_P := P^T * AP

My question regards what to do for a row map when constructing the Epetra_CrsMatrix for the second call to Multiply().  The docs indicate that the first matrix's row map must match that of the output argument.  Here, I'm transposing the first matrix (P), so intuitively one might think that I should use P's column map.  This is what my code below attempts, but I'm pretty sure that what's happening is that the column map varies in size from processor to processor, and therefore inappropriate to use as a row map on construction of a Epetra_CrsMatrix.  (When running on multiple processors, Multiply() outputs a message like "MatrixMatrix::Multiply: ERROR, dimensions of result C must match dimensions of op(A) * op(B). C has 4 rows, should have at least 6".)

My question then is twofold:

1. When the first matrix argument to Multiply() is rectangular and transposed, what is the appropriate row map to use to construct the output matrix?
2. How does one construct such a row map?

I'd also be interested in any comments about performance considerations here--though I'm guessing the basic answer is, avoid such explicit constructions.  A previous effort (involving a slightly different coarse operator, which didn't require any matrix multiplications involving the fine operator) completed my tests considerably faster (10-20x, I'd guess) on a single core than the current draft implementation does.

I'd very much appreciate any advice.

Regards,
Nate

  if (_P.get() == NULL) {
    constructProlongationOperator();
  }



  int maxRowSize = _P->MaxNumEntries();



  // compute A * P
  Epetra_CrsMatrix AP(::Copy, _finePartitionMap, maxRowSize);



  int err = EpetraExt::MatrixMatrix::Multiply(*A, false, *_P, false, AP);



  if (err != 0) {
    ...
  }



  Teuchos::RCP<Epetra_CrsMatrix> PT_A_P = Teuchos::rcp( new Epetra_CrsMatrix(::Copy, _P->ColMap(), maxRowSize) );



  // compute P^T * A * P
  err = EpetraExt::MatrixMatrix::Multiply(*_P, true, AP, false, *PT_A_P);
  if (err != 0) {
    ...
  }

----- Nathan V. Roberts, PhD. ---- Argonne Leadership Computing Facility ---- nvroberts at anl.gov<mailto:nvroberts at anl.gov> -----

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://software.sandia.gov/pipermail/trilinos-users/attachments/20141006/6bffbf32/attachment.html>