Actual source code: veccomp0.h
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, 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: static PetscErrorCode __SUF__(VecDot_Comp)(Vec a,Vec b,PetscScalar *z)
23: {
24: PetscScalar sum = 0.0,work;
25: PetscInt i;
26: Vec_Comp *as = (Vec_Comp*)a->data,*bs = (Vec_Comp*)b->data;
28: PetscFunctionBegin;
29: SlepcValidVecComp(a,1);
30: SlepcValidVecComp(b,2);
31: if (as->x[0]->ops->dot_local) {
32: for (i=0;i<as->n->n;i++) {
33: PetscUseTypeMethod(as->x[i],dot_local,bs->x[i],&work);
34: sum += work;
35: }
36: #if defined(__WITH_MPI__)
37: work = sum;
38: PetscCallMPI(MPIU_Allreduce(&work,&sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)a)));
39: #endif
40: } else {
41: for (i=0;i<as->n->n;i++) {
42: PetscCall(VecDot(as->x[i],bs->x[i],&work));
43: sum += work;
44: }
45: }
46: *z = sum;
47: PetscFunctionReturn(PETSC_SUCCESS);
48: }
50: static PetscErrorCode __SUF__(VecMDot_Comp)(Vec a,PetscInt n,const Vec b[],PetscScalar *z)
51: {
52: PetscScalar *work,*work0,*r;
53: Vec_Comp *as = (Vec_Comp*)a->data;
54: Vec *bx;
55: PetscInt i,j;
57: PetscFunctionBegin;
58: SlepcValidVecComp(a,1);
59: SlepcValidVecsComp(b,n,3);
61: if (as->n->n == 0) {
62: *z = 0;
63: PetscFunctionReturn(PETSC_SUCCESS);
64: }
66: PetscCall(PetscMalloc2(n,&work0,n,&bx));
68: #if defined(__WITH_MPI__)
69: if (as->x[0]->ops->mdot_local) {
70: r = work0;
71: work = z;
72: } else
73: #endif
74: {
75: r = z;
76: 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) PetscUseTypeMethod(as->x[j],mdot_local,n,bx,work);
84: else PetscCall(VecMDot(as->x[j],n,bx,work));
85: for (i=0;i<n;i++) r[i] += work[i];
86: }
88: /* If def(__WITH_MPI__) and exists mdot_local */
89: #if defined(__WITH_MPI__)
90: if (as->x[0]->ops->mdot_local) {
91: /* z[i] <- Allreduce(work[i]) */
92: PetscCallMPI(MPIU_Allreduce(r,z,n,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)a)));
93: }
94: #endif
96: PetscCall(PetscFree2(work0,bx));
97: PetscFunctionReturn(PETSC_SUCCESS);
98: }
100: static PetscErrorCode __SUF__(VecTDot_Comp)(Vec a,Vec b,PetscScalar *z)
101: {
102: PetscScalar sum = 0.0,work;
103: PetscInt i;
104: Vec_Comp *as = (Vec_Comp*)a->data,*bs = (Vec_Comp*)b->data;
106: PetscFunctionBegin;
107: SlepcValidVecComp(a,1);
108: SlepcValidVecComp(b,2);
109: if (as->x[0]->ops->tdot_local) {
110: for (i=0;i<as->n->n;i++) {
111: PetscUseTypeMethod(as->x[i],tdot_local,bs->x[i],&work);
112: sum += work;
113: }
114: #if defined(__WITH_MPI__)
115: work = sum;
116: PetscCallMPI(MPIU_Allreduce(&work,&sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)a)));
117: #endif
118: } else {
119: for (i=0;i<as->n->n;i++) {
120: PetscCall(VecTDot(as->x[i],bs->x[i],&work));
121: sum += work;
122: }
123: }
124: *z = sum;
125: PetscFunctionReturn(PETSC_SUCCESS);
126: }
128: static PetscErrorCode __SUF__(VecMTDot_Comp)(Vec a,PetscInt n,const Vec b[],PetscScalar *z)
129: {
130: PetscScalar *work,*work0,*r;
131: Vec_Comp *as = (Vec_Comp*)a->data;
132: Vec *bx;
133: PetscInt i,j;
135: PetscFunctionBegin;
136: SlepcValidVecComp(a,1);
137: SlepcValidVecsComp(b,n,3);
139: if (as->n->n == 0) {
140: *z = 0;
141: PetscFunctionReturn(PETSC_SUCCESS);
142: }
144: PetscCall(PetscMalloc2(n,&work0,n,&bx));
146: #if defined(__WITH_MPI__)
147: if (as->x[0]->ops->mtdot_local) {
148: r = work0;
149: work = z;
150: } else
151: #endif
152: {
153: r = z;
154: work = work0;
155: }
157: /* z[i] <- a.x' * b[i].x */
158: for (i=0;i<n;i++) r[i] = 0.0;
159: for (j=0;j<as->n->n;j++) {
160: for (i=0;i<n;i++) bx[i] = ((Vec_Comp*)b[i]->data)->x[j];
161: if (as->x[0]->ops->mtdot_local) PetscUseTypeMethod(as->x[j],mtdot_local,n,bx,work);
162: else PetscCall(VecMTDot(as->x[j],n,bx,work));
163: for (i=0;i<n;i++) r[i] += work[i];
164: }
166: /* If def(__WITH_MPI__) and exists mtdot_local */
167: #if defined(__WITH_MPI__)
168: if (as->x[0]->ops->mtdot_local) {
169: /* z[i] <- Allreduce(work[i]) */
170: PetscCallMPI(MPIU_Allreduce(r,z,n,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)a)));
171: }
172: #endif
174: PetscCall(PetscFree2(work0,bx));
175: PetscFunctionReturn(PETSC_SUCCESS);
176: }
178: static PetscErrorCode __SUF__(VecNorm_Comp)(Vec a,NormType t,PetscReal *norm)
179: {
180: PetscReal work[3],s=0.0;
181: Vec_Comp *as = (Vec_Comp*)a->data;
182: PetscInt i;
184: PetscFunctionBegin;
185: SlepcValidVecComp(a,1);
186: /* Initialize norm */
187: switch (t) {
188: case NORM_1:
189: case NORM_INFINITY:
190: *norm = 0.0;
191: break;
192: case NORM_2:
193: case NORM_FROBENIUS:
194: *norm = 1.0;
195: s = 0.0;
196: break;
197: case NORM_1_AND_2:
198: norm[0] = 0.0;
199: norm[1] = 1.0;
200: s = 0.0;
201: break;
202: }
203: for (i=0;i<as->n->n;i++) {
204: if (as->x[0]->ops->norm_local) PetscUseTypeMethod(as->x[i],norm_local,t,work);
205: else PetscCall(VecNorm(as->x[i],t,work));
206: /* norm+= work */
207: switch (t) {
208: case NORM_1:
209: *norm += *work;
210: break;
211: case NORM_2:
212: case NORM_FROBENIUS:
213: AddNorm2(norm,&s,*work);
214: break;
215: case NORM_1_AND_2:
216: norm[0] += work[0];
217: AddNorm2(&norm[1],&s,work[1]);
218: break;
219: case NORM_INFINITY:
220: *norm = PetscMax(*norm,*work);
221: break;
222: }
223: }
225: /* If def(__WITH_MPI__) and exists norm_local */
226: #if defined(__WITH_MPI__)
227: if (as->x[0]->ops->norm_local) {
228: PetscReal work0[3];
229: /* norm <- Allreduce(work) */
230: switch (t) {
231: case NORM_1:
232: work[0] = *norm;
233: PetscCallMPI(MPIU_Allreduce(work,norm,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)a)));
234: break;
235: case NORM_2:
236: case NORM_FROBENIUS:
237: work[0] = *norm;
238: work[1] = s;
239: PetscCallMPI(MPIU_Allreduce(work,work0,1,MPIU_NORM2,MPIU_NORM2_SUM,PetscObjectComm((PetscObject)a)));
240: *norm = GetNorm2(work0[0],work0[1]);
241: break;
242: case NORM_1_AND_2:
243: work[0] = norm[0];
244: work[1] = norm[1];
245: work[2] = s;
246: PetscCallMPI(MPIU_Allreduce(work,work0,1,MPIU_NORM1_AND_2,MPIU_NORM2_SUM,PetscObjectComm((PetscObject)a)));
247: norm[0] = work0[0];
248: norm[1] = GetNorm2(work0[1],work0[2]);
249: break;
250: case NORM_INFINITY:
251: work[0] = *norm;
252: PetscCallMPI(MPIU_Allreduce(work,norm,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)a)));
253: break;
254: }
255: }
256: #else
257: /* Norm correction */
258: switch (t) {
259: case NORM_1:
260: case NORM_INFINITY:
261: break;
262: case NORM_2:
263: case NORM_FROBENIUS:
264: *norm = GetNorm2(*norm,s);
265: break;
266: case NORM_1_AND_2:
267: norm[1] = GetNorm2(norm[1],s);
268: break;
269: }
270: #endif
271: PetscFunctionReturn(PETSC_SUCCESS);
272: }
274: PetscErrorCode __SUF__(VecDotNorm2_Comp)(Vec v,Vec w,PetscScalar *dp,PetscScalar *nm)
275: {
276: PetscScalar dp0=0.0,nm0=0.0,dp1=0.0,nm1=0.0;
277: const PetscScalar *vx,*wx;
278: Vec_Comp *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data;
279: PetscInt i,n;
280: PetscBool t0,t1;
281: #if defined(__WITH_MPI__)
282: PetscScalar work[4];
283: #endif
285: PetscFunctionBegin;
286: /* Compute recursively the local part */
287: PetscCall(PetscObjectTypeCompare((PetscObject)v,VECCOMP,&t0));
288: PetscCall(PetscObjectTypeCompare((PetscObject)w,VECCOMP,&t1));
289: if (t0 && t1) {
290: SlepcValidVecComp(v,1);
291: SlepcValidVecComp(w,2);
292: for (i=0;i<vs->n->n;i++) {
293: PetscCall(VecDotNorm2_Comp_Seq(vs->x[i],ws->x[i],&dp1,&nm1));
294: dp0 += dp1;
295: nm0 += nm1;
296: }
297: } else if (!t0 && !t1) {
298: PetscCall(VecGetLocalSize(v,&n));
299: PetscCall(VecGetArrayRead(v,&vx));
300: PetscCall(VecGetArrayRead(w,&wx));
301: for (i=0;i<n;i++) {
302: dp0 += vx[i]*PetscConj(wx[i]);
303: nm0 += wx[i]*PetscConj(wx[i]);
304: }
305: PetscCall(VecRestoreArrayRead(v,&vx));
306: PetscCall(VecRestoreArrayRead(w,&wx));
307: } else SETERRQ(PetscObjectComm((PetscObject)v),PETSC_ERR_ARG_INCOMP,"Incompatible vector types");
309: #if defined(__WITH_MPI__)
310: /* [dp, nm] <- Allreduce([dp0, nm0]) */
311: work[0] = dp0;
312: work[1] = nm0;
313: PetscCallMPI(MPIU_Allreduce(work,&work[2],2,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)v)));
314: *dp = work[2];
315: *nm = work[3];
316: #else
317: *dp = dp0;
318: *nm = nm0;
319: #endif
320: PetscFunctionReturn(PETSC_SUCCESS);
321: }