static char help[] = "Load data into subvector.\n\n";

#include <petscvec.h>
#include <petscviewerhdf5.h>

int main(int argc,char **argv)
{
  Vec            testvec, loadvec, writevec;
  IS             loadis;
  PetscMPIInt    size,rank;
  PetscInt       low,high;
  PetscInt       issize, isstart;
  PetscInt       skipuntil = 50;
  PetscViewer    viewer; 
  PetscScalar    sum;
  PetscBool      usehdf5 = PETSC_TRUE;
  PetscBool      usebinary = PETSC_FALSE;
  PetscErrorCode ierr;

  PetscFunctionBegin;

  ierr = PetscInitialize(&argc, &argv, (char*) 0, help);if (ierr) return ierr;
  ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);CHKERRQ(ierr);
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRQ(ierr);

  ierr = PetscOptionsGetBool(NULL,NULL,"-binary",&usebinary,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsGetBool(NULL,NULL,"-hdf5",&usehdf5,NULL);CHKERRQ(ierr);

  ierr = VecCreate(PETSC_COMM_WORLD, &testvec);CHKERRQ(ierr); 
  ierr = VecSetSizes(testvec, PETSC_DECIDE,100);CHKERRQ(ierr);
  ierr = VecSetUp(testvec);CHKERRQ(ierr);
  ierr = VecSet(testvec, (PetscScalar) 1);CHKERRQ(ierr); 

  ierr = VecCreate(PETSC_COMM_WORLD, &writevec);CHKERRQ(ierr); 
  ierr = VecSetSizes(writevec, PETSC_DECIDE,50);CHKERRQ(ierr);
  ierr = VecSetUp(writevec);CHKERRQ(ierr);
  ierr = VecSet(writevec, (PetscScalar) 2);CHKERRQ(ierr); 
  ierr = PetscObjectSetName((PetscObject)writevec,"temp"); CHKERRQ(ierr);

  if (usebinary) {
    ierr = PetscPrintf(PETSC_COMM_WORLD,"writing vector in binary to vector.dat ...\n");CHKERRQ(ierr);
    ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"vector.dat",FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);
  } else if (usehdf5) {
    ierr = PetscPrintf(PETSC_COMM_WORLD,"writing vector in hdf5 to vector.dat ...\n");CHKERRQ(ierr);
    ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD,"vector.dat",FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);
  }
  ierr = VecView(writevec,viewer);CHKERRQ(ierr);
  ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);

  ierr = VecGetOwnershipRange(testvec, &low, &high);CHKERRQ(ierr);
  if( low>=skipuntil ) {
    isstart = low;
    issize = high - low;
  } else if( low<=skipuntil && high>=skipuntil ) {
    isstart = skipuntil;
    issize = high - skipuntil;
  } else {
    isstart = low;
    issize  = 0;
  }
  ierr = ISCreateStride(PETSC_COMM_WORLD, issize, isstart, 1, &loadis);CHKERRQ(ierr);

  ierr = VecGetSubVector(testvec, loadis, &loadvec);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject)loadvec,"temp"); CHKERRQ(ierr);

  if (usebinary) {
    ierr = PetscPrintf(PETSC_COMM_WORLD,"reading vector in binary from vector.dat ...\n");CHKERRQ(ierr);
    ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"vector.dat",FILE_MODE_READ,&viewer);CHKERRQ(ierr);
  } else if (usehdf5) {
    ierr = PetscPrintf(PETSC_COMM_WORLD,"reading vector in hdf5 from vector.dat ...\n");CHKERRQ(ierr);
    ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD,"vector.dat",FILE_MODE_READ,&viewer);CHKERRQ(ierr);
  }
  ierr = VecLoad(loadvec, viewer);CHKERRQ(ierr); 
  ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);

  ierr = VecRestoreSubVector(testvec, loadis, &loadvec);CHKERRQ(ierr);
  ierr = VecDestroy(&loadvec); CHKERRQ(ierr);  
 
  ierr = VecSum(testvec, &sum); CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,"VecSum on testvec is : %e\n",sum);CHKERRQ(ierr);

  ierr = VecDestroy(&testvec); CHKERRQ(ierr);   
  ierr = PetscFinalize();
  return ierr;
}