Actual source code: ex15.c
petsc-3.6.1 2015-08-06
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 PetscScalar dmda_i_val[] = { 1.10, 2.3006, 2.32444, 3.44006, 66.9009 };
13: const PetscScalar dmda_j_val[] = { 0.0, 0.25, 0.5, 0.75 };
14: const PetscScalar dmda_k_val[] = { 0.0, 1.1, 2.2, 3.3, 4.4, 5.5 };
18: PetscErrorCode MyVecDump(const char fname[],PetscBool skippheader,PetscBool usempiio,Vec x)
19: {
20: MPI_Comm comm;
21: PetscViewer viewer;
22: PetscBool ismpiio,isskip;
24:
26: PetscObjectGetComm((PetscObject)x,&comm);
28: PetscViewerCreate(comm,&viewer);
29: PetscViewerSetType(viewer,PETSCVIEWERBINARY);
30: if (skippheader) { PetscViewerBinarySetSkipHeader(viewer,PETSC_TRUE); }
31: PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
32: if (usempiio) { PetscViewerBinarySetUseMPIIO(viewer,PETSC_TRUE); }
33: PetscViewerFileSetName(viewer,fname);
35: VecView(x,viewer);
37: PetscViewerBinaryGetUseMPIIO(viewer,&ismpiio);
38: if (ismpiio) { PetscPrintf(comm,"*** PetscViewer[write] using MPI-IO ***\n"); }
39: PetscViewerBinaryGetSkipHeader(viewer,&isskip);
40: if (isskip) { PetscPrintf(comm,"*** PetscViewer[write] skipping header ***\n"); }
42: PetscViewerDestroy(&viewer);
43: return(0);
44: }
48: PetscErrorCode MyVecLoad(const char fname[],PetscBool skippheader,PetscBool usempiio,Vec x)
49: {
50: MPI_Comm comm;
51: PetscViewer viewer;
52: PetscBool ismpiio,isskip;
56: PetscObjectGetComm((PetscObject)x,&comm);
58: PetscViewerCreate(comm,&viewer);
59: PetscViewerSetType(viewer,PETSCVIEWERBINARY);
60: if (skippheader) { PetscViewerBinarySetSkipHeader(viewer,PETSC_TRUE); }
61: PetscViewerFileSetMode(viewer,FILE_MODE_READ);
62: if (usempiio) { PetscViewerBinarySetUseMPIIO(viewer,PETSC_TRUE); }
63: PetscViewerFileSetName(viewer,fname);
65: VecLoad(x,viewer);
67: PetscViewerBinaryGetSkipHeader(viewer,&isskip);
68: if (isskip) { PetscPrintf(comm,"*** PetscViewer[load] skipping header ***\n"); }
69: PetscViewerBinaryGetUseMPIIO(viewer,&ismpiio);
70: if (ismpiio) { PetscPrintf(comm,"*** PetscViewer[load] using MPI-IO ***\n"); }
72: PetscViewerDestroy(&viewer);
73: return(0);
74: }
78: PetscErrorCode DMDAVecGenerateEntries(DM dm,Vec a)
79: {
80: PetscScalar ****LA_v;
81: PetscInt i,j,k,si,sj,sk,ni,nj,nk,M,N;
85: DMDAGetInfo(dm,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
86: DMDAGetCorners(dm,&si,&sj,&sk,&ni,&nj,&nk);
87: DMDAVecGetArrayDOF(dm,a,&LA_v);
88: for (k=sk; k<sk+nk; k++) {
89: for (j=sj; j<sj+nj; j++) {
90: for (i=si; i<si+ni; i++) {
91: PetscScalar test_value_s;
92:
93: 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));
94: LA_v[k][j][i][0] = 3.0 * test_value_s;
95: LA_v[k][j][i][1] = 3.0 * test_value_s + 1.0;
96: LA_v[k][j][i][2] = 3.0 * test_value_s + 2.0;
97: }
98: }
99: }
100: DMDAVecRestoreArrayDOF(dm,a,&LA_v);
101: return(0);
102: }
106: PetscErrorCode HeaderlessBinaryReadCheck(DM dm,const char name[])
107: {
109: int fdes;
110: PetscScalar buffer[DMDA_I*DMDA_J*DMDA_K*3];
111: PetscInt len,d,i,j,k,M,N;
112: PetscMPIInt rank;
113: PetscBool dataverified = PETSC_TRUE;
116: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
117: DMDAGetInfo(dm,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
118: len = DMDA_I*DMDA_J*DMDA_K*3;
119: if (!rank) {
120: PetscBinaryOpen(name,FILE_MODE_READ,&fdes);
121: PetscBinaryRead(fdes,buffer,len,PETSC_SCALAR);
122: PetscBinaryClose(fdes);
124: for (k=0; k<DMDA_K; k++) {
125: for (j=0; j<DMDA_J; j++) {
126: for (i=0; i<DMDA_I; i++) {
127: for (d=0; d<3; d++) {
128: PetscScalar v,test_value_s,test_value;
129: PetscInt index;
131: 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));
132: test_value = 3.0 * test_value_s + (PetscScalar)d;
134: index = 3*(i + j*M + k*M*N) + d;
135: v = PetscAbsScalar(test_value-buffer[index]);
136: #if defined(PETSC_USE_COMPLEX)
137: if ((PetscRealPart(v) > 1.0e-10) || (PetscImaginaryPart(v) > 1.0e-10)) {
138: 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);
139: dataverified = PETSC_FALSE;
140: }
141: #else
142: if (PetscRealPart(v) > 1.0e-10) {
143: 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);
144: dataverified = PETSC_FALSE;
145: }
146: #endif
147: }
148: }
149: }
150: }
151: if (dataverified) {
152: PetscPrintf(PETSC_COMM_SELF,"Headerless read of data verified for: %s\n",name);
153: }
154: }
155: return(0);
156: }
160: PetscErrorCode VecCompare(Vec a,Vec b)
161: {
162: PetscInt locmin[2],locmax[2];
163: PetscReal min[2],max[2];
164: Vec ref;
168: VecMin(a,&locmin[0],&min[0]);
169: VecMax(a,&locmax[0],&max[0]);
171: VecMin(b,&locmin[1],&min[1]);
172: VecMax(b,&locmax[1],&max[1]);
174: PetscPrintf(PETSC_COMM_WORLD,"VecCompare\n");
175: PetscPrintf(PETSC_COMM_WORLD," min(a) = %+1.2e [loc %D]\n",(double)min[0],locmin[0]);
176: PetscPrintf(PETSC_COMM_WORLD," max(a) = %+1.2e [loc %D]\n",(double)max[0],locmax[0]);
178: PetscPrintf(PETSC_COMM_WORLD," min(b) = %+1.2e [loc %D]\n",(double)min[1],locmin[1]);
179: PetscPrintf(PETSC_COMM_WORLD," max(b) = %+1.2e [loc %D]\n",(double)max[1],locmax[1]);
181: VecDuplicate(a,&ref);
182: VecCopy(a,ref);
183: VecAXPY(ref,-1.0,b);
184: VecMin(ref,&locmin[0],&min[0]);
185: if (PetscAbsReal(min[0]) > 1.0e-10) {
186: PetscPrintf(PETSC_COMM_WORLD," ERROR: min(a-b) > 1.0e-10\n");
187: PetscPrintf(PETSC_COMM_WORLD," min(a-b) = %+1.10e\n",(double)PetscAbsReal(min[0]));
188: } else {
189: PetscPrintf(PETSC_COMM_WORLD," min(a-b) < 1.0e-10\n");
190: }
191: VecDestroy(&ref);
192: return(0);
193: }
197: PetscErrorCode TestDMDAVec(PetscBool usempiio)
198: {
199: DM dm;
200: Vec x_ref,x_test;
201: PetscBool skipheader = PETSC_TRUE;
205: if (!usempiio) { PetscPrintf(PETSC_COMM_WORLD,"%s\n",__FUNCT__); }
206: else { PetscPrintf(PETSC_COMM_WORLD,"%s [using mpi-io]\n",__FUNCT__); }
207: DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,
208: DMDA_STENCIL_BOX,DMDA_I,DMDA_J,DMDA_K,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,
209: 3,2,NULL,NULL,NULL,&dm);
211: DMCreateGlobalVector(dm,&x_ref);
212: DMDAVecGenerateEntries(dm,x_ref);
214: if (!usempiio) {
215: MyVecDump("dmda.pbvec",skipheader,PETSC_FALSE,x_ref);
216: } else {
217: MyVecDump("dmda-mpiio.pbvec",skipheader,PETSC_TRUE,x_ref);
218: }
220: DMCreateGlobalVector(dm,&x_test);
222: if (!usempiio) {
223: MyVecLoad("dmda.pbvec",skipheader,usempiio,x_test);
224: } else {
225: MyVecLoad("dmda-mpiio.pbvec",skipheader,usempiio,x_test);
226: }
228: VecCompare(x_ref,x_test);
230: if (!usempiio) {
231: HeaderlessBinaryReadCheck(dm,"dmda.pbvec");
232: } else {
233: HeaderlessBinaryReadCheck(dm,"dmda-mpiio.pbvec");
234: }
235: VecDestroy(&x_ref);
236: VecDestroy(&x_test);
237: DMDestroy(&dm);
238: return(0);
239: }
243: int main(int argc,char **args)
244: {
246: PetscBool usempiio = PETSC_FALSE;
248: PetscInitialize(&argc,&args,(char*)0,help);
250: PetscOptionsGetBool(NULL,"-usempiio",&usempiio,NULL);
251: if (!usempiio) {
252: TestDMDAVec(PETSC_FALSE);
253: } else {
254: #if defined(PETSC_HAVE_MPIIO)
255: TestDMDAVec(PETSC_TRUE);
256: #else
257: PetscPrintf(PETSC_COMM_WORLD,"Warning: Executing TestDMDAVec(PETSC_TRUE) requires a working MPI-2 implementation\n");
258: #endif
259: }
260: PetscFinalize();
261: return 0;
262: }