Actual source code: veccomp0.h
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
11: #include <petsc/private/vecimpl.h>
13: #if defined(__WITH_MPI__)
15: #else
17: #endif
22: PetscErrorCode __SUF__(VecDot_Comp)(Vec a,Vec b,PetscScalar *z)
23: {
24: PetscScalar sum = 0.0,work;
25: PetscInt i;
27: Vec_Comp *as = (Vec_Comp*)a->data,*bs = (Vec_Comp*)b->data;
30: SlepcValidVecComp(a,1);
31: SlepcValidVecComp(b,2);
32: if (as->x[0]->ops->dot_local) {
33: for (i=0,sum=0.0;i<as->n->n;i++) {
34: as->x[i]->ops->dot_local(as->x[i],bs->x[i],&work);
35: sum += work;
36: }
37: #if defined(__WITH_MPI__)
38: work = sum;
39: MPI_Allreduce(&work,&sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)a));
40: #endif
41: } else {
42: for (i=0,sum=0.0;i<as->n->n;i++) {
43: VecDot(as->x[i],bs->x[i],&work);
44: sum += work;
45: }
46: }
47: *z = sum;
48: return(0);
49: }
51: PetscErrorCode __SUF__(VecMDot_Comp)(Vec a,PetscInt n,const Vec b[],PetscScalar *z)
52: {
53: PetscScalar *work,*work0,*r;
55: Vec_Comp *as = (Vec_Comp*)a->data;
56: Vec *bx;
57: PetscInt i,j;
60: SlepcValidVecComp(a,1);
61: for (i=0;i<n;i++) SlepcValidVecComp(b[i],3);
63: if (as->n->n == 0) {
64: *z = 0;
65: return(0);
66: }
68: PetscMalloc2(n,&work0,n,&bx);
70: #if defined(__WITH_MPI__)
71: if (as->x[0]->ops->mdot_local) {
72: r = work0; work = z;
73: } else
74: #endif
75: {
76: r = z; work = work0;
77: }
79: /* z[i] <- a.x' * b[i].x */
80: for (i=0;i<n;i++) r[i] = 0.0;
81: for (j=0;j<as->n->n;j++) {
82: for (i=0;i<n;i++) bx[i] = ((Vec_Comp*)b[i]->data)->x[j];
83: if (as->x[0]->ops->mdot_local) {
84: as->x[j]->ops->mdot_local(as->x[j],n,bx,work);
85: } else {
86: VecMDot(as->x[j],n,bx,work);
87: }
88: for (i=0;i<n;i++) r[i] += work[i];
89: }
91: /* If def(__WITH_MPI__) and exists mdot_local */
92: #if defined(__WITH_MPI__)
93: if (as->x[0]->ops->mdot_local) {
94: /* z[i] <- Allreduce(work[i]) */
95: MPI_Allreduce(r,z,n,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)a));
96: }
97: #endif
99: PetscFree2(work0,bx);
100: return(0);
101: }
103: PetscErrorCode __SUF__(VecTDot_Comp)(Vec a,Vec b,PetscScalar *z)
104: {
105: PetscScalar sum = 0.0,work;
106: PetscInt i;
108: Vec_Comp *as = (Vec_Comp*)a->data,*bs = (Vec_Comp*)b->data;
111: SlepcValidVecComp(a,1);
112: SlepcValidVecComp(b,2);
113: if (as->x[0]->ops->tdot_local) {
114: for (i=0,sum=0.0;i<as->n->n;i++) {
115: as->x[i]->ops->tdot_local(as->x[i],bs->x[i],&work);
116: sum += work;
117: }
118: #if defined(__WITH_MPI__)
119: work = sum;
120: MPI_Allreduce(&work,&sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)a));
121: #endif
122: } else {
123: for (i=0,sum=0.0;i<as->n->n;i++) {
124: VecTDot(as->x[i],bs->x[i],&work);
125: sum += work;
126: }
127: }
128: *z = sum;
129: return(0);
130: }
132: PetscErrorCode __SUF__(VecMTDot_Comp)(Vec a,PetscInt n,const Vec b[],PetscScalar *z)
133: {
134: PetscScalar *work,*work0,*r;
136: Vec_Comp *as = (Vec_Comp*)a->data;
137: Vec *bx;
138: PetscInt i,j;
141: SlepcValidVecComp(a,1);
142: for (i=0;i<n;i++) SlepcValidVecComp(b[i],3);
144: if (as->n->n == 0) {
145: *z = 0;
146: return(0);
147: }
149: PetscMalloc2(n,&work0,n,&bx);
151: #if defined(__WITH_MPI__)
152: if (as->x[0]->ops->mtdot_local) {
153: r = work0; work = z;
154: } else
155: #endif
156: {
157: r = z; work = work0;
158: }
160: /* z[i] <- a.x' * b[i].x */
161: for (i=0;i<n;i++) r[i] = 0.0;
162: for (j=0;j<as->n->n;j++) {
163: for (i=0;i<n;i++) bx[i] = ((Vec_Comp*)b[i]->data)->x[j];
164: if (as->x[0]->ops->mtdot_local) {
165: as->x[j]->ops->mtdot_local(as->x[j],n,bx,work);
166: } else {
167: VecMTDot(as->x[j],n,bx,work);
168: }
169: for (i=0;i<n;i++) r[i] += work[i];
170: }
172: /* If def(__WITH_MPI__) and exists mtdot_local */
173: #if defined(__WITH_MPI__)
174: if (as->x[0]->ops->mtdot_local) {
175: /* z[i] <- Allreduce(work[i]) */
176: MPI_Allreduce(r,z,n,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)a));
177: }
178: #endif
180: PetscFree2(work0,bx);
181: return(0);
182: }
184: PetscErrorCode __SUF__(VecNorm_Comp)(Vec a,NormType t,PetscReal *norm)
185: {
186: PetscReal work[3],s=0.0;
188: Vec_Comp *as = (Vec_Comp*)a->data;
189: PetscInt i;
192: SlepcValidVecComp(a,1);
193: /* Initialize norm */
194: switch (t) {
195: case NORM_1: case NORM_INFINITY: *norm = 0.0; break;
196: case NORM_2: case NORM_FROBENIUS: *norm = 1.0; s = 0.0; break;
197: case NORM_1_AND_2: norm[0] = 0.0; norm[1] = 1.0; s = 0.0; break;
198: }
199: for (i=0;i<as->n->n;i++) {
200: if (as->x[0]->ops->norm_local) {
201: as->x[0]->ops->norm_local(as->x[i],t,work);
202: } else {
203: VecNorm(as->x[i],t,work);
204: }
205: /* norm+= work */
206: switch (t) {
207: case NORM_1: *norm+= *work; break;
208: case NORM_2: case NORM_FROBENIUS: AddNorm2(norm,&s,*work); break;
209: case NORM_1_AND_2: norm[0]+= work[0]; AddNorm2(&norm[1],&s,work[1]); break;
210: case NORM_INFINITY: *norm = PetscMax(*norm,*work); break;
211: }
212: }
214: /* If def(__WITH_MPI__) and exists norm_local */
215: #if defined(__WITH_MPI__)
216: if (as->x[0]->ops->norm_local) {
217: PetscReal work0[3];
218: /* norm <- Allreduce(work) */
219: switch (t) {
220: case NORM_1:
221: work[0] = *norm;
222: MPI_Allreduce(work,norm,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)a));
223: break;
224: case NORM_2: case NORM_FROBENIUS:
225: work[0] = *norm; work[1] = s;
226: MPI_Allreduce(work,work0,1,MPIU_NORM2,MPIU_NORM2_SUM,PetscObjectComm((PetscObject)a));
227: *norm = GetNorm2(work0[0],work0[1]);
228: break;
229: case NORM_1_AND_2:
230: work[0] = norm[0]; work[1] = norm[1]; work[2] = s;
231: MPI_Allreduce(work,work0,1,MPIU_NORM1_AND_2,MPIU_NORM2_SUM,PetscObjectComm((PetscObject)a));
232: norm[0] = work0[0];
233: norm[1] = GetNorm2(work0[1],work0[2]);
234: break;
235: case NORM_INFINITY:
236: work[0] = *norm;
237: MPI_Allreduce(work,norm,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)a));
238: break;
239: }
240: }
241: #else
242: /* Norm correction */
243: switch (t) {
244: case NORM_2: case NORM_FROBENIUS: *norm = GetNorm2(*norm,s); break;
245: case NORM_1_AND_2: norm[1] = GetNorm2(norm[1],s); break;
246: default: ;
247: }
248: #endif
249: return(0);
250: }
252: PetscErrorCode __SUF__(VecDotNorm2_Comp)(Vec v,Vec w,PetscScalar *dp,PetscScalar *nm)
253: {
254: PetscErrorCode ierr;
255: PetscScalar dp0,nm0,dp1,nm1;
256: const PetscScalar *vx,*wx;
257: Vec_Comp *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data;
258: PetscInt i,n;
259: PetscBool t0,t1;
260: #if defined(__WITH_MPI__)
261: PetscScalar work[4];
262: #endif
265: /* Compute recursively the local part */
266: dp0 = nm0 = 0.0;
267: PetscObjectTypeCompare((PetscObject)v,VECCOMP,&t0);
268: PetscObjectTypeCompare((PetscObject)w,VECCOMP,&t1);
269: if (t0 && t1) {
270: SlepcValidVecComp(v,1);
271: SlepcValidVecComp(w,2);
272: for (i=0;i<vs->n->n;i++) {
273: VecDotNorm2_Comp_Seq(vs->x[i],ws->x[i],&dp1,&nm1);
274: dp0 += dp1;
275: nm0 += nm1;
276: }
277: } else if (!t0 && !t1) {
278: VecGetLocalSize(v,&n);
279: VecGetArrayRead(v,&vx);
280: VecGetArrayRead(w,&wx);
281: for (i=0;i<n;i++) {
282: dp0 += vx[i]*PetscConj(wx[i]);
283: nm0 += wx[i]*PetscConj(wx[i]);
284: }
285: VecRestoreArrayRead(v,&vx);
286: VecRestoreArrayRead(w,&wx);
287: } else SETERRQ(PetscObjectComm((PetscObject)v),PETSC_ERR_ARG_INCOMP,"Incompatible vector types");
289: #if defined(__WITH_MPI__)
290: /* [dp, nm] <- Allreduce([dp0, nm0]) */
291: work[0] = dp0; work[1] = nm0;
292: MPI_Allreduce(work,&work[2],2,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)v));
293: *dp = work[2]; *nm = work[3];
294: #else
295: *dp = dp0; *nm = nm0;
296: #endif
297: return(0);
298: }