1 /*
2 Implements the sequential vectors.
3 */
4
5 #include <../src/vec/vec/impls/dvecimpl.h> /*I "petscvec.h" I*/
6 #include <../src/vec/vec/impls/mpi/pvecimpl.h> /* For VecView_MPI_HDF5 */
7 #include <petsc/private/glvisviewerimpl.h>
8 #include <petsc/private/glvisvecimpl.h>
9 #include <petscblaslapack.h>
10
VecPointwiseApply_Seq(Vec win,Vec xin,Vec yin,PetscScalar (* const func)(PetscScalar,PetscScalar))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 */
15
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 }
27
MaxRealPart(PetscScalar x,PetscScalar y)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);
32
33 return PetscMax(rx, ry);
34 }
35
VecPointwiseMax_Seq(Vec win,Vec xin,Vec yin)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 }
42
MinRealPart(PetscScalar x,PetscScalar y)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);
47
48 return PetscMin(rx, ry);
49 }
50
VecPointwiseMin_Seq(Vec win,Vec xin,Vec yin)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 }
57
MaxAbs(PetscScalar x,PetscScalar y)58 static PetscScalar MaxAbs(PetscScalar x, PetscScalar y)
59 {
60 return PetscMax(PetscAbsScalar(x), PetscAbsScalar(y));
61 }
62
VecPointwiseMaxAbs_Seq(Vec win,Vec xin,Vec yin)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 }
69
70 #include <../src/vec/vec/impls/seq/ftn-kernels/fxtimesy.h>
71
VecPointwiseMult_Seq(Vec win,Vec xin,Vec yin)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 */
76
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 }
98
ScalDiv(PetscScalar x,PetscScalar y)99 static PetscScalar ScalDiv(PetscScalar x, PetscScalar y)
100 {
101 return y == 0.0 ? 0.0 : x / y;
102 }
103
VecPointwiseDivide_Seq(Vec win,Vec xin,Vec yin)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 }
110
VecSetRandom_Seq(Vec xin,PetscRandom r)111 PetscErrorCode VecSetRandom_Seq(Vec xin, PetscRandom r)
112 {
113 PetscScalar *xx;
114
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 }
121
VecGetSize_Seq(Vec vin,PetscInt * size)122 PetscErrorCode VecGetSize_Seq(Vec vin, PetscInt *size)
123 {
124 PetscFunctionBegin;
125 *size = vin->map->n;
126 PetscFunctionReturn(PETSC_SUCCESS);
127 }
128
VecConjugate_Seq(Vec xin)129 PetscErrorCode VecConjugate_Seq(Vec xin)
130 {
131 const PetscInt n = xin->map->n;
132 PetscScalar *x;
133
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 }
140
VecResetArray_Seq(Vec vin)141 PetscErrorCode VecResetArray_Seq(Vec vin)
142 {
143 Vec_Seq *v = (Vec_Seq *)vin->data;
144
145 PetscFunctionBegin;
146 v->array = v->unplacedarray;
147 v->unplacedarray = NULL;
148 PetscFunctionReturn(PETSC_SUCCESS);
149 }
150
VecCopy_Seq(Vec xin,Vec yin)151 PetscErrorCode VecCopy_Seq(Vec xin, Vec yin)
152 {
153 PetscFunctionBegin;
154 if (xin != yin) {
155 const PetscScalar *xa;
156 PetscScalar *ya;
157
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 }
166
VecSwap_Seq(Vec xin,Vec yin)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;
174
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 }
184
VecNorm_Seq(Vec xin,NormType type,PetscReal * z)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;
190
191 PetscFunctionBegin;
192 if (n) {
193 const PetscScalar *xx;
194 const PetscBLASInt one = 1;
195 PetscBLASInt bn = 0;
196
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]);
210
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 }
236
VecView_Seq_ASCII(Vec xin,PetscViewer viewer)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;
243
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;
274
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;
295
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 }
341
342 #include <petscdraw.h>
VecView_Seq_Draw_LG(Vec xin,PetscViewer v)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;
355
356 PetscFunctionBegin;
357 PetscCall(PetscViewerDrawGetDraw(v, 0, &draw));
358 PetscCall(PetscDrawIsNull(draw, &isnull));
359 if (isnull) PetscFunctionReturn(PETSC_SUCCESS);
360
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]);
380
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 }
389
VecView_Seq_Draw(Vec xin,PetscViewer v)390 static PetscErrorCode VecView_Seq_Draw(Vec xin, PetscViewer v)
391 {
392 PetscDraw draw;
393 PetscBool isnull;
394
395 PetscFunctionBegin;
396 PetscCall(PetscViewerDrawGetDraw(v, 0, &draw));
397 PetscCall(PetscDrawIsNull(draw, &isnull));
398 if (isnull) PetscFunctionReturn(PETSC_SUCCESS);
399
400 PetscCall(VecView_Seq_Draw_LG(xin, v));
401 PetscFunctionReturn(PETSC_SUCCESS);
402 }
403
VecView_Seq_Binary(Vec xin,PetscViewer viewer)404 static PetscErrorCode VecView_Seq_Binary(Vec xin, PetscViewer viewer)
405 {
406 return VecView_Binary(xin, viewer);
407 }
408
409 #if defined(PETSC_HAVE_MATLAB)
410 #include <petscmatlab.h>
411 #include <mat.h> /* MATLAB include file */
VecView_Seq_Matlab(Vec vec,PetscViewer viewer)412 PetscErrorCode VecView_Seq_Matlab(Vec vec, PetscViewer viewer)
413 {
414 PetscInt n;
415 const PetscScalar *array;
416
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
426
VecView_Seq(Vec xin,PetscViewer viewer)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
443
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
462
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 }
488
VecGetValues_Seq(Vec xin,PetscInt ni,const PetscInt ix[],PetscScalar y[])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;
493
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 }
507
VecSetValues_Seq(Vec xin,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode m)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;
512
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 }
534
VecSetValuesBlocked_Seq(Vec xin,PetscInt ni,const PetscInt ix[],const PetscScalar yin[],InsertMode m)535 PetscErrorCode VecSetValuesBlocked_Seq(Vec xin, PetscInt ni, const PetscInt ix[], const PetscScalar yin[], InsertMode m)
536 {
537 PetscScalar *xx;
538 PetscInt bs;
539
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];
546
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 }
560
VecResetPreallocationCOO_Seq(Vec x)561 static PetscErrorCode VecResetPreallocationCOO_Seq(Vec x)
562 {
563 Vec_Seq *vs = (Vec_Seq *)x->data;
564
565 PetscFunctionBegin;
566 if (vs) {
567 PetscCall(PetscFree(vs->jmap1)); /* Destroy old stuff */
568 PetscCall(PetscFree(vs->perm1));
569 }
570 PetscFunctionReturn(PETSC_SUCCESS);
571 }
572
VecSetPreallocationCOO_Seq(Vec x,PetscCount coo_n,const PetscInt coo_i[])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;
579
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;
591
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);
595
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));
600
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 }
608
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 }
616
VecSetValuesCOO_Seq(Vec x,const PetscScalar coo_v[],InsertMode imode)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;
623
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 }
635
VecDestroy_Seq(Vec v)636 PetscErrorCode VecDestroy_Seq(Vec v)
637 {
638 Vec_Seq *vs = (Vec_Seq *)v->data;
639
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 }
649
VecSetOption_Seq(Vec v,VecOption op,PetscBool flag)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 }
656
657 // duplicate w to v. v is half-baked, potentially already with arrays allocated.
VecDuplicate_Seq_Private(Vec w,Vec v)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));
664
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 }
673
VecDuplicate_Seq(Vec win,Vec * V)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 }
681
VecReplaceArray_Default_GEMV_Error(Vec v,const PetscScalar * a)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 }
690
VecDuplicateVecs_Seq_GEMV(Vec w,PetscInt m,Vec * V[])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"
695
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 }
713
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 };
806
807 /*
808 Create a VECSEQ with the given layout and array
809
810 Input Parameter:
811 + map - the layout
812 - array - the array on host
813
814 Output Parameter:
815 . V - The vector object
816 */
VecCreateSeqWithLayoutAndArray_Private(PetscLayout map,const PetscScalar array[],Vec * V)817 PetscErrorCode VecCreateSeqWithLayoutAndArray_Private(PetscLayout map, const PetscScalar array[], Vec *V)
818 {
819 PetscMPIInt size;
820
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 }
828
829 /*
830 This is called by VecCreate_Seq() (i.e. VecCreateSeq()) and VecCreateSeqWithArray()
831 */
VecCreate_Seq_Private(Vec v,const PetscScalar array[])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.
837
838 PetscFunctionBegin;
839 PetscCall(PetscNew(&s));
840 v->ops[0] = DvOps;
841
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));
844
845 // allocate multiple vectors together
846 if (mdot_use_gemv || maxpy_use_gemv) v->ops[0].duplicatevecs = VecDuplicateVecs_Seq_GEMV;
847
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;
855
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;
861
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 }
870
871 /*@
872 VecCreateSeqWithArray - Creates a standard,sequential array-style vector,
873 where the user provides the array space to store the vector values.
874
875 Collective
876
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.
882
883 Output Parameter:
884 . V - the vector
885
886 Level: intermediate
887
888 Notes:
889 Use `VecDuplicate()` or `VecDuplicateVecs(`) to form additional vectors of the
890 same type as an existing vector.
891
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.
894
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.
897
898 .seealso: `VecCreateMPIWithArray()`, `VecCreate()`, `VecDuplicate()`, `VecDuplicateVecs()`,
899 `VecCreateGhost()`, `VecCreateSeq()`, `VecPlaceArray()`
900 @*/
VecCreateSeqWithArray(MPI_Comm comm,PetscInt bs,PetscInt n,const PetscScalar array[],Vec * V)901 PetscErrorCode VecCreateSeqWithArray(MPI_Comm comm, PetscInt bs, PetscInt n, const PetscScalar array[], Vec *V)
902 {
903 PetscMPIInt size;
904
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 }
914