Actual source code: rvector.c
petsc-3.3-p7 2013-05-11
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) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths"); \
11: if ((x)->map->n != (y)->map->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");
16: /*@
17: VecMaxPointwiseDivide - Computes the maximum of the componentwise division max = max_i abs(x_i/y_i).
19: Logically Collective on Vec
21: Input Parameters:
22: . x, y - the vectors
24: Output Parameter:
25: . max - the result
27: Level: advanced
29: Notes: x and y may be the same vector
30: if a particular y_i is zero, it is treated as 1 in the above formula
32: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs()
33: @*/
34: PetscErrorCode VecMaxPointwiseDivide(Vec x,Vec y,PetscReal *max)
35: {
47: (*x->ops->maxpointwisedivide)(x,y,max);
48: return(0);
49: }
53: /*@
54: VecDot - Computes the vector dot product.
56: Collective on Vec
58: Input Parameters:
59: . x, y - the vectors
61: Output Parameter:
62: . val - the dot product
64: Performance Issues:
65: $ per-processor memory bandwidth
66: $ interprocessor latency
67: $ work load inbalance that causes certain processes to arrive much earlier than others
69: Notes for Users of Complex Numbers:
70: For complex vectors, VecDot() computes
71: $ val = (x,y) = y^H x,
72: where y^H denotes the conjugate transpose of y. Note that this corresponds to the usual "mathematicians" complex
73: inner product where the SECOND argument gets the complex conjugate. Since the BLASdot() complex conjugates the first
74: first argument we call the BLASdot() with the arguments reversed.
76: Use VecTDot() for the indefinite form
77: $ val = (x,y) = y^T x,
78: where y^T denotes the transpose of y.
80: Level: intermediate
82: Concepts: inner product
83: Concepts: vector^inner product
85: .seealso: VecMDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd()
86: @*/
87: PetscErrorCode VecDot(Vec x,Vec y,PetscScalar *val)
88: {
100: PetscLogEventBarrierBegin(VEC_DotBarrier,x,y,0,0,((PetscObject)x)->comm);
101: (*x->ops->dot)(x,y,val);
102: PetscLogEventBarrierEnd(VEC_DotBarrier,x,y,0,0,((PetscObject)x)->comm);
103: return(0);
104: }
108: /*@
109: VecNorm - Computes the vector norm.
111: Collective on Vec
113: Input Parameters:
114: + x - the vector
115: - type - one of NORM_1, NORM_2, NORM_INFINITY. Also available
116: NORM_1_AND_2, which computes both norms and stores them
117: in a two element array.
119: Output Parameter:
120: . val - the norm
122: Notes:
123: $ NORM_1 denotes sum_i |x_i|
124: $ NORM_2 denotes sqrt(sum_i (x_i)^2)
125: $ NORM_INFINITY denotes max_i |x_i|
127: Level: intermediate
129: Performance Issues:
130: $ per-processor memory bandwidth
131: $ interprocessor latency
132: $ work load inbalance that causes certain processes to arrive much earlier than others
134: Compile Option:
135: PETSC_HAVE_SLOW_BLAS_NORM2 will cause a C (loop unrolled) version of the norm to be used, rather
136: than the BLAS. This should probably only be used when one is using the FORTRAN BLAS routines
137: (as opposed to vendor provided) because the FORTRAN BLAS NRM2() routine is very slow.
139: Concepts: norm
140: Concepts: vector^norm
142: .seealso: VecDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecNormAvailable(),
143: VecNormBegin(), VecNormEnd()
145: @*/
146: PetscErrorCode VecNorm(Vec x,NormType type,PetscReal *val)
147: {
148: PetscBool flg;
155: if (((PetscObject)x)->precision != sizeof(PetscScalar)) SETERRQ(((PetscObject)x)->comm,PETSC_ERR_SUP,"Wrong precision of input argument");
157: /*
158: * Cached data?
159: */
160: if (type!=NORM_1_AND_2) {
161: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[type],*val,flg);
162: if (flg) return(0);
163: }
165: PetscLogEventBarrierBegin(VEC_NormBarrier,x,0,0,0,((PetscObject)x)->comm);
166: (*x->ops->norm)(x,type,val);
167: PetscLogEventBarrierEnd(VEC_NormBarrier,x,0,0,0,((PetscObject)x)->comm);
169: if (type!=NORM_1_AND_2) {
170: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[type],*val);
171: }
172: return(0);
173: }
177: /*@
178: VecNormAvailable - Returns the vector norm if it is already known.
180: Not Collective
182: Input Parameters:
183: + x - the vector
184: - type - one of NORM_1, NORM_2, NORM_INFINITY. Also available
185: NORM_1_AND_2, which computes both norms and stores them
186: in a two element array.
188: Output Parameter:
189: + available - PETSC_TRUE if the val returned is valid
190: - val - the norm
192: Notes:
193: $ NORM_1 denotes sum_i |x_i|
194: $ NORM_2 denotes sqrt(sum_i (x_i)^2)
195: $ NORM_INFINITY denotes max_i |x_i|
197: Level: intermediate
199: Performance Issues:
200: $ per-processor memory bandwidth
201: $ interprocessor latency
202: $ work load inbalance that causes certain processes to arrive much earlier than others
204: Compile Option:
205: PETSC_HAVE_SLOW_BLAS_NORM2 will cause a C (loop unrolled) version of the norm to be used, rather
206: than the BLAS. This should probably only be used when one is using the FORTRAN BLAS routines
207: (as opposed to vendor provided) because the FORTRAN BLAS NRM2() routine is very slow.
209: Concepts: norm
210: Concepts: vector^norm
212: .seealso: VecDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecNorm()
213: VecNormBegin(), VecNormEnd()
215: @*/
216: PetscErrorCode VecNormAvailable(Vec x,NormType type,PetscBool *available,PetscReal *val)
217: {
225: *available = PETSC_FALSE;
226: if (type!=NORM_1_AND_2) {
227: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[type],*val,*available);
228: }
229: return(0);
230: }
234: /*@
235: VecNormalize - Normalizes a vector by 2-norm.
237: Collective on Vec
239: Input Parameters:
240: + x - the vector
242: Output Parameter:
243: . x - the normalized vector
244: - val - the vector norm before normalization
246: Level: intermediate
248: Concepts: vector^normalizing
249: Concepts: normalizing^vector
251: @*/
252: PetscErrorCode VecNormalize(Vec x,PetscReal *val)
253: {
255: PetscReal norm;
260: PetscLogEventBegin(VEC_Normalize,x,0,0,0);
261: VecNorm(x,NORM_2,&norm);
262: if (norm == 0.0) {
263: PetscInfo(x,"Vector of zero norm can not be normalized; Returning only the zero norm\n");
264: } else if (norm != 1.0) {
265: PetscScalar tmp = 1.0/norm;
266: VecScale(x,tmp);
267: }
268: if (val) *val = norm;
269: PetscLogEventEnd(VEC_Normalize,x,0,0,0);
270: return(0);
271: }
275: /*@C
276: VecMax - Determines the maximum vector component and its location.
278: Collective on Vec
280: Input Parameter:
281: . x - the vector
283: Output Parameters:
284: + val - the maximum component
285: - p - the location of val (pass PETSC_NULL if you don't want this)
287: Notes:
288: Returns the value PETSC_MIN_REAL and p = -1 if the vector is of length 0.
290: Returns the smallest index with the maximum value
291: Level: intermediate
293: Concepts: maximum^of vector
294: Concepts: vector^maximum value
296: .seealso: VecNorm(), VecMin()
297: @*/
298: PetscErrorCode VecMax(Vec x,PetscInt *p,PetscReal *val)
299: {
306: PetscLogEventBegin(VEC_Max,x,0,0,0);
307: (*x->ops->max)(x,p,val);
308: PetscLogEventEnd(VEC_Max,x,0,0,0);
309: return(0);
310: }
314: /*@
315: VecMin - Determines the minimum vector component and its location.
317: Collective on Vec
319: Input Parameters:
320: . x - the vector
322: Output Parameter:
323: + val - the minimum component
324: - p - the location of val (pass PETSC_NULL if you don't want this location)
326: Level: intermediate
328: Notes:
329: Returns the value PETSC_MAX_REAL and p = -1 if the vector is of length 0.
331: This returns the smallest index with the minumum value
333: Concepts: minimum^of vector
334: Concepts: vector^minimum entry
336: .seealso: VecMax()
337: @*/
338: PetscErrorCode VecMin(Vec x,PetscInt *p,PetscReal *val)
339: {
346: PetscLogEventBegin(VEC_Min,x,0,0,0);
347: (*x->ops->min)(x,p,val);
348: PetscLogEventEnd(VEC_Min,x,0,0,0);
349: return(0);
350: }
354: /*@
355: VecTDot - Computes an indefinite vector dot product. That is, this
356: routine does NOT use the complex conjugate.
358: Collective on Vec
360: Input Parameters:
361: . x, y - the vectors
363: Output Parameter:
364: . val - the dot product
366: Notes for Users of Complex Numbers:
367: For complex vectors, VecTDot() computes the indefinite form
368: $ val = (x,y) = y^T x,
369: where y^T denotes the transpose of y.
371: Use VecDot() for the inner product
372: $ val = (x,y) = y^H x,
373: where y^H denotes the conjugate transpose of y.
375: Level: intermediate
377: Concepts: inner product^non-Hermitian
378: Concepts: vector^inner product
379: Concepts: non-Hermitian inner product
381: .seealso: VecDot(), VecMTDot()
382: @*/
383: PetscErrorCode VecTDot(Vec x,Vec y,PetscScalar *val)
384: {
396: PetscLogEventBegin(VEC_TDot,x,y,0,0);
397: (*x->ops->tdot)(x,y,val);
398: PetscLogEventEnd(VEC_TDot,x,y,0,0);
399: return(0);
400: }
404: /*@
405: VecScale - Scales a vector.
407: Not collective on Vec
409: Input Parameters:
410: + x - the vector
411: - alpha - the scalar
413: Output Parameter:
414: . x - the scaled vector
416: Note:
417: For a vector with n components, VecScale() computes
418: $ x[i] = alpha * x[i], for i=1,...,n.
420: Level: intermediate
422: Concepts: vector^scaling
423: Concepts: scaling^vector
425: @*/
426: PetscErrorCode VecScale (Vec x, PetscScalar alpha)
427: {
428: PetscReal norms[4] = {0.0,0.0,0.0, 0.0};
429: PetscBool flgs[4];
431: PetscInt i;
436: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
437: PetscLogEventBegin(VEC_Scale,x,0,0,0);
438: if (alpha != (PetscScalar)1.0) {
439: /* get current stashed norms */
440: for (i=0; i<4; i++) {
441: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],norms[i],flgs[i]);
442: }
443: (*x->ops->scale)(x,alpha);
444: PetscObjectStateIncrease((PetscObject)x);
445: /* put the scaled stashed norms back into the Vec */
446: for (i=0; i<4; i++) {
447: if (flgs[i]) {
448: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[i],PetscAbsScalar(alpha)*norms[i]);
449: }
450: }
451: }
452: PetscLogEventEnd(VEC_Scale,x,0,0,0);
453: return(0);
454: }
458: /*@
459: VecSet - Sets all components of a vector to a single scalar value.
461: Logically Collective on Vec
463: Input Parameters:
464: + x - the vector
465: - alpha - the scalar
467: Output Parameter:
468: . x - the vector
470: Note:
471: For a vector of dimension n, VecSet() computes
472: $ x[i] = alpha, for i=1,...,n,
473: so that all vector entries then equal the identical
474: scalar value, alpha. Use the more general routine
475: VecSetValues() to set different vector entries.
477: You CANNOT call this after you have called VecSetValues() but before you call
478: VecAssemblyBegin/End().
480: Level: beginner
482: .seealso VecSetValues(), VecSetValuesBlocked(), VecSetRandom()
484: Concepts: vector^setting to constant
486: @*/
487: PetscErrorCode VecSet(Vec x,PetscScalar alpha)
488: {
489: PetscReal val;
495: 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()");
498: PetscLogEventBegin(VEC_Set,x,0,0,0);
499: (*x->ops->set)(x,alpha);
500: PetscLogEventEnd(VEC_Set,x,0,0,0);
501: PetscObjectStateIncrease((PetscObject)x);
503: /* norms can be simply set */
504: val = PetscAbsScalar(alpha);
505: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_1],x->map->N * val);
506: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_INFINITY],val);
507: val = PetscSqrtReal((double)x->map->N) * val;
508: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_2],val);
509: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_FROBENIUS],val);
510: return(0);
511: }
516: /*@
517: VecAXPY - Computes y = alpha x + y.
519: Logically Collective on Vec
521: Input Parameters:
522: + alpha - the scalar
523: - x, y - the vectors
525: Output Parameter:
526: . y - output vector
528: Level: intermediate
530: Notes: x and y MUST be different vectors
532: Concepts: vector^BLAS
533: Concepts: BLAS
535: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY()
536: @*/
537: PetscErrorCode VecAXPY(Vec y,PetscScalar alpha,Vec x)
538: {
548: if (x == y) SETERRQ(((PetscObject)x)->comm,PETSC_ERR_ARG_IDN,"x and y cannot be the same vector");
551: PetscLogEventBegin(VEC_AXPY,x,y,0,0);
552: (*y->ops->axpy)(y,alpha,x);
553: PetscLogEventEnd(VEC_AXPY,x,y,0,0);
554: PetscObjectStateIncrease((PetscObject)y);
555: return(0);
556: }
560: /*@
561: VecAXPBY - Computes y = alpha x + beta y.
563: Logically Collective on Vec
565: Input Parameters:
566: + alpha,beta - the scalars
567: - x, y - the vectors
569: Output Parameter:
570: . y - output vector
572: Level: intermediate
574: Notes: x and y MUST be different vectors
576: Concepts: BLAS
577: Concepts: vector^BLAS
579: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPY()
580: @*/
581: PetscErrorCode VecAXPBY(Vec y,PetscScalar alpha,PetscScalar beta,Vec x)
582: {
592: if (x == y) SETERRQ(((PetscObject)x)->comm,PETSC_ERR_ARG_IDN,"x and y cannot be the same vector");
596: PetscLogEventBegin(VEC_AXPY,x,y,0,0);
597: (*y->ops->axpby)(y,alpha,beta,x);
598: PetscLogEventEnd(VEC_AXPY,x,y,0,0);
599: PetscObjectStateIncrease((PetscObject)y);
600: return(0);
601: }
605: /*@
606: VecAXPBYPCZ - Computes z = alpha x + beta y + gamma z
608: Logically Collective on Vec
610: Input Parameters:
611: + alpha,beta, gamma - the scalars
612: - x, y, z - the vectors
614: Output Parameter:
615: . z - output vector
617: Level: intermediate
619: Notes: x, y and z must be different vectors
621: Developer Note: alpha = 1 or gamma = 1 or gamma = 0.0 are handled as special cases
623: Concepts: BLAS
624: Concepts: vector^BLAS
626: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPY()
627: @*/
628: PetscErrorCode VecAXPBYPCZ(Vec z,PetscScalar alpha,PetscScalar beta,PetscScalar gamma,Vec x,Vec y)
629: {
643: if (x == y || x == z) SETERRQ(((PetscObject)x)->comm,PETSC_ERR_ARG_IDN,"x, y, and z must be different vectors");
644: if (y == z) SETERRQ(((PetscObject)y)->comm,PETSC_ERR_ARG_IDN,"x, y, and z must be different vectors");
649: PetscLogEventBegin(VEC_AXPBYPCZ,x,y,z,0);
650: (*y->ops->axpbypcz)(z,alpha,beta,gamma,x,y);
651: PetscLogEventEnd(VEC_AXPBYPCZ,x,y,z,0);
652: PetscObjectStateIncrease((PetscObject)z);
653: return(0);
654: }
658: /*@
659: VecAYPX - Computes y = x + alpha y.
661: Logically Collective on Vec
663: Input Parameters:
664: + alpha - the scalar
665: - x, y - the vectors
667: Output Parameter:
668: . y - output vector
670: Level: intermediate
672: Notes: x and y MUST be different vectors
674: Concepts: vector^BLAS
675: Concepts: BLAS
677: .seealso: VecAXPY(), VecWAXPY()
678: @*/
679: PetscErrorCode VecAYPX(Vec y,PetscScalar alpha,Vec x)
680: {
688: if (x == y) SETERRQ(((PetscObject)x)->comm,PETSC_ERR_ARG_IDN,"x and y must be different vectors");
691: PetscLogEventBegin(VEC_AYPX,x,y,0,0);
692: (*y->ops->aypx)(y,alpha,x);
693: PetscLogEventEnd(VEC_AYPX,x,y,0,0);
694: PetscObjectStateIncrease((PetscObject)y);
695: return(0);
696: }
701: /*@
702: VecWAXPY - Computes w = alpha x + y.
704: Logically Collective on Vec
706: Input Parameters:
707: + alpha - the scalar
708: - x, y - the vectors
710: Output Parameter:
711: . w - the result
713: Level: intermediate
715: Notes: w cannot be either x or y, but x and y can be the same
717: Concepts: vector^BLAS
718: Concepts: BLAS
720: .seealso: VecAXPY(), VecAYPX(), VecAXPBY()
721: @*/
722: PetscErrorCode VecWAXPY(Vec w,PetscScalar alpha,Vec x,Vec y)
723: {
737: if (w == y) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Result vector w cannot be same as input vector y, suggest VecAXPY()");
738: if (w == x) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Result vector w cannot be same as input vector x, suggest VecAYPX()");
741: PetscLogEventBegin(VEC_WAXPY,x,y,w,0);
742: (*w->ops->waxpy)(w,alpha,x,y);
743: PetscLogEventEnd(VEC_WAXPY,x,y,w,0);
744: PetscObjectStateIncrease((PetscObject)w);
745: return(0);
746: }
751: /*@
752: VecSetValues - Inserts or adds values into certain locations of a vector.
754: Not Collective
756: Input Parameters:
757: + x - vector to insert in
758: . ni - number of elements to add
759: . ix - indices where to add
760: . y - array of values
761: - iora - either INSERT_VALUES or ADD_VALUES, where
762: ADD_VALUES adds values to any existing entries, and
763: INSERT_VALUES replaces existing entries with new values
765: Notes:
766: VecSetValues() sets x[ix[i]] = y[i], for i=0,...,ni-1.
768: Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
769: options cannot be mixed without intervening calls to the assembly
770: routines.
772: These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd()
773: MUST be called after all calls to VecSetValues() have been completed.
775: VecSetValues() uses 0-based indices in Fortran as well as in C.
777: If you call VecSetOption(x, VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE),
778: negative indices may be passed in ix. These rows are
779: simply ignored. This allows easily inserting element load matrices
780: with homogeneous Dirchlet boundary conditions that you don't want represented
781: in the vector.
783: Level: beginner
785: Concepts: vector^setting values
787: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesLocal(),
788: VecSetValue(), VecSetValuesBlocked(), InsertMode, INSERT_VALUES, ADD_VALUES, VecGetValues()
789: @*/
790: PetscErrorCode VecSetValues(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
791: {
799: PetscLogEventBegin(VEC_SetValues,x,0,0,0);
800: (*x->ops->setvalues)(x,ni,ix,y,iora);
801: PetscLogEventEnd(VEC_SetValues,x,0,0,0);
802: PetscObjectStateIncrease((PetscObject)x);
803: return(0);
804: }
808: /*@
809: VecGetValues - Gets values from certain locations of a vector. Currently
810: can only get values on the same processor
812: Not Collective
813:
814: Input Parameters:
815: + x - vector to get values from
816: . ni - number of elements to get
817: - ix - indices where to get them from (in global 1d numbering)
819: Output Parameter:
820: . y - array of values
822: Notes:
823: The user provides the allocated array y; it is NOT allocated in this routine
825: VecGetValues() gets y[i] = x[ix[i]], for i=0,...,ni-1.
827: VecAssemblyBegin() and VecAssemblyEnd() MUST be called before calling this
829: VecGetValues() uses 0-based indices in Fortran as well as in C.
831: If you call VecSetOption(x, VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE),
832: negative indices may be passed in ix. These rows are
833: simply ignored.
835: Level: beginner
837: Concepts: vector^getting values
839: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecGetValuesLocal(),
840: VecGetValuesBlocked(), InsertMode, INSERT_VALUES, ADD_VALUES, VecSetValues()
841: @*/
842: PetscErrorCode VecGetValues(Vec x,PetscInt ni,const PetscInt ix[],PetscScalar y[])
843: {
851: (*x->ops->getvalues)(x,ni,ix,y);
852: return(0);
853: }
857: /*@
858: VecSetValuesBlocked - Inserts or adds blocks of values into certain locations of a vector.
860: Not Collective
862: Input Parameters:
863: + x - vector to insert in
864: . ni - number of blocks to add
865: . ix - indices where to add in block count, rather than element count
866: . y - array of values
867: - iora - either INSERT_VALUES or ADD_VALUES, where
868: ADD_VALUES adds values to any existing entries, and
869: INSERT_VALUES replaces existing entries with new values
871: Notes:
872: VecSetValuesBlocked() sets x[bs*ix[i]+j] = y[bs*i+j],
873: for j=0,...,bs, for i=0,...,ni-1. where bs was set with VecSetBlockSize().
875: Calls to VecSetValuesBlocked() with the INSERT_VALUES and ADD_VALUES
876: options cannot be mixed without intervening calls to the assembly
877: routines.
879: These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd()
880: MUST be called after all calls to VecSetValuesBlocked() have been completed.
882: VecSetValuesBlocked() uses 0-based indices in Fortran as well as in C.
884: Negative indices may be passed in ix, these rows are
885: simply ignored. This allows easily inserting element load matrices
886: with homogeneous Dirchlet boundary conditions that you don't want represented
887: in the vector.
889: Level: intermediate
891: Concepts: vector^setting values blocked
893: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesBlockedLocal(),
894: VecSetValues()
895: @*/
896: PetscErrorCode VecSetValuesBlocked(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
897: {
905: PetscLogEventBegin(VEC_SetValues,x,0,0,0);
906: (*x->ops->setvaluesblocked)(x,ni,ix,y,iora);
907: PetscLogEventEnd(VEC_SetValues,x,0,0,0);
908: PetscObjectStateIncrease((PetscObject)x);
909: return(0);
910: }
915: /*@
916: VecSetValuesLocal - Inserts or adds values into certain locations of a vector,
917: using a local ordering of the nodes.
919: Not Collective
921: Input Parameters:
922: + x - vector to insert in
923: . ni - number of elements to add
924: . ix - indices where to add
925: . y - array of values
926: - iora - either INSERT_VALUES or ADD_VALUES, where
927: ADD_VALUES adds values to any existing entries, and
928: INSERT_VALUES replaces existing entries with new values
930: Level: intermediate
932: Notes:
933: VecSetValuesLocal() sets x[ix[i]] = y[i], for i=0,...,ni-1.
935: Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
936: options cannot be mixed without intervening calls to the assembly
937: routines.
939: These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd()
940: MUST be called after all calls to VecSetValuesLocal() have been completed.
942: VecSetValuesLocal() uses 0-based indices in Fortran as well as in C.
944: Concepts: vector^setting values with local numbering
946: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetLocalToGlobalMapping(),
947: VecSetValuesBlockedLocal()
948: @*/
949: PetscErrorCode VecSetValuesLocal(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
950: {
952: PetscInt lixp[128],*lix = lixp;
960: PetscLogEventBegin(VEC_SetValues,x,0,0,0);
961: if (!x->ops->setvalueslocal) {
962: if (!x->map->mapping) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Local to global never set with VecSetLocalToGlobalMapping()");
963: if (ni > 128) {
964: PetscMalloc(ni*sizeof(PetscInt),&lix);
965: }
966: ISLocalToGlobalMappingApply(x->map->mapping,ni,(PetscInt*)ix,lix);
967: (*x->ops->setvalues)(x,ni,lix,y,iora);
968: if (ni > 128) {
969: PetscFree(lix);
970: }
971: } else {
972: (*x->ops->setvalueslocal)(x,ni,ix,y,iora);
973: }
974: PetscLogEventEnd(VEC_SetValues,x,0,0,0);
975: PetscObjectStateIncrease((PetscObject)x);
976: return(0);
977: }
981: /*@
982: VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
983: using a local ordering of the nodes.
985: Not Collective
987: Input Parameters:
988: + x - vector to insert in
989: . ni - number of blocks to add
990: . ix - indices where to add in block count, not element count
991: . y - array of values
992: - iora - either INSERT_VALUES or ADD_VALUES, where
993: ADD_VALUES adds values to any existing entries, and
994: INSERT_VALUES replaces existing entries with new values
996: Level: intermediate
998: Notes:
999: VecSetValuesBlockedLocal() sets x[bs*ix[i]+j] = y[bs*i+j],
1000: for j=0,..bs-1, for i=0,...,ni-1, where bs has been set with VecSetBlockSize().
1002: Calls to VecSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES
1003: options cannot be mixed without intervening calls to the assembly
1004: routines.
1006: These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd()
1007: MUST be called after all calls to VecSetValuesBlockedLocal() have been completed.
1009: VecSetValuesBlockedLocal() uses 0-based indices in Fortran as well as in C.
1012: Concepts: vector^setting values blocked with local numbering
1014: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetValuesBlocked(),
1015: VecSetLocalToGlobalMappingBlock()
1016: @*/
1017: PetscErrorCode VecSetValuesBlockedLocal(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
1018: {
1020: PetscInt lixp[128],*lix = lixp;
1027: if (!x->map->bmapping) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Local to global never set with VecSetLocalToGlobalMappingBlock()");
1028: if (ni > 128) {
1029: PetscMalloc(ni*sizeof(PetscInt),&lix);
1030: }
1032: PetscLogEventBegin(VEC_SetValues,x,0,0,0);
1033: ISLocalToGlobalMappingApply(x->map->bmapping,ni,(PetscInt*)ix,lix);
1034: (*x->ops->setvaluesblocked)(x,ni,lix,y,iora);
1035: PetscLogEventEnd(VEC_SetValues,x,0,0,0);
1036: if (ni > 128) {
1037: PetscFree(lix);
1038: }
1039: PetscObjectStateIncrease((PetscObject)x);
1040: return(0);
1041: }
1045: /*@
1046: VecMTDot - Computes indefinite vector multiple dot products.
1047: That is, it does NOT use the complex conjugate.
1049: Collective on Vec
1051: Input Parameters:
1052: + x - one vector
1053: . nv - number of vectors
1054: - y - array of vectors. Note that vectors are pointers
1056: Output Parameter:
1057: . val - array of the dot products
1059: Notes for Users of Complex Numbers:
1060: For complex vectors, VecMTDot() computes the indefinite form
1061: $ val = (x,y) = y^T x,
1062: where y^T denotes the transpose of y.
1064: Use VecMDot() for the inner product
1065: $ val = (x,y) = y^H x,
1066: where y^H denotes the conjugate transpose of y.
1068: Level: intermediate
1070: Concepts: inner product^multiple
1071: Concepts: vector^multiple inner products
1073: .seealso: VecMDot(), VecTDot()
1074: @*/
1075: PetscErrorCode VecMTDot(Vec x,PetscInt nv,const Vec y[],PetscScalar val[])
1076: {
1089: PetscLogEventBegin(VEC_MTDot,x,*y,0,0);
1090: (*x->ops->mtdot)(x,nv,y,val);
1091: PetscLogEventEnd(VEC_MTDot,x,*y,0,0);
1092: return(0);
1093: }
1097: /*@
1098: VecMDot - Computes vector multiple dot products.
1100: Collective on Vec
1102: Input Parameters:
1103: + x - one vector
1104: . nv - number of vectors
1105: - y - array of vectors.
1107: Output Parameter:
1108: . val - array of the dot products (does not allocate the array)
1110: Notes for Users of Complex Numbers:
1111: For complex vectors, VecMDot() computes
1112: $ val = (x,y) = y^H x,
1113: where y^H denotes the conjugate transpose of y.
1115: Use VecMTDot() for the indefinite form
1116: $ val = (x,y) = y^T x,
1117: where y^T denotes the transpose of y.
1119: Level: intermediate
1121: Concepts: inner product^multiple
1122: Concepts: vector^multiple inner products
1124: .seealso: VecMTDot(), VecDot()
1125: @*/
1126: PetscErrorCode VecMDot(Vec x,PetscInt nv,const Vec y[],PetscScalar val[])
1127: {
1132: if (!nv) return(0);
1133: if (nv < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of vectors (given %D) cannot be negative",nv);
1142: PetscLogEventBarrierBegin(VEC_MDotBarrier,x,*y,0,0,((PetscObject)x)->comm);
1143: (*x->ops->mdot)(x,nv,y,val);
1144: PetscLogEventBarrierEnd(VEC_MDotBarrier,x,*y,0,0,((PetscObject)x)->comm);
1145: return(0);
1146: }
1150: /*@
1151: VecMAXPY - Computes y = y + sum alpha[j] x[j]
1153: Logically Collective on Vec
1155: Input Parameters:
1156: + nv - number of scalars and x-vectors
1157: . alpha - array of scalars
1158: . y - one vector
1159: - x - array of vectors
1161: Level: intermediate
1163: Notes: y cannot be any of the x vectors
1165: Concepts: BLAS
1167: .seealso: VecAXPY(), VecWAXPY(), VecAYPX()
1168: @*/
1169: PetscErrorCode VecMAXPY(Vec y,PetscInt nv,const PetscScalar alpha[],Vec x[])
1170: {
1172: PetscInt i;
1176: if (!nv) return(0);
1177: if (nv < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of vectors (given %D) cannot be negative",nv);
1185: for (i=0; i<nv; i++) {
1187: }
1189: PetscLogEventBegin(VEC_MAXPY,*x,y,0,0);
1190: (*y->ops->maxpy)(y,nv,alpha,x);
1191: PetscLogEventEnd(VEC_MAXPY,*x,y,0,0);
1192: PetscObjectStateIncrease((PetscObject)y);
1193: return(0);
1194: }
1198: /*@
1199: VecGetSubVector - Gets a vector representing part of another vector
1201: Collective on IS (and Vec if nonlocal entries are needed)
1203: Input Arguments:
1204: + X - vector from which to extract a subvector
1205: - is - index set representing portion of X to extract
1207: Output Arguments:
1208: . Y - subvector corresponding to is
1210: Level: advanced
1212: Notes:
1213: The subvector Y should be returned with VecRestoreSubVector().
1215: This function may return a subvector without making a copy, therefore it is not safe to use the original vector while
1216: modifying the subvector. Other non-overlapping subvectors can still be obtained from X using this function.
1218: .seealso: MatGetSubMatrix()
1219: @*/
1220: PetscErrorCode VecGetSubVector(Vec X,IS is,Vec *Y)
1221: {
1223: Vec Z;
1224: PetscInt state;
1230: if (X->ops->getsubvector) {
1231: (*X->ops->getsubvector)(X,is,&Z);
1232: } else { /* Default implementation currently does no caching */
1233: PetscInt gstart,gend,start;
1234: PetscBool contiguous,gcontiguous;
1235: VecGetOwnershipRange(X,&gstart,&gend);
1236: ISContiguousLocal(is,gstart,gend,&start,&contiguous);
1237: MPI_Allreduce(&contiguous,&gcontiguous,1,MPI_INT,MPI_LAND,((PetscObject)is)->comm);
1238: if (gcontiguous) { /* We can do a no-copy implementation */
1239: PetscInt n,N;
1240: PetscScalar *x;
1241: PetscMPIInt size;
1242: ISGetLocalSize(is,&n);
1243: VecGetArray(X,&x);
1244: MPI_Comm_size(((PetscObject)X)->comm,&size);
1245: if (size == 1) {
1246: VecCreateSeqWithArray(((PetscObject)X)->comm,1,n,x+start,&Z);
1247: } else {
1248: ISGetSize(is,&N);
1249: VecCreateMPIWithArray(((PetscObject)X)->comm,1,n,N,x+start,&Z);
1250: }
1251: VecRestoreArray(X,&x);
1252: } else { /* Have to create a scatter and do a copy */
1253: VecScatter scatter;
1254: PetscInt n,N;
1255: ISGetLocalSize(is,&n);
1256: ISGetSize(is,&N);
1257: VecCreate(((PetscObject)is)->comm,&Z);
1258: VecSetSizes(Z,n,N);
1259: VecSetType(Z,((PetscObject)X)->type_name);
1260: VecScatterCreate(X,is,Z,PETSC_NULL,&scatter);
1261: VecScatterBegin(scatter,X,Z,INSERT_VALUES,SCATTER_FORWARD);
1262: VecScatterEnd(scatter,X,Z,INSERT_VALUES,SCATTER_FORWARD);
1263: VecScatterDestroy(&scatter);
1264: }
1265: }
1266: /* Record the state when the subvector was gotten so we know whether its values need to be put back */
1267: if (VecGetSubVectorSavedStateId < 0) {PetscObjectComposedDataRegister(&VecGetSubVectorSavedStateId);}
1268: PetscObjectStateQuery((PetscObject)Z,&state);
1269: PetscObjectComposedDataSetInt((PetscObject)Z,VecGetSubVectorSavedStateId,state);
1270: *Y = Z;
1271: return(0);
1272: }
1276: /*@
1277: VecRestoreSubVector - Restores a subvector extracted using VecGetSubVector()
1279: Collective on IS (and Vec if nonlocal entries need to be written)
1281: Input Arguments:
1282: + X - vector from which subvector was obtained
1283: . is - index set representing the subset of X
1284: - Y - subvector being restored
1286: Level: advanced
1288: .seealso: VecGetSubVector()
1289: @*/
1290: PetscErrorCode VecRestoreSubVector(Vec X,IS is,Vec *Y)
1291: {
1299: if (X->ops->restoresubvector) {
1300: (*X->ops->restoresubvector)(X,is,Y);
1301: } else {
1302: PetscInt savedstate=0,newstate;
1303: PetscBool valid;
1304: PetscObjectComposedDataGetInt((PetscObject)*Y,VecGetSubVectorSavedStateId,savedstate,valid);
1305: PetscObjectStateQuery((PetscObject)*Y,&newstate);
1306: if (valid && savedstate < newstate) {
1307: /* We might need to copy entries back, first check whether we have no-copy view */
1308: PetscInt gstart,gend,start;
1309: PetscBool contiguous,gcontiguous;
1310: VecGetOwnershipRange(X,&gstart,&gend);
1311: ISContiguousLocal(is,gstart,gend,&start,&contiguous);
1312: MPI_Allreduce(&contiguous,&gcontiguous,1,MPI_INT,MPI_LAND,((PetscObject)is)->comm);
1313: if (!gcontiguous) SETERRQ(((PetscObject)is)->comm,PETSC_ERR_SUP,"Unhandled case, values have been changed and need to be copied back into X");
1314: }
1315: VecDestroy(Y);
1316: }
1317: return(0);
1318: }
1320: /*MC
1321: VecGetArray - Returns a pointer to a contiguous array that contains this
1322: processor's portion of the vector data. For the standard PETSc
1323: vectors, VecGetArray() returns a pointer to the local data array and
1324: does not use any copies. If the underlying vector data is not stored
1325: in a contiquous array this routine will copy the data to a contiquous
1326: array and return a pointer to that. You MUST call VecRestoreArray()
1327: when you no longer need access to the array.
1329: Synopsis:
1330: PetscErrorCode VecGetArray(Vec x,PetscScalar *a[])
1332: Not Collective
1334: Input Parameter:
1335: . x - the vector
1337: Output Parameter:
1338: . a - location to put pointer to the array
1340: Fortran Note:
1341: This routine is used differently from Fortran 77
1342: $ Vec x
1343: $ PetscScalar x_array(1)
1344: $ PetscOffset i_x
1345: $ PetscErrorCode ierr
1346: $ call VecGetArray(x,x_array,i_x,ierr)
1347: $
1348: $ Access first local entry in vector with
1349: $ value = x_array(i_x + 1)
1350: $
1351: $ ...... other code
1352: $ call VecRestoreArray(x,x_array,i_x,ierr)
1353: For Fortran 90 see VecGetArrayF90()
1355: See the Fortran chapter of the users manual and
1356: petsc/src/snes/examples/tutorials/ex5f.F for details.
1358: Level: beginner
1360: Concepts: vector^accessing local values
1362: .seealso: VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(), VecGetArray2d()
1363: M*/
1368: /*@C
1369: VecGetArrays - Returns a pointer to the arrays in a set of vectors
1370: that were created by a call to VecDuplicateVecs(). You MUST call
1371: VecRestoreArrays() when you no longer need access to the array.
1373: Not Collective
1375: Input Parameter:
1376: + x - the vectors
1377: - n - the number of vectors
1379: Output Parameter:
1380: . a - location to put pointer to the array
1382: Fortran Note:
1383: This routine is not supported in Fortran.
1385: Level: intermediate
1387: .seealso: VecGetArray(), VecRestoreArrays()
1388: @*/
1389: PetscErrorCode VecGetArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1390: {
1392: PetscInt i;
1393: PetscScalar **q;
1399: if (n <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must get at least one array n = %D",n);
1400: PetscMalloc(n*sizeof(PetscScalar*),&q);
1401: for (i=0; i<n; ++i) {
1402: VecGetArray(x[i],&q[i]);
1403: }
1404: *a = q;
1405: return(0);
1406: }
1410: /*@C
1411: VecRestoreArrays - Restores a group of vectors after VecGetArrays()
1412: has been called.
1414: Not Collective
1416: Input Parameters:
1417: + x - the vector
1418: . n - the number of vectors
1419: - a - location of pointer to arrays obtained from VecGetArrays()
1421: Notes:
1422: For regular PETSc vectors this routine does not involve any copies. For
1423: any special vectors that do not store local vector data in a contiguous
1424: array, this routine will copy the data back into the underlying
1425: vector data structure from the arrays obtained with VecGetArrays().
1427: Fortran Note:
1428: This routine is not supported in Fortran.
1430: Level: intermediate
1432: .seealso: VecGetArrays(), VecRestoreArray()
1433: @*/
1434: PetscErrorCode VecRestoreArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1435: {
1437: PetscInt i;
1438: PetscScalar **q = *a;
1445: for(i=0;i<n;++i) {
1446: VecRestoreArray(x[i],&q[i]);
1447: }
1448: PetscFree(q);
1449: return(0);
1450: }
1452: /*MC
1453: VecRestoreArray - Restores a vector after VecGetArray() has been called.
1455: Synopsis:
1456: PetscErrorCode VecRestoreArray(Vec x,PetscScalar *a[])
1458: Not Collective
1460: Input Parameters:
1461: + x - the vector
1462: - a - location of pointer to array obtained from VecGetArray()
1464: Level: beginner
1466: Notes:
1467: For regular PETSc vectors this routine does not involve any copies. For
1468: any special vectors that do not store local vector data in a contiguous
1469: array, this routine will copy the data back into the underlying
1470: vector data structure from the array obtained with VecGetArray().
1472: This routine actually zeros out the a pointer. This is to prevent accidental
1473: us of the array after it has been restored. If you pass null for a it will
1474: not zero the array pointer a.
1476: Fortran Note:
1477: This routine is used differently from Fortran 77
1478: $ Vec x
1479: $ PetscScalar x_array(1)
1480: $ PetscOffset i_x
1481: $ PetscErrorCode ierr
1482: $ call VecGetArray(x,x_array,i_x,ierr)
1483: $
1484: $ Access first local entry in vector with
1485: $ value = x_array(i_x + 1)
1486: $
1487: $ ...... other code
1488: $ call VecRestoreArray(x,x_array,i_x,ierr)
1490: See the Fortran chapter of the users manual and
1491: petsc/src/snes/examples/tutorials/ex5f.F for details.
1492: For Fortran 90 see VecRestoreArrayF90()
1494: .seealso: VecGetArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(), VecRestoreArray2d()
1495: M*/
1499: /*@
1500: VecPlaceArray - Allows one to replace the array in a vector with an
1501: array provided by the user. This is useful to avoid copying an array
1502: into a vector.
1504: Not Collective
1506: Input Parameters:
1507: + vec - the vector
1508: - array - the array
1510: Notes:
1511: You can return to the original array with a call to VecResetArray()
1513: Level: developer
1515: .seealso: VecGetArray(), VecRestoreArray(), VecReplaceArray(), VecResetArray()
1517: @*/
1518: PetscErrorCode VecPlaceArray(Vec vec,const PetscScalar array[])
1519: {
1526: if (vec->ops->placearray) {
1527: (*vec->ops->placearray)(vec,array);
1528: } else SETERRQ(((PetscObject)vec)->comm,PETSC_ERR_SUP,"Cannot place array in this type of vector");
1529: PetscObjectStateIncrease((PetscObject)vec);
1530: return(0);
1531: }
1536: /*@C
1537: VecReplaceArray - Allows one to replace the array in a vector with an
1538: array provided by the user. This is useful to avoid copying an array
1539: into a vector.
1541: Not Collective
1543: Input Parameters:
1544: + vec - the vector
1545: - array - the array
1547: Notes:
1548: This permanently replaces the array and frees the memory associated
1549: with the old array.
1551: The memory passed in MUST be obtained with PetscMalloc() and CANNOT be
1552: freed by the user. It will be freed when the vector is destroy.
1554: Not supported from Fortran
1556: Level: developer
1558: .seealso: VecGetArray(), VecRestoreArray(), VecPlaceArray(), VecResetArray()
1560: @*/
1561: PetscErrorCode VecReplaceArray(Vec vec,const PetscScalar array[])
1562: {
1568: if (vec->ops->replacearray) {
1569: (*vec->ops->replacearray)(vec,array);
1570: } else SETERRQ(((PetscObject)vec)->comm,PETSC_ERR_SUP,"Cannot replace array in this type of vector");
1571: PetscObjectStateIncrease((PetscObject)vec);
1572: return(0);
1573: }
1575: /*MC
1576: VecDuplicateVecsF90 - Creates several vectors of the same type as an existing vector
1577: and makes them accessible via a Fortran90 pointer.
1579: Synopsis:
1580: VecDuplicateVecsF90(Vec x,PetscInt n,{Vec, pointer :: y(:)},integer ierr)
1582: Collective on Vec
1584: Input Parameters:
1585: + x - a vector to mimic
1586: - n - the number of vectors to obtain
1588: Output Parameters:
1589: + y - Fortran90 pointer to the array of vectors
1590: - ierr - error code
1592: Example of Usage:
1593: .vb
1594: Vec x
1595: Vec, pointer :: y(:)
1596: ....
1597: call VecDuplicateVecsF90(x,2,y,ierr)
1598: call VecSet(y(2),alpha,ierr)
1599: call VecSet(y(2),alpha,ierr)
1600: ....
1601: call VecDestroyVecsF90(2,y,ierr)
1602: .ve
1604: Notes:
1605: Not yet supported for all F90 compilers
1607: Use VecDestroyVecsF90() to free the space.
1609: Level: beginner
1611: .seealso: VecDestroyVecsF90(), VecDuplicateVecs()
1613: M*/
1615: /*MC
1616: VecRestoreArrayF90 - Restores a vector to a usable state after a call to
1617: VecGetArrayF90().
1619: Synopsis:
1620: VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
1622: Not collective
1624: Input Parameters:
1625: + x - vector
1626: - xx_v - the Fortran90 pointer to the array
1628: Output Parameter:
1629: . ierr - error code
1631: Example of Usage:
1632: .vb
1633: PetscScalar, pointer :: xx_v(:)
1634: ....
1635: call VecGetArrayF90(x,xx_v,ierr)
1636: a = xx_v(3)
1637: call VecRestoreArrayF90(x,xx_v,ierr)
1638: .ve
1639:
1640: Level: beginner
1642: .seealso: VecGetArrayF90(), VecGetArray(), VecRestoreArray(), UsingFortran
1644: M*/
1646: /*MC
1647: VecDestroyVecsF90 - Frees a block of vectors obtained with VecDuplicateVecsF90().
1649: Synopsis:
1650: VecDestroyVecsF90(PetscInt n,{Vec, pointer :: x(:)},PetscErrorCode ierr)
1652: Collective on Vec
1654: Input Parameters:
1655: + n - the number of vectors previously obtained
1656: - x - pointer to array of vector pointers
1658: Output Parameter:
1659: . ierr - error code
1661: Notes:
1662: Not yet supported for all F90 compilers
1664: Level: beginner
1666: .seealso: VecDestroyVecs(), VecDuplicateVecsF90()
1668: M*/
1670: /*MC
1671: VecGetArrayF90 - Accesses a vector array from Fortran90. For default PETSc
1672: vectors, VecGetArrayF90() returns a pointer to the local data array. Otherwise,
1673: this routine is implementation dependent. You MUST call VecRestoreArrayF90()
1674: when you no longer need access to the array.
1676: Synopsis:
1677: VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
1679: Not Collective
1681: Input Parameter:
1682: . x - vector
1684: Output Parameters:
1685: + xx_v - the Fortran90 pointer to the array
1686: - ierr - error code
1688: Example of Usage:
1689: .vb
1690: PetscScalar, pointer :: xx_v(:)
1691: ....
1692: call VecGetArrayF90(x,xx_v,ierr)
1693: a = xx_v(3)
1694: call VecRestoreArrayF90(x,xx_v,ierr)
1695: .ve
1697: Level: beginner
1699: .seealso: VecRestoreArrayF90(), VecGetArray(), VecRestoreArray(), UsingFortran
1701: M*/
1706: /*@C
1707: VecGetArray2d - Returns a pointer to a 2d contiguous array that contains this
1708: processor's portion of the vector data. You MUST call VecRestoreArray2d()
1709: when you no longer need access to the array.
1711: Not Collective
1713: Input Parameter:
1714: + x - the vector
1715: . m - first dimension of two dimensional array
1716: . n - second dimension of two dimensional array
1717: . mstart - first index you will use in first coordinate direction (often 0)
1718: - nstart - first index in the second coordinate direction (often 0)
1720: Output Parameter:
1721: . a - location to put pointer to the array
1723: Level: developer
1725: Notes:
1726: For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
1727: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
1728: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
1729: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().
1730:
1731: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
1733: Concepts: vector^accessing local values as 2d array
1735: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
1736: VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
1737: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
1738: @*/
1739: PetscErrorCode VecGetArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
1740: {
1742: PetscInt i,N;
1743: PetscScalar *aa;
1749: VecGetLocalSize(x,&N);
1750: 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);
1751: VecGetArray(x,&aa);
1753: PetscMalloc(m*sizeof(PetscScalar*),a);
1754: for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
1755: *a -= mstart;
1756: return(0);
1757: }
1761: /*@C
1762: VecRestoreArray2d - Restores a vector after VecGetArray2d() has been called.
1764: Not Collective
1766: Input Parameters:
1767: + x - the vector
1768: . m - first dimension of two dimensional array
1769: . n - second dimension of the two dimensional array
1770: . mstart - first index you will use in first coordinate direction (often 0)
1771: . nstart - first index in the second coordinate direction (often 0)
1772: - a - location of pointer to array obtained from VecGetArray2d()
1774: Level: developer
1776: Notes:
1777: For regular PETSc vectors this routine does not involve any copies. For
1778: any special vectors that do not store local vector data in a contiguous
1779: array, this routine will copy the data back into the underlying
1780: vector data structure from the array obtained with VecGetArray().
1782: This routine actually zeros out the a pointer.
1784: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
1785: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
1786: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
1787: @*/
1788: PetscErrorCode VecRestoreArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
1789: {
1791: void *dummy;
1797: dummy = (void*)(*a + mstart);
1798: PetscFree(dummy);
1799: VecRestoreArray(x,PETSC_NULL);
1800: return(0);
1801: }
1805: /*@C
1806: VecGetArray1d - Returns a pointer to a 1d contiguous array that contains this
1807: processor's portion of the vector data. You MUST call VecRestoreArray1d()
1808: when you no longer need access to the array.
1810: Not Collective
1812: Input Parameter:
1813: + x - the vector
1814: . m - first dimension of two dimensional array
1815: - mstart - first index you will use in first coordinate direction (often 0)
1817: Output Parameter:
1818: . a - location to put pointer to the array
1820: Level: developer
1822: Notes:
1823: For a vector obtained from DMCreateLocalVector() mstart are likely
1824: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
1825: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().
1826:
1827: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
1829: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
1830: VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
1831: VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
1832: @*/
1833: PetscErrorCode VecGetArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
1834: {
1836: PetscInt N;
1842: VecGetLocalSize(x,&N);
1843: if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
1844: VecGetArray(x,a);
1845: *a -= mstart;
1846: return(0);
1847: }
1851: /*@C
1852: VecRestoreArray1d - Restores a vector after VecGetArray1d() has been called.
1854: Not Collective
1856: Input Parameters:
1857: + x - the vector
1858: . m - first dimension of two dimensional array
1859: . mstart - first index you will use in first coordinate direction (often 0)
1860: - a - location of pointer to array obtained from VecGetArray21()
1862: Level: developer
1864: Notes:
1865: For regular PETSc vectors this routine does not involve any copies. For
1866: any special vectors that do not store local vector data in a contiguous
1867: array, this routine will copy the data back into the underlying
1868: vector data structure from the array obtained with VecGetArray1d().
1870: This routine actually zeros out the a pointer.
1872: Concepts: vector^accessing local values as 1d array
1874: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
1875: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
1876: VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
1877: @*/
1878: PetscErrorCode VecRestoreArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
1879: {
1885: VecRestoreArray(x,PETSC_NULL);
1886: return(0);
1887: }
1892: /*@C
1893: VecGetArray3d - Returns a pointer to a 3d contiguous array that contains this
1894: processor's portion of the vector data. You MUST call VecRestoreArray3d()
1895: when you no longer need access to the array.
1897: Not Collective
1899: Input Parameter:
1900: + x - the vector
1901: . m - first dimension of three dimensional array
1902: . n - second dimension of three dimensional array
1903: . p - third dimension of three dimensional array
1904: . mstart - first index you will use in first coordinate direction (often 0)
1905: . nstart - first index in the second coordinate direction (often 0)
1906: - pstart - first index in the third coordinate direction (often 0)
1908: Output Parameter:
1909: . a - location to put pointer to the array
1911: Level: developer
1913: Notes:
1914: For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
1915: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
1916: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
1917: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3d().
1918:
1919: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
1921: Concepts: vector^accessing local values as 3d array
1923: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
1924: VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
1925: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
1926: @*/
1927: PetscErrorCode VecGetArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
1928: {
1930: PetscInt i,N,j;
1931: PetscScalar *aa,**b;
1937: VecGetLocalSize(x,&N);
1938: 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);
1939: VecGetArray(x,&aa);
1941: PetscMalloc(m*sizeof(PetscScalar**)+m*n*sizeof(PetscScalar*),a);
1942: b = (PetscScalar **)((*a) + m);
1943: for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
1944: for (i=0; i<m; i++) {
1945: for (j=0; j<n; j++) {
1946: b[i*n+j] = aa + i*n*p + j*p - pstart;
1947: }
1948: }
1949: *a -= mstart;
1950: return(0);
1951: }
1955: /*@C
1956: VecRestoreArray3d - Restores a vector after VecGetArray3d() has been called.
1958: Not Collective
1960: Input Parameters:
1961: + x - the vector
1962: . m - first dimension of three dimensional array
1963: . n - second dimension of the three dimensional array
1964: . p - third dimension of the three dimensional array
1965: . mstart - first index you will use in first coordinate direction (often 0)
1966: . nstart - first index in the second coordinate direction (often 0)
1967: . pstart - first index in the third coordinate direction (often 0)
1968: - a - location of pointer to array obtained from VecGetArray3d()
1970: Level: developer
1972: Notes:
1973: For regular PETSc vectors this routine does not involve any copies. For
1974: any special vectors that do not store local vector data in a contiguous
1975: array, this routine will copy the data back into the underlying
1976: vector data structure from the array obtained with VecGetArray().
1978: This routine actually zeros out the a pointer.
1980: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
1981: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
1982: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
1983: @*/
1984: PetscErrorCode VecRestoreArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
1985: {
1987: void *dummy;
1993: dummy = (void*)(*a + mstart);
1994: PetscFree(dummy);
1995: VecRestoreArray(x,PETSC_NULL);
1996: return(0);
1997: }
2001: /*@C
2002: VecGetArray4d - Returns a pointer to a 4d contiguous array that contains this
2003: processor's portion of the vector data. You MUST call VecRestoreArray4d()
2004: when you no longer need access to the array.
2006: Not Collective
2008: Input Parameter:
2009: + x - the vector
2010: . m - first dimension of four dimensional array
2011: . n - second dimension of four dimensional array
2012: . p - third dimension of four dimensional array
2013: . q - fourth dimension of four dimensional array
2014: . mstart - first index you will use in first coordinate direction (often 0)
2015: . nstart - first index in the second coordinate direction (often 0)
2016: . pstart - first index in the third coordinate direction (often 0)
2017: - qstart - first index in the fourth coordinate direction (often 0)
2019: Output Parameter:
2020: . a - location to put pointer to the array
2022: Level: beginner
2024: Notes:
2025: For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
2026: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2027: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2028: the arguments from DMDAGet[Ghost}Corners() are reversed in the call to VecGetArray3d().
2029:
2030: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2032: Concepts: vector^accessing local values as 3d array
2034: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2035: VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2036: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2037: @*/
2038: PetscErrorCode VecGetArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
2039: {
2041: PetscInt i,N,j,k;
2042: PetscScalar *aa,***b,**c;
2048: VecGetLocalSize(x,&N);
2049: 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);
2050: VecGetArray(x,&aa);
2052: PetscMalloc(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p*sizeof(PetscScalar*),a);
2053: b = (PetscScalar ***)((*a) + m);
2054: c = (PetscScalar **)(b + m*n);
2055: for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
2056: for (i=0; i<m; i++) {
2057: for (j=0; j<n; j++) {
2058: b[i*n+j] = c + i*n*p + j*p - pstart;
2059: }
2060: }
2061: for (i=0; i<m; i++) {
2062: for (j=0; j<n; j++) {
2063: for (k=0; k<p; k++) {
2064: c[i*n*p+j*p+k] = aa + i*n*p*q + j*p*q + k*q - qstart;
2065: }
2066: }
2067: }
2068: *a -= mstart;
2069: return(0);
2070: }
2074: /*@C
2075: VecRestoreArray4d - Restores a vector after VecGetArray3d() has been called.
2077: Not Collective
2079: Input Parameters:
2080: + x - the vector
2081: . m - first dimension of four dimensional array
2082: . n - second dimension of the four dimensional array
2083: . p - third dimension of the four dimensional array
2084: . q - fourth dimension of the four dimensional array
2085: . mstart - first index you will use in first coordinate direction (often 0)
2086: . nstart - first index in the second coordinate direction (often 0)
2087: . pstart - first index in the third coordinate direction (often 0)
2088: . qstart - first index in the fourth coordinate direction (often 0)
2089: - a - location of pointer to array obtained from VecGetArray4d()
2091: Level: beginner
2093: Notes:
2094: For regular PETSc vectors this routine does not involve any copies. For
2095: any special vectors that do not store local vector data in a contiguous
2096: array, this routine will copy the data back into the underlying
2097: vector data structure from the array obtained with VecGetArray().
2099: This routine actually zeros out the a pointer.
2101: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2102: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2103: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
2104: @*/
2105: PetscErrorCode VecRestoreArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
2106: {
2108: void *dummy;
2114: dummy = (void*)(*a + mstart);
2115: PetscFree(dummy);
2116: VecRestoreArray(x,PETSC_NULL);
2117: return(0);
2118: }