[petsc-users] problem in calling KSPSetComputeRHS and KSPSolve multiple times to solve different RHS.

[Barry Smith]

On Sep 13, 2013, at 6:21 AM, Bishesh Khanal <bisheshkh at gmail.com> wrote:

> Dear all,
> I'm trying to reuse an operator for the linear system multiple times with different RHS. 
> But it seems that the RHS is not changed in the KSP object. Here is the code to illustrate what I'm doing:
> 
> Vec mX[3], mB[3];   //Solving Ax=b for 3 different RHS.
> ...
> ierr = KSPSetComputeOperators(mKsp,computeMatrix,this);CHKERRQ(ierr);
> ierr = KSPSetDM(mKsp,mDa);CHKERRQ(ierr);

   Move the line up above the previous line

> ierr = KSPSetFromOptions(mKsp);CHKERRQ(ierr);
> 
> for(int i = 0; i<3; ++i) {
>         mBComponent = i;   //sets different component for RHS.
>         ierr = KSPSetComputeRHS(mKsp,computeRhs,this);CHKERRQ(ierr);

   Put the above line in the setup code above, not in the loop

>         ierr = KSPSetDM(mKsp,mDa);CHKERRQ(ierr);

   Remove the line above

>         ierr = KSPSolve(mKsp,NULL,NULL);CHKERRQ(ierr);
>         ierr = KSPGetSolution(mKsp,&mX[i]);CHKERRQ(ierr);
>         ierr = KSPGetRhs(mKsp,&mB[i]);CHKERRQ(ierr);
>  }

   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.

>     
> 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.
> What am I doing wrong here ?
> 
> 
> Thanks,
> Bishesh
>