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