Actual source code: vinv.c
1: /*$Id: vinv.c,v 1.71 2001/09/11 16:31:37 bsmith Exp $*/
2: /*
3: Some useful vector utility functions.
4: */
5: #include src/vec/vecimpl.h
7: #undef __FUNCT__
9: /*@C
10: VecStrideNorm - Computes the norm of subvector of a vector defined
11: by a starting point and a stride.
13: Collective on Vec
15: Input Parameter:
16: + v - the vector
17: . start - starting point of the subvector (defined by a stride)
18: - ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY
20: Output Parameter:
21: . norm - the norm
23: Notes:
24: One must call VecSetBlockSize() before this routine to set the stride
25: information, or use a vector created from a multicomponent DA.
27: If x is the array representing the vector x then this computes the norm
28: of the array (x[start],x[start+stride],x[start+2*stride], ....)
30: This is useful for computing, say the norm of the pressure variable when
31: the pressure is stored (interlaced) with other variables, say density etc.
33: This will only work if the desire subvector is a stride subvector
35: Level: advanced
37: Concepts: norm^on stride of vector
38: Concepts: stride^norm
40: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
41: @*/
42: int VecStrideNorm(Vec v,int start,NormType ntype,PetscReal *nrm)
43: {
44: int i,n,ierr,bs;
45: PetscScalar *x;
46: PetscReal tnorm;
47: MPI_Comm comm;
51: VecGetLocalSize(v,&n);
52: VecGetArray(v,&x);
53: PetscObjectGetComm((PetscObject)v,&comm);
55: bs = v->bs;
56: if (start >= bs) {
57: SETERRQ2(1,"Start of stride subvector (%d) is too large for striden
58: Have you set the vector blocksize (%d) correctly with VecSetBlockSize()?",start,bs);
59: }
60: x += start;
62: if (ntype == NORM_2) {
63: PetscScalar sum = 0.0;
64: for (i=0; i<n; i+=bs) {
65: sum += x[i]*(PetscConj(x[i]));
66: }
67: tnorm = PetscRealPart(sum);
68: ierr = MPI_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPI_SUM,comm);
69: *nrm = sqrt(*nrm);
70: } else if (ntype == NORM_1) {
71: tnorm = 0.0;
72: for (i=0; i<n; i+=bs) {
73: tnorm += PetscAbsScalar(x[i]);
74: }
75: ierr = MPI_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPI_SUM,comm);
76: } else if (ntype == NORM_INFINITY) {
77: PetscReal tmp;
78: tnorm = 0.0;
80: for (i=0; i<n; i+=bs) {
81: if ((tmp = PetscAbsScalar(x[i])) > tnorm) tnorm = tmp;
82: /* check special case of tmp == NaN */
83: if (tmp != tmp) {tnorm = tmp; break;}
84: }
85: ierr = MPI_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPI_MAX,comm);
86: } else {
87: SETERRQ(1,"Unknown norm type");
88: }
90: VecRestoreArray(v,&x);
91: return(0);
92: }
94: #undef __FUNCT__
96: /*@C
97: VecStrideMax - Computes the maximum of subvector of a vector defined
98: by a starting point and a stride and optionally its location.
100: Collective on Vec
102: Input Parameter:
103: + v - the vector
104: - start - starting point of the subvector (defined by a stride)
106: Output Parameter:
107: + index - the location where the maximum occurred (not supported, pass PETSC_NULL,
108: if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
109: - nrm - the max
111: Notes:
112: One must call VecSetBlockSize() before this routine to set the stride
113: information, or use a vector created from a multicomponent DA.
115: If xa is the array representing the vector x, then this computes the max
116: of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)
118: This is useful for computing, say the maximum of the pressure variable when
119: the pressure is stored (interlaced) with other variables, e.g., density, etc.
120: This will only work if the desire subvector is a stride subvector.
122: Level: advanced
124: Concepts: maximum^on stride of vector
125: Concepts: stride^maximum
127: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
128: @*/
129: int VecStrideMax(Vec v,int start,int *idex,PetscReal *nrm)
130: {
131: int i,n,ierr,bs;
132: PetscScalar *x;
133: PetscReal max,tmp;
134: MPI_Comm comm;
138: if (idex) {
139: SETERRQ(1,"No support yet for returning index; send mail to petsc-maint@mcs.anl.gov asking for it");
140: }
141: VecGetLocalSize(v,&n);
142: VecGetArray(v,&x);
143: PetscObjectGetComm((PetscObject)v,&comm);
145: bs = v->bs;
146: if (start >= bs) {
147: SETERRQ2(1,"Start of stride subvector (%d) is too large for striden
148: Have you set the vector blocksize (%d) correctly with VecSetBlockSize()?",start,bs);
149: }
150: x += start;
152: if (!n) {
153: max = PETSC_MIN;
154: } else {
155: #if defined(PETSC_USE_COMPLEX)
156: max = PetscRealPart(x[0]);
157: #else
158: max = x[0];
159: #endif
160: for (i=bs; i<n; i+=bs) {
161: #if defined(PETSC_USE_COMPLEX)
162: if ((tmp = PetscRealPart(x[i])) > max) { max = tmp;}
163: #else
164: if ((tmp = x[i]) > max) { max = tmp; }
165: #endif
166: }
167: }
168: ierr = MPI_Allreduce(&max,nrm,1,MPIU_REAL,MPI_MAX,comm);
170: VecRestoreArray(v,&x);
171: return(0);
172: }
174: #undef __FUNCT__
176: /*@C
177: VecStrideMin - Computes the minimum of subvector of a vector defined
178: by a starting point and a stride and optionally its location.
180: Collective on Vec
182: Input Parameter:
183: + v - the vector
184: - start - starting point of the subvector (defined by a stride)
186: Output Parameter:
187: + idex - the location where the minimum occurred (not supported, pass PETSC_NULL,
188: if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
189: - nrm - the min
191: Level: advanced
193: Notes:
194: One must call VecSetBlockSize() before this routine to set the stride
195: information, or use a vector created from a multicomponent DA.
197: If xa is the array representing the vector x, then this computes the min
198: of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)
200: This is useful for computing, say the minimum of the pressure variable when
201: the pressure is stored (interlaced) with other variables, e.g., density, etc.
202: This will only work if the desire subvector is a stride subvector.
204: Concepts: minimum^on stride of vector
205: Concepts: stride^minimum
207: .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
208: @*/
209: int VecStrideMin(Vec v,int start,int *idex,PetscReal *nrm)
210: {
211: int i,n,ierr,bs;
212: PetscScalar *x;
213: PetscReal min,tmp;
214: MPI_Comm comm;
218: if (idex) {
219: SETERRQ(1,"No support yet for returning index; send mail to petsc-maint@mcs.anl.gov asking for it");
220: }
221: VecGetLocalSize(v,&n);
222: VecGetArray(v,&x);
223: PetscObjectGetComm((PetscObject)v,&comm);
225: bs = v->bs;
226: if (start >= bs) {
227: SETERRQ2(1,"Start of stride subvector (%d) is too large for striden
228: Have you set the vector blocksize (%d) correctly with VecSetBlockSize()?",start,bs);
229: }
230: x += start;
232: if (!n) {
233: min = PETSC_MAX;
234: } else {
235: #if defined(PETSC_USE_COMPLEX)
236: min = PetscRealPart(x[0]);
237: #else
238: min = x[0];
239: #endif
240: for (i=bs; i<n; i+=bs) {
241: #if defined(PETSC_USE_COMPLEX)
242: if ((tmp = PetscRealPart(x[i])) < min) { min = tmp;}
243: #else
244: if ((tmp = x[i]) < min) { min = tmp; }
245: #endif
246: }
247: }
248: ierr = MPI_Allreduce(&min,nrm,1,MPIU_REAL,MPI_MIN,comm);
250: VecRestoreArray(v,&x);
251: return(0);
252: }
254: #undef __FUNCT__
256: /*@
257: VecStrideGatherAll - Gathers all the single components from a multi-component vector into
258: seperate vectors.
260: Collective on Vec
262: Input Parameter:
263: + v - the vector
264: - addv - one of ADD_VALUES,SET_VALUES,MAX_VALUES
266: Output Parameter:
267: . s - the location where the subvectors are stored
269: Notes:
270: One must call VecSetBlockSize() before this routine to set the stride
271: information, or use a vector created from a multicomponent DA.
273: If x is the array representing the vector x then this gathers
274: the arrays (x[start],x[start+stride],x[start+2*stride], ....)
275: for start=0,1,2,...bs-1
277: The parallel layout of the vector and the subvector must be the same;
278: i.e., nlocal of v = stride*(nlocal of s)
280: Level: advanced
282: Concepts: gather^into strided vector
284: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
285: VecStrideScatterAll()
286: @*/
287: int VecStrideGatherAll(Vec v,Vec *s,InsertMode addv)
288: {
289: int i,n,ierr,bs,ns,j;
290: PetscScalar *x,**y;
295: VecGetLocalSize(v,&n);
296: VecGetLocalSize(*s,&ns);
297: VecGetArray(v,&x);
298: bs = v->bs;
300: PetscMalloc(bs*sizeof(PetscReal*),&y);
301: for (i=0; i<bs; i++) {
302: VecGetArray(s[i],&y[i]);
303: }
305: if (n != ns*bs) {
306: SETERRQ2(1,"Subvector length * blocksize %d not correct for gather from original vector %d",ns*bs,n);
307: }
308: n = n/bs;
310: if (addv == INSERT_VALUES) {
311: for (j=0; j<bs; j++) {
312: for (i=0; i<n; i++) {
313: y[j][i] = x[bs*i+j];
314: }
315: }
316: } else if (addv == ADD_VALUES) {
317: for (j=0; j<bs; j++) {
318: for (i=0; i<n; i++) {
319: y[j][i] += x[bs*i+j];
320: }
321: }
322: #if !defined(PETSC_USE_COMPLEX)
323: } else if (addv == MAX_VALUES) {
324: for (j=0; j<bs; j++) {
325: for (i=0; i<n; i++) {
326: y[j][i] = PetscMax(y[j][i],x[bs*i+j]);
327: }
328: }
329: #endif
330: } else {
331: SETERRQ(1,"Unknown insert type");
332: }
334: VecRestoreArray(v,&x);
335: for (i=0; i<bs; i++) {
336: VecRestoreArray(s[i],&y[i]);
337: }
338: PetscFree(y);
339: return(0);
340: }
341: #undef __FUNCT__
343: /*@
344: VecStrideScatterAll - Scatters all the single components from seperate vectors into
345: a multi-component vector.
347: Collective on Vec
349: Input Parameter:
350: + s - the location where the subvectors are stored
351: - addv - one of ADD_VALUES,SET_VALUES,MAX_VALUES
353: Output Parameter:
354: . v - the multicomponent vector
356: Notes:
357: One must call VecSetBlockSize() before this routine to set the stride
358: information, or use a vector created from a multicomponent DA.
360: The parallel layout of the vector and the subvector must be the same;
361: i.e., nlocal of v = stride*(nlocal of s)
363: Level: advanced
365: Concepts: scatter^into strided vector
367: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
368: VecStrideScatterAll()
369: @*/
370: int VecStrideScatterAll(Vec *s,Vec v,InsertMode addv)
371: {
372: int i,n,ierr,bs,ns,j;
373: PetscScalar *x,**y;
378: VecGetLocalSize(v,&n);
379: VecGetLocalSize(*s,&ns);
380: VecGetArray(v,&x);
381: bs = v->bs;
383: PetscMalloc(bs*sizeof(PetscReal*),&y);
384: for (i=0; i<bs; i++) {
385: VecGetArray(s[i],&y[i]);
386: }
388: if (n != ns*bs) {
389: SETERRQ2(1,"Subvector length * blocksize %d not correct for gather from original vector %d",ns*bs,n);
390: }
391: n = n/bs;
393: if (addv == INSERT_VALUES) {
394: for (j=0; j<bs; j++) {
395: for (i=0; i<n; i++) {
396: x[bs*i+j] = y[j][i];
397: }
398: }
399: } else if (addv == ADD_VALUES) {
400: for (j=0; j<bs; j++) {
401: for (i=0; i<n; i++) {
402: x[bs*i+j] += y[j][i];
403: }
404: }
405: #if !defined(PETSC_USE_COMPLEX)
406: } else if (addv == MAX_VALUES) {
407: for (j=0; j<bs; j++) {
408: for (i=0; i<n; i++) {
409: x[bs*i+j] = PetscMax(y[j][i],x[bs*i+j]);
410: }
411: }
412: #endif
413: } else {
414: SETERRQ(1,"Unknown insert type");
415: }
417: VecRestoreArray(v,&x);
418: for (i=0; i<bs; i++) {
419: VecRestoreArray(s[i],&y[i]);
420: }
421: PetscFree(y);
422: return(0);
423: }
425: #undef __FUNCT__
427: /*@
428: VecStrideGather - Gathers a single component from a multi-component vector into
429: another vector.
431: Collective on Vec
433: Input Parameter:
434: + v - the vector
435: . start - starting point of the subvector (defined by a stride)
436: - addv - one of ADD_VALUES,SET_VALUES,MAX_VALUES
438: Output Parameter:
439: . s - the location where the subvector is stored
441: Notes:
442: One must call VecSetBlockSize() before this routine to set the stride
443: information, or use a vector created from a multicomponent DA.
445: If x is the array representing the vector x then this gathers
446: the array (x[start],x[start+stride],x[start+2*stride], ....)
448: The parallel layout of the vector and the subvector must be the same;
449: i.e., nlocal of v = stride*(nlocal of s)
451: Level: advanced
453: Concepts: gather^into strided vector
455: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
456: VecStrideScatterAll()
457: @*/
458: int VecStrideGather(Vec v,int start,Vec s,InsertMode addv)
459: {
460: int i,n,ierr,bs,ns;
461: PetscScalar *x,*y;
466: VecGetLocalSize(v,&n);
467: VecGetLocalSize(s,&ns);
468: VecGetArray(v,&x);
469: VecGetArray(s,&y);
471: bs = v->bs;
472: if (start >= bs) {
473: SETERRQ2(1,"Start of stride subvector (%d) is too large for striden
474: Have you set the vector blocksize (%d) correctly with VecSetBlockSize()?",start,bs);
475: }
476: if (n != ns*bs) {
477: SETERRQ2(1,"Subvector length * blocksize %d not correct for gather from original vector %d",ns*bs,n);
478: }
479: x += start;
480: n = n/bs;
482: if (addv == INSERT_VALUES) {
483: for (i=0; i<n; i++) {
484: y[i] = x[bs*i];
485: }
486: } else if (addv == ADD_VALUES) {
487: for (i=0; i<n; i++) {
488: y[i] += x[bs*i];
489: }
490: #if !defined(PETSC_USE_COMPLEX)
491: } else if (addv == MAX_VALUES) {
492: for (i=0; i<n; i++) {
493: y[i] = PetscMax(y[i],x[bs*i]);
494: }
495: #endif
496: } else {
497: SETERRQ(1,"Unknown insert type");
498: }
500: VecRestoreArray(v,&x);
501: VecRestoreArray(s,&y);
502: return(0);
503: }
505: #undef __FUNCT__
507: /*@
508: VecStrideScatter - Scatters a single component from a vector into a multi-component vector.
510: Collective on Vec
512: Input Parameter:
513: + s - the single-component vector
514: . start - starting point of the subvector (defined by a stride)
515: - addv - one of ADD_VALUES,SET_VALUES,MAX_VALUES
517: Output Parameter:
518: . v - the location where the subvector is scattered (the multi-component vector)
520: Notes:
521: One must call VecSetBlockSize() on the multi-component vector before this
522: routine to set the stride information, or use a vector created from a multicomponent DA.
524: The parallel layout of the vector and the subvector must be the same;
525: i.e., nlocal of v = stride*(nlocal of s)
527: Level: advanced
529: Concepts: scatter^into strided vector
531: .seealso: VecStrideNorm(), VecStrideGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
532: VecStrideScatterAll()
533: @*/
534: int VecStrideScatter(Vec s,int start,Vec v,InsertMode addv)
535: {
536: int i,n,ierr,bs,ns;
537: PetscScalar *x,*y;
542: VecGetLocalSize(v,&n);
543: VecGetLocalSize(s,&ns);
544: VecGetArray(v,&x);
545: VecGetArray(s,&y);
547: bs = v->bs;
548: if (start >= bs) {
549: SETERRQ2(1,"Start of stride subvector (%d) is too large for striden
550: Have you set the vector blocksize (%d) correctly with VecSetBlockSize()?",start,bs);
551: }
552: if (n != ns*bs) {
553: SETERRQ2(1,"Subvector length * blocksize %d not correct for scatter to multicomponent vector %d",ns*bs,n);
554: }
555: x += start;
556: n = n/bs;
559: if (addv == INSERT_VALUES) {
560: for (i=0; i<n; i++) {
561: x[bs*i] = y[i];
562: }
563: } else if (addv == ADD_VALUES) {
564: for (i=0; i<n; i++) {
565: x[bs*i] += y[i];
566: }
567: #if !defined(PETSC_USE_COMPLEX)
568: } else if (addv == MAX_VALUES) {
569: for (i=0; i<n; i++) {
570: x[bs*i] = PetscMax(y[i],x[bs*i]);
571: }
572: #endif
573: } else {
574: SETERRQ(1,"Unknown insert type");
575: }
578: VecRestoreArray(v,&x);
579: VecRestoreArray(s,&y);
580: return(0);
581: }
583: #undef __FUNCT__
585: int VecReciprocal_Default(Vec v)
586: {
587: int i,n,ierr;
588: PetscScalar *x;
591: VecGetLocalSize(v,&n);
592: VecGetArrayFast(v,&x);
593: for (i=0; i<n; i++) {
594: if (x[i] != 0.0) x[i] = 1.0/x[i];
595: }
596: VecRestoreArrayFast(v,&x);
597: return(0);
598: }
600: #undef __FUNCT__
602: /*@
603: VecSqrt - Replaces each component of a vector by the square root of its magnitude.
605: Not collective
607: Input Parameter:
608: . v - The vector
610: Output Parameter:
611: . v - The vector square root
613: Level: beginner
615: Note: The actual function is sqrt(|x_i|)
617: .keywords: vector, sqrt, square root
618: @*/
619: int VecSqrt(Vec v)
620: {
621: PetscScalar *x;
622: int i, n;
623: int ierr;
627: VecGetLocalSize(v, &n);
628: VecGetArray(v, &x);
629: for(i = 0; i < n; i++) {
630: x[i] = sqrt(PetscAbsScalar(x[i]));
631: }
632: VecRestoreArray(v, &x);
633: return(0);
634: }
636: #undef __FUNCT__
638: /*@
639: VecSum - Computes the sum of all the components of a vector.
641: Collective on Vec
643: Input Parameter:
644: . v - the vector
646: Output Parameter:
647: . sum - the result
649: Level: beginner
651: Concepts: sum^of vector entries
653: .seealso: VecNorm()
654: @*/
655: int VecSum(Vec v,PetscScalar *sum)
656: {
657: int i,n,ierr;
658: PetscScalar *x,lsum = 0.0;
662: VecGetLocalSize(v,&n);
663: VecGetArray(v,&x);
664: for (i=0; i<n; i++) {
665: lsum += x[i];
666: }
667: MPI_Allreduce(&lsum,sum,1,MPIU_SCALAR,PetscSum_Op,v->comm);
668: VecRestoreArray(v,&x);
669: return(0);
670: }
672: #undef __FUNCT__
674: /*@
675: VecShift - Shifts all of the components of a vector by computing
676: x[i] = x[i] + shift.
678: Collective on Vec
680: Input Parameters:
681: + v - the vector
682: - shift - the shift
684: Output Parameter:
685: . v - the shifted vector
687: Level: intermediate
689: Concepts: vector^adding constant
691: @*/
692: int VecShift(const PetscScalar *shift,Vec v)
693: {
694: int i,n,ierr;
695: PetscScalar *x,lsum = *shift;
699: VecGetLocalSize(v,&n);
700: VecGetArray(v,&x);
701: for (i=0; i<n; i++) {
702: x[i] += lsum;
703: }
704: VecRestoreArray(v,&x);
705: return(0);
706: }
708: #undef __FUNCT__
710: /*@
711: VecAbs - Replaces every element in a vector with its absolute value.
713: Collective on Vec
715: Input Parameters:
716: . v - the vector
718: Level: intermediate
720: Concepts: vector^absolute value
722: @*/
723: int VecAbs(Vec v)
724: {
725: int i,n,ierr;
726: PetscScalar *x;
730: VecGetLocalSize(v,&n);
731: VecGetArray(v,&x);
732: for (i=0; i<n; i++) {
733: x[i] = PetscAbsScalar(x[i]);
734: }
735: VecRestoreArray(v,&x);
736: return(0);
737: }
739: #undef __FUNCT__
741: /*@
742: VecPermute - Permutes a vector in place using the given ordering.
744: Input Parameters:
745: + vec - The vector
746: . order - The ordering
747: - inv - The flag for inverting the permutation
749: Level: beginner
751: Note: This function does not yet support parallel Index Sets
753: .seealso: MatPermute()
754: .keywords: vec, permute
755: @*/
756: int VecPermute(Vec x, IS row, PetscTruth inv)
757: {
758: PetscScalar *array, *newArray;
759: int *idx;
760: int i;
761: int ierr;
764: ISGetIndices(row, &idx);
765: VecGetArray(x, &array);
766: PetscMalloc((x->n+1) * sizeof(PetscScalar), &newArray);
767: #ifdef PETSC_USE_BOPT_g
768: for(i = 0; i < x->n; i++) {
769: if ((idx[i] < 0) || (idx[i] >= x->n)) {
770: SETERRQ2(PETSC_ERR_ARG_CORRUPT, "Permutation index %d is out of bounds: %d", i, idx[i]);
771: }
772: }
773: #endif
774: if (inv == PETSC_FALSE) {
775: for(i = 0; i < x->n; i++) newArray[i] = array[idx[i]];
776: } else {
777: for(i = 0; i < x->n; i++) newArray[idx[i]] = array[i];
778: }
779: VecRestoreArray(x, &array);
780: ISRestoreIndices(row, &idx);
781: VecReplaceArray(x, newArray);
782: return(0);
783: }
785: #undef __FUNCT__
787: /*@
788: VecEqual - Compares two vectors.
790: Collective on Vec
792: Input Parameters:
793: + vec1 - the first matrix
794: - vec2 - the second matrix
796: Output Parameter:
797: . flg - PETSC_TRUE if the vectors are equal; PETSC_FALSE otherwise.
799: Level: intermediate
801: Concepts: equal^two vectors
802: Concepts: vector^equality
804: @*/
805: int VecEqual(Vec vec1,Vec vec2,PetscTruth *flg)
806: {
807: PetscScalar *v1,*v2;
808: int n1,n2,ierr;
809: PetscTruth flg1;
812: VecGetSize(vec1,&n1);
813: VecGetSize(vec2,&n2);
814: if (vec1 == vec2) {
815: flg1 = PETSC_TRUE;
816: } else if (n1 != n2) {
817: flg1 = PETSC_FALSE;
818: } else {
819: VecGetArray(vec1,&v1);
820: VecGetArray(vec2,&v2);
821: PetscMemcmp(v1,v2,n1*sizeof(PetscScalar),&flg1);
822: VecRestoreArray(vec1,&v1);
823: VecRestoreArray(vec2,&v2);
824: }
826: /* combine results from all processors */
827: MPI_Allreduce(&flg1,flg,1,MPI_INT,MPI_MIN,vec1->comm);
828:
830: return(0);
831: }