how to add a parallel MatMatSolve() function?

[Barry Smith]

   Some direct solver packages have support for solving directly with  
several right hand sides at the same time.
They could be a bit faster than solving one at a time; maybe 30%  
faster at most, not 10 times faster. What is more
important solving the problem you want to solve in a reasonable time  
or solving the problem a bit faster
after spending several weeks writing the much more complicated code?

   Barry

On Feb 29, 2008, at 8:01 PM, Yujie wrote:

> Dear Barry:
>
> 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.
> However, I am wondering whether it is slower than the method that  
> directly solves AX=B? thanks again.
>
> Regards,
> Yujie
>
> On 2/29/08, Barry Smith <bsmith at mcs.anl.gov> wrote:
>
>    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
>
>

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20080229/118818f9/attachment.htm>