static char help[]="Test for x = B f";

#include <petscdmda.h>

int main(int argc,char **argv) {
  PetscErrorCode ierr =0;
  PetscInitialize(&argc,&argv,(char*)0,help);
  int Mx = 5, Nx = 4; // domain size
  int bs = 1;		// depth of the boundary
  int ii,jj;
  DM dax,daf;		// domain and interior point da
  Vec x,f;		// vectors for each domain
  Mat B;		// matrix x = B f
  DMDACreate2d(PETSC_COMM_WORLD , DM_BOUNDARY_GHOSTED , DM_BOUNDARY_GHOSTED , DMDA_STENCIL_BOX ,
      Mx      , Nx      , PETSC_DECIDE , PETSC_DECIDE , 1 , 0 , 0 , 0 , &dax); 
  DMDACreate2d(PETSC_COMM_WORLD , DM_BOUNDARY_NONE    , DM_BOUNDARY_NONE    , DMDA_STENCIL_BOX ,
      Mx-2*bs , Nx-2*bs , PETSC_DECIDE , PETSC_DECIDE , 1 , 0 , 0 , 0 , &daf); 

  // Creating the vectors
  DMCreateGlobalVector(dax,&x); 
  DMCreateGlobalVector(daf,&f); 
  VecSet(f,0.); 
  VecSet(x,0.); 

  // Creating the matrix
  DMDALocalInfo infox,infof; 
  DMDAGetLocalInfo(dax,&infox); 
  DMDAGetLocalInfo(daf,&infof); 
  PetscInt nrows = infox.xm * infox.ym * infox.dof,
	   ncols = infof.xm * infof.ym * infof.dof; 
  MatCreate(PETSC_COMM_WORLD,&B); 
  MatSetSizes(B,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE); 
  MatSetFromOptions(B); 
  MatSeqAIJSetPreallocation(B,1,NULL); 
  MatMPIAIJSetPreallocation(B,1,NULL,1,NULL); 

  ISLocalToGlobalMapping xmap,fmap; 
  DMGetLocalToGlobalMapping(dax,&xmap); 
  DMGetLocalToGlobalMapping(daf,&fmap); 
  MatSetLocalToGlobalMapping(B,xmap,fmap); 

  // Filling the matrix
  PetscInt xs = infox.xs,
	   xe = infox.xs+infox.xm,
	   ys = infox.ys,
	   ye = infox.ys+infox.ym; 

  // trim off the points in the x vector that are not in the f vector
  if (ys ==        0) ys+=bs; 
  if (ye == infox.my) ye-=bs; 
  if (xs ==        0) xs+=bs; 
  if (xe == infox.mx) xe-=bs;
  
  // loop over rows of B matrix 
  for (int j = ys; j < ye; j++) {
    jj = j - infox.gys;
    for (int i = xs; i < xe; i++) {
      ii = i - infox.gxs;
      for (int k = 0; k < infof.dof; k++) {
	int row = k+ ii    * infox.dof +  jj     * infox.dof * infox.gxm; 
        int col = k+(i-bs-infof.gxs)* infof.dof + (j-bs-infof.gys) * infof.dof * infof.gxm;
        int rank;
        MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
	PetscScalar v = rank+1;
	MatSetValuesLocal(B,1,&row,1,&col,&v,INSERT_VALUES);
        VecSetValuesLocal(x,1,&row,&v,ADD_VALUES);
        VecSetValuesLocal(f,1,&col,&v,ADD_VALUES);
      }
    }
  }
  VecAssemblyBegin(x);
  VecAssemblyEnd(x);
  
  MatAssemblyBegin ( B,MAT_FINAL_ASSEMBLY ); 
  MatAssemblyEnd   ( B,MAT_FINAL_ASSEMBLY ); 

  // Multiplying and viewing
  //MatMult(B,f,x); 
  VecView(x,PETSC_VIEWER_STDOUT_WORLD);
  VecView(f,PETSC_VIEWER_STDOUT_WORLD);
  MatView(B,PETSC_VIEWER_STDOUT_WORLD);
  PetscFinalize(); 

}