Actual source code: dvdutils.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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: {
47: dvdPCWrapper *dvdpc = (dvdPCWrapper*)d->improvex_precond_data;
50: /* Free local data */
51: PCDestroy(&dvdpc->pc);
52: PetscFree(d->improvex_precond_data);
53: return(0);
54: }
56: static PetscErrorCode dvd_static_precond_PC_0(dvdDashboard *d,PetscInt i,Vec x,Vec Px)
57: {
59: dvdPCWrapper *dvdpc = (dvdPCWrapper*)d->improvex_precond_data;
62: PCApply(dvdpc->pc,x,Px);
63: return(0);
64: }
66: /*
67: Create a trivial preconditioner
68: */
69: static PetscErrorCode dvd_precond_none(dvdDashboard *d,PetscInt i,Vec x,Vec Px)
70: {
74: VecCopy(x,Px);
75: return(0);
76: }
78: /*
79: Create a static preconditioner from a PC
80: */
81: PetscErrorCode dvd_static_precond_PC(dvdDashboard *d,dvdBlackboard *b,PC pc)
82: {
84: dvdPCWrapper *dvdpc;
85: Mat P;
86: PetscBool t0,t1,t2;
89: /* Setup the step */
90: if (b->state >= DVD_STATE_CONF) {
91: /* If the preconditioner is valid */
92: if (pc) {
93: PetscNewLog(d->eps,&dvdpc);
94: dvdpc->pc = pc;
95: PetscObjectReference((PetscObject)pc);
96: d->improvex_precond_data = dvdpc;
97: d->improvex_precond = dvd_static_precond_PC_0;
99: /* PC saves the matrix associated with the linear system, and it has to
100: be initialize to a valid matrix */
101: PCGetOperatorsSet(pc,NULL,&t0);
102: PetscObjectTypeCompare((PetscObject)pc,PCNONE,&t1);
103: PetscObjectTypeCompare((PetscObject)pc,PCSHELL,&t2);
104: if (t0 && !t1) {
105: PCGetOperators(pc,NULL,&P);
106: PetscObjectReference((PetscObject)P);
107: PCSetOperators(pc,P,P);
108: PCSetReusePreconditioner(pc,PETSC_TRUE);
109: MatDestroy(&P);
110: } else if (t2) {
111: PCSetOperators(pc,d->A,d->A);
112: PCSetReusePreconditioner(pc,PETSC_TRUE);
113: } else {
114: d->improvex_precond = dvd_precond_none;
115: }
117: EPSDavidsonFLAdd(&d->destroyList,dvd_improvex_precond_d);
119: /* Else, use no preconditioner */
120: } else d->improvex_precond = dvd_precond_none;
121: }
122: return(0);
123: }
125: static PetscErrorCode dvd_harm_d(dvdDashboard *d)
126: {
130: /* Free local data */
131: PetscFree(d->calcpairs_W_data);
132: return(0);
133: }
135: static PetscErrorCode dvd_harm_transf(dvdHarmonic *dvdh,PetscScalar t)
136: {
138: switch (dvdh->mode) {
139: case DVD_HARM_RR: /* harmonic RR */
140: dvdh->Wa = 1.0; dvdh->Wb = t; dvdh->Pa = 0.0; dvdh->Pb = -1.0;
141: break;
142: case DVD_HARM_RRR: /* relative harmonic RR */
143: dvdh->Wa = 1.0; dvdh->Wb = t; dvdh->Pa = 1.0; dvdh->Pb = 0.0;
144: break;
145: case DVD_HARM_REIGS: /* rightmost eigenvalues */
146: dvdh->Wa = 1.0; dvdh->Wb = t; dvdh->Pa = 1.0; dvdh->Pb = -PetscConj(t);
147: break;
148: case DVD_HARM_LEIGS: /* largest eigenvalues */
149: dvdh->Wa = 0.0; dvdh->Wb = 1.0; dvdh->Pa = 1.0; dvdh->Pb = 0.0;
150: break;
151: case DVD_HARM_NONE:
152: default:
153: SETERRQ(PETSC_COMM_SELF,1,"Harmonic type not supported");
154: }
156: /* Check the transformation does not change the sign of the imaginary part */
157: #if !defined(PETSC_USE_COMPLEX)
158: if (dvdh->Pb*dvdh->Wa - dvdh->Wb*dvdh->Pa < 0.0) {
159: dvdh->Pa *= -1.0;
160: dvdh->Pb *= -1.0;
161: }
162: #endif
163: return(0);
164: }
166: static PetscErrorCode dvd_harm_updateW(dvdDashboard *d)
167: {
168: dvdHarmonic *data = (dvdHarmonic*)d->calcpairs_W_data;
170: PetscInt l,k;
171: BV BX = d->BX?d->BX:d->eps->V;
174: /* Update the target if it is necessary */
175: if (!data->withTarget) {
176: dvd_harm_transf(data,d->eigr[0]);
177: }
179: /* W(i) <- Wa*AV(i) - Wb*BV(i) */
180: BVGetActiveColumns(d->eps->V,&l,&k);
181: if (k != l+d->V_new_s) SETERRQ(PETSC_COMM_SELF,1,"Consistency broken");
182: BVSetActiveColumns(d->W,l+d->V_new_s,l+d->V_new_e);
183: BVSetActiveColumns(d->AX,l+d->V_new_s,l+d->V_new_e);
184: BVSetActiveColumns(BX,l+d->V_new_s,l+d->V_new_e);
185: BVCopy(d->AX,d->W);
186: BVScale(d->W,data->Wa);
187: BVMult(d->W,-data->Wb,1.0,BX,NULL);
188: BVSetActiveColumns(d->W,l,k);
189: BVSetActiveColumns(d->AX,l,k);
190: BVSetActiveColumns(BX,l,k);
191: return(0);
192: }
194: static PetscErrorCode dvd_harm_proj(dvdDashboard *d)
195: {
197: dvdHarmonic *data = (dvdHarmonic*)d->calcpairs_W_data;
198: PetscInt i,j,l0,l,k,ld;
199: PetscScalar h,g,*H,*G;
202: BVGetActiveColumns(d->eps->V,&l0,&k);
203: l = l0 + d->V_new_s;
204: k = l0 + d->V_new_e;
205: MatGetSize(d->H,&ld,NULL);
206: MatDenseGetArray(d->H,&H);
207: MatDenseGetArray(d->G,&G);
208: /* [H G] <- [Pa*H - Pb*G, Wa*H - Wb*G] */
209: /* Right part */
210: for (i=l;i<k;i++) {
211: for (j=l0;j<k;j++) {
212: h = H[ld*i+j];
213: g = G[ld*i+j];
214: H[ld*i+j] = data->Pa*h - data->Pb*g;
215: G[ld*i+j] = data->Wa*h - data->Wb*g;
216: }
217: }
218: /* Left part */
219: for (i=l0;i<l;i++) {
220: for (j=l;j<k;j++) {
221: h = H[ld*i+j];
222: g = G[ld*i+j];
223: H[ld*i+j] = data->Pa*h - data->Pb*g;
224: G[ld*i+j] = data->Wa*h - data->Wb*g;
225: }
226: }
227: MatDenseRestoreArray(d->H,&H);
228: MatDenseRestoreArray(d->G,&G);
229: return(0);
230: }
232: PetscErrorCode dvd_harm_updateproj(dvdDashboard *d)
233: {
235: dvdHarmonic *data = (dvdHarmonic*)d->calcpairs_W_data;
236: PetscInt i,j,l,k,ld;
237: PetscScalar h,g,*H,*G;
240: BVGetActiveColumns(d->eps->V,&l,&k);
241: k = l + d->V_tra_s;
242: MatGetSize(d->H,&ld,NULL);
243: MatDenseGetArray(d->H,&H);
244: MatDenseGetArray(d->G,&G);
245: /* [H G] <- [Pa*H - Pb*G, Wa*H - Wb*G] */
246: /* Right part */
247: for (i=l;i<k;i++) {
248: for (j=0;j<l;j++) {
249: h = H[ld*i+j];
250: g = G[ld*i+j];
251: H[ld*i+j] = data->Pa*h - data->Pb*g;
252: G[ld*i+j] = data->Wa*h - data->Wb*g;
253: }
254: }
255: /* Lower triangular part */
256: for (i=0;i<l;i++) {
257: for (j=l;j<k;j++) {
258: h = H[ld*i+j];
259: g = G[ld*i+j];
260: H[ld*i+j] = data->Pa*h - data->Pb*g;
261: G[ld*i+j] = data->Wa*h - data->Wb*g;
262: }
263: }
264: MatDenseRestoreArray(d->H,&H);
265: MatDenseRestoreArray(d->G,&G);
266: return(0);
267: }
269: static PetscErrorCode dvd_harm_backtrans(dvdHarmonic *data,PetscScalar *ar,PetscScalar *ai)
270: {
271: PetscScalar xr;
272: #if !defined(PETSC_USE_COMPLEX)
273: PetscScalar xi, k;
274: #endif
278: xr = *ar;
279: #if !defined(PETSC_USE_COMPLEX)
281: xi = *ai;
283: if (xi != 0.0) {
284: k = (data->Pa - data->Wa*xr)*(data->Pa - data->Wa*xr) + data->Wa*data->Wa*xi*xi;
285: *ar = (data->Pb*data->Pa - (data->Pb*data->Wa + data->Wb*data->Pa)*xr + data->Wb*data->Wa*(xr*xr + xi*xi))/k;
286: *ai = (data->Pb*data->Wa - data->Wb*data->Pa)*xi/k;
287: } else
288: #endif
289: *ar = (data->Pb - data->Wb*xr) / (data->Pa - data->Wa*xr);
290: return(0);
291: }
293: static PetscErrorCode dvd_harm_eig_backtrans(dvdDashboard *d,PetscScalar ar,PetscScalar ai,PetscScalar *br,PetscScalar *bi)
294: {
295: dvdHarmonic *data = (dvdHarmonic*)d->calcpairs_W_data;
299: dvd_harm_backtrans(data,&ar,&ai);
300: *br = ar;
301: *bi = ai;
302: return(0);
303: }
305: static PetscErrorCode dvd_harm_eigs_trans(dvdDashboard *d)
306: {
307: dvdHarmonic *data = (dvdHarmonic*)d->calcpairs_W_data;
308: PetscInt i,l,k;
312: BVGetActiveColumns(d->eps->V,&l,&k);
313: for (i=0;i<k-l;i++) {
314: dvd_harm_backtrans(data,&d->eigr[i],&d->eigi[i]);
315: }
316: return(0);
317: }
319: PetscErrorCode dvd_harm_conf(dvdDashboard *d,dvdBlackboard *b,HarmType_t mode,PetscBool fixedTarget,PetscScalar t)
320: {
322: dvdHarmonic *dvdh;
325: /* Set the problem to GNHEP:
326: d->G maybe is upper triangular due to biorthogonality of V and W */
327: d->sEP = d->sA = d->sB = 0;
329: /* Setup the step */
330: if (b->state >= DVD_STATE_CONF) {
331: PetscNewLog(d->eps,&dvdh);
332: dvdh->withTarget = fixedTarget;
333: dvdh->mode = mode;
334: if (fixedTarget) dvd_harm_transf(dvdh, t);
335: d->calcpairs_W_data = dvdh;
336: d->calcpairs_W = dvd_harm_updateW;
337: d->calcpairs_proj_trans = dvd_harm_proj;
338: d->calcpairs_eigs_trans = dvd_harm_eigs_trans;
339: d->calcpairs_eig_backtrans = dvd_harm_eig_backtrans;
341: EPSDavidsonFLAdd(&d->destroyList,dvd_harm_d);
342: }
343: return(0);
344: }