Actual source code: pdvec.c

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

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

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

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

 30:   PetscFunctionBegin;
 31:   PetscCall(PetscLogObjectState((PetscObject)v, "Length=%" PetscInt_FMT, v->map->N));
 32:   if (!x) PetscFunctionReturn(PETSC_SUCCESS);
 33:   PetscCall(PetscFree(x->array_allocated));

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

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

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

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

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

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

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

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

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

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

227:       PetscCall(VecGetLocalSize(xin, &nLen));
228:       PetscCall(PetscMPIIntCast(nLen, &n));
229:       PetscCall(VecGetBlockSize(xin, &bs));
230:       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);

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

255:       PetscCall(VecGetLocalSize(xin, &nLen));
256:       PetscCall(PetscMPIIntCast(nLen, &n));
257:       PetscCall(VecGetBlockSize(xin, &bs));
258:       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);

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

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

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

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

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

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

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

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

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

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

457:   PetscCall(PetscDrawCheckResizedWindow(draw));
458:   PetscCall(PetscDrawClear(draw));

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

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

480:   PetscCall(PetscDrawFlush(draw));
481:   PetscCall(PetscDrawPause(draw));
482:   PetscCall(PetscDrawSave(draw));
483:   PetscFunctionReturn(PETSC_SUCCESS);
484: }

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

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

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

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

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

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

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

533:   PetscFunctionBegin;
534:   PetscCall(PetscObjectGetName((PetscObject)xin, &vecname));

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

840:   PetscFunctionBegin;
841:   if (PetscDefined(USE_DEBUG)) {
842:     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");
843:     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");
844:   }
845:   PetscCall(VecGetArray(xin, &xx));
846:   xin->stash.insertmode = addv;
847:   for (PetscInt i = 0; i < ni; ++i) {
848:     PetscInt row;

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

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

874:   PetscFunctionBegin;
875:   PetscCall(VecGetArray(xin, &xx));
876:   if (PetscDefined(USE_DEBUG)) {
877:     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");
878:     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");
879:   }
880:   xin->stash.insertmode = addv;

882:   if (addv == INSERT_VALUES) {
883:     for (i = 0; i < ni; i++) {
884:       if ((row = bs * ix[i]) >= start && row < end) {
885:         for (j = 0; j < bs; j++) xx[row - start + j] = y[j];
886:       } else if (!xin->stash.donotstash) {
887:         if (ix[i] < 0) {
888:           y += bs;
889:           continue;
890:         }
891:         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);
892:         PetscCall(VecStashValuesBlocked_Private(&xin->bstash, ix[i], y));
893:       }
894:       y += bs;
895:     }
896:   } else {
897:     for (i = 0; i < ni; i++) {
898:       if ((row = bs * ix[i]) >= start && row < end) {
899:         for (j = 0; j < bs; j++) xx[row - start + j] += y[j];
900:       } else if (!xin->stash.donotstash) {
901:         if (ix[i] < 0) {
902:           y += bs;
903:           continue;
904:         }
905:         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);
906:         PetscCall(VecStashValuesBlocked_Private(&xin->bstash, ix[i], y));
907:       }
908:       y += bs;
909:     }
910:   }
911:   PetscCall(VecRestoreArray(xin, &xx));
912:   PetscFunctionReturn(PETSC_SUCCESS);
913: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1163:   PetscCall(PetscFree(i1));
1164:   PetscCall(PetscFree(offsets));
1165:   PetscCall(PetscFree2(sendto, nentries));

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

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

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

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

1204:   PetscCall(PetscFree(i2));

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

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

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

1237:   PetscFunctionBegin;
1238:   PetscCall(VecGetLocalSize(x, &m));
1239:   PetscCall(VecGetArray(x, &a));

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

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

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

1259:   PetscCall(VecRestoreArray(x, &a));
1260:   PetscFunctionReturn(PETSC_SUCCESS);
1261: }