[petsc-users] successive linear solves

[Umut Tabak]

Dear all,

I was trying to solve two linear systems in the same ksp context while 
setting the operators with KSPSetOperators. But I am getting a Memory 
corruption error. What  I try to do is to set two linear operators and 
rhs vectors and get the solutions, say

for solA, xa,

(Ka-lambs*Ma)*xa = -lambs*KsaT*modeS (I extract lambs from a vector of 
eigenvalues and modeS from a matrix of eigenvectors)

and for solS, xs

(Ks-lamba*Ms)*xs = -Ksa*modeP (similarly lamba is extracted from a 
vector and the modeP also from a matrix of vectors as columns)

I am using MatDuplicate functions to update the operator matrices at 
each iteration, there is not a problem in the first solve, at least no 
errors, but on the second one I am getting the Memory corruption error 
mentioned above. Also one more question is on where to use 
KSPSetFromOptions on the iterations. If I use KSPSetFromOptions on the 
same iteration twice(before solve)I get the same error from this 
function(this was not also reasonable to me). Section of the code is 
below(a bit crappy however should give the idea on what I tried to do). 
I am sure this is a simple problem but could not find the reason.

   KSP ksp;
   KSPCreate(PETSC_COMM_WORLD, &ksp);
   for(int j=0; j!=nR; j++)
   {
     //KSPCreate(PETSC_COMM_WORLD, &kspA);
     // duplicate the matrix
     MatDuplicate(Ka, MAT_COPY_VALUES, &OpA);
     //
     MatAXPY(OpA, -eSolverS.getRoot(j), Ma, SAME_NONZERO_PATTERN);
     MatGetColumnVector(eigVecsStr, modeS, j);
     // create the rhs for the operation
     MatMult(KsaT, modeS, rhsA);
     VecScale(rhsA, -eSolverS.getRoot(j));
     cout << "Solving the linear system for pressure side corrections" 
<< endl;
     KSPSetOperators(ksp, OpA, OpA, DIFFERENT_NONZERO_PATTERN);
     //KSPGetPC(kspA, &pcA);
     //PCSetType(pcA, PCJACOBI);
     //KSPSetTolerances(ksp, 1e-7, 
PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
     //KSPSetFromOptions(kspA);
     // solve
     KSPSolve(ksp,rhsA,solA);
     corA.assignSubVector(solA, 1, ra, j);
     //
     MatDuplicate(Ks, MAT_COPY_VALUES, &OpS);
     MatAXPY(OpS, -eSolverA.getRoot(j), Ms, SAME_NONZERO_PATTERN);
     MatGetColumnVector(eigVecsAco, modeA, j);
     // create the rhs for the operation
     MatMult(Ksa, modeA, rhsS);
     VecScale(rhsS, -1.0);
     cout << "Solving the linear system for structure side corrections" 
<< endl;
     //KSPCreate(PETSC_COMM_WORLD, &kspS);
     KSPSetOperators(ksp, OpS, OpS, DIFFERENT_NONZERO_PATTERN);
     //KSPGetPC(kspS, &pcS);
     //PCSetType(pcS, PCJACOBI);
     //KSPSetTolerances(kspS, 1e-7, 
PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
     //KSPSetFromOptions(kspS);
     // solve
     KSPSolve(ksp,rhsS,solS);
     //int r,m;
     //VecGetSize(rhsS, &r);
     //VecGetSize(solS, &m);
     //cout << r << " " << m << endl;
     corS.assignSubVector(solS, 1, rs, j);
     cout << "Solved  the linear system for structure corrections" << endl;
   }
   corA.finalize(); corS.finalize();

Best,
Umut