Actual source code: vinv.c
petsc-3.5.4 2015-05-23
2: /*
3: Some useful vector utility functions.
4: */
5: #include <../src/vec/vec/impls/mpi/pvecimpl.h> /*I "petscvec.h" I*/
6: extern MPI_Op VecMax_Local_Op;
7: extern MPI_Op VecMin_Local_Op;
11: /*@
12: VecStrideSet - Sets a subvector of a vector defined
13: by a starting point and a stride with a given value
15: Logically Collective on Vec
17: Input Parameter:
18: + v - the vector
19: . start - starting point of the subvector (defined by a stride)
20: - s - value to multiply each subvector entry by
22: Notes:
23: One must call VecSetBlockSize() before this routine to set the stride
24: information, or use a vector created from a multicomponent DMDA.
26: This will only work if the desire subvector is a stride subvector
28: Level: advanced
30: Concepts: scale^on stride of vector
31: Concepts: stride^scale
33: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideScale()
34: @*/
35: PetscErrorCode VecStrideSet(Vec v,PetscInt start,PetscScalar s)
36: {
38: PetscInt i,n,bs;
39: PetscScalar *x;
46: VecGetLocalSize(v,&n);
47: VecGetArray(v,&x);
48: VecGetBlockSize(v,&bs);
49: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
50: else if (start >= bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Start of stride subvector (%D) is too large for stride\n Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
51: x += start;
53: for (i=0; i<n; i+=bs) x[i] = s;
54: x -= start;
56: VecRestoreArray(v,&x);
57: return(0);
58: }
62: /*@
63: VecStrideScale - Scales a subvector of a vector defined
64: by a starting point and a stride.
66: Logically Collective on Vec
68: Input Parameter:
69: + v - the vector
70: . start - starting point of the subvector (defined by a stride)
71: - scale - value to multiply each subvector entry by
73: Notes:
74: One must call VecSetBlockSize() before this routine to set the stride
75: information, or use a vector created from a multicomponent DMDA.
77: This will only work if the desire subvector is a stride subvector
79: Level: advanced
81: Concepts: scale^on stride of vector
82: Concepts: stride^scale
84: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideScale()
85: @*/
86: PetscErrorCode VecStrideScale(Vec v,PetscInt start,PetscScalar scale)
87: {
89: PetscInt i,n,bs;
90: PetscScalar *x;
97: VecGetLocalSize(v,&n);
98: VecGetArray(v,&x);
99: VecGetBlockSize(v,&bs);
100: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
101: else if (start >= bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Start of stride subvector (%D) is too large for stride\n Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
102: x += start;
104: for (i=0; i<n; i+=bs) x[i] *= scale;
105: x -= start;
107: VecRestoreArray(v,&x);
108: return(0);
109: }
113: /*@
114: VecStrideNorm - Computes the norm of subvector of a vector defined
115: by a starting point and a stride.
117: Collective on Vec
119: Input Parameter:
120: + v - the vector
121: . start - starting point of the subvector (defined by a stride)
122: - ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY
124: Output Parameter:
125: . norm - the norm
127: Notes:
128: One must call VecSetBlockSize() before this routine to set the stride
129: information, or use a vector created from a multicomponent DMDA.
131: If x is the array representing the vector x then this computes the norm
132: of the array (x[start],x[start+stride],x[start+2*stride], ....)
134: This is useful for computing, say the norm of the pressure variable when
135: the pressure is stored (interlaced) with other variables, say density etc.
137: This will only work if the desire subvector is a stride subvector
139: Level: advanced
141: Concepts: norm^on stride of vector
142: Concepts: stride^norm
144: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
145: @*/
146: PetscErrorCode VecStrideNorm(Vec v,PetscInt start,NormType ntype,PetscReal *nrm)
147: {
149: PetscInt i,n,bs;
150: PetscScalar *x;
151: PetscReal tnorm;
152: MPI_Comm comm;
157: VecGetLocalSize(v,&n);
158: VecGetArray(v,&x);
159: PetscObjectGetComm((PetscObject)v,&comm);
161: VecGetBlockSize(v,&bs);
162: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
163: else if (start >= bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Start of stride subvector (%D) is too large for stride\n Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
164: x += start;
166: if (ntype == NORM_2) {
167: PetscScalar sum = 0.0;
168: for (i=0; i<n; i+=bs) sum += x[i]*(PetscConj(x[i]));
169: tnorm = PetscRealPart(sum);
170: MPI_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_SUM,comm);
171: *nrm = PetscSqrtReal(*nrm);
172: } else if (ntype == NORM_1) {
173: tnorm = 0.0;
174: for (i=0; i<n; i+=bs) tnorm += PetscAbsScalar(x[i]);
175: MPI_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_SUM,comm);
176: } else if (ntype == NORM_INFINITY) {
177: PetscReal tmp;
178: tnorm = 0.0;
180: for (i=0; i<n; i+=bs) {
181: if ((tmp = PetscAbsScalar(x[i])) > tnorm) tnorm = tmp;
182: /* check special case of tmp == NaN */
183: if (tmp != tmp) {tnorm = tmp; break;}
184: }
185: MPI_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_MAX,comm);
186: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
187: VecRestoreArray(v,&x);
188: return(0);
189: }
193: /*@
194: VecStrideMax - Computes the maximum of subvector of a vector defined
195: by a starting point and a stride and optionally its location.
197: Collective on Vec
199: Input Parameter:
200: + v - the vector
201: - start - starting point of the subvector (defined by a stride)
203: Output Parameter:
204: + index - the location where the maximum occurred (pass NULL if not required)
205: - nrm - the max
207: Notes:
208: One must call VecSetBlockSize() before this routine to set the stride
209: information, or use a vector created from a multicomponent DMDA.
211: If xa is the array representing the vector x, then this computes the max
212: of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)
214: This is useful for computing, say the maximum of the pressure variable when
215: the pressure is stored (interlaced) with other variables, e.g., density, etc.
216: This will only work if the desire subvector is a stride subvector.
218: Level: advanced
220: Concepts: maximum^on stride of vector
221: Concepts: stride^maximum
223: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
224: @*/
225: PetscErrorCode VecStrideMax(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
226: {
228: PetscInt i,n,bs,id;
229: PetscScalar *x;
230: PetscReal max,tmp;
231: MPI_Comm comm;
237: VecGetLocalSize(v,&n);
238: VecGetArray(v,&x);
239: PetscObjectGetComm((PetscObject)v,&comm);
241: VecGetBlockSize(v,&bs);
242: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
243: else if (start >= bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Start of stride subvector (%D) is too large for stride\n Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
244: x += start;
246: id = -1;
247: if (!n) max = PETSC_MIN_REAL;
248: else {
249: id = 0;
250: max = PetscRealPart(x[0]);
251: for (i=bs; i<n; i+=bs) {
252: if ((tmp = PetscRealPart(x[i])) > max) { max = tmp; id = i;}
253: }
254: }
255: VecRestoreArray(v,&x);
257: if (!idex) {
258: MPI_Allreduce(&max,nrm,1,MPIU_REAL,MPIU_MAX,comm);
259: } else {
260: PetscReal in[2],out[2];
261: PetscInt rstart;
263: VecGetOwnershipRange(v,&rstart,NULL);
264: in[0] = max;
265: in[1] = rstart+id+start;
266: MPI_Allreduce(in,out,2,MPIU_REAL,VecMax_Local_Op,PetscObjectComm((PetscObject)v));
267: *nrm = out[0];
268: *idex = (PetscInt)out[1];
269: }
270: return(0);
271: }
275: /*@
276: VecStrideMin - Computes the minimum of subvector of a vector defined
277: by a starting point and a stride and optionally its location.
279: Collective on Vec
281: Input Parameter:
282: + v - the vector
283: - start - starting point of the subvector (defined by a stride)
285: Output Parameter:
286: + idex - the location where the minimum occurred. (pass NULL if not required)
287: - nrm - the min
289: Level: advanced
291: Notes:
292: One must call VecSetBlockSize() before this routine to set the stride
293: information, or use a vector created from a multicomponent DMDA.
295: If xa is the array representing the vector x, then this computes the min
296: of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)
298: This is useful for computing, say the minimum of the pressure variable when
299: the pressure is stored (interlaced) with other variables, e.g., density, etc.
300: This will only work if the desire subvector is a stride subvector.
302: Concepts: minimum^on stride of vector
303: Concepts: stride^minimum
305: .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
306: @*/
307: PetscErrorCode VecStrideMin(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
308: {
310: PetscInt i,n,bs,id;
311: PetscScalar *x;
312: PetscReal min,tmp;
313: MPI_Comm comm;
319: VecGetLocalSize(v,&n);
320: VecGetArray(v,&x);
321: PetscObjectGetComm((PetscObject)v,&comm);
323: VecGetBlockSize(v,&bs);
324: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
325: else if (start >= bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Start of stride subvector (%D) is too large for stride\nHave you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
326: x += start;
328: id = -1;
329: if (!n) min = PETSC_MAX_REAL;
330: else {
331: id = 0;
332: min = PetscRealPart(x[0]);
333: for (i=bs; i<n; i+=bs) {
334: if ((tmp = PetscRealPart(x[i])) < min) { min = tmp; id = i;}
335: }
336: }
337: VecRestoreArray(v,&x);
339: if (!idex) {
340: MPI_Allreduce(&min,nrm,1,MPIU_REAL,MPIU_MIN,comm);
341: } else {
342: PetscReal in[2],out[2];
343: PetscInt rstart;
345: VecGetOwnershipRange(v,&rstart,NULL);
346: in[0] = min;
347: in[1] = rstart+id;
348: MPI_Allreduce(in,out,2,MPIU_REAL,VecMin_Local_Op,PetscObjectComm((PetscObject)v));
349: *nrm = out[0];
350: *idex = (PetscInt)out[1];
351: }
352: return(0);
353: }
357: /*@
358: VecStrideScaleAll - Scales the subvectors of a vector defined
359: by a starting point and a stride.
361: Logically Collective on Vec
363: Input Parameter:
364: + v - the vector
365: - scales - values to multiply each subvector entry by
367: Notes:
368: One must call VecSetBlockSize() before this routine to set the stride
369: information, or use a vector created from a multicomponent DMDA.
372: Level: advanced
374: Concepts: scale^on stride of vector
375: Concepts: stride^scale
377: .seealso: VecNorm(), VecStrideScale(), VecScale(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
378: @*/
379: PetscErrorCode VecStrideScaleAll(Vec v,const PetscScalar *scales)
380: {
382: PetscInt i,j,n,bs;
383: PetscScalar *x;
388: VecGetLocalSize(v,&n);
389: VecGetArray(v,&x);
390: VecGetBlockSize(v,&bs);
392: /* need to provide optimized code for each bs */
393: for (i=0; i<n; i+=bs) {
394: for (j=0; j<bs; j++) x[i+j] *= scales[j];
395: }
396: VecRestoreArray(v,&x);
397: return(0);
398: }
402: /*@
403: VecStrideNormAll - Computes the norms subvectors of a vector defined
404: by a starting point and a stride.
406: Collective on Vec
408: Input Parameter:
409: + v - the vector
410: - ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY
412: Output Parameter:
413: . nrm - the norms
415: Notes:
416: One must call VecSetBlockSize() before this routine to set the stride
417: information, or use a vector created from a multicomponent DMDA.
419: If x is the array representing the vector x then this computes the norm
420: of the array (x[start],x[start+stride],x[start+2*stride], ....)
422: This is useful for computing, say the norm of the pressure variable when
423: the pressure is stored (interlaced) with other variables, say density etc.
425: This will only work if the desire subvector is a stride subvector
427: Level: advanced
429: Concepts: norm^on stride of vector
430: Concepts: stride^norm
432: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
433: @*/
434: PetscErrorCode VecStrideNormAll(Vec v,NormType ntype,PetscReal nrm[])
435: {
437: PetscInt i,j,n,bs;
438: PetscScalar *x;
439: PetscReal tnorm[128];
440: MPI_Comm comm;
445: VecGetLocalSize(v,&n);
446: VecGetArray(v,&x);
447: PetscObjectGetComm((PetscObject)v,&comm);
449: VecGetBlockSize(v,&bs);
450: if (bs > 128) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");
452: if (ntype == NORM_2) {
453: PetscScalar sum[128];
454: for (j=0; j<bs; j++) sum[j] = 0.0;
455: for (i=0; i<n; i+=bs) {
456: for (j=0; j<bs; j++) sum[j] += x[i+j]*(PetscConj(x[i+j]));
457: }
458: for (j=0; j<bs; j++) tnorm[j] = PetscRealPart(sum[j]);
460: MPI_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_SUM,comm);
461: for (j=0; j<bs; j++) nrm[j] = PetscSqrtReal(nrm[j]);
462: } else if (ntype == NORM_1) {
463: for (j=0; j<bs; j++) tnorm[j] = 0.0;
465: for (i=0; i<n; i+=bs) {
466: for (j=0; j<bs; j++) tnorm[j] += PetscAbsScalar(x[i+j]);
467: }
469: MPI_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_SUM,comm);
470: } else if (ntype == NORM_INFINITY) {
471: PetscReal tmp;
472: for (j=0; j<bs; j++) tnorm[j] = 0.0;
474: for (i=0; i<n; i+=bs) {
475: for (j=0; j<bs; j++) {
476: if ((tmp = PetscAbsScalar(x[i+j])) > tnorm[j]) tnorm[j] = tmp;
477: /* check special case of tmp == NaN */
478: if (tmp != tmp) {tnorm[j] = tmp; break;}
479: }
480: }
481: MPI_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_MAX,comm);
482: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
483: VecRestoreArray(v,&x);
484: return(0);
485: }
489: /*@
490: VecStrideMaxAll - Computes the maximums of subvectors of a vector defined
491: by a starting point and a stride and optionally its location.
493: Collective on Vec
495: Input Parameter:
496: . v - the vector
498: Output Parameter:
499: + index - the location where the maximum occurred (not supported, pass NULL,
500: if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
501: - nrm - the maximums
503: Notes:
504: One must call VecSetBlockSize() before this routine to set the stride
505: information, or use a vector created from a multicomponent DMDA.
507: This is useful for computing, say the maximum of the pressure variable when
508: the pressure is stored (interlaced) with other variables, e.g., density, etc.
509: This will only work if the desire subvector is a stride subvector.
511: Level: advanced
513: Concepts: maximum^on stride of vector
514: Concepts: stride^maximum
516: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
517: @*/
518: PetscErrorCode VecStrideMaxAll(Vec v,PetscInt idex[],PetscReal nrm[])
519: {
521: PetscInt i,j,n,bs;
522: PetscScalar *x;
523: PetscReal max[128],tmp;
524: MPI_Comm comm;
529: 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");
530: VecGetLocalSize(v,&n);
531: VecGetArray(v,&x);
532: PetscObjectGetComm((PetscObject)v,&comm);
534: VecGetBlockSize(v,&bs);
535: if (bs > 128) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");
537: if (!n) {
538: for (j=0; j<bs; j++) max[j] = PETSC_MIN_REAL;
539: } else {
540: for (j=0; j<bs; j++) max[j] = PetscRealPart(x[j]);
542: for (i=bs; i<n; i+=bs) {
543: for (j=0; j<bs; j++) {
544: if ((tmp = PetscRealPart(x[i+j])) > max[j]) max[j] = tmp;
545: }
546: }
547: }
548: MPI_Allreduce(max,nrm,bs,MPIU_REAL,MPIU_MAX,comm);
550: VecRestoreArray(v,&x);
551: return(0);
552: }
556: /*@
557: VecStrideMinAll - Computes the minimum of subvector of a vector defined
558: by a starting point and a stride and optionally its location.
560: Collective on Vec
562: Input Parameter:
563: . v - the vector
565: Output Parameter:
566: + idex - the location where the minimum occurred (not supported, pass NULL,
567: if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
568: - nrm - the minimums
570: Level: advanced
572: Notes:
573: One must call VecSetBlockSize() before this routine to set the stride
574: information, or use a vector created from a multicomponent DMDA.
576: This is useful for computing, say the minimum of the pressure variable when
577: the pressure is stored (interlaced) with other variables, e.g., density, etc.
578: This will only work if the desire subvector is a stride subvector.
580: Concepts: minimum^on stride of vector
581: Concepts: stride^minimum
583: .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
584: @*/
585: PetscErrorCode VecStrideMinAll(Vec v,PetscInt idex[],PetscReal nrm[])
586: {
588: PetscInt i,n,bs,j;
589: PetscScalar *x;
590: PetscReal min[128],tmp;
591: MPI_Comm comm;
596: 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");
597: VecGetLocalSize(v,&n);
598: VecGetArray(v,&x);
599: PetscObjectGetComm((PetscObject)v,&comm);
601: VecGetBlockSize(v,&bs);
602: if (bs > 128) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");
604: if (!n) {
605: for (j=0; j<bs; j++) min[j] = PETSC_MAX_REAL;
606: } else {
607: for (j=0; j<bs; j++) min[j] = PetscRealPart(x[j]);
609: for (i=bs; i<n; i+=bs) {
610: for (j=0; j<bs; j++) {
611: if ((tmp = PetscRealPart(x[i+j])) < min[j]) min[j] = tmp;
612: }
613: }
614: }
615: MPI_Allreduce(min,nrm,bs,MPIU_REAL,MPIU_MIN,comm);
617: VecRestoreArray(v,&x);
618: return(0);
619: }
621: /*----------------------------------------------------------------------------------------------*/
624: /*@
625: VecStrideGatherAll - Gathers all the single components from a multi-component vector into
626: separate vectors.
628: Collective on Vec
630: Input Parameter:
631: + v - the vector
632: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
634: Output Parameter:
635: . s - the location where the subvectors are stored
637: Notes:
638: One must call VecSetBlockSize() before this routine to set the stride
639: information, or use a vector created from a multicomponent DMDA.
641: If x is the array representing the vector x then this gathers
642: the arrays (x[start],x[start+stride],x[start+2*stride], ....)
643: for start=0,1,2,...bs-1
645: The parallel layout of the vector and the subvector must be the same;
646: i.e., nlocal of v = stride*(nlocal of s)
648: Not optimized; could be easily
650: Level: advanced
652: Concepts: gather^into strided vector
654: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
655: VecStrideScatterAll()
656: @*/
657: PetscErrorCode VecStrideGatherAll(Vec v,Vec s[],InsertMode addv)
658: {
660: PetscInt i,n,n2,bs,j,k,*bss = NULL,nv,jj,nvc;
661: PetscScalar *x,**y;
667: VecGetLocalSize(v,&n);
668: VecGetLocalSize(s[0],&n2);
669: VecGetArray(v,&x);
670: VecGetBlockSize(v,&bs);
671: if (bs < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Input vector does not have a valid blocksize set");
673: PetscMalloc2(bs,&y,bs,&bss);
674: nv = 0;
675: nvc = 0;
676: for (i=0; i<bs; i++) {
677: VecGetBlockSize(s[i],&bss[i]);
678: if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1 Re: [PETSC #8241] VecStrideGatherAll */
679: VecGetArray(s[i],&y[i]);
680: nvc += bss[i];
681: nv++;
682: if (nvc > bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of subvectors in subvectors > number of vectors in main vector");
683: if (nvc == bs) break;
684: }
686: n = n/bs;
688: jj = 0;
689: if (addv == INSERT_VALUES) {
690: for (j=0; j<nv; j++) {
691: for (k=0; k<bss[j]; k++) {
692: for (i=0; i<n; i++) y[j][i*bss[j] + k] = x[bs*i+jj+k];
693: }
694: jj += bss[j];
695: }
696: } else if (addv == ADD_VALUES) {
697: for (j=0; j<nv; j++) {
698: for (k=0; k<bss[j]; k++) {
699: for (i=0; i<n; i++) y[j][i*bss[j] + k] += x[bs*i+jj+k];
700: }
701: jj += bss[j];
702: }
703: #if !defined(PETSC_USE_COMPLEX)
704: } else if (addv == MAX_VALUES) {
705: for (j=0; j<nv; j++) {
706: for (k=0; k<bss[j]; k++) {
707: for (i=0; i<n; i++) y[j][i*bss[j] + k] = PetscMax(y[j][i*bss[j] + k],x[bs*i+jj+k]);
708: }
709: jj += bss[j];
710: }
711: #endif
712: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
714: VecRestoreArray(v,&x);
715: for (i=0; i<nv; i++) {
716: VecRestoreArray(s[i],&y[i]);
717: }
719: PetscFree2(y,bss);
720: return(0);
721: }
725: /*@
726: VecStrideScatterAll - Scatters all the single components from separate vectors into
727: a multi-component vector.
729: Collective on Vec
731: Input Parameter:
732: + s - the location where the subvectors are stored
733: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
735: Output Parameter:
736: . v - the multicomponent vector
738: Notes:
739: One must call VecSetBlockSize() before this routine to set the stride
740: information, or use a vector created from a multicomponent DMDA.
742: The parallel layout of the vector and the subvector must be the same;
743: i.e., nlocal of v = stride*(nlocal of s)
745: Not optimized; could be easily
747: Level: advanced
749: Concepts: scatter^into strided vector
751: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
752: VecStrideScatterAll()
753: @*/
754: PetscErrorCode VecStrideScatterAll(Vec s[],Vec v,InsertMode addv)
755: {
757: PetscInt i,n,n2,bs,j,jj,k,*bss = NULL,nv,nvc;
758: PetscScalar *x,**y;
764: VecGetLocalSize(v,&n);
765: VecGetLocalSize(s[0],&n2);
766: VecGetArray(v,&x);
767: VecGetBlockSize(v,&bs);
768: if (bs < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Input vector does not have a valid blocksize set");
770: PetscMalloc2(bs,&y,bs,&bss);
771: nv = 0;
772: nvc = 0;
773: for (i=0; i<bs; i++) {
774: VecGetBlockSize(s[i],&bss[i]);
775: if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1 Re: [PETSC #8241] VecStrideGatherAll */
776: VecGetArray(s[i],&y[i]);
777: nvc += bss[i];
778: nv++;
779: if (nvc > bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of subvectors in subvectors > number of vectors in main vector");
780: if (nvc == bs) break;
781: }
783: n = n/bs;
785: jj = 0;
786: if (addv == INSERT_VALUES) {
787: for (j=0; j<nv; j++) {
788: for (k=0; k<bss[j]; k++) {
789: for (i=0; i<n; i++) x[bs*i+jj+k] = y[j][i*bss[j] + k];
790: }
791: jj += bss[j];
792: }
793: } else if (addv == ADD_VALUES) {
794: for (j=0; j<nv; j++) {
795: for (k=0; k<bss[j]; k++) {
796: for (i=0; i<n; i++) x[bs*i+jj+k] += y[j][i*bss[j] + k];
797: }
798: jj += bss[j];
799: }
800: #if !defined(PETSC_USE_COMPLEX)
801: } else if (addv == MAX_VALUES) {
802: for (j=0; j<nv; j++) {
803: for (k=0; k<bss[j]; k++) {
804: for (i=0; i<n; i++) x[bs*i+jj+k] = PetscMax(x[bs*i+jj+k],y[j][i*bss[j] + k]);
805: }
806: jj += bss[j];
807: }
808: #endif
809: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
811: VecRestoreArray(v,&x);
812: for (i=0; i<nv; i++) {
813: VecRestoreArray(s[i],&y[i]);
814: }
815: PetscFree2(y,bss);
816: return(0);
817: }
821: /*@
822: VecStrideGather - Gathers a single component from a multi-component vector into
823: another vector.
825: Collective on Vec
827: Input Parameter:
828: + v - the vector
829: . start - starting point of the subvector (defined by a stride)
830: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
832: Output Parameter:
833: . s - the location where the subvector is stored
835: Notes:
836: One must call VecSetBlockSize() before this routine to set the stride
837: information, or use a vector created from a multicomponent DMDA.
839: If x is the array representing the vector x then this gathers
840: the array (x[start],x[start+stride],x[start+2*stride], ....)
842: The parallel layout of the vector and the subvector must be the same;
843: i.e., nlocal of v = stride*(nlocal of s)
845: Not optimized; could be easily
847: Level: advanced
849: Concepts: gather^into strided vector
851: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
852: VecStrideScatterAll()
853: @*/
854: PetscErrorCode VecStrideGather(Vec v,PetscInt start,Vec s,InsertMode addv)
855: {
861: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
862: 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);
863: if (!v->ops->stridegather) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
864: (*v->ops->stridegather)(v,start,s,addv);
865: return(0);
866: }
870: /*@
871: VecStrideScatter - Scatters a single component from a vector into a multi-component vector.
873: Collective on Vec
875: Input Parameter:
876: + s - the single-component vector
877: . start - starting point of the subvector (defined by a stride)
878: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
880: Output Parameter:
881: . v - the location where the subvector is scattered (the multi-component vector)
883: Notes:
884: One must call VecSetBlockSize() on the multi-component vector before this
885: routine to set the stride information, or use a vector created from a multicomponent DMDA.
887: The parallel layout of the vector and the subvector must be the same;
888: i.e., nlocal of v = stride*(nlocal of s)
890: Not optimized; could be easily
892: Level: advanced
894: Concepts: scatter^into strided vector
896: .seealso: VecStrideNorm(), VecStrideGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
897: VecStrideScatterAll(), VecStrideSubSetScatter(), VecStrideSubSetGather()
898: @*/
899: PetscErrorCode VecStrideScatter(Vec s,PetscInt start,Vec v,InsertMode addv)
900: {
906: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
907: 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);
908: if (!v->ops->stridescatter) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
909: (*v->ops->stridescatter)(s,start,v,addv);
910: return(0);
911: }
915: /*@
916: VecStrideSubSetGather - Gathers a subset of components from a multi-component vector into
917: another vector.
919: Collective on Vec
921: Input Parameter:
922: + v - the vector
923: . nidx - the number of indices
924: . idxv - the indices of the components 0 <= idxv[0] ...idxv[nidx-1] < bs(v), they need not be sorted
925: . 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
926: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
928: Output Parameter:
929: . s - the location where the subvector is stored
931: Notes:
932: One must call VecSetBlockSize() on both vectors before this routine to set the stride
933: information, or use a vector created from a multicomponent DMDA.
936: The parallel layout of the vector and the subvector must be the same;
938: Not optimized; could be easily
940: Level: advanced
942: Concepts: gather^into strided vector
944: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideGather(), VecStrideSubSetScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
945: VecStrideScatterAll()
946: @*/
947: PetscErrorCode VecStrideSubSetGather(Vec v,PetscInt nidx,const PetscInt idxv[],const PetscInt idxs[],Vec s,InsertMode addv)
948: {
954: if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
955: if (!v->ops->stridesubsetgather) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
956: (*v->ops->stridesubsetgather)(v,nidx,idxv,idxs,s,addv);
957: return(0);
958: }
962: /*@
963: VecStrideSubSetScatter - Scatters components from a vector into a subset of components of a multi-component vector.
965: Collective on Vec
967: Input Parameter:
968: + s - the smaller-component vector
969: . nidx - the number of indices in idx
970: . 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
971: . idxv - the indices of the components in the larger-component vector, 0 <= idx[0] ...idx[nidx-1] < bs(v) they need not be sorted
972: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
974: Output Parameter:
975: . v - the location where the subvector is into scattered (the multi-component vector)
977: Notes:
978: One must call VecSetBlockSize() on the vectors before this
979: routine to set the stride information, or use a vector created from a multicomponent DMDA.
981: The parallel layout of the vector and the subvector must be the same;
983: Not optimized; could be easily
985: Level: advanced
987: Concepts: scatter^into strided vector
989: .seealso: VecStrideNorm(), VecStrideGather(), VecStrideGather(), VecStrideSubSetGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
990: VecStrideScatterAll()
991: @*/
992: PetscErrorCode VecStrideSubSetScatter(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)
993: {
999: if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
1000: if (!v->ops->stridesubsetscatter) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
1001: (*v->ops->stridesubsetscatter)(s,nidx,idxs,idxv,v,addv);
1002: return(0);
1003: }
1007: PetscErrorCode VecStrideGather_Default(Vec v,PetscInt start,Vec s,InsertMode addv)
1008: {
1010: PetscInt i,n,bs,ns;
1011: const PetscScalar *x;
1012: PetscScalar *y;
1015: VecGetLocalSize(v,&n);
1016: VecGetLocalSize(s,&ns);
1017: VecGetArrayRead(v,&x);
1018: VecGetArray(s,&y);
1020: bs = v->map->bs;
1021: 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);
1022: x += start;
1023: n = n/bs;
1025: if (addv == INSERT_VALUES) {
1026: for (i=0; i<n; i++) y[i] = x[bs*i];
1027: } else if (addv == ADD_VALUES) {
1028: for (i=0; i<n; i++) y[i] += x[bs*i];
1029: #if !defined(PETSC_USE_COMPLEX)
1030: } else if (addv == MAX_VALUES) {
1031: for (i=0; i<n; i++) y[i] = PetscMax(y[i],x[bs*i]);
1032: #endif
1033: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
1035: VecRestoreArrayRead(v,&x);
1036: VecRestoreArray(s,&y);
1037: return(0);
1038: }
1042: PetscErrorCode VecStrideScatter_Default(Vec s,PetscInt start,Vec v,InsertMode addv)
1043: {
1044: PetscErrorCode ierr;
1045: PetscInt i,n,bs,ns;
1046: PetscScalar *x;
1047: const PetscScalar *y;
1050: VecGetLocalSize(v,&n);
1051: VecGetLocalSize(s,&ns);
1052: VecGetArray(v,&x);
1053: VecGetArrayRead(s,&y);
1055: bs = v->map->bs;
1056: 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);
1057: x += start;
1058: n = n/bs;
1060: if (addv == INSERT_VALUES) {
1061: for (i=0; i<n; i++) x[bs*i] = y[i];
1062: } else if (addv == ADD_VALUES) {
1063: for (i=0; i<n; i++) x[bs*i] += y[i];
1064: #if !defined(PETSC_USE_COMPLEX)
1065: } else if (addv == MAX_VALUES) {
1066: for (i=0; i<n; i++) x[bs*i] = PetscMax(y[i],x[bs*i]);
1067: #endif
1068: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
1070: VecRestoreArray(v,&x);
1071: VecRestoreArrayRead(s,&y);
1072: return(0);
1073: }
1077: PetscErrorCode VecStrideSubSetGather_Default(Vec v,PetscInt nidx,const PetscInt idxv[],const PetscInt idxs[],Vec s,InsertMode addv)
1078: {
1079: PetscErrorCode ierr;
1080: PetscInt i,j,n,bs,bss,ns;
1081: const PetscScalar *x;
1082: PetscScalar *y;
1085: VecGetLocalSize(v,&n);
1086: VecGetLocalSize(s,&ns);
1087: VecGetArrayRead(v,&x);
1088: VecGetArray(s,&y);
1090: bs = v->map->bs;
1091: bss = s->map->bs;
1092: n = n/bs;
1094: #if defined(PETSC_DEBUG)
1095: if (n != ns/bss) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible layout of vectors");
1096: for (j=0; j<nidx; j++) {
1097: if (idxv[j] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is negative",j,idxv[j]);
1098: if (idxv[j] >= bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is greater than or equal to vector blocksize %D",j,idxv[j],bs);
1099: }
1100: if (!idxs && bss != nidx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide idxs when not gathering into all locations");
1101: #endif
1103: if (addv == INSERT_VALUES) {
1104: if (!idxs) {
1105: for (i=0; i<n; i++) {
1106: for (j=0; j<bss; j++) y[bss*i+j] = x[bs*i+idxv[j]];
1107: }
1108: } else {
1109: for (i=0; i<n; i++) {
1110: for (j=0; j<bss; j++) y[bss*i+idxs[j]] = x[bs*i+idxv[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++) y[bss*i+j] += x[bs*i+idxv[j]];
1117: }
1118: } else {
1119: for (i=0; i<n; i++) {
1120: for (j=0; j<bss; j++) y[bss*i+idxs[j]] += x[bs*i+idxv[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++) y[bss*i+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++) y[bss*i+idxs[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: VecRestoreArrayRead(v,&x);
1138: VecRestoreArray(s,&y);
1139: return(0);
1140: }
1144: PetscErrorCode VecStrideSubSetScatter_Default(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)
1145: {
1146: PetscErrorCode ierr;
1147: PetscInt j,i,n,bs,ns,bss;
1148: PetscScalar *x;
1149: const PetscScalar *y;
1152: VecGetLocalSize(v,&n);
1153: VecGetLocalSize(s,&ns);
1154: VecGetArray(v,&x);
1155: VecGetArrayRead(s,&y);
1157: bs = v->map->bs;
1158: bss = s->map->bs;
1159: n = n/bs;
1161: #if defined(PETSC_DEBUG)
1162: if (n != ns/bss) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible layout of vectors");
1163: for (j=0; j<bss; j++) {
1164: if (idx[j] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is negative",j,idx[j]);
1165: if (idx[j] >= bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is greater than or equal to vector blocksize %D",j,idx[j],bs);
1166: }
1167: if (!idxs && bss != nidx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide idxs when not scattering from all locations");
1168: #endif
1170: if (addv == INSERT_VALUES) {
1171: if (!idxs) {
1172: for (i=0; i<n; i++) {
1173: for (j=0; j<bss; j++) x[bs*i + idxv[j]] = y[bss*i+j];
1174: }
1175: } else {
1176: for (i=0; i<n; i++) {
1177: for (j=0; j<bss; j++) x[bs*i + idxv[j]] = y[bss*i+idxs[j]];
1178: }
1179: }
1180: } else if (addv == ADD_VALUES) {
1181: if (!idxs) {
1182: for (i=0; i<n; i++) {
1183: for (j=0; j<bss; j++) x[bs*i + idxv[j]] += y[bss*i+j];
1184: }
1185: } else {
1186: for (i=0; i<n; i++) {
1187: for (j=0; j<bss; j++) x[bs*i + idxv[j]] += y[bss*i+idxs[j]];
1188: }
1189: }
1190: #if !defined(PETSC_USE_COMPLEX)
1191: } else if (addv == MAX_VALUES) {
1192: if (!idxs) {
1193: for (i=0; i<n; i++) {
1194: for (j=0; j<bss; j++) x[bs*i + idxv[j]] = PetscMax(y[bss*i+j],x[bs*i + idxv[j]]);
1195: }
1196: } else {
1197: for (i=0; i<n; i++) {
1198: for (j=0; j<bss; j++) x[bs*i + idxv[j]] = PetscMax(y[bss*i+idxs[j]],x[bs*i + idxv[j]]);
1199: }
1200: }
1201: #endif
1202: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
1204: VecRestoreArray(v,&x);
1205: VecRestoreArrayRead(s,&y);
1206: return(0);
1207: }
1211: PetscErrorCode VecReciprocal_Default(Vec v)
1212: {
1214: PetscInt i,n;
1215: PetscScalar *x;
1218: VecGetLocalSize(v,&n);
1219: VecGetArray(v,&x);
1220: for (i=0; i<n; i++) {
1221: if (x[i] != (PetscScalar)0.0) x[i] = (PetscScalar)1.0/x[i];
1222: }
1223: VecRestoreArray(v,&x);
1224: return(0);
1225: }
1229: /*@
1230: VecExp - Replaces each component of a vector by e^x_i
1232: Not collective
1234: Input Parameter:
1235: . v - The vector
1237: Output Parameter:
1238: . v - The vector of exponents
1240: Level: beginner
1242: .seealso: VecLog(), VecAbs(), VecSqrtAbs(), VecReciprocal()
1244: .keywords: vector, sqrt, square root
1245: @*/
1246: PetscErrorCode VecExp(Vec v)
1247: {
1248: PetscScalar *x;
1249: PetscInt i, n;
1254: if (v->ops->exp) {
1255: (*v->ops->exp)(v);
1256: } else {
1257: VecGetLocalSize(v, &n);
1258: VecGetArray(v, &x);
1259: for (i = 0; i < n; i++) x[i] = PetscExpScalar(x[i]);
1260: VecRestoreArray(v, &x);
1261: }
1262: return(0);
1263: }
1267: /*@
1268: VecLog - Replaces each component of a vector by log(x_i), the natural logarithm
1270: Not collective
1272: Input Parameter:
1273: . v - The vector
1275: Output Parameter:
1276: . v - The vector of logs
1278: Level: beginner
1280: .seealso: VecExp(), VecAbs(), VecSqrtAbs(), VecReciprocal()
1282: .keywords: vector, sqrt, square root
1283: @*/
1284: PetscErrorCode VecLog(Vec v)
1285: {
1286: PetscScalar *x;
1287: PetscInt i, n;
1292: if (v->ops->log) {
1293: (*v->ops->log)(v);
1294: } else {
1295: VecGetLocalSize(v, &n);
1296: VecGetArray(v, &x);
1297: for (i = 0; i < n; i++) x[i] = PetscLogScalar(x[i]);
1298: VecRestoreArray(v, &x);
1299: }
1300: return(0);
1301: }
1305: /*@
1306: VecSqrtAbs - Replaces each component of a vector by the square root of its magnitude.
1308: Not collective
1310: Input Parameter:
1311: . v - The vector
1313: Output Parameter:
1314: . v - The vector square root
1316: Level: beginner
1318: Note: The actual function is sqrt(|x_i|)
1320: .seealso: VecLog(), VecExp(), VecReciprocal(), VecAbs()
1322: .keywords: vector, sqrt, square root
1323: @*/
1324: PetscErrorCode VecSqrtAbs(Vec v)
1325: {
1326: PetscScalar *x;
1327: PetscInt i, n;
1332: if (v->ops->sqrt) {
1333: (*v->ops->sqrt)(v);
1334: } else {
1335: VecGetLocalSize(v, &n);
1336: VecGetArray(v, &x);
1337: for (i = 0; i < n; i++) x[i] = PetscSqrtReal(PetscAbsScalar(x[i]));
1338: VecRestoreArray(v, &x);
1339: }
1340: return(0);
1341: }
1345: /*@
1346: VecDotNorm2 - computes the inner product of two vectors and the 2-norm squared of the second vector
1348: Collective on Vec
1350: Input Parameter:
1351: + s - first vector
1352: - t - second vector
1354: Output Parameter:
1355: + dp - s'conj(t)
1356: - nm - t'conj(t)
1358: Level: advanced
1360: Notes: conj(x) is the complex conjugate of x when x is complex
1363: .seealso: VecDot(), VecNorm(), VecDotBegin(), VecNormBegin(), VecDotEnd(), VecNormEnd()
1365: .keywords: vector, sqrt, square root
1366: @*/
1367: PetscErrorCode VecDotNorm2(Vec s,Vec t,PetscScalar *dp, PetscReal *nm)
1368: {
1369: PetscScalar *sx, *tx, dpx = 0.0, nmx = 0.0,work[2],sum[2];
1370: PetscInt i, n;
1381: if (s->map->N != t->map->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
1382: if (s->map->n != t->map->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");
1384: PetscLogEventBarrierBegin(VEC_DotNormBarrier,s,t,0,0,PetscObjectComm((PetscObject)s));
1385: if (s->ops->dotnorm2) {
1386: (*s->ops->dotnorm2)(s,t,dp,&dpx);
1387: *nm = PetscRealPart(dpx);
1388: } else {
1389: VecGetLocalSize(s, &n);
1390: VecGetArray(s, &sx);
1391: VecGetArray(t, &tx);
1393: for (i = 0; i<n; i++) {
1394: dpx += sx[i]*PetscConj(tx[i]);
1395: nmx += tx[i]*PetscConj(tx[i]);
1396: }
1397: work[0] = dpx;
1398: work[1] = nmx;
1400: MPI_Allreduce(work,sum,2,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)s));
1401: *dp = sum[0];
1402: *nm = PetscRealPart(sum[1]);
1404: VecRestoreArray(t, &tx);
1405: VecRestoreArray(s, &sx);
1406: PetscLogFlops(4.0*n);
1407: }
1408: PetscLogEventBarrierEnd(VEC_DotNormBarrier,s,t,0,0,PetscObjectComm((PetscObject)s));
1409: return(0);
1410: }
1414: /*@
1415: VecSum - Computes the sum of all the components of a vector.
1417: Collective on Vec
1419: Input Parameter:
1420: . v - the vector
1422: Output Parameter:
1423: . sum - the result
1425: Level: beginner
1427: Concepts: sum^of vector entries
1429: .seealso: VecNorm()
1430: @*/
1431: PetscErrorCode VecSum(Vec v,PetscScalar *sum)
1432: {
1434: PetscInt i,n;
1435: PetscScalar *x,lsum = 0.0;
1440: VecGetLocalSize(v,&n);
1441: VecGetArray(v,&x);
1442: for (i=0; i<n; i++) lsum += x[i];
1443: MPI_Allreduce(&lsum,sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)v));
1444: VecRestoreArray(v,&x);
1445: return(0);
1446: }
1450: /*@
1451: VecShift - Shifts all of the components of a vector by computing
1452: x[i] = x[i] + shift.
1454: Logically Collective on Vec
1456: Input Parameters:
1457: + v - the vector
1458: - shift - the shift
1460: Output Parameter:
1461: . v - the shifted vector
1463: Level: intermediate
1465: Concepts: vector^adding constant
1467: @*/
1468: PetscErrorCode VecShift(Vec v,PetscScalar shift)
1469: {
1471: PetscInt i,n;
1472: PetscScalar *x;
1477: if (v->ops->shift) {
1478: (*v->ops->shift)(v);
1479: } else {
1480: VecGetLocalSize(v,&n);
1481: VecGetArray(v,&x);
1482: for (i=0; i<n; i++) x[i] += shift;
1483: VecRestoreArray(v,&x);
1484: }
1485: return(0);
1486: }
1490: /*@
1491: VecAbs - Replaces every element in a vector with its absolute value.
1493: Logically Collective on Vec
1495: Input Parameters:
1496: . v - the vector
1498: Level: intermediate
1500: Concepts: vector^absolute value
1502: @*/
1503: PetscErrorCode VecAbs(Vec v)
1504: {
1506: PetscInt i,n;
1507: PetscScalar *x;
1511: if (v->ops->abs) {
1512: (*v->ops->abs)(v);
1513: } else {
1514: VecGetLocalSize(v,&n);
1515: VecGetArray(v,&x);
1516: for (i=0; i<n; i++) x[i] = PetscAbsScalar(x[i]);
1517: VecRestoreArray(v,&x);
1518: }
1519: return(0);
1520: }
1524: /*@
1525: VecPermute - Permutes a vector in place using the given ordering.
1527: Input Parameters:
1528: + vec - The vector
1529: . order - The ordering
1530: - inv - The flag for inverting the permutation
1532: Level: beginner
1534: Note: This function does not yet support parallel Index Sets with non-local permutations
1536: .seealso: MatPermute()
1537: .keywords: vec, permute
1538: @*/
1539: PetscErrorCode VecPermute(Vec x, IS row, PetscBool inv)
1540: {
1541: PetscScalar *array, *newArray;
1542: const PetscInt *idx;
1543: PetscInt i,rstart,rend;
1547: VecGetOwnershipRange(x,&rstart,&rend);
1548: ISGetIndices(row, &idx);
1549: VecGetArray(x, &array);
1550: PetscMalloc1(x->map->n, &newArray);
1551: #if defined(PETSC_USE_DEBUG)
1552: for (i = 0; i < x->map->n; i++) {
1553: 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]);
1554: }
1555: #endif
1556: if (!inv) {
1557: for (i = 0; i < x->map->n; i++) newArray[i] = array[idx[i]-rstart];
1558: } else {
1559: for (i = 0; i < x->map->n; i++) newArray[idx[i]-rstart] = array[i];
1560: }
1561: VecRestoreArray(x, &array);
1562: ISRestoreIndices(row, &idx);
1563: VecReplaceArray(x, newArray);
1564: return(0);
1565: }
1569: /*@
1570: VecEqual - Compares two vectors. Returns true if the two vectors are either pointing to the same memory buffer,
1571: or if the two vectors have the same local and global layout as well as bitwise equality of all entries.
1572: Does NOT take round-off errors into account.
1574: Collective on Vec
1576: Input Parameters:
1577: + vec1 - the first vector
1578: - vec2 - the second vector
1580: Output Parameter:
1581: . flg - PETSC_TRUE if the vectors are equal; PETSC_FALSE otherwise.
1583: Level: intermediate
1585: Concepts: equal^two vectors
1586: Concepts: vector^equality
1588: @*/
1589: PetscErrorCode VecEqual(Vec vec1,Vec vec2,PetscBool *flg)
1590: {
1591: const PetscScalar *v1,*v2;
1592: PetscErrorCode ierr;
1593: PetscInt n1,n2,N1,N2;
1594: PetscBool flg1;
1600: if (vec1 == vec2) *flg = PETSC_TRUE;
1601: else {
1602: VecGetSize(vec1,&N1);
1603: VecGetSize(vec2,&N2);
1604: if (N1 != N2) flg1 = PETSC_FALSE;
1605: else {
1606: VecGetLocalSize(vec1,&n1);
1607: VecGetLocalSize(vec2,&n2);
1608: if (n1 != n2) flg1 = PETSC_FALSE;
1609: else {
1610: VecGetArrayRead(vec1,&v1);
1611: VecGetArrayRead(vec2,&v2);
1612: PetscMemcmp(v1,v2,n1*sizeof(PetscScalar),&flg1);
1613: VecRestoreArrayRead(vec1,&v1);
1614: VecRestoreArrayRead(vec2,&v2);
1615: }
1616: }
1617: /* combine results from all processors */
1618: MPI_Allreduce(&flg1,flg,1,MPIU_BOOL,MPI_MIN,PetscObjectComm((PetscObject)vec1));
1619: }
1620: return(0);
1621: }
1625: /*@
1626: VecUniqueEntries - Compute the number of unique entries, and those entries
1628: Collective on Vec
1630: Input Parameter:
1631: . vec - the vector
1633: Output Parameters:
1634: + n - The number of unique entries
1635: - e - The entries
1637: Level: intermediate
1639: @*/
1640: PetscErrorCode VecUniqueEntries(Vec vec, PetscInt *n, PetscScalar **e)
1641: {
1642: PetscScalar *v, *tmp, *vals;
1643: PetscMPIInt *N, *displs, l;
1644: PetscInt ng, m, i, j, p;
1645: PetscMPIInt size;
1651: MPI_Comm_size(PetscObjectComm((PetscObject) vec), &size);
1652: VecGetLocalSize(vec, &m);
1653: VecGetArray(vec, &v);
1654: PetscMalloc2(m,&tmp,size,&N);
1655: for (i = 0, j = 0, l = 0; i < m; ++i) {
1656: /* Can speed this up with sorting */
1657: for (j = 0; j < l; ++j) {
1658: if (v[i] == tmp[j]) break;
1659: }
1660: if (j == l) {
1661: tmp[j] = v[i];
1662: ++l;
1663: }
1664: }
1665: /* Gather serial results */
1666: MPI_Allgather(&l, 1, MPI_INT, N, 1, MPI_INT, PetscObjectComm((PetscObject) vec));
1667: for (p = 0, ng = 0; p < size; ++p) {
1668: ng += N[p];
1669: }
1670: PetscMalloc2(ng,&vals,size+1,&displs);
1671: for (p = 1, displs[0] = 0; p <= size; ++p) {
1672: displs[p] = displs[p-1] + N[p-1];
1673: }
1674: MPI_Allgatherv(tmp, l, MPIU_SCALAR, vals, N, displs, MPIU_SCALAR, PetscObjectComm((PetscObject) vec));
1675: /* Find unique entries */
1676: #ifdef PETSC_USE_COMPLEX
1677: SETERRQ(PetscObjectComm((PetscObject) vec), PETSC_ERR_SUP, "Does not work with complex numbers");
1678: #else
1679: *n = displs[size];
1680: PetscSortRemoveDupsReal(n, (PetscReal *) vals);
1681: #endif
1682: if (e) {
1684: PetscMalloc1((*n), e);
1685: for (i = 0; i < *n; ++i) {
1686: (*e)[i] = vals[i];
1687: }
1688: }
1689: PetscFree2(vals,displs);
1690: PetscFree2(tmp,N);
1691: return(0);
1692: }