Barry,
I've attached a small example that replicates this behavior.
John
On Wed, Apr 13, 2016 at 9:25 AM, Barry Smith <[email protected]> wrote:
>
> John,
>
> We'll need the full code to debug this. It doesn't see this should
> happen.
>
> Barry
>
> > On Apr 13, 2016, at 7:11 AM, John Mousel <[email protected]> wrote:
> >
> > I have a colleague who is trying to run some simple test cases. He
> noticed if he provides the solution vector as the exact solution with a
> non-zero guess tagged, that the solution vector is zeroed upon return after
> 0 iterations. The system is just a 5x5 matrix corresponding to a heat
> equation with Dirichlet BCs. I've copied the KSPView below. He's using
> petsc-3.6.3.
> >
> > John
> >
> > solution before kspsolve:Vec Object: 1 MPI processes
> > type: mpi
> > Process [0]
> > 300
> > 300
> > 300
> > 300
> > 300
> > KSP Object: 1 MPI processes
> > type: bcgs
> > maximum iterations=100
> > tolerances: relative=1e-15, absolute=1e-50, divergence=10000
> > right preconditioning
> > using nonzero initial guess
> > using UNPRECONDITIONED norm type for convergence test
> > PC Object: 1 MPI processes
> > type: ilu
> > ILU: out-of-place factorization
> > 0 levels of fill
> > tolerance for zero pivot 2.22045e-14
> > matrix ordering: natural
> > factor fill ratio given 1, needed 1
> > Factored matrix follows:
> > Mat Object: 1 MPI processes
> > type: seqaij
> > rows=5, cols=5
> > package used to perform factorization: petsc
> > total: nonzeros=25, allocated nonzeros=25
> > total number of mallocs used during MatSetValues calls =0
> > using I-node routines: found 1 nodes, limit used is 5
> > linear system matrix = precond matrix:
> > Mat Object: 1 MPI processes
> > type: seqaij
> > rows=5, cols=5
> > total: nonzeros=25, allocated nonzeros=25
> > total number of mallocs used during MatSetValues calls =0
> > using I-node routines: found 1 nodes, limit used is 5
> > solution after kspsolve:Vec Object: 1 MPI processes
> > type: mpi
> > Process [0]
> > 0
> > 0
> > 0
> > 0
> > 0
>
>
static char help[] = "Solves a tridiagonal linear system with KSP.\n\n";
#include <petscksp.h>
#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **args)
{
Vec x, b; /* RHS, exact solution */
Mat A; /* linear system matrix */
KSP ksp; /* linear solver context */
PC pc; /* preconditioner context */
PetscInt i,n = 10,col[3],its;
PetscMPIInt size;
PetscScalar value[3];
PetscReal final_resid = 0.;
PetscInitialize(&argc,&args,(char*)0,help);
MPI_Comm_size(PETSC_COMM_WORLD,&size);
PetscPrintf(PETSC_COMM_WORLD,"\nargc after PetscInitialize = %D\n", argc);
PetscPrintf(PETSC_COMM_WORLD,"\nargv PetscInitialize = %s\n", args[0]);
VecCreate(PETSC_COMM_WORLD,&x);
VecSetSizes(x,PETSC_DECIDE,n);
VecSetFromOptions(x);
VecDuplicate(x,&b);
MatCreate(PETSC_COMM_WORLD,&A);
MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);
MatSetFromOptions(A);
MatSetUp(A);
value[0] = 0;
value[1] = 1;
value[2] = 0;
for (i=2; i < n-1; ++i)
{
col[0] = i-1;
col[1] = i;
col[2] = i+1;
MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
}
for (i = 0; i < 2; ++i)
{
col[0] = i;
col[1] = i + 1;
value[0] = 1.0;
value[1] = 0.0;
MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
}
i = n - 1;
col[0] = i - 2;
col[1] = i;
value[0] = 0;
value[1] = 1.;
MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
PetscPrintf(PETSC_COMM_WORLD,"\nView matrix _A values\n");
MatView(A, PETSC_VIEWER_STDOUT_WORLD);
KSPCreate(PETSC_COMM_WORLD,&ksp);
KSPSetOperators(ksp,A,A);
KSPSetType(ksp, KSPBCGS);
KSPSetTolerances(ksp,1e-15,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
KSPSetNormType(ksp, KSP_NORM_UNPRECONDITIONED);
KSPGetPC(ksp,&pc);
PCSetType(pc,PCILU);
KSPSetFromOptions(ksp);
const PetscScalar b_init = 1.;
VecSet(b,b_init);
PetscPrintf(PETSC_COMM_WORLD,"\nView vector _b values\n");
VecView(b, PETSC_VIEWER_STDOUT_WORLD);
const PetscScalar x_init = 1.;
VecSet(x, x_init);
KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
PetscPrintf(PETSC_COMM_WORLD,"\nView vector _x values before KSPSolve\n");
VecView(x, PETSC_VIEWER_STDOUT_WORLD);
KSPSolve(ksp,b,x);
KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);
KSPGetResidualNorm(ksp, &final_resid);
PetscPrintf(PETSC_COMM_WORLD,"norm residual %f\n", (double) final_resid);
KSPGetIterationNumber(ksp,&its);
PetscPrintf(PETSC_COMM_WORLD,"Iterations %D\n",its);
PetscPrintf(PETSC_COMM_WORLD,"\nView vector _x values after KSPSolve\n");
VecView(x,PETSC_VIEWER_STDOUT_WORLD);
VecDestroy(&x);
VecDestroy(&b); MatDestroy(&A);
KSPDestroy(&ksp);
PetscFinalize();
return 0;
}
