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