Actual source code: vinv.c
petsc-3.7.5 2017-01-01
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 set for each entry in that subvector
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;
45: VecLocked(v,1);
46:
47: VecGetLocalSize(v,&n);
48: VecGetArray(v,&x);
49: VecGetBlockSize(v,&bs);
50: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
51: 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);
52: x += start;
54: for (i=0; i<n; i+=bs) x[i] = s;
55: x -= start;
57: VecRestoreArray(v,&x);
58: return(0);
59: }
63: /*@
64: VecStrideScale - Scales a subvector of a vector defined
65: by a starting point and a stride.
67: Logically Collective on Vec
69: Input Parameter:
70: + v - the vector
71: . start - starting point of the subvector (defined by a stride)
72: - scale - value to multiply each subvector entry by
74: Notes:
75: One must call VecSetBlockSize() before this routine to set the stride
76: information, or use a vector created from a multicomponent DMDA.
78: This will only work if the desire subvector is a stride subvector
80: Level: advanced
82: Concepts: scale^on stride of vector
83: Concepts: stride^scale
85: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideScale()
86: @*/
87: PetscErrorCode VecStrideScale(Vec v,PetscInt start,PetscScalar scale)
88: {
90: PetscInt i,n,bs;
91: PetscScalar *x;
97: VecLocked(v,1);
98:
99: VecGetLocalSize(v,&n);
100: VecGetArray(v,&x);
101: VecGetBlockSize(v,&bs);
102: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
103: 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);
104: x += start;
106: for (i=0; i<n; i+=bs) x[i] *= scale;
107: x -= start;
109: VecRestoreArray(v,&x);
110: return(0);
111: }
115: /*@
116: VecStrideNorm - Computes the norm of subvector of a vector defined
117: by a starting point and a stride.
119: Collective on Vec
121: Input Parameter:
122: + v - the vector
123: . start - starting point of the subvector (defined by a stride)
124: - ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY
126: Output Parameter:
127: . norm - the norm
129: Notes:
130: One must call VecSetBlockSize() before this routine to set the stride
131: information, or use a vector created from a multicomponent DMDA.
133: If x is the array representing the vector x then this computes the norm
134: of the array (x[start],x[start+stride],x[start+2*stride], ....)
136: This is useful for computing, say the norm of the pressure variable when
137: the pressure is stored (interlaced) with other variables, say density etc.
139: This will only work if the desire subvector is a stride subvector
141: Level: advanced
143: Concepts: norm^on stride of vector
144: Concepts: stride^norm
146: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
147: @*/
148: PetscErrorCode VecStrideNorm(Vec v,PetscInt start,NormType ntype,PetscReal *nrm)
149: {
150: PetscErrorCode ierr;
151: PetscInt i,n,bs;
152: const PetscScalar *x;
153: PetscReal tnorm;
154: MPI_Comm comm;
159: VecGetLocalSize(v,&n);
160: VecGetArrayRead(v,&x);
161: PetscObjectGetComm((PetscObject)v,&comm);
163: VecGetBlockSize(v,&bs);
164: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
165: 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);
166: x += start;
168: if (ntype == NORM_2) {
169: PetscScalar sum = 0.0;
170: for (i=0; i<n; i+=bs) sum += x[i]*(PetscConj(x[i]));
171: tnorm = PetscRealPart(sum);
172: MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_SUM,comm);
173: *nrm = PetscSqrtReal(*nrm);
174: } else if (ntype == NORM_1) {
175: tnorm = 0.0;
176: for (i=0; i<n; i+=bs) tnorm += PetscAbsScalar(x[i]);
177: MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_SUM,comm);
178: } else if (ntype == NORM_INFINITY) {
179: PetscReal tmp;
180: tnorm = 0.0;
182: for (i=0; i<n; i+=bs) {
183: if ((tmp = PetscAbsScalar(x[i])) > tnorm) tnorm = tmp;
184: /* check special case of tmp == NaN */
185: if (tmp != tmp) {tnorm = tmp; break;}
186: }
187: MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_MAX,comm);
188: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
189: VecRestoreArrayRead(v,&x);
190: return(0);
191: }
195: /*@
196: VecStrideMax - Computes the maximum of subvector of a vector defined
197: by a starting point and a stride and optionally its location.
199: Collective on Vec
201: Input Parameter:
202: + v - the vector
203: - start - starting point of the subvector (defined by a stride)
205: Output Parameter:
206: + index - the location where the maximum occurred (pass NULL if not required)
207: - nrm - the maximum value in the subvector
209: Notes:
210: One must call VecSetBlockSize() before this routine to set the stride
211: information, or use a vector created from a multicomponent DMDA.
213: If xa is the array representing the vector x, then this computes the max
214: of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)
216: This is useful for computing, say the maximum of the pressure variable when
217: the pressure is stored (interlaced) with other variables, e.g., density, etc.
218: This will only work if the desire subvector is a stride subvector.
220: Level: advanced
222: Concepts: maximum^on stride of vector
223: Concepts: stride^maximum
225: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
226: @*/
227: PetscErrorCode VecStrideMax(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
228: {
229: PetscErrorCode ierr;
230: PetscInt i,n,bs,id;
231: const PetscScalar *x;
232: PetscReal max,tmp;
233: MPI_Comm comm;
239: VecGetLocalSize(v,&n);
240: VecGetArrayRead(v,&x);
241: PetscObjectGetComm((PetscObject)v,&comm);
243: VecGetBlockSize(v,&bs);
244: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
245: 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);
246: x += start;
248: id = -1;
249: if (!n) max = PETSC_MIN_REAL;
250: else {
251: id = 0;
252: max = PetscRealPart(x[0]);
253: for (i=bs; i<n; i+=bs) {
254: if ((tmp = PetscRealPart(x[i])) > max) { max = tmp; id = i;}
255: }
256: }
257: VecRestoreArrayRead(v,&x);
259: if (!idex) {
260: MPIU_Allreduce(&max,nrm,1,MPIU_REAL,MPIU_MAX,comm);
261: } else {
262: PetscReal in[2],out[2];
263: PetscInt rstart;
265: VecGetOwnershipRange(v,&rstart,NULL);
266: in[0] = max;
267: in[1] = rstart+id+start;
268: MPIU_Allreduce(in,out,2,MPIU_REAL,VecMax_Local_Op,PetscObjectComm((PetscObject)v));
269: *nrm = out[0];
270: *idex = (PetscInt)out[1];
271: }
272: return(0);
273: }
277: /*@
278: VecStrideMin - Computes the minimum of subvector of a vector defined
279: by a starting point and a stride and optionally its location.
281: Collective on Vec
283: Input Parameter:
284: + v - the vector
285: - start - starting point of the subvector (defined by a stride)
287: Output Parameter:
288: + idex - the location where the minimum occurred. (pass NULL if not required)
289: - nrm - the minimum value in the subvector
291: Level: advanced
293: Notes:
294: One must call VecSetBlockSize() before this routine to set the stride
295: information, or use a vector created from a multicomponent DMDA.
297: If xa is the array representing the vector x, then this computes the min
298: of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)
300: This is useful for computing, say the minimum of the pressure variable when
301: the pressure is stored (interlaced) with other variables, e.g., density, etc.
302: This will only work if the desire subvector is a stride subvector.
304: Concepts: minimum^on stride of vector
305: Concepts: stride^minimum
307: .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
308: @*/
309: PetscErrorCode VecStrideMin(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
310: {
311: PetscErrorCode ierr;
312: PetscInt i,n,bs,id;
313: const PetscScalar *x;
314: PetscReal min,tmp;
315: MPI_Comm comm;
321: VecGetLocalSize(v,&n);
322: VecGetArrayRead(v,&x);
323: PetscObjectGetComm((PetscObject)v,&comm);
325: VecGetBlockSize(v,&bs);
326: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
327: 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);
328: x += start;
330: id = -1;
331: if (!n) min = PETSC_MAX_REAL;
332: else {
333: id = 0;
334: min = PetscRealPart(x[0]);
335: for (i=bs; i<n; i+=bs) {
336: if ((tmp = PetscRealPart(x[i])) < min) { min = tmp; id = i;}
337: }
338: }
339: VecRestoreArrayRead(v,&x);
341: if (!idex) {
342: MPIU_Allreduce(&min,nrm,1,MPIU_REAL,MPIU_MIN,comm);
343: } else {
344: PetscReal in[2],out[2];
345: PetscInt rstart;
347: VecGetOwnershipRange(v,&rstart,NULL);
348: in[0] = min;
349: in[1] = rstart+id;
350: MPIU_Allreduce(in,out,2,MPIU_REAL,VecMin_Local_Op,PetscObjectComm((PetscObject)v));
351: *nrm = out[0];
352: *idex = (PetscInt)out[1];
353: }
354: return(0);
355: }
359: /*@
360: VecStrideScaleAll - Scales the subvectors of a vector defined
361: by a starting point and a stride.
363: Logically Collective on Vec
365: Input Parameter:
366: + v - the vector
367: - scales - values to multiply each subvector entry by
369: Notes:
370: One must call VecSetBlockSize() before this routine to set the stride
371: information, or use a vector created from a multicomponent DMDA.
373: The dimension of scales must be the same as the vector block size
376: Level: advanced
378: Concepts: scale^on stride of vector
379: Concepts: stride^scale
381: .seealso: VecNorm(), VecStrideScale(), VecScale(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
382: @*/
383: PetscErrorCode VecStrideScaleAll(Vec v,const PetscScalar *scales)
384: {
386: PetscInt i,j,n,bs;
387: PetscScalar *x;
392: VecLocked(v,1);
393:
394: VecGetLocalSize(v,&n);
395: VecGetArray(v,&x);
396: VecGetBlockSize(v,&bs);
398: /* need to provide optimized code for each bs */
399: for (i=0; i<n; i+=bs) {
400: for (j=0; j<bs; j++) x[i+j] *= scales[j];
401: }
402: VecRestoreArray(v,&x);
403: return(0);
404: }
408: /*@
409: VecStrideNormAll - Computes the norms of subvectors of a vector defined
410: by a starting point and a stride.
412: Collective on Vec
414: Input Parameter:
415: + v - the vector
416: - ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY
418: Output Parameter:
419: . nrm - the norms
421: Notes:
422: One must call VecSetBlockSize() before this routine to set the stride
423: information, or use a vector created from a multicomponent DMDA.
425: If x is the array representing the vector x then this computes the norm
426: of the array (x[start],x[start+stride],x[start+2*stride], ....) for each start < stride
428: The dimension of nrm must be the same as the vector block size
430: This will only work if the desire subvector is a stride subvector
432: Level: advanced
434: Concepts: norm^on stride of vector
435: Concepts: stride^norm
437: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
438: @*/
439: PetscErrorCode VecStrideNormAll(Vec v,NormType ntype,PetscReal nrm[])
440: {
441: PetscErrorCode ierr;
442: PetscInt i,j,n,bs;
443: const PetscScalar *x;
444: PetscReal tnorm[128];
445: MPI_Comm comm;
450: VecGetLocalSize(v,&n);
451: VecGetArrayRead(v,&x);
452: PetscObjectGetComm((PetscObject)v,&comm);
454: VecGetBlockSize(v,&bs);
455: if (bs > 128) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");
457: if (ntype == NORM_2) {
458: PetscScalar sum[128];
459: for (j=0; j<bs; j++) sum[j] = 0.0;
460: for (i=0; i<n; i+=bs) {
461: for (j=0; j<bs; j++) sum[j] += x[i+j]*(PetscConj(x[i+j]));
462: }
463: for (j=0; j<bs; j++) tnorm[j] = PetscRealPart(sum[j]);
465: MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_SUM,comm);
466: for (j=0; j<bs; j++) nrm[j] = PetscSqrtReal(nrm[j]);
467: } else if (ntype == NORM_1) {
468: for (j=0; j<bs; j++) tnorm[j] = 0.0;
470: for (i=0; i<n; i+=bs) {
471: for (j=0; j<bs; j++) tnorm[j] += PetscAbsScalar(x[i+j]);
472: }
474: MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_SUM,comm);
475: } else if (ntype == NORM_INFINITY) {
476: PetscReal tmp;
477: for (j=0; j<bs; j++) tnorm[j] = 0.0;
479: for (i=0; i<n; i+=bs) {
480: for (j=0; j<bs; j++) {
481: if ((tmp = PetscAbsScalar(x[i+j])) > tnorm[j]) tnorm[j] = tmp;
482: /* check special case of tmp == NaN */
483: if (tmp != tmp) {tnorm[j] = tmp; break;}
484: }
485: }
486: MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_MAX,comm);
487: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
488: VecRestoreArrayRead(v,&x);
489: return(0);
490: }
494: /*@
495: VecStrideMaxAll - Computes the maximums of subvectors of a vector defined
496: by a starting point and a stride and optionally its location.
498: Collective on Vec
500: Input Parameter:
501: . v - the vector
503: Output Parameter:
504: + index - the location where the maximum occurred (not supported, pass NULL,
505: if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
506: - nrm - the maximum values of each subvector
508: Notes:
509: One must call VecSetBlockSize() before this routine to set the stride
510: information, or use a vector created from a multicomponent DMDA.
512: The dimension of nrm must be the same as the vector block size
514: Level: advanced
516: Concepts: maximum^on stride of vector
517: Concepts: stride^maximum
519: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
520: @*/
521: PetscErrorCode VecStrideMaxAll(Vec v,PetscInt idex[],PetscReal nrm[])
522: {
523: PetscErrorCode ierr;
524: PetscInt i,j,n,bs;
525: const PetscScalar *x;
526: PetscReal max[128],tmp;
527: MPI_Comm comm;
532: 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");
533: VecGetLocalSize(v,&n);
534: VecGetArrayRead(v,&x);
535: PetscObjectGetComm((PetscObject)v,&comm);
537: VecGetBlockSize(v,&bs);
538: if (bs > 128) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");
540: if (!n) {
541: for (j=0; j<bs; j++) max[j] = PETSC_MIN_REAL;
542: } else {
543: for (j=0; j<bs; j++) max[j] = PetscRealPart(x[j]);
545: for (i=bs; i<n; i+=bs) {
546: for (j=0; j<bs; j++) {
547: if ((tmp = PetscRealPart(x[i+j])) > max[j]) max[j] = tmp;
548: }
549: }
550: }
551: MPIU_Allreduce(max,nrm,bs,MPIU_REAL,MPIU_MAX,comm);
553: VecRestoreArrayRead(v,&x);
554: return(0);
555: }
559: /*@
560: VecStrideMinAll - Computes the minimum of subvector of a vector defined
561: by a starting point and a stride and optionally its location.
563: Collective on Vec
565: Input Parameter:
566: . v - the vector
568: Output Parameter:
569: + idex - the location where the minimum occurred (not supported, pass NULL,
570: if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
571: - nrm - the minimums of each subvector
573: Level: advanced
575: Notes:
576: One must call VecSetBlockSize() before this routine to set the stride
577: information, or use a vector created from a multicomponent DMDA.
579: The dimension of nrm must be the same as the vector block size
581: Concepts: minimum^on stride of vector
582: Concepts: stride^minimum
584: .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
585: @*/
586: PetscErrorCode VecStrideMinAll(Vec v,PetscInt idex[],PetscReal nrm[])
587: {
588: PetscErrorCode ierr;
589: PetscInt i,n,bs,j;
590: const PetscScalar *x;
591: PetscReal min[128],tmp;
592: MPI_Comm comm;
597: 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");
598: VecGetLocalSize(v,&n);
599: VecGetArrayRead(v,&x);
600: PetscObjectGetComm((PetscObject)v,&comm);
602: VecGetBlockSize(v,&bs);
603: if (bs > 128) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");
605: if (!n) {
606: for (j=0; j<bs; j++) min[j] = PETSC_MAX_REAL;
607: } else {
608: for (j=0; j<bs; j++) min[j] = PetscRealPart(x[j]);
610: for (i=bs; i<n; i+=bs) {
611: for (j=0; j<bs; j++) {
612: if ((tmp = PetscRealPart(x[i+j])) < min[j]) min[j] = tmp;
613: }
614: }
615: }
616: MPIU_Allreduce(min,nrm,bs,MPIU_REAL,MPIU_MIN,comm);
618: VecRestoreArrayRead(v,&x);
619: return(0);
620: }
622: /*----------------------------------------------------------------------------------------------*/
625: /*@
626: VecStrideGatherAll - Gathers all the single components from a multi-component vector into
627: separate vectors.
629: Collective on Vec
631: Input Parameter:
632: + v - the vector
633: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
635: Output Parameter:
636: . s - the location where the subvectors are stored
638: Notes:
639: One must call VecSetBlockSize() before this routine to set the stride
640: information, or use a vector created from a multicomponent DMDA.
642: If x is the array representing the vector x then this gathers
643: the arrays (x[start],x[start+stride],x[start+2*stride], ....)
644: for start=0,1,2,...bs-1
646: The parallel layout of the vector and the subvector must be the same;
647: i.e., nlocal of v = stride*(nlocal of s)
649: Not optimized; could be easily
651: Level: advanced
653: Concepts: gather^into strided vector
655: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
656: VecStrideScatterAll()
657: @*/
658: PetscErrorCode VecStrideGatherAll(Vec v,Vec s[],InsertMode addv)
659: {
660: PetscErrorCode ierr;
661: PetscInt i,n,n2,bs,j,k,*bss = NULL,nv,jj,nvc;
662: PetscScalar **y;
663: const PetscScalar *x;
669: VecGetLocalSize(v,&n);
670: VecGetLocalSize(s[0],&n2);
671: VecGetArrayRead(v,&x);
672: VecGetBlockSize(v,&bs);
673: if (bs < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Input vector does not have a valid blocksize set");
675: PetscMalloc2(bs,&y,bs,&bss);
676: nv = 0;
677: nvc = 0;
678: for (i=0; i<bs; i++) {
679: VecGetBlockSize(s[i],&bss[i]);
680: if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1 Re: [PETSC #8241] VecStrideGatherAll */
681: VecGetArray(s[i],&y[i]);
682: nvc += bss[i];
683: nv++;
684: if (nvc > bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of subvectors in subvectors > number of vectors in main vector");
685: if (nvc == bs) break;
686: }
688: n = n/bs;
690: jj = 0;
691: if (addv == INSERT_VALUES) {
692: for (j=0; j<nv; j++) {
693: for (k=0; k<bss[j]; k++) {
694: for (i=0; i<n; i++) y[j][i*bss[j] + k] = x[bs*i+jj+k];
695: }
696: jj += bss[j];
697: }
698: } else if (addv == ADD_VALUES) {
699: for (j=0; j<nv; j++) {
700: for (k=0; k<bss[j]; k++) {
701: for (i=0; i<n; i++) y[j][i*bss[j] + k] += x[bs*i+jj+k];
702: }
703: jj += bss[j];
704: }
705: #if !defined(PETSC_USE_COMPLEX)
706: } else if (addv == MAX_VALUES) {
707: for (j=0; j<nv; j++) {
708: for (k=0; k<bss[j]; k++) {
709: for (i=0; i<n; i++) y[j][i*bss[j] + k] = PetscMax(y[j][i*bss[j] + k],x[bs*i+jj+k]);
710: }
711: jj += bss[j];
712: }
713: #endif
714: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
716: VecRestoreArrayRead(v,&x);
717: for (i=0; i<nv; i++) {
718: VecRestoreArray(s[i],&y[i]);
719: }
721: PetscFree2(y,bss);
722: return(0);
723: }
727: /*@
728: VecStrideScatterAll - Scatters all the single components from separate vectors into
729: a multi-component vector.
731: Collective on Vec
733: Input Parameter:
734: + s - the location where the subvectors are stored
735: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
737: Output Parameter:
738: . v - the multicomponent vector
740: Notes:
741: One must call VecSetBlockSize() before this routine to set the stride
742: information, or use a vector created from a multicomponent DMDA.
744: The parallel layout of the vector and the subvector must be the same;
745: i.e., nlocal of v = stride*(nlocal of s)
747: Not optimized; could be easily
749: Level: advanced
751: Concepts: scatter^into strided vector
753: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
754: VecStrideScatterAll()
755: @*/
756: PetscErrorCode VecStrideScatterAll(Vec s[],Vec v,InsertMode addv)
757: {
758: PetscErrorCode ierr;
759: PetscInt i,n,n2,bs,j,jj,k,*bss = NULL,nv,nvc;
760: PetscScalar *x,**y;
766: VecGetLocalSize(v,&n);
767: VecGetLocalSize(s[0],&n2);
768: VecGetArray(v,&x);
769: VecGetBlockSize(v,&bs);
770: if (bs < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Input vector does not have a valid blocksize set");
772: PetscMalloc2(bs,&y,bs,&bss);
773: nv = 0;
774: nvc = 0;
775: for (i=0; i<bs; i++) {
776: VecGetBlockSize(s[i],&bss[i]);
777: if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1 Re: [PETSC #8241] VecStrideGatherAll */
778: VecGetArray(s[i],&y[i]);
779: nvc += bss[i];
780: nv++;
781: if (nvc > bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of subvectors in subvectors > number of vectors in main vector");
782: if (nvc == bs) break;
783: }
785: n = n/bs;
787: jj = 0;
788: if (addv == INSERT_VALUES) {
789: for (j=0; j<nv; j++) {
790: for (k=0; k<bss[j]; k++) {
791: for (i=0; i<n; i++) x[bs*i+jj+k] = y[j][i*bss[j] + k];
792: }
793: jj += bss[j];
794: }
795: } else if (addv == ADD_VALUES) {
796: for (j=0; j<nv; j++) {
797: for (k=0; k<bss[j]; k++) {
798: for (i=0; i<n; i++) x[bs*i+jj+k] += y[j][i*bss[j] + k];
799: }
800: jj += bss[j];
801: }
802: #if !defined(PETSC_USE_COMPLEX)
803: } else if (addv == MAX_VALUES) {
804: for (j=0; j<nv; j++) {
805: for (k=0; k<bss[j]; k++) {
806: for (i=0; i<n; i++) x[bs*i+jj+k] = PetscMax(x[bs*i+jj+k],y[j][i*bss[j] + k]);
807: }
808: jj += bss[j];
809: }
810: #endif
811: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
813: VecRestoreArray(v,&x);
814: for (i=0; i<nv; i++) {
815: VecRestoreArray(s[i],&y[i]);
816: }
817: PetscFree2(y,bss);
818: return(0);
819: }
823: /*@
824: VecStrideGather - Gathers a single component from a multi-component vector into
825: another vector.
827: Collective on Vec
829: Input Parameter:
830: + v - the vector
831: . start - starting point of the subvector (defined by a stride)
832: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
834: Output Parameter:
835: . s - the location where the subvector is stored
837: Notes:
838: One must call VecSetBlockSize() before this routine to set the stride
839: information, or use a vector created from a multicomponent DMDA.
841: If x is the array representing the vector x then this gathers
842: the array (x[start],x[start+stride],x[start+2*stride], ....)
844: The parallel layout of the vector and the subvector must be the same;
845: i.e., nlocal of v = stride*(nlocal of s)
847: Not optimized; could be easily
849: Level: advanced
851: Concepts: gather^into strided vector
853: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
854: VecStrideScatterAll()
855: @*/
856: PetscErrorCode VecStrideGather(Vec v,PetscInt start,Vec s,InsertMode addv)
857: {
863: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
864: 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);
865: if (!v->ops->stridegather) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
866: (*v->ops->stridegather)(v,start,s,addv);
867: return(0);
868: }
872: /*@
873: VecStrideScatter - Scatters a single component from a vector into a multi-component vector.
875: Collective on Vec
877: Input Parameter:
878: + s - the single-component vector
879: . start - starting point of the subvector (defined by a stride)
880: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
882: Output Parameter:
883: . v - the location where the subvector is scattered (the multi-component vector)
885: Notes:
886: One must call VecSetBlockSize() on the multi-component vector before this
887: routine to set the stride information, or use a vector created from a multicomponent DMDA.
889: The parallel layout of the vector and the subvector must be the same;
890: i.e., nlocal of v = stride*(nlocal of s)
892: Not optimized; could be easily
894: Level: advanced
896: Concepts: scatter^into strided vector
898: .seealso: VecStrideNorm(), VecStrideGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
899: VecStrideScatterAll(), VecStrideSubSetScatter(), VecStrideSubSetGather()
900: @*/
901: PetscErrorCode VecStrideScatter(Vec s,PetscInt start,Vec v,InsertMode addv)
902: {
908: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
909: 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);
910: if (!v->ops->stridescatter) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
911: (*v->ops->stridescatter)(s,start,v,addv);
912: return(0);
913: }
917: /*@
918: VecStrideSubSetGather - Gathers a subset of components from a multi-component vector into
919: another vector.
921: Collective on Vec
923: Input Parameter:
924: + v - the vector
925: . nidx - the number of indices
926: . idxv - the indices of the components 0 <= idxv[0] ...idxv[nidx-1] < bs(v), they need not be sorted
927: . 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
928: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
930: Output Parameter:
931: . s - the location where the subvector is stored
933: Notes:
934: One must call VecSetBlockSize() on both vectors before this routine to set the stride
935: information, or use a vector created from a multicomponent DMDA.
938: The parallel layout of the vector and the subvector must be the same;
940: Not optimized; could be easily
942: Level: advanced
944: Concepts: gather^into strided vector
946: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideGather(), VecStrideSubSetScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
947: VecStrideScatterAll()
948: @*/
949: PetscErrorCode VecStrideSubSetGather(Vec v,PetscInt nidx,const PetscInt idxv[],const PetscInt idxs[],Vec s,InsertMode addv)
950: {
956: if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
957: if (!v->ops->stridesubsetgather) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
958: (*v->ops->stridesubsetgather)(v,nidx,idxv,idxs,s,addv);
959: return(0);
960: }
964: /*@
965: VecStrideSubSetScatter - Scatters components from a vector into a subset of components of a multi-component vector.
967: Collective on Vec
969: Input Parameter:
970: + s - the smaller-component vector
971: . nidx - the number of indices in idx
972: . 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
973: . idxv - the indices of the components in the larger-component vector, 0 <= idx[0] ...idx[nidx-1] < bs(v) they need not be sorted
974: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
976: Output Parameter:
977: . v - the location where the subvector is into scattered (the multi-component vector)
979: Notes:
980: One must call VecSetBlockSize() on the vectors before this
981: routine to set the stride information, or use a vector created from a multicomponent DMDA.
983: The parallel layout of the vector and the subvector must be the same;
985: Not optimized; could be easily
987: Level: advanced
989: Concepts: scatter^into strided vector
991: .seealso: VecStrideNorm(), VecStrideGather(), VecStrideGather(), VecStrideSubSetGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
992: VecStrideScatterAll()
993: @*/
994: PetscErrorCode VecStrideSubSetScatter(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)
995: {
1001: if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
1002: if (!v->ops->stridesubsetscatter) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
1003: (*v->ops->stridesubsetscatter)(s,nidx,idxs,idxv,v,addv);
1004: return(0);
1005: }
1009: PetscErrorCode VecStrideGather_Default(Vec v,PetscInt start,Vec s,InsertMode addv)
1010: {
1012: PetscInt i,n,bs,ns;
1013: const PetscScalar *x;
1014: PetscScalar *y;
1017: VecGetLocalSize(v,&n);
1018: VecGetLocalSize(s,&ns);
1019: VecGetArrayRead(v,&x);
1020: VecGetArray(s,&y);
1022: bs = v->map->bs;
1023: 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);
1024: x += start;
1025: n = n/bs;
1027: if (addv == INSERT_VALUES) {
1028: for (i=0; i<n; i++) y[i] = x[bs*i];
1029: } else if (addv == ADD_VALUES) {
1030: for (i=0; i<n; i++) y[i] += x[bs*i];
1031: #if !defined(PETSC_USE_COMPLEX)
1032: } else if (addv == MAX_VALUES) {
1033: for (i=0; i<n; i++) y[i] = PetscMax(y[i],x[bs*i]);
1034: #endif
1035: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
1037: VecRestoreArrayRead(v,&x);
1038: VecRestoreArray(s,&y);
1039: return(0);
1040: }
1044: PetscErrorCode VecStrideScatter_Default(Vec s,PetscInt start,Vec v,InsertMode addv)
1045: {
1046: PetscErrorCode ierr;
1047: PetscInt i,n,bs,ns;
1048: PetscScalar *x;
1049: const PetscScalar *y;
1052: VecGetLocalSize(v,&n);
1053: VecGetLocalSize(s,&ns);
1054: VecGetArray(v,&x);
1055: VecGetArrayRead(s,&y);
1057: bs = v->map->bs;
1058: 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);
1059: x += start;
1060: n = n/bs;
1062: if (addv == INSERT_VALUES) {
1063: for (i=0; i<n; i++) x[bs*i] = y[i];
1064: } else if (addv == ADD_VALUES) {
1065: for (i=0; i<n; i++) x[bs*i] += y[i];
1066: #if !defined(PETSC_USE_COMPLEX)
1067: } else if (addv == MAX_VALUES) {
1068: for (i=0; i<n; i++) x[bs*i] = PetscMax(y[i],x[bs*i]);
1069: #endif
1070: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
1072: VecRestoreArray(v,&x);
1073: VecRestoreArrayRead(s,&y);
1074: return(0);
1075: }
1079: PetscErrorCode VecStrideSubSetGather_Default(Vec v,PetscInt nidx,const PetscInt idxv[],const PetscInt idxs[],Vec s,InsertMode addv)
1080: {
1081: PetscErrorCode ierr;
1082: PetscInt i,j,n,bs,bss,ns;
1083: const PetscScalar *x;
1084: PetscScalar *y;
1087: VecGetLocalSize(v,&n);
1088: VecGetLocalSize(s,&ns);
1089: VecGetArrayRead(v,&x);
1090: VecGetArray(s,&y);
1092: bs = v->map->bs;
1093: bss = s->map->bs;
1094: n = n/bs;
1096: #if defined(PETSC_DEBUG)
1097: if (n != ns/bss) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible layout of vectors");
1098: for (j=0; j<nidx; j++) {
1099: if (idxv[j] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is negative",j,idxv[j]);
1100: 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);
1101: }
1102: if (!idxs && bss != nidx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide idxs when not gathering into all locations");
1103: #endif
1105: if (addv == INSERT_VALUES) {
1106: if (!idxs) {
1107: for (i=0; i<n; i++) {
1108: for (j=0; j<bss; j++) y[bss*i+j] = x[bs*i+idxv[j]];
1109: }
1110: } else {
1111: for (i=0; i<n; i++) {
1112: for (j=0; j<bss; j++) y[bss*i+idxs[j]] = x[bs*i+idxv[j]];
1113: }
1114: }
1115: } else if (addv == ADD_VALUES) {
1116: if (!idxs) {
1117: for (i=0; i<n; i++) {
1118: for (j=0; j<bss; j++) y[bss*i+j] += x[bs*i+idxv[j]];
1119: }
1120: } else {
1121: for (i=0; i<n; i++) {
1122: for (j=0; j<bss; j++) y[bss*i+idxs[j]] += x[bs*i+idxv[j]];
1123: }
1124: }
1125: #if !defined(PETSC_USE_COMPLEX)
1126: } else if (addv == MAX_VALUES) {
1127: if (!idxs) {
1128: for (i=0; i<n; i++) {
1129: for (j=0; j<bss; j++) y[bss*i+j] = PetscMax(y[bss*i+j],x[bs*i+idxv[j]]);
1130: }
1131: } else {
1132: for (i=0; i<n; i++) {
1133: for (j=0; j<bss; j++) y[bss*i+idxs[j]] = PetscMax(y[bss*i+idxs[j]],x[bs*i+idxv[j]]);
1134: }
1135: }
1136: #endif
1137: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
1139: VecRestoreArrayRead(v,&x);
1140: VecRestoreArray(s,&y);
1141: return(0);
1142: }
1146: PetscErrorCode VecStrideSubSetScatter_Default(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)
1147: {
1148: PetscErrorCode ierr;
1149: PetscInt j,i,n,bs,ns,bss;
1150: PetscScalar *x;
1151: const PetscScalar *y;
1154: VecGetLocalSize(v,&n);
1155: VecGetLocalSize(s,&ns);
1156: VecGetArray(v,&x);
1157: VecGetArrayRead(s,&y);
1159: bs = v->map->bs;
1160: bss = s->map->bs;
1161: n = n/bs;
1163: #if defined(PETSC_DEBUG)
1164: if (n != ns/bss) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible layout of vectors");
1165: for (j=0; j<bss; j++) {
1166: if (idx[j] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is negative",j,idx[j]);
1167: 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);
1168: }
1169: if (!idxs && bss != nidx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide idxs when not scattering from all locations");
1170: #endif
1172: if (addv == INSERT_VALUES) {
1173: if (!idxs) {
1174: for (i=0; i<n; i++) {
1175: for (j=0; j<bss; j++) x[bs*i + idxv[j]] = y[bss*i+j];
1176: }
1177: } else {
1178: for (i=0; i<n; i++) {
1179: for (j=0; j<bss; j++) x[bs*i + idxv[j]] = y[bss*i+idxs[j]];
1180: }
1181: }
1182: } else if (addv == ADD_VALUES) {
1183: if (!idxs) {
1184: for (i=0; i<n; i++) {
1185: for (j=0; j<bss; j++) x[bs*i + idxv[j]] += y[bss*i+j];
1186: }
1187: } else {
1188: for (i=0; i<n; i++) {
1189: for (j=0; j<bss; j++) x[bs*i + idxv[j]] += y[bss*i+idxs[j]];
1190: }
1191: }
1192: #if !defined(PETSC_USE_COMPLEX)
1193: } else if (addv == MAX_VALUES) {
1194: if (!idxs) {
1195: for (i=0; i<n; i++) {
1196: for (j=0; j<bss; j++) x[bs*i + idxv[j]] = PetscMax(y[bss*i+j],x[bs*i + idxv[j]]);
1197: }
1198: } else {
1199: for (i=0; i<n; i++) {
1200: for (j=0; j<bss; j++) x[bs*i + idxv[j]] = PetscMax(y[bss*i+idxs[j]],x[bs*i + idxv[j]]);
1201: }
1202: }
1203: #endif
1204: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
1206: VecRestoreArray(v,&x);
1207: VecRestoreArrayRead(s,&y);
1208: return(0);
1209: }
1213: PetscErrorCode VecReciprocal_Default(Vec v)
1214: {
1216: PetscInt i,n;
1217: PetscScalar *x;
1220: VecGetLocalSize(v,&n);
1221: VecGetArray(v,&x);
1222: for (i=0; i<n; i++) {
1223: if (x[i] != (PetscScalar)0.0) x[i] = (PetscScalar)1.0/x[i];
1224: }
1225: VecRestoreArray(v,&x);
1226: return(0);
1227: }
1231: /*@
1232: VecExp - Replaces each component of a vector by e^x_i
1234: Not collective
1236: Input Parameter:
1237: . v - The vector
1239: Output Parameter:
1240: . v - The vector of exponents
1242: Level: beginner
1244: .seealso: VecLog(), VecAbs(), VecSqrtAbs(), VecReciprocal()
1246: .keywords: vector, sqrt, square root
1247: @*/
1248: PetscErrorCode VecExp(Vec v)
1249: {
1250: PetscScalar *x;
1251: PetscInt i, n;
1256: if (v->ops->exp) {
1257: (*v->ops->exp)(v);
1258: } else {
1259: VecGetLocalSize(v, &n);
1260: VecGetArray(v, &x);
1261: for (i = 0; i < n; i++) x[i] = PetscExpScalar(x[i]);
1262: VecRestoreArray(v, &x);
1263: }
1264: return(0);
1265: }
1269: /*@
1270: VecLog - Replaces each component of a vector by log(x_i), the natural logarithm
1272: Not collective
1274: Input Parameter:
1275: . v - The vector
1277: Output Parameter:
1278: . v - The vector of logs
1280: Level: beginner
1282: .seealso: VecExp(), VecAbs(), VecSqrtAbs(), VecReciprocal()
1284: .keywords: vector, sqrt, square root
1285: @*/
1286: PetscErrorCode VecLog(Vec v)
1287: {
1288: PetscScalar *x;
1289: PetscInt i, n;
1294: if (v->ops->log) {
1295: (*v->ops->log)(v);
1296: } else {
1297: VecGetLocalSize(v, &n);
1298: VecGetArray(v, &x);
1299: for (i = 0; i < n; i++) x[i] = PetscLogScalar(x[i]);
1300: VecRestoreArray(v, &x);
1301: }
1302: return(0);
1303: }
1307: /*@
1308: VecSqrtAbs - Replaces each component of a vector by the square root of its magnitude.
1310: Not collective
1312: Input Parameter:
1313: . v - The vector
1315: Output Parameter:
1316: . v - The vector square root
1318: Level: beginner
1320: Note: The actual function is sqrt(|x_i|)
1322: .seealso: VecLog(), VecExp(), VecReciprocal(), VecAbs()
1324: .keywords: vector, sqrt, square root
1325: @*/
1326: PetscErrorCode VecSqrtAbs(Vec v)
1327: {
1328: PetscScalar *x;
1329: PetscInt i, n;
1334: if (v->ops->sqrt) {
1335: (*v->ops->sqrt)(v);
1336: } else {
1337: VecGetLocalSize(v, &n);
1338: VecGetArray(v, &x);
1339: for (i = 0; i < n; i++) x[i] = PetscSqrtReal(PetscAbsScalar(x[i]));
1340: VecRestoreArray(v, &x);
1341: }
1342: return(0);
1343: }
1347: /*@
1348: VecDotNorm2 - computes the inner product of two vectors and the 2-norm squared of the second vector
1350: Collective on Vec
1352: Input Parameter:
1353: + s - first vector
1354: - t - second vector
1356: Output Parameter:
1357: + dp - s'conj(t)
1358: - nm - t'conj(t)
1360: Level: advanced
1362: Notes: conj(x) is the complex conjugate of x when x is complex
1365: .seealso: VecDot(), VecNorm(), VecDotBegin(), VecNormBegin(), VecDotEnd(), VecNormEnd()
1367: .keywords: vector, sqrt, square root
1368: @*/
1369: PetscErrorCode VecDotNorm2(Vec s,Vec t,PetscScalar *dp, PetscReal *nm)
1370: {
1371: const PetscScalar *sx, *tx;
1372: PetscScalar dpx = 0.0, nmx = 0.0,work[2],sum[2];
1373: PetscInt i, n;
1374: PetscErrorCode ierr;
1384: if (s->map->N != t->map->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
1385: if (s->map->n != t->map->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");
1387: PetscLogEventBarrierBegin(VEC_DotNormBarrier,s,t,0,0,PetscObjectComm((PetscObject)s));
1388: if (s->ops->dotnorm2) {
1389: (*s->ops->dotnorm2)(s,t,dp,&dpx);
1390: *nm = PetscRealPart(dpx);
1391: } else {
1392: VecGetLocalSize(s, &n);
1393: VecGetArrayRead(s, &sx);
1394: VecGetArrayRead(t, &tx);
1396: for (i = 0; i<n; i++) {
1397: dpx += sx[i]*PetscConj(tx[i]);
1398: nmx += tx[i]*PetscConj(tx[i]);
1399: }
1400: work[0] = dpx;
1401: work[1] = nmx;
1403: MPIU_Allreduce(work,sum,2,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)s));
1404: *dp = sum[0];
1405: *nm = PetscRealPart(sum[1]);
1407: VecRestoreArrayRead(t, &tx);
1408: VecRestoreArrayRead(s, &sx);
1409: PetscLogFlops(4.0*n);
1410: }
1411: PetscLogEventBarrierEnd(VEC_DotNormBarrier,s,t,0,0,PetscObjectComm((PetscObject)s));
1412: return(0);
1413: }
1417: /*@
1418: VecSum - Computes the sum of all the components of a vector.
1420: Collective on Vec
1422: Input Parameter:
1423: . v - the vector
1425: Output Parameter:
1426: . sum - the result
1428: Level: beginner
1430: Concepts: sum^of vector entries
1432: .seealso: VecNorm()
1433: @*/
1434: PetscErrorCode VecSum(Vec v,PetscScalar *sum)
1435: {
1436: PetscErrorCode ierr;
1437: PetscInt i,n;
1438: const PetscScalar *x;
1439: PetscScalar lsum = 0.0;
1444: VecGetLocalSize(v,&n);
1445: VecGetArrayRead(v,&x);
1446: for (i=0; i<n; i++) lsum += x[i];
1447: MPIU_Allreduce(&lsum,sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)v));
1448: VecRestoreArrayRead(v,&x);
1449: return(0);
1450: }
1454: /*@
1455: VecShift - Shifts all of the components of a vector by computing
1456: x[i] = x[i] + shift.
1458: Logically Collective on Vec
1460: Input Parameters:
1461: + v - the vector
1462: - shift - the shift
1464: Output Parameter:
1465: . v - the shifted vector
1467: Level: intermediate
1469: Concepts: vector^adding constant
1471: @*/
1472: PetscErrorCode VecShift(Vec v,PetscScalar shift)
1473: {
1475: PetscInt i,n;
1476: PetscScalar *x;
1481: VecLocked(v,1);
1483: if (v->ops->shift) {
1484: (*v->ops->shift)(v,shift);
1485: } else {
1486: VecGetLocalSize(v,&n);
1487: VecGetArray(v,&x);
1488: for (i=0; i<n; i++) x[i] += shift;
1489: VecRestoreArray(v,&x);
1490: }
1491: return(0);
1492: }
1496: /*@
1497: VecAbs - Replaces every element in a vector with its absolute value.
1499: Logically Collective on Vec
1501: Input Parameters:
1502: . v - the vector
1504: Level: intermediate
1506: Concepts: vector^absolute value
1508: @*/
1509: PetscErrorCode VecAbs(Vec v)
1510: {
1512: PetscInt i,n;
1513: PetscScalar *x;
1517: VecLocked(v,1);
1518:
1519: if (v->ops->abs) {
1520: (*v->ops->abs)(v);
1521: } else {
1522: VecGetLocalSize(v,&n);
1523: VecGetArray(v,&x);
1524: for (i=0; i<n; i++) x[i] = PetscAbsScalar(x[i]);
1525: VecRestoreArray(v,&x);
1526: }
1527: return(0);
1528: }
1532: /*@
1533: VecPermute - Permutes a vector in place using the given ordering.
1535: Input Parameters:
1536: + vec - The vector
1537: . order - The ordering
1538: - inv - The flag for inverting the permutation
1540: Level: beginner
1542: Note: This function does not yet support parallel Index Sets with non-local permutations
1544: .seealso: MatPermute()
1545: .keywords: vec, permute
1546: @*/
1547: PetscErrorCode VecPermute(Vec x, IS row, PetscBool inv)
1548: {
1549: PetscScalar *array, *newArray;
1550: const PetscInt *idx;
1551: PetscInt i,rstart,rend;
1555: VecLocked(x,1);
1556:
1557: VecGetOwnershipRange(x,&rstart,&rend);
1558: ISGetIndices(row, &idx);
1559: VecGetArray(x, &array);
1560: PetscMalloc1(x->map->n, &newArray);
1561: #if defined(PETSC_USE_DEBUG)
1562: for (i = 0; i < x->map->n; i++) {
1563: 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]);
1564: }
1565: #endif
1566: if (!inv) {
1567: for (i = 0; i < x->map->n; i++) newArray[i] = array[idx[i]-rstart];
1568: } else {
1569: for (i = 0; i < x->map->n; i++) newArray[idx[i]-rstart] = array[i];
1570: }
1571: VecRestoreArray(x, &array);
1572: ISRestoreIndices(row, &idx);
1573: VecReplaceArray(x, newArray);
1574: return(0);
1575: }
1579: /*@
1580: VecEqual - Compares two vectors. Returns true if the two vectors are either pointing to the same memory buffer,
1581: or if the two vectors have the same local and global layout as well as bitwise equality of all entries.
1582: Does NOT take round-off errors into account.
1584: Collective on Vec
1586: Input Parameters:
1587: + vec1 - the first vector
1588: - vec2 - the second vector
1590: Output Parameter:
1591: . flg - PETSC_TRUE if the vectors are equal; PETSC_FALSE otherwise.
1593: Level: intermediate
1595: Concepts: equal^two vectors
1596: Concepts: vector^equality
1598: @*/
1599: PetscErrorCode VecEqual(Vec vec1,Vec vec2,PetscBool *flg)
1600: {
1601: const PetscScalar *v1,*v2;
1602: PetscErrorCode ierr;
1603: PetscInt n1,n2,N1,N2;
1604: PetscBool flg1;
1610: if (vec1 == vec2) *flg = PETSC_TRUE;
1611: else {
1612: VecGetSize(vec1,&N1);
1613: VecGetSize(vec2,&N2);
1614: if (N1 != N2) flg1 = PETSC_FALSE;
1615: else {
1616: VecGetLocalSize(vec1,&n1);
1617: VecGetLocalSize(vec2,&n2);
1618: if (n1 != n2) flg1 = PETSC_FALSE;
1619: else {
1620: VecGetArrayRead(vec1,&v1);
1621: VecGetArrayRead(vec2,&v2);
1622: PetscMemcmp(v1,v2,n1*sizeof(PetscScalar),&flg1);
1623: VecRestoreArrayRead(vec1,&v1);
1624: VecRestoreArrayRead(vec2,&v2);
1625: }
1626: }
1627: /* combine results from all processors */
1628: MPIU_Allreduce(&flg1,flg,1,MPIU_BOOL,MPI_MIN,PetscObjectComm((PetscObject)vec1));
1629: }
1630: return(0);
1631: }
1635: /*@
1636: VecUniqueEntries - Compute the number of unique entries, and those entries
1638: Collective on Vec
1640: Input Parameter:
1641: . vec - the vector
1643: Output Parameters:
1644: + n - The number of unique entries
1645: - e - The entries
1647: Level: intermediate
1649: @*/
1650: PetscErrorCode VecUniqueEntries(Vec vec, PetscInt *n, PetscScalar **e)
1651: {
1652: const PetscScalar *v;
1653: PetscScalar *tmp, *vals;
1654: PetscMPIInt *N, *displs, l;
1655: PetscInt ng, m, i, j, p;
1656: PetscMPIInt size;
1657: PetscErrorCode ierr;
1662: MPI_Comm_size(PetscObjectComm((PetscObject) vec), &size);
1663: VecGetLocalSize(vec, &m);
1664: VecGetArrayRead(vec, &v);
1665: PetscMalloc2(m,&tmp,size,&N);
1666: for (i = 0, j = 0, l = 0; i < m; ++i) {
1667: /* Can speed this up with sorting */
1668: for (j = 0; j < l; ++j) {
1669: if (v[i] == tmp[j]) break;
1670: }
1671: if (j == l) {
1672: tmp[j] = v[i];
1673: ++l;
1674: }
1675: }
1676: VecRestoreArrayRead(vec, &v);
1677: /* Gather serial results */
1678: MPI_Allgather(&l, 1, MPI_INT, N, 1, MPI_INT, PetscObjectComm((PetscObject) vec));
1679: for (p = 0, ng = 0; p < size; ++p) {
1680: ng += N[p];
1681: }
1682: PetscMalloc2(ng,&vals,size+1,&displs);
1683: for (p = 1, displs[0] = 0; p <= size; ++p) {
1684: displs[p] = displs[p-1] + N[p-1];
1685: }
1686: MPI_Allgatherv(tmp, l, MPIU_SCALAR, vals, N, displs, MPIU_SCALAR, PetscObjectComm((PetscObject) vec));
1687: /* Find unique entries */
1688: #ifdef PETSC_USE_COMPLEX
1689: SETERRQ(PetscObjectComm((PetscObject) vec), PETSC_ERR_SUP, "Does not work with complex numbers");
1690: #else
1691: *n = displs[size];
1692: PetscSortRemoveDupsReal(n, (PetscReal *) vals);
1693: if (e) {
1695: PetscMalloc1(*n, e);
1696: for (i = 0; i < *n; ++i) {
1697: (*e)[i] = vals[i];
1698: }
1699: }
1700: PetscFree2(vals,displs);
1701: PetscFree2(tmp,N);
1702: return(0);
1703: #endif
1704: }