Actual source code: vinv.c
petsc-main 2021-04-15
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 Parameter:
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
27: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideScale()
28: @*/
29: PetscErrorCode VecStrideSet(Vec v,PetscInt start,PetscScalar s)
30: {
32: PetscInt i,n,bs;
33: PetscScalar *x;
39: VecSetErrorIfLocked(v,1);
41: VecGetLocalSize(v,&n);
42: VecGetArray(v,&x);
43: VecGetBlockSize(v,&bs);
44: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
45: 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);
46: x += start;
47: for (i=0; i<n; i+=bs) x[i] = s;
48: x -= start;
49: VecRestoreArray(v,&x);
50: return(0);
51: }
53: /*@
54: VecStrideScale - Scales a subvector of a vector defined
55: by a starting point and a stride.
57: Logically Collective on Vec
59: Input Parameter:
60: + v - the vector
61: . start - starting point of the subvector (defined by a stride)
62: - scale - value to multiply each subvector entry by
64: Notes:
65: One must call VecSetBlockSize() before this routine to set the stride
66: information, or use a vector created from a multicomponent DMDA.
68: This will only work if the desire subvector is a stride subvector
70: Level: advanced
73: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideScale()
74: @*/
75: PetscErrorCode VecStrideScale(Vec v,PetscInt start,PetscScalar scale)
76: {
78: PetscInt i,n,bs;
79: PetscScalar *x;
85: VecSetErrorIfLocked(v,1);
87: VecGetLocalSize(v,&n);
88: VecGetArray(v,&x);
89: VecGetBlockSize(v,&bs);
90: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
91: 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);
92: x += start;
93: for (i=0; i<n; i+=bs) x[i] *= scale;
94: x -= start;
95: VecRestoreArray(v,&x);
96: return(0);
97: }
99: /*@
100: VecStrideNorm - Computes the norm of subvector of a vector defined
101: by a starting point and a stride.
103: Collective on Vec
105: Input Parameter:
106: + v - the vector
107: . start - starting point of the subvector (defined by a stride)
108: - ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY
110: Output Parameter:
111: . norm - the norm
113: Notes:
114: One must call VecSetBlockSize() before this routine to set the stride
115: information, or use a vector created from a multicomponent DMDA.
117: If x is the array representing the vector x then this computes the norm
118: of the array (x[start],x[start+stride],x[start+2*stride], ....)
120: This is useful for computing, say the norm of the pressure variable when
121: the pressure is stored (interlaced) with other variables, say density etc.
123: This will only work if the desire subvector is a stride subvector
125: Level: advanced
128: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
129: @*/
130: PetscErrorCode VecStrideNorm(Vec v,PetscInt start,NormType ntype,PetscReal *nrm)
131: {
132: PetscErrorCode ierr;
133: PetscInt i,n,bs;
134: const PetscScalar *x;
135: PetscReal tnorm;
136: MPI_Comm comm;
141: VecGetLocalSize(v,&n);
142: VecGetArrayRead(v,&x);
143: PetscObjectGetComm((PetscObject)v,&comm);
145: VecGetBlockSize(v,&bs);
146: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
147: 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);
148: x += start;
150: if (ntype == NORM_2) {
151: PetscScalar sum = 0.0;
152: for (i=0; i<n; i+=bs) sum += x[i]*(PetscConj(x[i]));
153: tnorm = PetscRealPart(sum);
154: MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_SUM,comm);
155: *nrm = PetscSqrtReal(*nrm);
156: } else if (ntype == NORM_1) {
157: tnorm = 0.0;
158: for (i=0; i<n; i+=bs) tnorm += PetscAbsScalar(x[i]);
159: MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_SUM,comm);
160: } else if (ntype == NORM_INFINITY) {
161: PetscReal tmp;
162: tnorm = 0.0;
164: for (i=0; i<n; i+=bs) {
165: if ((tmp = PetscAbsScalar(x[i])) > tnorm) tnorm = tmp;
166: /* check special case of tmp == NaN */
167: if (tmp != tmp) {tnorm = tmp; break;}
168: }
169: MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_MAX,comm);
170: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
171: VecRestoreArrayRead(v,&x);
172: return(0);
173: }
175: /*@
176: VecStrideMax - Computes the maximum of subvector of a vector defined
177: by a starting point and a stride and optionally its location.
179: Collective on Vec
181: Input Parameter:
182: + v - the vector
183: - start - starting point of the subvector (defined by a stride)
185: Output Parameter:
186: + index - the location where the maximum occurred (pass NULL if not required)
187: - nrm - the maximum value in the subvector
189: Notes:
190: One must call VecSetBlockSize() before this routine to set the stride
191: information, or use a vector created from a multicomponent DMDA.
193: If xa is the array representing the vector x, then this computes the max
194: of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)
196: This is useful for computing, say the maximum of the pressure variable when
197: the pressure is stored (interlaced) with other variables, e.g., density, etc.
198: This will only work if the desire subvector is a stride subvector.
200: Level: advanced
203: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
204: @*/
205: PetscErrorCode VecStrideMax(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
206: {
207: PetscErrorCode ierr;
208: PetscInt i,n,bs,id;
209: const PetscScalar *x;
210: PetscReal max,tmp;
211: MPI_Comm comm;
217: VecGetLocalSize(v,&n);
218: VecGetArrayRead(v,&x);
219: PetscObjectGetComm((PetscObject)v,&comm);
221: VecGetBlockSize(v,&bs);
222: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
223: 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);
224: x += start;
226: id = -1;
227: if (!n) max = PETSC_MIN_REAL;
228: else {
229: id = 0;
230: max = PetscRealPart(x[0]);
231: for (i=bs; i<n; i+=bs) {
232: if ((tmp = PetscRealPart(x[i])) > max) { max = tmp; id = i;}
233: }
234: }
235: VecRestoreArrayRead(v,&x);
237: if (!idex) {
238: MPIU_Allreduce(&max,nrm,1,MPIU_REAL,MPIU_MAX,comm);
239: } else {
240: PetscReal in[2],out[2];
241: PetscInt rstart;
243: VecGetOwnershipRange(v,&rstart,NULL);
244: in[0] = max;
245: in[1] = rstart+id+start;
246: MPIU_Allreduce(in,out,2,MPIU_REAL,MPIU_MAXINDEX_OP,PetscObjectComm((PetscObject)v));
247: *nrm = out[0];
248: *idex = (PetscInt)out[1];
249: }
250: return(0);
251: }
253: /*@
254: VecStrideMin - Computes the minimum of subvector of a vector defined
255: by a starting point and a stride and optionally its location.
257: Collective on Vec
259: Input Parameter:
260: + v - the vector
261: - start - starting point of the subvector (defined by a stride)
263: Output Parameter:
264: + idex - the location where the minimum occurred. (pass NULL if not required)
265: - nrm - the minimum value in the subvector
267: Level: advanced
269: Notes:
270: One must call VecSetBlockSize() before this routine to set the stride
271: information, or use a vector created from a multicomponent DMDA.
273: If xa is the array representing the vector x, then this computes the min
274: of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)
276: This is useful for computing, say the minimum of the pressure variable when
277: the pressure is stored (interlaced) with other variables, e.g., density, etc.
278: This will only work if the desire subvector is a stride subvector.
281: .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
282: @*/
283: PetscErrorCode VecStrideMin(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
284: {
285: PetscErrorCode ierr;
286: PetscInt i,n,bs,id;
287: const PetscScalar *x;
288: PetscReal min,tmp;
289: MPI_Comm comm;
295: VecGetLocalSize(v,&n);
296: VecGetArrayRead(v,&x);
297: PetscObjectGetComm((PetscObject)v,&comm);
299: VecGetBlockSize(v,&bs);
300: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
301: 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);
302: x += start;
304: id = -1;
305: if (!n) min = PETSC_MAX_REAL;
306: else {
307: id = 0;
308: min = PetscRealPart(x[0]);
309: for (i=bs; i<n; i+=bs) {
310: if ((tmp = PetscRealPart(x[i])) < min) { min = tmp; id = i;}
311: }
312: }
313: VecRestoreArrayRead(v,&x);
315: if (!idex) {
316: MPIU_Allreduce(&min,nrm,1,MPIU_REAL,MPIU_MIN,comm);
317: } else {
318: PetscReal in[2],out[2];
319: PetscInt rstart;
321: VecGetOwnershipRange(v,&rstart,NULL);
322: in[0] = min;
323: in[1] = rstart+id;
324: MPIU_Allreduce(in,out,2,MPIU_REAL,MPIU_MININDEX_OP,PetscObjectComm((PetscObject)v));
325: *nrm = out[0];
326: *idex = (PetscInt)out[1];
327: }
328: return(0);
329: }
331: /*@
332: VecStrideScaleAll - Scales the subvectors of a vector defined
333: by a starting point and a stride.
335: Logically Collective on Vec
337: Input Parameter:
338: + v - the vector
339: - scales - values to multiply each subvector entry by
341: Notes:
342: One must call VecSetBlockSize() before this routine to set the stride
343: information, or use a vector created from a multicomponent DMDA.
345: The dimension of scales must be the same as the vector block size
348: Level: advanced
351: .seealso: VecNorm(), VecStrideScale(), VecScale(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
352: @*/
353: PetscErrorCode VecStrideScaleAll(Vec v,const PetscScalar *scales)
354: {
356: PetscInt i,j,n,bs;
357: PetscScalar *x;
362: VecSetErrorIfLocked(v,1);
364: VecGetLocalSize(v,&n);
365: VecGetArray(v,&x);
366: VecGetBlockSize(v,&bs);
368: /* need to provide optimized code for each bs */
369: for (i=0; i<n; i+=bs) {
370: for (j=0; j<bs; j++) x[i+j] *= scales[j];
371: }
372: VecRestoreArray(v,&x);
373: return(0);
374: }
376: /*@
377: VecStrideNormAll - Computes the norms of subvectors of a vector defined
378: by a starting point and a stride.
380: Collective on Vec
382: Input Parameter:
383: + v - the vector
384: - ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY
386: Output Parameter:
387: . nrm - the norms
389: Notes:
390: One must call VecSetBlockSize() before this routine to set the stride
391: information, or use a vector created from a multicomponent DMDA.
393: If x is the array representing the vector x then this computes the norm
394: of the array (x[start],x[start+stride],x[start+2*stride], ....) for each start < stride
396: The dimension of nrm must be the same as the vector block size
398: This will only work if the desire subvector is a stride subvector
400: Level: advanced
403: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
404: @*/
405: PetscErrorCode VecStrideNormAll(Vec v,NormType ntype,PetscReal nrm[])
406: {
407: PetscErrorCode ierr;
408: PetscInt i,j,n,bs;
409: const PetscScalar *x;
410: PetscReal tnorm[128];
411: MPI_Comm comm;
416: VecGetLocalSize(v,&n);
417: VecGetArrayRead(v,&x);
418: PetscObjectGetComm((PetscObject)v,&comm);
420: VecGetBlockSize(v,&bs);
421: if (bs > 128) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");
423: if (ntype == NORM_2) {
424: PetscScalar sum[128];
425: for (j=0; j<bs; j++) sum[j] = 0.0;
426: for (i=0; i<n; i+=bs) {
427: for (j=0; j<bs; j++) sum[j] += x[i+j]*(PetscConj(x[i+j]));
428: }
429: for (j=0; j<bs; j++) tnorm[j] = PetscRealPart(sum[j]);
431: MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_SUM,comm);
432: for (j=0; j<bs; j++) nrm[j] = PetscSqrtReal(nrm[j]);
433: } else if (ntype == NORM_1) {
434: for (j=0; j<bs; j++) tnorm[j] = 0.0;
436: for (i=0; i<n; i+=bs) {
437: for (j=0; j<bs; j++) tnorm[j] += PetscAbsScalar(x[i+j]);
438: }
440: MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_SUM,comm);
441: } else if (ntype == NORM_INFINITY) {
442: PetscReal tmp;
443: for (j=0; j<bs; j++) tnorm[j] = 0.0;
445: for (i=0; i<n; i+=bs) {
446: for (j=0; j<bs; j++) {
447: if ((tmp = PetscAbsScalar(x[i+j])) > tnorm[j]) tnorm[j] = tmp;
448: /* check special case of tmp == NaN */
449: if (tmp != tmp) {tnorm[j] = tmp; break;}
450: }
451: }
452: MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_MAX,comm);
453: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
454: VecRestoreArrayRead(v,&x);
455: return(0);
456: }
458: /*@
459: VecStrideMaxAll - Computes the maximums of subvectors of a vector defined
460: by a starting point and a stride and optionally its location.
462: Collective on Vec
464: Input Parameter:
465: . v - the vector
467: Output Parameter:
468: + index - the location where the maximum occurred (not supported, pass NULL,
469: if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
470: - nrm - the maximum values of each subvector
472: Notes:
473: One must call VecSetBlockSize() before this routine to set the stride
474: information, or use a vector created from a multicomponent DMDA.
476: The dimension of nrm must be the same as the vector block size
478: Level: advanced
481: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
482: @*/
483: PetscErrorCode VecStrideMaxAll(Vec v,PetscInt idex[],PetscReal nrm[])
484: {
485: PetscErrorCode ierr;
486: PetscInt i,j,n,bs;
487: const PetscScalar *x;
488: PetscReal max[128],tmp;
489: MPI_Comm comm;
494: 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");
495: VecGetLocalSize(v,&n);
496: VecGetArrayRead(v,&x);
497: PetscObjectGetComm((PetscObject)v,&comm);
499: VecGetBlockSize(v,&bs);
500: if (bs > 128) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");
502: if (!n) {
503: for (j=0; j<bs; j++) max[j] = PETSC_MIN_REAL;
504: } else {
505: for (j=0; j<bs; j++) max[j] = PetscRealPart(x[j]);
507: for (i=bs; i<n; i+=bs) {
508: for (j=0; j<bs; j++) {
509: if ((tmp = PetscRealPart(x[i+j])) > max[j]) max[j] = tmp;
510: }
511: }
512: }
513: MPIU_Allreduce(max,nrm,bs,MPIU_REAL,MPIU_MAX,comm);
515: VecRestoreArrayRead(v,&x);
516: return(0);
517: }
519: /*@
520: VecStrideMinAll - Computes the minimum of subvector of a vector defined
521: by a starting point and a stride and optionally its location.
523: Collective on Vec
525: Input Parameter:
526: . v - the vector
528: Output Parameter:
529: + idex - the location where the minimum occurred (not supported, pass NULL,
530: if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
531: - nrm - the minimums of each subvector
533: Level: advanced
535: Notes:
536: One must call VecSetBlockSize() before this routine to set the stride
537: information, or use a vector created from a multicomponent DMDA.
539: The dimension of nrm must be the same as the vector block size
542: .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
543: @*/
544: PetscErrorCode VecStrideMinAll(Vec v,PetscInt idex[],PetscReal nrm[])
545: {
546: PetscErrorCode ierr;
547: PetscInt i,n,bs,j;
548: const PetscScalar *x;
549: PetscReal min[128],tmp;
550: MPI_Comm comm;
555: 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");
556: VecGetLocalSize(v,&n);
557: VecGetArrayRead(v,&x);
558: PetscObjectGetComm((PetscObject)v,&comm);
560: VecGetBlockSize(v,&bs);
561: if (bs > 128) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");
563: if (!n) {
564: for (j=0; j<bs; j++) min[j] = PETSC_MAX_REAL;
565: } else {
566: for (j=0; j<bs; j++) min[j] = PetscRealPart(x[j]);
568: for (i=bs; i<n; i+=bs) {
569: for (j=0; j<bs; j++) {
570: if ((tmp = PetscRealPart(x[i+j])) < min[j]) min[j] = tmp;
571: }
572: }
573: }
574: MPIU_Allreduce(min,nrm,bs,MPIU_REAL,MPIU_MIN,comm);
576: VecRestoreArrayRead(v,&x);
577: return(0);
578: }
580: /*----------------------------------------------------------------------------------------------*/
581: /*@
582: VecStrideGatherAll - Gathers all the single components from a multi-component vector into
583: separate vectors.
585: Collective on Vec
587: Input Parameter:
588: + v - the vector
589: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
591: Output Parameter:
592: . s - the location where the subvectors are stored
594: Notes:
595: One must call VecSetBlockSize() before this routine to set the stride
596: information, or use a vector created from a multicomponent DMDA.
598: If x is the array representing the vector x then this gathers
599: the arrays (x[start],x[start+stride],x[start+2*stride], ....)
600: for start=0,1,2,...bs-1
602: The parallel layout of the vector and the subvector must be the same;
603: i.e., nlocal of v = stride*(nlocal of s)
605: Not optimized; could be easily
607: Level: advanced
609: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
610: VecStrideScatterAll()
611: @*/
612: PetscErrorCode VecStrideGatherAll(Vec v,Vec s[],InsertMode addv)
613: {
614: PetscErrorCode ierr;
615: PetscInt i,n,n2,bs,j,k,*bss = NULL,nv,jj,nvc;
616: PetscScalar **y;
617: const PetscScalar *x;
623: VecGetLocalSize(v,&n);
624: VecGetLocalSize(s[0],&n2);
625: VecGetArrayRead(v,&x);
626: VecGetBlockSize(v,&bs);
627: if (bs <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Input vector does not have a valid blocksize set");
629: PetscMalloc2(bs,&y,bs,&bss);
630: nv = 0;
631: nvc = 0;
632: for (i=0; i<bs; i++) {
633: VecGetBlockSize(s[i],&bss[i]);
634: if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1 Re: [PETSC #8241] VecStrideGatherAll */
635: VecGetArray(s[i],&y[i]);
636: nvc += bss[i];
637: nv++;
638: if (nvc > bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of subvectors in subvectors > number of vectors in main vector");
639: if (nvc == bs) break;
640: }
642: n = n/bs;
644: jj = 0;
645: if (addv == INSERT_VALUES) {
646: for (j=0; j<nv; j++) {
647: for (k=0; k<bss[j]; k++) {
648: for (i=0; i<n; i++) y[j][i*bss[j] + k] = x[bs*i+jj+k];
649: }
650: jj += bss[j];
651: }
652: } else if (addv == ADD_VALUES) {
653: for (j=0; j<nv; j++) {
654: for (k=0; k<bss[j]; k++) {
655: for (i=0; i<n; i++) y[j][i*bss[j] + k] += x[bs*i+jj+k];
656: }
657: jj += bss[j];
658: }
659: #if !defined(PETSC_USE_COMPLEX)
660: } else if (addv == MAX_VALUES) {
661: for (j=0; j<nv; j++) {
662: for (k=0; k<bss[j]; k++) {
663: for (i=0; i<n; i++) y[j][i*bss[j] + k] = PetscMax(y[j][i*bss[j] + k],x[bs*i+jj+k]);
664: }
665: jj += bss[j];
666: }
667: #endif
668: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
670: VecRestoreArrayRead(v,&x);
671: for (i=0; i<nv; i++) {
672: VecRestoreArray(s[i],&y[i]);
673: }
675: PetscFree2(y,bss);
676: return(0);
677: }
679: /*@
680: VecStrideScatterAll - Scatters all the single components from separate vectors into
681: a multi-component vector.
683: Collective on Vec
685: Input Parameter:
686: + s - the location where the subvectors are stored
687: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
689: Output Parameter:
690: . v - the multicomponent vector
692: Notes:
693: One must call VecSetBlockSize() before this routine to set the stride
694: information, or use a vector created from a multicomponent DMDA.
696: The parallel layout of the vector and the subvector must be the same;
697: i.e., nlocal of v = stride*(nlocal of s)
699: Not optimized; could be easily
701: Level: advanced
703: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
704: VecStrideScatterAll()
705: @*/
706: PetscErrorCode VecStrideScatterAll(Vec s[],Vec v,InsertMode addv)
707: {
708: PetscErrorCode ierr;
709: PetscInt i,n,n2,bs,j,jj,k,*bss = NULL,nv,nvc;
710: PetscScalar *x;
711: PetscScalar const **y;
717: VecGetLocalSize(v,&n);
718: VecGetLocalSize(s[0],&n2);
719: VecGetArray(v,&x);
720: VecGetBlockSize(v,&bs);
721: if (bs <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Input vector does not have a valid blocksize set");
723: PetscMalloc2(bs,(PetscScalar***)&y,bs,&bss);
724: nv = 0;
725: nvc = 0;
726: for (i=0; i<bs; i++) {
727: VecGetBlockSize(s[i],&bss[i]);
728: if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1 Re: [PETSC #8241] VecStrideGatherAll */
729: VecGetArrayRead(s[i],&y[i]);
730: nvc += bss[i];
731: nv++;
732: if (nvc > bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of subvectors in subvectors > number of vectors in main vector");
733: if (nvc == bs) break;
734: }
736: n = n/bs;
738: jj = 0;
739: if (addv == INSERT_VALUES) {
740: for (j=0; j<nv; j++) {
741: for (k=0; k<bss[j]; k++) {
742: for (i=0; i<n; i++) x[bs*i+jj+k] = y[j][i*bss[j] + k];
743: }
744: jj += bss[j];
745: }
746: } else if (addv == ADD_VALUES) {
747: for (j=0; j<nv; j++) {
748: for (k=0; k<bss[j]; k++) {
749: for (i=0; i<n; i++) x[bs*i+jj+k] += y[j][i*bss[j] + k];
750: }
751: jj += bss[j];
752: }
753: #if !defined(PETSC_USE_COMPLEX)
754: } else if (addv == MAX_VALUES) {
755: for (j=0; j<nv; j++) {
756: for (k=0; k<bss[j]; k++) {
757: for (i=0; i<n; i++) x[bs*i+jj+k] = PetscMax(x[bs*i+jj+k],y[j][i*bss[j] + k]);
758: }
759: jj += bss[j];
760: }
761: #endif
762: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
764: VecRestoreArray(v,&x);
765: for (i=0; i<nv; i++) {
766: VecRestoreArrayRead(s[i],&y[i]);
767: }
768: PetscFree2(*(PetscScalar***)&y,bss);
769: return(0);
770: }
772: /*@
773: VecStrideGather - Gathers a single component from a multi-component vector into
774: another vector.
776: Collective on Vec
778: Input Parameter:
779: + v - the vector
780: . start - starting point of the subvector (defined by a stride)
781: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
783: Output Parameter:
784: . s - the location where the subvector is stored
786: Notes:
787: One must call VecSetBlockSize() before this routine to set the stride
788: information, or use a vector created from a multicomponent DMDA.
790: If x is the array representing the vector x then this gathers
791: the array (x[start],x[start+stride],x[start+2*stride], ....)
793: The parallel layout of the vector and the subvector must be the same;
794: i.e., nlocal of v = stride*(nlocal of s)
796: Not optimized; could be easily
798: Level: advanced
800: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
801: VecStrideScatterAll()
802: @*/
803: PetscErrorCode VecStrideGather(Vec v,PetscInt start,Vec s,InsertMode addv)
804: {
810: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
811: 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);
812: if (!v->ops->stridegather) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
813: (*v->ops->stridegather)(v,start,s,addv);
814: return(0);
815: }
817: /*@
818: VecStrideScatter - Scatters a single component from a vector into a multi-component vector.
820: Collective on Vec
822: Input Parameter:
823: + s - the single-component vector
824: . start - starting point of the subvector (defined by a stride)
825: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
827: Output Parameter:
828: . v - the location where the subvector is scattered (the multi-component vector)
830: Notes:
831: One must call VecSetBlockSize() on the multi-component vector before this
832: routine to set the stride information, or use a vector created from a multicomponent DMDA.
834: The parallel layout of the vector and the subvector must be the same;
835: i.e., nlocal of v = stride*(nlocal of s)
837: Not optimized; could be easily
839: Level: advanced
841: .seealso: VecStrideNorm(), VecStrideGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
842: VecStrideScatterAll(), VecStrideSubSetScatter(), VecStrideSubSetGather()
843: @*/
844: PetscErrorCode VecStrideScatter(Vec s,PetscInt start,Vec v,InsertMode addv)
845: {
851: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
852: 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);
853: if (!v->ops->stridescatter) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
854: (*v->ops->stridescatter)(s,start,v,addv);
855: return(0);
856: }
858: /*@
859: VecStrideSubSetGather - Gathers a subset of components from a multi-component vector into
860: another vector.
862: Collective on Vec
864: Input Parameter:
865: + v - the vector
866: . nidx - the number of indices
867: . idxv - the indices of the components 0 <= idxv[0] ...idxv[nidx-1] < bs(v), they need not be sorted
868: . 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
869: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
871: Output Parameter:
872: . s - the location where the subvector is stored
874: Notes:
875: One must call VecSetBlockSize() on both vectors before this routine to set the stride
876: information, or use a vector created from a multicomponent DMDA.
879: The parallel layout of the vector and the subvector must be the same;
881: Not optimized; could be easily
883: Level: advanced
885: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideGather(), VecStrideSubSetScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
886: VecStrideScatterAll()
887: @*/
888: PetscErrorCode VecStrideSubSetGather(Vec v,PetscInt nidx,const PetscInt idxv[],const PetscInt idxs[],Vec s,InsertMode addv)
889: {
895: if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
896: if (!v->ops->stridesubsetgather) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
897: (*v->ops->stridesubsetgather)(v,nidx,idxv,idxs,s,addv);
898: return(0);
899: }
901: /*@
902: VecStrideSubSetScatter - Scatters components from a vector into a subset of components of a multi-component vector.
904: Collective on Vec
906: Input Parameter:
907: + s - the smaller-component vector
908: . nidx - the number of indices in idx
909: . 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
910: . idxv - the indices of the components in the larger-component vector, 0 <= idx[0] ...idx[nidx-1] < bs(v) they need not be sorted
911: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
913: Output Parameter:
914: . v - the location where the subvector is into scattered (the multi-component vector)
916: Notes:
917: One must call VecSetBlockSize() on the vectors before this
918: routine to set the stride information, or use a vector created from a multicomponent DMDA.
920: The parallel layout of the vector and the subvector must be the same;
922: Not optimized; could be easily
924: Level: advanced
926: .seealso: VecStrideNorm(), VecStrideGather(), VecStrideGather(), VecStrideSubSetGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
927: VecStrideScatterAll()
928: @*/
929: PetscErrorCode VecStrideSubSetScatter(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)
930: {
936: if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
937: if (!v->ops->stridesubsetscatter) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
938: (*v->ops->stridesubsetscatter)(s,nidx,idxs,idxv,v,addv);
939: return(0);
940: }
942: PetscErrorCode VecStrideGather_Default(Vec v,PetscInt start,Vec s,InsertMode addv)
943: {
945: PetscInt i,n,bs,ns;
946: const PetscScalar *x;
947: PetscScalar *y;
950: VecGetLocalSize(v,&n);
951: VecGetLocalSize(s,&ns);
952: VecGetArrayRead(v,&x);
953: VecGetArray(s,&y);
955: bs = v->map->bs;
956: 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);
957: x += start;
958: n = n/bs;
960: if (addv == INSERT_VALUES) {
961: for (i=0; i<n; i++) y[i] = x[bs*i];
962: } else if (addv == ADD_VALUES) {
963: for (i=0; i<n; i++) y[i] += x[bs*i];
964: #if !defined(PETSC_USE_COMPLEX)
965: } else if (addv == MAX_VALUES) {
966: for (i=0; i<n; i++) y[i] = PetscMax(y[i],x[bs*i]);
967: #endif
968: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
970: VecRestoreArrayRead(v,&x);
971: VecRestoreArray(s,&y);
972: return(0);
973: }
975: PetscErrorCode VecStrideScatter_Default(Vec s,PetscInt start,Vec v,InsertMode addv)
976: {
977: PetscErrorCode ierr;
978: PetscInt i,n,bs,ns;
979: PetscScalar *x;
980: const PetscScalar *y;
983: VecGetLocalSize(v,&n);
984: VecGetLocalSize(s,&ns);
985: VecGetArray(v,&x);
986: VecGetArrayRead(s,&y);
988: bs = v->map->bs;
989: 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);
990: x += start;
991: n = n/bs;
993: if (addv == INSERT_VALUES) {
994: for (i=0; i<n; i++) x[bs*i] = y[i];
995: } else if (addv == ADD_VALUES) {
996: for (i=0; i<n; i++) x[bs*i] += y[i];
997: #if !defined(PETSC_USE_COMPLEX)
998: } else if (addv == MAX_VALUES) {
999: for (i=0; i<n; i++) x[bs*i] = PetscMax(y[i],x[bs*i]);
1000: #endif
1001: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
1003: VecRestoreArray(v,&x);
1004: VecRestoreArrayRead(s,&y);
1005: return(0);
1006: }
1008: PetscErrorCode VecStrideSubSetGather_Default(Vec v,PetscInt nidx,const PetscInt idxv[],const PetscInt idxs[],Vec s,InsertMode addv)
1009: {
1010: PetscErrorCode ierr;
1011: PetscInt i,j,n,bs,bss,ns;
1012: const PetscScalar *x;
1013: PetscScalar *y;
1016: VecGetLocalSize(v,&n);
1017: VecGetLocalSize(s,&ns);
1018: VecGetArrayRead(v,&x);
1019: VecGetArray(s,&y);
1021: bs = v->map->bs;
1022: bss = s->map->bs;
1023: n = n/bs;
1025: if (PetscDefined(USE_DEBUG)) {
1026: if (n != ns/bss) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible layout of vectors");
1027: for (j=0; j<nidx; j++) {
1028: if (idxv[j] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is negative",j,idxv[j]);
1029: 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);
1030: }
1031: if (!idxs && bss != nidx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide idxs when not gathering into all locations");
1032: }
1034: if (addv == INSERT_VALUES) {
1035: if (!idxs) {
1036: for (i=0; i<n; i++) {
1037: for (j=0; j<bss; j++) y[bss*i+j] = x[bs*i+idxv[j]];
1038: }
1039: } else {
1040: for (i=0; i<n; i++) {
1041: for (j=0; j<bss; j++) y[bss*i+idxs[j]] = x[bs*i+idxv[j]];
1042: }
1043: }
1044: } else if (addv == ADD_VALUES) {
1045: if (!idxs) {
1046: for (i=0; i<n; i++) {
1047: for (j=0; j<bss; j++) y[bss*i+j] += x[bs*i+idxv[j]];
1048: }
1049: } else {
1050: for (i=0; i<n; i++) {
1051: for (j=0; j<bss; j++) y[bss*i+idxs[j]] += x[bs*i+idxv[j]];
1052: }
1053: }
1054: #if !defined(PETSC_USE_COMPLEX)
1055: } else if (addv == MAX_VALUES) {
1056: if (!idxs) {
1057: for (i=0; i<n; i++) {
1058: for (j=0; j<bss; j++) y[bss*i+j] = PetscMax(y[bss*i+j],x[bs*i+idxv[j]]);
1059: }
1060: } else {
1061: for (i=0; i<n; i++) {
1062: for (j=0; j<bss; j++) y[bss*i+idxs[j]] = PetscMax(y[bss*i+idxs[j]],x[bs*i+idxv[j]]);
1063: }
1064: }
1065: #endif
1066: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
1068: VecRestoreArrayRead(v,&x);
1069: VecRestoreArray(s,&y);
1070: return(0);
1071: }
1073: PetscErrorCode VecStrideSubSetScatter_Default(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)
1074: {
1075: PetscErrorCode ierr;
1076: PetscInt j,i,n,bs,ns,bss;
1077: PetscScalar *x;
1078: const PetscScalar *y;
1081: VecGetLocalSize(v,&n);
1082: VecGetLocalSize(s,&ns);
1083: VecGetArray(v,&x);
1084: VecGetArrayRead(s,&y);
1086: bs = v->map->bs;
1087: bss = s->map->bs;
1088: n = n/bs;
1090: if (PetscDefined(USE_DEBUG)) {
1091: if (n != ns/bss) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible layout of vectors");
1092: for (j=0; j<bss; j++) {
1093: if (idxs) {
1094: if (idxs[j] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is negative",j,idxs[j]);
1095: 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);
1096: }
1097: }
1098: if (!idxs && bss != nidx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide idxs when not scattering from all locations");
1099: }
1101: if (addv == INSERT_VALUES) {
1102: if (!idxs) {
1103: for (i=0; i<n; i++) {
1104: for (j=0; j<bss; j++) x[bs*i + idxv[j]] = y[bss*i+j];
1105: }
1106: } else {
1107: for (i=0; i<n; i++) {
1108: for (j=0; j<bss; j++) x[bs*i + idxv[j]] = y[bss*i+idxs[j]];
1109: }
1110: }
1111: } else if (addv == ADD_VALUES) {
1112: if (!idxs) {
1113: for (i=0; i<n; i++) {
1114: for (j=0; j<bss; j++) x[bs*i + idxv[j]] += y[bss*i+j];
1115: }
1116: } else {
1117: for (i=0; i<n; i++) {
1118: for (j=0; j<bss; j++) x[bs*i + idxv[j]] += y[bss*i+idxs[j]];
1119: }
1120: }
1121: #if !defined(PETSC_USE_COMPLEX)
1122: } else if (addv == MAX_VALUES) {
1123: if (!idxs) {
1124: for (i=0; i<n; i++) {
1125: for (j=0; j<bss; j++) x[bs*i + idxv[j]] = PetscMax(y[bss*i+j],x[bs*i + idxv[j]]);
1126: }
1127: } else {
1128: for (i=0; i<n; i++) {
1129: for (j=0; j<bss; j++) x[bs*i + idxv[j]] = PetscMax(y[bss*i+idxs[j]],x[bs*i + idxv[j]]);
1130: }
1131: }
1132: #endif
1133: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
1135: VecRestoreArray(v,&x);
1136: VecRestoreArrayRead(s,&y);
1137: return(0);
1138: }
1140: PetscErrorCode VecReciprocal_Default(Vec v)
1141: {
1143: PetscInt i,n;
1144: PetscScalar *x;
1147: VecGetLocalSize(v,&n);
1148: VecGetArray(v,&x);
1149: for (i=0; i<n; i++) {
1150: if (x[i] != (PetscScalar)0.0) x[i] = (PetscScalar)1.0/x[i];
1151: }
1152: VecRestoreArray(v,&x);
1153: return(0);
1154: }
1156: /*@
1157: VecExp - Replaces each component of a vector by e^x_i
1159: Not collective
1161: Input Parameter:
1162: . v - The vector
1164: Output Parameter:
1165: . v - The vector of exponents
1167: Level: beginner
1169: .seealso: VecLog(), VecAbs(), VecSqrtAbs(), VecReciprocal()
1171: @*/
1172: PetscErrorCode VecExp(Vec v)
1173: {
1174: PetscScalar *x;
1175: PetscInt i, n;
1180: if (v->ops->exp) {
1181: (*v->ops->exp)(v);
1182: } else {
1183: VecGetLocalSize(v, &n);
1184: VecGetArray(v, &x);
1185: for (i = 0; i < n; i++) x[i] = PetscExpScalar(x[i]);
1186: VecRestoreArray(v, &x);
1187: }
1188: return(0);
1189: }
1191: /*@
1192: VecLog - Replaces each component of a vector by log(x_i), the natural logarithm
1194: Not collective
1196: Input Parameter:
1197: . v - The vector
1199: Output Parameter:
1200: . v - The vector of logs
1202: Level: beginner
1204: .seealso: VecExp(), VecAbs(), VecSqrtAbs(), VecReciprocal()
1206: @*/
1207: PetscErrorCode VecLog(Vec v)
1208: {
1209: PetscScalar *x;
1210: PetscInt i, n;
1215: if (v->ops->log) {
1216: (*v->ops->log)(v);
1217: } else {
1218: VecGetLocalSize(v, &n);
1219: VecGetArray(v, &x);
1220: for (i = 0; i < n; i++) x[i] = PetscLogScalar(x[i]);
1221: VecRestoreArray(v, &x);
1222: }
1223: return(0);
1224: }
1226: /*@
1227: VecSqrtAbs - Replaces each component of a vector by the square root of its magnitude.
1229: Not collective
1231: Input Parameter:
1232: . v - The vector
1234: Output Parameter:
1235: . v - The vector square root
1237: Level: beginner
1239: Note: The actual function is sqrt(|x_i|)
1241: .seealso: VecLog(), VecExp(), VecReciprocal(), VecAbs()
1243: @*/
1244: PetscErrorCode VecSqrtAbs(Vec v)
1245: {
1246: PetscScalar *x;
1247: PetscInt i, n;
1252: if (v->ops->sqrt) {
1253: (*v->ops->sqrt)(v);
1254: } else {
1255: VecGetLocalSize(v, &n);
1256: VecGetArray(v, &x);
1257: for (i = 0; i < n; i++) x[i] = PetscSqrtReal(PetscAbsScalar(x[i]));
1258: VecRestoreArray(v, &x);
1259: }
1260: return(0);
1261: }
1263: /*@
1264: VecDotNorm2 - computes the inner product of two vectors and the 2-norm squared of the second vector
1266: Collective on Vec
1268: Input Parameter:
1269: + s - first vector
1270: - t - second vector
1272: Output Parameter:
1273: + dp - s'conj(t)
1274: - nm - t'conj(t)
1276: Level: advanced
1278: Notes:
1279: conj(x) is the complex conjugate of x when x is complex
1282: .seealso: VecDot(), VecNorm(), VecDotBegin(), VecNormBegin(), VecDotEnd(), VecNormEnd()
1284: @*/
1285: PetscErrorCode VecDotNorm2(Vec s,Vec t,PetscScalar *dp, PetscReal *nm)
1286: {
1287: const PetscScalar *sx, *tx;
1288: PetscScalar dpx = 0.0, nmx = 0.0,work[2],sum[2];
1289: PetscInt i, n;
1290: PetscErrorCode ierr;
1300: if (s->map->N != t->map->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
1301: if (s->map->n != t->map->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");
1303: PetscLogEventBegin(VEC_DotNorm2,s,t,0,0);
1304: if (s->ops->dotnorm2) {
1305: (*s->ops->dotnorm2)(s,t,dp,&dpx);
1306: *nm = PetscRealPart(dpx);
1307: } else {
1308: VecGetLocalSize(s, &n);
1309: VecGetArrayRead(s, &sx);
1310: VecGetArrayRead(t, &tx);
1312: for (i = 0; i<n; i++) {
1313: dpx += sx[i]*PetscConj(tx[i]);
1314: nmx += tx[i]*PetscConj(tx[i]);
1315: }
1316: work[0] = dpx;
1317: work[1] = nmx;
1319: MPIU_Allreduce(work,sum,2,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)s));
1320: *dp = sum[0];
1321: *nm = PetscRealPart(sum[1]);
1323: VecRestoreArrayRead(t, &tx);
1324: VecRestoreArrayRead(s, &sx);
1325: PetscLogFlops(4.0*n);
1326: }
1327: PetscLogEventEnd(VEC_DotNorm2,s,t,0,0);
1328: return(0);
1329: }
1331: /*@
1332: VecSum - Computes the sum of all the components of a vector.
1334: Collective on Vec
1336: Input Parameter:
1337: . v - the vector
1339: Output Parameter:
1340: . sum - the result
1342: Level: beginner
1344: .seealso: VecNorm()
1345: @*/
1346: PetscErrorCode VecSum(Vec v,PetscScalar *sum)
1347: {
1348: PetscErrorCode ierr;
1349: PetscInt i,n;
1350: const PetscScalar *x;
1351: PetscScalar lsum = 0.0;
1356: VecGetLocalSize(v,&n);
1357: VecGetArrayRead(v,&x);
1358: for (i=0; i<n; i++) lsum += x[i];
1359: MPIU_Allreduce(&lsum,sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)v));
1360: VecRestoreArrayRead(v,&x);
1361: return(0);
1362: }
1364: /*@
1365: VecImaginaryPart - Replaces a complex vector with its imginary part
1367: Collective on Vec
1369: Input Parameter:
1370: . v - the vector
1372: Level: beginner
1374: .seealso: VecNorm(), VecRealPart()
1375: @*/
1376: PetscErrorCode VecImaginaryPart(Vec v)
1377: {
1378: PetscErrorCode ierr;
1379: PetscInt i,n;
1380: PetscScalar *x;
1384: VecGetLocalSize(v,&n);
1385: VecGetArray(v,&x);
1386: for (i=0; i<n; i++) x[i] = PetscImaginaryPart(x[i]);
1387: VecRestoreArray(v,&x);
1388: return(0);
1389: }
1391: /*@
1392: VecRealPart - Replaces a complex vector with its real part
1394: Collective on Vec
1396: Input Parameter:
1397: . v - the vector
1399: Level: beginner
1401: .seealso: VecNorm(), VecImaginaryPart()
1402: @*/
1403: PetscErrorCode VecRealPart(Vec v)
1404: {
1405: PetscErrorCode ierr;
1406: PetscInt i,n;
1407: PetscScalar *x;
1411: VecGetLocalSize(v,&n);
1412: VecGetArray(v,&x);
1413: for (i=0; i<n; i++) x[i] = PetscRealPart(x[i]);
1414: VecRestoreArray(v,&x);
1415: return(0);
1416: }
1418: /*@
1419: VecShift - Shifts all of the components of a vector by computing
1420: x[i] = x[i] + shift.
1422: Logically Collective on Vec
1424: Input Parameters:
1425: + v - the vector
1426: - shift - the shift
1428: Output Parameter:
1429: . v - the shifted vector
1431: Level: intermediate
1433: @*/
1434: PetscErrorCode VecShift(Vec v,PetscScalar shift)
1435: {
1437: PetscInt i,n;
1438: PetscScalar *x;
1443: VecSetErrorIfLocked(v,1);
1444: if (shift == 0.0) return(0);
1446: if (v->ops->shift) {
1447: (*v->ops->shift)(v,shift);
1448: } else {
1449: VecGetLocalSize(v,&n);
1450: VecGetArray(v,&x);
1451: for (i=0; i<n; i++) x[i] += shift;
1452: VecRestoreArray(v,&x);
1453: }
1454: return(0);
1455: }
1457: /*@
1458: VecAbs - Replaces every element in a vector with its absolute value.
1460: Logically Collective on Vec
1462: Input Parameters:
1463: . v - the vector
1465: Level: intermediate
1467: @*/
1468: PetscErrorCode VecAbs(Vec v)
1469: {
1471: PetscInt i,n;
1472: PetscScalar *x;
1476: VecSetErrorIfLocked(v,1);
1478: if (v->ops->abs) {
1479: (*v->ops->abs)(v);
1480: } else {
1481: VecGetLocalSize(v,&n);
1482: VecGetArray(v,&x);
1483: for (i=0; i<n; i++) x[i] = PetscAbsScalar(x[i]);
1484: VecRestoreArray(v,&x);
1485: }
1486: return(0);
1487: }
1489: /*@
1490: VecPermute - Permutes a vector in place using the given ordering.
1492: Input Parameters:
1493: + vec - The vector
1494: . order - The ordering
1495: - inv - The flag for inverting the permutation
1497: Level: beginner
1499: Note: This function does not yet support parallel Index Sets with non-local permutations
1501: .seealso: MatPermute()
1502: @*/
1503: PetscErrorCode VecPermute(Vec x, IS row, PetscBool inv)
1504: {
1505: const PetscScalar *array;
1506: PetscScalar *newArray;
1507: const PetscInt *idx;
1508: PetscInt i,rstart,rend;
1509: PetscErrorCode ierr;
1512: VecSetErrorIfLocked(x,1);
1514: VecGetOwnershipRange(x,&rstart,&rend);
1515: ISGetIndices(row, &idx);
1516: VecGetArrayRead(x, &array);
1517: PetscMalloc1(x->map->n, &newArray);
1518: if (PetscDefined(USE_DEBUG)) {
1519: for (i = 0; i < x->map->n; i++) {
1520: 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]);
1521: }
1522: }
1523: if (!inv) {
1524: for (i = 0; i < x->map->n; i++) newArray[i] = array[idx[i]-rstart];
1525: } else {
1526: for (i = 0; i < x->map->n; i++) newArray[idx[i]-rstart] = array[i];
1527: }
1528: VecRestoreArrayRead(x, &array);
1529: ISRestoreIndices(row, &idx);
1530: VecReplaceArray(x, newArray);
1531: return(0);
1532: }
1534: /*@
1535: VecEqual - Compares two vectors. Returns true if the two vectors are either pointing to the same memory buffer,
1536: or if the two vectors have the same local and global layout as well as bitwise equality of all entries.
1537: Does NOT take round-off errors into account.
1539: Collective on Vec
1541: Input Parameters:
1542: + vec1 - the first vector
1543: - vec2 - the second vector
1545: Output Parameter:
1546: . flg - PETSC_TRUE if the vectors are equal; PETSC_FALSE otherwise.
1548: Level: intermediate
1551: @*/
1552: PetscErrorCode VecEqual(Vec vec1,Vec vec2,PetscBool *flg)
1553: {
1554: const PetscScalar *v1,*v2;
1555: PetscErrorCode ierr;
1556: PetscInt n1,n2,N1,N2;
1557: PetscBool flg1;
1563: if (vec1 == vec2) *flg = PETSC_TRUE;
1564: else {
1565: VecGetSize(vec1,&N1);
1566: VecGetSize(vec2,&N2);
1567: if (N1 != N2) flg1 = PETSC_FALSE;
1568: else {
1569: VecGetLocalSize(vec1,&n1);
1570: VecGetLocalSize(vec2,&n2);
1571: if (n1 != n2) flg1 = PETSC_FALSE;
1572: else {
1573: VecGetArrayRead(vec1,&v1);
1574: VecGetArrayRead(vec2,&v2);
1575: PetscArraycmp(v1,v2,n1,&flg1);
1576: VecRestoreArrayRead(vec1,&v1);
1577: VecRestoreArrayRead(vec2,&v2);
1578: }
1579: }
1580: /* combine results from all processors */
1581: MPIU_Allreduce(&flg1,flg,1,MPIU_BOOL,MPI_MIN,PetscObjectComm((PetscObject)vec1));
1582: }
1583: return(0);
1584: }
1586: /*@
1587: VecUniqueEntries - Compute the number of unique entries, and those entries
1589: Collective on Vec
1591: Input Parameter:
1592: . vec - the vector
1594: Output Parameters:
1595: + n - The number of unique entries
1596: - e - The entries
1598: Level: intermediate
1600: @*/
1601: PetscErrorCode VecUniqueEntries(Vec vec, PetscInt *n, PetscScalar **e)
1602: {
1603: const PetscScalar *v;
1604: PetscScalar *tmp, *vals;
1605: PetscMPIInt *N, *displs, l;
1606: PetscInt ng, m, i, j, p;
1607: PetscMPIInt size;
1608: PetscErrorCode ierr;
1613: MPI_Comm_size(PetscObjectComm((PetscObject) vec), &size);
1614: VecGetLocalSize(vec, &m);
1615: VecGetArrayRead(vec, &v);
1616: PetscMalloc2(m,&tmp,size,&N);
1617: for (i = 0, j = 0, l = 0; i < m; ++i) {
1618: /* Can speed this up with sorting */
1619: for (j = 0; j < l; ++j) {
1620: if (v[i] == tmp[j]) break;
1621: }
1622: if (j == l) {
1623: tmp[j] = v[i];
1624: ++l;
1625: }
1626: }
1627: VecRestoreArrayRead(vec, &v);
1628: /* Gather serial results */
1629: MPI_Allgather(&l, 1, MPI_INT, N, 1, MPI_INT, PetscObjectComm((PetscObject) vec));
1630: for (p = 0, ng = 0; p < size; ++p) {
1631: ng += N[p];
1632: }
1633: PetscMalloc2(ng,&vals,size+1,&displs);
1634: for (p = 1, displs[0] = 0; p <= size; ++p) {
1635: displs[p] = displs[p-1] + N[p-1];
1636: }
1637: MPI_Allgatherv(tmp, l, MPIU_SCALAR, vals, N, displs, MPIU_SCALAR, PetscObjectComm((PetscObject) vec));
1638: /* Find unique entries */
1639: #ifdef PETSC_USE_COMPLEX
1640: SETERRQ(PetscObjectComm((PetscObject) vec), PETSC_ERR_SUP, "Does not work with complex numbers");
1641: #else
1642: *n = displs[size];
1643: PetscSortRemoveDupsReal(n, (PetscReal *) vals);
1644: if (e) {
1646: PetscMalloc1(*n, e);
1647: for (i = 0; i < *n; ++i) {
1648: (*e)[i] = vals[i];
1649: }
1650: }
1651: PetscFree2(vals,displs);
1652: PetscFree2(tmp,N);
1653: return(0);
1654: #endif
1655: }