Actual source code: pvec2.c
2: /*
3: Code for some of the parallel vector primatives.
4: */
5: #include <../src/vec/vec/impls/mpi/pvecimpl.h>
6: #include <petscblaslapack.h>
8: PetscErrorCode VecMDot_MPI(Vec xin,PetscInt nv,const Vec y[],PetscScalar *z)
9: {
10: PetscScalar awork[128],*work = awork;
14: if (nv > 128) {
15: PetscMalloc1(nv,&work);
16: }
17: VecMDot_Seq(xin,nv,y,work);
18: MPIU_Allreduce(work,z,nv,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)xin));
19: if (nv > 128) {
20: PetscFree(work);
21: }
22: return(0);
23: }
25: PetscErrorCode VecMTDot_MPI(Vec xin,PetscInt nv,const Vec y[],PetscScalar *z)
26: {
27: PetscScalar awork[128],*work = awork;
31: if (nv > 128) {
32: PetscMalloc1(nv,&work);
33: }
34: VecMTDot_Seq(xin,nv,y,work);
35: MPIU_Allreduce(work,z,nv,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)xin));
36: if (nv > 128) {
37: PetscFree(work);
38: }
39: return(0);
40: }
42: #include <../src/vec/vec/impls/seq/ftn-kernels/fnorm.h>
43: PetscErrorCode VecNorm_MPI(Vec xin,NormType type,PetscReal *z)
44: {
45: PetscReal sum,work = 0.0;
46: const PetscScalar *xx;
47: PetscErrorCode ierr;
48: PetscInt n = xin->map->n;
49: PetscBLASInt one = 1,bn = 0;
52: PetscBLASIntCast(n,&bn);
53: if (type == NORM_2 || type == NORM_FROBENIUS) {
54: VecGetArrayRead(xin,&xx);
55: work = PetscRealPart(BLASdot_(&bn,xx,&one,xx,&one));
56: VecRestoreArrayRead(xin,&xx);
57: MPIU_Allreduce(&work,&sum,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)xin));
58: *z = PetscSqrtReal(sum);
59: PetscLogFlops(2.0*xin->map->n);
60: } else if (type == NORM_1) {
61: /* Find the local part */
62: VecNorm_Seq(xin,NORM_1,&work);
63: /* Find the global max */
64: MPIU_Allreduce(&work,z,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)xin));
65: } else if (type == NORM_INFINITY) {
66: /* Find the local max */
67: VecNorm_Seq(xin,NORM_INFINITY,&work);
68: /* Find the global max */
69: MPIU_Allreduce(&work,z,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)xin));
70: } else if (type == NORM_1_AND_2) {
71: PetscReal temp[2];
72: VecNorm_Seq(xin,NORM_1,temp);
73: VecNorm_Seq(xin,NORM_2,temp+1);
74: temp[1] = temp[1]*temp[1];
75: MPIU_Allreduce(temp,z,2,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)xin));
76: z[1] = PetscSqrtReal(z[1]);
77: }
78: return(0);
79: }
81: PetscErrorCode VecMax_MPI(Vec xin,PetscInt *idx,PetscReal *z)
82: {
84: PetscReal work;
87: /* Find the local max */
88: VecMax_Seq(xin,idx,&work);
90: /* Find the global max */
91: if (!idx) {
92: MPIU_Allreduce(&work,z,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)xin));
93: } else {
94: PetscReal work2[2],z2[2];
95: PetscInt rstart;
96: rstart = xin->map->rstart;
97: work2[0] = work;
98: work2[1] = *idx + rstart;
99: MPIU_Allreduce(work2,z2,2,MPIU_REAL,MPIU_MAXINDEX_OP,PetscObjectComm((PetscObject)xin));
100: *z = z2[0];
101: *idx = (PetscInt)z2[1];
102: }
103: return(0);
104: }
106: PetscErrorCode VecMin_MPI(Vec xin,PetscInt *idx,PetscReal *z)
107: {
109: PetscReal work;
112: /* Find the local Min */
113: VecMin_Seq(xin,idx,&work);
115: /* Find the global Min */
116: if (!idx) {
117: MPIU_Allreduce(&work,z,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)xin));
118: } else {
119: PetscReal work2[2],z2[2];
120: PetscInt rstart;
122: VecGetOwnershipRange(xin,&rstart,NULL);
123: work2[0] = work;
124: work2[1] = *idx + rstart;
125: MPIU_Allreduce(work2,z2,2,MPIU_REAL,MPIU_MININDEX_OP,PetscObjectComm((PetscObject)xin));
126: *z = z2[0];
127: *idx = (PetscInt)z2[1];
128: }
129: return(0);
130: }