static char help[] = "Performs Matrix Vector product\n";

#include <stdio.h>
#include <stdlib.h>
#include "mpi.h"
#include <math.h>
#include "petsc.h"
#include "petscvec.h"
#include "petscmat.h"
#include "petscviewer.h"

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **argv)
{
    PetscErrorCode ierr;
    Mat AA;
    Vec bb,xx;
    PetscScalar val,*bb_v,*xx_v;
    PetscInt indi,indj;
    PetscViewer view_out;

    int n;
    int i,j,istart,iend;
    int nlocal;
    double *fx;
    int npes,myrank;

    /* Get information about the communicator */
    PetscInitialize(&argc,&argv,(char *)0,help);
    ierr=MPI_Comm_rank(PETSC_COMM_WORLD,&myrank);
    ierr=MPI_Comm_size(PETSC_COMM_WORLD,&npes);

    /* Set dimensions */
    nlocal=100000;
    n=nlocal*npes;

    /* Allocate memory */
    ierr=MatCreateMPIAIJ(PETSC_COMM_WORLD,nlocal,nlocal,n,n,3,PETSC_NULL,3,PETSC_NULL,&AA);
    ierr=VecCreateMPI(PETSC_COMM_WORLD,nlocal,n,&bb);
    ierr=VecSet(bb,0);
    ierr=VecCreateMPI(PETSC_COMM_WORLD,nlocal,n,&xx);
    ierr=VecSet(xx,0);

    fx=(double *)malloc(n*sizeof(double));

    /* Set matrix AA */
    if(myrank==0){
      indi=0;
      indj=0;
      val=2;
      ierr=MatSetValues(AA,1,&indi,1,&indj,&val,INSERT_VALUES);
      indi=0;
      indj=1;
      val=-3;
      ierr=MatSetValues(AA,1,&indi,1,&indj,&val,INSERT_VALUES);
      for(i=1;i<nlocal;i++){
         indi=i;
         indj=i-1;
         val=-1;
         ierr=MatSetValues(AA,1,&indi,1,&indj,&val,INSERT_VALUES);
         indi=i;
         indj=i;
         val=2;
         ierr=MatSetValues(AA,1,&indi,1,&indj,&val,INSERT_VALUES);
         indi=i;
         indj=i+1;
         val=-3;
         ierr=MatSetValues(AA,1,&indi,1,&indj,&val,INSERT_VALUES);
      }
    }
    else if(myrank==npes-1){
      for(i=0;i<nlocal-1;i++){
         indi=myrank*nlocal+i;
         indj=myrank*nlocal+i-1;
         val=-1;
         ierr=MatSetValues(AA,1,&indi,1,&indj,&val,INSERT_VALUES);
         indi=myrank*nlocal+i;
         indj=myrank*nlocal+i;
         val=2;
         ierr=MatSetValues(AA,1,&indi,1,&indj,&val,INSERT_VALUES);
         indi=myrank*nlocal+i;
         indj=myrank*nlocal+i+1;
         val=-3;
         ierr=MatSetValues(AA,1,&indi,1,&indj,&val,INSERT_VALUES);
      }
      indi=myrank*nlocal+nlocal-1;
      indj=myrank*nlocal+nlocal-2;
      val=-1;
      ierr=MatSetValues(AA,1,&indi,1,&indj,&val,INSERT_VALUES);
      indi=myrank*nlocal+nlocal-1;
      indj=myrank*nlocal+nlocal-1;
      val=2;
      ierr=MatSetValues(AA,1,&indi,1,&indj,&val,INSERT_VALUES);
    }
    else{
      for(i=0;i<nlocal;i++){
         indi=myrank*nlocal+i;
         indj=myrank*nlocal+i-1;
         val=-1;
         ierr=MatSetValues(AA,1,&indi,1,&indj,&val,INSERT_VALUES);
         indi=myrank*nlocal+i;
         indj=myrank*nlocal+i;
         val=2;
         ierr=MatSetValues(AA,1,&indi,1,&indj,&val,INSERT_VALUES);
         indi=myrank*nlocal+i;
         indj=myrank*nlocal+i+1;
         val=-3;
         ierr=MatSetValues(AA,1,&indi,1,&indj,&val,INSERT_VALUES);
      }
    }
    ierr=MatAssemblyBegin(AA,MAT_FINAL_ASSEMBLY);
    ierr=MatAssemblyEnd(AA,MAT_FINAL_ASSEMBLY);

    /* Set vector b */
    ierr=VecSet(bb,1.e0);

    ierr=PetscFinalize();
    return 0;
}