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 = 20, Nx = 10; // domain size
  int bs = 1;		// depth of the boundary
  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,1.); 
  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; 
  if (ys ==        0) ys+=bs; 
  if (ye == infox.my) ye-=bs; 
  if (xs ==        0) xs+=bs; 
  if (xe == infox.mx) xe-=bs; 

  for (int j = ys; j < ye; j++) {
    for (int i = xs; i < xe; i++) {
      for (int k = 0; k < infof.dof; k++) {
	int row = k+ i    * infox.dof +  j     * infox.dof * infox.mx; 
	int col = k+(i-bs)* infof.dof + (j-bs) * infof.dof * infof.mx; 
	const PetscScalar one = 1.; 
	MatSetValue(B,row,col,one,INSERT_VALUES); 
      }
    }
  }
  MatAssemblyBegin ( B,MAT_FINAL_ASSEMBLY ); 
  MatAssemblyEnd   ( B,MAT_FINAL_ASSEMBLY ); 

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

}