#include <petscmat.h>
#include <petscvec.h>
#include <petscksp.h>

int main(int argc, char** argv) {
  Mat A00, A;
  Vec b1, b2, x1, x2, b, x;
  PetscErrorCode ierr;

  PetscInitialize(&argc, &argv, 0, 0);

  ierr = VecCreateSeq(PETSC_COMM_SELF, 10, &b1); CHKERRQ(ierr);
  ierr = VecDuplicate(b1, &b2); CHKERRQ(ierr);
  ierr = VecDuplicate(b1, &x1); CHKERRQ(ierr);
  ierr = VecDuplicate(b1, &x2); CHKERRQ(ierr);
  ierr = VecSet(b1, 1.0); CHKERRQ(ierr);
  ierr = VecSet(b2, 1.0); CHKERRQ(ierr);
  ierr = VecSet(x1, 0.0); CHKERRQ(ierr);
  ierr = VecSet(x2, 0.0); CHKERRQ(ierr);

  ierr = MatCreateSeqAIJ(PETSC_COMM_SELF, 10, 10, 1, PETSC_NULL, &A00); CHKERRQ(ierr);
  ierr = MatDiagonalSet(A00, b1, INSERT_VALUES); CHKERRQ(ierr);

  Mat mata[] = {A00, NULL, NULL, A00};
  ierr = MatCreateNest(PETSC_COMM_SELF, 2, NULL, 2, NULL, mata, &A); CHKERRQ(ierr);
  Vec vecx[] = {x1, x2};
  ierr = VecCreateNest(PETSC_COMM_SELF, 2, NULL, vecx, &x); CHKERRQ(ierr);
  Vec vecb[] = {b1, b2};
  ierr = VecCreateNest(PETSC_COMM_SELF, 2, NULL, vecb, &b); CHKERRQ(ierr);

  KSP solver;
  ierr = KSPCreate(PETSC_COMM_SELF, &solver); CHKERRQ(ierr);
  ierr = KSPSetOperators(solver, A, A); CHKERRQ(ierr);
  ierr = KSPSetFromOptions(solver); CHKERRQ(ierr);
  ierr = KSPSetUp(solver); CHKERRQ(ierr); // comment this line causes an error
  ierr = KSPSolve(solver, b, x); CHKERRQ(ierr);
  ierr = KSPDestroy(&solver); CHKERRQ(ierr);

  PetscFinalize();

  return 0;
}