Actual source code: ex20.c
2: static char help[] = "Tests SDALocalToLocalxxx().\n\n";
4: #include petscda.h
6: /*
7: For testing purposes this example also creates a
8: DA context. Actually codes using SDA routines will probably
9: not also work with DA contexts.
10: */
14: int main(int argc,char **argv)
15: {
16: PetscMPIInt size,rank;
17: PetscInt M=8,dof=1,stencil_width=1,i,start,end,P=5,N = 6,m=PETSC_DECIDE,n=PETSC_DECIDE,p=PETSC_DECIDE,pt = 0,st = 0;
19: PetscTruth flg2,flg3,flg;
20: DAPeriodicType periodic = DA_NONPERIODIC;
21: DAStencilType stencil_type = DA_STENCIL_STAR;
22: DA da;
23: SDA sda;
24: Vec local,global,local_copy;
25: PetscScalar value,*in,*out;
26: PetscReal norm,work;
27: PetscViewer viewer;
28: char filename[PETSC_MAX_PATH_LEN];
29: FILE *file;
32: PetscInitialize(&argc,&argv,(char*)0,help);
34: PetscOptionsGetInt(PETSC_NULL,"-M",&M,PETSC_NULL);
35: PetscOptionsGetInt(PETSC_NULL,"-N",&N,PETSC_NULL);
36: PetscOptionsGetInt(PETSC_NULL,"-P",&P,PETSC_NULL);
37: PetscOptionsGetInt(PETSC_NULL,"-dof",&dof,PETSC_NULL);
38: PetscOptionsGetInt(PETSC_NULL,"-stencil_width",&stencil_width,PETSC_NULL);
39: PetscOptionsGetInt(PETSC_NULL,"-periodic",&pt,PETSC_NULL);
40: periodic = (DAPeriodicType) pt;
41: PetscOptionsGetInt(PETSC_NULL,"-stencil_type",&st,PETSC_NULL);
42: stencil_type = (DAStencilType) st;
44: PetscOptionsHasName(PETSC_NULL,"-1d",&flg2);
45: PetscOptionsHasName(PETSC_NULL,"-2d",&flg2);
46: PetscOptionsHasName(PETSC_NULL,"-3d",&flg3);
47: if (flg2) {
48: DACreate2d(PETSC_COMM_WORLD,periodic,stencil_type,M,N,m,n,dof,stencil_width,0,0,&da);
49: SDACreate2d(PETSC_COMM_WORLD,periodic,stencil_type,M,N,m,n,dof,stencil_width,0,0,&sda);
50: } else if (flg3) {
51: DACreate3d(PETSC_COMM_WORLD,periodic,stencil_type,M,N,P,m,n,p,dof,stencil_width,0,0,0,&da);
52: SDACreate3d(PETSC_COMM_WORLD,periodic,stencil_type,M,N,P,m,n,p,dof,stencil_width,0,0,0,&sda);
53: }
54: else {
55: DACreate1d(PETSC_COMM_WORLD,periodic,M,dof,stencil_width,PETSC_NULL,&da);
56: SDACreate1d(PETSC_COMM_WORLD,periodic,M,dof,stencil_width,PETSC_NULL,&sda);
57: }
59: DACreateGlobalVector(da,&global);
60: DACreateLocalVector(da,&local);
61: VecDuplicate(local,&local_copy);
62:
63: /* zero out vectors so that ghostpoints are zero */
64: value = 0;
65: VecSet(local,value);
66: VecSet(local_copy,value);
68: VecGetOwnershipRange(global,&start,&end);
69: for (i=start; i<end; i++) {
70: value = i + 1;
71: VecSetValues(global,1,&i,&value,INSERT_VALUES);
72: }
73: VecAssemblyBegin(global);
74: VecAssemblyEnd(global);
76: DAGlobalToLocalBegin(da,global,INSERT_VALUES,local);
77: DAGlobalToLocalEnd(da,global,INSERT_VALUES,local);
79: flg = PETSC_FALSE;
80: PetscOptionsGetTruth(PETSC_NULL,"-same_array",&flg,PETSC_NULL);
81: if (flg) {
82: /* test the case where the input and output array is the same */
83: VecCopy(local,local_copy);
84: VecGetArray(local_copy,&in);
85: VecRestoreArray(local_copy,PETSC_NULL);
86: SDALocalToLocalBegin(sda,in,INSERT_VALUES,in);
87: SDALocalToLocalEnd(sda,in,INSERT_VALUES,in);
88: } else {
89: VecGetArray(local,&out);
90: VecRestoreArray(local,PETSC_NULL);
91: VecGetArray(local_copy,&in);
92: VecRestoreArray(local_copy,PETSC_NULL);
93: SDALocalToLocalBegin(sda,out,INSERT_VALUES,in);
94: SDALocalToLocalEnd(sda,out,INSERT_VALUES,in);
95: }
97: flg = PETSC_FALSE;
98: PetscOptionsGetTruth(PETSC_NULL,"-save",&flg,PETSC_NULL);
99: if (flg) {
100: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
101: sprintf(filename,"local.%d",rank);
102: PetscViewerASCIIOpen(PETSC_COMM_SELF,filename,&viewer);
103: PetscViewerASCIIGetPointer(viewer,&file);
104: VecView(local,viewer);
105: fprintf(file,"Vector with correct ghost points\n");
106: VecView(local_copy,viewer);
107: PetscViewerDestroy(viewer);
108: }
110: VecAXPY(local_copy,-1.0,local);
111: VecNorm(local_copy,NORM_MAX,&work);
112: MPI_Allreduce(&work,&norm,1,MPIU_REAL,MPI_MAX,PETSC_COMM_WORLD);
113: if (norm != 0.0) {
114: MPI_Comm_size(PETSC_COMM_WORLD,&size);
115: PetscPrintf(PETSC_COMM_WORLD,"Norm of difference %G should be zero\n",norm);
116: PetscPrintf(PETSC_COMM_WORLD," Number of processors %d\n",size);
117: PetscPrintf(PETSC_COMM_WORLD," M,N,P,dof %D %D %D %D\n",M,N,P,dof);
118: PetscPrintf(PETSC_COMM_WORLD," stencil_width %D stencil_type %d periodic %d\n",stencil_width,(int)stencil_type,(int)periodic);
119: PetscPrintf(PETSC_COMM_WORLD," dimension %d\n",1 + (int) flg2 + (int) flg3);
120: }
121:
122: DADestroy(da);
123: SDADestroy(sda);
124: VecDestroy(local_copy);
125: VecDestroy(local);
126: VecDestroy(global);
127: PetscFinalize();
128: return 0;
129: }