Actual source code: bvec2.c

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

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

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

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

 28: static PetscScalar MaxRealPart(PetscScalar x, PetscScalar y)
 29: {
 30:   // use temporaries to avoid reevaluating side-effects
 31:   const PetscReal rx = PetscRealPart(x), ry = PetscRealPart(y);

 33:   return PetscMax(rx, ry);
 34: }

 36: PetscErrorCode VecPointwiseMax_Seq(Vec win, Vec xin, Vec yin)
 37: {
 38:   PetscFunctionBegin;
 39:   PetscCall(VecPointwiseApply_Seq(win, xin, yin, MaxRealPart));
 40:   PetscFunctionReturn(PETSC_SUCCESS);
 41: }

 43: static PetscScalar MinRealPart(PetscScalar x, PetscScalar y)
 44: {
 45:   // use temporaries to avoid reevaluating side-effects
 46:   const PetscReal rx = PetscRealPart(x), ry = PetscRealPart(y);

 48:   return PetscMin(rx, ry);
 49: }

 51: PetscErrorCode VecPointwiseMin_Seq(Vec win, Vec xin, Vec yin)
 52: {
 53:   PetscFunctionBegin;
 54:   PetscCall(VecPointwiseApply_Seq(win, xin, yin, MinRealPart));
 55:   PetscFunctionReturn(PETSC_SUCCESS);
 56: }

 58: static PetscScalar MaxAbs(PetscScalar x, PetscScalar y)
 59: {
 60:   return (PetscScalar)PetscMax(PetscAbsScalar(x), PetscAbsScalar(y));
 61: }

 63: PetscErrorCode VecPointwiseMaxAbs_Seq(Vec win, Vec xin, Vec yin)
 64: {
 65:   PetscFunctionBegin;
 66:   PetscCall(VecPointwiseApply_Seq(win, xin, yin, MaxAbs));
 67:   PetscFunctionReturn(PETSC_SUCCESS);
 68: }

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

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

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

 99: static PetscScalar ScalDiv(PetscScalar x, PetscScalar y)
100: {
101:   return y == 0.0 ? 0.0 : x / y;
102: }

104: PetscErrorCode VecPointwiseDivide_Seq(Vec win, Vec xin, Vec yin)
105: {
106:   PetscFunctionBegin;
107:   PetscCall(VecPointwiseApply_Seq(win, xin, yin, ScalDiv));
108:   PetscFunctionReturn(PETSC_SUCCESS);
109: }

111: PetscErrorCode VecSetRandom_Seq(Vec xin, PetscRandom r)
112: {
113:   PetscScalar *xx;

115:   PetscFunctionBegin;
116:   PetscCall(VecGetArrayWrite(xin, &xx));
117:   PetscCall(PetscRandomGetValues(r, xin->map->n, xx));
118:   PetscCall(VecRestoreArrayWrite(xin, &xx));
119:   PetscFunctionReturn(PETSC_SUCCESS);
120: }

122: PetscErrorCode VecGetSize_Seq(Vec vin, PetscInt *size)
123: {
124:   PetscFunctionBegin;
125:   *size = vin->map->n;
126:   PetscFunctionReturn(PETSC_SUCCESS);
127: }

129: PetscErrorCode VecConjugate_Seq(Vec xin)
130: {
131:   PetscFunctionBegin;
132:   if (PetscDefined(USE_COMPLEX)) {
133:     const PetscInt n = xin->map->n;
134:     PetscScalar   *x;

136:     PetscCall(VecGetArray(xin, &x));
137:     for (PetscInt i = 0; i < n; ++i) x[i] = PetscConj(x[i]);
138:     PetscCall(VecRestoreArray(xin, &x));
139:   }
140:   PetscFunctionReturn(PETSC_SUCCESS);
141: }

143: PetscErrorCode VecResetArray_Seq(Vec vin)
144: {
145:   Vec_Seq *v = (Vec_Seq *)vin->data;

147:   PetscFunctionBegin;
148:   v->array         = v->unplacedarray;
149:   v->unplacedarray = NULL;
150:   PetscFunctionReturn(PETSC_SUCCESS);
151: }

153: PetscErrorCode VecCopy_Seq(Vec xin, Vec yin)
154: {
155:   PetscFunctionBegin;
156:   if (xin != yin) {
157:     const PetscScalar *xa;
158:     PetscScalar       *ya;

160:     PetscCall(VecGetArrayRead(xin, &xa));
161:     PetscCall(VecGetArray(yin, &ya));
162:     PetscCall(PetscArraycpy(ya, xa, xin->map->n));
163:     PetscCall(VecRestoreArrayRead(xin, &xa));
164:     PetscCall(VecRestoreArray(yin, &ya));
165:   }
166:   PetscFunctionReturn(PETSC_SUCCESS);
167: }

169: PetscErrorCode VecSwap_Seq(Vec xin, Vec yin)
170: {
171:   PetscFunctionBegin;
172:   if (xin != yin) {
173:     const PetscBLASInt one = 1;
174:     PetscScalar       *ya, *xa;
175:     PetscBLASInt       bn;

177:     PetscCall(PetscBLASIntCast(xin->map->n, &bn));
178:     PetscCall(VecGetArray(xin, &xa));
179:     PetscCall(VecGetArray(yin, &ya));
180:     PetscCallBLAS("BLASswap", BLASswap_(&bn, xa, &one, ya, &one));
181:     PetscCall(VecRestoreArray(xin, &xa));
182:     PetscCall(VecRestoreArray(yin, &ya));
183:   }
184:   PetscFunctionReturn(PETSC_SUCCESS);
185: }

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

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

195:   PetscFunctionBegin;
196:   if (n) {
197:     const PetscScalar *xx;
198:     const PetscBLASInt one = 1;
199:     PetscBLASInt       bn  = 0;

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

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

241: static PetscErrorCode VecView_Seq_ASCII(Vec xin, PetscViewer viewer)
242: {
243:   PetscInt           i, n = xin->map->n;
244:   const char        *name;
245:   PetscViewerFormat  format;
246:   const PetscScalar *xv;

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

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

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

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

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

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

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

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

442:   PetscFunctionBegin;
443:   PetscCall(PetscViewerDrawGetDraw(v, 0, &draw));
444:   PetscCall(PetscDrawIsNull(draw, &isnull));
445:   if (isnull) PetscFunctionReturn(PETSC_SUCCESS);

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

465:     PetscCall(PetscDrawLGAddPoints(lg, n, &xx, &yy));
466:     PetscCall(PetscDrawLGDraw(lg));
467:     PetscCall(PetscDrawLGSave(lg));
468:   }
469:   PetscCall(VecRestoreArrayRead(xin, &xv));
470:   PetscCall(PetscFree2(xx, yy));
471:   PetscFunctionReturn(PETSC_SUCCESS);
472: }

474: static PetscErrorCode VecView_Seq_Draw(Vec xin, PetscViewer v)
475: {
476:   PetscDraw draw;
477:   PetscBool isnull;

479:   PetscFunctionBegin;
480:   PetscCall(PetscViewerDrawGetDraw(v, 0, &draw));
481:   PetscCall(PetscDrawIsNull(draw, &isnull));
482:   if (isnull) PetscFunctionReturn(PETSC_SUCCESS);

484:   PetscCall(VecView_Seq_Draw_LG(xin, v));
485:   PetscFunctionReturn(PETSC_SUCCESS);
486: }

488: static PetscErrorCode VecView_Seq_Binary(Vec xin, PetscViewer viewer)
489: {
490:   return VecView_Binary(xin, viewer);
491: }

493: #if defined(PETSC_HAVE_MATLAB)
494: #include <petscmatlab.h>
495:   #include <mat.h> /* MATLAB include file */
496: PetscErrorCode VecView_Seq_Matlab(Vec vec, PetscViewer viewer)
497: {
498:   PetscInt           n;
499:   const PetscScalar *array;

501:   PetscFunctionBegin;
502:   PetscCall(VecGetLocalSize(vec, &n));
503:   PetscCall(PetscObjectName((PetscObject)vec));
504:   PetscCall(VecGetArrayRead(vec, &array));
505:   PetscCall(PetscViewerMatlabPutArray(viewer, n, 1, array, ((PetscObject)vec)->name));
506:   PetscCall(VecRestoreArrayRead(vec, &array));
507:   PetscFunctionReturn(PETSC_SUCCESS);
508: }
509: #endif

511: PetscErrorCode VecView_Seq(Vec xin, PetscViewer viewer)
512: {
513:   PetscBool isdraw, iascii, issocket, isbinary;
514: #if defined(PETSC_HAVE_MATHEMATICA)
515:   PetscBool ismathematica;
516: #endif
517: #if defined(PETSC_HAVE_MATLAB)
518:   PetscBool ismatlab;
519: #endif
520: #if defined(PETSC_HAVE_HDF5)
521:   PetscBool ishdf5;
522: #endif
523:   PetscBool isglvis;
524: #if defined(PETSC_HAVE_ADIOS)
525:   PetscBool isadios;
526: #endif

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

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

573: PetscErrorCode VecGetValues_Seq(Vec xin, PetscInt ni, const PetscInt ix[], PetscScalar y[])
574: {
575:   const PetscBool    ignorenegidx = xin->stash.ignorenegidx;
576:   const PetscScalar *xx;

578:   PetscFunctionBegin;
579:   PetscCall(VecGetArrayRead(xin, &xx));
580:   for (PetscInt i = 0; i < ni; ++i) {
581:     if (ignorenegidx && (ix[i] < 0)) continue;
582:     if (PetscDefined(USE_DEBUG)) {
583:       PetscCheck(ix[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Out of range index value %" PetscInt_FMT " cannot be negative", ix[i]);
584:       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);
585:     }
586:     y[i] = xx[ix[i]];
587:   }
588:   PetscCall(VecRestoreArrayRead(xin, &xx));
589:   PetscFunctionReturn(PETSC_SUCCESS);
590: }

592: PetscErrorCode VecSetValues_Seq(Vec xin, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode m)
593: {
594:   const PetscBool ignorenegidx = xin->stash.ignorenegidx;
595:   PetscScalar    *xx;

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

618: PetscErrorCode VecSetValuesBlocked_Seq(Vec xin, PetscInt ni, const PetscInt ix[], const PetscScalar yin[], InsertMode m)
619: {
620:   PetscScalar *xx;
621:   PetscInt     bs;

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

630:     if (start < 0) continue;
631:     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);
632:     for (PetscInt j = 0; j < bs; j++) {
633:       if (m == INSERT_VALUES) {
634:         xx[start + j] = yin[j];
635:       } else {
636:         xx[start + j] += yin[j];
637:       }
638:     }
639:   }
640:   PetscCall(VecRestoreArray(xin, &xx));
641:   PetscFunctionReturn(PETSC_SUCCESS);
642: }

644: static PetscErrorCode VecResetPreallocationCOO_Seq(Vec x)
645: {
646:   Vec_Seq *vs = (Vec_Seq *)x->data;

648:   PetscFunctionBegin;
649:   if (vs) {
650:     PetscCall(PetscFree(vs->jmap1)); /* Destroy old stuff */
651:     PetscCall(PetscFree(vs->perm1));
652:   }
653:   PetscFunctionReturn(PETSC_SUCCESS);
654: }

656: PetscErrorCode VecSetPreallocationCOO_Seq(Vec x, PetscCount coo_n, const PetscInt coo_i[])
657: {
658:   PetscInt    m, *i;
659:   PetscCount  k, nneg;
660:   PetscCount *perm1, *jmap1;
661:   Vec_Seq    *vs = (Vec_Seq *)x->data;

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

675:   PetscCall(VecGetLocalSize(x, &m));
676:   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");
677:   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);

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

684:   if (nneg) { /* Discard leading negative indices */
685:     PetscCount *perm1_new;
686:     PetscCall(PetscMalloc1(coo_n - nneg, &perm1_new));
687:     PetscCall(PetscArraycpy(perm1_new, perm1 + nneg, coo_n - nneg));
688:     PetscCall(PetscFree(perm1));
689:     perm1 = perm1_new;
690:   }

692:   /* Record COO fields */
693:   vs->coo_n = coo_n;
694:   vs->tot1  = coo_n - nneg;
695:   vs->jmap1 = jmap1; /* [m+1] */
696:   vs->perm1 = perm1; /* [tot] */
697:   PetscFunctionReturn(PETSC_SUCCESS);
698: }

700: PetscErrorCode VecSetValuesCOO_Seq(Vec x, const PetscScalar coo_v[], InsertMode imode)
701: {
702:   Vec_Seq          *vs    = (Vec_Seq *)x->data;
703:   const PetscCount *perm1 = vs->perm1, *jmap1 = vs->jmap1;
704:   PetscScalar      *xv;
705:   PetscInt          m;

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

719: PetscErrorCode VecDestroy_Seq(Vec v)
720: {
721:   Vec_Seq *vs = (Vec_Seq *)v->data;

723:   PetscFunctionBegin;
724:   PetscCall(PetscLogObjectState((PetscObject)v, "Length=%" PetscInt_FMT, v->map->n));
725:   if (vs) PetscCall(PetscFree(vs->array_allocated));
726:   PetscCall(VecResetPreallocationCOO_Seq(v));
727:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscMatlabEnginePut_C", NULL));
728:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscMatlabEngineGet_C", NULL));
729:   PetscCall(PetscFree(v->data));
730:   PetscFunctionReturn(PETSC_SUCCESS);
731: }

733: PetscErrorCode VecSetOption_Seq(Vec v, VecOption op, PetscBool flag)
734: {
735:   PetscFunctionBegin;
736:   if (op == VEC_IGNORE_NEGATIVE_INDICES) v->stash.ignorenegidx = flag;
737:   PetscFunctionReturn(PETSC_SUCCESS);
738: }

740: PetscErrorCode VecDuplicate_Seq(Vec win, Vec *V)
741: {
742:   PetscFunctionBegin;
743:   PetscCall(VecCreateWithLayout_Private(win->map, V));
744:   PetscCall(VecSetType(*V, ((PetscObject)win)->type_name));
745:   PetscCall(PetscObjectListDuplicate(((PetscObject)win)->olist, &((PetscObject)(*V))->olist));
746:   PetscCall(PetscFunctionListDuplicate(((PetscObject)win)->qlist, &((PetscObject)(*V))->qlist));

748:   (*V)->ops->view          = win->ops->view;
749:   (*V)->stash.ignorenegidx = win->stash.ignorenegidx;
750:   PetscFunctionReturn(PETSC_SUCCESS);
751: }

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

845: /*
846:       This is called by VecCreate_Seq() (i.e. VecCreateSeq()) and VecCreateSeqWithArray()
847: */
848: PetscErrorCode VecCreate_Seq_Private(Vec v, const PetscScalar array[])
849: {
850:   Vec_Seq *s;

852:   PetscFunctionBegin;
853:   PetscCall(PetscNew(&s));
854:   v->ops[0] = DvOps;

856:   v->data            = (void *)s;
857:   v->petscnative     = PETSC_TRUE;
858:   s->array           = (PetscScalar *)array;
859:   s->array_allocated = NULL;
860:   if (array) v->offloadmask = PETSC_OFFLOAD_CPU;

862:   PetscCall(PetscLayoutSetUp(v->map));
863:   PetscCall(PetscObjectChangeTypeName((PetscObject)v, VECSEQ));
864: #if defined(PETSC_HAVE_MATLAB)
865:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscMatlabEnginePut_C", VecMatlabEnginePut_Default));
866:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscMatlabEngineGet_C", VecMatlabEngineGet_Default));
867: #endif
868:   PetscFunctionReturn(PETSC_SUCCESS);
869: }

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

875:   Collective

877:   Input Parameters:
878: + comm  - the communicator, should be `PETSC_COMM_SELF`
879: . bs    - the block size
880: . n     - the vector length
881: - array - memory where the vector elements are to be stored.

883:   Output Parameter:
884: . V - the vector

886:   Level: intermediate

888:   Notes:
889:   Use `VecDuplicate()` or `VecDuplicateVecs(`) to form additional vectors of the
890:   same type as an existing vector.

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

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

898: .seealso: `VecCreateMPIWithArray()`, `VecCreate()`, `VecDuplicate()`, `VecDuplicateVecs()`,
899:           `VecCreateGhost()`, `VecCreateSeq()`, `VecPlaceArray()`
900: @*/
901: PetscErrorCode VecCreateSeqWithArray(MPI_Comm comm, PetscInt bs, PetscInt n, const PetscScalar array[], Vec *V)
902: {
903:   PetscMPIInt size;

905:   PetscFunctionBegin;
906:   PetscCall(VecCreate(comm, V));
907:   PetscCall(VecSetSizes(*V, n, n));
908:   PetscCall(VecSetBlockSize(*V, bs));
909:   PetscCallMPI(MPI_Comm_size(comm, &size));
910:   PetscCheck(size <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot create VECSEQ on more than one process");
911:   PetscCall(VecCreate_Seq_Private(*V, array));
912:   PetscFunctionReturn(PETSC_SUCCESS);
913: }