Actual source code: bvec2.c
petsc-3.7.5 2017-01-01
2: /*
3: Implements the sequential vectors.
4: */
6: #include <../src/vec/vec/impls/dvecimpl.h> /*I "petscvec.h" I*/
7: #include <../src/vec/vec/impls/mpi/pvecimpl.h> /* For VecView_MPI_HDF5 */
8: #include <petscblaslapack.h>
10: #if defined(PETSC_HAVE_HDF5)
11: extern PetscErrorCode VecView_MPI_HDF5(Vec,PetscViewer);
12: #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: }
38: PetscErrorCode VecPointwiseMin_Seq(Vec win,Vec xin,Vec yin)
39: {
41: PetscInt n = win->map->n,i;
42: PetscScalar *ww,*xx,*yy; /* cannot make xx or yy const since might be ww */
45: VecGetArrayRead(xin,(const PetscScalar**)&xx);
46: VecGetArrayRead(yin,(const PetscScalar**)&yy);
47: VecGetArray(win,&ww);
49: for (i=0; i<n; i++) ww[i] = PetscMin(PetscRealPart(xx[i]),PetscRealPart(yy[i]));
51: VecRestoreArrayRead(xin,(const PetscScalar**)&xx);
52: VecRestoreArrayRead(yin,(const PetscScalar**)&yy);
53: VecRestoreArray(win,&ww);
54: PetscLogFlops(n);
55: return(0);
56: }
60: PetscErrorCode VecPointwiseMaxAbs_Seq(Vec win,Vec xin,Vec yin)
61: {
63: PetscInt n = win->map->n,i;
64: PetscScalar *ww,*xx,*yy; /* cannot make xx or yy const since might be ww */
67: VecGetArrayRead(xin,(const PetscScalar**)&xx);
68: VecGetArrayRead(yin,(const PetscScalar**)&yy);
69: VecGetArray(win,&ww);
71: for (i=0; i<n; i++) ww[i] = PetscMax(PetscAbsScalar(xx[i]),PetscAbsScalar(yy[i]));
73: PetscLogFlops(n);
74: VecRestoreArrayRead(xin,(const PetscScalar**)&xx);
75: VecRestoreArrayRead(yin,(const PetscScalar**)&yy);
76: VecRestoreArray(win,&ww);
77: return(0);
78: }
80: #include <../src/vec/vec/impls/seq/ftn-kernels/fxtimesy.h>
84: PetscErrorCode VecPointwiseMult_Seq(Vec win,Vec xin,Vec yin)
85: {
87: PetscInt n = win->map->n,i;
88: PetscScalar *ww,*xx,*yy; /* cannot make xx or yy const since might be ww */
91: VecGetArrayRead(xin,(const PetscScalar**)&xx);
92: VecGetArrayRead(yin,(const PetscScalar**)&yy);
93: VecGetArray(win,&ww);
94: if (ww == xx) {
95: for (i=0; i<n; i++) ww[i] *= yy[i];
96: } else if (ww == yy) {
97: for (i=0; i<n; i++) ww[i] *= xx[i];
98: } else {
99: #if defined(PETSC_USE_FORTRAN_KERNEL_XTIMESY)
100: fortranxtimesy_(xx,yy,ww,&n);
101: #else
102: for (i=0; i<n; i++) ww[i] = xx[i] * yy[i];
103: #endif
104: }
105: VecRestoreArrayRead(xin,(const PetscScalar**)&xx);
106: VecRestoreArrayRead(yin,(const PetscScalar**)&yy);
107: VecRestoreArray(win,&ww);
108: PetscLogFlops(n);
109: return(0);
110: }
114: PetscErrorCode VecPointwiseDivide_Seq(Vec win,Vec xin,Vec yin)
115: {
117: PetscInt n = win->map->n,i;
118: PetscScalar *ww,*xx,*yy; /* cannot make xx or yy const since might be ww */
121: VecGetArrayRead(xin,(const PetscScalar**)&xx);
122: VecGetArrayRead(yin,(const PetscScalar**)&yy);
123: VecGetArray(win,&ww);
125: for (i=0; i<n; i++) ww[i] = xx[i] / yy[i];
127: PetscLogFlops(n);
128: VecRestoreArrayRead(xin,(const PetscScalar**)&xx);
129: VecRestoreArrayRead(yin,(const PetscScalar**)&yy);
130: VecRestoreArray(win,&ww);
131: return(0);
132: }
136: PetscErrorCode VecSetRandom_Seq(Vec xin,PetscRandom r)
137: {
139: PetscInt n = xin->map->n,i;
140: PetscScalar *xx;
143: VecGetArray(xin,&xx);
144: for (i=0; i<n; i++) {PetscRandomGetValue(r,&xx[i]);}
145: VecRestoreArray(xin,&xx);
146: return(0);
147: }
151: PetscErrorCode VecGetSize_Seq(Vec vin,PetscInt *size)
152: {
154: *size = vin->map->n;
155: return(0);
156: }
162: PetscErrorCode VecConjugate_Seq(Vec xin)
163: {
164: PetscScalar *x;
165: PetscInt n = xin->map->n;
169: VecGetArray(xin,&x);
170: while (n-->0) {
171: *x = PetscConj(*x);
172: x++;
173: }
174: VecRestoreArray(xin,&x);
175: return(0);
176: }
180: PetscErrorCode VecResetArray_Seq(Vec vin)
181: {
182: Vec_Seq *v = (Vec_Seq*)vin->data;
185: v->array = v->unplacedarray;
186: v->unplacedarray = 0;
187: return(0);
188: }
192: PetscErrorCode VecCopy_Seq(Vec xin,Vec yin)
193: {
194: PetscScalar *ya;
195: const PetscScalar *xa;
196: PetscErrorCode ierr;
199: if (xin != yin) {
200: VecGetArrayRead(xin,&xa);
201: VecGetArray(yin,&ya);
202: PetscMemcpy(ya,xa,xin->map->n*sizeof(PetscScalar));
203: VecRestoreArrayRead(xin,&xa);
204: VecRestoreArray(yin,&ya);
205: }
206: return(0);
207: }
211: PetscErrorCode VecSwap_Seq(Vec xin,Vec yin)
212: {
213: PetscScalar *ya, *xa;
215: PetscBLASInt one = 1,bn;
218: if (xin != yin) {
219: PetscBLASIntCast(xin->map->n,&bn);
220: VecGetArray(xin,&xa);
221: VecGetArray(yin,&ya);
222: PetscStackCallBLAS("BLASswap",BLASswap_(&bn,xa,&one,ya,&one));
223: VecRestoreArray(xin,&xa);
224: VecRestoreArray(yin,&ya);
225: }
226: return(0);
227: }
229: #include <../src/vec/vec/impls/seq/ftn-kernels/fnorm.h>
233: PetscErrorCode VecNorm_Seq(Vec xin,NormType type,PetscReal *z)
234: {
235: const PetscScalar *xx;
236: PetscErrorCode ierr;
237: PetscInt n = xin->map->n;
238: PetscBLASInt one = 1, bn;
241: PetscBLASIntCast(n,&bn);
242: if (type == NORM_2 || type == NORM_FROBENIUS) {
243: VecGetArrayRead(xin,&xx);
244: *z = PetscRealPart(BLASdot_(&bn,xx,&one,xx,&one));
245: *z = PetscSqrtReal(*z);
246: VecRestoreArrayRead(xin,&xx);
247: PetscLogFlops(PetscMax(2.0*n-1,0.0));
248: } else if (type == NORM_INFINITY) {
249: PetscInt i;
250: PetscReal max = 0.0,tmp;
252: VecGetArrayRead(xin,&xx);
253: for (i=0; i<n; i++) {
254: if ((tmp = PetscAbsScalar(*xx)) > max) max = tmp;
255: /* check special case of tmp == NaN */
256: if (tmp != tmp) {max = tmp; break;}
257: xx++;
258: }
259: VecRestoreArrayRead(xin,&xx);
260: *z = max;
261: } else if (type == NORM_1) {
262: #if defined(PETSC_USE_COMPLEX)
263: PetscReal tmp = 0.0;
264: PetscInt i;
265: #endif
266: VecGetArrayRead(xin,&xx);
267: #if defined(PETSC_USE_COMPLEX)
268: /* BLASasum() returns the nonstandard 1 norm of the 1 norm of the complex entries so we provide a custom loop instead */
269: for (i=0; i<n; i++) {
270: tmp += PetscAbsScalar(xx[i]);
271: }
272: *z = tmp;
273: #else
274: PetscStackCallBLAS("BLASasum",*z = BLASasum_(&bn,xx,&one));
275: #endif
276: VecRestoreArrayRead(xin,&xx);
277: PetscLogFlops(PetscMax(n-1.0,0.0));
278: } else if (type == NORM_1_AND_2) {
279: VecNorm_Seq(xin,NORM_1,z);
280: VecNorm_Seq(xin,NORM_2,z+1);
281: }
282: return(0);
283: }
287: PetscErrorCode VecView_Seq_ASCII(Vec xin,PetscViewer viewer)
288: {
289: PetscErrorCode ierr;
290: PetscInt i,n = xin->map->n;
291: const char *name;
292: PetscViewerFormat format;
293: const PetscScalar *xv;
296: VecGetArrayRead(xin,&xv);
297: PetscViewerGetFormat(viewer,&format);
298: if (format == PETSC_VIEWER_ASCII_MATLAB) {
299: PetscObjectGetName((PetscObject)xin,&name);
300: PetscViewerASCIIPrintf(viewer,"%s = [\n",name);
301: for (i=0; i<n; i++) {
302: #if defined(PETSC_USE_COMPLEX)
303: if (PetscImaginaryPart(xv[i]) > 0.0) {
304: PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei\n",(double)PetscRealPart(xv[i]),(double)PetscImaginaryPart(xv[i]));
305: } else if (PetscImaginaryPart(xv[i]) < 0.0) {
306: PetscViewerASCIIPrintf(viewer,"%18.16e - %18.16ei\n",(double)PetscRealPart(xv[i]),-(double)PetscImaginaryPart(xv[i]));
307: } else {
308: PetscViewerASCIIPrintf(viewer,"%18.16e\n",(double)PetscRealPart(xv[i]));
309: }
310: #else
311: PetscViewerASCIIPrintf(viewer,"%18.16e\n",(double) xv[i]);
312: #endif
313: }
314: PetscViewerASCIIPrintf(viewer,"];\n");
315: } else if (format == PETSC_VIEWER_ASCII_SYMMODU) {
316: for (i=0; i<n; i++) {
317: #if defined(PETSC_USE_COMPLEX)
318: PetscViewerASCIIPrintf(viewer,"%18.16e %18.16e\n",(double)PetscRealPart(xv[i]),(double)PetscImaginaryPart(xv[i]));
319: #else
320: PetscViewerASCIIPrintf(viewer,"%18.16e\n",(double)xv[i]);
321: #endif
322: }
323: } else if (format == PETSC_VIEWER_ASCII_VTK || format == PETSC_VIEWER_ASCII_VTK_CELL) {
324: /*
325: state 0: No header has been output
326: state 1: Only POINT_DATA has been output
327: state 2: Only CELL_DATA has been output
328: state 3: Output both, POINT_DATA last
329: state 4: Output both, CELL_DATA last
330: */
331: static PetscInt stateId = -1;
332: int outputState = 0;
333: PetscBool hasState;
334: int doOutput = 0;
335: PetscInt bs, b;
337: if (stateId < 0) {
338: PetscObjectComposedDataRegister(&stateId);
339: }
340: PetscObjectComposedDataGetInt((PetscObject) viewer, stateId, outputState, hasState);
341: if (!hasState) outputState = 0;
342: PetscObjectGetName((PetscObject) xin, &name);
343: VecGetBlockSize(xin, &bs);
344: 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);
345: if (format == PETSC_VIEWER_ASCII_VTK) {
346: if (outputState == 0) {
347: outputState = 1;
348: doOutput = 1;
349: } else if (outputState == 1) doOutput = 0;
350: else if (outputState == 2) {
351: outputState = 3;
352: doOutput = 1;
353: } else if (outputState == 3) doOutput = 0;
354: else if (outputState == 4) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Tried to output POINT_DATA again after intervening CELL_DATA");
356: if (doOutput) {
357: PetscViewerASCIIPrintf(viewer, "POINT_DATA %d\n", n/bs);
358: }
359: } else {
360: if (outputState == 0) {
361: outputState = 2;
362: doOutput = 1;
363: } else if (outputState == 1) {
364: outputState = 4;
365: doOutput = 1;
366: } else if (outputState == 2) doOutput = 0;
367: else if (outputState == 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Tried to output CELL_DATA again after intervening POINT_DATA");
368: else if (outputState == 4) doOutput = 0;
370: if (doOutput) {
371: PetscViewerASCIIPrintf(viewer, "CELL_DATA %d\n", n);
372: }
373: }
374: PetscObjectComposedDataSetInt((PetscObject) viewer, stateId, outputState);
375: if (name) {
376: if (bs == 3) {
377: PetscViewerASCIIPrintf(viewer, "VECTORS %s double\n", name);
378: } else {
379: PetscViewerASCIIPrintf(viewer, "SCALARS %s double %d\n", name, bs);
380: }
381: } else {
382: PetscViewerASCIIPrintf(viewer, "SCALARS scalars double %d\n", bs);
383: }
384: if (bs != 3) {
385: PetscViewerASCIIPrintf(viewer, "LOOKUP_TABLE default\n");
386: }
387: for (i=0; i<n/bs; i++) {
388: for (b=0; b<bs; b++) {
389: if (b > 0) {
390: PetscViewerASCIIPrintf(viewer," ");
391: }
392: #if !defined(PETSC_USE_COMPLEX)
393: PetscViewerASCIIPrintf(viewer,"%g",(double)xv[i*bs+b]);
394: #endif
395: }
396: PetscViewerASCIIPrintf(viewer,"\n");
397: }
398: } else if (format == PETSC_VIEWER_ASCII_VTK_COORDS) {
399: PetscInt bs, b;
401: VecGetBlockSize(xin, &bs);
402: 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);
403: for (i=0; i<n/bs; i++) {
404: for (b=0; b<bs; b++) {
405: if (b > 0) {
406: PetscViewerASCIIPrintf(viewer," ");
407: }
408: #if !defined(PETSC_USE_COMPLEX)
409: PetscViewerASCIIPrintf(viewer,"%g",(double)xv[i*bs+b]);
410: #endif
411: }
412: for (b=bs; b<3; b++) {
413: PetscViewerASCIIPrintf(viewer," 0.0");
414: }
415: PetscViewerASCIIPrintf(viewer,"\n");
416: }
417: } else if (format == PETSC_VIEWER_ASCII_PCICE) {
418: PetscInt bs, b;
420: VecGetBlockSize(xin, &bs);
421: 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);
422: PetscViewerASCIIPrintf(viewer,"%D\n", xin->map->N/bs);
423: for (i=0; i<n/bs; i++) {
424: PetscViewerASCIIPrintf(viewer,"%7D ", i+1);
425: for (b=0; b<bs; b++) {
426: if (b > 0) {
427: PetscViewerASCIIPrintf(viewer," ");
428: }
429: #if !defined(PETSC_USE_COMPLEX)
430: PetscViewerASCIIPrintf(viewer,"% 12.5E",(double)xv[i*bs+b]);
431: #endif
432: }
433: PetscViewerASCIIPrintf(viewer,"\n");
434: }
435: } else if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
436: /* No info */
437: } else {
438: for (i=0; i<n; i++) {
439: if (format == PETSC_VIEWER_ASCII_INDEX) {
440: PetscViewerASCIIPrintf(viewer,"%D: ",i);
441: }
442: #if defined(PETSC_USE_COMPLEX)
443: if (PetscImaginaryPart(xv[i]) > 0.0) {
444: PetscViewerASCIIPrintf(viewer,"%g + %g i\n",(double)PetscRealPart(xv[i]),(double)PetscImaginaryPart(xv[i]));
445: } else if (PetscImaginaryPart(xv[i]) < 0.0) {
446: PetscViewerASCIIPrintf(viewer,"%g - %g i\n",(double)PetscRealPart(xv[i]),-(double)PetscImaginaryPart(xv[i]));
447: } else {
448: PetscViewerASCIIPrintf(viewer,"%g\n",(double)PetscRealPart(xv[i]));
449: }
450: #else
451: PetscViewerASCIIPrintf(viewer,"%g\n",(double)xv[i]);
452: #endif
453: }
454: }
455: PetscViewerFlush(viewer);
456: VecRestoreArrayRead(xin,&xv);
457: return(0);
458: }
460: #include <petscdraw.h>
463: PetscErrorCode VecView_Seq_Draw_LG(Vec xin,PetscViewer v)
464: {
465: PetscDraw draw;
466: PetscBool isnull;
467: PetscDrawLG lg;
468: PetscErrorCode ierr;
469: PetscInt i,c,bs = PetscAbs(xin->map->bs),n = xin->map->n/bs;
470: const PetscScalar *xv;
471: PetscReal *xx,*yy;
472: int colors[] = {PETSC_DRAW_RED};
475: PetscViewerDrawGetDraw(v,0,&draw);
476: PetscDrawIsNull(draw,&isnull);
477: if (isnull) return(0);
479: PetscMalloc2(n,&xx,n,&yy);
480: VecGetArrayRead(xin,&xv);
481: for (c=0; c<bs; c++) {
482: PetscViewerDrawGetDrawLG(v,c,&lg);
483: PetscDrawLGReset(lg);
484: PetscDrawLGSetDimension(lg,1);
485: PetscDrawLGSetColors(lg,colors);
486: for (i=0; i<n; i++) {
487: xx[i] = (PetscReal)i;
488: yy[i] = PetscRealPart(xv[c + i*bs]);
489: }
490: PetscDrawLGAddPoints(lg,n,&xx,&yy);
491: PetscDrawLGDraw(lg);
492: PetscDrawLGSave(lg);
493: }
494: VecRestoreArrayRead(xin,&xv);
495: PetscFree2(xx,yy);
496: return(0);
497: }
501: PetscErrorCode VecView_Seq_Draw(Vec xin,PetscViewer v)
502: {
503: PetscErrorCode ierr;
504: PetscDraw draw;
505: PetscBool isnull;
508: PetscViewerDrawGetDraw(v,0,&draw);
509: PetscDrawIsNull(draw,&isnull);
510: if (isnull) return(0);
511: PetscViewerPushFormat(v,PETSC_VIEWER_DRAW_LG);
512: VecView_Seq_Draw_LG(xin,v);
513: PetscViewerPopFormat(v);
514: return(0);
515: }
519: PetscErrorCode VecView_Seq_Binary(Vec xin,PetscViewer viewer)
520: {
521: PetscErrorCode ierr;
522: int fdes;
523: PetscInt n = xin->map->n,classid=VEC_FILE_CLASSID;
524: FILE *file;
525: const PetscScalar *xv;
526: #if defined(PETSC_HAVE_MPIIO)
527: PetscBool isMPIIO;
528: #endif
529: PetscBool skipHeader;
530: PetscViewerFormat format;
533: /* Write vector header */
534: PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
535: if (!skipHeader) {
536: PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);
537: PetscViewerBinaryWrite(viewer,&n,1,PETSC_INT,PETSC_FALSE);
538: }
540: /* Write vector contents */
541: #if defined(PETSC_HAVE_MPIIO)
542: PetscViewerBinaryGetUseMPIIO(viewer,&isMPIIO);
543: if (!isMPIIO) {
544: #endif
545: PetscViewerBinaryGetDescriptor(viewer,&fdes);
546: VecGetArrayRead(xin,&xv);
547: PetscBinaryWrite(fdes,(void*)xv,n,PETSC_SCALAR,PETSC_FALSE);
548: VecRestoreArrayRead(xin,&xv);
549: PetscViewerGetFormat(viewer,&format);
550: if (format == PETSC_VIEWER_BINARY_MATLAB) {
551: MPI_Comm comm;
552: FILE *info;
553: const char *name;
555: PetscObjectGetName((PetscObject)xin,&name);
556: PetscObjectGetComm((PetscObject)viewer,&comm);
557: PetscViewerBinaryGetInfoPointer(viewer,&info);
558: PetscFPrintf(comm,info,"#--- begin code written by PetscViewerBinary for MATLAB format ---#\n");
559: PetscFPrintf(comm,info,"#$$ Set.%s = PetscBinaryRead(fd);\n",name);
560: PetscFPrintf(comm,info,"#--- end code written by PetscViewerBinary for MATLAB format ---#\n\n");
561: }
562: #if defined(PETSC_HAVE_MPIIO)
563: } else {
564: MPI_Offset off;
565: MPI_File mfdes;
566: PetscMPIInt lsize;
568: PetscMPIIntCast(n,&lsize);
569: PetscViewerBinaryGetMPIIODescriptor(viewer,&mfdes);
570: PetscViewerBinaryGetMPIIOOffset(viewer,&off);
571: MPI_File_set_view(mfdes,off,MPIU_SCALAR,MPIU_SCALAR,(char*)"native",MPI_INFO_NULL);
572: VecGetArrayRead(xin,&xv);
573: MPIU_File_write_all(mfdes,(void*)xv,lsize,MPIU_SCALAR,MPI_STATUS_IGNORE);
574: VecRestoreArrayRead(xin,&xv);
575: PetscViewerBinaryAddMPIIOOffset(viewer,n*sizeof(PetscScalar));
576: }
577: #endif
579: PetscViewerBinaryGetInfoPointer(viewer,&file);
580: if (file) {
581: if (((PetscObject)xin)->prefix) {
582: PetscFPrintf(PETSC_COMM_SELF,file,"-%svecload_block_size %D\n",((PetscObject)xin)->prefix,PetscAbs(xin->map->bs));
583: } else {
584: PetscFPrintf(PETSC_COMM_SELF,file,"-vecload_block_size %D\n",PetscAbs(xin->map->bs));
585: }
586: }
587: return(0);
588: }
590: #if defined(PETSC_HAVE_MATLAB_ENGINE)
591: #include <petscmatlab.h>
592: #include <mat.h> /* MATLAB include file */
595: PetscErrorCode VecView_Seq_Matlab(Vec vec,PetscViewer viewer)
596: {
597: PetscErrorCode ierr;
598: PetscInt n;
599: const PetscScalar *array;
602: VecGetLocalSize(vec,&n);
603: PetscObjectName((PetscObject)vec);
604: VecGetArrayRead(vec,&array);
605: PetscViewerMatlabPutArray(viewer,n,1,array,((PetscObject)vec)->name);
606: VecRestoreArrayRead(vec,&array);
607: return(0);
608: }
609: #endif
613: PETSC_EXTERN PetscErrorCode VecView_Seq(Vec xin,PetscViewer viewer)
614: {
616: PetscBool isdraw,iascii,issocket,isbinary;
617: #if defined(PETSC_HAVE_MATHEMATICA)
618: PetscBool ismathematica;
619: #endif
620: #if defined(PETSC_HAVE_MATLAB_ENGINE)
621: PetscBool ismatlab;
622: #endif
623: #if defined(PETSC_HAVE_HDF5)
624: PetscBool ishdf5;
625: #endif
628: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
629: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
630: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSOCKET,&issocket);
631: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
632: #if defined(PETSC_HAVE_MATHEMATICA)
633: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERMATHEMATICA,&ismathematica);
634: #endif
635: #if defined(PETSC_HAVE_HDF5)
636: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
637: #endif
638: #if defined(PETSC_HAVE_MATLAB_ENGINE)
639: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERMATLAB,&ismatlab);
640: #endif
642: if (isdraw) {
643: VecView_Seq_Draw(xin,viewer);
644: } else if (iascii) {
645: VecView_Seq_ASCII(xin,viewer);
646: } else if (isbinary) {
647: VecView_Seq_Binary(xin,viewer);
648: #if defined(PETSC_HAVE_MATHEMATICA)
649: } else if (ismathematica) {
650: PetscViewerMathematicaPutVector(viewer,xin);
651: #endif
652: #if defined(PETSC_HAVE_HDF5)
653: } else if (ishdf5) {
654: VecView_MPI_HDF5(xin,viewer); /* Reusing VecView_MPI_HDF5 ... don't want code duplication*/
655: #endif
656: #if defined(PETSC_HAVE_MATLAB_ENGINE)
657: } else if (ismatlab) {
658: VecView_Seq_Matlab(xin,viewer);
659: #endif
660: }
661: return(0);
662: }
666: PetscErrorCode VecGetValues_Seq(Vec xin,PetscInt ni,const PetscInt ix[],PetscScalar y[])
667: {
668: const PetscScalar *xx;
669: PetscInt i;
670: PetscErrorCode ierr;
673: VecGetArrayRead(xin,&xx);
674: for (i=0; i<ni; i++) {
675: if (xin->stash.ignorenegidx && ix[i] < 0) continue;
676: #if defined(PETSC_USE_DEBUG)
677: if (ix[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D cannot be negative",ix[i]);
678: 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);
679: #endif
680: y[i] = xx[ix[i]];
681: }
682: VecRestoreArrayRead(xin,&xx);
683: return(0);
684: }
688: PetscErrorCode VecSetValues_Seq(Vec xin,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode m)
689: {
690: PetscScalar *xx;
691: PetscInt i;
695: VecGetArray(xin,&xx);
696: if (m == INSERT_VALUES) {
697: for (i=0; i<ni; i++) {
698: if (xin->stash.ignorenegidx && ix[i] < 0) continue;
699: #if defined(PETSC_USE_DEBUG)
700: if (ix[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D cannot be negative",ix[i]);
701: 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);
702: #endif
703: xx[ix[i]] = y[i];
704: }
705: } else {
706: for (i=0; i<ni; i++) {
707: if (xin->stash.ignorenegidx && ix[i] < 0) continue;
708: #if defined(PETSC_USE_DEBUG)
709: if (ix[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D cannot be negative",ix[i]);
710: 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);
711: #endif
712: xx[ix[i]] += y[i];
713: }
714: }
715: VecRestoreArray(xin,&xx);
716: return(0);
717: }
721: PetscErrorCode VecSetValuesBlocked_Seq(Vec xin,PetscInt ni,const PetscInt ix[],const PetscScalar yin[],InsertMode m)
722: {
723: PetscScalar *xx,*y = (PetscScalar*)yin;
724: PetscInt i,bs,start,j;
727: /*
728: For optimization could treat bs = 2, 3, 4, 5 as special cases with loop unrolling
729: */
731: VecGetBlockSize(xin,&bs);
732: VecGetArray(xin,&xx);
733: if (m == INSERT_VALUES) {
734: for (i=0; i<ni; i++) {
735: start = bs*ix[i];
736: if (start < 0) continue;
737: #if defined(PETSC_USE_DEBUG)
738: if (start >= xin->map->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D maximum %D",start,xin->map->n);
739: #endif
740: for (j=0; j<bs; j++) xx[start+j] = y[j];
741: y += bs;
742: }
743: } else {
744: for (i=0; i<ni; i++) {
745: start = bs*ix[i];
746: if (start < 0) continue;
747: #if defined(PETSC_USE_DEBUG)
748: if (start >= xin->map->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D maximum %D",start,xin->map->n);
749: #endif
750: for (j=0; j<bs; j++) xx[start+j] += y[j];
751: y += bs;
752: }
753: }
754: VecRestoreArray(xin,&xx);
755: return(0);
756: }
760: PetscErrorCode VecDestroy_Seq(Vec v)
761: {
762: Vec_Seq *vs = (Vec_Seq*)v->data;
766: #if defined(PETSC_USE_LOG)
767: PetscLogObjectState((PetscObject)v,"Length=%D",v->map->n);
768: #endif
769: PetscFree(vs->array_allocated);
770: PetscFree(v->data);
771: return(0);
772: }
776: PetscErrorCode VecSetOption_Seq(Vec v,VecOption op,PetscBool flag)
777: {
779: if (op == VEC_IGNORE_NEGATIVE_INDICES) v->stash.ignorenegidx = flag;
780: return(0);
781: }
785: PetscErrorCode VecDuplicate_Seq(Vec win,Vec *V)
786: {
790: VecCreate(PetscObjectComm((PetscObject)win),V);
791: PetscObjectSetPrecision((PetscObject)*V,((PetscObject)win)->precision);
792: VecSetSizes(*V,win->map->n,win->map->n);
793: VecSetType(*V,((PetscObject)win)->type_name);
794: PetscLayoutReference(win->map,&(*V)->map);
795: PetscObjectListDuplicate(((PetscObject)win)->olist,&((PetscObject)(*V))->olist);
796: PetscFunctionListDuplicate(((PetscObject)win)->qlist,&((PetscObject)(*V))->qlist);
798: (*V)->ops->view = win->ops->view;
799: (*V)->stash.ignorenegidx = win->stash.ignorenegidx;
800: return(0);
801: }
803: static struct _VecOps DvOps = {VecDuplicate_Seq, /* 1 */
804: VecDuplicateVecs_Default,
805: VecDestroyVecs_Default,
806: VecDot_Seq,
807: VecMDot_Seq,
808: VecNorm_Seq,
809: VecTDot_Seq,
810: VecMTDot_Seq,
811: VecScale_Seq,
812: VecCopy_Seq, /* 10 */
813: VecSet_Seq,
814: VecSwap_Seq,
815: VecAXPY_Seq,
816: VecAXPBY_Seq,
817: VecMAXPY_Seq,
818: VecAYPX_Seq,
819: VecWAXPY_Seq,
820: VecAXPBYPCZ_Seq,
821: VecPointwiseMult_Seq,
822: VecPointwiseDivide_Seq,
823: VecSetValues_Seq, /* 20 */
824: 0,0,
825: 0,
826: VecGetSize_Seq,
827: VecGetSize_Seq,
828: 0,
829: VecMax_Seq,
830: VecMin_Seq,
831: VecSetRandom_Seq,
832: VecSetOption_Seq, /* 30 */
833: VecSetValuesBlocked_Seq,
834: VecDestroy_Seq,
835: VecView_Seq,
836: VecPlaceArray_Seq,
837: VecReplaceArray_Seq,
838: VecDot_Seq,
839: VecTDot_Seq,
840: VecNorm_Seq,
841: VecMDot_Seq,
842: VecMTDot_Seq, /* 40 */
843: VecLoad_Default,
844: VecReciprocal_Default,
845: VecConjugate_Seq,
846: 0,
847: 0,
848: VecResetArray_Seq,
849: 0,
850: VecMaxPointwiseDivide_Seq,
851: VecPointwiseMax_Seq,
852: VecPointwiseMaxAbs_Seq,
853: VecPointwiseMin_Seq,
854: VecGetValues_Seq,
855: 0,
856: 0,
857: 0,
858: 0,
859: 0,
860: 0,
861: VecStrideGather_Default,
862: VecStrideScatter_Default,
863: 0,
864: 0,
865: 0,
866: 0,
867: 0,
868: VecStrideSubSetGather_Default,
869: VecStrideSubSetScatter_Default,
870: 0,
871: 0
872: };
875: /*
876: This is called by VecCreate_Seq() (i.e. VecCreateSeq()) and VecCreateSeqWithArray()
877: */
880: PetscErrorCode VecCreate_Seq_Private(Vec v,const PetscScalar array[])
881: {
882: Vec_Seq *s;
886: PetscNewLog(v,&s);
887: PetscMemcpy(v->ops,&DvOps,sizeof(DvOps));
889: v->data = (void*)s;
890: v->petscnative = PETSC_TRUE;
891: s->array = (PetscScalar*)array;
892: s->array_allocated = 0;
894: PetscLayoutSetUp(v->map);
895: PetscObjectChangeTypeName((PetscObject)v,VECSEQ);
896: #if defined(PETSC_HAVE_MATLAB_ENGINE)
897: PetscObjectComposeFunction((PetscObject)v,"PetscMatlabEnginePut_C",VecMatlabEnginePut_Default);
898: PetscObjectComposeFunction((PetscObject)v,"PetscMatlabEngineGet_C",VecMatlabEngineGet_Default);
899: #endif
900: #if defined(PETSC_USE_MIXED_PRECISION)
901: ((PetscObject)v)->precision = (PetscPrecision)sizeof(PetscReal);
902: #endif
903: return(0);
904: }
908: /*@C
909: VecCreateSeqWithArray - Creates a standard,sequential array-style vector,
910: where the user provides the array space to store the vector values.
912: Collective on MPI_Comm
914: Input Parameter:
915: + comm - the communicator, should be PETSC_COMM_SELF
916: . bs - the block size
917: . n - the vector length
918: - array - memory where the vector elements are to be stored.
920: Output Parameter:
921: . V - the vector
923: Notes:
924: Use VecDuplicate() or VecDuplicateVecs() to form additional vectors of the
925: same type as an existing vector.
927: If the user-provided array is NULL, then VecPlaceArray() can be used
928: at a later stage to SET the array for storing the vector values.
930: PETSc does NOT free the array when the vector is destroyed via VecDestroy().
931: The user should not free the array until the vector is destroyed.
933: Level: intermediate
935: Concepts: vectors^creating with array
937: .seealso: VecCreateMPIWithArray(), VecCreate(), VecDuplicate(), VecDuplicateVecs(),
938: VecCreateGhost(), VecCreateSeq(), VecPlaceArray()
939: @*/
940: PetscErrorCode VecCreateSeqWithArray(MPI_Comm comm,PetscInt bs,PetscInt n,const PetscScalar array[],Vec *V)
941: {
943: PetscMPIInt size;
946: VecCreate(comm,V);
947: VecSetSizes(*V,n,n);
948: VecSetBlockSize(*V,bs);
949: MPI_Comm_size(comm,&size);
950: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Cannot create VECSEQ on more than one process");
951: VecCreate_Seq_Private(*V,array);
952: return(0);
953: }