#include <petscmat.h>
#include <petscvec.h>
#include <petscksp.h>
#include <vector>
#include <fstream>
#include <iostream>

PetscErrorCode read_vec(Vec& out, const char* file)
{
	PetscErrorCode ierr = 0;
    ierr = VecCreate(PETSC_COMM_WORLD, &out); CHKERRQ(ierr);
    ierr = VecSetFromOptions(out); CHKERRQ(ierr);

	PetscViewer reader;
    ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &reader); CHKERRQ(ierr);
    ierr = VecLoad(out, reader); CHKERRQ(ierr);
	ierr = PetscViewerDestroy(&reader); CHKERRQ(ierr);
	return ierr;
}

PetscErrorCode read_mat(Mat& out, const char* file)
{
	PetscErrorCode ierr = 0;
    ierr = MatCreate(PETSC_COMM_WORLD, &out); CHKERRQ(ierr);
    ierr = MatSetFromOptions(out); CHKERRQ(ierr);

	PetscViewer reader;
    ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &reader); CHKERRQ(ierr);
    ierr = MatLoad(out, reader); CHKERRQ(ierr);
	ierr = PetscViewerDestroy(&reader); CHKERRQ(ierr);
	return ierr;
}

int main(int argc, char* argv[])
{
	PetscErrorCode ierr = 0;
	ierr = PetscInitialize(&argc, &argv, PETSC_NULL, PETSC_NULL); CHKERRQ(ierr);

	Mat A_mat;
	ierr = read_mat(A_mat, "a.petsc_bin"); CHKERRQ(ierr);

	Vec b_vec, x0_vec;
	ierr = read_vec(b_vec, "b.petsc_bin"); CHKERRQ(ierr);
	ierr = read_vec(x0_vec, "x0.petsc_bin"); CHKERRQ(ierr);

	KSP ksp;
	ierr = KSPCreate(PETSC_COMM_WORLD, &ksp); CHKERRQ(ierr);
	PC pc;
	ierr = KSPGetPC(ksp, &pc); CHKERRQ(ierr);
	ierr = PCSetType(pc, PCNONE); CHKERRQ(ierr);

	PetscReal rtol, abstol, dtol;
	PetscInt  max_iter;
	ierr = KSPGetTolerances(ksp, &rtol, &abstol, &dtol, &max_iter); CHKERRQ(ierr);

	rtol = 1E-8;
	abstol = 0;
	max_iter = 10;

	ierr = KSPSetTolerances(ksp, rtol, abstol, dtol, max_iter); CHKERRQ(ierr);

	ierr = KSPSetOperators(ksp, A_mat, A_mat, SAME_PRECONDITIONER); CHKERRQ(ierr);

	ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); CHKERRQ(ierr);

	ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr);

	ierr = KSPSolve(ksp, b_vec, x0_vec); CHKERRQ(ierr);

	ierr = KSPDestroy(&ksp); CHKERRQ(ierr);
	ierr = MatDestroy(&A_mat); CHKERRQ(ierr);
	ierr = VecDestroy(&b_vec); CHKERRQ(ierr);
	ierr = VecDestroy(&x0_vec); CHKERRQ(ierr);

	ierr = PetscFinalize(); CHKERRQ(ierr);

	return ierr;
}