how to add a parallel MatMatSolve() function?

Barry Smith bsmith at
Fri Feb 29 19:46:19 CST 2008

   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 
  if it crashes or produces incorrect


#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;

   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- 
   ierr = VecCreateMPIWithArray(A- 
   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);

#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)

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

    Most users should usually employ the simplified KSP interface for  
linear solvers
    instead of working directly with matrix algebra routines such as  
    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;

   if (X == B) SETERRQ(PETSC_ERR_ARG_IDN,"X and B must be different  
   if (!A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored  
   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", 
     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);

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

More information about the petsc-users mailing list