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;
}