Thanks a lot:). <br><br>Regards,<br>Yujie<br><br><div><span class="gmail_quote">On 3/1/08, <b class="gmail_sendername">Barry Smith</b> &lt;<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>&gt; 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">
<div><div><br></div>&nbsp;&nbsp;The code I provided works for any LU solver; the way the PETSc code is written you can customize<div>a routine for any specific matrix format, like the PETSc SuperLU_Dist format. You are certainly</div>
<div>free to try to write a custom one for SuperLU_dist() (see MatMatSolve_SeqAIJ() for how to do this).</div><div>It is your time, not mine. Personally I&#39;d rather have the computer run a few more minutes then spend&nbsp;</div>
<div>my time looking at code :-)</div><div><br></div><span class="sg"><div>&nbsp;&nbsp;Barry</div></span><div><span class="e" id="q_1186874f84c8804f_2"><div><br><div><div>On Feb 29, 2008, at 9:35 PM, Yujie wrote:</div><br><blockquote type="cite">
Dear Barry:<br><br>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. <br>how about your advice? <br>thanks a lot.<br><br>Regards,<br>
 Yujie<br><br><div><span class="gmail_quote">On 3/1/08, <b class="gmail_sendername">Barry Smith</b> &lt;<a href="mailto:bsmith@mcs.anl.gov" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)">bsmith@mcs.anl.gov</a>&gt; 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">
<div><div><br></div>&nbsp;&nbsp;Some direct solver packages have support for solving directly with several right hand sides at the same time.<div>They could be a bit faster than solving one at a time; maybe 30% faster at most, not 10 times faster. What is more</div>
<div>important solving the problem you want to solve in a reasonable time or solving the problem a bit faster</div><div>after spending several weeks writing the much more complicated code?</div><div><br></div><span><div> &nbsp;&nbsp;Barry</div>
</span><div><span><div><br><div><div>On Feb 29, 2008, at 8:01 PM, Yujie wrote:</div><br><blockquote type="cite">Dear Barry:<br><br>Thank you for your help. I check the codes roughly,&nbsp;the&nbsp;method&nbsp;in&nbsp;the&nbsp;codes&nbsp;is&nbsp;to&nbsp;use&nbsp;MatSolve()&nbsp;to&nbsp;solve&nbsp;AX=B&nbsp;in&nbsp;a&nbsp;loop.&nbsp;I&nbsp;also&nbsp;consider&nbsp;such&nbsp;a&nbsp;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> &lt;<a href="mailto:bsmith@mcs.anl.gov" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)">bsmith@mcs.anl.gov</a>&gt; 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>&nbsp;&nbsp; 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 &quot;make lib shared&quot; in that directory. Please let us know at <a href="mailto:petsc-maint@mcs.anl.gov" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)">petsc-maint@mcs.anl.gov</a><br>
 &nbsp;&nbsp;if it crashes or produces incorrect<br> results.<br><br>&nbsp;&nbsp;&nbsp;&nbsp;Barry<br><br><br><br> #undef __FUNCT__<br> #define __FUNCT__ &quot;MatMatSolve_Basic&quot;<br> PetscErrorCode PETSCMAT_DLLEXPORT MatMatSolve_Basic(Mat A,Mat B,Mat X)<br>
 {<br>&nbsp;&nbsp; PetscErrorCode ierr;<br>&nbsp;&nbsp; Vec&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;b,x;<br>&nbsp;&nbsp; PetscInt&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; m,N,i;<br>&nbsp;&nbsp; PetscScalar&nbsp;&nbsp;&nbsp;&nbsp;*bb,*xx;<br><br>&nbsp;&nbsp; PetscFunctionBegin;<br>&nbsp;&nbsp; ierr = MatGetArray(B,&amp;bb);CHKERRQ(ierr);<br>&nbsp;&nbsp; ierr = MatGetArray(X,&amp;xx);CHKERRQ(ierr);<br>
 &nbsp;&nbsp; ierr = MatGetLocalSize(B,&amp;m,PETSC_NULL);CHKERRQ(ierr);&nbsp;&nbsp;/* number<br> local rows */<br>&nbsp;&nbsp; ierr = MatGetSize(B,PETSC_NULL,&amp;N);CHKERRQ(ierr);&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; /* total<br> columns in dense matrix */<br>&nbsp;&nbsp; ierr = VecCreateMPIWithArray(A-<br>
 &nbsp;&nbsp;&gt;hdr.comm,m,PETSC_DETERMINE,PETSC_NULL,&amp;b);CHKERRQ(ierr);<br>&nbsp;&nbsp; ierr = VecCreateMPIWithArray(A-<br>&nbsp;&nbsp;&gt;hdr.comm,m,PETSC_DETERMINE,PETSC_NULL,&amp;x);CHKERRQ(ierr);<br>&nbsp;&nbsp; for (i=0; i&lt;N; i++) {<br>&nbsp;&nbsp;&nbsp;&nbsp; ierr = VecPlaceArray(b,bb + i*m);CHKERRQ(ierr);<br>
 &nbsp;&nbsp;&nbsp;&nbsp; ierr = VecPlaceArray(x,xx + i*m);CHKERRQ(ierr);<br>&nbsp;&nbsp;&nbsp;&nbsp; ierr = MatSolve(A,b,x);CHKERRQ(ierr);<br>&nbsp;&nbsp; }<br>&nbsp;&nbsp; ierr = VecDestroy(b);CHKERRQ(ierr);<br>&nbsp;&nbsp; ierr = VecDestroy(x);CHKERRQ(ierr);<br>&nbsp;&nbsp; PetscFunctionReturn(0);<br>
 }<br><br> #undef __FUNCT__<br> #define __FUNCT__ &quot;MatMatSolve&quot;<br> /*@<br>&nbsp;&nbsp;&nbsp;&nbsp;MatMatSolve - Solves A X = B, given a factored matrix.<br><br>&nbsp;&nbsp;&nbsp;&nbsp;Collective on Mat<br><br>&nbsp;&nbsp;&nbsp;&nbsp;Input Parameters:<br> +&nbsp;&nbsp;mat - the factored matrix<br>
 -&nbsp;&nbsp;B - the right-hand-side matrix&nbsp;&nbsp;(dense matrix)<br><br>&nbsp;&nbsp;&nbsp;&nbsp;Output Parameter:<br> .&nbsp;&nbsp;B - the result matrix (dense matrix)<br><br>&nbsp;&nbsp;&nbsp;&nbsp;Notes:<br>&nbsp;&nbsp;&nbsp;&nbsp;The matrices b and x cannot be the same.&nbsp;&nbsp;I.e., one cannot<br>&nbsp;&nbsp;&nbsp;&nbsp;call MatMatSolve(A,x,x).<br>
<br>&nbsp;&nbsp;&nbsp;&nbsp;Notes:<br>&nbsp;&nbsp;&nbsp;&nbsp;Most users should usually employ the simplified KSP interface for<br> linear solvers<br>&nbsp;&nbsp;&nbsp;&nbsp;instead of working directly with matrix algebra routines such as<br> this.<br>&nbsp;&nbsp;&nbsp;&nbsp;See, e.g., KSPCreate(). However KSP can only solve for one vector<br>
 (column of X)<br>&nbsp;&nbsp;&nbsp;&nbsp;at a time.<br><br>&nbsp;&nbsp;&nbsp;&nbsp;Level: developer<br><br>&nbsp;&nbsp;&nbsp;&nbsp;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>&nbsp;&nbsp; PetscErrorCode ierr;<br><br>&nbsp;&nbsp; PetscFunctionBegin;<br>&nbsp;&nbsp; PetscValidHeaderSpecific(A,MAT_COOKIE,1);<br>&nbsp;&nbsp; PetscValidType(A,1);<br>&nbsp;&nbsp; PetscValidHeaderSpecific(B,MAT_COOKIE,2);<br>
 &nbsp;&nbsp; PetscValidHeaderSpecific(X,MAT_COOKIE,3);<br>&nbsp;&nbsp; PetscCheckSameComm(A,1,B,2);<br>&nbsp;&nbsp; PetscCheckSameComm(A,1,X,3);<br>&nbsp;&nbsp; if (X == B) SETERRQ(PETSC_ERR_ARG_IDN,&quot;X and B must be different<br> matrices&quot;);<br>&nbsp;&nbsp; if (!A-&gt;factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,&quot;Unfactored<br>
 matrix&quot;);<br>&nbsp;&nbsp; if (A-&gt;cmap.N != X-&gt;rmap.N) SETERRQ2(PETSC_ERR_ARG_SIZ,&quot;Mat A,Mat<br> X: global dim %D %D&quot;,A-&gt;cmap.N,X-&gt;rmap.N);<br>&nbsp;&nbsp; if (A-&gt;rmap.N != B-&gt;rmap.N) SETERRQ2(PETSC_ERR_ARG_SIZ,&quot;Mat A,Mat<br>
 B: global dim %D %D&quot;,A-&gt;rmap.N,B-&gt;rmap.N);<br>&nbsp;&nbsp; if (A-&gt;rmap.n != B-&gt;rmap.n) SETERRQ2(PETSC_ERR_ARG_SIZ,&quot;Mat A,Mat<br> B: local dim %D %D&quot;,A-&gt;rmap.n,B-&gt;rmap.n);<br>&nbsp;&nbsp; if (!A-&gt;rmap.N &amp;&amp; !A-&gt;cmap.N) PetscFunctionReturn(0);<br>
 &nbsp;&nbsp; ierr = MatPreallocated(A);CHKERRQ(ierr);<br><br>&nbsp;&nbsp; ierr = PetscLogEventBegin(MAT_MatSolve,A,B,X,0);CHKERRQ(ierr);<br>&nbsp;&nbsp; if (!A-&gt;ops-&gt;matsolve) {<br>&nbsp;&nbsp;&nbsp;&nbsp; ierr = PetscInfo1(A,&quot;Mat type %s using basic MatMatSolve&quot;,<br>
 ((PetscObject)A)-&gt;type_name);CHKERRQ(ierr);<br>&nbsp;&nbsp;&nbsp;&nbsp; ierr = MatMatSolve_Basic(A,B,X);CHKERRQ(ierr);<br>&nbsp;&nbsp; } else {<br>&nbsp;&nbsp;&nbsp;&nbsp; ierr = (*A-&gt;ops-&gt;matsolve)(A,B,X);CHKERRQ(ierr);<br>&nbsp;&nbsp; }<br>&nbsp;&nbsp; ierr = PetscLogEventEnd(MAT_MatSolve,A,B,X,0);CHKERRQ(ierr);<br>
 &nbsp;&nbsp; ierr = PetscObjectStateIncrease((PetscObject)X);CHKERRQ(ierr);<br>&nbsp;&nbsp; PetscFunctionReturn(0);<br><br>}<br><br> On Feb 29, 2008, at 6:46 PM, Yujie wrote:<br><br> &gt; Hi, everyone<br> &gt;<br> &gt; I am considering to add a parallel MatMatSolve() into PETSc based on<br>
 &gt; SuperLu_DIST or Spooles. If I want to use it like current sequential<br> &gt; MatMatSolve() in the application codes, how to do it? Could you give<br> &gt; me some examples about how to add a new function?<br> &gt; thanks a lot.<br>
 &gt;<br> &gt; Regards,<br> &gt; Yujie<br><br></blockquote></div><br></blockquote></div><br></div></span></div></div></blockquote></div><br></blockquote></div><br></div></span></div></div></blockquote></div><br>