Actual source code: pdvec.c


  2: /*
  3:      Code for some of the parallel vector primitives.
  4: */
  5: #include <../src/vec/vec/impls/mpi/pvecimpl.h>
  6: #include <petsc/private/viewerhdf5impl.h>
  7: #include <petsc/private/glvisviewerimpl.h>
  8: #include <petsc/private/glvisvecimpl.h>
  9: #include <petscsf.h>

 11: static PetscErrorCode VecResetPreallocationCOO_MPI(Vec v)
 12: {
 13:   Vec_MPI *vmpi = (Vec_MPI *)v->data;

 15:   PetscFunctionBegin;
 16:   if (vmpi) {
 17:     PetscCall(PetscFree(vmpi->jmap1));
 18:     PetscCall(PetscFree(vmpi->perm1));
 19:     PetscCall(PetscFree(vmpi->Cperm));
 20:     PetscCall(PetscFree4(vmpi->imap2, vmpi->jmap2, vmpi->sendbuf, vmpi->recvbuf));
 21:     PetscCall(PetscFree(vmpi->perm2));
 22:     PetscCall(PetscSFDestroy(&vmpi->coo_sf));
 23:   }
 24:   PetscFunctionReturn(PETSC_SUCCESS);
 25: }

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

 31:   PetscFunctionBegin;
 32: #if defined(PETSC_USE_LOG)
 33:   PetscCall(PetscLogObjectState((PetscObject)v, "Length=%" PetscInt_FMT, v->map->N));
 34: #endif
 35:   if (!x) PetscFunctionReturn(PETSC_SUCCESS);
 36:   PetscCall(PetscFree(x->array_allocated));

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

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

 49:   PetscCall(VecResetPreallocationCOO_MPI(v));
 50:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscMatlabEnginePut_C", NULL));
 51:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscMatlabEngineGet_C", NULL));
 52:   PetscCall(PetscFree(v->data));
 53:   PetscFunctionReturn(PETSC_SUCCESS);
 54: }

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

 66:   PetscFunctionBegin;
 67:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)xin), &size));
 68:   PetscCall(PetscViewerGetFormat(viewer, &format));
 69:   if (format == PETSC_VIEWER_LOAD_BALANCE) {
 70:     PetscInt nmax = 0, nmin = xin->map->n, navg;
 71:     for (i = 0; i < (PetscInt)size; i++) {
 72:       nmax = PetscMax(nmax, xin->map->range[i + 1] - xin->map->range[i]);
 73:       nmin = PetscMin(nmin, xin->map->range[i + 1] - xin->map->range[i]);
 74:     }
 75:     navg = xin->map->N / size;
 76:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Load Balance - Local vector size Min %" PetscInt_FMT "  avg %" PetscInt_FMT "  max %" PetscInt_FMT "\n", nmin, navg, nmax));
 77:     PetscFunctionReturn(PETSC_SUCCESS);
 78:   }

 80:   PetscCall(VecGetArrayRead(xin, &xarray));
 81:   /* determine maximum message to arrive */
 82:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)xin), &rank));
 83:   PetscCallMPI(MPI_Reduce(&work, &len, 1, MPIU_INT, MPI_MAX, 0, PetscObjectComm((PetscObject)xin)));
 84:   if (format == PETSC_VIEWER_ASCII_GLVIS) rank = 0, len = 0; /* no parallel distributed write support from GLVis */
 85:   if (rank == 0) {
 86:     PetscCall(PetscMalloc1(len, &values));
 87:     /*
 88:         MATLAB format and ASCII format are very similar except
 89:         MATLAB uses %18.16e format while ASCII uses %g
 90:     */
 91:     if (format == PETSC_VIEWER_ASCII_MATLAB) {
 92:       PetscCall(PetscObjectGetName((PetscObject)xin, &name));
 93:       PetscCall(PetscViewerASCIIPrintf(viewer, "%s = [\n", name));
 94:       for (i = 0; i < xin->map->n; i++) {
 95: #if defined(PETSC_USE_COMPLEX)
 96:         if (PetscImaginaryPart(xarray[i]) > 0.0) {
 97:           PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e + %18.16ei\n", (double)PetscRealPart(xarray[i]), (double)PetscImaginaryPart(xarray[i])));
 98:         } else if (PetscImaginaryPart(xarray[i]) < 0.0) {
 99:           PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e - %18.16ei\n", (double)PetscRealPart(xarray[i]), -(double)PetscImaginaryPart(xarray[i])));
100:         } else {
101:           PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e\n", (double)PetscRealPart(xarray[i])));
102:         }
103: #else
104:         PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e\n", (double)xarray[i]));
105: #endif
106:       }
107:       /* receive and print messages */
108:       for (j = 1; j < size; j++) {
109:         PetscCallMPI(MPI_Recv(values, (PetscMPIInt)len, MPIU_SCALAR, j, tag, PetscObjectComm((PetscObject)xin), &status));
110:         PetscCallMPI(MPI_Get_count(&status, MPIU_SCALAR, &n));
111:         for (i = 0; i < n; i++) {
112: #if defined(PETSC_USE_COMPLEX)
113:           if (PetscImaginaryPart(values[i]) > 0.0) {
114:             PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e + %18.16ei\n", (double)PetscRealPart(values[i]), (double)PetscImaginaryPart(values[i])));
115:           } else if (PetscImaginaryPart(values[i]) < 0.0) {
116:             PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e - %18.16ei\n", (double)PetscRealPart(values[i]), -(double)PetscImaginaryPart(values[i])));
117:           } else {
118:             PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e\n", (double)PetscRealPart(values[i])));
119:           }
120: #else
121:           PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e\n", (double)values[i]));
122: #endif
123:         }
124:       }
125:       PetscCall(PetscViewerASCIIPrintf(viewer, "];\n"));

127:     } else if (format == PETSC_VIEWER_ASCII_SYMMODU) {
128:       for (i = 0; i < xin->map->n; i++) {
129: #if defined(PETSC_USE_COMPLEX)
130:         PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e %18.16e\n", (double)PetscRealPart(xarray[i]), (double)PetscImaginaryPart(xarray[i])));
131: #else
132:         PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e\n", (double)xarray[i]));
133: #endif
134:       }
135:       /* receive and print messages */
136:       for (j = 1; j < size; j++) {
137:         PetscCallMPI(MPI_Recv(values, (PetscMPIInt)len, MPIU_SCALAR, j, tag, PetscObjectComm((PetscObject)xin), &status));
138:         PetscCallMPI(MPI_Get_count(&status, MPIU_SCALAR, &n));
139:         for (i = 0; i < n; i++) {
140: #if defined(PETSC_USE_COMPLEX)
141:           PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e %18.16e\n", (double)PetscRealPart(values[i]), (double)PetscImaginaryPart(values[i])));
142: #else
143:           PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e\n", (double)values[i]));
144: #endif
145:         }
146:       }
147:     } else if (format == PETSC_VIEWER_ASCII_VTK_DEPRECATED || format == PETSC_VIEWER_ASCII_VTK_CELL_DEPRECATED) {
148:       /*
149:         state 0: No header has been output
150:         state 1: Only POINT_DATA has been output
151:         state 2: Only CELL_DATA has been output
152:         state 3: Output both, POINT_DATA last
153:         state 4: Output both, CELL_DATA last
154:       */
155:       static PetscInt stateId     = -1;
156:       int             outputState = 0;
157:       int             doOutput    = 0;
158:       PetscBool       hasState;
159:       PetscInt        bs, b;

161:       if (stateId < 0) PetscCall(PetscObjectComposedDataRegister(&stateId));
162:       PetscCall(PetscObjectComposedDataGetInt((PetscObject)viewer, stateId, outputState, hasState));
163:       if (!hasState) outputState = 0;

165:       PetscCall(PetscObjectGetName((PetscObject)xin, &name));
166:       PetscCall(VecGetLocalSize(xin, &nLen));
167:       PetscCall(PetscMPIIntCast(nLen, &n));
168:       PetscCall(VecGetBlockSize(xin, &bs));
169:       if (format == PETSC_VIEWER_ASCII_VTK_DEPRECATED) {
170:         if (outputState == 0) {
171:           outputState = 1;
172:           doOutput    = 1;
173:         } else if (outputState == 1) doOutput = 0;
174:         else if (outputState == 2) {
175:           outputState = 3;
176:           doOutput    = 1;
177:         } else if (outputState == 3) doOutput = 0;
178:         else PetscCheck(outputState != 4, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Tried to output POINT_DATA again after intervening CELL_DATA");

180:         if (doOutput) PetscCall(PetscViewerASCIIPrintf(viewer, "POINT_DATA %" PetscInt_FMT "\n", xin->map->N / bs));
181:       } else {
182:         if (outputState == 0) {
183:           outputState = 2;
184:           doOutput    = 1;
185:         } else if (outputState == 1) {
186:           outputState = 4;
187:           doOutput    = 1;
188:         } else if (outputState == 2) {
189:           doOutput = 0;
190:         } else {
191:           PetscCheck(outputState != 3, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Tried to output CELL_DATA again after intervening POINT_DATA");
192:           if (outputState == 4) doOutput = 0;
193:         }

195:         if (doOutput) PetscCall(PetscViewerASCIIPrintf(viewer, "CELL_DATA %" PetscInt_FMT "\n", xin->map->N / bs));
196:       }
197:       PetscCall(PetscObjectComposedDataSetInt((PetscObject)viewer, stateId, outputState));
198:       if (name) {
199:         if (bs == 3) {
200:           PetscCall(PetscViewerASCIIPrintf(viewer, "VECTORS %s double\n", name));
201:         } else {
202:           PetscCall(PetscViewerASCIIPrintf(viewer, "SCALARS %s double %" PetscInt_FMT "\n", name, bs));
203:         }
204:       } else {
205:         PetscCall(PetscViewerASCIIPrintf(viewer, "SCALARS scalars double %" PetscInt_FMT "\n", bs));
206:       }
207:       if (bs != 3) PetscCall(PetscViewerASCIIPrintf(viewer, "LOOKUP_TABLE default\n"));
208:       for (i = 0; i < n / bs; i++) {
209:         for (b = 0; b < bs; b++) {
210:           if (b > 0) PetscCall(PetscViewerASCIIPrintf(viewer, " "));
211:           PetscCall(PetscViewerASCIIPrintf(viewer, "%g", (double)PetscRealPart(xarray[i * bs + b])));
212:         }
213:         PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
214:       }
215:       for (j = 1; j < size; j++) {
216:         PetscCallMPI(MPI_Recv(values, (PetscMPIInt)len, MPIU_SCALAR, j, tag, PetscObjectComm((PetscObject)xin), &status));
217:         PetscCallMPI(MPI_Get_count(&status, MPIU_SCALAR, &n));
218:         for (i = 0; i < n / bs; i++) {
219:           for (b = 0; b < bs; b++) {
220:             if (b > 0) PetscCall(PetscViewerASCIIPrintf(viewer, " "));
221:             PetscCall(PetscViewerASCIIPrintf(viewer, "%g", (double)PetscRealPart(values[i * bs + b])));
222:           }
223:           PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
224:         }
225:       }
226:     } else if (format == PETSC_VIEWER_ASCII_VTK_COORDS_DEPRECATED) {
227:       PetscInt bs, b;

229:       PetscCall(VecGetLocalSize(xin, &nLen));
230:       PetscCall(PetscMPIIntCast(nLen, &n));
231:       PetscCall(VecGetBlockSize(xin, &bs));
232:       PetscCheck(bs >= 1 && bs <= 3, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "VTK can only handle 3D objects, but vector dimension is %" PetscInt_FMT, bs);

234:       for (i = 0; i < n / bs; i++) {
235:         for (b = 0; b < bs; b++) {
236:           if (b > 0) PetscCall(PetscViewerASCIIPrintf(viewer, " "));
237:           PetscCall(PetscViewerASCIIPrintf(viewer, "%g", (double)PetscRealPart(xarray[i * bs + b])));
238:         }
239:         for (b = bs; b < 3; b++) PetscCall(PetscViewerASCIIPrintf(viewer, " 0.0"));
240:         PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
241:       }
242:       for (j = 1; j < size; j++) {
243:         PetscCallMPI(MPI_Recv(values, (PetscMPIInt)len, MPIU_SCALAR, j, tag, PetscObjectComm((PetscObject)xin), &status));
244:         PetscCallMPI(MPI_Get_count(&status, MPIU_SCALAR, &n));
245:         for (i = 0; i < n / bs; i++) {
246:           for (b = 0; b < bs; b++) {
247:             if (b > 0) PetscCall(PetscViewerASCIIPrintf(viewer, " "));
248:             PetscCall(PetscViewerASCIIPrintf(viewer, "%g", (double)PetscRealPart(values[i * bs + b])));
249:           }
250:           for (b = bs; b < 3; b++) PetscCall(PetscViewerASCIIPrintf(viewer, " 0.0"));
251:           PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
252:         }
253:       }
254:     } else if (format == PETSC_VIEWER_ASCII_PCICE) {
255:       PetscInt bs, b, vertexCount = 1;

257:       PetscCall(VecGetLocalSize(xin, &nLen));
258:       PetscCall(PetscMPIIntCast(nLen, &n));
259:       PetscCall(VecGetBlockSize(xin, &bs));
260:       PetscCheck(bs >= 1 && bs <= 3, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "PCICE can only handle up to 3D objects, but vector dimension is %" PetscInt_FMT, bs);

262:       PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT "\n", xin->map->N / bs));
263:       for (i = 0; i < n / bs; i++) {
264:         PetscCall(PetscViewerASCIIPrintf(viewer, "%7" PetscInt_FMT "   ", vertexCount++));
265:         for (b = 0; b < bs; b++) {
266:           if (b > 0) PetscCall(PetscViewerASCIIPrintf(viewer, " "));
267: #if !defined(PETSC_USE_COMPLEX)
268:           PetscCall(PetscViewerASCIIPrintf(viewer, "% 12.5E", (double)xarray[i * bs + b]));
269: #endif
270:         }
271:         PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
272:       }
273:       for (j = 1; j < size; j++) {
274:         PetscCallMPI(MPI_Recv(values, (PetscMPIInt)len, MPIU_SCALAR, j, tag, PetscObjectComm((PetscObject)xin), &status));
275:         PetscCallMPI(MPI_Get_count(&status, MPIU_SCALAR, &n));
276:         for (i = 0; i < n / bs; i++) {
277:           PetscCall(PetscViewerASCIIPrintf(viewer, "%7" PetscInt_FMT "   ", vertexCount++));
278:           for (b = 0; b < bs; b++) {
279:             if (b > 0) PetscCall(PetscViewerASCIIPrintf(viewer, " "));
280: #if !defined(PETSC_USE_COMPLEX)
281:             PetscCall(PetscViewerASCIIPrintf(viewer, "% 12.5E", (double)values[i * bs + b]));
282: #endif
283:           }
284:           PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
285:         }
286:       }
287:     } else if (format == PETSC_VIEWER_ASCII_GLVIS) {
288:       /* GLVis ASCII visualization/dump: this function mimics mfem::GridFunction::Save() */
289:       const PetscScalar      *array;
290:       PetscInt                i, n, vdim, ordering = 1; /* mfem::FiniteElementSpace::Ordering::byVDIM */
291:       PetscContainer          glvis_container;
292:       PetscViewerGLVisVecInfo glvis_vec_info;
293:       PetscViewerGLVisInfo    glvis_info;

295:       /* mfem::FiniteElementSpace::Save() */
296:       PetscCall(VecGetBlockSize(xin, &vdim));
297:       PetscCall(PetscViewerASCIIPrintf(viewer, "FiniteElementSpace\n"));
298:       PetscCall(PetscObjectQuery((PetscObject)xin, "_glvis_info_container", (PetscObject *)&glvis_container));
299:       PetscCheck(glvis_container, PetscObjectComm((PetscObject)xin), PETSC_ERR_PLIB, "Missing GLVis container");
300:       PetscCall(PetscContainerGetPointer(glvis_container, (void **)&glvis_vec_info));
301:       PetscCall(PetscViewerASCIIPrintf(viewer, "%s\n", glvis_vec_info->fec_type));
302:       PetscCall(PetscViewerASCIIPrintf(viewer, "VDim: %" PetscInt_FMT "\n", vdim));
303:       PetscCall(PetscViewerASCIIPrintf(viewer, "Ordering: %" PetscInt_FMT "\n", ordering));
304:       PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
305:       /* mfem::Vector::Print() */
306:       PetscCall(PetscObjectQuery((PetscObject)viewer, "_glvis_info_container", (PetscObject *)&glvis_container));
307:       PetscCheck(glvis_container, PetscObjectComm((PetscObject)viewer), PETSC_ERR_PLIB, "Missing GLVis container");
308:       PetscCall(PetscContainerGetPointer(glvis_container, (void **)&glvis_info));
309:       if (glvis_info->enabled) {
310:         PetscCall(VecGetLocalSize(xin, &n));
311:         PetscCall(VecGetArrayRead(xin, &array));
312:         for (i = 0; i < n; i++) {
313:           PetscCall(PetscViewerASCIIPrintf(viewer, glvis_info->fmt, (double)PetscRealPart(array[i])));
314:           PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
315:         }
316:         PetscCall(VecRestoreArrayRead(xin, &array));
317:       }
318:     } else if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
319:       /* No info */
320:     } else {
321:       if (format != PETSC_VIEWER_ASCII_COMMON) PetscCall(PetscViewerASCIIPrintf(viewer, "Process [%d]\n", rank));
322:       cnt = 0;
323:       for (i = 0; i < xin->map->n; i++) {
324:         if (format == PETSC_VIEWER_ASCII_INDEX) PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT ": ", cnt++));
325: #if defined(PETSC_USE_COMPLEX)
326:         if (PetscImaginaryPart(xarray[i]) > 0.0) {
327:           PetscCall(PetscViewerASCIIPrintf(viewer, "%g + %g i\n", (double)PetscRealPart(xarray[i]), (double)PetscImaginaryPart(xarray[i])));
328:         } else if (PetscImaginaryPart(xarray[i]) < 0.0) {
329:           PetscCall(PetscViewerASCIIPrintf(viewer, "%g - %g i\n", (double)PetscRealPart(xarray[i]), -(double)PetscImaginaryPart(xarray[i])));
330:         } else {
331:           PetscCall(PetscViewerASCIIPrintf(viewer, "%g\n", (double)PetscRealPart(xarray[i])));
332:         }
333: #else
334:         PetscCall(PetscViewerASCIIPrintf(viewer, "%g\n", (double)xarray[i]));
335: #endif
336:       }
337:       /* receive and print messages */
338:       for (j = 1; j < size; j++) {
339:         PetscCallMPI(MPI_Recv(values, (PetscMPIInt)len, MPIU_SCALAR, j, tag, PetscObjectComm((PetscObject)xin), &status));
340:         PetscCallMPI(MPI_Get_count(&status, MPIU_SCALAR, &n));
341:         if (format != PETSC_VIEWER_ASCII_COMMON) PetscCall(PetscViewerASCIIPrintf(viewer, "Process [%d]\n", j));
342:         for (i = 0; i < n; i++) {
343:           if (format == PETSC_VIEWER_ASCII_INDEX) PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT ": ", cnt++));
344: #if defined(PETSC_USE_COMPLEX)
345:           if (PetscImaginaryPart(values[i]) > 0.0) {
346:             PetscCall(PetscViewerASCIIPrintf(viewer, "%g + %g i\n", (double)PetscRealPart(values[i]), (double)PetscImaginaryPart(values[i])));
347:           } else if (PetscImaginaryPart(values[i]) < 0.0) {
348:             PetscCall(PetscViewerASCIIPrintf(viewer, "%g - %g i\n", (double)PetscRealPart(values[i]), -(double)PetscImaginaryPart(values[i])));
349:           } else {
350:             PetscCall(PetscViewerASCIIPrintf(viewer, "%g\n", (double)PetscRealPart(values[i])));
351:           }
352: #else
353:           PetscCall(PetscViewerASCIIPrintf(viewer, "%g\n", (double)values[i]));
354: #endif
355:         }
356:       }
357:     }
358:     PetscCall(PetscFree(values));
359:   } else {
360:     if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
361:       /* Rank 0 is not trying to receive anything, so don't send anything */
362:     } else {
363:       if (format == PETSC_VIEWER_ASCII_MATLAB || format == PETSC_VIEWER_ASCII_VTK_DEPRECATED || format == PETSC_VIEWER_ASCII_VTK_CELL_DEPRECATED) {
364:         /* this may be a collective operation so make sure everyone calls it */
365:         PetscCall(PetscObjectGetName((PetscObject)xin, &name));
366:       }
367:       PetscCallMPI(MPI_Send((void *)xarray, xin->map->n, MPIU_SCALAR, 0, tag, PetscObjectComm((PetscObject)xin)));
368:     }
369:   }
370:   PetscCall(PetscViewerFlush(viewer));
371:   PetscCall(VecRestoreArrayRead(xin, &xarray));
372:   PetscFunctionReturn(PETSC_SUCCESS);
373: }

375: PetscErrorCode VecView_MPI_Binary(Vec xin, PetscViewer viewer)
376: {
377:   return VecView_Binary(xin, viewer);
378: }

380: #include <petscdraw.h>
381: PetscErrorCode VecView_MPI_Draw_LG(Vec xin, PetscViewer viewer)
382: {
383:   PetscDraw          draw;
384:   PetscBool          isnull;
385:   PetscDrawLG        lg;
386:   PetscMPIInt        i, size, rank, n, N, *lens = NULL, *disp = NULL;
387:   PetscReal         *values, *xx = NULL, *yy = NULL;
388:   const PetscScalar *xarray;
389:   int                colors[] = {PETSC_DRAW_RED};

391:   PetscFunctionBegin;
392:   PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
393:   PetscCall(PetscDrawIsNull(draw, &isnull));
394:   if (isnull) PetscFunctionReturn(PETSC_SUCCESS);
395:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)xin), &rank));
396:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)xin), &size));
397:   PetscCall(PetscMPIIntCast(xin->map->n, &n));
398:   PetscCall(PetscMPIIntCast(xin->map->N, &N));

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

421:   PetscCall(PetscViewerDrawGetDrawLG(viewer, 0, &lg));
422:   PetscCall(PetscDrawLGReset(lg));
423:   PetscCall(PetscDrawLGSetDimension(lg, 1));
424:   PetscCall(PetscDrawLGSetColors(lg, colors));
425:   if (rank == 0) {
426:     PetscCall(PetscDrawLGAddPoints(lg, N, &xx, &yy));
427:     PetscCall(PetscFree2(xx, yy));
428:   }
429:   PetscCall(PetscDrawLGDraw(lg));
430:   PetscCall(PetscDrawLGSave(lg));
431:   PetscFunctionReturn(PETSC_SUCCESS);
432: }

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

445:   PetscFunctionBegin;
446:   PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
447:   PetscCall(PetscDrawIsNull(draw, &isnull));
448:   if (isnull) PetscFunctionReturn(PETSC_SUCCESS);
449:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)xin), &size));
450:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)xin), &rank));

452:   PetscCall(VecMin(xin, NULL, &min));
453:   PetscCall(VecMax(xin, NULL, &max));
454:   if (min == max) {
455:     min -= 1.e-5;
456:     max += 1.e-5;
457:   }

459:   PetscCall(PetscDrawCheckResizedWindow(draw));
460:   PetscCall(PetscDrawClear(draw));

462:   PetscCall(PetscDrawAxisCreate(draw, &axis));
463:   PetscCall(PetscDrawAxisSetLimits(axis, 0.0, (PetscReal)xin->map->N, min, max));
464:   PetscCall(PetscDrawAxisDraw(axis));
465:   PetscCall(PetscDrawAxisDestroy(&axis));

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

482:   PetscCall(PetscDrawFlush(draw));
483:   PetscCall(PetscDrawPause(draw));
484:   PetscCall(PetscDrawSave(draw));
485:   PetscFunctionReturn(PETSC_SUCCESS);
486: }

488: #if defined(PETSC_HAVE_MATLAB)
489: PetscErrorCode VecView_MPI_Matlab(Vec xin, PetscViewer viewer)
490: {
491:   PetscMPIInt        rank, size, *lens;
492:   PetscInt           i, N = xin->map->N;
493:   const PetscScalar *xarray;
494:   PetscScalar       *xx;

496:   PetscFunctionBegin;
497:   PetscCall(VecGetArrayRead(xin, &xarray));
498:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)xin), &rank));
499:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)xin), &size));
500:   if (rank == 0) {
501:     PetscCall(PetscMalloc1(N, &xx));
502:     PetscCall(PetscMalloc1(size, &lens));
503:     for (i = 0; i < size; i++) lens[i] = xin->map->range[i + 1] - xin->map->range[i];

505:     PetscCallMPI(MPI_Gatherv((void *)xarray, xin->map->n, MPIU_SCALAR, xx, lens, xin->map->range, MPIU_SCALAR, 0, PetscObjectComm((PetscObject)xin)));
506:     PetscCall(PetscFree(lens));

508:     PetscCall(PetscObjectName((PetscObject)xin));
509:     PetscCall(PetscViewerMatlabPutArray(viewer, N, 1, xx, ((PetscObject)xin)->name));

511:     PetscCall(PetscFree(xx));
512:   } else {
513:     PetscCallMPI(MPI_Gatherv((void *)xarray, xin->map->n, MPIU_SCALAR, 0, 0, 0, MPIU_SCALAR, 0, PetscObjectComm((PetscObject)xin)));
514:   }
515:   PetscCall(VecRestoreArrayRead(xin, &xarray));
516:   PetscFunctionReturn(PETSC_SUCCESS);
517: }
518: #endif

520: #if defined(PETSC_HAVE_ADIOS)
521:   #include <adios.h>
522:   #include <adios_read.h>
523: #include <petsc/private/vieweradiosimpl.h>
524: #include <petsc/private/viewerimpl.h>

526: PetscErrorCode VecView_MPI_ADIOS(Vec xin, PetscViewer viewer)
527: {
528:   PetscViewer_ADIOS *adios = (PetscViewer_ADIOS *)viewer->data;
529:   const char        *vecname;
530:   int64_t            id;
531:   PetscInt           n, N, rstart;
532:   const PetscScalar *array;
533:   char               nglobalname[16], nlocalname[16], coffset[16];

535:   PetscFunctionBegin;
536:   PetscCall(PetscObjectGetName((PetscObject)xin, &vecname));

538:   PetscCall(VecGetLocalSize(xin, &n));
539:   PetscCall(VecGetSize(xin, &N));
540:   PetscCall(VecGetOwnershipRange(xin, &rstart, NULL));

542:   PetscCall(PetscSNPrintf(nlocalname, PETSC_STATIC_ARRAY_LENGTH(nlocalname), "%" PetscInt_FMT, n));
543:   PetscCall(PetscSNPrintf(nglobalname, PETSC_STATIC_ARRAY_LENGTH(nglobalname), "%" PetscInt_FMT, N));
544:   PetscCall(PetscSNPrintf(coffset, PETSC_STATIC_ARRAY_LENGTH(coffset), "%" PetscInt_FMT, rstart));
545:   id = adios_define_var(Petsc_adios_group, vecname, "", adios_double, nlocalname, nglobalname, coffset);
546:   PetscCall(VecGetArrayRead(xin, &array));
547:   PetscCallExternal(adios_write_byid, adios->adios_handle, id, array);
548:   PetscCall(VecRestoreArrayRead(xin, &array));
549:   PetscFunctionReturn(PETSC_SUCCESS);
550: }
551: #endif

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

575:   PetscFunctionBegin;
576:   PetscCall(PetscViewerHDF5OpenGroup(viewer, NULL, &file_id, &group));
577:   PetscCall(PetscViewerHDF5IsTimestepping(viewer, &timestepping));
578:   if (timestepping) PetscCall(PetscViewerHDF5GetTimestep(viewer, &timestep));
579:   PetscCall(PetscViewerHDF5GetBaseDimension2(viewer, &dim2));
580:   PetscCall(PetscViewerHDF5GetSPOutput(viewer, &spoutput));

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

603:   maxDims[dim]   = dims[dim];
604:   chunkDims[dim] = PetscMax(1, dims[dim]);
605:   chunksize *= chunkDims[dim];
606:   ++dim;
607:   if (bs > 1 || dim2) {
608:     dims[dim]      = bs;
609:     maxDims[dim]   = dims[dim];
610:     chunkDims[dim] = PetscMax(1, dims[dim]);
611:     chunksize *= chunkDims[dim];
612:     ++dim;
613:   }
614:   #if defined(PETSC_USE_COMPLEX)
615:   dims[dim]      = 2;
616:   maxDims[dim]   = dims[dim];
617:   chunkDims[dim] = PetscMax(1, dims[dim]);
618:   chunksize *= chunkDims[dim];
619:   /* hdf5 chunks must be less than 4GB */
620:   if (chunksize > PETSC_HDF5_MAX_CHUNKSIZE / 64) {
621:     if (bs > 1 || dim2) {
622:       if (chunkDims[dim - 2] > (hsize_t)PetscSqrtReal((PetscReal)(PETSC_HDF5_MAX_CHUNKSIZE / 128))) chunkDims[dim - 2] = (hsize_t)PetscSqrtReal((PetscReal)(PETSC_HDF5_MAX_CHUNKSIZE / 128));
623:       if (chunkDims[dim - 1] > (hsize_t)PetscSqrtReal((PetscReal)(PETSC_HDF5_MAX_CHUNKSIZE / 128))) chunkDims[dim - 1] = (hsize_t)PetscSqrtReal((PetscReal)(PETSC_HDF5_MAX_CHUNKSIZE / 128));
624:     } else {
625:       chunkDims[dim - 1] = PETSC_HDF5_MAX_CHUNKSIZE / 128;
626:     }
627:   }
628:   ++dim;
629:   #else
630:   /* hdf5 chunks must be less than 4GB */
631:   if (chunksize > PETSC_HDF5_MAX_CHUNKSIZE / 64) {
632:     if (bs > 1 || dim2) {
633:       if (chunkDims[dim - 2] > (hsize_t)PetscSqrtReal((PetscReal)(PETSC_HDF5_MAX_CHUNKSIZE / 64))) chunkDims[dim - 2] = (hsize_t)PetscSqrtReal((PetscReal)(PETSC_HDF5_MAX_CHUNKSIZE / 64));
634:       if (chunkDims[dim - 1] > (hsize_t)PetscSqrtReal((PetscReal)(PETSC_HDF5_MAX_CHUNKSIZE / 64))) chunkDims[dim - 1] = (hsize_t)PetscSqrtReal((PetscReal)(PETSC_HDF5_MAX_CHUNKSIZE / 64));
635:     } else {
636:       chunkDims[dim - 1] = PETSC_HDF5_MAX_CHUNKSIZE / 64;
637:     }
638:   }
639:   #endif

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

643:   #if defined(PETSC_USE_REAL_SINGLE)
644:   memscalartype  = H5T_NATIVE_FLOAT;
645:   filescalartype = H5T_NATIVE_FLOAT;
646:   #elif defined(PETSC_USE_REAL___FLOAT128)
647:     #error "HDF5 output with 128 bit floats not supported."
648:   #elif defined(PETSC_USE_REAL___FP16)
649:     #error "HDF5 output with 16 bit floats not supported."
650:   #else
651:   memscalartype = H5T_NATIVE_DOUBLE;
652:   if (spoutput == PETSC_TRUE) filescalartype = H5T_NATIVE_FLOAT;
653:   else filescalartype = H5T_NATIVE_DOUBLE;
654:   #endif

656:   /* Create the dataset with default properties and close filespace */
657:   PetscCall(PetscObjectGetName((PetscObject)xin, &vecname));
658:   if (H5Lexists(group, vecname, H5P_DEFAULT) < 1) {
659:     /* Create chunk */
660:     PetscCallHDF5Return(chunkspace, H5Pcreate, (H5P_DATASET_CREATE));
661:     PetscCallHDF5(H5Pset_chunk, (chunkspace, dim, chunkDims));

663:     PetscCallHDF5Return(dset_id, H5Dcreate2, (group, vecname, filescalartype, filespace, H5P_DEFAULT, chunkspace, H5P_DEFAULT));
664:     PetscCallHDF5(H5Pclose, (chunkspace));
665:   } else {
666:     PetscCallHDF5Return(dset_id, H5Dopen2, (group, vecname, H5P_DEFAULT));
667:     PetscCallHDF5(H5Dset_extent, (dset_id, dims));
668:   }
669:   PetscCallHDF5(H5Sclose, (filespace));

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

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

719:   PetscCall(VecGetArrayRead(xin, &x));
720:   PetscCallHDF5(H5Dwrite, (dset_id, memscalartype, memspace, filespace, hdf5->dxpl_id, x));
721:   PetscCallHDF5(H5Fflush, (file_id, H5F_SCOPE_GLOBAL));
722:   PetscCall(VecRestoreArrayRead(xin, &x));

724:   /* Close/release resources */
725:   PetscCallHDF5(H5Gclose, (group));
726:   PetscCallHDF5(H5Sclose, (filespace));
727:   PetscCallHDF5(H5Sclose, (memspace));
728:   PetscCallHDF5(H5Dclose, (dset_id));

730:   #if defined(PETSC_USE_COMPLEX)
731:   {
732:     PetscBool tru = PETSC_TRUE;
733:     PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)xin, "complex", PETSC_BOOL, &tru));
734:   }
735:   #endif
736:   if (timestepping) PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)xin, "timestepping", PETSC_BOOL, &timestepping));
737:   PetscCall(PetscInfo(xin, "Wrote Vec object with name %s\n", vecname));
738:   PetscFunctionReturn(PETSC_SUCCESS);
739: }
740: #endif

742: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec xin, PetscViewer viewer)
743: {
744:   PetscBool iascii, isbinary, isdraw;
745: #if defined(PETSC_HAVE_MATHEMATICA)
746:   PetscBool ismathematica;
747: #endif
748: #if defined(PETSC_HAVE_HDF5)
749:   PetscBool ishdf5;
750: #endif
751: #if defined(PETSC_HAVE_MATLAB)
752:   PetscBool ismatlab;
753: #endif
754: #if defined(PETSC_HAVE_ADIOS)
755:   PetscBool isadios;
756: #endif
757:   PetscBool isglvis;

759:   PetscFunctionBegin;
760:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
761:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
762:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
763: #if defined(PETSC_HAVE_MATHEMATICA)
764:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERMATHEMATICA, &ismathematica));
765: #endif
766: #if defined(PETSC_HAVE_HDF5)
767:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
768: #endif
769: #if defined(PETSC_HAVE_MATLAB)
770:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERMATLAB, &ismatlab));
771: #endif
772:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERGLVIS, &isglvis));
773: #if defined(PETSC_HAVE_ADIOS)
774:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERADIOS, &isadios));
775: #endif
776:   if (iascii) {
777:     PetscCall(VecView_MPI_ASCII(xin, viewer));
778:   } else if (isbinary) {
779:     PetscCall(VecView_MPI_Binary(xin, viewer));
780:   } else if (isdraw) {
781:     PetscViewerFormat format;
782:     PetscCall(PetscViewerGetFormat(viewer, &format));
783:     if (format == PETSC_VIEWER_DRAW_LG) {
784:       PetscCall(VecView_MPI_Draw_LG(xin, viewer));
785:     } else {
786:       PetscCall(VecView_MPI_Draw(xin, viewer));
787:     }
788: #if defined(PETSC_HAVE_MATHEMATICA)
789:   } else if (ismathematica) {
790:     PetscCall(PetscViewerMathematicaPutVector(viewer, xin));
791: #endif
792: #if defined(PETSC_HAVE_HDF5)
793:   } else if (ishdf5) {
794:     PetscCall(VecView_MPI_HDF5(xin, viewer));
795: #endif
796: #if defined(PETSC_HAVE_ADIOS)
797:   } else if (isadios) {
798:     PetscCall(VecView_MPI_ADIOS(xin, viewer));
799: #endif
800: #if defined(PETSC_HAVE_MATLAB)
801:   } else if (ismatlab) {
802:     PetscCall(VecView_MPI_Matlab(xin, viewer));
803: #endif
804:   } else if (isglvis) PetscCall(VecView_GLVis(xin, viewer));
805:   PetscFunctionReturn(PETSC_SUCCESS);
806: }

808: PetscErrorCode VecGetSize_MPI(Vec xin, PetscInt *N)
809: {
810:   PetscFunctionBegin;
811:   *N = xin->map->N;
812:   PetscFunctionReturn(PETSC_SUCCESS);
813: }

815: PetscErrorCode VecGetValues_MPI(Vec xin, PetscInt ni, const PetscInt ix[], PetscScalar y[])
816: {
817:   const PetscScalar *xx;
818:   const PetscInt     start = xin->map->range[xin->stash.rank];

820:   PetscFunctionBegin;
821:   PetscCall(VecGetArrayRead(xin, &xx));
822:   for (PetscInt i = 0; i < ni; i++) {
823:     if (xin->stash.ignorenegidx && ix[i] < 0) continue;
824:     const PetscInt tmp = ix[i] - start;

826:     PetscCheck(tmp >= 0 && tmp < xin->map->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Can only get local values, trying %" PetscInt_FMT, ix[i]);
827:     y[i] = xx[tmp];
828:   }
829:   PetscCall(VecRestoreArrayRead(xin, &xx));
830:   PetscFunctionReturn(PETSC_SUCCESS);
831: }

833: PetscErrorCode VecSetValues_MPI(Vec xin, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode addv)
834: {
835:   const PetscBool   ignorenegidx = xin->stash.ignorenegidx;
836:   const PetscBool   donotstash   = xin->stash.donotstash;
837:   const PetscMPIInt rank         = xin->stash.rank;
838:   const PetscInt   *owners       = xin->map->range;
839:   const PetscInt    start = owners[rank], end = owners[rank + 1];
840:   PetscScalar      *xx;

842:   PetscFunctionBegin;
843:   if (PetscDefined(USE_DEBUG)) {
844:     PetscCheck(xin->stash.insertmode != INSERT_VALUES || addv != ADD_VALUES, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "You have already inserted values; you cannot now add");
845:     PetscCheck(xin->stash.insertmode != ADD_VALUES || addv != INSERT_VALUES, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "You have already added values; you cannot now insert");
846:   }
847:   PetscCall(VecGetArray(xin, &xx));
848:   xin->stash.insertmode = addv;
849:   for (PetscInt i = 0; i < ni; ++i) {
850:     PetscInt row;

852:     if (ignorenegidx && ix[i] < 0) continue;
853:     PetscCheck(ix[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Out of range index value %" PetscInt_FMT " cannot be negative", ix[i]);
854:     if ((row = ix[i]) >= start && row < end) {
855:       if (addv == INSERT_VALUES) {
856:         xx[row - start] = y[i];
857:       } else {
858:         xx[row - start] += y[i];
859:       }
860:     } else if (!donotstash) {
861:       PetscCheck(ix[i] < xin->map->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Out of range index value %" PetscInt_FMT " maximum %" PetscInt_FMT, ix[i], xin->map->N);
862:       PetscCall(VecStashValue_Private(&xin->stash, row, y[i]));
863:     }
864:   }
865:   PetscCall(VecRestoreArray(xin, &xx));
866:   PetscFunctionReturn(PETSC_SUCCESS);
867: }

869: PetscErrorCode VecSetValuesBlocked_MPI(Vec xin, PetscInt ni, const PetscInt ix[], const PetscScalar yin[], InsertMode addv)
870: {
871:   PetscMPIInt  rank   = xin->stash.rank;
872:   PetscInt    *owners = xin->map->range, start = owners[rank];
873:   PetscInt     end = owners[rank + 1], i, row, bs = PetscAbs(xin->map->bs), j;
874:   PetscScalar *xx, *y = (PetscScalar *)yin;

876:   PetscFunctionBegin;
877:   PetscCall(VecGetArray(xin, &xx));
878:   if (PetscDefined(USE_DEBUG)) {
879:     PetscCheck(xin->stash.insertmode != INSERT_VALUES || addv != ADD_VALUES, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "You have already inserted values; you cannot now add");
880:     PetscCheck(xin->stash.insertmode != ADD_VALUES || addv != INSERT_VALUES, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "You have already added values; you cannot now insert");
881:   }
882:   xin->stash.insertmode = addv;

884:   if (addv == INSERT_VALUES) {
885:     for (i = 0; i < ni; i++) {
886:       if ((row = bs * ix[i]) >= start && row < end) {
887:         for (j = 0; j < bs; j++) xx[row - start + j] = y[j];
888:       } else if (!xin->stash.donotstash) {
889:         if (ix[i] < 0) {
890:           y += bs;
891:           continue;
892:         }
893:         PetscCheck(ix[i] < xin->map->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Out of range index value %" PetscInt_FMT " max %" PetscInt_FMT, ix[i], xin->map->N);
894:         PetscCall(VecStashValuesBlocked_Private(&xin->bstash, ix[i], y));
895:       }
896:       y += bs;
897:     }
898:   } else {
899:     for (i = 0; i < ni; i++) {
900:       if ((row = bs * ix[i]) >= start && row < end) {
901:         for (j = 0; j < bs; j++) xx[row - start + j] += y[j];
902:       } else if (!xin->stash.donotstash) {
903:         if (ix[i] < 0) {
904:           y += bs;
905:           continue;
906:         }
907:         PetscCheck(ix[i] <= xin->map->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Out of range index value %" PetscInt_FMT " max %" PetscInt_FMT, ix[i], xin->map->N);
908:         PetscCall(VecStashValuesBlocked_Private(&xin->bstash, ix[i], y));
909:       }
910:       y += bs;
911:     }
912:   }
913:   PetscCall(VecRestoreArray(xin, &xx));
914:   PetscFunctionReturn(PETSC_SUCCESS);
915: }

917: /*
918:    Since nsends or nreceives may be zero we add 1 in certain mallocs
919: to make sure we never malloc an empty one.
920: */
921: PetscErrorCode VecAssemblyBegin_MPI(Vec xin)
922: {
923:   PetscInt   *owners = xin->map->range, *bowners, i, bs, nstash, reallocs;
924:   PetscMPIInt size;
925:   InsertMode  addv;
926:   MPI_Comm    comm;

928:   PetscFunctionBegin;
929:   PetscCall(PetscObjectGetComm((PetscObject)xin, &comm));
930:   if (xin->stash.donotstash) PetscFunctionReturn(PETSC_SUCCESS);

932:   PetscCall(MPIU_Allreduce((PetscEnum *)&xin->stash.insertmode, (PetscEnum *)&addv, 1, MPIU_ENUM, MPI_BOR, comm));
933:   PetscCheck(addv != (ADD_VALUES | INSERT_VALUES), comm, PETSC_ERR_ARG_NOTSAMETYPE, "Some processors inserted values while others added");
934:   xin->stash.insertmode  = addv; /* in case this processor had no cache */
935:   xin->bstash.insertmode = addv; /* Block stash implicitly tracks InsertMode of scalar stash */

937:   PetscCall(VecGetBlockSize(xin, &bs));
938:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)xin), &size));
939:   if (!xin->bstash.bowners && xin->map->bs != -1) {
940:     PetscCall(PetscMalloc1(size + 1, &bowners));
941:     for (i = 0; i < size + 1; i++) bowners[i] = owners[i] / bs;
942:     xin->bstash.bowners = bowners;
943:   } else bowners = xin->bstash.bowners;

945:   PetscCall(VecStashScatterBegin_Private(&xin->stash, owners));
946:   PetscCall(VecStashScatterBegin_Private(&xin->bstash, bowners));
947:   PetscCall(VecStashGetInfo_Private(&xin->stash, &nstash, &reallocs));
948:   PetscCall(PetscInfo(xin, "Stash has %" PetscInt_FMT " entries, uses %" PetscInt_FMT " mallocs.\n", nstash, reallocs));
949:   PetscCall(VecStashGetInfo_Private(&xin->bstash, &nstash, &reallocs));
950:   PetscCall(PetscInfo(xin, "Block-Stash has %" PetscInt_FMT " entries, uses %" PetscInt_FMT " mallocs.\n", nstash, reallocs));
951:   PetscFunctionReturn(PETSC_SUCCESS);
952: }

954: PetscErrorCode VecAssemblyEnd_MPI(Vec vec)
955: {
956:   PetscInt     base, i, j, *row, flg, bs;
957:   PetscMPIInt  n;
958:   PetscScalar *val, *vv, *array, *xarray;

960:   PetscFunctionBegin;
961:   if (!vec->stash.donotstash) {
962:     PetscCall(VecGetArray(vec, &xarray));
963:     PetscCall(VecGetBlockSize(vec, &bs));
964:     base = vec->map->range[vec->stash.rank];

966:     /* Process the stash */
967:     while (1) {
968:       PetscCall(VecStashScatterGetMesg_Private(&vec->stash, &n, &row, &val, &flg));
969:       if (!flg) break;
970:       if (vec->stash.insertmode == ADD_VALUES) {
971:         for (i = 0; i < n; i++) xarray[row[i] - base] += val[i];
972:       } else if (vec->stash.insertmode == INSERT_VALUES) {
973:         for (i = 0; i < n; i++) xarray[row[i] - base] = val[i];
974:       } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_CORRUPT, "Insert mode is not set correctly; corrupted vector");
975:     }
976:     PetscCall(VecStashScatterEnd_Private(&vec->stash));

978:     /* now process the block-stash */
979:     while (1) {
980:       PetscCall(VecStashScatterGetMesg_Private(&vec->bstash, &n, &row, &val, &flg));
981:       if (!flg) break;
982:       for (i = 0; i < n; i++) {
983:         array = xarray + row[i] * bs - base;
984:         vv    = val + i * bs;
985:         if (vec->stash.insertmode == ADD_VALUES) {
986:           for (j = 0; j < bs; j++) array[j] += vv[j];
987:         } else if (vec->stash.insertmode == INSERT_VALUES) {
988:           for (j = 0; j < bs; j++) array[j] = vv[j];
989:         } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_CORRUPT, "Insert mode is not set correctly; corrupted vector");
990:       }
991:     }
992:     PetscCall(VecStashScatterEnd_Private(&vec->bstash));
993:     PetscCall(VecRestoreArray(vec, &xarray));
994:   }
995:   vec->stash.insertmode = NOT_SET_VALUES;
996:   PetscFunctionReturn(PETSC_SUCCESS);
997: }

999: PetscErrorCode VecSetPreallocationCOO_MPI(Vec x, PetscCount coo_n, const PetscInt coo_i[])
1000: {
1001:   PetscInt    m, M, rstart, rend;
1002:   Vec_MPI    *vmpi = (Vec_MPI *)x->data;
1003:   PetscCount  k, p, q, rem; /* Loop variables over coo arrays */
1004:   PetscMPIInt size;
1005:   MPI_Comm    comm;

1007:   PetscFunctionBegin;
1008:   PetscCall(PetscObjectGetComm((PetscObject)x, &comm));
1009:   PetscCallMPI(MPI_Comm_size(comm, &size));
1010:   PetscCall(VecResetPreallocationCOO_MPI(x));

1012:   PetscCall(PetscLayoutSetUp(x->map));
1013:   PetscCall(VecGetOwnershipRange(x, &rstart, &rend));
1014:   PetscCall(VecGetLocalSize(x, &m));
1015:   PetscCall(VecGetSize(x, &M));

1017:   /* ---------------------------------------------------------------------------*/
1018:   /* Sort COOs along with a permutation array, so that negative indices come    */
1019:   /* first, then local ones, then remote ones.                                  */
1020:   /* ---------------------------------------------------------------------------*/
1021:   PetscCount n1 = coo_n, nneg, *perm;
1022:   PetscInt  *i1; /* Copy of input COOs along with a permutation array */
1023:   PetscCall(PetscMalloc1(n1, &i1));
1024:   PetscCall(PetscMalloc1(n1, &perm));
1025:   PetscCall(PetscArraycpy(i1, coo_i, n1)); /* Make a copy since we'll modify it */
1026:   for (k = 0; k < n1; k++) perm[k] = k;

1028:   /* Manipulate i1[] so that entries with negative indices will have the smallest
1029:      index, local entries will have greater but negative indices, and remote entries
1030:      will have positive indices.
1031:   */
1032:   for (k = 0; k < n1; k++) {
1033:     if (i1[k] < 0) {
1034:       if (x->stash.ignorenegidx) i1[k] = PETSC_MIN_INT; /* e.g., -2^31, minimal to move them ahead */
1035:       else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Found a negative index in VecSetPreallocateCOO() but VEC_IGNORE_NEGATIVE_INDICES was not set");
1036:     } else if (i1[k] >= rstart && i1[k] < rend) {
1037:       i1[k] -= PETSC_MAX_INT; /* e.g., minus 2^31-1 to shift local rows to range of [-PETSC_MAX_INT, -1] */
1038:     } else {
1039:       PetscCheck(i1[k] < M, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Found index %" PetscInt_FMT " in VecSetPreallocateCOO() larger than the global size %" PetscInt_FMT "", i1[k], M);
1040:       if (x->stash.donotstash) i1[k] = PETSC_MIN_INT; /* Ignore off-proc indices as if they were negative */
1041:     }
1042:   }

1044:   /* Sort the indices, after that, [0,nneg) have ignored entries, [nneg,rem) have local entries and [rem,n1) have remote entries */
1045:   PetscCall(PetscSortIntWithCountArray(n1, i1, perm));
1046:   for (k = 0; k < n1; k++) {
1047:     if (i1[k] > PETSC_MIN_INT) break;
1048:   } /* Advance k to the first entry we need to take care of */
1049:   nneg = k;
1050:   PetscCall(PetscSortedIntUpperBound(i1, nneg, n1, rend - 1 - PETSC_MAX_INT, &rem)); /* rem is upper bound of the last local row */
1051:   for (k = nneg; k < rem; k++) i1[k] += PETSC_MAX_INT;                               /* Revert indices of local entries */

1053:   /* ---------------------------------------------------------------------------*/
1054:   /*           Build stuff for local entries                                    */
1055:   /* ---------------------------------------------------------------------------*/
1056:   PetscCount tot1, *jmap1, *perm1;
1057:   PetscCall(PetscCalloc1(m + 1, &jmap1));
1058:   for (k = nneg; k < rem; k++) jmap1[i1[k] - rstart + 1]++; /* Count repeats of each local entry */
1059:   for (k = 0; k < m; k++) jmap1[k + 1] += jmap1[k];         /* Transform jmap1[] to CSR-like data structure */
1060:   tot1 = jmap1[m];
1061:   PetscAssert(tot1 == rem - nneg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected errors in VecSetPreallocationCOO_MPI");
1062:   PetscCall(PetscMalloc1(tot1, &perm1));
1063:   PetscCall(PetscArraycpy(perm1, perm + nneg, tot1));

1065:   /* ---------------------------------------------------------------------------*/
1066:   /*        Record the permutation array for filling the send buffer            */
1067:   /* ---------------------------------------------------------------------------*/
1068:   PetscCount *Cperm;
1069:   PetscCall(PetscMalloc1(n1 - rem, &Cperm));
1070:   PetscCall(PetscArraycpy(Cperm, perm + rem, n1 - rem));
1071:   PetscCall(PetscFree(perm));

1073:   /* ---------------------------------------------------------------------------*/
1074:   /*           Send remote entries to their owner                                  */
1075:   /* ---------------------------------------------------------------------------*/
1076:   /* Find which entries should be sent to which remote ranks*/
1077:   PetscInt        nsend = 0; /* Number of MPI ranks to send data to */
1078:   PetscMPIInt    *sendto;    /* [nsend], storing remote ranks */
1079:   PetscInt       *nentries;  /* [nsend], storing number of entries sent to remote ranks; Assume PetscInt is big enough for this count, and error if not */
1080:   const PetscInt *ranges;
1081:   PetscInt        maxNsend = size >= 128 ? 128 : size; /* Assume max 128 neighbors; realloc when needed */

1083:   PetscCall(PetscLayoutGetRanges(x->map, &ranges));
1084:   PetscCall(PetscMalloc2(maxNsend, &sendto, maxNsend, &nentries));
1085:   for (k = rem; k < n1;) {
1086:     PetscMPIInt owner;
1087:     PetscInt    firstRow, lastRow;

1089:     /* Locate a row range */
1090:     firstRow = i1[k]; /* first row of this owner */
1091:     PetscCall(PetscLayoutFindOwner(x->map, firstRow, &owner));
1092:     lastRow = ranges[owner + 1] - 1; /* last row of this owner */

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

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

1103:       PetscCall(PetscMalloc2(maxNsend2, &sendto2, maxNsend2, &nentries2));
1104:       PetscCall(PetscArraycpy(sendto2, sendto, maxNsend));
1105:       PetscCall(PetscArraycpy(nentries2, nentries2, maxNsend + 1));
1106:       PetscCall(PetscFree2(sendto, nentries2));
1107:       sendto   = sendto2;
1108:       nentries = nentries2;
1109:       maxNsend = maxNsend2;
1110:     }
1111:     sendto[nsend]   = owner;
1112:     nentries[nsend] = p - k;
1113:     PetscCall(PetscCountCast(p - k, &nentries[nsend]));
1114:     nsend++;
1115:     k = p;
1116:   }

1118:   /* Build 1st SF to know offsets on remote to send data */
1119:   PetscSF      sf1;
1120:   PetscInt     nroots = 1, nroots2 = 0;
1121:   PetscInt     nleaves = nsend, nleaves2 = 0;
1122:   PetscInt    *offsets;
1123:   PetscSFNode *iremote;

1125:   PetscCall(PetscSFCreate(comm, &sf1));
1126:   PetscCall(PetscMalloc1(nsend, &iremote));
1127:   PetscCall(PetscMalloc1(nsend, &offsets));
1128:   for (k = 0; k < nsend; k++) {
1129:     iremote[k].rank  = sendto[k];
1130:     iremote[k].index = 0;
1131:     nleaves2 += nentries[k];
1132:     PetscCheck(nleaves2 >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of SF leaves is too large for PetscInt");
1133:   }
1134:   PetscCall(PetscSFSetGraph(sf1, nroots, nleaves, NULL, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
1135:   PetscCall(PetscSFFetchAndOpWithMemTypeBegin(sf1, MPIU_INT, PETSC_MEMTYPE_HOST, &nroots2 /*rootdata*/, PETSC_MEMTYPE_HOST, nentries /*leafdata*/, PETSC_MEMTYPE_HOST, offsets /*leafupdate*/, MPI_SUM));
1136:   PetscCall(PetscSFFetchAndOpEnd(sf1, MPIU_INT, &nroots2, nentries, offsets, MPI_SUM)); /* Would nroots2 overflow, we check offsets[] below */
1137:   PetscCall(PetscSFDestroy(&sf1));
1138:   PetscAssert(nleaves2 == n1 - rem, PETSC_COMM_SELF, PETSC_ERR_PLIB, "nleaves2 %" PetscInt_FMT " != number of remote entries %" PetscCount_FMT "", nleaves2, n1 - rem);

1140:   /* Build 2nd SF to send remote COOs to their owner */
1141:   PetscSF sf2;
1142:   nroots  = nroots2;
1143:   nleaves = nleaves2;
1144:   PetscCall(PetscSFCreate(comm, &sf2));
1145:   PetscCall(PetscSFSetFromOptions(sf2));
1146:   PetscCall(PetscMalloc1(nleaves, &iremote));
1147:   p = 0;
1148:   for (k = 0; k < nsend; k++) {
1149:     PetscCheck(offsets[k] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of SF roots is too large for PetscInt");
1150:     for (q = 0; q < nentries[k]; q++, p++) {
1151:       iremote[p].rank  = sendto[k];
1152:       iremote[p].index = offsets[k] + q;
1153:     }
1154:   }
1155:   PetscCall(PetscSFSetGraph(sf2, nroots, nleaves, NULL, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));

1157:   /* Send the remote COOs to their owner */
1158:   PetscInt    n2 = nroots, *i2; /* Buffers for received COOs from other ranks, along with a permutation array */
1159:   PetscCount *perm2;
1160:   PetscCall(PetscMalloc1(n2, &i2));
1161:   PetscCall(PetscMalloc1(n2, &perm2));
1162:   PetscCall(PetscSFReduceWithMemTypeBegin(sf2, MPIU_INT, PETSC_MEMTYPE_HOST, i1 + rem, PETSC_MEMTYPE_HOST, i2, MPI_REPLACE));
1163:   PetscCall(PetscSFReduceEnd(sf2, MPIU_INT, i1 + rem, i2, MPI_REPLACE));

1165:   PetscCall(PetscFree(i1));
1166:   PetscCall(PetscFree(offsets));
1167:   PetscCall(PetscFree2(sendto, nentries));

1169:   /* ---------------------------------------------------------------*/
1170:   /* Sort received COOs along with a permutation array            */
1171:   /* ---------------------------------------------------------------*/
1172:   PetscCount  *imap2;
1173:   PetscCount  *jmap2, nnz2;
1174:   PetscScalar *sendbuf, *recvbuf;
1175:   PetscInt     old;
1176:   PetscCount   sendlen = n1 - rem, recvlen = n2;

1178:   for (k = 0; k < n2; k++) perm2[k] = k;
1179:   PetscCall(PetscSortIntWithCountArray(n2, i2, perm2));

1181:   /* nnz2 will be # of unique entries in the recvbuf */
1182:   nnz2 = n2;
1183:   for (k = 1; k < n2; k++) {
1184:     if (i2[k] == i2[k - 1]) nnz2--;
1185:   }

1187:   /* Build imap2[] and jmap2[] for each unique entry */
1188:   PetscCall(PetscMalloc4(nnz2, &imap2, nnz2 + 1, &jmap2, sendlen, &sendbuf, recvlen, &recvbuf));
1189:   p        = -1;
1190:   old      = -1;
1191:   jmap2[0] = 0;
1192:   jmap2++;
1193:   for (k = 0; k < n2; k++) {
1194:     if (i2[k] != old) { /* Meet a new entry */
1195:       p++;
1196:       imap2[p] = i2[k] - rstart;
1197:       jmap2[p] = 1;
1198:       old      = i2[k];
1199:     } else {
1200:       jmap2[p]++;
1201:     }
1202:   }
1203:   jmap2--;
1204:   for (k = 0; k < nnz2; k++) jmap2[k + 1] += jmap2[k];

1206:   PetscCall(PetscFree(i2));

1208:   vmpi->coo_n = coo_n;
1209:   vmpi->tot1  = tot1;
1210:   vmpi->jmap1 = jmap1;
1211:   vmpi->perm1 = perm1;
1212:   vmpi->nnz2  = nnz2;
1213:   vmpi->imap2 = imap2;
1214:   vmpi->jmap2 = jmap2;
1215:   vmpi->perm2 = perm2;

1217:   vmpi->Cperm   = Cperm;
1218:   vmpi->sendbuf = sendbuf;
1219:   vmpi->recvbuf = recvbuf;
1220:   vmpi->sendlen = sendlen;
1221:   vmpi->recvlen = recvlen;
1222:   vmpi->coo_sf  = sf2;
1223:   PetscFunctionReturn(PETSC_SUCCESS);
1224: }

1226: PetscErrorCode VecSetValuesCOO_MPI(Vec x, const PetscScalar v[], InsertMode imode)
1227: {
1228:   Vec_MPI          *vmpi = (Vec_MPI *)x->data;
1229:   PetscInt          m;
1230:   PetscScalar      *a, *sendbuf = vmpi->sendbuf, *recvbuf = vmpi->recvbuf;
1231:   const PetscCount *jmap1 = vmpi->jmap1;
1232:   const PetscCount *perm1 = vmpi->perm1;
1233:   const PetscCount *imap2 = vmpi->imap2;
1234:   const PetscCount *jmap2 = vmpi->jmap2;
1235:   const PetscCount *perm2 = vmpi->perm2;
1236:   const PetscCount *Cperm = vmpi->Cperm;
1237:   const PetscCount  nnz2  = vmpi->nnz2;

1239:   PetscFunctionBegin;
1240:   PetscCall(VecGetLocalSize(x, &m));
1241:   PetscCall(VecGetArray(x, &a));

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

1246:   /* Send remote entries to their owner and overlap the communication with local computation */
1247:   PetscCall(PetscSFReduceWithMemTypeBegin(vmpi->coo_sf, MPIU_SCALAR, PETSC_MEMTYPE_HOST, sendbuf, PETSC_MEMTYPE_HOST, recvbuf, MPI_REPLACE));
1248:   /* Add local entries to A and B */
1249:   for (PetscInt i = 0; i < m; i++) { /* All entries in a[] are either zero'ed or added with a value (i.e., initialized) */
1250:     PetscScalar sum = 0.0;           /* Do partial summation first to improve numerical stability */
1251:     for (PetscCount k = jmap1[i]; k < jmap1[i + 1]; k++) sum += v[perm1[k]];
1252:     a[i] = (imode == INSERT_VALUES ? 0.0 : a[i]) + sum;
1253:   }
1254:   PetscCall(PetscSFReduceEnd(vmpi->coo_sf, MPIU_SCALAR, sendbuf, recvbuf, MPI_REPLACE));

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

1261:   PetscCall(VecRestoreArray(x, &a));
1262:   PetscFunctionReturn(PETSC_SUCCESS);
1263: }