Actual source code: pdvec.c
2: /*
3: Code for some of the parallel vector primitives.
4: */
5: #include <../src/vec/vec/impls/mpi/pvecimpl.h>
6: #include <petsc/private/viewerhdf5impl.h>
7: #include <petsc/private/glvisviewerimpl.h>
8: #include <petsc/private/glvisvecimpl.h>
9: #include <petscsf.h>
11: static PetscErrorCode VecResetPreallocationCOO_MPI(Vec v)
12: {
13: Vec_MPI *vmpi = (Vec_MPI *)v->data;
15: PetscFunctionBegin;
16: if (vmpi) {
17: PetscCall(PetscFree(vmpi->jmap1));
18: PetscCall(PetscFree(vmpi->perm1));
19: PetscCall(PetscFree(vmpi->Cperm));
20: PetscCall(PetscFree4(vmpi->imap2, vmpi->jmap2, vmpi->sendbuf, vmpi->recvbuf));
21: PetscCall(PetscFree(vmpi->perm2));
22: PetscCall(PetscSFDestroy(&vmpi->coo_sf));
23: }
24: PetscFunctionReturn(PETSC_SUCCESS);
25: }
27: PetscErrorCode VecDestroy_MPI(Vec v)
28: {
29: Vec_MPI *x = (Vec_MPI *)v->data;
31: PetscFunctionBegin;
32: #if defined(PETSC_USE_LOG)
33: PetscCall(PetscLogObjectState((PetscObject)v, "Length=%" PetscInt_FMT, v->map->N));
34: #endif
35: if (!x) PetscFunctionReturn(PETSC_SUCCESS);
36: PetscCall(PetscFree(x->array_allocated));
38: /* Destroy local representation of vector if it exists */
39: if (x->localrep) {
40: PetscCall(VecDestroy(&x->localrep));
41: PetscCall(VecScatterDestroy(&x->localupdate));
42: }
43: PetscCall(VecAssemblyReset_MPI(v));
45: /* Destroy the stashes: note the order - so that the tags are freed properly */
46: PetscCall(VecStashDestroy_Private(&v->bstash));
47: PetscCall(VecStashDestroy_Private(&v->stash));
49: PetscCall(VecResetPreallocationCOO_MPI(v));
50: PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscMatlabEnginePut_C", NULL));
51: PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscMatlabEngineGet_C", NULL));
52: PetscCall(PetscFree(v->data));
53: PetscFunctionReturn(PETSC_SUCCESS);
54: }
56: PetscErrorCode VecView_MPI_ASCII(Vec xin, PetscViewer viewer)
57: {
58: PetscInt i, work = xin->map->n, cnt, len, nLen;
59: PetscMPIInt j, n = 0, size, rank, tag = ((PetscObject)viewer)->tag;
60: MPI_Status status;
61: PetscScalar *values;
62: const PetscScalar *xarray;
63: const char *name;
64: PetscViewerFormat format;
66: PetscFunctionBegin;
67: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)xin), &size));
68: PetscCall(PetscViewerGetFormat(viewer, &format));
69: if (format == PETSC_VIEWER_LOAD_BALANCE) {
70: PetscInt nmax = 0, nmin = xin->map->n, navg;
71: for (i = 0; i < (PetscInt)size; i++) {
72: nmax = PetscMax(nmax, xin->map->range[i + 1] - xin->map->range[i]);
73: nmin = PetscMin(nmin, xin->map->range[i + 1] - xin->map->range[i]);
74: }
75: navg = xin->map->N / size;
76: PetscCall(PetscViewerASCIIPrintf(viewer, " Load Balance - Local vector size Min %" PetscInt_FMT " avg %" PetscInt_FMT " max %" PetscInt_FMT "\n", nmin, navg, nmax));
77: PetscFunctionReturn(PETSC_SUCCESS);
78: }
80: PetscCall(VecGetArrayRead(xin, &xarray));
81: /* determine maximum message to arrive */
82: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)xin), &rank));
83: PetscCallMPI(MPI_Reduce(&work, &len, 1, MPIU_INT, MPI_MAX, 0, PetscObjectComm((PetscObject)xin)));
84: if (format == PETSC_VIEWER_ASCII_GLVIS) rank = 0, len = 0; /* no parallel distributed write support from GLVis */
85: if (rank == 0) {
86: PetscCall(PetscMalloc1(len, &values));
87: /*
88: MATLAB format and ASCII format are very similar except
89: MATLAB uses %18.16e format while ASCII uses %g
90: */
91: if (format == PETSC_VIEWER_ASCII_MATLAB) {
92: PetscCall(PetscObjectGetName((PetscObject)xin, &name));
93: PetscCall(PetscViewerASCIIPrintf(viewer, "%s = [\n", name));
94: for (i = 0; i < xin->map->n; i++) {
95: #if defined(PETSC_USE_COMPLEX)
96: if (PetscImaginaryPart(xarray[i]) > 0.0) {
97: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e + %18.16ei\n", (double)PetscRealPart(xarray[i]), (double)PetscImaginaryPart(xarray[i])));
98: } else if (PetscImaginaryPart(xarray[i]) < 0.0) {
99: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e - %18.16ei\n", (double)PetscRealPart(xarray[i]), -(double)PetscImaginaryPart(xarray[i])));
100: } else {
101: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e\n", (double)PetscRealPart(xarray[i])));
102: }
103: #else
104: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e\n", (double)xarray[i]));
105: #endif
106: }
107: /* receive and print messages */
108: for (j = 1; j < size; j++) {
109: PetscCallMPI(MPI_Recv(values, (PetscMPIInt)len, MPIU_SCALAR, j, tag, PetscObjectComm((PetscObject)xin), &status));
110: PetscCallMPI(MPI_Get_count(&status, MPIU_SCALAR, &n));
111: for (i = 0; i < n; i++) {
112: #if defined(PETSC_USE_COMPLEX)
113: if (PetscImaginaryPart(values[i]) > 0.0) {
114: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e + %18.16ei\n", (double)PetscRealPart(values[i]), (double)PetscImaginaryPart(values[i])));
115: } else if (PetscImaginaryPart(values[i]) < 0.0) {
116: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e - %18.16ei\n", (double)PetscRealPart(values[i]), -(double)PetscImaginaryPart(values[i])));
117: } else {
118: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e\n", (double)PetscRealPart(values[i])));
119: }
120: #else
121: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e\n", (double)values[i]));
122: #endif
123: }
124: }
125: PetscCall(PetscViewerASCIIPrintf(viewer, "];\n"));
127: } else if (format == PETSC_VIEWER_ASCII_SYMMODU) {
128: for (i = 0; i < xin->map->n; i++) {
129: #if defined(PETSC_USE_COMPLEX)
130: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e %18.16e\n", (double)PetscRealPart(xarray[i]), (double)PetscImaginaryPart(xarray[i])));
131: #else
132: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e\n", (double)xarray[i]));
133: #endif
134: }
135: /* receive and print messages */
136: for (j = 1; j < size; j++) {
137: PetscCallMPI(MPI_Recv(values, (PetscMPIInt)len, MPIU_SCALAR, j, tag, PetscObjectComm((PetscObject)xin), &status));
138: PetscCallMPI(MPI_Get_count(&status, MPIU_SCALAR, &n));
139: for (i = 0; i < n; i++) {
140: #if defined(PETSC_USE_COMPLEX)
141: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e %18.16e\n", (double)PetscRealPart(values[i]), (double)PetscImaginaryPart(values[i])));
142: #else
143: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e\n", (double)values[i]));
144: #endif
145: }
146: }
147: } else if (format == PETSC_VIEWER_ASCII_VTK_DEPRECATED || format == PETSC_VIEWER_ASCII_VTK_CELL_DEPRECATED) {
148: /*
149: state 0: No header has been output
150: state 1: Only POINT_DATA has been output
151: state 2: Only CELL_DATA has been output
152: state 3: Output both, POINT_DATA last
153: state 4: Output both, CELL_DATA last
154: */
155: static PetscInt stateId = -1;
156: int outputState = 0;
157: int doOutput = 0;
158: PetscBool hasState;
159: PetscInt bs, b;
161: if (stateId < 0) PetscCall(PetscObjectComposedDataRegister(&stateId));
162: PetscCall(PetscObjectComposedDataGetInt((PetscObject)viewer, stateId, outputState, hasState));
163: if (!hasState) outputState = 0;
165: PetscCall(PetscObjectGetName((PetscObject)xin, &name));
166: PetscCall(VecGetLocalSize(xin, &nLen));
167: PetscCall(PetscMPIIntCast(nLen, &n));
168: PetscCall(VecGetBlockSize(xin, &bs));
169: if (format == PETSC_VIEWER_ASCII_VTK_DEPRECATED) {
170: if (outputState == 0) {
171: outputState = 1;
172: doOutput = 1;
173: } else if (outputState == 1) doOutput = 0;
174: else if (outputState == 2) {
175: outputState = 3;
176: doOutput = 1;
177: } else if (outputState == 3) doOutput = 0;
178: else PetscCheck(outputState != 4, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Tried to output POINT_DATA again after intervening CELL_DATA");
180: if (doOutput) PetscCall(PetscViewerASCIIPrintf(viewer, "POINT_DATA %" PetscInt_FMT "\n", xin->map->N / bs));
181: } else {
182: if (outputState == 0) {
183: outputState = 2;
184: doOutput = 1;
185: } else if (outputState == 1) {
186: outputState = 4;
187: doOutput = 1;
188: } else if (outputState == 2) {
189: doOutput = 0;
190: } else {
191: PetscCheck(outputState != 3, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Tried to output CELL_DATA again after intervening POINT_DATA");
192: if (outputState == 4) doOutput = 0;
193: }
195: if (doOutput) PetscCall(PetscViewerASCIIPrintf(viewer, "CELL_DATA %" PetscInt_FMT "\n", xin->map->N / bs));
196: }
197: PetscCall(PetscObjectComposedDataSetInt((PetscObject)viewer, stateId, outputState));
198: if (name) {
199: if (bs == 3) {
200: PetscCall(PetscViewerASCIIPrintf(viewer, "VECTORS %s double\n", name));
201: } else {
202: PetscCall(PetscViewerASCIIPrintf(viewer, "SCALARS %s double %" PetscInt_FMT "\n", name, bs));
203: }
204: } else {
205: PetscCall(PetscViewerASCIIPrintf(viewer, "SCALARS scalars double %" PetscInt_FMT "\n", bs));
206: }
207: if (bs != 3) PetscCall(PetscViewerASCIIPrintf(viewer, "LOOKUP_TABLE default\n"));
208: for (i = 0; i < n / bs; i++) {
209: for (b = 0; b < bs; b++) {
210: if (b > 0) PetscCall(PetscViewerASCIIPrintf(viewer, " "));
211: PetscCall(PetscViewerASCIIPrintf(viewer, "%g", (double)PetscRealPart(xarray[i * bs + b])));
212: }
213: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
214: }
215: for (j = 1; j < size; j++) {
216: PetscCallMPI(MPI_Recv(values, (PetscMPIInt)len, MPIU_SCALAR, j, tag, PetscObjectComm((PetscObject)xin), &status));
217: PetscCallMPI(MPI_Get_count(&status, MPIU_SCALAR, &n));
218: for (i = 0; i < n / bs; i++) {
219: for (b = 0; b < bs; b++) {
220: if (b > 0) PetscCall(PetscViewerASCIIPrintf(viewer, " "));
221: PetscCall(PetscViewerASCIIPrintf(viewer, "%g", (double)PetscRealPart(values[i * bs + b])));
222: }
223: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
224: }
225: }
226: } else if (format == PETSC_VIEWER_ASCII_VTK_COORDS_DEPRECATED) {
227: PetscInt bs, b;
229: PetscCall(VecGetLocalSize(xin, &nLen));
230: PetscCall(PetscMPIIntCast(nLen, &n));
231: PetscCall(VecGetBlockSize(xin, &bs));
232: PetscCheck(bs >= 1 && bs <= 3, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "VTK can only handle 3D objects, but vector dimension is %" PetscInt_FMT, bs);
234: for (i = 0; i < n / bs; i++) {
235: for (b = 0; b < bs; b++) {
236: if (b > 0) PetscCall(PetscViewerASCIIPrintf(viewer, " "));
237: PetscCall(PetscViewerASCIIPrintf(viewer, "%g", (double)PetscRealPart(xarray[i * bs + b])));
238: }
239: for (b = bs; b < 3; b++) PetscCall(PetscViewerASCIIPrintf(viewer, " 0.0"));
240: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
241: }
242: for (j = 1; j < size; j++) {
243: PetscCallMPI(MPI_Recv(values, (PetscMPIInt)len, MPIU_SCALAR, j, tag, PetscObjectComm((PetscObject)xin), &status));
244: PetscCallMPI(MPI_Get_count(&status, MPIU_SCALAR, &n));
245: for (i = 0; i < n / bs; i++) {
246: for (b = 0; b < bs; b++) {
247: if (b > 0) PetscCall(PetscViewerASCIIPrintf(viewer, " "));
248: PetscCall(PetscViewerASCIIPrintf(viewer, "%g", (double)PetscRealPart(values[i * bs + b])));
249: }
250: for (b = bs; b < 3; b++) PetscCall(PetscViewerASCIIPrintf(viewer, " 0.0"));
251: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
252: }
253: }
254: } else if (format == PETSC_VIEWER_ASCII_PCICE) {
255: PetscInt bs, b, vertexCount = 1;
257: PetscCall(VecGetLocalSize(xin, &nLen));
258: PetscCall(PetscMPIIntCast(nLen, &n));
259: PetscCall(VecGetBlockSize(xin, &bs));
260: PetscCheck(bs >= 1 && bs <= 3, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "PCICE can only handle up to 3D objects, but vector dimension is %" PetscInt_FMT, bs);
262: PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT "\n", xin->map->N / bs));
263: for (i = 0; i < n / bs; i++) {
264: PetscCall(PetscViewerASCIIPrintf(viewer, "%7" PetscInt_FMT " ", vertexCount++));
265: for (b = 0; b < bs; b++) {
266: if (b > 0) PetscCall(PetscViewerASCIIPrintf(viewer, " "));
267: #if !defined(PETSC_USE_COMPLEX)
268: PetscCall(PetscViewerASCIIPrintf(viewer, "% 12.5E", (double)xarray[i * bs + b]));
269: #endif
270: }
271: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
272: }
273: for (j = 1; j < size; j++) {
274: PetscCallMPI(MPI_Recv(values, (PetscMPIInt)len, MPIU_SCALAR, j, tag, PetscObjectComm((PetscObject)xin), &status));
275: PetscCallMPI(MPI_Get_count(&status, MPIU_SCALAR, &n));
276: for (i = 0; i < n / bs; i++) {
277: PetscCall(PetscViewerASCIIPrintf(viewer, "%7" PetscInt_FMT " ", vertexCount++));
278: for (b = 0; b < bs; b++) {
279: if (b > 0) PetscCall(PetscViewerASCIIPrintf(viewer, " "));
280: #if !defined(PETSC_USE_COMPLEX)
281: PetscCall(PetscViewerASCIIPrintf(viewer, "% 12.5E", (double)values[i * bs + b]));
282: #endif
283: }
284: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
285: }
286: }
287: } else if (format == PETSC_VIEWER_ASCII_GLVIS) {
288: /* GLVis ASCII visualization/dump: this function mimics mfem::GridFunction::Save() */
289: const PetscScalar *array;
290: PetscInt i, n, vdim, ordering = 1; /* mfem::FiniteElementSpace::Ordering::byVDIM */
291: PetscContainer glvis_container;
292: PetscViewerGLVisVecInfo glvis_vec_info;
293: PetscViewerGLVisInfo glvis_info;
295: /* mfem::FiniteElementSpace::Save() */
296: PetscCall(VecGetBlockSize(xin, &vdim));
297: PetscCall(PetscViewerASCIIPrintf(viewer, "FiniteElementSpace\n"));
298: PetscCall(PetscObjectQuery((PetscObject)xin, "_glvis_info_container", (PetscObject *)&glvis_container));
299: PetscCheck(glvis_container, PetscObjectComm((PetscObject)xin), PETSC_ERR_PLIB, "Missing GLVis container");
300: PetscCall(PetscContainerGetPointer(glvis_container, (void **)&glvis_vec_info));
301: PetscCall(PetscViewerASCIIPrintf(viewer, "%s\n", glvis_vec_info->fec_type));
302: PetscCall(PetscViewerASCIIPrintf(viewer, "VDim: %" PetscInt_FMT "\n", vdim));
303: PetscCall(PetscViewerASCIIPrintf(viewer, "Ordering: %" PetscInt_FMT "\n", ordering));
304: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
305: /* mfem::Vector::Print() */
306: PetscCall(PetscObjectQuery((PetscObject)viewer, "_glvis_info_container", (PetscObject *)&glvis_container));
307: PetscCheck(glvis_container, PetscObjectComm((PetscObject)viewer), PETSC_ERR_PLIB, "Missing GLVis container");
308: PetscCall(PetscContainerGetPointer(glvis_container, (void **)&glvis_info));
309: if (glvis_info->enabled) {
310: PetscCall(VecGetLocalSize(xin, &n));
311: PetscCall(VecGetArrayRead(xin, &array));
312: for (i = 0; i < n; i++) {
313: PetscCall(PetscViewerASCIIPrintf(viewer, glvis_info->fmt, (double)PetscRealPart(array[i])));
314: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
315: }
316: PetscCall(VecRestoreArrayRead(xin, &array));
317: }
318: } else if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
319: /* No info */
320: } else {
321: if (format != PETSC_VIEWER_ASCII_COMMON) PetscCall(PetscViewerASCIIPrintf(viewer, "Process [%d]\n", rank));
322: cnt = 0;
323: for (i = 0; i < xin->map->n; i++) {
324: if (format == PETSC_VIEWER_ASCII_INDEX) PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT ": ", cnt++));
325: #if defined(PETSC_USE_COMPLEX)
326: if (PetscImaginaryPart(xarray[i]) > 0.0) {
327: PetscCall(PetscViewerASCIIPrintf(viewer, "%g + %g i\n", (double)PetscRealPart(xarray[i]), (double)PetscImaginaryPart(xarray[i])));
328: } else if (PetscImaginaryPart(xarray[i]) < 0.0) {
329: PetscCall(PetscViewerASCIIPrintf(viewer, "%g - %g i\n", (double)PetscRealPart(xarray[i]), -(double)PetscImaginaryPart(xarray[i])));
330: } else {
331: PetscCall(PetscViewerASCIIPrintf(viewer, "%g\n", (double)PetscRealPart(xarray[i])));
332: }
333: #else
334: PetscCall(PetscViewerASCIIPrintf(viewer, "%g\n", (double)xarray[i]));
335: #endif
336: }
337: /* receive and print messages */
338: for (j = 1; j < size; j++) {
339: PetscCallMPI(MPI_Recv(values, (PetscMPIInt)len, MPIU_SCALAR, j, tag, PetscObjectComm((PetscObject)xin), &status));
340: PetscCallMPI(MPI_Get_count(&status, MPIU_SCALAR, &n));
341: if (format != PETSC_VIEWER_ASCII_COMMON) PetscCall(PetscViewerASCIIPrintf(viewer, "Process [%d]\n", j));
342: for (i = 0; i < n; i++) {
343: if (format == PETSC_VIEWER_ASCII_INDEX) PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT ": ", cnt++));
344: #if defined(PETSC_USE_COMPLEX)
345: if (PetscImaginaryPart(values[i]) > 0.0) {
346: PetscCall(PetscViewerASCIIPrintf(viewer, "%g + %g i\n", (double)PetscRealPart(values[i]), (double)PetscImaginaryPart(values[i])));
347: } else if (PetscImaginaryPart(values[i]) < 0.0) {
348: PetscCall(PetscViewerASCIIPrintf(viewer, "%g - %g i\n", (double)PetscRealPart(values[i]), -(double)PetscImaginaryPart(values[i])));
349: } else {
350: PetscCall(PetscViewerASCIIPrintf(viewer, "%g\n", (double)PetscRealPart(values[i])));
351: }
352: #else
353: PetscCall(PetscViewerASCIIPrintf(viewer, "%g\n", (double)values[i]));
354: #endif
355: }
356: }
357: }
358: PetscCall(PetscFree(values));
359: } else {
360: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
361: /* Rank 0 is not trying to receive anything, so don't send anything */
362: } else {
363: if (format == PETSC_VIEWER_ASCII_MATLAB || format == PETSC_VIEWER_ASCII_VTK_DEPRECATED || format == PETSC_VIEWER_ASCII_VTK_CELL_DEPRECATED) {
364: /* this may be a collective operation so make sure everyone calls it */
365: PetscCall(PetscObjectGetName((PetscObject)xin, &name));
366: }
367: PetscCallMPI(MPI_Send((void *)xarray, xin->map->n, MPIU_SCALAR, 0, tag, PetscObjectComm((PetscObject)xin)));
368: }
369: }
370: PetscCall(PetscViewerFlush(viewer));
371: PetscCall(VecRestoreArrayRead(xin, &xarray));
372: PetscFunctionReturn(PETSC_SUCCESS);
373: }
375: PetscErrorCode VecView_MPI_Binary(Vec xin, PetscViewer viewer)
376: {
377: return VecView_Binary(xin, viewer);
378: }
380: #include <petscdraw.h>
381: PetscErrorCode VecView_MPI_Draw_LG(Vec xin, PetscViewer viewer)
382: {
383: PetscDraw draw;
384: PetscBool isnull;
385: PetscDrawLG lg;
386: PetscMPIInt i, size, rank, n, N, *lens = NULL, *disp = NULL;
387: PetscReal *values, *xx = NULL, *yy = NULL;
388: const PetscScalar *xarray;
389: int colors[] = {PETSC_DRAW_RED};
391: PetscFunctionBegin;
392: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
393: PetscCall(PetscDrawIsNull(draw, &isnull));
394: if (isnull) PetscFunctionReturn(PETSC_SUCCESS);
395: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)xin), &rank));
396: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)xin), &size));
397: PetscCall(PetscMPIIntCast(xin->map->n, &n));
398: PetscCall(PetscMPIIntCast(xin->map->N, &N));
400: PetscCall(VecGetArrayRead(xin, &xarray));
401: #if defined(PETSC_USE_COMPLEX)
402: PetscCall(PetscMalloc1(n + 1, &values));
403: for (i = 0; i < n; i++) values[i] = PetscRealPart(xarray[i]);
404: #else
405: values = (PetscReal *)xarray;
406: #endif
407: if (rank == 0) {
408: PetscCall(PetscMalloc2(N, &xx, N, &yy));
409: for (i = 0; i < N; i++) xx[i] = (PetscReal)i;
410: PetscCall(PetscMalloc2(size, &lens, size, &disp));
411: for (i = 0; i < size; i++) lens[i] = (PetscMPIInt)xin->map->range[i + 1] - (PetscMPIInt)xin->map->range[i];
412: for (i = 0; i < size; i++) disp[i] = (PetscMPIInt)xin->map->range[i];
413: }
414: PetscCallMPI(MPI_Gatherv(values, n, MPIU_REAL, yy, lens, disp, MPIU_REAL, 0, PetscObjectComm((PetscObject)xin)));
415: PetscCall(PetscFree2(lens, disp));
416: #if defined(PETSC_USE_COMPLEX)
417: PetscCall(PetscFree(values));
418: #endif
419: PetscCall(VecRestoreArrayRead(xin, &xarray));
421: PetscCall(PetscViewerDrawGetDrawLG(viewer, 0, &lg));
422: PetscCall(PetscDrawLGReset(lg));
423: PetscCall(PetscDrawLGSetDimension(lg, 1));
424: PetscCall(PetscDrawLGSetColors(lg, colors));
425: if (rank == 0) {
426: PetscCall(PetscDrawLGAddPoints(lg, N, &xx, &yy));
427: PetscCall(PetscFree2(xx, yy));
428: }
429: PetscCall(PetscDrawLGDraw(lg));
430: PetscCall(PetscDrawLGSave(lg));
431: PetscFunctionReturn(PETSC_SUCCESS);
432: }
434: PetscErrorCode VecView_MPI_Draw(Vec xin, PetscViewer viewer)
435: {
436: PetscMPIInt rank, size, tag = ((PetscObject)viewer)->tag;
437: PetscInt i, start, end;
438: MPI_Status status;
439: PetscReal min, max, tmp = 0.0;
440: PetscDraw draw;
441: PetscBool isnull;
442: PetscDrawAxis axis;
443: const PetscScalar *xarray;
445: PetscFunctionBegin;
446: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
447: PetscCall(PetscDrawIsNull(draw, &isnull));
448: if (isnull) PetscFunctionReturn(PETSC_SUCCESS);
449: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)xin), &size));
450: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)xin), &rank));
452: PetscCall(VecMin(xin, NULL, &min));
453: PetscCall(VecMax(xin, NULL, &max));
454: if (min == max) {
455: min -= 1.e-5;
456: max += 1.e-5;
457: }
459: PetscCall(PetscDrawCheckResizedWindow(draw));
460: PetscCall(PetscDrawClear(draw));
462: PetscCall(PetscDrawAxisCreate(draw, &axis));
463: PetscCall(PetscDrawAxisSetLimits(axis, 0.0, (PetscReal)xin->map->N, min, max));
464: PetscCall(PetscDrawAxisDraw(axis));
465: PetscCall(PetscDrawAxisDestroy(&axis));
467: /* draw local part of vector */
468: PetscCall(VecGetArrayRead(xin, &xarray));
469: PetscCall(VecGetOwnershipRange(xin, &start, &end));
470: if (rank < size - 1) { /* send value to right */
471: PetscCallMPI(MPI_Send((void *)&xarray[xin->map->n - 1], 1, MPIU_REAL, rank + 1, tag, PetscObjectComm((PetscObject)xin)));
472: }
473: if (rank) { /* receive value from right */
474: PetscCallMPI(MPI_Recv(&tmp, 1, MPIU_REAL, rank - 1, tag, PetscObjectComm((PetscObject)xin), &status));
475: }
476: PetscDrawCollectiveBegin(draw);
477: if (rank) PetscCall(PetscDrawLine(draw, (PetscReal)start - 1, tmp, (PetscReal)start, PetscRealPart(xarray[0]), PETSC_DRAW_RED));
478: for (i = 1; i < xin->map->n; i++) PetscCall(PetscDrawLine(draw, (PetscReal)(i - 1 + start), PetscRealPart(xarray[i - 1]), (PetscReal)(i + start), PetscRealPart(xarray[i]), PETSC_DRAW_RED));
479: PetscDrawCollectiveEnd(draw);
480: PetscCall(VecRestoreArrayRead(xin, &xarray));
482: PetscCall(PetscDrawFlush(draw));
483: PetscCall(PetscDrawPause(draw));
484: PetscCall(PetscDrawSave(draw));
485: PetscFunctionReturn(PETSC_SUCCESS);
486: }
488: #if defined(PETSC_HAVE_MATLAB)
489: PetscErrorCode VecView_MPI_Matlab(Vec xin, PetscViewer viewer)
490: {
491: PetscMPIInt rank, size, *lens;
492: PetscInt i, N = xin->map->N;
493: const PetscScalar *xarray;
494: PetscScalar *xx;
496: PetscFunctionBegin;
497: PetscCall(VecGetArrayRead(xin, &xarray));
498: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)xin), &rank));
499: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)xin), &size));
500: if (rank == 0) {
501: PetscCall(PetscMalloc1(N, &xx));
502: PetscCall(PetscMalloc1(size, &lens));
503: for (i = 0; i < size; i++) lens[i] = xin->map->range[i + 1] - xin->map->range[i];
505: PetscCallMPI(MPI_Gatherv((void *)xarray, xin->map->n, MPIU_SCALAR, xx, lens, xin->map->range, MPIU_SCALAR, 0, PetscObjectComm((PetscObject)xin)));
506: PetscCall(PetscFree(lens));
508: PetscCall(PetscObjectName((PetscObject)xin));
509: PetscCall(PetscViewerMatlabPutArray(viewer, N, 1, xx, ((PetscObject)xin)->name));
511: PetscCall(PetscFree(xx));
512: } else {
513: PetscCallMPI(MPI_Gatherv((void *)xarray, xin->map->n, MPIU_SCALAR, 0, 0, 0, MPIU_SCALAR, 0, PetscObjectComm((PetscObject)xin)));
514: }
515: PetscCall(VecRestoreArrayRead(xin, &xarray));
516: PetscFunctionReturn(PETSC_SUCCESS);
517: }
518: #endif
520: #if defined(PETSC_HAVE_ADIOS)
521: #include <adios.h>
522: #include <adios_read.h>
523: #include <petsc/private/vieweradiosimpl.h>
524: #include <petsc/private/viewerimpl.h>
526: PetscErrorCode VecView_MPI_ADIOS(Vec xin, PetscViewer viewer)
527: {
528: PetscViewer_ADIOS *adios = (PetscViewer_ADIOS *)viewer->data;
529: const char *vecname;
530: int64_t id;
531: PetscInt n, N, rstart;
532: const PetscScalar *array;
533: char nglobalname[16], nlocalname[16], coffset[16];
535: PetscFunctionBegin;
536: PetscCall(PetscObjectGetName((PetscObject)xin, &vecname));
538: PetscCall(VecGetLocalSize(xin, &n));
539: PetscCall(VecGetSize(xin, &N));
540: PetscCall(VecGetOwnershipRange(xin, &rstart, NULL));
542: PetscCall(PetscSNPrintf(nlocalname, PETSC_STATIC_ARRAY_LENGTH(nlocalname), "%" PetscInt_FMT, n));
543: PetscCall(PetscSNPrintf(nglobalname, PETSC_STATIC_ARRAY_LENGTH(nglobalname), "%" PetscInt_FMT, N));
544: PetscCall(PetscSNPrintf(coffset, PETSC_STATIC_ARRAY_LENGTH(coffset), "%" PetscInt_FMT, rstart));
545: id = adios_define_var(Petsc_adios_group, vecname, "", adios_double, nlocalname, nglobalname, coffset);
546: PetscCall(VecGetArrayRead(xin, &array));
547: PetscCallExternal(adios_write_byid, adios->adios_handle, id, array);
548: PetscCall(VecRestoreArrayRead(xin, &array));
549: PetscFunctionReturn(PETSC_SUCCESS);
550: }
551: #endif
553: #if defined(PETSC_HAVE_HDF5)
554: PetscErrorCode VecView_MPI_HDF5(Vec xin, PetscViewer viewer)
555: {
556: PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5 *)viewer->data;
557: /* TODO: It looks like we can remove the H5Sclose(filespace) and H5Dget_space(dset_id). Why do we do this? */
558: hid_t filespace; /* file dataspace identifier */
559: hid_t chunkspace; /* chunk dataset property identifier */
560: hid_t dset_id; /* dataset identifier */
561: hid_t memspace; /* memory dataspace identifier */
562: hid_t file_id;
563: hid_t group;
564: hid_t memscalartype; /* scalar type for mem (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */
565: hid_t filescalartype; /* scalar type for file (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */
566: PetscInt bs = PetscAbs(xin->map->bs);
567: hsize_t dim;
568: hsize_t maxDims[4], dims[4], chunkDims[4], count[4], offset[4];
569: PetscBool timestepping, dim2, spoutput;
570: PetscInt timestep = PETSC_MIN_INT, low;
571: hsize_t chunksize;
572: const PetscScalar *x;
573: const char *vecname;
575: PetscFunctionBegin;
576: PetscCall(PetscViewerHDF5OpenGroup(viewer, NULL, &file_id, &group));
577: PetscCall(PetscViewerHDF5IsTimestepping(viewer, ×tepping));
578: if (timestepping) PetscCall(PetscViewerHDF5GetTimestep(viewer, ×tep));
579: PetscCall(PetscViewerHDF5GetBaseDimension2(viewer, &dim2));
580: PetscCall(PetscViewerHDF5GetSPOutput(viewer, &spoutput));
582: /* Create the dataspace for the dataset.
583: *
584: * dims - holds the current dimensions of the dataset
585: *
586: * maxDims - holds the maximum dimensions of the dataset (unlimited
587: * for the number of time steps with the current dimensions for the
588: * other dimensions; so only additional time steps can be added).
589: *
590: * chunkDims - holds the size of a single time step (required to
591: * permit extending dataset).
592: */
593: dim = 0;
594: chunksize = 1;
595: if (timestep >= 0) {
596: dims[dim] = timestep + 1;
597: maxDims[dim] = H5S_UNLIMITED;
598: chunkDims[dim] = 1;
599: ++dim;
600: }
601: PetscCall(PetscHDF5IntCast(xin->map->N / bs, dims + dim));
603: maxDims[dim] = dims[dim];
604: chunkDims[dim] = PetscMax(1, dims[dim]);
605: chunksize *= chunkDims[dim];
606: ++dim;
607: if (bs > 1 || dim2) {
608: dims[dim] = bs;
609: maxDims[dim] = dims[dim];
610: chunkDims[dim] = PetscMax(1, dims[dim]);
611: chunksize *= chunkDims[dim];
612: ++dim;
613: }
614: #if defined(PETSC_USE_COMPLEX)
615: dims[dim] = 2;
616: maxDims[dim] = dims[dim];
617: chunkDims[dim] = PetscMax(1, dims[dim]);
618: chunksize *= chunkDims[dim];
619: /* hdf5 chunks must be less than 4GB */
620: if (chunksize > PETSC_HDF5_MAX_CHUNKSIZE / 64) {
621: if (bs > 1 || dim2) {
622: if (chunkDims[dim - 2] > (hsize_t)PetscSqrtReal((PetscReal)(PETSC_HDF5_MAX_CHUNKSIZE / 128))) chunkDims[dim - 2] = (hsize_t)PetscSqrtReal((PetscReal)(PETSC_HDF5_MAX_CHUNKSIZE / 128));
623: if (chunkDims[dim - 1] > (hsize_t)PetscSqrtReal((PetscReal)(PETSC_HDF5_MAX_CHUNKSIZE / 128))) chunkDims[dim - 1] = (hsize_t)PetscSqrtReal((PetscReal)(PETSC_HDF5_MAX_CHUNKSIZE / 128));
624: } else {
625: chunkDims[dim - 1] = PETSC_HDF5_MAX_CHUNKSIZE / 128;
626: }
627: }
628: ++dim;
629: #else
630: /* hdf5 chunks must be less than 4GB */
631: if (chunksize > PETSC_HDF5_MAX_CHUNKSIZE / 64) {
632: if (bs > 1 || dim2) {
633: if (chunkDims[dim - 2] > (hsize_t)PetscSqrtReal((PetscReal)(PETSC_HDF5_MAX_CHUNKSIZE / 64))) chunkDims[dim - 2] = (hsize_t)PetscSqrtReal((PetscReal)(PETSC_HDF5_MAX_CHUNKSIZE / 64));
634: if (chunkDims[dim - 1] > (hsize_t)PetscSqrtReal((PetscReal)(PETSC_HDF5_MAX_CHUNKSIZE / 64))) chunkDims[dim - 1] = (hsize_t)PetscSqrtReal((PetscReal)(PETSC_HDF5_MAX_CHUNKSIZE / 64));
635: } else {
636: chunkDims[dim - 1] = PETSC_HDF5_MAX_CHUNKSIZE / 64;
637: }
638: }
639: #endif
641: PetscCallHDF5Return(filespace, H5Screate_simple, (dim, dims, maxDims));
643: #if defined(PETSC_USE_REAL_SINGLE)
644: memscalartype = H5T_NATIVE_FLOAT;
645: filescalartype = H5T_NATIVE_FLOAT;
646: #elif defined(PETSC_USE_REAL___FLOAT128)
647: #error "HDF5 output with 128 bit floats not supported."
648: #elif defined(PETSC_USE_REAL___FP16)
649: #error "HDF5 output with 16 bit floats not supported."
650: #else
651: memscalartype = H5T_NATIVE_DOUBLE;
652: if (spoutput == PETSC_TRUE) filescalartype = H5T_NATIVE_FLOAT;
653: else filescalartype = H5T_NATIVE_DOUBLE;
654: #endif
656: /* Create the dataset with default properties and close filespace */
657: PetscCall(PetscObjectGetName((PetscObject)xin, &vecname));
658: if (H5Lexists(group, vecname, H5P_DEFAULT) < 1) {
659: /* Create chunk */
660: PetscCallHDF5Return(chunkspace, H5Pcreate, (H5P_DATASET_CREATE));
661: PetscCallHDF5(H5Pset_chunk, (chunkspace, dim, chunkDims));
663: PetscCallHDF5Return(dset_id, H5Dcreate2, (group, vecname, filescalartype, filespace, H5P_DEFAULT, chunkspace, H5P_DEFAULT));
664: PetscCallHDF5(H5Pclose, (chunkspace));
665: } else {
666: PetscCallHDF5Return(dset_id, H5Dopen2, (group, vecname, H5P_DEFAULT));
667: PetscCallHDF5(H5Dset_extent, (dset_id, dims));
668: }
669: PetscCallHDF5(H5Sclose, (filespace));
671: /* Each process defines a dataset and writes it to the hyperslab in the file */
672: dim = 0;
673: if (timestep >= 0) {
674: count[dim] = 1;
675: ++dim;
676: }
677: PetscCall(PetscHDF5IntCast(xin->map->n / bs, count + dim));
678: ++dim;
679: if (bs > 1 || dim2) {
680: count[dim] = bs;
681: ++dim;
682: }
683: #if defined(PETSC_USE_COMPLEX)
684: count[dim] = 2;
685: ++dim;
686: #endif
687: if (xin->map->n > 0 || H5_VERSION_GE(1, 10, 0)) {
688: PetscCallHDF5Return(memspace, H5Screate_simple, (dim, count, NULL));
689: } else {
690: /* Can't create dataspace with zero for any dimension, so create null dataspace. */
691: PetscCallHDF5Return(memspace, H5Screate, (H5S_NULL));
692: }
694: /* Select hyperslab in the file */
695: PetscCall(VecGetOwnershipRange(xin, &low, NULL));
696: dim = 0;
697: if (timestep >= 0) {
698: offset[dim] = timestep;
699: ++dim;
700: }
701: PetscCall(PetscHDF5IntCast(low / bs, offset + dim));
702: ++dim;
703: if (bs > 1 || dim2) {
704: offset[dim] = 0;
705: ++dim;
706: }
707: #if defined(PETSC_USE_COMPLEX)
708: offset[dim] = 0;
709: ++dim;
710: #endif
711: if (xin->map->n > 0 || H5_VERSION_GE(1, 10, 0)) {
712: PetscCallHDF5Return(filespace, H5Dget_space, (dset_id));
713: PetscCallHDF5(H5Sselect_hyperslab, (filespace, H5S_SELECT_SET, offset, NULL, count, NULL));
714: } else {
715: /* Create null filespace to match null memspace. */
716: PetscCallHDF5Return(filespace, H5Screate, (H5S_NULL));
717: }
719: PetscCall(VecGetArrayRead(xin, &x));
720: PetscCallHDF5(H5Dwrite, (dset_id, memscalartype, memspace, filespace, hdf5->dxpl_id, x));
721: PetscCallHDF5(H5Fflush, (file_id, H5F_SCOPE_GLOBAL));
722: PetscCall(VecRestoreArrayRead(xin, &x));
724: /* Close/release resources */
725: PetscCallHDF5(H5Gclose, (group));
726: PetscCallHDF5(H5Sclose, (filespace));
727: PetscCallHDF5(H5Sclose, (memspace));
728: PetscCallHDF5(H5Dclose, (dset_id));
730: #if defined(PETSC_USE_COMPLEX)
731: {
732: PetscBool tru = PETSC_TRUE;
733: PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)xin, "complex", PETSC_BOOL, &tru));
734: }
735: #endif
736: if (timestepping) PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)xin, "timestepping", PETSC_BOOL, ×tepping));
737: PetscCall(PetscInfo(xin, "Wrote Vec object with name %s\n", vecname));
738: PetscFunctionReturn(PETSC_SUCCESS);
739: }
740: #endif
742: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec xin, PetscViewer viewer)
743: {
744: PetscBool iascii, isbinary, isdraw;
745: #if defined(PETSC_HAVE_MATHEMATICA)
746: PetscBool ismathematica;
747: #endif
748: #if defined(PETSC_HAVE_HDF5)
749: PetscBool ishdf5;
750: #endif
751: #if defined(PETSC_HAVE_MATLAB)
752: PetscBool ismatlab;
753: #endif
754: #if defined(PETSC_HAVE_ADIOS)
755: PetscBool isadios;
756: #endif
757: PetscBool isglvis;
759: PetscFunctionBegin;
760: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
761: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
762: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
763: #if defined(PETSC_HAVE_MATHEMATICA)
764: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERMATHEMATICA, &ismathematica));
765: #endif
766: #if defined(PETSC_HAVE_HDF5)
767: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
768: #endif
769: #if defined(PETSC_HAVE_MATLAB)
770: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERMATLAB, &ismatlab));
771: #endif
772: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERGLVIS, &isglvis));
773: #if defined(PETSC_HAVE_ADIOS)
774: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERADIOS, &isadios));
775: #endif
776: if (iascii) {
777: PetscCall(VecView_MPI_ASCII(xin, viewer));
778: } else if (isbinary) {
779: PetscCall(VecView_MPI_Binary(xin, viewer));
780: } else if (isdraw) {
781: PetscViewerFormat format;
782: PetscCall(PetscViewerGetFormat(viewer, &format));
783: if (format == PETSC_VIEWER_DRAW_LG) {
784: PetscCall(VecView_MPI_Draw_LG(xin, viewer));
785: } else {
786: PetscCall(VecView_MPI_Draw(xin, viewer));
787: }
788: #if defined(PETSC_HAVE_MATHEMATICA)
789: } else if (ismathematica) {
790: PetscCall(PetscViewerMathematicaPutVector(viewer, xin));
791: #endif
792: #if defined(PETSC_HAVE_HDF5)
793: } else if (ishdf5) {
794: PetscCall(VecView_MPI_HDF5(xin, viewer));
795: #endif
796: #if defined(PETSC_HAVE_ADIOS)
797: } else if (isadios) {
798: PetscCall(VecView_MPI_ADIOS(xin, viewer));
799: #endif
800: #if defined(PETSC_HAVE_MATLAB)
801: } else if (ismatlab) {
802: PetscCall(VecView_MPI_Matlab(xin, viewer));
803: #endif
804: } else if (isglvis) PetscCall(VecView_GLVis(xin, viewer));
805: PetscFunctionReturn(PETSC_SUCCESS);
806: }
808: PetscErrorCode VecGetSize_MPI(Vec xin, PetscInt *N)
809: {
810: PetscFunctionBegin;
811: *N = xin->map->N;
812: PetscFunctionReturn(PETSC_SUCCESS);
813: }
815: PetscErrorCode VecGetValues_MPI(Vec xin, PetscInt ni, const PetscInt ix[], PetscScalar y[])
816: {
817: const PetscScalar *xx;
818: const PetscInt start = xin->map->range[xin->stash.rank];
820: PetscFunctionBegin;
821: PetscCall(VecGetArrayRead(xin, &xx));
822: for (PetscInt i = 0; i < ni; i++) {
823: if (xin->stash.ignorenegidx && ix[i] < 0) continue;
824: const PetscInt tmp = ix[i] - start;
826: PetscCheck(tmp >= 0 && tmp < xin->map->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Can only get local values, trying %" PetscInt_FMT, ix[i]);
827: y[i] = xx[tmp];
828: }
829: PetscCall(VecRestoreArrayRead(xin, &xx));
830: PetscFunctionReturn(PETSC_SUCCESS);
831: }
833: PetscErrorCode VecSetValues_MPI(Vec xin, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode addv)
834: {
835: const PetscBool ignorenegidx = xin->stash.ignorenegidx;
836: const PetscBool donotstash = xin->stash.donotstash;
837: const PetscMPIInt rank = xin->stash.rank;
838: const PetscInt *owners = xin->map->range;
839: const PetscInt start = owners[rank], end = owners[rank + 1];
840: PetscScalar *xx;
842: PetscFunctionBegin;
843: if (PetscDefined(USE_DEBUG)) {
844: PetscCheck(xin->stash.insertmode != INSERT_VALUES || addv != ADD_VALUES, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "You have already inserted values; you cannot now add");
845: PetscCheck(xin->stash.insertmode != ADD_VALUES || addv != INSERT_VALUES, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "You have already added values; you cannot now insert");
846: }
847: PetscCall(VecGetArray(xin, &xx));
848: xin->stash.insertmode = addv;
849: for (PetscInt i = 0; i < ni; ++i) {
850: PetscInt row;
852: if (ignorenegidx && ix[i] < 0) continue;
853: PetscCheck(ix[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Out of range index value %" PetscInt_FMT " cannot be negative", ix[i]);
854: if ((row = ix[i]) >= start && row < end) {
855: if (addv == INSERT_VALUES) {
856: xx[row - start] = y[i];
857: } else {
858: xx[row - start] += y[i];
859: }
860: } else if (!donotstash) {
861: PetscCheck(ix[i] < xin->map->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Out of range index value %" PetscInt_FMT " maximum %" PetscInt_FMT, ix[i], xin->map->N);
862: PetscCall(VecStashValue_Private(&xin->stash, row, y[i]));
863: }
864: }
865: PetscCall(VecRestoreArray(xin, &xx));
866: PetscFunctionReturn(PETSC_SUCCESS);
867: }
869: PetscErrorCode VecSetValuesBlocked_MPI(Vec xin, PetscInt ni, const PetscInt ix[], const PetscScalar yin[], InsertMode addv)
870: {
871: PetscMPIInt rank = xin->stash.rank;
872: PetscInt *owners = xin->map->range, start = owners[rank];
873: PetscInt end = owners[rank + 1], i, row, bs = PetscAbs(xin->map->bs), j;
874: PetscScalar *xx, *y = (PetscScalar *)yin;
876: PetscFunctionBegin;
877: PetscCall(VecGetArray(xin, &xx));
878: if (PetscDefined(USE_DEBUG)) {
879: PetscCheck(xin->stash.insertmode != INSERT_VALUES || addv != ADD_VALUES, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "You have already inserted values; you cannot now add");
880: PetscCheck(xin->stash.insertmode != ADD_VALUES || addv != INSERT_VALUES, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "You have already added values; you cannot now insert");
881: }
882: xin->stash.insertmode = addv;
884: if (addv == INSERT_VALUES) {
885: for (i = 0; i < ni; i++) {
886: if ((row = bs * ix[i]) >= start && row < end) {
887: for (j = 0; j < bs; j++) xx[row - start + j] = y[j];
888: } else if (!xin->stash.donotstash) {
889: if (ix[i] < 0) {
890: y += bs;
891: continue;
892: }
893: PetscCheck(ix[i] < xin->map->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Out of range index value %" PetscInt_FMT " max %" PetscInt_FMT, ix[i], xin->map->N);
894: PetscCall(VecStashValuesBlocked_Private(&xin->bstash, ix[i], y));
895: }
896: y += bs;
897: }
898: } else {
899: for (i = 0; i < ni; i++) {
900: if ((row = bs * ix[i]) >= start && row < end) {
901: for (j = 0; j < bs; j++) xx[row - start + j] += y[j];
902: } else if (!xin->stash.donotstash) {
903: if (ix[i] < 0) {
904: y += bs;
905: continue;
906: }
907: PetscCheck(ix[i] <= xin->map->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Out of range index value %" PetscInt_FMT " max %" PetscInt_FMT, ix[i], xin->map->N);
908: PetscCall(VecStashValuesBlocked_Private(&xin->bstash, ix[i], y));
909: }
910: y += bs;
911: }
912: }
913: PetscCall(VecRestoreArray(xin, &xx));
914: PetscFunctionReturn(PETSC_SUCCESS);
915: }
917: /*
918: Since nsends or nreceives may be zero we add 1 in certain mallocs
919: to make sure we never malloc an empty one.
920: */
921: PetscErrorCode VecAssemblyBegin_MPI(Vec xin)
922: {
923: PetscInt *owners = xin->map->range, *bowners, i, bs, nstash, reallocs;
924: PetscMPIInt size;
925: InsertMode addv;
926: MPI_Comm comm;
928: PetscFunctionBegin;
929: PetscCall(PetscObjectGetComm((PetscObject)xin, &comm));
930: if (xin->stash.donotstash) PetscFunctionReturn(PETSC_SUCCESS);
932: PetscCall(MPIU_Allreduce((PetscEnum *)&xin->stash.insertmode, (PetscEnum *)&addv, 1, MPIU_ENUM, MPI_BOR, comm));
933: PetscCheck(addv != (ADD_VALUES | INSERT_VALUES), comm, PETSC_ERR_ARG_NOTSAMETYPE, "Some processors inserted values while others added");
934: xin->stash.insertmode = addv; /* in case this processor had no cache */
935: xin->bstash.insertmode = addv; /* Block stash implicitly tracks InsertMode of scalar stash */
937: PetscCall(VecGetBlockSize(xin, &bs));
938: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)xin), &size));
939: if (!xin->bstash.bowners && xin->map->bs != -1) {
940: PetscCall(PetscMalloc1(size + 1, &bowners));
941: for (i = 0; i < size + 1; i++) bowners[i] = owners[i] / bs;
942: xin->bstash.bowners = bowners;
943: } else bowners = xin->bstash.bowners;
945: PetscCall(VecStashScatterBegin_Private(&xin->stash, owners));
946: PetscCall(VecStashScatterBegin_Private(&xin->bstash, bowners));
947: PetscCall(VecStashGetInfo_Private(&xin->stash, &nstash, &reallocs));
948: PetscCall(PetscInfo(xin, "Stash has %" PetscInt_FMT " entries, uses %" PetscInt_FMT " mallocs.\n", nstash, reallocs));
949: PetscCall(VecStashGetInfo_Private(&xin->bstash, &nstash, &reallocs));
950: PetscCall(PetscInfo(xin, "Block-Stash has %" PetscInt_FMT " entries, uses %" PetscInt_FMT " mallocs.\n", nstash, reallocs));
951: PetscFunctionReturn(PETSC_SUCCESS);
952: }
954: PetscErrorCode VecAssemblyEnd_MPI(Vec vec)
955: {
956: PetscInt base, i, j, *row, flg, bs;
957: PetscMPIInt n;
958: PetscScalar *val, *vv, *array, *xarray;
960: PetscFunctionBegin;
961: if (!vec->stash.donotstash) {
962: PetscCall(VecGetArray(vec, &xarray));
963: PetscCall(VecGetBlockSize(vec, &bs));
964: base = vec->map->range[vec->stash.rank];
966: /* Process the stash */
967: while (1) {
968: PetscCall(VecStashScatterGetMesg_Private(&vec->stash, &n, &row, &val, &flg));
969: if (!flg) break;
970: if (vec->stash.insertmode == ADD_VALUES) {
971: for (i = 0; i < n; i++) xarray[row[i] - base] += val[i];
972: } else if (vec->stash.insertmode == INSERT_VALUES) {
973: for (i = 0; i < n; i++) xarray[row[i] - base] = val[i];
974: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_CORRUPT, "Insert mode is not set correctly; corrupted vector");
975: }
976: PetscCall(VecStashScatterEnd_Private(&vec->stash));
978: /* now process the block-stash */
979: while (1) {
980: PetscCall(VecStashScatterGetMesg_Private(&vec->bstash, &n, &row, &val, &flg));
981: if (!flg) break;
982: for (i = 0; i < n; i++) {
983: array = xarray + row[i] * bs - base;
984: vv = val + i * bs;
985: if (vec->stash.insertmode == ADD_VALUES) {
986: for (j = 0; j < bs; j++) array[j] += vv[j];
987: } else if (vec->stash.insertmode == INSERT_VALUES) {
988: for (j = 0; j < bs; j++) array[j] = vv[j];
989: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_CORRUPT, "Insert mode is not set correctly; corrupted vector");
990: }
991: }
992: PetscCall(VecStashScatterEnd_Private(&vec->bstash));
993: PetscCall(VecRestoreArray(vec, &xarray));
994: }
995: vec->stash.insertmode = NOT_SET_VALUES;
996: PetscFunctionReturn(PETSC_SUCCESS);
997: }
999: PetscErrorCode VecSetPreallocationCOO_MPI(Vec x, PetscCount coo_n, const PetscInt coo_i[])
1000: {
1001: PetscInt m, M, rstart, rend;
1002: Vec_MPI *vmpi = (Vec_MPI *)x->data;
1003: PetscCount k, p, q, rem; /* Loop variables over coo arrays */
1004: PetscMPIInt size;
1005: MPI_Comm comm;
1007: PetscFunctionBegin;
1008: PetscCall(PetscObjectGetComm((PetscObject)x, &comm));
1009: PetscCallMPI(MPI_Comm_size(comm, &size));
1010: PetscCall(VecResetPreallocationCOO_MPI(x));
1012: PetscCall(PetscLayoutSetUp(x->map));
1013: PetscCall(VecGetOwnershipRange(x, &rstart, &rend));
1014: PetscCall(VecGetLocalSize(x, &m));
1015: PetscCall(VecGetSize(x, &M));
1017: /* ---------------------------------------------------------------------------*/
1018: /* Sort COOs along with a permutation array, so that negative indices come */
1019: /* first, then local ones, then remote ones. */
1020: /* ---------------------------------------------------------------------------*/
1021: PetscCount n1 = coo_n, nneg, *perm;
1022: PetscInt *i1; /* Copy of input COOs along with a permutation array */
1023: PetscCall(PetscMalloc1(n1, &i1));
1024: PetscCall(PetscMalloc1(n1, &perm));
1025: PetscCall(PetscArraycpy(i1, coo_i, n1)); /* Make a copy since we'll modify it */
1026: for (k = 0; k < n1; k++) perm[k] = k;
1028: /* Manipulate i1[] so that entries with negative indices will have the smallest
1029: index, local entries will have greater but negative indices, and remote entries
1030: will have positive indices.
1031: */
1032: for (k = 0; k < n1; k++) {
1033: if (i1[k] < 0) {
1034: if (x->stash.ignorenegidx) i1[k] = PETSC_MIN_INT; /* e.g., -2^31, minimal to move them ahead */
1035: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Found a negative index in VecSetPreallocateCOO() but VEC_IGNORE_NEGATIVE_INDICES was not set");
1036: } else if (i1[k] >= rstart && i1[k] < rend) {
1037: i1[k] -= PETSC_MAX_INT; /* e.g., minus 2^31-1 to shift local rows to range of [-PETSC_MAX_INT, -1] */
1038: } else {
1039: PetscCheck(i1[k] < M, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Found index %" PetscInt_FMT " in VecSetPreallocateCOO() larger than the global size %" PetscInt_FMT "", i1[k], M);
1040: if (x->stash.donotstash) i1[k] = PETSC_MIN_INT; /* Ignore off-proc indices as if they were negative */
1041: }
1042: }
1044: /* Sort the indices, after that, [0,nneg) have ignored entries, [nneg,rem) have local entries and [rem,n1) have remote entries */
1045: PetscCall(PetscSortIntWithCountArray(n1, i1, perm));
1046: for (k = 0; k < n1; k++) {
1047: if (i1[k] > PETSC_MIN_INT) break;
1048: } /* Advance k to the first entry we need to take care of */
1049: nneg = k;
1050: PetscCall(PetscSortedIntUpperBound(i1, nneg, n1, rend - 1 - PETSC_MAX_INT, &rem)); /* rem is upper bound of the last local row */
1051: for (k = nneg; k < rem; k++) i1[k] += PETSC_MAX_INT; /* Revert indices of local entries */
1053: /* ---------------------------------------------------------------------------*/
1054: /* Build stuff for local entries */
1055: /* ---------------------------------------------------------------------------*/
1056: PetscCount tot1, *jmap1, *perm1;
1057: PetscCall(PetscCalloc1(m + 1, &jmap1));
1058: for (k = nneg; k < rem; k++) jmap1[i1[k] - rstart + 1]++; /* Count repeats of each local entry */
1059: for (k = 0; k < m; k++) jmap1[k + 1] += jmap1[k]; /* Transform jmap1[] to CSR-like data structure */
1060: tot1 = jmap1[m];
1061: PetscAssert(tot1 == rem - nneg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected errors in VecSetPreallocationCOO_MPI");
1062: PetscCall(PetscMalloc1(tot1, &perm1));
1063: PetscCall(PetscArraycpy(perm1, perm + nneg, tot1));
1065: /* ---------------------------------------------------------------------------*/
1066: /* Record the permutation array for filling the send buffer */
1067: /* ---------------------------------------------------------------------------*/
1068: PetscCount *Cperm;
1069: PetscCall(PetscMalloc1(n1 - rem, &Cperm));
1070: PetscCall(PetscArraycpy(Cperm, perm + rem, n1 - rem));
1071: PetscCall(PetscFree(perm));
1073: /* ---------------------------------------------------------------------------*/
1074: /* Send remote entries to their owner */
1075: /* ---------------------------------------------------------------------------*/
1076: /* Find which entries should be sent to which remote ranks*/
1077: PetscInt nsend = 0; /* Number of MPI ranks to send data to */
1078: PetscMPIInt *sendto; /* [nsend], storing remote ranks */
1079: PetscInt *nentries; /* [nsend], storing number of entries sent to remote ranks; Assume PetscInt is big enough for this count, and error if not */
1080: const PetscInt *ranges;
1081: PetscInt maxNsend = size >= 128 ? 128 : size; /* Assume max 128 neighbors; realloc when needed */
1083: PetscCall(PetscLayoutGetRanges(x->map, &ranges));
1084: PetscCall(PetscMalloc2(maxNsend, &sendto, maxNsend, &nentries));
1085: for (k = rem; k < n1;) {
1086: PetscMPIInt owner;
1087: PetscInt firstRow, lastRow;
1089: /* Locate a row range */
1090: firstRow = i1[k]; /* first row of this owner */
1091: PetscCall(PetscLayoutFindOwner(x->map, firstRow, &owner));
1092: lastRow = ranges[owner + 1] - 1; /* last row of this owner */
1094: /* Find the first index 'p' in [k,n) with i[p] belonging to next owner */
1095: PetscCall(PetscSortedIntUpperBound(i1, k, n1, lastRow, &p));
1097: /* All entries in [k,p) belong to this remote owner */
1098: if (nsend >= maxNsend) { /* Double the remote ranks arrays if not long enough */
1099: PetscMPIInt *sendto2;
1100: PetscInt *nentries2;
1101: PetscInt maxNsend2 = (maxNsend <= size / 2) ? maxNsend * 2 : size;
1103: PetscCall(PetscMalloc2(maxNsend2, &sendto2, maxNsend2, &nentries2));
1104: PetscCall(PetscArraycpy(sendto2, sendto, maxNsend));
1105: PetscCall(PetscArraycpy(nentries2, nentries2, maxNsend + 1));
1106: PetscCall(PetscFree2(sendto, nentries2));
1107: sendto = sendto2;
1108: nentries = nentries2;
1109: maxNsend = maxNsend2;
1110: }
1111: sendto[nsend] = owner;
1112: nentries[nsend] = p - k;
1113: PetscCall(PetscCountCast(p - k, &nentries[nsend]));
1114: nsend++;
1115: k = p;
1116: }
1118: /* Build 1st SF to know offsets on remote to send data */
1119: PetscSF sf1;
1120: PetscInt nroots = 1, nroots2 = 0;
1121: PetscInt nleaves = nsend, nleaves2 = 0;
1122: PetscInt *offsets;
1123: PetscSFNode *iremote;
1125: PetscCall(PetscSFCreate(comm, &sf1));
1126: PetscCall(PetscMalloc1(nsend, &iremote));
1127: PetscCall(PetscMalloc1(nsend, &offsets));
1128: for (k = 0; k < nsend; k++) {
1129: iremote[k].rank = sendto[k];
1130: iremote[k].index = 0;
1131: nleaves2 += nentries[k];
1132: PetscCheck(nleaves2 >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of SF leaves is too large for PetscInt");
1133: }
1134: PetscCall(PetscSFSetGraph(sf1, nroots, nleaves, NULL, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
1135: PetscCall(PetscSFFetchAndOpWithMemTypeBegin(sf1, MPIU_INT, PETSC_MEMTYPE_HOST, &nroots2 /*rootdata*/, PETSC_MEMTYPE_HOST, nentries /*leafdata*/, PETSC_MEMTYPE_HOST, offsets /*leafupdate*/, MPI_SUM));
1136: PetscCall(PetscSFFetchAndOpEnd(sf1, MPIU_INT, &nroots2, nentries, offsets, MPI_SUM)); /* Would nroots2 overflow, we check offsets[] below */
1137: PetscCall(PetscSFDestroy(&sf1));
1138: PetscAssert(nleaves2 == n1 - rem, PETSC_COMM_SELF, PETSC_ERR_PLIB, "nleaves2 %" PetscInt_FMT " != number of remote entries %" PetscCount_FMT "", nleaves2, n1 - rem);
1140: /* Build 2nd SF to send remote COOs to their owner */
1141: PetscSF sf2;
1142: nroots = nroots2;
1143: nleaves = nleaves2;
1144: PetscCall(PetscSFCreate(comm, &sf2));
1145: PetscCall(PetscSFSetFromOptions(sf2));
1146: PetscCall(PetscMalloc1(nleaves, &iremote));
1147: p = 0;
1148: for (k = 0; k < nsend; k++) {
1149: PetscCheck(offsets[k] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of SF roots is too large for PetscInt");
1150: for (q = 0; q < nentries[k]; q++, p++) {
1151: iremote[p].rank = sendto[k];
1152: iremote[p].index = offsets[k] + q;
1153: }
1154: }
1155: PetscCall(PetscSFSetGraph(sf2, nroots, nleaves, NULL, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
1157: /* Send the remote COOs to their owner */
1158: PetscInt n2 = nroots, *i2; /* Buffers for received COOs from other ranks, along with a permutation array */
1159: PetscCount *perm2;
1160: PetscCall(PetscMalloc1(n2, &i2));
1161: PetscCall(PetscMalloc1(n2, &perm2));
1162: PetscCall(PetscSFReduceWithMemTypeBegin(sf2, MPIU_INT, PETSC_MEMTYPE_HOST, i1 + rem, PETSC_MEMTYPE_HOST, i2, MPI_REPLACE));
1163: PetscCall(PetscSFReduceEnd(sf2, MPIU_INT, i1 + rem, i2, MPI_REPLACE));
1165: PetscCall(PetscFree(i1));
1166: PetscCall(PetscFree(offsets));
1167: PetscCall(PetscFree2(sendto, nentries));
1169: /* ---------------------------------------------------------------*/
1170: /* Sort received COOs along with a permutation array */
1171: /* ---------------------------------------------------------------*/
1172: PetscCount *imap2;
1173: PetscCount *jmap2, nnz2;
1174: PetscScalar *sendbuf, *recvbuf;
1175: PetscInt old;
1176: PetscCount sendlen = n1 - rem, recvlen = n2;
1178: for (k = 0; k < n2; k++) perm2[k] = k;
1179: PetscCall(PetscSortIntWithCountArray(n2, i2, perm2));
1181: /* nnz2 will be # of unique entries in the recvbuf */
1182: nnz2 = n2;
1183: for (k = 1; k < n2; k++) {
1184: if (i2[k] == i2[k - 1]) nnz2--;
1185: }
1187: /* Build imap2[] and jmap2[] for each unique entry */
1188: PetscCall(PetscMalloc4(nnz2, &imap2, nnz2 + 1, &jmap2, sendlen, &sendbuf, recvlen, &recvbuf));
1189: p = -1;
1190: old = -1;
1191: jmap2[0] = 0;
1192: jmap2++;
1193: for (k = 0; k < n2; k++) {
1194: if (i2[k] != old) { /* Meet a new entry */
1195: p++;
1196: imap2[p] = i2[k] - rstart;
1197: jmap2[p] = 1;
1198: old = i2[k];
1199: } else {
1200: jmap2[p]++;
1201: }
1202: }
1203: jmap2--;
1204: for (k = 0; k < nnz2; k++) jmap2[k + 1] += jmap2[k];
1206: PetscCall(PetscFree(i2));
1208: vmpi->coo_n = coo_n;
1209: vmpi->tot1 = tot1;
1210: vmpi->jmap1 = jmap1;
1211: vmpi->perm1 = perm1;
1212: vmpi->nnz2 = nnz2;
1213: vmpi->imap2 = imap2;
1214: vmpi->jmap2 = jmap2;
1215: vmpi->perm2 = perm2;
1217: vmpi->Cperm = Cperm;
1218: vmpi->sendbuf = sendbuf;
1219: vmpi->recvbuf = recvbuf;
1220: vmpi->sendlen = sendlen;
1221: vmpi->recvlen = recvlen;
1222: vmpi->coo_sf = sf2;
1223: PetscFunctionReturn(PETSC_SUCCESS);
1224: }
1226: PetscErrorCode VecSetValuesCOO_MPI(Vec x, const PetscScalar v[], InsertMode imode)
1227: {
1228: Vec_MPI *vmpi = (Vec_MPI *)x->data;
1229: PetscInt m;
1230: PetscScalar *a, *sendbuf = vmpi->sendbuf, *recvbuf = vmpi->recvbuf;
1231: const PetscCount *jmap1 = vmpi->jmap1;
1232: const PetscCount *perm1 = vmpi->perm1;
1233: const PetscCount *imap2 = vmpi->imap2;
1234: const PetscCount *jmap2 = vmpi->jmap2;
1235: const PetscCount *perm2 = vmpi->perm2;
1236: const PetscCount *Cperm = vmpi->Cperm;
1237: const PetscCount nnz2 = vmpi->nnz2;
1239: PetscFunctionBegin;
1240: PetscCall(VecGetLocalSize(x, &m));
1241: PetscCall(VecGetArray(x, &a));
1243: /* Pack entries to be sent to remote */
1244: for (PetscInt i = 0; i < vmpi->sendlen; i++) sendbuf[i] = v[Cperm[i]];
1246: /* Send remote entries to their owner and overlap the communication with local computation */
1247: PetscCall(PetscSFReduceWithMemTypeBegin(vmpi->coo_sf, MPIU_SCALAR, PETSC_MEMTYPE_HOST, sendbuf, PETSC_MEMTYPE_HOST, recvbuf, MPI_REPLACE));
1248: /* Add local entries to A and B */
1249: for (PetscInt i = 0; i < m; i++) { /* All entries in a[] are either zero'ed or added with a value (i.e., initialized) */
1250: PetscScalar sum = 0.0; /* Do partial summation first to improve numerical stability */
1251: for (PetscCount k = jmap1[i]; k < jmap1[i + 1]; k++) sum += v[perm1[k]];
1252: a[i] = (imode == INSERT_VALUES ? 0.0 : a[i]) + sum;
1253: }
1254: PetscCall(PetscSFReduceEnd(vmpi->coo_sf, MPIU_SCALAR, sendbuf, recvbuf, MPI_REPLACE));
1256: /* Add received remote entries to A and B */
1257: for (PetscInt i = 0; i < nnz2; i++) {
1258: for (PetscCount k = jmap2[i]; k < jmap2[i + 1]; k++) a[imap2[i]] += recvbuf[perm2[k]];
1259: }
1261: PetscCall(VecRestoreArray(x, &a));
1262: PetscFunctionReturn(PETSC_SUCCESS);
1263: }