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