Actual source code: ex15.c
2: static char help[] = "Tests VecView() functionality with DMDA objects when using:"\
3: "(i) a PetscViewer binary with MPI-IO support; and (ii) when the binary header is skipped.\n\n";
5: #include <petscdm.h>
6: #include <petscdmda.h>
8: #define DMDA_I 5
9: #define DMDA_J 4
10: #define DMDA_K 6
12: const PetscReal dmda_i_val[] = { 1.10, 2.3006, 2.32444, 3.44006, 66.9009 };
13: const PetscReal dmda_j_val[] = { 0.0, 0.25, 0.5, 0.75 };
14: const PetscReal dmda_k_val[] = { 0.0, 1.1, 2.2, 3.3, 4.4, 5.5 };
16: PetscErrorCode MyVecDump(const char fname[],PetscBool skippheader,PetscBool usempiio,Vec x)
17: {
18: MPI_Comm comm;
19: PetscViewer viewer;
20: PetscBool ismpiio,isskip;
23: PetscObjectGetComm((PetscObject)x,&comm);
25: PetscViewerCreate(comm,&viewer);
26: PetscViewerSetType(viewer,PETSCVIEWERBINARY);
27: if (skippheader) PetscViewerBinarySetSkipHeader(viewer,PETSC_TRUE);
28: PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
29: if (usempiio) PetscViewerBinarySetUseMPIIO(viewer,PETSC_TRUE);
30: PetscViewerFileSetName(viewer,fname);
32: VecView(x,viewer);
34: PetscViewerBinaryGetUseMPIIO(viewer,&ismpiio);
35: if (ismpiio) PetscPrintf(comm,"*** PetscViewer[write] using MPI-IO ***\n");
36: PetscViewerBinaryGetSkipHeader(viewer,&isskip);
37: if (isskip) PetscPrintf(comm,"*** PetscViewer[write] skipping header ***\n");
39: PetscViewerDestroy(&viewer);
40: return 0;
41: }
43: PetscErrorCode MyVecLoad(const char fname[],PetscBool skippheader,PetscBool usempiio,Vec x)
44: {
45: MPI_Comm comm;
46: PetscViewer viewer;
47: PetscBool ismpiio,isskip;
50: PetscObjectGetComm((PetscObject)x,&comm);
52: PetscViewerCreate(comm,&viewer);
53: PetscViewerSetType(viewer,PETSCVIEWERBINARY);
54: if (skippheader) PetscViewerBinarySetSkipHeader(viewer,PETSC_TRUE);
55: PetscViewerFileSetMode(viewer,FILE_MODE_READ);
56: if (usempiio) PetscViewerBinarySetUseMPIIO(viewer,PETSC_TRUE);
57: PetscViewerFileSetName(viewer,fname);
59: VecLoad(x,viewer);
61: PetscViewerBinaryGetSkipHeader(viewer,&isskip);
62: if (isskip) PetscPrintf(comm,"*** PetscViewer[load] skipping header ***\n");
63: PetscViewerBinaryGetUseMPIIO(viewer,&ismpiio);
64: if (ismpiio) PetscPrintf(comm,"*** PetscViewer[load] using MPI-IO ***\n");
66: PetscViewerDestroy(&viewer);
67: return 0;
68: }
70: PetscErrorCode DMDAVecGenerateEntries(DM dm,Vec a)
71: {
72: PetscScalar ****LA_v;
73: PetscInt i,j,k,l,si,sj,sk,ni,nj,nk,M,N,dof;
76: DMDAGetInfo(dm,NULL,&M,&N,NULL,NULL,NULL,NULL,&dof,NULL,NULL,NULL,NULL,NULL);
77: DMDAGetCorners(dm,&si,&sj,&sk,&ni,&nj,&nk);
78: DMDAVecGetArrayDOF(dm,a,&LA_v);
79: for (k=sk; k<sk+nk; k++) {
80: for (j=sj; j<sj+nj; j++) {
81: for (i=si; i<si+ni; i++) {
82: PetscScalar test_value_s;
84: test_value_s = dmda_i_val[i]*((PetscScalar)i) + dmda_j_val[j]*((PetscScalar)(i+j*M)) + dmda_k_val[k]*((PetscScalar)(i + j*M + k*M*N));
85: for (l=0; l<dof; l++) {
86: LA_v[k][j][i][l] = (PetscScalar)dof * test_value_s + (PetscScalar)l;
87: }
88: }
89: }
90: }
91: DMDAVecRestoreArrayDOF(dm,a,&LA_v);
92: return 0;
93: }
95: PetscErrorCode HeaderlessBinaryReadCheck(DM dm,const char name[])
96: {
97: int fdes;
98: PetscScalar buffer[DMDA_I*DMDA_J*DMDA_K*10];
99: PetscInt len,d,i,j,k,M,N,dof;
100: PetscMPIInt rank;
101: PetscBool dataverified = PETSC_TRUE;
104: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
105: DMDAGetInfo(dm,NULL,&M,&N,NULL,NULL,NULL,NULL,&dof,NULL,NULL,NULL,NULL,NULL);
106: len = DMDA_I*DMDA_J*DMDA_K*dof;
107: if (rank == 0) {
108: PetscBinaryOpen(name,FILE_MODE_READ,&fdes);
109: PetscBinaryRead(fdes,buffer,len,NULL,PETSC_SCALAR);
110: PetscBinaryClose(fdes);
112: for (k=0; k<DMDA_K; k++) {
113: for (j=0; j<DMDA_J; j++) {
114: for (i=0; i<DMDA_I; i++) {
115: for (d=0; d<dof; d++) {
116: PetscScalar v,test_value_s,test_value;
117: PetscInt index;
119: test_value_s = dmda_i_val[i]*((PetscScalar)i) + dmda_j_val[j]*((PetscScalar)(i+j*M)) + dmda_k_val[k]*((PetscScalar)(i + j*M + k*M*N));
120: test_value = (PetscScalar)dof * test_value_s + (PetscScalar)d;
122: index = dof*(i + j*M + k*M*N) + d;
123: v = PetscAbsScalar(test_value-buffer[index]);
124: #if defined(PETSC_USE_COMPLEX)
125: if ((PetscRealPart(v) > 1.0e-10) || (PetscImaginaryPart(v) > 1.0e-10)) {
126: PetscPrintf(PETSC_COMM_SELF,"ERROR: Difference > 1.0e-10 occurred (delta = (%+1.12e,%+1.12e) [loc %D,%D,%D(%D)])\n",(double)PetscRealPart(test_value),(double)PetscImaginaryPart(test_value),i,j,k,d);
127: dataverified = PETSC_FALSE;
128: }
129: #else
130: if (PetscRealPart(v) > 1.0e-10) {
131: PetscPrintf(PETSC_COMM_SELF,"ERROR: Difference > 1.0e-10 occurred (delta = %+1.12e [loc %D,%D,%D(%D)])\n",(double)PetscRealPart(test_value),i,j,k,d);
132: dataverified = PETSC_FALSE;
133: }
134: #endif
135: }
136: }
137: }
138: }
139: if (dataverified) {
140: PetscPrintf(PETSC_COMM_SELF,"Headerless read of data verified for: %s\n",name);
141: }
142: }
143: return 0;
144: }
146: PetscErrorCode VecCompare(Vec a,Vec b)
147: {
148: PetscInt locmin[2],locmax[2];
149: PetscReal min[2],max[2];
150: Vec ref;
153: VecMin(a,&locmin[0],&min[0]);
154: VecMax(a,&locmax[0],&max[0]);
156: VecMin(b,&locmin[1],&min[1]);
157: VecMax(b,&locmax[1],&max[1]);
159: PetscPrintf(PETSC_COMM_WORLD,"VecCompare\n");
160: PetscPrintf(PETSC_COMM_WORLD," min(a) = %+1.2e [loc %D]\n",(double)min[0],locmin[0]);
161: PetscPrintf(PETSC_COMM_WORLD," max(a) = %+1.2e [loc %D]\n",(double)max[0],locmax[0]);
163: PetscPrintf(PETSC_COMM_WORLD," min(b) = %+1.2e [loc %D]\n",(double)min[1],locmin[1]);
164: PetscPrintf(PETSC_COMM_WORLD," max(b) = %+1.2e [loc %D]\n",(double)max[1],locmax[1]);
166: VecDuplicate(a,&ref);
167: VecCopy(a,ref);
168: VecAXPY(ref,-1.0,b);
169: VecMin(ref,&locmin[0],&min[0]);
170: if (PetscAbsReal(min[0]) > 1.0e-10) {
171: PetscPrintf(PETSC_COMM_WORLD," ERROR: min(a-b) > 1.0e-10\n");
172: PetscPrintf(PETSC_COMM_WORLD," min(a-b) = %+1.10e\n",(double)PetscAbsReal(min[0]));
173: } else {
174: PetscPrintf(PETSC_COMM_WORLD," min(a-b) < 1.0e-10\n");
175: }
176: VecDestroy(&ref);
177: return 0;
178: }
180: PetscErrorCode TestDMDAVec(PetscBool usempiio)
181: {
182: DM dm;
183: Vec x_ref,x_test;
184: PetscBool skipheader = PETSC_TRUE;
188: if (!usempiio) PetscPrintf(PETSC_COMM_WORLD,"%s\n",PETSC_FUNCTION_NAME);
189: else PetscPrintf(PETSC_COMM_WORLD,"%s [using mpi-io]\n",PETSC_FUNCTION_NAME);
190: DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,DMDA_I,DMDA_J,DMDA_K,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,
191: 3,2,NULL,NULL,NULL,&dm);
192: DMSetFromOptions(dm);
193: DMSetUp(dm);
195: DMCreateGlobalVector(dm,&x_ref);
196: DMDAVecGenerateEntries(dm,x_ref);
198: if (!usempiio) {
199: MyVecDump("dmda.pbvec",skipheader,PETSC_FALSE,x_ref);
200: } else {
201: MyVecDump("dmda-mpiio.pbvec",skipheader,PETSC_TRUE,x_ref);
202: }
204: DMCreateGlobalVector(dm,&x_test);
206: if (!usempiio) {
207: MyVecLoad("dmda.pbvec",skipheader,usempiio,x_test);
208: } else {
209: MyVecLoad("dmda-mpiio.pbvec",skipheader,usempiio,x_test);
210: }
212: VecCompare(x_ref,x_test);
214: if (!usempiio) {
215: HeaderlessBinaryReadCheck(dm,"dmda.pbvec");
216: } else {
217: HeaderlessBinaryReadCheck(dm,"dmda-mpiio.pbvec");
218: }
219: VecDestroy(&x_ref);
220: VecDestroy(&x_test);
221: DMDestroy(&dm);
222: return 0;
223: }
225: int main(int argc,char **args)
226: {
227: PetscBool usempiio = PETSC_FALSE;
229: PetscInitialize(&argc,&args,(char*)0,help);
230: PetscOptionsGetBool(NULL,NULL,"-usempiio",&usempiio,NULL);
231: if (!usempiio) {
232: TestDMDAVec(PETSC_FALSE);
233: } else {
234: #if defined(PETSC_HAVE_MPIIO)
235: TestDMDAVec(PETSC_TRUE);
236: #else
237: PetscPrintf(PETSC_COMM_WORLD,"Warning: Executing TestDMDAVec(PETSC_TRUE) requires a working MPI-2 implementation\n");
238: #endif
239: }
240: PetscFinalize();
241: return 0;
242: }
244: /*TEST
246: test:
248: test:
249: suffix: 2
250: nsize: 12
252: test:
253: suffix: 3
254: nsize: 12
255: requires: defined(PETSC_HAVE_MPIIO)
256: args: -usempiio
258: TEST*/