# how to add a parallel MatMatSolve() function?

Barry Smith bsmith at mcs.anl.gov
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 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

@*/
PetscErrorCode PETSCMAT_DLLEXPORT MatMatSolve(Mat A,Mat B,Mat X)
{
PetscErrorCode ierr;

PetscFunctionBegin;
PetscValidType(A,1);
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