Line data    Source code

       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
       5             : 
       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"
      12             : 
      13             :    Step: improve the eigenvectors X with GD2
      14             : */
      15             : 
      16             : #include "davidson.h"
      17             : 
      18             : typedef struct {
      19             :   PetscInt size_X;
      20             : } dvdImprovex_gd2;
      21             : 
      22          17 : static PetscErrorCode dvd_improvex_gd2_d(dvdDashboard *d)
      23             : {
      24          17 :   dvdImprovex_gd2 *data = (dvdImprovex_gd2*)d->improveX_data;
      25             : 
      26          17 :   PetscFunctionBegin;
      27             :   /* Free local data and objects */
      28          17 :   PetscCall(PetscFree(data));
      29          17 :   PetscFunctionReturn(PETSC_SUCCESS);
      30             : }
      31             : 
      32        2274 : static PetscErrorCode dvd_improvex_gd2_gen(dvdDashboard *d,PetscInt r_s,PetscInt r_e,PetscInt *size_D)
      33             : {
      34        2274 :   dvdImprovex_gd2 *data = (dvdImprovex_gd2*)d->improveX_data;
      35        2274 :   PetscInt        i,j,n,s,ld,lv,kv,max_size_D;
      36        2274 :   PetscInt        oldnpreconv = d->npreconv;
      37        2274 :   PetscScalar     *pX,*b;
      38        2274 :   Vec             *Ax,*Bx,v,*x;
      39        2274 :   Mat             M;
      40        2274 :   BV              X;
      41             : 
      42        2274 :   PetscFunctionBegin;
      43             :   /* Compute the number of pairs to improve */
      44        2274 :   PetscCall(BVGetActiveColumns(d->eps->V,&lv,&kv));
      45        2274 :   max_size_D = d->eps->ncv-kv;
      46        2274 :   n = PetscMin(PetscMin(data->size_X*2,max_size_D),(r_e-r_s)*2)/2;
      47             : 
      48             :   /* Quick exit */
      49        2274 :   if (max_size_D == 0 || r_e-r_s <= 0 || n == 0) {
      50           0 :    *size_D = 0;
      51           0 :     PetscFunctionReturn(PETSC_SUCCESS);
      52             :   }
      53             : 
      54        2274 :   PetscCall(BVDuplicateResize(d->eps->V,4,&X));
      55        2274 :   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,4,2,NULL,&M));
      56             : 
      57             :   /* Compute the eigenvectors of the selected pairs */
      58        4548 :   for (i=r_s;i<r_s+n; i++) PetscCall(DSVectors(d->eps->ds,DS_MAT_X,&i,NULL));
      59        2274 :   PetscCall(DSGetArray(d->eps->ds,DS_MAT_X,&pX));
      60        2274 :   PetscCall(DSGetLeadingDimension(d->eps->ds,&ld));
      61             : 
      62        2274 :   PetscCall(SlepcVecPoolGetVecs(d->auxV,n,&Ax));
      63        2274 :   PetscCall(SlepcVecPoolGetVecs(d->auxV,n,&Bx));
      64             : 
      65             :   /* Bx <- B*X(i) */
      66        2274 :   if (d->BX) {
      67             :     /* Compute the norms of the eigenvectors */
      68         146 :     if (d->correctXnorm) PetscCall(dvd_improvex_compute_X(d,r_s,r_s+n,Bx,pX,ld));
      69             :     else {
      70           0 :       for (i=0;i<n;i++) d->nX[r_s+i] = 1.0;
      71             :     }
      72         292 :     for (i=0;i<n;i++) PetscCall(BVMultVec(d->BX,1.0,0.0,Bx[i],&pX[ld*(r_s+i)]));
      73        2128 :   } else if (d->B) {
      74           0 :     PetscCall(SlepcVecPoolGetVecs(d->auxV,1,&x));
      75           0 :     for (i=0;i<n;i++) {
      76             :       /* auxV(0) <- X(i) */
      77           0 :       PetscCall(dvd_improvex_compute_X(d,r_s+i,r_s+i+1,x,pX,ld));
      78             :       /* Bx(i) <- B*auxV(0) */
      79           0 :       PetscCall(MatMult(d->B,x[0],Bx[i]));
      80             :     }
      81           0 :     PetscCall(SlepcVecPoolRestoreVecs(d->auxV,1,&x));
      82             :   } else {
      83             :     /* Bx <- X */
      84        2128 :     PetscCall(dvd_improvex_compute_X(d,r_s,r_s+n,Bx,pX,ld));
      85             :   }
      86             : 
      87             :   /* Ax <- A*X(i) */
      88        4548 :   for (i=0;i<n;i++) PetscCall(BVMultVec(d->AX,1.0,0.0,Ax[i],&pX[ld*(i+r_s)]));
      89             : 
      90        2274 :   PetscCall(DSRestoreArray(d->eps->ds,DS_MAT_X,&pX));
      91             : 
      92        4493 :   for (i=0,s=0;i<n;i+=s) {
      93             : #if !defined(PETSC_USE_COMPLEX)
      94        2274 :     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           0 :       PetscCall(MatDenseGetArrayWrite(M,&b));
     100           0 :       b[0] = b[5] = 1.0/d->nX[r_s+i];
     101           0 :       b[2] = b[7] = -d->eigr[r_s+i]/d->nX[r_s+i];
     102           0 :       b[6] = -(b[3] = d->eigi[r_s+i]/d->nX[r_s+i]);
     103           0 :       b[1] = b[4] = 0.0;
     104           0 :       PetscCall(MatDenseRestoreArrayWrite(M,&b));
     105           0 :       PetscCall(BVInsertVec(X,0,Ax[i]));
     106           0 :       PetscCall(BVInsertVec(X,1,Ax[i+1]));
     107           0 :       PetscCall(BVInsertVec(X,2,Bx[i]));
     108           0 :       PetscCall(BVInsertVec(X,3,Bx[i+1]));
     109           0 :       PetscCall(BVSetActiveColumns(X,0,4));
     110           0 :       PetscCall(BVMultInPlace(X,M,0,2));
     111           0 :       PetscCall(BVCopyVec(X,0,Ax[i]));
     112           0 :       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        2274 :       PetscCall(MatDenseGetArrayWrite(M,&b));
     120        2274 :       b[0] = 1.0/d->nX[r_s+i];
     121        2274 :       b[1] = -d->eigr[r_s+i]/d->nX[r_s+i];
     122        2274 :       b[4] = PetscConj(d->eigr[r_s+i]/d->nX[r_s+i]);
     123        2274 :       b[5] = 1.0/d->nX[r_s+i];
     124        2274 :       PetscCall(MatDenseRestoreArrayWrite(M,&b));
     125        2274 :       PetscCall(BVInsertVec(X,0,Ax[i]));
     126        2274 :       PetscCall(BVInsertVec(X,1,Bx[i]));
     127        2274 :       PetscCall(BVSetActiveColumns(X,0,2));
     128        2274 :       PetscCall(BVMultInPlace(X,M,0,2));
     129        2274 :       PetscCall(BVCopyVec(X,0,Ax[i]));
     130        2274 :       PetscCall(BVCopyVec(X,1,Bx[i]));
     131             :       s = 1;
     132             :     }
     133        4548 :     for (j=0;j<s;j++) d->nX[r_s+i+j] = 1.0;
     134             : 
     135             :     /* Ax = R <- P*(Ax - eig_i*Bx) */
     136        2274 :     PetscCall(d->calcpairs_proj_res(d,r_s+i,r_s+i+s,&Ax[i]));
     137             : 
     138             :     /* Check if the first eigenpairs are converged */
     139        2274 :     if (i == 0) {
     140        2274 :       PetscCall(d->preTestConv(d,0,r_s+s,r_s+s,&d->npreconv));
     141        2274 :       if (d->npreconv > oldnpreconv) break;
     142             :     }
     143             :   }
     144             : 
     145             :   /* D <- K*[Ax Bx] */
     146        2274 :   if (d->npreconv <= oldnpreconv) {
     147        4438 :     for (i=0;i<n;i++) {
     148        2219 :       PetscCall(BVGetColumn(d->eps->V,kv+i,&v));
     149        2219 :       PetscCall(d->improvex_precond(d,r_s+i,Ax[i],v));
     150        2219 :       PetscCall(BVRestoreColumn(d->eps->V,kv+i,&v));
     151             :     }
     152        4438 :     for (i=n;i<n*2;i++) {
     153        2219 :       PetscCall(BVGetColumn(d->eps->V,kv+i,&v));
     154        2219 :       PetscCall(d->improvex_precond(d,r_s+i-n,Bx[i-n],v));
     155        2219 :       PetscCall(BVRestoreColumn(d->eps->V,kv+i,&v));
     156             :     }
     157        2219 :     *size_D = 2*n;
     158             : #if !defined(PETSC_USE_COMPLEX)
     159        2219 :     if (d->eigi[r_s] != 0.0) {
     160             :       s = 4;
     161             :     } else
     162             : #endif
     163             :     {
     164        2219 :       s = 2;
     165             :     }
     166             :     /* Prevent that short vectors are discarded in the orthogonalization */
     167        6657 :     for (i=0; i<s && i<*size_D; i++) {
     168        4438 :       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          55 :   } else *size_D = 0;
     171             : 
     172        2274 :   PetscCall(SlepcVecPoolRestoreVecs(d->auxV,n,&Bx));
     173        2274 :   PetscCall(SlepcVecPoolRestoreVecs(d->auxV,n,&Ax));
     174        2274 :   PetscCall(BVDestroy(&X));
     175        2274 :   PetscCall(MatDestroy(&M));
     176        2274 :   PetscFunctionReturn(PETSC_SUCCESS);
     177             : }
     178             : 
     179          51 : PetscErrorCode dvd_improvex_gd2(dvdDashboard *d,dvdBlackboard *b,KSP ksp,PetscInt max_bs)
     180             : {
     181          51 :   dvdImprovex_gd2 *data;
     182          51 :   PC              pc;
     183             : 
     184          51 :   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          51 :   if (!DVD_IS(d->sEP, DVD_EP_HERMITIAN)) {
     190          12 :     max_bs++;
     191          12 :     b->max_size_P = PetscMax(b->max_size_P,2);
     192             :   } else
     193             : #endif
     194             :   {
     195          39 :     b->max_size_P = PetscMax(b->max_size_P,1);
     196             :   }
     197          51 :   b->max_size_X = PetscMax(b->max_size_X,max_bs);
     198             : 
     199             :   /* Setup the preconditioner */
     200          51 :   if (ksp) {
     201          51 :     PetscCall(KSPGetPC(ksp,&pc));
     202          51 :     PetscCall(dvd_static_precond_PC(d,b,pc));
     203           0 :   } else PetscCall(dvd_static_precond_PC(d,b,NULL));
     204             : 
     205             :   /* Setup the step */
     206          51 :   if (b->state >= DVD_STATE_CONF) {
     207          17 :     PetscCall(PetscNew(&data));
     208          17 :     d->improveX_data = data;
     209          17 :     data->size_X = b->max_size_X;
     210          17 :     d->improveX = dvd_improvex_gd2_gen;
     211             : 
     212          17 :     PetscCall(EPSDavidsonFLAdd(&d->destroyList,dvd_improvex_gd2_d));
     213             :   }
     214          51 :   PetscFunctionReturn(PETSC_SUCCESS);
     215             : }