#include "petscksp.h"
#include <mpi.h>
#include<iostream>
#include <fstream>
#include <stdio.h>
#include <stdlib.h>
#include <iomanip>
#include <vector>
#include <cmath>
#include <string>
#include <sstream>
#include <limits>


static char help[] = "Solves a tridiagonal linear system with KSP.\n\n";

using namespace std;



int Rank,np;
int PetscPrintError(const char error[],...){
        if (Rank==0) cerr << "PETSc Error ... exiting" << endl;
        exit(1);
        return 0;
}




int main(int argc, char *argv[]) {

        // Initialize mpi
        MPI_Init(&argc,&argv);
        MPI_Comm_rank(MPI_COMM_WORLD, &Rank);
        MPI_Comm_size(MPI_COMM_WORLD, &np);



	KSP ksp;
	Mat A;
	Vec x, b;
	int n, its;
	PetscErrorCode ierr;
        PetscScalar value;
	PetscMPIInt np;
	int i,j; 

        PetscInitialize(&argc,&argv,(char *)0,help);
        PetscErrorPrintf = PetscPrintError;

	ierr = MPI_Comm_size(PETSC_COMM_WORLD,&np);CHKERRQ(ierr);
	if (np != 1) SETERRQ(PETSC_COMM_SELF,1,"This is a uniprocessor example only!");
	ierr = PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);CHKERRQ(ierr);


	/* linear solver context */
	/* matrix */
	/* solution, RHS vectors */
	/* problem dimension, number of iterations */
	n=2;

	cout<<"n="<<n<<endl;


	ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
	ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);
	ierr = MatSetFromOptions(A);CHKERRQ(ierr);


	value=1.; i=0; j=0;
	ierr = MatSetValues(A,1,&i,1,&j,&value,INSERT_VALUES);CHKERRQ(ierr);
	value=1.; i=0; j=1;
	ierr = MatSetValues(A,1,&i,1,&j,&value,INSERT_VALUES);CHKERRQ(ierr);
	value=1.; i=1; j=0;
	ierr = MatSetValues(A,1,&i,1,&j,&value,INSERT_VALUES);CHKERRQ(ierr);
	value=-1.; i=1; j=1;
	ierr = MatSetValues(A,1,&i,1,&j,&value,INSERT_VALUES);CHKERRQ(ierr);

	ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
	ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

	MatView(A,PETSC_VIEWER_STDOUT_WORLD);     

	exit(0);

	/* (code to assemble matrix not shown) */
	VecCreate(PETSC_COMM_WORLD,&x);
	VecSetSizes(x,PETSC_DECIDE, n);
	VecSetFromOptions(x);
	VecDuplicate(x,&b);
	/* (code to assemble RHS vector not shown)*/
	KSPCreate(PETSC_COMM_WORLD, &ksp);
	//KSPSetOperators(ksp, A, A, DIFFERENT_NONZERO_PATTERN);
	KSPSetOperators(ksp, A, A);
	KSPSetFromOptions(ksp);


	KSPSolve(ksp, b, x);
	KSPDestroy(&ksp);

}