Based on your comments I have simplified and improved the code slightly it is now
PetscErrorCode MatMatMultNumeric_MPIDense_MPIAIJ(Mat A,Mat B,Mat C)
{
PetscErrorCode ierr;
Mat At,Bt,Ct;
PetscFunctionBegin;
ierr = MatTranspose(A,MAT_INITIAL_MATRIX,&At);CHKERRQ(ierr);
ierr = MatTranspose(B,MAT_INITIAL_MATRIX,&Bt);CHKERRQ(ierr);
ierr = MatMatMult(Bt,At,MAT_INITIAL_MATRIX,1.0,&Ct);CHKERRQ(ierr);
ierr = MatDestroy(At);CHKERRQ(ierr);
ierr = MatDestroy(Bt);CHKERRQ(ierr);
ierr = MatTranspose(Ct,MAT_REUSE_MATRIX,&C);CHKERRQ(ierr);
ierr = MatDestroy(Ct);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "MatMatMultSymbolic_MPIDense_MPIAIJ"
PetscErrorCode MatMatMultSymbolic_MPIDense_MPIAIJ(Mat A,Mat B,PetscReal fill,Mat *C)
{
PetscErrorCode ierr;
PetscInt m=A->rmap.n,n=B->cmap.n;
Mat Cmat;
PetscFunctionBegin;
if (A->cmap.n != B->rmap.n) SETERRQ2(PETSC_ERR_ARG_SIZ,"A->cmap.n %d != B->rmap.n %d\n",A->cmap.n,B->rmap.n);
ierr = MatCreate(A->hdr.comm,&Cmat);CHKERRQ(ierr);
ierr = MatSetSizes(Cmat,m,n,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
ierr = MatSetType(Cmat,MATMPIDENSE);CHKERRQ(ierr);
ierr = MatMPIDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr);
ierr = MatAssemblyBegin(Cmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(Cmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
*C = Cmat;
PetscFunctionReturn(0);
}
As you note there was no need to call the assembly after the transpose. It is, needed, in general where I put it
in the symbolic portion of the code. (You could avoid it by directly setting some flags in the Cmat object
but it is better to do the right thing and call the assembly (which is free in terms of time anyways compare to the
numerical calculations).
Barry
On Mar 3, 2008, at 10:59 AM, Yujie wrote:
Dear Barry:
Why do you call these two functions in your code? Because I have realized the same method with yours to calculate MPIDENSE_MPIAIJ. However, I didn't call MatAssemblyBegin and MatAssemblyEnd after I call MatTranspose for matrix C.
In application, sometimes I don't know what do the functions I call in PETSc do? So, I don't know whether I should call these two functions although you give me these conditions. Regarding this example, which condition should it belong to? thanks a lot.
Regards,
Yujie
On 3/3/08, Barry Smith <bsmith@xxxxxxxxxxx> wrote:
MatAssembly does
1) moves any values set on one process that belongs on another
2) finalizes the sparse matrix data structure (for example takes any
holes out of the CSR format)
doesn't need to do this for dense matrices
3) sets up whatever VecScatters that are need to do matrix vector
products
4) marks the matrix as assembled and ready to be used
Barry
On Mar 2, 2008, at 11:47 PM, Yujie wrote:
> Hi, everyone
>
> I just check the codes added by Barry for MPIDENSE_MPIAIJ. The codes
> is as follows. I am wondering what context one may call
> MatAssemblyBegin and MatAssemblyEnd?
> From the webpage, I can find some information about two functions:
> "Begins assembling the matrix. This routine should be called after
> completing all calls to MatSetValues()."
> I mean that I don't know how to judge whether I should call two
> functions. Could you give me some advice? thanks a lot.
>
> 4723 */
> 4724 PetscErrorCode MatMatMultNumeric_MPIDense_MPIAIJ(Mat A,Mat
> B,Mat C)
> 4725 {
> 4726 PetscErrorCode ierr;
> 4727 Mat_MPIDense *sub_c = (Mat_MPIDense*)C->data;
> 4728 Mat At,Bt,Ct;
> 4729
> 4730 PetscFunctionBegin;
> 4731 ierr = MatTranspose(A,MAT_INITIAL_MATRIX,&At);CHKERRQ(ierr);
> 4732 ierr = MatTranspose(B,MAT_INITIAL_MATRIX,&Bt);CHKERRQ(ierr);
> 4733 ierr = MatMatMult(Bt,At,MAT_INITIAL_MATRIX,
> 1.0,&Ct);CHKERRQ(ierr);
> 4734 ierr = MatDestroy(At);CHKERRQ(ierr);
> 4735 ierr = MatDestroy(Bt);CHKERRQ(ierr);
> 4736 ierr = MatTranspose(Ct,MAT_REUSE_MATRIX,&C);CHKERRQ(ierr);
> 4737 ierr = MatDestroy(Ct);CHKERRQ(ierr);
> 4738
> 4739 C->assembled = PETSC_TRUE;
> 4740 ierr = MatAssemblyBegin(sub_c-
> >A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
> 4741 ierr = MatAssemblyEnd(sub_c->A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
> 4742 if (!C->was_assembled) {
> 4743 ierr = MatSetUpMultiply_MPIDense(C);CHKERRQ(ierr);
> 4744 }
> 4745 PetscFunctionReturn(0);
> 4746 }
>
>
>
> Regards,
> Yujie