Actual source code: vinv.c
petsc-master 2021-01-21
2: /*
3: Some useful vector utility functions.
4: */
5: #include <../src/vec/vec/impls/mpi/pvecimpl.h>
6: extern MPI_Op MPIU_MAXINDEX_OP;
7: extern MPI_Op MPIU_MININDEX_OP;
9: /*@
10: VecStrideSet - Sets a subvector of a vector defined
11: by a starting point and a stride with a given value
13: Logically Collective on Vec
15: Input Parameter:
16: + v - the vector
17: . start - starting point of the subvector (defined by a stride)
18: - s - value to set for each entry in that subvector
20: Notes:
21: One must call VecSetBlockSize() before this routine to set the stride
22: information, or use a vector created from a multicomponent DMDA.
24: This will only work if the desire subvector is a stride subvector
26: Level: advanced
29: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideScale()
30: @*/
31: PetscErrorCode VecStrideSet(Vec v,PetscInt start,PetscScalar s)
32: {
34: PetscInt i,n,bs;
35: PetscScalar *x;
41: VecSetErrorIfLocked(v,1);
43: VecGetLocalSize(v,&n);
44: VecGetArray(v,&x);
45: VecGetBlockSize(v,&bs);
46: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
47: 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);
48: x += start;
49: for (i=0; i<n; i+=bs) x[i] = s;
50: x -= start;
51: VecRestoreArray(v,&x);
52: return(0);
53: }
55: /*@
56: VecStrideScale - Scales a subvector of a vector defined
57: by a starting point and a stride.
59: Logically Collective on Vec
61: Input Parameter:
62: + v - the vector
63: . start - starting point of the subvector (defined by a stride)
64: - scale - value to multiply each subvector entry by
66: Notes:
67: One must call VecSetBlockSize() before this routine to set the stride
68: information, or use a vector created from a multicomponent DMDA.
70: This will only work if the desire subvector is a stride subvector
72: Level: advanced
75: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideScale()
76: @*/
77: PetscErrorCode VecStrideScale(Vec v,PetscInt start,PetscScalar scale)
78: {
80: PetscInt i,n,bs;
81: PetscScalar *x;
87: VecSetErrorIfLocked(v,1);
89: VecGetLocalSize(v,&n);
90: VecGetArray(v,&x);
91: VecGetBlockSize(v,&bs);
92: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
93: 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);
94: x += start;
95: for (i=0; i<n; i+=bs) x[i] *= scale;
96: x -= start;
97: VecRestoreArray(v,&x);
98: return(0);
99: }
101: /*@
102: VecStrideNorm - Computes the norm of subvector of a vector defined
103: by a starting point and a stride.
105: Collective on Vec
107: Input Parameter:
108: + v - the vector
109: . start - starting point of the subvector (defined by a stride)
110: - ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY
112: Output Parameter:
113: . norm - the norm
115: Notes:
116: One must call VecSetBlockSize() before this routine to set the stride
117: information, or use a vector created from a multicomponent DMDA.
119: If x is the array representing the vector x then this computes the norm
120: of the array (x[start],x[start+stride],x[start+2*stride], ....)
122: This is useful for computing, say the norm of the pressure variable when
123: the pressure is stored (interlaced) with other variables, say density etc.
125: This will only work if the desire subvector is a stride subvector
127: Level: advanced
130: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
131: @*/
132: PetscErrorCode VecStrideNorm(Vec v,PetscInt start,NormType ntype,PetscReal *nrm)
133: {
134: PetscErrorCode ierr;
135: PetscInt i,n,bs;
136: const PetscScalar *x;
137: PetscReal tnorm;
138: MPI_Comm comm;
143: VecGetLocalSize(v,&n);
144: VecGetArrayRead(v,&x);
145: PetscObjectGetComm((PetscObject)v,&comm);
147: VecGetBlockSize(v,&bs);
148: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
149: 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);
150: x += start;
152: if (ntype == NORM_2) {
153: PetscScalar sum = 0.0;
154: for (i=0; i<n; i+=bs) sum += x[i]*(PetscConj(x[i]));
155: tnorm = PetscRealPart(sum);
156: MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_SUM,comm);
157: *nrm = PetscSqrtReal(*nrm);
158: } else if (ntype == NORM_1) {
159: tnorm = 0.0;
160: for (i=0; i<n; i+=bs) tnorm += PetscAbsScalar(x[i]);
161: MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_SUM,comm);
162: } else if (ntype == NORM_INFINITY) {
163: PetscReal tmp;
164: tnorm = 0.0;
166: for (i=0; i<n; i+=bs) {
167: if ((tmp = PetscAbsScalar(x[i])) > tnorm) tnorm = tmp;
168: /* check special case of tmp == NaN */
169: if (tmp != tmp) {tnorm = tmp; break;}
170: }
171: MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_MAX,comm);
172: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
173: VecRestoreArrayRead(v,&x);
174: return(0);
175: }
177: /*@
178: VecStrideMax - Computes the maximum of subvector of a vector defined
179: by a starting point and a stride and optionally its location.
181: Collective on Vec
183: Input Parameter:
184: + v - the vector
185: - start - starting point of the subvector (defined by a stride)
187: Output Parameter:
188: + index - the location where the maximum occurred (pass NULL if not required)
189: - nrm - the maximum value in the subvector
191: Notes:
192: One must call VecSetBlockSize() before this routine to set the stride
193: information, or use a vector created from a multicomponent DMDA.
195: If xa is the array representing the vector x, then this computes the max
196: of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)
198: This is useful for computing, say the maximum of the pressure variable when
199: the pressure is stored (interlaced) with other variables, e.g., density, etc.
200: This will only work if the desire subvector is a stride subvector.
202: Level: advanced
205: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
206: @*/
207: PetscErrorCode VecStrideMax(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
208: {
209: PetscErrorCode ierr;
210: PetscInt i,n,bs,id;
211: const PetscScalar *x;
212: PetscReal max,tmp;
213: MPI_Comm comm;
219: VecGetLocalSize(v,&n);
220: VecGetArrayRead(v,&x);
221: PetscObjectGetComm((PetscObject)v,&comm);
223: VecGetBlockSize(v,&bs);
224: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
225: 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);
226: x += start;
228: id = -1;
229: if (!n) max = PETSC_MIN_REAL;
230: else {
231: id = 0;
232: max = PetscRealPart(x[0]);
233: for (i=bs; i<n; i+=bs) {
234: if ((tmp = PetscRealPart(x[i])) > max) { max = tmp; id = i;}
235: }
236: }
237: VecRestoreArrayRead(v,&x);
239: if (!idex) {
240: MPIU_Allreduce(&max,nrm,1,MPIU_REAL,MPIU_MAX,comm);
241: } else {
242: PetscReal in[2],out[2];
243: PetscInt rstart;
245: VecGetOwnershipRange(v,&rstart,NULL);
246: in[0] = max;
247: in[1] = rstart+id+start;
248: MPIU_Allreduce(in,out,2,MPIU_REAL,MPIU_MAXINDEX_OP,PetscObjectComm((PetscObject)v));
249: *nrm = out[0];
250: *idex = (PetscInt)out[1];
251: }
252: return(0);
253: }
255: /*@
256: VecStrideMin - Computes the minimum of subvector of a vector defined
257: by a starting point and a stride and optionally its location.
259: Collective on Vec
261: Input Parameter:
262: + v - the vector
263: - start - starting point of the subvector (defined by a stride)
265: Output Parameter:
266: + idex - the location where the minimum occurred. (pass NULL if not required)
267: - nrm - the minimum value in the subvector
269: Level: advanced
271: Notes:
272: One must call VecSetBlockSize() before this routine to set the stride
273: information, or use a vector created from a multicomponent DMDA.
275: If xa is the array representing the vector x, then this computes the min
276: of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)
278: This is useful for computing, say the minimum of the pressure variable when
279: the pressure is stored (interlaced) with other variables, e.g., density, etc.
280: This will only work if the desire subvector is a stride subvector.
283: .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
284: @*/
285: PetscErrorCode VecStrideMin(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
286: {
287: PetscErrorCode ierr;
288: PetscInt i,n,bs,id;
289: const PetscScalar *x;
290: PetscReal min,tmp;
291: MPI_Comm comm;
297: VecGetLocalSize(v,&n);
298: VecGetArrayRead(v,&x);
299: PetscObjectGetComm((PetscObject)v,&comm);
301: VecGetBlockSize(v,&bs);
302: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
303: 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);
304: x += start;
306: id = -1;
307: if (!n) min = PETSC_MAX_REAL;
308: else {
309: id = 0;
310: min = PetscRealPart(x[0]);
311: for (i=bs; i<n; i+=bs) {
312: if ((tmp = PetscRealPart(x[i])) < min) { min = tmp; id = i;}
313: }
314: }
315: VecRestoreArrayRead(v,&x);
317: if (!idex) {
318: MPIU_Allreduce(&min,nrm,1,MPIU_REAL,MPIU_MIN,comm);
319: } else {
320: PetscReal in[2],out[2];
321: PetscInt rstart;
323: VecGetOwnershipRange(v,&rstart,NULL);
324: in[0] = min;
325: in[1] = rstart+id;
326: MPIU_Allreduce(in,out,2,MPIU_REAL,MPIU_MININDEX_OP,PetscObjectComm((PetscObject)v));
327: *nrm = out[0];
328: *idex = (PetscInt)out[1];
329: }
330: return(0);
331: }
333: /*@
334: VecStrideScaleAll - Scales the subvectors of a vector defined
335: by a starting point and a stride.
337: Logically Collective on Vec
339: Input Parameter:
340: + v - the vector
341: - scales - values to multiply each subvector entry by
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: The dimension of scales must be the same as the vector block size
350: Level: advanced
353: .seealso: VecNorm(), VecStrideScale(), VecScale(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
354: @*/
355: PetscErrorCode VecStrideScaleAll(Vec v,const PetscScalar *scales)
356: {
358: PetscInt i,j,n,bs;
359: PetscScalar *x;
364: VecSetErrorIfLocked(v,1);
366: VecGetLocalSize(v,&n);
367: VecGetArray(v,&x);
368: VecGetBlockSize(v,&bs);
370: /* need to provide optimized code for each bs */
371: for (i=0; i<n; i+=bs) {
372: for (j=0; j<bs; j++) x[i+j] *= scales[j];
373: }
374: VecRestoreArray(v,&x);
375: return(0);
376: }
378: /*@
379: VecStrideNormAll - Computes the norms of subvectors of a vector defined
380: by a starting point and a stride.
382: Collective on Vec
384: Input Parameter:
385: + v - the vector
386: - ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY
388: Output Parameter:
389: . nrm - the norms
391: Notes:
392: One must call VecSetBlockSize() before this routine to set the stride
393: information, or use a vector created from a multicomponent DMDA.
395: If x is the array representing the vector x then this computes the norm
396: of the array (x[start],x[start+stride],x[start+2*stride], ....) for each start < stride
398: The dimension of nrm must be the same as the vector block size
400: This will only work if the desire subvector is a stride subvector
402: Level: advanced
405: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
406: @*/
407: PetscErrorCode VecStrideNormAll(Vec v,NormType ntype,PetscReal nrm[])
408: {
409: PetscErrorCode ierr;
410: PetscInt i,j,n,bs;
411: const PetscScalar *x;
412: PetscReal tnorm[128];
413: MPI_Comm comm;
418: VecGetLocalSize(v,&n);
419: VecGetArrayRead(v,&x);
420: PetscObjectGetComm((PetscObject)v,&comm);
422: VecGetBlockSize(v,&bs);
423: if (bs > 128) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");
425: if (ntype == NORM_2) {
426: PetscScalar sum[128];
427: for (j=0; j<bs; j++) sum[j] = 0.0;
428: for (i=0; i<n; i+=bs) {
429: for (j=0; j<bs; j++) sum[j] += x[i+j]*(PetscConj(x[i+j]));
430: }
431: for (j=0; j<bs; j++) tnorm[j] = PetscRealPart(sum[j]);
433: MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_SUM,comm);
434: for (j=0; j<bs; j++) nrm[j] = PetscSqrtReal(nrm[j]);
435: } else if (ntype == NORM_1) {
436: for (j=0; j<bs; j++) tnorm[j] = 0.0;
438: for (i=0; i<n; i+=bs) {
439: for (j=0; j<bs; j++) tnorm[j] += PetscAbsScalar(x[i+j]);
440: }
442: MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_SUM,comm);
443: } else if (ntype == NORM_INFINITY) {
444: PetscReal tmp;
445: for (j=0; j<bs; j++) tnorm[j] = 0.0;
447: for (i=0; i<n; i+=bs) {
448: for (j=0; j<bs; j++) {
449: if ((tmp = PetscAbsScalar(x[i+j])) > tnorm[j]) tnorm[j] = tmp;
450: /* check special case of tmp == NaN */
451: if (tmp != tmp) {tnorm[j] = tmp; break;}
452: }
453: }
454: MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_MAX,comm);
455: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
456: VecRestoreArrayRead(v,&x);
457: return(0);
458: }
460: /*@
461: VecStrideMaxAll - Computes the maximums of subvectors of a vector defined
462: by a starting point and a stride and optionally its location.
464: Collective on Vec
466: Input Parameter:
467: . v - the vector
469: Output Parameter:
470: + index - the location where the maximum occurred (not supported, pass NULL,
471: if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
472: - nrm - the maximum values of each subvector
474: Notes:
475: One must call VecSetBlockSize() before this routine to set the stride
476: information, or use a vector created from a multicomponent DMDA.
478: The dimension of nrm must be the same as the vector block size
480: Level: advanced
483: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
484: @*/
485: PetscErrorCode VecStrideMaxAll(Vec v,PetscInt idex[],PetscReal nrm[])
486: {
487: PetscErrorCode ierr;
488: PetscInt i,j,n,bs;
489: const PetscScalar *x;
490: PetscReal max[128],tmp;
491: MPI_Comm comm;
496: 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");
497: VecGetLocalSize(v,&n);
498: VecGetArrayRead(v,&x);
499: PetscObjectGetComm((PetscObject)v,&comm);
501: VecGetBlockSize(v,&bs);
502: if (bs > 128) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");
504: if (!n) {
505: for (j=0; j<bs; j++) max[j] = PETSC_MIN_REAL;
506: } else {
507: for (j=0; j<bs; j++) max[j] = PetscRealPart(x[j]);
509: for (i=bs; i<n; i+=bs) {
510: for (j=0; j<bs; j++) {
511: if ((tmp = PetscRealPart(x[i+j])) > max[j]) max[j] = tmp;
512: }
513: }
514: }
515: MPIU_Allreduce(max,nrm,bs,MPIU_REAL,MPIU_MAX,comm);
517: VecRestoreArrayRead(v,&x);
518: return(0);
519: }
521: /*@
522: VecStrideMinAll - Computes the minimum of subvector of a vector defined
523: by a starting point and a stride and optionally its location.
525: Collective on Vec
527: Input Parameter:
528: . v - the vector
530: Output Parameter:
531: + idex - the location where the minimum occurred (not supported, pass NULL,
532: if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
533: - nrm - the minimums of each subvector
535: Level: advanced
537: Notes:
538: One must call VecSetBlockSize() before this routine to set the stride
539: information, or use a vector created from a multicomponent DMDA.
541: The dimension of nrm must be the same as the vector block size
544: .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
545: @*/
546: PetscErrorCode VecStrideMinAll(Vec v,PetscInt idex[],PetscReal nrm[])
547: {
548: PetscErrorCode ierr;
549: PetscInt i,n,bs,j;
550: const PetscScalar *x;
551: PetscReal min[128],tmp;
552: MPI_Comm comm;
557: 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");
558: VecGetLocalSize(v,&n);
559: VecGetArrayRead(v,&x);
560: PetscObjectGetComm((PetscObject)v,&comm);
562: VecGetBlockSize(v,&bs);
563: if (bs > 128) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");
565: if (!n) {
566: for (j=0; j<bs; j++) min[j] = PETSC_MAX_REAL;
567: } else {
568: for (j=0; j<bs; j++) min[j] = PetscRealPart(x[j]);
570: for (i=bs; i<n; i+=bs) {
571: for (j=0; j<bs; j++) {
572: if ((tmp = PetscRealPart(x[i+j])) < min[j]) min[j] = tmp;
573: }
574: }
575: }
576: MPIU_Allreduce(min,nrm,bs,MPIU_REAL,MPIU_MIN,comm);
578: VecRestoreArrayRead(v,&x);
579: return(0);
580: }
582: /*----------------------------------------------------------------------------------------------*/
583: /*@
584: VecStrideGatherAll - Gathers all the single components from a multi-component vector into
585: separate vectors.
587: Collective on Vec
589: Input Parameter:
590: + v - the vector
591: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
593: Output Parameter:
594: . s - the location where the subvectors are stored
596: Notes:
597: One must call VecSetBlockSize() before this routine to set the stride
598: information, or use a vector created from a multicomponent DMDA.
600: If x is the array representing the vector x then this gathers
601: the arrays (x[start],x[start+stride],x[start+2*stride], ....)
602: for start=0,1,2,...bs-1
604: The parallel layout of the vector and the subvector must be the same;
605: i.e., nlocal of v = stride*(nlocal of s)
607: Not optimized; could be easily
609: Level: advanced
611: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
612: VecStrideScatterAll()
613: @*/
614: PetscErrorCode VecStrideGatherAll(Vec v,Vec s[],InsertMode addv)
615: {
616: PetscErrorCode ierr;
617: PetscInt i,n,n2,bs,j,k,*bss = NULL,nv,jj,nvc;
618: PetscScalar **y;
619: const PetscScalar *x;
625: VecGetLocalSize(v,&n);
626: VecGetLocalSize(s[0],&n2);
627: VecGetArrayRead(v,&x);
628: VecGetBlockSize(v,&bs);
629: if (bs <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Input vector does not have a valid blocksize set");
631: PetscMalloc2(bs,&y,bs,&bss);
632: nv = 0;
633: nvc = 0;
634: for (i=0; i<bs; i++) {
635: VecGetBlockSize(s[i],&bss[i]);
636: if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1 Re: [PETSC #8241] VecStrideGatherAll */
637: VecGetArray(s[i],&y[i]);
638: nvc += bss[i];
639: nv++;
640: if (nvc > bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of subvectors in subvectors > number of vectors in main vector");
641: if (nvc == bs) break;
642: }
644: n = n/bs;
646: jj = 0;
647: if (addv == INSERT_VALUES) {
648: for (j=0; j<nv; j++) {
649: for (k=0; k<bss[j]; k++) {
650: for (i=0; i<n; i++) y[j][i*bss[j] + k] = x[bs*i+jj+k];
651: }
652: jj += bss[j];
653: }
654: } else if (addv == ADD_VALUES) {
655: for (j=0; j<nv; j++) {
656: for (k=0; k<bss[j]; k++) {
657: for (i=0; i<n; i++) y[j][i*bss[j] + k] += x[bs*i+jj+k];
658: }
659: jj += bss[j];
660: }
661: #if !defined(PETSC_USE_COMPLEX)
662: } else if (addv == MAX_VALUES) {
663: for (j=0; j<nv; j++) {
664: for (k=0; k<bss[j]; k++) {
665: for (i=0; i<n; i++) y[j][i*bss[j] + k] = PetscMax(y[j][i*bss[j] + k],x[bs*i+jj+k]);
666: }
667: jj += bss[j];
668: }
669: #endif
670: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
672: VecRestoreArrayRead(v,&x);
673: for (i=0; i<nv; i++) {
674: VecRestoreArray(s[i],&y[i]);
675: }
677: PetscFree2(y,bss);
678: return(0);
679: }
681: /*@
682: VecStrideScatterAll - Scatters all the single components from separate vectors into
683: a multi-component vector.
685: Collective on Vec
687: Input Parameter:
688: + s - the location where the subvectors are stored
689: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
691: Output Parameter:
692: . v - the multicomponent vector
694: Notes:
695: One must call VecSetBlockSize() before this routine to set the stride
696: information, or use a vector created from a multicomponent DMDA.
698: The parallel layout of the vector and the subvector must be the same;
699: i.e., nlocal of v = stride*(nlocal of s)
701: Not optimized; could be easily
703: Level: advanced
705: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
706: VecStrideScatterAll()
707: @*/
708: PetscErrorCode VecStrideScatterAll(Vec s[],Vec v,InsertMode addv)
709: {
710: PetscErrorCode ierr;
711: PetscInt i,n,n2,bs,j,jj,k,*bss = NULL,nv,nvc;
712: PetscScalar *x;
713: PetscScalar const **y;
719: VecGetLocalSize(v,&n);
720: VecGetLocalSize(s[0],&n2);
721: VecGetArray(v,&x);
722: VecGetBlockSize(v,&bs);
723: if (bs <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Input vector does not have a valid blocksize set");
725: PetscMalloc2(bs,(PetscScalar***)&y,bs,&bss);
726: nv = 0;
727: nvc = 0;
728: for (i=0; i<bs; i++) {
729: VecGetBlockSize(s[i],&bss[i]);
730: if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1 Re: [PETSC #8241] VecStrideGatherAll */
731: VecGetArrayRead(s[i],&y[i]);
732: nvc += bss[i];
733: nv++;
734: if (nvc > bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of subvectors in subvectors > number of vectors in main vector");
735: if (nvc == bs) break;
736: }
738: n = n/bs;
740: jj = 0;
741: if (addv == INSERT_VALUES) {
742: for (j=0; j<nv; j++) {
743: for (k=0; k<bss[j]; k++) {
744: for (i=0; i<n; i++) x[bs*i+jj+k] = y[j][i*bss[j] + k];
745: }
746: jj += bss[j];
747: }
748: } else if (addv == ADD_VALUES) {
749: for (j=0; j<nv; j++) {
750: for (k=0; k<bss[j]; k++) {
751: for (i=0; i<n; i++) x[bs*i+jj+k] += y[j][i*bss[j] + k];
752: }
753: jj += bss[j];
754: }
755: #if !defined(PETSC_USE_COMPLEX)
756: } else if (addv == MAX_VALUES) {
757: for (j=0; j<nv; j++) {
758: for (k=0; k<bss[j]; k++) {
759: for (i=0; i<n; i++) x[bs*i+jj+k] = PetscMax(x[bs*i+jj+k],y[j][i*bss[j] + k]);
760: }
761: jj += bss[j];
762: }
763: #endif
764: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
766: VecRestoreArray(v,&x);
767: for (i=0; i<nv; i++) {
768: VecRestoreArrayRead(s[i],&y[i]);
769: }
770: PetscFree2(*(PetscScalar***)&y,bss);
771: return(0);
772: }
774: /*@
775: VecStrideGather - Gathers a single component from a multi-component vector into
776: another vector.
778: Collective on Vec
780: Input Parameter:
781: + v - the vector
782: . start - starting point of the subvector (defined by a stride)
783: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
785: Output Parameter:
786: . s - the location where the subvector is stored
788: Notes:
789: One must call VecSetBlockSize() before this routine to set the stride
790: information, or use a vector created from a multicomponent DMDA.
792: If x is the array representing the vector x then this gathers
793: the array (x[start],x[start+stride],x[start+2*stride], ....)
795: The parallel layout of the vector and the subvector must be the same;
796: i.e., nlocal of v = stride*(nlocal of s)
798: Not optimized; could be easily
800: Level: advanced
802: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
803: VecStrideScatterAll()
804: @*/
805: PetscErrorCode VecStrideGather(Vec v,PetscInt start,Vec s,InsertMode addv)
806: {
812: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
813: 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);
814: if (!v->ops->stridegather) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
815: (*v->ops->stridegather)(v,start,s,addv);
816: return(0);
817: }
819: /*@
820: VecStrideScatter - Scatters a single component from a vector into a multi-component vector.
822: Collective on Vec
824: Input Parameter:
825: + s - the single-component vector
826: . start - starting point of the subvector (defined by a stride)
827: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
829: Output Parameter:
830: . v - the location where the subvector is scattered (the multi-component vector)
832: Notes:
833: One must call VecSetBlockSize() on the multi-component vector before this
834: routine to set the stride information, or use a vector created from a multicomponent DMDA.
836: The parallel layout of the vector and the subvector must be the same;
837: i.e., nlocal of v = stride*(nlocal of s)
839: Not optimized; could be easily
841: Level: advanced
843: .seealso: VecStrideNorm(), VecStrideGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
844: VecStrideScatterAll(), VecStrideSubSetScatter(), VecStrideSubSetGather()
845: @*/
846: PetscErrorCode VecStrideScatter(Vec s,PetscInt start,Vec v,InsertMode addv)
847: {
853: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
854: 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);
855: if (!v->ops->stridescatter) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
856: (*v->ops->stridescatter)(s,start,v,addv);
857: return(0);
858: }
860: /*@
861: VecStrideSubSetGather - Gathers a subset of components from a multi-component vector into
862: another vector.
864: Collective on Vec
866: Input Parameter:
867: + v - the vector
868: . nidx - the number of indices
869: . idxv - the indices of the components 0 <= idxv[0] ...idxv[nidx-1] < bs(v), they need not be sorted
870: . 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
871: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
873: Output Parameter:
874: . s - the location where the subvector is stored
876: Notes:
877: One must call VecSetBlockSize() on both vectors before this routine to set the stride
878: information, or use a vector created from a multicomponent DMDA.
881: The parallel layout of the vector and the subvector must be the same;
883: Not optimized; could be easily
885: Level: advanced
887: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideGather(), VecStrideSubSetScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
888: VecStrideScatterAll()
889: @*/
890: PetscErrorCode VecStrideSubSetGather(Vec v,PetscInt nidx,const PetscInt idxv[],const PetscInt idxs[],Vec s,InsertMode addv)
891: {
897: if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
898: if (!v->ops->stridesubsetgather) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
899: (*v->ops->stridesubsetgather)(v,nidx,idxv,idxs,s,addv);
900: return(0);
901: }
903: /*@
904: VecStrideSubSetScatter - Scatters components from a vector into a subset of components of a multi-component vector.
906: Collective on Vec
908: Input Parameter:
909: + s - the smaller-component vector
910: . nidx - the number of indices in idx
911: . 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
912: . idxv - the indices of the components in the larger-component vector, 0 <= idx[0] ...idx[nidx-1] < bs(v) they need not be sorted
913: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
915: Output Parameter:
916: . v - the location where the subvector is into scattered (the multi-component vector)
918: Notes:
919: One must call VecSetBlockSize() on the vectors before this
920: routine to set the stride information, or use a vector created from a multicomponent DMDA.
922: The parallel layout of the vector and the subvector must be the same;
924: Not optimized; could be easily
926: Level: advanced
928: .seealso: VecStrideNorm(), VecStrideGather(), VecStrideGather(), VecStrideSubSetGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
929: VecStrideScatterAll()
930: @*/
931: PetscErrorCode VecStrideSubSetScatter(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)
932: {
938: if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
939: if (!v->ops->stridesubsetscatter) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
940: (*v->ops->stridesubsetscatter)(s,nidx,idxs,idxv,v,addv);
941: return(0);
942: }
944: PetscErrorCode VecStrideGather_Default(Vec v,PetscInt start,Vec s,InsertMode addv)
945: {
947: PetscInt i,n,bs,ns;
948: const PetscScalar *x;
949: PetscScalar *y;
952: VecGetLocalSize(v,&n);
953: VecGetLocalSize(s,&ns);
954: VecGetArrayRead(v,&x);
955: VecGetArray(s,&y);
957: bs = v->map->bs;
958: 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);
959: x += start;
960: n = n/bs;
962: if (addv == INSERT_VALUES) {
963: for (i=0; i<n; i++) y[i] = x[bs*i];
964: } else if (addv == ADD_VALUES) {
965: for (i=0; i<n; i++) y[i] += x[bs*i];
966: #if !defined(PETSC_USE_COMPLEX)
967: } else if (addv == MAX_VALUES) {
968: for (i=0; i<n; i++) y[i] = PetscMax(y[i],x[bs*i]);
969: #endif
970: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
972: VecRestoreArrayRead(v,&x);
973: VecRestoreArray(s,&y);
974: return(0);
975: }
977: PetscErrorCode VecStrideScatter_Default(Vec s,PetscInt start,Vec v,InsertMode addv)
978: {
979: PetscErrorCode ierr;
980: PetscInt i,n,bs,ns;
981: PetscScalar *x;
982: const PetscScalar *y;
985: VecGetLocalSize(v,&n);
986: VecGetLocalSize(s,&ns);
987: VecGetArray(v,&x);
988: VecGetArrayRead(s,&y);
990: bs = v->map->bs;
991: 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);
992: x += start;
993: n = n/bs;
995: if (addv == INSERT_VALUES) {
996: for (i=0; i<n; i++) x[bs*i] = y[i];
997: } else if (addv == ADD_VALUES) {
998: for (i=0; i<n; i++) x[bs*i] += y[i];
999: #if !defined(PETSC_USE_COMPLEX)
1000: } else if (addv == MAX_VALUES) {
1001: for (i=0; i<n; i++) x[bs*i] = PetscMax(y[i],x[bs*i]);
1002: #endif
1003: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
1005: VecRestoreArray(v,&x);
1006: VecRestoreArrayRead(s,&y);
1007: return(0);
1008: }
1010: PetscErrorCode VecStrideSubSetGather_Default(Vec v,PetscInt nidx,const PetscInt idxv[],const PetscInt idxs[],Vec s,InsertMode addv)
1011: {
1012: PetscErrorCode ierr;
1013: PetscInt i,j,n,bs,bss,ns;
1014: const PetscScalar *x;
1015: PetscScalar *y;
1018: VecGetLocalSize(v,&n);
1019: VecGetLocalSize(s,&ns);
1020: VecGetArrayRead(v,&x);
1021: VecGetArray(s,&y);
1023: bs = v->map->bs;
1024: bss = s->map->bs;
1025: n = n/bs;
1027: if (PetscDefined(USE_DEBUG)) {
1028: if (n != ns/bss) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible layout of vectors");
1029: for (j=0; j<nidx; j++) {
1030: if (idxv[j] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is negative",j,idxv[j]);
1031: 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);
1032: }
1033: if (!idxs && bss != nidx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide idxs when not gathering into all locations");
1034: }
1036: if (addv == INSERT_VALUES) {
1037: if (!idxs) {
1038: for (i=0; i<n; i++) {
1039: for (j=0; j<bss; j++) y[bss*i+j] = x[bs*i+idxv[j]];
1040: }
1041: } else {
1042: for (i=0; i<n; i++) {
1043: for (j=0; j<bss; j++) y[bss*i+idxs[j]] = x[bs*i+idxv[j]];
1044: }
1045: }
1046: } else if (addv == ADD_VALUES) {
1047: if (!idxs) {
1048: for (i=0; i<n; i++) {
1049: for (j=0; j<bss; j++) y[bss*i+j] += x[bs*i+idxv[j]];
1050: }
1051: } else {
1052: for (i=0; i<n; i++) {
1053: for (j=0; j<bss; j++) y[bss*i+idxs[j]] += x[bs*i+idxv[j]];
1054: }
1055: }
1056: #if !defined(PETSC_USE_COMPLEX)
1057: } else if (addv == MAX_VALUES) {
1058: if (!idxs) {
1059: for (i=0; i<n; i++) {
1060: for (j=0; j<bss; j++) y[bss*i+j] = PetscMax(y[bss*i+j],x[bs*i+idxv[j]]);
1061: }
1062: } else {
1063: for (i=0; i<n; i++) {
1064: for (j=0; j<bss; j++) y[bss*i+idxs[j]] = PetscMax(y[bss*i+idxs[j]],x[bs*i+idxv[j]]);
1065: }
1066: }
1067: #endif
1068: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
1070: VecRestoreArrayRead(v,&x);
1071: VecRestoreArray(s,&y);
1072: return(0);
1073: }
1075: PetscErrorCode VecStrideSubSetScatter_Default(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)
1076: {
1077: PetscErrorCode ierr;
1078: PetscInt j,i,n,bs,ns,bss;
1079: PetscScalar *x;
1080: const PetscScalar *y;
1083: VecGetLocalSize(v,&n);
1084: VecGetLocalSize(s,&ns);
1085: VecGetArray(v,&x);
1086: VecGetArrayRead(s,&y);
1088: bs = v->map->bs;
1089: bss = s->map->bs;
1090: n = n/bs;
1092: if (PetscDefined(USE_DEBUG)) {
1093: if (n != ns/bss) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible layout of vectors");
1094: for (j=0; j<bss; j++) {
1095: if (idxs) {
1096: if (idxs[j] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is negative",j,idxs[j]);
1097: 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);
1098: }
1099: }
1100: if (!idxs && bss != nidx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide idxs when not scattering from all locations");
1101: }
1103: if (addv == INSERT_VALUES) {
1104: if (!idxs) {
1105: for (i=0; i<n; i++) {
1106: for (j=0; j<bss; j++) x[bs*i + idxv[j]] = y[bss*i+j];
1107: }
1108: } else {
1109: for (i=0; i<n; i++) {
1110: for (j=0; j<bss; j++) x[bs*i + idxv[j]] = y[bss*i+idxs[j]];
1111: }
1112: }
1113: } else if (addv == ADD_VALUES) {
1114: if (!idxs) {
1115: for (i=0; i<n; i++) {
1116: for (j=0; j<bss; j++) x[bs*i + idxv[j]] += y[bss*i+j];
1117: }
1118: } else {
1119: for (i=0; i<n; i++) {
1120: for (j=0; j<bss; j++) x[bs*i + idxv[j]] += y[bss*i+idxs[j]];
1121: }
1122: }
1123: #if !defined(PETSC_USE_COMPLEX)
1124: } else if (addv == MAX_VALUES) {
1125: if (!idxs) {
1126: for (i=0; i<n; i++) {
1127: for (j=0; j<bss; j++) x[bs*i + idxv[j]] = PetscMax(y[bss*i+j],x[bs*i + idxv[j]]);
1128: }
1129: } else {
1130: for (i=0; i<n; i++) {
1131: for (j=0; j<bss; j++) x[bs*i + idxv[j]] = PetscMax(y[bss*i+idxs[j]],x[bs*i + idxv[j]]);
1132: }
1133: }
1134: #endif
1135: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
1137: VecRestoreArray(v,&x);
1138: VecRestoreArrayRead(s,&y);
1139: return(0);
1140: }
1142: PetscErrorCode VecReciprocal_Default(Vec v)
1143: {
1145: PetscInt i,n;
1146: PetscScalar *x;
1149: VecGetLocalSize(v,&n);
1150: VecGetArray(v,&x);
1151: for (i=0; i<n; i++) {
1152: if (x[i] != (PetscScalar)0.0) x[i] = (PetscScalar)1.0/x[i];
1153: }
1154: VecRestoreArray(v,&x);
1155: return(0);
1156: }
1158: /*@
1159: VecExp - Replaces each component of a vector by e^x_i
1161: Not collective
1163: Input Parameter:
1164: . v - The vector
1166: Output Parameter:
1167: . v - The vector of exponents
1169: Level: beginner
1171: .seealso: VecLog(), VecAbs(), VecSqrtAbs(), VecReciprocal()
1173: @*/
1174: PetscErrorCode VecExp(Vec v)
1175: {
1176: PetscScalar *x;
1177: PetscInt i, n;
1182: if (v->ops->exp) {
1183: (*v->ops->exp)(v);
1184: } else {
1185: VecGetLocalSize(v, &n);
1186: VecGetArray(v, &x);
1187: for (i = 0; i < n; i++) x[i] = PetscExpScalar(x[i]);
1188: VecRestoreArray(v, &x);
1189: }
1190: return(0);
1191: }
1193: /*@
1194: VecLog - Replaces each component of a vector by log(x_i), the natural logarithm
1196: Not collective
1198: Input Parameter:
1199: . v - The vector
1201: Output Parameter:
1202: . v - The vector of logs
1204: Level: beginner
1206: .seealso: VecExp(), VecAbs(), VecSqrtAbs(), VecReciprocal()
1208: @*/
1209: PetscErrorCode VecLog(Vec v)
1210: {
1211: PetscScalar *x;
1212: PetscInt i, n;
1217: if (v->ops->log) {
1218: (*v->ops->log)(v);
1219: } else {
1220: VecGetLocalSize(v, &n);
1221: VecGetArray(v, &x);
1222: for (i = 0; i < n; i++) x[i] = PetscLogScalar(x[i]);
1223: VecRestoreArray(v, &x);
1224: }
1225: return(0);
1226: }
1228: /*@
1229: VecSqrtAbs - Replaces each component of a vector by the square root of its magnitude.
1231: Not collective
1233: Input Parameter:
1234: . v - The vector
1236: Output Parameter:
1237: . v - The vector square root
1239: Level: beginner
1241: Note: The actual function is sqrt(|x_i|)
1243: .seealso: VecLog(), VecExp(), VecReciprocal(), VecAbs()
1245: @*/
1246: PetscErrorCode VecSqrtAbs(Vec v)
1247: {
1248: PetscScalar *x;
1249: PetscInt i, n;
1254: if (v->ops->sqrt) {
1255: (*v->ops->sqrt)(v);
1256: } else {
1257: VecGetLocalSize(v, &n);
1258: VecGetArray(v, &x);
1259: for (i = 0; i < n; i++) x[i] = PetscSqrtReal(PetscAbsScalar(x[i]));
1260: VecRestoreArray(v, &x);
1261: }
1262: return(0);
1263: }
1265: /*@
1266: VecDotNorm2 - computes the inner product of two vectors and the 2-norm squared of the second vector
1268: Collective on Vec
1270: Input Parameter:
1271: + s - first vector
1272: - t - second vector
1274: Output Parameter:
1275: + dp - s'conj(t)
1276: - nm - t'conj(t)
1278: Level: advanced
1280: Notes:
1281: conj(x) is the complex conjugate of x when x is complex
1284: .seealso: VecDot(), VecNorm(), VecDotBegin(), VecNormBegin(), VecDotEnd(), VecNormEnd()
1286: @*/
1287: PetscErrorCode VecDotNorm2(Vec s,Vec t,PetscScalar *dp, PetscReal *nm)
1288: {
1289: const PetscScalar *sx, *tx;
1290: PetscScalar dpx = 0.0, nmx = 0.0,work[2],sum[2];
1291: PetscInt i, n;
1292: PetscErrorCode ierr;
1302: if (s->map->N != t->map->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
1303: if (s->map->n != t->map->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");
1305: PetscLogEventBegin(VEC_DotNorm2,s,t,0,0);
1306: if (s->ops->dotnorm2) {
1307: (*s->ops->dotnorm2)(s,t,dp,&dpx);
1308: *nm = PetscRealPart(dpx);
1309: } else {
1310: VecGetLocalSize(s, &n);
1311: VecGetArrayRead(s, &sx);
1312: VecGetArrayRead(t, &tx);
1314: for (i = 0; i<n; i++) {
1315: dpx += sx[i]*PetscConj(tx[i]);
1316: nmx += tx[i]*PetscConj(tx[i]);
1317: }
1318: work[0] = dpx;
1319: work[1] = nmx;
1321: MPIU_Allreduce(work,sum,2,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)s));
1322: *dp = sum[0];
1323: *nm = PetscRealPart(sum[1]);
1325: VecRestoreArrayRead(t, &tx);
1326: VecRestoreArrayRead(s, &sx);
1327: PetscLogFlops(4.0*n);
1328: }
1329: PetscLogEventEnd(VEC_DotNorm2,s,t,0,0);
1330: return(0);
1331: }
1333: /*@
1334: VecSum - Computes the sum of all the components of a vector.
1336: Collective on Vec
1338: Input Parameter:
1339: . v - the vector
1341: Output Parameter:
1342: . sum - the result
1344: Level: beginner
1346: .seealso: VecNorm()
1347: @*/
1348: PetscErrorCode VecSum(Vec v,PetscScalar *sum)
1349: {
1350: PetscErrorCode ierr;
1351: PetscInt i,n;
1352: const PetscScalar *x;
1353: PetscScalar lsum = 0.0;
1358: VecGetLocalSize(v,&n);
1359: VecGetArrayRead(v,&x);
1360: for (i=0; i<n; i++) lsum += x[i];
1361: MPIU_Allreduce(&lsum,sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)v));
1362: VecRestoreArrayRead(v,&x);
1363: return(0);
1364: }
1366: /*@
1367: VecImaginaryPart - Replaces a complex vector with its imginary part
1369: Collective on Vec
1371: Input Parameter:
1372: . v - the vector
1374: Level: beginner
1376: .seealso: VecNorm(), VecRealPart()
1377: @*/
1378: PetscErrorCode VecImaginaryPart(Vec v)
1379: {
1380: PetscErrorCode ierr;
1381: PetscInt i,n;
1382: PetscScalar *x;
1386: VecGetLocalSize(v,&n);
1387: VecGetArray(v,&x);
1388: for (i=0; i<n; i++) x[i] = PetscImaginaryPart(x[i]);
1389: VecRestoreArray(v,&x);
1390: return(0);
1391: }
1393: /*@
1394: VecRealPart - Replaces a complex vector with its real part
1396: Collective on Vec
1398: Input Parameter:
1399: . v - the vector
1401: Level: beginner
1403: .seealso: VecNorm(), VecImaginaryPart()
1404: @*/
1405: PetscErrorCode VecRealPart(Vec v)
1406: {
1407: PetscErrorCode ierr;
1408: PetscInt i,n;
1409: PetscScalar *x;
1413: VecGetLocalSize(v,&n);
1414: VecGetArray(v,&x);
1415: for (i=0; i<n; i++) x[i] = PetscRealPart(x[i]);
1416: VecRestoreArray(v,&x);
1417: return(0);
1418: }
1420: /*@
1421: VecShift - Shifts all of the components of a vector by computing
1422: x[i] = x[i] + shift.
1424: Logically Collective on Vec
1426: Input Parameters:
1427: + v - the vector
1428: - shift - the shift
1430: Output Parameter:
1431: . v - the shifted vector
1433: Level: intermediate
1435: @*/
1436: PetscErrorCode VecShift(Vec v,PetscScalar shift)
1437: {
1439: PetscInt i,n;
1440: PetscScalar *x;
1445: VecSetErrorIfLocked(v,1);
1446: if (shift == 0.0) return(0);
1448: if (v->ops->shift) {
1449: (*v->ops->shift)(v,shift);
1450: } else {
1451: VecGetLocalSize(v,&n);
1452: VecGetArray(v,&x);
1453: for (i=0; i<n; i++) x[i] += shift;
1454: VecRestoreArray(v,&x);
1455: }
1456: return(0);
1457: }
1459: /*@
1460: VecAbs - Replaces every element in a vector with its absolute value.
1462: Logically Collective on Vec
1464: Input Parameters:
1465: . v - the vector
1467: Level: intermediate
1469: @*/
1470: PetscErrorCode VecAbs(Vec v)
1471: {
1473: PetscInt i,n;
1474: PetscScalar *x;
1478: VecSetErrorIfLocked(v,1);
1480: if (v->ops->abs) {
1481: (*v->ops->abs)(v);
1482: } else {
1483: VecGetLocalSize(v,&n);
1484: VecGetArray(v,&x);
1485: for (i=0; i<n; i++) x[i] = PetscAbsScalar(x[i]);
1486: VecRestoreArray(v,&x);
1487: }
1488: return(0);
1489: }
1491: /*@
1492: VecPermute - Permutes a vector in place using the given ordering.
1494: Input Parameters:
1495: + vec - The vector
1496: . order - The ordering
1497: - inv - The flag for inverting the permutation
1499: Level: beginner
1501: Note: This function does not yet support parallel Index Sets with non-local permutations
1503: .seealso: MatPermute()
1504: @*/
1505: PetscErrorCode VecPermute(Vec x, IS row, PetscBool inv)
1506: {
1507: const PetscScalar *array;
1508: PetscScalar *newArray;
1509: const PetscInt *idx;
1510: PetscInt i,rstart,rend;
1511: PetscErrorCode ierr;
1514: VecSetErrorIfLocked(x,1);
1516: VecGetOwnershipRange(x,&rstart,&rend);
1517: ISGetIndices(row, &idx);
1518: VecGetArrayRead(x, &array);
1519: PetscMalloc1(x->map->n, &newArray);
1520: if (PetscDefined(USE_DEBUG)) {
1521: for (i = 0; i < x->map->n; i++) {
1522: 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]);
1523: }
1524: }
1525: if (!inv) {
1526: for (i = 0; i < x->map->n; i++) newArray[i] = array[idx[i]-rstart];
1527: } else {
1528: for (i = 0; i < x->map->n; i++) newArray[idx[i]-rstart] = array[i];
1529: }
1530: VecRestoreArrayRead(x, &array);
1531: ISRestoreIndices(row, &idx);
1532: VecReplaceArray(x, newArray);
1533: return(0);
1534: }
1536: /*@
1537: VecEqual - Compares two vectors. Returns true if the two vectors are either pointing to the same memory buffer,
1538: or if the two vectors have the same local and global layout as well as bitwise equality of all entries.
1539: Does NOT take round-off errors into account.
1541: Collective on Vec
1543: Input Parameters:
1544: + vec1 - the first vector
1545: - vec2 - the second vector
1547: Output Parameter:
1548: . flg - PETSC_TRUE if the vectors are equal; PETSC_FALSE otherwise.
1550: Level: intermediate
1553: @*/
1554: PetscErrorCode VecEqual(Vec vec1,Vec vec2,PetscBool *flg)
1555: {
1556: const PetscScalar *v1,*v2;
1557: PetscErrorCode ierr;
1558: PetscInt n1,n2,N1,N2;
1559: PetscBool flg1;
1565: if (vec1 == vec2) *flg = PETSC_TRUE;
1566: else {
1567: VecGetSize(vec1,&N1);
1568: VecGetSize(vec2,&N2);
1569: if (N1 != N2) flg1 = PETSC_FALSE;
1570: else {
1571: VecGetLocalSize(vec1,&n1);
1572: VecGetLocalSize(vec2,&n2);
1573: if (n1 != n2) flg1 = PETSC_FALSE;
1574: else {
1575: VecGetArrayRead(vec1,&v1);
1576: VecGetArrayRead(vec2,&v2);
1577: PetscArraycmp(v1,v2,n1,&flg1);
1578: VecRestoreArrayRead(vec1,&v1);
1579: VecRestoreArrayRead(vec2,&v2);
1580: }
1581: }
1582: /* combine results from all processors */
1583: MPIU_Allreduce(&flg1,flg,1,MPIU_BOOL,MPI_MIN,PetscObjectComm((PetscObject)vec1));
1584: }
1585: return(0);
1586: }
1588: /*@
1589: VecUniqueEntries - Compute the number of unique entries, and those entries
1591: Collective on Vec
1593: Input Parameter:
1594: . vec - the vector
1596: Output Parameters:
1597: + n - The number of unique entries
1598: - e - The entries
1600: Level: intermediate
1602: @*/
1603: PetscErrorCode VecUniqueEntries(Vec vec, PetscInt *n, PetscScalar **e)
1604: {
1605: const PetscScalar *v;
1606: PetscScalar *tmp, *vals;
1607: PetscMPIInt *N, *displs, l;
1608: PetscInt ng, m, i, j, p;
1609: PetscMPIInt size;
1610: PetscErrorCode ierr;
1615: MPI_Comm_size(PetscObjectComm((PetscObject) vec), &size);CHKERRMPI(ierr);
1616: VecGetLocalSize(vec, &m);
1617: VecGetArrayRead(vec, &v);
1618: PetscMalloc2(m,&tmp,size,&N);
1619: for (i = 0, j = 0, l = 0; i < m; ++i) {
1620: /* Can speed this up with sorting */
1621: for (j = 0; j < l; ++j) {
1622: if (v[i] == tmp[j]) break;
1623: }
1624: if (j == l) {
1625: tmp[j] = v[i];
1626: ++l;
1627: }
1628: }
1629: VecRestoreArrayRead(vec, &v);
1630: /* Gather serial results */
1631: MPI_Allgather(&l, 1, MPI_INT, N, 1, MPI_INT, PetscObjectComm((PetscObject) vec));CHKERRMPI(ierr);
1632: for (p = 0, ng = 0; p < size; ++p) {
1633: ng += N[p];
1634: }
1635: PetscMalloc2(ng,&vals,size+1,&displs);
1636: for (p = 1, displs[0] = 0; p <= size; ++p) {
1637: displs[p] = displs[p-1] + N[p-1];
1638: }
1639: MPI_Allgatherv(tmp, l, MPIU_SCALAR, vals, N, displs, MPIU_SCALAR, PetscObjectComm((PetscObject) vec));CHKERRMPI(ierr);
1640: /* Find unique entries */
1641: #ifdef PETSC_USE_COMPLEX
1642: SETERRQ(PetscObjectComm((PetscObject) vec), PETSC_ERR_SUP, "Does not work with complex numbers");
1643: #else
1644: *n = displs[size];
1645: PetscSortRemoveDupsReal(n, (PetscReal *) vals);
1646: if (e) {
1648: PetscMalloc1(*n, e);
1649: for (i = 0; i < *n; ++i) {
1650: (*e)[i] = vals[i];
1651: }
1652: }
1653: PetscFree2(vals,displs);
1654: PetscFree2(tmp,N);
1655: return(0);
1656: #endif
1657: }