Actual source code: dvdutils.c

slepc-3.21.1 2024-04-26
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:   PetscFunctionBegin;
 49:   /* Free local data */
 50:   PetscCall(PCDestroy(&dvdpc->pc));
 51:   PetscCall(PetscFree(d->improvex_precond_data));
 52:   PetscFunctionReturn(PETSC_SUCCESS);
 53: }

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

 59:   PetscFunctionBegin;
 60:   PetscCall(PCApply(dvdpc->pc,x,Px));
 61:   PetscFunctionReturn(PETSC_SUCCESS);
 62: }

 64: /*
 65:   Create a trivial preconditioner
 66: */
 67: static PetscErrorCode dvd_precond_none(dvdDashboard *d,PetscInt i,Vec x,Vec Px)
 68: {
 69:   PetscFunctionBegin;
 70:   PetscCall(VecCopy(x,Px));
 71:   PetscFunctionReturn(PETSC_SUCCESS);
 72: }

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

 83:   PetscFunctionBegin;
 84:   /* Setup the step */
 85:   if (b->state >= DVD_STATE_CONF) {
 86:     /* If the preconditioner is valid */
 87:     if (pc) {
 88:       PetscCall(PetscNew(&dvdpc));
 89:       dvdpc->pc = pc;
 90:       PetscCall(PetscObjectReference((PetscObject)pc));
 91:       d->improvex_precond_data = dvdpc;
 92:       d->improvex_precond = dvd_static_precond_PC_0;

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

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

114:     /* Else, use no preconditioner */
115:     } else d->improvex_precond = dvd_precond_none;
116:   }
117:   PetscFunctionReturn(PETSC_SUCCESS);
118: }

120: static PetscErrorCode dvd_harm_d(dvdDashboard *d)
121: {
122:   PetscFunctionBegin;
123:   /* Free local data */
124:   PetscCall(PetscFree(d->calcpairs_W_data));
125:   PetscFunctionReturn(PETSC_SUCCESS);
126: }

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

149:   /* Check the transformation does not change the sign of the imaginary part */
150: #if !defined(PETSC_USE_COMPLEX)
151:   if (dvdh->Pb*dvdh->Wa - dvdh->Wb*dvdh->Pa < 0.0) {
152:     dvdh->Pa *= -1.0;
153:     dvdh->Pb *= -1.0;
154:   }
155: #endif
156:   PetscFunctionReturn(PETSC_SUCCESS);
157: }

159: static PetscErrorCode dvd_harm_updateW(dvdDashboard *d)
160: {
161:   dvdHarmonic    *data = (dvdHarmonic*)d->calcpairs_W_data;
162:   PetscInt       l,k;
163:   BV             BX = d->BX?d->BX:d->eps->V;

165:   PetscFunctionBegin;
166:   /* Update the target if it is necessary */
167:   if (!data->withTarget) PetscCall(dvd_harm_transf(data,d->eigr[0]));

169:   /* W(i) <- Wa*AV(i) - Wb*BV(i) */
170:   PetscCall(BVGetActiveColumns(d->eps->V,&l,&k));
171:   PetscAssert(k==l+d->V_new_s,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Consistency broken");
172:   PetscCall(BVSetActiveColumns(d->W,l+d->V_new_s,l+d->V_new_e));
173:   PetscCall(BVSetActiveColumns(d->AX,l+d->V_new_s,l+d->V_new_e));
174:   PetscCall(BVSetActiveColumns(BX,l+d->V_new_s,l+d->V_new_e));
175:   PetscCall(BVCopy(d->AX,d->W));
176:   PetscCall(BVScale(d->W,data->Wa));
177:   PetscCall(BVMult(d->W,-data->Wb,1.0,BX,NULL));
178:   PetscCall(BVSetActiveColumns(d->W,l,k));
179:   PetscCall(BVSetActiveColumns(d->AX,l,k));
180:   PetscCall(BVSetActiveColumns(BX,l,k));
181:   PetscFunctionReturn(PETSC_SUCCESS);
182: }

184: static PetscErrorCode dvd_harm_proj(dvdDashboard *d)
185: {
186:   dvdHarmonic    *data = (dvdHarmonic*)d->calcpairs_W_data;
187:   PetscInt       i,j,l0,l,k,ld;
188:   PetscScalar    h,g,*H,*G;

190:   PetscFunctionBegin;
191:   PetscCall(BVGetActiveColumns(d->eps->V,&l0,&k));
192:   l = l0 + d->V_new_s;
193:   k = l0 + d->V_new_e;
194:   PetscCall(MatGetSize(d->H,&ld,NULL));
195:   PetscCall(MatDenseGetArray(d->H,&H));
196:   PetscCall(MatDenseGetArray(d->G,&G));
197:   /* [H G] <- [Pa*H - Pb*G, Wa*H - Wb*G] */
198:   /* Right part */
199:   for (i=l;i<k;i++) {
200:     for (j=l0;j<k;j++) {
201:       h = H[ld*i+j];
202:       g = G[ld*i+j];
203:       H[ld*i+j] = data->Pa*h - data->Pb*g;
204:       G[ld*i+j] = data->Wa*h - data->Wb*g;
205:     }
206:   }
207:   /* Left part */
208:   for (i=l0;i<l;i++) {
209:     for (j=l;j<k;j++) {
210:       h = H[ld*i+j];
211:       g = G[ld*i+j];
212:       H[ld*i+j] = data->Pa*h - data->Pb*g;
213:       G[ld*i+j] = data->Wa*h - data->Wb*g;
214:     }
215:   }
216:   PetscCall(MatDenseRestoreArray(d->H,&H));
217:   PetscCall(MatDenseRestoreArray(d->G,&G));
218:   PetscFunctionReturn(PETSC_SUCCESS);
219: }

221: PetscErrorCode dvd_harm_updateproj(dvdDashboard *d)
222: {
223:   dvdHarmonic    *data = (dvdHarmonic*)d->calcpairs_W_data;
224:   PetscInt       i,j,l,k,ld;
225:   PetscScalar    h,g,*H,*G;

227:   PetscFunctionBegin;
228:   PetscCall(BVGetActiveColumns(d->eps->V,&l,&k));
229:   k = l + d->V_tra_s;
230:   PetscCall(MatGetSize(d->H,&ld,NULL));
231:   PetscCall(MatDenseGetArray(d->H,&H));
232:   PetscCall(MatDenseGetArray(d->G,&G));
233:   /* [H G] <- [Pa*H - Pb*G, Wa*H - Wb*G] */
234:   /* Right part */
235:   for (i=l;i<k;i++) {
236:     for (j=0;j<l;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:   /* Lower triangular part */
244:   for (i=0;i<l;i++) {
245:     for (j=l;j<k;j++) {
246:       h = H[ld*i+j];
247:       g = G[ld*i+j];
248:       H[ld*i+j] = data->Pa*h - data->Pb*g;
249:       G[ld*i+j] = data->Wa*h - data->Wb*g;
250:     }
251:   }
252:   PetscCall(MatDenseRestoreArray(d->H,&H));
253:   PetscCall(MatDenseRestoreArray(d->G,&G));
254:   PetscFunctionReturn(PETSC_SUCCESS);
255: }

257: static PetscErrorCode dvd_harm_backtrans(dvdHarmonic *data,PetscScalar *ar,PetscScalar *ai)
258: {
259:   PetscScalar xr;
260: #if !defined(PETSC_USE_COMPLEX)
261:   PetscScalar xi, k;
262: #endif

264:   PetscFunctionBegin;
265:   xr = *ar;
266: #if !defined(PETSC_USE_COMPLEX)
267:   xi = *ai;
268:   if (PetscUnlikely(xi != 0.0)) {
269:     k = (data->Pa - data->Wa*xr)*(data->Pa - data->Wa*xr) + data->Wa*data->Wa*xi*xi;
270:     *ar = (data->Pb*data->Pa - (data->Pb*data->Wa + data->Wb*data->Pa)*xr + data->Wb*data->Wa*(xr*xr + xi*xi))/k;
271:     *ai = (data->Pb*data->Wa - data->Wb*data->Pa)*xi/k;
272:   } else
273: #endif
274:     *ar = (data->Pb - data->Wb*xr) / (data->Pa - data->Wa*xr);
275:   PetscFunctionReturn(PETSC_SUCCESS);
276: }

278: static PetscErrorCode dvd_harm_eig_backtrans(dvdDashboard *d,PetscScalar ar,PetscScalar ai,PetscScalar *br,PetscScalar *bi)
279: {
280:   dvdHarmonic    *data = (dvdHarmonic*)d->calcpairs_W_data;

282:   PetscFunctionBegin;
283:   PetscCall(dvd_harm_backtrans(data,&ar,&ai));
284:   *br = ar;
285:   *bi = ai;
286:   PetscFunctionReturn(PETSC_SUCCESS);
287: }

289: static PetscErrorCode dvd_harm_eigs_trans(dvdDashboard *d)
290: {
291:   dvdHarmonic    *data = (dvdHarmonic*)d->calcpairs_W_data;
292:   PetscInt       i,l,k;

294:   PetscFunctionBegin;
295:   PetscCall(BVGetActiveColumns(d->eps->V,&l,&k));
296:   for (i=0;i<k-l;i++) PetscCall(dvd_harm_backtrans(data,&d->eigr[i],&d->eigi[i]));
297:   PetscFunctionReturn(PETSC_SUCCESS);
298: }

300: PetscErrorCode dvd_harm_conf(dvdDashboard *d,dvdBlackboard *b,HarmType_t mode,PetscBool fixedTarget,PetscScalar t)
301: {
302:   dvdHarmonic    *dvdh;

304:   PetscFunctionBegin;
305:   /* Set the problem to GNHEP:
306:      d->G maybe is upper triangular due to biorthogonality of V and W */
307:   d->sEP = d->sA = d->sB = 0;

309:   /* Setup the step */
310:   if (b->state >= DVD_STATE_CONF) {
311:     PetscCall(PetscNew(&dvdh));
312:     dvdh->withTarget = fixedTarget;
313:     dvdh->mode = mode;
314:     if (fixedTarget) PetscCall(dvd_harm_transf(dvdh, t));
315:     d->calcpairs_W_data = dvdh;
316:     d->calcpairs_W = dvd_harm_updateW;
317:     d->calcpairs_proj_trans = dvd_harm_proj;
318:     d->calcpairs_eigs_trans = dvd_harm_eigs_trans;
319:     d->calcpairs_eig_backtrans = dvd_harm_eig_backtrans;

321:     PetscCall(EPSDavidsonFLAdd(&d->destroyList,dvd_harm_d));
322:   }
323:   PetscFunctionReturn(PETSC_SUCCESS);
324: }