LCOV - code coverage report
Current view: top level - eps/impls/davidson - dvdgd2.c (source / functions) Hit Total Coverage
Test: coverage.info Lines: 88 111 79.3 %
Date: 2019-07-19 01:26:51 Functions: 3 3 100.0 %

          Line data    Source code
       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-2019, 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             :   PetscErrorCode  ierr;
      25          17 :   dvdImprovex_gd2 *data = (dvdImprovex_gd2*)d->improveX_data;
      26             : 
      27          17 :   PetscFunctionBegin;
      28             :   /* Free local data and objects */
      29          17 :   ierr = PetscFree(data);CHKERRQ(ierr);
      30          17 :   PetscFunctionReturn(0);
      31             : }
      32             : 
      33        2515 : static PetscErrorCode dvd_improvex_gd2_gen(dvdDashboard *d,PetscInt r_s,PetscInt r_e,PetscInt *size_D)
      34             : {
      35        2515 :   dvdImprovex_gd2 *data = (dvdImprovex_gd2*)d->improveX_data;
      36             :   PetscErrorCode  ierr;
      37             :   PetscInt        i,j,n,s,ld,lv,kv,max_size_D;
      38        2515 :   PetscInt        oldnpreconv = d->npreconv;
      39             :   PetscScalar     *pX,*b;
      40             :   Vec             *Ax,*Bx,v,*x;
      41             :   Mat             M;
      42             :   BV              X;
      43             : 
      44        2515 :   PetscFunctionBegin;
      45             :   /* Compute the number of pairs to improve */
      46        2515 :   ierr = BVGetActiveColumns(d->eps->V,&lv,&kv);CHKERRQ(ierr);
      47        2515 :   max_size_D = d->eps->ncv-kv;
      48        2515 :   n = PetscMin(PetscMin(data->size_X*2,max_size_D),(r_e-r_s)*2)/2;
      49             : 
      50             :   /* Quick exit */
      51        2515 :   if (max_size_D == 0 || r_e-r_s <= 0 || n == 0) {
      52           0 :    *size_D = 0;
      53           0 :     PetscFunctionReturn(0);
      54             :   }
      55             : 
      56        2515 :   ierr = BVDuplicateResize(d->eps->V,4,&X);CHKERRQ(ierr);
      57        2515 :   ierr = MatCreateSeqDense(PETSC_COMM_SELF,4,2,NULL,&M);CHKERRQ(ierr);
      58             : 
      59             :   /* Compute the eigenvectors of the selected pairs */
      60        5030 :   for (i=r_s;i<r_s+n; i++) {
      61        2515 :     ierr = DSVectors(d->eps->ds,DS_MAT_X,&i,NULL);CHKERRQ(ierr);
      62             :   }
      63        2515 :   ierr = DSGetArray(d->eps->ds,DS_MAT_X,&pX);CHKERRQ(ierr);
      64        2515 :   ierr = DSGetLeadingDimension(d->eps->ds,&ld);CHKERRQ(ierr);
      65             : 
      66        2515 :   ierr = SlepcVecPoolGetVecs(d->auxV,n,&Ax);CHKERRQ(ierr);
      67        2515 :   ierr = SlepcVecPoolGetVecs(d->auxV,n,&Bx);CHKERRQ(ierr);
      68             : 
      69             :   /* Bx <- B*X(i) */
      70        2515 :   if (d->BX) {
      71             :     /* Compute the norms of the eigenvectors */
      72         187 :     if (d->correctXnorm) {
      73         187 :       ierr = dvd_improvex_compute_X(d,r_s,r_s+n,Bx,pX,ld);CHKERRQ(ierr);
      74             :     } else {
      75           0 :       for (i=0;i<n;i++) d->nX[r_s+i] = 1.0;
      76             :     }
      77         374 :     for (i=0;i<n;i++) {
      78         187 :       ierr = BVMultVec(d->BX,1.0,0.0,Bx[i],&pX[ld*(r_s+i)]);CHKERRQ(ierr);
      79             :     }
      80        2328 :   } else if (d->B) {
      81           0 :     ierr = SlepcVecPoolGetVecs(d->auxV,1,&x);CHKERRQ(ierr);
      82           0 :     for (i=0;i<n;i++) {
      83             :       /* auxV(0) <- X(i) */
      84           0 :       ierr = dvd_improvex_compute_X(d,r_s+i,r_s+i+1,x,pX,ld);CHKERRQ(ierr);
      85             :       /* Bx(i) <- B*auxV(0) */
      86           0 :       ierr = MatMult(d->B,x[0],Bx[i]);CHKERRQ(ierr);
      87             :     }
      88           0 :     ierr = SlepcVecPoolRestoreVecs(d->auxV,1,&x);CHKERRQ(ierr);
      89             :   } else {
      90             :     /* Bx <- X */
      91        2328 :     ierr = dvd_improvex_compute_X(d,r_s,r_s+n,Bx,pX,ld);CHKERRQ(ierr);
      92             :   }
      93             : 
      94             :   /* Ax <- A*X(i) */
      95        5030 :   for (i=0;i<n;i++) {
      96        2515 :     ierr = BVMultVec(d->AX,1.0,0.0,Ax[i],&pX[ld*(i+r_s)]);CHKERRQ(ierr);
      97             :   }
      98             : 
      99        2515 :   ierr = DSRestoreArray(d->eps->ds,DS_MAT_X,&pX);CHKERRQ(ierr);
     100             : 
     101        4975 :   for (i=0,s=0;i<n;i+=s) {
     102             : #if !defined(PETSC_USE_COMPLEX)
     103        2515 :     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           0 :       ierr = MatDenseGetArray(M,&b);CHKERRQ(ierr);
     109           0 :       b[0] = b[5] = 1.0/d->nX[r_s+i];
     110           0 :       b[2] = b[7] = -d->eigr[r_s+i]/d->nX[r_s+i];
     111           0 :       b[6] = -(b[3] = d->eigi[r_s+i]/d->nX[r_s+i]);
     112           0 :       b[1] = b[4] = 0.0;
     113           0 :       ierr = MatDenseRestoreArray(M,&b);CHKERRQ(ierr);
     114           0 :       ierr = BVInsertVec(X,0,Ax[i]);CHKERRQ(ierr);
     115           0 :       ierr = BVInsertVec(X,1,Ax[i+1]);CHKERRQ(ierr);
     116           0 :       ierr = BVInsertVec(X,2,Bx[i]);CHKERRQ(ierr);
     117           0 :       ierr = BVInsertVec(X,3,Bx[i+1]);CHKERRQ(ierr);
     118           0 :       ierr = BVSetActiveColumns(X,0,4);CHKERRQ(ierr);
     119           0 :       ierr = BVMultInPlace(X,M,0,2);CHKERRQ(ierr);
     120           0 :       ierr = BVCopyVec(X,0,Ax[i]);CHKERRQ(ierr);
     121           0 :       ierr = BVCopyVec(X,1,Ax[i+1]);CHKERRQ(ierr);
     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        2515 :       ierr = MatDenseGetArray(M,&b);CHKERRQ(ierr);
     129        2515 :       b[0] = 1.0/d->nX[r_s+i];
     130        2515 :       b[1] = -d->eigr[r_s+i]/d->nX[r_s+i];
     131        2515 :       b[4] = PetscConj(d->eigr[r_s+i]/d->nX[r_s+i]);
     132        2515 :       b[5] = 1.0/d->nX[r_s+i];
     133        2515 :       ierr = MatDenseRestoreArray(M,&b);CHKERRQ(ierr);
     134        2515 :       ierr = BVInsertVec(X,0,Ax[i]);CHKERRQ(ierr);
     135        2515 :       ierr = BVInsertVec(X,1,Bx[i]);CHKERRQ(ierr);
     136        2515 :       ierr = BVSetActiveColumns(X,0,2);CHKERRQ(ierr);
     137        2515 :       ierr = BVMultInPlace(X,M,0,2);CHKERRQ(ierr);
     138        2515 :       ierr = BVCopyVec(X,0,Ax[i]);CHKERRQ(ierr);
     139        2515 :       ierr = BVCopyVec(X,1,Bx[i]);CHKERRQ(ierr);
     140             :       s = 1;
     141             :     }
     142        2515 :     for (j=0;j<s;j++) d->nX[r_s+i+j] = 1.0;
     143             : 
     144             :     /* Ax = R <- P*(Ax - eig_i*Bx) */
     145        2515 :     ierr = d->calcpairs_proj_res(d,r_s+i,r_s+i+s,&Ax[i]);CHKERRQ(ierr);
     146             : 
     147             :     /* Check if the first eigenpairs are converged */
     148        2515 :     if (i == 0) {
     149        2515 :       ierr = d->preTestConv(d,0,r_s+s,r_s+s,&d->npreconv);CHKERRQ(ierr);
     150        2515 :       if (d->npreconv > oldnpreconv) break;
     151             :     }
     152             :   }
     153             : 
     154             :   /* D <- K*[Ax Bx] */
     155        2515 :   if (d->npreconv <= oldnpreconv) {
     156        4920 :     for (i=0;i<n;i++) {
     157        2460 :       ierr = BVGetColumn(d->eps->V,kv+i,&v);CHKERRQ(ierr);
     158        2460 :       ierr = d->improvex_precond(d,r_s+i,Ax[i],v);CHKERRQ(ierr);
     159        2460 :       ierr = BVRestoreColumn(d->eps->V,kv+i,&v);CHKERRQ(ierr);
     160             :     }
     161        4920 :     for (i=n;i<n*2;i++) {
     162        2460 :       ierr = BVGetColumn(d->eps->V,kv+i,&v);CHKERRQ(ierr);
     163        2460 :       ierr = d->improvex_precond(d,r_s+i-n,Bx[i-n],v);CHKERRQ(ierr);
     164        2460 :       ierr = BVRestoreColumn(d->eps->V,kv+i,&v);CHKERRQ(ierr);
     165             :     }
     166        2460 :     *size_D = 2*n;
     167             : #if !defined(PETSC_USE_COMPLEX)
     168        2460 :     if (d->eigi[r_s] != 0.0) {
     169             :       s = 4;
     170             :     } else
     171             : #endif
     172             :     {
     173        2457 :       s = 2;
     174             :     }
     175             :     /* Prevent that short vectors are discarded in the orthogonalization */
     176        7380 :     for (i=0; i<s && i<*size_D; i++) {
     177        4920 :       if (d->eps->errest[d->nconv+r_s+i] > PETSC_MACHINE_EPSILON && d->eps->errest[d->nconv+r_s+i] < PETSC_MAX_REAL) {
     178        4920 :         ierr = BVScaleColumn(d->eps->V,i+kv,1.0/d->eps->errest[d->nconv+r_s+i]);CHKERRQ(ierr);
     179             :       }
     180             :     }
     181          55 :   } else *size_D = 0;
     182             : 
     183        2515 :   ierr = SlepcVecPoolRestoreVecs(d->auxV,n,&Bx);CHKERRQ(ierr);
     184        2515 :   ierr = SlepcVecPoolRestoreVecs(d->auxV,n,&Ax);CHKERRQ(ierr);
     185        2515 :   ierr = BVDestroy(&X);CHKERRQ(ierr);
     186        2515 :   ierr = MatDestroy(&M);CHKERRQ(ierr);
     187        2515 :   PetscFunctionReturn(0);
     188             : }
     189             : 
     190          51 : PetscErrorCode dvd_improvex_gd2(dvdDashboard *d,dvdBlackboard *b,KSP ksp,PetscInt max_bs)
     191             : {
     192             :   PetscErrorCode  ierr;
     193             :   dvdImprovex_gd2 *data;
     194             :   PC              pc;
     195             : 
     196          51 :   PetscFunctionBegin;
     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          51 :   if (!DVD_IS(d->sEP, DVD_EP_HERMITIAN)) {
     202          12 :     max_bs++;
     203          12 :     b->max_size_P = PetscMax(b->max_size_P,2);
     204             :   } else
     205             : #endif
     206             :   {
     207          39 :     b->max_size_P = PetscMax(b->max_size_P,1);
     208             :   }
     209          51 :   b->max_size_X = PetscMax(b->max_size_X,max_bs);
     210             : 
     211             :   /* Setup the preconditioner */
     212          51 :   if (ksp) {
     213          51 :     ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
     214          51 :     ierr = dvd_static_precond_PC(d,b,pc);CHKERRQ(ierr);
     215             :   } else {
     216           0 :     ierr = dvd_static_precond_PC(d,b,0);CHKERRQ(ierr);
     217             :   }
     218             : 
     219             :   /* Setup the step */
     220          51 :   if (b->state >= DVD_STATE_CONF) {
     221          17 :     ierr = PetscNewLog(d->eps,&data);CHKERRQ(ierr);
     222          17 :     d->improveX_data = data;
     223          17 :     data->size_X = b->max_size_X;
     224          17 :     d->improveX = dvd_improvex_gd2_gen;
     225             : 
     226          17 :     ierr = EPSDavidsonFLAdd(&d->destroyList,dvd_improvex_gd2_d);CHKERRQ(ierr);
     227             :   }
     228          51 :   PetscFunctionReturn(0);
     229             : }
     230             : 

Generated by: LCOV version 1.13