Actual source code: dvdutils.c
slepc-3.18.0 2022-10-01
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: return 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: return 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: return 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: PetscNew(&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: return 0;
114: }
116: static PetscErrorCode dvd_harm_d(dvdDashboard *d)
117: {
118: /* Free local data */
119: PetscFree(d->calcpairs_W_data);
120: return 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: return 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: return 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: return 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: return 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: return 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: return 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: return 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: PetscNew(&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: return 0;
311: }