Actual source code: vecio.c
2: /*
3: This file contains simple binary input routines for vectors. The
4: analogous output routines are within each vector implementation's
5: VecView (with viewer types PETSCVIEWERBINARY)
6: */
8: #include <petscsys.h>
9: #include <petscvec.h>
10: #include <petsc/private/vecimpl.h>
11: #include <petsc/private/viewerimpl.h>
12: #include <petsclayouthdf5.h>
14: PetscErrorCode VecView_Binary(Vec vec,PetscViewer viewer)
15: {
16: PetscBool skipHeader;
17: PetscLayout map;
18: PetscInt tr[2],n,s,N;
19: const PetscScalar *xarray;
22: PetscViewerSetUp(viewer);
23: PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
25: VecGetLayout(vec,&map);
26: PetscLayoutGetLocalSize(map,&n);
27: PetscLayoutGetRange(map,&s,NULL);
28: PetscLayoutGetSize(map,&N);
30: tr[0] = VEC_FILE_CLASSID; tr[1] = N;
31: if (!skipHeader) PetscViewerBinaryWrite(viewer,tr,2,PETSC_INT);
33: VecGetArrayRead(vec,&xarray);
34: PetscViewerBinaryWriteAll(viewer,xarray,n,s,N,PETSC_SCALAR);
35: VecRestoreArrayRead(vec,&xarray);
37: { /* write to the viewer's .info file */
38: FILE *info;
39: PetscMPIInt rank;
40: PetscViewerFormat format;
41: const char *name,*pre;
43: PetscViewerBinaryGetInfoPointer(viewer,&info);
44: MPI_Comm_rank(PetscObjectComm((PetscObject)vec),&rank);
46: PetscViewerGetFormat(viewer,&format);
47: if (format == PETSC_VIEWER_BINARY_MATLAB) {
48: PetscObjectGetName((PetscObject)vec,&name);
49: if (rank == 0 && info) {
50: PetscFPrintf(PETSC_COMM_SELF,info,"#--- begin code written by PetscViewerBinary for MATLAB format ---#\n");
51: PetscFPrintf(PETSC_COMM_SELF,info,"#$$ Set.%s = PetscBinaryRead(fd);\n",name);
52: PetscFPrintf(PETSC_COMM_SELF,info,"#--- end code written by PetscViewerBinary for MATLAB format ---#\n\n");
53: }
54: }
56: PetscObjectGetOptionsPrefix((PetscObject)vec,&pre);
57: if (rank == 0 && info) PetscFPrintf(PETSC_COMM_SELF,info,"-%svecload_block_size %" PetscInt_FMT "\n",pre?pre:"",PetscAbs(vec->map->bs));
58: }
59: return 0;
60: }
62: PetscErrorCode VecLoad_Binary(Vec vec, PetscViewer viewer)
63: {
64: PetscBool skipHeader,flg;
65: PetscInt tr[2],rows,N,n,s,bs;
66: PetscScalar *array;
67: PetscLayout map;
69: PetscViewerSetUp(viewer);
70: PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
72: VecGetLayout(vec,&map);
73: PetscLayoutGetSize(map,&N);
75: /* read vector header */
76: if (!skipHeader) {
77: PetscViewerBinaryRead(viewer,tr,2,NULL,PETSC_INT);
81: rows = tr[1];
82: } else {
84: rows = N;
85: }
87: /* set vector size, blocksize, and type if not already set; block size first so that local sizes will be compatible. */
88: PetscLayoutGetBlockSize(map,&bs);
89: PetscOptionsGetInt(((PetscObject)viewer)->options,((PetscObject)vec)->prefix,"-vecload_block_size",&bs,&flg);
90: if (flg) VecSetBlockSize(vec,bs);
91: PetscLayoutGetLocalSize(map,&n);
92: if (N < 0) VecSetSizes(vec,n,rows);
93: VecSetUp(vec);
95: /* get vector sizes and check global size */
96: VecGetSize(vec,&N);
97: VecGetLocalSize(vec,&n);
98: VecGetOwnershipRange(vec,&s,NULL);
101: /* read vector values */
102: VecGetArray(vec,&array);
103: PetscViewerBinaryReadAll(viewer,array,n,s,N,PETSC_SCALAR);
104: VecRestoreArray(vec,&array);
105: return 0;
106: }
108: #if defined(PETSC_HAVE_HDF5)
109: PetscErrorCode VecLoad_HDF5(Vec xin, PetscViewer viewer)
110: {
111: hid_t scalartype; /* scalar type (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */
112: PetscScalar *x,*arr;
113: const char *vecname;
116: #if defined(PETSC_USE_REAL_SINGLE)
117: scalartype = H5T_NATIVE_FLOAT;
118: #elif defined(PETSC_USE_REAL___FLOAT128)
119: #error "HDF5 output with 128 bit floats not supported."
120: #elif defined(PETSC_USE_REAL___FP16)
121: #error "HDF5 output with 16 bit floats not supported."
122: #else
123: scalartype = H5T_NATIVE_DOUBLE;
124: #endif
125: PetscObjectGetName((PetscObject)xin, &vecname);
126: PetscViewerHDF5Load(viewer,vecname,xin->map,scalartype,(void**)&x);
127: VecSetUp(xin); /* VecSetSizes might have not been called so ensure ops->create has been called */
128: if (!xin->ops->replacearray) {
129: VecGetArray(xin,&arr);
130: PetscArraycpy(arr,x,xin->map->n);
131: PetscFree(x);
132: VecRestoreArray(xin,&arr);
133: } else {
134: VecReplaceArray(xin,x);
135: }
136: return 0;
137: }
138: #endif
140: #if defined(PETSC_HAVE_ADIOS)
141: #include <adios.h>
142: #include <adios_read.h>
143: #include <petsc/private/vieweradiosimpl.h>
144: #include <petsc/private/viewerimpl.h>
146: PetscErrorCode VecLoad_ADIOS(Vec xin, PetscViewer viewer)
147: {
148: PetscViewer_ADIOS *adios = (PetscViewer_ADIOS*)viewer->data;
149: PetscScalar *x;
150: PetscInt Nfile,N,rstart,n;
151: uint64_t N_t,rstart_t;
152: const char *vecname;
153: ADIOS_VARINFO *v;
154: ADIOS_SELECTION *sel;
156: PetscObjectGetName((PetscObject) xin, &vecname);
158: v = adios_inq_var(adios->adios_fp, vecname);
160: Nfile = (PetscInt) v->dims[0];
162: /* Set Vec sizes,blocksize,and type if not already set */
163: if ((xin)->map->n < 0 && (xin)->map->N < 0) VecSetSizes(xin, PETSC_DECIDE, Nfile);
164: /* If sizes and type already set,check if the vector global size is correct */
165: VecGetSize(xin, &N);
166: VecGetLocalSize(xin, &n);
169: VecGetOwnershipRange(xin,&rstart,NULL);
170: rstart_t = rstart;
171: N_t = n;
172: sel = adios_selection_boundingbox (v->ndim, &rstart_t, &N_t);
173: VecGetArray(xin,&x);
174: adios_schedule_read(adios->adios_fp, sel, vecname, 0, 1, x);
175: adios_perform_reads (adios->adios_fp, 1);
176: VecRestoreArray(xin,&x);
177: adios_selection_delete(sel);
179: return 0;
180: }
181: #endif
183: PetscErrorCode VecLoad_Default(Vec newvec, PetscViewer viewer)
184: {
185: PetscBool isbinary;
186: #if defined(PETSC_HAVE_HDF5)
187: PetscBool ishdf5;
188: #endif
189: #if defined(PETSC_HAVE_ADIOS)
190: PetscBool isadios;
191: #endif
193: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
194: #if defined(PETSC_HAVE_HDF5)
195: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
196: #endif
197: #if defined(PETSC_HAVE_ADIOS)
198: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERADIOS,&isadios);
199: #endif
201: #if defined(PETSC_HAVE_HDF5)
202: if (ishdf5) {
203: if (!((PetscObject)newvec)->name) {
204: PetscLogEventEnd(VEC_Load,viewer,0,0,0);
205: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Since HDF5 format gives ASCII name for each object in file; must use VecLoad() after setting name of Vec with PetscObjectSetName()");
206: }
207: VecLoad_HDF5(newvec, viewer);
208: } else
209: #endif
210: #if defined(PETSC_HAVE_ADIOS)
211: if (isadios) {
212: if (!((PetscObject)newvec)->name) {
213: PetscLogEventEnd(VEC_Load,viewer,0,0,0);
214: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Since ADIOS format gives ASCII name for each object in file; must use VecLoad() after setting name of Vec with PetscObjectSetName()");
215: }
216: VecLoad_ADIOS(newvec, viewer);
217: } else
218: #endif
219: {
220: VecLoad_Binary(newvec, viewer);
221: }
222: return 0;
223: }
225: /*@
226: VecChop - Set all values in the vector with an absolute value less than the tolerance to zero
228: Input Parameters:
229: + v - The vector
230: - tol - The zero tolerance
232: Output Parameters:
233: . v - The chopped vector
235: Level: intermediate
237: .seealso: VecCreate(), VecSet()
238: @*/
239: PetscErrorCode VecChop(Vec v, PetscReal tol)
240: {
241: PetscScalar *a;
242: PetscInt n, i;
244: VecGetLocalSize(v, &n);
245: VecGetArray(v, &a);
246: for (i = 0; i < n; ++i) {
247: if (PetscAbsScalar(a[i]) < tol) a[i] = 0.0;
248: }
249: VecRestoreArray(v, &a);
250: return 0;
251: }