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