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: }