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

[Pierre Jolivet via petsc-dev]

> On 20 Sep 2019, at 7:36 AM, Jed Brown  > wrote:
> 
> Pierre Jolivet via petsc-dev  > 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, ) 
>> but the result is not correct in at least two cases:
> 
> Do you want MATOP_MAT_MULT and MATOP_TRANSPOSE_MAT_MULT?

Ah, OK, MATMAT => MatMatMat, so you are right.
I’ll make sure MatHasOperation_Transpose and MatHasOperation_Nest return the 
correct value then.

>> 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
>>  
>> 
>>  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,);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.

Perfect, that’s all I need!

Thanks,
Pierre

>> 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.

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

[Jed Brown via petsc-dev]

Pierre Jolivet via petsc-dev  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, ) 
> 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
>  
> 
>  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,);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.