static char help[] = "Tests VecLoad.\n\n";

#include "petsc.h"

int main(int argc,char **argv)
{
  PetscErrorCode      ierr;
  PetscInt            nx=3,ny=4,mode;
  PetscInt            ndofU=1,ndofV=1;
  PetscInt            i,j,c;
  PetscInt            xs,xm,ys,ym;
  DM                  daU,daV;
  Vec                 U,V;
  PetscReal           ***v,***u;
  PetscViewer         viewer;
  
  ierr = PetscInitialize(&argc,&argv,(char*)0,help);CHKERRQ(ierr);
  ierr = PetscOptionsGetInt(PETSC_NULL,"-nx",&nx,PETSC_NULL);CHKERRQ(ierr);
  ierr = PetscOptionsGetInt(PETSC_NULL,"-ny",&ny,PETSC_NULL);CHKERRQ(ierr);
  ierr = PetscOptionsGetInt(PETSC_NULL,"-ndofU",&ndofU,PETSC_NULL);CHKERRQ(ierr);
  ierr = PetscOptionsGetInt(PETSC_NULL,"-ndofV",&ndofV,PETSC_NULL);CHKERRQ(ierr);
  ierr = PetscOptionsGetInt(PETSC_NULL,"-mode",&mode,PETSC_NULL);CHKERRQ(ierr);

  
  if (mode == 0) {
    ierr = DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,
						DMDA_STENCIL_BOX,nx,ny,PETSC_DECIDE,PETSC_DECIDE,ndofU,1,
						PETSC_NULL,PETSC_NULL,&daU);CHKERRQ(ierr);    
    ierr = DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,
						DMDA_STENCIL_BOX,nx,ny,PETSC_DECIDE,PETSC_DECIDE,ndofV,1,
						PETSC_NULL,PETSC_NULL,&daV);CHKERRQ(ierr);    
		ierr = DMDAGetCorners(daU,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);CHKERRQ(ierr);

    ierr = DMGetGlobalVector(daU,&U);CHKERRQ(ierr);
    ierr = PetscObjectSetName((PetscObject) U,"UVec");CHKERRQ(ierr);
    ierr = DMGetGlobalVector(daV,&V);CHKERRQ(ierr);
    ierr = PetscObjectSetName((PetscObject) V,"VVec");CHKERRQ(ierr);

    ierr = DMDAVecGetArrayDOF(daU,U,&u);CHKERRQ(ierr);
    ierr = DMDAVecGetArrayDOF(daV,V,&v);CHKERRQ(ierr);
    for (j = ys; j < ys + ym; j++) {
      for (i = xs; i < xs + xm; i++) {
        for (c = 0; c < ndofU; c++) {
          u[j][i][c] = c + 10*i + j;
        }
        for (c = 0; c < ndofV; c++) {
          v[j][i][c] = -100*c -10*i - j;
        }
      }
    }
    ierr = DMDAVecRestoreArrayDOF(daV,V,&v);CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayDOF(daU,U,&u);CHKERRQ(ierr);
    ierr = VecView(U,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
    ierr = VecView(V,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, "Writing vectors to file\n");CHKERRQ(ierr);
    ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"vector.dat",FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);
    ierr = DMView(daU,viewer);CHKERRQ(ierr);
    ierr = DMView(daV,viewer);CHKERRQ(ierr);
    ierr = VecView(U, viewer);CHKERRQ(ierr);
    ierr = VecView(V, viewer);CHKERRQ(ierr);

    ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
  } else {
    ierr = PetscPrintf(PETSC_COMM_WORLD, "Reading vectors from file\n");CHKERRQ(ierr);
    ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"vector.dat",FILE_MODE_READ,&viewer);CHKERRQ(ierr);

    ierr = DMCreate(PETSC_COMM_WORLD,&daU);CHKERRQ(ierr);
    ierr = DMLoad(daU,viewer);CHKERRQ(ierr);
    ierr = DMDAGetInfo(daU,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
                       &ndofU,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
    ierr = DMView(daU,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
    
    ierr = DMCreate(PETSC_COMM_WORLD,&daV);CHKERRQ(ierr);
    ierr = DMLoad(daV,viewer);CHKERRQ(ierr);
    ierr = DMDAGetInfo(daV,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
                       &ndofV,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
    ierr = DMView(daV,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
    
    ierr = VecCreate(PETSC_COMM_WORLD,&U);CHKERRQ(ierr);
    ierr = VecSetBlockSize(U,ndofU);CHKERRQ(ierr);
    ierr = PetscObjectSetName((PetscObject) U,"UVec");CHKERRQ(ierr);
    ierr = VecCreate(PETSC_COMM_WORLD,&V);CHKERRQ(ierr);
    ierr = VecSetBlockSize(V,ndofV);CHKERRQ(ierr);
    ierr = PetscObjectSetName((PetscObject) V,"VVec");CHKERRQ(ierr);
    
    ierr = VecLoad(U,viewer);CHKERRQ(ierr);
    ierr = VecLoad(V,viewer);CHKERRQ(ierr);
    ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
    ierr = VecView(U,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
    ierr = VecView(V,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  }
  ierr = DMRestoreGlobalVector(daV,&V);CHKERRQ(ierr);
  ierr = DMRestoreGlobalVector(daU,&U);CHKERRQ(ierr);
  ierr = DMDestroy(&daU);CHKERRQ(ierr);
  ierr = DMDestroy(&daV);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return(0);
}