Actual source code: ex20.c

  1: /*$Id: ex1.c,v 1.19 2001/08/10 03:34:58 bsmith Exp $*/

  3: static char help[] = "Tests SDALocalToLocalxxx().\n\n";

 5:  #include petscda.h
 6:  #include petscsys.h

  8: /*
  9:          For testing purposes this example also creates a 
 10:    DA context. Actually codes using SDA routines will probably 
 11:    not also work with DA contexts.
 12: */

 16: int main(int argc,char **argv)
 17: {
 18:   int            size,rank,M=8,ierr,dof=1,stencil_width=1,i,start,end,P=5;
 19:   int            N = 6,m=PETSC_DECIDE,n=PETSC_DECIDE,p=PETSC_DECIDE;
 20:   PetscTruth     flg2,flg3,flg;
 21:   DAPeriodicType periodic = DA_NONPERIODIC;
 22:   DAStencilType  stencil_type = DA_STENCIL_STAR;
 23:   DA             da;
 24:   SDA            sda;
 25:   Vec            local,global,local_copy;
 26:   PetscScalar    value,mone = -1.0,*in,*out;
 27:   PetscReal      norm,work;
 28:   PetscViewer    viewer;
 29:   char           filename[64];
 30:   FILE           *file;


 33:   PetscInitialize(&argc,&argv,(char*)0,help);

 35:   PetscOptionsGetInt(PETSC_NULL,"-M",&M,PETSC_NULL);
 36:   PetscOptionsGetInt(PETSC_NULL,"-N",&N,PETSC_NULL);
 37:   PetscOptionsGetInt(PETSC_NULL,"-P",&P,PETSC_NULL);
 38:   PetscOptionsGetInt(PETSC_NULL,"-dof",&dof,PETSC_NULL);
 39:   PetscOptionsGetInt(PETSC_NULL,"-stencil_width",&stencil_width,PETSC_NULL);
 40:   PetscOptionsGetInt(PETSC_NULL,"-periodic",(int*)&periodic,PETSC_NULL);
 41:   PetscOptionsGetInt(PETSC_NULL,"-stencil_type",(int*)&stencil_type,PETSC_NULL);

 43:   PetscOptionsHasName(PETSC_NULL,"-1d",&flg2);
 44:   PetscOptionsHasName(PETSC_NULL,"-2d",&flg2);
 45:   PetscOptionsHasName(PETSC_NULL,"-3d",&flg3);
 46:   if (flg2) {
 47:     DACreate2d(PETSC_COMM_WORLD,periodic,stencil_type,M,N,m,n,dof,stencil_width,0,0,&da);
 48: 
 49:     SDACreate2d(PETSC_COMM_WORLD,periodic,stencil_type,M,N,m,n,dof,stencil_width,0,0,&sda);
 50: 
 51:   } else if (flg3) {
 52:     DACreate3d(PETSC_COMM_WORLD,periodic,stencil_type,M,N,P,m,n,p,dof,stencil_width,
 53:                       0,0,0,&da);
 54:     SDACreate3d(PETSC_COMM_WORLD,periodic,stencil_type,M,N,P,m,n,p,dof,stencil_width,
 55:                       0,0,0,&sda);
 56:   }
 57:   else {
 58:     DACreate1d(PETSC_COMM_WORLD,periodic,M,dof,stencil_width,PETSC_NULL,&da);
 59:     SDACreate1d(PETSC_COMM_WORLD,periodic,M,dof,stencil_width,PETSC_NULL,&sda);
 60:   }

 62:   DACreateGlobalVector(da,&global);
 63:   DACreateLocalVector(da,&local);
 64:   VecDuplicate(local,&local_copy);

 66: 
 67:   /* zero out vectors so that ghostpoints are zero */
 68:   value = 0;
 69:   VecSet(&value,local);
 70:   VecSet(&value,local_copy);

 72:   VecGetOwnershipRange(global,&start,&end);
 73:   for (i=start; i<end; i++) {
 74:     value = i + 1;
 75:     VecSetValues(global,1,&i,&value,INSERT_VALUES);
 76:   }
 77:   VecAssemblyBegin(global);
 78:   VecAssemblyEnd(global);

 80:   DAGlobalToLocalBegin(da,global,INSERT_VALUES,local);
 81:   DAGlobalToLocalEnd(da,global,INSERT_VALUES,local);


 84:   PetscOptionsHasName(PETSC_NULL,"-same_array",&flg);
 85:   if (flg) {
 86:     /* test the case where the input and output array is the same */
 87:     VecCopy(local,local_copy);
 88:     VecGetArray(local_copy,&in);
 89:     VecRestoreArray(local_copy,PETSC_NULL);
 90:     SDALocalToLocalBegin(sda,in,INSERT_VALUES,in);
 91:     SDALocalToLocalEnd(sda,in,INSERT_VALUES,in);
 92:   } else {
 93:     VecGetArray(local,&out);
 94:     VecRestoreArray(local,PETSC_NULL);
 95:     VecGetArray(local_copy,&in);
 96:     VecRestoreArray(local_copy,PETSC_NULL);
 97:     SDALocalToLocalBegin(sda,out,INSERT_VALUES,in);
 98:     SDALocalToLocalEnd(sda,out,INSERT_VALUES,in);
 99:   }

101:   PetscOptionsHasName(PETSC_NULL,"-save",&flg);
102:   if (flg) {
103:     MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
104:     sprintf(filename,"local.%d",rank);
105:     PetscViewerASCIIOpen(PETSC_COMM_SELF,filename,&viewer);
106:     PetscViewerASCIIGetPointer(viewer,&file);
107:     VecView(local,viewer);
108:     fprintf(file,"Vector with correct ghost points\n");
109:     VecView(local_copy,viewer);
110:     PetscViewerDestroy(viewer);
111:   }

113:   VecAXPY(&mone,local,local_copy);
114:   VecNorm(local_copy,NORM_MAX,&work);
115:   MPI_Allreduce(&work,&norm,1,MPIU_REAL,MPI_MAX,PETSC_COMM_WORLD);
116:   if (norm != 0.0) {
117:     MPI_Comm_size(PETSC_COMM_WORLD,&size);
118:     PetscPrintf(PETSC_COMM_WORLD,"Norm of difference %g should be zero\n",norm);
119:     PetscPrintf(PETSC_COMM_WORLD,"  Number of processors %d\n",size);
120:     PetscPrintf(PETSC_COMM_WORLD,"  M,N,P,dof %d %d %d %d\n",M,N,P,dof);
121:     PetscPrintf(PETSC_COMM_WORLD,"  stencil_width %d stencil_type %d periodic %d\n",stencil_width,stencil_type,periodic);
122:     PetscPrintf(PETSC_COMM_WORLD,"  dimension %d\n",1 + (int) flg2 + (int) flg3);
123:   }
124: 
125:   DADestroy(da);
126:   SDADestroy(sda);
127:   VecDestroy(local_copy);
128:   VecDestroy(local);
129:   VecDestroy(global);
130:   PetscFinalize();
131:   return 0;
132: }