Actual source code: bvec2.c
2: /*
3: Implements the sequential vectors.
4: */
6: #include <../src/vec/vec/impls/dvecimpl.h>
7: #include <../src/vec/vec/impls/mpi/pvecimpl.h>
8: #include <petsc/private/glvisviewerimpl.h>
9: #include <petsc/private/glvisvecimpl.h>
10: #include <petscblaslapack.h>
12: #if defined(PETSC_HAVE_HDF5)
13: extern PetscErrorCode VecView_MPI_HDF5(Vec,PetscViewer);
14: #endif
16: PetscErrorCode VecPointwiseMax_Seq(Vec win,Vec xin,Vec yin)
17: {
18: PetscInt n = win->map->n,i;
19: PetscScalar *ww,*xx,*yy; /* cannot make xx or yy const since might be ww */
21: VecGetArrayRead(xin,(const PetscScalar**)&xx);
22: VecGetArrayRead(yin,(const PetscScalar**)&yy);
23: VecGetArray(win,&ww);
25: for (i=0; i<n; i++) ww[i] = PetscMax(PetscRealPart(xx[i]),PetscRealPart(yy[i]));
27: VecRestoreArrayRead(xin,(const PetscScalar**)&xx);
28: VecRestoreArrayRead(yin,(const PetscScalar**)&yy);
29: VecRestoreArray(win,&ww);
30: PetscLogFlops(n);
31: return 0;
32: }
34: PetscErrorCode VecPointwiseMin_Seq(Vec win,Vec xin,Vec yin)
35: {
36: PetscInt n = win->map->n,i;
37: PetscScalar *ww,*xx,*yy; /* cannot make xx or yy const since might be ww */
39: VecGetArrayRead(xin,(const PetscScalar**)&xx);
40: VecGetArrayRead(yin,(const PetscScalar**)&yy);
41: VecGetArray(win,&ww);
43: for (i=0; i<n; i++) ww[i] = PetscMin(PetscRealPart(xx[i]),PetscRealPart(yy[i]));
45: VecRestoreArrayRead(xin,(const PetscScalar**)&xx);
46: VecRestoreArrayRead(yin,(const PetscScalar**)&yy);
47: VecRestoreArray(win,&ww);
48: PetscLogFlops(n);
49: return 0;
50: }
52: PetscErrorCode VecPointwiseMaxAbs_Seq(Vec win,Vec xin,Vec yin)
53: {
54: PetscInt n = win->map->n,i;
55: PetscScalar *ww,*xx,*yy; /* cannot make xx or yy const since might be ww */
57: VecGetArrayRead(xin,(const PetscScalar**)&xx);
58: VecGetArrayRead(yin,(const PetscScalar**)&yy);
59: VecGetArray(win,&ww);
61: for (i=0; i<n; i++) ww[i] = PetscMax(PetscAbsScalar(xx[i]),PetscAbsScalar(yy[i]));
63: PetscLogFlops(n);
64: VecRestoreArrayRead(xin,(const PetscScalar**)&xx);
65: VecRestoreArrayRead(yin,(const PetscScalar**)&yy);
66: VecRestoreArray(win,&ww);
67: return 0;
68: }
70: #include <../src/vec/vec/impls/seq/ftn-kernels/fxtimesy.h>
72: PetscErrorCode VecPointwiseMult_Seq(Vec win,Vec xin,Vec yin)
73: {
74: PetscInt n = win->map->n,i;
75: PetscScalar *ww,*xx,*yy; /* cannot make xx or yy const since might be ww */
77: VecGetArrayRead(xin,(const PetscScalar**)&xx);
78: VecGetArrayRead(yin,(const PetscScalar**)&yy);
79: VecGetArray(win,&ww);
80: if (ww == xx) {
81: for (i=0; i<n; i++) ww[i] *= yy[i];
82: } else if (ww == yy) {
83: for (i=0; i<n; i++) ww[i] *= xx[i];
84: } else {
85: #if defined(PETSC_USE_FORTRAN_KERNEL_XTIMESY)
86: fortranxtimesy_(xx,yy,ww,&n);
87: #else
88: for (i=0; i<n; i++) ww[i] = xx[i] * yy[i];
89: #endif
90: }
91: VecRestoreArrayRead(xin,(const PetscScalar**)&xx);
92: VecRestoreArrayRead(yin,(const PetscScalar**)&yy);
93: VecRestoreArray(win,&ww);
94: PetscLogFlops(n);
95: return 0;
96: }
98: PetscErrorCode VecPointwiseDivide_Seq(Vec win,Vec xin,Vec yin)
99: {
100: PetscInt n = win->map->n,i;
101: PetscScalar *ww,*xx,*yy; /* cannot make xx or yy const since might be ww */
103: VecGetArrayRead(xin,(const PetscScalar**)&xx);
104: VecGetArrayRead(yin,(const PetscScalar**)&yy);
105: VecGetArray(win,&ww);
107: for (i=0; i<n; i++) {
108: if (yy[i] != 0.0) ww[i] = xx[i] / yy[i];
109: else ww[i] = 0.0;
110: }
112: PetscLogFlops(n);
113: VecRestoreArrayRead(xin,(const PetscScalar**)&xx);
114: VecRestoreArrayRead(yin,(const PetscScalar**)&yy);
115: VecRestoreArray(win,&ww);
116: return 0;
117: }
119: PetscErrorCode VecSetRandom_Seq(Vec xin,PetscRandom r)
120: {
121: PetscInt n = xin->map->n,i;
122: PetscScalar *xx;
124: VecGetArrayWrite(xin,&xx);
125: for (i=0; i<n; i++) PetscRandomGetValue(r,&xx[i]);
126: VecRestoreArrayWrite(xin,&xx);
127: return 0;
128: }
130: PetscErrorCode VecGetSize_Seq(Vec vin,PetscInt *size)
131: {
132: *size = vin->map->n;
133: return 0;
134: }
136: PetscErrorCode VecConjugate_Seq(Vec xin)
137: {
138: PetscScalar *x;
139: PetscInt n = xin->map->n;
141: VecGetArray(xin,&x);
142: while (n-->0) {
143: *x = PetscConj(*x);
144: x++;
145: }
146: VecRestoreArray(xin,&x);
147: return 0;
148: }
150: PetscErrorCode VecResetArray_Seq(Vec vin)
151: {
152: Vec_Seq *v = (Vec_Seq*)vin->data;
154: v->array = v->unplacedarray;
155: v->unplacedarray = NULL;
156: return 0;
157: }
159: PetscErrorCode VecCopy_Seq(Vec xin,Vec yin)
160: {
161: PetscScalar *ya;
162: const PetscScalar *xa;
164: if (xin != yin) {
165: VecGetArrayRead(xin,&xa);
166: VecGetArray(yin,&ya);
167: PetscArraycpy(ya,xa,xin->map->n);
168: VecRestoreArrayRead(xin,&xa);
169: VecRestoreArray(yin,&ya);
170: }
171: return 0;
172: }
174: PetscErrorCode VecSwap_Seq(Vec xin,Vec yin)
175: {
176: PetscScalar *ya, *xa;
177: PetscBLASInt one = 1,bn;
179: if (xin != yin) {
180: PetscBLASIntCast(xin->map->n,&bn);
181: VecGetArray(xin,&xa);
182: VecGetArray(yin,&ya);
183: PetscStackCallBLAS("BLASswap",BLASswap_(&bn,xa,&one,ya,&one));
184: VecRestoreArray(xin,&xa);
185: VecRestoreArray(yin,&ya);
186: }
187: return 0;
188: }
190: #include <../src/vec/vec/impls/seq/ftn-kernels/fnorm.h>
192: PetscErrorCode VecNorm_Seq(Vec xin,NormType type,PetscReal *z)
193: {
194: const PetscScalar *xx;
195: PetscInt n = xin->map->n;
196: PetscBLASInt one = 1, bn = 0;
198: PetscBLASIntCast(n,&bn);
199: if (type == NORM_2 || type == NORM_FROBENIUS) {
200: VecGetArrayRead(xin,&xx);
201: #if defined(PETSC_USE_REAL___FP16)
202: PetscStackCallBLAS("BLASnrm2",*z = BLASnrm2_(&bn,xx,&one));
203: #else
204: PetscStackCallBLAS("BLASdot",*z = PetscRealPart(BLASdot_(&bn,xx,&one,xx,&one)));
205: *z = PetscSqrtReal(*z);
206: #endif
207: VecRestoreArrayRead(xin,&xx);
208: PetscLogFlops(PetscMax(2.0*n-1,0.0));
209: } else if (type == NORM_INFINITY) {
210: PetscInt i;
211: PetscReal max = 0.0,tmp;
213: VecGetArrayRead(xin,&xx);
214: for (i=0; i<n; i++) {
215: if ((tmp = PetscAbsScalar(*xx)) > max) max = tmp;
216: /* check special case of tmp == NaN */
217: if (tmp != tmp) {max = tmp; break;}
218: xx++;
219: }
220: VecRestoreArrayRead(xin,&xx);
221: *z = max;
222: } else if (type == NORM_1) {
223: #if defined(PETSC_USE_COMPLEX)
224: PetscReal tmp = 0.0;
225: PetscInt i;
226: #endif
227: VecGetArrayRead(xin,&xx);
228: #if defined(PETSC_USE_COMPLEX)
229: /* BLASasum() returns the nonstandard 1 norm of the 1 norm of the complex entries so we provide a custom loop instead */
230: for (i=0; i<n; i++) {
231: tmp += PetscAbsScalar(xx[i]);
232: }
233: *z = tmp;
234: #else
235: PetscStackCallBLAS("BLASasum",*z = BLASasum_(&bn,xx,&one));
236: #endif
237: VecRestoreArrayRead(xin,&xx);
238: PetscLogFlops(PetscMax(n-1.0,0.0));
239: } else if (type == NORM_1_AND_2) {
240: VecNorm_Seq(xin,NORM_1,z);
241: VecNorm_Seq(xin,NORM_2,z+1);
242: }
243: return 0;
244: }
246: PetscErrorCode VecView_Seq_ASCII(Vec xin,PetscViewer viewer)
247: {
248: PetscInt i,n = xin->map->n;
249: const char *name;
250: PetscViewerFormat format;
251: const PetscScalar *xv;
253: VecGetArrayRead(xin,&xv);
254: PetscViewerGetFormat(viewer,&format);
255: if (format == PETSC_VIEWER_ASCII_MATLAB) {
256: PetscObjectGetName((PetscObject)xin,&name);
257: PetscViewerASCIIPrintf(viewer,"%s = [\n",name);
258: for (i=0; i<n; i++) {
259: #if defined(PETSC_USE_COMPLEX)
260: if (PetscImaginaryPart(xv[i]) > 0.0) {
261: PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei\n",(double)PetscRealPart(xv[i]),(double)PetscImaginaryPart(xv[i]));
262: } else if (PetscImaginaryPart(xv[i]) < 0.0) {
263: PetscViewerASCIIPrintf(viewer,"%18.16e - %18.16ei\n",(double)PetscRealPart(xv[i]),-(double)PetscImaginaryPart(xv[i]));
264: } else {
265: PetscViewerASCIIPrintf(viewer,"%18.16e\n",(double)PetscRealPart(xv[i]));
266: }
267: #else
268: PetscViewerASCIIPrintf(viewer,"%18.16e\n",(double) xv[i]);
269: #endif
270: }
271: PetscViewerASCIIPrintf(viewer,"];\n");
272: } else if (format == PETSC_VIEWER_ASCII_SYMMODU) {
273: for (i=0; i<n; i++) {
274: #if defined(PETSC_USE_COMPLEX)
275: PetscViewerASCIIPrintf(viewer,"%18.16e %18.16e\n",(double)PetscRealPart(xv[i]),(double)PetscImaginaryPart(xv[i]));
276: #else
277: PetscViewerASCIIPrintf(viewer,"%18.16e\n",(double)xv[i]);
278: #endif
279: }
280: } else if (format == PETSC_VIEWER_ASCII_VTK_DEPRECATED || format == PETSC_VIEWER_ASCII_VTK_CELL_DEPRECATED) {
281: /*
282: state 0: No header has been output
283: state 1: Only POINT_DATA has been output
284: state 2: Only CELL_DATA has been output
285: state 3: Output both, POINT_DATA last
286: state 4: Output both, CELL_DATA last
287: */
288: static PetscInt stateId = -1;
289: int outputState = 0;
290: PetscBool hasState;
291: int doOutput = 0;
292: PetscInt bs, b;
294: if (stateId < 0) {
295: PetscObjectComposedDataRegister(&stateId);
296: }
297: PetscObjectComposedDataGetInt((PetscObject) viewer, stateId, outputState, hasState);
298: if (!hasState) outputState = 0;
299: PetscObjectGetName((PetscObject) xin, &name);
300: VecGetBlockSize(xin, &bs);
302: if (format == PETSC_VIEWER_ASCII_VTK_DEPRECATED) {
303: if (outputState == 0) {
304: outputState = 1;
305: doOutput = 1;
306: } else if (outputState == 1) doOutput = 0;
307: else if (outputState == 2) {
308: outputState = 3;
309: doOutput = 1;
310: } else if (outputState == 3) doOutput = 0;
313: if (doOutput) {
314: PetscViewerASCIIPrintf(viewer, "POINT_DATA %" PetscInt_FMT "\n", n/bs);
315: }
316: } else {
317: if (outputState == 0) {
318: outputState = 2;
319: doOutput = 1;
320: } else if (outputState == 1) {
321: outputState = 4;
322: doOutput = 1;
323: } else if (outputState == 2) doOutput = 0;
325: else if (outputState == 4) doOutput = 0;
327: if (doOutput) {
328: PetscViewerASCIIPrintf(viewer, "CELL_DATA %" PetscInt_FMT "\n", n);
329: }
330: }
331: PetscObjectComposedDataSetInt((PetscObject) viewer, stateId, outputState);
332: if (name) {
333: if (bs == 3) {
334: PetscViewerASCIIPrintf(viewer, "VECTORS %s double\n", name);
335: } else {
336: PetscViewerASCIIPrintf(viewer, "SCALARS %s double %" PetscInt_FMT "\n", name, bs);
337: }
338: } else {
339: PetscViewerASCIIPrintf(viewer, "SCALARS scalars double %" PetscInt_FMT "\n", bs);
340: }
341: if (bs != 3) {
342: PetscViewerASCIIPrintf(viewer, "LOOKUP_TABLE default\n");
343: }
344: for (i=0; i<n/bs; i++) {
345: for (b=0; b<bs; b++) {
346: if (b > 0) {
347: PetscViewerASCIIPrintf(viewer," ");
348: }
349: #if !defined(PETSC_USE_COMPLEX)
350: PetscViewerASCIIPrintf(viewer,"%g",(double)xv[i*bs+b]);
351: #endif
352: }
353: PetscViewerASCIIPrintf(viewer,"\n");
354: }
355: } else if (format == PETSC_VIEWER_ASCII_VTK_COORDS_DEPRECATED) {
356: PetscInt bs, b;
358: VecGetBlockSize(xin, &bs);
360: for (i=0; i<n/bs; i++) {
361: for (b=0; b<bs; b++) {
362: if (b > 0) {
363: PetscViewerASCIIPrintf(viewer," ");
364: }
365: #if !defined(PETSC_USE_COMPLEX)
366: PetscViewerASCIIPrintf(viewer,"%g",(double)xv[i*bs+b]);
367: #endif
368: }
369: for (b=bs; b<3; b++) {
370: PetscViewerASCIIPrintf(viewer," 0.0");
371: }
372: PetscViewerASCIIPrintf(viewer,"\n");
373: }
374: } else if (format == PETSC_VIEWER_ASCII_PCICE) {
375: PetscInt bs, b;
377: VecGetBlockSize(xin, &bs);
379: PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT "\n", xin->map->N/bs);
380: for (i=0; i<n/bs; i++) {
381: PetscViewerASCIIPrintf(viewer,"%7" PetscInt_FMT " ", i+1);
382: for (b=0; b<bs; b++) {
383: if (b > 0) {
384: PetscViewerASCIIPrintf(viewer," ");
385: }
386: #if !defined(PETSC_USE_COMPLEX)
387: PetscViewerASCIIPrintf(viewer,"% 12.5E",(double)xv[i*bs+b]);
388: #endif
389: }
390: PetscViewerASCIIPrintf(viewer,"\n");
391: }
392: } else if (format == PETSC_VIEWER_ASCII_GLVIS) {
393: /* GLVis ASCII visualization/dump: this function mimics mfem::GridFunction::Save() */
394: const PetscScalar *array;
395: PetscInt i,n,vdim, ordering = 1; /* mfem::FiniteElementSpace::Ordering::byVDIM */
396: PetscContainer glvis_container;
397: PetscViewerGLVisVecInfo glvis_vec_info;
398: PetscViewerGLVisInfo glvis_info;
400: /* mfem::FiniteElementSpace::Save() */
401: VecGetBlockSize(xin,&vdim);
402: PetscViewerASCIIPrintf(viewer,"FiniteElementSpace\n");
403: PetscObjectQuery((PetscObject)xin,"_glvis_info_container",(PetscObject*)&glvis_container);
405: PetscContainerGetPointer(glvis_container,(void**)&glvis_vec_info);
406: PetscViewerASCIIPrintf(viewer,"%s\n",glvis_vec_info->fec_type);
407: PetscViewerASCIIPrintf(viewer,"VDim: %" PetscInt_FMT "\n",vdim);
408: PetscViewerASCIIPrintf(viewer,"Ordering: %" PetscInt_FMT "\n",ordering);
409: PetscViewerASCIIPrintf(viewer,"\n");
410: /* mfem::Vector::Print() */
411: PetscObjectQuery((PetscObject)viewer,"_glvis_info_container",(PetscObject*)&glvis_container);
413: PetscContainerGetPointer(glvis_container,(void**)&glvis_info);
414: if (glvis_info->enabled) {
415: VecGetLocalSize(xin,&n);
416: VecGetArrayRead(xin,&array);
417: for (i=0;i<n;i++) {
418: PetscViewerASCIIPrintf(viewer,glvis_info->fmt,(double)PetscRealPart(array[i]));
419: PetscViewerASCIIPrintf(viewer,"\n");
420: }
421: VecRestoreArrayRead(xin,&array);
422: }
423: } else if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
424: /* No info */
425: } else {
426: for (i=0; i<n; i++) {
427: if (format == PETSC_VIEWER_ASCII_INDEX) {
428: PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT ": ",i);
429: }
430: #if defined(PETSC_USE_COMPLEX)
431: if (PetscImaginaryPart(xv[i]) > 0.0) {
432: PetscViewerASCIIPrintf(viewer,"%g + %g i\n",(double)PetscRealPart(xv[i]),(double)PetscImaginaryPart(xv[i]));
433: } else if (PetscImaginaryPart(xv[i]) < 0.0) {
434: PetscViewerASCIIPrintf(viewer,"%g - %g i\n",(double)PetscRealPart(xv[i]),-(double)PetscImaginaryPart(xv[i]));
435: } else {
436: PetscViewerASCIIPrintf(viewer,"%g\n",(double)PetscRealPart(xv[i]));
437: }
438: #else
439: PetscViewerASCIIPrintf(viewer,"%g\n",(double)xv[i]);
440: #endif
441: }
442: }
443: PetscViewerFlush(viewer);
444: VecRestoreArrayRead(xin,&xv);
445: return 0;
446: }
448: #include <petscdraw.h>
449: PetscErrorCode VecView_Seq_Draw_LG(Vec xin,PetscViewer v)
450: {
451: PetscDraw draw;
452: PetscBool isnull;
453: PetscDrawLG lg;
454: PetscInt i,c,bs = PetscAbs(xin->map->bs),n = xin->map->n/bs;
455: const PetscScalar *xv;
456: PetscReal *xx,*yy,xmin,xmax,h;
457: int colors[] = {PETSC_DRAW_RED};
458: PetscViewerFormat format;
459: PetscDrawAxis axis;
461: PetscViewerDrawGetDraw(v,0,&draw);
462: PetscDrawIsNull(draw,&isnull);
463: if (isnull) return 0;
465: PetscViewerGetFormat(v,&format);
466: PetscMalloc2(n,&xx,n,&yy);
467: VecGetArrayRead(xin,&xv);
468: for (c=0; c<bs; c++) {
469: PetscViewerDrawGetDrawLG(v,c,&lg);
470: PetscDrawLGReset(lg);
471: PetscDrawLGSetDimension(lg,1);
472: PetscDrawLGSetColors(lg,colors);
473: if (format == PETSC_VIEWER_DRAW_LG_XRANGE) {
474: PetscDrawLGGetAxis(lg,&axis);
475: PetscDrawAxisGetLimits(axis,&xmin,&xmax,NULL,NULL);
476: h = (xmax - xmin)/n;
477: for (i=0; i<n; i++) xx[i] = i*h + 0.5*h; /* cell center */
478: } else {
479: for (i=0; i<n; i++) xx[i] = (PetscReal)i;
480: }
481: for (i=0; i<n; i++) yy[i] = PetscRealPart(xv[c + i*bs]);
483: PetscDrawLGAddPoints(lg,n,&xx,&yy);
484: PetscDrawLGDraw(lg);
485: PetscDrawLGSave(lg);
486: }
487: VecRestoreArrayRead(xin,&xv);
488: PetscFree2(xx,yy);
489: return 0;
490: }
492: PetscErrorCode VecView_Seq_Draw(Vec xin,PetscViewer v)
493: {
494: PetscDraw draw;
495: PetscBool isnull;
497: PetscViewerDrawGetDraw(v,0,&draw);
498: PetscDrawIsNull(draw,&isnull);
499: if (isnull) return 0;
501: VecView_Seq_Draw_LG(xin,v);
502: return 0;
503: }
505: PetscErrorCode VecView_Seq_Binary(Vec xin,PetscViewer viewer)
506: {
507: return VecView_Binary(xin,viewer);
508: }
510: #if defined(PETSC_HAVE_MATLAB_ENGINE)
511: #include <petscmatlab.h>
512: #include <mat.h> /* MATLAB include file */
513: PetscErrorCode VecView_Seq_Matlab(Vec vec,PetscViewer viewer)
514: {
515: PetscInt n;
516: const PetscScalar *array;
518: VecGetLocalSize(vec,&n);
519: PetscObjectName((PetscObject)vec);
520: VecGetArrayRead(vec,&array);
521: PetscViewerMatlabPutArray(viewer,n,1,array,((PetscObject)vec)->name);
522: VecRestoreArrayRead(vec,&array);
523: return 0;
524: }
525: #endif
527: PETSC_EXTERN PetscErrorCode VecView_Seq(Vec xin,PetscViewer viewer)
528: {
529: PetscBool isdraw,iascii,issocket,isbinary;
530: #if defined(PETSC_HAVE_MATHEMATICA)
531: PetscBool ismathematica;
532: #endif
533: #if defined(PETSC_HAVE_MATLAB_ENGINE)
534: PetscBool ismatlab;
535: #endif
536: #if defined(PETSC_HAVE_HDF5)
537: PetscBool ishdf5;
538: #endif
539: PetscBool isglvis;
540: #if defined(PETSC_HAVE_ADIOS)
541: PetscBool isadios;
542: #endif
544: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
545: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
546: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSOCKET,&issocket);
547: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
548: #if defined(PETSC_HAVE_MATHEMATICA)
549: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERMATHEMATICA,&ismathematica);
550: #endif
551: #if defined(PETSC_HAVE_HDF5)
552: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
553: #endif
554: #if defined(PETSC_HAVE_MATLAB_ENGINE)
555: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERMATLAB,&ismatlab);
556: #endif
557: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERGLVIS,&isglvis);
558: #if defined(PETSC_HAVE_ADIOS)
559: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERADIOS,&isadios);
560: #endif
562: if (isdraw) {
563: VecView_Seq_Draw(xin,viewer);
564: } else if (iascii) {
565: VecView_Seq_ASCII(xin,viewer);
566: } else if (isbinary) {
567: VecView_Seq_Binary(xin,viewer);
568: #if defined(PETSC_HAVE_MATHEMATICA)
569: } else if (ismathematica) {
570: PetscViewerMathematicaPutVector(viewer,xin);
571: #endif
572: #if defined(PETSC_HAVE_HDF5)
573: } else if (ishdf5) {
574: VecView_MPI_HDF5(xin,viewer); /* Reusing VecView_MPI_HDF5 ... don't want code duplication*/
575: #endif
576: #if defined(PETSC_HAVE_ADIOS)
577: } else if (isadios) {
578: VecView_MPI_ADIOS(xin,viewer); /* Reusing VecView_MPI_ADIOS ... don't want code duplication*/
579: #endif
580: #if defined(PETSC_HAVE_MATLAB_ENGINE)
581: } else if (ismatlab) {
582: VecView_Seq_Matlab(xin,viewer);
583: #endif
584: } else if (isglvis) {
585: VecView_GLVis(xin,viewer);
586: }
587: return 0;
588: }
590: PetscErrorCode VecGetValues_Seq(Vec xin,PetscInt ni,const PetscInt ix[],PetscScalar y[])
591: {
592: const PetscScalar *xx;
593: PetscInt i;
595: VecGetArrayRead(xin,&xx);
596: for (i=0; i<ni; i++) {
597: if (xin->stash.ignorenegidx && ix[i] < 0) continue;
598: if (PetscDefined(USE_DEBUG)) {
601: }
602: y[i] = xx[ix[i]];
603: }
604: VecRestoreArrayRead(xin,&xx);
605: return 0;
606: }
608: PetscErrorCode VecSetValues_Seq(Vec xin,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode m)
609: {
610: PetscScalar *xx;
611: PetscInt i;
613: VecGetArray(xin,&xx);
614: if (m == INSERT_VALUES) {
615: for (i=0; i<ni; i++) {
616: if (xin->stash.ignorenegidx && ix[i] < 0) continue;
617: if (PetscDefined(USE_DEBUG)) {
620: }
621: xx[ix[i]] = y[i];
622: }
623: } else {
624: for (i=0; i<ni; i++) {
625: if (xin->stash.ignorenegidx && ix[i] < 0) continue;
626: if (PetscDefined(USE_DEBUG)) {
629: }
630: xx[ix[i]] += y[i];
631: }
632: }
633: VecRestoreArray(xin,&xx);
634: return 0;
635: }
637: PetscErrorCode VecSetValuesBlocked_Seq(Vec xin,PetscInt ni,const PetscInt ix[],const PetscScalar yin[],InsertMode m)
638: {
639: PetscScalar *xx,*y = (PetscScalar*)yin;
640: PetscInt i,bs,start,j;
642: /*
643: For optimization could treat bs = 2, 3, 4, 5 as special cases with loop unrolling
644: */
645: VecGetBlockSize(xin,&bs);
646: VecGetArray(xin,&xx);
647: if (m == INSERT_VALUES) {
648: for (i=0; i<ni; i++, y+=bs) {
649: start = bs*ix[i];
650: if (start < 0) continue;
652: for (j=0; j<bs; j++) xx[start+j] = y[j];
653: }
654: } else {
655: for (i=0; i<ni; i++, y+=bs) {
656: start = bs*ix[i];
657: if (start < 0) continue;
659: for (j=0; j<bs; j++) xx[start+j] += y[j];
660: }
661: }
662: VecRestoreArray(xin,&xx);
663: return 0;
664: }
666: PetscErrorCode VecDestroy_Seq(Vec v)
667: {
668: Vec_Seq *vs = (Vec_Seq*)v->data;
670: #if defined(PETSC_USE_LOG)
671: PetscLogObjectState((PetscObject)v,"Length=%" PetscInt_FMT,v->map->n);
672: #endif
673: if (vs) PetscFree(vs->array_allocated);
674: PetscFree(v->data);
675: return 0;
676: }
678: PetscErrorCode VecSetOption_Seq(Vec v,VecOption op,PetscBool flag)
679: {
680: if (op == VEC_IGNORE_NEGATIVE_INDICES) v->stash.ignorenegidx = flag;
681: return 0;
682: }
684: PetscErrorCode VecDuplicate_Seq(Vec win,Vec *V)
685: {
686: VecCreate(PetscObjectComm((PetscObject)win),V);
687: VecSetSizes(*V,win->map->n,win->map->n);
688: VecSetType(*V,((PetscObject)win)->type_name);
689: PetscLayoutReference(win->map,&(*V)->map);
690: PetscObjectListDuplicate(((PetscObject)win)->olist,&((PetscObject)(*V))->olist);
691: PetscFunctionListDuplicate(((PetscObject)win)->qlist,&((PetscObject)(*V))->qlist);
693: (*V)->ops->view = win->ops->view;
694: (*V)->stash.ignorenegidx = win->stash.ignorenegidx;
695: return 0;
696: }
698: static struct _VecOps DvOps = {
699: PetscDesignatedInitializer(duplicate,VecDuplicate_Seq), /* 1 */
700: PetscDesignatedInitializer(duplicatevecs,VecDuplicateVecs_Default),
701: PetscDesignatedInitializer(destroyvecs,VecDestroyVecs_Default),
702: PetscDesignatedInitializer(dot,VecDot_Seq),
703: PetscDesignatedInitializer(mdot,VecMDot_Seq),
704: PetscDesignatedInitializer(norm,VecNorm_Seq),
705: PetscDesignatedInitializer(tdot,VecTDot_Seq),
706: PetscDesignatedInitializer(mtdot,VecMTDot_Seq),
707: PetscDesignatedInitializer(scale,VecScale_Seq),
708: PetscDesignatedInitializer(copy,VecCopy_Seq), /* 10 */
709: PetscDesignatedInitializer(set,VecSet_Seq),
710: PetscDesignatedInitializer(swap,VecSwap_Seq),
711: PetscDesignatedInitializer(axpy,VecAXPY_Seq),
712: PetscDesignatedInitializer(axpby,VecAXPBY_Seq),
713: PetscDesignatedInitializer(maxpy,VecMAXPY_Seq),
714: PetscDesignatedInitializer(aypx,VecAYPX_Seq),
715: PetscDesignatedInitializer(waxpy,VecWAXPY_Seq),
716: PetscDesignatedInitializer(axpbypcz,VecAXPBYPCZ_Seq),
717: PetscDesignatedInitializer(pointwisemult,VecPointwiseMult_Seq),
718: PetscDesignatedInitializer(pointwisedivide,VecPointwiseDivide_Seq),
719: PetscDesignatedInitializer(setvalues,VecSetValues_Seq), /* 20 */
720: PetscDesignatedInitializer(assemblybegin,NULL),
721: PetscDesignatedInitializer(assemblyend,NULL),
722: PetscDesignatedInitializer(getarray,NULL),
723: PetscDesignatedInitializer(getsize,VecGetSize_Seq),
724: PetscDesignatedInitializer(getlocalsize,VecGetSize_Seq),
725: PetscDesignatedInitializer(restorearray,NULL),
726: PetscDesignatedInitializer(max,VecMax_Seq),
727: PetscDesignatedInitializer(min,VecMin_Seq),
728: PetscDesignatedInitializer(setrandom,VecSetRandom_Seq),
729: PetscDesignatedInitializer(setoption,VecSetOption_Seq), /* 30 */
730: PetscDesignatedInitializer(setvaluesblocked,VecSetValuesBlocked_Seq),
731: PetscDesignatedInitializer(destroy,VecDestroy_Seq),
732: PetscDesignatedInitializer(view,VecView_Seq),
733: PetscDesignatedInitializer(placearray,VecPlaceArray_Seq),
734: PetscDesignatedInitializer(replacearray,VecReplaceArray_Seq),
735: PetscDesignatedInitializer(dot_local,VecDot_Seq),
736: PetscDesignatedInitializer(tdot_local,VecTDot_Seq),
737: PetscDesignatedInitializer(norm_local,VecNorm_Seq),
738: PetscDesignatedInitializer(mdot_local,VecMDot_Seq),
739: PetscDesignatedInitializer(mtdot_local,VecMTDot_Seq), /* 40 */
740: PetscDesignatedInitializer(load,VecLoad_Default),
741: PetscDesignatedInitializer(reciprocal,VecReciprocal_Default),
742: PetscDesignatedInitializer(conjugate,VecConjugate_Seq),
743: PetscDesignatedInitializer(setlocaltoglobalmapping,NULL),
744: PetscDesignatedInitializer(setvalueslocal,NULL),
745: PetscDesignatedInitializer(resetarray,VecResetArray_Seq),
746: PetscDesignatedInitializer(setfromoptions,NULL),
747: PetscDesignatedInitializer(maxpointwisedivide,VecMaxPointwiseDivide_Seq),
748: PetscDesignatedInitializer(pointwisemax,VecPointwiseMax_Seq),
749: PetscDesignatedInitializer(pointwisemaxabs,VecPointwiseMaxAbs_Seq),
750: PetscDesignatedInitializer(pointwisemin,VecPointwiseMin_Seq),
751: PetscDesignatedInitializer(getvalues,VecGetValues_Seq),
752: PetscDesignatedInitializer(sqrt,NULL),
753: PetscDesignatedInitializer(abs,NULL),
754: PetscDesignatedInitializer(exp,NULL),
755: PetscDesignatedInitializer(log,NULL),
756: PetscDesignatedInitializer(shift,NULL),
757: PetscDesignatedInitializer(create,NULL),
758: PetscDesignatedInitializer(stridegather,VecStrideGather_Default),
759: PetscDesignatedInitializer(stridescatter,VecStrideScatter_Default),
760: PetscDesignatedInitializer(dotnorm2,NULL),
761: PetscDesignatedInitializer(getsubvector,NULL),
762: PetscDesignatedInitializer(restoresubvector,NULL),
763: PetscDesignatedInitializer(getarrayread,NULL),
764: PetscDesignatedInitializer(restorearrayread,NULL),
765: PetscDesignatedInitializer(stridesubsetgather,VecStrideSubSetGather_Default),
766: PetscDesignatedInitializer(stridesubsetscatter,VecStrideSubSetScatter_Default),
767: PetscDesignatedInitializer(viewnative,VecView_Seq),
768: PetscDesignatedInitializer(loadnative,NULL),
769: PetscDesignatedInitializer(getlocalvector,NULL),
770: };
772: /*
773: This is called by VecCreate_Seq() (i.e. VecCreateSeq()) and VecCreateSeqWithArray()
774: */
775: PetscErrorCode VecCreate_Seq_Private(Vec v,const PetscScalar array[])
776: {
777: Vec_Seq *s;
779: PetscNewLog(v,&s);
780: PetscMemcpy(v->ops,&DvOps,sizeof(DvOps));
782: v->data = (void*)s;
783: v->petscnative = PETSC_TRUE;
784: s->array = (PetscScalar*)array;
785: s->array_allocated = NULL;
786: if (array) v->offloadmask = PETSC_OFFLOAD_CPU;
788: PetscLayoutSetUp(v->map);
789: PetscObjectChangeTypeName((PetscObject)v,VECSEQ);
790: #if defined(PETSC_HAVE_MATLAB_ENGINE)
791: PetscObjectComposeFunction((PetscObject)v,"PetscMatlabEnginePut_C",VecMatlabEnginePut_Default);
792: PetscObjectComposeFunction((PetscObject)v,"PetscMatlabEngineGet_C",VecMatlabEngineGet_Default);
793: #endif
794: return 0;
795: }
797: /*@C
798: VecCreateSeqWithArray - Creates a standard,sequential array-style vector,
799: where the user provides the array space to store the vector values.
801: Collective
803: Input Parameters:
804: + comm - the communicator, should be PETSC_COMM_SELF
805: . bs - the block size
806: . n - the vector length
807: - array - memory where the vector elements are to be stored.
809: Output Parameter:
810: . V - the vector
812: Notes:
813: Use VecDuplicate() or VecDuplicateVecs() to form additional vectors of the
814: same type as an existing vector.
816: If the user-provided array is NULL, then VecPlaceArray() can be used
817: at a later stage to SET the array for storing the vector values.
819: PETSc does NOT free the array when the vector is destroyed via VecDestroy().
820: The user should not free the array until the vector is destroyed.
822: Level: intermediate
824: .seealso: VecCreateMPIWithArray(), VecCreate(), VecDuplicate(), VecDuplicateVecs(),
825: VecCreateGhost(), VecCreateSeq(), VecPlaceArray()
826: @*/
827: PetscErrorCode VecCreateSeqWithArray(MPI_Comm comm,PetscInt bs,PetscInt n,const PetscScalar array[],Vec *V)
828: {
829: PetscMPIInt size;
831: VecCreate(comm,V);
832: VecSetSizes(*V,n,n);
833: VecSetBlockSize(*V,bs);
834: MPI_Comm_size(comm,&size);
836: VecCreate_Seq_Private(*V,array);
837: return 0;
838: }