Actual source code: dvdgd2.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: Step: improve the eigenvectors X with GD2
14: */
16: #include "davidson.h"
18: typedef struct {
19: PetscInt size_X;
20: } dvdImprovex_gd2;
22: static PetscErrorCode dvd_improvex_gd2_d(dvdDashboard *d)
23: {
24: dvdImprovex_gd2 *data = (dvdImprovex_gd2*)d->improveX_data;
26: PetscFunctionBegin;
27: /* Free local data and objects */
28: PetscCall(PetscFree(data));
29: PetscFunctionReturn(PETSC_SUCCESS);
30: }
32: static PetscErrorCode dvd_improvex_gd2_gen(dvdDashboard *d,PetscInt r_s,PetscInt r_e,PetscInt *size_D)
33: {
34: dvdImprovex_gd2 *data = (dvdImprovex_gd2*)d->improveX_data;
35: PetscInt i,j,n,s,ld,lv,kv,max_size_D;
36: PetscInt oldnpreconv = d->npreconv;
37: PetscScalar *pX,*b;
38: Vec *Ax,*Bx,v,*x;
39: Mat M;
40: BV X;
42: PetscFunctionBegin;
43: /* Compute the number of pairs to improve */
44: PetscCall(BVGetActiveColumns(d->eps->V,&lv,&kv));
45: max_size_D = d->eps->ncv-kv;
46: n = PetscMin(PetscMin(data->size_X*2,max_size_D),(r_e-r_s)*2)/2;
48: /* Quick exit */
49: if (max_size_D == 0 || r_e-r_s <= 0 || n == 0) {
50: *size_D = 0;
51: PetscFunctionReturn(PETSC_SUCCESS);
52: }
54: PetscCall(BVDuplicateResize(d->eps->V,4,&X));
55: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,4,2,NULL,&M));
57: /* Compute the eigenvectors of the selected pairs */
58: for (i=r_s;i<r_s+n; i++) PetscCall(DSVectors(d->eps->ds,DS_MAT_X,&i,NULL));
59: PetscCall(DSGetArray(d->eps->ds,DS_MAT_X,&pX));
60: PetscCall(DSGetLeadingDimension(d->eps->ds,&ld));
62: PetscCall(SlepcVecPoolGetVecs(d->auxV,n,&Ax));
63: PetscCall(SlepcVecPoolGetVecs(d->auxV,n,&Bx));
65: /* Bx <- B*X(i) */
66: if (d->BX) {
67: /* Compute the norms of the eigenvectors */
68: if (d->correctXnorm) PetscCall(dvd_improvex_compute_X(d,r_s,r_s+n,Bx,pX,ld));
69: else {
70: for (i=0;i<n;i++) d->nX[r_s+i] = 1.0;
71: }
72: for (i=0;i<n;i++) PetscCall(BVMultVec(d->BX,1.0,0.0,Bx[i],&pX[ld*(r_s+i)]));
73: } else if (d->B) {
74: PetscCall(SlepcVecPoolGetVecs(d->auxV,1,&x));
75: for (i=0;i<n;i++) {
76: /* auxV(0) <- X(i) */
77: PetscCall(dvd_improvex_compute_X(d,r_s+i,r_s+i+1,x,pX,ld));
78: /* Bx(i) <- B*auxV(0) */
79: PetscCall(MatMult(d->B,x[0],Bx[i]));
80: }
81: PetscCall(SlepcVecPoolRestoreVecs(d->auxV,1,&x));
82: } else {
83: /* Bx <- X */
84: PetscCall(dvd_improvex_compute_X(d,r_s,r_s+n,Bx,pX,ld));
85: }
87: /* Ax <- A*X(i) */
88: for (i=0;i<n;i++) PetscCall(BVMultVec(d->AX,1.0,0.0,Ax[i],&pX[ld*(i+r_s)]));
90: PetscCall(DSRestoreArray(d->eps->ds,DS_MAT_X,&pX));
92: for (i=0,s=0;i<n;i+=s) {
93: #if !defined(PETSC_USE_COMPLEX)
94: if (d->eigi[r_s+i] != 0.0 && i+2<=n) {
95: /* [Ax_i Ax_i+1 Bx_i Bx_i+1]*= [ 1 0
96: 0 1
97: -eigr_i -eigi_i
98: eigi_i -eigr_i] */
99: PetscCall(MatDenseGetArrayWrite(M,&b));
100: b[0] = b[5] = 1.0/d->nX[r_s+i];
101: b[2] = b[7] = -d->eigr[r_s+i]/d->nX[r_s+i];
102: b[6] = -(b[3] = d->eigi[r_s+i]/d->nX[r_s+i]);
103: b[1] = b[4] = 0.0;
104: PetscCall(MatDenseRestoreArrayWrite(M,&b));
105: PetscCall(BVInsertVec(X,0,Ax[i]));
106: PetscCall(BVInsertVec(X,1,Ax[i+1]));
107: PetscCall(BVInsertVec(X,2,Bx[i]));
108: PetscCall(BVInsertVec(X,3,Bx[i+1]));
109: PetscCall(BVSetActiveColumns(X,0,4));
110: PetscCall(BVMultInPlace(X,M,0,2));
111: PetscCall(BVCopyVec(X,0,Ax[i]));
112: PetscCall(BVCopyVec(X,1,Ax[i+1]));
113: s = 2;
114: } else
115: #endif
116: {
117: /* [Ax_i Bx_i]*= [ 1/nX_i conj(eig_i/nX_i)
118: -eig_i/nX_i 1/nX_i ] */
119: PetscCall(MatDenseGetArrayWrite(M,&b));
120: b[0] = 1.0/d->nX[r_s+i];
121: b[1] = -d->eigr[r_s+i]/d->nX[r_s+i];
122: b[4] = PetscConj(d->eigr[r_s+i]/d->nX[r_s+i]);
123: b[5] = 1.0/d->nX[r_s+i];
124: PetscCall(MatDenseRestoreArrayWrite(M,&b));
125: PetscCall(BVInsertVec(X,0,Ax[i]));
126: PetscCall(BVInsertVec(X,1,Bx[i]));
127: PetscCall(BVSetActiveColumns(X,0,2));
128: PetscCall(BVMultInPlace(X,M,0,2));
129: PetscCall(BVCopyVec(X,0,Ax[i]));
130: PetscCall(BVCopyVec(X,1,Bx[i]));
131: s = 1;
132: }
133: for (j=0;j<s;j++) d->nX[r_s+i+j] = 1.0;
135: /* Ax = R <- P*(Ax - eig_i*Bx) */
136: PetscCall(d->calcpairs_proj_res(d,r_s+i,r_s+i+s,&Ax[i]));
138: /* Check if the first eigenpairs are converged */
139: if (i == 0) {
140: PetscCall(d->preTestConv(d,0,r_s+s,r_s+s,&d->npreconv));
141: if (d->npreconv > oldnpreconv) break;
142: }
143: }
145: /* D <- K*[Ax Bx] */
146: if (d->npreconv <= oldnpreconv) {
147: for (i=0;i<n;i++) {
148: PetscCall(BVGetColumn(d->eps->V,kv+i,&v));
149: PetscCall(d->improvex_precond(d,r_s+i,Ax[i],v));
150: PetscCall(BVRestoreColumn(d->eps->V,kv+i,&v));
151: }
152: for (i=n;i<n*2;i++) {
153: PetscCall(BVGetColumn(d->eps->V,kv+i,&v));
154: PetscCall(d->improvex_precond(d,r_s+i-n,Bx[i-n],v));
155: PetscCall(BVRestoreColumn(d->eps->V,kv+i,&v));
156: }
157: *size_D = 2*n;
158: #if !defined(PETSC_USE_COMPLEX)
159: if (d->eigi[r_s] != 0.0) {
160: s = 4;
161: } else
162: #endif
163: {
164: s = 2;
165: }
166: /* Prevent that short vectors are discarded in the orthogonalization */
167: for (i=0; i<s && i<*size_D; i++) {
168: if (d->eps->errest[d->nconv+r_s+i] > PETSC_MACHINE_EPSILON && d->eps->errest[d->nconv+r_s+i] < PETSC_MAX_REAL) PetscCall(BVScaleColumn(d->eps->V,i+kv,1.0/d->eps->errest[d->nconv+r_s+i]));
169: }
170: } else *size_D = 0;
172: PetscCall(SlepcVecPoolRestoreVecs(d->auxV,n,&Bx));
173: PetscCall(SlepcVecPoolRestoreVecs(d->auxV,n,&Ax));
174: PetscCall(BVDestroy(&X));
175: PetscCall(MatDestroy(&M));
176: PetscFunctionReturn(PETSC_SUCCESS);
177: }
179: PetscErrorCode dvd_improvex_gd2(dvdDashboard *d,dvdBlackboard *b,KSP ksp,PetscInt max_bs)
180: {
181: dvdImprovex_gd2 *data;
182: PC pc;
184: PetscFunctionBegin;
185: /* Setting configuration constrains */
186: /* If the arithmetic is real and the problem is not Hermitian, then
187: the block size is incremented in one */
188: #if !defined(PETSC_USE_COMPLEX)
189: if (!DVD_IS(d->sEP, DVD_EP_HERMITIAN)) {
190: max_bs++;
191: b->max_size_P = PetscMax(b->max_size_P,2);
192: } else
193: #endif
194: {
195: b->max_size_P = PetscMax(b->max_size_P,1);
196: }
197: b->max_size_X = PetscMax(b->max_size_X,max_bs);
199: /* Setup the preconditioner */
200: if (ksp) {
201: PetscCall(KSPGetPC(ksp,&pc));
202: PetscCall(dvd_static_precond_PC(d,b,pc));
203: } else PetscCall(dvd_static_precond_PC(d,b,NULL));
205: /* Setup the step */
206: if (b->state >= DVD_STATE_CONF) {
207: PetscCall(PetscNew(&data));
208: d->improveX_data = data;
209: data->size_X = b->max_size_X;
210: d->improveX = dvd_improvex_gd2_gen;
212: PetscCall(EPSDavidsonFLAdd(&d->destroyList,dvd_improvex_gd2_d));
213: }
214: PetscFunctionReturn(PETSC_SUCCESS);
215: }