Actual source code: bvec2.c
1: /*
2: Implements the sequential vectors.
3: */
5: #include <../src/vec/vec/impls/dvecimpl.h>
6: #include <../src/vec/vec/impls/mpi/pvecimpl.h>
7: #include <petsc/private/glvisviewerimpl.h>
8: #include <petsc/private/glvisvecimpl.h>
9: #include <petscblaslapack.h>
11: static PetscErrorCode VecPointwiseApply_Seq(Vec win, Vec xin, Vec yin, PetscScalar (*const func)(PetscScalar, PetscScalar))
12: {
13: const PetscInt n = win->map->n;
14: PetscScalar *ww, *xx, *yy; /* cannot make xx or yy const since might be ww */
16: PetscFunctionBegin;
17: PetscCall(VecGetArrayRead(xin, (const PetscScalar **)&xx));
18: PetscCall(VecGetArrayRead(yin, (const PetscScalar **)&yy));
19: PetscCall(VecGetArray(win, &ww));
20: for (PetscInt i = 0; i < n; ++i) ww[i] = func(xx[i], yy[i]);
21: PetscCall(VecRestoreArrayRead(xin, (const PetscScalar **)&xx));
22: PetscCall(VecRestoreArrayRead(yin, (const PetscScalar **)&yy));
23: PetscCall(VecRestoreArray(win, &ww));
24: PetscCall(PetscLogFlops(n));
25: PetscFunctionReturn(PETSC_SUCCESS);
26: }
28: static PetscScalar MaxRealPart(PetscScalar x, PetscScalar y)
29: {
30: // use temporaries to avoid reevaluating side-effects
31: const PetscReal rx = PetscRealPart(x), ry = PetscRealPart(y);
33: return PetscMax(rx, ry);
34: }
36: PetscErrorCode VecPointwiseMax_Seq(Vec win, Vec xin, Vec yin)
37: {
38: PetscFunctionBegin;
39: PetscCall(VecPointwiseApply_Seq(win, xin, yin, MaxRealPart));
40: PetscFunctionReturn(PETSC_SUCCESS);
41: }
43: static PetscScalar MinRealPart(PetscScalar x, PetscScalar y)
44: {
45: // use temporaries to avoid reevaluating side-effects
46: const PetscReal rx = PetscRealPart(x), ry = PetscRealPart(y);
48: return PetscMin(rx, ry);
49: }
51: PetscErrorCode VecPointwiseMin_Seq(Vec win, Vec xin, Vec yin)
52: {
53: PetscFunctionBegin;
54: PetscCall(VecPointwiseApply_Seq(win, xin, yin, MinRealPart));
55: PetscFunctionReturn(PETSC_SUCCESS);
56: }
58: static PetscScalar MaxAbs(PetscScalar x, PetscScalar y)
59: {
60: return (PetscScalar)PetscMax(PetscAbsScalar(x), PetscAbsScalar(y));
61: }
63: PetscErrorCode VecPointwiseMaxAbs_Seq(Vec win, Vec xin, Vec yin)
64: {
65: PetscFunctionBegin;
66: PetscCall(VecPointwiseApply_Seq(win, xin, yin, MaxAbs));
67: PetscFunctionReturn(PETSC_SUCCESS);
68: }
70: #include <../src/vec/vec/impls/seq/ftn-kernels/fxtimesy.h>
72: PetscErrorCode VecPointwiseMult_Seq(Vec win, Vec xin, Vec yin)
73: {
74: PetscInt n = win->map->n, i;
75: PetscScalar *ww, *xx, *yy; /* cannot make xx or yy const since might be ww */
77: PetscFunctionBegin;
78: PetscCall(VecGetArrayRead(xin, (const PetscScalar **)&xx));
79: PetscCall(VecGetArrayRead(yin, (const PetscScalar **)&yy));
80: PetscCall(VecGetArray(win, &ww));
81: if (ww == xx) {
82: for (i = 0; i < n; i++) ww[i] *= yy[i];
83: } else if (ww == yy) {
84: for (i = 0; i < n; i++) ww[i] *= xx[i];
85: } else {
86: #if defined(PETSC_USE_FORTRAN_KERNEL_XTIMESY)
87: fortranxtimesy_(xx, yy, ww, &n);
88: #else
89: for (i = 0; i < n; i++) ww[i] = xx[i] * yy[i];
90: #endif
91: }
92: PetscCall(VecRestoreArrayRead(xin, (const PetscScalar **)&xx));
93: PetscCall(VecRestoreArrayRead(yin, (const PetscScalar **)&yy));
94: PetscCall(VecRestoreArray(win, &ww));
95: PetscCall(PetscLogFlops(n));
96: PetscFunctionReturn(PETSC_SUCCESS);
97: }
99: static PetscScalar ScalDiv(PetscScalar x, PetscScalar y)
100: {
101: return y == 0.0 ? 0.0 : x / y;
102: }
104: PetscErrorCode VecPointwiseDivide_Seq(Vec win, Vec xin, Vec yin)
105: {
106: PetscFunctionBegin;
107: PetscCall(VecPointwiseApply_Seq(win, xin, yin, ScalDiv));
108: PetscFunctionReturn(PETSC_SUCCESS);
109: }
111: PetscErrorCode VecSetRandom_Seq(Vec xin, PetscRandom r)
112: {
113: PetscScalar *xx;
115: PetscFunctionBegin;
116: PetscCall(VecGetArrayWrite(xin, &xx));
117: PetscCall(PetscRandomGetValues(r, xin->map->n, xx));
118: PetscCall(VecRestoreArrayWrite(xin, &xx));
119: PetscFunctionReturn(PETSC_SUCCESS);
120: }
122: PetscErrorCode VecGetSize_Seq(Vec vin, PetscInt *size)
123: {
124: PetscFunctionBegin;
125: *size = vin->map->n;
126: PetscFunctionReturn(PETSC_SUCCESS);
127: }
129: PetscErrorCode VecConjugate_Seq(Vec xin)
130: {
131: PetscFunctionBegin;
132: if (PetscDefined(USE_COMPLEX)) {
133: const PetscInt n = xin->map->n;
134: PetscScalar *x;
136: PetscCall(VecGetArray(xin, &x));
137: for (PetscInt i = 0; i < n; ++i) x[i] = PetscConj(x[i]);
138: PetscCall(VecRestoreArray(xin, &x));
139: }
140: PetscFunctionReturn(PETSC_SUCCESS);
141: }
143: PetscErrorCode VecResetArray_Seq(Vec vin)
144: {
145: Vec_Seq *v = (Vec_Seq *)vin->data;
147: PetscFunctionBegin;
148: v->array = v->unplacedarray;
149: v->unplacedarray = NULL;
150: PetscFunctionReturn(PETSC_SUCCESS);
151: }
153: PetscErrorCode VecCopy_Seq(Vec xin, Vec yin)
154: {
155: PetscFunctionBegin;
156: if (xin != yin) {
157: const PetscScalar *xa;
158: PetscScalar *ya;
160: PetscCall(VecGetArrayRead(xin, &xa));
161: PetscCall(VecGetArray(yin, &ya));
162: PetscCall(PetscArraycpy(ya, xa, xin->map->n));
163: PetscCall(VecRestoreArrayRead(xin, &xa));
164: PetscCall(VecRestoreArray(yin, &ya));
165: }
166: PetscFunctionReturn(PETSC_SUCCESS);
167: }
169: PetscErrorCode VecSwap_Seq(Vec xin, Vec yin)
170: {
171: PetscFunctionBegin;
172: if (xin != yin) {
173: const PetscBLASInt one = 1;
174: PetscScalar *ya, *xa;
175: PetscBLASInt bn;
177: PetscCall(PetscBLASIntCast(xin->map->n, &bn));
178: PetscCall(VecGetArray(xin, &xa));
179: PetscCall(VecGetArray(yin, &ya));
180: PetscCallBLAS("BLASswap", BLASswap_(&bn, xa, &one, ya, &one));
181: PetscCall(VecRestoreArray(xin, &xa));
182: PetscCall(VecRestoreArray(yin, &ya));
183: }
184: PetscFunctionReturn(PETSC_SUCCESS);
185: }
187: #include <../src/vec/vec/impls/seq/ftn-kernels/fnorm.h>
189: PetscErrorCode VecNorm_Seq(Vec xin, NormType type, PetscReal *z)
190: {
191: // use a local variable to ensure compiler doesn't think z aliases any of the other arrays
192: PetscReal ztmp[] = {0.0, 0.0};
193: const PetscInt n = xin->map->n;
195: PetscFunctionBegin;
196: if (n) {
197: const PetscScalar *xx;
198: const PetscBLASInt one = 1;
199: PetscBLASInt bn = 0;
201: PetscCall(PetscBLASIntCast(n, &bn));
202: PetscCall(VecGetArrayRead(xin, &xx));
203: if (type == NORM_2 || type == NORM_FROBENIUS) {
204: NORM_1_AND_2_DOING_NORM_2:
205: if (PetscDefined(USE_REAL___FP16)) {
206: PetscCallBLAS("BLASnrm2", ztmp[type == NORM_1_AND_2] = BLASnrm2_(&bn, xx, &one));
207: } else {
208: PetscCallBLAS("BLASdot", ztmp[type == NORM_1_AND_2] = PetscSqrtReal(PetscRealPart(BLASdot_(&bn, xx, &one, xx, &one))));
209: }
210: PetscCall(PetscLogFlops(2.0 * n - 1));
211: } else if (type == NORM_INFINITY) {
212: for (PetscInt i = 0; i < n; ++i) {
213: const PetscReal tmp = PetscAbsScalar(xx[i]);
215: /* check special case of tmp == NaN */
216: if ((tmp > ztmp[0]) || (tmp != tmp)) {
217: ztmp[0] = tmp;
218: if (tmp != tmp) break;
219: }
220: }
221: } else if (type == NORM_1 || type == NORM_1_AND_2) {
222: if (PetscDefined(USE_COMPLEX)) {
223: // BLASasum() returns the nonstandard 1 norm of the 1 norm of the complex entries so we
224: // provide a custom loop instead
225: for (PetscInt i = 0; i < n; ++i) ztmp[0] += PetscAbsScalar(xx[i]);
226: } else {
227: PetscCallBLAS("BLASasum", ztmp[0] = BLASasum_(&bn, xx, &one));
228: }
229: PetscCall(PetscLogFlops(n - 1.0));
230: /* slight reshuffle so we can skip getting the array again (but still log the flops) if we
231: do norm2 after this */
232: if (type == NORM_1_AND_2) goto NORM_1_AND_2_DOING_NORM_2;
233: }
234: PetscCall(VecRestoreArrayRead(xin, &xx));
235: }
236: z[0] = ztmp[0];
237: if (type == NORM_1_AND_2) z[1] = ztmp[1];
238: PetscFunctionReturn(PETSC_SUCCESS);
239: }
241: static PetscErrorCode VecView_Seq_ASCII(Vec xin, PetscViewer viewer)
242: {
243: PetscInt i, n = xin->map->n;
244: const char *name;
245: PetscViewerFormat format;
246: const PetscScalar *xv;
248: PetscFunctionBegin;
249: PetscCall(VecGetArrayRead(xin, &xv));
250: PetscCall(PetscViewerGetFormat(viewer, &format));
251: if (format == PETSC_VIEWER_ASCII_MATLAB) {
252: PetscCall(PetscObjectGetName((PetscObject)xin, &name));
253: PetscCall(PetscViewerASCIIPrintf(viewer, "%s = [\n", name));
254: for (i = 0; i < n; i++) {
255: #if defined(PETSC_USE_COMPLEX)
256: if (PetscImaginaryPart(xv[i]) > 0.0) {
257: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e + %18.16ei\n", (double)PetscRealPart(xv[i]), (double)PetscImaginaryPart(xv[i])));
258: } else if (PetscImaginaryPart(xv[i]) < 0.0) {
259: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e - %18.16ei\n", (double)PetscRealPart(xv[i]), -(double)PetscImaginaryPart(xv[i])));
260: } else {
261: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e\n", (double)PetscRealPart(xv[i])));
262: }
263: #else
264: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e\n", (double)xv[i]));
265: #endif
266: }
267: PetscCall(PetscViewerASCIIPrintf(viewer, "];\n"));
268: } else if (format == PETSC_VIEWER_ASCII_SYMMODU) {
269: for (i = 0; i < n; i++) {
270: #if defined(PETSC_USE_COMPLEX)
271: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e %18.16e\n", (double)PetscRealPart(xv[i]), (double)PetscImaginaryPart(xv[i])));
272: #else
273: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e\n", (double)xv[i]));
274: #endif
275: }
276: } else if (format == PETSC_VIEWER_ASCII_VTK_DEPRECATED || format == PETSC_VIEWER_ASCII_VTK_CELL_DEPRECATED) {
277: /*
278: state 0: No header has been output
279: state 1: Only POINT_DATA has been output
280: state 2: Only CELL_DATA has been output
281: state 3: Output both, POINT_DATA last
282: state 4: Output both, CELL_DATA last
283: */
284: static PetscInt stateId = -1;
285: int outputState = 0;
286: PetscBool hasState;
287: int doOutput = 0;
288: PetscInt bs, b;
290: if (stateId < 0) PetscCall(PetscObjectComposedDataRegister(&stateId));
291: PetscCall(PetscObjectComposedDataGetInt((PetscObject)viewer, stateId, outputState, hasState));
292: if (!hasState) outputState = 0;
293: PetscCall(PetscObjectGetName((PetscObject)xin, &name));
294: PetscCall(VecGetBlockSize(xin, &bs));
295: PetscCheck(bs >= 1 && bs <= 3, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "VTK can only handle 3D objects, but vector dimension is %" PetscInt_FMT, bs);
296: if (format == PETSC_VIEWER_ASCII_VTK_DEPRECATED) {
297: if (outputState == 0) {
298: outputState = 1;
299: doOutput = 1;
300: } else if (outputState == 1) doOutput = 0;
301: else if (outputState == 2) {
302: outputState = 3;
303: doOutput = 1;
304: } else if (outputState == 3) doOutput = 0;
305: else PetscCheck(outputState != 4, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Tried to output POINT_DATA again after intervening CELL_DATA");
307: if (doOutput) PetscCall(PetscViewerASCIIPrintf(viewer, "POINT_DATA %" PetscInt_FMT "\n", n / bs));
308: } else {
309: if (outputState == 0) {
310: outputState = 2;
311: doOutput = 1;
312: } else if (outputState == 1) {
313: outputState = 4;
314: doOutput = 1;
315: } else if (outputState == 2) {
316: doOutput = 0;
317: } else {
318: PetscCheck(outputState != 3, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Tried to output CELL_DATA again after intervening POINT_DATA");
319: if (outputState == 4) doOutput = 0;
320: }
322: if (doOutput) PetscCall(PetscViewerASCIIPrintf(viewer, "CELL_DATA %" PetscInt_FMT "\n", n));
323: }
324: PetscCall(PetscObjectComposedDataSetInt((PetscObject)viewer, stateId, outputState));
325: if (name) {
326: if (bs == 3) {
327: PetscCall(PetscViewerASCIIPrintf(viewer, "VECTORS %s double\n", name));
328: } else {
329: PetscCall(PetscViewerASCIIPrintf(viewer, "SCALARS %s double %" PetscInt_FMT "\n", name, bs));
330: }
331: } else {
332: PetscCall(PetscViewerASCIIPrintf(viewer, "SCALARS scalars double %" PetscInt_FMT "\n", bs));
333: }
334: if (bs != 3) PetscCall(PetscViewerASCIIPrintf(viewer, "LOOKUP_TABLE default\n"));
335: for (i = 0; i < n / bs; i++) {
336: for (b = 0; b < bs; b++) {
337: if (b > 0) PetscCall(PetscViewerASCIIPrintf(viewer, " "));
338: #if !defined(PETSC_USE_COMPLEX)
339: PetscCall(PetscViewerASCIIPrintf(viewer, "%g", (double)xv[i * bs + b]));
340: #endif
341: }
342: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
343: }
344: } else if (format == PETSC_VIEWER_ASCII_VTK_COORDS_DEPRECATED) {
345: PetscInt bs, b;
347: PetscCall(VecGetBlockSize(xin, &bs));
348: PetscCheck(bs >= 1 && bs <= 3, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "VTK can only handle 3D objects, but vector dimension is %" PetscInt_FMT, bs);
349: for (i = 0; i < n / bs; i++) {
350: for (b = 0; b < bs; b++) {
351: if (b > 0) PetscCall(PetscViewerASCIIPrintf(viewer, " "));
352: #if !defined(PETSC_USE_COMPLEX)
353: PetscCall(PetscViewerASCIIPrintf(viewer, "%g", (double)xv[i * bs + b]));
354: #endif
355: }
356: for (b = bs; b < 3; b++) PetscCall(PetscViewerASCIIPrintf(viewer, " 0.0"));
357: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
358: }
359: } else if (format == PETSC_VIEWER_ASCII_PCICE) {
360: PetscInt bs, b;
362: PetscCall(VecGetBlockSize(xin, &bs));
363: PetscCheck(bs >= 1 && bs <= 3, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "PCICE can only handle up to 3D objects, but vector dimension is %" PetscInt_FMT, bs);
364: PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT "\n", xin->map->N / bs));
365: for (i = 0; i < n / bs; i++) {
366: PetscCall(PetscViewerASCIIPrintf(viewer, "%7" PetscInt_FMT " ", i + 1));
367: for (b = 0; b < bs; b++) {
368: if (b > 0) PetscCall(PetscViewerASCIIPrintf(viewer, " "));
369: #if !defined(PETSC_USE_COMPLEX)
370: PetscCall(PetscViewerASCIIPrintf(viewer, "% 12.5E", (double)xv[i * bs + b]));
371: #endif
372: }
373: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
374: }
375: } else if (format == PETSC_VIEWER_ASCII_GLVIS) {
376: /* GLVis ASCII visualization/dump: this function mimics mfem::GridFunction::Save() */
377: const PetscScalar *array;
378: PetscInt i, n, vdim, ordering = 1; /* mfem::FiniteElementSpace::Ordering::byVDIM */
379: PetscContainer glvis_container;
380: PetscViewerGLVisVecInfo glvis_vec_info;
381: PetscViewerGLVisInfo glvis_info;
383: /* mfem::FiniteElementSpace::Save() */
384: PetscCall(VecGetBlockSize(xin, &vdim));
385: PetscCall(PetscViewerASCIIPrintf(viewer, "FiniteElementSpace\n"));
386: PetscCall(PetscObjectQuery((PetscObject)xin, "_glvis_info_container", (PetscObject *)&glvis_container));
387: PetscCheck(glvis_container, PetscObjectComm((PetscObject)xin), PETSC_ERR_PLIB, "Missing GLVis container");
388: PetscCall(PetscContainerGetPointer(glvis_container, (void **)&glvis_vec_info));
389: PetscCall(PetscViewerASCIIPrintf(viewer, "%s\n", glvis_vec_info->fec_type));
390: PetscCall(PetscViewerASCIIPrintf(viewer, "VDim: %" PetscInt_FMT "\n", vdim));
391: PetscCall(PetscViewerASCIIPrintf(viewer, "Ordering: %" PetscInt_FMT "\n", ordering));
392: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
393: /* mfem::Vector::Print() */
394: PetscCall(PetscObjectQuery((PetscObject)viewer, "_glvis_info_container", (PetscObject *)&glvis_container));
395: PetscCheck(glvis_container, PetscObjectComm((PetscObject)viewer), PETSC_ERR_PLIB, "Missing GLVis container");
396: PetscCall(PetscContainerGetPointer(glvis_container, (void **)&glvis_info));
397: if (glvis_info->enabled) {
398: PetscCall(VecGetLocalSize(xin, &n));
399: PetscCall(VecGetArrayRead(xin, &array));
400: for (i = 0; i < n; i++) {
401: PetscCall(PetscViewerASCIIPrintf(viewer, glvis_info->fmt, (double)PetscRealPart(array[i])));
402: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
403: }
404: PetscCall(VecRestoreArrayRead(xin, &array));
405: }
406: } else if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
407: /* No info */
408: } else {
409: for (i = 0; i < n; i++) {
410: if (format == PETSC_VIEWER_ASCII_INDEX) PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT ": ", i));
411: #if defined(PETSC_USE_COMPLEX)
412: if (PetscImaginaryPart(xv[i]) > 0.0) {
413: PetscCall(PetscViewerASCIIPrintf(viewer, "%g + %g i\n", (double)PetscRealPart(xv[i]), (double)PetscImaginaryPart(xv[i])));
414: } else if (PetscImaginaryPart(xv[i]) < 0.0) {
415: PetscCall(PetscViewerASCIIPrintf(viewer, "%g - %g i\n", (double)PetscRealPart(xv[i]), -(double)PetscImaginaryPart(xv[i])));
416: } else {
417: PetscCall(PetscViewerASCIIPrintf(viewer, "%g\n", (double)PetscRealPart(xv[i])));
418: }
419: #else
420: PetscCall(PetscViewerASCIIPrintf(viewer, "%g\n", (double)xv[i]));
421: #endif
422: }
423: }
424: PetscCall(PetscViewerFlush(viewer));
425: PetscCall(VecRestoreArrayRead(xin, &xv));
426: PetscFunctionReturn(PETSC_SUCCESS);
427: }
429: #include <petscdraw.h>
430: static PetscErrorCode VecView_Seq_Draw_LG(Vec xin, PetscViewer v)
431: {
432: PetscDraw draw;
433: PetscBool isnull;
434: PetscDrawLG lg;
435: PetscInt i, c, bs = PetscAbs(xin->map->bs), n = xin->map->n / bs;
436: const PetscScalar *xv;
437: PetscReal *xx, *yy, xmin, xmax, h;
438: int colors[] = {PETSC_DRAW_RED};
439: PetscViewerFormat format;
440: PetscDrawAxis axis;
442: PetscFunctionBegin;
443: PetscCall(PetscViewerDrawGetDraw(v, 0, &draw));
444: PetscCall(PetscDrawIsNull(draw, &isnull));
445: if (isnull) PetscFunctionReturn(PETSC_SUCCESS);
447: PetscCall(PetscViewerGetFormat(v, &format));
448: PetscCall(PetscMalloc2(n, &xx, n, &yy));
449: PetscCall(VecGetArrayRead(xin, &xv));
450: for (c = 0; c < bs; c++) {
451: PetscCall(PetscViewerDrawGetDrawLG(v, c, &lg));
452: PetscCall(PetscDrawLGReset(lg));
453: PetscCall(PetscDrawLGSetDimension(lg, 1));
454: PetscCall(PetscDrawLGSetColors(lg, colors));
455: if (format == PETSC_VIEWER_DRAW_LG_XRANGE) {
456: PetscCall(PetscDrawLGGetAxis(lg, &axis));
457: PetscCall(PetscDrawAxisGetLimits(axis, &xmin, &xmax, NULL, NULL));
458: h = (xmax - xmin) / n;
459: for (i = 0; i < n; i++) xx[i] = i * h + 0.5 * h; /* cell center */
460: } else {
461: for (i = 0; i < n; i++) xx[i] = (PetscReal)i;
462: }
463: for (i = 0; i < n; i++) yy[i] = PetscRealPart(xv[c + i * bs]);
465: PetscCall(PetscDrawLGAddPoints(lg, n, &xx, &yy));
466: PetscCall(PetscDrawLGDraw(lg));
467: PetscCall(PetscDrawLGSave(lg));
468: }
469: PetscCall(VecRestoreArrayRead(xin, &xv));
470: PetscCall(PetscFree2(xx, yy));
471: PetscFunctionReturn(PETSC_SUCCESS);
472: }
474: static PetscErrorCode VecView_Seq_Draw(Vec xin, PetscViewer v)
475: {
476: PetscDraw draw;
477: PetscBool isnull;
479: PetscFunctionBegin;
480: PetscCall(PetscViewerDrawGetDraw(v, 0, &draw));
481: PetscCall(PetscDrawIsNull(draw, &isnull));
482: if (isnull) PetscFunctionReturn(PETSC_SUCCESS);
484: PetscCall(VecView_Seq_Draw_LG(xin, v));
485: PetscFunctionReturn(PETSC_SUCCESS);
486: }
488: static PetscErrorCode VecView_Seq_Binary(Vec xin, PetscViewer viewer)
489: {
490: return VecView_Binary(xin, viewer);
491: }
493: #if defined(PETSC_HAVE_MATLAB)
494: #include <petscmatlab.h>
495: #include <mat.h> /* MATLAB include file */
496: PetscErrorCode VecView_Seq_Matlab(Vec vec, PetscViewer viewer)
497: {
498: PetscInt n;
499: const PetscScalar *array;
501: PetscFunctionBegin;
502: PetscCall(VecGetLocalSize(vec, &n));
503: PetscCall(PetscObjectName((PetscObject)vec));
504: PetscCall(VecGetArrayRead(vec, &array));
505: PetscCall(PetscViewerMatlabPutArray(viewer, n, 1, array, ((PetscObject)vec)->name));
506: PetscCall(VecRestoreArrayRead(vec, &array));
507: PetscFunctionReturn(PETSC_SUCCESS);
508: }
509: #endif
511: PetscErrorCode VecView_Seq(Vec xin, PetscViewer viewer)
512: {
513: PetscBool isdraw, iascii, issocket, isbinary;
514: #if defined(PETSC_HAVE_MATHEMATICA)
515: PetscBool ismathematica;
516: #endif
517: #if defined(PETSC_HAVE_MATLAB)
518: PetscBool ismatlab;
519: #endif
520: #if defined(PETSC_HAVE_HDF5)
521: PetscBool ishdf5;
522: #endif
523: PetscBool isglvis;
524: #if defined(PETSC_HAVE_ADIOS)
525: PetscBool isadios;
526: #endif
528: PetscFunctionBegin;
529: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
530: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
531: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSOCKET, &issocket));
532: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
533: #if defined(PETSC_HAVE_MATHEMATICA)
534: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERMATHEMATICA, &ismathematica));
535: #endif
536: #if defined(PETSC_HAVE_HDF5)
537: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
538: #endif
539: #if defined(PETSC_HAVE_MATLAB)
540: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERMATLAB, &ismatlab));
541: #endif
542: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERGLVIS, &isglvis));
543: #if defined(PETSC_HAVE_ADIOS)
544: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERADIOS, &isadios));
545: #endif
547: if (isdraw) {
548: PetscCall(VecView_Seq_Draw(xin, viewer));
549: } else if (iascii) {
550: PetscCall(VecView_Seq_ASCII(xin, viewer));
551: } else if (isbinary) {
552: PetscCall(VecView_Seq_Binary(xin, viewer));
553: #if defined(PETSC_HAVE_MATHEMATICA)
554: } else if (ismathematica) {
555: PetscCall(PetscViewerMathematicaPutVector(viewer, xin));
556: #endif
557: #if defined(PETSC_HAVE_HDF5)
558: } else if (ishdf5) {
559: PetscCall(VecView_MPI_HDF5(xin, viewer)); /* Reusing VecView_MPI_HDF5 ... don't want code duplication*/
560: #endif
561: #if defined(PETSC_HAVE_ADIOS)
562: } else if (isadios) {
563: PetscCall(VecView_MPI_ADIOS(xin, viewer)); /* Reusing VecView_MPI_ADIOS ... don't want code duplication*/
564: #endif
565: #if defined(PETSC_HAVE_MATLAB)
566: } else if (ismatlab) {
567: PetscCall(VecView_Seq_Matlab(xin, viewer));
568: #endif
569: } else if (isglvis) PetscCall(VecView_GLVis(xin, viewer));
570: PetscFunctionReturn(PETSC_SUCCESS);
571: }
573: PetscErrorCode VecGetValues_Seq(Vec xin, PetscInt ni, const PetscInt ix[], PetscScalar y[])
574: {
575: const PetscBool ignorenegidx = xin->stash.ignorenegidx;
576: const PetscScalar *xx;
578: PetscFunctionBegin;
579: PetscCall(VecGetArrayRead(xin, &xx));
580: for (PetscInt i = 0; i < ni; ++i) {
581: if (ignorenegidx && (ix[i] < 0)) continue;
582: if (PetscDefined(USE_DEBUG)) {
583: PetscCheck(ix[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Out of range index value %" PetscInt_FMT " cannot be negative", ix[i]);
584: PetscCheck(ix[i] < xin->map->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Out of range index value %" PetscInt_FMT " to large maximum allowed %" PetscInt_FMT, ix[i], xin->map->n);
585: }
586: y[i] = xx[ix[i]];
587: }
588: PetscCall(VecRestoreArrayRead(xin, &xx));
589: PetscFunctionReturn(PETSC_SUCCESS);
590: }
592: PetscErrorCode VecSetValues_Seq(Vec xin, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode m)
593: {
594: const PetscBool ignorenegidx = xin->stash.ignorenegidx;
595: PetscScalar *xx;
597: PetscFunctionBegin;
598: // call to getarray (not e.g. getarraywrite() if m is INSERT_VALUES) is deliberate! If this
599: // is secretly a VECSEQCUDA it may have values currently on the device, in which case --
600: // unless we are replacing the entire array -- we need to copy them up
601: PetscCall(VecGetArray(xin, &xx));
602: for (PetscInt i = 0; i < ni; i++) {
603: if (ignorenegidx && (ix[i] < 0)) continue;
604: if (PetscDefined(USE_DEBUG)) {
605: PetscCheck(ix[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Out of range index value %" PetscInt_FMT " cannot be negative", ix[i]);
606: PetscCheck(ix[i] < xin->map->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Out of range index value %" PetscInt_FMT " maximum %" PetscInt_FMT, ix[i], xin->map->n);
607: }
608: if (m == INSERT_VALUES) {
609: xx[ix[i]] = y[i];
610: } else {
611: xx[ix[i]] += y[i];
612: }
613: }
614: PetscCall(VecRestoreArray(xin, &xx));
615: PetscFunctionReturn(PETSC_SUCCESS);
616: }
618: PetscErrorCode VecSetValuesBlocked_Seq(Vec xin, PetscInt ni, const PetscInt ix[], const PetscScalar yin[], InsertMode m)
619: {
620: PetscScalar *xx;
621: PetscInt bs;
623: /* For optimization could treat bs = 2, 3, 4, 5 as special cases with loop unrolling */
624: PetscFunctionBegin;
625: PetscCall(VecGetBlockSize(xin, &bs));
626: PetscCall(VecGetArray(xin, &xx));
627: for (PetscInt i = 0; i < ni; ++i, yin += bs) {
628: const PetscInt start = bs * ix[i];
630: if (start < 0) continue;
631: PetscCheck(start < xin->map->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Out of range index value %" PetscInt_FMT " maximum %" PetscInt_FMT, start, xin->map->n);
632: for (PetscInt j = 0; j < bs; j++) {
633: if (m == INSERT_VALUES) {
634: xx[start + j] = yin[j];
635: } else {
636: xx[start + j] += yin[j];
637: }
638: }
639: }
640: PetscCall(VecRestoreArray(xin, &xx));
641: PetscFunctionReturn(PETSC_SUCCESS);
642: }
644: static PetscErrorCode VecResetPreallocationCOO_Seq(Vec x)
645: {
646: Vec_Seq *vs = (Vec_Seq *)x->data;
648: PetscFunctionBegin;
649: if (vs) {
650: PetscCall(PetscFree(vs->jmap1)); /* Destroy old stuff */
651: PetscCall(PetscFree(vs->perm1));
652: }
653: PetscFunctionReturn(PETSC_SUCCESS);
654: }
656: PetscErrorCode VecSetPreallocationCOO_Seq(Vec x, PetscCount coo_n, const PetscInt coo_i[])
657: {
658: PetscInt m, *i;
659: PetscCount k, nneg;
660: PetscCount *perm1, *jmap1;
661: Vec_Seq *vs = (Vec_Seq *)x->data;
663: PetscFunctionBegin;
664: PetscCall(VecResetPreallocationCOO_Seq(x)); /* Destroy old stuff */
665: PetscCall(PetscMalloc1(coo_n, &i));
666: PetscCall(PetscArraycpy(i, coo_i, coo_n)); /* Make a copy since we'll modify it */
667: PetscCall(PetscMalloc1(coo_n, &perm1));
668: for (k = 0; k < coo_n; k++) perm1[k] = k;
669: PetscCall(PetscSortIntWithCountArray(coo_n, i, perm1));
670: for (k = 0; k < coo_n; k++) {
671: if (i[k] >= 0) break;
672: } /* Advance k to the first entry with a non-negative index */
673: nneg = k;
675: PetscCall(VecGetLocalSize(x, &m));
676: PetscCheck(!nneg || x->stash.ignorenegidx, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Found a negative index in VecSetPreallocateCOO() but VEC_IGNORE_NEGATIVE_INDICES was not set");
677: PetscCheck(!coo_n || i[coo_n - 1] < m, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Found index (%" PetscInt_FMT ") greater than the size of the vector (%" PetscInt_FMT ") in VecSetPreallocateCOO()", i[coo_n - 1], m);
679: PetscCall(PetscCalloc1(m + 1, &jmap1));
680: for (; k < coo_n; k++) jmap1[i[k] + 1]++; /* Count repeats of each entry */
681: for (k = 0; k < m; k++) jmap1[k + 1] += jmap1[k]; /* Transform jmap[] to CSR-like data structure */
682: PetscCall(PetscFree(i));
684: if (nneg) { /* Discard leading negative indices */
685: PetscCount *perm1_new;
686: PetscCall(PetscMalloc1(coo_n - nneg, &perm1_new));
687: PetscCall(PetscArraycpy(perm1_new, perm1 + nneg, coo_n - nneg));
688: PetscCall(PetscFree(perm1));
689: perm1 = perm1_new;
690: }
692: /* Record COO fields */
693: vs->coo_n = coo_n;
694: vs->tot1 = coo_n - nneg;
695: vs->jmap1 = jmap1; /* [m+1] */
696: vs->perm1 = perm1; /* [tot] */
697: PetscFunctionReturn(PETSC_SUCCESS);
698: }
700: PetscErrorCode VecSetValuesCOO_Seq(Vec x, const PetscScalar coo_v[], InsertMode imode)
701: {
702: Vec_Seq *vs = (Vec_Seq *)x->data;
703: const PetscCount *perm1 = vs->perm1, *jmap1 = vs->jmap1;
704: PetscScalar *xv;
705: PetscInt m;
707: PetscFunctionBegin;
708: PetscCall(VecGetLocalSize(x, &m));
709: PetscCall(VecGetArray(x, &xv));
710: for (PetscInt i = 0; i < m; i++) {
711: PetscScalar sum = 0.0;
712: for (PetscCount j = jmap1[i]; j < jmap1[i + 1]; j++) sum += coo_v[perm1[j]];
713: xv[i] = (imode == INSERT_VALUES ? 0.0 : xv[i]) + sum;
714: }
715: PetscCall(VecRestoreArray(x, &xv));
716: PetscFunctionReturn(PETSC_SUCCESS);
717: }
719: PetscErrorCode VecDestroy_Seq(Vec v)
720: {
721: Vec_Seq *vs = (Vec_Seq *)v->data;
723: PetscFunctionBegin;
724: PetscCall(PetscLogObjectState((PetscObject)v, "Length=%" PetscInt_FMT, v->map->n));
725: if (vs) PetscCall(PetscFree(vs->array_allocated));
726: PetscCall(VecResetPreallocationCOO_Seq(v));
727: PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscMatlabEnginePut_C", NULL));
728: PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscMatlabEngineGet_C", NULL));
729: PetscCall(PetscFree(v->data));
730: PetscFunctionReturn(PETSC_SUCCESS);
731: }
733: PetscErrorCode VecSetOption_Seq(Vec v, VecOption op, PetscBool flag)
734: {
735: PetscFunctionBegin;
736: if (op == VEC_IGNORE_NEGATIVE_INDICES) v->stash.ignorenegidx = flag;
737: PetscFunctionReturn(PETSC_SUCCESS);
738: }
740: PetscErrorCode VecDuplicate_Seq(Vec win, Vec *V)
741: {
742: PetscFunctionBegin;
743: PetscCall(VecCreateWithLayout_Private(win->map, V));
744: PetscCall(VecSetType(*V, ((PetscObject)win)->type_name));
745: PetscCall(PetscObjectListDuplicate(((PetscObject)win)->olist, &((PetscObject)(*V))->olist));
746: PetscCall(PetscFunctionListDuplicate(((PetscObject)win)->qlist, &((PetscObject)(*V))->qlist));
748: (*V)->ops->view = win->ops->view;
749: (*V)->stash.ignorenegidx = win->stash.ignorenegidx;
750: PetscFunctionReturn(PETSC_SUCCESS);
751: }
753: static const struct _VecOps DvOps = {
754: PetscDesignatedInitializer(duplicate, VecDuplicate_Seq), /* 1 */
755: PetscDesignatedInitializer(duplicatevecs, VecDuplicateVecs_Default),
756: PetscDesignatedInitializer(destroyvecs, VecDestroyVecs_Default),
757: PetscDesignatedInitializer(dot, VecDot_Seq),
758: PetscDesignatedInitializer(mdot, VecMDot_Seq),
759: PetscDesignatedInitializer(norm, VecNorm_Seq),
760: PetscDesignatedInitializer(tdot, VecTDot_Seq),
761: PetscDesignatedInitializer(mtdot, VecMTDot_Seq),
762: PetscDesignatedInitializer(scale, VecScale_Seq),
763: PetscDesignatedInitializer(copy, VecCopy_Seq), /* 10 */
764: PetscDesignatedInitializer(set, VecSet_Seq),
765: PetscDesignatedInitializer(swap, VecSwap_Seq),
766: PetscDesignatedInitializer(axpy, VecAXPY_Seq),
767: PetscDesignatedInitializer(axpby, VecAXPBY_Seq),
768: PetscDesignatedInitializer(maxpy, VecMAXPY_Seq),
769: PetscDesignatedInitializer(aypx, VecAYPX_Seq),
770: PetscDesignatedInitializer(waxpy, VecWAXPY_Seq),
771: PetscDesignatedInitializer(axpbypcz, VecAXPBYPCZ_Seq),
772: PetscDesignatedInitializer(pointwisemult, VecPointwiseMult_Seq),
773: PetscDesignatedInitializer(pointwisedivide, VecPointwiseDivide_Seq),
774: PetscDesignatedInitializer(setvalues, VecSetValues_Seq), /* 20 */
775: PetscDesignatedInitializer(assemblybegin, NULL),
776: PetscDesignatedInitializer(assemblyend, NULL),
777: PetscDesignatedInitializer(getarray, NULL),
778: PetscDesignatedInitializer(getsize, VecGetSize_Seq),
779: PetscDesignatedInitializer(getlocalsize, VecGetSize_Seq),
780: PetscDesignatedInitializer(restorearray, NULL),
781: PetscDesignatedInitializer(max, VecMax_Seq),
782: PetscDesignatedInitializer(min, VecMin_Seq),
783: PetscDesignatedInitializer(setrandom, VecSetRandom_Seq),
784: PetscDesignatedInitializer(setoption, VecSetOption_Seq), /* 30 */
785: PetscDesignatedInitializer(setvaluesblocked, VecSetValuesBlocked_Seq),
786: PetscDesignatedInitializer(destroy, VecDestroy_Seq),
787: PetscDesignatedInitializer(view, VecView_Seq),
788: PetscDesignatedInitializer(placearray, VecPlaceArray_Seq),
789: PetscDesignatedInitializer(replacearray, VecReplaceArray_Seq),
790: PetscDesignatedInitializer(dot_local, VecDot_Seq),
791: PetscDesignatedInitializer(tdot_local, VecTDot_Seq),
792: PetscDesignatedInitializer(norm_local, VecNorm_Seq),
793: PetscDesignatedInitializer(mdot_local, VecMDot_Seq),
794: PetscDesignatedInitializer(mtdot_local, VecMTDot_Seq), /* 40 */
795: PetscDesignatedInitializer(load, VecLoad_Default),
796: PetscDesignatedInitializer(reciprocal, VecReciprocal_Default),
797: PetscDesignatedInitializer(conjugate, VecConjugate_Seq),
798: PetscDesignatedInitializer(setlocaltoglobalmapping, NULL),
799: PetscDesignatedInitializer(setvalueslocal, NULL),
800: PetscDesignatedInitializer(resetarray, VecResetArray_Seq),
801: PetscDesignatedInitializer(setfromoptions, NULL),
802: PetscDesignatedInitializer(maxpointwisedivide, VecMaxPointwiseDivide_Seq),
803: PetscDesignatedInitializer(pointwisemax, VecPointwiseMax_Seq),
804: PetscDesignatedInitializer(pointwisemaxabs, VecPointwiseMaxAbs_Seq),
805: PetscDesignatedInitializer(pointwisemin, VecPointwiseMin_Seq),
806: PetscDesignatedInitializer(getvalues, VecGetValues_Seq),
807: PetscDesignatedInitializer(sqrt, NULL),
808: PetscDesignatedInitializer(abs, NULL),
809: PetscDesignatedInitializer(exp, NULL),
810: PetscDesignatedInitializer(log, NULL),
811: PetscDesignatedInitializer(shift, NULL),
812: PetscDesignatedInitializer(create, NULL),
813: PetscDesignatedInitializer(stridegather, VecStrideGather_Default),
814: PetscDesignatedInitializer(stridescatter, VecStrideScatter_Default),
815: PetscDesignatedInitializer(dotnorm2, NULL),
816: PetscDesignatedInitializer(getsubvector, NULL),
817: PetscDesignatedInitializer(restoresubvector, NULL),
818: PetscDesignatedInitializer(getarrayread, NULL),
819: PetscDesignatedInitializer(restorearrayread, NULL),
820: PetscDesignatedInitializer(stridesubsetgather, VecStrideSubSetGather_Default),
821: PetscDesignatedInitializer(stridesubsetscatter, VecStrideSubSetScatter_Default),
822: PetscDesignatedInitializer(viewnative, VecView_Seq),
823: PetscDesignatedInitializer(loadnative, NULL),
824: PetscDesignatedInitializer(createlocalvector, NULL),
825: PetscDesignatedInitializer(getlocalvector, NULL),
826: PetscDesignatedInitializer(restorelocalvector, NULL),
827: PetscDesignatedInitializer(getlocalvectorread, NULL),
828: PetscDesignatedInitializer(restorelocalvectorread, NULL),
829: PetscDesignatedInitializer(bindtocpu, NULL),
830: PetscDesignatedInitializer(getarraywrite, NULL),
831: PetscDesignatedInitializer(restorearraywrite, NULL),
832: PetscDesignatedInitializer(getarrayandmemtype, NULL),
833: PetscDesignatedInitializer(restorearrayandmemtype, NULL),
834: PetscDesignatedInitializer(getarrayreadandmemtype, NULL),
835: PetscDesignatedInitializer(restorearrayreadandmemtype, NULL),
836: PetscDesignatedInitializer(getarraywriteandmemtype, NULL),
837: PetscDesignatedInitializer(restorearraywriteandmemtype, NULL),
838: PetscDesignatedInitializer(concatenate, NULL),
839: PetscDesignatedInitializer(sum, NULL),
840: PetscDesignatedInitializer(setpreallocationcoo, VecSetPreallocationCOO_Seq),
841: PetscDesignatedInitializer(setvaluescoo, VecSetValuesCOO_Seq),
842: PetscDesignatedInitializer(errorwnorm, NULL),
843: };
845: /*
846: This is called by VecCreate_Seq() (i.e. VecCreateSeq()) and VecCreateSeqWithArray()
847: */
848: PetscErrorCode VecCreate_Seq_Private(Vec v, const PetscScalar array[])
849: {
850: Vec_Seq *s;
852: PetscFunctionBegin;
853: PetscCall(PetscNew(&s));
854: v->ops[0] = DvOps;
856: v->data = (void *)s;
857: v->petscnative = PETSC_TRUE;
858: s->array = (PetscScalar *)array;
859: s->array_allocated = NULL;
860: if (array) v->offloadmask = PETSC_OFFLOAD_CPU;
862: PetscCall(PetscLayoutSetUp(v->map));
863: PetscCall(PetscObjectChangeTypeName((PetscObject)v, VECSEQ));
864: #if defined(PETSC_HAVE_MATLAB)
865: PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscMatlabEnginePut_C", VecMatlabEnginePut_Default));
866: PetscCall(PetscObjectComposeFunction((PetscObject)v, "PetscMatlabEngineGet_C", VecMatlabEngineGet_Default));
867: #endif
868: PetscFunctionReturn(PETSC_SUCCESS);
869: }
871: /*@C
872: VecCreateSeqWithArray - Creates a standard,sequential array-style vector,
873: where the user provides the array space to store the vector values.
875: Collective
877: Input Parameters:
878: + comm - the communicator, should be `PETSC_COMM_SELF`
879: . bs - the block size
880: . n - the vector length
881: - array - memory where the vector elements are to be stored.
883: Output Parameter:
884: . V - the vector
886: Level: intermediate
888: Notes:
889: Use `VecDuplicate()` or `VecDuplicateVecs(`) to form additional vectors of the
890: same type as an existing vector.
892: If the user-provided array is` NULL`, then `VecPlaceArray()` can be used
893: at a later stage to SET the array for storing the vector values.
895: PETSc does NOT free the array when the vector is destroyed via `VecDestroy()`.
896: The user should not free the array until the vector is destroyed.
898: .seealso: `VecCreateMPIWithArray()`, `VecCreate()`, `VecDuplicate()`, `VecDuplicateVecs()`,
899: `VecCreateGhost()`, `VecCreateSeq()`, `VecPlaceArray()`
900: @*/
901: PetscErrorCode VecCreateSeqWithArray(MPI_Comm comm, PetscInt bs, PetscInt n, const PetscScalar array[], Vec *V)
902: {
903: PetscMPIInt size;
905: PetscFunctionBegin;
906: PetscCall(VecCreate(comm, V));
907: PetscCall(VecSetSizes(*V, n, n));
908: PetscCall(VecSetBlockSize(*V, bs));
909: PetscCallMPI(MPI_Comm_size(comm, &size));
910: PetscCheck(size <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot create VECSEQ on more than one process");
911: PetscCall(VecCreate_Seq_Private(*V, array));
912: PetscFunctionReturn(PETSC_SUCCESS);
913: }