[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: about Dense Matrix Re: several problems about codes Re: how to add a parallel MatMatSolve() function?
- To: "Matthew Knepley" <knepley@xxxxxxxxx>
- Subject: Re: about Dense Matrix Re: several problems about codes Re: how to add a parallel MatMatSolve() function?
- From: Yujie <recrusader@xxxxxxxxx>
- Date: Sat, 1 Mar 2008 11:12:12 -0800
- Cc: petsc-users@xxxxxxxxxxx
- Dkim-signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:received:received:message-id:date:from:to:subject:cc:in-reply-to:mime-version:content-type:references; bh=xJMrWXVbEF8HTp/J+UYVKAUS+t57wr+cWJdFnwBXpEU=; b=VNm+cnzXodptW3VKcgGjzVbAQe8Wsgl7I810/1uC7oFId4u0fieaO90VUSH3Formy/T6UoguH79a2SE23HTKQt2NLfJOcE3dqAWUIA3SOpuM1FJdyxOqulSO4iq/6nGeQ+XqVh/dM067vJgt6P2qod3/wORG8xbe35yLfhfjT/8=
- Domainkey-signature: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=message-id:date:from:to:subject:cc:in-reply-to:mime-version:content-type:references; b=gmuJE6qE6Xb61iJcs+TNcMl0DaLOTTbxbGI0VSjfBu1CAJ2nm2IgOLwim0BoYZkj5EY3tQF+CunyHkRl2IQVkiXqhb0TJCCviZyLCZ9Fp4dhPVVyGt/z9uIOGsdK4ExsJSgJwmorJ8zuyF07C2fQZaUU/uDZ3JoK+QjuY5MNjW4=
- In-reply-to: <a9f269830803011041s5921bd53jf84d59fa88dae016@mail.gmail.com>
- References: <7ff0ee010803011038q6d0fb498ne8d65d97f5711992@mail.gmail.com> <a9f269830803011041s5921bd53jf84d59fa88dae016@mail.gmail.com>
- Reply-to: petsc-users@xxxxxxxxxxx
- Sender: owner-petsc-users@xxxxxxxxxxx
However, you have realized A*B(A is AIJ, B is dense matrix) to my knowledge.
PETSc is the Portable, Extensible Toolkit for Scientific computation, not the Sparse Portable, Extensible Toolkit for Scientific computation :).
It is easy to realize A*B(A is dense matrix, B is AIJ)?
Thanks,
Yujie
On 3/1/08, Matthew Knepley <knepley@xxxxxxxxx> wrote:
If you want dense matrices, I would suggest PLAPACK/FLAME instead of PETSc. We
specialize in sparse matrices.
Thanks,
Matt
On Sat, Mar 1, 2008 at 12:38 PM, Yujie <recrusader@xxxxxxxxx> wrote:
> Dear Barry:
>
> I just find MatMatMult() doesn't support Dense Matrix. I want to do A*B (A
> is dense matrix, B is AIJ). I think a possible method is to convert A to
> AIJ, however, it will take much more memory than using dense matrix, I
> think. Because AIJ needs to save positions and corresponding values. Do you
> have any good idea? Thanks a lot.
>
> Regards,
> Yujie
>
> On 3/1/08, Barry Smith <bsmith@xxxxxxxxxxx> wrote:
> >
> >
> >
> >
> > On Mar 1, 2008, at 2:01 AM, Yujie wrote:
> >
> > Dear Barry:
> >
> > I can't understand several problems as follows about codes. Could you give
> me some explanations? thanks a lot.
> >
> > 1. About MatGetArray()
> > I check the codes of PETSc. I find only matrix types, SeqAij and MPIDense,
> support this fnction. If I want to use SuperLU_dist or Spooles. I need to
> convert matrix A into corresponding type
> > before I call MatMatSolve(). If it is, MatGetArray() doesn't work, I
> think.
> >
> > MatGetArray() is called on B and X, it is not called on A. B and A must
> be dense or the whole thing doesn't make
> > sense, since X will always end up being dense.
> >
> >
> >
> >
> > 2. A->hdr.comm should be A->comm?
> >
> > Yes, for your version of PETSc
> >
> >
> >
> > 3.
> > 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);
> > }
> >
> > should be
> >
> > for (i=0; i<N; i++) {
> > ierr = VecPlaceArray(b,bb + i*m);CHKERRQ(ierr);
> > ierr = MatSolve(A,b,x);CHKERRQ(ierr);
> > ierr = VecPlaceArray(xx + i*m,x);CHKERRQ(ierr);
> > }
> > MatRestoreArray(B,&xx);
> >
> > First replace b with bb+i*m, and solve Ax=b, final replace xx+i*m with x.
> >
> > No, you misunderstand what VecPlaceArray() does. The code above will solve
> the columns in the wrong place.
> >
> >
> > Yes, I did forget the MatRestoreArray(), thanks for pointing that out.
> >
> >
> > Barry
> >
> >
> >
> >
> > After finishing the loop, the code should call MatRestoreArray() to
> restore matrix B.
> >
> > Regards,
> > Yujie
> >
> >
> >
> > On 3/1/08, Barry Smith <bsmith@xxxxxxxxxxx> 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@xxxxxxxxxxx
> > > 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
> > >
> > >
> >
> >
> >
>
>
--
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which
their experiments lead.
-- Norbert Wiener