Actual source code: dvdutils.c
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: */
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: }