Dear Barry:<br><br>Thank you for your help. I check the codes roughly, the method in the codes is to use MatSolve() to solve AX=B in a loop. I also consider such a method.<br> However, I am wondering whether it is slower than the method that directly solves AX=B? thanks again.<br>
<br>Regards,<br>Yujie<br><br><div><span class="gmail_quote">On 2/29/08, <b class="gmail_sendername">Barry Smith</b> <<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>> wrote:</span><blockquote class="gmail_quote" style="margin-top: 0; margin-right: 0; margin-bottom: 0; margin-left: 0; margin-left: 0.80ex; border-left-color: #cccccc; border-left-width: 1px; border-left-style: solid; padding-left: 1ex">
<br> Please edit the file src/mat/interface/matrix.c and remove the<br> function MatMatSolve(). Replace it with the following two functions,<br> then run "make lib shared" in that directory. Please let us know at <a href="mailto:petsc-maint@mcs.anl.gov">petsc-maint@mcs.anl.gov</a><br>
if it crashes or produces incorrect<br> results.<br><br> Barry<br><br><br><br> #undef __FUNCT__<br> #define __FUNCT__ "MatMatSolve_Basic"<br> PetscErrorCode PETSCMAT_DLLEXPORT MatMatSolve_Basic(Mat A,Mat B,Mat X)<br>
{<br> PetscErrorCode ierr;<br> Vec b,x;<br> PetscInt m,N,i;<br> PetscScalar *bb,*xx;<br><br> PetscFunctionBegin;<br> ierr = MatGetArray(B,&bb);CHKERRQ(ierr);<br> ierr = MatGetArray(X,&xx);CHKERRQ(ierr);<br>
ierr = MatGetLocalSize(B,&m,PETSC_NULL);CHKERRQ(ierr); /* number<br> local rows */<br> ierr = MatGetSize(B,PETSC_NULL,&N);CHKERRQ(ierr); /* total<br> columns in dense matrix */<br> ierr = VecCreateMPIWithArray(A-<br>
>hdr.comm,m,PETSC_DETERMINE,PETSC_NULL,&b);CHKERRQ(ierr);<br> ierr = VecCreateMPIWithArray(A-<br> >hdr.comm,m,PETSC_DETERMINE,PETSC_NULL,&x);CHKERRQ(ierr);<br> for (i=0; i<N; i++) {<br> ierr = VecPlaceArray(b,bb + i*m);CHKERRQ(ierr);<br>
ierr = VecPlaceArray(x,xx + i*m);CHKERRQ(ierr);<br> ierr = MatSolve(A,b,x);CHKERRQ(ierr);<br> }<br> ierr = VecDestroy(b);CHKERRQ(ierr);<br> ierr = VecDestroy(x);CHKERRQ(ierr);<br> PetscFunctionReturn(0);<br>
}<br><br> #undef __FUNCT__<br> #define __FUNCT__ "MatMatSolve"<br> /*@<br> MatMatSolve - Solves A X = B, given a factored matrix.<br><br> Collective on Mat<br><br> Input Parameters:<br> + mat - the factored matrix<br>
- B - the right-hand-side matrix (dense matrix)<br><br> Output Parameter:<br> . B - the result matrix (dense matrix)<br><br> Notes:<br> The matrices b and x cannot be the same. I.e., one cannot<br> call MatMatSolve(A,x,x).<br>
<br> Notes:<br> Most users should usually employ the simplified KSP interface for<br> linear solvers<br> instead of working directly with matrix algebra routines such as<br> this.<br> See, e.g., KSPCreate(). However KSP can only solve for one vector<br>
(column of X)<br> at a time.<br><br> Level: developer<br><br> Concepts: matrices^triangular solves<br><br> .seealso: MatMatSolveAdd(), MatMatSolveTranspose(),<br> MatMatSolveTransposeAdd(), MatLUFactor(), MatCholeskyFactor()<br>
@*/<br> PetscErrorCode PETSCMAT_DLLEXPORT MatMatSolve(Mat A,Mat B,Mat X)<br> {<br> PetscErrorCode ierr;<br><br> PetscFunctionBegin;<br> PetscValidHeaderSpecific(A,MAT_COOKIE,1);<br> PetscValidType(A,1);<br> PetscValidHeaderSpecific(B,MAT_COOKIE,2);<br>
PetscValidHeaderSpecific(X,MAT_COOKIE,3);<br> PetscCheckSameComm(A,1,B,2);<br> PetscCheckSameComm(A,1,X,3);<br> if (X == B) SETERRQ(PETSC_ERR_ARG_IDN,"X and B must be different<br> matrices");<br> if (!A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored<br>
matrix");<br> if (A->cmap.N != X->rmap.N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat A,Mat<br> X: global dim %D %D",A->cmap.N,X->rmap.N);<br> if (A->rmap.N != B->rmap.N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat A,Mat<br>
B: global dim %D %D",A->rmap.N,B->rmap.N);<br> if (A->rmap.n != B->rmap.n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat A,Mat<br> B: local dim %D %D",A->rmap.n,B->rmap.n);<br> if (!A->rmap.N && !A->cmap.N) PetscFunctionReturn(0);<br>
ierr = MatPreallocated(A);CHKERRQ(ierr);<br><br> ierr = PetscLogEventBegin(MAT_MatSolve,A,B,X,0);CHKERRQ(ierr);<br> if (!A->ops->matsolve) {<br> ierr = PetscInfo1(A,"Mat type %s using basic MatMatSolve",<br>
((PetscObject)A)->type_name);CHKERRQ(ierr);<br> ierr = MatMatSolve_Basic(A,B,X);CHKERRQ(ierr);<br> } else {<br> ierr = (*A->ops->matsolve)(A,B,X);CHKERRQ(ierr);<br> }<br> ierr = PetscLogEventEnd(MAT_MatSolve,A,B,X,0);CHKERRQ(ierr);<br>
ierr = PetscObjectStateIncrease((PetscObject)X);CHKERRQ(ierr);<br> PetscFunctionReturn(0);<br><br>}<br><br> On Feb 29, 2008, at 6:46 PM, Yujie wrote:<br><br> > Hi, everyone<br> ><br> > I am considering to add a parallel MatMatSolve() into PETSc based on<br>
> SuperLu_DIST or Spooles. If I want to use it like current sequential<br> > MatMatSolve() in the application codes, how to do it? Could you give<br> > me some examples about how to add a new function?<br> > thanks a lot.<br>
><br> > Regards,<br> > Yujie<br><br></blockquote></div><br>