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: {
31: PetscInt i,n,bs;
32: PetscScalar *x;
36: VecGetLocalSize(v,&n);
37: VecGetBlockSize(v,&bs);
38: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
39: else if (start >= bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Start of stride subvector (%D) is too large for stride %D",start,bs);
40: VecGetArray(v,&x);
41: for (i=start; i<n; i+=bs) x[i] = s;
42: VecRestoreArray(v,&x);
43: return(0);
44: }
46: /*@
47: VecStrideScale - Scales a subvector of a vector defined
48: by a starting point and a stride.
50: Logically Collective on Vec
52: Input Parameters:
53: + v - the vector
54: . start - starting point of the subvector (defined by a stride)
55: - scale - value to multiply each subvector entry by
57: Notes:
58: One must call VecSetBlockSize() before this routine to set the stride
59: information, or use a vector created from a multicomponent DMDA.
61: This will only work if the desire subvector is a stride subvector
63: Level: advanced
65: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideScale()
66: @*/
67: PetscErrorCode VecStrideScale(Vec v,PetscInt start,PetscScalar scale)
68: {
70: PetscInt i,n,bs;
71: PetscScalar *x;
75: VecGetLocalSize(v,&n);
76: VecGetBlockSize(v,&bs);
77: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
78: else if (start >= bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Start of stride subvector (%D) is too large for stride %D",start,bs);
79: VecGetArray(v,&x);
80: for (i=start; i<n; i+=bs) x[i] *= scale;
81: VecRestoreArray(v,&x);
82: return(0);
83: }
85: /*@
86: VecStrideNorm - Computes the norm of subvector of a vector defined
87: by a starting point and a stride.
89: Collective on Vec
91: Input Parameters:
92: + v - the vector
93: . start - starting point of the subvector (defined by a stride)
94: - ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY
96: Output Parameter:
97: . norm - the norm
99: Notes:
100: One must call VecSetBlockSize() before this routine to set the stride
101: information, or use a vector created from a multicomponent DMDA.
103: If x is the array representing the vector x then this computes the norm
104: of the array (x[start],x[start+stride],x[start+2*stride], ....)
106: This is useful for computing, say the norm of the pressure variable when
107: the pressure is stored (interlaced) with other variables, say density etc.
109: This will only work if the desire subvector is a stride subvector
111: Level: advanced
113: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
114: @*/
115: PetscErrorCode VecStrideNorm(Vec v,PetscInt start,NormType ntype,PetscReal *nrm)
116: {
117: PetscErrorCode ierr;
118: PetscInt i,n,bs;
119: const PetscScalar *x;
120: PetscReal tnorm;
126: VecGetLocalSize(v,&n);
127: VecGetBlockSize(v,&bs);
128: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
129: else if (start >= bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Start of stride subvector (%D) is too large for stride %D",start,bs);
130: VecGetArrayRead(v,&x);
131: if (ntype == NORM_2) {
132: PetscScalar sum = 0.0;
133: for (i=start; i<n; i+=bs) sum += x[i]*(PetscConj(x[i]));
134: tnorm = PetscRealPart(sum);
135: MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)v));
136: *nrm = PetscSqrtReal(*nrm);
137: } else if (ntype == NORM_1) {
138: tnorm = 0.0;
139: for (i=start; i<n; i+=bs) tnorm += PetscAbsScalar(x[i]);
140: MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)v));
141: } else if (ntype == NORM_INFINITY) {
142: tnorm = 0.0;
143: for (i=start; i<n; i+=bs) {
144: if (PetscAbsScalar(x[i]) > tnorm) tnorm = PetscAbsScalar(x[i]);
145: }
146: MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)v));
147: } else SETERRQ(PetscObjectComm((PetscObject)v),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
148: VecRestoreArrayRead(v,&x);
149: return(0);
150: }
152: /*@
153: VecStrideMax - Computes the maximum of subvector of a vector defined
154: by a starting point and a stride and optionally its location.
156: Collective on Vec
158: Input Parameters:
159: + v - the vector
160: - start - starting point of the subvector (defined by a stride)
162: Output Parameters:
163: + idex - the location where the maximum occurred (pass NULL if not required)
164: - nrm - the maximum value in the subvector
166: Notes:
167: One must call VecSetBlockSize() before this routine to set the stride
168: information, or use a vector created from a multicomponent DMDA.
170: If xa is the array representing the vector x, then this computes the max
171: of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)
173: This is useful for computing, say the maximum of the pressure variable when
174: the pressure is stored (interlaced) with other variables, e.g., density, etc.
175: This will only work if the desire subvector is a stride subvector.
177: Level: advanced
179: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
180: @*/
181: PetscErrorCode VecStrideMax(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
182: {
183: PetscErrorCode ierr;
184: PetscInt i,n,bs,id = -1;
185: const PetscScalar *x;
186: PetscReal max = PETSC_MIN_REAL;
191: VecGetLocalSize(v,&n);
192: VecGetBlockSize(v,&bs);
193: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
194: else if (start >= bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Start of stride subvector (%D) is too large for stride %D",start,bs);
195: VecGetArrayRead(v,&x);
196: for (i=start; i<n; i+=bs) {
197: if (PetscRealPart(x[i]) > max) { max = PetscRealPart(x[i]); id = i;}
198: }
199: VecRestoreArrayRead(v,&x);
200: #if defined(PETSC_HAVE_MPIUNI)
201: *nrm = max;
202: if (idex) *idex = id;
203: #else
204: if (!idex) {
205: MPIU_Allreduce(&max,nrm,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)v));
206: } else {
207: struct { PetscReal v; PetscInt i; } in,out;
208: PetscInt rstart;
210: VecGetOwnershipRange(v,&rstart,NULL);
211: in.v = max;
212: in.i = rstart+id;
213: MPIU_Allreduce(&in,&out,1,MPIU_REAL_INT,MPIU_MAXLOC,PetscObjectComm((PetscObject)v));
214: *nrm = out.v;
215: *idex = out.i;
216: }
217: #endif
218: return(0);
219: }
221: /*@
222: VecStrideMin - Computes the minimum of subvector of a vector defined
223: by a starting point and a stride and optionally its location.
225: Collective on Vec
227: Input Parameters:
228: + v - the vector
229: - start - starting point of the subvector (defined by a stride)
231: Output Parameters:
232: + idex - the location where the minimum occurred. (pass NULL if not required)
233: - nrm - the minimum value in the subvector
235: Level: advanced
237: Notes:
238: One must call VecSetBlockSize() before this routine to set the stride
239: information, or use a vector created from a multicomponent DMDA.
241: If xa is the array representing the vector x, then this computes the min
242: of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)
244: This is useful for computing, say the minimum of the pressure variable when
245: the pressure is stored (interlaced) with other variables, e.g., density, etc.
246: This will only work if the desire subvector is a stride subvector.
248: .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
249: @*/
250: PetscErrorCode VecStrideMin(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
251: {
252: PetscErrorCode ierr;
253: PetscInt i,n,bs,id = -1;
254: const PetscScalar *x;
255: PetscReal min = PETSC_MAX_REAL;
260: VecGetLocalSize(v,&n);
261: VecGetBlockSize(v,&bs);
262: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
263: else if (start >= bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Start of stride subvector (%D) is too large for stride %D",start,bs);
264: VecGetArrayRead(v,&x);
265: for (i=start; i<n; i+=bs) {
266: if (PetscRealPart(x[i]) < min) { min = PetscRealPart(x[i]); id = i;}
267: }
268: VecRestoreArrayRead(v,&x);
269: #if defined(PETSC_HAVE_MPIUNI)
270: *nrm = min;
271: if (idex) *idex = id;
272: #else
273: if (!idex) {
274: MPIU_Allreduce(&min,nrm,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)v));
275: } else {
276: struct { PetscReal v; PetscInt i; } in,out;
277: PetscInt rstart;
279: VecGetOwnershipRange(v,&rstart,NULL);
280: in.v = min;
281: in.i = rstart+id;
282: MPIU_Allreduce(&in,&out,1,MPIU_REAL_INT,MPIU_MINLOC,PetscObjectComm((PetscObject)v));
283: *nrm = out.v;
284: *idex = out.i;
285: }
286: #endif
287: return(0);
288: }
290: /*@
291: VecStrideScaleAll - Scales the subvectors of a vector defined
292: by a starting point and a stride.
294: Logically Collective on Vec
296: Input Parameters:
297: + v - the vector
298: - scales - values to multiply each subvector entry by
300: Notes:
301: One must call VecSetBlockSize() before this routine to set the stride
302: information, or use a vector created from a multicomponent DMDA.
304: The dimension of scales must be the same as the vector block size
306: Level: advanced
308: .seealso: VecNorm(), VecStrideScale(), VecScale(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
309: @*/
310: PetscErrorCode VecStrideScaleAll(Vec v,const PetscScalar *scales)
311: {
313: PetscInt i,j,n,bs;
314: PetscScalar *x;
319: VecGetLocalSize(v,&n);
320: VecGetBlockSize(v,&bs);
321: VecGetArray(v,&x);
322: /* need to provide optimized code for each bs */
323: for (i=0; i<n; i+=bs) {
324: for (j=0; j<bs; j++) x[i+j] *= scales[j];
325: }
326: VecRestoreArray(v,&x);
327: return(0);
328: }
330: /*@
331: VecStrideNormAll - Computes the norms of subvectors of a vector defined
332: by a starting point and a stride.
334: Collective on Vec
336: Input Parameters:
337: + v - the vector
338: - ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY
340: Output Parameter:
341: . nrm - the norms
343: Notes:
344: One must call VecSetBlockSize() before this routine to set the stride
345: information, or use a vector created from a multicomponent DMDA.
347: If x is the array representing the vector x then this computes the norm
348: of the array (x[start],x[start+stride],x[start+2*stride], ....) for each start < stride
350: The dimension of nrm must be the same as the vector block size
352: This will only work if the desire subvector is a stride subvector
354: Level: advanced
356: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
357: @*/
358: PetscErrorCode VecStrideNormAll(Vec v,NormType ntype,PetscReal nrm[])
359: {
360: PetscErrorCode ierr;
361: PetscInt i,j,n,bs;
362: const PetscScalar *x;
363: PetscReal tnorm[128];
364: MPI_Comm comm;
370: VecGetLocalSize(v,&n);
371: VecGetArrayRead(v,&x);
372: PetscObjectGetComm((PetscObject)v,&comm);
374: VecGetBlockSize(v,&bs);
375: if (bs > 128) SETERRQ(comm,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");
377: if (ntype == NORM_2) {
378: PetscScalar sum[128];
379: for (j=0; j<bs; j++) sum[j] = 0.0;
380: for (i=0; i<n; i+=bs) {
381: for (j=0; j<bs; j++) sum[j] += x[i+j]*(PetscConj(x[i+j]));
382: }
383: for (j=0; j<bs; j++) tnorm[j] = PetscRealPart(sum[j]);
385: MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_SUM,comm);
386: for (j=0; j<bs; j++) nrm[j] = PetscSqrtReal(nrm[j]);
387: } else if (ntype == NORM_1) {
388: for (j=0; j<bs; j++) tnorm[j] = 0.0;
390: for (i=0; i<n; i+=bs) {
391: for (j=0; j<bs; j++) tnorm[j] += PetscAbsScalar(x[i+j]);
392: }
394: MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_SUM,comm);
395: } else if (ntype == NORM_INFINITY) {
396: PetscReal tmp;
397: for (j=0; j<bs; j++) tnorm[j] = 0.0;
399: for (i=0; i<n; i+=bs) {
400: for (j=0; j<bs; j++) {
401: if ((tmp = PetscAbsScalar(x[i+j])) > tnorm[j]) tnorm[j] = tmp;
402: /* check special case of tmp == NaN */
403: if (tmp != tmp) {tnorm[j] = tmp; break;}
404: }
405: }
406: MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_MAX,comm);
407: } else SETERRQ(comm,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
408: VecRestoreArrayRead(v,&x);
409: return(0);
410: }
412: /*@
413: VecStrideMaxAll - Computes the maximums of subvectors of a vector defined
414: by a starting point and a stride and optionally its location.
416: Collective on Vec
418: Input Parameter:
419: . v - the vector
421: Output Parameters:
422: + index - the location where the maximum occurred (not supported, pass NULL,
423: if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
424: - nrm - the maximum values of each subvector
426: Notes:
427: One must call VecSetBlockSize() before this routine to set the stride
428: information, or use a vector created from a multicomponent DMDA.
430: The dimension of nrm must be the same as the vector block size
432: Level: advanced
434: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
435: @*/
436: PetscErrorCode VecStrideMaxAll(Vec v,PetscInt idex[],PetscReal nrm[])
437: {
438: PetscErrorCode ierr;
439: PetscInt i,j,n,bs;
440: const PetscScalar *x;
441: PetscReal max[128],tmp;
442: MPI_Comm comm;
447: 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");
448: VecGetLocalSize(v,&n);
449: VecGetArrayRead(v,&x);
450: PetscObjectGetComm((PetscObject)v,&comm);
452: VecGetBlockSize(v,&bs);
453: if (bs > 128) SETERRQ(comm,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");
455: if (!n) {
456: for (j=0; j<bs; j++) max[j] = PETSC_MIN_REAL;
457: } else {
458: for (j=0; j<bs; j++) max[j] = PetscRealPart(x[j]);
460: for (i=bs; i<n; i+=bs) {
461: for (j=0; j<bs; j++) {
462: if ((tmp = PetscRealPart(x[i+j])) > max[j]) max[j] = tmp;
463: }
464: }
465: }
466: MPIU_Allreduce(max,nrm,bs,MPIU_REAL,MPIU_MAX,comm);
468: VecRestoreArrayRead(v,&x);
469: return(0);
470: }
472: /*@
473: VecStrideMinAll - Computes the minimum of subvector of a vector defined
474: by a starting point and a stride and optionally its location.
476: Collective on Vec
478: Input Parameter:
479: . v - the vector
481: Output Parameters:
482: + idex - the location where the minimum occurred (not supported, pass NULL,
483: if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
484: - nrm - the minimums of each subvector
486: Level: advanced
488: Notes:
489: One must call VecSetBlockSize() before this routine to set the stride
490: information, or use a vector created from a multicomponent DMDA.
492: The dimension of nrm must be the same as the vector block size
494: .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
495: @*/
496: PetscErrorCode VecStrideMinAll(Vec v,PetscInt idex[],PetscReal nrm[])
497: {
498: PetscErrorCode ierr;
499: PetscInt i,n,bs,j;
500: const PetscScalar *x;
501: PetscReal min[128],tmp;
502: MPI_Comm comm;
507: 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");
508: VecGetLocalSize(v,&n);
509: VecGetArrayRead(v,&x);
510: PetscObjectGetComm((PetscObject)v,&comm);
512: VecGetBlockSize(v,&bs);
513: if (bs > 128) SETERRQ(comm,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");
515: if (!n) {
516: for (j=0; j<bs; j++) min[j] = PETSC_MAX_REAL;
517: } else {
518: for (j=0; j<bs; j++) min[j] = PetscRealPart(x[j]);
520: for (i=bs; i<n; i+=bs) {
521: for (j=0; j<bs; j++) {
522: if ((tmp = PetscRealPart(x[i+j])) < min[j]) min[j] = tmp;
523: }
524: }
525: }
526: MPIU_Allreduce(min,nrm,bs,MPIU_REAL,MPIU_MIN,comm);
528: VecRestoreArrayRead(v,&x);
529: return(0);
530: }
532: /*----------------------------------------------------------------------------------------------*/
533: /*@
534: VecStrideGatherAll - Gathers all the single components from a multi-component vector into
535: separate vectors.
537: Collective on Vec
539: Input Parameters:
540: + v - the vector
541: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
543: Output Parameter:
544: . s - the location where the subvectors are stored
546: Notes:
547: One must call VecSetBlockSize() before this routine to set the stride
548: information, or use a vector created from a multicomponent DMDA.
550: If x is the array representing the vector x then this gathers
551: the arrays (x[start],x[start+stride],x[start+2*stride], ....)
552: for start=0,1,2,...bs-1
554: The parallel layout of the vector and the subvector must be the same;
555: i.e., nlocal of v = stride*(nlocal of s)
557: Not optimized; could be easily
559: Level: advanced
561: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
562: VecStrideScatterAll()
563: @*/
564: PetscErrorCode VecStrideGatherAll(Vec v,Vec s[],InsertMode addv)
565: {
566: PetscErrorCode ierr;
567: PetscInt i,n,n2,bs,j,k,*bss = NULL,nv,jj,nvc;
568: PetscScalar **y;
569: const PetscScalar *x;
575: VecGetLocalSize(v,&n);
576: VecGetLocalSize(s[0],&n2);
577: VecGetArrayRead(v,&x);
578: VecGetBlockSize(v,&bs);
579: if (bs <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Input vector does not have a valid blocksize set");
581: PetscMalloc2(bs,&y,bs,&bss);
582: nv = 0;
583: nvc = 0;
584: for (i=0; i<bs; i++) {
585: VecGetBlockSize(s[i],&bss[i]);
586: if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1 Re: [PETSC #8241] VecStrideGatherAll */
587: VecGetArray(s[i],&y[i]);
588: nvc += bss[i];
589: nv++;
590: if (nvc > bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of subvectors in subvectors > number of vectors in main vector");
591: if (nvc == bs) break;
592: }
594: n = n/bs;
596: jj = 0;
597: if (addv == INSERT_VALUES) {
598: for (j=0; j<nv; j++) {
599: for (k=0; k<bss[j]; k++) {
600: for (i=0; i<n; i++) y[j][i*bss[j] + k] = x[bs*i+jj+k];
601: }
602: jj += bss[j];
603: }
604: } else if (addv == ADD_VALUES) {
605: for (j=0; j<nv; j++) {
606: for (k=0; k<bss[j]; k++) {
607: for (i=0; i<n; i++) y[j][i*bss[j] + k] += x[bs*i+jj+k];
608: }
609: jj += bss[j];
610: }
611: #if !defined(PETSC_USE_COMPLEX)
612: } else if (addv == MAX_VALUES) {
613: for (j=0; j<nv; j++) {
614: for (k=0; k<bss[j]; k++) {
615: for (i=0; i<n; i++) y[j][i*bss[j] + k] = PetscMax(y[j][i*bss[j] + k],x[bs*i+jj+k]);
616: }
617: jj += bss[j];
618: }
619: #endif
620: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
622: VecRestoreArrayRead(v,&x);
623: for (i=0; i<nv; i++) {
624: VecRestoreArray(s[i],&y[i]);
625: }
627: PetscFree2(y,bss);
628: return(0);
629: }
631: /*@
632: VecStrideScatterAll - Scatters all the single components from separate vectors into
633: a multi-component vector.
635: Collective on Vec
637: Input Parameters:
638: + s - the location where the subvectors are stored
639: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
641: Output Parameter:
642: . v - the multicomponent vector
644: Notes:
645: One must call VecSetBlockSize() before this routine to set the stride
646: information, or use a vector created from a multicomponent DMDA.
648: The parallel layout of the vector and the subvector must be the same;
649: i.e., nlocal of v = stride*(nlocal of s)
651: Not optimized; could be easily
653: Level: advanced
655: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
656: VecStrideScatterAll()
657: @*/
658: PetscErrorCode VecStrideScatterAll(Vec s[],Vec v,InsertMode addv)
659: {
660: PetscErrorCode ierr;
661: PetscInt i,n,n2,bs,j,jj,k,*bss = NULL,nv,nvc;
662: PetscScalar *x;
663: PetscScalar const **y;
669: VecGetLocalSize(v,&n);
670: VecGetLocalSize(s[0],&n2);
671: VecGetArray(v,&x);
672: VecGetBlockSize(v,&bs);
673: if (bs <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Input vector does not have a valid blocksize set");
675: PetscMalloc2(bs,(PetscScalar***)&y,bs,&bss);
676: nv = 0;
677: nvc = 0;
678: for (i=0; i<bs; i++) {
679: VecGetBlockSize(s[i],&bss[i]);
680: if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1 Re: [PETSC #8241] VecStrideGatherAll */
681: VecGetArrayRead(s[i],&y[i]);
682: nvc += bss[i];
683: nv++;
684: if (nvc > bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of subvectors in subvectors > number of vectors in main vector");
685: if (nvc == bs) break;
686: }
688: n = n/bs;
690: jj = 0;
691: if (addv == INSERT_VALUES) {
692: for (j=0; j<nv; j++) {
693: for (k=0; k<bss[j]; k++) {
694: for (i=0; i<n; i++) x[bs*i+jj+k] = y[j][i*bss[j] + k];
695: }
696: jj += bss[j];
697: }
698: } else if (addv == ADD_VALUES) {
699: for (j=0; j<nv; j++) {
700: for (k=0; k<bss[j]; k++) {
701: for (i=0; i<n; i++) x[bs*i+jj+k] += y[j][i*bss[j] + k];
702: }
703: jj += bss[j];
704: }
705: #if !defined(PETSC_USE_COMPLEX)
706: } else if (addv == MAX_VALUES) {
707: for (j=0; j<nv; j++) {
708: for (k=0; k<bss[j]; k++) {
709: for (i=0; i<n; i++) x[bs*i+jj+k] = PetscMax(x[bs*i+jj+k],y[j][i*bss[j] + k]);
710: }
711: jj += bss[j];
712: }
713: #endif
714: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
716: VecRestoreArray(v,&x);
717: for (i=0; i<nv; i++) {
718: VecRestoreArrayRead(s[i],&y[i]);
719: }
720: PetscFree2(*(PetscScalar***)&y,bss);
721: return(0);
722: }
724: /*@
725: VecStrideGather - Gathers a single component from a multi-component vector into
726: another vector.
728: Collective on Vec
730: Input Parameters:
731: + v - the vector
732: . start - starting point of the subvector (defined by a stride)
733: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
735: Output Parameter:
736: . s - the location where the subvector is stored
738: Notes:
739: One must call VecSetBlockSize() before this routine to set the stride
740: information, or use a vector created from a multicomponent DMDA.
742: If x is the array representing the vector x then this gathers
743: the array (x[start],x[start+stride],x[start+2*stride], ....)
745: The parallel layout of the vector and the subvector must be the same;
746: i.e., nlocal of v = stride*(nlocal of s)
748: Not optimized; could be easily
750: Level: advanced
752: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
753: VecStrideScatterAll()
754: @*/
755: PetscErrorCode VecStrideGather(Vec v,PetscInt start,Vec s,InsertMode addv)
756: {
762: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
763: 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);
764: if (!v->ops->stridegather) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
765: (*v->ops->stridegather)(v,start,s,addv);
766: return(0);
767: }
769: /*@
770: VecStrideScatter - Scatters a single component from a vector into a multi-component vector.
772: Collective on Vec
774: Input Parameters:
775: + s - the single-component vector
776: . start - starting point of the subvector (defined by a stride)
777: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
779: Output Parameter:
780: . v - the location where the subvector is scattered (the multi-component vector)
782: Notes:
783: One must call VecSetBlockSize() on the multi-component vector before this
784: routine to set the stride information, or use a vector created from a multicomponent DMDA.
786: The parallel layout of the vector and the subvector must be the same;
787: i.e., nlocal of v = stride*(nlocal of s)
789: Not optimized; could be easily
791: Level: advanced
793: .seealso: VecStrideNorm(), VecStrideGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
794: VecStrideScatterAll(), VecStrideSubSetScatter(), VecStrideSubSetGather()
795: @*/
796: PetscErrorCode VecStrideScatter(Vec s,PetscInt start,Vec v,InsertMode addv)
797: {
803: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
804: 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);
805: if (!v->ops->stridescatter) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
806: (*v->ops->stridescatter)(s,start,v,addv);
807: return(0);
808: }
810: /*@
811: VecStrideSubSetGather - Gathers a subset of components from a multi-component vector into
812: another vector.
814: Collective on Vec
816: Input Parameters:
817: + v - the vector
818: . nidx - the number of indices
819: . idxv - the indices of the components 0 <= idxv[0] ...idxv[nidx-1] < bs(v), they need not be sorted
820: . 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
821: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
823: Output Parameter:
824: . s - the location where the subvector is stored
826: Notes:
827: One must call VecSetBlockSize() on both vectors before this routine to set the stride
828: information, or use a vector created from a multicomponent DMDA.
830: The parallel layout of the vector and the subvector must be the same;
832: Not optimized; could be easily
834: Level: advanced
836: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideGather(), VecStrideSubSetScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
837: VecStrideScatterAll()
838: @*/
839: PetscErrorCode VecStrideSubSetGather(Vec v,PetscInt nidx,const PetscInt idxv[],const PetscInt idxs[],Vec s,InsertMode addv)
840: {
846: if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
847: if (!v->ops->stridesubsetgather) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
848: (*v->ops->stridesubsetgather)(v,nidx,idxv,idxs,s,addv);
849: return(0);
850: }
852: /*@
853: VecStrideSubSetScatter - Scatters components from a vector into a subset of components of a multi-component vector.
855: Collective on Vec
857: Input Parameters:
858: + s - the smaller-component vector
859: . nidx - the number of indices in idx
860: . 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
861: . idxv - the indices of the components in the larger-component vector, 0 <= idx[0] ...idx[nidx-1] < bs(v) they need not be sorted
862: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
864: Output Parameter:
865: . v - the location where the subvector is into scattered (the multi-component vector)
867: Notes:
868: One must call VecSetBlockSize() on the vectors before this
869: routine to set the stride information, or use a vector created from a multicomponent DMDA.
871: The parallel layout of the vector and the subvector must be the same;
873: Not optimized; could be easily
875: Level: advanced
877: .seealso: VecStrideNorm(), VecStrideGather(), VecStrideGather(), VecStrideSubSetGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
878: VecStrideScatterAll()
879: @*/
880: PetscErrorCode VecStrideSubSetScatter(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)
881: {
887: if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
888: if (!v->ops->stridesubsetscatter) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
889: (*v->ops->stridesubsetscatter)(s,nidx,idxs,idxv,v,addv);
890: return(0);
891: }
893: PetscErrorCode VecStrideGather_Default(Vec v,PetscInt start,Vec s,InsertMode addv)
894: {
896: PetscInt i,n,bs,ns;
897: const PetscScalar *x;
898: PetscScalar *y;
901: VecGetLocalSize(v,&n);
902: VecGetLocalSize(s,&ns);
903: VecGetArrayRead(v,&x);
904: VecGetArray(s,&y);
906: bs = v->map->bs;
907: 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);
908: x += start;
909: n = n/bs;
911: if (addv == INSERT_VALUES) {
912: for (i=0; i<n; i++) y[i] = x[bs*i];
913: } else if (addv == ADD_VALUES) {
914: for (i=0; i<n; i++) y[i] += x[bs*i];
915: #if !defined(PETSC_USE_COMPLEX)
916: } else if (addv == MAX_VALUES) {
917: for (i=0; i<n; i++) y[i] = PetscMax(y[i],x[bs*i]);
918: #endif
919: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
921: VecRestoreArrayRead(v,&x);
922: VecRestoreArray(s,&y);
923: return(0);
924: }
926: PetscErrorCode VecStrideScatter_Default(Vec s,PetscInt start,Vec v,InsertMode addv)
927: {
928: PetscErrorCode ierr;
929: PetscInt i,n,bs,ns;
930: PetscScalar *x;
931: const PetscScalar *y;
934: VecGetLocalSize(v,&n);
935: VecGetLocalSize(s,&ns);
936: VecGetArray(v,&x);
937: VecGetArrayRead(s,&y);
939: bs = v->map->bs;
940: 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);
941: x += start;
942: n = n/bs;
944: if (addv == INSERT_VALUES) {
945: for (i=0; i<n; i++) x[bs*i] = y[i];
946: } else if (addv == ADD_VALUES) {
947: for (i=0; i<n; i++) x[bs*i] += y[i];
948: #if !defined(PETSC_USE_COMPLEX)
949: } else if (addv == MAX_VALUES) {
950: for (i=0; i<n; i++) x[bs*i] = PetscMax(y[i],x[bs*i]);
951: #endif
952: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
954: VecRestoreArray(v,&x);
955: VecRestoreArrayRead(s,&y);
956: return(0);
957: }
959: PetscErrorCode VecStrideSubSetGather_Default(Vec v,PetscInt nidx,const PetscInt idxv[],const PetscInt idxs[],Vec s,InsertMode addv)
960: {
961: PetscErrorCode ierr;
962: PetscInt i,j,n,bs,bss,ns;
963: const PetscScalar *x;
964: PetscScalar *y;
967: VecGetLocalSize(v,&n);
968: VecGetLocalSize(s,&ns);
969: VecGetArrayRead(v,&x);
970: VecGetArray(s,&y);
972: bs = v->map->bs;
973: bss = s->map->bs;
974: n = n/bs;
976: if (PetscDefined(USE_DEBUG)) {
977: if (n != ns/bss) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible layout of vectors");
978: for (j=0; j<nidx; j++) {
979: if (idxv[j] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is negative",j,idxv[j]);
980: if (idxv[j] >= bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is greater than or equal to vector blocksize %D",j,idxv[j],bs);
981: }
982: if (!idxs && bss != nidx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide idxs when not gathering into all locations");
983: }
985: if (addv == INSERT_VALUES) {
986: if (!idxs) {
987: for (i=0; i<n; i++) {
988: for (j=0; j<bss; j++) y[bss*i+j] = x[bs*i+idxv[j]];
989: }
990: } else {
991: for (i=0; i<n; i++) {
992: for (j=0; j<bss; j++) y[bss*i+idxs[j]] = x[bs*i+idxv[j]];
993: }
994: }
995: } else if (addv == ADD_VALUES) {
996: if (!idxs) {
997: for (i=0; i<n; i++) {
998: for (j=0; j<bss; j++) y[bss*i+j] += x[bs*i+idxv[j]];
999: }
1000: } else {
1001: for (i=0; i<n; i++) {
1002: for (j=0; j<bss; j++) y[bss*i+idxs[j]] += x[bs*i+idxv[j]];
1003: }
1004: }
1005: #if !defined(PETSC_USE_COMPLEX)
1006: } else if (addv == MAX_VALUES) {
1007: if (!idxs) {
1008: for (i=0; i<n; i++) {
1009: for (j=0; j<bss; j++) y[bss*i+j] = PetscMax(y[bss*i+j],x[bs*i+idxv[j]]);
1010: }
1011: } else {
1012: for (i=0; i<n; i++) {
1013: for (j=0; j<bss; j++) y[bss*i+idxs[j]] = PetscMax(y[bss*i+idxs[j]],x[bs*i+idxv[j]]);
1014: }
1015: }
1016: #endif
1017: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
1019: VecRestoreArrayRead(v,&x);
1020: VecRestoreArray(s,&y);
1021: return(0);
1022: }
1024: PetscErrorCode VecStrideSubSetScatter_Default(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)
1025: {
1026: PetscErrorCode ierr;
1027: PetscInt j,i,n,bs,ns,bss;
1028: PetscScalar *x;
1029: const PetscScalar *y;
1032: VecGetLocalSize(v,&n);
1033: VecGetLocalSize(s,&ns);
1034: VecGetArray(v,&x);
1035: VecGetArrayRead(s,&y);
1037: bs = v->map->bs;
1038: bss = s->map->bs;
1039: n = n/bs;
1041: if (PetscDefined(USE_DEBUG)) {
1042: if (n != ns/bss) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible layout of vectors");
1043: for (j=0; j<bss; j++) {
1044: if (idxs) {
1045: if (idxs[j] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is negative",j,idxs[j]);
1046: if (idxs[j] >= bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is greater than or equal to vector blocksize %D",j,idxs[j],bs);
1047: }
1048: }
1049: if (!idxs && bss != nidx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide idxs when not scattering from all locations");
1050: }
1052: if (addv == INSERT_VALUES) {
1053: if (!idxs) {
1054: for (i=0; i<n; i++) {
1055: for (j=0; j<bss; j++) x[bs*i + idxv[j]] = y[bss*i+j];
1056: }
1057: } else {
1058: for (i=0; i<n; i++) {
1059: for (j=0; j<bss; j++) x[bs*i + idxv[j]] = y[bss*i+idxs[j]];
1060: }
1061: }
1062: } else if (addv == ADD_VALUES) {
1063: if (!idxs) {
1064: for (i=0; i<n; i++) {
1065: for (j=0; j<bss; j++) x[bs*i + idxv[j]] += y[bss*i+j];
1066: }
1067: } else {
1068: for (i=0; i<n; i++) {
1069: for (j=0; j<bss; j++) x[bs*i + idxv[j]] += y[bss*i+idxs[j]];
1070: }
1071: }
1072: #if !defined(PETSC_USE_COMPLEX)
1073: } else if (addv == MAX_VALUES) {
1074: if (!idxs) {
1075: for (i=0; i<n; i++) {
1076: for (j=0; j<bss; j++) x[bs*i + idxv[j]] = PetscMax(y[bss*i+j],x[bs*i + idxv[j]]);
1077: }
1078: } else {
1079: for (i=0; i<n; i++) {
1080: for (j=0; j<bss; j++) x[bs*i + idxv[j]] = PetscMax(y[bss*i+idxs[j]],x[bs*i + idxv[j]]);
1081: }
1082: }
1083: #endif
1084: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
1086: VecRestoreArray(v,&x);
1087: VecRestoreArrayRead(s,&y);
1088: return(0);
1089: }
1091: PetscErrorCode VecReciprocal_Default(Vec v)
1092: {
1094: PetscInt i,n;
1095: PetscScalar *x;
1098: VecGetLocalSize(v,&n);
1099: VecGetArray(v,&x);
1100: for (i=0; i<n; i++) {
1101: if (x[i] != (PetscScalar)0.0) x[i] = (PetscScalar)1.0/x[i];
1102: }
1103: VecRestoreArray(v,&x);
1104: return(0);
1105: }
1107: /*@
1108: VecExp - Replaces each component of a vector by e^x_i
1110: Not collective
1112: Input Parameter:
1113: . v - The vector
1115: Output Parameter:
1116: . v - The vector of exponents
1118: Level: beginner
1120: .seealso: VecLog(), VecAbs(), VecSqrtAbs(), VecReciprocal()
1122: @*/
1123: PetscErrorCode VecExp(Vec v)
1124: {
1125: PetscScalar *x;
1126: PetscInt i, n;
1131: if (v->ops->exp) {
1132: (*v->ops->exp)(v);
1133: } else {
1134: VecGetLocalSize(v, &n);
1135: VecGetArray(v, &x);
1136: for (i = 0; i < n; i++) x[i] = PetscExpScalar(x[i]);
1137: VecRestoreArray(v, &x);
1138: }
1139: return(0);
1140: }
1142: /*@
1143: VecLog - Replaces each component of a vector by log(x_i), the natural logarithm
1145: Not collective
1147: Input Parameter:
1148: . v - The vector
1150: Output Parameter:
1151: . v - The vector of logs
1153: Level: beginner
1155: .seealso: VecExp(), VecAbs(), VecSqrtAbs(), VecReciprocal()
1157: @*/
1158: PetscErrorCode VecLog(Vec v)
1159: {
1160: PetscScalar *x;
1161: PetscInt i, n;
1166: if (v->ops->log) {
1167: (*v->ops->log)(v);
1168: } else {
1169: VecGetLocalSize(v, &n);
1170: VecGetArray(v, &x);
1171: for (i = 0; i < n; i++) x[i] = PetscLogScalar(x[i]);
1172: VecRestoreArray(v, &x);
1173: }
1174: return(0);
1175: }
1177: /*@
1178: VecSqrtAbs - Replaces each component of a vector by the square root of its magnitude.
1180: Not collective
1182: Input Parameter:
1183: . v - The vector
1185: Output Parameter:
1186: . v - The vector square root
1188: Level: beginner
1190: Note: The actual function is sqrt(|x_i|)
1192: .seealso: VecLog(), VecExp(), VecReciprocal(), VecAbs()
1194: @*/
1195: PetscErrorCode VecSqrtAbs(Vec v)
1196: {
1197: PetscScalar *x;
1198: PetscInt i, n;
1203: if (v->ops->sqrt) {
1204: (*v->ops->sqrt)(v);
1205: } else {
1206: VecGetLocalSize(v, &n);
1207: VecGetArray(v, &x);
1208: for (i = 0; i < n; i++) x[i] = PetscSqrtReal(PetscAbsScalar(x[i]));
1209: VecRestoreArray(v, &x);
1210: }
1211: return(0);
1212: }
1214: /*@
1215: VecDotNorm2 - computes the inner product of two vectors and the 2-norm squared of the second vector
1217: Collective on Vec
1219: Input Parameters:
1220: + s - first vector
1221: - t - second vector
1223: Output Parameters:
1224: + dp - s'conj(t)
1225: - nm - t'conj(t)
1227: Level: advanced
1229: Notes:
1230: conj(x) is the complex conjugate of x when x is complex
1232: .seealso: VecDot(), VecNorm(), VecDotBegin(), VecNormBegin(), VecDotEnd(), VecNormEnd()
1234: @*/
1235: PetscErrorCode VecDotNorm2(Vec s,Vec t,PetscScalar *dp, PetscReal *nm)
1236: {
1237: const PetscScalar *sx, *tx;
1238: PetscScalar dpx = 0.0, nmx = 0.0,work[2],sum[2];
1239: PetscInt i, n;
1240: PetscErrorCode ierr;
1250: if (s->map->N != t->map->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
1251: if (s->map->n != t->map->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");
1253: PetscLogEventBegin(VEC_DotNorm2,s,t,0,0);
1254: if (s->ops->dotnorm2) {
1255: (*s->ops->dotnorm2)(s,t,dp,&dpx);
1256: *nm = PetscRealPart(dpx);
1257: } else {
1258: VecGetLocalSize(s, &n);
1259: VecGetArrayRead(s, &sx);
1260: VecGetArrayRead(t, &tx);
1262: for (i = 0; i<n; i++) {
1263: dpx += sx[i]*PetscConj(tx[i]);
1264: nmx += tx[i]*PetscConj(tx[i]);
1265: }
1266: work[0] = dpx;
1267: work[1] = nmx;
1269: MPIU_Allreduce(work,sum,2,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)s));
1270: *dp = sum[0];
1271: *nm = PetscRealPart(sum[1]);
1273: VecRestoreArrayRead(t, &tx);
1274: VecRestoreArrayRead(s, &sx);
1275: PetscLogFlops(4.0*n);
1276: }
1277: PetscLogEventEnd(VEC_DotNorm2,s,t,0,0);
1278: return(0);
1279: }
1281: /*@
1282: VecSum - Computes the sum of all the components of a vector.
1284: Collective on Vec
1286: Input Parameter:
1287: . v - the vector
1289: Output Parameter:
1290: . sum - the result
1292: Level: beginner
1294: .seealso: VecMean(), VecNorm()
1295: @*/
1296: PetscErrorCode VecSum(Vec v,PetscScalar *sum)
1297: {
1298: PetscErrorCode ierr;
1299: PetscInt i,n;
1300: const PetscScalar *x;
1305: *sum = 0.0;
1306: if (v->ops->sum) {
1307: (*v->ops->sum)(v,sum);
1308: } else {
1309: VecGetLocalSize(v,&n);
1310: VecGetArrayRead(v,&x);
1311: for (i=0; i<n; i++) *sum += x[i];
1312: VecRestoreArrayRead(v,&x);
1313: }
1314: MPIU_Allreduce(MPI_IN_PLACE,sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)v));
1315: return(0);
1316: }
1318: /*@
1319: VecMean - Computes the arithmetic mean of all the components of a vector.
1321: Collective on Vec
1323: Input Parameter:
1324: . v - the vector
1326: Output Parameter:
1327: . mean - the result
1329: Level: beginner
1331: .seealso: VecSum(), VecNorm()
1332: @*/
1333: PetscErrorCode VecMean(Vec v,PetscScalar *mean)
1334: {
1335: PetscErrorCode ierr;
1336: PetscInt n;
1341: VecGetSize(v,&n);
1342: VecSum(v,mean);
1343: *mean /= n;
1344: return(0);
1345: }
1347: /*@
1348: VecImaginaryPart - Replaces a complex vector with its imginary part
1350: Collective on Vec
1352: Input Parameter:
1353: . v - the vector
1355: Level: beginner
1357: .seealso: VecNorm(), VecRealPart()
1358: @*/
1359: PetscErrorCode VecImaginaryPart(Vec v)
1360: {
1361: PetscErrorCode ierr;
1362: PetscInt i,n;
1363: PetscScalar *x;
1367: VecGetLocalSize(v,&n);
1368: VecGetArray(v,&x);
1369: for (i=0; i<n; i++) x[i] = PetscImaginaryPart(x[i]);
1370: VecRestoreArray(v,&x);
1371: return(0);
1372: }
1374: /*@
1375: VecRealPart - Replaces a complex vector with its real part
1377: Collective on Vec
1379: Input Parameter:
1380: . v - the vector
1382: Level: beginner
1384: .seealso: VecNorm(), VecImaginaryPart()
1385: @*/
1386: PetscErrorCode VecRealPart(Vec v)
1387: {
1388: PetscErrorCode ierr;
1389: PetscInt i,n;
1390: PetscScalar *x;
1394: VecGetLocalSize(v,&n);
1395: VecGetArray(v,&x);
1396: for (i=0; i<n; i++) x[i] = PetscRealPart(x[i]);
1397: VecRestoreArray(v,&x);
1398: return(0);
1399: }
1401: /*@
1402: VecShift - Shifts all of the components of a vector by computing
1403: x[i] = x[i] + shift.
1405: Logically Collective on Vec
1407: Input Parameters:
1408: + v - the vector
1409: - shift - the shift
1411: Level: intermediate
1413: @*/
1414: PetscErrorCode VecShift(Vec v,PetscScalar shift)
1415: {
1417: PetscInt i,n;
1418: PetscScalar *x;
1423: VecSetErrorIfLocked(v,1);
1424: if (shift == 0.0) return(0);
1426: if (v->ops->shift) {
1427: (*v->ops->shift)(v,shift);
1428: } else {
1429: VecGetLocalSize(v,&n);
1430: VecGetArray(v,&x);
1431: for (i=0; i<n; i++) x[i] += shift;
1432: VecRestoreArray(v,&x);
1433: }
1434: return(0);
1435: }
1437: /*@
1438: VecAbs - Replaces every element in a vector with its absolute value.
1440: Logically Collective on Vec
1442: Input Parameters:
1443: . v - the vector
1445: Level: intermediate
1447: @*/
1448: PetscErrorCode VecAbs(Vec v)
1449: {
1451: PetscInt i,n;
1452: PetscScalar *x;
1456: VecSetErrorIfLocked(v,1);
1458: if (v->ops->abs) {
1459: (*v->ops->abs)(v);
1460: } else {
1461: VecGetLocalSize(v,&n);
1462: VecGetArray(v,&x);
1463: for (i=0; i<n; i++) x[i] = PetscAbsScalar(x[i]);
1464: VecRestoreArray(v,&x);
1465: }
1466: return(0);
1467: }
1469: /*@
1470: VecPermute - Permutes a vector in place using the given ordering.
1472: Input Parameters:
1473: + vec - The vector
1474: . order - The ordering
1475: - inv - The flag for inverting the permutation
1477: Level: beginner
1479: Note: This function does not yet support parallel Index Sets with non-local permutations
1481: .seealso: MatPermute()
1482: @*/
1483: PetscErrorCode VecPermute(Vec x, IS row, PetscBool inv)
1484: {
1485: const PetscScalar *array;
1486: PetscScalar *newArray;
1487: const PetscInt *idx;
1488: PetscInt i,rstart,rend;
1489: PetscErrorCode ierr;
1494: VecSetErrorIfLocked(x,1);
1495: VecGetOwnershipRange(x,&rstart,&rend);
1496: ISGetIndices(row, &idx);
1497: VecGetArrayRead(x, &array);
1498: PetscMalloc1(x->map->n, &newArray);
1499: if (PetscDefined(USE_DEBUG)) {
1500: for (i = 0; i < x->map->n; i++) {
1501: 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]);
1502: }
1503: }
1504: if (!inv) {
1505: for (i = 0; i < x->map->n; i++) newArray[i] = array[idx[i]-rstart];
1506: } else {
1507: for (i = 0; i < x->map->n; i++) newArray[idx[i]-rstart] = array[i];
1508: }
1509: VecRestoreArrayRead(x, &array);
1510: ISRestoreIndices(row, &idx);
1511: VecReplaceArray(x, newArray);
1512: return(0);
1513: }
1515: /*@
1516: VecEqual - Compares two vectors. Returns true if the two vectors are either pointing to the same memory buffer,
1517: or if the two vectors have the same local and global layout as well as bitwise equality of all entries.
1518: Does NOT take round-off errors into account.
1520: Collective on Vec
1522: Input Parameters:
1523: + vec1 - the first vector
1524: - vec2 - the second vector
1526: Output Parameter:
1527: . flg - PETSC_TRUE if the vectors are equal; PETSC_FALSE otherwise.
1529: Level: intermediate
1530: @*/
1531: PetscErrorCode VecEqual(Vec vec1,Vec vec2,PetscBool *flg)
1532: {
1533: const PetscScalar *v1,*v2;
1534: PetscErrorCode ierr;
1535: PetscInt n1,n2,N1,N2;
1536: PetscBool flg1;
1542: if (vec1 == vec2) *flg = PETSC_TRUE;
1543: else {
1544: VecGetSize(vec1,&N1);
1545: VecGetSize(vec2,&N2);
1546: if (N1 != N2) flg1 = PETSC_FALSE;
1547: else {
1548: VecGetLocalSize(vec1,&n1);
1549: VecGetLocalSize(vec2,&n2);
1550: if (n1 != n2) flg1 = PETSC_FALSE;
1551: else {
1552: VecGetArrayRead(vec1,&v1);
1553: VecGetArrayRead(vec2,&v2);
1554: PetscArraycmp(v1,v2,n1,&flg1);
1555: VecRestoreArrayRead(vec1,&v1);
1556: VecRestoreArrayRead(vec2,&v2);
1557: }
1558: }
1559: /* combine results from all processors */
1560: MPIU_Allreduce(&flg1,flg,1,MPIU_BOOL,MPI_MIN,PetscObjectComm((PetscObject)vec1));
1561: }
1562: return(0);
1563: }
1565: /*@
1566: VecUniqueEntries - Compute the number of unique entries, and those entries
1568: Collective on Vec
1570: Input Parameter:
1571: . vec - the vector
1573: Output Parameters:
1574: + n - The number of unique entries
1575: - e - The entries
1577: Level: intermediate
1579: @*/
1580: PetscErrorCode VecUniqueEntries(Vec vec, PetscInt *n, PetscScalar **e)
1581: {
1582: const PetscScalar *v;
1583: PetscScalar *tmp, *vals;
1584: PetscMPIInt *N, *displs, l;
1585: PetscInt ng, m, i, j, p;
1586: PetscMPIInt size;
1587: PetscErrorCode ierr;
1592: MPI_Comm_size(PetscObjectComm((PetscObject) vec), &size);
1593: VecGetLocalSize(vec, &m);
1594: VecGetArrayRead(vec, &v);
1595: PetscMalloc2(m,&tmp,size,&N);
1596: for (i = 0, j = 0, l = 0; i < m; ++i) {
1597: /* Can speed this up with sorting */
1598: for (j = 0; j < l; ++j) {
1599: if (v[i] == tmp[j]) break;
1600: }
1601: if (j == l) {
1602: tmp[j] = v[i];
1603: ++l;
1604: }
1605: }
1606: VecRestoreArrayRead(vec, &v);
1607: /* Gather serial results */
1608: MPI_Allgather(&l, 1, MPI_INT, N, 1, MPI_INT, PetscObjectComm((PetscObject) vec));
1609: for (p = 0, ng = 0; p < size; ++p) {
1610: ng += N[p];
1611: }
1612: PetscMalloc2(ng,&vals,size+1,&displs);
1613: for (p = 1, displs[0] = 0; p <= size; ++p) {
1614: displs[p] = displs[p-1] + N[p-1];
1615: }
1616: MPI_Allgatherv(tmp, l, MPIU_SCALAR, vals, N, displs, MPIU_SCALAR, PetscObjectComm((PetscObject) vec));
1617: /* Find unique entries */
1618: #ifdef PETSC_USE_COMPLEX
1619: SETERRQ(PetscObjectComm((PetscObject) vec), PETSC_ERR_SUP, "Does not work with complex numbers");
1620: #else
1621: *n = displs[size];
1622: PetscSortRemoveDupsReal(n, (PetscReal *) vals);
1623: if (e) {
1625: PetscMalloc1(*n, e);
1626: for (i = 0; i < *n; ++i) {
1627: (*e)[i] = vals[i];
1628: }
1629: }
1630: PetscFree2(vals,displs);
1631: PetscFree2(tmp,N);
1632: return(0);
1633: #endif
1634: }