Actual source code: vector.c

  1: /*$Id: vector.c,v 1.238 2001/09/11 16:31:48 bsmith Exp $*/
  2: /*
  3:      Provides the interface functions for all vector operations.
  4:    These are the vector functions the user calls.
  5: */
 6:  #include src/vec/vecimpl.h

  8: /* Logging support */
  9: int VEC_COOKIE;
 10: int VEC_View, VEC_Max, VEC_Min, VEC_DotBarrier, VEC_Dot, VEC_MDotBarrier, VEC_MDot, VEC_TDot, VEC_MTDot, VEC_NormBarrier;
 11: int VEC_Norm, VEC_Scale, VEC_Copy, VEC_Set, VEC_AXPY, VEC_AYPX, VEC_WAXPY, VEC_MAXPY, VEC_Swap, VEC_AssemblyBegin;
 12: int VEC_AssemblyEnd, VEC_PointwiseMult, VEC_SetValues, VEC_Load, VEC_ScatterBarrier, VEC_ScatterBegin, VEC_ScatterEnd;
 13: int VEC_SetRandom, VEC_ReduceArithmetic, VEC_ReduceBarrier, VEC_ReduceCommunication;

 15: #undef __FUNCT__  
 17: /*
 18:   VecSetTypeFromOptions_Private - Sets the type of vector from user options. Defaults to a PETSc sequential vector on one
 19:   processor and a PETSc MPI vector on more than one processor.

 21:   Collective on Vec

 23:   Input Parameter:
 24: . vec - The vector

 26:   Level: intermediate

 28: .keywords: Vec, set, options, database, type
 29: .seealso: VecSetFromOptions(), VecSetType()
 30: */
 31: static int VecSetTypeFromOptions_Private(Vec vec)
 32: {
 33:   PetscTruth opt;
 34:   char      *defaultType;
 35:   char       typeName[256];
 36:   int        numProcs;
 37:   int        ierr;

 40:   if (vec->type_name != PETSC_NULL) {
 41:     defaultType = vec->type_name;
 42:   } else {
 43:     MPI_Comm_size(vec->comm, &numProcs);
 44:     if (numProcs > 1) {
 45:       defaultType = VECMPI;
 46:     } else {
 47:       defaultType = VECSEQ;
 48:     }
 49:   }

 51:   if (!VecRegisterAllCalled) {
 52:     VecRegisterAll(PETSC_NULL);
 53:   }
 54:   PetscOptionsList("-vec_type", "Vector type"," VecSetType", VecList, defaultType, typeName, 256, &opt);
 55: 
 56:   if (opt == PETSC_TRUE) {
 57:     VecSetType(vec, typeName);
 58:   } else {
 59:     VecSetType(vec, defaultType);
 60:   }
 61:   return(0);
 62: }

 64: #undef __FUNCT__  
 66: /*@C
 67:   VecSetFromOptions - Configures the vector from the options database.

 69:   Collective on Vec

 71:   Input Parameter:
 72: . vec - The vector

 74:   Notes:  To see all options, run your program with the -help option, or consult the users manual.
 75:           Must be called after VecCreate() but before the vector is used.

 77:   Level: beginner

 79:   Concepts: vectors^setting options
 80:   Concepts: vectors^setting type

 82: .keywords: Vec, set, options, database
 83: .seealso: VecCreate(), VecPrintHelp(), VechSetOptionsPrefix()
 84: @*/
 85: int VecSetFromOptions(Vec vec)
 86: {
 87:   PetscTruth opt;
 88:   int        ierr;


 93:   PetscOptionsBegin(vec->comm, vec->prefix, "Vector options", "Vec");

 95:   /* Handle generic vector options */
 96:   PetscOptionsHasName(PETSC_NULL, "-help", &opt);
 97:   if (opt == PETSC_TRUE) {
 98:     VecPrintHelp(vec);
 99:   }

101:   /* Handle vector type options */
102:   VecSetTypeFromOptions_Private(vec);

104:   /* Handle specific vector options */
105:   if (vec->ops->setfromoptions != PETSC_NULL) {
106:     (*vec->ops->setfromoptions)(vec);
107:   }
108:   PetscOptionsEnd();

110: #if defined(__cplusplus) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) && defined(PETSC_HAVE_CXX_NAMESPACE)
111:   VecESISetFromOptions(vec);
112: #endif

114:   VecViewFromOptions(vec, vec->name);
115:   return(0);
116: }

118: #undef __FUNCT__  
120: /*@
121:   VecPrintHelp - Prints some options for the Vec.

123:   Input Parameter:
124: . vec - The vector

126:   Options Database Keys:
127: $  -help, -h

129:   Level: intermediate

131: .keywords: Vec, help
132: .seealso: VecSetFromOptions()
133: @*/
134: int VecPrintHelp(Vec vec)
135: {
138:   return(0);
139: }

141: #undef __FUNCT__  
143: /*@
144:   VecSetSizes - Sets the local and global sizes, and checks to determine compatibility

146:   Collective on Vec

148:   Input Parameters:
149: + v - the vector
150: . n - the local size (or PETSC_DECIDE to have it set)
151: - N - the global size (or PETSC_DECIDE)

153:   Notes:
154:   n and N cannot be both PETSC_DECIDE
155:   If one processor calls this with N of PETSC_DECIDE then all processors must, otherwise the program will hang.

157:   Level: intermediate

159: .seealso: VecGetSize(), PetscSplitOwnership()
160: @*/
161: int VecSetSizes(Vec v, int n, int N)
162: {

167:   v->n = n;
168:   v->N = N;
169:   PetscSplitOwnership(v->comm, &v->n, &v->N);
170:   return(0);
171: }

173: #undef __FUNCT__  
175: /*@
176:    VecSetBlockSize - Sets the blocksize for future calls to VecSetValuesBlocked()
177:    and VecSetValuesBlockedLocal().

179:    Collective on Vec

181:    Input Parameter:
182: +  v - the vector
183: -  bs - the blocksize

185:    Notes:
186:    All vectors obtained by VecDuplicate() inherit the same blocksize.

188:    Level: advanced

190: .seealso: VecSetValuesBlocked(), VecSetLocalToGlobalMappingBlocked(), VecGetBlockSize()

192:   Concepts: block size^vectors
193: @*/
194: int VecSetBlockSize(Vec v,int bs)
195: {
198:   if (bs <= 0) bs = 1;
199:   if (bs == v->bs) return(0);
200:   if (v->bs != -1) SETERRQ2(PETSC_ERR_ARG_WRONGSTATE,"Cannot reset blocksize. Current size %d new %d",v->bs,bs);
201:   if (v->N != -1 && v->N % bs) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Vector length not divisible by blocksize %d %d",v->N,bs);
202:   if (v->n != -1 && v->n % bs) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Local vector length not divisible by blocksize %d %dn
203:    Try setting blocksize before setting the vector type",v->n,bs);
204: 
205:   v->bs        = bs;
206:   v->bstash.bs = bs; /* use the same blocksize for the vec's block-stash */
207:   return(0);
208: }

210: #undef __FUNCT__  
212: /*@
213:    VecGetBlockSize - Gets the blocksize for the vector, i.e. what is used for VecSetValuesBlocked()
214:    and VecSetValuesBlockedLocal().

216:    Collective on Vec

218:    Input Parameter:
219: .  v - the vector

221:    Output Parameter:
222: .  bs - the blocksize

224:    Notes:
225:    All vectors obtained by VecDuplicate() inherit the same blocksize.

227:    Level: advanced

229: .seealso: VecSetValuesBlocked(), VecSetLocalToGlobalMappingBlocked(), VecSetBlockSize()

231:    Concepts: vector^block size
232:    Concepts: block^vector

234: @*/
235: int VecGetBlockSize(Vec v,int *bs)
236: {
239:   *bs = v->bs;
240:   return(0);
241: }

243: #undef __FUNCT__  
245: /*@
246:    VecValid - Checks whether a vector object is valid.

248:    Not Collective

250:    Input Parameter:
251: .  v - the object to check

253:    Output Parameter:
254:    flg - flag indicating vector status, either
255:    PETSC_TRUE if vector is valid, or PETSC_FALSE otherwise.

257:    Level: developer

259: @*/
260: int VecValid(Vec v,PetscTruth *flg)
261: {
264:   if (!v)                           *flg = PETSC_FALSE;
265:   else if (v->cookie != VEC_COOKIE) *flg = PETSC_FALSE;
266:   else                              *flg = PETSC_TRUE;
267:   return(0);
268: }

270: #undef __FUNCT__  
272: /*@
273:    VecDot - Computes the vector dot product.

275:    Collective on Vec

277:    Input Parameters:
278: .  x, y - the vectors

280:    Output Parameter:
281: .  alpha - the dot product

283:    Performance Issues:
284: +    per-processor memory bandwidth
285: .    interprocessor latency
286: -    work load inbalance that causes certain processes to arrive much earlier than
287:      others

289:    Notes for Users of Complex Numbers:
290:    For complex vectors, VecDot() computes 
291: $     val = (x,y) = y^H x,
292:    where y^H denotes the conjugate transpose of y.

294:    Use VecTDot() for the indefinite form
295: $     val = (x,y) = y^T x,
296:    where y^T denotes the transpose of y.

298:    Level: intermediate

300:    Concepts: inner product
301:    Concepts: vector^inner product

303: .seealso: VecMDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd()
304: @*/
305: int VecDot(Vec x,Vec y,PetscScalar *val)
306: {

317:   if (x->N != y->N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
318:   if (x->n != y->n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

320:   PetscLogEventBarrierBegin(VEC_DotBarrier,x,y,0,0,x->comm);
321:   (*x->ops->dot)(x,y,val);
322:   PetscLogEventBarrierEnd(VEC_DotBarrier,x,y,0,0,x->comm);
323:   /*
324:      The next block is for incremental debugging
325:   */
326:   if (PetscCompare) {
327:     int flag;
328:     MPI_Comm_compare(PETSC_COMM_WORLD,x->comm,&flag);
329:     if (flag != MPI_UNEQUAL) {
330:       PetscCompareScalar(*val);
331:     }
332:   }
333:   return(0);
334: }

336: #undef __FUNCT__  
338: /*@
339:    VecNorm  - Computes the vector norm.

341:    Collective on Vec

343:    Input Parameters:
344: +  x - the vector
345: -  type - one of NORM_1, NORM_2, NORM_INFINITY.  Also available
346:           NORM_1_AND_2, which computes both norms and stores them
347:           in a two element array.

349:    Output Parameter:
350: .  val - the norm 

352:    Notes:
353: $     NORM_1 denotes sum_i |x_i|
354: $     NORM_2 denotes sqrt(sum_i (x_i)^2)
355: $     NORM_INFINITY denotes max_i |x_i|

357:    Level: intermediate

359:    Performance Issues:
360: +    per-processor memory bandwidth
361: .    interprocessor latency
362: -    work load inbalance that causes certain processes to arrive much earlier than
363:      others

365:    Compile Option:
366:    PETSC_HAVE_SLOW_NRM2 will cause a C (loop unrolled) version of the norm to be used, rather
367:  than the BLAS. This should probably only be used when one is using the FORTRAN BLAS routines 
368:  (as opposed to vendor provided) because the FORTRAN BLAS NRM2() routine is very slow. 

370:    Concepts: norm
371:    Concepts: vector^norm

373: .seealso: VecDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), 
374:           VecNormBegin(), VecNormEnd()

376: @*/
377: int VecNorm(Vec x,NormType type,PetscReal *val)
378: {

384:   PetscLogEventBarrierBegin(VEC_NormBarrier,x,0,0,0,x->comm);
385:   (*x->ops->norm)(x,type,val);
386:   PetscLogEventBarrierEnd(VEC_NormBarrier,x,0,0,0,x->comm);
387:   /*
388:      The next block is for incremental debugging
389:   */
390:   if (PetscCompare) {
391:     int flag;
392:     MPI_Comm_compare(PETSC_COMM_WORLD,x->comm,&flag);
393:     if (flag != MPI_UNEQUAL) {
394:       PetscCompareDouble(*val);
395:     }
396:   }
397:   return(0);
398: }

400: #undef __FUNCT__  
402: /*@C
403:    VecMax - Determines the maximum vector component and its location.

405:    Collective on Vec

407:    Input Parameter:
408: .  x - the vector

410:    Output Parameters:
411: +  val - the maximum component
412: -  p - the location of val

414:    Notes:
415:    Returns the value PETSC_MIN and p = -1 if the vector is of length 0.

417:    Level: intermediate

419:    Concepts: maximum^of vector
420:    Concepts: vector^maximum value

422: .seealso: VecNorm(), VecMin()
423: @*/
424: int VecMax(Vec x,int *p,PetscReal *val)
425: {

432:   PetscLogEventBegin(VEC_Max,x,0,0,0);
433:   (*x->ops->max)(x,p,val);
434:   PetscLogEventEnd(VEC_Max,x,0,0,0);
435:   return(0);
436: }

438: #undef __FUNCT__  
440: /*@
441:    VecMin - Determines the minimum vector component and its location.

443:    Collective on Vec

445:    Input Parameters:
446: .  x - the vector

448:    Output Parameter:
449: +  val - the minimum component
450: -  p - the location of val

452:    Level: intermediate

454:    Notes:
455:    Returns the value PETSC_MAX and p = -1 if the vector is of length 0.

457:    Concepts: minimum^of vector
458:    Concepts: vector^minimum entry

460: .seealso: VecMax()
461: @*/
462: int VecMin(Vec x,int *p,PetscReal *val)
463: {

470:   PetscLogEventBegin(VEC_Min,x,0,0,0);
471:   (*x->ops->min)(x,p,val);
472:   PetscLogEventEnd(VEC_Min,x,0,0,0);
473:   return(0);
474: }

476: #undef __FUNCT__  
478: /*@
479:    VecTDot - Computes an indefinite vector dot product. That is, this
480:    routine does NOT use the complex conjugate.

482:    Collective on Vec

484:    Input Parameters:
485: .  x, y - the vectors

487:    Output Parameter:
488: .  val - the dot product

490:    Notes for Users of Complex Numbers:
491:    For complex vectors, VecTDot() computes the indefinite form
492: $     val = (x,y) = y^T x,
493:    where y^T denotes the transpose of y.

495:    Use VecDot() for the inner product
496: $     val = (x,y) = y^H x,
497:    where y^H denotes the conjugate transpose of y.

499:    Level: intermediate

501:    Concepts: inner product^non-Hermitian
502:    Concepts: vector^inner product
503:    Concepts: non-Hermitian inner product

505: .seealso: VecDot(), VecMTDot()
506: @*/
507: int VecTDot(Vec x,Vec y,PetscScalar *val)
508: {

519:   if (x->N != y->N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
520:   if (x->n != y->n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

522:   PetscLogEventBegin(VEC_TDot,x,y,0,0);
523:   (*x->ops->tdot)(x,y,val);
524:   PetscLogEventEnd(VEC_TDot,x,y,0,0);
525:   return(0);
526: }

528: #undef __FUNCT__  
530: /*@
531:    VecScale - Scales a vector. 

533:    Collective on Vec

535:    Input Parameters:
536: +  x - the vector
537: -  alpha - the scalar

539:    Output Parameter:
540: .  x - the scaled vector

542:    Note:
543:    For a vector with n components, VecScale() computes 
544: $      x[i] = alpha * x[i], for i=1,...,n.

546:    Level: intermediate

548:    Concepts: vector^scaling
549:    Concepts: scaling^vector

551: @*/
552: int VecScale (const PetscScalar *alpha,Vec x)
553: {

560:   PetscLogEventBegin(VEC_Scale,x,0,0,0);
561:   (*x->ops->scale)(alpha,x);
562:   PetscLogEventEnd(VEC_Scale,x,0,0,0);
563:   return(0);
564: }

566: #undef __FUNCT__  
568: /*@
569:    VecCopy - Copies a vector. 

571:    Collective on Vec

573:    Input Parameter:
574: .  x - the vector

576:    Output Parameter:
577: .  y - the copy

579:    Notes:
580:    For default parallel PETSc vectors, both x and y must be distributed in
581:    the same manner; local copies are done.

583:    Level: beginner

585: .seealso: VecDuplicate()
586: @*/
587: int VecCopy(Vec x,Vec y)
588: {

597:   if (x->N != y->N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
598:   if (x->n != y->n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

600:   PetscLogEventBegin(VEC_Copy,x,y,0,0);
601:   (*x->ops->copy)(x,y);
602:   PetscLogEventEnd(VEC_Copy,x,y,0,0);
603:   return(0);
604: }

606: #undef __FUNCT__  
608: /*@
609:    VecSet - Sets all components of a vector to a single scalar value. 

611:    Collective on Vec

613:    Input Parameters:
614: +  alpha - the scalar
615: -  x  - the vector

617:    Output Parameter:
618: .  x  - the vector

620:    Note:
621:    For a vector of dimension n, VecSet() computes
622: $     x[i] = alpha, for i=1,...,n,
623:    so that all vector entries then equal the identical
624:    scalar value, alpha.  Use the more general routine
625:    VecSetValues() to set different vector entries.

627:    Level: beginner

629: .seealso VecSetValues(), VecSetValuesBlocked(), VecSetRandom()

631:    Concepts: vector^setting to constant

633: @*/
634: int VecSet(const PetscScalar *alpha,Vec x)
635: {


643:   PetscLogEventBegin(VEC_Set,x,0,0,0);
644:   (*x->ops->set)(alpha,x);
645:   PetscLogEventEnd(VEC_Set,x,0,0,0);
646:   return(0);
647: }

649: #undef __FUNCT__  
651: /*@C
652:    VecSetRandom - Sets all components of a vector to random numbers.

654:    Collective on Vec

656:    Input Parameters:
657: +  rctx - the random number context, formed by PetscRandomCreate(), or PETSC_NULL and
658:           it will create one internally.
659: -  x  - the vector

661:    Output Parameter:
662: .  x  - the vector

664:    Example of Usage:
665: .vb
666:      PetscRandomCreate(PETSC_COMM_WORLD,RANDOM_DEFAULT,&rctx);
667:      VecSetRandom(rctx,x);
668:      PetscRandomDestroy(rctx);
669: .ve

671:    Level: intermediate

673:    Concepts: vector^setting to random
674:    Concepts: random^vector

676: .seealso: VecSet(), VecSetValues(), PetscRandomCreate(), PetscRandomDestroy()
677: @*/
678: int VecSetRandom(PetscRandom rctx,Vec x)
679: {
680:   int         ierr;
681:   PetscRandom randObj = PETSC_NULL;


688:   if (!rctx) {
689:     MPI_Comm    comm;
690:     PetscObjectGetComm((PetscObject)x,&comm);
691:     PetscRandomCreate(comm,RANDOM_DEFAULT,&randObj);
692:     rctx = randObj;
693:   }

695:   PetscLogEventBegin(VEC_SetRandom,x,rctx,0,0);
696:   (*x->ops->setrandom)(rctx,x);
697:   PetscLogEventEnd(VEC_SetRandom,x,rctx,0,0);
698: 
699:   if (randObj) {
700:     PetscRandomDestroy(randObj);
701:   }
702:   return(0);
703: }

705: #undef __FUNCT__  
707: /*@
708:    VecAXPY - Computes y = alpha x + y. 

710:    Collective on Vec

712:    Input Parameters:
713: +  alpha - the scalar
714: -  x, y  - the vectors

716:    Output Parameter:
717: .  y - output vector

719:    Level: intermediate

721:    Concepts: vector^BLAS
722:    Concepts: BLAS

724: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY()
725: @*/
726: int VecAXPY(const PetscScalar *alpha,Vec x,Vec y)
727: {

738:   if (x->N != y->N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
739:   if (x->n != y->n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

741:   PetscLogEventBegin(VEC_AXPY,x,y,0,0);
742:   (*x->ops->axpy)(alpha,x,y);
743:   PetscLogEventEnd(VEC_AXPY,x,y,0,0);
744:   return(0);
745: }

747: #undef __FUNCT__  
749: /*@
750:    VecAXPBY - Computes y = alpha x + beta y. 

752:    Collective on Vec

754:    Input Parameters:
755: +  alpha,beta - the scalars
756: -  x, y  - the vectors

758:    Output Parameter:
759: .  y - output vector

761:    Level: intermediate

763:    Concepts: BLAS
764:    Concepts: vector^BLAS

766: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPY()
767: @*/
768: int VecAXPBY(const PetscScalar *alpha,const PetscScalar *beta,Vec x,Vec y)
769: {

781:   if (x->N != y->N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
782:   if (x->n != y->n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

784:   PetscLogEventBegin(VEC_AXPY,x,y,0,0);
785:   (*x->ops->axpby)(alpha,beta,x,y);
786:   PetscLogEventEnd(VEC_AXPY,x,y,0,0);
787:   return(0);
788: }

790: #undef __FUNCT__  
792: /*@
793:    VecAYPX - Computes y = x + alpha y.

795:    Collective on Vec

797:    Input Parameters:
798: +  alpha - the scalar
799: -  x, y  - the vectors

801:    Output Parameter:
802: .  y - output vector

804:    Level: intermediate

806:    Concepts: vector^BLAS
807:    Concepts: BLAS

809: .seealso: VecAXPY(), VecWAXPY()
810: @*/
811: int VecAYPX(const PetscScalar *alpha,Vec x,Vec y)
812: {

823:   if (x->N != y->N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
824:   if (x->n != y->n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

826:   PetscLogEventBegin(VEC_AYPX,x,y,0,0);
827:    (*x->ops->aypx)(alpha,x,y);
828:   PetscLogEventEnd(VEC_AYPX,x,y,0,0);
829:   return(0);
830: }

832: #undef __FUNCT__  
834: /*@
835:    VecSwap - Swaps the vectors x and y.

837:    Collective on Vec

839:    Input Parameters:
840: .  x, y  - the vectors

842:    Level: advanced

844:    Concepts: vector^swapping values

846: @*/
847: int VecSwap(Vec x,Vec y)
848: {

858:   if (x->N != y->N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
859:   if (x->n != y->n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

861:   PetscLogEventBegin(VEC_Swap,x,y,0,0);
862:   (*x->ops->swap)(x,y);
863:   PetscLogEventEnd(VEC_Swap,x,y,0,0);
864:   return(0);
865: }

867: #undef __FUNCT__  
869: /*@
870:    VecWAXPY - Computes w = alpha x + y.

872:    Collective on Vec

874:    Input Parameters:
875: +  alpha - the scalar
876: -  x, y  - the vectors

878:    Output Parameter:
879: .  w - the result

881:    Level: intermediate

883:    Concepts: vector^BLAS
884:    Concepts: BLAS

886: .seealso: VecAXPY(), VecAYPX()
887: @*/
888: int VecWAXPY(const PetscScalar *alpha,Vec x,Vec y,Vec w)
889: {

902:   if (x->N != y->N || x->N != w->N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
903:   if (x->n != y->n || x->n != w->n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

905:   PetscLogEventBegin(VEC_WAXPY,x,y,w,0);
906:    (*x->ops->waxpy)(alpha,x,y,w);
907:   PetscLogEventEnd(VEC_WAXPY,x,y,w,0);
908:   return(0);
909: }

911: #undef __FUNCT__  
913: /*@
914:    VecPointwiseMult - Computes the componentwise multiplication w = x*y.

916:    Collective on Vec

918:    Input Parameters:
919: .  x, y  - the vectors

921:    Output Parameter:
922: .  w - the result

924:    Level: advanced

926:    Notes: any subset of the x, y, and w may be the same vector.

928:    Concepts: vector^pointwise multiply

930: .seealso: VecPointwiseDivide()
931: @*/
932: int VecPointwiseMult(Vec x,Vec y,Vec w)
933: {

945:   if (x->N != y->N || x->N != w->N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
946:   if (x->n != y->n || x->n != w->n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

948:   PetscLogEventBegin(VEC_PointwiseMult,x,y,w,0);
949:   (*x->ops->pointwisemult)(x,y,w);
950:   PetscLogEventEnd(VEC_PointwiseMult,x,y,w,0);
951:   return(0);
952: }

954: #undef __FUNCT__  
956: /*@
957:    VecPointwiseDivide - Computes the componentwise division w = x/y.

959:    Collective on Vec

961:    Input Parameters:
962: .  x, y  - the vectors

964:    Output Parameter:
965: .  w - the result

967:    Level: advanced

969:    Notes: any subset of the x, y, and w may be the same vector.

971:    Concepts: vector^pointwise divide

973: .seealso: VecPointwiseMult()
974: @*/
975: int VecPointwiseDivide(Vec x,Vec y,Vec w)
976: {

988:   if (x->N != y->N || x->N != w->N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
989:   if (x->n != y->n || x->n != w->n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

991:   (*x->ops->pointwisedivide)(x,y,w);
992:   return(0);
993: }

995: #undef __FUNCT__  
997: /*@
998:    VecMaxPointwiseDivide - Computes the maximum of the componentwise division w = abs(x/y).

1000:    Collective on Vec

1002:    Input Parameters:
1003: .  x, y  - the vectors

1005:    Output Parameter:
1006: .  w - the result

1008:    Level: advanced

1010:    Notes: any subset of the x, y, and w may be the same vector.

1012: .seealso: VecPointwiseDivide(), VecPointwiseMult()
1013: @*/
1014: int VecMaxPointwiseDivide(Vec x,Vec y,PetscReal *max)
1015: {

1025:   if (x->N != y->N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
1026:   if (x->n != y->n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

1028:   (*x->ops->maxpointwisedivide)(x,y,max);
1029:   return(0);
1030: }

1032: #undef __FUNCT__  
1034: /*@C
1035:    VecDuplicate - Creates a new vector of the same type as an existing vector.

1037:    Collective on Vec

1039:    Input Parameters:
1040: .  v - a vector to mimic

1042:    Output Parameter:
1043: .  newv - location to put new vector

1045:    Notes:
1046:    VecDuplicate() does not copy the vector, but rather allocates storage
1047:    for the new vector.  Use VecCopy() to copy a vector.

1049:    Use VecDestroy() to free the space. Use VecDuplicateVecs() to get several
1050:    vectors. 

1052:    Level: beginner

1054: .seealso: VecDestroy(), VecDuplicateVecs(), VecCreate(), VecCopy()
1055: @*/
1056: int VecDuplicate(Vec x,Vec *newv)
1057: {

1064:   (*x->ops->duplicate)(x,newv);
1065:   return(0);
1066: }

1068: #undef __FUNCT__  
1070: /*@C
1071:    VecDestroy - Destroys a vector.

1073:    Collective on Vec

1075:    Input Parameters:
1076: .  v  - the vector

1078:    Level: beginner

1080: .seealso: VecDuplicate(), VecDestroyVecs()
1081: @*/
1082: int VecDestroy(Vec v)
1083: {

1088:   if (--v->refct > 0) return(0);
1089:   /* destroy the internal part */
1090:   if (v->ops->destroy) {
1091:     (*v->ops->destroy)(v);
1092:   }
1093:   /* destroy the external/common part */
1094:   if (v->mapping) {
1095:     ISLocalToGlobalMappingDestroy(v->mapping);
1096:   }
1097:   if (v->bmapping) {
1098:     ISLocalToGlobalMappingDestroy(v->bmapping);
1099:   }
1100:   if (v->map) {
1101:     PetscMapDestroy(v->map);
1102:   }
1103:   PetscLogObjectDestroy(v);
1104:   PetscHeaderDestroy(v);
1105:   return(0);
1106: }

1108: #undef __FUNCT__  
1110: /*@C
1111:    VecDuplicateVecs - Creates several vectors of the same type as an existing vector.

1113:    Collective on Vec

1115:    Input Parameters:
1116: +  m - the number of vectors to obtain
1117: -  v - a vector to mimic

1119:    Output Parameter:
1120: .  V - location to put pointer to array of vectors

1122:    Notes:
1123:    Use VecDestroyVecs() to free the space. Use VecDuplicate() to form a single
1124:    vector.

1126:    Fortran Note:
1127:    The Fortran interface is slightly different from that given below, it 
1128:    requires one to pass in V a Vec (integer) array of size at least m.
1129:    See the Fortran chapter of the users manual and petsc/src/vec/examples for details.

1131:    Level: intermediate

1133: .seealso:  VecDestroyVecs(), VecDuplicate(), VecCreate(), VecDuplicateVecsF90()
1134: @*/
1135: int VecDuplicateVecs(Vec v,int m,Vec *V[])
1136: {

1143:   (*v->ops->duplicatevecs)(v, m,V);
1144:   return(0);
1145: }

1147: #undef __FUNCT__  
1149: /*@C
1150:    VecDestroyVecs - Frees a block of vectors obtained with VecDuplicateVecs().

1152:    Collective on Vec

1154:    Input Parameters:
1155: +  vv - pointer to array of vector pointers
1156: -  m - the number of vectors previously obtained

1158:    Fortran Note:
1159:    The Fortran interface is slightly different from that given below.
1160:    See the Fortran chapter of the users manual and 
1161:    petsc/src/vec/examples for details.

1163:    Level: intermediate

1165: .seealso: VecDuplicateVecs(), VecDestroyVecsF90()
1166: @*/
1167: int VecDestroyVecs(const Vec vv[],int m)
1168: {

1172:   if (!vv) SETERRQ(PETSC_ERR_ARG_BADPTR,"Null vectors");
1175:   (*(*vv)->ops->destroyvecs)(vv,m);
1176:   return(0);
1177: }

1179: #undef __FUNCT__  
1181: /*@
1182:    VecSetValues - Inserts or adds values into certain locations of a vector. 

1184:    Input Parameters:
1185:    Not Collective

1187: +  x - vector to insert in
1188: .  ni - number of elements to add
1189: .  ix - indices where to add
1190: .  y - array of values
1191: -  iora - either INSERT_VALUES or ADD_VALUES, where
1192:    ADD_VALUES adds values to any existing entries, and
1193:    INSERT_VALUES replaces existing entries with new values

1195:    Notes: 
1196:    VecSetValues() sets x[ix[i]] = y[i], for i=0,...,ni-1.

1198:    Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES 
1199:    options cannot be mixed without intervening calls to the assembly
1200:    routines.

1202:    These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd() 
1203:    MUST be called after all calls to VecSetValues() have been completed.

1205:    VecSetValues() uses 0-based indices in Fortran as well as in C.

1207:    Negative indices may be passed in ix, these rows are 
1208:    simply ignored. This allows easily inserting element load matrices
1209:    with homogeneous Dirchlet boundary conditions that you don't want represented
1210:    in the vector.

1212:    Level: beginner

1214:    Concepts: vector^setting values

1216: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesLocal(),
1217:            VecSetValue(), VecSetValuesBlocked()
1218: @*/
1219: int VecSetValues(Vec x,int ni,const int ix[],const PetscScalar y[],InsertMode iora)
1220: {

1228:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
1229:   (*x->ops->setvalues)(x,ni,ix,y,iora);
1230:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
1231:   return(0);
1232: }

1234: #undef __FUNCT__  
1236: /*@
1237:    VecSetValuesBlocked - Inserts or adds blocks of values into certain locations of a vector. 

1239:    Not Collective

1241:    Input Parameters:
1242: +  x - vector to insert in
1243: .  ni - number of blocks to add
1244: .  ix - indices where to add in block count, rather than element count
1245: .  y - array of values
1246: -  iora - either INSERT_VALUES or ADD_VALUES, where
1247:    ADD_VALUES adds values to any existing entries, and
1248:    INSERT_VALUES replaces existing entries with new values

1250:    Notes: 
1251:    VecSetValuesBlocked() sets x[ix[bs*i]+j] = y[bs*i+j], 
1252:    for j=0,...,bs, for i=0,...,ni-1. where bs was set with VecSetBlockSize().

1254:    Calls to VecSetValuesBlocked() with the INSERT_VALUES and ADD_VALUES 
1255:    options cannot be mixed without intervening calls to the assembly
1256:    routines.

1258:    These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd() 
1259:    MUST be called after all calls to VecSetValuesBlocked() have been completed.

1261:    VecSetValuesBlocked() uses 0-based indices in Fortran as well as in C.

1263:    Negative indices may be passed in ix, these rows are 
1264:    simply ignored. This allows easily inserting element load matrices
1265:    with homogeneous Dirchlet boundary conditions that you don't want represented
1266:    in the vector.

1268:    Level: intermediate

1270:    Concepts: vector^setting values blocked

1272: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesBlockedLocal(),
1273:            VecSetValues()
1274: @*/
1275: int VecSetValuesBlocked(Vec x,int ni,const int ix[],const PetscScalar y[],InsertMode iora)
1276: {

1284:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
1285:   (*x->ops->setvaluesblocked)(x,ni,ix,y,iora);
1286:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
1287:   return(0);
1288: }

1290: /*MC
1291:    VecSetValue - Set a single entry into a vector.

1293:    Synopsis:
1294:    int VecSetValue(Vec v,int row,PetscScalar value, InsertMode mode);

1296:    Not Collective

1298:    Input Parameters:
1299: +  v - the vector
1300: .  row - the row location of the entry
1301: .  value - the value to insert
1302: -  mode - either INSERT_VALUES or ADD_VALUES

1304:    Notes:
1305:    For efficiency one should use VecSetValues() and set several or 
1306:    many values simultaneously if possible.

1308:    Note that VecSetValue() does NOT return an error code (since this
1309:    is checked internally).

1311:    These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd() 
1312:    MUST be called after all calls to VecSetValues() have been completed.

1314:    VecSetValues() uses 0-based indices in Fortran as well as in C.

1316:    Level: beginner

1318: .seealso: VecSetValues(), VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesBlockedLocal()
1319: M*/

1321: #undef __FUNCT__  
1323: /*@
1324:    VecSetLocalToGlobalMapping - Sets a local numbering to global numbering used
1325:    by the routine VecSetValuesLocal() to allow users to insert vector entries
1326:    using a local (per-processor) numbering.

1328:    Collective on Vec

1330:    Input Parameters:
1331: +  x - vector
1332: -  mapping - mapping created with ISLocalToGlobalMappingCreate() or ISLocalToGlobalMappingCreateIS()

1334:    Notes: 
1335:    All vectors obtained with VecDuplicate() from this vector inherit the same mapping.

1337:    Level: intermediate

1339:    Concepts: vector^setting values with local numbering

1341: seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetValuesLocal(),
1342:            VecSetLocalToGlobalMappingBlocked(), VecSetValuesBlockedLocal()
1343: @*/
1344: int VecSetLocalToGlobalMapping(Vec x,ISLocalToGlobalMapping mapping)
1345: {

1351:   if (x->mapping) {
1352:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Mapping already set for vector");
1353:   }

1355:   if (x->ops->setlocaltoglobalmapping) {
1356:     (*x->ops->setlocaltoglobalmapping)(x,mapping);
1357:   } else {
1358:     x->mapping = mapping;
1359:     PetscObjectReference((PetscObject)mapping);
1360:   }
1361:   return(0);
1362: }

1364: #undef __FUNCT__  
1366: /*@
1367:    VecSetLocalToGlobalMappingBlock - Sets a local numbering to global numbering used
1368:    by the routine VecSetValuesBlockedLocal() to allow users to insert vector entries
1369:    using a local (per-processor) numbering.

1371:    Collective on Vec

1373:    Input Parameters:
1374: +  x - vector
1375: -  mapping - mapping created with ISLocalToGlobalMappingCreate() or ISLocalToGlobalMappingCreateIS()

1377:    Notes: 
1378:    All vectors obtained with VecDuplicate() from this vector inherit the same mapping.

1380:    Level: intermediate

1382:    Concepts: vector^setting values blocked with local numbering

1384: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetValuesLocal(),
1385:            VecSetLocalToGlobalMapping(), VecSetValuesBlockedLocal()
1386: @*/
1387: int VecSetLocalToGlobalMappingBlock(Vec x,ISLocalToGlobalMapping mapping)
1388: {

1394:   if (x->bmapping) {
1395:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Mapping already set for vector");
1396:   }
1397:   x->bmapping = mapping;
1398:   PetscObjectReference((PetscObject)mapping);
1399:   return(0);
1400: }

1402: #undef __FUNCT__  
1404: /*@
1405:    VecSetValuesLocal - Inserts or adds values into certain locations of a vector,
1406:    using a local ordering of the nodes. 

1408:    Not Collective

1410:    Input Parameters:
1411: +  x - vector to insert in
1412: .  ni - number of elements to add
1413: .  ix - indices where to add
1414: .  y - array of values
1415: -  iora - either INSERT_VALUES or ADD_VALUES, where
1416:    ADD_VALUES adds values to any existing entries, and
1417:    INSERT_VALUES replaces existing entries with new values

1419:    Level: intermediate

1421:    Notes: 
1422:    VecSetValuesLocal() sets x[ix[i]] = y[i], for i=0,...,ni-1.

1424:    Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES 
1425:    options cannot be mixed without intervening calls to the assembly
1426:    routines.

1428:    These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd() 
1429:    MUST be called after all calls to VecSetValuesLocal() have been completed.

1431:    VecSetValuesLocal() uses 0-based indices in Fortran as well as in C.

1433:    Concepts: vector^setting values with local numbering

1435: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetLocalToGlobalMapping(),
1436:            VecSetValuesBlockedLocal()
1437: @*/
1438: int VecSetValuesLocal(Vec x,int ni,const int ix[],const PetscScalar y[],InsertMode iora)
1439: {
1440:   int ierr,lixp[128],*lix = lixp;


1448:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
1449:   if (!x->ops->setvalueslocal) {
1450:     if (!x->mapping) {
1451:       SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Local to global never set with VecSetLocalToGlobalMapping()");
1452:     }
1453:     if (ni > 128) {
1454:       PetscMalloc(ni*sizeof(int),&lix);
1455:     }
1456:     ISLocalToGlobalMappingApply(x->mapping,ni,(int*)ix,lix);
1457:     (*x->ops->setvalues)(x,ni,lix,y,iora);
1458:     if (ni > 128) {
1459:       PetscFree(lix);
1460:     }
1461:   } else {
1462:     (*x->ops->setvalueslocal)(x,ni,ix,y,iora);
1463:   }
1464:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
1465:   return(0);
1466: }

1468: #undef __FUNCT__  
1470: /*@
1471:    VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
1472:    using a local ordering of the nodes. 

1474:    Not Collective

1476:    Input Parameters:
1477: +  x - vector to insert in
1478: .  ni - number of blocks to add
1479: .  ix - indices where to add in block count, not element count
1480: .  y - array of values
1481: -  iora - either INSERT_VALUES or ADD_VALUES, where
1482:    ADD_VALUES adds values to any existing entries, and
1483:    INSERT_VALUES replaces existing entries with new values

1485:    Level: intermediate

1487:    Notes: 
1488:    VecSetValuesBlockedLocal() sets x[bs*ix[i]+j] = y[bs*i+j], 
1489:    for j=0,..bs-1, for i=0,...,ni-1, where bs has been set with VecSetBlockSize().

1491:    Calls to VecSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES 
1492:    options cannot be mixed without intervening calls to the assembly
1493:    routines.

1495:    These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd() 
1496:    MUST be called after all calls to VecSetValuesBlockedLocal() have been completed.

1498:    VecSetValuesBlockedLocal() uses 0-based indices in Fortran as well as in C.


1501:    Concepts: vector^setting values blocked with local numbering

1503: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetValuesBlocked(), 
1504:            VecSetLocalToGlobalMappingBlocked()
1505: @*/
1506: int VecSetValuesBlockedLocal(Vec x,int ni,const int ix[],const PetscScalar y[],InsertMode iora)
1507: {
1508:   int ierr,lixp[128],*lix = lixp;

1515:   if (!x->bmapping) {
1516:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Local to global never set with VecSetLocalToGlobalMappingBlocked()");
1517:   }
1518:   if (ni > 128) {
1519:     PetscMalloc(ni*sizeof(int),&lix);
1520:   }

1522:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
1523:   ISLocalToGlobalMappingApply(x->bmapping,ni,(int*)ix,lix);
1524:   (*x->ops->setvaluesblocked)(x,ni,lix,y,iora);
1525:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
1526:   if (ni > 128) {
1527:     PetscFree(lix);
1528:   }
1529:   return(0);
1530: }

1532: #undef __FUNCT__  
1534: /*@
1535:    VecAssemblyBegin - Begins assembling the vector.  This routine should
1536:    be called after completing all calls to VecSetValues().

1538:    Collective on Vec

1540:    Input Parameter:
1541: .  vec - the vector

1543:    Level: beginner

1545:    Concepts: assembly^vectors

1547: .seealso: VecAssemblyEnd(), VecSetValues()
1548: @*/
1549: int VecAssemblyBegin(Vec vec)
1550: {
1551:   int        ierr;
1552:   PetscTruth flg;


1558:   PetscOptionsHasName(vec->prefix,"-vec_view_stash",&flg);
1559:   if (flg) {
1560:     VecStashView(vec,PETSC_VIEWER_STDOUT_(vec->comm));
1561:   }

1563:   PetscLogEventBegin(VEC_AssemblyBegin,vec,0,0,0);
1564:   if (vec->ops->assemblybegin) {
1565:     (*vec->ops->assemblybegin)(vec);
1566:   }
1567:   PetscLogEventEnd(VEC_AssemblyBegin,vec,0,0,0);
1568:   return(0);
1569: }

1571: #undef __FUNCT__  
1573: /*@
1574:    VecAssemblyEnd - Completes assembling the vector.  This routine should
1575:    be called after VecAssemblyBegin().

1577:    Collective on Vec

1579:    Input Parameter:
1580: .  vec - the vector

1582:    Options Database Keys:
1583: .  -vec_view - Prints vector in ASCII format
1584: .  -vec_view_matlab - Prints vector in Matlab format
1585: .  -vec_view_draw - Activates vector viewing using drawing tools
1586: .  -display <name> - Sets display name (default is host)
1587: .  -draw_pause <sec> - Sets number of seconds to pause after display
1588: .  -vec_view_socket - Activates vector viewing using a socket
1589: -  -vec_view_ams - Activates vector viewing using the ALICE Memory Snooper (AMS)
1590:  
1591:    Level: beginner

1593: .seealso: VecAssemblyBegin(), VecSetValues()
1594: @*/
1595: int VecAssemblyEnd(Vec vec)
1596: {
1597:   int        ierr;
1598:   PetscTruth flg;

1602:   PetscLogEventBegin(VEC_AssemblyEnd,vec,0,0,0);
1604:   if (vec->ops->assemblyend) {
1605:     (*vec->ops->assemblyend)(vec);
1606:   }
1607:   PetscLogEventEnd(VEC_AssemblyEnd,vec,0,0,0);
1608:   PetscOptionsHasName(vec->prefix,"-vec_view",&flg);
1609:   if (flg) {
1610:     VecView(vec,PETSC_VIEWER_STDOUT_(vec->comm));
1611:   }
1612:   PetscOptionsHasName(PETSC_NULL,"-vec_view_matlab",&flg);
1613:   if (flg) {
1614:     PetscViewerPushFormat(PETSC_VIEWER_STDOUT_(vec->comm),PETSC_VIEWER_ASCII_MATLAB);
1615:     VecView(vec,PETSC_VIEWER_STDOUT_(vec->comm));
1616:     PetscViewerPopFormat(PETSC_VIEWER_STDOUT_(vec->comm));
1617:   }
1618:   PetscOptionsHasName(PETSC_NULL,"-vec_view_draw",&flg);
1619:   if (flg) {
1620:     VecView(vec,PETSC_VIEWER_DRAW_(vec->comm));
1621:     PetscViewerFlush(PETSC_VIEWER_DRAW_(vec->comm));
1622:   }
1623:   PetscOptionsHasName(PETSC_NULL,"-vec_view_draw_lg",&flg);
1624:   if (flg) {
1625:     PetscViewerSetFormat(PETSC_VIEWER_DRAW_(vec->comm),PETSC_VIEWER_DRAW_LG);
1626:     VecView(vec,PETSC_VIEWER_DRAW_(vec->comm));
1627:     PetscViewerFlush(PETSC_VIEWER_DRAW_(vec->comm));
1628:   }
1629:   PetscOptionsHasName(PETSC_NULL,"-vec_view_socket",&flg);
1630:   if (flg) {
1631:     VecView(vec,PETSC_VIEWER_SOCKET_(vec->comm));
1632:     PetscViewerFlush(PETSC_VIEWER_SOCKET_(vec->comm));
1633:   }
1634:   PetscOptionsHasName(PETSC_NULL,"-vec_view_binary",&flg);
1635:   if (flg) {
1636:     VecView(vec,PETSC_VIEWER_BINARY_(vec->comm));
1637:     PetscViewerFlush(PETSC_VIEWER_BINARY_(vec->comm));
1638:   }
1639: #if defined(PETSC_HAVE_AMS)
1640:   PetscOptionsHasName(PETSC_NULL,"-vec_view_ams",&flg);
1641:   if (flg) {
1642:     VecView(vec,PETSC_VIEWER_AMS_(vec->comm));
1643:   }
1644: #endif
1645:   return(0);
1646: }


1649: #undef __FUNCT__  
1651: /*@C
1652:    VecMTDot - Computes indefinite vector multiple dot products. 
1653:    That is, it does NOT use the complex conjugate.

1655:    Collective on Vec

1657:    Input Parameters:
1658: +  nv - number of vectors
1659: .  x - one vector
1660: -  y - array of vectors.  Note that vectors are pointers

1662:    Output Parameter:
1663: .  val - array of the dot products

1665:    Notes for Users of Complex Numbers:
1666:    For complex vectors, VecMTDot() computes the indefinite form
1667: $      val = (x,y) = y^T x,
1668:    where y^T denotes the transpose of y.

1670:    Use VecMDot() for the inner product
1671: $      val = (x,y) = y^H x,
1672:    where y^H denotes the conjugate transpose of y.

1674:    Level: intermediate

1676:    Concepts: inner product^multiple
1677:    Concepts: vector^multiple inner products

1679: .seealso: VecMDot(), VecTDot()
1680: @*/
1681: int VecMTDot(int nv,Vec x,const Vec y[],PetscScalar *val)
1682: {

1693:   if (x->N != (*y)->N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
1694:   if (x->n != (*y)->n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

1696:   PetscLogEventBegin(VEC_MTDot,x,*y,0,0);
1697:   (*x->ops->mtdot)(nv,x,y,val);
1698:   PetscLogEventEnd(VEC_MTDot,x,*y,0,0);
1699:   return(0);
1700: }

1702: #undef __FUNCT__  
1704: /*@C
1705:    VecMDot - Computes vector multiple dot products. 

1707:    Collective on Vec

1709:    Input Parameters:
1710: +  nv - number of vectors
1711: .  x - one vector
1712: -  y - array of vectors. 

1714:    Output Parameter:
1715: .  val - array of the dot products

1717:    Notes for Users of Complex Numbers:
1718:    For complex vectors, VecMDot() computes 
1719: $     val = (x,y) = y^H x,
1720:    where y^H denotes the conjugate transpose of y.

1722:    Use VecMTDot() for the indefinite form
1723: $     val = (x,y) = y^T x,
1724:    where y^T denotes the transpose of y.

1726:    Level: intermediate

1728:    Concepts: inner product^multiple
1729:    Concepts: vector^multiple inner products

1731: .seealso: VecMTDot(), VecDot()
1732: @*/
1733: int VecMDot(int nv,Vec x,const Vec y[],PetscScalar *val)
1734: {

1745:   if (x->N != (*y)->N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
1746:   if (x->n != (*y)->n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

1748:   PetscLogEventBarrierBegin(VEC_MDotBarrier,x,*y,0,0,x->comm);
1749:   (*x->ops->mdot)(nv,x,y,val);
1750:   PetscLogEventBarrierEnd(VEC_MDotBarrier,x,*y,0,0,x->comm);
1751:   return(0);
1752: }

1754: #undef __FUNCT__  
1756: /*@C
1757:    VecMAXPY - Computes y = y + sum alpha[j] x[j]

1759:    Collective on Vec

1761:    Input Parameters:
1762: +  nv - number of scalars and x-vectors
1763: .  alpha - array of scalars
1764: .  y - one vector
1765: -  x - array of vectors

1767:    Level: intermediate

1769:    Concepts: BLAS

1771: .seealso: VecAXPY(), VecWAXPY(), VecAYPX()
1772: @*/
1773: int  VecMAXPY(int nv,const PetscScalar *alpha,Vec y,Vec *x)
1774: {

1785:   if (y->N != (*x)->N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
1786:   if (y->n != (*x)->n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

1788:   PetscLogEventBegin(VEC_MAXPY,*x,y,0,0);
1789:   (*y->ops->maxpy)(nv,alpha,y,x);
1790:   PetscLogEventEnd(VEC_MAXPY,*x,y,0,0);
1791:   return(0);
1792: }

1794: #undef __FUNCT__  
1796: /*@C
1797:    VecGetArray - Returns a pointer to a contiguous array that contains this 
1798:    processor's portion of the vector data. For the standard PETSc
1799:    vectors, VecGetArray() returns a pointer to the local data array and
1800:    does not use any copies. If the underlying vector data is not stored
1801:    in a contiquous array this routine will copy the data to a contiquous
1802:    array and return a pointer to that. You MUST call VecRestoreArray() 
1803:    when you no longer need access to the array.

1805:    Not Collective

1807:    Input Parameter:
1808: .  x - the vector

1810:    Output Parameter:
1811: .  a - location to put pointer to the array

1813:    Fortran Note:
1814:    This routine is used differently from Fortran 77
1815: $    Vec         x
1816: $    PetscScalar x_array(1)
1817: $    PetscOffset i_x
1818: $    int         ierr
1819: $       call VecGetArray(x,x_array,i_x,ierr)
1820: $
1821: $   Access first local entry in vector with
1822: $      value = x_array(i_x + 1)
1823: $
1824: $      ...... other code
1825: $       call VecRestoreArray(x,x_array,i_x,ierr)
1826:    For Fortran 90 see VecGetArrayF90()

1828:    See the Fortran chapter of the users manual and 
1829:    petsc/src/snes/examples/tutorials/ex5f.F for details.

1831:    Level: beginner

1833:    Concepts: vector^accessing local values

1835: .seealso: VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(), VecGetArray2d()
1836: @*/
1837: int VecGetArray(Vec x,PetscScalar *a[])
1838: {

1845:   (*x->ops->getarray)(x,a);
1846:   return(0);
1847: }


1850: #undef __FUNCT__  
1852: /*@C
1853:    VecGetArrays - Returns a pointer to the arrays in a set of vectors
1854:    that were created by a call to VecDuplicateVecs().  You MUST call
1855:    VecRestoreArrays() when you no longer need access to the array.

1857:    Not Collective

1859:    Input Parameter:
1860: +  x - the vectors
1861: -  n - the number of vectors

1863:    Output Parameter:
1864: .  a - location to put pointer to the array

1866:    Fortran Note:
1867:    This routine is not supported in Fortran.

1869:    Level: intermediate

1871: .seealso: VecGetArray(), VecRestoreArrays()
1872: @*/
1873: int VecGetArrays(const Vec x[],int n,PetscScalar **a[])
1874: {
1875:   int         i,ierr;
1876:   PetscScalar **q;

1881:   if (n <= 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Must get at least one array n = %d",n);
1882:   PetscMalloc(n*sizeof(PetscScalar*),&q);
1883:   for (i=0; i<n; ++i) {
1884:     VecGetArray(x[i],&q[i]);
1885:   }
1886:   *a = q;
1887:   return(0);
1888: }

1890: #undef __FUNCT__  
1892: /*@C
1893:    VecRestoreArrays - Restores a group of vectors after VecGetArrays()
1894:    has been called.

1896:    Not Collective

1898:    Input Parameters:
1899: +  x - the vector
1900: .  n - the number of vectors
1901: -  a - location of pointer to arrays obtained from VecGetArrays()

1903:    Notes:
1904:    For regular PETSc vectors this routine does not involve any copies. For
1905:    any special vectors that do not store local vector data in a contiguous
1906:    array, this routine will copy the data back into the underlying 
1907:    vector data structure from the arrays obtained with VecGetArrays().

1909:    Fortran Note:
1910:    This routine is not supported in Fortran.

1912:    Level: intermediate

1914: .seealso: VecGetArrays(), VecRestoreArray()
1915: @*/
1916: int VecRestoreArrays(const Vec x[],int n,PetscScalar **a[])
1917: {
1918:   int         i,ierr;
1919:   PetscScalar **q = *a;


1925:   for(i=0;i<n;++i) {
1926:     VecRestoreArray(x[i],&q[i]);
1927:   }
1928:   PetscFree(q);
1929:   return(0);
1930: }

1932: #undef __FUNCT__  
1934: /*@C
1935:    VecRestoreArray - Restores a vector after VecGetArray() has been called.

1937:    Not Collective

1939:    Input Parameters:
1940: +  x - the vector
1941: -  a - location of pointer to array obtained from VecGetArray()

1943:    Level: beginner

1945:    Notes:
1946:    For regular PETSc vectors this routine does not involve any copies. For
1947:    any special vectors that do not store local vector data in a contiguous
1948:    array, this routine will copy the data back into the underlying 
1949:    vector data structure from the array obtained with VecGetArray().

1951:    This routine actually zeros out the a pointer. This is to prevent accidental
1952:    us of the array after it has been restored. If you pass null for a it will 
1953:    not zero the array pointer a.

1955:    Fortran Note:
1956:    This routine is used differently from Fortran 77
1957: $    Vec         x
1958: $    PetscScalar x_array(1)
1959: $    PetscOffset i_x
1960: $    int         ierr
1961: $       call VecGetArray(x,x_array,i_x,ierr)
1962: $
1963: $   Access first local entry in vector with
1964: $      value = x_array(i_x + 1)
1965: $
1966: $      ...... other code
1967: $       call VecRestoreArray(x,x_array,i_x,ierr)

1969:    See the Fortran chapter of the users manual and 
1970:    petsc/src/snes/examples/tutorials/ex5f.F for details.
1971:    For Fortran 90 see VecRestoreArrayF90()

1973: .seealso: VecGetArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(), VecRestoreArray2d()
1974: @*/
1975: int VecRestoreArray(Vec x,PetscScalar *a[])
1976: {

1983:   if (x->ops->restorearray) {
1984:     (*x->ops->restorearray)(x,a);
1985:   }
1986:   return(0);
1987: }

1989: #undef  __FUNCT__
1991: /*@
1992:   VecViewFromOptions - This function visualizes the vector based upon user options.

1994:   Collective on Vec

1996:   Input Parameters:
1997: . vec   - The vector
1998: . title - The title

2000:   Level: intermediate

2002: .keywords: Vec, view, options, database
2003: .seealso: VecSetFromOptions(), VecView()
2004: @*/
2005: int VecViewFromOptions(Vec vec, char *title)
2006: {
2007:   PetscViewer viewer;
2008:   PetscDraw   draw;
2009:   PetscTruth  opt;
2010:   char       *titleStr;
2011:   char        typeName[1024];
2012:   char        fileName[PETSC_MAX_PATH_LEN];
2013:   int         len;
2014:   int         ierr;

2017:   PetscOptionsHasName(vec->prefix, "-vec_view", &opt);
2018:   if (opt == PETSC_TRUE) {
2019:     PetscOptionsGetString(vec->prefix, "-vec_view", typeName, 1024, &opt);
2020:     PetscStrlen(typeName, &len);
2021:     if (len > 0) {
2022:       PetscViewerCreate(vec->comm, &viewer);
2023:       PetscViewerSetType(viewer, typeName);
2024:       PetscOptionsGetString(vec->prefix, "-vec_view_file", fileName, 1024, &opt);
2025:       if (opt == PETSC_TRUE) {
2026:         PetscViewerSetFilename(viewer, fileName);
2027:       } else {
2028:         PetscViewerSetFilename(viewer, vec->name);
2029:       }
2030:       VecView(vec, viewer);
2031:       PetscViewerFlush(viewer);
2032:       PetscViewerDestroy(viewer);
2033:     } else {
2034:       VecView(vec, PETSC_VIEWER_STDOUT_(vec->comm));
2035:     }
2036:   }
2037:   PetscOptionsHasName(vec->prefix, "-vec_view_draw", &opt);
2038:   if (opt == PETSC_TRUE) {
2039:     PetscViewerDrawOpen(vec->comm, 0, 0, 0, 0, 300, 300, &viewer);
2040:     PetscViewerDrawGetDraw(viewer, 0, &draw);
2041:     if (title != PETSC_NULL) {
2042:       titleStr = title;
2043:     } else {
2044:       PetscObjectName((PetscObject) vec);                                                          CHKERRQ(ierr) ;
2045:       titleStr = vec->name;
2046:     }
2047:     PetscDrawSetTitle(draw, titleStr);
2048:     VecView(vec, viewer);
2049:     PetscViewerFlush(viewer);
2050:     PetscDrawPause(draw);
2051:     PetscViewerDestroy(viewer);
2052:   }
2053:   return(0);
2054: }

2056: #undef __FUNCT__  
2058: /*@C
2059:    VecView - Views a vector object. 

2061:    Collective on Vec

2063:    Input Parameters:
2064: +  v - the vector
2065: -  viewer - an optional visualization context

2067:    Notes:
2068:    The available visualization contexts include
2069: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
2070: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
2071:          output where only the first processor opens
2072:          the file.  All other processors send their 
2073:          data to the first processor to print. 

2075:    You can change the format the vector is printed using the 
2076:    option PetscViewerSetFormat().

2078:    The user can open alternative visualization contexts with
2079: +    PetscViewerASCIIOpen() - Outputs vector to a specified file
2080: .    PetscViewerBinaryOpen() - Outputs vector in binary to a
2081:          specified file; corresponding input uses VecLoad()
2082: .    PetscViewerDrawOpen() - Outputs vector to an X window display
2083: -    PetscViewerSocketOpen() - Outputs vector to Socket viewer

2085:    The user can call PetscViewerSetFormat() to specify the output
2086:    format of ASCII printed objects (when using PETSC_VIEWER_STDOUT_SELF,
2087:    PETSC_VIEWER_STDOUT_WORLD and PetscViewerASCIIOpen).  Available formats include
2088: +    PETSC_VIEWER_ASCII_DEFAULT - default, prints vector contents
2089: .    PETSC_VIEWER_ASCII_MATLAB - prints vector contents in Matlab format
2090: .    PETSC_VIEWER_ASCII_INDEX - prints vector contents, including indices of vector elements
2091: -    PETSC_VIEWER_ASCII_COMMON - prints vector contents, using a 
2092:          format common among all vector types

2094:    Level: beginner

2096:    Concepts: vector^printing
2097:    Concepts: vector^saving to disk

2099: .seealso: PetscViewerASCIIOpen(), PetscViewerDrawOpen(), PetscDrawLGCreate(),
2100:           PetscViewerSocketOpen(), PetscViewerBinaryOpen(), VecLoad(), PetscViewerCreate(),
2101:           PetscRealView(), PetscScalarView(), PetscIntView()
2102: @*/
2103: int VecView(Vec vec,PetscViewer viewer)
2104: {
2105:   int               ierr;
2106:   PetscViewerFormat format;

2111:   if (!viewer) viewer = PETSC_VIEWER_STDOUT_(vec->comm);
2114:   if (vec->stash.n || vec->bstash.n) SETERRQ(1,"Must call VecAssemblyBegin/End() before viewing this vector");

2116:   /*
2117:      Check if default viewer has been overridden, but user request it anyways
2118:   */
2119:   PetscViewerGetFormat(viewer,&format);
2120:   if (vec->ops->viewnative && format == PETSC_VIEWER_NATIVE) {
2121:     ierr   = PetscViewerPopFormat(viewer);
2122:     (*vec->ops->viewnative)(vec,viewer);
2123:     ierr   = PetscViewerPushFormat(viewer,PETSC_VIEWER_NATIVE);
2124:   } else {
2125:     (*vec->ops->view)(vec,viewer);
2126:   }
2127:   return(0);
2128: }

2130: #undef __FUNCT__  
2132: /*@
2133:    VecGetSize - Returns the global number of elements of the vector.

2135:    Not Collective

2137:    Input Parameter:
2138: .  x - the vector

2140:    Output Parameters:
2141: .  size - the global length of the vector

2143:    Level: beginner

2145:    Concepts: vector^local size

2147: .seealso: VecGetLocalSize()
2148: @*/
2149: int VecGetSize(Vec x,int *size)
2150: {

2157:   (*x->ops->getsize)(x,size);
2158:   return(0);
2159: }

2161: #undef __FUNCT__  
2163: /*@
2164:    VecGetLocalSize - Returns the number of elements of the vector stored 
2165:    in local memory. This routine may be implementation dependent, so use 
2166:    with care.

2168:    Not Collective

2170:    Input Parameter:
2171: .  x - the vector

2173:    Output Parameter:
2174: .  size - the length of the local piece of the vector

2176:    Level: beginner

2178:    Concepts: vector^size

2180: .seealso: VecGetSize()
2181: @*/
2182: int VecGetLocalSize(Vec x,int *size)
2183: {

2190:   (*x->ops->getlocalsize)(x,size);
2191:   return(0);
2192: }

2194: #undef __FUNCT__  
2196: /*@C
2197:    VecGetOwnershipRange - Returns the range of indices owned by 
2198:    this processor, assuming that the vectors are laid out with the
2199:    first n1 elements on the first processor, next n2 elements on the
2200:    second, etc.  For certain parallel layouts this range may not be 
2201:    well defined. 

2203:    Not Collective

2205:    Input Parameter:
2206: .  x - the vector

2208:    Output Parameters:
2209: +  low - the first local element, pass in PETSC_NULL if not interested
2210: -  high - one more than the last local element, pass in PETSC_NULL if not interested

2212:    Note:
2213:    The high argument is one more than the last element stored locally.

2215:    Fortran: PETSC_NULL_INTEGER should be used instead of PETSC_NULL

2217:    Level: beginner

2219:    Concepts: ownership^of vectors
2220:    Concepts: vector^ownership of elements

2222: @*/
2223: int VecGetOwnershipRange(Vec x,int *low,int *high)
2224: {
2225:   int      ierr;

2232:   PetscMapGetLocalRange(x->map,low,high);
2233:   return(0);
2234: }

2236: #undef __FUNCT__  
2238: /*@C
2239:    VecGetPetscMap - Returns the map associated with the vector

2241:    Not Collective

2243:    Input Parameter:
2244: .  x - the vector

2246:    Output Parameters:
2247: .  map - the map

2249:    Level: developer

2251: @*/
2252: int VecGetPetscMap(Vec x,PetscMap *map)
2253: {
2257:   *map = x->map;
2258:   return(0);
2259: }

2261: #undef __FUNCT__  
2263: /*@
2264:    VecSetOption - Sets an option for controling a vector's behavior.

2266:    Collective on Vec

2268:    Input Parameter:
2269: +  x - the vector
2270: -  op - the option

2272:    Supported Options:
2273: +     VEC_IGNORE_OFF_PROC_ENTRIES, which causes VecSetValues() to ignore 
2274:       entries destined to be stored on a seperate processor. This can be used
2275:       to eliminate the global reduction in the VecAssemblyXXXX() if you know 
2276:       that you have only used VecSetValues() to set local elements
2277: -   VEC_TREAT_OFF_PROC_ENTRIES restores the treatment of off processor entries.

2279:    Level: intermediate

2281: @*/
2282: int VecSetOption(Vec x,VecOption op)
2283: {

2289:   if (x->ops->setoption) {
2290:     (*x->ops->setoption)(x,op);
2291:   }
2292:   return(0);
2293: }

2295: #undef __FUNCT__  
2297: /* Default routines for obtaining and releasing; */
2298: /* may be used by any implementation */
2299: int VecDuplicateVecs_Default(Vec w,int m,Vec *V[])
2300: {
2301:   int  i,ierr;

2306:   if (m <= 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"m must be > 0: m = %d",m);
2307:   PetscMalloc(m*sizeof(Vec*),V);
2308:   for (i=0; i<m; i++) {VecDuplicate(w,*V+i);}
2309:   return(0);
2310: }

2312: #undef __FUNCT__  
2314: int VecDestroyVecs_Default(const Vec v[], int m)
2315: {
2316:   int i,ierr;

2320:   if (m <= 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"m must be > 0: m = %d",m);
2321:   for (i=0; i<m; i++) {VecDestroy(v[i]);}
2322:   PetscFree((Vec*)v);
2323:   return(0);
2324: }

2326: #undef __FUNCT__  
2328: /*@
2329:    VecPlaceArray - Allows one to replace the array in a vector with an
2330:    array provided by the user. This is useful to avoid copying an array
2331:    into a vector.

2333:    Not Collective

2335:    Input Parameters:
2336: +  vec - the vector
2337: -  array - the array

2339:    Notes:
2340:    You can return to the original array with a call to VecResetArray()

2342:    Level: developer

2344: .seealso: VecGetArray(), VecRestoreArray(), VecReplaceArray(), VecResetArray()

2346: @*/
2347: int VecPlaceArray(Vec vec,const PetscScalar array[])
2348: {

2354:   if (vec->ops->placearray) {
2355:     (*vec->ops->placearray)(vec,array);
2356:   } else {
2357:     SETERRQ(1,"Cannot place array in this type of vector");
2358:   }
2359:   return(0);
2360: }

2362: #undef __FUNCT__  
2364: /*@
2365:    VecResetArray - Resets a vector to use its default memory. Call this 
2366:    after the use of VecPlaceArray().

2368:    Not Collective

2370:    Input Parameters:
2371: .  vec - the vector

2373:    Level: developer

2375: .seealso: VecGetArray(), VecRestoreArray(), VecReplaceArray(), VecPlaceArray()

2377: @*/
2378: int VecResetArray(Vec vec)
2379: {

2385:   if (vec->ops->resetarray) {
2386:     (*vec->ops->resetarray)(vec);
2387:   } else {
2388:     SETERRQ(1,"Cannot reset array in this type of vector");
2389:   }
2390:   return(0);
2391: }

2393: #undef __FUNCT__  
2395: /*@C
2396:    VecReplaceArray - Allows one to replace the array in a vector with an
2397:    array provided by the user. This is useful to avoid copying an array
2398:    into a vector.

2400:    Not Collective

2402:    Input Parameters:
2403: +  vec - the vector
2404: -  array - the array

2406:    Notes:
2407:    This permanently replaces the array and frees the memory associated
2408:    with the old array.

2410:    The memory passed in MUST be obtained with PetscMalloc() and CANNOT be
2411:    freed by the user. It will be freed when the vector is destroy. 

2413:    Not supported from Fortran

2415:    Level: developer

2417: .seealso: VecGetArray(), VecRestoreArray(), VecPlaceArray(), VecResetArray()

2419: @*/
2420: int VecReplaceArray(Vec vec,const PetscScalar array[])
2421: {

2427:   if (vec->ops->replacearray) {
2428:     (*vec->ops->replacearray)(vec,array);
2429:   } else {
2430:     SETERRQ(1,"Cannot replace array in this type of vector");
2431:   }
2432:   return(0);
2433: }

2435: /*MC
2436:     VecDuplicateVecsF90 - Creates several vectors of the same type as an existing vector
2437:     and makes them accessible via a Fortran90 pointer.

2439:     Synopsis:
2440:     VecDuplicateVecsF90(Vec x,int n,{Vec, pointer :: y(:)},integer ierr)

2442:     Collective on Vec

2444:     Input Parameters:
2445: +   x - a vector to mimic
2446: -   n - the number of vectors to obtain

2448:     Output Parameters:
2449: +   y - Fortran90 pointer to the array of vectors
2450: -   ierr - error code

2452:     Example of Usage: 
2453: .vb
2454:     Vec x
2455:     Vec, pointer :: y(:)
2456:     ....
2457:     call VecDuplicateVecsF90(x,2,y,ierr)
2458:     call VecSet(alpha,y(2),ierr)
2459:     call VecSet(alpha,y(2),ierr)
2460:     ....
2461:     call VecDestroyVecsF90(y,2,ierr)
2462: .ve

2464:     Notes:
2465:     Not yet supported for all F90 compilers

2467:     Use VecDestroyVecsF90() to free the space.

2469:     Level: beginner

2471: .seealso:  VecDestroyVecsF90(), VecDuplicateVecs()

2473: M*/

2475: /*MC
2476:     VecRestoreArrayF90 - Restores a vector to a usable state after a call to
2477:     VecGetArrayF90().

2479:     Synopsis:
2480:     VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2482:     Not collective

2484:     Input Parameters:
2485: +   x - vector
2486: -   xx_v - the Fortran90 pointer to the array

2488:     Output Parameter:
2489: .   ierr - error code

2491:     Example of Usage: 
2492: .vb
2493:     PetscScalar, pointer :: xx_v(:)
2494:     ....
2495:     call VecGetArrayF90(x,xx_v,ierr)
2496:     a = xx_v(3)
2497:     call VecRestoreArrayF90(x,xx_v,ierr)
2498: .ve
2499:    
2500:     Notes:
2501:     Not yet supported for all F90 compilers

2503:     Level: beginner

2505: .seealso:  VecGetArrayF90(), VecGetArray(), VecRestoreArray()

2507: M*/

2509: /*MC
2510:     VecDestroyVecsF90 - Frees a block of vectors obtained with VecDuplicateVecsF90().

2512:     Synopsis:
2513:     VecDestroyVecsF90({Vec, pointer :: x(:)},integer n,integer ierr)

2515:     Input Parameters:
2516: +   x - pointer to array of vector pointers
2517: -   n - the number of vectors previously obtained

2519:     Output Parameter:
2520: .   ierr - error code

2522:     Notes:
2523:     Not yet supported for all F90 compilers

2525:     Level: beginner

2527: .seealso:  VecDestroyVecs(), VecDuplicateVecsF90()

2529: M*/

2531: /*MC
2532:     VecGetArrayF90 - Accesses a vector array from Fortran90. For default PETSc
2533:     vectors, VecGetArrayF90() returns a pointer to the local data array. Otherwise,
2534:     this routine is implementation dependent. You MUST call VecRestoreArrayF90() 
2535:     when you no longer need access to the array.

2537:     Synopsis:
2538:     VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

2540:     Not Collective 

2542:     Input Parameter:
2543: .   x - vector

2545:     Output Parameters:
2546: +   xx_v - the Fortran90 pointer to the array
2547: -   ierr - error code

2549:     Example of Usage: 
2550: .vb
2551:     PetscScalar, pointer :: xx_v(:)
2552:     ....
2553:     call VecGetArrayF90(x,xx_v,ierr)
2554:     a = xx_v(3)
2555:     call VecRestoreArrayF90(x,xx_v,ierr)
2556: .ve

2558:     Notes:
2559:     Not yet supported for all F90 compilers

2561:     Level: beginner

2563: .seealso:  VecRestoreArrayF90(), VecGetArray(), VecRestoreArray()

2565: M*/

2567: #undef __FUNCT__  
2569: /*@C 
2570:   VecLoadIntoVector - Loads a vector that has been stored in binary format
2571:   with VecView().

2573:   Collective on PetscViewer 

2575:   Input Parameters:
2576: + viewer - binary file viewer, obtained from PetscViewerBinaryOpen()
2577: - vec - vector to contain files values (must be of correct length)

2579:   Level: intermediate

2581:   Notes:
2582:   The input file must contain the full global vector, as
2583:   written by the routine VecView().

2585:   Use VecLoad() to create the vector as the values are read in

2587:   Notes for advanced users:
2588:   Most users should not need to know the details of the binary storage
2589:   format, since VecLoad() and VecView() completely hide these details.
2590:   But for anyone who's interested, the standard binary matrix storage
2591:   format is
2592: .vb
2593:      int    VEC_FILE_COOKIE
2594:      int    number of rows
2595:      PetscScalar *values of all nonzeros
2596: .ve

2598:    Note for Cray users, the int's stored in the binary file are 32 bit
2599: integers; not 64 as they are represented in the memory, so if you
2600: write your own routines to read/write these binary files from the Cray
2601: you need to adjust the integer sizes that you read in, see
2602: PetscReadBinary() and PetscWriteBinary() to see how this may be
2603: done.

2605:    In addition, PETSc automatically does the byte swapping for
2606: machines that store the bytes reversed, e.g.  DEC alpha, freebsd,
2607: linux, nt and the paragon; thus if you write your own binary
2608: read/write routines you have to swap the bytes; see PetscReadBinary()
2609: and PetscWriteBinary() to see how this may be done.

2611:    Concepts: vector^loading from file

2613: .seealso: PetscViewerBinaryOpen(), VecView(), MatLoad(), VecLoad() 
2614: @*/
2615: int VecLoadIntoVector(PetscViewer viewer,Vec vec)
2616: {

2623:   if (!vec->ops->loadintovector) {
2624:     SETERRQ(1,"Vector does not support load");
2625:   }
2626:   (*vec->ops->loadintovector)(viewer,vec);
2627:   return(0);
2628: }

2630: #undef __FUNCT__  
2632: /*@
2633:    VecReciprocal - Replaces each component of a vector by its reciprocal.

2635:    Collective on Vec

2637:    Input Parameter:
2638: .  v - the vector 

2640:    Output Parameter:
2641: .  v - the vector reciprocal

2643:    Level: intermediate

2645:    Concepts: vector^reciprocal

2647: @*/
2648: int VecReciprocal(Vec vec)
2649: {
2650:   int    ierr;

2655:   if (!vec->ops->reciprocal) {
2656:     SETERRQ(1,"Vector does not support reciprocal operation");
2657:   }
2658:   (*vec->ops->reciprocal)(vec);
2659:   return(0);
2660: }

2662: #undef __FUNCT__  
2664: int VecSetOperation(Vec vec,VecOperation op, void (*f)(void))
2665: {
2668:   /* save the native version of the viewer */
2669:   if (op == VECOP_VIEW && !vec->ops->viewnative) {
2670:     vec->ops->viewnative = vec->ops->view;
2671:   }
2672:   (((void(**)(void))vec->ops)[(int)op]) = f;
2673:   return(0);
2674: }

2676: #undef __FUNCT__  
2678: /*@
2679:    VecSetStashInitialSize - sets the sizes of the vec-stash, that is
2680:    used during the assembly process to store values that belong to 
2681:    other processors.

2683:    Collective on Vec

2685:    Input Parameters:
2686: +  vec   - the vector
2687: .  size  - the initial size of the stash.
2688: -  bsize - the initial size of the block-stash(if used).

2690:    Options Database Keys:
2691: +   -vecstash_initial_size <size> or <size0,size1,...sizep-1>
2692: -   -vecstash_block_initial_size <bsize> or <bsize0,bsize1,...bsizep-1>

2694:    Level: intermediate

2696:    Notes: 
2697:      The block-stash is used for values set with VecSetValuesBlocked() while
2698:      the stash is used for values set with VecSetValues()

2700:      Run with the option -log_info and look for output of the form
2701:      VecAssemblyBegin_MPIXXX:Stash has MM entries, uses nn mallocs.
2702:      to determine the appropriate value, MM, to use for size and 
2703:      VecAssemblyBegin_MPIXXX:Block-Stash has BMM entries, uses nn mallocs.
2704:      to determine the value, BMM to use for bsize

2706:    Concepts: vector^stash
2707:    Concepts: stash^vector

2709: .seealso: VecSetBlockSize(), VecSetValues(), VecSetValuesBlocked(), VecStashView()

2711: @*/
2712: int VecSetStashInitialSize(Vec vec,int size,int bsize)
2713: {

2718:   VecStashSetInitialSize_Private(&vec->stash,size);
2719:   VecStashSetInitialSize_Private(&vec->bstash,bsize);
2720:   return(0);
2721: }

2723: #undef __FUNCT__  
2725: /*@
2726:    VecStashView - Prints the entries in the vector stash and block stash.

2728:    Collective on Vec

2730:    Input Parameters:
2731: +  vec   - the vector
2732: -  viewer - the viewer

2734:    Level: advanced

2736:    Concepts: vector^stash
2737:    Concepts: stash^vector

2739: .seealso: VecSetBlockSize(), VecSetValues(), VecSetValuesBlocked()

2741: @*/
2742: int VecStashView(Vec v,PetscViewer viewer)
2743: {
2744:   int          ierr,rank,i,j;
2745:   PetscTruth   match;
2746:   VecStash     *s;
2747:   PetscScalar  val;


2754:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&match);
2755:   if (!match) SETERRQ1(1,"Stash viewer only works with ASCII viewer not %sn",((PetscObject)v)->type_name);
2756:   PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
2757:   MPI_Comm_rank(v->comm,&rank);
2758:   s = &v->bstash;

2760:   /* print block stash */
2761:   PetscViewerASCIISynchronizedPrintf(viewer,"[%d]Vector Block stash size %d block size %dn",rank,s->n,s->bs);
2762:   for (i=0; i<s->n; i++) {
2763:     PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %d ",rank,s->idx[i]);
2764:     for (j=0; j<s->bs; j++) {
2765:       val = s->array[i*s->bs+j];
2766: #if defined(PETSC_USE_COMPLEX)
2767:       PetscViewerASCIISynchronizedPrintf(viewer,"(%18.16e %18.16e) ",PetscRealPart(val),PetscImaginaryPart(val));
2768: #else
2769:       PetscViewerASCIISynchronizedPrintf(viewer,"%18.16e ",val);
2770: #endif
2771:     }
2772:     PetscViewerASCIISynchronizedPrintf(viewer,"n");
2773:   }
2774:   PetscViewerFlush(viewer);

2776:   s = &v->stash;

2778:   /* print basic stash */
2779:   PetscViewerASCIISynchronizedPrintf(viewer,"[%d]Vector stash size %dn",rank,s->n);
2780:   for (i=0; i<s->n; i++) {
2781:     val = s->array[i];
2782: #if defined(PETSC_USE_COMPLEX)
2783:       PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %d (%18.16e %18.16e) ",rank,s->idx[i],PetscRealPart(val),PetscImaginaryPart(val));
2784: #else
2785:     PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %d %18.16en",rank,s->idx[i],val);
2786: #endif
2787:   }
2788:   PetscViewerFlush(viewer);

2790:   PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
2791:   return(0);
2792: }

2794: #undef __FUNCT__  
2796: /*@C
2797:    VecGetArray2d - Returns a pointer to a 2d contiguous array that contains this 
2798:    processor's portion of the vector data.  You MUST call VecRestoreArray2d() 
2799:    when you no longer need access to the array.

2801:    Not Collective

2803:    Input Parameter:
2804: +  x - the vector
2805: .  m - first dimension of two dimensional array
2806: .  n - second dimension of two dimensional array
2807: .  mstart - first index you will use in first coordinate direction (often 0)
2808: -  nstart - first index in the second coordinate direction (often 0)

2810:    Output Parameter:
2811: .  a - location to put pointer to the array

2813:    Level: beginner

2815:   Notes:
2816:    For a vector obtained from DACreateLocalVector() mstart and nstart are likely
2817:    obtained from the corner indices obtained from DAGetGhostCorners() while for
2818:    DACreateGlobalVector() they are the corner indices from DAGetCorners(). In both cases
2819:    the arguments from DAGet[Ghost}Corners() are reversed in the call to VecGetArray2d().
2820:    
2821:    For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.

2823:    Concepts: vector^accessing local values as 2d array

2825: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2826:           VecRestoreArray2d(), DAVecGetarray(), DAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2827:           VecGetarray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2828: @*/
2829: int VecGetArray2d(Vec x,int m,int n,int mstart,int nstart,PetscScalar **a[])
2830: {
2831:   int         i,ierr,N;
2832:   PetscScalar *aa;

2838:   VecGetLocalSize(x,&N);
2839:   if (m*n != N) SETERRQ3(1,"Local array size %d does not match 2d array dimensions %d by %d",N,m,n);
2840:   VecGetArray(x,&aa);

2842:   PetscMalloc(m*sizeof(PetscScalar*),a);
2843:   for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
2844:   *a -= mstart;
2845:   return(0);
2846: }

2848: #undef __FUNCT__  
2850: /*@C
2851:    VecRestoreArray2d - Restores a vector after VecGetArray2d() has been called.

2853:    Not Collective

2855:    Input Parameters:
2856: +  x - the vector
2857: .  m - first dimension of two dimensional array
2858: .  n - second dimension of the two dimensional array
2859: .  mstart - first index you will use in first coordinate direction (often 0)
2860: .  nstart - first index in the second coordinate direction (often 0)
2861: -  a - location of pointer to array obtained from VecGetArray2d()

2863:    Level: beginner

2865:    Notes:
2866:    For regular PETSc vectors this routine does not involve any copies. For
2867:    any special vectors that do not store local vector data in a contiguous
2868:    array, this routine will copy the data back into the underlying 
2869:    vector data structure from the array obtained with VecGetArray().

2871:    This routine actually zeros out the a pointer. 

2873: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2874:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DAVecGetArray(), DAVecRestoreArray()
2875:           VecGetarray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2876: @*/
2877: int VecRestoreArray2d(Vec x,int m,int n,int mstart,int nstart,PetscScalar **a[])
2878: {

2884:   PetscFree(*a + mstart);
2885:   VecRestoreArray(x,PETSC_NULL);
2886:   return(0);
2887: }

2889: #undef __FUNCT__  
2891: /*@C
2892:    VecGetArray1d - Returns a pointer to a 1d contiguous array that contains this 
2893:    processor's portion of the vector data.  You MUST call VecRestoreArray1d() 
2894:    when you no longer need access to the array.

2896:    Not Collective

2898:    Input Parameter:
2899: +  x - the vector
2900: .  m - first dimension of two dimensional array
2901: -  mstart - first index you will use in first coordinate direction (often 0)

2903:    Output Parameter:
2904: .  a - location to put pointer to the array

2906:    Level: beginner

2908:   Notes:
2909:    For a vector obtained from DACreateLocalVector() mstart are likely
2910:    obtained from the corner indices obtained from DAGetGhostCorners() while for
2911:    DACreateGlobalVector() they are the corner indices from DAGetCorners(). 
2912:    
2913:    For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.

2915: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2916:           VecRestoreArray2d(), DAVecGetarray(), DAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2917:           VecGetarray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2918: @*/
2919: int VecGetArray1d(Vec x,int m,int mstart,PetscScalar *a[])
2920: {
2921:   int ierr,N;

2927:   VecGetLocalSize(x,&N);
2928:   if (m != N) SETERRQ2(1,"Local array size %d does not match 1d array dimensions %d",N,m);
2929:   VecGetArray(x,a);
2930:   *a  -= mstart;
2931:   return(0);
2932: }

2934: #undef __FUNCT__  
2936: /*@C
2937:    VecRestoreArray1d - Restores a vector after VecGetArray1d() has been called.

2939:    Not Collective

2941:    Input Parameters:
2942: +  x - the vector
2943: .  m - first dimension of two dimensional array
2944: .  mstart - first index you will use in first coordinate direction (often 0)
2945: -  a - location of pointer to array obtained from VecGetArray21()

2947:    Level: beginner

2949:    Notes:
2950:    For regular PETSc vectors this routine does not involve any copies. For
2951:    any special vectors that do not store local vector data in a contiguous
2952:    array, this routine will copy the data back into the underlying 
2953:    vector data structure from the array obtained with VecGetArray1d().

2955:    This routine actually zeros out the a pointer. 

2957:    Concepts: vector^accessing local values as 1d array

2959: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2960:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DAVecGetArray(), DAVecRestoreArray()
2961:           VecGetarray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
2962: @*/
2963: int VecRestoreArray1d(Vec x,int m,int mstart,PetscScalar *a[])
2964: {

2970:   VecRestoreArray(x,PETSC_NULL);
2971:   return(0);
2972: }

2974: #undef __FUNCT__  
2976: /*@C
2977:    VecConjugate - Conjugates a vector.

2979:    Collective on Vec

2981:    Input Parameters:
2982: .  x - the vector

2984:    Level: intermediate

2986:    Concepts: vector^conjugate

2988: @*/
2989: int VecConjugate(Vec x)
2990: {
2991: #ifdef PETSC_USE_COMPLEX

2997:   (*x->ops->conjugate)(x);
2998:   return(0);
2999: #else
3000:   return(0);
3001: #endif
3002: }

3004: #undef __FUNCT__  
3006: /*@C
3007:    VecGetArray3d - Returns a pointer to a 3d contiguous array that contains this 
3008:    processor's portion of the vector data.  You MUST call VecRestoreArray3d() 
3009:    when you no longer need access to the array.

3011:    Not Collective

3013:    Input Parameter:
3014: +  x - the vector
3015: .  m - first dimension of three dimensional array
3016: .  n - second dimension of three dimensional array
3017: .  p - third dimension of three dimensional array
3018: .  mstart - first index you will use in first coordinate direction (often 0)
3019: .  nstart - first index in the second coordinate direction (often 0)
3020: -  pstart - first index in the third coordinate direction (often 0)

3022:    Output Parameter:
3023: .  a - location to put pointer to the array

3025:    Level: beginner

3027:   Notes:
3028:    For a vector obtained from DACreateLocalVector() mstart, nstart, and pstart are likely
3029:    obtained from the corner indices obtained from DAGetGhostCorners() while for
3030:    DACreateGlobalVector() they are the corner indices from DAGetCorners(). In both cases
3031:    the arguments from DAGet[Ghost}Corners() are reversed in the call to VecGetArray3d().
3032:    
3033:    For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.

3035:    Concepts: vector^accessing local values as 3d array

3037: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
3038:           VecRestoreArray2d(), DAVecGetarray(), DAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
3039:           VecGetarray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
3040: @*/
3041: int VecGetArray3d(Vec x,int m,int n,int p,int mstart,int nstart,int pstart,PetscScalar ***a[])
3042: {
3043:   int         i,ierr,N,j;
3044:   PetscScalar *aa,**b;

3050:   VecGetLocalSize(x,&N);
3051:   if (m*n*p != N) SETERRQ4(1,"Local array size %d does not match 3d array dimensions %d by %d by %d",N,m,n,p);
3052:   VecGetArray(x,&aa);

3054:   PetscMalloc(m*sizeof(PetscScalar**)+m*n*sizeof(PetscScalar*),a);
3055:   b    = (PetscScalar **)((*a) + m);
3056:   for (i=0; i<m; i++)   (*a)[i] = b + i*n - nstart;
3057:   for (i=0; i<m; i++) {
3058:     for (j=0; j<n; j++) {
3059:       b[i*n+j] = aa + i*n*p + j*p - pstart;
3060:     }
3061:   }
3062:   *a -= mstart;
3063:   return(0);
3064: }

3066: #undef __FUNCT__  
3068: /*@C
3069:    VecRestoreArray3d - Restores a vector after VecGetArray3d() has been called.

3071:    Not Collective

3073:    Input Parameters:
3074: +  x - the vector
3075: .  m - first dimension of three dimensional array
3076: .  n - second dimension of the three dimensional array
3077: .  p - third dimension of the three dimensional array
3078: .  mstart - first index you will use in first coordinate direction (often 0)
3079: .  nstart - first index in the second coordinate direction (often 0)
3080: .  pstart - first index in the third coordinate direction (often 0)
3081: -  a - location of pointer to array obtained from VecGetArray3d()

3083:    Level: beginner

3085:    Notes:
3086:    For regular PETSc vectors this routine does not involve any copies. For
3087:    any special vectors that do not store local vector data in a contiguous
3088:    array, this routine will copy the data back into the underlying 
3089:    vector data structure from the array obtained with VecGetArray().

3091:    This routine actually zeros out the a pointer. 

3093: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
3094:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DAVecGetArray(), DAVecRestoreArray()
3095:           VecGetarray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
3096: @*/
3097: int VecRestoreArray3d(Vec x,int m,int n,int p,int mstart,int nstart,int pstart,PetscScalar ***a[])
3098: {

3104:   PetscFree(*a + mstart);
3105:   VecRestoreArray(x,PETSC_NULL);
3106:   return(0);
3107: }