Actual source code: vinv.c
2: /*
3: Some useful vector utility functions.
4: */
5: #include <../src/vec/vec/impls/mpi/pvecimpl.h>
7: /*@
8: VecStrideSet - Sets a subvector of a vector defined
9: by a starting point and a stride with a given value
11: Logically Collective on Vec
13: Input Parameters:
14: + v - the vector
15: . start - starting point of the subvector (defined by a stride)
16: - s - value to set for each entry in that subvector
18: Notes:
19: One must call VecSetBlockSize() before this routine to set the stride
20: information, or use a vector created from a multicomponent DMDA.
22: This will only work if the desire subvector is a stride subvector
24: Level: advanced
26: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideScale()
27: @*/
28: PetscErrorCode VecStrideSet(Vec v,PetscInt start,PetscScalar s)
29: {
30: PetscInt i,n,bs;
31: PetscScalar *x;
34: VecGetLocalSize(v,&n);
35: VecGetBlockSize(v,&bs);
38: VecGetArray(v,&x);
39: for (i=start; i<n; i+=bs) x[i] = s;
40: VecRestoreArray(v,&x);
41: return 0;
42: }
44: /*@
45: VecStrideScale - Scales a subvector of a vector defined
46: by a starting point and a stride.
48: Logically Collective on Vec
50: Input Parameters:
51: + v - the vector
52: . start - starting point of the subvector (defined by a stride)
53: - scale - value to multiply each subvector entry by
55: Notes:
56: One must call VecSetBlockSize() before this routine to set the stride
57: information, or use a vector created from a multicomponent DMDA.
59: This will only work if the desire subvector is a stride subvector
61: Level: advanced
63: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideScale()
64: @*/
65: PetscErrorCode VecStrideScale(Vec v,PetscInt start,PetscScalar scale)
66: {
67: PetscInt i,n,bs;
68: PetscScalar *x;
71: VecGetLocalSize(v,&n);
72: VecGetBlockSize(v,&bs);
75: VecGetArray(v,&x);
76: for (i=start; i<n; i+=bs) x[i] *= scale;
77: VecRestoreArray(v,&x);
78: return 0;
79: }
81: /*@
82: VecStrideNorm - Computes the norm of subvector of a vector defined
83: by a starting point and a stride.
85: Collective on Vec
87: Input Parameters:
88: + v - the vector
89: . start - starting point of the subvector (defined by a stride)
90: - ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY
92: Output Parameter:
93: . norm - the norm
95: Notes:
96: One must call VecSetBlockSize() before this routine to set the stride
97: information, or use a vector created from a multicomponent DMDA.
99: If x is the array representing the vector x then this computes the norm
100: of the array (x[start],x[start+stride],x[start+2*stride], ....)
102: This is useful for computing, say the norm of the pressure variable when
103: the pressure is stored (interlaced) with other variables, say density etc.
105: This will only work if the desire subvector is a stride subvector
107: Level: advanced
109: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
110: @*/
111: PetscErrorCode VecStrideNorm(Vec v,PetscInt start,NormType ntype,PetscReal *nrm)
112: {
113: PetscInt i,n,bs;
114: const PetscScalar *x;
115: PetscReal tnorm;
120: VecGetLocalSize(v,&n);
121: VecGetBlockSize(v,&bs);
124: VecGetArrayRead(v,&x);
125: if (ntype == NORM_2) {
126: PetscScalar sum = 0.0;
127: for (i=start; i<n; i+=bs) sum += x[i]*(PetscConj(x[i]));
128: tnorm = PetscRealPart(sum);
129: MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)v));
130: *nrm = PetscSqrtReal(*nrm);
131: } else if (ntype == NORM_1) {
132: tnorm = 0.0;
133: for (i=start; i<n; i+=bs) tnorm += PetscAbsScalar(x[i]);
134: MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)v));
135: } else if (ntype == NORM_INFINITY) {
136: tnorm = 0.0;
137: for (i=start; i<n; i+=bs) {
138: if (PetscAbsScalar(x[i]) > tnorm) tnorm = PetscAbsScalar(x[i]);
139: }
140: MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)v));
141: } else SETERRQ(PetscObjectComm((PetscObject)v),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
142: VecRestoreArrayRead(v,&x);
143: return 0;
144: }
146: /*@
147: VecStrideMax - Computes the maximum of subvector of a vector defined
148: by a starting point and a stride and optionally its location.
150: Collective on Vec
152: Input Parameters:
153: + v - the vector
154: - start - starting point of the subvector (defined by a stride)
156: Output Parameters:
157: + idex - the location where the maximum occurred (pass NULL if not required)
158: - nrm - the maximum value in the subvector
160: Notes:
161: One must call VecSetBlockSize() before this routine to set the stride
162: information, or use a vector created from a multicomponent DMDA.
164: If xa is the array representing the vector x, then this computes the max
165: of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)
167: This is useful for computing, say the maximum of the pressure variable when
168: the pressure is stored (interlaced) with other variables, e.g., density, etc.
169: This will only work if the desire subvector is a stride subvector.
171: Level: advanced
173: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
174: @*/
175: PetscErrorCode VecStrideMax(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
176: {
177: PetscInt i,n,bs,id = -1;
178: const PetscScalar *x;
179: PetscReal max = PETSC_MIN_REAL;
183: VecGetLocalSize(v,&n);
184: VecGetBlockSize(v,&bs);
187: VecGetArrayRead(v,&x);
188: for (i=start; i<n; i+=bs) {
189: if (PetscRealPart(x[i]) > max) { max = PetscRealPart(x[i]); id = i;}
190: }
191: VecRestoreArrayRead(v,&x);
192: #if defined(PETSC_HAVE_MPIUNI)
193: *nrm = max;
194: if (idex) *idex = id;
195: #else
196: if (!idex) {
197: MPIU_Allreduce(&max,nrm,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)v));
198: } else {
199: struct { PetscReal v; PetscInt i; } in,out;
200: PetscInt rstart;
202: VecGetOwnershipRange(v,&rstart,NULL);
203: in.v = max;
204: in.i = rstart+id;
205: MPIU_Allreduce(&in,&out,1,MPIU_REAL_INT,MPIU_MAXLOC,PetscObjectComm((PetscObject)v));
206: *nrm = out.v;
207: *idex = out.i;
208: }
209: #endif
210: return 0;
211: }
213: /*@
214: VecStrideMin - Computes the minimum of subvector of a vector defined
215: by a starting point and a stride and optionally its location.
217: Collective on Vec
219: Input Parameters:
220: + v - the vector
221: - start - starting point of the subvector (defined by a stride)
223: Output Parameters:
224: + idex - the location where the minimum occurred. (pass NULL if not required)
225: - nrm - the minimum value in the subvector
227: Level: advanced
229: Notes:
230: One must call VecSetBlockSize() before this routine to set the stride
231: information, or use a vector created from a multicomponent DMDA.
233: If xa is the array representing the vector x, then this computes the min
234: of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)
236: This is useful for computing, say the minimum of the pressure variable when
237: the pressure is stored (interlaced) with other variables, e.g., density, etc.
238: This will only work if the desire subvector is a stride subvector.
240: .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
241: @*/
242: PetscErrorCode VecStrideMin(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
243: {
244: PetscInt i,n,bs,id = -1;
245: const PetscScalar *x;
246: PetscReal min = PETSC_MAX_REAL;
250: VecGetLocalSize(v,&n);
251: VecGetBlockSize(v,&bs);
254: VecGetArrayRead(v,&x);
255: for (i=start; i<n; i+=bs) {
256: if (PetscRealPart(x[i]) < min) { min = PetscRealPart(x[i]); id = i;}
257: }
258: VecRestoreArrayRead(v,&x);
259: #if defined(PETSC_HAVE_MPIUNI)
260: *nrm = min;
261: if (idex) *idex = id;
262: #else
263: if (!idex) {
264: MPIU_Allreduce(&min,nrm,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)v));
265: } else {
266: struct { PetscReal v; PetscInt i; } in,out;
267: PetscInt rstart;
269: VecGetOwnershipRange(v,&rstart,NULL);
270: in.v = min;
271: in.i = rstart+id;
272: MPIU_Allreduce(&in,&out,1,MPIU_REAL_INT,MPIU_MINLOC,PetscObjectComm((PetscObject)v));
273: *nrm = out.v;
274: *idex = out.i;
275: }
276: #endif
277: return 0;
278: }
280: /*@
281: VecStrideScaleAll - Scales the subvectors of a vector defined
282: by a starting point and a stride.
284: Logically Collective on Vec
286: Input Parameters:
287: + v - the vector
288: - scales - values to multiply each subvector entry by
290: Notes:
291: One must call VecSetBlockSize() before this routine to set the stride
292: information, or use a vector created from a multicomponent DMDA.
294: The dimension of scales must be the same as the vector block size
296: Level: advanced
298: .seealso: VecNorm(), VecStrideScale(), VecScale(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
299: @*/
300: PetscErrorCode VecStrideScaleAll(Vec v,const PetscScalar *scales)
301: {
302: PetscInt i,j,n,bs;
303: PetscScalar *x;
307: VecGetLocalSize(v,&n);
308: VecGetBlockSize(v,&bs);
309: VecGetArray(v,&x);
310: /* need to provide optimized code for each bs */
311: for (i=0; i<n; i+=bs) {
312: for (j=0; j<bs; j++) x[i+j] *= scales[j];
313: }
314: VecRestoreArray(v,&x);
315: return 0;
316: }
318: /*@
319: VecStrideNormAll - Computes the norms of subvectors of a vector defined
320: by a starting point and a stride.
322: Collective on Vec
324: Input Parameters:
325: + v - the vector
326: - ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY
328: Output Parameter:
329: . nrm - the norms
331: Notes:
332: One must call VecSetBlockSize() before this routine to set the stride
333: information, or use a vector created from a multicomponent DMDA.
335: If x is the array representing the vector x then this computes the norm
336: of the array (x[start],x[start+stride],x[start+2*stride], ....) for each start < stride
338: The dimension of nrm must be the same as the vector block size
340: This will only work if the desire subvector is a stride subvector
342: Level: advanced
344: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
345: @*/
346: PetscErrorCode VecStrideNormAll(Vec v,NormType ntype,PetscReal nrm[])
347: {
348: PetscInt i,j,n,bs;
349: const PetscScalar *x;
350: PetscReal tnorm[128];
351: MPI_Comm comm;
356: VecGetLocalSize(v,&n);
357: VecGetArrayRead(v,&x);
358: PetscObjectGetComm((PetscObject)v,&comm);
360: VecGetBlockSize(v,&bs);
363: if (ntype == NORM_2) {
364: PetscScalar sum[128];
365: for (j=0; j<bs; j++) sum[j] = 0.0;
366: for (i=0; i<n; i+=bs) {
367: for (j=0; j<bs; j++) sum[j] += x[i+j]*(PetscConj(x[i+j]));
368: }
369: for (j=0; j<bs; j++) tnorm[j] = PetscRealPart(sum[j]);
371: MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_SUM,comm);
372: for (j=0; j<bs; j++) nrm[j] = PetscSqrtReal(nrm[j]);
373: } else if (ntype == NORM_1) {
374: for (j=0; j<bs; j++) tnorm[j] = 0.0;
376: for (i=0; i<n; i+=bs) {
377: for (j=0; j<bs; j++) tnorm[j] += PetscAbsScalar(x[i+j]);
378: }
380: MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_SUM,comm);
381: } else if (ntype == NORM_INFINITY) {
382: PetscReal tmp;
383: for (j=0; j<bs; j++) tnorm[j] = 0.0;
385: for (i=0; i<n; i+=bs) {
386: for (j=0; j<bs; j++) {
387: if ((tmp = PetscAbsScalar(x[i+j])) > tnorm[j]) tnorm[j] = tmp;
388: /* check special case of tmp == NaN */
389: if (tmp != tmp) {tnorm[j] = tmp; break;}
390: }
391: }
392: MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_MAX,comm);
393: } else SETERRQ(comm,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
394: VecRestoreArrayRead(v,&x);
395: return 0;
396: }
398: /*@
399: VecStrideMaxAll - Computes the maximums of subvectors of a vector defined
400: by a starting point and a stride and optionally its location.
402: Collective on Vec
404: Input Parameter:
405: . v - the vector
407: Output Parameters:
408: + index - the location where the maximum occurred (not supported, pass NULL,
409: if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
410: - nrm - the maximum values of each subvector
412: Notes:
413: One must call VecSetBlockSize() before this routine to set the stride
414: information, or use a vector created from a multicomponent DMDA.
416: The dimension of nrm must be the same as the vector block size
418: Level: advanced
420: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
421: @*/
422: PetscErrorCode VecStrideMaxAll(Vec v,PetscInt idex[],PetscReal nrm[])
423: {
424: PetscInt i,j,n,bs;
425: const PetscScalar *x;
426: PetscReal max[128],tmp;
427: MPI_Comm comm;
432: VecGetLocalSize(v,&n);
433: VecGetArrayRead(v,&x);
434: PetscObjectGetComm((PetscObject)v,&comm);
436: VecGetBlockSize(v,&bs);
439: if (!n) {
440: for (j=0; j<bs; j++) max[j] = PETSC_MIN_REAL;
441: } else {
442: for (j=0; j<bs; j++) max[j] = PetscRealPart(x[j]);
444: for (i=bs; i<n; i+=bs) {
445: for (j=0; j<bs; j++) {
446: if ((tmp = PetscRealPart(x[i+j])) > max[j]) max[j] = tmp;
447: }
448: }
449: }
450: MPIU_Allreduce(max,nrm,bs,MPIU_REAL,MPIU_MAX,comm);
452: VecRestoreArrayRead(v,&x);
453: return 0;
454: }
456: /*@
457: VecStrideMinAll - Computes the minimum of subvector of a vector defined
458: by a starting point and a stride and optionally its location.
460: Collective on Vec
462: Input Parameter:
463: . v - the vector
465: Output Parameters:
466: + idex - the location where the minimum occurred (not supported, pass NULL,
467: if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
468: - nrm - the minimums of each subvector
470: Level: advanced
472: Notes:
473: One must call VecSetBlockSize() before this routine to set the stride
474: information, or use a vector created from a multicomponent DMDA.
476: The dimension of nrm must be the same as the vector block size
478: .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
479: @*/
480: PetscErrorCode VecStrideMinAll(Vec v,PetscInt idex[],PetscReal nrm[])
481: {
482: PetscInt i,n,bs,j;
483: const PetscScalar *x;
484: PetscReal min[128],tmp;
485: MPI_Comm comm;
490: VecGetLocalSize(v,&n);
491: VecGetArrayRead(v,&x);
492: PetscObjectGetComm((PetscObject)v,&comm);
494: VecGetBlockSize(v,&bs);
497: if (!n) {
498: for (j=0; j<bs; j++) min[j] = PETSC_MAX_REAL;
499: } else {
500: for (j=0; j<bs; j++) min[j] = PetscRealPart(x[j]);
502: for (i=bs; i<n; i+=bs) {
503: for (j=0; j<bs; j++) {
504: if ((tmp = PetscRealPart(x[i+j])) < min[j]) min[j] = tmp;
505: }
506: }
507: }
508: MPIU_Allreduce(min,nrm,bs,MPIU_REAL,MPIU_MIN,comm);
510: VecRestoreArrayRead(v,&x);
511: return 0;
512: }
514: /*----------------------------------------------------------------------------------------------*/
515: /*@
516: VecStrideGatherAll - Gathers all the single components from a multi-component vector into
517: separate vectors.
519: Collective on Vec
521: Input Parameters:
522: + v - the vector
523: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
525: Output Parameter:
526: . s - the location where the subvectors are stored
528: Notes:
529: One must call VecSetBlockSize() before this routine to set the stride
530: information, or use a vector created from a multicomponent DMDA.
532: If x is the array representing the vector x then this gathers
533: the arrays (x[start],x[start+stride],x[start+2*stride], ....)
534: for start=0,1,2,...bs-1
536: The parallel layout of the vector and the subvector must be the same;
537: i.e., nlocal of v = stride*(nlocal of s)
539: Not optimized; could be easily
541: Level: advanced
543: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
544: VecStrideScatterAll()
545: @*/
546: PetscErrorCode VecStrideGatherAll(Vec v,Vec s[],InsertMode addv)
547: {
548: PetscInt i,n,n2,bs,j,k,*bss = NULL,nv,jj,nvc;
549: PetscScalar **y;
550: const PetscScalar *x;
555: VecGetLocalSize(v,&n);
556: VecGetLocalSize(s[0],&n2);
557: VecGetArrayRead(v,&x);
558: VecGetBlockSize(v,&bs);
561: PetscMalloc2(bs,&y,bs,&bss);
562: nv = 0;
563: nvc = 0;
564: for (i=0; i<bs; i++) {
565: VecGetBlockSize(s[i],&bss[i]);
566: if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1 Re: [PETSC #8241] VecStrideGatherAll */
567: VecGetArray(s[i],&y[i]);
568: nvc += bss[i];
569: nv++;
571: if (nvc == bs) break;
572: }
574: n = n/bs;
576: jj = 0;
577: if (addv == INSERT_VALUES) {
578: for (j=0; j<nv; j++) {
579: for (k=0; k<bss[j]; k++) {
580: for (i=0; i<n; i++) y[j][i*bss[j] + k] = x[bs*i+jj+k];
581: }
582: jj += bss[j];
583: }
584: } else if (addv == ADD_VALUES) {
585: for (j=0; j<nv; j++) {
586: for (k=0; k<bss[j]; k++) {
587: for (i=0; i<n; i++) y[j][i*bss[j] + k] += x[bs*i+jj+k];
588: }
589: jj += bss[j];
590: }
591: #if !defined(PETSC_USE_COMPLEX)
592: } else if (addv == MAX_VALUES) {
593: for (j=0; j<nv; j++) {
594: for (k=0; k<bss[j]; k++) {
595: for (i=0; i<n; i++) y[j][i*bss[j] + k] = PetscMax(y[j][i*bss[j] + k],x[bs*i+jj+k]);
596: }
597: jj += bss[j];
598: }
599: #endif
600: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
602: VecRestoreArrayRead(v,&x);
603: for (i=0; i<nv; i++) {
604: VecRestoreArray(s[i],&y[i]);
605: }
607: PetscFree2(y,bss);
608: return 0;
609: }
611: /*@
612: VecStrideScatterAll - Scatters all the single components from separate vectors into
613: a multi-component vector.
615: Collective on Vec
617: Input Parameters:
618: + s - the location where the subvectors are stored
619: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
621: Output Parameter:
622: . v - the multicomponent vector
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: The parallel layout of the vector and the subvector must be the same;
629: i.e., nlocal of v = stride*(nlocal of s)
631: Not optimized; could be easily
633: Level: advanced
635: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
636: VecStrideScatterAll()
637: @*/
638: PetscErrorCode VecStrideScatterAll(Vec s[],Vec v,InsertMode addv)
639: {
640: PetscInt i,n,n2,bs,j,jj,k,*bss = NULL,nv,nvc;
641: PetscScalar *x;
642: PetscScalar const **y;
647: VecGetLocalSize(v,&n);
648: VecGetLocalSize(s[0],&n2);
649: VecGetArray(v,&x);
650: VecGetBlockSize(v,&bs);
653: PetscMalloc2(bs,(PetscScalar***)&y,bs,&bss);
654: nv = 0;
655: nvc = 0;
656: for (i=0; i<bs; i++) {
657: VecGetBlockSize(s[i],&bss[i]);
658: if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1 Re: [PETSC #8241] VecStrideGatherAll */
659: VecGetArrayRead(s[i],&y[i]);
660: nvc += bss[i];
661: nv++;
663: if (nvc == bs) break;
664: }
666: n = n/bs;
668: jj = 0;
669: if (addv == INSERT_VALUES) {
670: for (j=0; j<nv; j++) {
671: for (k=0; k<bss[j]; k++) {
672: for (i=0; i<n; i++) x[bs*i+jj+k] = y[j][i*bss[j] + k];
673: }
674: jj += bss[j];
675: }
676: } else if (addv == ADD_VALUES) {
677: for (j=0; j<nv; j++) {
678: for (k=0; k<bss[j]; k++) {
679: for (i=0; i<n; i++) x[bs*i+jj+k] += y[j][i*bss[j] + k];
680: }
681: jj += bss[j];
682: }
683: #if !defined(PETSC_USE_COMPLEX)
684: } else if (addv == MAX_VALUES) {
685: for (j=0; j<nv; j++) {
686: for (k=0; k<bss[j]; k++) {
687: for (i=0; i<n; i++) x[bs*i+jj+k] = PetscMax(x[bs*i+jj+k],y[j][i*bss[j] + k]);
688: }
689: jj += bss[j];
690: }
691: #endif
692: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
694: VecRestoreArray(v,&x);
695: for (i=0; i<nv; i++) {
696: VecRestoreArrayRead(s[i],&y[i]);
697: }
698: PetscFree2(*(PetscScalar***)&y,bss);
699: return 0;
700: }
702: /*@
703: VecStrideGather - Gathers a single component from a multi-component vector into
704: another vector.
706: Collective on Vec
708: Input Parameters:
709: + v - the vector
710: . start - starting point of the subvector (defined by a stride)
711: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
713: Output Parameter:
714: . s - the location where the subvector is stored
716: Notes:
717: One must call VecSetBlockSize() before this routine to set the stride
718: information, or use a vector created from a multicomponent DMDA.
720: If x is the array representing the vector x then this gathers
721: the array (x[start],x[start+stride],x[start+2*stride], ....)
723: The parallel layout of the vector and the subvector must be the same;
724: i.e., nlocal of v = stride*(nlocal of s)
726: Not optimized; could be easily
728: Level: advanced
730: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
731: VecStrideScatterAll()
732: @*/
733: PetscErrorCode VecStrideGather(Vec v,PetscInt start,Vec s,InsertMode addv)
734: {
740: (*v->ops->stridegather)(v,start,s,addv);
741: return 0;
742: }
744: /*@
745: VecStrideScatter - Scatters a single component from a vector into a multi-component vector.
747: Collective on Vec
749: Input Parameters:
750: + s - the single-component vector
751: . start - starting point of the subvector (defined by a stride)
752: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
754: Output Parameter:
755: . v - the location where the subvector is scattered (the multi-component vector)
757: Notes:
758: One must call VecSetBlockSize() on the multi-component vector before this
759: routine to set the stride information, or use a vector created from a multicomponent DMDA.
761: The parallel layout of the vector and the subvector must be the same;
762: i.e., nlocal of v = stride*(nlocal of s)
764: Not optimized; could be easily
766: Level: advanced
768: .seealso: VecStrideNorm(), VecStrideGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
769: VecStrideScatterAll(), VecStrideSubSetScatter(), VecStrideSubSetGather()
770: @*/
771: PetscErrorCode VecStrideScatter(Vec s,PetscInt start,Vec v,InsertMode addv)
772: {
778: (*v->ops->stridescatter)(s,start,v,addv);
779: return 0;
780: }
782: /*@
783: VecStrideSubSetGather - Gathers a subset of components from a multi-component vector into
784: another vector.
786: Collective on Vec
788: Input Parameters:
789: + v - the vector
790: . nidx - the number of indices
791: . idxv - the indices of the components 0 <= idxv[0] ...idxv[nidx-1] < bs(v), they need not be sorted
792: . idxs - the indices of the components 0 <= idxs[0] ...idxs[nidx-1] < bs(s), they need not be sorted, may be null if nidx == bs(s) or is PETSC_DETERMINE
793: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
795: Output Parameter:
796: . s - the location where the subvector is stored
798: Notes:
799: One must call VecSetBlockSize() on both vectors before this routine to set the stride
800: information, or use a vector created from a multicomponent DMDA.
802: The parallel layout of the vector and the subvector must be the same;
804: Not optimized; could be easily
806: Level: advanced
808: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideGather(), VecStrideSubSetScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
809: VecStrideScatterAll()
810: @*/
811: PetscErrorCode VecStrideSubSetGather(Vec v,PetscInt nidx,const PetscInt idxv[],const PetscInt idxs[],Vec s,InsertMode addv)
812: {
815: if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
817: (*v->ops->stridesubsetgather)(v,nidx,idxv,idxs,s,addv);
818: return 0;
819: }
821: /*@
822: VecStrideSubSetScatter - Scatters components from a vector into a subset of components of a multi-component vector.
824: Collective on Vec
826: Input Parameters:
827: + s - the smaller-component vector
828: . nidx - the number of indices in idx
829: . idxs - the indices of the components in the smaller-component vector, 0 <= idxs[0] ...idxs[nidx-1] < bs(s) they need not be sorted, may be null if nidx == bs(s) or is PETSC_DETERMINE
830: . idxv - the indices of the components in the larger-component vector, 0 <= idx[0] ...idx[nidx-1] < bs(v) they need not be sorted
831: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
833: Output Parameter:
834: . v - the location where the subvector is into scattered (the multi-component vector)
836: Notes:
837: One must call VecSetBlockSize() on the vectors before this
838: routine to set the stride information, or use a vector created from a multicomponent DMDA.
840: The parallel layout of the vector and the subvector must be the same;
842: Not optimized; could be easily
844: Level: advanced
846: .seealso: VecStrideNorm(), VecStrideGather(), VecStrideGather(), VecStrideSubSetGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
847: VecStrideScatterAll()
848: @*/
849: PetscErrorCode VecStrideSubSetScatter(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)
850: {
853: if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
855: (*v->ops->stridesubsetscatter)(s,nidx,idxs,idxv,v,addv);
856: return 0;
857: }
859: PetscErrorCode VecStrideGather_Default(Vec v,PetscInt start,Vec s,InsertMode addv)
860: {
861: PetscInt i,n,bs,ns;
862: const PetscScalar *x;
863: PetscScalar *y;
865: VecGetLocalSize(v,&n);
866: VecGetLocalSize(s,&ns);
867: VecGetArrayRead(v,&x);
868: VecGetArray(s,&y);
870: bs = v->map->bs;
872: x += start;
873: n = n/bs;
875: if (addv == INSERT_VALUES) {
876: for (i=0; i<n; i++) y[i] = x[bs*i];
877: } else if (addv == ADD_VALUES) {
878: for (i=0; i<n; i++) y[i] += x[bs*i];
879: #if !defined(PETSC_USE_COMPLEX)
880: } else if (addv == MAX_VALUES) {
881: for (i=0; i<n; i++) y[i] = PetscMax(y[i],x[bs*i]);
882: #endif
883: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
885: VecRestoreArrayRead(v,&x);
886: VecRestoreArray(s,&y);
887: return 0;
888: }
890: PetscErrorCode VecStrideScatter_Default(Vec s,PetscInt start,Vec v,InsertMode addv)
891: {
892: PetscInt i,n,bs,ns;
893: PetscScalar *x;
894: const PetscScalar *y;
896: VecGetLocalSize(v,&n);
897: VecGetLocalSize(s,&ns);
898: VecGetArray(v,&x);
899: VecGetArrayRead(s,&y);
901: bs = v->map->bs;
903: x += start;
904: n = n/bs;
906: if (addv == INSERT_VALUES) {
907: for (i=0; i<n; i++) x[bs*i] = y[i];
908: } else if (addv == ADD_VALUES) {
909: for (i=0; i<n; i++) x[bs*i] += y[i];
910: #if !defined(PETSC_USE_COMPLEX)
911: } else if (addv == MAX_VALUES) {
912: for (i=0; i<n; i++) x[bs*i] = PetscMax(y[i],x[bs*i]);
913: #endif
914: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
916: VecRestoreArray(v,&x);
917: VecRestoreArrayRead(s,&y);
918: return 0;
919: }
921: PetscErrorCode VecStrideSubSetGather_Default(Vec v,PetscInt nidx,const PetscInt idxv[],const PetscInt idxs[],Vec s,InsertMode addv)
922: {
923: PetscInt i,j,n,bs,bss,ns;
924: const PetscScalar *x;
925: PetscScalar *y;
927: VecGetLocalSize(v,&n);
928: VecGetLocalSize(s,&ns);
929: VecGetArrayRead(v,&x);
930: VecGetArray(s,&y);
932: bs = v->map->bs;
933: bss = s->map->bs;
934: n = n/bs;
936: if (PetscDefined(USE_DEBUG)) {
938: for (j=0; j<nidx; j++) {
941: }
943: }
945: if (addv == INSERT_VALUES) {
946: if (!idxs) {
947: for (i=0; i<n; i++) {
948: for (j=0; j<bss; j++) y[bss*i+j] = x[bs*i+idxv[j]];
949: }
950: } else {
951: for (i=0; i<n; i++) {
952: for (j=0; j<bss; j++) y[bss*i+idxs[j]] = x[bs*i+idxv[j]];
953: }
954: }
955: } else if (addv == ADD_VALUES) {
956: if (!idxs) {
957: for (i=0; i<n; i++) {
958: for (j=0; j<bss; j++) y[bss*i+j] += x[bs*i+idxv[j]];
959: }
960: } else {
961: for (i=0; i<n; i++) {
962: for (j=0; j<bss; j++) y[bss*i+idxs[j]] += x[bs*i+idxv[j]];
963: }
964: }
965: #if !defined(PETSC_USE_COMPLEX)
966: } else if (addv == MAX_VALUES) {
967: if (!idxs) {
968: for (i=0; i<n; i++) {
969: for (j=0; j<bss; j++) y[bss*i+j] = PetscMax(y[bss*i+j],x[bs*i+idxv[j]]);
970: }
971: } else {
972: for (i=0; i<n; i++) {
973: for (j=0; j<bss; j++) y[bss*i+idxs[j]] = PetscMax(y[bss*i+idxs[j]],x[bs*i+idxv[j]]);
974: }
975: }
976: #endif
977: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
979: VecRestoreArrayRead(v,&x);
980: VecRestoreArray(s,&y);
981: return 0;
982: }
984: PetscErrorCode VecStrideSubSetScatter_Default(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)
985: {
986: PetscInt j,i,n,bs,ns,bss;
987: PetscScalar *x;
988: const PetscScalar *y;
990: VecGetLocalSize(v,&n);
991: VecGetLocalSize(s,&ns);
992: VecGetArray(v,&x);
993: VecGetArrayRead(s,&y);
995: bs = v->map->bs;
996: bss = s->map->bs;
997: n = n/bs;
999: if (PetscDefined(USE_DEBUG)) {
1001: for (j=0; j<bss; j++) {
1002: if (idxs) {
1005: }
1006: }
1008: }
1010: if (addv == INSERT_VALUES) {
1011: if (!idxs) {
1012: for (i=0; i<n; i++) {
1013: for (j=0; j<bss; j++) x[bs*i + idxv[j]] = y[bss*i+j];
1014: }
1015: } else {
1016: for (i=0; i<n; i++) {
1017: for (j=0; j<bss; j++) x[bs*i + idxv[j]] = y[bss*i+idxs[j]];
1018: }
1019: }
1020: } else if (addv == ADD_VALUES) {
1021: if (!idxs) {
1022: for (i=0; i<n; i++) {
1023: for (j=0; j<bss; j++) x[bs*i + idxv[j]] += y[bss*i+j];
1024: }
1025: } else {
1026: for (i=0; i<n; i++) {
1027: for (j=0; j<bss; j++) x[bs*i + idxv[j]] += y[bss*i+idxs[j]];
1028: }
1029: }
1030: #if !defined(PETSC_USE_COMPLEX)
1031: } else if (addv == MAX_VALUES) {
1032: if (!idxs) {
1033: for (i=0; i<n; i++) {
1034: for (j=0; j<bss; j++) x[bs*i + idxv[j]] = PetscMax(y[bss*i+j],x[bs*i + idxv[j]]);
1035: }
1036: } else {
1037: for (i=0; i<n; i++) {
1038: for (j=0; j<bss; j++) x[bs*i + idxv[j]] = PetscMax(y[bss*i+idxs[j]],x[bs*i + idxv[j]]);
1039: }
1040: }
1041: #endif
1042: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
1044: VecRestoreArray(v,&x);
1045: VecRestoreArrayRead(s,&y);
1046: return 0;
1047: }
1049: PetscErrorCode VecReciprocal_Default(Vec v)
1050: {
1051: PetscInt i,n;
1052: PetscScalar *x;
1054: VecGetLocalSize(v,&n);
1055: VecGetArray(v,&x);
1056: for (i=0; i<n; i++) {
1057: if (x[i] != (PetscScalar)0.0) x[i] = (PetscScalar)1.0/x[i];
1058: }
1059: VecRestoreArray(v,&x);
1060: return 0;
1061: }
1063: /*@
1064: VecExp - Replaces each component of a vector by e^x_i
1066: Not collective
1068: Input Parameter:
1069: . v - The vector
1071: Output Parameter:
1072: . v - The vector of exponents
1074: Level: beginner
1076: .seealso: VecLog(), VecAbs(), VecSqrtAbs(), VecReciprocal()
1078: @*/
1079: PetscErrorCode VecExp(Vec v)
1080: {
1081: PetscScalar *x;
1082: PetscInt i, n;
1085: if (v->ops->exp) {
1086: (*v->ops->exp)(v);
1087: } else {
1088: VecGetLocalSize(v, &n);
1089: VecGetArray(v, &x);
1090: for (i = 0; i < n; i++) x[i] = PetscExpScalar(x[i]);
1091: VecRestoreArray(v, &x);
1092: }
1093: return 0;
1094: }
1096: /*@
1097: VecLog - Replaces each component of a vector by log(x_i), the natural logarithm
1099: Not collective
1101: Input Parameter:
1102: . v - The vector
1104: Output Parameter:
1105: . v - The vector of logs
1107: Level: beginner
1109: .seealso: VecExp(), VecAbs(), VecSqrtAbs(), VecReciprocal()
1111: @*/
1112: PetscErrorCode VecLog(Vec v)
1113: {
1114: PetscScalar *x;
1115: PetscInt i, n;
1118: if (v->ops->log) {
1119: (*v->ops->log)(v);
1120: } else {
1121: VecGetLocalSize(v, &n);
1122: VecGetArray(v, &x);
1123: for (i = 0; i < n; i++) x[i] = PetscLogScalar(x[i]);
1124: VecRestoreArray(v, &x);
1125: }
1126: return 0;
1127: }
1129: /*@
1130: VecSqrtAbs - Replaces each component of a vector by the square root of its magnitude.
1132: Not collective
1134: Input Parameter:
1135: . v - The vector
1137: Output Parameter:
1138: . v - The vector square root
1140: Level: beginner
1142: Note: The actual function is sqrt(|x_i|)
1144: .seealso: VecLog(), VecExp(), VecReciprocal(), VecAbs()
1146: @*/
1147: PetscErrorCode VecSqrtAbs(Vec v)
1148: {
1149: PetscScalar *x;
1150: PetscInt i, n;
1153: if (v->ops->sqrt) {
1154: (*v->ops->sqrt)(v);
1155: } else {
1156: VecGetLocalSize(v, &n);
1157: VecGetArray(v, &x);
1158: for (i = 0; i < n; i++) x[i] = PetscSqrtReal(PetscAbsScalar(x[i]));
1159: VecRestoreArray(v, &x);
1160: }
1161: return 0;
1162: }
1164: /*@
1165: VecDotNorm2 - computes the inner product of two vectors and the 2-norm squared of the second vector
1167: Collective on Vec
1169: Input Parameters:
1170: + s - first vector
1171: - t - second vector
1173: Output Parameters:
1174: + dp - s'conj(t)
1175: - nm - t'conj(t)
1177: Level: advanced
1179: Notes:
1180: conj(x) is the complex conjugate of x when x is complex
1182: .seealso: VecDot(), VecNorm(), VecDotBegin(), VecNormBegin(), VecDotEnd(), VecNormEnd()
1184: @*/
1185: PetscErrorCode VecDotNorm2(Vec s,Vec t,PetscScalar *dp, PetscReal *nm)
1186: {
1187: const PetscScalar *sx, *tx;
1188: PetscScalar dpx = 0.0, nmx = 0.0,work[2],sum[2];
1189: PetscInt i, n;
1201: PetscLogEventBegin(VEC_DotNorm2,s,t,0,0);
1202: if (s->ops->dotnorm2) {
1203: (*s->ops->dotnorm2)(s,t,dp,&dpx);
1204: *nm = PetscRealPart(dpx);
1205: } else {
1206: VecGetLocalSize(s, &n);
1207: VecGetArrayRead(s, &sx);
1208: VecGetArrayRead(t, &tx);
1210: for (i = 0; i<n; i++) {
1211: dpx += sx[i]*PetscConj(tx[i]);
1212: nmx += tx[i]*PetscConj(tx[i]);
1213: }
1214: work[0] = dpx;
1215: work[1] = nmx;
1217: MPIU_Allreduce(work,sum,2,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)s));
1218: *dp = sum[0];
1219: *nm = PetscRealPart(sum[1]);
1221: VecRestoreArrayRead(t, &tx);
1222: VecRestoreArrayRead(s, &sx);
1223: PetscLogFlops(4.0*n);
1224: }
1225: PetscLogEventEnd(VEC_DotNorm2,s,t,0,0);
1226: return 0;
1227: }
1229: /*@
1230: VecSum - Computes the sum of all the components of a vector.
1232: Collective on Vec
1234: Input Parameter:
1235: . v - the vector
1237: Output Parameter:
1238: . sum - the result
1240: Level: beginner
1242: .seealso: VecMean(), VecNorm()
1243: @*/
1244: PetscErrorCode VecSum(Vec v,PetscScalar *sum)
1245: {
1246: PetscInt i,n;
1247: const PetscScalar *x;
1251: *sum = 0.0;
1252: if (v->ops->sum) {
1253: (*v->ops->sum)(v,sum);
1254: } else {
1255: VecGetLocalSize(v,&n);
1256: VecGetArrayRead(v,&x);
1257: for (i=0; i<n; i++) *sum += x[i];
1258: VecRestoreArrayRead(v,&x);
1259: }
1260: MPIU_Allreduce(MPI_IN_PLACE,sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)v));
1261: return 0;
1262: }
1264: /*@
1265: VecMean - Computes the arithmetic mean of all the components of a vector.
1267: Collective on Vec
1269: Input Parameter:
1270: . v - the vector
1272: Output Parameter:
1273: . mean - the result
1275: Level: beginner
1277: .seealso: VecSum(), VecNorm()
1278: @*/
1279: PetscErrorCode VecMean(Vec v,PetscScalar *mean)
1280: {
1281: PetscInt n;
1285: VecGetSize(v,&n);
1286: VecSum(v,mean);
1287: *mean /= n;
1288: return 0;
1289: }
1291: /*@
1292: VecImaginaryPart - Replaces a complex vector with its imginary part
1294: Collective on Vec
1296: Input Parameter:
1297: . v - the vector
1299: Level: beginner
1301: .seealso: VecNorm(), VecRealPart()
1302: @*/
1303: PetscErrorCode VecImaginaryPart(Vec v)
1304: {
1305: PetscInt i,n;
1306: PetscScalar *x;
1309: VecGetLocalSize(v,&n);
1310: VecGetArray(v,&x);
1311: for (i=0; i<n; i++) x[i] = PetscImaginaryPart(x[i]);
1312: VecRestoreArray(v,&x);
1313: return 0;
1314: }
1316: /*@
1317: VecRealPart - Replaces a complex vector with its real part
1319: Collective on Vec
1321: Input Parameter:
1322: . v - the vector
1324: Level: beginner
1326: .seealso: VecNorm(), VecImaginaryPart()
1327: @*/
1328: PetscErrorCode VecRealPart(Vec v)
1329: {
1330: PetscInt i,n;
1331: PetscScalar *x;
1334: VecGetLocalSize(v,&n);
1335: VecGetArray(v,&x);
1336: for (i=0; i<n; i++) x[i] = PetscRealPart(x[i]);
1337: VecRestoreArray(v,&x);
1338: return 0;
1339: }
1341: /*@
1342: VecShift - Shifts all of the components of a vector by computing
1343: x[i] = x[i] + shift.
1345: Logically Collective on Vec
1347: Input Parameters:
1348: + v - the vector
1349: - shift - the shift
1351: Level: intermediate
1353: @*/
1354: PetscErrorCode VecShift(Vec v,PetscScalar shift)
1355: {
1356: PetscInt i,n;
1357: PetscScalar *x;
1361: VecSetErrorIfLocked(v,1);
1362: if (shift == 0.0) return 0;
1364: if (v->ops->shift) {
1365: (*v->ops->shift)(v,shift);
1366: } else {
1367: VecGetLocalSize(v,&n);
1368: VecGetArray(v,&x);
1369: for (i=0; i<n; i++) x[i] += shift;
1370: VecRestoreArray(v,&x);
1371: }
1372: return 0;
1373: }
1375: /*@
1376: VecAbs - Replaces every element in a vector with its absolute value.
1378: Logically Collective on Vec
1380: Input Parameters:
1381: . v - the vector
1383: Level: intermediate
1385: @*/
1386: PetscErrorCode VecAbs(Vec v)
1387: {
1388: PetscInt i,n;
1389: PetscScalar *x;
1392: VecSetErrorIfLocked(v,1);
1394: if (v->ops->abs) {
1395: (*v->ops->abs)(v);
1396: } else {
1397: VecGetLocalSize(v,&n);
1398: VecGetArray(v,&x);
1399: for (i=0; i<n; i++) x[i] = PetscAbsScalar(x[i]);
1400: VecRestoreArray(v,&x);
1401: }
1402: return 0;
1403: }
1405: /*@
1406: VecPermute - Permutes a vector in place using the given ordering.
1408: Input Parameters:
1409: + vec - The vector
1410: . order - The ordering
1411: - inv - The flag for inverting the permutation
1413: Level: beginner
1415: Note: This function does not yet support parallel Index Sets with non-local permutations
1417: .seealso: MatPermute()
1418: @*/
1419: PetscErrorCode VecPermute(Vec x, IS row, PetscBool inv)
1420: {
1421: const PetscScalar *array;
1422: PetscScalar *newArray;
1423: const PetscInt *idx;
1424: PetscInt i,rstart,rend;
1428: VecSetErrorIfLocked(x,1);
1429: VecGetOwnershipRange(x,&rstart,&rend);
1430: ISGetIndices(row, &idx);
1431: VecGetArrayRead(x, &array);
1432: PetscMalloc1(x->map->n, &newArray);
1433: if (PetscDefined(USE_DEBUG)) {
1434: for (i = 0; i < x->map->n; i++) {
1436: }
1437: }
1438: if (!inv) {
1439: for (i = 0; i < x->map->n; i++) newArray[i] = array[idx[i]-rstart];
1440: } else {
1441: for (i = 0; i < x->map->n; i++) newArray[idx[i]-rstart] = array[i];
1442: }
1443: VecRestoreArrayRead(x, &array);
1444: ISRestoreIndices(row, &idx);
1445: VecReplaceArray(x, newArray);
1446: return 0;
1447: }
1449: /*@
1450: VecEqual - Compares two vectors. Returns true if the two vectors are either pointing to the same memory buffer,
1451: or if the two vectors have the same local and global layout as well as bitwise equality of all entries.
1452: Does NOT take round-off errors into account.
1454: Collective on Vec
1456: Input Parameters:
1457: + vec1 - the first vector
1458: - vec2 - the second vector
1460: Output Parameter:
1461: . flg - PETSC_TRUE if the vectors are equal; PETSC_FALSE otherwise.
1463: Level: intermediate
1464: @*/
1465: PetscErrorCode VecEqual(Vec vec1,Vec vec2,PetscBool *flg)
1466: {
1467: const PetscScalar *v1,*v2;
1468: PetscInt n1,n2,N1,N2;
1469: PetscBool flg1;
1474: if (vec1 == vec2) *flg = PETSC_TRUE;
1475: else {
1476: VecGetSize(vec1,&N1);
1477: VecGetSize(vec2,&N2);
1478: if (N1 != N2) flg1 = PETSC_FALSE;
1479: else {
1480: VecGetLocalSize(vec1,&n1);
1481: VecGetLocalSize(vec2,&n2);
1482: if (n1 != n2) flg1 = PETSC_FALSE;
1483: else {
1484: VecGetArrayRead(vec1,&v1);
1485: VecGetArrayRead(vec2,&v2);
1486: PetscArraycmp(v1,v2,n1,&flg1);
1487: VecRestoreArrayRead(vec1,&v1);
1488: VecRestoreArrayRead(vec2,&v2);
1489: }
1490: }
1491: /* combine results from all processors */
1492: MPIU_Allreduce(&flg1,flg,1,MPIU_BOOL,MPI_MIN,PetscObjectComm((PetscObject)vec1));
1493: }
1494: return 0;
1495: }
1497: /*@
1498: VecUniqueEntries - Compute the number of unique entries, and those entries
1500: Collective on Vec
1502: Input Parameter:
1503: . vec - the vector
1505: Output Parameters:
1506: + n - The number of unique entries
1507: - e - The entries
1509: Level: intermediate
1511: @*/
1512: PetscErrorCode VecUniqueEntries(Vec vec, PetscInt *n, PetscScalar **e)
1513: {
1514: const PetscScalar *v;
1515: PetscScalar *tmp, *vals;
1516: PetscMPIInt *N, *displs, l;
1517: PetscInt ng, m, i, j, p;
1518: PetscMPIInt size;
1522: MPI_Comm_size(PetscObjectComm((PetscObject) vec), &size);
1523: VecGetLocalSize(vec, &m);
1524: VecGetArrayRead(vec, &v);
1525: PetscMalloc2(m,&tmp,size,&N);
1526: for (i = 0, j = 0, l = 0; i < m; ++i) {
1527: /* Can speed this up with sorting */
1528: for (j = 0; j < l; ++j) {
1529: if (v[i] == tmp[j]) break;
1530: }
1531: if (j == l) {
1532: tmp[j] = v[i];
1533: ++l;
1534: }
1535: }
1536: VecRestoreArrayRead(vec, &v);
1537: /* Gather serial results */
1538: MPI_Allgather(&l, 1, MPI_INT, N, 1, MPI_INT, PetscObjectComm((PetscObject) vec));
1539: for (p = 0, ng = 0; p < size; ++p) {
1540: ng += N[p];
1541: }
1542: PetscMalloc2(ng,&vals,size+1,&displs);
1543: for (p = 1, displs[0] = 0; p <= size; ++p) {
1544: displs[p] = displs[p-1] + N[p-1];
1545: }
1546: MPI_Allgatherv(tmp, l, MPIU_SCALAR, vals, N, displs, MPIU_SCALAR, PetscObjectComm((PetscObject) vec));
1547: /* Find unique entries */
1548: #ifdef PETSC_USE_COMPLEX
1549: SETERRQ(PetscObjectComm((PetscObject) vec), PETSC_ERR_SUP, "Does not work with complex numbers");
1550: #else
1551: *n = displs[size];
1552: PetscSortRemoveDupsReal(n, (PetscReal *) vals);
1553: if (e) {
1555: PetscMalloc1(*n, e);
1556: for (i = 0; i < *n; ++i) {
1557: (*e)[i] = vals[i];
1558: }
1559: }
1560: PetscFree2(vals,displs);
1561: PetscFree2(tmp,N);
1562: return 0;
1563: #endif
1564: }