#include <petsc.h>
static char help[] = "";

int main(int argc,char **args)
{
  PetscInt  n = 5;
  Mat       array[4],N;
  Vec       v[2],view,x,b;
  KSP       ksp;
  PetscBool flg = PETSC_FALSE;

  PetscCall(PetscInitialize(&argc,&args,NULL,help));
  PetscCall(MatCreateConstantDiagonal(PETSC_COMM_WORLD,n,n,n,n,-1.0,array));
  PetscCall(MatConvert(array[0],MATDENSE,MAT_INPLACE_MATRIX,array));
  PetscCall(MatSetRandom(array[0],NULL));
  PetscCall(MatCreateConstantDiagonal(PETSC_COMM_WORLD,n,n,n,n,-1.0,array+3));
  array[1] = array[2] = NULL;
  PetscCall(MatCreateNest(PETSC_COMM_WORLD,2,NULL,2,NULL,array,&N));
  PetscCall(MatNestSetVecType(N,VECNEST));
  PetscCall(KSPCreate(PETSC_COMM_WORLD,&ksp));
  PetscCall(KSPSetFromOptions(ksp));
  PetscCall(KSPSetOperators(ksp,N,N));
  PetscCall(MatCreateVecs(N,v+1,v));
  PetscCall(MatCreateVecs(array[0],&x,&b));
  PetscCall(PetscOptionsGetBool(NULL,NULL,"-nest_subvec",&flg,NULL));
  if (flg) {
    PetscCall(VecNestGetSubVec(v[0],1,&view));
    PetscCall(VecSet(view,0.0));
  } else PetscCall(VecSet(v[0],0.0));
  PetscCall(VecNestGetSubVec(v[0],0,&view));
  PetscCall(VecSetRandom(b,NULL));
  PetscCall(VecCopy(b,view));
  PetscCall(KSPSolve(ksp,v[0],v[1]));
  PetscCall(VecDestroy(v+1));
  PetscCall(VecDestroy(v));
  PetscCall(KSPDestroy(&ksp));
  PetscCall(VecDestroy(&x));
  PetscCall(VecDestroy(&b));
  PetscCall(MatDestroy(array));
  PetscCall(MatDestroy(array+3));
  PetscCall(MatDestroy(&N));
  PetscCall(PetscFinalize());
  return 0;
}