Actual source code: veccomp0.h
slepc-3.21.1 2024-04-26
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,sum=0.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: PetscCall(MPIU_Allreduce(&work,&sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)a)));
39: #endif
40: } else {
41: for (i=0,sum=0.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: for (i=0;i<n;i++) SlepcValidVecComp(b[i],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; work = z;
71: } else
72: #endif
73: {
74: r = z; work = work0;
75: }
77: /* z[i] <- a.x' * b[i].x */
78: for (i=0;i<n;i++) r[i] = 0.0;
79: for (j=0;j<as->n->n;j++) {
80: for (i=0;i<n;i++) bx[i] = ((Vec_Comp*)b[i]->data)->x[j];
81: if (as->x[0]->ops->mdot_local) PetscUseTypeMethod(as->x[j],mdot_local,n,bx,work);
82: else PetscCall(VecMDot(as->x[j],n,bx,work));
83: for (i=0;i<n;i++) r[i] += work[i];
84: }
86: /* If def(__WITH_MPI__) and exists mdot_local */
87: #if defined(__WITH_MPI__)
88: if (as->x[0]->ops->mdot_local) {
89: /* z[i] <- Allreduce(work[i]) */
90: PetscCall(MPIU_Allreduce(r,z,n,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)a)));
91: }
92: #endif
94: PetscCall(PetscFree2(work0,bx));
95: PetscFunctionReturn(PETSC_SUCCESS);
96: }
98: static PetscErrorCode __SUF__(VecTDot_Comp)(Vec a,Vec b,PetscScalar *z)
99: {
100: PetscScalar sum = 0.0,work;
101: PetscInt i;
102: Vec_Comp *as = (Vec_Comp*)a->data,*bs = (Vec_Comp*)b->data;
104: PetscFunctionBegin;
105: SlepcValidVecComp(a,1);
106: SlepcValidVecComp(b,2);
107: if (as->x[0]->ops->tdot_local) {
108: for (i=0,sum=0.0;i<as->n->n;i++) {
109: PetscUseTypeMethod(as->x[i],tdot_local,bs->x[i],&work);
110: sum += work;
111: }
112: #if defined(__WITH_MPI__)
113: work = sum;
114: PetscCall(MPIU_Allreduce(&work,&sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)a)));
115: #endif
116: } else {
117: for (i=0,sum=0.0;i<as->n->n;i++) {
118: PetscCall(VecTDot(as->x[i],bs->x[i],&work));
119: sum += work;
120: }
121: }
122: *z = sum;
123: PetscFunctionReturn(PETSC_SUCCESS);
124: }
126: static PetscErrorCode __SUF__(VecMTDot_Comp)(Vec a,PetscInt n,const Vec b[],PetscScalar *z)
127: {
128: PetscScalar *work,*work0,*r;
129: Vec_Comp *as = (Vec_Comp*)a->data;
130: Vec *bx;
131: PetscInt i,j;
133: PetscFunctionBegin;
134: SlepcValidVecComp(a,1);
135: for (i=0;i<n;i++) SlepcValidVecComp(b[i],3);
137: if (as->n->n == 0) {
138: *z = 0;
139: PetscFunctionReturn(PETSC_SUCCESS);
140: }
142: PetscCall(PetscMalloc2(n,&work0,n,&bx));
144: #if defined(__WITH_MPI__)
145: if (as->x[0]->ops->mtdot_local) {
146: r = work0; work = z;
147: } else
148: #endif
149: {
150: r = z; work = work0;
151: }
153: /* z[i] <- a.x' * b[i].x */
154: for (i=0;i<n;i++) r[i] = 0.0;
155: for (j=0;j<as->n->n;j++) {
156: for (i=0;i<n;i++) bx[i] = ((Vec_Comp*)b[i]->data)->x[j];
157: if (as->x[0]->ops->mtdot_local) PetscUseTypeMethod(as->x[j],mtdot_local,n,bx,work);
158: else PetscCall(VecMTDot(as->x[j],n,bx,work));
159: for (i=0;i<n;i++) r[i] += work[i];
160: }
162: /* If def(__WITH_MPI__) and exists mtdot_local */
163: #if defined(__WITH_MPI__)
164: if (as->x[0]->ops->mtdot_local) {
165: /* z[i] <- Allreduce(work[i]) */
166: PetscCall(MPIU_Allreduce(r,z,n,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)a)));
167: }
168: #endif
170: PetscCall(PetscFree2(work0,bx));
171: PetscFunctionReturn(PETSC_SUCCESS);
172: }
174: static PetscErrorCode __SUF__(VecNorm_Comp)(Vec a,NormType t,PetscReal *norm)
175: {
176: PetscReal work[3],s=0.0;
177: Vec_Comp *as = (Vec_Comp*)a->data;
178: PetscInt i;
180: PetscFunctionBegin;
181: SlepcValidVecComp(a,1);
182: /* Initialize norm */
183: switch (t) {
184: case NORM_1: case NORM_INFINITY: *norm = 0.0; break;
185: case NORM_2: case NORM_FROBENIUS: *norm = 1.0; s = 0.0; break;
186: case NORM_1_AND_2: norm[0] = 0.0; norm[1] = 1.0; s = 0.0; break;
187: }
188: for (i=0;i<as->n->n;i++) {
189: if (as->x[0]->ops->norm_local) PetscUseTypeMethod(as->x[i],norm_local,t,work);
190: else PetscCall(VecNorm(as->x[i],t,work));
191: /* norm+= work */
192: switch (t) {
193: case NORM_1: *norm+= *work; break;
194: case NORM_2: case NORM_FROBENIUS: AddNorm2(norm,&s,*work); break;
195: case NORM_1_AND_2: norm[0]+= work[0]; AddNorm2(&norm[1],&s,work[1]); break;
196: case NORM_INFINITY: *norm = PetscMax(*norm,*work); break;
197: }
198: }
200: /* If def(__WITH_MPI__) and exists norm_local */
201: #if defined(__WITH_MPI__)
202: if (as->x[0]->ops->norm_local) {
203: PetscReal work0[3];
204: /* norm <- Allreduce(work) */
205: switch (t) {
206: case NORM_1:
207: work[0] = *norm;
208: PetscCall(MPIU_Allreduce(work,norm,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)a)));
209: break;
210: case NORM_2: case NORM_FROBENIUS:
211: work[0] = *norm; work[1] = s;
212: PetscCall(MPIU_Allreduce(work,work0,1,MPIU_NORM2,MPIU_NORM2_SUM,PetscObjectComm((PetscObject)a)));
213: *norm = GetNorm2(work0[0],work0[1]);
214: break;
215: case NORM_1_AND_2:
216: work[0] = norm[0]; work[1] = norm[1]; work[2] = s;
217: PetscCall(MPIU_Allreduce(work,work0,1,MPIU_NORM1_AND_2,MPIU_NORM2_SUM,PetscObjectComm((PetscObject)a)));
218: norm[0] = work0[0];
219: norm[1] = GetNorm2(work0[1],work0[2]);
220: break;
221: case NORM_INFINITY:
222: work[0] = *norm;
223: PetscCall(MPIU_Allreduce(work,norm,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)a)));
224: break;
225: }
226: }
227: #else
228: /* Norm correction */
229: switch (t) {
230: case NORM_2: case NORM_FROBENIUS: *norm = GetNorm2(*norm,s); break;
231: case NORM_1_AND_2: norm[1] = GetNorm2(norm[1],s); break;
232: default: ;
233: }
234: #endif
235: PetscFunctionReturn(PETSC_SUCCESS);
236: }
238: PetscErrorCode __SUF__(VecDotNorm2_Comp)(Vec v,Vec w,PetscScalar *dp,PetscScalar *nm)
239: {
240: PetscScalar dp0=0.0,nm0=0.0,dp1=0.0,nm1=0.0;
241: const PetscScalar *vx,*wx;
242: Vec_Comp *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data;
243: PetscInt i,n;
244: PetscBool t0,t1;
245: #if defined(__WITH_MPI__)
246: PetscScalar work[4];
247: #endif
249: PetscFunctionBegin;
250: /* Compute recursively the local part */
251: PetscCall(PetscObjectTypeCompare((PetscObject)v,VECCOMP,&t0));
252: PetscCall(PetscObjectTypeCompare((PetscObject)w,VECCOMP,&t1));
253: if (t0 && t1) {
254: SlepcValidVecComp(v,1);
255: SlepcValidVecComp(w,2);
256: for (i=0;i<vs->n->n;i++) {
257: PetscCall(VecDotNorm2_Comp_Seq(vs->x[i],ws->x[i],&dp1,&nm1));
258: dp0 += dp1;
259: nm0 += nm1;
260: }
261: } else if (!t0 && !t1) {
262: PetscCall(VecGetLocalSize(v,&n));
263: PetscCall(VecGetArrayRead(v,&vx));
264: PetscCall(VecGetArrayRead(w,&wx));
265: for (i=0;i<n;i++) {
266: dp0 += vx[i]*PetscConj(wx[i]);
267: nm0 += wx[i]*PetscConj(wx[i]);
268: }
269: PetscCall(VecRestoreArrayRead(v,&vx));
270: PetscCall(VecRestoreArrayRead(w,&wx));
271: } else SETERRQ(PetscObjectComm((PetscObject)v),PETSC_ERR_ARG_INCOMP,"Incompatible vector types");
273: #if defined(__WITH_MPI__)
274: /* [dp, nm] <- Allreduce([dp0, nm0]) */
275: work[0] = dp0; work[1] = nm0;
276: PetscCall(MPIU_Allreduce(work,&work[2],2,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)v)));
277: *dp = work[2]; *nm = work[3];
278: #else
279: *dp = dp0; *nm = nm0;
280: #endif
281: PetscFunctionReturn(PETSC_SUCCESS);
282: }