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: "krylovschur"
      12             : 
      13             :    Method: Two-sided Arnoldi with Krylov-Schur restart (for left eigenvectors)
      14             : 
      15             :    References:
      16             : 
      17             :        [1] I.N. Zwaan and M.E. Hochstenbach, "Krylov-Schur-type restarts
      18             :            for the two-sided Arnoldi method", SIAM J. Matrix Anal. Appl.
      19             :            38(2):297-321, 2017.
      20             : 
      21             : */
      22             : 
      23             : #include <slepc/private/epsimpl.h>
      24             : #include "krylovschur.h"
      25             : #include <slepcblaslapack.h>
      26             : 
      27          82 : static PetscErrorCode EPSTwoSidedRQUpdate1(EPS eps,Mat M,PetscInt nv,PetscReal beta,PetscReal betat)
      28             : {
      29          82 :   PetscScalar       *T,*S,*A,*w;
      30          82 :   const PetscScalar *pM;
      31          82 :   Vec               u;
      32          82 :   PetscInt          ld,ncv=eps->ncv,i,l,nnv;
      33          82 :   PetscBLASInt      info,n_,ncv_,*p,one=1;
      34             : 
      35          82 :   PetscFunctionBegin;
      36          82 :   PetscCall(DSGetLeadingDimension(eps->ds,&ld));
      37          82 :   PetscCall(PetscMalloc3(nv,&p,ncv*ncv,&A,ncv,&w));
      38          82 :   PetscCall(BVGetActiveColumns(eps->V,&l,&nnv));
      39          82 :   PetscCall(BVSetActiveColumns(eps->V,0,nv));
      40          82 :   PetscCall(BVSetActiveColumns(eps->W,0,nv));
      41          82 :   PetscCall(BVGetColumn(eps->V,nv,&u));
      42          82 :   PetscCall(BVDotVec(eps->W,u,w));
      43          82 :   PetscCall(BVRestoreColumn(eps->V,nv,&u));
      44          82 :   PetscCall(MatDenseGetArrayRead(M,&pM));
      45          82 :   PetscCall(PetscArraycpy(A,pM,ncv*ncv));
      46          82 :   PetscCall(MatDenseRestoreArrayRead(M,&pM));
      47          82 :   PetscCall(PetscBLASIntCast(nv,&n_));
      48          82 :   PetscCall(PetscBLASIntCast(ncv,&ncv_));
      49          82 :   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
      50          82 :   PetscCallBLAS("LAPACKgetrf",LAPACKgetrf_(&n_,&n_,A,&ncv_,p,&info));
      51          82 :   SlepcCheckLapackInfo("getrf",info);
      52          82 :   PetscCall(PetscLogFlops(2.0*n_*n_*n_/3.0));
      53          82 :   PetscCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&n_,&one,A,&ncv_,p,w,&ncv_,&info));
      54          82 :   SlepcCheckLapackInfo("getrs",info);
      55          82 :   PetscCall(PetscLogFlops(2.0*n_*n_-n_));
      56          82 :   PetscCall(BVMultColumn(eps->V,-1.0,1.0,nv,w));
      57          82 :   PetscCall(DSGetArray(eps->ds,DS_MAT_A,&S));
      58        1650 :   for (i=0;i<nv;i++) S[(nv-1)*ld+i] += beta*w[i];
      59          82 :   PetscCall(DSRestoreArray(eps->ds,DS_MAT_A,&S));
      60          82 :   PetscCall(BVGetColumn(eps->W,nv,&u));
      61          82 :   PetscCall(BVDotVec(eps->V,u,w));
      62          82 :   PetscCall(BVRestoreColumn(eps->W,nv,&u));
      63          82 :   PetscCallBLAS("LAPACKgetrs",LAPACKgetrs_("C",&n_,&one,A,&ncv_,p,w,&ncv_,&info));
      64          82 :   PetscCall(PetscFPTrapPop());
      65          82 :   PetscCall(BVMultColumn(eps->W,-1.0,1.0,nv,w));
      66          82 :   PetscCall(DSGetArray(eps->ds,DS_MAT_B,&T));
      67        1650 :   for (i=0;i<nv;i++) T[(nv-1)*ld+i] += betat*w[i];
      68          82 :   PetscCall(DSRestoreArray(eps->ds,DS_MAT_B,&T));
      69          82 :   PetscCall(PetscFree3(p,A,w));
      70          82 :   PetscCall(BVSetActiveColumns(eps->V,l,nnv));
      71          82 :   PetscCall(BVSetActiveColumns(eps->W,l,nnv));
      72          82 :   PetscFunctionReturn(PETSC_SUCCESS);
      73             : }
      74             : 
      75          69 : static PetscErrorCode EPSTwoSidedRQUpdate2(EPS eps,Mat M,PetscInt k)
      76             : {
      77          69 :   PetscScalar    *Q,*pM,*w,zero=0.0,sone=1.0,*c,*A;
      78          69 :   PetscBLASInt   n_,ncv_,ld_;
      79          69 :   PetscReal      norm;
      80          69 :   PetscInt       l,nv,ncv=eps->ncv,ld,i,j;
      81             : 
      82          69 :   PetscFunctionBegin;
      83          69 :   PetscCall(DSGetLeadingDimension(eps->ds,&ld));
      84          69 :   PetscCall(BVGetActiveColumns(eps->V,&l,&nv));
      85          69 :   PetscCall(BVSetActiveColumns(eps->V,0,nv));
      86          69 :   PetscCall(BVSetActiveColumns(eps->W,0,nv));
      87          69 :   PetscCall(PetscMalloc2(ncv*ncv,&w,ncv,&c));
      88             :   /* u = u - V*V'*u */
      89          69 :   PetscCall(BVOrthogonalizeColumn(eps->V,k,c,&norm,NULL));
      90          69 :   PetscCall(BVScaleColumn(eps->V,k,1.0/norm));
      91          69 :   PetscCall(DSGetArray(eps->ds,DS_MAT_A,&A));
      92             :   /* H = H + V'*u*b' */
      93         714 :   for (j=l;j<k;j++) {
      94        6885 :     for (i=0;i<k;i++) A[i+j*ld] += c[i]*A[k+j*ld];
      95         645 :     A[k+j*ld] *= norm;
      96             :   }
      97          69 :   PetscCall(DSRestoreArray(eps->ds,DS_MAT_A,&A));
      98          69 :   PetscCall(BVOrthogonalizeColumn(eps->W,k,c,&norm,NULL));
      99          69 :   PetscCall(BVScaleColumn(eps->W,k,1.0/norm));
     100          69 :   PetscCall(DSGetArray(eps->ds,DS_MAT_B,&A));
     101             :   /* H = H + V'*u*b' */
     102         714 :   for (j=l;j<k;j++) {
     103        6885 :     for (i=0;i<k;i++) A[i+j*ld] += c[i]*A[k+j*ld];
     104         645 :     A[k+j*ld] *= norm;
     105             :   }
     106          69 :   PetscCall(DSRestoreArray(eps->ds,DS_MAT_B,&A));
     107             : 
     108             :   /* M = Q'*M*Q */
     109          69 :   PetscCall(MatDenseGetArray(M,&pM));
     110          69 :   PetscCall(PetscBLASIntCast(ncv,&ncv_));
     111          69 :   PetscCall(PetscBLASIntCast(nv,&n_));
     112          69 :   PetscCall(PetscBLASIntCast(ld,&ld_));
     113          69 :   PetscCall(DSGetArray(eps->ds,DS_MAT_Q,&Q));
     114          69 :   PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&n_,&n_,&sone,pM,&ncv_,Q,&ld_,&zero,w,&ncv_));
     115          69 :   PetscCall(DSRestoreArray(eps->ds,DS_MAT_Q,&Q));
     116          69 :   PetscCall(DSGetArray(eps->ds,DS_MAT_Z,&Q));
     117          69 :   PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&n_,&n_,&n_,&sone,Q,&ld_,w,&ncv_,&zero,pM,&ncv_));
     118          69 :   PetscCall(DSRestoreArray(eps->ds,DS_MAT_Z,&Q));
     119          69 :   PetscCall(MatDenseRestoreArray(M,&pM));
     120          69 :   PetscCall(PetscFree2(w,c));
     121          69 :   PetscCall(BVSetActiveColumns(eps->V,l,nv));
     122          69 :   PetscCall(BVSetActiveColumns(eps->W,l,nv));
     123          69 :   PetscFunctionReturn(PETSC_SUCCESS);
     124             : }
     125             : 
     126          13 : PetscErrorCode EPSSolve_KrylovSchur_TwoSided(EPS eps)
     127             : {
     128          13 :   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
     129          13 :   Mat             M,U,Op,OpHT,S,T;
     130          13 :   PetscReal       norm,norm2,beta,betat;
     131          13 :   PetscInt        ld,l,nv,nvt,k,nconv,dsn,dsk;
     132          13 :   PetscBool       breakdownt,breakdown,breakdownl;
     133             : 
     134          13 :   PetscFunctionBegin;
     135          13 :   PetscCall(DSGetLeadingDimension(eps->ds,&ld));
     136          13 :   PetscCall(EPSGetStartVector(eps,0,NULL));
     137          13 :   PetscCall(EPSGetLeftStartVector(eps,0,NULL));
     138          13 :   l = 0;
     139          13 :   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,eps->ncv,eps->ncv,NULL,&M));
     140             : 
     141          13 :   PetscCall(STGetOperator(eps->st,&Op));
     142          13 :   PetscCall(MatCreateHermitianTranspose(Op,&OpHT));
     143             : 
     144             :   /* Restart loop */
     145          95 :   while (eps->reason == EPS_CONVERGED_ITERATING) {
     146          82 :     eps->its++;
     147             : 
     148             :     /* Compute an nv-step Arnoldi factorization for Op */
     149          82 :     nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
     150          82 :     PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
     151          82 :     PetscCall(DSGetMat(eps->ds,DS_MAT_A,&S));
     152          82 :     PetscCall(BVMatArnoldi(eps->V,Op,S,eps->nconv+l,&nv,&beta,&breakdown));
     153          82 :     PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&S));
     154             : 
     155             :     /* Compute an nv-step Arnoldi factorization for Op' */
     156          82 :     nvt = nv;
     157          82 :     PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
     158          82 :     PetscCall(DSGetMat(eps->ds,DS_MAT_B,&T));
     159          82 :     PetscCall(BVMatArnoldi(eps->W,OpHT,T,eps->nconv+l,&nvt,&betat,&breakdownt));
     160          82 :     PetscCall(DSRestoreMat(eps->ds,DS_MAT_B,&T));
     161             : 
     162             :     /* Make sure both factorizations have the same length */
     163          82 :     nv = PetscMin(nv,nvt);
     164          82 :     PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
     165          82 :     if (l==0) PetscCall(DSSetState(eps->ds,DS_STATE_INTERMEDIATE));
     166          69 :     else PetscCall(DSSetState(eps->ds,DS_STATE_RAW));
     167          82 :     breakdown = (breakdown || breakdownt)? PETSC_TRUE: PETSC_FALSE;
     168             : 
     169             :     /* Update M, modify Rayleigh quotients S and T */
     170          82 :     PetscCall(BVSetActiveColumns(eps->V,eps->nconv+l,nv));
     171          82 :     PetscCall(BVSetActiveColumns(eps->W,eps->nconv+l,nv));
     172          82 :     PetscCall(BVMatProject(eps->V,NULL,eps->W,M));
     173             : 
     174          82 :     PetscCall(EPSTwoSidedRQUpdate1(eps,M,nv,beta,betat));
     175             : 
     176             :     /* Solve projected problem */
     177          82 :     PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
     178          82 :     PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
     179          82 :     PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
     180          82 :     PetscCall(DSUpdateExtraRow(eps->ds));
     181             : 
     182             :     /* Check convergence */
     183          82 :     PetscCall(BVNormColumn(eps->V,nv,NORM_2,&norm));
     184          82 :     PetscCall(BVNormColumn(eps->W,nv,NORM_2,&norm2));
     185          82 :     PetscCall(EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,beta*norm,betat*norm2,1.0,&k));
     186          82 :     PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx));
     187          82 :     nconv = k;
     188             : 
     189             :     /* Update l */
     190          82 :     if (eps->reason != EPS_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
     191             :     else {
     192          69 :       l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
     193          69 :       PetscCall(DSGetTruncateSize(eps->ds,k,nv,&l));
     194             :     }
     195          82 :     if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
     196          82 :     if (l) PetscCall(PetscInfo(eps,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
     197             : 
     198             :     /* Update the corresponding vectors V(:,idx) = V*Q(:,idx) */
     199          82 :     PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
     200          82 :     PetscCall(BVSetActiveColumns(eps->W,eps->nconv,nv));
     201          82 :     PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&U));
     202          82 :     PetscCall(BVMultInPlace(eps->V,U,eps->nconv,k+l));
     203          82 :     PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&U));
     204          82 :     PetscCall(DSGetMat(eps->ds,DS_MAT_Z,&U));
     205          82 :     PetscCall(BVMultInPlace(eps->W,U,eps->nconv,k+l));
     206          82 :     PetscCall(DSRestoreMat(eps->ds,DS_MAT_Z,&U));
     207          82 :     if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) {
     208          69 :       PetscCall(BVCopyColumn(eps->V,nv,k+l));
     209          69 :       PetscCall(BVCopyColumn(eps->W,nv,k+l));
     210             :     }
     211             : 
     212          82 :     if (eps->reason == EPS_CONVERGED_ITERATING) {
     213          69 :       if (breakdown || k==nv) {
     214             :         /* Start a new Arnoldi factorization */
     215           0 :         PetscCall(PetscInfo(eps,"Breakdown in Krylov-Schur method (it=%" PetscInt_FMT " norm=%g)\n",eps->its,(double)beta));
     216           0 :         if (k<eps->nev) {
     217           0 :           PetscCall(EPSGetStartVector(eps,k,&breakdown));
     218           0 :           PetscCall(EPSGetLeftStartVector(eps,k,&breakdownl));
     219           0 :           if (breakdown || breakdownl) {
     220           0 :             eps->reason = EPS_DIVERGED_BREAKDOWN;
     221           0 :             PetscCall(PetscInfo(eps,"Unable to generate more start vectors\n"));
     222             :           }
     223             :         }
     224             :       } else {
     225          69 :         PetscCall(DSGetDimensions(eps->ds,&dsn,NULL,&dsk,NULL));
     226          69 :         PetscCall(DSSetDimensions(eps->ds,dsn,k,dsk));
     227          69 :         PetscCall(DSTruncate(eps->ds,k+l,PETSC_FALSE));
     228             :       }
     229          69 :       PetscCall(EPSTwoSidedRQUpdate2(eps,M,k+l));
     230             :     }
     231          82 :     eps->nconv = k;
     232          95 :     PetscCall(EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv));
     233             :   }
     234             : 
     235          13 :   PetscCall(STRestoreOperator(eps->st,&Op));
     236          13 :   PetscCall(MatDestroy(&OpHT));
     237             : 
     238          13 :   PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
     239          13 :   PetscCall(MatDestroy(&M));
     240          13 :   PetscFunctionReturn(PETSC_SUCCESS);
     241             : }