Actual source code: bvec2.c
petsc-3.5.4 2015-05-23
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>
9: #include <petscthreadcomm.h>
11: #if defined(PETSC_HAVE_HDF5)
12: extern PetscErrorCode VecView_MPI_HDF5(Vec,PetscViewer);
13: #endif
15: #if defined(PETSC_THREADCOMM_ACTIVE)
16: PetscErrorCode VecPointwiseMax_kernel(PetscInt thread_id,Vec win,Vec xin,Vec yin)
17: {
18: PetscErrorCode ierr;
19: PetscInt *trstarts=win->map->trstarts;
20: PetscInt i;
21: const PetscScalar *xx,*yy;
22: PetscScalar *ww;
24: VecGetArrayRead(xin,&xx);
25: VecGetArrayRead(yin,&yy);
26: VecGetArray(win,&ww);
28: for (i=trstarts[thread_id]; i < trstarts[thread_id+1];i++) ww[i] = PetscMax(PetscRealPart(xx[i]),PetscRealPart(yy[i]));
30: VecRestoreArrayRead(xin,&xx);
31: VecRestoreArrayRead(yin,&yy);
32: VecRestoreArray(win,&ww);
33: return 0;
34: }
38: PetscErrorCode VecPointwiseMax_Seq(Vec win,Vec xin,Vec yin)
39: {
43: PetscThreadCommRunKernel3(PetscObjectComm((PetscObject)win),(PetscThreadKernel)VecPointwiseMax_kernel,win,xin,yin);
44: PetscLogFlops(win->map->n);
45: return(0);
46: }
47: #else
50: PetscErrorCode VecPointwiseMax_Seq(Vec win,Vec xin,Vec yin)
51: {
53: PetscInt n = win->map->n,i;
54: 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(PetscRealPart(xx[i]),PetscRealPart(yy[i]));
63: VecRestoreArrayRead(xin,(const PetscScalar**)&xx);
64: VecRestoreArrayRead(yin,(const PetscScalar**)&yy);
65: VecRestoreArray(win,&ww);
66: PetscLogFlops(n);
67: return(0);
68: }
69: #endif
71: #if defined(PETSC_THREADCOMM_ACTIVE)
72: PetscErrorCode VecPointwiseMin_kernel(PetscInt thread_id,Vec win,Vec xin,Vec yin)
73: {
74: PetscErrorCode ierr;
75: PetscInt *trstarts=win->map->trstarts;
76: PetscInt i;
77: const PetscScalar *xx,*yy;
78: PetscScalar *ww;
80: VecGetArrayRead(xin,&xx);
81: VecGetArrayRead(yin,&yy);
82: VecGetArray(win,&ww);
84: for (i=trstarts[thread_id]; i < trstarts[thread_id+1];i++) ww[i] = PetscMin(PetscRealPart(xx[i]),PetscRealPart(yy[i]));
86: VecRestoreArrayRead(xin,&xx);
87: VecRestoreArrayRead(yin,&yy);
88: VecRestoreArray(win,&ww);
89: return 0;
90: }
94: PetscErrorCode VecPointwiseMin_Seq(Vec win,Vec xin,Vec yin)
95: {
99: PetscThreadCommRunKernel3(PetscObjectComm((PetscObject)win),(PetscThreadKernel)VecPointwiseMin_kernel,win,xin,yin);
100: PetscLogFlops(win->map->n);
101: return(0);
102: }
103: #else
106: PetscErrorCode VecPointwiseMin_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++) ww[i] = PetscMin(PetscRealPart(xx[i]),PetscRealPart(yy[i]));
119: VecRestoreArrayRead(xin,(const PetscScalar**)&xx);
120: VecRestoreArrayRead(yin,(const PetscScalar**)&yy);
121: VecRestoreArray(win,&ww);
122: PetscLogFlops(n);
123: return(0);
124: }
125: #endif
127: #if defined(PETSC_THREADCOMM_ACTIVE)
128: PetscErrorCode VecPointwiseMaxAbs_kernel(PetscInt thread_id,Vec win,Vec xin,Vec yin)
129: {
130: PetscErrorCode ierr;
131: PetscInt *trstarts=win->map->trstarts;
132: PetscInt i;
133: const PetscScalar *xx,*yy;
134: PetscScalar *ww;
136: VecGetArrayRead(xin,&xx);
137: VecGetArrayRead(yin,&yy);
138: VecGetArray(win,&ww);
140: for (i=trstarts[thread_id]; i < trstarts[thread_id+1];i++) ww[i] = PetscMax(PetscAbsScalar(xx[i]),PetscAbsScalar(yy[i]));
142: VecRestoreArrayRead(xin,&xx);
143: VecRestoreArrayRead(yin,&yy);
144: VecRestoreArray(win,&ww);
145: return 0;
146: }
150: PetscErrorCode VecPointwiseMaxAbs_Seq(Vec win,Vec xin,Vec yin)
151: {
155: PetscThreadCommRunKernel3(PetscObjectComm((PetscObject)win),(PetscThreadKernel)VecPointwiseMaxAbs_kernel,win,xin,yin);
156: PetscLogFlops(win->map->n);
157: return(0);
158: }
159: #else
162: PetscErrorCode VecPointwiseMaxAbs_Seq(Vec win,Vec xin,Vec yin)
163: {
165: PetscInt n = win->map->n,i;
166: PetscScalar *ww,*xx,*yy; /* cannot make xx or yy const since might be ww */
169: VecGetArrayRead(xin,(const PetscScalar**)&xx);
170: VecGetArrayRead(yin,(const PetscScalar**)&yy);
171: VecGetArray(win,&ww);
173: for (i=0; i<n; i++) ww[i] = PetscMax(PetscAbsScalar(xx[i]),PetscAbsScalar(yy[i]));
175: PetscLogFlops(n);
176: VecRestoreArrayRead(xin,(const PetscScalar**)&xx);
177: VecRestoreArrayRead(yin,(const PetscScalar**)&yy);
178: VecRestoreArray(win,&ww);
179: return(0);
180: }
181: #endif
183: #include <../src/vec/vec/impls/seq/ftn-kernels/fxtimesy.h>
184: #if defined(PETSC_THREADCOMM_ACTIVE)
185: PetscErrorCode VecPointwiseMult_kernel(PetscInt thread_id,Vec win,Vec xin,Vec yin)
186: {
187: PetscErrorCode ierr;
188: PetscInt *trstarts=win->map->trstarts;
189: PetscInt i;
190: const PetscScalar *xx,*yy;
191: PetscScalar *ww;
193: VecGetArrayRead(xin,&xx);
194: VecGetArrayRead(yin,&yy);
195: VecGetArray(win,&ww);
196: if (ww == xx) {
197: for (i=trstarts[thread_id]; i < trstarts[thread_id+1]; i++) ww[i] *= yy[i];
198: } else if (ww == yy) {
199: for (i=trstarts[thread_id]; i < trstarts[thread_id+1]; i++) ww[i] *= xx[i];
200: } else {
201: #if defined(PETSC_USE_FORTRAN_KERNEL_XTIMESY)
202: PetscInt start,n;
203: start = trstarts[thread_id];
204: n = trstarts[thread_id+1] - trstarts[thread_id];
205: fortranxtimesy_(xx+start,yy+start,ww+start,&n);
206: }
207: #else
208: for (i=trstarts[thread_id]; i < trstarts[thread_id+1]; i++) ww[i] = xx[i] * yy[i];
209: }
210: #endif
211: VecRestoreArrayRead(xin,&xx);
212: VecRestoreArrayRead(yin,&yy);
213: VecRestoreArray(win,&ww);
214: return 0;
215: }
219: PetscErrorCode VecPointwiseMult_Seq(Vec win,Vec xin,Vec yin)
220: {
224: PetscThreadCommRunKernel3(PetscObjectComm((PetscObject)win),(PetscThreadKernel)VecPointwiseMult_kernel,win,xin,yin);
225: PetscLogFlops(win->map->n);
226: return(0);
227: }
228: #else
231: PetscErrorCode VecPointwiseMult_Seq(Vec win,Vec xin,Vec yin)
232: {
234: PetscInt n = win->map->n,i;
235: PetscScalar *ww,*xx,*yy; /* cannot make xx or yy const since might be ww */
238: VecGetArrayRead(xin,(const PetscScalar**)&xx);
239: VecGetArrayRead(yin,(const PetscScalar**)&yy);
240: VecGetArray(win,&ww);
241: if (ww == xx) {
242: for (i=0; i<n; i++) ww[i] *= yy[i];
243: } else if (ww == yy) {
244: for (i=0; i<n; i++) ww[i] *= xx[i];
245: } else {
246: #if defined(PETSC_USE_FORTRAN_KERNEL_XTIMESY)
247: fortranxtimesy_(xx,yy,ww,&n);
248: #else
249: for (i=0; i<n; i++) ww[i] = xx[i] * yy[i];
250: #endif
251: }
252: VecRestoreArrayRead(xin,(const PetscScalar**)&xx);
253: VecRestoreArrayRead(yin,(const PetscScalar**)&yy);
254: VecRestoreArray(win,&ww);
255: PetscLogFlops(n);
256: return(0);
257: }
258: #endif
260: #if defined(PETSC_THREADCOMM_ACTIVE)
261: PetscErrorCode VecPointwiseDivide_kernel(PetscInt thread_id,Vec win,Vec xin,Vec yin)
262: {
263: PetscErrorCode ierr;
264: PetscInt *trstarts=win->map->trstarts;
265: PetscInt i;
266: const PetscScalar *xx,*yy;
267: PetscScalar *ww;
269: VecGetArrayRead(xin,&xx);
270: VecGetArrayRead(yin,&yy);
271: VecGetArray(win,&ww);
273: for (i=trstarts[thread_id]; i < trstarts[thread_id+1]; i++) ww[i] = xx[i] / yy[i];
275: VecRestoreArrayRead(xin,&xx);
276: VecRestoreArrayRead(yin,&yy);
277: VecRestoreArray(win,&ww);
278: return 0;
279: }
283: PetscErrorCode VecPointwiseDivide_Seq(Vec win,Vec xin,Vec yin)
284: {
288: PetscThreadCommRunKernel3(PetscObjectComm((PetscObject)win),(PetscThreadKernel)VecPointwiseDivide_kernel,win,xin,yin);
289: PetscLogFlops(win->map->n);
290: return(0);
291: }
292: #else
295: PetscErrorCode VecPointwiseDivide_Seq(Vec win,Vec xin,Vec yin)
296: {
298: PetscInt n = win->map->n,i;
299: PetscScalar *ww,*xx,*yy; /* cannot make xx or yy const since might be ww */
302: VecGetArrayRead(xin,(const PetscScalar**)&xx);
303: VecGetArrayRead(yin,(const PetscScalar**)&yy);
304: VecGetArray(win,&ww);
306: for (i=0; i<n; i++) ww[i] = xx[i] / yy[i];
308: PetscLogFlops(n);
309: VecRestoreArrayRead(xin,(const PetscScalar**)&xx);
310: VecRestoreArrayRead(yin,(const PetscScalar**)&yy);
311: VecRestoreArray(win,&ww);
312: return(0);
313: }
314: #endif
316: #if defined(PETSC_THREADCOMM_ACTIVE)
317: PetscErrorCode VecSetRandom_kernel(PetscInt thread_id,Vec xin, PetscRandom r)
318: {
320: PetscInt *trstarts=xin->map->trstarts;
321: PetscInt i;
322: PetscScalar *xx;
324: VecGetArray(xin,&xx);
325: for (i=trstarts[thread_id]; i < trstarts[thread_id+1]; i++) {
326: PetscRandomGetValue(r,&xx[i]);
327: }
328: VecRestoreArray(xin,&xx);
329: return 0;
330: }
334: PetscErrorCode VecSetRandom_Seq(Vec xin,PetscRandom r)
335: {
339: PetscThreadCommRunKernel2(PetscObjectComm((PetscObject)xin),(PetscThreadKernel)VecSetRandom_kernel,xin,r);
340: return(0);
341: }
342: #else
345: PetscErrorCode VecSetRandom_Seq(Vec xin,PetscRandom r)
346: {
348: PetscInt n = xin->map->n,i;
349: PetscScalar *xx;
352: VecGetArray(xin,&xx);
353: for (i=0; i<n; i++) {PetscRandomGetValue(r,&xx[i]);}
354: VecRestoreArray(xin,&xx);
355: return(0);
356: }
357: #endif
361: PetscErrorCode VecGetSize_Seq(Vec vin,PetscInt *size)
362: {
364: *size = vin->map->n;
365: return(0);
366: }
369: #if defined(PETSC_THREADCOMM_ACTIVE)
370: PetscErrorCode VecConjugate_kernel(PetscInt thread_id,Vec xin)
371: {
373: PetscScalar *x;
374: PetscInt *trstarts=xin->map->trstarts,i;
376: VecGetArray(xin,&x);
377: for (i=trstarts[thread_id]; i < trstarts[thread_id+1]; i++) x[i] = PetscConj(x[i]);
378: VecRestoreArray(xin,&x);
379: return 0;
380: }
384: PetscErrorCode VecConjugate_Seq(Vec xin)
385: {
389: PetscThreadCommRunKernel1(PetscObjectComm((PetscObject)xin),(PetscThreadKernel)VecConjugate_kernel,xin);
390: return(0);
391: }
392: #else
395: PetscErrorCode VecConjugate_Seq(Vec xin)
396: {
397: PetscScalar *x;
398: PetscInt n = xin->map->n;
402: VecGetArray(xin,&x);
403: while (n-->0) {
404: *x = PetscConj(*x);
405: x++;
406: }
407: VecRestoreArray(xin,&x);
408: return(0);
409: }
410: #endif
414: PetscErrorCode VecResetArray_Seq(Vec vin)
415: {
416: Vec_Seq *v = (Vec_Seq*)vin->data;
419: v->array = v->unplacedarray;
420: v->unplacedarray = 0;
421: return(0);
422: }
424: #if defined(PETSC_THREADCOMM_ACTIVE)
425: PetscErrorCode VecCopy_kernel(PetscInt thread_id,Vec xin,Vec yin)
426: {
427: PetscErrorCode ierr;
428: PetscInt *trstarts=yin->map->trstarts;
429: PetscScalar *ya;
430: const PetscScalar *xa;
431: PetscInt start,end,n;
433: VecGetArrayRead(xin,&xa);
434: VecGetArray(yin,&ya);
435: start = trstarts[thread_id];
436: end = trstarts[thread_id+1];
437: n = end-start;
438: PetscMemcpy(ya+start,xa+start,n*sizeof(PetscScalar));
439: VecRestoreArrayRead(xin,&xa);
440: VecRestoreArray(yin,&ya);
441: return 0;
442: }
446: PetscErrorCode VecCopy_Seq(Vec xin,Vec yin)
447: {
451: if (xin != yin) {
452: PetscThreadCommRunKernel2(PetscObjectComm((PetscObject)xin),(PetscThreadKernel)VecCopy_kernel,xin,yin);
453: }
454: return(0);
455: }
456: #else
459: PetscErrorCode VecCopy_Seq(Vec xin,Vec yin)
460: {
461: PetscScalar *ya;
462: const PetscScalar *xa;
463: PetscErrorCode ierr;
466: if (xin != yin) {
467: VecGetArrayRead(xin,&xa);
468: VecGetArray(yin,&ya);
469: PetscMemcpy(ya,xa,xin->map->n*sizeof(PetscScalar));
470: VecRestoreArrayRead(xin,&xa);
471: VecRestoreArray(yin,&ya);
472: }
473: return(0);
474: }
475: #endif
477: #if defined(PETSC_THREADCOMM_ACTIVE)
478: PetscErrorCode VecSwap_kernel(PetscInt thread_id,Vec xin,Vec yin)
479: {
480: PetscErrorCode ierr;
481: PetscInt *trstarts=xin->map->trstarts;
482: PetscInt start,end,n;
483: PetscBLASInt one=1,bn;
484: PetscScalar *xa,*ya;
486: VecGetArray(xin,&xa);
487: VecGetArray(yin,&ya);
488: start = trstarts[thread_id];
489: end = trstarts[thread_id+1];
490: n = end-start;
491: PetscBLASIntCast(n,&bn);
492: PetscStackCallBLAS("BLASswap",BLASswap_(&bn,xa+start,&one,ya+start,&one));
493: VecRestoreArray(xin,&xa);
494: VecRestoreArray(yin,&ya);
495: return 0;
496: }
500: PetscErrorCode VecSwap_Seq(Vec xin,Vec yin)
501: {
505: if (xin != yin) {
506: PetscThreadCommRunKernel2(PetscObjectComm((PetscObject)xin),(PetscThreadKernel)VecSwap_kernel,xin,yin);
507: }
508: return(0);
509: }
510: #else
514: PetscErrorCode VecSwap_Seq(Vec xin,Vec yin)
515: {
516: PetscScalar *ya, *xa;
518: PetscBLASInt one = 1,bn;
521: if (xin != yin) {
522: PetscBLASIntCast(xin->map->n,&bn);
523: VecGetArray(xin,&xa);
524: VecGetArray(yin,&ya);
525: PetscStackCallBLAS("BLASswap",BLASswap_(&bn,xa,&one,ya,&one));
526: VecRestoreArray(xin,&xa);
527: VecRestoreArray(yin,&ya);
528: }
529: return(0);
530: }
531: #endif
533: #include <../src/vec/vec/impls/seq/ftn-kernels/fnorm.h>
534: #if defined(PETSC_THREADCOMM_ACTIVE)
535: PetscErrorCode VecNorm_kernel(PetscInt thread_id,Vec xin,NormType* type_p,PetscThreadCommReduction red)
536: {
538: PetscInt *trstarts=xin->map->trstarts;
539: PetscInt start,end,n;
540: PetscBLASInt one = 1,bn;
541: const PetscScalar *xx;
542: PetscReal z_loc;
543: NormType type=*type_p;
545: start = trstarts[thread_id];
546: end = trstarts[thread_id+1];
547: n = end - start;
548: PetscBLASIntCast(n,&bn);
549: VecGetArrayRead(xin,&xx);
550: if (type == NORM_2 || type == NORM_FROBENIUS) {
551: z_loc = PetscRealPart(BLASdot_(&bn,xx+start,&one,xx+start,&one));
552: PetscThreadReductionKernelPost(thread_id,red,(void*)&z_loc);
553: } else if (type == NORM_INFINITY) {
554: PetscInt i;
555: PetscReal max=0.0,tmp;
556: for (i=trstarts[thread_id]; i < trstarts[thread_id+1]; i++) {
557: if ((tmp = PetscAbsScalar(xx[i])) > max) max = tmp;
558: /* check special case of tmp == NaN */
559: if (tmp != tmp) {max = tmp; break;}
560: }
561: PetscThreadReductionKernelPost(thread_id,red,(void*)&max);
562: } else if (type == NORM_1) {
563: PetscStackCallBLAS("BLASasum",z_loc = BLASasum_(&bn,xx+start,&one));
564: PetscThreadReductionKernelPost(thread_id,red,(void*)&z_loc);
565: }
566: VecRestoreArrayRead(xin,&xx);
567: return 0;
568: }
572: PetscErrorCode VecNorm_Seq(Vec xin,NormType type,PetscReal *z)
573: {
574: PetscErrorCode ierr;
575: PetscThreadCommReduction red;
578: if (type == NORM_2 || type == NORM_FROBENIUS) {
579: PetscThreadReductionBegin(PetscObjectComm((PetscObject)xin),THREADCOMM_SUM,PETSC_REAL,1,&red);
580: PetscThreadCommRunKernel3(PetscObjectComm((PetscObject)xin),(PetscThreadKernel)VecNorm_kernel,xin,(void*)&type,red);
581: PetscThreadReductionEnd(red,z);
582: *z = PetscSqrtReal(*z);
583: PetscLogFlops(PetscMax(2*xin->map->n-1.0,0.0));
584: } else if (type == NORM_INFINITY) {
585: PetscThreadReductionBegin(PetscObjectComm((PetscObject)xin),THREADCOMM_MAX,PETSC_REAL,1,&red);
586: PetscThreadCommRunKernel3(PetscObjectComm((PetscObject)xin),(PetscThreadKernel)VecNorm_kernel,xin,(void*)&type,red);
587: PetscThreadReductionEnd(red,z);
588: } else if (type == NORM_1) {
589: PetscThreadReductionBegin(PetscObjectComm((PetscObject)xin),THREADCOMM_SUM,PETSC_REAL,1,&red);
590: PetscThreadCommRunKernel3(PetscObjectComm((PetscObject)xin),(PetscThreadKernel)VecNorm_kernel,xin,(void*)&type,red);
591: PetscThreadReductionEnd(red,z);
592: PetscLogFlops(PetscMax(xin->map->n-1.0,0.0));
593: } else if (type == NORM_1_AND_2) {
594: VecNorm_Seq(xin,NORM_1,z);
595: VecNorm_Seq(xin,NORM_2,z+1);
596: }
597: return(0);
598: }
599: #else
603: PetscErrorCode VecNorm_Seq(Vec xin,NormType type,PetscReal *z)
604: {
605: const PetscScalar *xx;
606: PetscErrorCode ierr;
607: PetscInt n = xin->map->n;
608: PetscBLASInt one = 1, bn;
611: PetscBLASIntCast(n,&bn);
612: if (type == NORM_2 || type == NORM_FROBENIUS) {
613: VecGetArrayRead(xin,&xx);
614: *z = PetscRealPart(BLASdot_(&bn,xx,&one,xx,&one));
615: *z = PetscSqrtReal(*z);
616: VecRestoreArrayRead(xin,&xx);
617: PetscLogFlops(PetscMax(2.0*n-1,0.0));
618: } else if (type == NORM_INFINITY) {
619: PetscInt i;
620: PetscReal max = 0.0,tmp;
622: VecGetArrayRead(xin,&xx);
623: for (i=0; i<n; i++) {
624: if ((tmp = PetscAbsScalar(*xx)) > max) max = tmp;
625: /* check special case of tmp == NaN */
626: if (tmp != tmp) {max = tmp; break;}
627: xx++;
628: }
629: VecRestoreArrayRead(xin,&xx);
630: *z = max;
631: } else if (type == NORM_1) {
632: VecGetArrayRead(xin,&xx);
633: PetscStackCallBLAS("BLASasum",*z = BLASasum_(&bn,xx,&one));
634: VecRestoreArrayRead(xin,&xx);
635: PetscLogFlops(PetscMax(n-1.0,0.0));
636: } else if (type == NORM_1_AND_2) {
637: VecNorm_Seq(xin,NORM_1,z);
638: VecNorm_Seq(xin,NORM_2,z+1);
639: }
640: return(0);
641: }
642: #endif
646: PetscErrorCode VecView_Seq_ASCII(Vec xin,PetscViewer viewer)
647: {
648: PetscErrorCode ierr;
649: PetscInt i,n = xin->map->n;
650: const char *name;
651: PetscViewerFormat format;
652: const PetscScalar *xv;
655: VecGetArrayRead(xin,&xv);
656: PetscViewerGetFormat(viewer,&format);
657: if (format == PETSC_VIEWER_ASCII_MATLAB) {
658: PetscObjectGetName((PetscObject)xin,&name);
659: PetscViewerASCIIPrintf(viewer,"%s = [\n",name);
660: for (i=0; i<n; i++) {
661: #if defined(PETSC_USE_COMPLEX)
662: if (PetscImaginaryPart(xv[i]) > 0.0) {
663: PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei\n",(double)PetscRealPart(xv[i]),(double)PetscImaginaryPart(xv[i]));
664: } else if (PetscImaginaryPart(xv[i]) < 0.0) {
665: PetscViewerASCIIPrintf(viewer,"%18.16e - %18.16ei\n",(double)PetscRealPart(xv[i]),-(double)PetscImaginaryPart(xv[i]));
666: } else {
667: PetscViewerASCIIPrintf(viewer,"%18.16e\n",(double)PetscRealPart(xv[i]));
668: }
669: #else
670: PetscViewerASCIIPrintf(viewer,"%18.16e\n",(double) xv[i]);
671: #endif
672: }
673: PetscViewerASCIIPrintf(viewer,"];\n");
674: } else if (format == PETSC_VIEWER_ASCII_SYMMODU) {
675: for (i=0; i<n; i++) {
676: #if defined(PETSC_USE_COMPLEX)
677: PetscViewerASCIIPrintf(viewer,"%18.16e %18.16e\n",(double)PetscRealPart(xv[i]),(double)PetscImaginaryPart(xv[i]));
678: #else
679: PetscViewerASCIIPrintf(viewer,"%18.16e\n",(double)xv[i]);
680: #endif
681: }
682: } else if (format == PETSC_VIEWER_ASCII_VTK || format == PETSC_VIEWER_ASCII_VTK_CELL) {
683: /*
684: state 0: No header has been output
685: state 1: Only POINT_DATA has been output
686: state 2: Only CELL_DATA has been output
687: state 3: Output both, POINT_DATA last
688: state 4: Output both, CELL_DATA last
689: */
690: static PetscInt stateId = -1;
691: int outputState = 0;
692: PetscBool hasState;
693: int doOutput = 0;
694: PetscInt bs, b;
696: if (stateId < 0) {
697: PetscObjectComposedDataRegister(&stateId);
698: }
699: PetscObjectComposedDataGetInt((PetscObject) viewer, stateId, outputState, hasState);
700: if (!hasState) outputState = 0;
701: PetscObjectGetName((PetscObject) xin, &name);
702: VecGetBlockSize(xin, &bs);
703: 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);
704: if (format == PETSC_VIEWER_ASCII_VTK) {
705: if (outputState == 0) {
706: outputState = 1;
707: doOutput = 1;
708: } else if (outputState == 1) doOutput = 0;
709: else if (outputState == 2) {
710: outputState = 3;
711: doOutput = 1;
712: } else if (outputState == 3) doOutput = 0;
713: else if (outputState == 4) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Tried to output POINT_DATA again after intervening CELL_DATA");
715: if (doOutput) {
716: PetscViewerASCIIPrintf(viewer, "POINT_DATA %d\n", n/bs);
717: }
718: } else {
719: if (outputState == 0) {
720: outputState = 2;
721: doOutput = 1;
722: } else if (outputState == 1) {
723: outputState = 4;
724: doOutput = 1;
725: } else if (outputState == 2) doOutput = 0;
726: else if (outputState == 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Tried to output CELL_DATA again after intervening POINT_DATA");
727: else if (outputState == 4) doOutput = 0;
729: if (doOutput) {
730: PetscViewerASCIIPrintf(viewer, "CELL_DATA %d\n", n);
731: }
732: }
733: PetscObjectComposedDataSetInt((PetscObject) viewer, stateId, outputState);
734: if (name) {
735: if (bs == 3) {
736: PetscViewerASCIIPrintf(viewer, "VECTORS %s double\n", name);
737: } else {
738: PetscViewerASCIIPrintf(viewer, "SCALARS %s double %d\n", name, bs);
739: }
740: } else {
741: PetscViewerASCIIPrintf(viewer, "SCALARS scalars double %d\n", bs);
742: }
743: if (bs != 3) {
744: PetscViewerASCIIPrintf(viewer, "LOOKUP_TABLE default\n");
745: }
746: for (i=0; i<n/bs; i++) {
747: for (b=0; b<bs; b++) {
748: if (b > 0) {
749: PetscViewerASCIIPrintf(viewer," ");
750: }
751: #if !defined(PETSC_USE_COMPLEX)
752: PetscViewerASCIIPrintf(viewer,"%g",(double)xv[i*bs+b]);
753: #endif
754: }
755: PetscViewerASCIIPrintf(viewer,"\n");
756: }
757: } else if (format == PETSC_VIEWER_ASCII_VTK_COORDS) {
758: PetscInt bs, b;
760: VecGetBlockSize(xin, &bs);
761: 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);
762: for (i=0; i<n/bs; i++) {
763: for (b=0; b<bs; b++) {
764: if (b > 0) {
765: PetscViewerASCIIPrintf(viewer," ");
766: }
767: #if !defined(PETSC_USE_COMPLEX)
768: PetscViewerASCIIPrintf(viewer,"%g",(double)xv[i*bs+b]);
769: #endif
770: }
771: for (b=bs; b<3; b++) {
772: PetscViewerASCIIPrintf(viewer," 0.0");
773: }
774: PetscViewerASCIIPrintf(viewer,"\n");
775: }
776: } else if (format == PETSC_VIEWER_ASCII_PCICE) {
777: PetscInt bs, b;
779: VecGetBlockSize(xin, &bs);
780: 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);
781: PetscViewerASCIIPrintf(viewer,"%D\n", xin->map->N/bs);
782: for (i=0; i<n/bs; i++) {
783: PetscViewerASCIIPrintf(viewer,"%7D ", i+1);
784: for (b=0; b<bs; b++) {
785: if (b > 0) {
786: PetscViewerASCIIPrintf(viewer," ");
787: }
788: #if !defined(PETSC_USE_COMPLEX)
789: PetscViewerASCIIPrintf(viewer,"% 12.5E",(double)xv[i*bs+b]);
790: #endif
791: }
792: PetscViewerASCIIPrintf(viewer,"\n");
793: }
794: } else if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
795: /* No info */
796: } else {
797: for (i=0; i<n; i++) {
798: if (format == PETSC_VIEWER_ASCII_INDEX) {
799: PetscViewerASCIIPrintf(viewer,"%D: ",i);
800: }
801: #if defined(PETSC_USE_COMPLEX)
802: if (PetscImaginaryPart(xv[i]) > 0.0) {
803: PetscViewerASCIIPrintf(viewer,"%g + %g i\n",(double)PetscRealPart(xv[i]),(double)PetscImaginaryPart(xv[i]));
804: } else if (PetscImaginaryPart(xv[i]) < 0.0) {
805: PetscViewerASCIIPrintf(viewer,"%g - %g i\n",(double)PetscRealPart(xv[i]),-(double)PetscImaginaryPart(xv[i]));
806: } else {
807: PetscViewerASCIIPrintf(viewer,"%g\n",(double)PetscRealPart(xv[i]));
808: }
809: #else
810: PetscViewerASCIIPrintf(viewer,"%g\n",(double)xv[i]);
811: #endif
812: }
813: }
814: PetscViewerFlush(viewer);
815: VecRestoreArrayRead(xin,&xv);
816: return(0);
817: }
819: #include <petscdraw.h>
822: PetscErrorCode VecView_Seq_Draw_LG(Vec xin,PetscViewer v)
823: {
824: PetscErrorCode ierr;
825: PetscInt i,c,bs = PetscAbs(xin->map->bs),n = xin->map->n/bs;
826: PetscDraw win;
827: PetscReal *xx;
828: PetscDrawLG lg;
829: const PetscScalar *xv;
830: PetscReal *yy;
833: PetscMalloc1(n,&xx);
834: PetscMalloc1(n,&yy);
835: VecGetArrayRead(xin,&xv);
836: for (c=0; c<bs; c++) {
837: PetscViewerDrawGetDrawLG(v,c,&lg);
838: PetscDrawLGGetDraw(lg,&win);
839: PetscDrawCheckResizedWindow(win);
840: PetscDrawLGReset(lg);
841: for (i=0; i<n; i++) {
842: xx[i] = (PetscReal) i;
843: yy[i] = PetscRealPart(xv[c + i*bs]);
844: }
845: PetscDrawLGAddPoints(lg,n,&xx,&yy);
846: PetscDrawLGDraw(lg);
847: PetscDrawSynchronizedFlush(win);
848: }
849: VecRestoreArrayRead(xin,&xv);
850: PetscFree(yy);
851: PetscFree(xx);
852: return(0);
853: }
857: PetscErrorCode VecView_Seq_Draw(Vec xin,PetscViewer v)
858: {
859: PetscErrorCode ierr;
860: PetscDraw draw;
861: PetscBool isnull;
862: PetscViewerFormat format;
865: PetscViewerDrawGetDraw(v,0,&draw);
866: PetscDrawIsNull(draw,&isnull); if (isnull) return(0);
868: PetscViewerGetFormat(v,&format);
869: /*
870: Currently it only supports drawing to a line graph */
871: if (format != PETSC_VIEWER_DRAW_LG) {
872: PetscViewerPushFormat(v,PETSC_VIEWER_DRAW_LG);
873: }
874: VecView_Seq_Draw_LG(xin,v);
875: if (format != PETSC_VIEWER_DRAW_LG) {
876: PetscViewerPopFormat(v);
877: }
878: return(0);
879: }
883: PetscErrorCode VecView_Seq_Binary(Vec xin,PetscViewer viewer)
884: {
885: PetscErrorCode ierr;
886: int fdes;
887: PetscInt n = xin->map->n,classid=VEC_FILE_CLASSID;
888: FILE *file;
889: const PetscScalar *xv;
890: #if defined(PETSC_HAVE_MPIIO)
891: PetscBool isMPIIO;
892: #endif
893: PetscBool skipHeader;
894: PetscViewerFormat format;
897: /* Write vector header */
898: PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
899: if (!skipHeader) {
900: PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);
901: PetscViewerBinaryWrite(viewer,&n,1,PETSC_INT,PETSC_FALSE);
902: }
904: /* Write vector contents */
905: #if defined(PETSC_HAVE_MPIIO)
906: PetscViewerBinaryGetMPIIO(viewer,&isMPIIO);
907: if (!isMPIIO) {
908: #endif
909: PetscViewerBinaryGetDescriptor(viewer,&fdes);
910: VecGetArrayRead(xin,&xv);
911: PetscBinaryWrite(fdes,(void*)xv,n,PETSC_SCALAR,PETSC_FALSE);
912: VecRestoreArrayRead(xin,&xv);
913: PetscViewerGetFormat(viewer,&format);
914: if (format == PETSC_VIEWER_BINARY_MATLAB) {
915: MPI_Comm comm;
916: FILE *info;
917: const char *name;
919: PetscObjectGetName((PetscObject)xin,&name);
920: PetscObjectGetComm((PetscObject)viewer,&comm);
921: PetscViewerBinaryGetInfoPointer(viewer,&info);
922: PetscFPrintf(comm,info,"#--- begin code written by PetscViewerBinary for MATLAB format ---#\n");
923: PetscFPrintf(comm,info,"#$$ Set.%s = PetscBinaryRead(fd);\n",name);
924: PetscFPrintf(comm,info,"#--- end code written by PetscViewerBinary for MATLAB format ---#\n\n");
925: }
926: #if defined(PETSC_HAVE_MPIIO)
927: } else {
928: MPI_Offset off;
929: MPI_File mfdes;
930: PetscMPIInt lsize;
932: PetscMPIIntCast(n,&lsize);
933: PetscViewerBinaryGetMPIIODescriptor(viewer,&mfdes);
934: PetscViewerBinaryGetMPIIOOffset(viewer,&off);
935: MPI_File_set_view(mfdes,off,MPIU_SCALAR,MPIU_SCALAR,(char*)"native",MPI_INFO_NULL);
936: VecGetArrayRead(xin,&xv);
937: MPIU_File_write_all(mfdes,(void*)xv,lsize,MPIU_SCALAR,MPI_STATUS_IGNORE);
938: VecRestoreArrayRead(xin,&xv);
939: PetscViewerBinaryAddMPIIOOffset(viewer,n*sizeof(PetscScalar));
940: }
941: #endif
943: PetscViewerBinaryGetInfoPointer(viewer,&file);
944: if (file) {
945: if (((PetscObject)xin)->prefix) {
946: PetscFPrintf(PETSC_COMM_SELF,file,"-%svecload_block_size %D\n",((PetscObject)xin)->prefix,PetscAbs(xin->map->bs));
947: } else {
948: PetscFPrintf(PETSC_COMM_SELF,file,"-vecload_block_size %D\n",PetscAbs(xin->map->bs));
949: }
950: }
951: return(0);
952: }
954: #if defined(PETSC_HAVE_MATLAB_ENGINE)
955: #include <petscmatlab.h>
956: #include <mat.h> /* MATLAB include file */
959: PetscErrorCode VecView_Seq_Matlab(Vec vec,PetscViewer viewer)
960: {
961: PetscErrorCode ierr;
962: PetscInt n;
963: const PetscScalar *array;
966: VecGetLocalSize(vec,&n);
967: PetscObjectName((PetscObject)vec);
968: VecGetArrayRead(vec,&array);
969: PetscViewerMatlabPutArray(viewer,n,1,array,((PetscObject)vec)->name);
970: VecRestoreArrayRead(vec,&array);
971: return(0);
972: }
973: #endif
977: PetscErrorCode VecView_Seq(Vec xin,PetscViewer viewer)
978: {
980: PetscBool isdraw,iascii,issocket,isbinary;
981: #if defined(PETSC_HAVE_MATHEMATICA)
982: PetscBool ismathematica;
983: #endif
984: #if defined(PETSC_HAVE_MATLAB_ENGINE)
985: PetscBool ismatlab;
986: #endif
987: #if defined(PETSC_HAVE_HDF5)
988: PetscBool ishdf5;
989: #endif
992: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
993: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
994: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSOCKET,&issocket);
995: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
996: #if defined(PETSC_HAVE_MATHEMATICA)
997: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERMATHEMATICA,&ismathematica);
998: #endif
999: #if defined(PETSC_HAVE_HDF5)
1000: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
1001: #endif
1002: #if defined(PETSC_HAVE_MATLAB_ENGINE)
1003: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERMATLAB,&ismatlab);
1004: #endif
1006: if (isdraw) {
1007: VecView_Seq_Draw(xin,viewer);
1008: } else if (iascii) {
1009: VecView_Seq_ASCII(xin,viewer);
1010: } else if (isbinary) {
1011: VecView_Seq_Binary(xin,viewer);
1012: #if defined(PETSC_HAVE_MATHEMATICA)
1013: } else if (ismathematica) {
1014: PetscViewerMathematicaPutVector(viewer,xin);
1015: #endif
1016: #if defined(PETSC_HAVE_HDF5)
1017: } else if (ishdf5) {
1018: VecView_MPI_HDF5(xin,viewer); /* Reusing VecView_MPI_HDF5 ... don't want code duplication*/
1019: #endif
1020: #if defined(PETSC_HAVE_MATLAB_ENGINE)
1021: } else if (ismatlab) {
1022: VecView_Seq_Matlab(xin,viewer);
1023: #endif
1024: }
1025: return(0);
1026: }
1030: PetscErrorCode VecGetValues_Seq(Vec xin,PetscInt ni,const PetscInt ix[],PetscScalar y[])
1031: {
1032: const PetscScalar *xx;
1033: PetscInt i;
1034: PetscErrorCode ierr;
1037: VecGetArrayRead(xin,&xx);
1038: for (i=0; i<ni; i++) {
1039: if (xin->stash.ignorenegidx && ix[i] < 0) continue;
1040: #if defined(PETSC_USE_DEBUG)
1041: if (ix[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D cannot be negative",ix[i]);
1042: 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);
1043: #endif
1044: y[i] = xx[ix[i]];
1045: }
1046: VecRestoreArrayRead(xin,&xx);
1047: return(0);
1048: }
1052: PetscErrorCode VecSetValues_Seq(Vec xin,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode m)
1053: {
1054: PetscScalar *xx;
1055: PetscInt i;
1059: VecGetArray(xin,&xx);
1060: if (m == INSERT_VALUES) {
1061: for (i=0; i<ni; i++) {
1062: if (xin->stash.ignorenegidx && ix[i] < 0) continue;
1063: #if defined(PETSC_USE_DEBUG)
1064: if (ix[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D cannot be negative",ix[i]);
1065: 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);
1066: #endif
1067: xx[ix[i]] = y[i];
1068: }
1069: } else {
1070: for (i=0; i<ni; i++) {
1071: if (xin->stash.ignorenegidx && ix[i] < 0) continue;
1072: #if defined(PETSC_USE_DEBUG)
1073: if (ix[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D cannot be negative",ix[i]);
1074: 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);
1075: #endif
1076: xx[ix[i]] += y[i];
1077: }
1078: }
1079: VecRestoreArray(xin,&xx);
1080: return(0);
1081: }
1085: PetscErrorCode VecSetValuesBlocked_Seq(Vec xin,PetscInt ni,const PetscInt ix[],const PetscScalar yin[],InsertMode m)
1086: {
1087: PetscScalar *xx,*y = (PetscScalar*)yin;
1088: PetscInt i,bs,start,j;
1091: /*
1092: For optimization could treat bs = 2, 3, 4, 5 as special cases with loop unrolling
1093: */
1095: VecGetBlockSize(xin,&bs);
1096: VecGetArray(xin,&xx);
1097: if (m == INSERT_VALUES) {
1098: for (i=0; i<ni; i++) {
1099: start = bs*ix[i];
1100: if (start < 0) continue;
1101: #if defined(PETSC_USE_DEBUG)
1102: 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);
1103: #endif
1104: for (j=0; j<bs; j++) xx[start+j] = y[j];
1105: y += bs;
1106: }
1107: } else {
1108: for (i=0; i<ni; i++) {
1109: start = bs*ix[i];
1110: if (start < 0) continue;
1111: #if defined(PETSC_USE_DEBUG)
1112: 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);
1113: #endif
1114: for (j=0; j<bs; j++) xx[start+j] += y[j];
1115: y += bs;
1116: }
1117: }
1118: VecRestoreArray(xin,&xx);
1119: return(0);
1120: }
1124: PetscErrorCode VecDestroy_Seq(Vec v)
1125: {
1126: Vec_Seq *vs = (Vec_Seq*)v->data;
1130: #if defined(PETSC_USE_LOG)
1131: PetscLogObjectState((PetscObject)v,"Length=%D",v->map->n);
1132: #endif
1133: PetscFree(vs->array_allocated);
1134: PetscFree(v->data);
1135: return(0);
1136: }
1140: PetscErrorCode VecSetOption_Seq(Vec v,VecOption op,PetscBool flag)
1141: {
1143: if (op == VEC_IGNORE_NEGATIVE_INDICES) v->stash.ignorenegidx = flag;
1144: return(0);
1145: }
1149: PetscErrorCode VecDuplicate_Seq(Vec win,Vec *V)
1150: {
1154: VecCreate(PetscObjectComm((PetscObject)win),V);
1155: PetscObjectSetPrecision((PetscObject)*V,((PetscObject)win)->precision);
1156: VecSetSizes(*V,win->map->n,win->map->n);
1157: VecSetType(*V,((PetscObject)win)->type_name);
1158: PetscLayoutReference(win->map,&(*V)->map);
1159: PetscObjectListDuplicate(((PetscObject)win)->olist,&((PetscObject)(*V))->olist);
1160: PetscFunctionListDuplicate(((PetscObject)win)->qlist,&((PetscObject)(*V))->qlist);
1162: (*V)->ops->view = win->ops->view;
1163: (*V)->stash.ignorenegidx = win->stash.ignorenegidx;
1164: return(0);
1165: }
1167: static struct _VecOps DvOps = {VecDuplicate_Seq, /* 1 */
1168: VecDuplicateVecs_Default,
1169: VecDestroyVecs_Default,
1170: VecDot_Seq,
1171: VecMDot_Seq,
1172: VecNorm_Seq,
1173: VecTDot_Seq,
1174: VecMTDot_Seq,
1175: VecScale_Seq,
1176: VecCopy_Seq, /* 10 */
1177: VecSet_Seq,
1178: VecSwap_Seq,
1179: VecAXPY_Seq,
1180: VecAXPBY_Seq,
1181: VecMAXPY_Seq,
1182: VecAYPX_Seq,
1183: VecWAXPY_Seq,
1184: VecAXPBYPCZ_Seq,
1185: VecPointwiseMult_Seq,
1186: VecPointwiseDivide_Seq,
1187: VecSetValues_Seq, /* 20 */
1188: 0,0,
1189: 0,
1190: VecGetSize_Seq,
1191: VecGetSize_Seq,
1192: 0,
1193: VecMax_Seq,
1194: VecMin_Seq,
1195: VecSetRandom_Seq,
1196: VecSetOption_Seq, /* 30 */
1197: VecSetValuesBlocked_Seq,
1198: VecDestroy_Seq,
1199: VecView_Seq,
1200: VecPlaceArray_Seq,
1201: VecReplaceArray_Seq,
1202: VecDot_Seq,
1203: VecTDot_Seq,
1204: VecNorm_Seq,
1205: VecMDot_Seq,
1206: VecMTDot_Seq, /* 40 */
1207: VecLoad_Default,
1208: VecReciprocal_Default,
1209: VecConjugate_Seq,
1210: 0,
1211: 0,
1212: VecResetArray_Seq,
1213: 0,
1214: VecMaxPointwiseDivide_Seq,
1215: VecPointwiseMax_Seq,
1216: VecPointwiseMaxAbs_Seq,
1217: VecPointwiseMin_Seq,
1218: VecGetValues_Seq,
1219: 0,
1220: 0,
1221: 0,
1222: 0,
1223: 0,
1224: 0,
1225: VecStrideGather_Default,
1226: VecStrideScatter_Default,
1227: 0,
1228: 0,
1229: 0,
1230: 0,
1231: 0,
1232: VecStrideSubSetGather_Default,
1233: VecStrideSubSetScatter_Default,
1234: 0,
1235: 0
1236: };
1239: /*
1240: This is called by VecCreate_Seq() (i.e. VecCreateSeq()) and VecCreateSeqWithArray()
1241: */
1244: PetscErrorCode VecCreate_Seq_Private(Vec v,const PetscScalar array[])
1245: {
1246: Vec_Seq *s;
1250: PetscNewLog(v,&s);
1251: PetscMemcpy(v->ops,&DvOps,sizeof(DvOps));
1253: v->data = (void*)s;
1254: v->petscnative = PETSC_TRUE;
1255: s->array = (PetscScalar*)array;
1256: s->array_allocated = 0;
1258: PetscLayoutSetUp(v->map);
1259: PetscObjectChangeTypeName((PetscObject)v,VECSEQ);
1260: #if defined(PETSC_HAVE_MATLAB_ENGINE)
1261: PetscObjectComposeFunction((PetscObject)v,"PetscMatlabEnginePut_C",VecMatlabEnginePut_Default);
1262: PetscObjectComposeFunction((PetscObject)v,"PetscMatlabEngineGet_C",VecMatlabEngineGet_Default);
1263: #endif
1264: #if defined(PETSC_USE_MIXED_PRECISION)
1265: ((PetscObject)v)->precision = (PetscPrecision)sizeof(PetscReal);
1266: #endif
1267: return(0);
1268: }
1272: /*@C
1273: VecCreateSeqWithArray - Creates a standard,sequential array-style vector,
1274: where the user provides the array space to store the vector values.
1276: Collective on MPI_Comm
1278: Input Parameter:
1279: + comm - the communicator, should be PETSC_COMM_SELF
1280: . bs - the block size
1281: . n - the vector length
1282: - array - memory where the vector elements are to be stored.
1284: Output Parameter:
1285: . V - the vector
1287: Notes:
1288: Use VecDuplicate() or VecDuplicateVecs() to form additional vectors of the
1289: same type as an existing vector.
1291: If the user-provided array is NULL, then VecPlaceArray() can be used
1292: at a later stage to SET the array for storing the vector values.
1294: PETSc does NOT free the array when the vector is destroyed via VecDestroy().
1295: The user should not free the array until the vector is destroyed.
1297: Level: intermediate
1299: Concepts: vectors^creating with array
1301: .seealso: VecCreateMPIWithArray(), VecCreate(), VecDuplicate(), VecDuplicateVecs(),
1302: VecCreateGhost(), VecCreateSeq(), VecPlaceArray()
1303: @*/
1304: PetscErrorCode VecCreateSeqWithArray(MPI_Comm comm,PetscInt bs,PetscInt n,const PetscScalar array[],Vec *V)
1305: {
1307: PetscMPIInt size;
1310: VecCreate(comm,V);
1311: VecSetSizes(*V,n,n);
1312: VecSetBlockSize(*V,bs);
1313: MPI_Comm_size(comm,&size);
1314: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Cannot create VECSEQ on more than one process");
1315: VecCreate_Seq_Private(*V,array);
1316: return(0);
1317: }