how to add a parallel MatMatSolve() function?

[Barry Smith]

   Please edit the file src/mat/interface/matrix.c and remove the  
function MatMatSolve(). Replace it with the following two functions,
then run "make lib shared" in that directory. Please let us know at petsc-maint at mcs.anl.gov 
  if it crashes or produces incorrect
results.

    Barry



#undef __FUNCT__
#define __FUNCT__ "MatMatSolve_Basic"
PetscErrorCode PETSCMAT_DLLEXPORT MatMatSolve_Basic(Mat A,Mat B,Mat X)
{
   PetscErrorCode ierr;
   Vec            b,x;
   PetscInt       m,N,i;
   PetscScalar    *bb,*xx;

   PetscFunctionBegin;
   ierr = MatGetArray(B,&bb);CHKERRQ(ierr);
   ierr = MatGetArray(X,&xx);CHKERRQ(ierr);
   ierr = MatGetLocalSize(B,&m,PETSC_NULL);CHKERRQ(ierr);  /* number  
local rows */
   ierr = MatGetSize(B,PETSC_NULL,&N);CHKERRQ(ierr);       /* total  
columns in dense matrix */
   ierr = VecCreateMPIWithArray(A- 
 >hdr.comm,m,PETSC_DETERMINE,PETSC_NULL,&b);CHKERRQ(ierr);
   ierr = VecCreateMPIWithArray(A- 
 >hdr.comm,m,PETSC_DETERMINE,PETSC_NULL,&x);CHKERRQ(ierr);
   for (i=0; i<N; i++) {
     ierr = VecPlaceArray(b,bb + i*m);CHKERRQ(ierr);
     ierr = VecPlaceArray(x,xx + i*m);CHKERRQ(ierr);
     ierr = MatSolve(A,b,x);CHKERRQ(ierr);
   }
   ierr = VecDestroy(b);CHKERRQ(ierr);
   ierr = VecDestroy(x);CHKERRQ(ierr);
   PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "MatMatSolve"
/*@
    MatMatSolve - Solves A X = B, given a factored matrix.

    Collective on Mat

    Input Parameters:
+  mat - the factored matrix
-  B - the right-hand-side matrix  (dense matrix)

    Output Parameter:
.  B - the result matrix (dense matrix)

    Notes:
    The matrices b and x cannot be the same.  I.e., one cannot
    call MatMatSolve(A,x,x).

    Notes:
    Most users should usually employ the simplified KSP interface for  
linear solvers
    instead of working directly with matrix algebra routines such as  
this.
    See, e.g., KSPCreate(). However KSP can only solve for one vector  
(column of X)
    at a time.

    Level: developer

    Concepts: matrices^triangular solves

.seealso: MatMatSolveAdd(), MatMatSolveTranspose(),  
MatMatSolveTransposeAdd(), MatLUFactor(), MatCholeskyFactor()
@*/
PetscErrorCode PETSCMAT_DLLEXPORT MatMatSolve(Mat A,Mat B,Mat X)
{
   PetscErrorCode ierr;

   PetscFunctionBegin;
   PetscValidHeaderSpecific(A,MAT_COOKIE,1);
   PetscValidType(A,1);
   PetscValidHeaderSpecific(B,MAT_COOKIE,2);
   PetscValidHeaderSpecific(X,MAT_COOKIE,3);
   PetscCheckSameComm(A,1,B,2);
   PetscCheckSameComm(A,1,X,3);
   if (X == B) SETERRQ(PETSC_ERR_ARG_IDN,"X and B must be different  
matrices");
   if (!A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored  
matrix");
   if (A->cmap.N != X->rmap.N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat A,Mat  
X: global dim %D %D",A->cmap.N,X->rmap.N);
   if (A->rmap.N != B->rmap.N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat A,Mat  
B: global dim %D %D",A->rmap.N,B->rmap.N);
   if (A->rmap.n != B->rmap.n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat A,Mat  
B: local dim %D %D",A->rmap.n,B->rmap.n);
   if (!A->rmap.N && !A->cmap.N) PetscFunctionReturn(0);
   ierr = MatPreallocated(A);CHKERRQ(ierr);

   ierr = PetscLogEventBegin(MAT_MatSolve,A,B,X,0);CHKERRQ(ierr);
   if (!A->ops->matsolve) {
     ierr = PetscInfo1(A,"Mat type %s using basic MatMatSolve", 
((PetscObject)A)->type_name);CHKERRQ(ierr);
     ierr = MatMatSolve_Basic(A,B,X);CHKERRQ(ierr);
   } else {
     ierr = (*A->ops->matsolve)(A,B,X);CHKERRQ(ierr);
   }
   ierr = PetscLogEventEnd(MAT_MatSolve,A,B,X,0);CHKERRQ(ierr);
   ierr = PetscObjectStateIncrease((PetscObject)X);CHKERRQ(ierr);
   PetscFunctionReturn(0);
}

On Feb 29, 2008, at 6:46 PM, Yujie wrote:

> Hi, everyone
>
> I am considering to add a parallel MatMatSolve() into PETSc based on  
> SuperLu_DIST or Spooles. If I want to use it like current sequential  
> MatMatSolve() in the application codes, how to do it? Could you give  
> me some examples about how to add a new function?
> thanks a lot.
>
> Regards,
> Yujie