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;
12: if (nv > 128) {
13: PetscMalloc1(nv,&work);
14: }
15: VecMDot_Seq(xin,nv,y,work);
16: MPIU_Allreduce(work,z,nv,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)xin));
17: if (nv > 128) {
18: PetscFree(work);
19: }
20: return 0;
21: }
23: PetscErrorCode VecMTDot_MPI(Vec xin,PetscInt nv,const Vec y[],PetscScalar *z)
24: {
25: PetscScalar awork[128],*work = awork;
27: if (nv > 128) {
28: PetscMalloc1(nv,&work);
29: }
30: VecMTDot_Seq(xin,nv,y,work);
31: MPIU_Allreduce(work,z,nv,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)xin));
32: if (nv > 128) {
33: PetscFree(work);
34: }
35: return 0;
36: }
38: #include <../src/vec/vec/impls/seq/ftn-kernels/fnorm.h>
39: PetscErrorCode VecNorm_MPI(Vec xin,NormType type,PetscReal *z)
40: {
41: PetscReal sum,work = 0.0;
42: const PetscScalar *xx;
43: PetscInt n = xin->map->n;
44: PetscBLASInt one = 1,bn = 0;
46: PetscBLASIntCast(n,&bn);
47: if (type == NORM_2 || type == NORM_FROBENIUS) {
48: VecGetArrayRead(xin,&xx);
49: work = PetscRealPart(BLASdot_(&bn,xx,&one,xx,&one));
50: VecRestoreArrayRead(xin,&xx);
51: MPIU_Allreduce(&work,&sum,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)xin));
52: *z = PetscSqrtReal(sum);
53: PetscLogFlops(2.0*xin->map->n);
54: } else if (type == NORM_1) {
55: /* Find the local part */
56: VecNorm_Seq(xin,NORM_1,&work);
57: /* Find the global max */
58: MPIU_Allreduce(&work,z,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)xin));
59: } else if (type == NORM_INFINITY) {
60: /* Find the local max */
61: VecNorm_Seq(xin,NORM_INFINITY,&work);
62: /* Find the global max */
63: MPIU_Allreduce(&work,z,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)xin));
64: } else if (type == NORM_1_AND_2) {
65: PetscReal temp[2];
66: VecNorm_Seq(xin,NORM_1,temp);
67: VecNorm_Seq(xin,NORM_2,temp+1);
68: temp[1] = temp[1]*temp[1];
69: MPIU_Allreduce(temp,z,2,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)xin));
70: z[1] = PetscSqrtReal(z[1]);
71: }
72: return 0;
73: }
75: PetscErrorCode VecMax_MPI(Vec xin,PetscInt *idx,PetscReal *z)
76: {
77: PetscReal work;
79: /* Find the local max */
80: VecMax_Seq(xin,idx,&work);
81: #if defined(PETSC_HAVE_MPIUNI)
82: *z = work;
83: #else
84: /* Find the global max */
85: if (!idx) {
86: MPIU_Allreduce(&work,z,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)xin));
87: } else {
88: struct { PetscReal v; PetscInt i; } in,out;
90: in.v = work;
91: in.i = *idx + xin->map->rstart;
92: MPIU_Allreduce(&in,&out,1,MPIU_REAL_INT,MPIU_MAXLOC,PetscObjectComm((PetscObject)xin));
93: *z = out.v;
94: *idx = out.i;
95: }
96: #endif
97: return 0;
98: }
100: PetscErrorCode VecMin_MPI(Vec xin,PetscInt *idx,PetscReal *z)
101: {
102: PetscReal work;
104: /* Find the local Min */
105: VecMin_Seq(xin,idx,&work);
106: #if defined(PETSC_HAVE_MPIUNI)
107: *z = work;
108: #else
109: /* Find the global Min */
110: if (!idx) {
111: MPIU_Allreduce(&work,z,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)xin));
112: } else {
113: struct { PetscReal v; PetscInt i; } in,out;
115: in.v = work;
116: in.i = *idx + xin->map->rstart;
117: MPIU_Allreduce(&in,&out,1,MPIU_REAL_INT,MPIU_MINLOC,PetscObjectComm((PetscObject)xin));
118: *z = out.v;
119: *idx = out.i;
120: }
121: #endif
122: return 0;
123: }