#include "petscda.h"
#include "petsc.h"

const PetscInt BLOCK = 50;
const PetscInt NVECS = 100;

int main(int argc, char **args)
{
  PetscInitialize(&argc, &args, PETSC_NULL, PETSC_NULL);

  PetscInt m, n, p, is, js, ks, i, j, k, nv;
  DA das[NVECS];

  // Create the distributed arrays:
  for(nv = 0; nv < NVECS; nv++) {
    
    DACreate3d(PETSC_COMM_WORLD, DA_NONPERIODIC, DA_STENCIL_BOX,
	       BLOCK, BLOCK, BLOCK,
	       PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE,
	       1, 1,
	       PETSC_NULL, PETSC_NULL, PETSC_NULL,
	       &das[nv]);
  }

  DAGetCorners(das[0], &ks, &js, &is, &p, &n, &m);

  // Create a bunch of vectors:
  Vec globals[NVECS], locals[NVECS];
  PetscReal ***ldata, ***gdata;

  for(i = 0; i < NVECS; i++) {
    DACreateGlobalVector(das[i], &globals[i]);
    DACreateLocalVector(das[i], &locals[i]);
  }

  // Populate the vectors:
  for(nv = 0; nv < NVECS; nv++) {
    
    DAVecGetArray(das[nv], globals[nv], &gdata);

    for(i = is; i < is + m; i++)
      for(j = js; j < js + n; j++)
	for(k = ks; k < ks + p; k++) {
	  
	  // Do some math:
	  gdata[i][j][k] = nv*nv*nv + 100.0*nv;
      }

    DAVecRestoreArray(das[nv], globals[nv], &gdata);
  }

  // Scatter the vectors:
  for(nv = 0; nv < NVECS; nv++) {
    DAGlobalToLocalBegin(das[nv], globals[nv], INSERT_VALUES, locals[nv]);
    DAGlobalToLocalEnd(das[nv], globals[nv], INSERT_VALUES, locals[nv]);
  }

  // CLEAN UP:
  for(i = 0; i < NVECS; i++) {
    VecDestroy(globals[i]);
    VecDestroy(locals[i]);
    DADestroy(das[i]);
  }

  PetscFinalize();
  return 0;
}