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: }