[petsc-dev] How to check that MatMatMult is available

[Jed Brown]

Pierre Jolivet via petsc-dev <petsc-dev at mcs.anl.gov> writes:

> Hello,
> Given a Mat A, I’d like to know if there is an implementation available for doing C=A*B
> I was previously using MatHasOperation(A, MATOP_MATMAT_MULT, &hasMatMatMult) but the result is not correct in at least two cases:

Do you want MATOP_MAT_MULT and MATOP_TRANSPOSE_MAT_MULT?

> 1) A is a MATTRANSPOSE and the underlying Mat B=A^T has a MatTransposeMatMult implementation (there is currently no MATOP_MATTRANSPOSEMAT_MULT)
> 2) A is a MATNEST. This could be fixed in MatHasOperation_Nest, by checking MATOP_MATMAT_MULT of all matrices in the MATNEST, but this would be incorrect again if there is a single MATTRANPOSE in there
> What is then the proper way to check that I can indeed call MatMatMult(A,…)?
> Do I need to copy/paste all this https://www.mcs.anl.gov/petsc/petsc-current/src/mat/interface/matrix.c.html#line9801 <https://www.mcs.anl.gov/petsc/petsc-current/src/mat/interface/matrix.c.html#line9801> in user code?

Unfortunately, I don't think there is any common interface for the
string handling, though it would make sense to add one because code of
this sort is copied many times:

    /* dispatch based on the type of A and B from their PetscObject's PetscFunctionLists. */
    char multname[256];
    ierr = PetscStrncpy(multname,"MatMatMult_",sizeof(multname));CHKERRQ(ierr);
    ierr = PetscStrlcat(multname,((PetscObject)A)->type_name,sizeof(multname));CHKERRQ(ierr);
    ierr = PetscStrlcat(multname,"_",sizeof(multname));CHKERRQ(ierr);
    ierr = PetscStrlcat(multname,((PetscObject)B)->type_name,sizeof(multname));CHKERRQ(ierr);
    ierr = PetscStrlcat(multname,"_C",sizeof(multname));CHKERRQ(ierr); /* e.g., multname = "MatMatMult_seqdense_seqaij_C" */
    ierr = PetscObjectQueryFunction((PetscObject)B,multname,&mult);CHKERRQ(ierr);
    if (!mult) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"MatMatMult requires A, %s, to be compatible with B, %s",((PetscObject)A)->type_name,((PetscObject)B)->type_name);

> Thanks,
> Pierre
>
> PS: in my case, C and B are always of type MATDENSE. Should we handle
> this in MatMatMult and never error out for such a simple case. 

I would say yes.

> Indeed, one can just loop on the columns of B and C by doing multiple
> MatMult. This is what I’m currently doing in user code when
> hasMatMatMult == PETSC_FALSE.