static char help[] = "Solves a tridiagonal linear system with KSP.\n\n";
/*
 * resolution d'un systeme linaire
*/

#include "petscvec.h"
#include "petscmat.h"
#include "petscksp.h"
#include "petsctime.h"

int main(int argc,char **args)
{
  KSP ksp;
  PC  pc; 
  PetscErrorCode ierr;
  Mat mat;
  Vec x,rhs,sol;  //a calculer, rhs, solution exacte
  PetscReal norm,u_exact;
  PetscInt i,Nx,Ny,column[3],start,end,istart,iend,debug,its;
  PetscScalar neg_one=-1.0,value[3];
  PetscLogDouble v1,v2,elapsed_time;


  PetscInitialize(&argc,&args,(char*)0,help);

  //Taille du probleme
  // Attention Nx=Ny sinon pas de sens
  PetscOptionsGetInt(PETSC_NULL,"-nx",&Nx,PETSC_NULL);
  PetscOptionsGetInt(PETSC_NULL,"-ny",&Ny,PETSC_NULL);
  PetscOptionsGetInt(PETSC_NULL,"-debug",&debug,PETSC_NULL);

/*
  - Creation de la matrice du probleme
*/
  MatCreate(PETSC_COMM_WORLD,&mat);
  MatSetSizes(mat,PETSC_DECIDE,PETSC_DECIDE,Nx,Ny);
  MatSetType(mat,MATAIJ); 
  //MatSetType(mat,MATAIJMUMPS);   // apres aijmumps
  MatSetFromOptions(mat);

  /* on remplit le vecteur suivant sont decoupage */
  MatGetOwnershipRange(mat,&start,&end);
  istart = start; if (start == 0) istart = 1;
  iend = end; if (iend == Nx) iend = Nx-1;
  value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
  
  for (i=istart; i<iend; i++) {
    column[0] = i-1; column[1] = i; column[2] = i+1;
    MatSetValues(mat,1,&i,3,column,value,INSERT_VALUES);
  }
  //On assemble les bords
  //1ere ligne :
  if (start == 0) {
    i=0;
    column[0]=0; column[1]=1;
    value[0]=1; value[1]=0;
    MatSetValues(mat,1,&i,2,column,value,INSERT_VALUES); 
  }
  //derniere ligne
  //attention erreur classique : Ca commence a 0, pas 1 !!!
  //Ca decale tout !!!
  if (end == Nx) {
    i=Nx-1;
    column[0]=Nx-2; column[1]=Nx-1;
    value[0]=0; value[1]=1;
    MatSetValues(mat,1,&i,2,column,value,INSERT_VALUES); 
  }

  MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
  MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);

/*
  affichage de la matrice
*/
  if (debug) MatView(mat,PETSC_VIEWER_STDOUT_WORLD);

/*
  - Creation du rhs
*/
  VecCreate(PETSC_COMM_WORLD,&rhs);
  VecSetSizes(rhs,PETSC_DECIDE,Nx);
  VecSetType(rhs,VECMPI);
  VecSetFromOptions(rhs);

  /* on remplit le vecteur suivant sont decoupage */
  VecGetOwnershipRange(rhs,&start,&end);
  if (start == 0)  VecSetValue(rhs,start,0,INSERT_VALUES);
  if (end == Nx)  VecSetValue(rhs,end-1,1,INSERT_VALUES);

  VecAssemblyBegin(rhs);
  VecAssemblyEnd(rhs);

/*
  View the vector
*/
  if (debug) VecView(rhs,PETSC_VIEWER_STDOUT_WORLD);

/*
  - On resoud
*/
  VecDuplicate(rhs,&x);
  VecDuplicate(rhs,&sol);

  KSPCreate(PETSC_COMM_WORLD,&ksp);

  KSPSetOperators(ksp,mat,mat,DIFFERENT_NONZERO_PATTERN);

/* Le preconditionneur */
/*  KSPGetPC(ksp,&pc); */
/*  PCSetType(pc,PCJACOBI); */

  KSPSetTolerances(ksp,1.e-12,PETSC_DEFAULT,PETSC_DEFAULT,100000);
  KSPSetFromOptions(ksp);

  PetscGetTime(&v1);
  KSPSolve(ksp,rhs,x);
  PetscGetTime(&v2);
  elapsed_time = v2 - v1;
  PetscPrintf(PETSC_COMM_WORLD,"\nTemps de resolution = %G s\n",elapsed_time); 

/*
  - On affiche la solution
*/
  if (debug) VecView(x,PETSC_VIEWER_STDOUT_WORLD);

/*
  - Calcul de l'erreur
*/
  /* solution exact */

  VecGetOwnershipRange(sol,&start,&end);
  for (i=start; i<end; i++) {
    u_exact=(i*1.0)/(Nx-1);
    VecSetValue(sol,i,u_exact,INSERT_VALUES);
  }

  VecAssemblyBegin(sol);
  VecAssemblyEnd(sol);

  if (start == 0) {
    PetscInt  val[1];
    PetscReal  res[1];
    val[0]=3;
    VecGetValues(x,1,val,res);
    u_exact=res[0]-(3*1.0)/(Nx-1);
    printf("erreur = %12.7e\n",u_exact);
  }


  VecAXPY(x,neg_one,sol);
  VecNorm(x,NORM_2,&norm);
  KSPGetIterationNumber(ksp,&its);
  PetscPrintf(PETSC_COMM_WORLD,"Norme de l'erreur %A, Iterations %D\n",
    norm,its);

  PetscFree(mat);
  PetscFree(x);
  PetscFree(rhs);
  PetscFree(sol);
  
  PetscFinalize();
  return 0;
}
 
