Actual source code: veccomp0.h
slepc-3.18.3 2023-03-24
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: 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: SlepcValidVecComp(a,1);
29: SlepcValidVecComp(b,2);
30: if (as->x[0]->ops->dot_local) {
31: for (i=0,sum=0.0;i<as->n->n;i++) {
32: PetscUseTypeMethod(as->x[i],dot_local,bs->x[i],&work);
33: sum += work;
34: }
35: #if defined(__WITH_MPI__)
36: work = sum;
37: MPIU_Allreduce(&work,&sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)a));
38: #endif
39: } else {
40: for (i=0,sum=0.0;i<as->n->n;i++) {
41: VecDot(as->x[i],bs->x[i],&work);
42: sum += work;
43: }
44: }
45: *z = sum;
46: return 0;
47: }
49: PetscErrorCode __SUF__(VecMDot_Comp)(Vec a,PetscInt n,const Vec b[],PetscScalar *z)
50: {
51: PetscScalar *work,*work0,*r;
52: Vec_Comp *as = (Vec_Comp*)a->data;
53: Vec *bx;
54: PetscInt i,j;
56: SlepcValidVecComp(a,1);
57: for (i=0;i<n;i++) SlepcValidVecComp(b[i],3);
59: if (as->n->n == 0) {
60: *z = 0;
61: return 0;
62: }
64: PetscMalloc2(n,&work0,n,&bx);
66: #if defined(__WITH_MPI__)
67: if (as->x[0]->ops->mdot_local) {
68: r = work0; work = z;
69: } else
70: #endif
71: {
72: r = z; work = work0;
73: }
75: /* z[i] <- a.x' * b[i].x */
76: for (i=0;i<n;i++) r[i] = 0.0;
77: for (j=0;j<as->n->n;j++) {
78: for (i=0;i<n;i++) bx[i] = ((Vec_Comp*)b[i]->data)->x[j];
79: if (as->x[0]->ops->mdot_local) PetscUseTypeMethod(as->x[j],mdot_local,n,bx,work);
80: else VecMDot(as->x[j],n,bx,work);
81: for (i=0;i<n;i++) r[i] += work[i];
82: }
84: /* If def(__WITH_MPI__) and exists mdot_local */
85: #if defined(__WITH_MPI__)
86: if (as->x[0]->ops->mdot_local) {
87: /* z[i] <- Allreduce(work[i]) */
88: MPIU_Allreduce(r,z,n,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)a));
89: }
90: #endif
92: PetscFree2(work0,bx);
93: return 0;
94: }
96: PetscErrorCode __SUF__(VecTDot_Comp)(Vec a,Vec b,PetscScalar *z)
97: {
98: PetscScalar sum = 0.0,work;
99: PetscInt i;
100: Vec_Comp *as = (Vec_Comp*)a->data,*bs = (Vec_Comp*)b->data;
102: SlepcValidVecComp(a,1);
103: SlepcValidVecComp(b,2);
104: if (as->x[0]->ops->tdot_local) {
105: for (i=0,sum=0.0;i<as->n->n;i++) {
106: PetscUseTypeMethod(as->x[i],tdot_local,bs->x[i],&work);
107: sum += work;
108: }
109: #if defined(__WITH_MPI__)
110: work = sum;
111: MPIU_Allreduce(&work,&sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)a));
112: #endif
113: } else {
114: for (i=0,sum=0.0;i<as->n->n;i++) {
115: VecTDot(as->x[i],bs->x[i],&work);
116: sum += work;
117: }
118: }
119: *z = sum;
120: return 0;
121: }
123: PetscErrorCode __SUF__(VecMTDot_Comp)(Vec a,PetscInt n,const Vec b[],PetscScalar *z)
124: {
125: PetscScalar *work,*work0,*r;
126: Vec_Comp *as = (Vec_Comp*)a->data;
127: Vec *bx;
128: PetscInt i,j;
130: SlepcValidVecComp(a,1);
131: for (i=0;i<n;i++) SlepcValidVecComp(b[i],3);
133: if (as->n->n == 0) {
134: *z = 0;
135: return 0;
136: }
138: PetscMalloc2(n,&work0,n,&bx);
140: #if defined(__WITH_MPI__)
141: if (as->x[0]->ops->mtdot_local) {
142: r = work0; work = z;
143: } else
144: #endif
145: {
146: r = z; work = work0;
147: }
149: /* z[i] <- a.x' * b[i].x */
150: for (i=0;i<n;i++) r[i] = 0.0;
151: for (j=0;j<as->n->n;j++) {
152: for (i=0;i<n;i++) bx[i] = ((Vec_Comp*)b[i]->data)->x[j];
153: if (as->x[0]->ops->mtdot_local) PetscUseTypeMethod(as->x[j],mtdot_local,n,bx,work);
154: else VecMTDot(as->x[j],n,bx,work);
155: for (i=0;i<n;i++) r[i] += work[i];
156: }
158: /* If def(__WITH_MPI__) and exists mtdot_local */
159: #if defined(__WITH_MPI__)
160: if (as->x[0]->ops->mtdot_local) {
161: /* z[i] <- Allreduce(work[i]) */
162: MPIU_Allreduce(r,z,n,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)a));
163: }
164: #endif
166: PetscFree2(work0,bx);
167: return 0;
168: }
170: PetscErrorCode __SUF__(VecNorm_Comp)(Vec a,NormType t,PetscReal *norm)
171: {
172: PetscReal work[3],s=0.0;
173: Vec_Comp *as = (Vec_Comp*)a->data;
174: PetscInt i;
176: SlepcValidVecComp(a,1);
177: /* Initialize norm */
178: switch (t) {
179: case NORM_1: case NORM_INFINITY: *norm = 0.0; break;
180: case NORM_2: case NORM_FROBENIUS: *norm = 1.0; s = 0.0; break;
181: case NORM_1_AND_2: norm[0] = 0.0; norm[1] = 1.0; s = 0.0; break;
182: }
183: for (i=0;i<as->n->n;i++) {
184: if (as->x[0]->ops->norm_local) PetscUseTypeMethod(as->x[i],norm_local,t,work);
185: else VecNorm(as->x[i],t,work);
186: /* norm+= work */
187: switch (t) {
188: case NORM_1: *norm+= *work; break;
189: case NORM_2: case NORM_FROBENIUS: AddNorm2(norm,&s,*work); break;
190: case NORM_1_AND_2: norm[0]+= work[0]; AddNorm2(&norm[1],&s,work[1]); break;
191: case NORM_INFINITY: *norm = PetscMax(*norm,*work); break;
192: }
193: }
195: /* If def(__WITH_MPI__) and exists norm_local */
196: #if defined(__WITH_MPI__)
197: if (as->x[0]->ops->norm_local) {
198: PetscReal work0[3];
199: /* norm <- Allreduce(work) */
200: switch (t) {
201: case NORM_1:
202: work[0] = *norm;
203: MPIU_Allreduce(work,norm,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)a));
204: break;
205: case NORM_2: case NORM_FROBENIUS:
206: work[0] = *norm; work[1] = s;
207: MPIU_Allreduce(work,work0,1,MPIU_NORM2,MPIU_NORM2_SUM,PetscObjectComm((PetscObject)a));
208: *norm = GetNorm2(work0[0],work0[1]);
209: break;
210: case NORM_1_AND_2:
211: work[0] = norm[0]; work[1] = norm[1]; work[2] = s;
212: MPIU_Allreduce(work,work0,1,MPIU_NORM1_AND_2,MPIU_NORM2_SUM,PetscObjectComm((PetscObject)a));
213: norm[0] = work0[0];
214: norm[1] = GetNorm2(work0[1],work0[2]);
215: break;
216: case NORM_INFINITY:
217: work[0] = *norm;
218: MPIU_Allreduce(work,norm,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)a));
219: break;
220: }
221: }
222: #else
223: /* Norm correction */
224: switch (t) {
225: case NORM_2: case NORM_FROBENIUS: *norm = GetNorm2(*norm,s); break;
226: case NORM_1_AND_2: norm[1] = GetNorm2(norm[1],s); break;
227: default: ;
228: }
229: #endif
230: return 0;
231: }
233: PetscErrorCode __SUF__(VecDotNorm2_Comp)(Vec v,Vec w,PetscScalar *dp,PetscScalar *nm)
234: {
235: PetscScalar dp0=0.0,nm0=0.0,dp1=0.0,nm1=0.0;
236: const PetscScalar *vx,*wx;
237: Vec_Comp *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data;
238: PetscInt i,n;
239: PetscBool t0,t1;
240: #if defined(__WITH_MPI__)
241: PetscScalar work[4];
242: #endif
244: /* Compute recursively the local part */
245: PetscObjectTypeCompare((PetscObject)v,VECCOMP,&t0);
246: PetscObjectTypeCompare((PetscObject)w,VECCOMP,&t1);
247: if (t0 && t1) {
248: SlepcValidVecComp(v,1);
249: SlepcValidVecComp(w,2);
250: for (i=0;i<vs->n->n;i++) {
251: VecDotNorm2_Comp_Seq(vs->x[i],ws->x[i],&dp1,&nm1);
252: dp0 += dp1;
253: nm0 += nm1;
254: }
255: } else if (!t0 && !t1) {
256: VecGetLocalSize(v,&n);
257: VecGetArrayRead(v,&vx);
258: VecGetArrayRead(w,&wx);
259: for (i=0;i<n;i++) {
260: dp0 += vx[i]*PetscConj(wx[i]);
261: nm0 += wx[i]*PetscConj(wx[i]);
262: }
263: VecRestoreArrayRead(v,&vx);
264: VecRestoreArrayRead(w,&wx);
265: } else SETERRQ(PetscObjectComm((PetscObject)v),PETSC_ERR_ARG_INCOMP,"Incompatible vector types");
267: #if defined(__WITH_MPI__)
268: /* [dp, nm] <- Allreduce([dp0, nm0]) */
269: work[0] = dp0; work[1] = nm0;
270: MPIU_Allreduce(work,&work[2],2,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)v));
271: *dp = work[2]; *nm = work[3];
272: #else
273: *dp = dp0; *nm = nm0;
274: #endif
275: return 0;
276: }