However, you have realized A*B(A is AIJ, B is dense matrix) to my knowledge. <br>PETSc is the Portable, Extensible Toolkit for Scientific computation, not the Sparse Portable, Extensible Toolkit for Scientific computation :). <br>
It is easy to realize A*B(A is dense matrix, B is AIJ)?<br> <br>Thanks,<br>Yujie<br><br><div><span class="gmail_quote">On 3/1/08, <b class="gmail_sendername">Matthew Knepley</b> <<a href="mailto:knepley@gmail.com">knepley@gmail.com</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">
If you want dense matrices, I would suggest PLAPACK/FLAME instead of PETSc. We<br> specialize in sparse matrices.<br><br> Thanks,<br><br> Matt<br><br><br> On Sat, Mar 1, 2008 at 12:38 PM, Yujie <<a href="mailto:recrusader@gmail.com">recrusader@gmail.com</a>> wrote:<br>
> Dear Barry:<br> ><br> > I just find MatMatMult() doesn't support Dense Matrix. I want to do A*B (A<br> > is dense matrix, B is AIJ). I think a possible method is to convert A to<br> > AIJ, however, it will take much more memory than using dense matrix, I<br>
> think. Because AIJ needs to save positions and corresponding values. Do you<br> > have any good idea? Thanks a lot.<br> ><br> > Regards,<br> > Yujie<br> ><br> > On 3/1/08, Barry Smith <<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>> wrote:<br>
> ><br> > ><br> > ><br> > ><br> > > On Mar 1, 2008, at 2:01 AM, Yujie wrote:<br> > ><br> > > Dear Barry:<br> > ><br> > > I can't understand several problems as follows about codes. Could you give<br>
> me some explanations? thanks a lot.<br> > ><br> > > 1. About MatGetArray()<br> > > I check the codes of PETSc. I find only matrix types, SeqAij and MPIDense,<br> > support this fnction. If I want to use SuperLU_dist or Spooles. I need to<br>
> convert matrix A into corresponding type<br> > > before I call MatMatSolve(). If it is, MatGetArray() doesn't work, I<br> > think.<br> > ><br> > > MatGetArray() is called on B and X, it is not called on A. B and A must<br>
> be dense or the whole thing doesn't make<br> > > sense, since X will always end up being dense.<br> > ><br> > ><br> > ><br> > ><br> > > 2. A->hdr.comm should be A->comm?<br>
> ><br> > > Yes, for your version of PETSc<br> > ><br> > ><br> > ><br> > > 3.<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> > ><br> > > should be<br> > ><br> > > for (i=0; i<N; i++) {<br> > > ierr = VecPlaceArray(b,bb + i*m);CHKERRQ(ierr);<br>
> > ierr = MatSolve(A,b,x);CHKERRQ(ierr);<br> > > ierr = VecPlaceArray(xx + i*m,x);CHKERRQ(ierr);<br> > > }<br> > > MatRestoreArray(B,&xx);<br> > ><br> > > First replace b with bb+i*m, and solve Ax=b, final replace xx+i*m with x.<br>
> ><br> > > No, you misunderstand what VecPlaceArray() does. The code above will solve<br> > the columns in the wrong place.<br> > ><br> > ><br> > > Yes, I did forget the MatRestoreArray(), thanks for pointing that out.<br>
> ><br> > ><br> > > Barry<br> > ><br> > ><br> > ><br> > ><br> > > After finishing the loop, the code should call MatRestoreArray() to<br> > restore matrix B.<br> > ><br>
> > Regards,<br> > > Yujie<br> > ><br> > ><br> > ><br> > > On 3/1/08, Barry Smith <<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>> wrote:<br> > > ><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<br> > <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> > > ><br> > ><br> > ><br> > ><br> ><br> ><br>
<br><br><br><br>--<br> What most experimenters take for granted before they begin their<br> experiments is infinitely more interesting than any results to which<br> their experiments lead.<br> -- Norbert Wiener<br></blockquote>
</div><br>