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 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: const PetscInt n = xin->map->n;
132: PetscScalar *x;
134: PetscFunctionBegin;
135: PetscCall(VecGetArray(xin, &x));
136: for (PetscInt i = 0; i < n; ++i) x[i] = PetscConj(x[i]);
137: PetscCall(VecRestoreArray(xin, &x));
138: PetscFunctionReturn(PETSC_SUCCESS);
139: }
141: PetscErrorCode VecResetArray_Seq(Vec vin)
142: {
143: Vec_Seq *v = (Vec_Seq *)vin->data;
145: PetscFunctionBegin;
146: v->array = v->unplacedarray;
147: v->unplacedarray = NULL;
148: PetscFunctionReturn(PETSC_SUCCESS);
149: }
151: PetscErrorCode VecCopy_Seq(Vec xin, Vec yin)
152: {
153: PetscFunctionBegin;
154: if (xin != yin) {
155: const PetscScalar *xa;
156: PetscScalar *ya;
158: PetscCall(VecGetArrayRead(xin, &xa));
159: PetscCall(VecGetArrayWrite(yin, &ya));
160: PetscCall(PetscArraycpy(ya, xa, xin->map->n));
161: PetscCall(VecRestoreArrayRead(xin, &xa));
162: PetscCall(VecRestoreArrayWrite(yin, &ya));
163: }
164: PetscFunctionReturn(PETSC_SUCCESS);
165: }
167: PetscErrorCode VecSwap_Seq(Vec xin, Vec yin)
168: {
169: PetscFunctionBegin;
170: if (xin != yin) {
171: const PetscBLASInt one = 1;
172: PetscScalar *ya, *xa;
173: PetscBLASInt bn;
175: PetscCall(PetscBLASIntCast(xin->map->n, &bn));
176: PetscCall(VecGetArray(xin, &xa));
177: PetscCall(VecGetArray(yin, &ya));
178: PetscCallBLAS("BLASswap", BLASswap_(&bn, xa, &one, ya, &one));
179: PetscCall(VecRestoreArray(xin, &xa));
180: PetscCall(VecRestoreArray(yin, &ya));
181: }
182: PetscFunctionReturn(PETSC_SUCCESS);
183: }
185: PetscErrorCode VecNorm_Seq(Vec xin, NormType type, PetscReal *z)
186: {
187: // use a local variable to ensure compiler doesn't think z aliases any of the other arrays
188: PetscReal ztmp[] = {0.0, 0.0};
189: const PetscInt n = xin->map->n;
191: PetscFunctionBegin;
192: if (n) {
193: const PetscScalar *xx;
194: const PetscBLASInt one = 1;
195: PetscBLASInt bn = 0;
197: PetscCall(PetscBLASIntCast(n, &bn));
198: PetscCall(VecGetArrayRead(xin, &xx));
199: if (type == NORM_2 || type == NORM_FROBENIUS) {
200: NORM_1_AND_2_DOING_NORM_2:
201: if (PetscDefined(USE_REAL___FP16)) {
202: PetscCallBLAS("BLASnrm2", ztmp[type == NORM_1_AND_2] = BLASnrm2_(&bn, xx, &one));
203: } else {
204: PetscCallBLAS("BLASdot", ztmp[type == NORM_1_AND_2] = PetscSqrtReal(PetscRealPart(BLASdot_(&bn, xx, &one, xx, &one))));
205: }
206: PetscCall(PetscLogFlops(2.0 * n - 1));
207: } else if (type == NORM_INFINITY) {
208: for (PetscInt i = 0; i < n; ++i) {
209: const PetscReal tmp = PetscAbsScalar(xx[i]);
211: /* check special case of tmp == NaN */
212: if ((tmp > ztmp[0]) || (tmp != tmp)) {
213: ztmp[0] = tmp;
214: if (tmp != tmp) break;
215: }
216: }
217: } else if (type == NORM_1 || type == NORM_1_AND_2) {
218: if (PetscDefined(USE_COMPLEX)) {
219: // BLASasum() returns the nonstandard 1 norm of the 1 norm of the complex entries so we
220: // provide a custom loop instead
221: for (PetscInt i = 0; i < n; ++i) ztmp[0] += PetscAbsScalar(xx[i]);
222: } else {
223: PetscCallBLAS("BLASasum", ztmp[0] = BLASasum_(&bn, xx, &one));
224: }
225: PetscCall(PetscLogFlops(n - 1.0));
226: /* slight reshuffle so we can skip getting the array again (but still log the flops) if we
227: do norm2 after this */
228: if (type == NORM_1_AND_2) goto NORM_1_AND_2_DOING_NORM_2;
229: }
230: PetscCall(VecRestoreArrayRead(xin, &xx));
231: }
232: z[0] = ztmp[0];
233: if (type == NORM_1_AND_2) z[1] = ztmp[1];
234: PetscFunctionReturn(PETSC_SUCCESS);
235: }
237: static PetscErrorCode VecView_Seq_ASCII(Vec xin, PetscViewer viewer)
238: {
239: PetscInt i, n = xin->map->n;
240: const char *name;
241: PetscViewerFormat format;
242: const PetscScalar *xv;
244: PetscFunctionBegin;
245: PetscCall(VecGetArrayRead(xin, &xv));
246: PetscCall(PetscViewerGetFormat(viewer, &format));
247: if (format == PETSC_VIEWER_ASCII_MATLAB) {
248: PetscCall(PetscObjectGetName((PetscObject)xin, &name));
249: PetscCall(PetscViewerASCIIPrintf(viewer, "%s = [\n", name));
250: for (i = 0; i < n; i++) {
251: #if defined(PETSC_USE_COMPLEX)
252: if (PetscImaginaryPart(xv[i]) > 0.0) {
253: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e + %18.16ei\n", (double)PetscRealPart(xv[i]), (double)PetscImaginaryPart(xv[i])));
254: } else if (PetscImaginaryPart(xv[i]) < 0.0) {
255: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e - %18.16ei\n", (double)PetscRealPart(xv[i]), -(double)PetscImaginaryPart(xv[i])));
256: } else {
257: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e\n", (double)PetscRealPart(xv[i])));
258: }
259: #else
260: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e\n", (double)xv[i]));
261: #endif
262: }
263: PetscCall(PetscViewerASCIIPrintf(viewer, "];\n"));
264: } else if (format == PETSC_VIEWER_ASCII_SYMMODU) {
265: for (i = 0; i < n; i++) {
266: #if defined(PETSC_USE_COMPLEX)
267: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e %18.16e\n", (double)PetscRealPart(xv[i]), (double)PetscImaginaryPart(xv[i])));
268: #else
269: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e\n", (double)xv[i]));
270: #endif
271: }
272: } else if (format == PETSC_VIEWER_ASCII_PCICE) {
273: PetscInt bs, b;
275: PetscCall(VecGetBlockSize(xin, &bs));
276: 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);
277: PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT "\n", xin->map->N / bs));
278: for (i = 0; i < n / bs; i++) {
279: PetscCall(PetscViewerASCIIPrintf(viewer, "%7" PetscInt_FMT " ", i + 1));
280: for (b = 0; b < bs; b++) {
281: if (b > 0) PetscCall(PetscViewerASCIIPrintf(viewer, " "));
282: #if !defined(PETSC_USE_COMPLEX)
283: PetscCall(PetscViewerASCIIPrintf(viewer, "% 12.5E", (double)xv[i * bs + b]));
284: #endif
285: }
286: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
287: }
288: } else if (format == PETSC_VIEWER_ASCII_GLVIS) {
289: /* GLVis ASCII visualization/dump: this function mimics mfem::GridFunction::Save() */
290: const PetscScalar *array;
291: PetscInt i, n, vdim, ordering = 1; /* mfem::FiniteElementSpace::Ordering::byVDIM */
292: PetscContainer glvis_container;
293: PetscViewerGLVisVecInfo glvis_vec_info;
294: PetscViewerGLVisInfo glvis_info;
296: /* mfem::FiniteElementSpace::Save() */
297: PetscCall(VecGetBlockSize(xin, &vdim));
298: PetscCall(PetscViewerASCIIPrintf(viewer, "FiniteElementSpace\n"));
299: PetscCall(PetscObjectQuery((PetscObject)xin, "_glvis_info_container", (PetscObject *)&glvis_container));
300: PetscCheck(glvis_container, PetscObjectComm((PetscObject)xin), PETSC_ERR_PLIB, "Missing GLVis container");
301: PetscCall(PetscContainerGetPointer(glvis_container, &glvis_vec_info));
302: PetscCall(PetscViewerASCIIPrintf(viewer, "%s\n", glvis_vec_info->fec_type));
303: PetscCall(PetscViewerASCIIPrintf(viewer, "VDim: %" PetscInt_FMT "\n", vdim));
304: PetscCall(PetscViewerASCIIPrintf(viewer, "Ordering: %" PetscInt_FMT "\n", ordering));
305: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
306: /* mfem::Vector::Print() */
307: PetscCall(PetscObjectQuery((PetscObject)viewer, "_glvis_info_container", (PetscObject *)&glvis_container));
308: PetscCheck(glvis_container, PetscObjectComm((PetscObject)viewer), PETSC_ERR_PLIB, "Missing GLVis container");
309: PetscCall(PetscContainerGetPointer(glvis_container, &glvis_info));
310: if (glvis_info->enabled) {
311: PetscCall(VecGetLocalSize(xin, &n));
312: PetscCall(VecGetArrayRead(xin, &array));
313: for (i = 0; i < n; i++) {
314: PetscCall(PetscViewerASCIIPrintf(viewer, glvis_info->fmt, (double)PetscRealPart(array[i])));
315: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
316: }
317: PetscCall(VecRestoreArrayRead(xin, &array));
318: }
319: } else if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
320: /* No info */
321: } else {
322: for (i = 0; i < n; i++) {
323: if (format == PETSC_VIEWER_ASCII_INDEX) PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT ": ", i));
324: #if defined(PETSC_USE_COMPLEX)
325: if (PetscImaginaryPart(xv[i]) > 0.0) {
326: PetscCall(PetscViewerASCIIPrintf(viewer, "%g + %g i\n", (double)PetscRealPart(xv[i]), (double)PetscImaginaryPart(xv[i])));
327: } else if (PetscImaginaryPart(xv[i]) < 0.0) {
328: PetscCall(PetscViewerASCIIPrintf(viewer, "%g - %g i\n", (double)PetscRealPart(xv[i]), -(double)PetscImaginaryPart(xv[i])));
329: } else {
330: PetscCall(PetscViewerASCIIPrintf(viewer, "%g\n", (double)PetscRealPart(xv[i])));
331: }
332: #else
333: PetscCall(PetscViewerASCIIPrintf(viewer, "%g\n", (double)xv[i]));
334: #endif
335: }
336: }
337: PetscCall(PetscViewerFlush(viewer));
338: PetscCall(VecRestoreArrayRead(xin, &xv));
339: PetscFunctionReturn(PETSC_SUCCESS);
340: }
342: #include <petscdraw.h>
343: static PetscErrorCode VecView_Seq_Draw_LG(Vec xin, PetscViewer v)
344: {
345: PetscDraw draw;
346: PetscBool isnull;
347: PetscDrawLG lg;
348: PetscInt i, c, bs = xin->map->bs, n = xin->map->n / bs;
349: const PetscScalar *xv;
350: PetscReal *xx, *yy, xmin, xmax, h;
351: int colors[] = {PETSC_DRAW_RED};
352: PetscViewerFormat format;
353: PetscDrawAxis axis;
354: const char *name;
356: PetscFunctionBegin;
357: PetscCall(PetscViewerDrawGetDraw(v, 0, &draw));
358: PetscCall(PetscDrawIsNull(draw, &isnull));
359: if (isnull) PetscFunctionReturn(PETSC_SUCCESS);
361: PetscCall(PetscObjectGetName((PetscObject)xin, &name));
362: PetscCall(PetscDrawSetTitle(draw, name));
363: PetscCall(PetscViewerGetFormat(v, &format));
364: PetscCall(PetscMalloc2(n, &xx, n, &yy));
365: PetscCall(VecGetArrayRead(xin, &xv));
366: for (c = 0; c < bs; c++) {
367: PetscCall(PetscViewerDrawGetDrawLG(v, c, &lg));
368: PetscCall(PetscDrawLGReset(lg));
369: PetscCall(PetscDrawLGSetDimension(lg, 1));
370: PetscCall(PetscDrawLGSetColors(lg, colors));
371: if (format == PETSC_VIEWER_DRAW_LG_XRANGE) {
372: PetscCall(PetscDrawLGGetAxis(lg, &axis));
373: PetscCall(PetscDrawAxisGetLimits(axis, &xmin, &xmax, NULL, NULL));
374: h = (xmax - xmin) / n;
375: for (i = 0; i < n; i++) xx[i] = i * h + 0.5 * h; /* cell center */
376: } else {
377: for (i = 0; i < n; i++) xx[i] = (PetscReal)i;
378: }
379: for (i = 0; i < n; i++) yy[i] = PetscRealPart(xv[c + i * bs]);
381: PetscCall(PetscDrawLGAddPoints(lg, n, &xx, &yy));
382: PetscCall(PetscDrawLGDraw(lg));
383: PetscCall(PetscDrawLGSave(lg));
384: }
385: PetscCall(VecRestoreArrayRead(xin, &xv));
386: PetscCall(PetscFree2(xx, yy));
387: PetscFunctionReturn(PETSC_SUCCESS);
388: }
390: static PetscErrorCode VecView_Seq_Draw(Vec xin, PetscViewer v)
391: {
392: PetscDraw draw;
393: PetscBool isnull;
395: PetscFunctionBegin;
396: PetscCall(PetscViewerDrawGetDraw(v, 0, &draw));
397: PetscCall(PetscDrawIsNull(draw, &isnull));
398: if (isnull) PetscFunctionReturn(PETSC_SUCCESS);
400: PetscCall(VecView_Seq_Draw_LG(xin, v));
401: PetscFunctionReturn(PETSC_SUCCESS);
402: }
404: static PetscErrorCode VecView_Seq_Binary(Vec xin, PetscViewer viewer)
405: {
406: return VecView_Binary(xin, viewer);
407: }
409: #if defined(PETSC_HAVE_MATLAB)
410: #include <petscmatlab.h>
411: #include <mat.h> /* MATLAB include file */
412: PetscErrorCode VecView_Seq_Matlab(Vec vec, PetscViewer viewer)
413: {
414: PetscInt n;
415: const PetscScalar *array;
417: PetscFunctionBegin;
418: PetscCall(VecGetLocalSize(vec, &n));
419: PetscCall(PetscObjectName((PetscObject)vec));
420: PetscCall(VecGetArrayRead(vec, &array));
421: PetscCall(PetscViewerMatlabPutArray(viewer, n, 1, array, ((PetscObject)vec)->name));
422: PetscCall(VecRestoreArrayRead(vec, &array));
423: PetscFunctionReturn(PETSC_SUCCESS);
424: }
425: #endif
427: PetscErrorCode VecView_Seq(Vec xin, PetscViewer viewer)
428: {
429: PetscBool isdraw, isascii, issocket, isbinary;
430: #if defined(PETSC_HAVE_MATHEMATICA)
431: PetscBool ismathematica;
432: #endif
433: #if defined(PETSC_HAVE_MATLAB)
434: PetscBool ismatlab;
435: #endif
436: #if defined(PETSC_HAVE_HDF5)
437: PetscBool ishdf5;
438: #endif
439: PetscBool isglvis;
440: #if defined(PETSC_HAVE_ADIOS)
441: PetscBool isadios;
442: #endif
444: PetscFunctionBegin;
445: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
446: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
447: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSOCKET, &issocket));
448: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
449: #if defined(PETSC_HAVE_MATHEMATICA)
450: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERMATHEMATICA, &ismathematica));
451: #endif
452: #if defined(PETSC_HAVE_HDF5)
453: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
454: #endif
455: #if defined(PETSC_HAVE_MATLAB)
456: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERMATLAB, &ismatlab));
457: #endif
458: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERGLVIS, &isglvis));
459: #if defined(PETSC_HAVE_ADIOS)
460: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERADIOS, &isadios));
461: #endif
463: if (isdraw) {
464: PetscCall(VecView_Seq_Draw(xin, viewer));
465: } else if (isascii) {
466: PetscCall(VecView_Seq_ASCII(xin, viewer));
467: } else if (isbinary) {
468: PetscCall(VecView_Seq_Binary(xin, viewer));
469: #if defined(PETSC_HAVE_MATHEMATICA)
470: } else if (ismathematica) {
471: PetscCall(PetscViewerMathematicaPutVector(viewer, xin));
472: #endif
473: #if defined(PETSC_HAVE_HDF5)
474: } else if (ishdf5) {
475: PetscCall(VecView_MPI_HDF5(xin, viewer)); /* Reusing VecView_MPI_HDF5 ... don't want code duplication*/
476: #endif
477: #if defined(PETSC_HAVE_ADIOS)
478: } else if (isadios) {
479: PetscCall(VecView_MPI_ADIOS(xin, viewer)); /* Reusing VecView_MPI_ADIOS ... don't want code duplication*/
480: #endif
481: #if defined(PETSC_HAVE_MATLAB)
482: } else if (ismatlab) {
483: PetscCall(VecView_Seq_Matlab(xin, viewer));
484: #endif
485: } else if (isglvis) PetscCall(VecView_GLVis(xin, viewer));
486: PetscFunctionReturn(PETSC_SUCCESS);
487: }
489: PetscErrorCode VecGetValues_Seq(Vec xin, PetscInt ni, const PetscInt ix[], PetscScalar y[])
490: {
491: const PetscBool ignorenegidx = xin->stash.ignorenegidx;
492: const PetscScalar *xx;
494: PetscFunctionBegin;
495: PetscCall(VecGetArrayRead(xin, &xx));
496: for (PetscInt i = 0; i < ni; ++i) {
497: if (ignorenegidx && (ix[i] < 0)) continue;
498: if (PetscDefined(USE_DEBUG)) {
499: PetscCheck(ix[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Out of range index value %" PetscInt_FMT " cannot be negative", ix[i]);
500: PetscCheck(ix[i] < xin->map->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Out of range index value %" PetscInt_FMT ", should be less than %" PetscInt_FMT, ix[i], xin->map->n);
501: }
502: y[i] = xx[ix[i]];
503: }
504: PetscCall(VecRestoreArrayRead(xin, &xx));
505: PetscFunctionReturn(PETSC_SUCCESS);
506: }
508: PetscErrorCode VecSetValues_Seq(Vec xin, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode m)
509: {
510: const PetscBool ignorenegidx = xin->stash.ignorenegidx;
511: PetscScalar *xx;
513: PetscFunctionBegin;
514: // call to getarray (not e.g. getarraywrite() if m is INSERT_VALUES) is deliberate! If this
515: // is secretly a VECSEQCUDA it may have values currently on the device, in which case --
516: // unless we are replacing the entire array -- we need to copy them up
517: PetscCall(VecGetArray(xin, &xx));
518: for (PetscInt i = 0; i < ni; i++) {
519: if (ignorenegidx && (ix[i] < 0)) continue;
520: PetscScalar yv = y ? y[i] : 0.0;
521: if (PetscDefined(USE_DEBUG)) {
522: PetscCheck(ix[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Out of range index value %" PetscInt_FMT " cannot be negative", ix[i]);
523: PetscCheck(ix[i] < xin->map->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Out of range index value %" PetscInt_FMT ", should be less than %" PetscInt_FMT, ix[i], xin->map->n);
524: }
525: if (m == INSERT_VALUES) {
526: xx[ix[i]] = yv;
527: } else {
528: xx[ix[i]] += yv;
529: }
530: }
531: PetscCall(VecRestoreArray(xin, &xx));
532: PetscFunctionReturn(PETSC_SUCCESS);
533: }
535: PetscErrorCode VecSetValuesBlocked_Seq(Vec xin, PetscInt ni, const PetscInt ix[], const PetscScalar yin[], InsertMode m)
536: {
537: PetscScalar *xx;
538: PetscInt bs;
540: /* For optimization could treat bs = 2, 3, 4, 5 as special cases with loop unrolling */
541: PetscFunctionBegin;
542: PetscCall(VecGetBlockSize(xin, &bs));
543: PetscCall(VecGetArray(xin, &xx));
544: for (PetscInt i = 0; i < ni; ++i, yin += bs) {
545: const PetscInt start = bs * ix[i];
547: if (start < 0) continue;
548: PetscCheck(start < xin->map->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Out of range index value %" PetscInt_FMT ", should be less than %" PetscInt_FMT, start, xin->map->n);
549: for (PetscInt j = 0; j < bs; j++) {
550: if (m == INSERT_VALUES) {
551: xx[start + j] = yin ? yin[j] : 0.0;
552: } else {
553: xx[start + j] += yin ? yin[j] : 0.0;
554: }
555: }
556: }
557: PetscCall(VecRestoreArray(xin, &xx));
558: PetscFunctionReturn(PETSC_SUCCESS);
559: }
561: static PetscErrorCode VecResetPreallocationCOO_Seq(Vec x)
562: {
563: Vec_Seq *vs = (Vec_Seq *)x->data;
565: PetscFunctionBegin;
566: if (vs) {
567: PetscCall(PetscFree(vs->jmap1)); /* Destroy old stuff */
568: PetscCall(PetscFree(vs->perm1));
569: }
570: PetscFunctionReturn(PETSC_SUCCESS);
571: }
573: PetscErrorCode VecSetPreallocationCOO_Seq(Vec x, PetscCount coo_n, const PetscInt coo_i[])
574: {
575: PetscInt m, *i;
576: PetscCount k, nneg;
577: PetscCount *perm1, *jmap1;
578: Vec_Seq *vs = (Vec_Seq *)x->data;
580: PetscFunctionBegin;
581: PetscCall(VecResetPreallocationCOO_Seq(x)); /* Destroy old stuff */
582: PetscCall(PetscMalloc1(coo_n, &i));
583: PetscCall(PetscArraycpy(i, coo_i, coo_n)); /* Make a copy since we'll modify it */
584: PetscCall(PetscMalloc1(coo_n, &perm1));
585: for (k = 0; k < coo_n; k++) perm1[k] = k;
586: PetscCall(PetscSortIntWithCountArray(coo_n, i, perm1));
587: for (k = 0; k < coo_n; k++) {
588: if (i[k] >= 0) break;
589: } /* Advance k to the first entry with a non-negative index */
590: nneg = k;
592: PetscCall(VecGetLocalSize(x, &m));
593: 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");
594: 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);
596: PetscCall(PetscCalloc1(m + 1, &jmap1));
597: for (; k < coo_n; k++) jmap1[i[k] + 1]++; /* Count repeats of each entry */
598: for (k = 0; k < m; k++) jmap1[k + 1] += jmap1[k]; /* Transform jmap[] to CSR-like data structure */
599: PetscCall(PetscFree(i));
601: if (nneg) { /* Discard leading negative indices */
602: PetscCount *perm1_new;
603: PetscCall(PetscMalloc1(coo_n - nneg, &perm1_new));
604: PetscCall(PetscArraycpy(perm1_new, perm1 + nneg, coo_n - nneg));
605: PetscCall(PetscFree(perm1));
606: perm1 = perm1_new;
607: }
609: /* Record COO fields */
610: vs->coo_n = coo_n;
611: vs->tot1 = coo_n - nneg;
612: vs->jmap1 = jmap1; /* [m+1] */
613: vs->perm1 = perm1; /* [tot] */
614: PetscFunctionReturn(PETSC_SUCCESS);
615: }
617: PetscErrorCode VecSetValuesCOO_Seq(Vec x, const PetscScalar coo_v[], InsertMode imode)
618: {
619: Vec_Seq *vs = (Vec_Seq *)x->data;
620: const PetscCount *perm1 = vs->perm1, *jmap1 = vs->jmap1;
621: PetscScalar *xv;
622: PetscInt m;
624: PetscFunctionBegin;
625: PetscCall(VecGetLocalSize(x, &m));
626: PetscCall(VecGetArray(x, &xv));
627: for (PetscInt i = 0; i < m; i++) {
628: PetscScalar sum = 0.0;
629: for (PetscCount j = jmap1[i]; j < jmap1[i + 1]; j++) sum += coo_v[perm1[j]];
630: xv[i] = (imode == INSERT_VALUES ? 0.0 : xv[i]) + sum;
631: }
632: PetscCall(VecRestoreArray(x, &xv));
633: PetscFunctionReturn(PETSC_SUCCESS);
634: }
636: PetscErrorCode VecDestroy_Seq(Vec v)
637: {
638: Vec_Seq *vs = (Vec_Seq *)v->data;
640: PetscFunctionBegin;
641: PetscCall(PetscLogObjectState((PetscObject)v, "Length=%" PetscInt_FMT, v->map->n));
642: if (vs) PetscCall(PetscShmgetDeallocateArray((void **)&vs->array_allocated));
643: PetscCall(VecResetPreallocationCOO_Seq(v));
644: PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscMatlabEnginePut_C", NULL));
645: PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscMatlabEngineGet_C", NULL));
646: PetscCall(PetscFree(v->data));
647: PetscFunctionReturn(PETSC_SUCCESS);
648: }
650: PetscErrorCode VecSetOption_Seq(Vec v, VecOption op, PetscBool flag)
651: {
652: PetscFunctionBegin;
653: if (op == VEC_IGNORE_NEGATIVE_INDICES) v->stash.ignorenegidx = flag;
654: PetscFunctionReturn(PETSC_SUCCESS);
655: }
657: // duplicate w to v. v is half-baked, potentially already with arrays allocated.
658: static PetscErrorCode VecDuplicate_Seq_Private(Vec w, Vec v)
659: {
660: PetscFunctionBegin;
661: PetscCall(VecSetType(v, ((PetscObject)w)->type_name));
662: PetscCall(PetscObjectListDuplicate(((PetscObject)w)->olist, &((PetscObject)v)->olist));
663: PetscCall(PetscFunctionListDuplicate(((PetscObject)w)->qlist, &((PetscObject)v)->qlist));
665: // Vec ops are not necessarily fully set by VecSetType(), e.g., see DMCreateGlobalVector_DA, so we copy w's to v
666: v->ops[0] = w->ops[0];
667: #if defined(PETSC_HAVE_DEVICE)
668: v->boundtocpu = w->boundtocpu;
669: v->bindingpropagates = w->bindingpropagates;
670: #endif
671: PetscFunctionReturn(PETSC_SUCCESS);
672: }
674: PetscErrorCode VecDuplicate_Seq(Vec win, Vec *V)
675: {
676: PetscFunctionBegin;
677: PetscCall(VecCreateWithLayout_Private(win->map, V));
678: PetscCall(VecDuplicate_Seq_Private(win, *V));
679: PetscFunctionReturn(PETSC_SUCCESS);
680: }
682: PetscErrorCode VecReplaceArray_Default_GEMV_Error(Vec v, const PetscScalar *a)
683: {
684: PetscFunctionBegin;
685: PetscCheck(PETSC_FALSE, PetscObjectComm((PetscObject)v), PETSC_ERR_SUP, "VecReplaceArray() is not supported on the first Vec obtained from VecDuplicateVecs(). \
686: You could either 1) use -vec_mdot_use_gemv 0 -vec_maxpy_use_gemv 0 to turn off an optimization to allow your current code to work or 2) use VecDuplicate() to duplicate the vector.");
687: (void)a;
688: PetscFunctionReturn(PETSC_SUCCESS);
689: }
691: static PetscErrorCode VecDuplicateVecs_Seq_GEMV(Vec w, PetscInt m, Vec *V[])
692: {
693: PetscScalar *array;
694: PetscInt64 lda; // use 64-bit as we will do "m * lda"
696: PetscFunctionBegin;
697: PetscCall(PetscMalloc1(m, V));
698: VecGetLocalSizeAligned(w, 64, &lda); // get in lda the 64-bytes aligned local size
699: PetscCall(PetscCalloc1(m * lda, &array));
700: for (PetscInt i = 0; i < m; i++) {
701: Vec v;
702: PetscCall(VecCreateSeqWithLayoutAndArray_Private(w->map, PetscSafePointerPlusOffset(array, i * lda), &v));
703: PetscCall(VecDuplicate_Seq_Private(w, v));
704: (*V)[i] = v;
705: }
706: // so when the first vector is destroyed it will destroy the array
707: if (m) ((Vec_Seq *)(*V)[0]->data)->array_allocated = array;
708: // disable replacearray of the first vector, as freeing its memory also frees others in the group.
709: // But replacearray of others is ok, as they don't own their array.
710: if (m > 1) (*V)[0]->ops->replacearray = VecReplaceArray_Default_GEMV_Error;
711: PetscFunctionReturn(PETSC_SUCCESS);
712: }
714: static struct _VecOps DvOps = {
715: PetscDesignatedInitializer(duplicate, VecDuplicate_Seq), /* 1 */
716: PetscDesignatedInitializer(duplicatevecs, VecDuplicateVecs_Default),
717: PetscDesignatedInitializer(destroyvecs, VecDestroyVecs_Default),
718: PetscDesignatedInitializer(dot, VecDot_Seq),
719: PetscDesignatedInitializer(mdot, VecMDot_Seq),
720: PetscDesignatedInitializer(norm, VecNorm_Seq),
721: PetscDesignatedInitializer(tdot, VecTDot_Seq),
722: PetscDesignatedInitializer(mtdot, VecMTDot_Seq),
723: PetscDesignatedInitializer(scale, VecScale_Seq),
724: PetscDesignatedInitializer(copy, VecCopy_Seq), /* 10 */
725: PetscDesignatedInitializer(set, VecSet_Seq),
726: PetscDesignatedInitializer(swap, VecSwap_Seq),
727: PetscDesignatedInitializer(axpy, VecAXPY_Seq),
728: PetscDesignatedInitializer(axpby, VecAXPBY_Seq),
729: PetscDesignatedInitializer(maxpy, VecMAXPY_Seq),
730: PetscDesignatedInitializer(aypx, VecAYPX_Seq),
731: PetscDesignatedInitializer(waxpy, VecWAXPY_Seq),
732: PetscDesignatedInitializer(axpbypcz, VecAXPBYPCZ_Seq),
733: PetscDesignatedInitializer(pointwisemult, VecPointwiseMult_Seq),
734: PetscDesignatedInitializer(pointwisedivide, VecPointwiseDivide_Seq),
735: PetscDesignatedInitializer(setvalues, VecSetValues_Seq), /* 20 */
736: PetscDesignatedInitializer(assemblybegin, NULL),
737: PetscDesignatedInitializer(assemblyend, NULL),
738: PetscDesignatedInitializer(getarray, NULL),
739: PetscDesignatedInitializer(getsize, VecGetSize_Seq),
740: PetscDesignatedInitializer(getlocalsize, VecGetSize_Seq),
741: PetscDesignatedInitializer(restorearray, NULL),
742: PetscDesignatedInitializer(max, VecMax_Seq),
743: PetscDesignatedInitializer(min, VecMin_Seq),
744: PetscDesignatedInitializer(setrandom, VecSetRandom_Seq),
745: PetscDesignatedInitializer(setoption, VecSetOption_Seq), /* 30 */
746: PetscDesignatedInitializer(setvaluesblocked, VecSetValuesBlocked_Seq),
747: PetscDesignatedInitializer(destroy, VecDestroy_Seq),
748: PetscDesignatedInitializer(view, VecView_Seq),
749: PetscDesignatedInitializer(placearray, VecPlaceArray_Seq),
750: PetscDesignatedInitializer(replacearray, VecReplaceArray_Seq),
751: PetscDesignatedInitializer(dot_local, VecDot_Seq),
752: PetscDesignatedInitializer(tdot_local, VecTDot_Seq),
753: PetscDesignatedInitializer(norm_local, VecNorm_Seq),
754: PetscDesignatedInitializer(mdot_local, VecMDot_Seq),
755: PetscDesignatedInitializer(mtdot_local, VecMTDot_Seq), /* 40 */
756: PetscDesignatedInitializer(load, VecLoad_Default),
757: PetscDesignatedInitializer(reciprocal, VecReciprocal_Default),
758: PetscDesignatedInitializer(conjugate, VecConjugate_Seq),
759: PetscDesignatedInitializer(setlocaltoglobalmapping, NULL),
760: PetscDesignatedInitializer(getlocaltoglobalmapping, NULL),
761: PetscDesignatedInitializer(resetarray, VecResetArray_Seq),
762: PetscDesignatedInitializer(setfromoptions, NULL),
763: PetscDesignatedInitializer(maxpointwisedivide, VecMaxPointwiseDivide_Seq),
764: PetscDesignatedInitializer(pointwisemax, VecPointwiseMax_Seq),
765: PetscDesignatedInitializer(pointwisemaxabs, VecPointwiseMaxAbs_Seq),
766: PetscDesignatedInitializer(pointwisemin, VecPointwiseMin_Seq),
767: PetscDesignatedInitializer(getvalues, VecGetValues_Seq),
768: PetscDesignatedInitializer(sqrt, NULL),
769: PetscDesignatedInitializer(abs, NULL),
770: PetscDesignatedInitializer(exp, NULL),
771: PetscDesignatedInitializer(log, NULL),
772: PetscDesignatedInitializer(shift, NULL),
773: PetscDesignatedInitializer(create, NULL),
774: PetscDesignatedInitializer(stridegather, VecStrideGather_Default),
775: PetscDesignatedInitializer(stridescatter, VecStrideScatter_Default),
776: PetscDesignatedInitializer(dotnorm2, NULL),
777: PetscDesignatedInitializer(getsubvector, NULL),
778: PetscDesignatedInitializer(restoresubvector, NULL),
779: PetscDesignatedInitializer(getarrayread, NULL),
780: PetscDesignatedInitializer(restorearrayread, NULL),
781: PetscDesignatedInitializer(stridesubsetgather, VecStrideSubSetGather_Default),
782: PetscDesignatedInitializer(stridesubsetscatter, VecStrideSubSetScatter_Default),
783: PetscDesignatedInitializer(viewnative, VecView_Seq),
784: PetscDesignatedInitializer(loadnative, NULL),
785: PetscDesignatedInitializer(createlocalvector, NULL),
786: PetscDesignatedInitializer(getlocalvector, NULL),
787: PetscDesignatedInitializer(restorelocalvector, NULL),
788: PetscDesignatedInitializer(getlocalvectorread, NULL),
789: PetscDesignatedInitializer(restorelocalvectorread, NULL),
790: PetscDesignatedInitializer(bindtocpu, NULL),
791: PetscDesignatedInitializer(getarraywrite, NULL),
792: PetscDesignatedInitializer(restorearraywrite, NULL),
793: PetscDesignatedInitializer(getarrayandmemtype, NULL),
794: PetscDesignatedInitializer(restorearrayandmemtype, NULL),
795: PetscDesignatedInitializer(getarrayreadandmemtype, NULL),
796: PetscDesignatedInitializer(restorearrayreadandmemtype, NULL),
797: PetscDesignatedInitializer(getarraywriteandmemtype, NULL),
798: PetscDesignatedInitializer(restorearraywriteandmemtype, NULL),
799: PetscDesignatedInitializer(concatenate, NULL),
800: PetscDesignatedInitializer(sum, NULL),
801: PetscDesignatedInitializer(setpreallocationcoo, VecSetPreallocationCOO_Seq),
802: PetscDesignatedInitializer(setvaluescoo, VecSetValuesCOO_Seq),
803: PetscDesignatedInitializer(errorwnorm, NULL),
804: PetscDesignatedInitializer(maxpby, NULL),
805: };
807: /*
808: Create a VECSEQ with the given layout and array
810: Input Parameter:
811: + map - the layout
812: - array - the array on host
814: Output Parameter:
815: . V - The vector object
816: */
817: PetscErrorCode VecCreateSeqWithLayoutAndArray_Private(PetscLayout map, const PetscScalar array[], Vec *V)
818: {
819: PetscMPIInt size;
821: PetscFunctionBegin;
822: PetscCall(VecCreateWithLayout_Private(map, V));
823: PetscCallMPI(MPI_Comm_size(map->comm, &size));
824: PetscCheck(size == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot create VECSEQ on more than one process");
825: PetscCall(VecCreate_Seq_Private(*V, array));
826: PetscFunctionReturn(PETSC_SUCCESS);
827: }
829: /*
830: This is called by VecCreate_Seq() (i.e. VecCreateSeq()) and VecCreateSeqWithArray()
831: */
832: PetscErrorCode VecCreate_Seq_Private(Vec v, const PetscScalar array[])
833: {
834: Vec_Seq *s;
835: PetscBool mdot_use_gemv = PETSC_TRUE;
836: PetscBool maxpy_use_gemv = PETSC_FALSE; // default is false as we saw bad performance with vendors' GEMV with tall skinny matrices.
838: PetscFunctionBegin;
839: PetscCall(PetscNew(&s));
840: v->ops[0] = DvOps;
842: PetscCall(PetscOptionsGetBool(NULL, NULL, "-vec_mdot_use_gemv", &mdot_use_gemv, NULL));
843: PetscCall(PetscOptionsGetBool(NULL, NULL, "-vec_maxpy_use_gemv", &maxpy_use_gemv, NULL));
845: // allocate multiple vectors together
846: if (mdot_use_gemv || maxpy_use_gemv) v->ops[0].duplicatevecs = VecDuplicateVecs_Seq_GEMV;
848: if (mdot_use_gemv) {
849: v->ops[0].mdot = VecMDot_Seq_GEMV;
850: v->ops[0].mdot_local = VecMDot_Seq_GEMV;
851: v->ops[0].mtdot = VecMTDot_Seq_GEMV;
852: v->ops[0].mtdot_local = VecMTDot_Seq_GEMV;
853: }
854: if (maxpy_use_gemv) v->ops[0].maxpy = VecMAXPY_Seq_GEMV;
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: /*@
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: }