Pierre Jolivet via petsc-dev <petsc-dev@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.