Actual source code: bvec2.c
1: /*$Id: bvec2.c,v 1.202 2001/09/12 03:26:24 bsmith Exp $*/
2: /*
3: Implements the sequential vectors.
4: */
6: #include src/vec/vecimpl.h
7: #include src/vec/impls/dvecimpl.h
8: #include src/inline/dot.h
9: #include petscblaslapack.h
10: #if defined(PETSC_HAVE_AMS)
11: EXTERN int PetscViewerAMSGetAMSComm(PetscViewer,AMS_Comm *);
12: #endif
14: #undef __FUNCT__
16: int VecNorm_Seq(Vec xin,NormType type,PetscReal* z)
17: {
18: PetscScalar *xx;
19: int n=xin->n,ierr,one = 1;
22: if (type == NORM_2) {
23: VecGetArrayFast(xin,&xx);
24: /*
25: This is because the Fortran BLAS 1 Norm is very slow!
26: */
27: #if defined(PETSC_HAVE_SLOW_NRM2)
28: #if defined(PETSC_USE_FORTRAN_KERNEL_NORM)
29: fortrannormsqr_(xx,&n,z);
30: *z = sqrt(*z);
31: #elif defined(PETSC_USE_UNROLLED_NORM)
32: {
33: PetscReal work = 0.0;
34: switch (n & 0x3) {
35: case 3: work += PetscRealPart(xx[0]*PetscConj(xx[0])); xx++;
36: case 2: work += PetscRealPart(xx[0]*PetscConj(xx[0])); xx++;
37: case 1: work += PetscRealPart(xx[0]*PetscConj(xx[0])); xx++; n -= 4;
38: }
39: while (n>0) {
40: work += PetscRealPart(xx[0]*PetscConj(xx[0])+xx[1]*PetscConj(xx[1])+
41: xx[2]*PetscConj(xx[2])+xx[3]*PetscConj(xx[3]));
42: xx += 4; n -= 4;
43: }
44: *z = sqrt(work);}
45: #else
46: {
47: int i;
48: PetscScalar sum=0.0;
49: for (i=0; i<n; i++) {
50: sum += (xx[i])*(PetscConj(xx[i]));
51: }
52: *z = sqrt(PetscRealPart(sum));
53: }
54: #endif
55: #else
56: *z = BLnrm2_(&n,xx,&one);
57: #endif
58: VecRestoreArrayFast(xin,&xx);
59: PetscLogFlops(2*n-1);
60: } else if (type == NORM_INFINITY) {
61: int i;
62: PetscReal max = 0.0,tmp;
64: VecGetArrayFast(xin,&xx);
65: for (i=0; i<n; i++) {
66: if ((tmp = PetscAbsScalar(*xx)) > max) max = tmp;
67: /* check special case of tmp == NaN */
68: if (tmp != tmp) {max = tmp; break;}
69: xx++;
70: }
71: VecRestoreArrayFast(xin,&xx);
72: *z = max;
73: } else if (type == NORM_1) {
74: VecGetArrayFast(xin,&xx);
75: *z = BLasum_(&n,xx,&one);
76: VecRestoreArrayFast(xin,&xx);
77: PetscLogFlops(n-1);
78: } else if (type == NORM_1_AND_2) {
79: VecNorm_Seq(xin,NORM_1,z);
80: VecNorm_Seq(xin,NORM_2,z+1);
81: }
82: return(0);
83: }
85: #include petscviewer.h
86: #include petscsys.h
88: #undef __FUNCT__
90: int VecView_Seq_File(Vec xin,PetscViewer viewer)
91: {
92: Vec_Seq *x = (Vec_Seq *)xin->data;
93: int i,n = xin->n,ierr;
94: char *name;
95: PetscViewerFormat format;
98: PetscViewerGetFormat(viewer,&format);
99: if (format == PETSC_VIEWER_ASCII_MATLAB) {
100: PetscObjectGetName((PetscObject)xin,&name);
101: PetscViewerASCIIPrintf(viewer,"%s = [n",name);
102: for (i=0; i<n; i++) {
103: #if defined(PETSC_USE_COMPLEX)
104: if (PetscImaginaryPart(x->array[i]) > 0.0) {
105: PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ein",PetscRealPart(x->array[i]),PetscImaginaryPart(x->array[i]));
106: } else if (PetscImaginaryPart(x->array[i]) < 0.0) {
107: PetscViewerASCIIPrintf(viewer,"%18.16e - %18.16ein",PetscRealPart(x->array[i]),-PetscImaginaryPart(x->array[i]));
108: } else {
109: PetscViewerASCIIPrintf(viewer,"%18.16en",PetscRealPart(x->array[i]));
110: }
111: #else
112: PetscViewerASCIIPrintf(viewer,"%18.16en",x->array[i]);
113: #endif
114: }
115: PetscViewerASCIIPrintf(viewer,"];n");
116: } else if (format == PETSC_VIEWER_ASCII_SYMMODU) {
117: for (i=0; i<n; i++) {
118: #if defined(PETSC_USE_COMPLEX)
119: PetscViewerASCIIPrintf(viewer,"%18.16e %18.16en",PetscRealPart(x->array[i]),PetscImaginaryPart(x->array[i]));
120: #else
121: PetscViewerASCIIPrintf(viewer,"%18.16en",x->array[i]);
122: #endif
123: }
124: } else {
125: for (i=0; i<n; i++) {
126: if (format == PETSC_VIEWER_ASCII_INDEX) {
127: PetscViewerASCIIPrintf(viewer,"%d: ",i);
128: }
129: #if defined(PETSC_USE_COMPLEX)
130: if (PetscImaginaryPart(x->array[i]) > 0.0) {
131: PetscViewerASCIIPrintf(viewer,"%g + %g in",PetscRealPart(x->array[i]),PetscImaginaryPart(x->array[i]));
132: } else if (PetscImaginaryPart(x->array[i]) < 0.0) {
133: PetscViewerASCIIPrintf(viewer,"%g - %g in",PetscRealPart(x->array[i]),-PetscImaginaryPart(x->array[i]));
134: } else {
135: PetscViewerASCIIPrintf(viewer,"%gn",PetscRealPart(x->array[i]));
136: }
137: #else
138: PetscViewerASCIIPrintf(viewer,"%gn",x->array[i]);
139: #endif
140: }
141: }
142: PetscViewerFlush(viewer);
143: return(0);
144: }
146: #undef __FUNCT__
148: static int VecView_Seq_Draw_LG(Vec xin,PetscViewer v)
149: {
150: Vec_Seq *x = (Vec_Seq *)xin->data;
151: int i,n = xin->n,ierr;
152: PetscDraw win;
153: PetscReal *xx;
154: PetscDrawLG lg;
157: PetscViewerDrawGetDrawLG(v,0,&lg);
158: PetscDrawLGGetDraw(lg,&win);
159: PetscDrawCheckResizedWindow(win);
160: PetscDrawLGReset(lg);
161: PetscMalloc((n+1)*sizeof(PetscReal),&xx);
162: for (i=0; i<n; i++) {
163: xx[i] = (PetscReal) i;
164: }
165: #if !defined(PETSC_USE_COMPLEX)
166: PetscDrawLGAddPoints(lg,n,&xx,&x->array);
167: #else
168: {
169: PetscReal *yy;
170: PetscMalloc((n+1)*sizeof(PetscReal),&yy);
171: for (i=0; i<n; i++) {
172: yy[i] = PetscRealPart(x->array[i]);
173: }
174: PetscDrawLGAddPoints(lg,n,&xx,&yy);
175: PetscFree(yy);
176: }
177: #endif
178: PetscFree(xx);
179: PetscDrawLGDraw(lg);
180: PetscDrawSynchronizedFlush(win);
181: return(0);
182: }
184: #undef __FUNCT__
186: static int VecView_Seq_Draw(Vec xin,PetscViewer v)
187: {
188: int ierr;
189: PetscDraw draw;
190: PetscTruth isnull;
191: PetscViewerFormat format;
194: PetscViewerDrawGetDraw(v,0,&draw);
195: PetscDrawIsNull(draw,&isnull); if (isnull) return(0);
196:
197: PetscViewerGetFormat(v,&format);
198: /*
199: Currently it only supports drawing to a line graph */
200: if (format != PETSC_VIEWER_DRAW_LG) {
201: PetscViewerPushFormat(v,PETSC_VIEWER_DRAW_LG);
202: }
203: VecView_Seq_Draw_LG(xin,v);
204: if (format != PETSC_VIEWER_DRAW_LG) {
205: PetscViewerPopFormat(v);
206: }
208: return(0);
209: }
211: #undef __FUNCT__
213: static int VecView_Seq_Binary(Vec xin,PetscViewer viewer)
214: {
215: Vec_Seq *x = (Vec_Seq *)xin->data;
216: int ierr,fdes,n = xin->n,cookie=VEC_FILE_COOKIE;
217: FILE *file;
220: ierr = PetscViewerBinaryGetDescriptor(viewer,&fdes);
221: /* Write vector header */
222: PetscBinaryWrite(fdes,&cookie,1,PETSC_INT,0);
223: PetscBinaryWrite(fdes,&n,1,PETSC_INT,0);
225: /* Write vector contents */
226: PetscBinaryWrite(fdes,x->array,n,PETSC_SCALAR,0);
228: PetscViewerBinaryGetInfoPointer(viewer,&file);
229: if (file && xin->bs > 1) {
230: fprintf(file,"-vecload_block_size %dn",xin->bs);
231: }
232: return(0);
233: }
236: #undef __FUNCT__
238: int VecView_Seq(Vec xin,PetscViewer viewer)
239: {
240: Vec_Seq *x = (Vec_Seq *)xin->data;
241: int ierr;
242: PetscTruth isdraw,isascii,issocket,isbinary,ismathematica;
245: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);
246: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);
247: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_SOCKET,&issocket);
248: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_BINARY,&isbinary);
249: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_MATHEMATICA,&ismathematica);
250: if (isdraw){
251: VecView_Seq_Draw(xin,viewer);
252: } else if (isascii){
253: VecView_Seq_File(xin,viewer);
254: } else if (issocket) {
255: PetscViewerSocketPutScalar(viewer,xin->n,1,x->array);
256: } else if (isbinary) {
257: VecView_Seq_Binary(xin,viewer);
258: } else if (ismathematica) {
259: PetscViewerMathematicaPutVector(viewer,xin);
260: } else {
261: SETERRQ1(1,"Viewer type %s not supported by this vector object",((PetscObject)viewer)->type_name);
262: }
263: return(0);
264: }
266: #undef __FUNCT__
268: int VecSetValues_Seq(Vec xin,int ni,const int ix[],const PetscScalar y[],InsertMode m)
269: {
270: Vec_Seq *x = (Vec_Seq *)xin->data;
271: PetscScalar *xx = x->array;
272: int i;
275: if (m == INSERT_VALUES) {
276: for (i=0; i<ni; i++) {
277: if (ix[i] < 0) continue;
278: #if defined(PETSC_USE_BOPT_g)
279: if (ix[i] >= xin->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %d maximum %d",ix[i],xin->n);
280: #endif
281: xx[ix[i]] = y[i];
282: }
283: } else {
284: for (i=0; i<ni; i++) {
285: if (ix[i] < 0) continue;
286: #if defined(PETSC_USE_BOPT_g)
287: if (ix[i] >= xin->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %d maximum %d",ix[i],xin->n);
288: #endif
289: xx[ix[i]] += y[i];
290: }
291: }
292: return(0);
293: }
295: #undef __FUNCT__
297: int VecSetValuesBlocked_Seq(Vec xin,int ni,const int ix[],const PetscScalar yin[],InsertMode m)
298: {
299: Vec_Seq *x = (Vec_Seq *)xin->data;
300: PetscScalar *xx = x->array,*y = (PetscScalar*)yin;
301: int i,bs = xin->bs,start,j;
303: /*
304: For optimization could treat bs = 2, 3, 4, 5 as special cases with loop unrolling
305: */
307: if (m == INSERT_VALUES) {
308: for (i=0; i<ni; i++) {
309: start = bs*ix[i];
310: if (start < 0) continue;
311: #if defined(PETSC_USE_BOPT_g)
312: if (start >= xin->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %d maximum %d",start,xin->n);
313: #endif
314: for (j=0; j<bs; j++) {
315: xx[start+j] = y[j];
316: }
317: y += bs;
318: }
319: } else {
320: for (i=0; i<ni; i++) {
321: start = bs*ix[i];
322: if (start < 0) continue;
323: #if defined(PETSC_USE_BOPT_g)
324: if (start >= xin->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %d maximum %d",start,xin->n);
325: #endif
326: for (j=0; j<bs; j++) {
327: xx[start+j] += y[j];
328: }
329: y += bs;
330: }
331: }
332: return(0);
333: }
336: #undef __FUNCT__
338: int VecDestroy_Seq(Vec v)
339: {
340: Vec_Seq *vs = (Vec_Seq*)v->data;
341: int ierr;
345: /* if memory was published with AMS then destroy it */
346: PetscObjectDepublish(v);
348: #if defined(PETSC_USE_LOG)
349: PetscLogObjectState((PetscObject)v,"Length=%d",v->n);
350: #endif
351: if (vs->array_allocated) {PetscFree(vs->array_allocated);}
352: PetscFree(vs);
354: return(0);
355: }
357: #undef __FUNCT__
359: static int VecPublish_Seq(PetscObject obj)
360: {
361: #if defined(PETSC_HAVE_AMS)
362: Vec v = (Vec) obj;
363: Vec_Seq *s = (Vec_Seq*)v->data;
364: int ierr,(*f)(AMS_Memory,char *,Vec);
365: #endif
369: #if defined(PETSC_HAVE_AMS)
370: /* if it is already published then return */
371: if (v->amem >=0) return(0);
373: /* if array in vector was not allocated (for example PCSetUp_BJacobi_Singleblock()) then
374: cannot AMS publish the object*/
375: if (!s->array) return(0);
377: PetscObjectPublishBaseBegin(obj);
378: AMS_Memory_add_field((AMS_Memory)v->amem,"values",s->array,v->n,AMS_DOUBLE,AMS_READ,
379: AMS_DISTRIBUTED,AMS_REDUCT_UNDEF);
381: /* if the vector knows its "layout" let it set it*/
382: PetscObjectQueryFunction(obj,"AMSSetFieldBlock_C",(void (**)(void))&f);
383: if (f) {
384: (*f)((AMS_Memory)v->amem,"values",v);
385: }
386: PetscObjectPublishBaseEnd(obj);
387: #endif
389: return(0);
390: }
392: static struct _VecOps DvOps = {VecDuplicate_Seq,
393: VecDuplicateVecs_Default,
394: VecDestroyVecs_Default,
395: VecDot_Seq,
396: VecMDot_Seq,
397: VecNorm_Seq,
398: VecTDot_Seq,
399: VecMTDot_Seq,
400: VecScale_Seq,
401: VecCopy_Seq,
402: VecSet_Seq,
403: VecSwap_Seq,
404: VecAXPY_Seq,
405: VecAXPBY_Seq,
406: VecMAXPY_Seq,
407: VecAYPX_Seq,
408: VecWAXPY_Seq,
409: VecPointwiseMult_Seq,
410: VecPointwiseDivide_Seq,
411: VecSetValues_Seq,0,0,
412: VecGetArray_Seq,
413: VecGetSize_Seq,
414: VecGetSize_Seq,
415: VecRestoreArray_Seq,
416: VecMax_Seq,
417: VecMin_Seq,
418: VecSetRandom_Seq,0,
419: VecSetValuesBlocked_Seq,
420: VecDestroy_Seq,
421: VecView_Seq,
422: VecPlaceArray_Seq,
423: VecReplaceArray_Seq,
424: VecDot_Seq,
425: VecTDot_Seq,
426: VecNorm_Seq,
427: VecLoadIntoVector_Default,
428: VecReciprocal_Default,
429: 0, /* VecViewNative */
430: VecConjugate_Seq,
431: 0,
432: 0,
433: VecResetArray_Seq,
434: 0,
435: VecMaxPointwiseDivide_Seq};
438: /*
439: This is called by VecCreate_Seq() (i.e. VecCreateSeq()) and VecCreateSeqWithArray()
440: */
441: #undef __FUNCT__
443: static int VecCreate_Seq_Private(Vec v,const PetscScalar array[])
444: {
445: Vec_Seq *s;
446: int ierr;
449: PetscMemcpy(v->ops,&DvOps,sizeof(DvOps));
450: PetscNew(Vec_Seq,&s);
451: PetscMemzero(s,sizeof(Vec_Seq));
452: v->data = (void*)s;
453: v->bops->publish = VecPublish_Seq;
454: v->n = PetscMax(v->n,v->N);
455: v->N = PetscMax(v->n,v->N);
456: v->bs = -1;
457: v->petscnative = PETSC_TRUE;
458: s->array = (PetscScalar *)array;
459: s->array_allocated = 0;
460: if (!v->map) {
461: PetscMapCreateMPI(v->comm,v->n,v->N,&v->map);
462: }
463: PetscObjectChangeTypeName((PetscObject)v,VECSEQ);
464: #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE)
465: PetscObjectComposeFunctionDynamic((PetscObject)v,"PetscMatlabEnginePut_C","VecMatlabEnginePut_Default",VecMatlabEnginePut_Default);
466: PetscObjectComposeFunctionDynamic((PetscObject)v,"PetscMatlabEngineGet_C","VecMatlabEngineGet_Default",VecMatlabEngineGet_Default);
467: #endif
468: PetscPublishAll(v);
469: return(0);
470: }
472: #undef __FUNCT__
474: /*@C
475: VecCreateSeqWithArray - Creates a standard,sequential array-style vector,
476: where the user provides the array space to store the vector values.
478: Collective on MPI_Comm
480: Input Parameter:
481: + comm - the communicator, should be PETSC_COMM_SELF
482: . n - the vector length
483: - array - memory where the vector elements are to be stored.
485: Output Parameter:
486: . V - the vector
488: Notes:
489: Use VecDuplicate() or VecDuplicateVecs() to form additional vectors of the
490: same type as an existing vector.
492: If the user-provided array is PETSC_NULL, then VecPlaceArray() can be used
493: at a later stage to SET the array for storing the vector values.
495: PETSc does NOT free the array when the vector is destroyed via VecDestroy().
496: The user should not free the array until the vector is destroyed.
498: Level: intermediate
500: Concepts: vectors^creating with array
502: .seealso: VecCreateMPIWithArray(), VecCreate(), VecDuplicate(), VecDuplicateVecs(),
503: VecCreateGhost(), VecCreateSeq(), VecPlaceArray()
504: @*/
505: int VecCreateSeqWithArray(MPI_Comm comm,int n,const PetscScalar array[],Vec *V)
506: {
507: int ierr;
510: VecCreate(comm,V);
511: VecSetSizes(*V,n,n);
512: VecCreate_Seq_Private(*V,array);
513: return(0);
514: }
516: EXTERN_C_BEGIN
517: #undef __FUNCT__
519: int VecCreate_Seq(Vec V)
520: {
521: Vec_Seq *s;
522: PetscScalar *array;
523: int ierr,n = PetscMax(V->n,V->N);
526: PetscMalloc((n+1)*sizeof(PetscScalar),&array);
527: PetscMemzero(array,n*sizeof(PetscScalar));
528: VecCreate_Seq_Private(V,array);
529: s = (Vec_Seq*)V->data;
530: s->array_allocated = array;
531: VecSetSerializeType(V,VEC_SER_SEQ_BINARY);
532: return(0);
533: }
535: #undef __FUNCT__
537: int VecSerialize_Seq(MPI_Comm comm, Vec *vec, PetscViewer viewer, PetscTruth store)
538: {
539: Vec v;
540: Vec_Seq *x;
541: PetscScalar *array;
542: int fd;
543: int vars;
544: int ierr;
547: PetscViewerBinaryGetDescriptor(viewer, &fd);
548: if (store) {
549: v = *vec;
550: x = (Vec_Seq *) v->data;
551: PetscBinaryWrite(fd, &v->n, 1, PETSC_INT, 0);
552: PetscBinaryWrite(fd, x->array, v->n, PETSC_SCALAR, 0);
553: } else {
554: PetscBinaryRead(fd, &vars, 1, PETSC_INT);
555: VecCreate(comm, &v);
556: VecSetSizes(v, vars, vars);
557: PetscMalloc(vars * sizeof(PetscScalar), &array);
558: PetscBinaryRead(fd, array, vars, PETSC_SCALAR);
559: VecCreate_Seq_Private(v, array);
560: ((Vec_Seq *) v->data)->array_allocated = array;
562: VecAssemblyBegin(v);
563: VecAssemblyEnd(v);
564: *vec = v;
565: }
567: return(0);
568: }
569: EXTERN_C_END
572: #undef __FUNCT__
574: int VecDuplicate_Seq(Vec win,Vec *V)
575: {
576: int ierr;
579: VecCreateSeq(win->comm,win->n,V);
580: if (win->mapping) {
581: (*V)->mapping = win->mapping;
582: PetscObjectReference((PetscObject)win->mapping);
583: }
584: if (win->bmapping) {
585: (*V)->bmapping = win->bmapping;
586: PetscObjectReference((PetscObject)win->bmapping);
587: }
588: (*V)->bs = win->bs;
589: PetscOListDuplicate(win->olist,&(*V)->olist);
590: PetscFListDuplicate(win->qlist,&(*V)->qlist);
591: return(0);
592: }