Actual source code: vinv.c
2: /*
3: Some useful vector utility functions.
4: */
5: #include <../src/vec/vec/impls/mpi/pvecimpl.h>
7: /*@
8: VecStrideSet - Sets a subvector of a vector defined
9: by a starting point and a stride with a given value
11: Logically Collective on Vec
13: Input Parameters:
14: + v - the vector
15: . start - starting point of the subvector (defined by a stride)
16: - s - value to set for each entry in that subvector
18: Notes:
19: One must call VecSetBlockSize() before this routine to set the stride
20: information, or use a vector created from a multicomponent DMDA.
22: This will only work if the desire subvector is a stride subvector
24: Level: advanced
26: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideScale()
27: @*/
28: PetscErrorCode VecStrideSet(Vec v,PetscInt start,PetscScalar s)
29: {
31: PetscInt i,n,bs;
32: PetscScalar *x;
38: VecSetErrorIfLocked(v,1);
40: VecGetLocalSize(v,&n);
41: VecGetArray(v,&x);
42: VecGetBlockSize(v,&bs);
43: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
44: 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);
45: x += start;
46: for (i=0; i<n; i+=bs) x[i] = s;
47: x -= start;
48: VecRestoreArray(v,&x);
49: return(0);
50: }
52: /*@
53: VecStrideScale - Scales a subvector of a vector defined
54: by a starting point and a stride.
56: Logically Collective on Vec
58: Input Parameters:
59: + v - the vector
60: . start - starting point of the subvector (defined by a stride)
61: - scale - value to multiply each subvector entry by
63: Notes:
64: One must call VecSetBlockSize() before this routine to set the stride
65: information, or use a vector created from a multicomponent DMDA.
67: This will only work if the desire subvector is a stride subvector
69: Level: advanced
71: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideScale()
72: @*/
73: PetscErrorCode VecStrideScale(Vec v,PetscInt start,PetscScalar scale)
74: {
76: PetscInt i,n,bs;
77: PetscScalar *x;
83: VecSetErrorIfLocked(v,1);
85: VecGetLocalSize(v,&n);
86: VecGetArray(v,&x);
87: VecGetBlockSize(v,&bs);
88: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
89: 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);
90: x += start;
91: for (i=0; i<n; i+=bs) x[i] *= scale;
92: x -= start;
93: VecRestoreArray(v,&x);
94: return(0);
95: }
97: /*@
98: VecStrideNorm - Computes the norm of subvector of a vector defined
99: by a starting point and a stride.
101: Collective on Vec
103: Input Parameters:
104: + v - the vector
105: . start - starting point of the subvector (defined by a stride)
106: - ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY
108: Output Parameter:
109: . norm - the norm
111: Notes:
112: One must call VecSetBlockSize() before this routine to set the stride
113: information, or use a vector created from a multicomponent DMDA.
115: If x is the array representing the vector x then this computes the norm
116: of the array (x[start],x[start+stride],x[start+2*stride], ....)
118: This is useful for computing, say the norm of the pressure variable when
119: the pressure is stored (interlaced) with other variables, say density etc.
121: This will only work if the desire subvector is a stride subvector
123: Level: advanced
125: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
126: @*/
127: PetscErrorCode VecStrideNorm(Vec v,PetscInt start,NormType ntype,PetscReal *nrm)
128: {
129: PetscErrorCode ierr;
130: PetscInt i,n,bs;
131: const PetscScalar *x;
132: PetscReal tnorm;
133: MPI_Comm comm;
138: VecGetLocalSize(v,&n);
139: VecGetArrayRead(v,&x);
140: PetscObjectGetComm((PetscObject)v,&comm);
142: VecGetBlockSize(v,&bs);
143: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
144: 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);
145: x += start;
147: if (ntype == NORM_2) {
148: PetscScalar sum = 0.0;
149: for (i=0; i<n; i+=bs) sum += x[i]*(PetscConj(x[i]));
150: tnorm = PetscRealPart(sum);
151: MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_SUM,comm);
152: *nrm = PetscSqrtReal(*nrm);
153: } else if (ntype == NORM_1) {
154: tnorm = 0.0;
155: for (i=0; i<n; i+=bs) tnorm += PetscAbsScalar(x[i]);
156: MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_SUM,comm);
157: } else if (ntype == NORM_INFINITY) {
158: PetscReal tmp;
159: tnorm = 0.0;
161: for (i=0; i<n; i+=bs) {
162: if ((tmp = PetscAbsScalar(x[i])) > tnorm) tnorm = tmp;
163: /* check special case of tmp == NaN */
164: if (tmp != tmp) {tnorm = tmp; break;}
165: }
166: MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_MAX,comm);
167: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
168: VecRestoreArrayRead(v,&x);
169: return(0);
170: }
172: /*@
173: VecStrideMax - Computes the maximum of subvector of a vector defined
174: by a starting point and a stride and optionally its location.
176: Collective on Vec
178: Input Parameters:
179: + v - the vector
180: - start - starting point of the subvector (defined by a stride)
182: Output Parameters:
183: + idex - the location where the maximum occurred (pass NULL if not required)
184: - nrm - the maximum value in the subvector
186: Notes:
187: One must call VecSetBlockSize() before this routine to set the stride
188: information, or use a vector created from a multicomponent DMDA.
190: If xa is the array representing the vector x, then this computes the max
191: of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)
193: This is useful for computing, say the maximum of the pressure variable when
194: the pressure is stored (interlaced) with other variables, e.g., density, etc.
195: This will only work if the desire subvector is a stride subvector.
197: Level: advanced
199: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
200: @*/
201: PetscErrorCode VecStrideMax(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
202: {
203: PetscErrorCode ierr;
204: PetscInt i,n,bs,id;
205: const PetscScalar *x;
206: PetscReal max,tmp;
212: VecGetLocalSize(v,&n);
213: VecGetArrayRead(v,&x);
215: VecGetBlockSize(v,&bs);
216: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
217: 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);
218: x += start;
220: id = -1;
221: if (!n) max = PETSC_MIN_REAL;
222: else {
223: id = 0;
224: max = PetscRealPart(x[0]);
225: for (i=bs; i<n; i+=bs) {
226: if ((tmp = PetscRealPart(x[i])) > max) { max = tmp; id = i;}
227: }
228: }
229: VecRestoreArrayRead(v,&x);
231: if (!idex) {
232: MPIU_Allreduce(&max,nrm,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)v));
233: } else {
234: struct { PetscReal v; PetscInt i; } in,out;
235: PetscInt rstart;
237: VecGetOwnershipRange(v,&rstart,NULL);
238: in.v = max;
239: in.i = rstart+id+start;
240: MPIU_Allreduce(&in,&out,1,MPIU_REAL_INT,MPIU_MAXLOC,PetscObjectComm((PetscObject)v));
241: *nrm = out.v;
242: *idex = out.i;
243: }
244: return(0);
245: }
247: /*@
248: VecStrideMin - Computes the minimum of subvector of a vector defined
249: by a starting point and a stride and optionally its location.
251: Collective on Vec
253: Input Parameters:
254: + v - the vector
255: - start - starting point of the subvector (defined by a stride)
257: Output Parameters:
258: + idex - the location where the minimum occurred. (pass NULL if not required)
259: - nrm - the minimum value in the subvector
261: Level: advanced
263: Notes:
264: One must call VecSetBlockSize() before this routine to set the stride
265: information, or use a vector created from a multicomponent DMDA.
267: If xa is the array representing the vector x, then this computes the min
268: of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)
270: This is useful for computing, say the minimum of the pressure variable when
271: the pressure is stored (interlaced) with other variables, e.g., density, etc.
272: This will only work if the desire subvector is a stride subvector.
274: .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
275: @*/
276: PetscErrorCode VecStrideMin(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
277: {
278: PetscErrorCode ierr;
279: PetscInt i,n,bs,id;
280: const PetscScalar *x;
281: PetscReal min,tmp;
282: MPI_Comm comm;
288: VecGetLocalSize(v,&n);
289: VecGetArrayRead(v,&x);
290: PetscObjectGetComm((PetscObject)v,&comm);
292: VecGetBlockSize(v,&bs);
293: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
294: 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);
295: x += start;
297: id = -1;
298: if (!n) min = PETSC_MAX_REAL;
299: else {
300: id = 0;
301: min = PetscRealPart(x[0]);
302: for (i=bs; i<n; i+=bs) {
303: if ((tmp = PetscRealPart(x[i])) < min) { min = tmp; id = i;}
304: }
305: }
306: VecRestoreArrayRead(v,&x);
308: if (!idex) {
309: MPIU_Allreduce(&min,nrm,1,MPIU_REAL,MPIU_MIN,comm);
310: } else {
311: struct { PetscReal v; PetscInt i; } in,out;
312: PetscInt rstart;
314: VecGetOwnershipRange(v,&rstart,NULL);
315: in.v = min;
316: in.i = rstart+id+start;
317: MPIU_Allreduce(&in,&out,1,MPIU_REAL_INT,MPIU_MINLOC,PetscObjectComm((PetscObject)v));
318: *nrm = out.v;
319: *idex = out.i;
320: }
321: return(0);
322: }
324: /*@
325: VecStrideScaleAll - Scales the subvectors of a vector defined
326: by a starting point and a stride.
328: Logically Collective on Vec
330: Input Parameters:
331: + v - the vector
332: - scales - values to multiply each subvector entry by
334: Notes:
335: One must call VecSetBlockSize() before this routine to set the stride
336: information, or use a vector created from a multicomponent DMDA.
338: The dimension of scales must be the same as the vector block size
340: Level: advanced
342: .seealso: VecNorm(), VecStrideScale(), VecScale(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
343: @*/
344: PetscErrorCode VecStrideScaleAll(Vec v,const PetscScalar *scales)
345: {
347: PetscInt i,j,n,bs;
348: PetscScalar *x;
353: VecSetErrorIfLocked(v,1);
355: VecGetLocalSize(v,&n);
356: VecGetArray(v,&x);
357: VecGetBlockSize(v,&bs);
359: /* need to provide optimized code for each bs */
360: for (i=0; i<n; i+=bs) {
361: for (j=0; j<bs; j++) x[i+j] *= scales[j];
362: }
363: VecRestoreArray(v,&x);
364: return(0);
365: }
367: /*@
368: VecStrideNormAll - Computes the norms of subvectors of a vector defined
369: by a starting point and a stride.
371: Collective on Vec
373: Input Parameters:
374: + v - the vector
375: - ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY
377: Output Parameter:
378: . nrm - the norms
380: Notes:
381: One must call VecSetBlockSize() before this routine to set the stride
382: information, or use a vector created from a multicomponent DMDA.
384: If x is the array representing the vector x then this computes the norm
385: of the array (x[start],x[start+stride],x[start+2*stride], ....) for each start < stride
387: The dimension of nrm must be the same as the vector block size
389: This will only work if the desire subvector is a stride subvector
391: Level: advanced
393: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
394: @*/
395: PetscErrorCode VecStrideNormAll(Vec v,NormType ntype,PetscReal nrm[])
396: {
397: PetscErrorCode ierr;
398: PetscInt i,j,n,bs;
399: const PetscScalar *x;
400: PetscReal tnorm[128];
401: MPI_Comm comm;
406: VecGetLocalSize(v,&n);
407: VecGetArrayRead(v,&x);
408: PetscObjectGetComm((PetscObject)v,&comm);
410: VecGetBlockSize(v,&bs);
411: if (bs > 128) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");
413: if (ntype == NORM_2) {
414: PetscScalar sum[128];
415: for (j=0; j<bs; j++) sum[j] = 0.0;
416: for (i=0; i<n; i+=bs) {
417: for (j=0; j<bs; j++) sum[j] += x[i+j]*(PetscConj(x[i+j]));
418: }
419: for (j=0; j<bs; j++) tnorm[j] = PetscRealPart(sum[j]);
421: MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_SUM,comm);
422: for (j=0; j<bs; j++) nrm[j] = PetscSqrtReal(nrm[j]);
423: } else if (ntype == NORM_1) {
424: for (j=0; j<bs; j++) tnorm[j] = 0.0;
426: for (i=0; i<n; i+=bs) {
427: for (j=0; j<bs; j++) tnorm[j] += PetscAbsScalar(x[i+j]);
428: }
430: MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_SUM,comm);
431: } else if (ntype == NORM_INFINITY) {
432: PetscReal tmp;
433: for (j=0; j<bs; j++) tnorm[j] = 0.0;
435: for (i=0; i<n; i+=bs) {
436: for (j=0; j<bs; j++) {
437: if ((tmp = PetscAbsScalar(x[i+j])) > tnorm[j]) tnorm[j] = tmp;
438: /* check special case of tmp == NaN */
439: if (tmp != tmp) {tnorm[j] = tmp; break;}
440: }
441: }
442: MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_MAX,comm);
443: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
444: VecRestoreArrayRead(v,&x);
445: return(0);
446: }
448: /*@
449: VecStrideMaxAll - Computes the maximums of subvectors of a vector defined
450: by a starting point and a stride and optionally its location.
452: Collective on Vec
454: Input Parameter:
455: . v - the vector
457: Output Parameters:
458: + index - the location where the maximum occurred (not supported, pass NULL,
459: if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
460: - nrm - the maximum values of each subvector
462: Notes:
463: One must call VecSetBlockSize() before this routine to set the stride
464: information, or use a vector created from a multicomponent DMDA.
466: The dimension of nrm must be the same as the vector block size
468: Level: advanced
470: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
471: @*/
472: PetscErrorCode VecStrideMaxAll(Vec v,PetscInt idex[],PetscReal nrm[])
473: {
474: PetscErrorCode ierr;
475: PetscInt i,j,n,bs;
476: const PetscScalar *x;
477: PetscReal max[128],tmp;
478: MPI_Comm comm;
483: 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");
484: VecGetLocalSize(v,&n);
485: VecGetArrayRead(v,&x);
486: PetscObjectGetComm((PetscObject)v,&comm);
488: VecGetBlockSize(v,&bs);
489: if (bs > 128) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");
491: if (!n) {
492: for (j=0; j<bs; j++) max[j] = PETSC_MIN_REAL;
493: } else {
494: for (j=0; j<bs; j++) max[j] = PetscRealPart(x[j]);
496: for (i=bs; i<n; i+=bs) {
497: for (j=0; j<bs; j++) {
498: if ((tmp = PetscRealPart(x[i+j])) > max[j]) max[j] = tmp;
499: }
500: }
501: }
502: MPIU_Allreduce(max,nrm,bs,MPIU_REAL,MPIU_MAX,comm);
504: VecRestoreArrayRead(v,&x);
505: return(0);
506: }
508: /*@
509: VecStrideMinAll - Computes the minimum of subvector of a vector defined
510: by a starting point and a stride and optionally its location.
512: Collective on Vec
514: Input Parameter:
515: . v - the vector
517: Output Parameters:
518: + idex - the location where the minimum occurred (not supported, pass NULL,
519: if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
520: - nrm - the minimums of each subvector
522: Level: advanced
524: Notes:
525: One must call VecSetBlockSize() before this routine to set the stride
526: information, or use a vector created from a multicomponent DMDA.
528: The dimension of nrm must be the same as the vector block size
530: .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
531: @*/
532: PetscErrorCode VecStrideMinAll(Vec v,PetscInt idex[],PetscReal nrm[])
533: {
534: PetscErrorCode ierr;
535: PetscInt i,n,bs,j;
536: const PetscScalar *x;
537: PetscReal min[128],tmp;
538: MPI_Comm comm;
543: 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");
544: VecGetLocalSize(v,&n);
545: VecGetArrayRead(v,&x);
546: PetscObjectGetComm((PetscObject)v,&comm);
548: VecGetBlockSize(v,&bs);
549: if (bs > 128) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");
551: if (!n) {
552: for (j=0; j<bs; j++) min[j] = PETSC_MAX_REAL;
553: } else {
554: for (j=0; j<bs; j++) min[j] = PetscRealPart(x[j]);
556: for (i=bs; i<n; i+=bs) {
557: for (j=0; j<bs; j++) {
558: if ((tmp = PetscRealPart(x[i+j])) < min[j]) min[j] = tmp;
559: }
560: }
561: }
562: MPIU_Allreduce(min,nrm,bs,MPIU_REAL,MPIU_MIN,comm);
564: VecRestoreArrayRead(v,&x);
565: return(0);
566: }
568: /*----------------------------------------------------------------------------------------------*/
569: /*@
570: VecStrideGatherAll - Gathers all the single components from a multi-component vector into
571: separate vectors.
573: Collective on Vec
575: Input Parameters:
576: + v - the vector
577: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
579: Output Parameter:
580: . s - the location where the subvectors are stored
582: Notes:
583: One must call VecSetBlockSize() before this routine to set the stride
584: information, or use a vector created from a multicomponent DMDA.
586: If x is the array representing the vector x then this gathers
587: the arrays (x[start],x[start+stride],x[start+2*stride], ....)
588: for start=0,1,2,...bs-1
590: The parallel layout of the vector and the subvector must be the same;
591: i.e., nlocal of v = stride*(nlocal of s)
593: Not optimized; could be easily
595: Level: advanced
597: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
598: VecStrideScatterAll()
599: @*/
600: PetscErrorCode VecStrideGatherAll(Vec v,Vec s[],InsertMode addv)
601: {
602: PetscErrorCode ierr;
603: PetscInt i,n,n2,bs,j,k,*bss = NULL,nv,jj,nvc;
604: PetscScalar **y;
605: const PetscScalar *x;
611: VecGetLocalSize(v,&n);
612: VecGetLocalSize(s[0],&n2);
613: VecGetArrayRead(v,&x);
614: VecGetBlockSize(v,&bs);
615: if (bs <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Input vector does not have a valid blocksize set");
617: PetscMalloc2(bs,&y,bs,&bss);
618: nv = 0;
619: nvc = 0;
620: for (i=0; i<bs; i++) {
621: VecGetBlockSize(s[i],&bss[i]);
622: if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1 Re: [PETSC #8241] VecStrideGatherAll */
623: VecGetArray(s[i],&y[i]);
624: nvc += bss[i];
625: nv++;
626: if (nvc > bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of subvectors in subvectors > number of vectors in main vector");
627: if (nvc == bs) break;
628: }
630: n = n/bs;
632: jj = 0;
633: if (addv == INSERT_VALUES) {
634: for (j=0; j<nv; j++) {
635: for (k=0; k<bss[j]; k++) {
636: for (i=0; i<n; i++) y[j][i*bss[j] + k] = x[bs*i+jj+k];
637: }
638: jj += bss[j];
639: }
640: } else if (addv == ADD_VALUES) {
641: for (j=0; j<nv; j++) {
642: for (k=0; k<bss[j]; k++) {
643: for (i=0; i<n; i++) y[j][i*bss[j] + k] += x[bs*i+jj+k];
644: }
645: jj += bss[j];
646: }
647: #if !defined(PETSC_USE_COMPLEX)
648: } else if (addv == MAX_VALUES) {
649: for (j=0; j<nv; j++) {
650: for (k=0; k<bss[j]; k++) {
651: for (i=0; i<n; i++) y[j][i*bss[j] + k] = PetscMax(y[j][i*bss[j] + k],x[bs*i+jj+k]);
652: }
653: jj += bss[j];
654: }
655: #endif
656: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
658: VecRestoreArrayRead(v,&x);
659: for (i=0; i<nv; i++) {
660: VecRestoreArray(s[i],&y[i]);
661: }
663: PetscFree2(y,bss);
664: return(0);
665: }
667: /*@
668: VecStrideScatterAll - Scatters all the single components from separate vectors into
669: a multi-component vector.
671: Collective on Vec
673: Input Parameters:
674: + s - the location where the subvectors are stored
675: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
677: Output Parameter:
678: . v - the multicomponent vector
680: Notes:
681: One must call VecSetBlockSize() before this routine to set the stride
682: information, or use a vector created from a multicomponent DMDA.
684: The parallel layout of the vector and the subvector must be the same;
685: i.e., nlocal of v = stride*(nlocal of s)
687: Not optimized; could be easily
689: Level: advanced
691: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
692: VecStrideScatterAll()
693: @*/
694: PetscErrorCode VecStrideScatterAll(Vec s[],Vec v,InsertMode addv)
695: {
696: PetscErrorCode ierr;
697: PetscInt i,n,n2,bs,j,jj,k,*bss = NULL,nv,nvc;
698: PetscScalar *x;
699: PetscScalar const **y;
705: VecGetLocalSize(v,&n);
706: VecGetLocalSize(s[0],&n2);
707: VecGetArray(v,&x);
708: VecGetBlockSize(v,&bs);
709: if (bs <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Input vector does not have a valid blocksize set");
711: PetscMalloc2(bs,(PetscScalar***)&y,bs,&bss);
712: nv = 0;
713: nvc = 0;
714: for (i=0; i<bs; i++) {
715: VecGetBlockSize(s[i],&bss[i]);
716: if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1 Re: [PETSC #8241] VecStrideGatherAll */
717: VecGetArrayRead(s[i],&y[i]);
718: nvc += bss[i];
719: nv++;
720: if (nvc > bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of subvectors in subvectors > number of vectors in main vector");
721: if (nvc == bs) break;
722: }
724: n = n/bs;
726: jj = 0;
727: if (addv == INSERT_VALUES) {
728: for (j=0; j<nv; j++) {
729: for (k=0; k<bss[j]; k++) {
730: for (i=0; i<n; i++) x[bs*i+jj+k] = y[j][i*bss[j] + k];
731: }
732: jj += bss[j];
733: }
734: } else if (addv == ADD_VALUES) {
735: for (j=0; j<nv; j++) {
736: for (k=0; k<bss[j]; k++) {
737: for (i=0; i<n; i++) x[bs*i+jj+k] += y[j][i*bss[j] + k];
738: }
739: jj += bss[j];
740: }
741: #if !defined(PETSC_USE_COMPLEX)
742: } else if (addv == MAX_VALUES) {
743: for (j=0; j<nv; j++) {
744: for (k=0; k<bss[j]; k++) {
745: for (i=0; i<n; i++) x[bs*i+jj+k] = PetscMax(x[bs*i+jj+k],y[j][i*bss[j] + k]);
746: }
747: jj += bss[j];
748: }
749: #endif
750: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
752: VecRestoreArray(v,&x);
753: for (i=0; i<nv; i++) {
754: VecRestoreArrayRead(s[i],&y[i]);
755: }
756: PetscFree2(*(PetscScalar***)&y,bss);
757: return(0);
758: }
760: /*@
761: VecStrideGather - Gathers a single component from a multi-component vector into
762: another vector.
764: Collective on Vec
766: Input Parameters:
767: + v - the vector
768: . start - starting point of the subvector (defined by a stride)
769: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
771: Output Parameter:
772: . s - the location where the subvector is stored
774: Notes:
775: One must call VecSetBlockSize() before this routine to set the stride
776: information, or use a vector created from a multicomponent DMDA.
778: If x is the array representing the vector x then this gathers
779: the array (x[start],x[start+stride],x[start+2*stride], ....)
781: The parallel layout of the vector and the subvector must be the same;
782: i.e., nlocal of v = stride*(nlocal of s)
784: Not optimized; could be easily
786: Level: advanced
788: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
789: VecStrideScatterAll()
790: @*/
791: PetscErrorCode VecStrideGather(Vec v,PetscInt start,Vec s,InsertMode addv)
792: {
798: if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
799: 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);
800: if (!v->ops->stridegather) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
801: (*v->ops->stridegather)(v,start,s,addv);
802: return(0);
803: }
805: /*@
806: VecStrideScatter - Scatters a single component from a vector into a multi-component vector.
808: Collective on Vec
810: Input Parameters:
811: + s - the single-component vector
812: . start - starting point of the subvector (defined by a stride)
813: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
815: Output Parameter:
816: . v - the location where the subvector is scattered (the multi-component vector)
818: Notes:
819: One must call VecSetBlockSize() on the multi-component vector before this
820: routine to set the stride information, or use a vector created from a multicomponent DMDA.
822: The parallel layout of the vector and the subvector must be the same;
823: i.e., nlocal of v = stride*(nlocal of s)
825: Not optimized; could be easily
827: Level: advanced
829: .seealso: VecStrideNorm(), VecStrideGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
830: VecStrideScatterAll(), VecStrideSubSetScatter(), VecStrideSubSetGather()
831: @*/
832: PetscErrorCode VecStrideScatter(Vec s,PetscInt start,Vec v,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->stridescatter) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
842: (*v->ops->stridescatter)(s,start,v,addv);
843: return(0);
844: }
846: /*@
847: VecStrideSubSetGather - Gathers a subset of components from a multi-component vector into
848: another vector.
850: Collective on Vec
852: Input Parameters:
853: + v - the vector
854: . nidx - the number of indices
855: . idxv - the indices of the components 0 <= idxv[0] ...idxv[nidx-1] < bs(v), they need not be sorted
856: . 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
857: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
859: Output Parameter:
860: . s - the location where the subvector is stored
862: Notes:
863: One must call VecSetBlockSize() on both vectors before this routine to set the stride
864: information, or use a vector created from a multicomponent DMDA.
866: The parallel layout of the vector and the subvector must be the same;
868: Not optimized; could be easily
870: Level: advanced
872: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideGather(), VecStrideSubSetScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
873: VecStrideScatterAll()
874: @*/
875: PetscErrorCode VecStrideSubSetGather(Vec v,PetscInt nidx,const PetscInt idxv[],const PetscInt idxs[],Vec s,InsertMode addv)
876: {
882: if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
883: if (!v->ops->stridesubsetgather) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
884: (*v->ops->stridesubsetgather)(v,nidx,idxv,idxs,s,addv);
885: return(0);
886: }
888: /*@
889: VecStrideSubSetScatter - Scatters components from a vector into a subset of components of a multi-component vector.
891: Collective on Vec
893: Input Parameters:
894: + s - the smaller-component vector
895: . nidx - the number of indices in idx
896: . 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
897: . idxv - the indices of the components in the larger-component vector, 0 <= idx[0] ...idx[nidx-1] < bs(v) they need not be sorted
898: - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
900: Output Parameter:
901: . v - the location where the subvector is into scattered (the multi-component vector)
903: Notes:
904: One must call VecSetBlockSize() on the vectors before this
905: routine to set the stride information, or use a vector created from a multicomponent DMDA.
907: The parallel layout of the vector and the subvector must be the same;
909: Not optimized; could be easily
911: Level: advanced
913: .seealso: VecStrideNorm(), VecStrideGather(), VecStrideGather(), VecStrideSubSetGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
914: VecStrideScatterAll()
915: @*/
916: PetscErrorCode VecStrideSubSetScatter(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)
917: {
923: if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
924: if (!v->ops->stridesubsetscatter) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
925: (*v->ops->stridesubsetscatter)(s,nidx,idxs,idxv,v,addv);
926: return(0);
927: }
929: PetscErrorCode VecStrideGather_Default(Vec v,PetscInt start,Vec s,InsertMode addv)
930: {
932: PetscInt i,n,bs,ns;
933: const PetscScalar *x;
934: PetscScalar *y;
937: VecGetLocalSize(v,&n);
938: VecGetLocalSize(s,&ns);
939: VecGetArrayRead(v,&x);
940: VecGetArray(s,&y);
942: bs = v->map->bs;
943: 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);
944: x += start;
945: n = n/bs;
947: if (addv == INSERT_VALUES) {
948: for (i=0; i<n; i++) y[i] = x[bs*i];
949: } else if (addv == ADD_VALUES) {
950: for (i=0; i<n; i++) y[i] += x[bs*i];
951: #if !defined(PETSC_USE_COMPLEX)
952: } else if (addv == MAX_VALUES) {
953: for (i=0; i<n; i++) y[i] = PetscMax(y[i],x[bs*i]);
954: #endif
955: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
957: VecRestoreArrayRead(v,&x);
958: VecRestoreArray(s,&y);
959: return(0);
960: }
962: PetscErrorCode VecStrideScatter_Default(Vec s,PetscInt start,Vec v,InsertMode addv)
963: {
964: PetscErrorCode ierr;
965: PetscInt i,n,bs,ns;
966: PetscScalar *x;
967: const PetscScalar *y;
970: VecGetLocalSize(v,&n);
971: VecGetLocalSize(s,&ns);
972: VecGetArray(v,&x);
973: VecGetArrayRead(s,&y);
975: bs = v->map->bs;
976: 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);
977: x += start;
978: n = n/bs;
980: if (addv == INSERT_VALUES) {
981: for (i=0; i<n; i++) x[bs*i] = y[i];
982: } else if (addv == ADD_VALUES) {
983: for (i=0; i<n; i++) x[bs*i] += y[i];
984: #if !defined(PETSC_USE_COMPLEX)
985: } else if (addv == MAX_VALUES) {
986: for (i=0; i<n; i++) x[bs*i] = PetscMax(y[i],x[bs*i]);
987: #endif
988: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
990: VecRestoreArray(v,&x);
991: VecRestoreArrayRead(s,&y);
992: return(0);
993: }
995: PetscErrorCode VecStrideSubSetGather_Default(Vec v,PetscInt nidx,const PetscInt idxv[],const PetscInt idxs[],Vec s,InsertMode addv)
996: {
997: PetscErrorCode ierr;
998: PetscInt i,j,n,bs,bss,ns;
999: const PetscScalar *x;
1000: PetscScalar *y;
1003: VecGetLocalSize(v,&n);
1004: VecGetLocalSize(s,&ns);
1005: VecGetArrayRead(v,&x);
1006: VecGetArray(s,&y);
1008: bs = v->map->bs;
1009: bss = s->map->bs;
1010: n = n/bs;
1012: if (PetscDefined(USE_DEBUG)) {
1013: if (n != ns/bss) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible layout of vectors");
1014: for (j=0; j<nidx; j++) {
1015: if (idxv[j] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is negative",j,idxv[j]);
1016: if (idxv[j] >= bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is greater than or equal to vector blocksize %D",j,idxv[j],bs);
1017: }
1018: if (!idxs && bss != nidx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide idxs when not gathering into all locations");
1019: }
1021: if (addv == INSERT_VALUES) {
1022: if (!idxs) {
1023: for (i=0; i<n; i++) {
1024: for (j=0; j<bss; j++) y[bss*i+j] = x[bs*i+idxv[j]];
1025: }
1026: } else {
1027: for (i=0; i<n; i++) {
1028: for (j=0; j<bss; j++) y[bss*i+idxs[j]] = x[bs*i+idxv[j]];
1029: }
1030: }
1031: } else if (addv == ADD_VALUES) {
1032: if (!idxs) {
1033: for (i=0; i<n; i++) {
1034: for (j=0; j<bss; j++) y[bss*i+j] += x[bs*i+idxv[j]];
1035: }
1036: } else {
1037: for (i=0; i<n; i++) {
1038: for (j=0; j<bss; j++) y[bss*i+idxs[j]] += x[bs*i+idxv[j]];
1039: }
1040: }
1041: #if !defined(PETSC_USE_COMPLEX)
1042: } else if (addv == MAX_VALUES) {
1043: if (!idxs) {
1044: for (i=0; i<n; i++) {
1045: for (j=0; j<bss; j++) y[bss*i+j] = PetscMax(y[bss*i+j],x[bs*i+idxv[j]]);
1046: }
1047: } else {
1048: for (i=0; i<n; i++) {
1049: for (j=0; j<bss; j++) y[bss*i+idxs[j]] = PetscMax(y[bss*i+idxs[j]],x[bs*i+idxv[j]]);
1050: }
1051: }
1052: #endif
1053: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
1055: VecRestoreArrayRead(v,&x);
1056: VecRestoreArray(s,&y);
1057: return(0);
1058: }
1060: PetscErrorCode VecStrideSubSetScatter_Default(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)
1061: {
1062: PetscErrorCode ierr;
1063: PetscInt j,i,n,bs,ns,bss;
1064: PetscScalar *x;
1065: const PetscScalar *y;
1068: VecGetLocalSize(v,&n);
1069: VecGetLocalSize(s,&ns);
1070: VecGetArray(v,&x);
1071: VecGetArrayRead(s,&y);
1073: bs = v->map->bs;
1074: bss = s->map->bs;
1075: n = n/bs;
1077: if (PetscDefined(USE_DEBUG)) {
1078: if (n != ns/bss) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible layout of vectors");
1079: for (j=0; j<bss; j++) {
1080: if (idxs) {
1081: if (idxs[j] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is negative",j,idxs[j]);
1082: if (idxs[j] >= bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is greater than or equal to vector blocksize %D",j,idxs[j],bs);
1083: }
1084: }
1085: if (!idxs && bss != nidx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide idxs when not scattering from all locations");
1086: }
1088: if (addv == INSERT_VALUES) {
1089: if (!idxs) {
1090: for (i=0; i<n; i++) {
1091: for (j=0; j<bss; j++) x[bs*i + idxv[j]] = y[bss*i+j];
1092: }
1093: } else {
1094: for (i=0; i<n; i++) {
1095: for (j=0; j<bss; j++) x[bs*i + idxv[j]] = y[bss*i+idxs[j]];
1096: }
1097: }
1098: } else if (addv == ADD_VALUES) {
1099: if (!idxs) {
1100: for (i=0; i<n; i++) {
1101: for (j=0; j<bss; j++) x[bs*i + idxv[j]] += y[bss*i+j];
1102: }
1103: } else {
1104: for (i=0; i<n; i++) {
1105: for (j=0; j<bss; j++) x[bs*i + idxv[j]] += y[bss*i+idxs[j]];
1106: }
1107: }
1108: #if !defined(PETSC_USE_COMPLEX)
1109: } else if (addv == MAX_VALUES) {
1110: if (!idxs) {
1111: for (i=0; i<n; i++) {
1112: for (j=0; j<bss; j++) x[bs*i + idxv[j]] = PetscMax(y[bss*i+j],x[bs*i + idxv[j]]);
1113: }
1114: } else {
1115: for (i=0; i<n; i++) {
1116: for (j=0; j<bss; j++) x[bs*i + idxv[j]] = PetscMax(y[bss*i+idxs[j]],x[bs*i + idxv[j]]);
1117: }
1118: }
1119: #endif
1120: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
1122: VecRestoreArray(v,&x);
1123: VecRestoreArrayRead(s,&y);
1124: return(0);
1125: }
1127: PetscErrorCode VecReciprocal_Default(Vec v)
1128: {
1130: PetscInt i,n;
1131: PetscScalar *x;
1134: VecGetLocalSize(v,&n);
1135: VecGetArray(v,&x);
1136: for (i=0; i<n; i++) {
1137: if (x[i] != (PetscScalar)0.0) x[i] = (PetscScalar)1.0/x[i];
1138: }
1139: VecRestoreArray(v,&x);
1140: return(0);
1141: }
1143: /*@
1144: VecExp - Replaces each component of a vector by e^x_i
1146: Not collective
1148: Input Parameter:
1149: . v - The vector
1151: Output Parameter:
1152: . v - The vector of exponents
1154: Level: beginner
1156: .seealso: VecLog(), VecAbs(), VecSqrtAbs(), VecReciprocal()
1158: @*/
1159: PetscErrorCode VecExp(Vec v)
1160: {
1161: PetscScalar *x;
1162: PetscInt i, n;
1167: if (v->ops->exp) {
1168: (*v->ops->exp)(v);
1169: } else {
1170: VecGetLocalSize(v, &n);
1171: VecGetArray(v, &x);
1172: for (i = 0; i < n; i++) x[i] = PetscExpScalar(x[i]);
1173: VecRestoreArray(v, &x);
1174: }
1175: return(0);
1176: }
1178: /*@
1179: VecLog - Replaces each component of a vector by log(x_i), the natural logarithm
1181: Not collective
1183: Input Parameter:
1184: . v - The vector
1186: Output Parameter:
1187: . v - The vector of logs
1189: Level: beginner
1191: .seealso: VecExp(), VecAbs(), VecSqrtAbs(), VecReciprocal()
1193: @*/
1194: PetscErrorCode VecLog(Vec v)
1195: {
1196: PetscScalar *x;
1197: PetscInt i, n;
1202: if (v->ops->log) {
1203: (*v->ops->log)(v);
1204: } else {
1205: VecGetLocalSize(v, &n);
1206: VecGetArray(v, &x);
1207: for (i = 0; i < n; i++) x[i] = PetscLogScalar(x[i]);
1208: VecRestoreArray(v, &x);
1209: }
1210: return(0);
1211: }
1213: /*@
1214: VecSqrtAbs - Replaces each component of a vector by the square root of its magnitude.
1216: Not collective
1218: Input Parameter:
1219: . v - The vector
1221: Output Parameter:
1222: . v - The vector square root
1224: Level: beginner
1226: Note: The actual function is sqrt(|x_i|)
1228: .seealso: VecLog(), VecExp(), VecReciprocal(), VecAbs()
1230: @*/
1231: PetscErrorCode VecSqrtAbs(Vec v)
1232: {
1233: PetscScalar *x;
1234: PetscInt i, n;
1239: if (v->ops->sqrt) {
1240: (*v->ops->sqrt)(v);
1241: } else {
1242: VecGetLocalSize(v, &n);
1243: VecGetArray(v, &x);
1244: for (i = 0; i < n; i++) x[i] = PetscSqrtReal(PetscAbsScalar(x[i]));
1245: VecRestoreArray(v, &x);
1246: }
1247: return(0);
1248: }
1250: /*@
1251: VecDotNorm2 - computes the inner product of two vectors and the 2-norm squared of the second vector
1253: Collective on Vec
1255: Input Parameters:
1256: + s - first vector
1257: - t - second vector
1259: Output Parameters:
1260: + dp - s'conj(t)
1261: - nm - t'conj(t)
1263: Level: advanced
1265: Notes:
1266: conj(x) is the complex conjugate of x when x is complex
1268: .seealso: VecDot(), VecNorm(), VecDotBegin(), VecNormBegin(), VecDotEnd(), VecNormEnd()
1270: @*/
1271: PetscErrorCode VecDotNorm2(Vec s,Vec t,PetscScalar *dp, PetscReal *nm)
1272: {
1273: const PetscScalar *sx, *tx;
1274: PetscScalar dpx = 0.0, nmx = 0.0,work[2],sum[2];
1275: PetscInt i, n;
1276: PetscErrorCode ierr;
1286: if (s->map->N != t->map->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
1287: if (s->map->n != t->map->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");
1289: PetscLogEventBegin(VEC_DotNorm2,s,t,0,0);
1290: if (s->ops->dotnorm2) {
1291: (*s->ops->dotnorm2)(s,t,dp,&dpx);
1292: *nm = PetscRealPart(dpx);
1293: } else {
1294: VecGetLocalSize(s, &n);
1295: VecGetArrayRead(s, &sx);
1296: VecGetArrayRead(t, &tx);
1298: for (i = 0; i<n; i++) {
1299: dpx += sx[i]*PetscConj(tx[i]);
1300: nmx += tx[i]*PetscConj(tx[i]);
1301: }
1302: work[0] = dpx;
1303: work[1] = nmx;
1305: MPIU_Allreduce(work,sum,2,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)s));
1306: *dp = sum[0];
1307: *nm = PetscRealPart(sum[1]);
1309: VecRestoreArrayRead(t, &tx);
1310: VecRestoreArrayRead(s, &sx);
1311: PetscLogFlops(4.0*n);
1312: }
1313: PetscLogEventEnd(VEC_DotNorm2,s,t,0,0);
1314: return(0);
1315: }
1317: /*@
1318: VecSum - Computes the sum of all the components of a vector.
1320: Collective on Vec
1322: Input Parameter:
1323: . v - the vector
1325: Output Parameter:
1326: . sum - the result
1328: Level: beginner
1330: .seealso: VecMean(), VecNorm()
1331: @*/
1332: PetscErrorCode VecSum(Vec v,PetscScalar *sum)
1333: {
1334: PetscErrorCode ierr;
1335: PetscInt i,n;
1336: const PetscScalar *x;
1341: *sum = 0.0;
1342: if (v->ops->sum) {
1343: (*v->ops->sum)(v,sum);
1344: } else {
1345: VecGetLocalSize(v,&n);
1346: VecGetArrayRead(v,&x);
1347: for (i=0; i<n; i++) *sum += x[i];
1348: VecRestoreArrayRead(v,&x);
1349: }
1350: MPIU_Allreduce(MPI_IN_PLACE,sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)v));
1351: return(0);
1352: }
1354: /*@
1355: VecMean - Computes the arithmetic mean of all the components of a vector.
1357: Collective on Vec
1359: Input Parameter:
1360: . v - the vector
1362: Output Parameter:
1363: . mean - the result
1365: Level: beginner
1367: .seealso: VecSum(), VecNorm()
1368: @*/
1369: PetscErrorCode VecMean(Vec v,PetscScalar *mean)
1370: {
1371: PetscErrorCode ierr;
1372: PetscInt n;
1377: VecGetSize(v,&n);
1378: VecSum(v,mean);
1379: *mean /= n;
1380: return(0);
1381: }
1383: /*@
1384: VecImaginaryPart - Replaces a complex vector with its imginary part
1386: Collective on Vec
1388: Input Parameter:
1389: . v - the vector
1391: Level: beginner
1393: .seealso: VecNorm(), VecRealPart()
1394: @*/
1395: PetscErrorCode VecImaginaryPart(Vec v)
1396: {
1397: PetscErrorCode ierr;
1398: PetscInt i,n;
1399: PetscScalar *x;
1403: VecGetLocalSize(v,&n);
1404: VecGetArray(v,&x);
1405: for (i=0; i<n; i++) x[i] = PetscImaginaryPart(x[i]);
1406: VecRestoreArray(v,&x);
1407: return(0);
1408: }
1410: /*@
1411: VecRealPart - Replaces a complex vector with its real part
1413: Collective on Vec
1415: Input Parameter:
1416: . v - the vector
1418: Level: beginner
1420: .seealso: VecNorm(), VecImaginaryPart()
1421: @*/
1422: PetscErrorCode VecRealPart(Vec v)
1423: {
1424: PetscErrorCode ierr;
1425: PetscInt i,n;
1426: PetscScalar *x;
1430: VecGetLocalSize(v,&n);
1431: VecGetArray(v,&x);
1432: for (i=0; i<n; i++) x[i] = PetscRealPart(x[i]);
1433: VecRestoreArray(v,&x);
1434: return(0);
1435: }
1437: /*@
1438: VecShift - Shifts all of the components of a vector by computing
1439: x[i] = x[i] + shift.
1441: Logically Collective on Vec
1443: Input Parameters:
1444: + v - the vector
1445: - shift - the shift
1447: Level: intermediate
1449: @*/
1450: PetscErrorCode VecShift(Vec v,PetscScalar shift)
1451: {
1453: PetscInt i,n;
1454: PetscScalar *x;
1459: VecSetErrorIfLocked(v,1);
1460: if (shift == 0.0) return(0);
1462: if (v->ops->shift) {
1463: (*v->ops->shift)(v,shift);
1464: } else {
1465: VecGetLocalSize(v,&n);
1466: VecGetArray(v,&x);
1467: for (i=0; i<n; i++) x[i] += shift;
1468: VecRestoreArray(v,&x);
1469: }
1470: return(0);
1471: }
1473: /*@
1474: VecAbs - Replaces every element in a vector with its absolute value.
1476: Logically Collective on Vec
1478: Input Parameters:
1479: . v - the vector
1481: Level: intermediate
1483: @*/
1484: PetscErrorCode VecAbs(Vec v)
1485: {
1487: PetscInt i,n;
1488: PetscScalar *x;
1492: VecSetErrorIfLocked(v,1);
1494: if (v->ops->abs) {
1495: (*v->ops->abs)(v);
1496: } else {
1497: VecGetLocalSize(v,&n);
1498: VecGetArray(v,&x);
1499: for (i=0; i<n; i++) x[i] = PetscAbsScalar(x[i]);
1500: VecRestoreArray(v,&x);
1501: }
1502: return(0);
1503: }
1505: /*@
1506: VecPermute - Permutes a vector in place using the given ordering.
1508: Input Parameters:
1509: + vec - The vector
1510: . order - The ordering
1511: - inv - The flag for inverting the permutation
1513: Level: beginner
1515: Note: This function does not yet support parallel Index Sets with non-local permutations
1517: .seealso: MatPermute()
1518: @*/
1519: PetscErrorCode VecPermute(Vec x, IS row, PetscBool inv)
1520: {
1521: const PetscScalar *array;
1522: PetscScalar *newArray;
1523: const PetscInt *idx;
1524: PetscInt i,rstart,rend;
1525: PetscErrorCode ierr;
1530: VecSetErrorIfLocked(x,1);
1531: VecGetOwnershipRange(x,&rstart,&rend);
1532: ISGetIndices(row, &idx);
1533: VecGetArrayRead(x, &array);
1534: PetscMalloc1(x->map->n, &newArray);
1535: if (PetscDefined(USE_DEBUG)) {
1536: for (i = 0; i < x->map->n; i++) {
1537: 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]);
1538: }
1539: }
1540: if (!inv) {
1541: for (i = 0; i < x->map->n; i++) newArray[i] = array[idx[i]-rstart];
1542: } else {
1543: for (i = 0; i < x->map->n; i++) newArray[idx[i]-rstart] = array[i];
1544: }
1545: VecRestoreArrayRead(x, &array);
1546: ISRestoreIndices(row, &idx);
1547: VecReplaceArray(x, newArray);
1548: return(0);
1549: }
1551: /*@
1552: VecEqual - Compares two vectors. Returns true if the two vectors are either pointing to the same memory buffer,
1553: or if the two vectors have the same local and global layout as well as bitwise equality of all entries.
1554: Does NOT take round-off errors into account.
1556: Collective on Vec
1558: Input Parameters:
1559: + vec1 - the first vector
1560: - vec2 - the second vector
1562: Output Parameter:
1563: . flg - PETSC_TRUE if the vectors are equal; PETSC_FALSE otherwise.
1565: Level: intermediate
1566: @*/
1567: PetscErrorCode VecEqual(Vec vec1,Vec vec2,PetscBool *flg)
1568: {
1569: const PetscScalar *v1,*v2;
1570: PetscErrorCode ierr;
1571: PetscInt n1,n2,N1,N2;
1572: PetscBool flg1;
1578: if (vec1 == vec2) *flg = PETSC_TRUE;
1579: else {
1580: VecGetSize(vec1,&N1);
1581: VecGetSize(vec2,&N2);
1582: if (N1 != N2) flg1 = PETSC_FALSE;
1583: else {
1584: VecGetLocalSize(vec1,&n1);
1585: VecGetLocalSize(vec2,&n2);
1586: if (n1 != n2) flg1 = PETSC_FALSE;
1587: else {
1588: VecGetArrayRead(vec1,&v1);
1589: VecGetArrayRead(vec2,&v2);
1590: PetscArraycmp(v1,v2,n1,&flg1);
1591: VecRestoreArrayRead(vec1,&v1);
1592: VecRestoreArrayRead(vec2,&v2);
1593: }
1594: }
1595: /* combine results from all processors */
1596: MPIU_Allreduce(&flg1,flg,1,MPIU_BOOL,MPI_MIN,PetscObjectComm((PetscObject)vec1));
1597: }
1598: return(0);
1599: }
1601: /*@
1602: VecUniqueEntries - Compute the number of unique entries, and those entries
1604: Collective on Vec
1606: Input Parameter:
1607: . vec - the vector
1609: Output Parameters:
1610: + n - The number of unique entries
1611: - e - The entries
1613: Level: intermediate
1615: @*/
1616: PetscErrorCode VecUniqueEntries(Vec vec, PetscInt *n, PetscScalar **e)
1617: {
1618: const PetscScalar *v;
1619: PetscScalar *tmp, *vals;
1620: PetscMPIInt *N, *displs, l;
1621: PetscInt ng, m, i, j, p;
1622: PetscMPIInt size;
1623: PetscErrorCode ierr;
1628: MPI_Comm_size(PetscObjectComm((PetscObject) vec), &size);
1629: VecGetLocalSize(vec, &m);
1630: VecGetArrayRead(vec, &v);
1631: PetscMalloc2(m,&tmp,size,&N);
1632: for (i = 0, j = 0, l = 0; i < m; ++i) {
1633: /* Can speed this up with sorting */
1634: for (j = 0; j < l; ++j) {
1635: if (v[i] == tmp[j]) break;
1636: }
1637: if (j == l) {
1638: tmp[j] = v[i];
1639: ++l;
1640: }
1641: }
1642: VecRestoreArrayRead(vec, &v);
1643: /* Gather serial results */
1644: MPI_Allgather(&l, 1, MPI_INT, N, 1, MPI_INT, PetscObjectComm((PetscObject) vec));
1645: for (p = 0, ng = 0; p < size; ++p) {
1646: ng += N[p];
1647: }
1648: PetscMalloc2(ng,&vals,size+1,&displs);
1649: for (p = 1, displs[0] = 0; p <= size; ++p) {
1650: displs[p] = displs[p-1] + N[p-1];
1651: }
1652: MPI_Allgatherv(tmp, l, MPIU_SCALAR, vals, N, displs, MPIU_SCALAR, PetscObjectComm((PetscObject) vec));
1653: /* Find unique entries */
1654: #ifdef PETSC_USE_COMPLEX
1655: SETERRQ(PetscObjectComm((PetscObject) vec), PETSC_ERR_SUP, "Does not work with complex numbers");
1656: #else
1657: *n = displs[size];
1658: PetscSortRemoveDupsReal(n, (PetscReal *) vals);
1659: if (e) {
1661: PetscMalloc1(*n, e);
1662: for (i = 0; i < *n; ++i) {
1663: (*e)[i] = vals[i];
1664: }
1665: }
1666: PetscFree2(vals,displs);
1667: PetscFree2(tmp,N);
1668: return(0);
1669: #endif
1670: }