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: }