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: {
19: PetscInt n = win->map->n,i;
20: PetscScalar *ww,*xx,*yy; /* cannot make xx or yy const since might be ww */
23: VecGetArrayRead(xin,(const PetscScalar**)&xx);
24: VecGetArrayRead(yin,(const PetscScalar**)&yy);
25: VecGetArray(win,&ww);
27: for (i=0; i<n; i++) ww[i] = PetscMax(PetscRealPart(xx[i]),PetscRealPart(yy[i]));
29: VecRestoreArrayRead(xin,(const PetscScalar**)&xx);
30: VecRestoreArrayRead(yin,(const PetscScalar**)&yy);
31: VecRestoreArray(win,&ww);
32: PetscLogFlops(n);
33: return(0);
34: }
36: PetscErrorCode VecPointwiseMin_Seq(Vec win,Vec xin,Vec yin)
37: {
39: PetscInt n = win->map->n,i;
40: PetscScalar *ww,*xx,*yy; /* cannot make xx or yy const since might be ww */
43: VecGetArrayRead(xin,(const PetscScalar**)&xx);
44: VecGetArrayRead(yin,(const PetscScalar**)&yy);
45: VecGetArray(win,&ww);
47: for (i=0; i<n; i++) ww[i] = PetscMin(PetscRealPart(xx[i]),PetscRealPart(yy[i]));
49: VecRestoreArrayRead(xin,(const PetscScalar**)&xx);
50: VecRestoreArrayRead(yin,(const PetscScalar**)&yy);
51: VecRestoreArray(win,&ww);
52: PetscLogFlops(n);
53: return(0);
54: }
56: PetscErrorCode VecPointwiseMaxAbs_Seq(Vec win,Vec xin,Vec yin)
57: {
59: PetscInt n = win->map->n,i;
60: PetscScalar *ww,*xx,*yy; /* cannot make xx or yy const since might be ww */
63: VecGetArrayRead(xin,(const PetscScalar**)&xx);
64: VecGetArrayRead(yin,(const PetscScalar**)&yy);
65: VecGetArray(win,&ww);
67: for (i=0; i<n; i++) ww[i] = PetscMax(PetscAbsScalar(xx[i]),PetscAbsScalar(yy[i]));
69: PetscLogFlops(n);
70: VecRestoreArrayRead(xin,(const PetscScalar**)&xx);
71: VecRestoreArrayRead(yin,(const PetscScalar**)&yy);
72: VecRestoreArray(win,&ww);
73: return(0);
74: }
76: #include <../src/vec/vec/impls/seq/ftn-kernels/fxtimesy.h>
78: PetscErrorCode VecPointwiseMult_Seq(Vec win,Vec xin,Vec yin)
79: {
81: PetscInt n = win->map->n,i;
82: PetscScalar *ww,*xx,*yy; /* cannot make xx or yy const since might be ww */
85: VecGetArrayRead(xin,(const PetscScalar**)&xx);
86: VecGetArrayRead(yin,(const PetscScalar**)&yy);
87: VecGetArray(win,&ww);
88: if (ww == xx) {
89: for (i=0; i<n; i++) ww[i] *= yy[i];
90: } else if (ww == yy) {
91: for (i=0; i<n; i++) ww[i] *= xx[i];
92: } else {
93: #if defined(PETSC_USE_FORTRAN_KERNEL_XTIMESY)
94: fortranxtimesy_(xx,yy,ww,&n);
95: #else
96: for (i=0; i<n; i++) ww[i] = xx[i] * yy[i];
97: #endif
98: }
99: VecRestoreArrayRead(xin,(const PetscScalar**)&xx);
100: VecRestoreArrayRead(yin,(const PetscScalar**)&yy);
101: VecRestoreArray(win,&ww);
102: PetscLogFlops(n);
103: return(0);
104: }
106: PetscErrorCode VecPointwiseDivide_Seq(Vec win,Vec xin,Vec yin)
107: {
109: PetscInt n = win->map->n,i;
110: PetscScalar *ww,*xx,*yy; /* cannot make xx or yy const since might be ww */
113: VecGetArrayRead(xin,(const PetscScalar**)&xx);
114: VecGetArrayRead(yin,(const PetscScalar**)&yy);
115: VecGetArray(win,&ww);
117: for (i=0; i<n; i++) {
118: if (yy[i] != 0.0) ww[i] = xx[i] / yy[i];
119: else ww[i] = 0.0;
120: }
122: PetscLogFlops(n);
123: VecRestoreArrayRead(xin,(const PetscScalar**)&xx);
124: VecRestoreArrayRead(yin,(const PetscScalar**)&yy);
125: VecRestoreArray(win,&ww);
126: return(0);
127: }
129: PetscErrorCode VecSetRandom_Seq(Vec xin,PetscRandom r)
130: {
132: PetscInt n = xin->map->n,i;
133: PetscScalar *xx;
136: VecGetArrayWrite(xin,&xx);
137: for (i=0; i<n; i++) {PetscRandomGetValue(r,&xx[i]);}
138: VecRestoreArrayWrite(xin,&xx);
139: return(0);
140: }
142: PetscErrorCode VecGetSize_Seq(Vec vin,PetscInt *size)
143: {
145: *size = vin->map->n;
146: return(0);
147: }
149: PetscErrorCode VecConjugate_Seq(Vec xin)
150: {
151: PetscScalar *x;
152: PetscInt n = xin->map->n;
156: VecGetArray(xin,&x);
157: while (n-->0) {
158: *x = PetscConj(*x);
159: x++;
160: }
161: VecRestoreArray(xin,&x);
162: return(0);
163: }
165: PetscErrorCode VecResetArray_Seq(Vec vin)
166: {
167: Vec_Seq *v = (Vec_Seq*)vin->data;
170: v->array = v->unplacedarray;
171: v->unplacedarray = NULL;
172: return(0);
173: }
175: PetscErrorCode VecCopy_Seq(Vec xin,Vec yin)
176: {
177: PetscScalar *ya;
178: const PetscScalar *xa;
179: PetscErrorCode ierr;
182: if (xin != yin) {
183: VecGetArrayRead(xin,&xa);
184: VecGetArray(yin,&ya);
185: PetscArraycpy(ya,xa,xin->map->n);
186: VecRestoreArrayRead(xin,&xa);
187: VecRestoreArray(yin,&ya);
188: }
189: return(0);
190: }
192: PetscErrorCode VecSwap_Seq(Vec xin,Vec yin)
193: {
194: PetscScalar *ya, *xa;
196: PetscBLASInt one = 1,bn;
199: if (xin != yin) {
200: PetscBLASIntCast(xin->map->n,&bn);
201: VecGetArray(xin,&xa);
202: VecGetArray(yin,&ya);
203: PetscStackCallBLAS("BLASswap",BLASswap_(&bn,xa,&one,ya,&one));
204: VecRestoreArray(xin,&xa);
205: VecRestoreArray(yin,&ya);
206: }
207: return(0);
208: }
210: #include <../src/vec/vec/impls/seq/ftn-kernels/fnorm.h>
212: PetscErrorCode VecNorm_Seq(Vec xin,NormType type,PetscReal *z)
213: {
214: const PetscScalar *xx;
215: PetscErrorCode ierr;
216: PetscInt n = xin->map->n;
217: PetscBLASInt one = 1, bn = 0;
220: PetscBLASIntCast(n,&bn);
221: if (type == NORM_2 || type == NORM_FROBENIUS) {
222: VecGetArrayRead(xin,&xx);
223: #if defined(PETSC_USE_REAL___FP16)
224: PetscStackCallBLAS("BLASnrm2",*z = BLASnrm2_(&bn,xx,&one));
225: #else
226: PetscStackCallBLAS("BLASdot",*z = PetscRealPart(BLASdot_(&bn,xx,&one,xx,&one)));
227: *z = PetscSqrtReal(*z);
228: #endif
229: VecRestoreArrayRead(xin,&xx);
230: PetscLogFlops(PetscMax(2.0*n-1,0.0));
231: } else if (type == NORM_INFINITY) {
232: PetscInt i;
233: PetscReal max = 0.0,tmp;
235: VecGetArrayRead(xin,&xx);
236: for (i=0; i<n; i++) {
237: if ((tmp = PetscAbsScalar(*xx)) > max) max = tmp;
238: /* check special case of tmp == NaN */
239: if (tmp != tmp) {max = tmp; break;}
240: xx++;
241: }
242: VecRestoreArrayRead(xin,&xx);
243: *z = max;
244: } else if (type == NORM_1) {
245: #if defined(PETSC_USE_COMPLEX)
246: PetscReal tmp = 0.0;
247: PetscInt i;
248: #endif
249: VecGetArrayRead(xin,&xx);
250: #if defined(PETSC_USE_COMPLEX)
251: /* BLASasum() returns the nonstandard 1 norm of the 1 norm of the complex entries so we provide a custom loop instead */
252: for (i=0; i<n; i++) {
253: tmp += PetscAbsScalar(xx[i]);
254: }
255: *z = tmp;
256: #else
257: PetscStackCallBLAS("BLASasum",*z = BLASasum_(&bn,xx,&one));
258: #endif
259: VecRestoreArrayRead(xin,&xx);
260: PetscLogFlops(PetscMax(n-1.0,0.0));
261: } else if (type == NORM_1_AND_2) {
262: VecNorm_Seq(xin,NORM_1,z);
263: VecNorm_Seq(xin,NORM_2,z+1);
264: }
265: return(0);
266: }
268: PetscErrorCode VecView_Seq_ASCII(Vec xin,PetscViewer viewer)
269: {
270: PetscErrorCode ierr;
271: PetscInt i,n = xin->map->n;
272: const char *name;
273: PetscViewerFormat format;
274: const PetscScalar *xv;
277: VecGetArrayRead(xin,&xv);
278: PetscViewerGetFormat(viewer,&format);
279: if (format == PETSC_VIEWER_ASCII_MATLAB) {
280: PetscObjectGetName((PetscObject)xin,&name);
281: PetscViewerASCIIPrintf(viewer,"%s = [\n",name);
282: for (i=0; i<n; i++) {
283: #if defined(PETSC_USE_COMPLEX)
284: if (PetscImaginaryPart(xv[i]) > 0.0) {
285: PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei\n",(double)PetscRealPart(xv[i]),(double)PetscImaginaryPart(xv[i]));
286: } else if (PetscImaginaryPart(xv[i]) < 0.0) {
287: PetscViewerASCIIPrintf(viewer,"%18.16e - %18.16ei\n",(double)PetscRealPart(xv[i]),-(double)PetscImaginaryPart(xv[i]));
288: } else {
289: PetscViewerASCIIPrintf(viewer,"%18.16e\n",(double)PetscRealPart(xv[i]));
290: }
291: #else
292: PetscViewerASCIIPrintf(viewer,"%18.16e\n",(double) xv[i]);
293: #endif
294: }
295: PetscViewerASCIIPrintf(viewer,"];\n");
296: } else if (format == PETSC_VIEWER_ASCII_SYMMODU) {
297: for (i=0; i<n; i++) {
298: #if defined(PETSC_USE_COMPLEX)
299: PetscViewerASCIIPrintf(viewer,"%18.16e %18.16e\n",(double)PetscRealPart(xv[i]),(double)PetscImaginaryPart(xv[i]));
300: #else
301: PetscViewerASCIIPrintf(viewer,"%18.16e\n",(double)xv[i]);
302: #endif
303: }
304: } else if (format == PETSC_VIEWER_ASCII_VTK_DEPRECATED || format == PETSC_VIEWER_ASCII_VTK_CELL_DEPRECATED) {
305: /*
306: state 0: No header has been output
307: state 1: Only POINT_DATA has been output
308: state 2: Only CELL_DATA has been output
309: state 3: Output both, POINT_DATA last
310: state 4: Output both, CELL_DATA last
311: */
312: static PetscInt stateId = -1;
313: int outputState = 0;
314: PetscBool hasState;
315: int doOutput = 0;
316: PetscInt bs, b;
318: if (stateId < 0) {
319: PetscObjectComposedDataRegister(&stateId);
320: }
321: PetscObjectComposedDataGetInt((PetscObject) viewer, stateId, outputState, hasState);
322: if (!hasState) outputState = 0;
323: PetscObjectGetName((PetscObject) xin, &name);
324: VecGetBlockSize(xin, &bs);
325: if ((bs < 1) || (bs > 3)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "VTK can only handle 3D objects, but vector dimension is %d", bs);
326: if (format == PETSC_VIEWER_ASCII_VTK_DEPRECATED) {
327: if (outputState == 0) {
328: outputState = 1;
329: doOutput = 1;
330: } else if (outputState == 1) doOutput = 0;
331: else if (outputState == 2) {
332: outputState = 3;
333: doOutput = 1;
334: } else if (outputState == 3) doOutput = 0;
335: else if (outputState == 4) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Tried to output POINT_DATA again after intervening CELL_DATA");
337: if (doOutput) {
338: PetscViewerASCIIPrintf(viewer, "POINT_DATA %d\n", n/bs);
339: }
340: } else {
341: if (outputState == 0) {
342: outputState = 2;
343: doOutput = 1;
344: } else if (outputState == 1) {
345: outputState = 4;
346: doOutput = 1;
347: } else if (outputState == 2) doOutput = 0;
348: else if (outputState == 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Tried to output CELL_DATA again after intervening POINT_DATA");
349: else if (outputState == 4) doOutput = 0;
351: if (doOutput) {
352: PetscViewerASCIIPrintf(viewer, "CELL_DATA %d\n", n);
353: }
354: }
355: PetscObjectComposedDataSetInt((PetscObject) viewer, stateId, outputState);
356: if (name) {
357: if (bs == 3) {
358: PetscViewerASCIIPrintf(viewer, "VECTORS %s double\n", name);
359: } else {
360: PetscViewerASCIIPrintf(viewer, "SCALARS %s double %d\n", name, bs);
361: }
362: } else {
363: PetscViewerASCIIPrintf(viewer, "SCALARS scalars double %d\n", bs);
364: }
365: if (bs != 3) {
366: PetscViewerASCIIPrintf(viewer, "LOOKUP_TABLE default\n");
367: }
368: for (i=0; i<n/bs; i++) {
369: for (b=0; b<bs; b++) {
370: if (b > 0) {
371: PetscViewerASCIIPrintf(viewer," ");
372: }
373: #if !defined(PETSC_USE_COMPLEX)
374: PetscViewerASCIIPrintf(viewer,"%g",(double)xv[i*bs+b]);
375: #endif
376: }
377: PetscViewerASCIIPrintf(viewer,"\n");
378: }
379: } else if (format == PETSC_VIEWER_ASCII_VTK_COORDS_DEPRECATED) {
380: PetscInt bs, b;
382: VecGetBlockSize(xin, &bs);
383: if ((bs < 1) || (bs > 3)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "VTK can only handle 3D objects, but vector dimension is %d", bs);
384: for (i=0; i<n/bs; i++) {
385: for (b=0; b<bs; b++) {
386: if (b > 0) {
387: PetscViewerASCIIPrintf(viewer," ");
388: }
389: #if !defined(PETSC_USE_COMPLEX)
390: PetscViewerASCIIPrintf(viewer,"%g",(double)xv[i*bs+b]);
391: #endif
392: }
393: for (b=bs; b<3; b++) {
394: PetscViewerASCIIPrintf(viewer," 0.0");
395: }
396: PetscViewerASCIIPrintf(viewer,"\n");
397: }
398: } else if (format == PETSC_VIEWER_ASCII_PCICE) {
399: PetscInt bs, b;
401: VecGetBlockSize(xin, &bs);
402: if ((bs < 1) || (bs > 3)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "PCICE can only handle up to 3D objects, but vector dimension is %d", bs);
403: PetscViewerASCIIPrintf(viewer,"%D\n", xin->map->N/bs);
404: for (i=0; i<n/bs; i++) {
405: PetscViewerASCIIPrintf(viewer,"%7D ", i+1);
406: for (b=0; b<bs; b++) {
407: if (b > 0) {
408: PetscViewerASCIIPrintf(viewer," ");
409: }
410: #if !defined(PETSC_USE_COMPLEX)
411: PetscViewerASCIIPrintf(viewer,"% 12.5E",(double)xv[i*bs+b]);
412: #endif
413: }
414: PetscViewerASCIIPrintf(viewer,"\n");
415: }
416: } else if (format == PETSC_VIEWER_ASCII_GLVIS) {
417: /* GLVis ASCII visualization/dump: this function mimics mfem::GridFunction::Save() */
418: const PetscScalar *array;
419: PetscInt i,n,vdim, ordering = 1; /* mfem::FiniteElementSpace::Ordering::byVDIM */
420: PetscContainer glvis_container;
421: PetscViewerGLVisVecInfo glvis_vec_info;
422: PetscViewerGLVisInfo glvis_info;
423: PetscErrorCode ierr;
425: /* mfem::FiniteElementSpace::Save() */
426: VecGetBlockSize(xin,&vdim);
427: PetscViewerASCIIPrintf(viewer,"FiniteElementSpace\n");
428: PetscObjectQuery((PetscObject)xin,"_glvis_info_container",(PetscObject*)&glvis_container);
429: if (!glvis_container) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_PLIB,"Missing GLVis container");
430: PetscContainerGetPointer(glvis_container,(void**)&glvis_vec_info);
431: PetscViewerASCIIPrintf(viewer,"%s\n",glvis_vec_info->fec_type);
432: PetscViewerASCIIPrintf(viewer,"VDim: %d\n",vdim);
433: PetscViewerASCIIPrintf(viewer,"Ordering: %d\n",ordering);
434: PetscViewerASCIIPrintf(viewer,"\n");
435: /* mfem::Vector::Print() */
436: PetscObjectQuery((PetscObject)viewer,"_glvis_info_container",(PetscObject*)&glvis_container);
437: if (!glvis_container) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_PLIB,"Missing GLVis container");
438: PetscContainerGetPointer(glvis_container,(void**)&glvis_info);
439: if (glvis_info->enabled) {
440: VecGetLocalSize(xin,&n);
441: VecGetArrayRead(xin,&array);
442: for (i=0;i<n;i++) {
443: PetscViewerASCIIPrintf(viewer,glvis_info->fmt,(double)PetscRealPart(array[i]));
444: PetscViewerASCIIPrintf(viewer,"\n");
445: }
446: VecRestoreArrayRead(xin,&array);
447: }
448: } else if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
449: /* No info */
450: } else {
451: for (i=0; i<n; i++) {
452: if (format == PETSC_VIEWER_ASCII_INDEX) {
453: PetscViewerASCIIPrintf(viewer,"%D: ",i);
454: }
455: #if defined(PETSC_USE_COMPLEX)
456: if (PetscImaginaryPart(xv[i]) > 0.0) {
457: PetscViewerASCIIPrintf(viewer,"%g + %g i\n",(double)PetscRealPart(xv[i]),(double)PetscImaginaryPart(xv[i]));
458: } else if (PetscImaginaryPart(xv[i]) < 0.0) {
459: PetscViewerASCIIPrintf(viewer,"%g - %g i\n",(double)PetscRealPart(xv[i]),-(double)PetscImaginaryPart(xv[i]));
460: } else {
461: PetscViewerASCIIPrintf(viewer,"%g\n",(double)PetscRealPart(xv[i]));
462: }
463: #else
464: PetscViewerASCIIPrintf(viewer,"%g\n",(double)xv[i]);
465: #endif
466: }
467: }
468: PetscViewerFlush(viewer);
469: VecRestoreArrayRead(xin,&xv);
470: return(0);
471: }
473: #include <petscdraw.h>
474: PetscErrorCode VecView_Seq_Draw_LG(Vec xin,PetscViewer v)
475: {
476: PetscDraw draw;
477: PetscBool isnull;
478: PetscDrawLG lg;
479: PetscErrorCode ierr;
480: PetscInt i,c,bs = PetscAbs(xin->map->bs),n = xin->map->n/bs;
481: const PetscScalar *xv;
482: PetscReal *xx,*yy,xmin,xmax,h;
483: int colors[] = {PETSC_DRAW_RED};
484: PetscViewerFormat format;
485: PetscDrawAxis axis;
488: PetscViewerDrawGetDraw(v,0,&draw);
489: PetscDrawIsNull(draw,&isnull);
490: if (isnull) return(0);
492: PetscViewerGetFormat(v,&format);
493: PetscMalloc2(n,&xx,n,&yy);
494: VecGetArrayRead(xin,&xv);
495: for (c=0; c<bs; c++) {
496: PetscViewerDrawGetDrawLG(v,c,&lg);
497: PetscDrawLGReset(lg);
498: PetscDrawLGSetDimension(lg,1);
499: PetscDrawLGSetColors(lg,colors);
500: if (format == PETSC_VIEWER_DRAW_LG_XRANGE) {
501: PetscDrawLGGetAxis(lg,&axis);
502: PetscDrawAxisGetLimits(axis,&xmin,&xmax,NULL,NULL);
503: h = (xmax - xmin)/n;
504: for (i=0; i<n; i++) xx[i] = i*h + 0.5*h; /* cell center */
505: } else {
506: for (i=0; i<n; i++) xx[i] = (PetscReal)i;
507: }
508: for (i=0; i<n; i++) yy[i] = PetscRealPart(xv[c + i*bs]);
510: PetscDrawLGAddPoints(lg,n,&xx,&yy);
511: PetscDrawLGDraw(lg);
512: PetscDrawLGSave(lg);
513: }
514: VecRestoreArrayRead(xin,&xv);
515: PetscFree2(xx,yy);
516: return(0);
517: }
519: PetscErrorCode VecView_Seq_Draw(Vec xin,PetscViewer v)
520: {
521: PetscErrorCode ierr;
522: PetscDraw draw;
523: PetscBool isnull;
526: PetscViewerDrawGetDraw(v,0,&draw);
527: PetscDrawIsNull(draw,&isnull);
528: if (isnull) return(0);
530: VecView_Seq_Draw_LG(xin,v);
531: return(0);
532: }
534: PetscErrorCode VecView_Seq_Binary(Vec xin,PetscViewer viewer)
535: {
536: return VecView_Binary(xin,viewer);
537: }
539: #if defined(PETSC_HAVE_MATLAB_ENGINE)
540: #include <petscmatlab.h>
541: #include <mat.h> /* MATLAB include file */
542: PetscErrorCode VecView_Seq_Matlab(Vec vec,PetscViewer viewer)
543: {
544: PetscErrorCode ierr;
545: PetscInt n;
546: const PetscScalar *array;
549: VecGetLocalSize(vec,&n);
550: PetscObjectName((PetscObject)vec);
551: VecGetArrayRead(vec,&array);
552: PetscViewerMatlabPutArray(viewer,n,1,array,((PetscObject)vec)->name);
553: VecRestoreArrayRead(vec,&array);
554: return(0);
555: }
556: #endif
558: PETSC_EXTERN PetscErrorCode VecView_Seq(Vec xin,PetscViewer viewer)
559: {
561: PetscBool isdraw,iascii,issocket,isbinary;
562: #if defined(PETSC_HAVE_MATHEMATICA)
563: PetscBool ismathematica;
564: #endif
565: #if defined(PETSC_HAVE_MATLAB_ENGINE)
566: PetscBool ismatlab;
567: #endif
568: #if defined(PETSC_HAVE_HDF5)
569: PetscBool ishdf5;
570: #endif
571: PetscBool isglvis;
572: #if defined(PETSC_HAVE_ADIOS)
573: PetscBool isadios;
574: #endif
577: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
578: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
579: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSOCKET,&issocket);
580: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
581: #if defined(PETSC_HAVE_MATHEMATICA)
582: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERMATHEMATICA,&ismathematica);
583: #endif
584: #if defined(PETSC_HAVE_HDF5)
585: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
586: #endif
587: #if defined(PETSC_HAVE_MATLAB_ENGINE)
588: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERMATLAB,&ismatlab);
589: #endif
590: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERGLVIS,&isglvis);
591: #if defined(PETSC_HAVE_ADIOS)
592: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERADIOS,&isadios);
593: #endif
595: if (isdraw) {
596: VecView_Seq_Draw(xin,viewer);
597: } else if (iascii) {
598: VecView_Seq_ASCII(xin,viewer);
599: } else if (isbinary) {
600: VecView_Seq_Binary(xin,viewer);
601: #if defined(PETSC_HAVE_MATHEMATICA)
602: } else if (ismathematica) {
603: PetscViewerMathematicaPutVector(viewer,xin);
604: #endif
605: #if defined(PETSC_HAVE_HDF5)
606: } else if (ishdf5) {
607: VecView_MPI_HDF5(xin,viewer); /* Reusing VecView_MPI_HDF5 ... don't want code duplication*/
608: #endif
609: #if defined(PETSC_HAVE_ADIOS)
610: } else if (isadios) {
611: VecView_MPI_ADIOS(xin,viewer); /* Reusing VecView_MPI_ADIOS ... don't want code duplication*/
612: #endif
613: #if defined(PETSC_HAVE_MATLAB_ENGINE)
614: } else if (ismatlab) {
615: VecView_Seq_Matlab(xin,viewer);
616: #endif
617: } else if (isglvis) {
618: VecView_GLVis(xin,viewer);
619: }
620: return(0);
621: }
623: PetscErrorCode VecGetValues_Seq(Vec xin,PetscInt ni,const PetscInt ix[],PetscScalar y[])
624: {
625: const PetscScalar *xx;
626: PetscInt i;
627: PetscErrorCode ierr;
630: VecGetArrayRead(xin,&xx);
631: for (i=0; i<ni; i++) {
632: if (xin->stash.ignorenegidx && ix[i] < 0) continue;
633: if (PetscDefined(USE_DEBUG)) {
634: if (ix[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D cannot be negative",ix[i]);
635: if (ix[i] >= xin->map->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D to large maximum allowed %D",ix[i],xin->map->n);
636: }
637: y[i] = xx[ix[i]];
638: }
639: VecRestoreArrayRead(xin,&xx);
640: return(0);
641: }
643: PetscErrorCode VecSetValues_Seq(Vec xin,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode m)
644: {
645: PetscScalar *xx;
646: PetscInt i;
650: VecGetArray(xin,&xx);
651: if (m == INSERT_VALUES) {
652: for (i=0; i<ni; i++) {
653: if (xin->stash.ignorenegidx && ix[i] < 0) continue;
654: if (PetscDefined(USE_DEBUG)) {
655: if (ix[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D cannot be negative",ix[i]);
656: if (ix[i] >= xin->map->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D maximum %D",ix[i],xin->map->n);
657: }
658: xx[ix[i]] = y[i];
659: }
660: } else {
661: for (i=0; i<ni; i++) {
662: if (xin->stash.ignorenegidx && ix[i] < 0) continue;
663: if (PetscDefined(USE_DEBUG)) {
664: if (ix[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D cannot be negative",ix[i]);
665: if (ix[i] >= xin->map->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D maximum %D",ix[i],xin->map->n);
666: }
667: xx[ix[i]] += y[i];
668: }
669: }
670: VecRestoreArray(xin,&xx);
671: return(0);
672: }
674: PetscErrorCode VecSetValuesBlocked_Seq(Vec xin,PetscInt ni,const PetscInt ix[],const PetscScalar yin[],InsertMode m)
675: {
676: PetscScalar *xx,*y = (PetscScalar*)yin;
677: PetscInt i,bs,start,j;
680: /*
681: For optimization could treat bs = 2, 3, 4, 5 as special cases with loop unrolling
682: */
684: VecGetBlockSize(xin,&bs);
685: VecGetArray(xin,&xx);
686: if (m == INSERT_VALUES) {
687: for (i=0; i<ni; i++, y+=bs) {
688: start = bs*ix[i];
689: if (start < 0) continue;
690: if (PetscUnlikelyDebug(start >= xin->map->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D maximum %D",start,xin->map->n);
691: for (j=0; j<bs; j++) xx[start+j] = y[j];
692: }
693: } else {
694: for (i=0; i<ni; i++, y+=bs) {
695: start = bs*ix[i];
696: if (start < 0) continue;
697: if (PetscUnlikelyDebug(start >= xin->map->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D maximum %D",start,xin->map->n);
698: for (j=0; j<bs; j++) xx[start+j] += y[j];
699: }
700: }
701: VecRestoreArray(xin,&xx);
702: return(0);
703: }
705: PetscErrorCode VecDestroy_Seq(Vec v)
706: {
707: Vec_Seq *vs = (Vec_Seq*)v->data;
711: #if defined(PETSC_USE_LOG)
712: PetscLogObjectState((PetscObject)v,"Length=%D",v->map->n);
713: #endif
714: if (vs) { PetscFree(vs->array_allocated); }
715: PetscFree(v->data);
716: return(0);
717: }
719: PetscErrorCode VecSetOption_Seq(Vec v,VecOption op,PetscBool flag)
720: {
722: if (op == VEC_IGNORE_NEGATIVE_INDICES) v->stash.ignorenegidx = flag;
723: return(0);
724: }
726: PetscErrorCode VecDuplicate_Seq(Vec win,Vec *V)
727: {
731: VecCreate(PetscObjectComm((PetscObject)win),V);
732: VecSetSizes(*V,win->map->n,win->map->n);
733: VecSetType(*V,((PetscObject)win)->type_name);
734: PetscLayoutReference(win->map,&(*V)->map);
735: PetscObjectListDuplicate(((PetscObject)win)->olist,&((PetscObject)(*V))->olist);
736: PetscFunctionListDuplicate(((PetscObject)win)->qlist,&((PetscObject)(*V))->qlist);
738: (*V)->ops->view = win->ops->view;
739: (*V)->stash.ignorenegidx = win->stash.ignorenegidx;
740: return(0);
741: }
743: static struct _VecOps DvOps = {VecDuplicate_Seq, /* 1 */
744: VecDuplicateVecs_Default,
745: VecDestroyVecs_Default,
746: VecDot_Seq,
747: VecMDot_Seq,
748: VecNorm_Seq,
749: VecTDot_Seq,
750: VecMTDot_Seq,
751: VecScale_Seq,
752: VecCopy_Seq, /* 10 */
753: VecSet_Seq,
754: VecSwap_Seq,
755: VecAXPY_Seq,
756: VecAXPBY_Seq,
757: VecMAXPY_Seq,
758: VecAYPX_Seq,
759: VecWAXPY_Seq,
760: VecAXPBYPCZ_Seq,
761: VecPointwiseMult_Seq,
762: VecPointwiseDivide_Seq,
763: VecSetValues_Seq, /* 20 */
764: NULL,
765: NULL,
766: NULL,
767: VecGetSize_Seq,
768: VecGetSize_Seq,
769: NULL,
770: VecMax_Seq,
771: VecMin_Seq,
772: VecSetRandom_Seq,
773: VecSetOption_Seq, /* 30 */
774: VecSetValuesBlocked_Seq,
775: VecDestroy_Seq,
776: VecView_Seq,
777: VecPlaceArray_Seq,
778: VecReplaceArray_Seq,
779: VecDot_Seq,
780: VecTDot_Seq,
781: VecNorm_Seq,
782: VecMDot_Seq,
783: VecMTDot_Seq, /* 40 */
784: VecLoad_Default,
785: VecReciprocal_Default,
786: VecConjugate_Seq,
787: NULL,
788: NULL,
789: VecResetArray_Seq,
790: NULL,
791: VecMaxPointwiseDivide_Seq,
792: VecPointwiseMax_Seq,
793: VecPointwiseMaxAbs_Seq,
794: VecPointwiseMin_Seq,
795: VecGetValues_Seq,
796: NULL,
797: NULL,
798: NULL,
799: NULL,
800: NULL,
801: NULL,
802: VecStrideGather_Default,
803: VecStrideScatter_Default,
804: NULL,
805: NULL,
806: NULL,
807: NULL,
808: NULL,
809: VecStrideSubSetGather_Default,
810: VecStrideSubSetScatter_Default,
811: NULL,
812: NULL,
813: NULL
814: };
816: /*
817: This is called by VecCreate_Seq() (i.e. VecCreateSeq()) and VecCreateSeqWithArray()
818: */
819: PetscErrorCode VecCreate_Seq_Private(Vec v,const PetscScalar array[])
820: {
821: Vec_Seq *s;
825: PetscNewLog(v,&s);
826: PetscMemcpy(v->ops,&DvOps,sizeof(DvOps));
828: v->data = (void*)s;
829: v->petscnative = PETSC_TRUE;
830: s->array = (PetscScalar*)array;
831: s->array_allocated = NULL;
832: if (array) v->offloadmask = PETSC_OFFLOAD_CPU;
834: PetscLayoutSetUp(v->map);
835: PetscObjectChangeTypeName((PetscObject)v,VECSEQ);
836: #if defined(PETSC_HAVE_MATLAB_ENGINE)
837: PetscObjectComposeFunction((PetscObject)v,"PetscMatlabEnginePut_C",VecMatlabEnginePut_Default);
838: PetscObjectComposeFunction((PetscObject)v,"PetscMatlabEngineGet_C",VecMatlabEngineGet_Default);
839: #endif
840: return(0);
841: }
843: /*@C
844: VecCreateSeqWithArray - Creates a standard,sequential array-style vector,
845: where the user provides the array space to store the vector values.
847: Collective
849: Input Parameters:
850: + comm - the communicator, should be PETSC_COMM_SELF
851: . bs - the block size
852: . n - the vector length
853: - array - memory where the vector elements are to be stored.
855: Output Parameter:
856: . V - the vector
858: Notes:
859: Use VecDuplicate() or VecDuplicateVecs() to form additional vectors of the
860: same type as an existing vector.
862: If the user-provided array is NULL, then VecPlaceArray() can be used
863: at a later stage to SET the array for storing the vector values.
865: PETSc does NOT free the array when the vector is destroyed via VecDestroy().
866: The user should not free the array until the vector is destroyed.
868: Level: intermediate
870: .seealso: VecCreateMPIWithArray(), VecCreate(), VecDuplicate(), VecDuplicateVecs(),
871: VecCreateGhost(), VecCreateSeq(), VecPlaceArray()
872: @*/
873: PetscErrorCode VecCreateSeqWithArray(MPI_Comm comm,PetscInt bs,PetscInt n,const PetscScalar array[],Vec *V)
874: {
876: PetscMPIInt size;
879: VecCreate(comm,V);
880: VecSetSizes(*V,n,n);
881: VecSetBlockSize(*V,bs);
882: MPI_Comm_size(comm,&size);
883: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Cannot create VECSEQ on more than one process");
884: VecCreate_Seq_Private(*V,array);
885: return(0);
886: }