Actual source code: bvec2.c
petsc-3.3-p7 2013-05-11
2: /*
3: Implements the sequential vectors.
4: */
6: #include <petsc-private/vecimpl.h> /*I "petscvec.h" I*/
7: #include <../src/vec/vec/impls/dvecimpl.h>
8: #include <../src/vec/vec/impls/mpi/pvecimpl.h> /* For VecView_MPI_HDF5 */
9: #include <petscblaslapack.h>
11: #if defined(PETSC_HAVE_HDF5)
12: extern PetscErrorCode VecView_MPI_HDF5(Vec,PetscViewer);
13: #endif
17: PetscErrorCode VecPointwiseMax_Seq(Vec win,Vec xin,Vec yin)
18: {
20: PetscInt n = win->map->n,i;
21: PetscScalar *ww,*xx,*yy; /* cannot make xx or yy const since might be ww */
24: VecGetArrayRead(xin,(const PetscScalar**)&xx);
25: VecGetArrayRead(yin,(const PetscScalar**)&yy);
26: VecGetArray(win,&ww);
27: for (i=0; i<n; i++) {
28: ww[i] = PetscMax(PetscRealPart(xx[i]),PetscRealPart(yy[i]));
29: }
30: VecRestoreArrayRead(xin,(const PetscScalar**)&xx);
31: VecRestoreArrayRead(yin,(const PetscScalar**)&yy);
32: VecRestoreArray(win,&ww);
33: PetscLogFlops(n);
34: return(0);
35: }
39: PetscErrorCode VecPointwiseMin_Seq(Vec win,Vec xin,Vec yin)
40: {
42: PetscInt n = win->map->n,i;
43: PetscScalar *ww,*xx,*yy; /* cannot make xx or yy const since might be ww */
46: VecGetArrayRead(xin,(const PetscScalar**)&xx);
47: VecGetArrayRead(yin,(const PetscScalar**)&yy);
48: VecGetArray(win,&ww);
49: for (i=0; i<n; i++) {
50: ww[i] = PetscMin(PetscRealPart(xx[i]),PetscRealPart(yy[i]));
51: }
52: VecRestoreArrayRead(xin,(const PetscScalar**)&xx);
53: VecRestoreArrayRead(yin,(const PetscScalar**)&yy);
54: VecRestoreArray(win,&ww);
55: PetscLogFlops(n);
56: return(0);
57: }
61: PetscErrorCode VecPointwiseMaxAbs_Seq(Vec win,Vec xin,Vec yin)
62: {
64: PetscInt n = win->map->n,i;
65: PetscScalar *ww,*xx,*yy; /* cannot make xx or yy const since might be ww */
68: VecGetArrayRead(xin,(const PetscScalar**)&xx);
69: VecGetArrayRead(yin,(const PetscScalar**)&yy);
70: VecGetArray(win,&ww);
71: for (i=0; i<n; i++) {
72: ww[i] = PetscMax(PetscAbsScalar(xx[i]),PetscAbsScalar(yy[i]));
73: }
74: PetscLogFlops(n);
75: VecRestoreArrayRead(xin,(const PetscScalar**)&xx);
76: VecRestoreArrayRead(yin,(const PetscScalar**)&yy);
77: VecRestoreArray(win,&ww);
78: return(0);
79: }
81: #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);
124: for (i=0; i<n; i++) {
125: ww[i] = xx[i] / yy[i];
126: }
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: }
160: PetscErrorCode VecConjugate_Seq(Vec xin)
161: {
162: PetscScalar *x;
163: PetscInt n = xin->map->n;
167: VecGetArray(xin,&x);
168: while (n-->0) {
169: *x = PetscConj(*x);
170: x++;
171: }
172: VecRestoreArray(xin,&x);
173: return(0);
174: }
178: PetscErrorCode VecResetArray_Seq(Vec vin)
179: {
180: Vec_Seq *v = (Vec_Seq *)vin->data;
183: v->array = v->unplacedarray;
184: v->unplacedarray = 0;
185: return(0);
186: }
190: PetscErrorCode VecCopy_Seq(Vec xin,Vec yin)
191: {
192: PetscScalar *ya;
193: const PetscScalar *xa;
194: PetscErrorCode ierr;
197: if (xin != yin) {
198: VecGetArrayRead(xin,&xa);
199: VecGetArray(yin,&ya);
200: PetscMemcpy(ya,xa,xin->map->n*sizeof(PetscScalar));
201: VecRestoreArrayRead(xin,&xa);
202: VecRestoreArray(yin,&ya);
203: }
204: return(0);
205: }
209: PetscErrorCode VecSwap_Seq(Vec xin,Vec yin)
210: {
211: PetscScalar *ya, *xa;
213: PetscBLASInt one = 1,bn = PetscBLASIntCast(xin->map->n);
216: if (xin != yin) {
217: VecGetArray(xin,&xa);
218: VecGetArray(yin,&ya);
219: BLASswap_(&bn,xa,&one,ya,&one);
220: VecRestoreArray(xin,&xa);
221: VecRestoreArray(yin,&ya);
222: }
223: return(0);
224: }
226: #include <../src/vec/vec/impls/seq/ftn-kernels/fnorm.h>
229: PetscErrorCode VecNorm_Seq(Vec xin,NormType type,PetscReal* z)
230: {
231: const PetscScalar *xx;
232: PetscErrorCode ierr;
233: PetscInt n = xin->map->n;
234: PetscBLASInt one = 1, bn = PetscBLASIntCast(n);
237: if (type == NORM_2 || type == NORM_FROBENIUS) {
238: VecGetArrayRead(xin,&xx);
239: *z = PetscRealPart(BLASdot_(&bn,xx,&one,xx,&one));
240: *z = PetscSqrtReal(*z);
241: VecRestoreArrayRead(xin,&xx);
242: PetscLogFlops(PetscMax(2.0*n-1,0.0));
243: } else if (type == NORM_INFINITY) {
244: PetscInt i;
245: PetscReal max = 0.0,tmp;
247: VecGetArrayRead(xin,&xx);
248: for (i=0; i<n; i++) {
249: if ((tmp = PetscAbsScalar(*xx)) > max) max = tmp;
250: /* check special case of tmp == NaN */
251: if (tmp != tmp) {max = tmp; break;}
252: xx++;
253: }
254: VecRestoreArrayRead(xin,&xx);
255: *z = max;
256: } else if (type == NORM_1) {
257: VecGetArrayRead(xin,&xx);
258: *z = BLASasum_(&bn,xx,&one);
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: }
270: PetscErrorCode VecView_Seq_ASCII(Vec xin,PetscViewer viewer)
271: {
272: PetscErrorCode ierr;
273: PetscInt i,n = xin->map->n;
274: const char *name;
275: PetscViewerFormat format;
276: const PetscScalar *xv;
279: VecGetArrayRead(xin,&xv);
280: PetscViewerGetFormat(viewer,&format);
281: if (format == PETSC_VIEWER_ASCII_MATLAB) {
282: PetscObjectGetName((PetscObject)xin,&name);
283: PetscViewerASCIIPrintf(viewer,"%s = [\n",name);
284: for (i=0; i<n; i++) {
285: #if defined(PETSC_USE_COMPLEX)
286: if (PetscImaginaryPart(xv[i]) > 0.0) {
287: PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei\n",PetscRealPart(xv[i]),PetscImaginaryPart(xv[i]));
288: } else if (PetscImaginaryPart(xv[i]) < 0.0) {
289: PetscViewerASCIIPrintf(viewer,"%18.16e - %18.16ei\n",PetscRealPart(xv[i]),-PetscImaginaryPart(xv[i]));
290: } else {
291: PetscViewerASCIIPrintf(viewer,"%18.16e\n",PetscRealPart(xv[i]));
292: }
293: #else
294: PetscViewerASCIIPrintf(viewer,"%18.16e\n",(double) xv[i]);
295: #endif
296: }
297: PetscViewerASCIIPrintf(viewer,"];\n");
298: } else if (format == PETSC_VIEWER_ASCII_SYMMODU) {
299: for (i=0; i<n; i++) {
300: #if defined(PETSC_USE_COMPLEX)
301: PetscViewerASCIIPrintf(viewer,"%18.16e %18.16e\n",PetscRealPart(xv[i]),PetscImaginaryPart(xv[i]));
302: #else
303: PetscViewerASCIIPrintf(viewer,"%18.16e\n",xv[i]);
304: #endif
305: }
306: } else if (format == PETSC_VIEWER_ASCII_VTK || format == PETSC_VIEWER_ASCII_VTK_CELL) {
307: /*
308: state 0: No header has been output
309: state 1: Only POINT_DATA has been output
310: state 2: Only CELL_DATA has been output
311: state 3: Output both, POINT_DATA last
312: state 4: Output both, CELL_DATA last
313: */
314: static PetscInt stateId = -1;
315: int outputState = 0;
316: PetscBool hasState;
317: int doOutput = 0;
318: PetscInt bs, b;
320: if (stateId < 0) {
321: PetscObjectComposedDataRegister(&stateId);
322: }
323: PetscObjectComposedDataGetInt((PetscObject) viewer, stateId, outputState, hasState);
324: if (!hasState) {
325: outputState = 0;
326: }
327: PetscObjectGetName((PetscObject) xin, &name);
328: VecGetBlockSize(xin, &bs);
329: if ((bs < 1) || (bs > 3)) {
330: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "VTK can only handle 3D objects, but vector dimension is %d", bs);
331: }
332: if (format == PETSC_VIEWER_ASCII_VTK) {
333: if (outputState == 0) {
334: outputState = 1;
335: doOutput = 1;
336: } else if (outputState == 1) {
337: doOutput = 0;
338: } else if (outputState == 2) {
339: outputState = 3;
340: doOutput = 1;
341: } else if (outputState == 3) {
342: doOutput = 0;
343: } else if (outputState == 4) {
344: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Tried to output POINT_DATA again after intervening CELL_DATA");
345: }
346: if (doOutput) {
347: PetscViewerASCIIPrintf(viewer, "POINT_DATA %d\n", n/bs);
348: }
349: } else {
350: if (outputState == 0) {
351: outputState = 2;
352: doOutput = 1;
353: } else if (outputState == 1) {
354: outputState = 4;
355: doOutput = 1;
356: } else if (outputState == 2) {
357: doOutput = 0;
358: } else if (outputState == 3) {
359: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Tried to output CELL_DATA again after intervening POINT_DATA");
360: } else if (outputState == 4) {
361: doOutput = 0;
362: }
363: if (doOutput) {
364: PetscViewerASCIIPrintf(viewer, "CELL_DATA %d\n", n);
365: }
366: }
367: PetscObjectComposedDataSetInt((PetscObject) viewer, stateId, outputState);
368: if (name) {
369: if (bs == 3) {
370: PetscViewerASCIIPrintf(viewer, "VECTORS %s double\n", name);
371: } else {
372: PetscViewerASCIIPrintf(viewer, "SCALARS %s double %d\n", name, bs);
373: }
374: } else {
375: PetscViewerASCIIPrintf(viewer, "SCALARS scalars double %d\n", bs);
376: }
377: if (bs != 3) {
378: PetscViewerASCIIPrintf(viewer, "LOOKUP_TABLE default\n");
379: }
380: for (i=0; i<n/bs; i++) {
381: for (b=0; b<bs; b++) {
382: if (b > 0) {
383: PetscViewerASCIIPrintf(viewer," ");
384: }
385: #if !defined(PETSC_USE_COMPLEX)
386: PetscViewerASCIIPrintf(viewer,"%G",xv[i*bs+b]);
387: #endif
388: }
389: PetscViewerASCIIPrintf(viewer,"\n");
390: }
391: } else if (format == PETSC_VIEWER_ASCII_VTK_COORDS) {
392: PetscInt bs, b;
394: VecGetBlockSize(xin, &bs);
395: if ((bs < 1) || (bs > 3)) {
396: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "VTK can only handle 3D objects, but vector dimension is %d", bs);
397: }
398: for (i=0; i<n/bs; i++) {
399: for (b=0; b<bs; b++) {
400: if (b > 0) {
401: PetscViewerASCIIPrintf(viewer," ");
402: }
403: #if !defined(PETSC_USE_COMPLEX)
404: PetscViewerASCIIPrintf(viewer,"%G",xv[i*bs+b]);
405: #endif
406: }
407: for (b=bs; b<3; b++) {
408: PetscViewerASCIIPrintf(viewer," 0.0");
409: }
410: PetscViewerASCIIPrintf(viewer,"\n");
411: }
412: } else if (format == PETSC_VIEWER_ASCII_PCICE) {
413: PetscInt bs, b;
415: VecGetBlockSize(xin, &bs);
416: if ((bs < 1) || (bs > 3)) {
417: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "PCICE can only handle up to 3D objects, but vector dimension is %d", bs);
418: }
419: PetscViewerASCIIPrintf(viewer,"%D\n", xin->map->N/bs);
420: for (i=0; i<n/bs; i++) {
421: PetscViewerASCIIPrintf(viewer,"%7D ", i+1);
422: for (b=0; b<bs; b++) {
423: if (b > 0) {
424: PetscViewerASCIIPrintf(viewer," ");
425: }
426: #if !defined(PETSC_USE_COMPLEX)
427: PetscViewerASCIIPrintf(viewer,"% 12.5E",xv[i*bs+b]);
428: #endif
429: }
430: PetscViewerASCIIPrintf(viewer,"\n");
431: }
432: } else if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
433: return(0);
434: } else {
435: PetscObjectPrintClassNamePrefixType((PetscObject)xin,viewer,"Vector Object");
436: for (i=0; i<n; i++) {
437: if (format == PETSC_VIEWER_ASCII_INDEX) {
438: PetscViewerASCIIPrintf(viewer,"%D: ",i);
439: }
440: #if defined(PETSC_USE_COMPLEX)
441: if (PetscImaginaryPart(xv[i]) > 0.0) {
442: PetscViewerASCIIPrintf(viewer,"%G + %G i\n",PetscRealPart(xv[i]),PetscImaginaryPart(xv[i]));
443: } else if (PetscImaginaryPart(xv[i]) < 0.0) {
444: PetscViewerASCIIPrintf(viewer,"%G - %G i\n",PetscRealPart(xv[i]),-PetscImaginaryPart(xv[i]));
445: } else {
446: PetscViewerASCIIPrintf(viewer,"%G\n",PetscRealPart(xv[i]));
447: }
448: #else
449: PetscViewerASCIIPrintf(viewer,"%G\n",xv[i]);
450: #endif
451: }
452: }
453: PetscViewerFlush(viewer);
454: VecRestoreArrayRead(xin,&xv);
455: return(0);
456: }
460: PetscErrorCode VecView_Seq_Draw_LG(Vec xin,PetscViewer v)
461: {
462: PetscErrorCode ierr;
463: PetscInt i,n = xin->map->n;
464: PetscDraw win;
465: PetscReal *xx;
466: PetscDrawLG lg;
467: const PetscScalar *xv;
470: PetscViewerDrawGetDrawLG(v,0,&lg);
471: PetscDrawLGGetDraw(lg,&win);
472: PetscDrawCheckResizedWindow(win);
473: PetscDrawLGReset(lg);
474: PetscMalloc((n+1)*sizeof(PetscReal),&xx);
475: for (i=0; i<n; i++) {
476: xx[i] = (PetscReal) i;
477: }
478: VecGetArrayRead(xin,&xv);
479: #if !defined(PETSC_USE_COMPLEX)
480: PetscDrawLGAddPoints(lg,n,&xx,(PetscReal**)&xv);
481: #else
482: {
483: PetscReal *yy;
484: PetscMalloc((n+1)*sizeof(PetscReal),&yy);
485: for (i=0; i<n; i++) {
486: yy[i] = PetscRealPart(xv[i]);
487: }
488: PetscDrawLGAddPoints(lg,n,&xx,&yy);
489: PetscFree(yy);
490: }
491: #endif
492: VecRestoreArrayRead(xin,&xv);
493: PetscFree(xx);
494: PetscDrawLGDraw(lg);
495: PetscDrawSynchronizedFlush(win);
496: return(0);
497: }
501: PetscErrorCode VecView_Seq_Draw(Vec xin,PetscViewer v)
502: {
503: PetscErrorCode ierr;
504: PetscDraw draw;
505: PetscBool isnull;
506: PetscViewerFormat format;
509: PetscViewerDrawGetDraw(v,0,&draw);
510: PetscDrawIsNull(draw,&isnull); if (isnull) return(0);
511:
512: PetscViewerGetFormat(v,&format);
513: /*
514: Currently it only supports drawing to a line graph */
515: if (format != PETSC_VIEWER_DRAW_LG) {
516: PetscViewerPushFormat(v,PETSC_VIEWER_DRAW_LG);
517: }
518: VecView_Seq_Draw_LG(xin,v);
519: if (format != PETSC_VIEWER_DRAW_LG) {
520: PetscViewerPopFormat(v);
521: }
523: return(0);
524: }
528: PetscErrorCode VecView_Seq_Binary(Vec xin,PetscViewer viewer)
529: {
530: PetscErrorCode ierr;
531: int fdes;
532: PetscInt n = xin->map->n,classid=VEC_FILE_CLASSID;
533: FILE *file;
534: const PetscScalar *xv;
535: #if defined(PETSC_HAVE_MPIIO)
536: PetscBool isMPIIO;
537: #endif
538: PetscBool skipHeader;
542: /* Write vector header */
543: PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
544: if (!skipHeader) {
545: PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);
546: PetscViewerBinaryWrite(viewer,&n,1,PETSC_INT,PETSC_FALSE);
547: }
549: /* Write vector contents */
550: #if defined(PETSC_HAVE_MPIIO)
551: PetscViewerBinaryGetMPIIO(viewer,&isMPIIO);
552: if (!isMPIIO) {
553: #endif
554: PetscViewerBinaryGetDescriptor(viewer,&fdes);
555: VecGetArrayRead(xin,&xv);
556: PetscBinaryWrite(fdes,(void*)xv,n,PETSC_SCALAR,PETSC_FALSE);
557: VecRestoreArrayRead(xin,&xv);
558: #if defined(PETSC_HAVE_MPIIO)
559: } else {
560: MPI_Offset off;
561: MPI_File mfdes;
562: PetscMPIInt gsizes[1],lsizes[1],lstarts[1];
563: MPI_Datatype view;
565: gsizes[0] = PetscMPIIntCast(n);
566: lsizes[0] = PetscMPIIntCast(n);
567: lstarts[0] = 0;
568: MPI_Type_create_subarray(1,gsizes,lsizes,lstarts,MPI_ORDER_FORTRAN,MPIU_SCALAR,&view);
569: MPI_Type_commit(&view);
571: PetscViewerBinaryGetMPIIODescriptor(viewer,&mfdes);
572: PetscViewerBinaryGetMPIIOOffset(viewer,&off);
573: VecGetArrayRead(xin,&xv);
574: MPIU_File_write_all(mfdes,(void*)xv,lsizes[0],MPIU_SCALAR,MPI_STATUS_IGNORE);
575: VecRestoreArrayRead(xin,&xv);
576: PetscViewerBinaryAddMPIIOOffset(viewer,n*sizeof(PetscScalar));
577: MPI_Type_free(&view);
578: }
579: #endif
581: PetscViewerBinaryGetInfoPointer(viewer,&file);
582: if (file) {
583: if (((PetscObject)xin)->prefix) {
584: PetscFPrintf(PETSC_COMM_SELF,file,"-%s_vecload_block_size %D\n",((PetscObject)xin)->prefix,xin->map->bs);
585: } else {
586: PetscFPrintf(PETSC_COMM_SELF,file,"-vecload_block_size %D\n",xin->map->bs);
587: }
588: }
589: return(0);
590: }
592: #if defined(PETSC_HAVE_MATLAB_ENGINE)
593: #include <mat.h> /* MATLAB include file */
594: EXTERN_C_BEGIN
597: PetscErrorCode VecView_Seq_Matlab(Vec vec,PetscViewer viewer)
598: {
599: PetscErrorCode ierr;
600: PetscInt n;
601: const PetscScalar *array;
602:
604: VecGetLocalSize(vec,&n);
605: PetscObjectName((PetscObject)vec);
606: VecGetArrayRead(vec,&array);
607: PetscViewerMatlabPutArray(viewer,n,1,array,((PetscObject)vec)->name);
608: VecRestoreArrayRead(vec,&array);
609: return(0);
610: }
611: EXTERN_C_END
612: #endif
616: PetscErrorCode VecView_Seq(Vec xin,PetscViewer viewer)
617: {
619: PetscBool isdraw,iascii,issocket,isbinary;
620: #if defined(PETSC_HAVE_MATHEMATICA)
621: PetscBool ismathematica;
622: #endif
623: #if defined(PETSC_HAVE_MATLAB_ENGINE)
624: PetscBool ismatlab;
625: #endif
626: #if defined(PETSC_HAVE_HDF5)
627: PetscBool ishdf5;
628: #endif
631: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
632: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
633: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSOCKET,&issocket);
634: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
635: #if defined(PETSC_HAVE_MATHEMATICA)
636: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERMATHEMATICA,&ismathematica);
637: #endif
638: #if defined(PETSC_HAVE_HDF5)
639: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
640: #endif
641: #if defined(PETSC_HAVE_MATLAB_ENGINE)
642: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERMATLAB,&ismatlab);
643: #endif
645: if (isdraw){
646: VecView_Seq_Draw(xin,viewer);
647: } else if (iascii){
648: VecView_Seq_ASCII(xin,viewer);
649: } else if (isbinary) {
650: VecView_Seq_Binary(xin,viewer);
651: #if defined(PETSC_HAVE_MATHEMATICA)
652: } else if (ismathematica) {
653: PetscViewerMathematicaPutVector(viewer,xin);
654: #endif
655: #if defined(PETSC_HAVE_HDF5)
656: } else if (ishdf5) {
657: VecView_MPI_HDF5(xin,viewer); /* Reusing VecView_MPI_HDF5 ... don't want code duplication*/
658: #endif
659: #if defined(PETSC_HAVE_MATLAB_ENGINE)
660: } else if (ismatlab) {
661: VecView_Seq_Matlab(xin,viewer);
662: #endif
663: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported by this vector object",((PetscObject)viewer)->type_name);
664: return(0);
665: }
669: PetscErrorCode VecGetValues_Seq(Vec xin,PetscInt ni,const PetscInt ix[],PetscScalar y[])
670: {
671: const PetscScalar *xx;
672: PetscInt i;
673: PetscErrorCode ierr;
676: VecGetArrayRead(xin,&xx);
677: for (i=0; i<ni; i++) {
678: if (xin->stash.ignorenegidx && ix[i] < 0) continue;
679: #if defined(PETSC_USE_DEBUG)
680: if (ix[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D cannot be negative",ix[i]);
681: 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);
682: #endif
683: y[i] = xx[ix[i]];
684: }
685: VecRestoreArrayRead(xin,&xx);
686: return(0);
687: }
691: PetscErrorCode VecSetValues_Seq(Vec xin,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode m)
692: {
693: PetscScalar *xx;
694: PetscInt i;
698: VecGetArray(xin,&xx);
699: if (m == INSERT_VALUES) {
700: for (i=0; i<ni; i++) {
701: if (xin->stash.ignorenegidx && ix[i] < 0) continue;
702: #if defined(PETSC_USE_DEBUG)
703: if (ix[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D cannot be negative",ix[i]);
704: 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);
705: #endif
706: xx[ix[i]] = y[i];
707: }
708: } else {
709: for (i=0; i<ni; i++) {
710: if (xin->stash.ignorenegidx && ix[i] < 0) continue;
711: #if defined(PETSC_USE_DEBUG)
712: if (ix[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D cannot be negative",ix[i]);
713: 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);
714: #endif
715: xx[ix[i]] += y[i];
716: }
717: }
718: VecRestoreArray(xin,&xx);
719: return(0);
720: }
724: PetscErrorCode VecSetValuesBlocked_Seq(Vec xin,PetscInt ni,const PetscInt ix[],const PetscScalar yin[],InsertMode m)
725: {
726: PetscScalar *xx,*y = (PetscScalar*)yin;
727: PetscInt i,bs = xin->map->bs,start,j;
730: /*
731: For optimization could treat bs = 2, 3, 4, 5 as special cases with loop unrolling
732: */
734: VecGetArray(xin,&xx);
735: if (m == INSERT_VALUES) {
736: for (i=0; i<ni; i++) {
737: start = bs*ix[i];
738: if (start < 0) continue;
739: #if defined(PETSC_USE_DEBUG)
740: 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);
741: #endif
742: for (j=0; j<bs; j++) {
743: xx[start+j] = y[j];
744: }
745: y += bs;
746: }
747: } else {
748: for (i=0; i<ni; i++) {
749: start = bs*ix[i];
750: if (start < 0) continue;
751: #if defined(PETSC_USE_DEBUG)
752: 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);
753: #endif
754: for (j=0; j<bs; j++) {
755: xx[start+j] += y[j];
756: }
757: y += bs;
758: }
759: }
760: VecRestoreArray(xin,&xx);
761: return(0);
762: }
766: PetscErrorCode VecDestroy_Seq(Vec v)
767: {
768: Vec_Seq *vs = (Vec_Seq*)v->data;
772: PetscObjectDepublish(v);
774: #if defined(PETSC_USE_LOG)
775: PetscLogObjectState((PetscObject)v,"Length=%D",v->map->n);
776: #endif
777: PetscFree(vs->array_allocated);
778: PetscFree(v->data);
779: return(0);
780: }
784: PetscErrorCode VecSetOption_Seq(Vec v,VecOption op,PetscBool flag)
785: {
787: if (op == VEC_IGNORE_NEGATIVE_INDICES) {
788: v->stash.ignorenegidx = flag;
789: }
790: return(0);
791: }
795: PetscErrorCode VecDuplicate_Seq(Vec win,Vec *V)
796: {
800: VecCreate(((PetscObject)win)->comm,V);
801: PetscObjectSetPrecision((PetscObject)*V,((PetscObject)win)->precision);
802: VecSetSizes(*V,win->map->n,win->map->n);
803: VecSetType(*V,((PetscObject)win)->type_name);
804: PetscLayoutReference(win->map,&(*V)->map);
805: PetscOListDuplicate(((PetscObject)win)->olist,&((PetscObject)(*V))->olist);
806: PetscFListDuplicate(((PetscObject)win)->qlist,&((PetscObject)(*V))->qlist);
808: (*V)->stash.ignorenegidx = win->stash.ignorenegidx;
810: return(0);
811: }
813: static struct _VecOps DvOps = {VecDuplicate_Seq, /* 1 */
814: VecDuplicateVecs_Default,
815: VecDestroyVecs_Default,
816: VecDot_Seq,
817: VecMDot_Seq,
818: VecNorm_Seq,
819: VecTDot_Seq,
820: VecMTDot_Seq,
821: VecScale_Seq,
822: VecCopy_Seq, /* 10 */
823: VecSet_Seq,
824: VecSwap_Seq,
825: VecAXPY_Seq,
826: VecAXPBY_Seq,
827: VecMAXPY_Seq,
828: VecAYPX_Seq,
829: VecWAXPY_Seq,
830: VecAXPBYPCZ_Seq,
831: VecPointwiseMult_Seq,
832: VecPointwiseDivide_Seq,
833: VecSetValues_Seq, /* 20 */
834: 0,0,
835: 0,
836: VecGetSize_Seq,
837: VecGetSize_Seq,
838: 0,
839: VecMax_Seq,
840: VecMin_Seq,
841: VecSetRandom_Seq,
842: VecSetOption_Seq, /* 30 */
843: VecSetValuesBlocked_Seq,
844: VecDestroy_Seq,
845: VecView_Seq,
846: VecPlaceArray_Seq,
847: VecReplaceArray_Seq,
848: VecDot_Seq,
849: VecTDot_Seq,
850: VecNorm_Seq,
851: VecMDot_Seq,
852: VecMTDot_Seq, /* 40 */
853: VecLoad_Default,
854: VecReciprocal_Default,
855: VecConjugate_Seq,
856: 0,
857: 0,
858: VecResetArray_Seq,
859: 0,
860: VecMaxPointwiseDivide_Seq,
861: VecPointwiseMax_Seq,
862: VecPointwiseMaxAbs_Seq,
863: VecPointwiseMin_Seq,
864: VecGetValues_Seq,
865: 0,
866: 0,
867: 0,
868: 0,
869: 0,
870: 0,
871: VecStrideGather_Default,
872: VecStrideScatter_Default
873: };
876: /*
877: This is called by VecCreate_Seq() (i.e. VecCreateSeq()) and VecCreateSeqWithArray()
878: */
881: PetscErrorCode VecCreate_Seq_Private(Vec v,const PetscScalar array[])
882: {
883: Vec_Seq *s;
887: PetscNewLog(v,Vec_Seq,&s);
888: 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: PetscObjectComposeFunctionDynamic((PetscObject)v,"PetscMatlabEnginePut_C","VecMatlabEnginePut_Default",VecMatlabEnginePut_Default);
898: PetscObjectComposeFunctionDynamic((PetscObject)v,"PetscMatlabEngineGet_C","VecMatlabEngineGet_Default",VecMatlabEngineGet_Default);
899: #endif
900: #if defined(PETSC_USE_MIXED_PRECISION)
901: ((PetscObject)v)->precision = (PetscPrecision)sizeof(PetscScalar);
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: . n - the vector length
917: - array - memory where the vector elements are to be stored.
919: Output Parameter:
920: . V - the vector
922: Notes:
923: Use VecDuplicate() or VecDuplicateVecs() to form additional vectors of the
924: same type as an existing vector.
926: If the user-provided array is PETSC_NULL, then VecPlaceArray() can be used
927: at a later stage to SET the array for storing the vector values.
929: PETSc does NOT free the array when the vector is destroyed via VecDestroy().
930: The user should not free the array until the vector is destroyed.
932: Level: intermediate
934: Concepts: vectors^creating with array
936: .seealso: VecCreateMPIWithArray(), VecCreate(), VecDuplicate(), VecDuplicateVecs(),
937: VecCreateGhost(), VecCreateSeq(), VecPlaceArray()
938: @*/
939: PetscErrorCode VecCreateSeqWithArray(MPI_Comm comm,PetscInt bs,PetscInt n,const PetscScalar array[],Vec *V)
940: {
942: PetscMPIInt size;
945: VecCreate(comm,V);
946: VecSetSizes(*V,n,n);
947: VecSetBlockSize(*V,bs);
948: MPI_Comm_size(comm,&size);
949: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Cannot create VECSEQ on more than one process");
950: VecCreate_Seq_Private(*V,array);
951: return(0);
952: }