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:   if (vmpi) {
 16:     PetscFree(vmpi->jmap1);
 17:     PetscFree(vmpi->perm1);
 18:     PetscFree(vmpi->Cperm);
 19:     PetscFree4(vmpi->imap2, vmpi->jmap2, vmpi->sendbuf, vmpi->recvbuf);
 20:     PetscFree(vmpi->perm2);
 21:     PetscSFDestroy(&vmpi->coo_sf);
 22:   }
 23:   return 0;
 24: }

 26: PetscErrorCode VecDestroy_MPI(Vec v)
 27: {
 28:   Vec_MPI *x = (Vec_MPI *)v->data;

 30: #if defined(PETSC_USE_LOG)
 31:   PetscLogObjectState((PetscObject)v, "Length=%" PetscInt_FMT, v->map->N);
 32: #endif
 33:   if (!x) return 0;
 34:   PetscFree(x->array_allocated);

 36:   /* Destroy local representation of vector if it exists */
 37:   if (x->localrep) {
 38:     VecDestroy(&x->localrep);
 39:     VecScatterDestroy(&x->localupdate);
 40:   }
 41:   VecAssemblyReset_MPI(v);

 43:   /* Destroy the stashes: note the order - so that the tags are freed properly */
 44:   VecStashDestroy_Private(&v->bstash);
 45:   VecStashDestroy_Private(&v->stash);

 47:   VecResetPreallocationCOO_MPI(v);
 48:   PetscObjectComposeFunction((PetscObject)v, "PetscMatlabEnginePut_C", NULL);
 49:   PetscObjectComposeFunction((PetscObject)v, "PetscMatlabEngineGet_C", NULL);
 50:   PetscFree(v->data);
 51:   return 0;
 52: }

 54: PetscErrorCode VecView_MPI_ASCII(Vec xin, PetscViewer viewer)
 55: {
 56:   PetscInt           i, work = xin->map->n, cnt, len, nLen;
 57:   PetscMPIInt        j, n = 0, size, rank, tag = ((PetscObject)viewer)->tag;
 58:   MPI_Status         status;
 59:   PetscScalar       *values;
 60:   const PetscScalar *xarray;
 61:   const char        *name;
 62:   PetscViewerFormat  format;

 64:   MPI_Comm_size(PetscObjectComm((PetscObject)xin), &size);
 65:   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:     PetscViewerASCIIPrintf(viewer, "  Load Balance - Local vector size Min %" PetscInt_FMT "  avg %" PetscInt_FMT "  max %" PetscInt_FMT "\n", nmin, navg, nmax);
 74:     return 0;
 75:   }

 77:   VecGetArrayRead(xin, &xarray);
 78:   /* determine maximum message to arrive */
 79:   MPI_Comm_rank(PetscObjectComm((PetscObject)xin), &rank);
 80:   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:     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:       PetscObjectGetName((PetscObject)xin, &name);
 90:       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:           PetscViewerASCIIPrintf(viewer, "%18.16e + %18.16ei\n", (double)PetscRealPart(xarray[i]), (double)PetscImaginaryPart(xarray[i]));
 95:         } else if (PetscImaginaryPart(xarray[i]) < 0.0) {
 96:           PetscViewerASCIIPrintf(viewer, "%18.16e - %18.16ei\n", (double)PetscRealPart(xarray[i]), -(double)PetscImaginaryPart(xarray[i]));
 97:         } else {
 98:           PetscViewerASCIIPrintf(viewer, "%18.16e\n", (double)PetscRealPart(xarray[i]));
 99:         }
100: #else
101:         PetscViewerASCIIPrintf(viewer, "%18.16e\n", (double)xarray[i]);
102: #endif
103:       }
104:       /* receive and print messages */
105:       for (j = 1; j < size; j++) {
106:         MPI_Recv(values, (PetscMPIInt)len, MPIU_SCALAR, j, tag, PetscObjectComm((PetscObject)xin), &status);
107:         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:             PetscViewerASCIIPrintf(viewer, "%18.16e + %18.16ei\n", (double)PetscRealPart(values[i]), (double)PetscImaginaryPart(values[i]));
112:           } else if (PetscImaginaryPart(values[i]) < 0.0) {
113:             PetscViewerASCIIPrintf(viewer, "%18.16e - %18.16ei\n", (double)PetscRealPart(values[i]), -(double)PetscImaginaryPart(values[i]));
114:           } else {
115:             PetscViewerASCIIPrintf(viewer, "%18.16e\n", (double)PetscRealPart(values[i]));
116:           }
117: #else
118:           PetscViewerASCIIPrintf(viewer, "%18.16e\n", (double)values[i]);
119: #endif
120:         }
121:       }
122:       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:         PetscViewerASCIIPrintf(viewer, "%18.16e %18.16e\n", (double)PetscRealPart(xarray[i]), (double)PetscImaginaryPart(xarray[i]));
128: #else
129:         PetscViewerASCIIPrintf(viewer, "%18.16e\n", (double)xarray[i]);
130: #endif
131:       }
132:       /* receive and print messages */
133:       for (j = 1; j < size; j++) {
134:         MPI_Recv(values, (PetscMPIInt)len, MPIU_SCALAR, j, tag, PetscObjectComm((PetscObject)xin), &status);
135:         MPI_Get_count(&status, MPIU_SCALAR, &n);
136:         for (i = 0; i < n; i++) {
137: #if defined(PETSC_USE_COMPLEX)
138:           PetscViewerASCIIPrintf(viewer, "%18.16e %18.16e\n", (double)PetscRealPart(values[i]), (double)PetscImaginaryPart(values[i]));
139: #else
140:           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) PetscObjectComposedDataRegister(&stateId);
159:       PetscObjectComposedDataGetInt((PetscObject)viewer, stateId, outputState, hasState);
160:       if (!hasState) outputState = 0;

162:       PetscObjectGetName((PetscObject)xin, &name);
163:       VecGetLocalSize(xin, &nLen);
164:       PetscMPIIntCast(nLen, &n);
165:       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;

177:         if (doOutput) 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 {
189:           if (outputState == 4) doOutput = 0;
190:         }

192:         if (doOutput) PetscViewerASCIIPrintf(viewer, "CELL_DATA %" PetscInt_FMT "\n", xin->map->N / bs);
193:       }
194:       PetscObjectComposedDataSetInt((PetscObject)viewer, stateId, outputState);
195:       if (name) {
196:         if (bs == 3) {
197:           PetscViewerASCIIPrintf(viewer, "VECTORS %s double\n", name);
198:         } else {
199:           PetscViewerASCIIPrintf(viewer, "SCALARS %s double %" PetscInt_FMT "\n", name, bs);
200:         }
201:       } else {
202:         PetscViewerASCIIPrintf(viewer, "SCALARS scalars double %" PetscInt_FMT "\n", bs);
203:       }
204:       if (bs != 3) 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) PetscViewerASCIIPrintf(viewer, " ");
208:           PetscViewerASCIIPrintf(viewer, "%g", (double)PetscRealPart(xarray[i * bs + b]));
209:         }
210:         PetscViewerASCIIPrintf(viewer, "\n");
211:       }
212:       for (j = 1; j < size; j++) {
213:         MPI_Recv(values, (PetscMPIInt)len, MPIU_SCALAR, j, tag, PetscObjectComm((PetscObject)xin), &status);
214:         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) PetscViewerASCIIPrintf(viewer, " ");
218:             PetscViewerASCIIPrintf(viewer, "%g", (double)PetscRealPart(values[i * bs + b]));
219:           }
220:           PetscViewerASCIIPrintf(viewer, "\n");
221:         }
222:       }
223:     } else if (format == PETSC_VIEWER_ASCII_VTK_COORDS_DEPRECATED) {
224:       PetscInt bs, b;

226:       VecGetLocalSize(xin, &nLen);
227:       PetscMPIIntCast(nLen, &n);
228:       VecGetBlockSize(xin, &bs);

231:       for (i = 0; i < n / bs; i++) {
232:         for (b = 0; b < bs; b++) {
233:           if (b > 0) PetscViewerASCIIPrintf(viewer, " ");
234:           PetscViewerASCIIPrintf(viewer, "%g", (double)PetscRealPart(xarray[i * bs + b]));
235:         }
236:         for (b = bs; b < 3; b++) PetscViewerASCIIPrintf(viewer, " 0.0");
237:         PetscViewerASCIIPrintf(viewer, "\n");
238:       }
239:       for (j = 1; j < size; j++) {
240:         MPI_Recv(values, (PetscMPIInt)len, MPIU_SCALAR, j, tag, PetscObjectComm((PetscObject)xin), &status);
241:         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) PetscViewerASCIIPrintf(viewer, " ");
245:             PetscViewerASCIIPrintf(viewer, "%g", (double)PetscRealPart(values[i * bs + b]));
246:           }
247:           for (b = bs; b < 3; b++) PetscViewerASCIIPrintf(viewer, " 0.0");
248:           PetscViewerASCIIPrintf(viewer, "\n");
249:         }
250:       }
251:     } else if (format == PETSC_VIEWER_ASCII_PCICE) {
252:       PetscInt bs, b, vertexCount = 1;

254:       VecGetLocalSize(xin, &nLen);
255:       PetscMPIIntCast(nLen, &n);
256:       VecGetBlockSize(xin, &bs);

259:       PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT "\n", xin->map->N / bs);
260:       for (i = 0; i < n / bs; i++) {
261:         PetscViewerASCIIPrintf(viewer, "%7" PetscInt_FMT "   ", vertexCount++);
262:         for (b = 0; b < bs; b++) {
263:           if (b > 0) PetscViewerASCIIPrintf(viewer, " ");
264: #if !defined(PETSC_USE_COMPLEX)
265:           PetscViewerASCIIPrintf(viewer, "% 12.5E", (double)xarray[i * bs + b]);
266: #endif
267:         }
268:         PetscViewerASCIIPrintf(viewer, "\n");
269:       }
270:       for (j = 1; j < size; j++) {
271:         MPI_Recv(values, (PetscMPIInt)len, MPIU_SCALAR, j, tag, PetscObjectComm((PetscObject)xin), &status);
272:         MPI_Get_count(&status, MPIU_SCALAR, &n);
273:         for (i = 0; i < n / bs; i++) {
274:           PetscViewerASCIIPrintf(viewer, "%7" PetscInt_FMT "   ", vertexCount++);
275:           for (b = 0; b < bs; b++) {
276:             if (b > 0) PetscViewerASCIIPrintf(viewer, " ");
277: #if !defined(PETSC_USE_COMPLEX)
278:             PetscViewerASCIIPrintf(viewer, "% 12.5E", (double)values[i * bs + b]);
279: #endif
280:           }
281:           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:       VecGetBlockSize(xin, &vdim);
294:       PetscViewerASCIIPrintf(viewer, "FiniteElementSpace\n");
295:       PetscObjectQuery((PetscObject)xin, "_glvis_info_container", (PetscObject *)&glvis_container);
297:       PetscContainerGetPointer(glvis_container, (void **)&glvis_vec_info);
298:       PetscViewerASCIIPrintf(viewer, "%s\n", glvis_vec_info->fec_type);
299:       PetscViewerASCIIPrintf(viewer, "VDim: %" PetscInt_FMT "\n", vdim);
300:       PetscViewerASCIIPrintf(viewer, "Ordering: %" PetscInt_FMT "\n", ordering);
301:       PetscViewerASCIIPrintf(viewer, "\n");
302:       /* mfem::Vector::Print() */
303:       PetscObjectQuery((PetscObject)viewer, "_glvis_info_container", (PetscObject *)&glvis_container);
305:       PetscContainerGetPointer(glvis_container, (void **)&glvis_info);
306:       if (glvis_info->enabled) {
307:         VecGetLocalSize(xin, &n);
308:         VecGetArrayRead(xin, &array);
309:         for (i = 0; i < n; i++) {
310:           PetscViewerASCIIPrintf(viewer, glvis_info->fmt, (double)PetscRealPart(array[i]));
311:           PetscViewerASCIIPrintf(viewer, "\n");
312:         }
313:         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) 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) PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT ": ", cnt++);
322: #if defined(PETSC_USE_COMPLEX)
323:         if (PetscImaginaryPart(xarray[i]) > 0.0) {
324:           PetscViewerASCIIPrintf(viewer, "%g + %g i\n", (double)PetscRealPart(xarray[i]), (double)PetscImaginaryPart(xarray[i]));
325:         } else if (PetscImaginaryPart(xarray[i]) < 0.0) {
326:           PetscViewerASCIIPrintf(viewer, "%g - %g i\n", (double)PetscRealPart(xarray[i]), -(double)PetscImaginaryPart(xarray[i]));
327:         } else {
328:           PetscViewerASCIIPrintf(viewer, "%g\n", (double)PetscRealPart(xarray[i]));
329:         }
330: #else
331:         PetscViewerASCIIPrintf(viewer, "%g\n", (double)xarray[i]);
332: #endif
333:       }
334:       /* receive and print messages */
335:       for (j = 1; j < size; j++) {
336:         MPI_Recv(values, (PetscMPIInt)len, MPIU_SCALAR, j, tag, PetscObjectComm((PetscObject)xin), &status);
337:         MPI_Get_count(&status, MPIU_SCALAR, &n);
338:         if (format != PETSC_VIEWER_ASCII_COMMON) PetscViewerASCIIPrintf(viewer, "Process [%d]\n", j);
339:         for (i = 0; i < n; i++) {
340:           if (format == PETSC_VIEWER_ASCII_INDEX) PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT ": ", cnt++);
341: #if defined(PETSC_USE_COMPLEX)
342:           if (PetscImaginaryPart(values[i]) > 0.0) {
343:             PetscViewerASCIIPrintf(viewer, "%g + %g i\n", (double)PetscRealPart(values[i]), (double)PetscImaginaryPart(values[i]));
344:           } else if (PetscImaginaryPart(values[i]) < 0.0) {
345:             PetscViewerASCIIPrintf(viewer, "%g - %g i\n", (double)PetscRealPart(values[i]), -(double)PetscImaginaryPart(values[i]));
346:           } else {
347:             PetscViewerASCIIPrintf(viewer, "%g\n", (double)PetscRealPart(values[i]));
348:           }
349: #else
350:           PetscViewerASCIIPrintf(viewer, "%g\n", (double)values[i]);
351: #endif
352:         }
353:       }
354:     }
355:     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:         PetscObjectGetName((PetscObject)xin, &name);
363:       }
364:       MPI_Send((void *)xarray, xin->map->n, MPIU_SCALAR, 0, tag, PetscObjectComm((PetscObject)xin));
365:     }
366:   }
367:   PetscViewerFlush(viewer);
368:   VecRestoreArrayRead(xin, &xarray);
369:   return 0;
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:   PetscViewerDrawGetDraw(viewer, 0, &draw);
389:   PetscDrawIsNull(draw, &isnull);
390:   if (isnull) return 0;
391:   MPI_Comm_rank(PetscObjectComm((PetscObject)xin), &rank);
392:   MPI_Comm_size(PetscObjectComm((PetscObject)xin), &size);
393:   PetscMPIIntCast(xin->map->n, &n);
394:   PetscMPIIntCast(xin->map->N, &N);

396:   VecGetArrayRead(xin, &xarray);
397: #if defined(PETSC_USE_COMPLEX)
398:   PetscMalloc1(n + 1, &values);
399:   for (i = 0; i < n; i++) values[i] = PetscRealPart(xarray[i]);
400: #else
401:   values = (PetscReal *)xarray;
402: #endif
403:   if (rank == 0) {
404:     PetscMalloc2(N, &xx, N, &yy);
405:     for (i = 0; i < N; i++) xx[i] = (PetscReal)i;
406:     PetscMalloc2(size, &lens, size, &disp);
407:     for (i = 0; i < size; i++) lens[i] = (PetscMPIInt)xin->map->range[i + 1] - (PetscMPIInt)xin->map->range[i];
408:     for (i = 0; i < size; i++) disp[i] = (PetscMPIInt)xin->map->range[i];
409:   }
410:   MPI_Gatherv(values, n, MPIU_REAL, yy, lens, disp, MPIU_REAL, 0, PetscObjectComm((PetscObject)xin));
411:   PetscFree2(lens, disp);
412: #if defined(PETSC_USE_COMPLEX)
413:   PetscFree(values);
414: #endif
415:   VecRestoreArrayRead(xin, &xarray);

417:   PetscViewerDrawGetDrawLG(viewer, 0, &lg);
418:   PetscDrawLGReset(lg);
419:   PetscDrawLGSetDimension(lg, 1);
420:   PetscDrawLGSetColors(lg, colors);
421:   if (rank == 0) {
422:     PetscDrawLGAddPoints(lg, N, &xx, &yy);
423:     PetscFree2(xx, yy);
424:   }
425:   PetscDrawLGDraw(lg);
426:   PetscDrawLGSave(lg);
427:   return 0;
428: }

430: PetscErrorCode VecView_MPI_Draw(Vec xin, PetscViewer viewer)
431: {
432:   PetscMPIInt        rank, size, tag = ((PetscObject)viewer)->tag;
433:   PetscInt           i, start, end;
434:   MPI_Status         status;
435:   PetscReal          min, max, tmp = 0.0;
436:   PetscDraw          draw;
437:   PetscBool          isnull;
438:   PetscDrawAxis      axis;
439:   const PetscScalar *xarray;

441:   PetscViewerDrawGetDraw(viewer, 0, &draw);
442:   PetscDrawIsNull(draw, &isnull);
443:   if (isnull) return 0;
444:   MPI_Comm_size(PetscObjectComm((PetscObject)xin), &size);
445:   MPI_Comm_rank(PetscObjectComm((PetscObject)xin), &rank);

447:   VecMin(xin, NULL, &min);
448:   VecMax(xin, NULL, &max);
449:   if (min == max) {
450:     min -= 1.e-5;
451:     max += 1.e-5;
452:   }

454:   PetscDrawCheckResizedWindow(draw);
455:   PetscDrawClear(draw);

457:   PetscDrawAxisCreate(draw, &axis);
458:   PetscDrawAxisSetLimits(axis, 0.0, (PetscReal)xin->map->N, min, max);
459:   PetscDrawAxisDraw(axis);
460:   PetscDrawAxisDestroy(&axis);

462:   /* draw local part of vector */
463:   VecGetArrayRead(xin, &xarray);
464:   VecGetOwnershipRange(xin, &start, &end);
465:   if (rank < size - 1) { /* send value to right */
466:     MPI_Send((void *)&xarray[xin->map->n - 1], 1, MPIU_REAL, rank + 1, tag, PetscObjectComm((PetscObject)xin));
467:   }
468:   if (rank) { /* receive value from right */
469:     MPI_Recv(&tmp, 1, MPIU_REAL, rank - 1, tag, PetscObjectComm((PetscObject)xin), &status);
470:   }
471:   PetscDrawCollectiveBegin(draw);
472:   if (rank) PetscDrawLine(draw, (PetscReal)start - 1, tmp, (PetscReal)start, PetscRealPart(xarray[0]), PETSC_DRAW_RED);
473:   for (i = 1; i < xin->map->n; i++) PetscDrawLine(draw, (PetscReal)(i - 1 + start), PetscRealPart(xarray[i - 1]), (PetscReal)(i + start), PetscRealPart(xarray[i]), PETSC_DRAW_RED);
474:   PetscDrawCollectiveEnd(draw);
475:   VecRestoreArrayRead(xin, &xarray);

477:   PetscDrawFlush(draw);
478:   PetscDrawPause(draw);
479:   PetscDrawSave(draw);
480:   return 0;
481: }

483: #if defined(PETSC_HAVE_MATLAB)
484: PetscErrorCode VecView_MPI_Matlab(Vec xin, PetscViewer viewer)
485: {
486:   PetscMPIInt        rank, size, *lens;
487:   PetscInt           i, N = xin->map->N;
488:   const PetscScalar *xarray;
489:   PetscScalar       *xx;

491:   VecGetArrayRead(xin, &xarray);
492:   MPI_Comm_rank(PetscObjectComm((PetscObject)xin), &rank);
493:   MPI_Comm_size(PetscObjectComm((PetscObject)xin), &size);
494:   if (rank == 0) {
495:     PetscMalloc1(N, &xx);
496:     PetscMalloc1(size, &lens);
497:     for (i = 0; i < size; i++) lens[i] = xin->map->range[i + 1] - xin->map->range[i];

499:     MPI_Gatherv((void *)xarray, xin->map->n, MPIU_SCALAR, xx, lens, xin->map->range, MPIU_SCALAR, 0, PetscObjectComm((PetscObject)xin));
500:     PetscFree(lens);

502:     PetscObjectName((PetscObject)xin);
503:     PetscViewerMatlabPutArray(viewer, N, 1, xx, ((PetscObject)xin)->name);

505:     PetscFree(xx);
506:   } else {
507:     MPI_Gatherv((void *)xarray, xin->map->n, MPIU_SCALAR, 0, 0, 0, MPIU_SCALAR, 0, PetscObjectComm((PetscObject)xin));
508:   }
509:   VecRestoreArrayRead(xin, &xarray);
510:   return 0;
511: }
512: #endif

514: #if defined(PETSC_HAVE_ADIOS)
515:   #include <adios.h>
516:   #include <adios_read.h>
517: #include <petsc/private/vieweradiosimpl.h>
518: #include <petsc/private/viewerimpl.h>

520: PetscErrorCode VecView_MPI_ADIOS(Vec xin, PetscViewer viewer)
521: {
522:   PetscViewer_ADIOS *adios = (PetscViewer_ADIOS *)viewer->data;
523:   const char        *vecname;
524:   int64_t            id;
525:   PetscInt           n, N, rstart;
526:   const PetscScalar *array;
527:   char               nglobalname[16], nlocalname[16], coffset[16];

529:   PetscObjectGetName((PetscObject)xin, &vecname);

531:   VecGetLocalSize(xin, &n);
532:   VecGetSize(xin, &N);
533:   VecGetOwnershipRange(xin, &rstart, NULL);

535:   sprintf(nlocalname, "%d", (int)n);
536:   sprintf(nglobalname, "%d", (int)N);
537:   sprintf(coffset, "%d", (int)rstart);
538:   id = adios_define_var(Petsc_adios_group, vecname, "", adios_double, nlocalname, nglobalname, coffset);
539:   VecGetArrayRead(xin, &array);
540:   adios_write_byid(adios->adios_handle, id, array);
541:   VecRestoreArrayRead(xin, &array);
542:   return 0;
543: }
544: #endif

546: #if defined(PETSC_HAVE_HDF5)
547: PetscErrorCode VecView_MPI_HDF5(Vec xin, PetscViewer viewer)
548: {
549:   PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5 *)viewer->data;
550:   /* TODO: It looks like we can remove the H5Sclose(filespace) and H5Dget_space(dset_id). Why do we do this? */
551:   hid_t              filespace;  /* file dataspace identifier */
552:   hid_t              chunkspace; /* chunk dataset property identifier */
553:   hid_t              dset_id;    /* dataset identifier */
554:   hid_t              memspace;   /* memory dataspace identifier */
555:   hid_t              file_id;
556:   hid_t              group;
557:   hid_t              memscalartype;  /* scalar type for mem (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */
558:   hid_t              filescalartype; /* scalar type for file (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */
559:   PetscInt           bs = PetscAbs(xin->map->bs);
560:   hsize_t            dim;
561:   hsize_t            maxDims[4], dims[4], chunkDims[4], count[4], offset[4];
562:   PetscBool          timestepping, dim2, spoutput;
563:   PetscInt           timestep = PETSC_MIN_INT, low;
564:   hsize_t            chunksize;
565:   const PetscScalar *x;
566:   const char        *vecname;

568:   PetscViewerHDF5OpenGroup(viewer, NULL, &file_id, &group);
569:   PetscViewerHDF5IsTimestepping(viewer, &timestepping);
570:   if (timestepping) PetscViewerHDF5GetTimestep(viewer, &timestep);
571:   PetscViewerHDF5GetBaseDimension2(viewer, &dim2);
572:   PetscViewerHDF5GetSPOutput(viewer, &spoutput);

574:   /* Create the dataspace for the dataset.
575:    *
576:    * dims - holds the current dimensions of the dataset
577:    *
578:    * maxDims - holds the maximum dimensions of the dataset (unlimited
579:    * for the number of time steps with the current dimensions for the
580:    * other dimensions; so only additional time steps can be added).
581:    *
582:    * chunkDims - holds the size of a single time step (required to
583:    * permit extending dataset).
584:    */
585:   dim       = 0;
586:   chunksize = 1;
587:   if (timestep >= 0) {
588:     dims[dim]      = timestep + 1;
589:     maxDims[dim]   = H5S_UNLIMITED;
590:     chunkDims[dim] = 1;
591:     ++dim;
592:   }
593:   PetscHDF5IntCast(xin->map->N / bs, dims + dim);

595:   maxDims[dim]   = dims[dim];
596:   chunkDims[dim] = PetscMax(1, dims[dim]);
597:   chunksize *= chunkDims[dim];
598:   ++dim;
599:   if (bs > 1 || dim2) {
600:     dims[dim]      = bs;
601:     maxDims[dim]   = dims[dim];
602:     chunkDims[dim] = PetscMax(1, dims[dim]);
603:     chunksize *= chunkDims[dim];
604:     ++dim;
605:   }
606:   #if defined(PETSC_USE_COMPLEX)
607:   dims[dim]      = 2;
608:   maxDims[dim]   = dims[dim];
609:   chunkDims[dim] = PetscMax(1, dims[dim]);
610:   chunksize *= chunkDims[dim];
611:   /* hdf5 chunks must be less than 4GB */
612:   if (chunksize > PETSC_HDF5_MAX_CHUNKSIZE / 64) {
613:     if (bs > 1 || dim2) {
614:       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));
615:       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));
616:     } else {
617:       chunkDims[dim - 1] = PETSC_HDF5_MAX_CHUNKSIZE / 128;
618:     }
619:   }
620:   ++dim;
621:   #else
622:   /* hdf5 chunks must be less than 4GB */
623:   if (chunksize > PETSC_HDF5_MAX_CHUNKSIZE / 64) {
624:     if (bs > 1 || dim2) {
625:       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));
626:       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));
627:     } else {
628:       chunkDims[dim - 1] = PETSC_HDF5_MAX_CHUNKSIZE / 64;
629:     }
630:   }
631:   #endif

633:   PetscCallHDF5Return(filespace, H5Screate_simple, (dim, dims, maxDims));

635:   #if defined(PETSC_USE_REAL_SINGLE)
636:   memscalartype  = H5T_NATIVE_FLOAT;
637:   filescalartype = H5T_NATIVE_FLOAT;
638:   #elif defined(PETSC_USE_REAL___FLOAT128)
639:     #error "HDF5 output with 128 bit floats not supported."
640:   #elif defined(PETSC_USE_REAL___FP16)
641:     #error "HDF5 output with 16 bit floats not supported."
642:   #else
643:   memscalartype = H5T_NATIVE_DOUBLE;
644:   if (spoutput == PETSC_TRUE) filescalartype = H5T_NATIVE_FLOAT;
645:   else filescalartype = H5T_NATIVE_DOUBLE;
646:   #endif

648:   /* Create the dataset with default properties and close filespace */
649:   PetscObjectGetName((PetscObject)xin, &vecname);
650:   if (H5Lexists(group, vecname, H5P_DEFAULT) < 1) {
651:     /* Create chunk */
652:     PetscCallHDF5Return(chunkspace, H5Pcreate, (H5P_DATASET_CREATE));
653:     PetscCallHDF5(H5Pset_chunk, (chunkspace, dim, chunkDims));

655:     PetscCallHDF5Return(dset_id, H5Dcreate2, (group, vecname, filescalartype, filespace, H5P_DEFAULT, chunkspace, H5P_DEFAULT));
656:     PetscCallHDF5(H5Pclose, (chunkspace));
657:   } else {
658:     PetscCallHDF5Return(dset_id, H5Dopen2, (group, vecname, H5P_DEFAULT));
659:     PetscCallHDF5(H5Dset_extent, (dset_id, dims));
660:   }
661:   PetscCallHDF5(H5Sclose, (filespace));

663:   /* Each process defines a dataset and writes it to the hyperslab in the file */
664:   dim = 0;
665:   if (timestep >= 0) {
666:     count[dim] = 1;
667:     ++dim;
668:   }
669:   PetscHDF5IntCast(xin->map->n / bs, count + dim);
670:   ++dim;
671:   if (bs > 1 || dim2) {
672:     count[dim] = bs;
673:     ++dim;
674:   }
675:   #if defined(PETSC_USE_COMPLEX)
676:   count[dim] = 2;
677:   ++dim;
678:   #endif
679:   if (xin->map->n > 0 || H5_VERSION_GE(1, 10, 0)) {
680:     PetscCallHDF5Return(memspace, H5Screate_simple, (dim, count, NULL));
681:   } else {
682:     /* Can't create dataspace with zero for any dimension, so create null dataspace. */
683:     PetscCallHDF5Return(memspace, H5Screate, (H5S_NULL));
684:   }

686:   /* Select hyperslab in the file */
687:   VecGetOwnershipRange(xin, &low, NULL);
688:   dim = 0;
689:   if (timestep >= 0) {
690:     offset[dim] = timestep;
691:     ++dim;
692:   }
693:   PetscHDF5IntCast(low / bs, offset + dim);
694:   ++dim;
695:   if (bs > 1 || dim2) {
696:     offset[dim] = 0;
697:     ++dim;
698:   }
699:   #if defined(PETSC_USE_COMPLEX)
700:   offset[dim] = 0;
701:   ++dim;
702:   #endif
703:   if (xin->map->n > 0 || H5_VERSION_GE(1, 10, 0)) {
704:     PetscCallHDF5Return(filespace, H5Dget_space, (dset_id));
705:     PetscCallHDF5(H5Sselect_hyperslab, (filespace, H5S_SELECT_SET, offset, NULL, count, NULL));
706:   } else {
707:     /* Create null filespace to match null memspace. */
708:     PetscCallHDF5Return(filespace, H5Screate, (H5S_NULL));
709:   }

711:   VecGetArrayRead(xin, &x);
712:   PetscCallHDF5(H5Dwrite, (dset_id, memscalartype, memspace, filespace, hdf5->dxpl_id, x));
713:   PetscCallHDF5(H5Fflush, (file_id, H5F_SCOPE_GLOBAL));
714:   VecRestoreArrayRead(xin, &x);

716:   /* Close/release resources */
717:   PetscCallHDF5(H5Gclose, (group));
718:   PetscCallHDF5(H5Sclose, (filespace));
719:   PetscCallHDF5(H5Sclose, (memspace));
720:   PetscCallHDF5(H5Dclose, (dset_id));

722:   #if defined(PETSC_USE_COMPLEX)
723:   {
724:     PetscBool tru = PETSC_TRUE;
725:     PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)xin, "complex", PETSC_BOOL, &tru);
726:   }
727:   #endif
728:   if (timestepping) PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)xin, "timestepping", PETSC_BOOL, &timestepping);
729:   PetscInfo(xin, "Wrote Vec object with name %s\n", vecname);
730:   return 0;
731: }
732: #endif

734: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec xin, PetscViewer viewer)
735: {
736:   PetscBool iascii, isbinary, isdraw;
737: #if defined(PETSC_HAVE_MATHEMATICA)
738:   PetscBool ismathematica;
739: #endif
740: #if defined(PETSC_HAVE_HDF5)
741:   PetscBool ishdf5;
742: #endif
743: #if defined(PETSC_HAVE_MATLAB)
744:   PetscBool ismatlab;
745: #endif
746: #if defined(PETSC_HAVE_ADIOS)
747:   PetscBool isadios;
748: #endif
749:   PetscBool isglvis;

751:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
752:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary);
753:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);
754: #if defined(PETSC_HAVE_MATHEMATICA)
755:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERMATHEMATICA, &ismathematica);
756: #endif
757: #if defined(PETSC_HAVE_HDF5)
758:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5);
759: #endif
760: #if defined(PETSC_HAVE_MATLAB)
761:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERMATLAB, &ismatlab);
762: #endif
763:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERGLVIS, &isglvis);
764: #if defined(PETSC_HAVE_ADIOS)
765:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERADIOS, &isadios);
766: #endif
767:   if (iascii) {
768:     VecView_MPI_ASCII(xin, viewer);
769:   } else if (isbinary) {
770:     VecView_MPI_Binary(xin, viewer);
771:   } else if (isdraw) {
772:     PetscViewerFormat format;
773:     PetscViewerGetFormat(viewer, &format);
774:     if (format == PETSC_VIEWER_DRAW_LG) {
775:       VecView_MPI_Draw_LG(xin, viewer);
776:     } else {
777:       VecView_MPI_Draw(xin, viewer);
778:     }
779: #if defined(PETSC_HAVE_MATHEMATICA)
780:   } else if (ismathematica) {
781:     PetscViewerMathematicaPutVector(viewer, xin);
782: #endif
783: #if defined(PETSC_HAVE_HDF5)
784:   } else if (ishdf5) {
785:     VecView_MPI_HDF5(xin, viewer);
786: #endif
787: #if defined(PETSC_HAVE_ADIOS)
788:   } else if (isadios) {
789:     VecView_MPI_ADIOS(xin, viewer);
790: #endif
791: #if defined(PETSC_HAVE_MATLAB)
792:   } else if (ismatlab) {
793:     VecView_MPI_Matlab(xin, viewer);
794: #endif
795:   } else if (isglvis) VecView_GLVis(xin, viewer);
796:   return 0;
797: }

799: PetscErrorCode VecGetSize_MPI(Vec xin, PetscInt *N)
800: {
801:   *N = xin->map->N;
802:   return 0;
803: }

805: PetscErrorCode VecGetValues_MPI(Vec xin, PetscInt ni, const PetscInt ix[], PetscScalar y[])
806: {
807:   const PetscScalar *xx;
808:   PetscInt           i, tmp, start = xin->map->range[xin->stash.rank];

810:   VecGetArrayRead(xin, &xx);
811:   for (i = 0; i < ni; i++) {
812:     if (xin->stash.ignorenegidx && ix[i] < 0) continue;
813:     tmp = ix[i] - start;
815:     y[i] = xx[tmp];
816:   }
817:   VecRestoreArrayRead(xin, &xx);
818:   return 0;
819: }

821: PetscErrorCode VecSetValues_MPI(Vec xin, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode addv)
822: {
823:   PetscMPIInt  rank   = xin->stash.rank;
824:   PetscInt    *owners = xin->map->range, start = owners[rank];
825:   PetscInt     end = owners[rank + 1], i, row;
826:   PetscScalar *xx;

828:   if (PetscDefined(USE_DEBUG)) {
831:   }
832:   VecGetArray(xin, &xx);
833:   xin->stash.insertmode = addv;

835:   if (addv == INSERT_VALUES) {
836:     for (i = 0; i < ni; i++) {
837:       if (xin->stash.ignorenegidx && ix[i] < 0) continue;
839:       if ((row = ix[i]) >= start && row < end) {
840:         xx[row - start] = y[i];
841:       } else if (!xin->stash.donotstash) {
843:         VecStashValue_Private(&xin->stash, row, y[i]);
844:       }
845:     }
846:   } else {
847:     for (i = 0; i < ni; i++) {
848:       if (xin->stash.ignorenegidx && ix[i] < 0) continue;
850:       if ((row = ix[i]) >= start && row < end) {
851:         xx[row - start] += y[i];
852:       } else if (!xin->stash.donotstash) {
854:         VecStashValue_Private(&xin->stash, row, y[i]);
855:       }
856:     }
857:   }
858:   VecRestoreArray(xin, &xx);
859:   return 0;
860: }

862: PetscErrorCode VecSetValuesBlocked_MPI(Vec xin, PetscInt ni, const PetscInt ix[], const PetscScalar yin[], InsertMode addv)
863: {
864:   PetscMPIInt  rank   = xin->stash.rank;
865:   PetscInt    *owners = xin->map->range, start = owners[rank];
866:   PetscInt     end = owners[rank + 1], i, row, bs = PetscAbs(xin->map->bs), j;
867:   PetscScalar *xx, *y = (PetscScalar *)yin;

869:   VecGetArray(xin, &xx);
870:   if (PetscDefined(USE_DEBUG)) {
873:   }
874:   xin->stash.insertmode = addv;

876:   if (addv == INSERT_VALUES) {
877:     for (i = 0; i < ni; i++) {
878:       if ((row = bs * ix[i]) >= start && row < end) {
879:         for (j = 0; j < bs; j++) xx[row - start + j] = y[j];
880:       } else if (!xin->stash.donotstash) {
881:         if (ix[i] < 0) {
882:           y += bs;
883:           continue;
884:         }
886:         VecStashValuesBlocked_Private(&xin->bstash, ix[i], y);
887:       }
888:       y += bs;
889:     }
890:   } else {
891:     for (i = 0; i < ni; i++) {
892:       if ((row = bs * ix[i]) >= start && row < end) {
893:         for (j = 0; j < bs; j++) xx[row - start + j] += y[j];
894:       } else if (!xin->stash.donotstash) {
895:         if (ix[i] < 0) {
896:           y += bs;
897:           continue;
898:         }
900:         VecStashValuesBlocked_Private(&xin->bstash, ix[i], y);
901:       }
902:       y += bs;
903:     }
904:   }
905:   VecRestoreArray(xin, &xx);
906:   return 0;
907: }

909: /*
910:    Since nsends or nreceives may be zero we add 1 in certain mallocs
911: to make sure we never malloc an empty one.
912: */
913: PetscErrorCode VecAssemblyBegin_MPI(Vec xin)
914: {
915:   PetscInt   *owners = xin->map->range, *bowners, i, bs, nstash, reallocs;
916:   PetscMPIInt size;
917:   InsertMode  addv;
918:   MPI_Comm    comm;

920:   PetscObjectGetComm((PetscObject)xin, &comm);
921:   if (xin->stash.donotstash) return 0;

923:   MPIU_Allreduce((PetscEnum *)&xin->stash.insertmode, (PetscEnum *)&addv, 1, MPIU_ENUM, MPI_BOR, comm);
925:   xin->stash.insertmode  = addv; /* in case this processor had no cache */
926:   xin->bstash.insertmode = addv; /* Block stash implicitly tracks InsertMode of scalar stash */

928:   VecGetBlockSize(xin, &bs);
929:   MPI_Comm_size(PetscObjectComm((PetscObject)xin), &size);
930:   if (!xin->bstash.bowners && xin->map->bs != -1) {
931:     PetscMalloc1(size + 1, &bowners);
932:     for (i = 0; i < size + 1; i++) bowners[i] = owners[i] / bs;
933:     xin->bstash.bowners = bowners;
934:   } else bowners = xin->bstash.bowners;

936:   VecStashScatterBegin_Private(&xin->stash, owners);
937:   VecStashScatterBegin_Private(&xin->bstash, bowners);
938:   VecStashGetInfo_Private(&xin->stash, &nstash, &reallocs);
939:   PetscInfo(xin, "Stash has %" PetscInt_FMT " entries, uses %" PetscInt_FMT " mallocs.\n", nstash, reallocs);
940:   VecStashGetInfo_Private(&xin->bstash, &nstash, &reallocs);
941:   PetscInfo(xin, "Block-Stash has %" PetscInt_FMT " entries, uses %" PetscInt_FMT " mallocs.\n", nstash, reallocs);
942:   return 0;
943: }

945: PetscErrorCode VecAssemblyEnd_MPI(Vec vec)
946: {
947:   PetscInt     base, i, j, *row, flg, bs;
948:   PetscMPIInt  n;
949:   PetscScalar *val, *vv, *array, *xarray;

951:   if (!vec->stash.donotstash) {
952:     VecGetArray(vec, &xarray);
953:     VecGetBlockSize(vec, &bs);
954:     base = vec->map->range[vec->stash.rank];

956:     /* Process the stash */
957:     while (1) {
958:       VecStashScatterGetMesg_Private(&vec->stash, &n, &row, &val, &flg);
959:       if (!flg) break;
960:       if (vec->stash.insertmode == ADD_VALUES) {
961:         for (i = 0; i < n; i++) xarray[row[i] - base] += val[i];
962:       } else if (vec->stash.insertmode == INSERT_VALUES) {
963:         for (i = 0; i < n; i++) xarray[row[i] - base] = val[i];
964:       } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_CORRUPT, "Insert mode is not set correctly; corrupted vector");
965:     }
966:     VecStashScatterEnd_Private(&vec->stash);

968:     /* now process the block-stash */
969:     while (1) {
970:       VecStashScatterGetMesg_Private(&vec->bstash, &n, &row, &val, &flg);
971:       if (!flg) break;
972:       for (i = 0; i < n; i++) {
973:         array = xarray + row[i] * bs - base;
974:         vv    = val + i * bs;
975:         if (vec->stash.insertmode == ADD_VALUES) {
976:           for (j = 0; j < bs; j++) array[j] += vv[j];
977:         } else if (vec->stash.insertmode == INSERT_VALUES) {
978:           for (j = 0; j < bs; j++) array[j] = vv[j];
979:         } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_CORRUPT, "Insert mode is not set correctly; corrupted vector");
980:       }
981:     }
982:     VecStashScatterEnd_Private(&vec->bstash);
983:     VecRestoreArray(vec, &xarray);
984:   }
985:   vec->stash.insertmode = NOT_SET_VALUES;
986:   return 0;
987: }

989: PetscErrorCode VecSetPreallocationCOO_MPI(Vec x, PetscCount coo_n, const PetscInt coo_i[])
990: {
991:   PetscInt    m, M, rstart, rend;
992:   Vec_MPI    *vmpi = (Vec_MPI *)x->data;
993:   PetscCount  k, p, q, rem; /* Loop variables over coo arrays */
994:   PetscMPIInt size;
995:   MPI_Comm    comm;

997:   PetscObjectGetComm((PetscObject)x, &comm);
998:   MPI_Comm_size(comm, &size);
999:   VecResetPreallocationCOO_MPI(x);

1001:   PetscLayoutSetUp(x->map);
1002:   VecGetOwnershipRange(x, &rstart, &rend);
1003:   VecGetLocalSize(x, &m);
1004:   VecGetSize(x, &M);

1006:   /* ---------------------------------------------------------------------------*/
1007:   /* Sort COOs along with a permuation array, so that negative indices come     */
1008:   /* first, then local ones, then remote ones.                                  */
1009:   /* ---------------------------------------------------------------------------*/
1010:   PetscCount n1 = coo_n, nneg, *perm;
1011:   PetscInt  *i1; /* Copy of input COOs along with a permutation array */
1012:   PetscMalloc1(n1, &i1);
1013:   PetscMalloc1(n1, &perm);
1014:   PetscArraycpy(i1, coo_i, n1); /* Make a copy since we'll modify it */
1015:   for (k = 0; k < n1; k++) perm[k] = k;

1017:   /* Manipulate i1[] so that entries with negative indices will have the smallest
1018:      index, local entries will have greater but negative indices, and remote entries
1019:      will have positive indices.
1020:   */
1021:   for (k = 0; k < n1; k++) {
1022:     if (i1[k] < 0) {
1023:       if (x->stash.ignorenegidx) i1[k] = PETSC_MIN_INT; /* e.g., -2^31, minimal to move them ahead */
1024:       else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Found a negative index in VecSetPreallocateCOO() but VEC_IGNORE_NEGATIVE_INDICES was not set");
1025:     } else if (i1[k] >= rstart && i1[k] < rend) {
1026:       i1[k] -= PETSC_MAX_INT; /* e.g., minus 2^31-1 to shift local rows to range of [-PETSC_MAX_INT, -1] */
1027:     } else {
1029:       if (x->stash.donotstash) i1[k] = PETSC_MIN_INT; /* Ignore off-proc indices as if they were negative */
1030:     }
1031:   }

1033:   /* Sort the indices, after that, [0,nneg) have ignored entires, [nneg,rem) have local entries and [rem,n1) have remote entries */
1034:   PetscSortIntWithCountArray(n1, i1, perm);
1035:   for (k = 0; k < n1; k++) {
1036:     if (i1[k] > PETSC_MIN_INT) break;
1037:   } /* Advance k to the first entry we need to take care of */
1038:   nneg = k;
1039:   PetscSortedIntUpperBound(i1, nneg, n1, rend - 1 - PETSC_MAX_INT, &rem); /* rem is upper bound of the last local row */
1040:   for (k = nneg; k < rem; k++) i1[k] += PETSC_MAX_INT;                               /* Revert indices of local entries */

1042:   /* ---------------------------------------------------------------------------*/
1043:   /*           Build stuff for local entries                                    */
1044:   /* ---------------------------------------------------------------------------*/
1045:   PetscCount tot1, *jmap1, *perm1;
1046:   PetscCalloc1(m + 1, &jmap1);
1047:   for (k = nneg; k < rem; k++) jmap1[i1[k] - rstart + 1]++; /* Count repeats of each local entry */
1048:   for (k = 0; k < m; k++) jmap1[k + 1] += jmap1[k];         /* Transform jmap1[] to CSR-like data structure */
1049:   tot1 = jmap1[m];
1050:   PetscAssert(tot1 == rem - nneg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected errors in VecSetPreallocationCOO_MPI");
1051:   PetscMalloc1(tot1, &perm1);
1052:   PetscArraycpy(perm1, perm + nneg, tot1);

1054:   /* ---------------------------------------------------------------------------*/
1055:   /*        Record the permutation array for filling the send buffer            */
1056:   /* ---------------------------------------------------------------------------*/
1057:   PetscCount *Cperm;
1058:   PetscMalloc1(n1 - rem, &Cperm);
1059:   PetscArraycpy(Cperm, perm + rem, n1 - rem);
1060:   PetscFree(perm);

1062:   /* ---------------------------------------------------------------------------*/
1063:   /*           Send remote entries to their owner                                  */
1064:   /* ---------------------------------------------------------------------------*/
1065:   /* Find which entries should be sent to which remote ranks*/
1066:   PetscInt        nsend = 0; /* Number of MPI ranks to send data to */
1067:   PetscMPIInt    *sendto;    /* [nsend], storing remote ranks */
1068:   PetscInt       *nentries;  /* [nsend], storing number of entries sent to remote ranks; Assume PetscInt is big enough for this count, and error if not */
1069:   const PetscInt *ranges;
1070:   PetscInt        maxNsend = size >= 128 ? 128 : size; /* Assume max 128 neighbors; realloc when needed */

1072:   PetscLayoutGetRanges(x->map, &ranges);
1073:   PetscMalloc2(maxNsend, &sendto, maxNsend, &nentries);
1074:   for (k = rem; k < n1;) {
1075:     PetscMPIInt owner;
1076:     PetscInt    firstRow, lastRow;

1078:     /* Locate a row range */
1079:     firstRow = i1[k]; /* first row of this owner */
1080:     PetscLayoutFindOwner(x->map, firstRow, &owner);
1081:     lastRow = ranges[owner + 1] - 1; /* last row of this owner */

1083:     /* Find the first index 'p' in [k,n) with i[p] belonging to next owner */
1084:     PetscSortedIntUpperBound(i1, k, n1, lastRow, &p);

1086:     /* All entries in [k,p) belong to this remote owner */
1087:     if (nsend >= maxNsend) { /* Double the remote ranks arrays if not long enough */
1088:       PetscMPIInt *sendto2;
1089:       PetscInt    *nentries2;
1090:       PetscInt     maxNsend2 = (maxNsend <= size / 2) ? maxNsend * 2 : size;

1092:       PetscMalloc2(maxNsend2, &sendto2, maxNsend2, &nentries2);
1093:       PetscArraycpy(sendto2, sendto, maxNsend);
1094:       PetscArraycpy(nentries2, nentries2, maxNsend + 1);
1095:       PetscFree2(sendto, nentries2);
1096:       sendto   = sendto2;
1097:       nentries = nentries2;
1098:       maxNsend = maxNsend2;
1099:     }
1100:     sendto[nsend]   = owner;
1101:     nentries[nsend] = p - k;
1102:     PetscCountCast(p - k, &nentries[nsend]);
1103:     nsend++;
1104:     k = p;
1105:   }

1107:   /* Build 1st SF to know offsets on remote to send data */
1108:   PetscSF      sf1;
1109:   PetscInt     nroots = 1, nroots2 = 0;
1110:   PetscInt     nleaves = nsend, nleaves2 = 0;
1111:   PetscInt    *offsets;
1112:   PetscSFNode *iremote;

1114:   PetscSFCreate(comm, &sf1);
1115:   PetscMalloc1(nsend, &iremote);
1116:   PetscMalloc1(nsend, &offsets);
1117:   for (k = 0; k < nsend; k++) {
1118:     iremote[k].rank  = sendto[k];
1119:     iremote[k].index = 0;
1120:     nleaves2 += nentries[k];
1122:   }
1123:   PetscSFSetGraph(sf1, nroots, nleaves, NULL, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER);
1124:   PetscSFFetchAndOpWithMemTypeBegin(sf1, MPIU_INT, PETSC_MEMTYPE_HOST, &nroots2 /*rootdata*/, PETSC_MEMTYPE_HOST, nentries /*leafdata*/, PETSC_MEMTYPE_HOST, offsets /*leafupdate*/, MPI_SUM);
1125:   PetscSFFetchAndOpEnd(sf1, MPIU_INT, &nroots2, nentries, offsets, MPI_SUM); /* Would nroots2 overflow, we check offsets[] below */
1126:   PetscSFDestroy(&sf1);
1127:   PetscAssert(nleaves2 == n1 - rem, PETSC_COMM_SELF, PETSC_ERR_PLIB, "nleaves2 %" PetscInt_FMT " != number of remote entries %" PetscCount_FMT "", nleaves2, n1 - rem);

1129:   /* Build 2nd SF to send remote COOs to their owner */
1130:   PetscSF sf2;
1131:   nroots  = nroots2;
1132:   nleaves = nleaves2;
1133:   PetscSFCreate(comm, &sf2);
1134:   PetscSFSetFromOptions(sf2);
1135:   PetscMalloc1(nleaves, &iremote);
1136:   p = 0;
1137:   for (k = 0; k < nsend; k++) {
1139:     for (q = 0; q < nentries[k]; q++, p++) {
1140:       iremote[p].rank  = sendto[k];
1141:       iremote[p].index = offsets[k] + q;
1142:     }
1143:   }
1144:   PetscSFSetGraph(sf2, nroots, nleaves, NULL, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER);

1146:   /* Send the remote COOs to their owner */
1147:   PetscInt    n2 = nroots, *i2; /* Buffers for received COOs from other ranks, along with a permutation array */
1148:   PetscCount *perm2;
1149:   PetscMalloc1(n2, &i2);
1150:   PetscMalloc1(n2, &perm2);
1151:   PetscSFReduceWithMemTypeBegin(sf2, MPIU_INT, PETSC_MEMTYPE_HOST, i1 + rem, PETSC_MEMTYPE_HOST, i2, MPI_REPLACE);
1152:   PetscSFReduceEnd(sf2, MPIU_INT, i1 + rem, i2, MPI_REPLACE);

1154:   PetscFree(i1);
1155:   PetscFree(offsets);
1156:   PetscFree2(sendto, nentries);

1158:   /* ---------------------------------------------------------------*/
1159:   /* Sort received COOs along with a permutation array            */
1160:   /* ---------------------------------------------------------------*/
1161:   PetscCount  *imap2;
1162:   PetscCount  *jmap2, nnz2;
1163:   PetscScalar *sendbuf, *recvbuf;
1164:   PetscInt     old;
1165:   PetscCount   sendlen = n1 - rem, recvlen = n2;

1167:   for (k = 0; k < n2; k++) perm2[k] = k;
1168:   PetscSortIntWithCountArray(n2, i2, perm2);

1170:   /* nnz2 will be # of unique entries in the recvbuf */
1171:   nnz2 = n2;
1172:   for (k = 1; k < n2; k++) {
1173:     if (i2[k] == i2[k - 1]) nnz2--;
1174:   }

1176:   /* Build imap2[] and jmap2[] for each unique entry */
1177:   PetscMalloc4(nnz2, &imap2, nnz2 + 1, &jmap2, sendlen, &sendbuf, recvlen, &recvbuf);
1178:   p        = -1;
1179:   old      = -1;
1180:   jmap2[0] = 0;
1181:   jmap2++;
1182:   for (k = 0; k < n2; k++) {
1183:     if (i2[k] != old) { /* Meet a new entry */
1184:       p++;
1185:       imap2[p] = i2[k] - rstart;
1186:       jmap2[p] = 1;
1187:       old      = i2[k];
1188:     } else {
1189:       jmap2[p]++;
1190:     }
1191:   }
1192:   jmap2--;
1193:   for (k = 0; k < nnz2; k++) jmap2[k + 1] += jmap2[k];

1195:   PetscFree(i2);

1197:   vmpi->coo_n = coo_n;
1198:   vmpi->tot1  = tot1;
1199:   vmpi->jmap1 = jmap1;
1200:   vmpi->perm1 = perm1;
1201:   vmpi->nnz2  = nnz2;
1202:   vmpi->imap2 = imap2;
1203:   vmpi->jmap2 = jmap2;
1204:   vmpi->perm2 = perm2;

1206:   vmpi->Cperm   = Cperm;
1207:   vmpi->sendbuf = sendbuf;
1208:   vmpi->recvbuf = recvbuf;
1209:   vmpi->sendlen = sendlen;
1210:   vmpi->recvlen = recvlen;
1211:   vmpi->coo_sf  = sf2;
1212:   return 0;
1213: }

1215: PetscErrorCode VecSetValuesCOO_MPI(Vec x, const PetscScalar v[], InsertMode imode)
1216: {
1217:   Vec_MPI          *vmpi = (Vec_MPI *)x->data;
1218:   PetscInt          m;
1219:   PetscScalar      *a, *sendbuf = vmpi->sendbuf, *recvbuf = vmpi->recvbuf;
1220:   const PetscCount *jmap1 = vmpi->jmap1;
1221:   const PetscCount *perm1 = vmpi->perm1;
1222:   const PetscCount *imap2 = vmpi->imap2;
1223:   const PetscCount *jmap2 = vmpi->jmap2;
1224:   const PetscCount *perm2 = vmpi->perm2;
1225:   const PetscCount *Cperm = vmpi->Cperm;
1226:   const PetscCount  nnz2  = vmpi->nnz2;

1228:   VecGetLocalSize(x, &m);
1229:   VecGetArray(x, &a);

1231:   /* Pack entries to be sent to remote */
1232:   for (PetscInt i = 0; i < vmpi->sendlen; i++) sendbuf[i] = v[Cperm[i]];

1234:   /* Send remote entries to their owner and overlap the communication with local computation */
1235:   PetscSFReduceWithMemTypeBegin(vmpi->coo_sf, MPIU_SCALAR, PETSC_MEMTYPE_HOST, sendbuf, PETSC_MEMTYPE_HOST, recvbuf, MPI_REPLACE);
1236:   /* Add local entries to A and B */
1237:   for (PetscInt i = 0; i < m; i++) { /* All entries in a[] are either zero'ed or added with a value (i.e., initialized) */
1238:     PetscScalar sum = 0.0;           /* Do partial summation first to improve numerical stablility */
1239:     for (PetscCount k = jmap1[i]; k < jmap1[i + 1]; k++) sum += v[perm1[k]];
1240:     a[i] = (imode == INSERT_VALUES ? 0.0 : a[i]) + sum;
1241:   }
1242:   PetscSFReduceEnd(vmpi->coo_sf, MPIU_SCALAR, sendbuf, recvbuf, MPI_REPLACE);

1244:   /* Add received remote entries to A and B */
1245:   for (PetscInt i = 0; i < nnz2; i++) {
1246:     for (PetscCount k = jmap2[i]; k < jmap2[i + 1]; k++) a[imap2[i]] += recvbuf[perm2[k]];
1247:   }

1249:   VecRestoreArray(x, &a);
1250:   return 0;
1251: }