Actual source code: vinv.c
petsc-3.10.5 2019-03-28
2: /*
3: Some useful vector utility functions.
4: */
5: #include <../src/vec/vec/impls/mpi/pvecimpl.h>
6: extern MPI_Op MPIU_MAXINDEX_OP;
7: extern MPI_Op MPIU_MININDEX_OP;
9: /*@
10: VecStrideSet - Sets a subvector of a vector defined
11: by a starting point and a stride with a given value
13: Logically Collective on Vec
15: Input Parameter:
16: + v - the vector
17: . start - starting point of the subvector (defined by a stride)
18: - s - value to set for each entry in that subvector
20: Notes:
21: One must call VecSetBlockSize() before this routine to set the stride
22: information, or use a vector created from a multicomponent DMDA.
24: This will only work if the desire subvector is a stride subvector
26: Level: advanced
28: Concepts: scale^on stride of vector
29: Concepts: stride^scale
31: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideScale()
32: @*/
33: PetscErrorCode VecStrideSet(Vec v,PetscInt start,PetscScalar s)
34: {
36: PetscInt i,n,bs;
37: PetscScalar *x;
43: VecLocked(v,1);
44:
45: VecGetLocalSize(v,&n);
46: VecGetArray(v,&x);
47: VecGetBlockSize(v,&bs);
48: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
49: 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);
50: x += start;
52: for (i=0; i<n; i+=bs) x[i] = s;
53: x -= start;
55: VecRestoreArray(v,&x);
56: return(0);
57: }
59: /*@
60: VecStrideScale - Scales a subvector of a vector defined
61: by a starting point and a stride.
63: Logically Collective on Vec
65: Input Parameter:
66: + v - the vector
67: . start - starting point of the subvector (defined by a stride)
68: - scale - value to multiply each subvector entry by
70: Notes:
71: One must call VecSetBlockSize() before this routine to set the stride
72: information, or use a vector created from a multicomponent DMDA.
74: This will only work if the desire subvector is a stride subvector
76: Level: advanced
78: Concepts: scale^on stride of vector
79: Concepts: stride^scale
81: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideScale()
82: @*/
83: PetscErrorCode VecStrideScale(Vec v,PetscInt start,PetscScalar scale)
84: {
86: PetscInt i,n,bs;
87: PetscScalar *x;
93: VecLocked(v,1);
94:
95: VecGetLocalSize(v,&n);
96: VecGetArray(v,&x);
97: VecGetBlockSize(v,&bs);
98: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
99: 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);
100: x += start;
102: for (i=0; i<n; i+=bs) x[i] *= scale;
103: x -= start;
105: VecRestoreArray(v,&x);
106: return(0);
107: }
109: /*@
110: VecStrideNorm - Computes the norm of subvector of a vector defined
111: by a starting point and a stride.
113: Collective on Vec
115: Input Parameter:
116: + v - the vector
117: . start - starting point of the subvector (defined by a stride)
118: - ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY
120: Output Parameter:
121: . norm - the norm
123: Notes:
124: One must call VecSetBlockSize() before this routine to set the stride
125: information, or use a vector created from a multicomponent DMDA.
127: If x is the array representing the vector x then this computes the norm
128: of the array (x[start],x[start+stride],x[start+2*stride], ....)
130: This is useful for computing, say the norm of the pressure variable when
131: the pressure is stored (interlaced) with other variables, say density etc.
133: This will only work if the desire subvector is a stride subvector
135: Level: advanced
137: Concepts: norm^on stride of vector
138: Concepts: stride^norm
140: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
141: @*/
142: PetscErrorCode VecStrideNorm(Vec v,PetscInt start,NormType ntype,PetscReal *nrm)
143: {
144: PetscErrorCode ierr;
145: PetscInt i,n,bs;
146: const PetscScalar *x;
147: PetscReal tnorm;
148: MPI_Comm comm;
153: VecGetLocalSize(v,&n);
154: VecGetArrayRead(v,&x);
155: PetscObjectGetComm((PetscObject)v,&comm);
157: VecGetBlockSize(v,&bs);
158: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
159: 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);
160: x += start;
162: if (ntype == NORM_2) {
163: PetscScalar sum = 0.0;
164: for (i=0; i<n; i+=bs) sum += x[i]*(PetscConj(x[i]));
165: tnorm = PetscRealPart(sum);
166: MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_SUM,comm);
167: *nrm = PetscSqrtReal(*nrm);
168: } else if (ntype == NORM_1) {
169: tnorm = 0.0;
170: for (i=0; i<n; i+=bs) tnorm += PetscAbsScalar(x[i]);
171: MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_SUM,comm);
172: } else if (ntype == NORM_INFINITY) {
173: PetscReal tmp;
174: tnorm = 0.0;
176: for (i=0; i<n; i+=bs) {
177: if ((tmp = PetscAbsScalar(x[i])) > tnorm) tnorm = tmp;
178: /* check special case of tmp == NaN */
179: if (tmp != tmp) {tnorm = tmp; break;}
180: }
181: MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_MAX,comm);
182: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
183: VecRestoreArrayRead(v,&x);
184: return(0);
185: }
187: /*@
188: VecStrideMax - Computes the maximum of subvector of a vector defined
189: by a starting point and a stride and optionally its location.
191: Collective on Vec
193: Input Parameter:
194: + v - the vector
195: - start - starting point of the subvector (defined by a stride)
197: Output Parameter:
198: + index - the location where the maximum occurred (pass NULL if not required)
199: - nrm - the maximum value in the subvector
201: Notes:
202: One must call VecSetBlockSize() before this routine to set the stride
203: information, or use a vector created from a multicomponent DMDA.
205: If xa is the array representing the vector x, then this computes the max
206: of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)
208: This is useful for computing, say the maximum of the pressure variable when
209: the pressure is stored (interlaced) with other variables, e.g., density, etc.
210: This will only work if the desire subvector is a stride subvector.
212: Level: advanced
214: Concepts: maximum^on stride of vector
215: Concepts: stride^maximum
217: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
218: @*/
219: PetscErrorCode VecStrideMax(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
220: {
221: PetscErrorCode ierr;
222: PetscInt i,n,bs,id;
223: const PetscScalar *x;
224: PetscReal max,tmp;
225: MPI_Comm comm;
231: VecGetLocalSize(v,&n);
232: VecGetArrayRead(v,&x);
233: PetscObjectGetComm((PetscObject)v,&comm);
235: VecGetBlockSize(v,&bs);
236: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
237: 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);
238: x += start;
240: id = -1;
241: if (!n) max = PETSC_MIN_REAL;
242: else {
243: id = 0;
244: max = PetscRealPart(x[0]);
245: for (i=bs; i<n; i+=bs) {
246: if ((tmp = PetscRealPart(x[i])) > max) { max = tmp; id = i;}
247: }
248: }
249: VecRestoreArrayRead(v,&x);
251: if (!idex) {
252: MPIU_Allreduce(&max,nrm,1,MPIU_REAL,MPIU_MAX,comm);
253: } else {
254: PetscReal in[2],out[2];
255: PetscInt rstart;
257: VecGetOwnershipRange(v,&rstart,NULL);
258: in[0] = max;
259: in[1] = rstart+id+start;
260: MPIU_Allreduce(in,out,2,MPIU_REAL,MPIU_MAXINDEX_OP,PetscObjectComm((PetscObject)v));
261: *nrm = out[0];
262: *idex = (PetscInt)out[1];
263: }
264: return(0);
265: }
267: /*@
268: VecStrideMin - Computes the minimum of subvector of a vector defined
269: by a starting point and a stride and optionally its location.
271: Collective on Vec
273: Input Parameter:
274: + v - the vector
275: - start - starting point of the subvector (defined by a stride)
277: Output Parameter:
278: + idex - the location where the minimum occurred. (pass NULL if not required)
279: - nrm - the minimum value in the subvector
281: Level: advanced
283: Notes:
284: One must call VecSetBlockSize() before this routine to set the stride
285: information, or use a vector created from a multicomponent DMDA.
287: If xa is the array representing the vector x, then this computes the min
288: of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)
290: This is useful for computing, say the minimum of the pressure variable when
291: the pressure is stored (interlaced) with other variables, e.g., density, etc.
292: This will only work if the desire subvector is a stride subvector.
294: Concepts: minimum^on stride of vector
295: Concepts: stride^minimum
297: .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
298: @*/
299: PetscErrorCode VecStrideMin(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
300: {
301: PetscErrorCode ierr;
302: PetscInt i,n,bs,id;
303: const PetscScalar *x;
304: PetscReal min,tmp;
305: MPI_Comm comm;
311: VecGetLocalSize(v,&n);
312: VecGetArrayRead(v,&x);
313: PetscObjectGetComm((PetscObject)v,&comm);
315: VecGetBlockSize(v,&bs);
316: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
317: 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);
318: x += start;
320: id = -1;
321: if (!n) min = PETSC_MAX_REAL;
322: else {
323: id = 0;
324: min = PetscRealPart(x[0]);
325: for (i=bs; i<n; i+=bs) {
326: if ((tmp = PetscRealPart(x[i])) < min) { min = tmp; id = i;}
327: }
328: }
329: VecRestoreArrayRead(v,&x);
331: if (!idex) {
332: MPIU_Allreduce(&min,nrm,1,MPIU_REAL,MPIU_MIN,comm);
333: } else {
334: PetscReal in[2],out[2];
335: PetscInt rstart;
337: VecGetOwnershipRange(v,&rstart,NULL);
338: in[0] = min;
339: in[1] = rstart+id;
340: MPIU_Allreduce(in,out,2,MPIU_REAL,MPIU_MININDEX_OP,PetscObjectComm((PetscObject)v));
341: *nrm = out[0];
342: *idex = (PetscInt)out[1];
343: }
344: return(0);
345: }
347: /*@
348: VecStrideScaleAll - Scales the subvectors of a vector defined
349: by a starting point and a stride.
351: Logically Collective on Vec
353: Input Parameter:
354: + v - the vector
355: - scales - values to multiply each subvector entry by
357: Notes:
358: One must call VecSetBlockSize() before this routine to set the stride
359: information, or use a vector created from a multicomponent DMDA.
361: The dimension of scales must be the same as the vector block size
364: Level: advanced
366: Concepts: scale^on stride of vector
367: Concepts: stride^scale
369: .seealso: VecNorm(), VecStrideScale(), VecScale(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
370: @*/
371: PetscErrorCode VecStrideScaleAll(Vec v,const PetscScalar *scales)
372: {
374: PetscInt i,j,n,bs;
375: PetscScalar *x;
380: VecLocked(v,1);
381:
382: VecGetLocalSize(v,&n);
383: VecGetArray(v,&x);
384: VecGetBlockSize(v,&bs);
386: /* need to provide optimized code for each bs */
387: for (i=0; i<n; i+=bs) {
388: for (j=0; j<bs; j++) x[i+j] *= scales[j];
389: }
390: VecRestoreArray(v,&x);
391: return(0);
392: }
394: /*@
395: VecStrideNormAll - Computes the norms of subvectors of a vector defined
396: by a starting point and a stride.
398: Collective on Vec
400: Input Parameter:
401: + v - the vector
402: - ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY
404: Output Parameter:
405: . nrm - the norms
407: Notes:
408: One must call VecSetBlockSize() before this routine to set the stride
409: information, or use a vector created from a multicomponent DMDA.
411: If x is the array representing the vector x then this computes the norm
412: of the array (x[start],x[start+stride],x[start+2*stride], ....) for each start < stride
414: The dimension of nrm must be the same as the vector block size
416: This will only work if the desire subvector is a stride subvector
418: Level: advanced
420: Concepts: norm^on stride of vector
421: Concepts: stride^norm
423: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
424: @*/
425: PetscErrorCode VecStrideNormAll(Vec v,NormType ntype,PetscReal nrm[])
426: {
427: PetscErrorCode ierr;
428: PetscInt i,j,n,bs;
429: const PetscScalar *x;
430: PetscReal tnorm[128];
431: MPI_Comm comm;
436: VecGetLocalSize(v,&n);
437: VecGetArrayRead(v,&x);
438: PetscObjectGetComm((PetscObject)v,&comm);
440: VecGetBlockSize(v,&bs);
441: if (bs > 128) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");
443: if (ntype == NORM_2) {
444: PetscScalar sum[128];
445: for (j=0; j<bs; j++) sum[j] = 0.0;
446: for (i=0; i<n; i+=bs) {
447: for (j=0; j<bs; j++) sum[j] += x[i+j]*(PetscConj(x[i+j]));
448: }
449: for (j=0; j<bs; j++) tnorm[j] = PetscRealPart(sum[j]);
451: MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_SUM,comm);
452: for (j=0; j<bs; j++) nrm[j] = PetscSqrtReal(nrm[j]);
453: } else if (ntype == NORM_1) {
454: for (j=0; j<bs; j++) tnorm[j] = 0.0;
456: for (i=0; i<n; i+=bs) {
457: for (j=0; j<bs; j++) tnorm[j] += PetscAbsScalar(x[i+j]);
458: }
460: MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_SUM,comm);
461: } else if (ntype == NORM_INFINITY) {
462: PetscReal tmp;
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++) {
467: if ((tmp = PetscAbsScalar(x[i+j])) > tnorm[j]) tnorm[j] = tmp;
468: /* check special case of tmp == NaN */
469: if (tmp != tmp) {tnorm[j] = tmp; break;}
470: }
471: }
472: MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_MAX,comm);
473: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
474: VecRestoreArrayRead(v,&x);
475: return(0);
476: }
478: /*@
479: VecStrideMaxAll - Computes the maximums of subvectors of a vector defined
480: by a starting point and a stride and optionally its location.
482: Collective on Vec
484: Input Parameter:
485: . v - the vector
487: Output Parameter:
488: + index - the location where the maximum occurred (not supported, pass NULL,
489: if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
490: - nrm - the maximum values of each subvector
492: Notes:
493: One must call VecSetBlockSize() before this routine to set the stride
494: information, or use a vector created from a multicomponent DMDA.
496: The dimension of nrm must be the same as the vector block size
498: Level: advanced
500: Concepts: maximum^on stride of vector
501: Concepts: stride^maximum
503: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
504: @*/
505: PetscErrorCode VecStrideMaxAll(Vec v,PetscInt idex[],PetscReal nrm[])
506: {
507: PetscErrorCode ierr;
508: PetscInt i,j,n,bs;
509: const PetscScalar *x;
510: PetscReal max[128],tmp;
511: MPI_Comm comm;
516: 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");
517: VecGetLocalSize(v,&n);
518: VecGetArrayRead(v,&x);
519: PetscObjectGetComm((PetscObject)v,&comm);
521: VecGetBlockSize(v,&bs);
522: if (bs > 128) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");
524: if (!n) {
525: for (j=0; j<bs; j++) max[j] = PETSC_MIN_REAL;
526: } else {
527: for (j=0; j<bs; j++) max[j] = PetscRealPart(x[j]);
529: for (i=bs; i<n; i+=bs) {
530: for (j=0; j<bs; j++) {
531: if ((tmp = PetscRealPart(x[i+j])) > max[j]) max[j] = tmp;
532: }
533: }
534: }
535: MPIU_Allreduce(max,nrm,bs,MPIU_REAL,MPIU_MAX,comm);
537: VecRestoreArrayRead(v,&x);
538: return(0);
539: }
541: /*@
542: VecStrideMinAll - Computes the minimum of subvector of a vector defined
543: by a starting point and a stride and optionally its location.
545: Collective on Vec
547: Input Parameter:
548: . v - the vector
550: Output Parameter:
551: + idex - the location where the minimum occurred (not supported, pass NULL,
552: if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
553: - nrm - the minimums of each subvector
555: Level: advanced
557: Notes:
558: One must call VecSetBlockSize() before this routine to set the stride
559: information, or use a vector created from a multicomponent DMDA.
561: The dimension of nrm must be the same as the vector block size
563: Concepts: minimum^on stride of vector
564: Concepts: stride^minimum
566: .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
567: @*/
568: PetscErrorCode VecStrideMinAll(Vec v,PetscInt idex[],PetscReal nrm[])
569: {
570: PetscErrorCode ierr;
571: PetscInt i,n,bs,j;
572: const PetscScalar *x;
573: PetscReal min[128],tmp;
574: MPI_Comm comm;
579: 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");
580: VecGetLocalSize(v,&n);
581: VecGetArrayRead(v,&x);
582: PetscObjectGetComm((PetscObject)v,&comm);
584: VecGetBlockSize(v,&bs);
585: if (bs > 128) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");
587: if (!n) {
588: for (j=0; j<bs; j++) min[j] = PETSC_MAX_REAL;
589: } else {
590: for (j=0; j<bs; j++) min[j] = PetscRealPart(x[j]);
592: for (i=bs; i<n; i+=bs) {
593: for (j=0; j<bs; j++) {
594: if ((tmp = PetscRealPart(x[i+j])) < min[j]) min[j] = tmp;
595: }
596: }
597: }
598: MPIU_Allreduce(min,nrm,bs,MPIU_REAL,MPIU_MIN,comm);
600: VecRestoreArrayRead(v,&x);
601: return(0);
602: }
604: /*----------------------------------------------------------------------------------------------*/
605: /*@
606: VecStrideGatherAll - Gathers all the single components from a multi-component vector into
607: separate vectors.
609: Collective on Vec
611: Input Parameter:
612: + v - the vector
613: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
615: Output Parameter:
616: . s - the location where the subvectors are stored
618: Notes:
619: One must call VecSetBlockSize() before this routine to set the stride
620: information, or use a vector created from a multicomponent DMDA.
622: If x is the array representing the vector x then this gathers
623: the arrays (x[start],x[start+stride],x[start+2*stride], ....)
624: for start=0,1,2,...bs-1
626: The parallel layout of the vector and the subvector must be the same;
627: i.e., nlocal of v = stride*(nlocal of s)
629: Not optimized; could be easily
631: Level: advanced
633: Concepts: gather^into strided vector
635: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
636: VecStrideScatterAll()
637: @*/
638: PetscErrorCode VecStrideGatherAll(Vec v,Vec s[],InsertMode addv)
639: {
640: PetscErrorCode ierr;
641: PetscInt i,n,n2,bs,j,k,*bss = NULL,nv,jj,nvc;
642: PetscScalar **y;
643: const PetscScalar *x;
649: VecGetLocalSize(v,&n);
650: VecGetLocalSize(s[0],&n2);
651: VecGetArrayRead(v,&x);
652: VecGetBlockSize(v,&bs);
653: if (bs <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Input vector does not have a valid blocksize set");
655: PetscMalloc2(bs,&y,bs,&bss);
656: nv = 0;
657: nvc = 0;
658: for (i=0; i<bs; i++) {
659: VecGetBlockSize(s[i],&bss[i]);
660: if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1 Re: [PETSC #8241] VecStrideGatherAll */
661: VecGetArray(s[i],&y[i]);
662: nvc += bss[i];
663: nv++;
664: if (nvc > bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of subvectors in subvectors > number of vectors in main vector");
665: if (nvc == bs) break;
666: }
668: n = n/bs;
670: jj = 0;
671: if (addv == INSERT_VALUES) {
672: for (j=0; j<nv; j++) {
673: for (k=0; k<bss[j]; k++) {
674: for (i=0; i<n; i++) y[j][i*bss[j] + k] = x[bs*i+jj+k];
675: }
676: jj += bss[j];
677: }
678: } else if (addv == ADD_VALUES) {
679: for (j=0; j<nv; j++) {
680: for (k=0; k<bss[j]; k++) {
681: for (i=0; i<n; i++) y[j][i*bss[j] + k] += x[bs*i+jj+k];
682: }
683: jj += bss[j];
684: }
685: #if !defined(PETSC_USE_COMPLEX)
686: } else if (addv == MAX_VALUES) {
687: for (j=0; j<nv; j++) {
688: for (k=0; k<bss[j]; k++) {
689: for (i=0; i<n; i++) y[j][i*bss[j] + k] = PetscMax(y[j][i*bss[j] + k],x[bs*i+jj+k]);
690: }
691: jj += bss[j];
692: }
693: #endif
694: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
696: VecRestoreArrayRead(v,&x);
697: for (i=0; i<nv; i++) {
698: VecRestoreArray(s[i],&y[i]);
699: }
701: PetscFree2(y,bss);
702: return(0);
703: }
705: /*@
706: VecStrideScatterAll - Scatters all the single components from separate vectors into
707: a multi-component vector.
709: Collective on Vec
711: Input Parameter:
712: + s - the location where the subvectors are stored
713: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
715: Output Parameter:
716: . v - the multicomponent vector
718: Notes:
719: One must call VecSetBlockSize() before this routine to set the stride
720: information, or use a vector created from a multicomponent DMDA.
722: The parallel layout of the vector and the subvector must be the same;
723: i.e., nlocal of v = stride*(nlocal of s)
725: Not optimized; could be easily
727: Level: advanced
729: Concepts: scatter^into strided vector
731: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
732: VecStrideScatterAll()
733: @*/
734: PetscErrorCode VecStrideScatterAll(Vec s[],Vec v,InsertMode addv)
735: {
736: PetscErrorCode ierr;
737: PetscInt i,n,n2,bs,j,jj,k,*bss = NULL,nv,nvc;
738: PetscScalar *x,**y;
744: VecGetLocalSize(v,&n);
745: VecGetLocalSize(s[0],&n2);
746: VecGetArray(v,&x);
747: VecGetBlockSize(v,&bs);
748: if (bs <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Input vector does not have a valid blocksize set");
750: PetscMalloc2(bs,&y,bs,&bss);
751: nv = 0;
752: nvc = 0;
753: for (i=0; i<bs; i++) {
754: VecGetBlockSize(s[i],&bss[i]);
755: if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1 Re: [PETSC #8241] VecStrideGatherAll */
756: VecGetArray(s[i],&y[i]);
757: nvc += bss[i];
758: nv++;
759: if (nvc > bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of subvectors in subvectors > number of vectors in main vector");
760: if (nvc == bs) break;
761: }
763: n = n/bs;
765: jj = 0;
766: if (addv == INSERT_VALUES) {
767: for (j=0; j<nv; j++) {
768: for (k=0; k<bss[j]; k++) {
769: for (i=0; i<n; i++) x[bs*i+jj+k] = y[j][i*bss[j] + k];
770: }
771: jj += bss[j];
772: }
773: } else if (addv == ADD_VALUES) {
774: for (j=0; j<nv; j++) {
775: for (k=0; k<bss[j]; k++) {
776: for (i=0; i<n; i++) x[bs*i+jj+k] += y[j][i*bss[j] + k];
777: }
778: jj += bss[j];
779: }
780: #if !defined(PETSC_USE_COMPLEX)
781: } else if (addv == MAX_VALUES) {
782: for (j=0; j<nv; j++) {
783: for (k=0; k<bss[j]; k++) {
784: for (i=0; i<n; i++) x[bs*i+jj+k] = PetscMax(x[bs*i+jj+k],y[j][i*bss[j] + k]);
785: }
786: jj += bss[j];
787: }
788: #endif
789: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
791: VecRestoreArray(v,&x);
792: for (i=0; i<nv; i++) {
793: VecRestoreArray(s[i],&y[i]);
794: }
795: PetscFree2(y,bss);
796: return(0);
797: }
799: /*@
800: VecStrideGather - Gathers a single component from a multi-component vector into
801: another vector.
803: Collective on Vec
805: Input Parameter:
806: + v - the vector
807: . start - starting point of the subvector (defined by a stride)
808: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
810: Output Parameter:
811: . s - the location where the subvector is stored
813: Notes:
814: One must call VecSetBlockSize() before this routine to set the stride
815: information, or use a vector created from a multicomponent DMDA.
817: If x is the array representing the vector x then this gathers
818: the array (x[start],x[start+stride],x[start+2*stride], ....)
820: The parallel layout of the vector and the subvector must be the same;
821: i.e., nlocal of v = stride*(nlocal of s)
823: Not optimized; could be easily
825: Level: advanced
827: Concepts: gather^into strided vector
829: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
830: VecStrideScatterAll()
831: @*/
832: PetscErrorCode VecStrideGather(Vec v,PetscInt start,Vec s,InsertMode addv)
833: {
839: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
840: 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);
841: if (!v->ops->stridegather) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
842: (*v->ops->stridegather)(v,start,s,addv);
843: return(0);
844: }
846: /*@
847: VecStrideScatter - Scatters a single component from a vector into a multi-component vector.
849: Collective on Vec
851: Input Parameter:
852: + s - the single-component vector
853: . start - starting point of the subvector (defined by a stride)
854: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
856: Output Parameter:
857: . v - the location where the subvector is scattered (the multi-component vector)
859: Notes:
860: One must call VecSetBlockSize() on the multi-component vector before this
861: routine to set the stride information, or use a vector created from a multicomponent DMDA.
863: The parallel layout of the vector and the subvector must be the same;
864: i.e., nlocal of v = stride*(nlocal of s)
866: Not optimized; could be easily
868: Level: advanced
870: Concepts: scatter^into strided vector
872: .seealso: VecStrideNorm(), VecStrideGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
873: VecStrideScatterAll(), VecStrideSubSetScatter(), VecStrideSubSetGather()
874: @*/
875: PetscErrorCode VecStrideScatter(Vec s,PetscInt start,Vec v,InsertMode addv)
876: {
882: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
883: 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);
884: if (!v->ops->stridescatter) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
885: (*v->ops->stridescatter)(s,start,v,addv);
886: return(0);
887: }
889: /*@
890: VecStrideSubSetGather - Gathers a subset of components from a multi-component vector into
891: another vector.
893: Collective on Vec
895: Input Parameter:
896: + v - the vector
897: . nidx - the number of indices
898: . idxv - the indices of the components 0 <= idxv[0] ...idxv[nidx-1] < bs(v), they need not be sorted
899: . 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
900: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
902: Output Parameter:
903: . s - the location where the subvector is stored
905: Notes:
906: One must call VecSetBlockSize() on both vectors before this routine to set the stride
907: information, or use a vector created from a multicomponent DMDA.
910: The parallel layout of the vector and the subvector must be the same;
912: Not optimized; could be easily
914: Level: advanced
916: Concepts: gather^into strided vector
918: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideGather(), VecStrideSubSetScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
919: VecStrideScatterAll()
920: @*/
921: PetscErrorCode VecStrideSubSetGather(Vec v,PetscInt nidx,const PetscInt idxv[],const PetscInt idxs[],Vec s,InsertMode addv)
922: {
928: if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
929: if (!v->ops->stridesubsetgather) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
930: (*v->ops->stridesubsetgather)(v,nidx,idxv,idxs,s,addv);
931: return(0);
932: }
934: /*@
935: VecStrideSubSetScatter - Scatters components from a vector into a subset of components of a multi-component vector.
937: Collective on Vec
939: Input Parameter:
940: + s - the smaller-component vector
941: . nidx - the number of indices in idx
942: . 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
943: . idxv - the indices of the components in the larger-component vector, 0 <= idx[0] ...idx[nidx-1] < bs(v) they need not be sorted
944: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
946: Output Parameter:
947: . v - the location where the subvector is into scattered (the multi-component vector)
949: Notes:
950: One must call VecSetBlockSize() on the vectors before this
951: routine to set the stride information, or use a vector created from a multicomponent DMDA.
953: The parallel layout of the vector and the subvector must be the same;
955: Not optimized; could be easily
957: Level: advanced
959: Concepts: scatter^into strided vector
961: .seealso: VecStrideNorm(), VecStrideGather(), VecStrideGather(), VecStrideSubSetGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
962: VecStrideScatterAll()
963: @*/
964: PetscErrorCode VecStrideSubSetScatter(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)
965: {
971: if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
972: if (!v->ops->stridesubsetscatter) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
973: (*v->ops->stridesubsetscatter)(s,nidx,idxs,idxv,v,addv);
974: return(0);
975: }
977: PetscErrorCode VecStrideGather_Default(Vec v,PetscInt start,Vec s,InsertMode addv)
978: {
980: PetscInt i,n,bs,ns;
981: const PetscScalar *x;
982: PetscScalar *y;
985: VecGetLocalSize(v,&n);
986: VecGetLocalSize(s,&ns);
987: VecGetArrayRead(v,&x);
988: VecGetArray(s,&y);
990: bs = v->map->bs;
991: if (n != ns*bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Subvector length * blocksize %D not correct for gather from original vector %D",ns*bs,n);
992: x += start;
993: n = n/bs;
995: if (addv == INSERT_VALUES) {
996: for (i=0; i<n; i++) y[i] = x[bs*i];
997: } else if (addv == ADD_VALUES) {
998: for (i=0; i<n; i++) y[i] += x[bs*i];
999: #if !defined(PETSC_USE_COMPLEX)
1000: } else if (addv == MAX_VALUES) {
1001: for (i=0; i<n; i++) y[i] = PetscMax(y[i],x[bs*i]);
1002: #endif
1003: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
1005: VecRestoreArrayRead(v,&x);
1006: VecRestoreArray(s,&y);
1007: return(0);
1008: }
1010: PetscErrorCode VecStrideScatter_Default(Vec s,PetscInt start,Vec v,InsertMode addv)
1011: {
1012: PetscErrorCode ierr;
1013: PetscInt i,n,bs,ns;
1014: PetscScalar *x;
1015: const PetscScalar *y;
1018: VecGetLocalSize(v,&n);
1019: VecGetLocalSize(s,&ns);
1020: VecGetArray(v,&x);
1021: VecGetArrayRead(s,&y);
1023: bs = v->map->bs;
1024: 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);
1025: x += start;
1026: n = n/bs;
1028: if (addv == INSERT_VALUES) {
1029: for (i=0; i<n; i++) x[bs*i] = y[i];
1030: } else if (addv == ADD_VALUES) {
1031: for (i=0; i<n; i++) x[bs*i] += y[i];
1032: #if !defined(PETSC_USE_COMPLEX)
1033: } else if (addv == MAX_VALUES) {
1034: for (i=0; i<n; i++) x[bs*i] = PetscMax(y[i],x[bs*i]);
1035: #endif
1036: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
1038: VecRestoreArray(v,&x);
1039: VecRestoreArrayRead(s,&y);
1040: return(0);
1041: }
1043: PetscErrorCode VecStrideSubSetGather_Default(Vec v,PetscInt nidx,const PetscInt idxv[],const PetscInt idxs[],Vec s,InsertMode addv)
1044: {
1045: PetscErrorCode ierr;
1046: PetscInt i,j,n,bs,bss,ns;
1047: const PetscScalar *x;
1048: PetscScalar *y;
1051: VecGetLocalSize(v,&n);
1052: VecGetLocalSize(s,&ns);
1053: VecGetArrayRead(v,&x);
1054: VecGetArray(s,&y);
1056: bs = v->map->bs;
1057: bss = s->map->bs;
1058: n = n/bs;
1060: #if defined(PETSC_DEBUG)
1061: if (n != ns/bss) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible layout of vectors");
1062: for (j=0; j<nidx; j++) {
1063: if (idxv[j] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is negative",j,idxv[j]);
1064: 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);
1065: }
1066: if (!idxs && bss != nidx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide idxs when not gathering into all locations");
1067: #endif
1069: if (addv == INSERT_VALUES) {
1070: if (!idxs) {
1071: for (i=0; i<n; i++) {
1072: for (j=0; j<bss; j++) y[bss*i+j] = x[bs*i+idxv[j]];
1073: }
1074: } else {
1075: for (i=0; i<n; i++) {
1076: for (j=0; j<bss; j++) y[bss*i+idxs[j]] = x[bs*i+idxv[j]];
1077: }
1078: }
1079: } else if (addv == ADD_VALUES) {
1080: if (!idxs) {
1081: for (i=0; i<n; i++) {
1082: for (j=0; j<bss; j++) y[bss*i+j] += x[bs*i+idxv[j]];
1083: }
1084: } else {
1085: for (i=0; i<n; i++) {
1086: for (j=0; j<bss; j++) y[bss*i+idxs[j]] += x[bs*i+idxv[j]];
1087: }
1088: }
1089: #if !defined(PETSC_USE_COMPLEX)
1090: } else if (addv == MAX_VALUES) {
1091: if (!idxs) {
1092: for (i=0; i<n; i++) {
1093: for (j=0; j<bss; j++) y[bss*i+j] = PetscMax(y[bss*i+j],x[bs*i+idxv[j]]);
1094: }
1095: } else {
1096: for (i=0; i<n; i++) {
1097: for (j=0; j<bss; j++) y[bss*i+idxs[j]] = PetscMax(y[bss*i+idxs[j]],x[bs*i+idxv[j]]);
1098: }
1099: }
1100: #endif
1101: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
1103: VecRestoreArrayRead(v,&x);
1104: VecRestoreArray(s,&y);
1105: return(0);
1106: }
1108: PetscErrorCode VecStrideSubSetScatter_Default(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)
1109: {
1110: PetscErrorCode ierr;
1111: PetscInt j,i,n,bs,ns,bss;
1112: PetscScalar *x;
1113: const PetscScalar *y;
1116: VecGetLocalSize(v,&n);
1117: VecGetLocalSize(s,&ns);
1118: VecGetArray(v,&x);
1119: VecGetArrayRead(s,&y);
1121: bs = v->map->bs;
1122: bss = s->map->bs;
1123: n = n/bs;
1125: #if defined(PETSC_DEBUG)
1126: if (n != ns/bss) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible layout of vectors");
1127: for (j=0; j<bss; j++) {
1128: if (idx[j] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is negative",j,idx[j]);
1129: 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);
1130: }
1131: if (!idxs && bss != nidx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide idxs when not scattering from all locations");
1132: #endif
1134: if (addv == INSERT_VALUES) {
1135: if (!idxs) {
1136: for (i=0; i<n; i++) {
1137: for (j=0; j<bss; j++) x[bs*i + idxv[j]] = y[bss*i+j];
1138: }
1139: } else {
1140: for (i=0; i<n; i++) {
1141: for (j=0; j<bss; j++) x[bs*i + idxv[j]] = y[bss*i+idxs[j]];
1142: }
1143: }
1144: } else if (addv == ADD_VALUES) {
1145: if (!idxs) {
1146: for (i=0; i<n; i++) {
1147: for (j=0; j<bss; j++) x[bs*i + idxv[j]] += y[bss*i+j];
1148: }
1149: } else {
1150: for (i=0; i<n; i++) {
1151: for (j=0; j<bss; j++) x[bs*i + idxv[j]] += y[bss*i+idxs[j]];
1152: }
1153: }
1154: #if !defined(PETSC_USE_COMPLEX)
1155: } else if (addv == MAX_VALUES) {
1156: if (!idxs) {
1157: for (i=0; i<n; i++) {
1158: for (j=0; j<bss; j++) x[bs*i + idxv[j]] = PetscMax(y[bss*i+j],x[bs*i + idxv[j]]);
1159: }
1160: } else {
1161: for (i=0; i<n; i++) {
1162: for (j=0; j<bss; j++) x[bs*i + idxv[j]] = PetscMax(y[bss*i+idxs[j]],x[bs*i + idxv[j]]);
1163: }
1164: }
1165: #endif
1166: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
1168: VecRestoreArray(v,&x);
1169: VecRestoreArrayRead(s,&y);
1170: return(0);
1171: }
1173: PetscErrorCode VecReciprocal_Default(Vec v)
1174: {
1176: PetscInt i,n;
1177: PetscScalar *x;
1180: VecGetLocalSize(v,&n);
1181: VecGetArray(v,&x);
1182: for (i=0; i<n; i++) {
1183: if (x[i] != (PetscScalar)0.0) x[i] = (PetscScalar)1.0/x[i];
1184: }
1185: VecRestoreArray(v,&x);
1186: return(0);
1187: }
1189: /*@
1190: VecExp - Replaces each component of a vector by e^x_i
1192: Not collective
1194: Input Parameter:
1195: . v - The vector
1197: Output Parameter:
1198: . v - The vector of exponents
1200: Level: beginner
1202: .seealso: VecLog(), VecAbs(), VecSqrtAbs(), VecReciprocal()
1204: .keywords: vector, sqrt, square root
1205: @*/
1206: PetscErrorCode VecExp(Vec v)
1207: {
1208: PetscScalar *x;
1209: PetscInt i, n;
1214: if (v->ops->exp) {
1215: (*v->ops->exp)(v);
1216: } else {
1217: VecGetLocalSize(v, &n);
1218: VecGetArray(v, &x);
1219: for (i = 0; i < n; i++) x[i] = PetscExpScalar(x[i]);
1220: VecRestoreArray(v, &x);
1221: }
1222: return(0);
1223: }
1225: /*@
1226: VecLog - Replaces each component of a vector by log(x_i), the natural logarithm
1228: Not collective
1230: Input Parameter:
1231: . v - The vector
1233: Output Parameter:
1234: . v - The vector of logs
1236: Level: beginner
1238: .seealso: VecExp(), VecAbs(), VecSqrtAbs(), VecReciprocal()
1240: .keywords: vector, sqrt, square root
1241: @*/
1242: PetscErrorCode VecLog(Vec v)
1243: {
1244: PetscScalar *x;
1245: PetscInt i, n;
1250: if (v->ops->log) {
1251: (*v->ops->log)(v);
1252: } else {
1253: VecGetLocalSize(v, &n);
1254: VecGetArray(v, &x);
1255: for (i = 0; i < n; i++) x[i] = PetscLogScalar(x[i]);
1256: VecRestoreArray(v, &x);
1257: }
1258: return(0);
1259: }
1261: /*@
1262: VecSqrtAbs - Replaces each component of a vector by the square root of its magnitude.
1264: Not collective
1266: Input Parameter:
1267: . v - The vector
1269: Output Parameter:
1270: . v - The vector square root
1272: Level: beginner
1274: Note: The actual function is sqrt(|x_i|)
1276: .seealso: VecLog(), VecExp(), VecReciprocal(), VecAbs()
1278: .keywords: vector, sqrt, square root
1279: @*/
1280: PetscErrorCode VecSqrtAbs(Vec v)
1281: {
1282: PetscScalar *x;
1283: PetscInt i, n;
1288: if (v->ops->sqrt) {
1289: (*v->ops->sqrt)(v);
1290: } else {
1291: VecGetLocalSize(v, &n);
1292: VecGetArray(v, &x);
1293: for (i = 0; i < n; i++) x[i] = PetscSqrtReal(PetscAbsScalar(x[i]));
1294: VecRestoreArray(v, &x);
1295: }
1296: return(0);
1297: }
1299: /*@
1300: VecDotNorm2 - computes the inner product of two vectors and the 2-norm squared of the second vector
1302: Collective on Vec
1304: Input Parameter:
1305: + s - first vector
1306: - t - second vector
1308: Output Parameter:
1309: + dp - s'conj(t)
1310: - nm - t'conj(t)
1312: Level: advanced
1314: Notes:
1315: conj(x) is the complex conjugate of x when x is complex
1318: .seealso: VecDot(), VecNorm(), VecDotBegin(), VecNormBegin(), VecDotEnd(), VecNormEnd()
1320: .keywords: vector, sqrt, square root
1321: @*/
1322: PetscErrorCode VecDotNorm2(Vec s,Vec t,PetscScalar *dp, PetscReal *nm)
1323: {
1324: const PetscScalar *sx, *tx;
1325: PetscScalar dpx = 0.0, nmx = 0.0,work[2],sum[2];
1326: PetscInt i, n;
1327: PetscErrorCode ierr;
1337: if (s->map->N != t->map->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
1338: if (s->map->n != t->map->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");
1340: PetscLogEventBegin(VEC_DotNorm2,s,t,0,0);
1341: if (s->ops->dotnorm2) {
1342: (*s->ops->dotnorm2)(s,t,dp,&dpx);
1343: *nm = PetscRealPart(dpx);
1344: } else {
1345: VecGetLocalSize(s, &n);
1346: VecGetArrayRead(s, &sx);
1347: VecGetArrayRead(t, &tx);
1349: for (i = 0; i<n; i++) {
1350: dpx += sx[i]*PetscConj(tx[i]);
1351: nmx += tx[i]*PetscConj(tx[i]);
1352: }
1353: work[0] = dpx;
1354: work[1] = nmx;
1356: MPIU_Allreduce(work,sum,2,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)s));
1357: *dp = sum[0];
1358: *nm = PetscRealPart(sum[1]);
1360: VecRestoreArrayRead(t, &tx);
1361: VecRestoreArrayRead(s, &sx);
1362: PetscLogFlops(4.0*n);
1363: }
1364: PetscLogEventEnd(VEC_DotNorm2,s,t,0,0);
1365: return(0);
1366: }
1368: /*@
1369: VecSum - Computes the sum of all the components of a vector.
1371: Collective on Vec
1373: Input Parameter:
1374: . v - the vector
1376: Output Parameter:
1377: . sum - the result
1379: Level: beginner
1381: Concepts: sum^of vector entries
1383: .seealso: VecNorm()
1384: @*/
1385: PetscErrorCode VecSum(Vec v,PetscScalar *sum)
1386: {
1387: PetscErrorCode ierr;
1388: PetscInt i,n;
1389: const PetscScalar *x;
1390: PetscScalar lsum = 0.0;
1395: VecGetLocalSize(v,&n);
1396: VecGetArrayRead(v,&x);
1397: for (i=0; i<n; i++) lsum += x[i];
1398: MPIU_Allreduce(&lsum,sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)v));
1399: VecRestoreArrayRead(v,&x);
1400: return(0);
1401: }
1403: /*@
1404: VecShift - Shifts all of the components of a vector by computing
1405: x[i] = x[i] + shift.
1407: Logically Collective on Vec
1409: Input Parameters:
1410: + v - the vector
1411: - shift - the shift
1413: Output Parameter:
1414: . v - the shifted vector
1416: Level: intermediate
1418: Concepts: vector^adding constant
1420: @*/
1421: PetscErrorCode VecShift(Vec v,PetscScalar shift)
1422: {
1424: PetscInt i,n;
1425: PetscScalar *x;
1430: VecLocked(v,1);
1432: if (v->ops->shift) {
1433: (*v->ops->shift)(v,shift);
1434: } else {
1435: VecGetLocalSize(v,&n);
1436: VecGetArray(v,&x);
1437: for (i=0; i<n; i++) x[i] += shift;
1438: VecRestoreArray(v,&x);
1439: }
1440: return(0);
1441: }
1443: /*@
1444: VecAbs - Replaces every element in a vector with its absolute value.
1446: Logically Collective on Vec
1448: Input Parameters:
1449: . v - the vector
1451: Level: intermediate
1453: Concepts: vector^absolute value
1455: @*/
1456: PetscErrorCode VecAbs(Vec v)
1457: {
1459: PetscInt i,n;
1460: PetscScalar *x;
1464: VecLocked(v,1);
1465:
1466: if (v->ops->abs) {
1467: (*v->ops->abs)(v);
1468: } else {
1469: VecGetLocalSize(v,&n);
1470: VecGetArray(v,&x);
1471: for (i=0; i<n; i++) x[i] = PetscAbsScalar(x[i]);
1472: VecRestoreArray(v,&x);
1473: }
1474: return(0);
1475: }
1477: /*@
1478: VecPermute - Permutes a vector in place using the given ordering.
1480: Input Parameters:
1481: + vec - The vector
1482: . order - The ordering
1483: - inv - The flag for inverting the permutation
1485: Level: beginner
1487: Note: This function does not yet support parallel Index Sets with non-local permutations
1489: .seealso: MatPermute()
1490: .keywords: vec, permute
1491: @*/
1492: PetscErrorCode VecPermute(Vec x, IS row, PetscBool inv)
1493: {
1494: PetscScalar *array, *newArray;
1495: const PetscInt *idx;
1496: PetscInt i,rstart,rend;
1500: VecLocked(x,1);
1501:
1502: VecGetOwnershipRange(x,&rstart,&rend);
1503: ISGetIndices(row, &idx);
1504: VecGetArray(x, &array);
1505: PetscMalloc1(x->map->n, &newArray);
1506: #if defined(PETSC_USE_DEBUG)
1507: for (i = 0; i < x->map->n; i++) {
1508: 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]);
1509: }
1510: #endif
1511: if (!inv) {
1512: for (i = 0; i < x->map->n; i++) newArray[i] = array[idx[i]-rstart];
1513: } else {
1514: for (i = 0; i < x->map->n; i++) newArray[idx[i]-rstart] = array[i];
1515: }
1516: VecRestoreArray(x, &array);
1517: ISRestoreIndices(row, &idx);
1518: VecReplaceArray(x, newArray);
1519: return(0);
1520: }
1522: /*@
1523: VecEqual - Compares two vectors. Returns true if the two vectors are either pointing to the same memory buffer,
1524: or if the two vectors have the same local and global layout as well as bitwise equality of all entries.
1525: Does NOT take round-off errors into account.
1527: Collective on Vec
1529: Input Parameters:
1530: + vec1 - the first vector
1531: - vec2 - the second vector
1533: Output Parameter:
1534: . flg - PETSC_TRUE if the vectors are equal; PETSC_FALSE otherwise.
1536: Level: intermediate
1538: Concepts: equal^two vectors
1539: Concepts: vector^equality
1541: @*/
1542: PetscErrorCode VecEqual(Vec vec1,Vec vec2,PetscBool *flg)
1543: {
1544: const PetscScalar *v1,*v2;
1545: PetscErrorCode ierr;
1546: PetscInt n1,n2,N1,N2;
1547: PetscBool flg1;
1553: if (vec1 == vec2) *flg = PETSC_TRUE;
1554: else {
1555: VecGetSize(vec1,&N1);
1556: VecGetSize(vec2,&N2);
1557: if (N1 != N2) flg1 = PETSC_FALSE;
1558: else {
1559: VecGetLocalSize(vec1,&n1);
1560: VecGetLocalSize(vec2,&n2);
1561: if (n1 != n2) flg1 = PETSC_FALSE;
1562: else {
1563: VecGetArrayRead(vec1,&v1);
1564: VecGetArrayRead(vec2,&v2);
1565: PetscMemcmp(v1,v2,n1*sizeof(PetscScalar),&flg1);
1566: VecRestoreArrayRead(vec1,&v1);
1567: VecRestoreArrayRead(vec2,&v2);
1568: }
1569: }
1570: /* combine results from all processors */
1571: MPIU_Allreduce(&flg1,flg,1,MPIU_BOOL,MPI_MIN,PetscObjectComm((PetscObject)vec1));
1572: }
1573: return(0);
1574: }
1576: /*@
1577: VecUniqueEntries - Compute the number of unique entries, and those entries
1579: Collective on Vec
1581: Input Parameter:
1582: . vec - the vector
1584: Output Parameters:
1585: + n - The number of unique entries
1586: - e - The entries
1588: Level: intermediate
1590: @*/
1591: PetscErrorCode VecUniqueEntries(Vec vec, PetscInt *n, PetscScalar **e)
1592: {
1593: const PetscScalar *v;
1594: PetscScalar *tmp, *vals;
1595: PetscMPIInt *N, *displs, l;
1596: PetscInt ng, m, i, j, p;
1597: PetscMPIInt size;
1598: PetscErrorCode ierr;
1603: MPI_Comm_size(PetscObjectComm((PetscObject) vec), &size);
1604: VecGetLocalSize(vec, &m);
1605: VecGetArrayRead(vec, &v);
1606: PetscMalloc2(m,&tmp,size,&N);
1607: for (i = 0, j = 0, l = 0; i < m; ++i) {
1608: /* Can speed this up with sorting */
1609: for (j = 0; j < l; ++j) {
1610: if (v[i] == tmp[j]) break;
1611: }
1612: if (j == l) {
1613: tmp[j] = v[i];
1614: ++l;
1615: }
1616: }
1617: VecRestoreArrayRead(vec, &v);
1618: /* Gather serial results */
1619: MPI_Allgather(&l, 1, MPI_INT, N, 1, MPI_INT, PetscObjectComm((PetscObject) vec));
1620: for (p = 0, ng = 0; p < size; ++p) {
1621: ng += N[p];
1622: }
1623: PetscMalloc2(ng,&vals,size+1,&displs);
1624: for (p = 1, displs[0] = 0; p <= size; ++p) {
1625: displs[p] = displs[p-1] + N[p-1];
1626: }
1627: MPI_Allgatherv(tmp, l, MPIU_SCALAR, vals, N, displs, MPIU_SCALAR, PetscObjectComm((PetscObject) vec));
1628: /* Find unique entries */
1629: #ifdef PETSC_USE_COMPLEX
1630: SETERRQ(PetscObjectComm((PetscObject) vec), PETSC_ERR_SUP, "Does not work with complex numbers");
1631: #else
1632: *n = displs[size];
1633: PetscSortRemoveDupsReal(n, (PetscReal *) vals);
1634: if (e) {
1636: PetscMalloc1(*n, e);
1637: for (i = 0; i < *n; ++i) {
1638: (*e)[i] = vals[i];
1639: }
1640: }
1641: PetscFree2(vals,displs);
1642: PetscFree2(tmp,N);
1643: return(0);
1644: #endif
1645: }