Actual source code: vinv.c
petsc-3.3-p7 2013-05-11
2: /*
3: Some useful vector utility functions.
4: */
5: #include <petsc-private/vecimpl.h> /*I "petscvec.h" I*/
7: extern MPI_Op VecMax_Local_Op;
8: extern MPI_Op VecMin_Local_Op;
12: /*@
13: VecStrideSet - Sets a subvector of a vector defined
14: by a starting point and a stride with a given value
16: Logically Collective on Vec
18: Input Parameter:
19: + v - the vector
20: . start - starting point of the subvector (defined by a stride)
21: - s - value to multiply each subvector entry by
23: Notes:
24: One must call VecSetBlockSize() before this routine to set the stride
25: information, or use a vector created from a multicomponent DMDA.
27: This will only work if the desire subvector is a stride subvector
29: Level: advanced
31: Concepts: scale^on stride of vector
32: Concepts: stride^scale
34: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideScale()
35: @*/
36: PetscErrorCode VecStrideSet(Vec v,PetscInt start,PetscScalar s)
37: {
39: PetscInt i,n,bs;
40: PetscScalar *x;
47: VecGetLocalSize(v,&n);
48: VecGetArray(v,&x);
49: bs = v->map->bs;
50: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
51: else if (start >= bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Start of stride subvector (%D) is too large for stride\n Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
52: x += start;
54: for (i=0; i<n; i+=bs) {
55: x[i] = s;
56: }
57: x -= start;
59: VecRestoreArray(v,&x);
60: return(0);
61: }
65: /*@
66: VecStrideScale - Scales a subvector of a vector defined
67: by a starting point and a stride.
69: Logically Collective on Vec
71: Input Parameter:
72: + v - the vector
73: . start - starting point of the subvector (defined by a stride)
74: - scale - value to multiply each subvector entry by
76: Notes:
77: One must call VecSetBlockSize() before this routine to set the stride
78: information, or use a vector created from a multicomponent DMDA.
80: This will only work if the desire subvector is a stride subvector
82: Level: advanced
84: Concepts: scale^on stride of vector
85: Concepts: stride^scale
87: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideScale()
88: @*/
89: PetscErrorCode VecStrideScale(Vec v,PetscInt start,PetscScalar scale)
90: {
92: PetscInt i,n,bs;
93: PetscScalar *x;
100: VecGetLocalSize(v,&n);
101: VecGetArray(v,&x);
102: bs = v->map->bs;
103: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
104: else if (start >= bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Start of stride subvector (%D) is too large for stride\n Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
105: x += start;
107: for (i=0; i<n; i+=bs) {
108: x[i] *= scale;
109: }
110: x -= start;
112: VecRestoreArray(v,&x);
113: return(0);
114: }
118: /*@
119: VecStrideNorm - Computes the norm of subvector of a vector defined
120: by a starting point and a stride.
122: Collective on Vec
124: Input Parameter:
125: + v - the vector
126: . start - starting point of the subvector (defined by a stride)
127: - ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY
129: Output Parameter:
130: . norm - the norm
132: Notes:
133: One must call VecSetBlockSize() before this routine to set the stride
134: information, or use a vector created from a multicomponent DMDA.
136: If x is the array representing the vector x then this computes the norm
137: of the array (x[start],x[start+stride],x[start+2*stride], ....)
139: This is useful for computing, say the norm of the pressure variable when
140: the pressure is stored (interlaced) with other variables, say density etc.
142: This will only work if the desire subvector is a stride subvector
144: Level: advanced
146: Concepts: norm^on stride of vector
147: Concepts: stride^norm
149: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
150: @*/
151: PetscErrorCode VecStrideNorm(Vec v,PetscInt start,NormType ntype,PetscReal *nrm)
152: {
154: PetscInt i,n,bs;
155: PetscScalar *x;
156: PetscReal tnorm;
157: MPI_Comm comm;
162: VecGetLocalSize(v,&n);
163: VecGetArray(v,&x);
164: PetscObjectGetComm((PetscObject)v,&comm);
166: bs = v->map->bs;
167: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
168: else if (start >= bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Start of stride subvector (%D) is too large for stride\n Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
169: x += start;
171: if (ntype == NORM_2) {
172: PetscScalar sum = 0.0;
173: for (i=0; i<n; i+=bs) {
174: sum += x[i]*(PetscConj(x[i]));
175: }
176: tnorm = PetscRealPart(sum);
177: MPI_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_SUM,comm);
178: *nrm = PetscSqrtReal(*nrm);
179: } else if (ntype == NORM_1) {
180: tnorm = 0.0;
181: for (i=0; i<n; i+=bs) {
182: tnorm += PetscAbsScalar(x[i]);
183: }
184: MPI_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_SUM,comm);
185: } else if (ntype == NORM_INFINITY) {
186: PetscReal tmp;
187: tnorm = 0.0;
189: for (i=0; i<n; i+=bs) {
190: if ((tmp = PetscAbsScalar(x[i])) > tnorm) tnorm = tmp;
191: /* check special case of tmp == NaN */
192: if (tmp != tmp) {tnorm = tmp; break;}
193: }
194: MPI_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_MAX,comm);
195: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
196: VecRestoreArray(v,&x);
197: return(0);
198: }
202: /*@
203: VecStrideMax - Computes the maximum of subvector of a vector defined
204: by a starting point and a stride and optionally its location.
206: Collective on Vec
208: Input Parameter:
209: + v - the vector
210: - start - starting point of the subvector (defined by a stride)
212: Output Parameter:
213: + index - the location where the maximum occurred (pass PETSC_NULL if not required)
214: - nrm - the max
216: Notes:
217: One must call VecSetBlockSize() before this routine to set the stride
218: information, or use a vector created from a multicomponent DMDA.
220: If xa is the array representing the vector x, then this computes the max
221: of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)
223: This is useful for computing, say the maximum of the pressure variable when
224: the pressure is stored (interlaced) with other variables, e.g., density, etc.
225: This will only work if the desire subvector is a stride subvector.
227: Level: advanced
229: Concepts: maximum^on stride of vector
230: Concepts: stride^maximum
232: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
233: @*/
234: PetscErrorCode VecStrideMax(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
235: {
237: PetscInt i,n,bs,id;
238: PetscScalar *x;
239: PetscReal max,tmp;
240: MPI_Comm comm;
246: VecGetLocalSize(v,&n);
247: VecGetArray(v,&x);
248: PetscObjectGetComm((PetscObject)v,&comm);
250: bs = v->map->bs;
251: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
252: else if (start >= bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Start of stride subvector (%D) is too large for stride\n Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
253: x += start;
255: id = -1;
256: if (!n) {
257: max = PETSC_MIN_REAL;
258: } else {
259: id = 0;
260: #if defined(PETSC_USE_COMPLEX)
261: max = PetscRealPart(x[0]);
262: #else
263: max = x[0];
264: #endif
265: for (i=bs; i<n; i+=bs) {
266: #if defined(PETSC_USE_COMPLEX)
267: if ((tmp = PetscRealPart(x[i])) > max) { max = tmp; id = i;}
268: #else
269: if ((tmp = x[i]) > max) { max = tmp; id = i;}
270: #endif
271: }
272: }
273: VecRestoreArray(v,&x);
275: if (!idex) {
276: MPI_Allreduce(&max,nrm,1,MPIU_REAL,MPIU_MAX,comm);
277: } else {
278: PetscReal in[2],out[2];
279: PetscInt rstart;
281: VecGetOwnershipRange(v,&rstart,PETSC_NULL);
282: in[0] = max;
283: in[1] = rstart+id+start;
284: MPI_Allreduce(in,out,2,MPIU_REAL,VecMax_Local_Op,((PetscObject)v)->comm);
285: *nrm = out[0];
286: *idex = (PetscInt)out[1];
287: }
289: return(0);
290: }
294: /*@
295: VecStrideMin - Computes the minimum of subvector of a vector defined
296: by a starting point and a stride and optionally its location.
298: Collective on Vec
300: Input Parameter:
301: + v - the vector
302: - start - starting point of the subvector (defined by a stride)
304: Output Parameter:
305: + idex - the location where the minimum occurred. (pass PETSC_NULL if not required)
306: - nrm - the min
308: Level: advanced
310: Notes:
311: One must call VecSetBlockSize() before this routine to set the stride
312: information, or use a vector created from a multicomponent DMDA.
314: If xa is the array representing the vector x, then this computes the min
315: of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)
317: This is useful for computing, say the minimum of the pressure variable when
318: the pressure is stored (interlaced) with other variables, e.g., density, etc.
319: This will only work if the desire subvector is a stride subvector.
321: Concepts: minimum^on stride of vector
322: Concepts: stride^minimum
324: .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
325: @*/
326: PetscErrorCode VecStrideMin(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
327: {
329: PetscInt i,n,bs,id;
330: PetscScalar *x;
331: PetscReal min,tmp;
332: MPI_Comm comm;
338: VecGetLocalSize(v,&n);
339: VecGetArray(v,&x);
340: PetscObjectGetComm((PetscObject)v,&comm);
342: bs = v->map->bs;
343: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
344: else if (start >= bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Start of stride subvector (%D) is too large for stride\nHave you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
345: x += start;
347: id = -1;
348: if (!n) {
349: min = PETSC_MAX_REAL;
350: } else {
351: id = 0;
352: #if defined(PETSC_USE_COMPLEX)
353: min = PetscRealPart(x[0]);
354: #else
355: min = x[0];
356: #endif
357: for (i=bs; i<n; i+=bs) {
358: #if defined(PETSC_USE_COMPLEX)
359: if ((tmp = PetscRealPart(x[i])) < min) { min = tmp; id = i;}
360: #else
361: if ((tmp = x[i]) < min) { min = tmp; id = i;}
362: #endif
363: }
364: }
365: VecRestoreArray(v,&x);
367: if (!idex) {
368: MPI_Allreduce(&min,nrm,1,MPIU_REAL,MPIU_MIN,comm);
369: } else {
370: PetscReal in[2],out[2];
371: PetscInt rstart;
373: VecGetOwnershipRange(v,&rstart,PETSC_NULL);
374: in[0] = min;
375: in[1] = rstart+id;
376: MPI_Allreduce(in,out,2,MPIU_REAL,VecMin_Local_Op,((PetscObject)v)->comm);
377: *nrm = out[0];
378: *idex = (PetscInt)out[1];
379: }
381: return(0);
382: }
386: /*@
387: VecStrideScaleAll - Scales the subvectors of a vector defined
388: by a starting point and a stride.
390: Logically Collective on Vec
392: Input Parameter:
393: + v - the vector
394: - scales - values to multiply each subvector entry by
396: Notes:
397: One must call VecSetBlockSize() before this routine to set the stride
398: information, or use a vector created from a multicomponent DMDA.
401: Level: advanced
403: Concepts: scale^on stride of vector
404: Concepts: stride^scale
406: .seealso: VecNorm(), VecStrideScale(), VecScale(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
407: @*/
408: PetscErrorCode VecStrideScaleAll(Vec v,const PetscScalar *scales)
409: {
411: PetscInt i,j,n,bs;
412: PetscScalar *x;
417: VecGetLocalSize(v,&n);
418: VecGetArray(v,&x);
420: bs = v->map->bs;
422: /* need to provide optimized code for each bs */
423: for (i=0; i<n; i+=bs) {
424: for (j=0; j<bs; j++) {
425: x[i+j] *= scales[j];
426: }
427: }
428: VecRestoreArray(v,&x);
429: return(0);
430: }
434: /*@
435: VecStrideNormAll - Computes the norms subvectors of a vector defined
436: by a starting point and a stride.
438: Collective on Vec
440: Input Parameter:
441: + v - the vector
442: - ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY
444: Output Parameter:
445: . nrm - the norms
447: Notes:
448: One must call VecSetBlockSize() before this routine to set the stride
449: information, or use a vector created from a multicomponent DMDA.
451: If x is the array representing the vector x then this computes the norm
452: of the array (x[start],x[start+stride],x[start+2*stride], ....)
454: This is useful for computing, say the norm of the pressure variable when
455: the pressure is stored (interlaced) with other variables, say density etc.
457: This will only work if the desire subvector is a stride subvector
459: Level: advanced
461: Concepts: norm^on stride of vector
462: Concepts: stride^norm
464: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
465: @*/
466: PetscErrorCode VecStrideNormAll(Vec v,NormType ntype,PetscReal nrm[])
467: {
469: PetscInt i,j,n,bs;
470: PetscScalar *x;
471: PetscReal tnorm[128];
472: MPI_Comm comm;
477: VecGetLocalSize(v,&n);
478: VecGetArray(v,&x);
479: PetscObjectGetComm((PetscObject)v,&comm);
481: bs = v->map->bs;
482: if (bs > 128) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");
484: if (ntype == NORM_2) {
485: PetscScalar sum[128];
486: for (j=0; j<bs; j++) sum[j] = 0.0;
487: for (i=0; i<n; i+=bs) {
488: for (j=0; j<bs; j++) {
489: sum[j] += x[i+j]*(PetscConj(x[i+j]));
490: }
491: }
492: for (j=0; j<bs; j++) {
493: tnorm[j] = PetscRealPart(sum[j]);
494: }
495: MPI_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_SUM,comm);
496: for (j=0; j<bs; j++) {
497: nrm[j] = PetscSqrtReal(nrm[j]);
498: }
499: } else if (ntype == NORM_1) {
500: for (j=0; j<bs; j++) {
501: tnorm[j] = 0.0;
502: }
503: for (i=0; i<n; i+=bs) {
504: for (j=0; j<bs; j++) {
505: tnorm[j] += PetscAbsScalar(x[i+j]);
506: }
507: }
508: MPI_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_SUM,comm);
509: } else if (ntype == NORM_INFINITY) {
510: PetscReal tmp;
511: for (j=0; j<bs; j++) {
512: tnorm[j] = 0.0;
513: }
515: for (i=0; i<n; i+=bs) {
516: for (j=0; j<bs; j++) {
517: if ((tmp = PetscAbsScalar(x[i+j])) > tnorm[j]) tnorm[j] = tmp;
518: /* check special case of tmp == NaN */
519: if (tmp != tmp) {tnorm[j] = tmp; break;}
520: }
521: }
522: MPI_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_MAX,comm);
523: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
524: VecRestoreArray(v,&x);
525: return(0);
526: }
530: /*@
531: VecStrideMaxAll - Computes the maximums of subvectors of a vector defined
532: by a starting point and a stride and optionally its location.
534: Collective on Vec
536: Input Parameter:
537: . v - the vector
539: Output Parameter:
540: + index - the location where the maximum occurred (not supported, pass PETSC_NULL,
541: if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
542: - nrm - the maximums
544: Notes:
545: One must call VecSetBlockSize() before this routine to set the stride
546: information, or use a vector created from a multicomponent DMDA.
548: This is useful for computing, say the maximum of the pressure variable when
549: the pressure is stored (interlaced) with other variables, e.g., density, etc.
550: This will only work if the desire subvector is a stride subvector.
552: Level: advanced
554: Concepts: maximum^on stride of vector
555: Concepts: stride^maximum
557: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
558: @*/
559: PetscErrorCode VecStrideMaxAll(Vec v,PetscInt idex[],PetscReal nrm[])
560: {
562: PetscInt i,j,n,bs;
563: PetscScalar *x;
564: PetscReal max[128],tmp;
565: MPI_Comm comm;
570: if (idex) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support yet for returning index; send mail to petsc-maint@mcs.anl.gov asking for it");
571: VecGetLocalSize(v,&n);
572: VecGetArray(v,&x);
573: PetscObjectGetComm((PetscObject)v,&comm);
575: bs = v->map->bs;
576: if (bs > 128) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");
578: if (!n) {
579: for (j=0; j<bs; j++) {
580: max[j] = PETSC_MIN_REAL;
581: }
582: } else {
583: for (j=0; j<bs; j++) {
584: #if defined(PETSC_USE_COMPLEX)
585: max[j] = PetscRealPart(x[j]);
586: #else
587: max[j] = x[j];
588: #endif
589: }
590: for (i=bs; i<n; i+=bs) {
591: for (j=0; j<bs; j++) {
592: #if defined(PETSC_USE_COMPLEX)
593: if ((tmp = PetscRealPart(x[i+j])) > max[j]) { max[j] = tmp;}
594: #else
595: if ((tmp = x[i+j]) > max[j]) { max[j] = tmp; }
596: #endif
597: }
598: }
599: }
600: MPI_Allreduce(max,nrm,bs,MPIU_REAL,MPIU_MAX,comm);
602: VecRestoreArray(v,&x);
603: return(0);
604: }
608: /*@
609: VecStrideMinAll - Computes the minimum of subvector of a vector defined
610: by a starting point and a stride and optionally its location.
612: Collective on Vec
614: Input Parameter:
615: . v - the vector
617: Output Parameter:
618: + idex - the location where the minimum occurred (not supported, pass PETSC_NULL,
619: if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
620: - nrm - the minimums
622: Level: advanced
624: Notes:
625: One must call VecSetBlockSize() before this routine to set the stride
626: information, or use a vector created from a multicomponent DMDA.
628: This is useful for computing, say the minimum of the pressure variable when
629: the pressure is stored (interlaced) with other variables, e.g., density, etc.
630: This will only work if the desire subvector is a stride subvector.
632: Concepts: minimum^on stride of vector
633: Concepts: stride^minimum
635: .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
636: @*/
637: PetscErrorCode VecStrideMinAll(Vec v,PetscInt idex[],PetscReal nrm[])
638: {
640: PetscInt i,n,bs,j;
641: PetscScalar *x;
642: PetscReal min[128],tmp;
643: MPI_Comm comm;
648: if (idex) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support yet for returning index; send mail to petsc-maint@mcs.anl.gov asking for it");
649: VecGetLocalSize(v,&n);
650: VecGetArray(v,&x);
651: PetscObjectGetComm((PetscObject)v,&comm);
653: bs = v->map->bs;
654: if (bs > 128) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");
656: if (!n) {
657: for (j=0; j<bs; j++) {
658: min[j] = PETSC_MAX_REAL;
659: }
660: } else {
661: for (j=0; j<bs; j++) {
662: #if defined(PETSC_USE_COMPLEX)
663: min[j] = PetscRealPart(x[j]);
664: #else
665: min[j] = x[j];
666: #endif
667: }
668: for (i=bs; i<n; i+=bs) {
669: for (j=0; j<bs; j++) {
670: #if defined(PETSC_USE_COMPLEX)
671: if ((tmp = PetscRealPart(x[i+j])) < min[j]) { min[j] = tmp;}
672: #else
673: if ((tmp = x[i+j]) < min[j]) { min[j] = tmp; }
674: #endif
675: }
676: }
677: }
678: MPI_Allreduce(min,nrm,bs,MPIU_REAL,MPIU_MIN,comm);
680: VecRestoreArray(v,&x);
681: return(0);
682: }
684: /*----------------------------------------------------------------------------------------------*/
687: /*@
688: VecStrideGatherAll - Gathers all the single components from a multi-component vector into
689: separate vectors.
691: Collective on Vec
693: Input Parameter:
694: + v - the vector
695: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
697: Output Parameter:
698: . s - the location where the subvectors are stored
700: Notes:
701: One must call VecSetBlockSize() before this routine to set the stride
702: information, or use a vector created from a multicomponent DMDA.
704: If x is the array representing the vector x then this gathers
705: the arrays (x[start],x[start+stride],x[start+2*stride], ....)
706: for start=0,1,2,...bs-1
708: The parallel layout of the vector and the subvector must be the same;
709: i.e., nlocal of v = stride*(nlocal of s)
711: Not optimized; could be easily
713: Level: advanced
715: Concepts: gather^into strided vector
717: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
718: VecStrideScatterAll()
719: @*/
720: PetscErrorCode VecStrideGatherAll(Vec v,Vec s[],InsertMode addv)
721: {
723: PetscInt i,n,n2,bs,j,k,*bss = PETSC_NULL,nv,jj,nvc;
724: PetscScalar *x,**y;
730: VecGetLocalSize(v,&n);
731: VecGetLocalSize(s[0],&n2);
732: VecGetArray(v,&x);
733: bs = v->map->bs;
734: if (bs < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Input vector does not have a valid blocksize set");
736: PetscMalloc2(bs,PetscReal*,&y,bs,PetscInt,&bss);
737: nv = 0;
738: nvc = 0;
739: for (i=0; i<bs; i++) {
740: VecGetBlockSize(s[i],&bss[i]);
741: if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1 Re: [PETSC #8241] VecStrideGatherAll */
742: VecGetArray(s[i],&y[i]);
743: nvc += bss[i];
744: nv++;
745: if (nvc > bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of subvectors in subvectors > number of vectors in main vector");
746: if (nvc == bs) break;
747: }
749: n = n/bs;
751: jj = 0;
752: if (addv == INSERT_VALUES) {
753: for (j=0; j<nv; j++) {
754: for (k=0; k<bss[j]; k++) {
755: for (i=0; i<n; i++) {
756: y[j][i*bss[j] + k] = x[bs*i+jj+k];
757: }
758: }
759: jj += bss[j];
760: }
761: } else if (addv == ADD_VALUES) {
762: for (j=0; j<nv; j++) {
763: for (k=0; k<bss[j]; k++) {
764: for (i=0; i<n; i++) {
765: y[j][i*bss[j] + k] += x[bs*i+jj+k];
766: }
767: }
768: jj += bss[j];
769: }
770: #if !defined(PETSC_USE_COMPLEX)
771: } else if (addv == MAX_VALUES) {
772: for (j=0; j<nv; j++) {
773: for (k=0; k<bss[j]; k++) {
774: for (i=0; i<n; i++) {
775: y[j][i*bss[j] + k] = PetscMax(y[j][i*bss[j] + k],x[bs*i+jj+k]);
776: }
777: }
778: jj += bss[j];
779: }
780: #endif
781: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
783: VecRestoreArray(v,&x);
784: for (i=0; i<nv; i++) {
785: VecRestoreArray(s[i],&y[i]);
786: }
788: PetscFree2(y,bss);
789: return(0);
790: }
794: /*@
795: VecStrideScatterAll - Scatters all the single components from separate vectors into
796: a multi-component vector.
798: Collective on Vec
800: Input Parameter:
801: + s - the location where the subvectors are stored
802: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
804: Output Parameter:
805: . v - the multicomponent vector
807: Notes:
808: One must call VecSetBlockSize() before this routine to set the stride
809: information, or use a vector created from a multicomponent DMDA.
811: The parallel layout of the vector and the subvector must be the same;
812: i.e., nlocal of v = stride*(nlocal of s)
814: Not optimized; could be easily
816: Level: advanced
818: Concepts: scatter^into strided vector
820: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
821: VecStrideScatterAll()
822: @*/
823: PetscErrorCode VecStrideScatterAll(Vec s[],Vec v,InsertMode addv)
824: {
826: PetscInt i,n,n2,bs,j,jj,k,*bss = PETSC_NULL,nv,nvc;
827: PetscScalar *x,**y;
833: VecGetLocalSize(v,&n);
834: VecGetLocalSize(s[0],&n2);
835: VecGetArray(v,&x);
836: bs = v->map->bs;
837: if (bs < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Input vector does not have a valid blocksize set");
839: PetscMalloc2(bs,PetscScalar**,&y,bs,PetscInt,&bss);
840: nv = 0;
841: nvc = 0;
842: for (i=0; i<bs; i++) {
843: VecGetBlockSize(s[i],&bss[i]);
844: if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1 Re: [PETSC #8241] VecStrideGatherAll */
845: VecGetArray(s[i],&y[i]);
846: nvc += bss[i];
847: nv++;
848: if (nvc > bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of subvectors in subvectors > number of vectors in main vector");
849: if (nvc == bs) break;
850: }
852: n = n/bs;
854: jj = 0;
855: if (addv == INSERT_VALUES) {
856: for (j=0; j<nv; j++) {
857: for (k=0; k<bss[j]; k++) {
858: for (i=0; i<n; i++) {
859: x[bs*i+jj+k] = y[j][i*bss[j] + k];
860: }
861: }
862: jj += bss[j];
863: }
864: } else if (addv == ADD_VALUES) {
865: for (j=0; j<nv; j++) {
866: for (k=0; k<bss[j]; k++) {
867: for (i=0; i<n; i++) {
868: x[bs*i+jj+k] += y[j][i*bss[j] + k];
869: }
870: }
871: jj += bss[j];
872: }
873: #if !defined(PETSC_USE_COMPLEX)
874: } else if (addv == MAX_VALUES) {
875: for (j=0; j<nv; j++) {
876: for (k=0; k<bss[j]; k++) {
877: for (i=0; i<n; i++) {
878: x[bs*i+jj+k] = PetscMax(x[bs*i+jj+k],y[j][i*bss[j] + k]);
879: }
880: }
881: jj += bss[j];
882: }
883: #endif
884: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
886: VecRestoreArray(v,&x);
887: for (i=0; i<nv; i++) {
888: VecRestoreArray(s[i],&y[i]);
889: }
890: PetscFree2(y,bss);
891: return(0);
892: }
896: /*@
897: VecStrideGather - Gathers a single component from a multi-component vector into
898: another vector.
900: Collective on Vec
902: Input Parameter:
903: + v - the vector
904: . start - starting point of the subvector (defined by a stride)
905: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
907: Output Parameter:
908: . s - the location where the subvector is stored
910: Notes:
911: One must call VecSetBlockSize() before this routine to set the stride
912: information, or use a vector created from a multicomponent DMDA.
914: If x is the array representing the vector x then this gathers
915: the array (x[start],x[start+stride],x[start+2*stride], ....)
917: The parallel layout of the vector and the subvector must be the same;
918: i.e., nlocal of v = stride*(nlocal of s)
920: Not optimized; could be easily
922: Level: advanced
924: Concepts: gather^into strided vector
926: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
927: VecStrideScatterAll()
928: @*/
929: PetscErrorCode VecStrideGather(Vec v,PetscInt start,Vec s,InsertMode addv)
930: {
936: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
937: if (start >= v->map->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Start of stride subvector (%D) is too large for stride\n Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,v->map->bs);
938: if (!v->ops->stridegather) SETERRQ(((PetscObject)s)->comm,PETSC_ERR_SUP,"Not implemented for this Vec class");
939: (*v->ops->stridegather)(v,start,s,addv);
940: return(0);
941: }
945: /*@
946: VecStrideScatter - Scatters a single component from a vector into a multi-component vector.
948: Collective on Vec
950: Input Parameter:
951: + s - the single-component vector
952: . start - starting point of the subvector (defined by a stride)
953: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
955: Output Parameter:
956: . v - the location where the subvector is scattered (the multi-component vector)
958: Notes:
959: One must call VecSetBlockSize() on the multi-component vector before this
960: routine to set the stride information, or use a vector created from a multicomponent DMDA.
962: The parallel layout of the vector and the subvector must be the same;
963: i.e., nlocal of v = stride*(nlocal of s)
965: Not optimized; could be easily
967: Level: advanced
969: Concepts: scatter^into strided vector
971: .seealso: VecStrideNorm(), VecStrideGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
972: VecStrideScatterAll()
973: @*/
974: PetscErrorCode VecStrideScatter(Vec s,PetscInt start,Vec v,InsertMode addv)
975: {
981: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
982: if (start >= v->map->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Start of stride subvector (%D) is too large for stride\n Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,v->map->bs);
983: if (!v->ops->stridescatter) SETERRQ(((PetscObject)s)->comm,PETSC_ERR_SUP,"Not implemented for this Vec class");
984: (*v->ops->stridescatter)(s,start,v,addv);
985: return(0);
986: }
990: PetscErrorCode VecStrideGather_Default(Vec v,PetscInt start,Vec s,InsertMode addv)
991: {
993: PetscInt i,n,bs,ns;
994: PetscScalar *x,*y;
997: VecGetLocalSize(v,&n);
998: VecGetLocalSize(s,&ns);
999: VecGetArray(v,&x);
1000: VecGetArray(s,&y);
1002: bs = v->map->bs;
1003: if (n != ns*bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Subvector length * blocksize %D not correct for gather from original vector %D",ns*bs,n);
1004: x += start;
1005: n = n/bs;
1007: if (addv == INSERT_VALUES) {
1008: for (i=0; i<n; i++) {
1009: y[i] = x[bs*i];
1010: }
1011: } else if (addv == ADD_VALUES) {
1012: for (i=0; i<n; i++) {
1013: y[i] += x[bs*i];
1014: }
1015: #if !defined(PETSC_USE_COMPLEX)
1016: } else if (addv == MAX_VALUES) {
1017: for (i=0; i<n; i++) {
1018: y[i] = PetscMax(y[i],x[bs*i]);
1019: }
1020: #endif
1021: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
1023: VecRestoreArray(v,&x);
1024: VecRestoreArray(s,&y);
1025: return(0);
1026: }
1030: PetscErrorCode VecStrideScatter_Default(Vec s,PetscInt start,Vec v,InsertMode addv)
1031: {
1033: PetscInt i,n,bs,ns;
1034: PetscScalar *x,*y;
1037: VecGetLocalSize(v,&n);
1038: VecGetLocalSize(s,&ns);
1039: VecGetArray(v,&x);
1040: VecGetArray(s,&y);
1042: bs = v->map->bs;
1043: if (n != ns*bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Subvector length * blocksize %D not correct for scatter to multicomponent vector %D",ns*bs,n);
1044: x += start;
1045: n = n/bs;
1047: if (addv == INSERT_VALUES) {
1048: for (i=0; i<n; i++) {
1049: x[bs*i] = y[i];
1050: }
1051: } else if (addv == ADD_VALUES) {
1052: for (i=0; i<n; i++) {
1053: x[bs*i] += y[i];
1054: }
1055: #if !defined(PETSC_USE_COMPLEX)
1056: } else if (addv == MAX_VALUES) {
1057: for (i=0; i<n; i++) {
1058: x[bs*i] = PetscMax(y[i],x[bs*i]);
1059: }
1060: #endif
1061: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
1063: VecRestoreArray(v,&x);
1064: VecRestoreArray(s,&y);
1065: return(0);
1066: }
1070: PetscErrorCode VecReciprocal_Default(Vec v)
1071: {
1073: PetscInt i,n;
1074: PetscScalar *x;
1077: VecGetLocalSize(v,&n);
1078: VecGetArray(v,&x);
1079: for (i=0; i<n; i++) {
1080: if (x[i] != (PetscScalar)0.0) x[i] = (PetscScalar)1.0/x[i];
1081: }
1082: VecRestoreArray(v,&x);
1083: return(0);
1084: }
1088: /*@
1089: VecExp - Replaces each component of a vector by e^x_i
1091: Not collective
1093: Input Parameter:
1094: . v - The vector
1096: Output Parameter:
1097: . v - The vector of exponents
1099: Level: beginner
1101: .seealso: VecLog(), VecAbs(), VecSqrtAbs(), VecReciprocal()
1103: .keywords: vector, sqrt, square root
1104: @*/
1105: PetscErrorCode VecExp(Vec v)
1106: {
1107: PetscScalar *x;
1108: PetscInt i, n;
1113: if (v->ops->exp) {
1114: (*v->ops->exp)(v);
1115: } else {
1116: VecGetLocalSize(v, &n);
1117: VecGetArray(v, &x);
1118: for(i = 0; i < n; i++) {
1119: x[i] = PetscExpScalar(x[i]);
1120: }
1121: VecRestoreArray(v, &x);
1122: }
1123: return(0);
1124: }
1128: /*@
1129: VecLog - Replaces each component of a vector by log(x_i), the natural logarithm
1131: Not collective
1133: Input Parameter:
1134: . v - The vector
1136: Output Parameter:
1137: . v - The vector of logs
1139: Level: beginner
1141: .seealso: VecExp(), VecAbs(), VecSqrtAbs(), VecReciprocal()
1143: .keywords: vector, sqrt, square root
1144: @*/
1145: PetscErrorCode VecLog(Vec v)
1146: {
1147: PetscScalar *x;
1148: PetscInt i, n;
1153: if (v->ops->log) {
1154: (*v->ops->log)(v);
1155: } else {
1156: VecGetLocalSize(v, &n);
1157: VecGetArray(v, &x);
1158: for(i = 0; i < n; i++) {
1159: x[i] = PetscLogScalar(x[i]);
1160: }
1161: VecRestoreArray(v, &x);
1162: }
1163: return(0);
1164: }
1168: /*@
1169: VecSqrtAbs - Replaces each component of a vector by the square root of its magnitude.
1171: Not collective
1173: Input Parameter:
1174: . v - The vector
1176: Output Parameter:
1177: . v - The vector square root
1179: Level: beginner
1181: Note: The actual function is sqrt(|x_i|)
1183: .seealso: VecLog(), VecExp(), VecReciprocal(), VecAbs()
1185: .keywords: vector, sqrt, square root
1186: @*/
1187: PetscErrorCode VecSqrtAbs(Vec v)
1188: {
1189: PetscScalar *x;
1190: PetscInt i, n;
1195: if (v->ops->sqrt) {
1196: (*v->ops->sqrt)(v);
1197: } else {
1198: VecGetLocalSize(v, &n);
1199: VecGetArray(v, &x);
1200: for(i = 0; i < n; i++) {
1201: x[i] = PetscSqrtReal(PetscAbsScalar(x[i]));
1202: }
1203: VecRestoreArray(v, &x);
1204: }
1205: return(0);
1206: }
1210: /*@
1211: VecDotNorm2 - computes the inner product of two vectors and the 2-norm squared of the second vector
1213: Collective on Vec
1215: Input Parameter:
1216: + s - first vector
1217: - t - second vector
1219: Output Parameter:
1220: + dp - s't
1221: - nm - t't
1223: Level: advanced
1225: Developer Notes: Even though the second return argument is a norm and hence could be a PetscReal value it is returned as PetscScalar
1227: .seealso: VecDot(), VecNorm(), VecDotBegin(), VecNormBegin(), VecDotEnd(), VecNormEnd()
1229: .keywords: vector, sqrt, square root
1230: @*/
1231: PetscErrorCode VecDotNorm2(Vec s,Vec t,PetscScalar *dp, PetscScalar *nm)
1232: {
1233: PetscScalar *sx, *tx, dpx = 0.0, nmx = 0.0,work[2],sum[2];
1234: PetscInt i, n;
1245: if (s->map->N != t->map->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
1246: if (s->map->n != t->map->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");
1248: PetscLogEventBarrierBegin(VEC_DotNormBarrier,s,t,0,0,((PetscObject)s)->comm);
1249: if (s->ops->dotnorm2) {
1250: (*s->ops->dotnorm2)(s,t,dp,nm);
1251: } else {
1252: VecGetLocalSize(s, &n);
1253: VecGetArray(s, &sx);
1254: VecGetArray(t, &tx);
1256: for (i = 0; i<n; i++) {
1257: dpx += sx[i]*PetscConj(tx[i]);
1258: nmx += tx[i]*PetscConj(tx[i]);
1259: }
1260: work[0] = dpx;
1261: work[1] = nmx;
1262: MPI_Allreduce(&work,&sum,2,MPIU_SCALAR,MPIU_SUM,((PetscObject)s)->comm);
1263: *dp = sum[0];
1264: *nm = sum[1];
1266: VecRestoreArray(t, &tx);
1267: VecRestoreArray(s, &sx);
1268: PetscLogFlops(4.0*n);
1269: }
1270: PetscLogEventBarrierEnd(VEC_DotNormBarrier,s,t,0,0,((PetscObject)s)->comm);
1271: return(0);
1272: }
1276: /*@
1277: VecSum - Computes the sum of all the components of a vector.
1279: Collective on Vec
1281: Input Parameter:
1282: . v - the vector
1284: Output Parameter:
1285: . sum - the result
1287: Level: beginner
1289: Concepts: sum^of vector entries
1291: .seealso: VecNorm()
1292: @*/
1293: PetscErrorCode VecSum(Vec v,PetscScalar *sum)
1294: {
1296: PetscInt i,n;
1297: PetscScalar *x,lsum = 0.0;
1302: VecGetLocalSize(v,&n);
1303: VecGetArray(v,&x);
1304: for (i=0; i<n; i++) {
1305: lsum += x[i];
1306: }
1307: MPI_Allreduce(&lsum,sum,1,MPIU_SCALAR,MPIU_SUM,((PetscObject)v)->comm);
1308: VecRestoreArray(v,&x);
1309: return(0);
1310: }
1314: /*@
1315: VecShift - Shifts all of the components of a vector by computing
1316: x[i] = x[i] + shift.
1318: Logically Collective on Vec
1320: Input Parameters:
1321: + v - the vector
1322: - shift - the shift
1324: Output Parameter:
1325: . v - the shifted vector
1327: Level: intermediate
1329: Concepts: vector^adding constant
1331: @*/
1332: PetscErrorCode VecShift(Vec v,PetscScalar shift)
1333: {
1335: PetscInt i,n;
1336: PetscScalar *x;
1341: if (v->ops->shift) {
1342: (*v->ops->shift)(v);
1343: } else {
1344: VecGetLocalSize(v,&n);
1345: VecGetArray(v,&x);
1346: for (i=0; i<n; i++) {
1347: x[i] += shift;
1348: }
1349: VecRestoreArray(v,&x);
1350: }
1351: return(0);
1352: }
1356: /*@
1357: VecAbs - Replaces every element in a vector with its absolute value.
1359: Logically Collective on Vec
1361: Input Parameters:
1362: . v - the vector
1364: Level: intermediate
1366: Concepts: vector^absolute value
1368: @*/
1369: PetscErrorCode VecAbs(Vec v)
1370: {
1372: PetscInt i,n;
1373: PetscScalar *x;
1377: if (v->ops->abs) {
1378: (*v->ops->abs)(v);
1379: } else {
1380: VecGetLocalSize(v,&n);
1381: VecGetArray(v,&x);
1382: for (i=0; i<n; i++) {
1383: x[i] = PetscAbsScalar(x[i]);
1384: }
1385: VecRestoreArray(v,&x);
1386: }
1387: return(0);
1388: }
1392: /*@
1393: VecPermute - Permutes a vector in place using the given ordering.
1395: Input Parameters:
1396: + vec - The vector
1397: . order - The ordering
1398: - inv - The flag for inverting the permutation
1400: Level: beginner
1402: Note: This function does not yet support parallel Index Sets
1404: .seealso: MatPermute()
1405: .keywords: vec, permute
1406: @*/
1407: PetscErrorCode VecPermute(Vec x, IS row, PetscBool inv)
1408: {
1409: PetscScalar *array, *newArray;
1410: const PetscInt *idx;
1411: PetscInt i;
1415: ISGetIndices(row, &idx);
1416: VecGetArray(x, &array);
1417: PetscMalloc(x->map->n*sizeof(PetscScalar), &newArray);
1418: #ifdef PETSC_USE_DEBUG
1419: for(i = 0; i < x->map->n; i++) {
1420: if ((idx[i] < 0) || (idx[i] >= x->map->n)) {
1421: SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT, "Permutation index %D is out of bounds: %D", i, idx[i]);
1422: }
1423: }
1424: #endif
1425: if (!inv) {
1426: for(i = 0; i < x->map->n; i++) newArray[i] = array[idx[i]];
1427: } else {
1428: for(i = 0; i < x->map->n; i++) newArray[idx[i]] = array[i];
1429: }
1430: VecRestoreArray(x, &array);
1431: ISRestoreIndices(row, &idx);
1432: VecReplaceArray(x, newArray);
1433: return(0);
1434: }
1438: /*@
1439: VecEqual - Compares two vectors.
1441: Collective on Vec
1443: Input Parameters:
1444: + vec1 - the first vector
1445: - vec2 - the second vector
1447: Output Parameter:
1448: . flg - PETSC_TRUE if the vectors are equal; PETSC_FALSE otherwise.
1450: Level: intermediate
1452: Concepts: equal^two vectors
1453: Concepts: vector^equality
1455: @*/
1456: PetscErrorCode VecEqual(Vec vec1,Vec vec2,PetscBool *flg)
1457: {
1458: PetscScalar *v1,*v2;
1460: PetscInt n1,n2,N1,N2;
1461: PetscInt state1,state2;
1462: PetscBool flg1;
1463:
1468: if (vec1 == vec2) {
1469: *flg = PETSC_TRUE;
1470: } else {
1471: VecGetSize(vec1,&N1);
1472: VecGetSize(vec2,&N2);
1473: if (N1 != N2) {
1474: flg1 = PETSC_FALSE;
1475: } else {
1476: VecGetLocalSize(vec1,&n1);
1477: VecGetLocalSize(vec2,&n2);
1478: if (n1 != n2) {
1479: flg1 = PETSC_FALSE;
1480: } else {
1481: PetscObjectStateQuery((PetscObject) vec1,&state1);
1482: PetscObjectStateQuery((PetscObject) vec2,&state2);
1483: VecGetArray(vec1,&v1);
1484: VecGetArray(vec2,&v2);
1485: #if defined(PETSC_USE_COMPLEX)
1486: {
1487: PetscInt k;
1488: flg1 = PETSC_TRUE;
1489: for (k=0; k<n1; k++){
1490: if (PetscRealPart(v1[k]) != PetscRealPart(v2[k]) || PetscImaginaryPart(v1[k]) != PetscImaginaryPart(v2[k])){
1491: flg1 = PETSC_FALSE;
1492: break;
1493: }
1494: }
1495: }
1496: #else
1497: PetscMemcmp(v1,v2,n1*sizeof(PetscScalar),&flg1);
1498: #endif
1499: VecRestoreArray(vec1,&v1);
1500: VecRestoreArray(vec2,&v2);
1501: PetscObjectSetState((PetscObject) vec1,state1);
1502: PetscObjectSetState((PetscObject) vec2,state2);
1503: }
1504: }
1505: /* combine results from all processors */
1506: MPI_Allreduce(&flg1,flg,1,MPI_INT,MPI_MIN,((PetscObject)vec1)->comm);
1507: }
1508: return(0);
1509: }