how to add a parallel MatMatSolve() function?
Barry Smith
bsmith at mcs.anl.gov
Fri Feb 29 21:48:33 CST 2008
The code I provided works for any LU solver; the way the PETSc code
is written you can customize
a routine for any specific matrix format, like the PETSc SuperLU_Dist
format. You are certainly
free to try to write a custom one for SuperLU_dist() (see
MatMatSolve_SeqAIJ() for how to do this).
It is your time, not mine. Personally I'd rather have the computer run
a few more minutes then spend
my time looking at code :-)
Barry
On Feb 29, 2008, at 9:35 PM, Yujie wrote:
> Dear Barry:
>
> I have checked SuperLU_Dist codes. It looks like relative easy to
> write codes for AX=B based on MatSolve(). This is why I ask you the
> above problem.
> how about your advice?
> thanks a lot.
>
> Regards,
> Yujie
>
> On 3/1/08, Barry Smith <bsmith at mcs.anl.gov> wrote:
>
> 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/6a3cb797/attachment.htm>
More information about the petsc-users
mailing list