Actual source code: vinv.c
petsc-3.4.5 2014-06-29
2: /*
3: Some useful vector utility functions.
4: */
5: #include <../src/vec/vec/impls/mpi/pvecimpl.h> /*I "petscvec.h" I*/
6: extern MPI_Op VecMax_Local_Op;
7: extern MPI_Op VecMin_Local_Op;
11: /*@
12: VecStrideSet - Sets a subvector of a vector defined
13: by a starting point and a stride with a given value
15: Logically Collective on Vec
17: Input Parameter:
18: + v - the vector
19: . start - starting point of the subvector (defined by a stride)
20: - s - value to multiply each subvector entry by
22: Notes:
23: One must call VecSetBlockSize() before this routine to set the stride
24: information, or use a vector created from a multicomponent DMDA.
26: This will only work if the desire subvector is a stride subvector
28: Level: advanced
30: Concepts: scale^on stride of vector
31: Concepts: stride^scale
33: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideScale()
34: @*/
35: PetscErrorCode VecStrideSet(Vec v,PetscInt start,PetscScalar s)
36: {
38: PetscInt i,n,bs;
39: PetscScalar *x;
46: VecGetLocalSize(v,&n);
47: VecGetArray(v,&x);
48: bs = v->map->bs;
49: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
50: 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);
51: x += start;
53: for (i=0; i<n; i+=bs) x[i] = s;
54: x -= start;
56: VecRestoreArray(v,&x);
57: return(0);
58: }
62: /*@
63: VecStrideScale - Scales a subvector of a vector defined
64: by a starting point and a stride.
66: Logically Collective on Vec
68: Input Parameter:
69: + v - the vector
70: . start - starting point of the subvector (defined by a stride)
71: - scale - value to multiply each subvector entry by
73: Notes:
74: One must call VecSetBlockSize() before this routine to set the stride
75: information, or use a vector created from a multicomponent DMDA.
77: This will only work if the desire subvector is a stride subvector
79: Level: advanced
81: Concepts: scale^on stride of vector
82: Concepts: stride^scale
84: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideScale()
85: @*/
86: PetscErrorCode VecStrideScale(Vec v,PetscInt start,PetscScalar scale)
87: {
89: PetscInt i,n,bs;
90: PetscScalar *x;
97: VecGetLocalSize(v,&n);
98: VecGetArray(v,&x);
99: bs = v->map->bs;
100: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
101: 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);
102: x += start;
104: for (i=0; i<n; i+=bs) x[i] *= scale;
105: x -= start;
107: VecRestoreArray(v,&x);
108: return(0);
109: }
113: /*@
114: VecStrideNorm - Computes the norm of subvector of a vector defined
115: by a starting point and a stride.
117: Collective on Vec
119: Input Parameter:
120: + v - the vector
121: . start - starting point of the subvector (defined by a stride)
122: - ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY
124: Output Parameter:
125: . norm - the norm
127: Notes:
128: One must call VecSetBlockSize() before this routine to set the stride
129: information, or use a vector created from a multicomponent DMDA.
131: If x is the array representing the vector x then this computes the norm
132: of the array (x[start],x[start+stride],x[start+2*stride], ....)
134: This is useful for computing, say the norm of the pressure variable when
135: the pressure is stored (interlaced) with other variables, say density etc.
137: This will only work if the desire subvector is a stride subvector
139: Level: advanced
141: Concepts: norm^on stride of vector
142: Concepts: stride^norm
144: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
145: @*/
146: PetscErrorCode VecStrideNorm(Vec v,PetscInt start,NormType ntype,PetscReal *nrm)
147: {
149: PetscInt i,n,bs;
150: PetscScalar *x;
151: PetscReal tnorm;
152: MPI_Comm comm;
157: VecGetLocalSize(v,&n);
158: VecGetArray(v,&x);
159: PetscObjectGetComm((PetscObject)v,&comm);
161: bs = v->map->bs;
162: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
163: 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);
164: x += start;
166: if (ntype == NORM_2) {
167: PetscScalar sum = 0.0;
168: for (i=0; i<n; i+=bs) sum += x[i]*(PetscConj(x[i]));
169: tnorm = PetscRealPart(sum);
170: MPI_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_SUM,comm);
171: *nrm = PetscSqrtReal(*nrm);
172: } else if (ntype == NORM_1) {
173: tnorm = 0.0;
174: for (i=0; i<n; i+=bs) tnorm += PetscAbsScalar(x[i]);
175: MPI_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_SUM,comm);
176: } else if (ntype == NORM_INFINITY) {
177: PetscReal tmp;
178: tnorm = 0.0;
180: for (i=0; i<n; i+=bs) {
181: if ((tmp = PetscAbsScalar(x[i])) > tnorm) tnorm = tmp;
182: /* check special case of tmp == NaN */
183: if (tmp != tmp) {tnorm = tmp; break;}
184: }
185: MPI_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_MAX,comm);
186: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
187: VecRestoreArray(v,&x);
188: return(0);
189: }
193: /*@
194: VecStrideMax - Computes the maximum of subvector of a vector defined
195: by a starting point and a stride and optionally its location.
197: Collective on Vec
199: Input Parameter:
200: + v - the vector
201: - start - starting point of the subvector (defined by a stride)
203: Output Parameter:
204: + index - the location where the maximum occurred (pass NULL if not required)
205: - nrm - the max
207: Notes:
208: One must call VecSetBlockSize() before this routine to set the stride
209: information, or use a vector created from a multicomponent DMDA.
211: If xa is the array representing the vector x, then this computes the max
212: of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)
214: This is useful for computing, say the maximum of the pressure variable when
215: the pressure is stored (interlaced) with other variables, e.g., density, etc.
216: This will only work if the desire subvector is a stride subvector.
218: Level: advanced
220: Concepts: maximum^on stride of vector
221: Concepts: stride^maximum
223: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
224: @*/
225: PetscErrorCode VecStrideMax(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
226: {
228: PetscInt i,n,bs,id;
229: PetscScalar *x;
230: PetscReal max,tmp;
231: MPI_Comm comm;
237: VecGetLocalSize(v,&n);
238: VecGetArray(v,&x);
239: PetscObjectGetComm((PetscObject)v,&comm);
241: bs = v->map->bs;
242: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
243: 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);
244: x += start;
246: id = -1;
247: if (!n) max = PETSC_MIN_REAL;
248: else {
249: id = 0;
250: max = PetscRealPart(x[0]);
251: for (i=bs; i<n; i+=bs) {
252: if ((tmp = PetscRealPart(x[i])) > max) { max = tmp; id = i;}
253: }
254: }
255: VecRestoreArray(v,&x);
257: if (!idex) {
258: MPI_Allreduce(&max,nrm,1,MPIU_REAL,MPIU_MAX,comm);
259: } else {
260: PetscReal in[2],out[2];
261: PetscInt rstart;
263: VecGetOwnershipRange(v,&rstart,NULL);
264: in[0] = max;
265: in[1] = rstart+id+start;
266: MPI_Allreduce(in,out,2,MPIU_REAL,VecMax_Local_Op,PetscObjectComm((PetscObject)v));
267: *nrm = out[0];
268: *idex = (PetscInt)out[1];
269: }
270: return(0);
271: }
275: /*@
276: VecStrideMin - Computes the minimum of subvector of a vector defined
277: by a starting point and a stride and optionally its location.
279: Collective on Vec
281: Input Parameter:
282: + v - the vector
283: - start - starting point of the subvector (defined by a stride)
285: Output Parameter:
286: + idex - the location where the minimum occurred. (pass NULL if not required)
287: - nrm - the min
289: Level: advanced
291: Notes:
292: One must call VecSetBlockSize() before this routine to set the stride
293: information, or use a vector created from a multicomponent DMDA.
295: If xa is the array representing the vector x, then this computes the min
296: of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)
298: This is useful for computing, say the minimum of the pressure variable when
299: the pressure is stored (interlaced) with other variables, e.g., density, etc.
300: This will only work if the desire subvector is a stride subvector.
302: Concepts: minimum^on stride of vector
303: Concepts: stride^minimum
305: .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
306: @*/
307: PetscErrorCode VecStrideMin(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
308: {
310: PetscInt i,n,bs,id;
311: PetscScalar *x;
312: PetscReal min,tmp;
313: MPI_Comm comm;
319: VecGetLocalSize(v,&n);
320: VecGetArray(v,&x);
321: PetscObjectGetComm((PetscObject)v,&comm);
323: bs = v->map->bs;
324: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
325: 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);
326: x += start;
328: id = -1;
329: if (!n) min = PETSC_MAX_REAL;
330: else {
331: id = 0;
332: min = PetscRealPart(x[0]);
333: for (i=bs; i<n; i+=bs) {
334: if ((tmp = PetscRealPart(x[i])) < min) { min = tmp; id = i;}
335: }
336: }
337: VecRestoreArray(v,&x);
339: if (!idex) {
340: MPI_Allreduce(&min,nrm,1,MPIU_REAL,MPIU_MIN,comm);
341: } else {
342: PetscReal in[2],out[2];
343: PetscInt rstart;
345: VecGetOwnershipRange(v,&rstart,NULL);
346: in[0] = min;
347: in[1] = rstart+id;
348: MPI_Allreduce(in,out,2,MPIU_REAL,VecMin_Local_Op,PetscObjectComm((PetscObject)v));
349: *nrm = out[0];
350: *idex = (PetscInt)out[1];
351: }
352: return(0);
353: }
357: /*@
358: VecStrideScaleAll - Scales the subvectors of a vector defined
359: by a starting point and a stride.
361: Logically Collective on Vec
363: Input Parameter:
364: + v - the vector
365: - scales - values to multiply each subvector entry by
367: Notes:
368: One must call VecSetBlockSize() before this routine to set the stride
369: information, or use a vector created from a multicomponent DMDA.
372: Level: advanced
374: Concepts: scale^on stride of vector
375: Concepts: stride^scale
377: .seealso: VecNorm(), VecStrideScale(), VecScale(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
378: @*/
379: PetscErrorCode VecStrideScaleAll(Vec v,const PetscScalar *scales)
380: {
382: PetscInt i,j,n,bs;
383: PetscScalar *x;
388: VecGetLocalSize(v,&n);
389: VecGetArray(v,&x);
391: bs = v->map->bs;
393: /* need to provide optimized code for each bs */
394: for (i=0; i<n; i+=bs) {
395: for (j=0; j<bs; j++) x[i+j] *= scales[j];
396: }
397: VecRestoreArray(v,&x);
398: return(0);
399: }
403: /*@
404: VecStrideNormAll - Computes the norms subvectors of a vector defined
405: by a starting point and a stride.
407: Collective on Vec
409: Input Parameter:
410: + v - the vector
411: - ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY
413: Output Parameter:
414: . nrm - the norms
416: Notes:
417: One must call VecSetBlockSize() before this routine to set the stride
418: information, or use a vector created from a multicomponent DMDA.
420: If x is the array representing the vector x then this computes the norm
421: of the array (x[start],x[start+stride],x[start+2*stride], ....)
423: This is useful for computing, say the norm of the pressure variable when
424: the pressure is stored (interlaced) with other variables, say density etc.
426: This will only work if the desire subvector is a stride subvector
428: Level: advanced
430: Concepts: norm^on stride of vector
431: Concepts: stride^norm
433: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
434: @*/
435: PetscErrorCode VecStrideNormAll(Vec v,NormType ntype,PetscReal nrm[])
436: {
438: PetscInt i,j,n,bs;
439: PetscScalar *x;
440: PetscReal tnorm[128];
441: MPI_Comm comm;
446: VecGetLocalSize(v,&n);
447: VecGetArray(v,&x);
448: PetscObjectGetComm((PetscObject)v,&comm);
450: bs = v->map->bs;
451: if (bs > 128) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");
453: if (ntype == NORM_2) {
454: PetscScalar sum[128];
455: for (j=0; j<bs; j++) sum[j] = 0.0;
456: for (i=0; i<n; i+=bs) {
457: for (j=0; j<bs; j++) sum[j] += x[i+j]*(PetscConj(x[i+j]));
458: }
459: for (j=0; j<bs; j++) tnorm[j] = PetscRealPart(sum[j]);
461: MPI_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_SUM,comm);
462: for (j=0; j<bs; j++) nrm[j] = PetscSqrtReal(nrm[j]);
463: } else if (ntype == NORM_1) {
464: for (j=0; j<bs; j++) tnorm[j] = 0.0;
466: for (i=0; i<n; i+=bs) {
467: for (j=0; j<bs; j++) tnorm[j] += PetscAbsScalar(x[i+j]);
468: }
470: MPI_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_SUM,comm);
471: } else if (ntype == NORM_INFINITY) {
472: PetscReal tmp;
473: for (j=0; j<bs; j++) tnorm[j] = 0.0;
475: for (i=0; i<n; i+=bs) {
476: for (j=0; j<bs; j++) {
477: if ((tmp = PetscAbsScalar(x[i+j])) > tnorm[j]) tnorm[j] = tmp;
478: /* check special case of tmp == NaN */
479: if (tmp != tmp) {tnorm[j] = tmp; break;}
480: }
481: }
482: MPI_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_MAX,comm);
483: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
484: VecRestoreArray(v,&x);
485: return(0);
486: }
490: /*@
491: VecStrideMaxAll - Computes the maximums of subvectors of a vector defined
492: by a starting point and a stride and optionally its location.
494: Collective on Vec
496: Input Parameter:
497: . v - the vector
499: Output Parameter:
500: + index - the location where the maximum occurred (not supported, pass NULL,
501: if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
502: - nrm - the maximums
504: Notes:
505: One must call VecSetBlockSize() before this routine to set the stride
506: information, or use a vector created from a multicomponent DMDA.
508: This is useful for computing, say the maximum of the pressure variable when
509: the pressure is stored (interlaced) with other variables, e.g., density, etc.
510: This will only work if the desire subvector is a stride subvector.
512: Level: advanced
514: Concepts: maximum^on stride of vector
515: Concepts: stride^maximum
517: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
518: @*/
519: PetscErrorCode VecStrideMaxAll(Vec v,PetscInt idex[],PetscReal nrm[])
520: {
522: PetscInt i,j,n,bs;
523: PetscScalar *x;
524: PetscReal max[128],tmp;
525: MPI_Comm comm;
530: 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");
531: VecGetLocalSize(v,&n);
532: VecGetArray(v,&x);
533: PetscObjectGetComm((PetscObject)v,&comm);
535: bs = v->map->bs;
536: if (bs > 128) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");
538: if (!n) {
539: for (j=0; j<bs; j++) max[j] = PETSC_MIN_REAL;
540: } else {
541: for (j=0; j<bs; j++) max[j] = PetscRealPart(x[j]);
543: for (i=bs; i<n; i+=bs) {
544: for (j=0; j<bs; j++) {
545: if ((tmp = PetscRealPart(x[i+j])) > max[j]) max[j] = tmp;
546: }
547: }
548: }
549: MPI_Allreduce(max,nrm,bs,MPIU_REAL,MPIU_MAX,comm);
551: VecRestoreArray(v,&x);
552: return(0);
553: }
557: /*@
558: VecStrideMinAll - Computes the minimum of subvector of a vector defined
559: by a starting point and a stride and optionally its location.
561: Collective on Vec
563: Input Parameter:
564: . v - the vector
566: Output Parameter:
567: + idex - the location where the minimum occurred (not supported, pass NULL,
568: if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
569: - nrm - the minimums
571: Level: advanced
573: Notes:
574: One must call VecSetBlockSize() before this routine to set the stride
575: information, or use a vector created from a multicomponent DMDA.
577: This is useful for computing, say the minimum of the pressure variable when
578: the pressure is stored (interlaced) with other variables, e.g., density, etc.
579: This will only work if the desire subvector is a stride subvector.
581: Concepts: minimum^on stride of vector
582: Concepts: stride^minimum
584: .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
585: @*/
586: PetscErrorCode VecStrideMinAll(Vec v,PetscInt idex[],PetscReal nrm[])
587: {
589: PetscInt i,n,bs,j;
590: PetscScalar *x;
591: PetscReal min[128],tmp;
592: MPI_Comm comm;
597: 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");
598: VecGetLocalSize(v,&n);
599: VecGetArray(v,&x);
600: PetscObjectGetComm((PetscObject)v,&comm);
602: bs = v->map->bs;
603: if (bs > 128) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");
605: if (!n) {
606: for (j=0; j<bs; j++) min[j] = PETSC_MAX_REAL;
607: } else {
608: for (j=0; j<bs; j++) min[j] = PetscRealPart(x[j]);
610: for (i=bs; i<n; i+=bs) {
611: for (j=0; j<bs; j++) {
612: if ((tmp = PetscRealPart(x[i+j])) < min[j]) min[j] = tmp;
613: }
614: }
615: }
616: MPI_Allreduce(min,nrm,bs,MPIU_REAL,MPIU_MIN,comm);
618: VecRestoreArray(v,&x);
619: return(0);
620: }
622: /*----------------------------------------------------------------------------------------------*/
625: /*@
626: VecStrideGatherAll - Gathers all the single components from a multi-component vector into
627: separate vectors.
629: Collective on Vec
631: Input Parameter:
632: + v - the vector
633: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
635: Output Parameter:
636: . s - the location where the subvectors are stored
638: Notes:
639: One must call VecSetBlockSize() before this routine to set the stride
640: information, or use a vector created from a multicomponent DMDA.
642: If x is the array representing the vector x then this gathers
643: the arrays (x[start],x[start+stride],x[start+2*stride], ....)
644: for start=0,1,2,...bs-1
646: The parallel layout of the vector and the subvector must be the same;
647: i.e., nlocal of v = stride*(nlocal of s)
649: Not optimized; could be easily
651: Level: advanced
653: Concepts: gather^into strided vector
655: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
656: VecStrideScatterAll()
657: @*/
658: PetscErrorCode VecStrideGatherAll(Vec v,Vec s[],InsertMode addv)
659: {
661: PetscInt i,n,n2,bs,j,k,*bss = NULL,nv,jj,nvc;
662: PetscScalar *x,**y;
668: VecGetLocalSize(v,&n);
669: VecGetLocalSize(s[0],&n2);
670: VecGetArray(v,&x);
671: bs = v->map->bs;
672: if (bs < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Input vector does not have a valid blocksize set");
674: PetscMalloc2(bs,PetscReal*,&y,bs,PetscInt,&bss);
675: nv = 0;
676: nvc = 0;
677: for (i=0; i<bs; i++) {
678: VecGetBlockSize(s[i],&bss[i]);
679: if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1 Re: [PETSC #8241] VecStrideGatherAll */
680: VecGetArray(s[i],&y[i]);
681: nvc += bss[i];
682: nv++;
683: if (nvc > bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of subvectors in subvectors > number of vectors in main vector");
684: if (nvc == bs) break;
685: }
687: n = n/bs;
689: jj = 0;
690: if (addv == INSERT_VALUES) {
691: for (j=0; j<nv; j++) {
692: for (k=0; k<bss[j]; k++) {
693: for (i=0; i<n; i++) y[j][i*bss[j] + k] = x[bs*i+jj+k];
694: }
695: jj += bss[j];
696: }
697: } else if (addv == ADD_VALUES) {
698: for (j=0; j<nv; j++) {
699: for (k=0; k<bss[j]; k++) {
700: for (i=0; i<n; i++) y[j][i*bss[j] + k] += x[bs*i+jj+k];
701: }
702: jj += bss[j];
703: }
704: #if !defined(PETSC_USE_COMPLEX)
705: } else if (addv == MAX_VALUES) {
706: for (j=0; j<nv; j++) {
707: for (k=0; k<bss[j]; k++) {
708: for (i=0; i<n; i++) y[j][i*bss[j] + k] = PetscMax(y[j][i*bss[j] + k],x[bs*i+jj+k]);
709: }
710: jj += bss[j];
711: }
712: #endif
713: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
715: VecRestoreArray(v,&x);
716: for (i=0; i<nv; i++) {
717: VecRestoreArray(s[i],&y[i]);
718: }
720: PetscFree2(y,bss);
721: return(0);
722: }
726: /*@
727: VecStrideScatterAll - Scatters all the single components from separate vectors into
728: a multi-component vector.
730: Collective on Vec
732: Input Parameter:
733: + s - the location where the subvectors are stored
734: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
736: Output Parameter:
737: . v - the multicomponent vector
739: Notes:
740: One must call VecSetBlockSize() before this routine to set the stride
741: information, or use a vector created from a multicomponent DMDA.
743: The parallel layout of the vector and the subvector must be the same;
744: i.e., nlocal of v = stride*(nlocal of s)
746: Not optimized; could be easily
748: Level: advanced
750: Concepts: scatter^into strided vector
752: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
753: VecStrideScatterAll()
754: @*/
755: PetscErrorCode VecStrideScatterAll(Vec s[],Vec v,InsertMode addv)
756: {
758: PetscInt i,n,n2,bs,j,jj,k,*bss = NULL,nv,nvc;
759: PetscScalar *x,**y;
765: VecGetLocalSize(v,&n);
766: VecGetLocalSize(s[0],&n2);
767: VecGetArray(v,&x);
768: bs = v->map->bs;
769: if (bs < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Input vector does not have a valid blocksize set");
771: PetscMalloc2(bs,PetscScalar**,&y,bs,PetscInt,&bss);
772: nv = 0;
773: nvc = 0;
774: for (i=0; i<bs; i++) {
775: VecGetBlockSize(s[i],&bss[i]);
776: if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1 Re: [PETSC #8241] VecStrideGatherAll */
777: VecGetArray(s[i],&y[i]);
778: nvc += bss[i];
779: nv++;
780: if (nvc > bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of subvectors in subvectors > number of vectors in main vector");
781: if (nvc == bs) break;
782: }
784: n = n/bs;
786: jj = 0;
787: if (addv == INSERT_VALUES) {
788: for (j=0; j<nv; j++) {
789: for (k=0; k<bss[j]; k++) {
790: for (i=0; i<n; i++) x[bs*i+jj+k] = y[j][i*bss[j] + k];
791: }
792: jj += bss[j];
793: }
794: } else if (addv == ADD_VALUES) {
795: for (j=0; j<nv; j++) {
796: for (k=0; k<bss[j]; k++) {
797: for (i=0; i<n; i++) x[bs*i+jj+k] += y[j][i*bss[j] + k];
798: }
799: jj += bss[j];
800: }
801: #if !defined(PETSC_USE_COMPLEX)
802: } else if (addv == MAX_VALUES) {
803: for (j=0; j<nv; j++) {
804: for (k=0; k<bss[j]; k++) {
805: for (i=0; i<n; i++) x[bs*i+jj+k] = PetscMax(x[bs*i+jj+k],y[j][i*bss[j] + k]);
806: }
807: jj += bss[j];
808: }
809: #endif
810: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
812: VecRestoreArray(v,&x);
813: for (i=0; i<nv; i++) {
814: VecRestoreArray(s[i],&y[i]);
815: }
816: PetscFree2(y,bss);
817: return(0);
818: }
822: /*@
823: VecStrideGather - Gathers a single component from a multi-component vector into
824: another vector.
826: Collective on Vec
828: Input Parameter:
829: + v - the vector
830: . start - starting point of the subvector (defined by a stride)
831: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
833: Output Parameter:
834: . s - the location where the subvector is stored
836: Notes:
837: One must call VecSetBlockSize() before this routine to set the stride
838: information, or use a vector created from a multicomponent DMDA.
840: If x is the array representing the vector x then this gathers
841: the array (x[start],x[start+stride],x[start+2*stride], ....)
843: The parallel layout of the vector and the subvector must be the same;
844: i.e., nlocal of v = stride*(nlocal of s)
846: Not optimized; could be easily
848: Level: advanced
850: Concepts: gather^into strided vector
852: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
853: VecStrideScatterAll()
854: @*/
855: PetscErrorCode VecStrideGather(Vec v,PetscInt start,Vec s,InsertMode addv)
856: {
862: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
863: 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);
864: if (!v->ops->stridegather) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
865: (*v->ops->stridegather)(v,start,s,addv);
866: return(0);
867: }
871: /*@
872: VecStrideScatter - Scatters a single component from a vector into a multi-component vector.
874: Collective on Vec
876: Input Parameter:
877: + s - the single-component vector
878: . start - starting point of the subvector (defined by a stride)
879: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
881: Output Parameter:
882: . v - the location where the subvector is scattered (the multi-component vector)
884: Notes:
885: One must call VecSetBlockSize() on the multi-component vector before this
886: routine to set the stride information, or use a vector created from a multicomponent DMDA.
888: The parallel layout of the vector and the subvector must be the same;
889: i.e., nlocal of v = stride*(nlocal of s)
891: Not optimized; could be easily
893: Level: advanced
895: Concepts: scatter^into strided vector
897: .seealso: VecStrideNorm(), VecStrideGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
898: VecStrideScatterAll()
899: @*/
900: PetscErrorCode VecStrideScatter(Vec s,PetscInt start,Vec v,InsertMode addv)
901: {
907: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
908: 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);
909: if (!v->ops->stridescatter) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
910: (*v->ops->stridescatter)(s,start,v,addv);
911: return(0);
912: }
916: PetscErrorCode VecStrideGather_Default(Vec v,PetscInt start,Vec s,InsertMode addv)
917: {
919: PetscInt i,n,bs,ns;
920: PetscScalar *x,*y;
923: VecGetLocalSize(v,&n);
924: VecGetLocalSize(s,&ns);
925: VecGetArray(v,&x);
926: VecGetArray(s,&y);
928: bs = v->map->bs;
929: 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);
930: x += start;
931: n = n/bs;
933: if (addv == INSERT_VALUES) {
934: for (i=0; i<n; i++) y[i] = x[bs*i];
935: } else if (addv == ADD_VALUES) {
936: for (i=0; i<n; i++) y[i] += x[bs*i];
937: #if !defined(PETSC_USE_COMPLEX)
938: } else if (addv == MAX_VALUES) {
939: for (i=0; i<n; i++) y[i] = PetscMax(y[i],x[bs*i]);
940: #endif
941: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
943: VecRestoreArray(v,&x);
944: VecRestoreArray(s,&y);
945: return(0);
946: }
950: PetscErrorCode VecStrideScatter_Default(Vec s,PetscInt start,Vec v,InsertMode addv)
951: {
953: PetscInt i,n,bs,ns;
954: PetscScalar *x,*y;
957: VecGetLocalSize(v,&n);
958: VecGetLocalSize(s,&ns);
959: VecGetArray(v,&x);
960: VecGetArray(s,&y);
962: bs = v->map->bs;
963: 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);
964: x += start;
965: n = n/bs;
967: if (addv == INSERT_VALUES) {
968: for (i=0; i<n; i++) x[bs*i] = y[i];
969: } else if (addv == ADD_VALUES) {
970: for (i=0; i<n; i++) x[bs*i] += y[i];
971: #if !defined(PETSC_USE_COMPLEX)
972: } else if (addv == MAX_VALUES) {
973: for (i=0; i<n; i++) x[bs*i] = PetscMax(y[i],x[bs*i]);
974: #endif
975: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
977: VecRestoreArray(v,&x);
978: VecRestoreArray(s,&y);
979: return(0);
980: }
984: PetscErrorCode VecReciprocal_Default(Vec v)
985: {
987: PetscInt i,n;
988: PetscScalar *x;
991: VecGetLocalSize(v,&n);
992: VecGetArray(v,&x);
993: for (i=0; i<n; i++) {
994: if (x[i] != (PetscScalar)0.0) x[i] = (PetscScalar)1.0/x[i];
995: }
996: VecRestoreArray(v,&x);
997: return(0);
998: }
1002: /*@
1003: VecExp - Replaces each component of a vector by e^x_i
1005: Not collective
1007: Input Parameter:
1008: . v - The vector
1010: Output Parameter:
1011: . v - The vector of exponents
1013: Level: beginner
1015: .seealso: VecLog(), VecAbs(), VecSqrtAbs(), VecReciprocal()
1017: .keywords: vector, sqrt, square root
1018: @*/
1019: PetscErrorCode VecExp(Vec v)
1020: {
1021: PetscScalar *x;
1022: PetscInt i, n;
1027: if (v->ops->exp) {
1028: (*v->ops->exp)(v);
1029: } else {
1030: VecGetLocalSize(v, &n);
1031: VecGetArray(v, &x);
1032: for (i = 0; i < n; i++) x[i] = PetscExpScalar(x[i]);
1033: VecRestoreArray(v, &x);
1034: }
1035: return(0);
1036: }
1040: /*@
1041: VecLog - Replaces each component of a vector by log(x_i), the natural logarithm
1043: Not collective
1045: Input Parameter:
1046: . v - The vector
1048: Output Parameter:
1049: . v - The vector of logs
1051: Level: beginner
1053: .seealso: VecExp(), VecAbs(), VecSqrtAbs(), VecReciprocal()
1055: .keywords: vector, sqrt, square root
1056: @*/
1057: PetscErrorCode VecLog(Vec v)
1058: {
1059: PetscScalar *x;
1060: PetscInt i, n;
1065: if (v->ops->log) {
1066: (*v->ops->log)(v);
1067: } else {
1068: VecGetLocalSize(v, &n);
1069: VecGetArray(v, &x);
1070: for (i = 0; i < n; i++) x[i] = PetscLogScalar(x[i]);
1071: VecRestoreArray(v, &x);
1072: }
1073: return(0);
1074: }
1078: /*@
1079: VecSqrtAbs - Replaces each component of a vector by the square root of its magnitude.
1081: Not collective
1083: Input Parameter:
1084: . v - The vector
1086: Output Parameter:
1087: . v - The vector square root
1089: Level: beginner
1091: Note: The actual function is sqrt(|x_i|)
1093: .seealso: VecLog(), VecExp(), VecReciprocal(), VecAbs()
1095: .keywords: vector, sqrt, square root
1096: @*/
1097: PetscErrorCode VecSqrtAbs(Vec v)
1098: {
1099: PetscScalar *x;
1100: PetscInt i, n;
1105: if (v->ops->sqrt) {
1106: (*v->ops->sqrt)(v);
1107: } else {
1108: VecGetLocalSize(v, &n);
1109: VecGetArray(v, &x);
1110: for (i = 0; i < n; i++) x[i] = PetscSqrtReal(PetscAbsScalar(x[i]));
1111: VecRestoreArray(v, &x);
1112: }
1113: return(0);
1114: }
1118: /*@
1119: VecDotNorm2 - computes the inner product of two vectors and the 2-norm squared of the second vector
1121: Collective on Vec
1123: Input Parameter:
1124: + s - first vector
1125: - t - second vector
1127: Output Parameter:
1128: + dp - s'conj(t)
1129: - nm - t'conj(t)
1131: Level: advanced
1133: Notes: conj(x) is the complex conjugate of x when x is complex
1136: .seealso: VecDot(), VecNorm(), VecDotBegin(), VecNormBegin(), VecDotEnd(), VecNormEnd()
1138: .keywords: vector, sqrt, square root
1139: @*/
1140: PetscErrorCode VecDotNorm2(Vec s,Vec t,PetscScalar *dp, PetscReal *nm)
1141: {
1142: PetscScalar *sx, *tx, dpx = 0.0, nmx = 0.0,work[2],sum[2];
1143: PetscInt i, n;
1154: if (s->map->N != t->map->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
1155: if (s->map->n != t->map->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");
1157: PetscLogEventBarrierBegin(VEC_DotNormBarrier,s,t,0,0,PetscObjectComm((PetscObject)s));
1158: if (s->ops->dotnorm2) {
1159: (*s->ops->dotnorm2)(s,t,dp,&dpx);
1160: *nm = PetscRealPart(dpx);
1161: } else {
1162: VecGetLocalSize(s, &n);
1163: VecGetArray(s, &sx);
1164: VecGetArray(t, &tx);
1166: for (i = 0; i<n; i++) {
1167: dpx += sx[i]*PetscConj(tx[i]);
1168: nmx += tx[i]*PetscConj(tx[i]);
1169: }
1170: work[0] = dpx;
1171: work[1] = nmx;
1173: MPI_Allreduce(work,sum,2,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)s));
1174: *dp = sum[0];
1175: *nm = PetscRealPart(sum[1]);
1177: VecRestoreArray(t, &tx);
1178: VecRestoreArray(s, &sx);
1179: PetscLogFlops(4.0*n);
1180: }
1181: PetscLogEventBarrierEnd(VEC_DotNormBarrier,s,t,0,0,PetscObjectComm((PetscObject)s));
1182: return(0);
1183: }
1187: /*@
1188: VecSum - Computes the sum of all the components of a vector.
1190: Collective on Vec
1192: Input Parameter:
1193: . v - the vector
1195: Output Parameter:
1196: . sum - the result
1198: Level: beginner
1200: Concepts: sum^of vector entries
1202: .seealso: VecNorm()
1203: @*/
1204: PetscErrorCode VecSum(Vec v,PetscScalar *sum)
1205: {
1207: PetscInt i,n;
1208: PetscScalar *x,lsum = 0.0;
1213: VecGetLocalSize(v,&n);
1214: VecGetArray(v,&x);
1215: for (i=0; i<n; i++) lsum += x[i];
1216: MPI_Allreduce(&lsum,sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)v));
1217: VecRestoreArray(v,&x);
1218: return(0);
1219: }
1223: /*@
1224: VecShift - Shifts all of the components of a vector by computing
1225: x[i] = x[i] + shift.
1227: Logically Collective on Vec
1229: Input Parameters:
1230: + v - the vector
1231: - shift - the shift
1233: Output Parameter:
1234: . v - the shifted vector
1236: Level: intermediate
1238: Concepts: vector^adding constant
1240: @*/
1241: PetscErrorCode VecShift(Vec v,PetscScalar shift)
1242: {
1244: PetscInt i,n;
1245: PetscScalar *x;
1250: if (v->ops->shift) {
1251: (*v->ops->shift)(v);
1252: } else {
1253: VecGetLocalSize(v,&n);
1254: VecGetArray(v,&x);
1255: for (i=0; i<n; i++) x[i] += shift;
1256: VecRestoreArray(v,&x);
1257: }
1258: return(0);
1259: }
1263: /*@
1264: VecAbs - Replaces every element in a vector with its absolute value.
1266: Logically Collective on Vec
1268: Input Parameters:
1269: . v - the vector
1271: Level: intermediate
1273: Concepts: vector^absolute value
1275: @*/
1276: PetscErrorCode VecAbs(Vec v)
1277: {
1279: PetscInt i,n;
1280: PetscScalar *x;
1284: if (v->ops->abs) {
1285: (*v->ops->abs)(v);
1286: } else {
1287: VecGetLocalSize(v,&n);
1288: VecGetArray(v,&x);
1289: for (i=0; i<n; i++) x[i] = PetscAbsScalar(x[i]);
1290: VecRestoreArray(v,&x);
1291: }
1292: return(0);
1293: }
1297: /*@
1298: VecPermute - Permutes a vector in place using the given ordering.
1300: Input Parameters:
1301: + vec - The vector
1302: . order - The ordering
1303: - inv - The flag for inverting the permutation
1305: Level: beginner
1307: Note: This function does not yet support parallel Index Sets with non-local permutations
1309: .seealso: MatPermute()
1310: .keywords: vec, permute
1311: @*/
1312: PetscErrorCode VecPermute(Vec x, IS row, PetscBool inv)
1313: {
1314: PetscScalar *array, *newArray;
1315: const PetscInt *idx;
1316: PetscInt i,rstart,rend;
1320: VecGetOwnershipRange(x,&rstart,&rend);
1321: ISGetIndices(row, &idx);
1322: VecGetArray(x, &array);
1323: PetscMalloc(x->map->n*sizeof(PetscScalar), &newArray);
1324: #if defined(PETSC_USE_DEBUG)
1325: for (i = 0; i < x->map->n; i++) {
1326: if ((idx[i] < rstart) || (idx[i] >= rend)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT, "Permutation index %D is out of bounds: %D", i, idx[i]);
1327: }
1328: #endif
1329: if (!inv) {
1330: for (i = 0; i < x->map->n; i++) newArray[i] = array[idx[i]-rstart];
1331: } else {
1332: for (i = 0; i < x->map->n; i++) newArray[idx[i]-rstart] = array[i];
1333: }
1334: VecRestoreArray(x, &array);
1335: ISRestoreIndices(row, &idx);
1336: VecReplaceArray(x, newArray);
1337: return(0);
1338: }
1342: /*@
1343: VecEqual - Compares two vectors. Returns true if the two vectors are either pointing to the same memory buffer,
1344: or if the two vectors have the same local and global layout as well as bitwise equality of all entries.
1345: Does NOT take round-off errors into account.
1347: Collective on Vec
1349: Input Parameters:
1350: + vec1 - the first vector
1351: - vec2 - the second vector
1353: Output Parameter:
1354: . flg - PETSC_TRUE if the vectors are equal; PETSC_FALSE otherwise.
1356: Level: intermediate
1358: Concepts: equal^two vectors
1359: Concepts: vector^equality
1361: @*/
1362: PetscErrorCode VecEqual(Vec vec1,Vec vec2,PetscBool *flg)
1363: {
1364: PetscScalar *v1,*v2;
1366: PetscInt n1,n2,N1,N2;
1367: PetscInt state1,state2;
1368: PetscBool flg1;
1374: if (vec1 == vec2) *flg = PETSC_TRUE;
1375: else {
1376: VecGetSize(vec1,&N1);
1377: VecGetSize(vec2,&N2);
1378: if (N1 != N2) flg1 = PETSC_FALSE;
1379: else {
1380: VecGetLocalSize(vec1,&n1);
1381: VecGetLocalSize(vec2,&n2);
1382: if (n1 != n2) flg1 = PETSC_FALSE;
1383: else {
1384: PetscObjectStateQuery((PetscObject) vec1,&state1);
1385: PetscObjectStateQuery((PetscObject) vec2,&state2);
1386: VecGetArray(vec1,&v1);
1387: VecGetArray(vec2,&v2);
1388: #if defined(PETSC_USE_COMPLEX)
1389: {
1390: PetscInt k;
1391: flg1 = PETSC_TRUE;
1392: for (k=0; k<n1; k++) {
1393: if (PetscRealPart(v1[k]) != PetscRealPart(v2[k]) || PetscImaginaryPart(v1[k]) != PetscImaginaryPart(v2[k])) {
1394: flg1 = PETSC_FALSE;
1395: break;
1396: }
1397: }
1398: }
1399: #else
1400: PetscMemcmp(v1,v2,n1*sizeof(PetscScalar),&flg1);
1401: #endif
1402: VecRestoreArray(vec1,&v1);
1403: VecRestoreArray(vec2,&v2);
1404: PetscObjectSetState((PetscObject) vec1,state1);
1405: PetscObjectSetState((PetscObject) vec2,state2);
1406: }
1407: }
1408: /* combine results from all processors */
1409: MPI_Allreduce(&flg1,flg,1,MPIU_BOOL,MPI_MIN,PetscObjectComm((PetscObject)vec1));
1410: }
1411: return(0);
1412: }