This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

PETScBaseMatrix problem porting code from dolphin 1.0 to 1.4

0 votes

Hi,

I have defined a function to multiply two Petsc matrix AtB. Using dolphin 1.0 I write

/ /C=AtB
    void AtB(boost::shared_ptr<dolfin::PETScMatrix> mA,boost::shared_ptr<dolfin::PETScMatrix> mB,boost::shared_ptr<dolfin::PETScMatrix> mC) {
          boost::shared_ptr<Mat> A = mA->mat();
          boost::shared_ptr<Mat> B = mB->mat();
          boost::shared_ptr<Mat> C = mC->mat();

          MatMatMultTranspose(*A,*B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&(*C));

        }

this code run as expected.

With dolphin 1.4, I changed boost::shared_ptr to std::shared_ptr and MatMatMultTranspose to MatTransposeMatMult

  void AtB(std::shared_ptr<const dolfin::PETScMatrix> mA,std::shared_ptr<const dolfin::PETScMatrix> mB,std::shared_ptr<dolfin::PETScMatrix> mC) {

      Mat* pC=&(mC->mat());  //error: lvalue required as unary ‘&’ operand

      MatTransposeMatMult(mA->mat(),mB->mat(),MAT_INITIAL_MATRIX,PETSC_DEFAULT,pC);

    }

Now I obtain a unary ‘&’ operand error when I want to obtain a Mat* ith the mat() function. Note that MatTransposeMatMult want a Mat* for the C matrix.

In fact the definition of the mat() function have changed in PETScBaseMatrix

dolfin 1.4

class PETScBaseMatrix : .... {
  public:
    ......
    /// Return PETSc Mat pointer
    Mat mat() const
    { return _matA; }
    Mat* pmat() //const
        { return &_matA; }
  protected:

    // PETSc Mat pointer
    Mat _matA;
};

dolfin 1.0

class PETScBaseMatrix : .... {
  public:
    ......

   /// Return PETSc Mat pointer
    boost::shared_ptr<Mat> mat() const
    { return A; }

    /// Return informal string representation (pretty-print)
    virtual std::string str(bool verbose) const = 0;

  protected:

    // PETSc Mat pointer
    boost::shared_ptr<Mat> A;

};

To obtain a Mat* from PETScMatrix I propose to add this member function to PETScBaseMatrix

Mat* pmat() //const
        { return &_matA; }

Maybe this solution is too dangerous or too intrusive, or completely wrong !

I will prefer a solution that don't change the dolfin 1.4 class.

Please give me your expert opinion.

Sincerely yours

Pierre Joyot

asked Jun 5, 2014 by Pierre Joyot FEniCS Novice (170 points)

1 Answer

+1 vote

Here's a matmattransposemult that should work with newest dolfin:

std::shared_ptr<GenericMatrix> MatMatTransposeMult(GenericMatrix& A, GenericMatrix& B)
  {
    const dolfin::PETScMatrix* Ap = &as_type<const dolfin::PETScMatrix>(A);
    const dolfin::PETScMatrix* Bp = &as_type<const dolfin::PETScMatrix>(B);
    Mat CC;
    PetscErrorCode ierr = MatMatTransposeMult(Ap->mat(), Bp->mat(), MAT_INITIAL_MATRIX, PETSC_DEFAULT, &CC);
    dolfin::PETScMatrix CCC = PETScMatrix(CC);
    return CCC.copy();
  }

Or you could do it in Python as explained here

answered Jun 6, 2014 by mikael-mortensen FEniCS Expert (29,340 points)

Thank you,
this solution runs well

...