Actual source code: bvec2.c


  2: /*
  3:    Implements the sequential vectors.
  4: */

  6: #include <../src/vec/vec/impls/dvecimpl.h>
  7: #include <../src/vec/vec/impls/mpi/pvecimpl.h>
  8: #include <petsc/private/glvisviewerimpl.h>
  9: #include <petsc/private/glvisvecimpl.h>
 10: #include <petscblaslapack.h>

 12: #if defined(PETSC_HAVE_HDF5)
 13: extern PetscErrorCode VecView_MPI_HDF5(Vec, PetscViewer);
 14: #endif

 16: static PetscErrorCode VecPointwiseApply_Seq(Vec win, Vec xin, Vec yin, PetscScalar (*const func)(PetscScalar, PetscScalar))
 17: {
 18:   const PetscInt n = win->map->n;
 19:   PetscScalar   *ww, *xx, *yy; /* cannot make xx or yy const since might be ww */

 21:   PetscFunctionBegin;
 22:   PetscCall(VecGetArrayRead(xin, (const PetscScalar **)&xx));
 23:   PetscCall(VecGetArrayRead(yin, (const PetscScalar **)&yy));
 24:   PetscCall(VecGetArray(win, &ww));
 25:   for (PetscInt i = 0; i < n; ++i) ww[i] = func(xx[i], yy[i]);
 26:   PetscCall(VecRestoreArrayRead(xin, (const PetscScalar **)&xx));
 27:   PetscCall(VecRestoreArrayRead(yin, (const PetscScalar **)&yy));
 28:   PetscCall(VecRestoreArray(win, &ww));
 29:   PetscCall(PetscLogFlops(n));
 30:   PetscFunctionReturn(PETSC_SUCCESS);
 31: }

 33: static PetscScalar MaxRealPart(PetscScalar x, PetscScalar y)
 34: {
 35:   // use temporaries to avoid reevaluating side-effects
 36:   const PetscReal rx = PetscRealPart(x), ry = PetscRealPart(y);

 38:   return PetscMax(rx, ry);
 39: }

 41: PetscErrorCode VecPointwiseMax_Seq(Vec win, Vec xin, Vec yin)
 42: {
 43:   PetscFunctionBegin;
 44:   PetscCall(VecPointwiseApply_Seq(win, xin, yin, MaxRealPart));
 45:   PetscFunctionReturn(PETSC_SUCCESS);
 46: }

 48: static PetscScalar MinRealPart(PetscScalar x, PetscScalar y)
 49: {
 50:   // use temporaries to avoid reevaluating side-effects
 51:   const PetscReal rx = PetscRealPart(x), ry = PetscRealPart(y);

 53:   return PetscMin(rx, ry);
 54: }

 56: PetscErrorCode VecPointwiseMin_Seq(Vec win, Vec xin, Vec yin)
 57: {
 58:   PetscFunctionBegin;
 59:   PetscCall(VecPointwiseApply_Seq(win, xin, yin, MinRealPart));
 60:   PetscFunctionReturn(PETSC_SUCCESS);
 61: }

 63: static PetscScalar MaxAbs(PetscScalar x, PetscScalar y)
 64: {
 65:   return (PetscScalar)PetscMax(PetscAbsScalar(x), PetscAbsScalar(y));
 66: }

 68: PetscErrorCode VecPointwiseMaxAbs_Seq(Vec win, Vec xin, Vec yin)
 69: {
 70:   PetscFunctionBegin;
 71:   PetscCall(VecPointwiseApply_Seq(win, xin, yin, MaxAbs));
 72:   PetscFunctionReturn(PETSC_SUCCESS);
 73: }

 75: #include <../src/vec/vec/impls/seq/ftn-kernels/fxtimesy.h>

 77: PetscErrorCode VecPointwiseMult_Seq(Vec win, Vec xin, Vec yin)
 78: {
 79:   PetscInt     n = win->map->n, i;
 80:   PetscScalar *ww, *xx, *yy; /* cannot make xx or yy const since might be ww */

 82:   PetscFunctionBegin;
 83:   PetscCall(VecGetArrayRead(xin, (const PetscScalar **)&xx));
 84:   PetscCall(VecGetArrayRead(yin, (const PetscScalar **)&yy));
 85:   PetscCall(VecGetArray(win, &ww));
 86:   if (ww == xx) {
 87:     for (i = 0; i < n; i++) ww[i] *= yy[i];
 88:   } else if (ww == yy) {
 89:     for (i = 0; i < n; i++) ww[i] *= xx[i];
 90:   } else {
 91: #if defined(PETSC_USE_FORTRAN_KERNEL_XTIMESY)
 92:     fortranxtimesy_(xx, yy, ww, &n);
 93: #else
 94:     for (i = 0; i < n; i++) ww[i] = xx[i] * yy[i];
 95: #endif
 96:   }
 97:   PetscCall(VecRestoreArrayRead(xin, (const PetscScalar **)&xx));
 98:   PetscCall(VecRestoreArrayRead(yin, (const PetscScalar **)&yy));
 99:   PetscCall(VecRestoreArray(win, &ww));
100:   PetscCall(PetscLogFlops(n));
101:   PetscFunctionReturn(PETSC_SUCCESS);
102: }

104: static PetscScalar ScalDiv(PetscScalar x, PetscScalar y)
105: {
106:   return y == 0.0 ? 0.0 : x / y;
107: }

109: PetscErrorCode VecPointwiseDivide_Seq(Vec win, Vec xin, Vec yin)
110: {
111:   PetscFunctionBegin;
112:   PetscCall(VecPointwiseApply_Seq(win, xin, yin, ScalDiv));
113:   PetscFunctionReturn(PETSC_SUCCESS);
114: }

116: PetscErrorCode VecSetRandom_Seq(Vec xin, PetscRandom r)
117: {
118:   PetscScalar *xx;

120:   PetscFunctionBegin;
121:   PetscCall(VecGetArrayWrite(xin, &xx));
122:   PetscCall(PetscRandomGetValues(r, xin->map->n, xx));
123:   PetscCall(VecRestoreArrayWrite(xin, &xx));
124:   PetscFunctionReturn(PETSC_SUCCESS);
125: }

127: PetscErrorCode VecGetSize_Seq(Vec vin, PetscInt *size)
128: {
129:   PetscFunctionBegin;
130:   *size = vin->map->n;
131:   PetscFunctionReturn(PETSC_SUCCESS);
132: }

134: PetscErrorCode VecConjugate_Seq(Vec xin)
135: {
136:   PetscFunctionBegin;
137:   if (PetscDefined(USE_COMPLEX)) {
138:     const PetscInt n = xin->map->n;
139:     PetscScalar   *x;

141:     PetscCall(VecGetArray(xin, &x));
142:     for (PetscInt i = 0; i < n; ++i) x[i] = PetscConj(x[i]);
143:     PetscCall(VecRestoreArray(xin, &x));
144:   }
145:   PetscFunctionReturn(PETSC_SUCCESS);
146: }

148: PetscErrorCode VecResetArray_Seq(Vec vin)
149: {
150:   Vec_Seq *v = (Vec_Seq *)vin->data;

152:   PetscFunctionBegin;
153:   v->array         = v->unplacedarray;
154:   v->unplacedarray = NULL;
155:   PetscFunctionReturn(PETSC_SUCCESS);
156: }

158: PetscErrorCode VecCopy_Seq(Vec xin, Vec yin)
159: {
160:   PetscFunctionBegin;
161:   if (xin != yin) {
162:     const PetscScalar *xa;
163:     PetscScalar       *ya;

165:     PetscCall(VecGetArrayRead(xin, &xa));
166:     PetscCall(VecGetArray(yin, &ya));
167:     PetscCall(PetscArraycpy(ya, xa, xin->map->n));
168:     PetscCall(VecRestoreArrayRead(xin, &xa));
169:     PetscCall(VecRestoreArray(yin, &ya));
170:   }
171:   PetscFunctionReturn(PETSC_SUCCESS);
172: }

174: PetscErrorCode VecSwap_Seq(Vec xin, Vec yin)
175: {
176:   PetscFunctionBegin;
177:   if (xin != yin) {
178:     const PetscBLASInt one = 1;
179:     PetscScalar       *ya, *xa;
180:     PetscBLASInt       bn;

182:     PetscCall(PetscBLASIntCast(xin->map->n, &bn));
183:     PetscCall(VecGetArray(xin, &xa));
184:     PetscCall(VecGetArray(yin, &ya));
185:     PetscCallBLAS("BLASswap", BLASswap_(&bn, xa, &one, ya, &one));
186:     PetscCall(VecRestoreArray(xin, &xa));
187:     PetscCall(VecRestoreArray(yin, &ya));
188:   }
189:   PetscFunctionReturn(PETSC_SUCCESS);
190: }

192: #include <../src/vec/vec/impls/seq/ftn-kernels/fnorm.h>

194: PetscErrorCode VecNorm_Seq(Vec xin, NormType type, PetscReal *z)
195: {
196:   // use a local variable to ensure compiler doesn't think z aliases any of the other arrays
197:   PetscReal      ztmp[] = {0.0, 0.0};
198:   const PetscInt n      = xin->map->n;

200:   PetscFunctionBegin;
201:   if (n) {
202:     const PetscScalar *xx;
203:     const PetscBLASInt one = 1;
204:     PetscBLASInt       bn  = 0;

206:     PetscCall(PetscBLASIntCast(n, &bn));
207:     PetscCall(VecGetArrayRead(xin, &xx));
208:     if (type == NORM_2 || type == NORM_FROBENIUS) {
209:     NORM_1_AND_2_DOING_NORM_2:
210:       if (PetscDefined(USE_REAL___FP16)) {
211:         PetscCallBLAS("BLASnrm2", ztmp[type == NORM_1_AND_2] = BLASnrm2_(&bn, xx, &one));
212:       } else {
213:         PetscCallBLAS("BLASdot", ztmp[type == NORM_1_AND_2] = PetscSqrtReal(PetscRealPart(BLASdot_(&bn, xx, &one, xx, &one))));
214:       }
215:       PetscCall(PetscLogFlops(2.0 * n - 1));
216:     } else if (type == NORM_INFINITY) {
217:       for (PetscInt i = 0; i < n; ++i) {
218:         const PetscReal tmp = PetscAbsScalar(xx[i]);

220:         /* check special case of tmp == NaN */
221:         if ((tmp > ztmp[0]) || (tmp != tmp)) {
222:           ztmp[0] = tmp;
223:           if (tmp != tmp) break;
224:         }
225:       }
226:     } else if (type == NORM_1 || type == NORM_1_AND_2) {
227:       if (PetscDefined(USE_COMPLEX)) {
228:         // BLASasum() returns the nonstandard 1 norm of the 1 norm of the complex entries so we
229:         // provide a custom loop instead
230:         for (PetscInt i = 0; i < n; ++i) ztmp[0] += PetscAbsScalar(xx[i]);
231:       } else {
232:         PetscCallBLAS("BLASasum", ztmp[0] = BLASasum_(&bn, xx, &one));
233:       }
234:       PetscCall(PetscLogFlops(n - 1.0));
235:       /* slight reshuffle so we can skip getting the array again (but still log the flops) if we
236:          do norm2 after this */
237:       if (type == NORM_1_AND_2) goto NORM_1_AND_2_DOING_NORM_2;
238:     }
239:     PetscCall(VecRestoreArrayRead(xin, &xx));
240:   }
241:   z[0] = ztmp[0];
242:   if (type == NORM_1_AND_2) z[1] = ztmp[1];
243:   PetscFunctionReturn(PETSC_SUCCESS);
244: }

246: PetscErrorCode VecView_Seq_ASCII(Vec xin, PetscViewer viewer)
247: {
248:   PetscInt           i, n = xin->map->n;
249:   const char        *name;
250:   PetscViewerFormat  format;
251:   const PetscScalar *xv;

253:   PetscFunctionBegin;
254:   PetscCall(VecGetArrayRead(xin, &xv));
255:   PetscCall(PetscViewerGetFormat(viewer, &format));
256:   if (format == PETSC_VIEWER_ASCII_MATLAB) {
257:     PetscCall(PetscObjectGetName((PetscObject)xin, &name));
258:     PetscCall(PetscViewerASCIIPrintf(viewer, "%s = [\n", name));
259:     for (i = 0; i < n; i++) {
260: #if defined(PETSC_USE_COMPLEX)
261:       if (PetscImaginaryPart(xv[i]) > 0.0) {
262:         PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e + %18.16ei\n", (double)PetscRealPart(xv[i]), (double)PetscImaginaryPart(xv[i])));
263:       } else if (PetscImaginaryPart(xv[i]) < 0.0) {
264:         PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e - %18.16ei\n", (double)PetscRealPart(xv[i]), -(double)PetscImaginaryPart(xv[i])));
265:       } else {
266:         PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e\n", (double)PetscRealPart(xv[i])));
267:       }
268: #else
269:       PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e\n", (double)xv[i]));
270: #endif
271:     }
272:     PetscCall(PetscViewerASCIIPrintf(viewer, "];\n"));
273:   } else if (format == PETSC_VIEWER_ASCII_SYMMODU) {
274:     for (i = 0; i < n; i++) {
275: #if defined(PETSC_USE_COMPLEX)
276:       PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e %18.16e\n", (double)PetscRealPart(xv[i]), (double)PetscImaginaryPart(xv[i])));
277: #else
278:       PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e\n", (double)xv[i]));
279: #endif
280:     }
281:   } else if (format == PETSC_VIEWER_ASCII_VTK_DEPRECATED || format == PETSC_VIEWER_ASCII_VTK_CELL_DEPRECATED) {
282:     /*
283:        state 0: No header has been output
284:        state 1: Only POINT_DATA has been output
285:        state 2: Only CELL_DATA has been output
286:        state 3: Output both, POINT_DATA last
287:        state 4: Output both, CELL_DATA last
288:     */
289:     static PetscInt stateId     = -1;
290:     int             outputState = 0;
291:     PetscBool       hasState;
292:     int             doOutput = 0;
293:     PetscInt        bs, b;

295:     if (stateId < 0) PetscCall(PetscObjectComposedDataRegister(&stateId));
296:     PetscCall(PetscObjectComposedDataGetInt((PetscObject)viewer, stateId, outputState, hasState));
297:     if (!hasState) outputState = 0;
298:     PetscCall(PetscObjectGetName((PetscObject)xin, &name));
299:     PetscCall(VecGetBlockSize(xin, &bs));
300:     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);
301:     if (format == PETSC_VIEWER_ASCII_VTK_DEPRECATED) {
302:       if (outputState == 0) {
303:         outputState = 1;
304:         doOutput    = 1;
305:       } else if (outputState == 1) doOutput = 0;
306:       else if (outputState == 2) {
307:         outputState = 3;
308:         doOutput    = 1;
309:       } else if (outputState == 3) doOutput = 0;
310:       else PetscCheck(outputState != 4, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Tried to output POINT_DATA again after intervening CELL_DATA");

312:       if (doOutput) PetscCall(PetscViewerASCIIPrintf(viewer, "POINT_DATA %" PetscInt_FMT "\n", n / bs));
313:     } else {
314:       if (outputState == 0) {
315:         outputState = 2;
316:         doOutput    = 1;
317:       } else if (outputState == 1) {
318:         outputState = 4;
319:         doOutput    = 1;
320:       } else if (outputState == 2) {
321:         doOutput = 0;
322:       } else {
323:         PetscCheck(outputState != 3, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Tried to output CELL_DATA again after intervening POINT_DATA");
324:         if (outputState == 4) doOutput = 0;
325:       }

327:       if (doOutput) PetscCall(PetscViewerASCIIPrintf(viewer, "CELL_DATA %" PetscInt_FMT "\n", n));
328:     }
329:     PetscCall(PetscObjectComposedDataSetInt((PetscObject)viewer, stateId, outputState));
330:     if (name) {
331:       if (bs == 3) {
332:         PetscCall(PetscViewerASCIIPrintf(viewer, "VECTORS %s double\n", name));
333:       } else {
334:         PetscCall(PetscViewerASCIIPrintf(viewer, "SCALARS %s double %" PetscInt_FMT "\n", name, bs));
335:       }
336:     } else {
337:       PetscCall(PetscViewerASCIIPrintf(viewer, "SCALARS scalars double %" PetscInt_FMT "\n", bs));
338:     }
339:     if (bs != 3) PetscCall(PetscViewerASCIIPrintf(viewer, "LOOKUP_TABLE default\n"));
340:     for (i = 0; i < n / bs; i++) {
341:       for (b = 0; b < bs; b++) {
342:         if (b > 0) PetscCall(PetscViewerASCIIPrintf(viewer, " "));
343: #if !defined(PETSC_USE_COMPLEX)
344:         PetscCall(PetscViewerASCIIPrintf(viewer, "%g", (double)xv[i * bs + b]));
345: #endif
346:       }
347:       PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
348:     }
349:   } else if (format == PETSC_VIEWER_ASCII_VTK_COORDS_DEPRECATED) {
350:     PetscInt bs, b;

352:     PetscCall(VecGetBlockSize(xin, &bs));
353:     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);
354:     for (i = 0; i < n / bs; i++) {
355:       for (b = 0; b < bs; b++) {
356:         if (b > 0) PetscCall(PetscViewerASCIIPrintf(viewer, " "));
357: #if !defined(PETSC_USE_COMPLEX)
358:         PetscCall(PetscViewerASCIIPrintf(viewer, "%g", (double)xv[i * bs + b]));
359: #endif
360:       }
361:       for (b = bs; b < 3; b++) PetscCall(PetscViewerASCIIPrintf(viewer, " 0.0"));
362:       PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
363:     }
364:   } else if (format == PETSC_VIEWER_ASCII_PCICE) {
365:     PetscInt bs, b;

367:     PetscCall(VecGetBlockSize(xin, &bs));
368:     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);
369:     PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT "\n", xin->map->N / bs));
370:     for (i = 0; i < n / bs; i++) {
371:       PetscCall(PetscViewerASCIIPrintf(viewer, "%7" PetscInt_FMT "   ", i + 1));
372:       for (b = 0; b < bs; b++) {
373:         if (b > 0) PetscCall(PetscViewerASCIIPrintf(viewer, " "));
374: #if !defined(PETSC_USE_COMPLEX)
375:         PetscCall(PetscViewerASCIIPrintf(viewer, "% 12.5E", (double)xv[i * bs + b]));
376: #endif
377:       }
378:       PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
379:     }
380:   } else if (format == PETSC_VIEWER_ASCII_GLVIS) {
381:     /* GLVis ASCII visualization/dump: this function mimics mfem::GridFunction::Save() */
382:     const PetscScalar      *array;
383:     PetscInt                i, n, vdim, ordering = 1; /* mfem::FiniteElementSpace::Ordering::byVDIM */
384:     PetscContainer          glvis_container;
385:     PetscViewerGLVisVecInfo glvis_vec_info;
386:     PetscViewerGLVisInfo    glvis_info;

388:     /* mfem::FiniteElementSpace::Save() */
389:     PetscCall(VecGetBlockSize(xin, &vdim));
390:     PetscCall(PetscViewerASCIIPrintf(viewer, "FiniteElementSpace\n"));
391:     PetscCall(PetscObjectQuery((PetscObject)xin, "_glvis_info_container", (PetscObject *)&glvis_container));
392:     PetscCheck(glvis_container, PetscObjectComm((PetscObject)xin), PETSC_ERR_PLIB, "Missing GLVis container");
393:     PetscCall(PetscContainerGetPointer(glvis_container, (void **)&glvis_vec_info));
394:     PetscCall(PetscViewerASCIIPrintf(viewer, "%s\n", glvis_vec_info->fec_type));
395:     PetscCall(PetscViewerASCIIPrintf(viewer, "VDim: %" PetscInt_FMT "\n", vdim));
396:     PetscCall(PetscViewerASCIIPrintf(viewer, "Ordering: %" PetscInt_FMT "\n", ordering));
397:     PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
398:     /* mfem::Vector::Print() */
399:     PetscCall(PetscObjectQuery((PetscObject)viewer, "_glvis_info_container", (PetscObject *)&glvis_container));
400:     PetscCheck(glvis_container, PetscObjectComm((PetscObject)viewer), PETSC_ERR_PLIB, "Missing GLVis container");
401:     PetscCall(PetscContainerGetPointer(glvis_container, (void **)&glvis_info));
402:     if (glvis_info->enabled) {
403:       PetscCall(VecGetLocalSize(xin, &n));
404:       PetscCall(VecGetArrayRead(xin, &array));
405:       for (i = 0; i < n; i++) {
406:         PetscCall(PetscViewerASCIIPrintf(viewer, glvis_info->fmt, (double)PetscRealPart(array[i])));
407:         PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
408:       }
409:       PetscCall(VecRestoreArrayRead(xin, &array));
410:     }
411:   } else if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
412:     /* No info */
413:   } else {
414:     for (i = 0; i < n; i++) {
415:       if (format == PETSC_VIEWER_ASCII_INDEX) PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT ": ", i));
416: #if defined(PETSC_USE_COMPLEX)
417:       if (PetscImaginaryPart(xv[i]) > 0.0) {
418:         PetscCall(PetscViewerASCIIPrintf(viewer, "%g + %g i\n", (double)PetscRealPart(xv[i]), (double)PetscImaginaryPart(xv[i])));
419:       } else if (PetscImaginaryPart(xv[i]) < 0.0) {
420:         PetscCall(PetscViewerASCIIPrintf(viewer, "%g - %g i\n", (double)PetscRealPart(xv[i]), -(double)PetscImaginaryPart(xv[i])));
421:       } else {
422:         PetscCall(PetscViewerASCIIPrintf(viewer, "%g\n", (double)PetscRealPart(xv[i])));
423:       }
424: #else
425:       PetscCall(PetscViewerASCIIPrintf(viewer, "%g\n", (double)xv[i]));
426: #endif
427:     }
428:   }
429:   PetscCall(PetscViewerFlush(viewer));
430:   PetscCall(VecRestoreArrayRead(xin, &xv));
431:   PetscFunctionReturn(PETSC_SUCCESS);
432: }

434: #include <petscdraw.h>
435: PetscErrorCode VecView_Seq_Draw_LG(Vec xin, PetscViewer v)
436: {
437:   PetscDraw          draw;
438:   PetscBool          isnull;
439:   PetscDrawLG        lg;
440:   PetscInt           i, c, bs = PetscAbs(xin->map->bs), n = xin->map->n / bs;
441:   const PetscScalar *xv;
442:   PetscReal         *xx, *yy, xmin, xmax, h;
443:   int                colors[] = {PETSC_DRAW_RED};
444:   PetscViewerFormat  format;
445:   PetscDrawAxis      axis;

447:   PetscFunctionBegin;
448:   PetscCall(PetscViewerDrawGetDraw(v, 0, &draw));
449:   PetscCall(PetscDrawIsNull(draw, &isnull));
450:   if (isnull) PetscFunctionReturn(PETSC_SUCCESS);

452:   PetscCall(PetscViewerGetFormat(v, &format));
453:   PetscCall(PetscMalloc2(n, &xx, n, &yy));
454:   PetscCall(VecGetArrayRead(xin, &xv));
455:   for (c = 0; c < bs; c++) {
456:     PetscCall(PetscViewerDrawGetDrawLG(v, c, &lg));
457:     PetscCall(PetscDrawLGReset(lg));
458:     PetscCall(PetscDrawLGSetDimension(lg, 1));
459:     PetscCall(PetscDrawLGSetColors(lg, colors));
460:     if (format == PETSC_VIEWER_DRAW_LG_XRANGE) {
461:       PetscCall(PetscDrawLGGetAxis(lg, &axis));
462:       PetscCall(PetscDrawAxisGetLimits(axis, &xmin, &xmax, NULL, NULL));
463:       h = (xmax - xmin) / n;
464:       for (i = 0; i < n; i++) xx[i] = i * h + 0.5 * h; /* cell center */
465:     } else {
466:       for (i = 0; i < n; i++) xx[i] = (PetscReal)i;
467:     }
468:     for (i = 0; i < n; i++) yy[i] = PetscRealPart(xv[c + i * bs]);

470:     PetscCall(PetscDrawLGAddPoints(lg, n, &xx, &yy));
471:     PetscCall(PetscDrawLGDraw(lg));
472:     PetscCall(PetscDrawLGSave(lg));
473:   }
474:   PetscCall(VecRestoreArrayRead(xin, &xv));
475:   PetscCall(PetscFree2(xx, yy));
476:   PetscFunctionReturn(PETSC_SUCCESS);
477: }

479: PetscErrorCode VecView_Seq_Draw(Vec xin, PetscViewer v)
480: {
481:   PetscDraw draw;
482:   PetscBool isnull;

484:   PetscFunctionBegin;
485:   PetscCall(PetscViewerDrawGetDraw(v, 0, &draw));
486:   PetscCall(PetscDrawIsNull(draw, &isnull));
487:   if (isnull) PetscFunctionReturn(PETSC_SUCCESS);

489:   PetscCall(VecView_Seq_Draw_LG(xin, v));
490:   PetscFunctionReturn(PETSC_SUCCESS);
491: }

493: PetscErrorCode VecView_Seq_Binary(Vec xin, PetscViewer viewer)
494: {
495:   return VecView_Binary(xin, viewer);
496: }

498: #if defined(PETSC_HAVE_MATLAB)
499: #include <petscmatlab.h>
500:   #include <mat.h> /* MATLAB include file */
501: PetscErrorCode VecView_Seq_Matlab(Vec vec, PetscViewer viewer)
502: {
503:   PetscInt           n;
504:   const PetscScalar *array;

506:   PetscFunctionBegin;
507:   PetscCall(VecGetLocalSize(vec, &n));
508:   PetscCall(PetscObjectName((PetscObject)vec));
509:   PetscCall(VecGetArrayRead(vec, &array));
510:   PetscCall(PetscViewerMatlabPutArray(viewer, n, 1, array, ((PetscObject)vec)->name));
511:   PetscCall(VecRestoreArrayRead(vec, &array));
512:   PetscFunctionReturn(PETSC_SUCCESS);
513: }
514: #endif

516: PETSC_EXTERN PetscErrorCode VecView_Seq(Vec xin, PetscViewer viewer)
517: {
518:   PetscBool isdraw, iascii, issocket, isbinary;
519: #if defined(PETSC_HAVE_MATHEMATICA)
520:   PetscBool ismathematica;
521: #endif
522: #if defined(PETSC_HAVE_MATLAB)
523:   PetscBool ismatlab;
524: #endif
525: #if defined(PETSC_HAVE_HDF5)
526:   PetscBool ishdf5;
527: #endif
528:   PetscBool isglvis;
529: #if defined(PETSC_HAVE_ADIOS)
530:   PetscBool isadios;
531: #endif

533:   PetscFunctionBegin;
534:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
535:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
536:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSOCKET, &issocket));
537:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
538: #if defined(PETSC_HAVE_MATHEMATICA)
539:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERMATHEMATICA, &ismathematica));
540: #endif
541: #if defined(PETSC_HAVE_HDF5)
542:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
543: #endif
544: #if defined(PETSC_HAVE_MATLAB)
545:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERMATLAB, &ismatlab));
546: #endif
547:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERGLVIS, &isglvis));
548: #if defined(PETSC_HAVE_ADIOS)
549:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERADIOS, &isadios));
550: #endif

552:   if (isdraw) {
553:     PetscCall(VecView_Seq_Draw(xin, viewer));
554:   } else if (iascii) {
555:     PetscCall(VecView_Seq_ASCII(xin, viewer));
556:   } else if (isbinary) {
557:     PetscCall(VecView_Seq_Binary(xin, viewer));
558: #if defined(PETSC_HAVE_MATHEMATICA)
559:   } else if (ismathematica) {
560:     PetscCall(PetscViewerMathematicaPutVector(viewer, xin));
561: #endif
562: #if defined(PETSC_HAVE_HDF5)
563:   } else if (ishdf5) {
564:     PetscCall(VecView_MPI_HDF5(xin, viewer)); /* Reusing VecView_MPI_HDF5 ... don't want code duplication*/
565: #endif
566: #if defined(PETSC_HAVE_ADIOS)
567:   } else if (isadios) {
568:     PetscCall(VecView_MPI_ADIOS(xin, viewer)); /* Reusing VecView_MPI_ADIOS ... don't want code duplication*/
569: #endif
570: #if defined(PETSC_HAVE_MATLAB)
571:   } else if (ismatlab) {
572:     PetscCall(VecView_Seq_Matlab(xin, viewer));
573: #endif
574:   } else if (isglvis) PetscCall(VecView_GLVis(xin, viewer));
575:   PetscFunctionReturn(PETSC_SUCCESS);
576: }

578: PetscErrorCode VecGetValues_Seq(Vec xin, PetscInt ni, const PetscInt ix[], PetscScalar y[])
579: {
580:   const PetscBool    ignorenegidx = xin->stash.ignorenegidx;
581:   const PetscScalar *xx;

583:   PetscFunctionBegin;
584:   PetscCall(VecGetArrayRead(xin, &xx));
585:   for (PetscInt i = 0; i < ni; ++i) {
586:     if (ignorenegidx && (ix[i] < 0)) continue;
587:     if (PetscDefined(USE_DEBUG)) {
588:       PetscCheck(ix[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Out of range index value %" PetscInt_FMT " cannot be negative", ix[i]);
589:       PetscCheck(ix[i] < xin->map->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Out of range index value %" PetscInt_FMT " to large maximum allowed %" PetscInt_FMT, ix[i], xin->map->n);
590:     }
591:     y[i] = xx[ix[i]];
592:   }
593:   PetscCall(VecRestoreArrayRead(xin, &xx));
594:   PetscFunctionReturn(PETSC_SUCCESS);
595: }

597: PetscErrorCode VecSetValues_Seq(Vec xin, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode m)
598: {
599:   const PetscBool ignorenegidx = xin->stash.ignorenegidx;
600:   PetscScalar    *xx;

602:   PetscFunctionBegin;
603:   // call to getarray (not e.g. getarraywrite() if m is INSERT_VALUES) is deliberate! If this
604:   // is secretly a VECSEQCUDA it may have values currently on the device, in which case --
605:   // unless we are replacing the entire array -- we need to copy them up
606:   PetscCall(VecGetArray(xin, &xx));
607:   for (PetscInt i = 0; i < ni; i++) {
608:     if (ignorenegidx && (ix[i] < 0)) continue;
609:     if (PetscDefined(USE_DEBUG)) {
610:       PetscCheck(ix[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Out of range index value %" PetscInt_FMT " cannot be negative", ix[i]);
611:       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);
612:     }
613:     if (m == INSERT_VALUES) {
614:       xx[ix[i]] = y[i];
615:     } else {
616:       xx[ix[i]] += y[i];
617:     }
618:   }
619:   PetscCall(VecRestoreArray(xin, &xx));
620:   PetscFunctionReturn(PETSC_SUCCESS);
621: }

623: PetscErrorCode VecSetValuesBlocked_Seq(Vec xin, PetscInt ni, const PetscInt ix[], const PetscScalar yin[], InsertMode m)
624: {
625:   PetscScalar *xx;
626:   PetscInt     bs;

628:   /* For optimization could treat bs = 2, 3, 4, 5 as special cases with loop unrolling */
629:   PetscFunctionBegin;
630:   PetscCall(VecGetBlockSize(xin, &bs));
631:   PetscCall(VecGetArray(xin, &xx));
632:   for (PetscInt i = 0; i < ni; ++i, yin += bs) {
633:     const PetscInt start = bs * ix[i];

635:     if (start < 0) continue;
636:     PetscCheck(start < xin->map->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Out of range index value %" PetscInt_FMT " maximum %" PetscInt_FMT, start, xin->map->n);
637:     for (PetscInt j = 0; j < bs; j++) {
638:       if (m == INSERT_VALUES) {
639:         xx[start + j] = yin[j];
640:       } else {
641:         xx[start + j] += yin[j];
642:       }
643:     }
644:   }
645:   PetscCall(VecRestoreArray(xin, &xx));
646:   PetscFunctionReturn(PETSC_SUCCESS);
647: }

649: static PetscErrorCode VecResetPreallocationCOO_Seq(Vec x)
650: {
651:   Vec_Seq *vs = (Vec_Seq *)x->data;

653:   PetscFunctionBegin;
654:   if (vs) {
655:     PetscCall(PetscFree(vs->jmap1)); /* Destroy old stuff */
656:     PetscCall(PetscFree(vs->perm1));
657:   }
658:   PetscFunctionReturn(PETSC_SUCCESS);
659: }

661: PetscErrorCode VecSetPreallocationCOO_Seq(Vec x, PetscCount coo_n, const PetscInt coo_i[])
662: {
663:   PetscInt    m, *i;
664:   PetscCount  k, nneg;
665:   PetscCount *perm1, *jmap1;
666:   Vec_Seq    *vs = (Vec_Seq *)x->data;

668:   PetscFunctionBegin;
669:   PetscCall(VecResetPreallocationCOO_Seq(x)); /* Destroy old stuff */
670:   PetscCall(PetscMalloc1(coo_n, &i));
671:   PetscCall(PetscArraycpy(i, coo_i, coo_n)); /* Make a copy since we'll modify it */
672:   PetscCall(PetscMalloc1(coo_n, &perm1));
673:   for (k = 0; k < coo_n; k++) perm1[k] = k;
674:   PetscCall(PetscSortIntWithCountArray(coo_n, i, perm1));
675:   for (k = 0; k < coo_n; k++) {
676:     if (i[k] >= 0) break;
677:   } /* Advance k to the first entry with a non-negative index */
678:   nneg = k;

680:   PetscCall(VecGetLocalSize(x, &m));
681:   PetscCheck(!nneg || x->stash.ignorenegidx, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Found a negative index in VecSetPreallocateCOO() but VEC_IGNORE_NEGATIVE_INDICES was not set");
682:   PetscCheck(!coo_n || i[coo_n - 1] < m, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Found index (%" PetscInt_FMT ") greater than the size of the vector (%" PetscInt_FMT ") in VecSetPreallocateCOO()", i[coo_n - 1], m);

684:   PetscCall(PetscCalloc1(m + 1, &jmap1));
685:   for (; k < coo_n; k++) jmap1[i[k] + 1]++;         /* Count repeats of each entry */
686:   for (k = 0; k < m; k++) jmap1[k + 1] += jmap1[k]; /* Transform jmap[] to CSR-like data structure */
687:   PetscCall(PetscFree(i));

689:   if (nneg) { /* Discard leading negative indices */
690:     PetscCount *perm1_new;
691:     PetscCall(PetscMalloc1(coo_n - nneg, &perm1_new));
692:     PetscCall(PetscArraycpy(perm1_new, perm1 + nneg, coo_n - nneg));
693:     PetscCall(PetscFree(perm1));
694:     perm1 = perm1_new;
695:   }

697:   /* Record COO fields */
698:   vs->coo_n = coo_n;
699:   vs->tot1  = coo_n - nneg;
700:   vs->jmap1 = jmap1; /* [m+1] */
701:   vs->perm1 = perm1; /* [tot] */
702:   PetscFunctionReturn(PETSC_SUCCESS);
703: }

705: PetscErrorCode VecSetValuesCOO_Seq(Vec x, const PetscScalar coo_v[], InsertMode imode)
706: {
707:   Vec_Seq          *vs    = (Vec_Seq *)x->data;
708:   const PetscCount *perm1 = vs->perm1, *jmap1 = vs->jmap1;
709:   PetscScalar      *xv;
710:   PetscInt          m;

712:   PetscFunctionBegin;
713:   PetscCall(VecGetLocalSize(x, &m));
714:   PetscCall(VecGetArray(x, &xv));
715:   for (PetscInt i = 0; i < m; i++) {
716:     PetscScalar sum = 0.0;
717:     for (PetscCount j = jmap1[i]; j < jmap1[i + 1]; j++) sum += coo_v[perm1[j]];
718:     xv[i] = (imode == INSERT_VALUES ? 0.0 : xv[i]) + sum;
719:   }
720:   PetscCall(VecRestoreArray(x, &xv));
721:   PetscFunctionReturn(PETSC_SUCCESS);
722: }

724: PetscErrorCode VecDestroy_Seq(Vec v)
725: {
726:   Vec_Seq *vs = (Vec_Seq *)v->data;

728:   PetscFunctionBegin;
729: #if defined(PETSC_USE_LOG)
730:   PetscCall(PetscLogObjectState((PetscObject)v, "Length=%" PetscInt_FMT, v->map->n));
731: #endif
732:   if (vs) PetscCall(PetscFree(vs->array_allocated));
733:   PetscCall(VecResetPreallocationCOO_Seq(v));
734:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscMatlabEnginePut_C", NULL));
735:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscMatlabEngineGet_C", NULL));
736:   PetscCall(PetscFree(v->data));
737:   PetscFunctionReturn(PETSC_SUCCESS);
738: }

740: PetscErrorCode VecSetOption_Seq(Vec v, VecOption op, PetscBool flag)
741: {
742:   PetscFunctionBegin;
743:   if (op == VEC_IGNORE_NEGATIVE_INDICES) v->stash.ignorenegidx = flag;
744:   PetscFunctionReturn(PETSC_SUCCESS);
745: }

747: PetscErrorCode VecDuplicate_Seq(Vec win, Vec *V)
748: {
749:   PetscFunctionBegin;
750:   PetscCall(VecCreate(PetscObjectComm((PetscObject)win), V));
751:   PetscCall(VecSetSizes(*V, win->map->n, win->map->n));
752:   PetscCall(VecSetType(*V, ((PetscObject)win)->type_name));
753:   PetscCall(PetscLayoutReference(win->map, &(*V)->map));
754:   PetscCall(PetscObjectListDuplicate(((PetscObject)win)->olist, &((PetscObject)(*V))->olist));
755:   PetscCall(PetscFunctionListDuplicate(((PetscObject)win)->qlist, &((PetscObject)(*V))->qlist));

757:   (*V)->ops->view          = win->ops->view;
758:   (*V)->stash.ignorenegidx = win->stash.ignorenegidx;
759:   PetscFunctionReturn(PETSC_SUCCESS);
760: }

762: static struct _VecOps DvOps = {
763:   PetscDesignatedInitializer(duplicate, VecDuplicate_Seq), /* 1 */
764:   PetscDesignatedInitializer(duplicatevecs, VecDuplicateVecs_Default),
765:   PetscDesignatedInitializer(destroyvecs, VecDestroyVecs_Default),
766:   PetscDesignatedInitializer(dot, VecDot_Seq),
767:   PetscDesignatedInitializer(mdot, VecMDot_Seq),
768:   PetscDesignatedInitializer(norm, VecNorm_Seq),
769:   PetscDesignatedInitializer(tdot, VecTDot_Seq),
770:   PetscDesignatedInitializer(mtdot, VecMTDot_Seq),
771:   PetscDesignatedInitializer(scale, VecScale_Seq),
772:   PetscDesignatedInitializer(copy, VecCopy_Seq), /* 10 */
773:   PetscDesignatedInitializer(set, VecSet_Seq),
774:   PetscDesignatedInitializer(swap, VecSwap_Seq),
775:   PetscDesignatedInitializer(axpy, VecAXPY_Seq),
776:   PetscDesignatedInitializer(axpby, VecAXPBY_Seq),
777:   PetscDesignatedInitializer(maxpy, VecMAXPY_Seq),
778:   PetscDesignatedInitializer(aypx, VecAYPX_Seq),
779:   PetscDesignatedInitializer(waxpy, VecWAXPY_Seq),
780:   PetscDesignatedInitializer(axpbypcz, VecAXPBYPCZ_Seq),
781:   PetscDesignatedInitializer(pointwisemult, VecPointwiseMult_Seq),
782:   PetscDesignatedInitializer(pointwisedivide, VecPointwiseDivide_Seq),
783:   PetscDesignatedInitializer(setvalues, VecSetValues_Seq), /* 20 */
784:   PetscDesignatedInitializer(assemblybegin, NULL),
785:   PetscDesignatedInitializer(assemblyend, NULL),
786:   PetscDesignatedInitializer(getarray, NULL),
787:   PetscDesignatedInitializer(getsize, VecGetSize_Seq),
788:   PetscDesignatedInitializer(getlocalsize, VecGetSize_Seq),
789:   PetscDesignatedInitializer(restorearray, NULL),
790:   PetscDesignatedInitializer(max, VecMax_Seq),
791:   PetscDesignatedInitializer(min, VecMin_Seq),
792:   PetscDesignatedInitializer(setrandom, VecSetRandom_Seq),
793:   PetscDesignatedInitializer(setoption, VecSetOption_Seq), /* 30 */
794:   PetscDesignatedInitializer(setvaluesblocked, VecSetValuesBlocked_Seq),
795:   PetscDesignatedInitializer(destroy, VecDestroy_Seq),
796:   PetscDesignatedInitializer(view, VecView_Seq),
797:   PetscDesignatedInitializer(placearray, VecPlaceArray_Seq),
798:   PetscDesignatedInitializer(replacearray, VecReplaceArray_Seq),
799:   PetscDesignatedInitializer(dot_local, VecDot_Seq),
800:   PetscDesignatedInitializer(tdot_local, VecTDot_Seq),
801:   PetscDesignatedInitializer(norm_local, VecNorm_Seq),
802:   PetscDesignatedInitializer(mdot_local, VecMDot_Seq),
803:   PetscDesignatedInitializer(mtdot_local, VecMTDot_Seq), /* 40 */
804:   PetscDesignatedInitializer(load, VecLoad_Default),
805:   PetscDesignatedInitializer(reciprocal, VecReciprocal_Default),
806:   PetscDesignatedInitializer(conjugate, VecConjugate_Seq),
807:   PetscDesignatedInitializer(setlocaltoglobalmapping, NULL),
808:   PetscDesignatedInitializer(setvalueslocal, NULL),
809:   PetscDesignatedInitializer(resetarray, VecResetArray_Seq),
810:   PetscDesignatedInitializer(setfromoptions, NULL),
811:   PetscDesignatedInitializer(maxpointwisedivide, VecMaxPointwiseDivide_Seq),
812:   PetscDesignatedInitializer(pointwisemax, VecPointwiseMax_Seq),
813:   PetscDesignatedInitializer(pointwisemaxabs, VecPointwiseMaxAbs_Seq),
814:   PetscDesignatedInitializer(pointwisemin, VecPointwiseMin_Seq),
815:   PetscDesignatedInitializer(getvalues, VecGetValues_Seq),
816:   PetscDesignatedInitializer(sqrt, NULL),
817:   PetscDesignatedInitializer(abs, NULL),
818:   PetscDesignatedInitializer(exp, NULL),
819:   PetscDesignatedInitializer(log, NULL),
820:   PetscDesignatedInitializer(shift, NULL),
821:   PetscDesignatedInitializer(create, NULL),
822:   PetscDesignatedInitializer(stridegather, VecStrideGather_Default),
823:   PetscDesignatedInitializer(stridescatter, VecStrideScatter_Default),
824:   PetscDesignatedInitializer(dotnorm2, NULL),
825:   PetscDesignatedInitializer(getsubvector, NULL),
826:   PetscDesignatedInitializer(restoresubvector, NULL),
827:   PetscDesignatedInitializer(getarrayread, NULL),
828:   PetscDesignatedInitializer(restorearrayread, NULL),
829:   PetscDesignatedInitializer(stridesubsetgather, VecStrideSubSetGather_Default),
830:   PetscDesignatedInitializer(stridesubsetscatter, VecStrideSubSetScatter_Default),
831:   PetscDesignatedInitializer(viewnative, VecView_Seq),
832:   PetscDesignatedInitializer(loadnative, NULL),
833:   PetscDesignatedInitializer(createlocalvector, NULL),
834:   PetscDesignatedInitializer(getlocalvector, NULL),
835:   PetscDesignatedInitializer(restorelocalvector, NULL),
836:   PetscDesignatedInitializer(getlocalvectorread, NULL),
837:   PetscDesignatedInitializer(restorelocalvectorread, NULL),
838:   PetscDesignatedInitializer(bindtocpu, NULL),
839:   PetscDesignatedInitializer(getarraywrite, NULL),
840:   PetscDesignatedInitializer(restorearraywrite, NULL),
841:   PetscDesignatedInitializer(getarrayandmemtype, NULL),
842:   PetscDesignatedInitializer(restorearrayandmemtype, NULL),
843:   PetscDesignatedInitializer(getarrayreadandmemtype, NULL),
844:   PetscDesignatedInitializer(restorearrayreadandmemtype, NULL),
845:   PetscDesignatedInitializer(getarraywriteandmemtype, NULL),
846:   PetscDesignatedInitializer(restorearraywriteandmemtype, NULL),
847:   PetscDesignatedInitializer(concatenate, NULL),
848:   PetscDesignatedInitializer(sum, NULL),
849:   PetscDesignatedInitializer(setpreallocationcoo, VecSetPreallocationCOO_Seq),
850:   PetscDesignatedInitializer(setvaluescoo, VecSetValuesCOO_Seq),
851: };

853: /*
854:       This is called by VecCreate_Seq() (i.e. VecCreateSeq()) and VecCreateSeqWithArray()
855: */
856: PetscErrorCode VecCreate_Seq_Private(Vec v, const PetscScalar array[])
857: {
858:   Vec_Seq *s;

860:   PetscFunctionBegin;
861:   PetscCall(PetscNew(&s));
862:   PetscCall(PetscMemcpy(v->ops, &DvOps, sizeof(DvOps)));

864:   v->data            = (void *)s;
865:   v->petscnative     = PETSC_TRUE;
866:   s->array           = (PetscScalar *)array;
867:   s->array_allocated = NULL;
868:   if (array) v->offloadmask = PETSC_OFFLOAD_CPU;

870:   PetscCall(PetscLayoutSetUp(v->map));
871:   PetscCall(PetscObjectChangeTypeName((PetscObject)v, VECSEQ));
872: #if defined(PETSC_HAVE_MATLAB)
873:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscMatlabEnginePut_C", VecMatlabEnginePut_Default));
874:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscMatlabEngineGet_C", VecMatlabEngineGet_Default));
875: #endif
876:   PetscFunctionReturn(PETSC_SUCCESS);
877: }

879: /*@C
880:    VecCreateSeqWithArray - Creates a standard,sequential array-style vector,
881:    where the user provides the array space to store the vector values.

883:    Collective

885:    Input Parameters:
886: +  comm - the communicator, should be `PETSC_COMM_SELF`
887: .  bs - the block size
888: .  n - the vector length
889: -  array - memory where the vector elements are to be stored.

891:    Output Parameter:
892: .  V - the vector

894:    Level: intermediate

896:    Notes:
897:    Use `VecDuplicate()` or `VecDuplicateVecs(`) to form additional vectors of the
898:    same type as an existing vector.

900:    If the user-provided array is` NULL`, then `VecPlaceArray()` can be used
901:    at a later stage to SET the array for storing the vector values.

903:    PETSc does NOT free the array when the vector is destroyed via `VecDestroy()`.
904:    The user should not free the array until the vector is destroyed.

906: .seealso: `VecCreateMPIWithArray()`, `VecCreate()`, `VecDuplicate()`, `VecDuplicateVecs()`,
907:           `VecCreateGhost()`, `VecCreateSeq()`, `VecPlaceArray()`
908: @*/
909: PetscErrorCode VecCreateSeqWithArray(MPI_Comm comm, PetscInt bs, PetscInt n, const PetscScalar array[], Vec *V)
910: {
911:   PetscMPIInt size;

913:   PetscFunctionBegin;
914:   PetscCall(VecCreate(comm, V));
915:   PetscCall(VecSetSizes(*V, n, n));
916:   PetscCall(VecSetBlockSize(*V, bs));
917:   PetscCallMPI(MPI_Comm_size(comm, &size));
918:   PetscCheck(size <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot create VECSEQ on more than one process");
919:   PetscCall(VecCreate_Seq_Private(*V, array));
920:   PetscFunctionReturn(PETSC_SUCCESS);
921: }