Actual source code: ex3.c
petsc-3.5.4 2015-05-23
2: static char help[] = "Tests PetscViewerBinary VecView()/VecLoad() function correctly when binary header is skipped.\n\n";
4: /*T
5: Concepts: viewers^skipheader
6: T*/
8: #include <petscviewer.h>
9: #include <petscvec.h>
11: #define VEC_LEN 10
12: const PetscScalar test_values[] = { 0.311256, 88.068, 11.077444, 9953.62, 7.345, 64.8943, 3.1458, 6699.95, 0.00084, 0.0647 };
16: PetscErrorCode MyVecDump(const char fname[],PetscBool skippheader,PetscBool usempiio,Vec x)
17: {
18: MPI_Comm comm;
19: PetscViewer viewer;
20: PetscBool ismpiio,isskip;
24: PetscObjectGetComm((PetscObject)x,&comm);
26: PetscViewerCreate(comm,&viewer);
27: PetscViewerSetType(viewer,PETSCVIEWERBINARY);
28: if (skippheader) { PetscViewerBinarySetSkipHeader(viewer,PETSC_TRUE); }
29: PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
30: if (usempiio) { PetscViewerBinarySetMPIIO(viewer); }
31: PetscViewerFileSetName(viewer,fname);
33: VecView(x,viewer);
35: PetscViewerBinaryGetMPIIO(viewer,&ismpiio);
36: if (ismpiio) { PetscPrintf(comm,"*** PetscViewer[write] using MPI-IO ***\n"); }
37: PetscViewerBinaryGetSkipHeader(viewer,&isskip);
38: if (isskip) { PetscPrintf(comm,"*** PetscViewer[write] skipping header ***\n"); }
40: PetscViewerDestroy(&viewer);
41: return(0);
42: }
46: PetscErrorCode MyVecLoad(const char fname[],PetscBool skippheader,PetscBool usempiio,Vec x)
47: {
48: MPI_Comm comm;
49: PetscViewer viewer;
50: PetscBool ismpiio,isskip;
54: PetscObjectGetComm((PetscObject)x,&comm);
56: PetscViewerCreate(comm,&viewer);
57: PetscViewerSetType(viewer,PETSCVIEWERBINARY);
58: if (skippheader) { PetscViewerBinarySetSkipHeader(viewer,PETSC_TRUE); }
59: PetscViewerFileSetMode(viewer,FILE_MODE_READ);
60: if (usempiio) { PetscViewerBinarySetMPIIO(viewer); }
61: PetscViewerFileSetName(viewer,fname);
63: VecLoad(x,viewer);
65: PetscViewerBinaryGetSkipHeader(viewer,&isskip);
66: if (isskip) { PetscPrintf(comm,"*** PetscViewer[load] skipping header ***\n"); }
67: PetscViewerBinaryGetMPIIO(viewer,&ismpiio);
68: if (ismpiio) { PetscPrintf(comm,"*** PetscViewer[load] using MPI-IO ***\n"); }
70: PetscViewerDestroy(&viewer);
71: return(0);
72: }
76: PetscErrorCode VecFill(Vec x)
77: {
79: PetscInt i,s,e;
82: VecGetOwnershipRange(x,&s,&e);
83: for (i=s; i<e; i++) {
84: VecSetValue(x,i,test_values[i],INSERT_VALUES);
85: }
86: VecAssemblyBegin(x);
87: VecAssemblyEnd(x);
88: return(0);
89: }
93: PetscErrorCode VecCompare(Vec a,Vec b)
94: {
95: PetscInt locmin[2],locmax[2];
96: PetscReal min[2],max[2];
97: Vec ref;
101: VecMin(a,&locmin[0],&min[0]);
102: VecMax(a,&locmax[0],&max[0]);
104: VecMin(b,&locmin[1],&min[1]);
105: VecMax(b,&locmax[1],&max[1]);
107: PetscPrintf(PETSC_COMM_WORLD,"VecCompare\n");
108: PetscPrintf(PETSC_COMM_WORLD," min(a) = %+1.2e [loc %D]\n",(double)min[0],locmin[0]);
109: PetscPrintf(PETSC_COMM_WORLD," max(a) = %+1.2e [loc %D]\n",(double)max[0],locmax[0]);
111: PetscPrintf(PETSC_COMM_WORLD," min(b) = %+1.2e [loc %D]\n",(double)min[1],locmin[1]);
112: PetscPrintf(PETSC_COMM_WORLD," max(b) = %+1.2e [loc %D]\n",(double)max[1],locmax[1]);
114: VecDuplicate(a,&ref);
115: VecCopy(a,ref);
116: VecAXPY(ref,-1.0,b);
117: VecMin(ref,&locmin[0],&min[0]);
118: if (PetscAbsReal(min[0]) > 1.0e-10) {
119: PetscPrintf(PETSC_COMM_WORLD," ERROR: min(a-b) > 1.0e-10\n");
120: PetscPrintf(PETSC_COMM_WORLD," min(a-b) = %+1.10e\n",(double)PetscAbsReal(min[0]));
121: } else {
122: PetscPrintf(PETSC_COMM_WORLD," min(a-b) < 1.0e-10\n");
123: }
124: VecDestroy(&ref);
125: return(0);
126: }
130: PetscErrorCode HeaderlessBinaryRead(const char name[])
131: {
133: int fdes;
134: PetscScalar buffer[VEC_LEN];
135: PetscInt i;
136: PetscMPIInt rank;
137: PetscBool dataverified = PETSC_TRUE;
140: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
141: if (!rank) {
142: PetscBinaryOpen(name,FILE_MODE_READ,&fdes);
143: PetscBinaryRead(fdes,buffer,VEC_LEN,PETSC_SCALAR);
144: PetscBinaryClose(fdes);
146: for (i=0; i<VEC_LEN; i++) {
147: PetscScalar v;
148: v = PetscAbsScalar(test_values[i]-buffer[i]);
149: #if defined(PETSC_USE_COMPLEX)
150: if ((PetscRealPart(v) > 1.0e-10) || (PetscImaginaryPart(v) > 1.0e-10)) {
151: PetscPrintf(PETSC_COMM_SELF,"ERROR: Difference > 1.0e-10 occurred (delta = (%+1.12e,%+1.12e) [loc %D])\n",(double)PetscRealPart(buffer[i]),(double)PetscImaginaryPart(buffer[i]),i);
152: dataverified = PETSC_FALSE;
153: }
154: #else
155: if (PetscRealPart(v) > 1.0e-10) {
156: PetscPrintf(PETSC_COMM_SELF,"ERROR: Difference > 1.0e-10 occurred (delta = %+1.12e [loc %D])\n",(double)PetscRealPart(buffer[i]),i);
157: dataverified = PETSC_FALSE;
158: }
159: #endif
160: }
161: if (dataverified) {
162: PetscPrintf(PETSC_COMM_SELF,"Headerless read of data verified\n");
163: }
164: }
165: return(0);
166: }
170: PetscErrorCode TestBinary(void)
171: {
173: Vec x,y;
174: PetscBool skipheader = PETSC_TRUE;
175: PetscBool usempiio = PETSC_FALSE;
176:
178: VecCreate(PETSC_COMM_WORLD,&x);
179: VecSetSizes(x,PETSC_DECIDE,VEC_LEN);
180: VecSetFromOptions(x);
181: VecFill(x);
182: MyVecDump("xH.pbvec",skipheader,usempiio,x);
184: VecCreate(PETSC_COMM_WORLD,&y);
185: VecSetSizes(y,PETSC_DECIDE,VEC_LEN);
186: VecSetFromOptions(y);
188: MyVecLoad("xH.pbvec",skipheader,usempiio,y);
189: VecCompare(x,y);
191: VecDestroy(&y);
192: VecDestroy(&x);
194: HeaderlessBinaryRead("xH.pbvec");
195: return(0);
196: }
198: #if defined(PETSC_HAVE_MPIIO)
201: PetscErrorCode TestBinaryMPIIO(void)
202: {
204: Vec x,y;
205: PetscBool skipheader = PETSC_TRUE;
206: PetscBool usempiio = PETSC_TRUE;
209: VecCreate(PETSC_COMM_WORLD,&x);
210: VecSetSizes(x,PETSC_DECIDE,VEC_LEN);
211: VecSetFromOptions(x);
212: VecFill(x);
213: MyVecDump("xHmpi.pbvec",skipheader,usempiio,x);
215: VecCreate(PETSC_COMM_WORLD,&y);
216: VecSetSizes(y,PETSC_DECIDE,VEC_LEN);
217: VecSetFromOptions(y);
219: MyVecLoad("xHmpi.pbvec",skipheader,usempiio,y);
220: VecCompare(x,y);
222: VecDestroy(&y);
223: VecDestroy(&x);
225: HeaderlessBinaryRead("xHmpi.pbvec");
226: return(0);
227: }
228: #endif
232: int main(int argc,char **args)
233: {
235: PetscBool usempiio = PETSC_FALSE;
236:
237: PetscInitialize(&argc,&args,(char*)0,help);
239: PetscOptionsGetBool(NULL,"-usempiio",&usempiio,NULL);
240: if (!usempiio) {
241: TestBinary();
242: } else {
243: #if defined(PETSC_HAVE_MPIIO)
244: TestBinaryMPIIO();
245: #else
246: PetscPrintf(PETSC_COMM_WORLD,"Warning: Executing TestBinaryMPIIO() requires a working MPI-2 implementation\n");
247: #endif
248: }
249: PetscFinalize();
250: return 0;
251: }