<div dir="ltr"><br><div class="gmail_extra"><br><br><div class="gmail_quote">On Fri, Sep 13, 2013 at 2:47 PM, Barry Smith <span dir="ltr"><<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div class="im"><br>
On Sep 13, 2013, at 6:21 AM, Bishesh Khanal <<a href="mailto:bisheshkh@gmail.com">bisheshkh@gmail.com</a>> wrote:<br>
<br>
> Dear all,<br>
> I'm trying to reuse an operator for the linear system multiple times with different RHS.<br>
> But it seems that the RHS is not changed in the KSP object. Here is the code to illustrate what I'm doing:<br>
><br>
> Vec mX[3], mB[3]; //Solving Ax=b for 3 different RHS.<br>
> ...<br>
> ierr = KSPSetComputeOperators(mKsp,computeMatrix,this);CHKERRQ(ierr);<br>
> ierr = KSPSetDM(mKsp,mDa);CHKERRQ(ierr);<br>
<br>
</div> Move the line up above the previous line<br>
<div class="im"><br>
> ierr = KSPSetFromOptions(mKsp);CHKERRQ(ierr);<br>
><br>
> for(int i = 0; i<3; ++i) {<br>
> mBComponent = i; //sets different component for RHS.<br>
> ierr = KSPSetComputeRHS(mKsp,computeRhs,this);CHKERRQ(ierr);<br>
<br>
</div> Put the above line in the setup code above, not in the loop<br>
<br>
> ierr = KSPSetDM(mKsp,mDa);CHKERRQ(ierr);<br>
<br>
Remove the line above<br>
<div class="im"><br>
> ierr = KSPSolve(mKsp,NULL,NULL);CHKERRQ(ierr);<br>
> ierr = KSPGetSolution(mKsp,&mX[i]);CHKERRQ(ierr);<br>
> ierr = KSPGetRhs(mKsp,&mB[i]);CHKERRQ(ierr);<br>
> }<br>
<br>
</div> After you make this changes it should work. If it does not then run with -start_in_debugger noxterm and put a break point in computeRHS, make sure it is being called.<br></blockquote><div>Thanks, I did as you suggested but there is still a problem. <br>
I had kept KSPComputeRHS within the loop so that it gets called to update the RHS based on the value of mBComponent, using computeRhs function.<br></div><div>Now I kept it outside the loop in the setup code as you suggested and while doing debugging, I see that the program does go inside the computeRhs function three times even in this case. <br>
However when I see the results with the following code, it still has the same rhs values for all three cases.<br>ierr = PetscObjectSetName((PetscObject)mB[0],"rhsX");CHKERRQ(ierr);<br> ierr = VecView(mB[0],viewer);CHKERRQ(ierr);<br>
ierr = PetscObjectSetName((PetscObject)mB[1],"rhsY");CHKERRQ(ierr);<br> ierr = VecView(mB[1],viewer);CHKERRQ(ierr);<br> ierr = PetscObjectSetName((PetscObject)mB[2],"rhsZ");CHKERRQ(ierr);<br>
ierr = VecView(mB[2],viewer);CHKERRQ(ierr);<br><br>Now I am getting confused with two things:<br></div><div>1) It seems that computeRhs function will be called every time we call KspSolve. If this is the case, is it the same with KSPSetComputeOperators too ? Because I do not want the computeMatrix function to be called repeatedly to assemble the same matrix operator three times. Or there must be sth I'm missing here.<br>
<br></div><div>2) The computeRhs is called three times with different values of mBComponent, but I can't figure out why the rhs vectors in mB[] array are same. Here is my computeRhs function:<br><br>PetscErrorCode HarmCurvRegSolver::computeRhs(KSP ksp, Vec b, void *ctx)<br>
{<br> PetscFunctionBeginUser;<br> HarmCurvRegSolver *context = (HarmCurvRegSolver*)ctx;<br><br> PetscErrorCode ierr;<br> DM da;<br> ierr = KSPGetDM(ksp,&da);CHKERRQ(ierr);<br> DMDALocalInfo info;<br>
ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);<br><br> PetscScalar ***rhs;<br> ierr = DMDAVecGetArray(da,b,&rhs);CHKERRQ(ierr);<br><br> //The RHS should be provided by the model!<br> for (PetscInt k = info.zs; k<info.zs+<a href="http://info.zm">info.zm</a>; ++k) {<br>
for (PetscInt j = info.ys; j < info.ys+info.ym; ++j) {<br> for (PetscInt i = info.xs; i < info.xs+info.xm; ++i) {<br> rhs[k][j][i] = context->getProblem()->getRhsAtPosition(i,j,k,context->getCurrentRhsComponent()); <br>
//I've checked well the getRhsAtPosition(...) function! It does give different values of the 4th argument.<br>//and getCurrentRhsComponent() just returns the value of mBComponent.<br> }<br> }<br> }<br>
<br> ierr = DMDAVecRestoreArray(da,b,&rhs);CHKERRQ(ierr);<br> ierr = VecAssemblyBegin(b);CHKERRQ(ierr);<br> ierr = VecAssemblyEnd(b);CHKERRQ(ierr);<br><br> PetscFunctionReturn(0);<br><br>}<br><br></div><div>
<br></div><div> <br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<div class=""><div class="h5"><br>
><br>
> What I expect is to have three different solutions in mX array for the different Rhs computed by computeRhs function. However, the result when I see using VecView is that the Rhs mB[0], mB[1] and mB[2] are all identical (which should not have been the case). And mX also has three identical vectors as solution.<br>
> What am I doing wrong here ?<br>
><br>
><br>
> Thanks,<br>
> Bishesh<br>
><br>
<br>
</div></div></blockquote></div><br></div></div>