Actual source code: dvdutils.c

slepc-3.17.0 2022-03-31
Report Typos and Errors
  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: */
 10: /*
 11:    SLEPc eigensolver: "davidson"

 13:    Some utils
 14: */

 16: #include "davidson.h"

 18: typedef struct {
 19:   PC pc;
 20: } dvdPCWrapper;

 22: /*
 23:   Configure the harmonics.
 24:   switch (mode) {
 25:   DVD_HARM_RR:    harmonic RR
 26:   DVD_HARM_RRR:   relative harmonic RR
 27:   DVD_HARM_REIGS: rightmost eigenvalues
 28:   DVD_HARM_LEIGS: largest eigenvalues
 29:   }
 30:   fixedTarged, if true use the target instead of the best eigenvalue
 31:   target, the fixed target to be used
 32: */
 33: typedef struct {
 34:   PetscScalar Wa,Wb;       /* span{W} = span{Wa*AV - Wb*BV} */
 35:   PetscScalar Pa,Pb;       /* H=W'*(Pa*AV - Pb*BV), G=W'*(Wa*AV - Wb*BV) */
 36:   PetscBool   withTarget;
 37:   HarmType_t  mode;
 38: } dvdHarmonic;

 40: typedef struct {
 41:   Vec diagA, diagB;
 42: } dvdJacobiPrecond;

 44: static PetscErrorCode dvd_improvex_precond_d(dvdDashboard *d)
 45: {
 46:   dvdPCWrapper   *dvdpc = (dvdPCWrapper*)d->improvex_precond_data;

 48:   /* Free local data */
 49:   PCDestroy(&dvdpc->pc);
 50:   PetscFree(d->improvex_precond_data);
 51:   PetscFunctionReturn(0);
 52: }

 54: static PetscErrorCode dvd_static_precond_PC_0(dvdDashboard *d,PetscInt i,Vec x,Vec Px)
 55: {
 56:   dvdPCWrapper   *dvdpc = (dvdPCWrapper*)d->improvex_precond_data;

 58:   PCApply(dvdpc->pc,x,Px);
 59:   PetscFunctionReturn(0);
 60: }

 62: /*
 63:   Create a trivial preconditioner
 64: */
 65: static PetscErrorCode dvd_precond_none(dvdDashboard *d,PetscInt i,Vec x,Vec Px)
 66: {
 67:   VecCopy(x,Px);
 68:   PetscFunctionReturn(0);
 69: }

 71: /*
 72:   Create a static preconditioner from a PC
 73: */
 74: PetscErrorCode dvd_static_precond_PC(dvdDashboard *d,dvdBlackboard *b,PC pc)
 75: {
 76:   dvdPCWrapper   *dvdpc;
 77:   Mat            P;
 78:   PetscBool      t0,t1,t2;

 80:   /* Setup the step */
 81:   if (b->state >= DVD_STATE_CONF) {
 82:     /* If the preconditioner is valid */
 83:     if (pc) {
 84:       PetscNewLog(d->eps,&dvdpc);
 85:       dvdpc->pc = pc;
 86:       PetscObjectReference((PetscObject)pc);
 87:       d->improvex_precond_data = dvdpc;
 88:       d->improvex_precond = dvd_static_precond_PC_0;

 90:       /* PC saves the matrix associated with the linear system, and it has to
 91:          be initialize to a valid matrix */
 92:       PCGetOperatorsSet(pc,NULL,&t0);
 93:       PetscObjectTypeCompare((PetscObject)pc,PCNONE,&t1);
 94:       PetscObjectTypeCompare((PetscObject)pc,PCSHELL,&t2);
 95:       if (t0 && !t1) {
 96:         PCGetOperators(pc,NULL,&P);
 97:         PetscObjectReference((PetscObject)P);
 98:         PCSetOperators(pc,P,P);
 99:         PCSetReusePreconditioner(pc,PETSC_TRUE);
100:         MatDestroy(&P);
101:       } else if (t2) {
102:         PCSetOperators(pc,d->A,d->A);
103:         PCSetReusePreconditioner(pc,PETSC_TRUE);
104:       } else {
105:         d->improvex_precond = dvd_precond_none;
106:       }

108:       EPSDavidsonFLAdd(&d->destroyList,dvd_improvex_precond_d);

110:     /* Else, use no preconditioner */
111:     } else d->improvex_precond = dvd_precond_none;
112:   }
113:   PetscFunctionReturn(0);
114: }

116: static PetscErrorCode dvd_harm_d(dvdDashboard *d)
117: {
118:   /* Free local data */
119:   PetscFree(d->calcpairs_W_data);
120:   PetscFunctionReturn(0);
121: }

123: static PetscErrorCode dvd_harm_transf(dvdHarmonic *dvdh,PetscScalar t)
124: {
125:   switch (dvdh->mode) {
126:   case DVD_HARM_RR:    /* harmonic RR */
127:     dvdh->Wa = 1.0; dvdh->Wb = t;   dvdh->Pa = 0.0; dvdh->Pb = -1.0;
128:     break;
129:   case DVD_HARM_RRR:   /* relative harmonic RR */
130:     dvdh->Wa = 1.0; dvdh->Wb = t;   dvdh->Pa = 1.0; dvdh->Pb = 0.0;
131:     break;
132:   case DVD_HARM_REIGS: /* rightmost eigenvalues */
133:     dvdh->Wa = 1.0; dvdh->Wb = t;   dvdh->Pa = 1.0; dvdh->Pb = -PetscConj(t);
134:     break;
135:   case DVD_HARM_LEIGS: /* largest eigenvalues */
136:     dvdh->Wa = 0.0; dvdh->Wb = 1.0; dvdh->Pa = 1.0; dvdh->Pb = 0.0;
137:     break;
138:   case DVD_HARM_NONE:
139:   default:
140:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Harmonic type not supported");
141:   }

143:   /* Check the transformation does not change the sign of the imaginary part */
144: #if !defined(PETSC_USE_COMPLEX)
145:   if (dvdh->Pb*dvdh->Wa - dvdh->Wb*dvdh->Pa < 0.0) {
146:     dvdh->Pa *= -1.0;
147:     dvdh->Pb *= -1.0;
148:   }
149: #endif
150:   PetscFunctionReturn(0);
151: }

153: static PetscErrorCode dvd_harm_updateW(dvdDashboard *d)
154: {
155:   dvdHarmonic    *data = (dvdHarmonic*)d->calcpairs_W_data;
156:   PetscInt       l,k;
157:   BV             BX = d->BX?d->BX:d->eps->V;

159:   /* Update the target if it is necessary */
160:   if (!data->withTarget) dvd_harm_transf(data,d->eigr[0]);

162:   /* W(i) <- Wa*AV(i) - Wb*BV(i) */
163:   BVGetActiveColumns(d->eps->V,&l,&k);
164:   PetscAssert(k==l+d->V_new_s,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Consistency broken");
165:   BVSetActiveColumns(d->W,l+d->V_new_s,l+d->V_new_e);
166:   BVSetActiveColumns(d->AX,l+d->V_new_s,l+d->V_new_e);
167:   BVSetActiveColumns(BX,l+d->V_new_s,l+d->V_new_e);
168:   BVCopy(d->AX,d->W);
169:   BVScale(d->W,data->Wa);
170:   BVMult(d->W,-data->Wb,1.0,BX,NULL);
171:   BVSetActiveColumns(d->W,l,k);
172:   BVSetActiveColumns(d->AX,l,k);
173:   BVSetActiveColumns(BX,l,k);
174:   PetscFunctionReturn(0);
175: }

177: static PetscErrorCode dvd_harm_proj(dvdDashboard *d)
178: {
179:   dvdHarmonic    *data = (dvdHarmonic*)d->calcpairs_W_data;
180:   PetscInt       i,j,l0,l,k,ld;
181:   PetscScalar    h,g,*H,*G;

183:   BVGetActiveColumns(d->eps->V,&l0,&k);
184:   l = l0 + d->V_new_s;
185:   k = l0 + d->V_new_e;
186:   MatGetSize(d->H,&ld,NULL);
187:   MatDenseGetArray(d->H,&H);
188:   MatDenseGetArray(d->G,&G);
189:   /* [H G] <- [Pa*H - Pb*G, Wa*H - Wb*G] */
190:   /* Right part */
191:   for (i=l;i<k;i++) {
192:     for (j=l0;j<k;j++) {
193:       h = H[ld*i+j];
194:       g = G[ld*i+j];
195:       H[ld*i+j] = data->Pa*h - data->Pb*g;
196:       G[ld*i+j] = data->Wa*h - data->Wb*g;
197:     }
198:   }
199:   /* Left part */
200:   for (i=l0;i<l;i++) {
201:     for (j=l;j<k;j++) {
202:       h = H[ld*i+j];
203:       g = G[ld*i+j];
204:       H[ld*i+j] = data->Pa*h - data->Pb*g;
205:       G[ld*i+j] = data->Wa*h - data->Wb*g;
206:     }
207:   }
208:   MatDenseRestoreArray(d->H,&H);
209:   MatDenseRestoreArray(d->G,&G);
210:   PetscFunctionReturn(0);
211: }

213: PetscErrorCode dvd_harm_updateproj(dvdDashboard *d)
214: {
215:   dvdHarmonic    *data = (dvdHarmonic*)d->calcpairs_W_data;
216:   PetscInt       i,j,l,k,ld;
217:   PetscScalar    h,g,*H,*G;

219:   BVGetActiveColumns(d->eps->V,&l,&k);
220:   k = l + d->V_tra_s;
221:   MatGetSize(d->H,&ld,NULL);
222:   MatDenseGetArray(d->H,&H);
223:   MatDenseGetArray(d->G,&G);
224:   /* [H G] <- [Pa*H - Pb*G, Wa*H - Wb*G] */
225:   /* Right part */
226:   for (i=l;i<k;i++) {
227:     for (j=0;j<l;j++) {
228:       h = H[ld*i+j];
229:       g = G[ld*i+j];
230:       H[ld*i+j] = data->Pa*h - data->Pb*g;
231:       G[ld*i+j] = data->Wa*h - data->Wb*g;
232:     }
233:   }
234:   /* Lower triangular part */
235:   for (i=0;i<l;i++) {
236:     for (j=l;j<k;j++) {
237:       h = H[ld*i+j];
238:       g = G[ld*i+j];
239:       H[ld*i+j] = data->Pa*h - data->Pb*g;
240:       G[ld*i+j] = data->Wa*h - data->Wb*g;
241:     }
242:   }
243:   MatDenseRestoreArray(d->H,&H);
244:   MatDenseRestoreArray(d->G,&G);
245:   PetscFunctionReturn(0);
246: }

248: static PetscErrorCode dvd_harm_backtrans(dvdHarmonic *data,PetscScalar *ar,PetscScalar *ai)
249: {
250:   PetscScalar xr;
251: #if !defined(PETSC_USE_COMPLEX)
252:   PetscScalar xi, k;
253: #endif

255:   xr = *ar;
256: #if !defined(PETSC_USE_COMPLEX)
257:   xi = *ai;
258:   if (PetscUnlikely(xi != 0.0)) {
259:     k = (data->Pa - data->Wa*xr)*(data->Pa - data->Wa*xr) + data->Wa*data->Wa*xi*xi;
260:     *ar = (data->Pb*data->Pa - (data->Pb*data->Wa + data->Wb*data->Pa)*xr + data->Wb*data->Wa*(xr*xr + xi*xi))/k;
261:     *ai = (data->Pb*data->Wa - data->Wb*data->Pa)*xi/k;
262:   } else
263: #endif
264:     *ar = (data->Pb - data->Wb*xr) / (data->Pa - data->Wa*xr);
265:   PetscFunctionReturn(0);
266: }

268: static PetscErrorCode dvd_harm_eig_backtrans(dvdDashboard *d,PetscScalar ar,PetscScalar ai,PetscScalar *br,PetscScalar *bi)
269: {
270:   dvdHarmonic    *data = (dvdHarmonic*)d->calcpairs_W_data;

272:   dvd_harm_backtrans(data,&ar,&ai);
273:   *br = ar;
274:   *bi = ai;
275:   PetscFunctionReturn(0);
276: }

278: static PetscErrorCode dvd_harm_eigs_trans(dvdDashboard *d)
279: {
280:   dvdHarmonic    *data = (dvdHarmonic*)d->calcpairs_W_data;
281:   PetscInt       i,l,k;

283:   BVGetActiveColumns(d->eps->V,&l,&k);
284:   for (i=0;i<k-l;i++) dvd_harm_backtrans(data,&d->eigr[i],&d->eigi[i]);
285:   PetscFunctionReturn(0);
286: }

288: PetscErrorCode dvd_harm_conf(dvdDashboard *d,dvdBlackboard *b,HarmType_t mode,PetscBool fixedTarget,PetscScalar t)
289: {
290:   dvdHarmonic    *dvdh;

292:   /* Set the problem to GNHEP:
293:      d->G maybe is upper triangular due to biorthogonality of V and W */
294:   d->sEP = d->sA = d->sB = 0;

296:   /* Setup the step */
297:   if (b->state >= DVD_STATE_CONF) {
298:     PetscNewLog(d->eps,&dvdh);
299:     dvdh->withTarget = fixedTarget;
300:     dvdh->mode = mode;
301:     if (fixedTarget) dvd_harm_transf(dvdh, t);
302:     d->calcpairs_W_data = dvdh;
303:     d->calcpairs_W = dvd_harm_updateW;
304:     d->calcpairs_proj_trans = dvd_harm_proj;
305:     d->calcpairs_eigs_trans = dvd_harm_eigs_trans;
306:     d->calcpairs_eig_backtrans = dvd_harm_eig_backtrans;

308:     EPSDavidsonFLAdd(&d->destroyList,dvd_harm_d);
309:   }
310:   PetscFunctionReturn(0);
311: }