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: PetscErrorCode VecPointwiseMax_Seq(Vec win, Vec xin, Vec yin)
 17: {
 18:   PetscInt     n = win->map->n, i;
 19:   PetscScalar *ww, *xx, *yy; /* cannot make xx or yy const since might be ww */

 21:   VecGetArrayRead(xin, (const PetscScalar **)&xx);
 22:   VecGetArrayRead(yin, (const PetscScalar **)&yy);
 23:   VecGetArray(win, &ww);

 25:   for (i = 0; i < n; i++) ww[i] = PetscMax(PetscRealPart(xx[i]), PetscRealPart(yy[i]));

 27:   VecRestoreArrayRead(xin, (const PetscScalar **)&xx);
 28:   VecRestoreArrayRead(yin, (const PetscScalar **)&yy);
 29:   VecRestoreArray(win, &ww);
 30:   PetscLogFlops(n);
 31:   return 0;
 32: }

 34: PetscErrorCode VecPointwiseMin_Seq(Vec win, Vec xin, Vec yin)
 35: {
 36:   PetscInt     n = win->map->n, i;
 37:   PetscScalar *ww, *xx, *yy; /* cannot make xx or yy const since might be ww */

 39:   VecGetArrayRead(xin, (const PetscScalar **)&xx);
 40:   VecGetArrayRead(yin, (const PetscScalar **)&yy);
 41:   VecGetArray(win, &ww);

 43:   for (i = 0; i < n; i++) ww[i] = PetscMin(PetscRealPart(xx[i]), PetscRealPart(yy[i]));

 45:   VecRestoreArrayRead(xin, (const PetscScalar **)&xx);
 46:   VecRestoreArrayRead(yin, (const PetscScalar **)&yy);
 47:   VecRestoreArray(win, &ww);
 48:   PetscLogFlops(n);
 49:   return 0;
 50: }

 52: PetscErrorCode VecPointwiseMaxAbs_Seq(Vec win, Vec xin, Vec yin)
 53: {
 54:   PetscInt     n = win->map->n, i;
 55:   PetscScalar *ww, *xx, *yy; /* cannot make xx or yy const since might be ww */

 57:   VecGetArrayRead(xin, (const PetscScalar **)&xx);
 58:   VecGetArrayRead(yin, (const PetscScalar **)&yy);
 59:   VecGetArray(win, &ww);

 61:   for (i = 0; i < n; i++) ww[i] = PetscMax(PetscAbsScalar(xx[i]), PetscAbsScalar(yy[i]));

 63:   PetscLogFlops(n);
 64:   VecRestoreArrayRead(xin, (const PetscScalar **)&xx);
 65:   VecRestoreArrayRead(yin, (const PetscScalar **)&yy);
 66:   VecRestoreArray(win, &ww);
 67:   return 0;
 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:   VecGetArrayRead(xin, (const PetscScalar **)&xx);
 78:   VecGetArrayRead(yin, (const PetscScalar **)&yy);
 79:   VecGetArray(win, &ww);
 80:   if (ww == xx) {
 81:     for (i = 0; i < n; i++) ww[i] *= yy[i];
 82:   } else if (ww == yy) {
 83:     for (i = 0; i < n; i++) ww[i] *= xx[i];
 84:   } else {
 85: #if defined(PETSC_USE_FORTRAN_KERNEL_XTIMESY)
 86:     fortranxtimesy_(xx, yy, ww, &n);
 87: #else
 88:     for (i = 0; i < n; i++) ww[i] = xx[i] * yy[i];
 89: #endif
 90:   }
 91:   VecRestoreArrayRead(xin, (const PetscScalar **)&xx);
 92:   VecRestoreArrayRead(yin, (const PetscScalar **)&yy);
 93:   VecRestoreArray(win, &ww);
 94:   PetscLogFlops(n);
 95:   return 0;
 96: }

 98: PetscErrorCode VecPointwiseDivide_Seq(Vec win, Vec xin, Vec yin)
 99: {
100:   PetscInt     n = win->map->n, i;
101:   PetscScalar *ww, *xx, *yy; /* cannot make xx or yy const since might be ww */

103:   VecGetArrayRead(xin, (const PetscScalar **)&xx);
104:   VecGetArrayRead(yin, (const PetscScalar **)&yy);
105:   VecGetArray(win, &ww);

107:   for (i = 0; i < n; i++) {
108:     if (yy[i] != 0.0) ww[i] = xx[i] / yy[i];
109:     else ww[i] = 0.0;
110:   }

112:   PetscLogFlops(n);
113:   VecRestoreArrayRead(xin, (const PetscScalar **)&xx);
114:   VecRestoreArrayRead(yin, (const PetscScalar **)&yy);
115:   VecRestoreArray(win, &ww);
116:   return 0;
117: }

119: PetscErrorCode VecSetRandom_Seq(Vec xin, PetscRandom r)
120: {
121:   PetscInt     n = xin->map->n, i;
122:   PetscScalar *xx;

124:   VecGetArrayWrite(xin, &xx);
125:   for (i = 0; i < n; i++) PetscRandomGetValue(r, &xx[i]);
126:   VecRestoreArrayWrite(xin, &xx);
127:   return 0;
128: }

130: PetscErrorCode VecGetSize_Seq(Vec vin, PetscInt *size)
131: {
132:   *size = vin->map->n;
133:   return 0;
134: }

136: PetscErrorCode VecConjugate_Seq(Vec xin)
137: {
138:   PetscScalar *x;
139:   PetscInt     n = xin->map->n;

141:   VecGetArray(xin, &x);
142:   while (n-- > 0) {
143:     *x = PetscConj(*x);
144:     x++;
145:   }
146:   VecRestoreArray(xin, &x);
147:   return 0;
148: }

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

154:   v->array         = v->unplacedarray;
155:   v->unplacedarray = NULL;
156:   return 0;
157: }

159: PetscErrorCode VecCopy_Seq(Vec xin, Vec yin)
160: {
161:   PetscScalar       *ya;
162:   const PetscScalar *xa;

164:   if (xin != yin) {
165:     VecGetArrayRead(xin, &xa);
166:     VecGetArray(yin, &ya);
167:     PetscArraycpy(ya, xa, xin->map->n);
168:     VecRestoreArrayRead(xin, &xa);
169:     VecRestoreArray(yin, &ya);
170:   }
171:   return 0;
172: }

174: PetscErrorCode VecSwap_Seq(Vec xin, Vec yin)
175: {
176:   PetscScalar *ya, *xa;
177:   PetscBLASInt one = 1, bn;

179:   if (xin != yin) {
180:     PetscBLASIntCast(xin->map->n, &bn);
181:     VecGetArray(xin, &xa);
182:     VecGetArray(yin, &ya);
183:     PetscCallBLAS("BLASswap", BLASswap_(&bn, xa, &one, ya, &one));
184:     VecRestoreArray(xin, &xa);
185:     VecRestoreArray(yin, &ya);
186:   }
187:   return 0;
188: }

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

192: PetscErrorCode VecNorm_Seq(Vec xin, NormType type, PetscReal *z)
193: {
194:   const PetscScalar *xx;
195:   PetscInt           n   = xin->map->n;
196:   PetscBLASInt       one = 1, bn = 0;

198:   PetscBLASIntCast(n, &bn);
199:   if (type == NORM_2 || type == NORM_FROBENIUS) {
200:     VecGetArrayRead(xin, &xx);
201: #if defined(PETSC_USE_REAL___FP16)
202:     PetscCallBLAS("BLASnrm2", *z = BLASnrm2_(&bn, xx, &one));
203: #else
204:     PetscCallBLAS("BLASdot", *z = PetscRealPart(BLASdot_(&bn, xx, &one, xx, &one)));
205:     *z = PetscSqrtReal(*z);
206: #endif
207:     VecRestoreArrayRead(xin, &xx);
208:     PetscLogFlops(PetscMax(2.0 * n - 1, 0.0));
209:   } else if (type == NORM_INFINITY) {
210:     PetscInt  i;
211:     PetscReal max = 0.0, tmp;

213:     VecGetArrayRead(xin, &xx);
214:     for (i = 0; i < n; i++) {
215:       if ((tmp = PetscAbsScalar(*xx)) > max) max = tmp;
216:       /* check special case of tmp == NaN */
217:       if (tmp != tmp) {
218:         max = tmp;
219:         break;
220:       }
221:       xx++;
222:     }
223:     VecRestoreArrayRead(xin, &xx);
224:     *z = max;
225:   } else if (type == NORM_1) {
226: #if defined(PETSC_USE_COMPLEX)
227:     PetscReal tmp = 0.0;
228:     PetscInt  i;
229: #endif
230:     VecGetArrayRead(xin, &xx);
231: #if defined(PETSC_USE_COMPLEX)
232:     /* BLASasum() returns the nonstandard 1 norm of the 1 norm of the complex entries so we provide a custom loop instead */
233:     for (i = 0; i < n; i++) tmp += PetscAbsScalar(xx[i]);
234:     *z = tmp;
235: #else
236:     PetscCallBLAS("BLASasum", *z = BLASasum_(&bn, xx, &one));
237: #endif
238:     VecRestoreArrayRead(xin, &xx);
239:     PetscLogFlops(PetscMax(n - 1.0, 0.0));
240:   } else if (type == NORM_1_AND_2) {
241:     VecNorm_Seq(xin, NORM_1, z);
242:     VecNorm_Seq(xin, NORM_2, z + 1);
243:   }
244:   return 0;
245: }

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

254:   VecGetArrayRead(xin, &xv);
255:   PetscViewerGetFormat(viewer, &format);
256:   if (format == PETSC_VIEWER_ASCII_MATLAB) {
257:     PetscObjectGetName((PetscObject)xin, &name);
258:     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:         PetscViewerASCIIPrintf(viewer, "%18.16e + %18.16ei\n", (double)PetscRealPart(xv[i]), (double)PetscImaginaryPart(xv[i]));
263:       } else if (PetscImaginaryPart(xv[i]) < 0.0) {
264:         PetscViewerASCIIPrintf(viewer, "%18.16e - %18.16ei\n", (double)PetscRealPart(xv[i]), -(double)PetscImaginaryPart(xv[i]));
265:       } else {
266:         PetscViewerASCIIPrintf(viewer, "%18.16e\n", (double)PetscRealPart(xv[i]));
267:       }
268: #else
269:       PetscViewerASCIIPrintf(viewer, "%18.16e\n", (double)xv[i]);
270: #endif
271:     }
272:     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:       PetscViewerASCIIPrintf(viewer, "%18.16e %18.16e\n", (double)PetscRealPart(xv[i]), (double)PetscImaginaryPart(xv[i]));
277: #else
278:       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) PetscObjectComposedDataRegister(&stateId);
296:     PetscObjectComposedDataGetInt((PetscObject)viewer, stateId, outputState, hasState);
297:     if (!hasState) outputState = 0;
298:     PetscObjectGetName((PetscObject)xin, &name);
299:     VecGetBlockSize(xin, &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;

312:       if (doOutput) 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 {
324:         if (outputState == 4) doOutput = 0;
325:       }

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

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

367:     VecGetBlockSize(xin, &bs);
369:     PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT "\n", xin->map->N / bs);
370:     for (i = 0; i < n / bs; i++) {
371:       PetscViewerASCIIPrintf(viewer, "%7" PetscInt_FMT "   ", i + 1);
372:       for (b = 0; b < bs; b++) {
373:         if (b > 0) PetscViewerASCIIPrintf(viewer, " ");
374: #if !defined(PETSC_USE_COMPLEX)
375:         PetscViewerASCIIPrintf(viewer, "% 12.5E", (double)xv[i * bs + b]);
376: #endif
377:       }
378:       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:     VecGetBlockSize(xin, &vdim);
390:     PetscViewerASCIIPrintf(viewer, "FiniteElementSpace\n");
391:     PetscObjectQuery((PetscObject)xin, "_glvis_info_container", (PetscObject *)&glvis_container);
393:     PetscContainerGetPointer(glvis_container, (void **)&glvis_vec_info);
394:     PetscViewerASCIIPrintf(viewer, "%s\n", glvis_vec_info->fec_type);
395:     PetscViewerASCIIPrintf(viewer, "VDim: %" PetscInt_FMT "\n", vdim);
396:     PetscViewerASCIIPrintf(viewer, "Ordering: %" PetscInt_FMT "\n", ordering);
397:     PetscViewerASCIIPrintf(viewer, "\n");
398:     /* mfem::Vector::Print() */
399:     PetscObjectQuery((PetscObject)viewer, "_glvis_info_container", (PetscObject *)&glvis_container);
401:     PetscContainerGetPointer(glvis_container, (void **)&glvis_info);
402:     if (glvis_info->enabled) {
403:       VecGetLocalSize(xin, &n);
404:       VecGetArrayRead(xin, &array);
405:       for (i = 0; i < n; i++) {
406:         PetscViewerASCIIPrintf(viewer, glvis_info->fmt, (double)PetscRealPart(array[i]));
407:         PetscViewerASCIIPrintf(viewer, "\n");
408:       }
409:       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) PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT ": ", i);
416: #if defined(PETSC_USE_COMPLEX)
417:       if (PetscImaginaryPart(xv[i]) > 0.0) {
418:         PetscViewerASCIIPrintf(viewer, "%g + %g i\n", (double)PetscRealPart(xv[i]), (double)PetscImaginaryPart(xv[i]));
419:       } else if (PetscImaginaryPart(xv[i]) < 0.0) {
420:         PetscViewerASCIIPrintf(viewer, "%g - %g i\n", (double)PetscRealPart(xv[i]), -(double)PetscImaginaryPart(xv[i]));
421:       } else {
422:         PetscViewerASCIIPrintf(viewer, "%g\n", (double)PetscRealPart(xv[i]));
423:       }
424: #else
425:       PetscViewerASCIIPrintf(viewer, "%g\n", (double)xv[i]);
426: #endif
427:     }
428:   }
429:   PetscViewerFlush(viewer);
430:   VecRestoreArrayRead(xin, &xv);
431:   return 0;
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:   PetscViewerDrawGetDraw(v, 0, &draw);
448:   PetscDrawIsNull(draw, &isnull);
449:   if (isnull) return 0;

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

469:     PetscDrawLGAddPoints(lg, n, &xx, &yy);
470:     PetscDrawLGDraw(lg);
471:     PetscDrawLGSave(lg);
472:   }
473:   VecRestoreArrayRead(xin, &xv);
474:   PetscFree2(xx, yy);
475:   return 0;
476: }

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

483:   PetscViewerDrawGetDraw(v, 0, &draw);
484:   PetscDrawIsNull(draw, &isnull);
485:   if (isnull) return 0;

487:   VecView_Seq_Draw_LG(xin, v);
488:   return 0;
489: }

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

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

504:   VecGetLocalSize(vec, &n);
505:   PetscObjectName((PetscObject)vec);
506:   VecGetArrayRead(vec, &array);
507:   PetscViewerMatlabPutArray(viewer, n, 1, array, ((PetscObject)vec)->name);
508:   VecRestoreArrayRead(vec, &array);
509:   return 0;
510: }
511: #endif

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

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

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

574: PetscErrorCode VecGetValues_Seq(Vec xin, PetscInt ni, const PetscInt ix[], PetscScalar y[])
575: {
576:   const PetscScalar *xx;
577:   PetscInt           i;

579:   VecGetArrayRead(xin, &xx);
580:   for (i = 0; i < ni; i++) {
581:     if (xin->stash.ignorenegidx && ix[i] < 0) continue;
582:     if (PetscDefined(USE_DEBUG)) {
585:     }
586:     y[i] = xx[ix[i]];
587:   }
588:   VecRestoreArrayRead(xin, &xx);
589:   return 0;
590: }

592: PetscErrorCode VecSetValues_Seq(Vec xin, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode m)
593: {
594:   PetscScalar *xx;
595:   PetscInt     i;

597:   VecGetArray(xin, &xx);
598:   if (m == INSERT_VALUES) {
599:     for (i = 0; i < ni; i++) {
600:       if (xin->stash.ignorenegidx && ix[i] < 0) continue;
601:       if (PetscDefined(USE_DEBUG)) {
604:       }
605:       xx[ix[i]] = y[i];
606:     }
607:   } else {
608:     for (i = 0; i < ni; i++) {
609:       if (xin->stash.ignorenegidx && ix[i] < 0) continue;
610:       if (PetscDefined(USE_DEBUG)) {
613:       }
614:       xx[ix[i]] += y[i];
615:     }
616:   }
617:   VecRestoreArray(xin, &xx);
618:   return 0;
619: }

621: PetscErrorCode VecSetValuesBlocked_Seq(Vec xin, PetscInt ni, const PetscInt ix[], const PetscScalar yin[], InsertMode m)
622: {
623:   PetscScalar *xx, *y = (PetscScalar *)yin;
624:   PetscInt     i, bs, start, j;

626:   /*
627:        For optimization could treat bs = 2, 3, 4, 5 as special cases with loop unrolling
628:   */
629:   VecGetBlockSize(xin, &bs);
630:   VecGetArray(xin, &xx);
631:   if (m == INSERT_VALUES) {
632:     for (i = 0; i < ni; i++, y += bs) {
633:       start = bs * ix[i];
634:       if (start < 0) continue;
636:       for (j = 0; j < bs; j++) xx[start + j] = y[j];
637:     }
638:   } else {
639:     for (i = 0; i < ni; i++, y += bs) {
640:       start = bs * ix[i];
641:       if (start < 0) continue;
643:       for (j = 0; j < bs; j++) xx[start + j] += y[j];
644:     }
645:   }
646:   VecRestoreArray(xin, &xx);
647:   return 0;
648: }

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

654:   if (vs) {
655:     PetscFree(vs->jmap1); /* Destroy old stuff */
656:     PetscFree(vs->perm1);
657:   }
658:   return 0;
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:   VecResetPreallocationCOO_Seq(x); /* Destroy old stuff */
669:   PetscMalloc1(coo_n, &i);
670:   PetscArraycpy(i, coo_i, coo_n); /* Make a copy since we'll modify it */
671:   PetscMalloc1(coo_n, &perm1);
672:   for (k = 0; k < coo_n; k++) perm1[k] = k;
673:   PetscSortIntWithCountArray(coo_n, i, perm1);
674:   for (k = 0; k < coo_n; k++) {
675:     if (i[k] >= 0) break;
676:   } /* Advance k to the first entry with a non-negative index */
677:   nneg = k;

679:   VecGetLocalSize(x, &m);

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

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

696:   /* Record COO fields */
697:   vs->coo_n = coo_n;
698:   vs->tot1  = coo_n - nneg;
699:   vs->jmap1 = jmap1; /* [m+1] */
700:   vs->perm1 = perm1; /* [tot] */
701:   return 0;
702: }

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

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

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

726: #if defined(PETSC_USE_LOG)
727:   PetscLogObjectState((PetscObject)v, "Length=%" PetscInt_FMT, v->map->n);
728: #endif
729:   if (vs) PetscFree(vs->array_allocated);
730:   VecResetPreallocationCOO_Seq(v);
731:   PetscObjectComposeFunction((PetscObject)v, "PetscMatlabEnginePut_C", NULL);
732:   PetscObjectComposeFunction((PetscObject)v, "PetscMatlabEngineGet_C", NULL);
733:   PetscFree(v->data);
734:   return 0;
735: }

737: PetscErrorCode VecSetOption_Seq(Vec v, VecOption op, PetscBool flag)
738: {
739:   if (op == VEC_IGNORE_NEGATIVE_INDICES) v->stash.ignorenegidx = flag;
740:   return 0;
741: }

743: PetscErrorCode VecDuplicate_Seq(Vec win, Vec *V)
744: {
745:   VecCreate(PetscObjectComm((PetscObject)win), V);
746:   VecSetSizes(*V, win->map->n, win->map->n);
747:   VecSetType(*V, ((PetscObject)win)->type_name);
748:   PetscLayoutReference(win->map, &(*V)->map);
749:   PetscObjectListDuplicate(((PetscObject)win)->olist, &((PetscObject)(*V))->olist);
750:   PetscFunctionListDuplicate(((PetscObject)win)->qlist, &((PetscObject)(*V))->qlist);

752:   (*V)->ops->view          = win->ops->view;
753:   (*V)->stash.ignorenegidx = win->stash.ignorenegidx;
754:   return 0;
755: }

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

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

855:   PetscNew(&s);
856:   PetscMemcpy(v->ops, &DvOps, sizeof(DvOps));

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

864:   PetscLayoutSetUp(v->map);
865:   PetscObjectChangeTypeName((PetscObject)v, VECSEQ);
866: #if defined(PETSC_HAVE_MATLAB)
867:   PetscObjectComposeFunction((PetscObject)v, "PetscMatlabEnginePut_C", VecMatlabEnginePut_Default);
868:   PetscObjectComposeFunction((PetscObject)v, "PetscMatlabEngineGet_C", VecMatlabEngineGet_Default);
869: #endif
870:   return 0;
871: }

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

877:    Collective

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

885:    Output Parameter:
886: .  V - the vector

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:    Level: intermediate

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

907:   VecCreate(comm, V);
908:   VecSetSizes(*V, n, n);
909:   VecSetBlockSize(*V, bs);
910:   MPI_Comm_size(comm, &size);
912:   VecCreate_Seq_Private(*V, array);
913:   return 0;
914: }