[Trilinos-Users] Multiplying an Epetra_CrsMatrix by another Epetra_CrsMatrix

[Baker, Christopher G.]

Hello Mehmet,

There is functionality in package EpetraExt for multiplying two Epetra_CrsMatrix objects into a third. However, the effect you seek is best provided by the following methods on the Epetra_CrsMatrix:
Epetra_CrsMatrix::RightScale(Epetra_Vector &x)
Epetra_CrsMatrix::LeftScale(Epetra_Vector &y)
These allow the easy scaling, from the left or right (both, in your case), of a CrsMatrix by a diagonal matrix ( B^{-1/2}, in your case. )

You will need to put the values of B^{-1/2} into a Epetra_Vector. For example, if you have the data stored in two Epetra_CrsMatrix objects, A and B, the following code roughly describes what should be done:
given: RCP<Epetra_CrsMatrix> A, B
assumed: B is diagonal, with strictly positive entries
// declare a Epetra_Vector to hold the diagonal entries of B
RCP<Epetra_Vector> Bdiag = rcp( new Epetra_Vector( B->RowMap() ) );
int epetra_err = B->ExtractDiagonalCopy( *Bdiag ); assert( epetra_err == 0 );
// done with B matrix (it never should have been stored CRS in the first place)
B = Teuchos::null;
// take inverse-root of diagonal elements from B
double * diagPtr = Bdiag->Values();
const int numDiags = Bdiag->MyLength();
for (int i=0; i < numDiags; ++i) {
  diagPtr[i] = 1.0 / sqrt( diagPtr[i] );
}
diagPtr = NULL; // done with this dangerous low-level pointer that I created above
// compute A = B^{-1/2} * A * B^{-1/2}
epetra_err = A->RightScale( *Bdiag );  assert( epetra_err == 0 );
epetra_err = A->LeftScale( *Bdiag );   assert( epetra_err == 0 );

// now, A has been transformed to represent the standard eigenvalue problem
solve the problem...
get eigenvectors into RCP<Epetra_MultiVector> V ...

// perform change of coordinates on V, from standard back to generalized eigenvectors
// we want B^{1/2} * V
// Bdiag contains inverse square root, so we do this via diag(1 / Bdiag) * V
// if you're worried about taking the inverse of these entries twice,
// then you can save a copy after the square root, before the  inversion.
// the call below does V = diag(Bdiag) \ V
// it is documented under Epetra_MultiVector
epetra_err =  V->ReciprocalMultiply( 1.0, *Bdiag, *V, 0.0 ); assert( epetra_err == 0 );

Good luck!

Chris