Actual source code: rvector.c
petsc-3.4.5 2014-06-29
2: /*
3: Provides the interface functions for vector operations that have PetscScalar/PetscReal in the signature
4: These are the vector functions the user calls.
5: */
6: #include <petsc-private/vecimpl.h> /*I "petscvec.h" I*/
7: static PetscInt VecGetSubVectorSavedStateId = -1;
10: if ((x)->map->N != (y)->map->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths %d != %d", (x)->map->N, (y)->map->N); \
11: if ((x)->map->n != (y)->map->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths %d != %d", (x)->map->n, (y)->map->n);
15: PETSC_EXTERN PetscErrorCode VecValidValues(Vec vec,PetscInt argnum,PetscBool begin)
16: {
17: #if defined(PETSC_USE_DEBUG)
18: PetscErrorCode ierr;
19: PetscInt n,i;
20: const PetscScalar *x;
23: if (vec->petscnative || vec->ops->getarray) {
24: VecGetLocalSize(vec,&n);
25: VecGetArrayRead(vec,&x);
26: for (i=0; i<n; i++) {
27: if (begin) {
28: if (PetscIsInfOrNanScalar(x[i])) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FP,"Vec entry at local location %D is not-a-number or infinite at beginning of function: Parameter number %D",i,argnum);
29: } else {
30: if (PetscIsInfOrNanScalar(x[i])) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FP,"Vec entry at local location %D is not-a-number or infinite at end of function: Parameter number %D",i,argnum);
31: }
32: }
33: VecRestoreArrayRead(vec,&x);
34: }
35: return(0);
36: #else
37: return 0;
38: #endif
39: }
43: /*@
44: VecMaxPointwiseDivide - Computes the maximum of the componentwise division max = max_i abs(x_i/y_i).
46: Logically Collective on Vec
48: Input Parameters:
49: . x, y - the vectors
51: Output Parameter:
52: . max - the result
54: Level: advanced
56: Notes: x and y may be the same vector
57: if a particular y_i is zero, it is treated as 1 in the above formula
59: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs()
60: @*/
61: PetscErrorCode VecMaxPointwiseDivide(Vec x,Vec y,PetscReal *max)
62: {
74: (*x->ops->maxpointwisedivide)(x,y,max);
75: return(0);
76: }
80: /*@
81: VecDot - Computes the vector dot product.
83: Collective on Vec
85: Input Parameters:
86: . x, y - the vectors
88: Output Parameter:
89: . val - the dot product
91: Performance Issues:
92: $ per-processor memory bandwidth
93: $ interprocessor latency
94: $ work load inbalance that causes certain processes to arrive much earlier than others
96: Notes for Users of Complex Numbers:
97: For complex vectors, VecDot() computes
98: $ val = (x,y) = y^H x,
99: where y^H denotes the conjugate transpose of y. Note that this corresponds to the usual "mathematicians" complex
100: inner product where the SECOND argument gets the complex conjugate. Since the BLASdot() complex conjugates the first
101: first argument we call the BLASdot() with the arguments reversed.
103: Use VecTDot() for the indefinite form
104: $ val = (x,y) = y^T x,
105: where y^T denotes the transpose of y.
107: Level: intermediate
109: Concepts: inner product
110: Concepts: vector^inner product
112: .seealso: VecMDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecDotRealPart()
113: @*/
114: PetscErrorCode VecDot(Vec x,Vec y,PetscScalar *val)
115: {
127: PetscLogEventBarrierBegin(VEC_DotBarrier,x,y,0,0,PetscObjectComm((PetscObject)x));
128: (*x->ops->dot)(x,y,val);
129: PetscLogEventBarrierEnd(VEC_DotBarrier,x,y,0,0,PetscObjectComm((PetscObject)x));
130: return(0);
131: }
135: /*@
136: VecDotRealPart - Computes the real part of the vector dot product.
138: Collective on Vec
140: Input Parameters:
141: . x, y - the vectors
143: Output Parameter:
144: . val - the real part of the dot product;
146: Performance Issues:
147: $ per-processor memory bandwidth
148: $ interprocessor latency
149: $ work load inbalance that causes certain processes to arrive much earlier than others
151: Notes for Users of Complex Numbers:
152: See VecDot() for more details on the definition of the dot product for complex numbers
154: For real numbers this returns the same value as VecDot()
156: For complex numbers in C^n (that is a vector of n components with a complex number for each component) this is equal to the usual real dot product on the
157: the space R^{2n} (that is a vector of 2n components with the real or imaginary part of the complex numbers for components)
159: Developer Note: This is not currently optimized to compute only the real part of the dot product.
161: Level: intermediate
163: Concepts: inner product
164: Concepts: vector^inner product
166: .seealso: VecMDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecDot(), VecDotNorm2()
167: @*/
168: PetscErrorCode VecDotRealPart(Vec x,Vec y,PetscReal *val)
169: {
171: PetscScalar fdot;
174: VecDot(x,y,&fdot);
175: *val = PetscRealPart(fdot);
176: return(0);
177: }
181: /*@
182: VecNorm - Computes the vector norm.
184: Collective on Vec
186: Input Parameters:
187: + x - the vector
188: - type - one of NORM_1, NORM_2, NORM_INFINITY. Also available
189: NORM_1_AND_2, which computes both norms and stores them
190: in a two element array.
192: Output Parameter:
193: . val - the norm
195: Notes:
196: $ NORM_1 denotes sum_i |x_i|
197: $ NORM_2 denotes sqrt(sum_i (x_i)^2)
198: $ NORM_INFINITY denotes max_i |x_i|
200: Level: intermediate
202: Performance Issues:
203: $ per-processor memory bandwidth
204: $ interprocessor latency
205: $ work load inbalance that causes certain processes to arrive much earlier than others
207: Compile Option:
208: PETSC_HAVE_SLOW_BLAS_NORM2 will cause a C (loop unrolled) version of the norm to be used, rather
209: than the BLAS. This should probably only be used when one is using the FORTRAN BLAS routines
210: (as opposed to vendor provided) because the FORTRAN BLAS NRM2() routine is very slow.
212: Concepts: norm
213: Concepts: vector^norm
215: .seealso: VecDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecNormAvailable(),
216: VecNormBegin(), VecNormEnd()
218: @*/
219: PetscErrorCode VecNorm(Vec x,NormType type,PetscReal *val)
220: {
221: PetscBool flg;
228: if (((PetscObject)x)->precision != sizeof(PetscReal)) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Wrong precision of input argument");
230: /*
231: * Cached data?
232: */
233: if (type!=NORM_1_AND_2) {
234: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[type],*val,flg);
235: if (flg) return(0);
236: }
237: PetscLogEventBarrierBegin(VEC_NormBarrier,x,0,0,0,PetscObjectComm((PetscObject)x));
238: (*x->ops->norm)(x,type,val);
239: PetscLogEventBarrierEnd(VEC_NormBarrier,x,0,0,0,PetscObjectComm((PetscObject)x));
241: if (type!=NORM_1_AND_2) {
242: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[type],*val);
243: }
244: return(0);
245: }
249: /*@
250: VecNormAvailable - Returns the vector norm if it is already known.
252: Not Collective
254: Input Parameters:
255: + x - the vector
256: - type - one of NORM_1, NORM_2, NORM_INFINITY. Also available
257: NORM_1_AND_2, which computes both norms and stores them
258: in a two element array.
260: Output Parameter:
261: + available - PETSC_TRUE if the val returned is valid
262: - val - the norm
264: Notes:
265: $ NORM_1 denotes sum_i |x_i|
266: $ NORM_2 denotes sqrt(sum_i (x_i)^2)
267: $ NORM_INFINITY denotes max_i |x_i|
269: Level: intermediate
271: Performance Issues:
272: $ per-processor memory bandwidth
273: $ interprocessor latency
274: $ work load inbalance that causes certain processes to arrive much earlier than others
276: Compile Option:
277: PETSC_HAVE_SLOW_BLAS_NORM2 will cause a C (loop unrolled) version of the norm to be used, rather
278: than the BLAS. This should probably only be used when one is using the FORTRAN BLAS routines
279: (as opposed to vendor provided) because the FORTRAN BLAS NRM2() routine is very slow.
281: Concepts: norm
282: Concepts: vector^norm
284: .seealso: VecDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecNorm()
285: VecNormBegin(), VecNormEnd()
287: @*/
288: PetscErrorCode VecNormAvailable(Vec x,NormType type,PetscBool *available,PetscReal *val)
289: {
297: *available = PETSC_FALSE;
298: if (type!=NORM_1_AND_2) {
299: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[type],*val,*available);
300: }
301: return(0);
302: }
306: /*@
307: VecNormalize - Normalizes a vector by 2-norm.
309: Collective on Vec
311: Input Parameters:
312: + x - the vector
314: Output Parameter:
315: . x - the normalized vector
316: - val - the vector norm before normalization
318: Level: intermediate
320: Concepts: vector^normalizing
321: Concepts: normalizing^vector
323: @*/
324: PetscErrorCode VecNormalize(Vec x,PetscReal *val)
325: {
327: PetscReal norm;
332: PetscLogEventBegin(VEC_Normalize,x,0,0,0);
333: VecNorm(x,NORM_2,&norm);
334: if (norm == 0.0) {
335: PetscInfo(x,"Vector of zero norm can not be normalized; Returning only the zero norm\n");
336: } else if (norm != 1.0) {
337: PetscScalar tmp = 1.0/norm;
338: VecScale(x,tmp);
339: }
340: if (val) *val = norm;
341: PetscLogEventEnd(VEC_Normalize,x,0,0,0);
342: return(0);
343: }
347: /*@C
348: VecMax - Determines the maximum vector component and its location.
350: Collective on Vec
352: Input Parameter:
353: . x - the vector
355: Output Parameters:
356: + val - the maximum component
357: - p - the location of val (pass NULL if you don't want this)
359: Notes:
360: Returns the value PETSC_MIN_REAL and p = -1 if the vector is of length 0.
362: Returns the smallest index with the maximum value
363: Level: intermediate
365: Concepts: maximum^of vector
366: Concepts: vector^maximum value
368: .seealso: VecNorm(), VecMin()
369: @*/
370: PetscErrorCode VecMax(Vec x,PetscInt *p,PetscReal *val)
371: {
378: PetscLogEventBegin(VEC_Max,x,0,0,0);
379: (*x->ops->max)(x,p,val);
380: PetscLogEventEnd(VEC_Max,x,0,0,0);
381: return(0);
382: }
386: /*@
387: VecMin - Determines the minimum vector component and its location.
389: Collective on Vec
391: Input Parameters:
392: . x - the vector
394: Output Parameter:
395: + val - the minimum component
396: - p - the location of val (pass NULL if you don't want this location)
398: Level: intermediate
400: Notes:
401: Returns the value PETSC_MAX_REAL and p = -1 if the vector is of length 0.
403: This returns the smallest index with the minumum value
405: Concepts: minimum^of vector
406: Concepts: vector^minimum entry
408: .seealso: VecMax()
409: @*/
410: PetscErrorCode VecMin(Vec x,PetscInt *p,PetscReal *val)
411: {
418: PetscLogEventBegin(VEC_Min,x,0,0,0);
419: (*x->ops->min)(x,p,val);
420: PetscLogEventEnd(VEC_Min,x,0,0,0);
421: return(0);
422: }
426: /*@
427: VecTDot - Computes an indefinite vector dot product. That is, this
428: routine does NOT use the complex conjugate.
430: Collective on Vec
432: Input Parameters:
433: . x, y - the vectors
435: Output Parameter:
436: . val - the dot product
438: Notes for Users of Complex Numbers:
439: For complex vectors, VecTDot() computes the indefinite form
440: $ val = (x,y) = y^T x,
441: where y^T denotes the transpose of y.
443: Use VecDot() for the inner product
444: $ val = (x,y) = y^H x,
445: where y^H denotes the conjugate transpose of y.
447: Level: intermediate
449: Concepts: inner product^non-Hermitian
450: Concepts: vector^inner product
451: Concepts: non-Hermitian inner product
453: .seealso: VecDot(), VecMTDot()
454: @*/
455: PetscErrorCode VecTDot(Vec x,Vec y,PetscScalar *val)
456: {
468: PetscLogEventBegin(VEC_TDot,x,y,0,0);
469: (*x->ops->tdot)(x,y,val);
470: PetscLogEventEnd(VEC_TDot,x,y,0,0);
471: return(0);
472: }
476: /*@
477: VecScale - Scales a vector.
479: Not collective on Vec
481: Input Parameters:
482: + x - the vector
483: - alpha - the scalar
485: Output Parameter:
486: . x - the scaled vector
488: Note:
489: For a vector with n components, VecScale() computes
490: $ x[i] = alpha * x[i], for i=1,...,n.
492: Level: intermediate
494: Concepts: vector^scaling
495: Concepts: scaling^vector
497: @*/
498: PetscErrorCode VecScale(Vec x, PetscScalar alpha)
499: {
500: PetscReal norms[4] = {0.0,0.0,0.0, 0.0};
501: PetscBool flgs[4];
503: PetscInt i;
508: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
509: PetscLogEventBegin(VEC_Scale,x,0,0,0);
510: if (alpha != (PetscScalar)1.0) {
511: /* get current stashed norms */
512: for (i=0; i<4; i++) {
513: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],norms[i],flgs[i]);
514: }
515: (*x->ops->scale)(x,alpha);
516: PetscObjectStateIncrease((PetscObject)x);
517: /* put the scaled stashed norms back into the Vec */
518: for (i=0; i<4; i++) {
519: if (flgs[i]) {
520: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[i],PetscAbsScalar(alpha)*norms[i]);
521: }
522: }
523: }
524: PetscLogEventEnd(VEC_Scale,x,0,0,0);
525: return(0);
526: }
530: /*@
531: VecSet - Sets all components of a vector to a single scalar value.
533: Logically Collective on Vec
535: Input Parameters:
536: + x - the vector
537: - alpha - the scalar
539: Output Parameter:
540: . x - the vector
542: Note:
543: For a vector of dimension n, VecSet() computes
544: $ x[i] = alpha, for i=1,...,n,
545: so that all vector entries then equal the identical
546: scalar value, alpha. Use the more general routine
547: VecSetValues() to set different vector entries.
549: You CANNOT call this after you have called VecSetValues() but before you call
550: VecAssemblyBegin/End().
552: Level: beginner
554: .seealso VecSetValues(), VecSetValuesBlocked(), VecSetRandom()
556: Concepts: vector^setting to constant
558: @*/
559: PetscErrorCode VecSet(Vec x,PetscScalar alpha)
560: {
561: PetscReal val;
567: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"You cannot call this after you have called VecSetValues() but\n before you have called VecAssemblyBegin/End()");
570: PetscLogEventBegin(VEC_Set,x,0,0,0);
571: (*x->ops->set)(x,alpha);
572: PetscLogEventEnd(VEC_Set,x,0,0,0);
573: PetscObjectStateIncrease((PetscObject)x);
575: /* norms can be simply set */
576: val = PetscAbsScalar(alpha);
577: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_1],x->map->N * val);
578: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_INFINITY],val);
579: val = PetscSqrtReal((double)x->map->N) * val;
580: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_2],val);
581: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_FROBENIUS],val);
582: return(0);
583: }
588: /*@
589: VecAXPY - Computes y = alpha x + y.
591: Logically Collective on Vec
593: Input Parameters:
594: + alpha - the scalar
595: - x, y - the vectors
597: Output Parameter:
598: . y - output vector
600: Level: intermediate
602: Notes: x and y MUST be different vectors
604: Concepts: vector^BLAS
605: Concepts: BLAS
607: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY()
608: @*/
609: PetscErrorCode VecAXPY(Vec y,PetscScalar alpha,Vec x)
610: {
620: if (x == y) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x and y cannot be the same vector");
623: PetscLogEventBegin(VEC_AXPY,x,y,0,0);
624: (*y->ops->axpy)(y,alpha,x);
625: PetscLogEventEnd(VEC_AXPY,x,y,0,0);
626: PetscObjectStateIncrease((PetscObject)y);
627: return(0);
628: }
632: /*@
633: VecAXPBY - Computes y = alpha x + beta y.
635: Logically Collective on Vec
637: Input Parameters:
638: + alpha,beta - the scalars
639: - x, y - the vectors
641: Output Parameter:
642: . y - output vector
644: Level: intermediate
646: Notes: x and y MUST be different vectors
648: Concepts: BLAS
649: Concepts: vector^BLAS
651: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPY()
652: @*/
653: PetscErrorCode VecAXPBY(Vec y,PetscScalar alpha,PetscScalar beta,Vec x)
654: {
664: if (x == y) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x and y cannot be the same vector");
668: PetscLogEventBegin(VEC_AXPY,x,y,0,0);
669: (*y->ops->axpby)(y,alpha,beta,x);
670: PetscLogEventEnd(VEC_AXPY,x,y,0,0);
671: PetscObjectStateIncrease((PetscObject)y);
672: return(0);
673: }
677: /*@
678: VecAXPBYPCZ - Computes z = alpha x + beta y + gamma z
680: Logically Collective on Vec
682: Input Parameters:
683: + alpha,beta, gamma - the scalars
684: - x, y, z - the vectors
686: Output Parameter:
687: . z - output vector
689: Level: intermediate
691: Notes: x, y and z must be different vectors
693: Developer Note: alpha = 1 or gamma = 1 or gamma = 0.0 are handled as special cases
695: Concepts: BLAS
696: Concepts: vector^BLAS
698: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPY()
699: @*/
700: PetscErrorCode VecAXPBYPCZ(Vec z,PetscScalar alpha,PetscScalar beta,PetscScalar gamma,Vec x,Vec y)
701: {
715: if (x == y || x == z) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x, y, and z must be different vectors");
716: if (y == z) SETERRQ(PetscObjectComm((PetscObject)y),PETSC_ERR_ARG_IDN,"x, y, and z must be different vectors");
721: PetscLogEventBegin(VEC_AXPBYPCZ,x,y,z,0);
722: (*y->ops->axpbypcz)(z,alpha,beta,gamma,x,y);
723: PetscLogEventEnd(VEC_AXPBYPCZ,x,y,z,0);
724: PetscObjectStateIncrease((PetscObject)z);
725: return(0);
726: }
730: /*@
731: VecAYPX - Computes y = x + alpha y.
733: Logically Collective on Vec
735: Input Parameters:
736: + alpha - the scalar
737: - x, y - the vectors
739: Output Parameter:
740: . y - output vector
742: Level: intermediate
744: Notes: x and y MUST be different vectors
746: Concepts: vector^BLAS
747: Concepts: BLAS
749: .seealso: VecAXPY(), VecWAXPY()
750: @*/
751: PetscErrorCode VecAYPX(Vec y,PetscScalar alpha,Vec x)
752: {
760: if (x == y) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
763: PetscLogEventBegin(VEC_AYPX,x,y,0,0);
764: (*y->ops->aypx)(y,alpha,x);
765: PetscLogEventEnd(VEC_AYPX,x,y,0,0);
766: PetscObjectStateIncrease((PetscObject)y);
767: return(0);
768: }
773: /*@
774: VecWAXPY - Computes w = alpha x + y.
776: Logically Collective on Vec
778: Input Parameters:
779: + alpha - the scalar
780: - x, y - the vectors
782: Output Parameter:
783: . w - the result
785: Level: intermediate
787: Notes: w cannot be either x or y, but x and y can be the same
789: Concepts: vector^BLAS
790: Concepts: BLAS
792: .seealso: VecAXPY(), VecAYPX(), VecAXPBY()
793: @*/
794: PetscErrorCode VecWAXPY(Vec w,PetscScalar alpha,Vec x,Vec y)
795: {
809: if (w == y) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Result vector w cannot be same as input vector y, suggest VecAXPY()");
810: if (w == x) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Result vector w cannot be same as input vector x, suggest VecAYPX()");
813: PetscLogEventBegin(VEC_WAXPY,x,y,w,0);
814: (*w->ops->waxpy)(w,alpha,x,y);
815: PetscLogEventEnd(VEC_WAXPY,x,y,w,0);
816: PetscObjectStateIncrease((PetscObject)w);
817: return(0);
818: }
823: /*@
824: VecSetValues - Inserts or adds values into certain locations of a vector.
826: Not Collective
828: Input Parameters:
829: + x - vector to insert in
830: . ni - number of elements to add
831: . ix - indices where to add
832: . y - array of values
833: - iora - either INSERT_VALUES or ADD_VALUES, where
834: ADD_VALUES adds values to any existing entries, and
835: INSERT_VALUES replaces existing entries with new values
837: Notes:
838: VecSetValues() sets x[ix[i]] = y[i], for i=0,...,ni-1.
840: Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
841: options cannot be mixed without intervening calls to the assembly
842: routines.
844: These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd()
845: MUST be called after all calls to VecSetValues() have been completed.
847: VecSetValues() uses 0-based indices in Fortran as well as in C.
849: If you call VecSetOption(x, VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE),
850: negative indices may be passed in ix. These rows are
851: simply ignored. This allows easily inserting element load matrices
852: with homogeneous Dirchlet boundary conditions that you don't want represented
853: in the vector.
855: Level: beginner
857: Concepts: vector^setting values
859: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesLocal(),
860: VecSetValue(), VecSetValuesBlocked(), InsertMode, INSERT_VALUES, ADD_VALUES, VecGetValues()
861: @*/
862: PetscErrorCode VecSetValues(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
863: {
871: PetscLogEventBegin(VEC_SetValues,x,0,0,0);
872: (*x->ops->setvalues)(x,ni,ix,y,iora);
873: PetscLogEventEnd(VEC_SetValues,x,0,0,0);
874: PetscObjectStateIncrease((PetscObject)x);
875: return(0);
876: }
880: /*@
881: VecGetValues - Gets values from certain locations of a vector. Currently
882: can only get values on the same processor
884: Not Collective
886: Input Parameters:
887: + x - vector to get values from
888: . ni - number of elements to get
889: - ix - indices where to get them from (in global 1d numbering)
891: Output Parameter:
892: . y - array of values
894: Notes:
895: The user provides the allocated array y; it is NOT allocated in this routine
897: VecGetValues() gets y[i] = x[ix[i]], for i=0,...,ni-1.
899: VecAssemblyBegin() and VecAssemblyEnd() MUST be called before calling this
901: VecGetValues() uses 0-based indices in Fortran as well as in C.
903: If you call VecSetOption(x, VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE),
904: negative indices may be passed in ix. These rows are
905: simply ignored.
907: Level: beginner
909: Concepts: vector^getting values
911: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecGetValuesLocal(),
912: VecGetValuesBlocked(), InsertMode, INSERT_VALUES, ADD_VALUES, VecSetValues()
913: @*/
914: PetscErrorCode VecGetValues(Vec x,PetscInt ni,const PetscInt ix[],PetscScalar y[])
915: {
923: (*x->ops->getvalues)(x,ni,ix,y);
924: return(0);
925: }
929: /*@
930: VecSetValuesBlocked - Inserts or adds blocks of values into certain locations of a vector.
932: Not Collective
934: Input Parameters:
935: + x - vector to insert in
936: . ni - number of blocks to add
937: . ix - indices where to add in block count, rather than element count
938: . y - array of values
939: - iora - either INSERT_VALUES or ADD_VALUES, where
940: ADD_VALUES adds values to any existing entries, and
941: INSERT_VALUES replaces existing entries with new values
943: Notes:
944: VecSetValuesBlocked() sets x[bs*ix[i]+j] = y[bs*i+j],
945: for j=0,...,bs, for i=0,...,ni-1. where bs was set with VecSetBlockSize().
947: Calls to VecSetValuesBlocked() with the INSERT_VALUES and ADD_VALUES
948: options cannot be mixed without intervening calls to the assembly
949: routines.
951: These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd()
952: MUST be called after all calls to VecSetValuesBlocked() have been completed.
954: VecSetValuesBlocked() uses 0-based indices in Fortran as well as in C.
956: Negative indices may be passed in ix, these rows are
957: simply ignored. This allows easily inserting element load matrices
958: with homogeneous Dirchlet boundary conditions that you don't want represented
959: in the vector.
961: Level: intermediate
963: Concepts: vector^setting values blocked
965: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesBlockedLocal(),
966: VecSetValues()
967: @*/
968: PetscErrorCode VecSetValuesBlocked(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
969: {
977: PetscLogEventBegin(VEC_SetValues,x,0,0,0);
978: (*x->ops->setvaluesblocked)(x,ni,ix,y,iora);
979: PetscLogEventEnd(VEC_SetValues,x,0,0,0);
980: PetscObjectStateIncrease((PetscObject)x);
981: return(0);
982: }
987: /*@
988: VecSetValuesLocal - Inserts or adds values into certain locations of a vector,
989: using a local ordering of the nodes.
991: Not Collective
993: Input Parameters:
994: + x - vector to insert in
995: . ni - number of elements to add
996: . ix - indices where to add
997: . y - array of values
998: - iora - either INSERT_VALUES or ADD_VALUES, where
999: ADD_VALUES adds values to any existing entries, and
1000: INSERT_VALUES replaces existing entries with new values
1002: Level: intermediate
1004: Notes:
1005: VecSetValuesLocal() sets x[ix[i]] = y[i], for i=0,...,ni-1.
1007: Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
1008: options cannot be mixed without intervening calls to the assembly
1009: routines.
1011: These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd()
1012: MUST be called after all calls to VecSetValuesLocal() have been completed.
1014: VecSetValuesLocal() uses 0-based indices in Fortran as well as in C.
1016: Concepts: vector^setting values with local numbering
1018: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetLocalToGlobalMapping(),
1019: VecSetValuesBlockedLocal()
1020: @*/
1021: PetscErrorCode VecSetValuesLocal(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
1022: {
1024: PetscInt lixp[128],*lix = lixp;
1032: PetscLogEventBegin(VEC_SetValues,x,0,0,0);
1033: if (!x->ops->setvalueslocal) {
1034: if (!x->map->mapping) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Local to global never set with VecSetLocalToGlobalMapping()");
1035: if (ni > 128) {
1036: PetscMalloc(ni*sizeof(PetscInt),&lix);
1037: }
1038: ISLocalToGlobalMappingApply(x->map->mapping,ni,(PetscInt*)ix,lix);
1039: (*x->ops->setvalues)(x,ni,lix,y,iora);
1040: if (ni > 128) {
1041: PetscFree(lix);
1042: }
1043: } else {
1044: (*x->ops->setvalueslocal)(x,ni,ix,y,iora);
1045: }
1046: PetscLogEventEnd(VEC_SetValues,x,0,0,0);
1047: PetscObjectStateIncrease((PetscObject)x);
1048: return(0);
1049: }
1053: /*@
1054: VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
1055: using a local ordering of the nodes.
1057: Not Collective
1059: Input Parameters:
1060: + x - vector to insert in
1061: . ni - number of blocks to add
1062: . ix - indices where to add in block count, not element count
1063: . y - array of values
1064: - iora - either INSERT_VALUES or ADD_VALUES, where
1065: ADD_VALUES adds values to any existing entries, and
1066: INSERT_VALUES replaces existing entries with new values
1068: Level: intermediate
1070: Notes:
1071: VecSetValuesBlockedLocal() sets x[bs*ix[i]+j] = y[bs*i+j],
1072: for j=0,..bs-1, for i=0,...,ni-1, where bs has been set with VecSetBlockSize().
1074: Calls to VecSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES
1075: options cannot be mixed without intervening calls to the assembly
1076: routines.
1078: These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd()
1079: MUST be called after all calls to VecSetValuesBlockedLocal() have been completed.
1081: VecSetValuesBlockedLocal() uses 0-based indices in Fortran as well as in C.
1084: Concepts: vector^setting values blocked with local numbering
1086: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetValuesBlocked(),
1087: VecSetLocalToGlobalMappingBlock()
1088: @*/
1089: PetscErrorCode VecSetValuesBlockedLocal(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
1090: {
1092: PetscInt lixp[128],*lix = lixp;
1099: if (!x->map->bmapping) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Local to global never set with VecSetLocalToGlobalMappingBlock()");
1100: if (ni > 128) {
1101: PetscMalloc(ni*sizeof(PetscInt),&lix);
1102: }
1104: PetscLogEventBegin(VEC_SetValues,x,0,0,0);
1105: ISLocalToGlobalMappingApply(x->map->bmapping,ni,(PetscInt*)ix,lix);
1106: (*x->ops->setvaluesblocked)(x,ni,lix,y,iora);
1107: PetscLogEventEnd(VEC_SetValues,x,0,0,0);
1108: if (ni > 128) {
1109: PetscFree(lix);
1110: }
1111: PetscObjectStateIncrease((PetscObject)x);
1112: return(0);
1113: }
1117: /*@
1118: VecMTDot - Computes indefinite vector multiple dot products.
1119: That is, it does NOT use the complex conjugate.
1121: Collective on Vec
1123: Input Parameters:
1124: + x - one vector
1125: . nv - number of vectors
1126: - y - array of vectors. Note that vectors are pointers
1128: Output Parameter:
1129: . val - array of the dot products
1131: Notes for Users of Complex Numbers:
1132: For complex vectors, VecMTDot() computes the indefinite form
1133: $ val = (x,y) = y^T x,
1134: where y^T denotes the transpose of y.
1136: Use VecMDot() for the inner product
1137: $ val = (x,y) = y^H x,
1138: where y^H denotes the conjugate transpose of y.
1140: Level: intermediate
1142: Concepts: inner product^multiple
1143: Concepts: vector^multiple inner products
1145: .seealso: VecMDot(), VecTDot()
1146: @*/
1147: PetscErrorCode VecMTDot(Vec x,PetscInt nv,const Vec y[],PetscScalar val[])
1148: {
1161: PetscLogEventBegin(VEC_MTDot,x,*y,0,0);
1162: (*x->ops->mtdot)(x,nv,y,val);
1163: PetscLogEventEnd(VEC_MTDot,x,*y,0,0);
1164: return(0);
1165: }
1169: /*@
1170: VecMDot - Computes vector multiple dot products.
1172: Collective on Vec
1174: Input Parameters:
1175: + x - one vector
1176: . nv - number of vectors
1177: - y - array of vectors.
1179: Output Parameter:
1180: . val - array of the dot products (does not allocate the array)
1182: Notes for Users of Complex Numbers:
1183: For complex vectors, VecMDot() computes
1184: $ val = (x,y) = y^H x,
1185: where y^H denotes the conjugate transpose of y.
1187: Use VecMTDot() for the indefinite form
1188: $ val = (x,y) = y^T x,
1189: where y^T denotes the transpose of y.
1191: Level: intermediate
1193: Concepts: inner product^multiple
1194: Concepts: vector^multiple inner products
1196: .seealso: VecMTDot(), VecDot()
1197: @*/
1198: PetscErrorCode VecMDot(Vec x,PetscInt nv,const Vec y[],PetscScalar val[])
1199: {
1204: if (!nv) return(0);
1205: if (nv < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of vectors (given %D) cannot be negative",nv);
1214: PetscLogEventBarrierBegin(VEC_MDotBarrier,x,*y,0,0,PetscObjectComm((PetscObject)x));
1215: (*x->ops->mdot)(x,nv,y,val);
1216: PetscLogEventBarrierEnd(VEC_MDotBarrier,x,*y,0,0,PetscObjectComm((PetscObject)x));
1217: return(0);
1218: }
1222: /*@
1223: VecMAXPY - Computes y = y + sum alpha[j] x[j]
1225: Logically Collective on Vec
1227: Input Parameters:
1228: + nv - number of scalars and x-vectors
1229: . alpha - array of scalars
1230: . y - one vector
1231: - x - array of vectors
1233: Level: intermediate
1235: Notes: y cannot be any of the x vectors
1237: Concepts: BLAS
1239: .seealso: VecAXPY(), VecWAXPY(), VecAYPX()
1240: @*/
1241: PetscErrorCode VecMAXPY(Vec y,PetscInt nv,const PetscScalar alpha[],Vec x[])
1242: {
1244: PetscInt i;
1248: if (!nv) return(0);
1249: if (nv < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of vectors (given %D) cannot be negative",nv);
1259: PetscLogEventBegin(VEC_MAXPY,*x,y,0,0);
1260: (*y->ops->maxpy)(y,nv,alpha,x);
1261: PetscLogEventEnd(VEC_MAXPY,*x,y,0,0);
1262: PetscObjectStateIncrease((PetscObject)y);
1263: return(0);
1264: }
1268: /*@
1269: VecGetSubVector - Gets a vector representing part of another vector
1271: Collective on IS (and Vec if nonlocal entries are needed)
1273: Input Arguments:
1274: + X - vector from which to extract a subvector
1275: - is - index set representing portion of X to extract
1277: Output Arguments:
1278: . Y - subvector corresponding to is
1280: Level: advanced
1282: Notes:
1283: The subvector Y should be returned with VecRestoreSubVector().
1285: This function may return a subvector without making a copy, therefore it is not safe to use the original vector while
1286: modifying the subvector. Other non-overlapping subvectors can still be obtained from X using this function.
1288: .seealso: MatGetSubMatrix()
1289: @*/
1290: PetscErrorCode VecGetSubVector(Vec X,IS is,Vec *Y)
1291: {
1293: Vec Z;
1294: PetscInt state;
1300: if (X->ops->getsubvector) {
1301: (*X->ops->getsubvector)(X,is,&Z);
1302: } else { /* Default implementation currently does no caching */
1303: PetscInt gstart,gend,start;
1304: PetscBool contiguous,gcontiguous;
1305: VecGetOwnershipRange(X,&gstart,&gend);
1306: ISContiguousLocal(is,gstart,gend,&start,&contiguous);
1307: MPI_Allreduce(&contiguous,&gcontiguous,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)is));
1308: if (gcontiguous) { /* We can do a no-copy implementation */
1309: PetscInt n,N;
1310: PetscScalar *x;
1311: PetscMPIInt size;
1312: ISGetLocalSize(is,&n);
1313: VecGetArray(X,&x);
1314: MPI_Comm_size(PetscObjectComm((PetscObject)X),&size);
1315: if (size == 1) {
1316: VecCreateSeqWithArray(PetscObjectComm((PetscObject)X),1,n,x+start,&Z);
1317: } else {
1318: ISGetSize(is,&N);
1319: VecCreateMPIWithArray(PetscObjectComm((PetscObject)X),1,n,N,x+start,&Z);
1320: }
1321: VecRestoreArray(X,&x);
1322: } else { /* Have to create a scatter and do a copy */
1323: VecScatter scatter;
1324: PetscInt n,N;
1325: ISGetLocalSize(is,&n);
1326: ISGetSize(is,&N);
1327: VecCreate(PetscObjectComm((PetscObject)is),&Z);
1328: VecSetSizes(Z,n,N);
1329: VecSetType(Z,((PetscObject)X)->type_name);
1330: VecScatterCreate(X,is,Z,NULL,&scatter);
1331: VecScatterBegin(scatter,X,Z,INSERT_VALUES,SCATTER_FORWARD);
1332: VecScatterEnd(scatter,X,Z,INSERT_VALUES,SCATTER_FORWARD);
1333: VecScatterDestroy(&scatter);
1334: }
1335: }
1336: /* Record the state when the subvector was gotten so we know whether its values need to be put back */
1337: if (VecGetSubVectorSavedStateId < 0) {PetscObjectComposedDataRegister(&VecGetSubVectorSavedStateId);}
1338: PetscObjectStateQuery((PetscObject)Z,&state);
1339: PetscObjectComposedDataSetInt((PetscObject)Z,VecGetSubVectorSavedStateId,state);
1340: *Y = Z;
1341: return(0);
1342: }
1346: /*@
1347: VecRestoreSubVector - Restores a subvector extracted using VecGetSubVector()
1349: Collective on IS (and Vec if nonlocal entries need to be written)
1351: Input Arguments:
1352: + X - vector from which subvector was obtained
1353: . is - index set representing the subset of X
1354: - Y - subvector being restored
1356: Level: advanced
1358: .seealso: VecGetSubVector()
1359: @*/
1360: PetscErrorCode VecRestoreSubVector(Vec X,IS is,Vec *Y)
1361: {
1369: if (X->ops->restoresubvector) {
1370: (*X->ops->restoresubvector)(X,is,Y);
1371: } else {
1372: PetscInt savedstate=0,newstate;
1373: PetscBool valid;
1374: PetscObjectComposedDataGetInt((PetscObject)*Y,VecGetSubVectorSavedStateId,savedstate,valid);
1375: PetscObjectStateQuery((PetscObject)*Y,&newstate);
1376: if (valid && savedstate < newstate) {
1377: /* We might need to copy entries back, first check whether we have no-copy view */
1378: PetscInt gstart,gend,start;
1379: PetscBool contiguous,gcontiguous;
1380: VecGetOwnershipRange(X,&gstart,&gend);
1381: ISContiguousLocal(is,gstart,gend,&start,&contiguous);
1382: MPI_Allreduce(&contiguous,&gcontiguous,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)is));
1383: if (!gcontiguous) SETERRQ(PetscObjectComm((PetscObject)is),PETSC_ERR_SUP,"Unhandled case, values have been changed and need to be copied back into X");
1384: }
1385: VecDestroy(Y);
1386: }
1387: return(0);
1388: }
1392: /*@C
1393: VecGetArray - Returns a pointer to a contiguous array that contains this
1394: processor's portion of the vector data. For the standard PETSc
1395: vectors, VecGetArray() returns a pointer to the local data array and
1396: does not use any copies. If the underlying vector data is not stored
1397: in a contiquous array this routine will copy the data to a contiquous
1398: array and return a pointer to that. You MUST call VecRestoreArray()
1399: when you no longer need access to the array.
1401: Logically Collective on Vec
1403: Input Parameter:
1404: . x - the vector
1406: Output Parameter:
1407: . a - location to put pointer to the array
1409: Fortran Note:
1410: This routine is used differently from Fortran 77
1411: $ Vec x
1412: $ PetscScalar x_array(1)
1413: $ PetscOffset i_x
1414: $ PetscErrorCode ierr
1415: $ call VecGetArray(x,x_array,i_x,ierr)
1416: $
1417: $ Access first local entry in vector with
1418: $ value = x_array(i_x + 1)
1419: $
1420: $ ...... other code
1421: $ call VecRestoreArray(x,x_array,i_x,ierr)
1422: For Fortran 90 see VecGetArrayF90()
1424: See the Fortran chapter of the users manual and
1425: petsc/src/snes/examples/tutorials/ex5f.F for details.
1427: Level: beginner
1429: Concepts: vector^accessing local values
1431: .seealso: VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(), VecGetArray2d()
1432: @*/
1433: PetscErrorCode VecGetArray(Vec x,PetscScalar **a)
1434: {
1439: if (x->petscnative) {
1440: #if defined(PETSC_HAVE_CUSP)
1441: if (x->valid_GPU_array == PETSC_CUSP_GPU) {
1442: VecCUSPCopyFromGPU(x);
1443: }
1444: #endif
1445: *a = *((PetscScalar **)x->data);
1446: } else {
1447: (*x->ops->getarray)(x,a);
1448: }
1449: return(0);
1450: }
1454: /*@C
1455: VecGetArrayRead - Get read-only pointer to contiguous array containing this processor's portion of the vector data.
1457: Not Collective
1459: Input Parameters:
1460: . x - the vector
1462: Output Parameter:
1463: . a - the array
1465: Level: beginner
1467: Notes:
1468: The array must be returned using a matching call to VecRestoreArrayRead().
1470: Unlike VecGetArray(), this routine is not collective and preserves cached information like vector norms.
1472: Standard PETSc vectors use contiguous storage so that this routine does not perform a copy. Other vector
1473: implementations may require a copy, but must such implementations should cache the contiguous representation so that
1474: only one copy is performed when this routine is called multiple times in sequence.
1476: .seealso: VecGetArray(), VecRestoreArray()
1477: @*/
1478: PetscErrorCode VecGetArrayRead(Vec x,const PetscScalar **a)
1479: {
1484: if (x->petscnative) {
1485: #if defined(PETSC_HAVE_CUSP)
1486: if (x->valid_GPU_array == PETSC_CUSP_GPU) {
1487: VecCUSPCopyFromGPU(x);
1488: }
1489: #endif
1490: *a = *((PetscScalar **)x->data);
1491: } else if (x->ops->getarrayread) {
1492: (*x->ops->getarrayread)(x,a);
1493: } else {
1494: (*x->ops->getarray)(x,(PetscScalar**)a);
1495: }
1496: return(0);
1497: }
1501: /*@C
1502: VecGetArrays - Returns a pointer to the arrays in a set of vectors
1503: that were created by a call to VecDuplicateVecs(). You MUST call
1504: VecRestoreArrays() when you no longer need access to the array.
1506: Logically Collective on Vec
1508: Input Parameter:
1509: + x - the vectors
1510: - n - the number of vectors
1512: Output Parameter:
1513: . a - location to put pointer to the array
1515: Fortran Note:
1516: This routine is not supported in Fortran.
1518: Level: intermediate
1520: .seealso: VecGetArray(), VecRestoreArrays()
1521: @*/
1522: PetscErrorCode VecGetArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1523: {
1525: PetscInt i;
1526: PetscScalar **q;
1532: if (n <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must get at least one array n = %D",n);
1533: PetscMalloc(n*sizeof(PetscScalar*),&q);
1534: for (i=0; i<n; ++i) {
1535: VecGetArray(x[i],&q[i]);
1536: }
1537: *a = q;
1538: return(0);
1539: }
1543: /*@C
1544: VecRestoreArrays - Restores a group of vectors after VecGetArrays()
1545: has been called.
1547: Logically Collective on Vec
1549: Input Parameters:
1550: + x - the vector
1551: . n - the number of vectors
1552: - a - location of pointer to arrays obtained from VecGetArrays()
1554: Notes:
1555: For regular PETSc vectors this routine does not involve any copies. For
1556: any special vectors that do not store local vector data in a contiguous
1557: array, this routine will copy the data back into the underlying
1558: vector data structure from the arrays obtained with VecGetArrays().
1560: Fortran Note:
1561: This routine is not supported in Fortran.
1563: Level: intermediate
1565: .seealso: VecGetArrays(), VecRestoreArray()
1566: @*/
1567: PetscErrorCode VecRestoreArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1568: {
1570: PetscInt i;
1571: PetscScalar **q = *a;
1578: for (i=0; i<n; ++i) {
1579: VecRestoreArray(x[i],&q[i]);
1580: }
1581: PetscFree(q);
1582: return(0);
1583: }
1587: /*@C
1588: VecRestoreArray - Restores a vector after VecGetArray() has been called.
1590: Logically Collective on Vec
1592: Input Parameters:
1593: + x - the vector
1594: - a - location of pointer to array obtained from VecGetArray()
1596: Level: beginner
1598: Notes:
1599: For regular PETSc vectors this routine does not involve any copies. For
1600: any special vectors that do not store local vector data in a contiguous
1601: array, this routine will copy the data back into the underlying
1602: vector data structure from the array obtained with VecGetArray().
1604: This routine actually zeros out the a pointer. This is to prevent accidental
1605: us of the array after it has been restored. If you pass null for a it will
1606: not zero the array pointer a.
1608: Fortran Note:
1609: This routine is used differently from Fortran 77
1610: $ Vec x
1611: $ PetscScalar x_array(1)
1612: $ PetscOffset i_x
1613: $ PetscErrorCode ierr
1614: $ call VecGetArray(x,x_array,i_x,ierr)
1615: $
1616: $ Access first local entry in vector with
1617: $ value = x_array(i_x + 1)
1618: $
1619: $ ...... other code
1620: $ call VecRestoreArray(x,x_array,i_x,ierr)
1622: See the Fortran chapter of the users manual and
1623: petsc/src/snes/examples/tutorials/ex5f.F for details.
1624: For Fortran 90 see VecRestoreArrayF90()
1626: .seealso: VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(), VecRestoreArray2d()
1627: @*/
1628: PetscErrorCode VecRestoreArray(Vec x,PetscScalar **a)
1629: {
1634: if (x->petscnative) {
1635: #if defined(PETSC_HAVE_CUSP)
1636: x->valid_GPU_array = PETSC_CUSP_CPU;
1637: #endif
1638: } else {
1639: (*x->ops->restorearray)(x,a);
1640: }
1641: if (a) *a = NULL;
1642: PetscObjectStateIncrease((PetscObject)x);
1643: return(0);
1644: }
1648: /*@C
1649: VecRestoreArrayRead - Restore array obtained with VecGetArrayRead()
1651: Not Collective
1653: Input Parameters:
1654: + vec - the vector
1655: - array - the array
1657: Level: beginner
1659: .seealso: VecGetArray(), VecRestoreArray()
1660: @*/
1661: PetscErrorCode VecRestoreArrayRead(Vec x,const PetscScalar **a)
1662: {
1667: if (x->petscnative) {
1668: #if defined(PETSC_HAVE_CUSP)
1669: x->valid_GPU_array = PETSC_CUSP_CPU;
1670: #endif
1671: } else if (x->ops->restorearrayread) {
1672: (*x->ops->restorearrayread)(x,a);
1673: } else {
1674: (*x->ops->restorearray)(x,(PetscScalar**)a);
1675: }
1676: if (a) *a = NULL;
1677: return(0);
1678: }
1682: /*@
1683: VecPlaceArray - Allows one to replace the array in a vector with an
1684: array provided by the user. This is useful to avoid copying an array
1685: into a vector.
1687: Not Collective
1689: Input Parameters:
1690: + vec - the vector
1691: - array - the array
1693: Notes:
1694: You can return to the original array with a call to VecResetArray()
1696: Level: developer
1698: .seealso: VecGetArray(), VecRestoreArray(), VecReplaceArray(), VecResetArray()
1700: @*/
1701: PetscErrorCode VecPlaceArray(Vec vec,const PetscScalar array[])
1702: {
1709: if (vec->ops->placearray) {
1710: (*vec->ops->placearray)(vec,array);
1711: } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot place array in this type of vector");
1712: PetscObjectStateIncrease((PetscObject)vec);
1713: return(0);
1714: }
1718: /*@C
1719: VecReplaceArray - Allows one to replace the array in a vector with an
1720: array provided by the user. This is useful to avoid copying an array
1721: into a vector.
1723: Not Collective
1725: Input Parameters:
1726: + vec - the vector
1727: - array - the array
1729: Notes:
1730: This permanently replaces the array and frees the memory associated
1731: with the old array.
1733: The memory passed in MUST be obtained with PetscMalloc() and CANNOT be
1734: freed by the user. It will be freed when the vector is destroy.
1736: Not supported from Fortran
1738: Level: developer
1740: .seealso: VecGetArray(), VecRestoreArray(), VecPlaceArray(), VecResetArray()
1742: @*/
1743: PetscErrorCode VecReplaceArray(Vec vec,const PetscScalar array[])
1744: {
1750: if (vec->ops->replacearray) {
1751: (*vec->ops->replacearray)(vec,array);
1752: } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot replace array in this type of vector");
1753: PetscObjectStateIncrease((PetscObject)vec);
1754: return(0);
1755: }
1757: /*MC
1758: VecDuplicateVecsF90 - Creates several vectors of the same type as an existing vector
1759: and makes them accessible via a Fortran90 pointer.
1761: Synopsis:
1762: VecDuplicateVecsF90(Vec x,PetscInt n,{Vec, pointer :: y(:)},integer ierr)
1764: Collective on Vec
1766: Input Parameters:
1767: + x - a vector to mimic
1768: - n - the number of vectors to obtain
1770: Output Parameters:
1771: + y - Fortran90 pointer to the array of vectors
1772: - ierr - error code
1774: Example of Usage:
1775: .vb
1776: Vec x
1777: Vec, pointer :: y(:)
1778: ....
1779: call VecDuplicateVecsF90(x,2,y,ierr)
1780: call VecSet(y(2),alpha,ierr)
1781: call VecSet(y(2),alpha,ierr)
1782: ....
1783: call VecDestroyVecsF90(2,y,ierr)
1784: .ve
1786: Notes:
1787: Not yet supported for all F90 compilers
1789: Use VecDestroyVecsF90() to free the space.
1791: Level: beginner
1793: .seealso: VecDestroyVecsF90(), VecDuplicateVecs()
1795: M*/
1797: /*MC
1798: VecRestoreArrayF90 - Restores a vector to a usable state after a call to
1799: VecGetArrayF90().
1801: Synopsis:
1802: VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
1804: Logically Collective on Vec
1806: Input Parameters:
1807: + x - vector
1808: - xx_v - the Fortran90 pointer to the array
1810: Output Parameter:
1811: . ierr - error code
1813: Example of Usage:
1814: .vb
1815: PetscScalar, pointer :: xx_v(:)
1816: ....
1817: call VecGetArrayF90(x,xx_v,ierr)
1818: a = xx_v(3)
1819: call VecRestoreArrayF90(x,xx_v,ierr)
1820: .ve
1822: Level: beginner
1824: .seealso: VecGetArrayF90(), VecGetArray(), VecRestoreArray(), UsingFortran
1826: M*/
1828: /*MC
1829: VecDestroyVecsF90 - Frees a block of vectors obtained with VecDuplicateVecsF90().
1831: Synopsis:
1832: VecDestroyVecsF90(PetscInt n,{Vec, pointer :: x(:)},PetscErrorCode ierr)
1834: Collective on Vec
1836: Input Parameters:
1837: + n - the number of vectors previously obtained
1838: - x - pointer to array of vector pointers
1840: Output Parameter:
1841: . ierr - error code
1843: Notes:
1844: Not yet supported for all F90 compilers
1846: Level: beginner
1848: .seealso: VecDestroyVecs(), VecDuplicateVecsF90()
1850: M*/
1852: /*MC
1853: VecGetArrayF90 - Accesses a vector array from Fortran90. For default PETSc
1854: vectors, VecGetArrayF90() returns a pointer to the local data array. Otherwise,
1855: this routine is implementation dependent. You MUST call VecRestoreArrayF90()
1856: when you no longer need access to the array.
1858: Synopsis:
1859: VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
1861: Logically Collective on Vec
1863: Input Parameter:
1864: . x - vector
1866: Output Parameters:
1867: + xx_v - the Fortran90 pointer to the array
1868: - ierr - error code
1870: Example of Usage:
1871: .vb
1872: PetscScalar, pointer :: xx_v(:)
1873: ....
1874: call VecGetArrayF90(x,xx_v,ierr)
1875: a = xx_v(3)
1876: call VecRestoreArrayF90(x,xx_v,ierr)
1877: .ve
1879: Level: beginner
1881: .seealso: VecRestoreArrayF90(), VecGetArray(), VecRestoreArray(), UsingFortran
1883: M*/
1888: /*@C
1889: VecGetArray2d - Returns a pointer to a 2d contiguous array that contains this
1890: processor's portion of the vector data. You MUST call VecRestoreArray2d()
1891: when you no longer need access to the array.
1893: Logically Collective
1895: Input Parameter:
1896: + x - the vector
1897: . m - first dimension of two dimensional array
1898: . n - second dimension of two dimensional array
1899: . mstart - first index you will use in first coordinate direction (often 0)
1900: - nstart - first index in the second coordinate direction (often 0)
1902: Output Parameter:
1903: . a - location to put pointer to the array
1905: Level: developer
1907: Notes:
1908: For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
1909: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
1910: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
1911: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().
1913: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
1915: Concepts: vector^accessing local values as 2d array
1917: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
1918: VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
1919: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
1920: @*/
1921: PetscErrorCode VecGetArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
1922: {
1924: PetscInt i,N;
1925: PetscScalar *aa;
1931: VecGetLocalSize(x,&N);
1932: if (m*n != N) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local array size %D does not match 2d array dimensions %D by %D",N,m,n);
1933: VecGetArray(x,&aa);
1935: PetscMalloc(m*sizeof(PetscScalar*),a);
1936: for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
1937: *a -= mstart;
1938: return(0);
1939: }
1943: /*@C
1944: VecRestoreArray2d - Restores a vector after VecGetArray2d() has been called.
1946: Logically Collective
1948: Input Parameters:
1949: + x - the vector
1950: . m - first dimension of two dimensional array
1951: . n - second dimension of the two dimensional array
1952: . mstart - first index you will use in first coordinate direction (often 0)
1953: . nstart - first index in the second coordinate direction (often 0)
1954: - a - location of pointer to array obtained from VecGetArray2d()
1956: Level: developer
1958: Notes:
1959: For regular PETSc vectors this routine does not involve any copies. For
1960: any special vectors that do not store local vector data in a contiguous
1961: array, this routine will copy the data back into the underlying
1962: vector data structure from the array obtained with VecGetArray().
1964: This routine actually zeros out the a pointer.
1966: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
1967: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
1968: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
1969: @*/
1970: PetscErrorCode VecRestoreArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
1971: {
1973: void *dummy;
1979: dummy = (void*)(*a + mstart);
1980: PetscFree(dummy);
1981: VecRestoreArray(x,NULL);
1982: return(0);
1983: }
1987: /*@C
1988: VecGetArray1d - Returns a pointer to a 1d contiguous array that contains this
1989: processor's portion of the vector data. You MUST call VecRestoreArray1d()
1990: when you no longer need access to the array.
1992: Logically Collective
1994: Input Parameter:
1995: + x - the vector
1996: . m - first dimension of two dimensional array
1997: - mstart - first index you will use in first coordinate direction (often 0)
1999: Output Parameter:
2000: . a - location to put pointer to the array
2002: Level: developer
2004: Notes:
2005: For a vector obtained from DMCreateLocalVector() mstart are likely
2006: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2007: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().
2009: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2011: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2012: VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2013: VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2014: @*/
2015: PetscErrorCode VecGetArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2016: {
2018: PetscInt N;
2024: VecGetLocalSize(x,&N);
2025: if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
2026: VecGetArray(x,a);
2027: *a -= mstart;
2028: return(0);
2029: }
2033: /*@C
2034: VecRestoreArray1d - Restores a vector after VecGetArray1d() has been called.
2036: Logically Collective
2038: Input Parameters:
2039: + x - the vector
2040: . m - first dimension of two dimensional array
2041: . mstart - first index you will use in first coordinate direction (often 0)
2042: - a - location of pointer to array obtained from VecGetArray21()
2044: Level: developer
2046: Notes:
2047: For regular PETSc vectors this routine does not involve any copies. For
2048: any special vectors that do not store local vector data in a contiguous
2049: array, this routine will copy the data back into the underlying
2050: vector data structure from the array obtained with VecGetArray1d().
2052: This routine actually zeros out the a pointer.
2054: Concepts: vector^accessing local values as 1d array
2056: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2057: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2058: VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
2059: @*/
2060: PetscErrorCode VecRestoreArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2061: {
2067: VecRestoreArray(x,NULL);
2068: return(0);
2069: }
2074: /*@C
2075: VecGetArray3d - Returns a pointer to a 3d contiguous array that contains this
2076: processor's portion of the vector data. You MUST call VecRestoreArray3d()
2077: when you no longer need access to the array.
2079: Logically Collective
2081: Input Parameter:
2082: + x - the vector
2083: . m - first dimension of three dimensional array
2084: . n - second dimension of three dimensional array
2085: . p - third dimension of three dimensional array
2086: . mstart - first index you will use in first coordinate direction (often 0)
2087: . nstart - first index in the second coordinate direction (often 0)
2088: - pstart - first index in the third coordinate direction (often 0)
2090: Output Parameter:
2091: . a - location to put pointer to the array
2093: Level: developer
2095: Notes:
2096: For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
2097: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2098: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2099: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3d().
2101: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2103: Concepts: vector^accessing local values as 3d array
2105: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2106: VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2107: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2108: @*/
2109: PetscErrorCode VecGetArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
2110: {
2112: PetscInt i,N,j;
2113: PetscScalar *aa,**b;
2119: VecGetLocalSize(x,&N);
2120: if (m*n*p != N) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local array size %D does not match 3d array dimensions %D by %D by %D",N,m,n,p);
2121: VecGetArray(x,&aa);
2123: PetscMalloc(m*sizeof(PetscScalar**)+m*n*sizeof(PetscScalar*),a);
2124: b = (PetscScalar**)((*a) + m);
2125: for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
2126: for (i=0; i<m; i++)
2127: for (j=0; j<n; j++)
2128: b[i*n+j] = aa + i*n*p + j*p - pstart;
2130: *a -= mstart;
2131: return(0);
2132: }
2136: /*@C
2137: VecRestoreArray3d - Restores a vector after VecGetArray3d() has been called.
2139: Logically Collective
2141: Input Parameters:
2142: + x - the vector
2143: . m - first dimension of three dimensional array
2144: . n - second dimension of the three dimensional array
2145: . p - third dimension of the three dimensional array
2146: . mstart - first index you will use in first coordinate direction (often 0)
2147: . nstart - first index in the second coordinate direction (often 0)
2148: . pstart - first index in the third coordinate direction (often 0)
2149: - a - location of pointer to array obtained from VecGetArray3d()
2151: Level: developer
2153: Notes:
2154: For regular PETSc vectors this routine does not involve any copies. For
2155: any special vectors that do not store local vector data in a contiguous
2156: array, this routine will copy the data back into the underlying
2157: vector data structure from the array obtained with VecGetArray().
2159: This routine actually zeros out the a pointer.
2161: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2162: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2163: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
2164: @*/
2165: PetscErrorCode VecRestoreArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
2166: {
2168: void *dummy;
2174: dummy = (void*)(*a + mstart);
2175: PetscFree(dummy);
2176: VecRestoreArray(x,NULL);
2177: return(0);
2178: }
2182: /*@C
2183: VecGetArray4d - Returns a pointer to a 4d contiguous array that contains this
2184: processor's portion of the vector data. You MUST call VecRestoreArray4d()
2185: when you no longer need access to the array.
2187: Logically Collective
2189: Input Parameter:
2190: + x - the vector
2191: . m - first dimension of four dimensional array
2192: . n - second dimension of four dimensional array
2193: . p - third dimension of four dimensional array
2194: . q - fourth dimension of four dimensional array
2195: . mstart - first index you will use in first coordinate direction (often 0)
2196: . nstart - first index in the second coordinate direction (often 0)
2197: . pstart - first index in the third coordinate direction (often 0)
2198: - qstart - first index in the fourth coordinate direction (often 0)
2200: Output Parameter:
2201: . a - location to put pointer to the array
2203: Level: beginner
2205: Notes:
2206: For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
2207: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2208: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2209: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3d().
2211: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2213: Concepts: vector^accessing local values as 3d array
2215: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2216: VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2217: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2218: @*/
2219: PetscErrorCode VecGetArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
2220: {
2222: PetscInt i,N,j,k;
2223: PetscScalar *aa,***b,**c;
2229: VecGetLocalSize(x,&N);
2230: if (m*n*p*q != N) SETERRQ5(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local array size %D does not match 4d array dimensions %D by %D by %D by %D",N,m,n,p,q);
2231: VecGetArray(x,&aa);
2233: PetscMalloc(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p*sizeof(PetscScalar*),a);
2234: b = (PetscScalar***)((*a) + m);
2235: c = (PetscScalar**)(b + m*n);
2236: for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
2237: for (i=0; i<m; i++)
2238: for (j=0; j<n; j++)
2239: b[i*n+j] = c + i*n*p + j*p - pstart;
2240: for (i=0; i<m; i++)
2241: for (j=0; j<n; j++)
2242: for (k=0; k<p; k++)
2243: c[i*n*p+j*p+k] = aa + i*n*p*q + j*p*q + k*q - qstart;
2244: *a -= mstart;
2245: return(0);
2246: }
2250: /*@C
2251: VecRestoreArray4d - Restores a vector after VecGetArray3d() has been called.
2253: Logically Collective
2255: Input Parameters:
2256: + x - the vector
2257: . m - first dimension of four dimensional array
2258: . n - second dimension of the four dimensional array
2259: . p - third dimension of the four dimensional array
2260: . q - fourth dimension of the four dimensional array
2261: . mstart - first index you will use in first coordinate direction (often 0)
2262: . nstart - first index in the second coordinate direction (often 0)
2263: . pstart - first index in the third coordinate direction (often 0)
2264: . qstart - first index in the fourth coordinate direction (often 0)
2265: - a - location of pointer to array obtained from VecGetArray4d()
2267: Level: beginner
2269: Notes:
2270: For regular PETSc vectors this routine does not involve any copies. For
2271: any special vectors that do not store local vector data in a contiguous
2272: array, this routine will copy the data back into the underlying
2273: vector data structure from the array obtained with VecGetArray().
2275: This routine actually zeros out the a pointer.
2277: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2278: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2279: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
2280: @*/
2281: PetscErrorCode VecRestoreArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
2282: {
2284: void *dummy;
2290: dummy = (void*)(*a + mstart);
2291: PetscFree(dummy);
2292: VecRestoreArray(x,NULL);
2293: return(0);
2294: }