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: Krylov-Schur for symmetric-indefinite eigenproblems
      14             : */
      15             : #include <slepc/private/epsimpl.h>
      16             : #include "krylovschur.h"
      17             : 
      18          13 : PetscErrorCode EPSSolve_KrylovSchur_Indefinite(EPS eps)
      19             : {
      20          13 :   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
      21          13 :   PetscInt        k,l,ld,nv,t,nconv=0;
      22          13 :   Mat             U,D;
      23          13 :   Vec             vomega,w=eps->work[0];
      24          13 :   PetscReal       *a,*b,beta,beta1=1.0,*omega;
      25          13 :   PetscBool       breakdown=PETSC_FALSE,symmlost=PETSC_FALSE;
      26             : 
      27          13 :   PetscFunctionBegin;
      28          13 :   PetscCall(DSGetLeadingDimension(eps->ds,&ld));
      29             : 
      30             :   /* Get the starting Lanczos vector */
      31          13 :   PetscCall(EPSGetStartVector(eps,0,NULL));
      32             : 
      33             :   /* Extract sigma[0] from BV, computed during normalization */
      34          13 :   PetscCall(DSSetDimensions(eps->ds,1,PETSC_DEFAULT,PETSC_DEFAULT));
      35          13 :   PetscCall(BVSetActiveColumns(eps->V,0,1));
      36          13 :   PetscCall(DSGetMatAndColumn(eps->ds,DS_MAT_D,0,&D,&vomega));
      37          13 :   PetscCall(BVGetSignature(eps->V,vomega));
      38          13 :   PetscCall(DSRestoreMatAndColumn(eps->ds,DS_MAT_D,0,&D,&vomega));
      39          13 :   l = 0;
      40             : 
      41             :   /* Restart loop */
      42          13 :   while (eps->reason == EPS_CONVERGED_ITERATING) {
      43          61 :     eps->its++;
      44             : 
      45             :     /* Compute an nv-step Lanczos factorization */
      46          61 :     nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
      47          61 :     PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
      48          61 :     PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&a));
      49          61 :     b = a + ld;
      50          61 :     PetscCall(DSGetArrayReal(eps->ds,DS_MAT_D,&omega));
      51          61 :     PetscCall(EPSPseudoLanczos(eps,a,b,omega,eps->nconv+l,&nv,&breakdown,&symmlost,NULL,w));
      52          61 :     if (symmlost) {
      53           0 :       eps->reason = EPS_DIVERGED_SYMMETRY_LOST;
      54           0 :       if (nv==eps->nconv+l+1) { eps->nconv = nconv; break; }
      55             :     }
      56          61 :     beta = b[nv-1];
      57          61 :     PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&a));
      58          61 :     PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_D,&omega));
      59          61 :     PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
      60          61 :     PetscCall(DSSetState(eps->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
      61          61 :     PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
      62             : 
      63             :     /* Solve projected problem */
      64          61 :     PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
      65          61 :     PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
      66          61 :     PetscCall(DSUpdateExtraRow(eps->ds));
      67          61 :     PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
      68             : 
      69             :     /* Check convergence */
      70          61 :     PetscCall(DSGetDimensions(eps->ds,NULL,NULL,NULL,&t));
      71             : #if 0
      72             :     /* take into account also left residual */
      73             :     PetscCall(BVGetColumn(eps->V,nv,&u));
      74             :     PetscCall(VecNorm(u,NORM_2,&beta1));
      75             :     PetscCall(BVRestoreColumn(eps->V,nv,&u));
      76             :     PetscCall(VecNorm(w,NORM_2,&beta2));  /* w contains B*V[nv] */
      77             :     beta1 = PetscMax(beta1,beta2);
      78             : #endif
      79          61 :     PetscCall(EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,t-eps->nconv,beta*beta1,0.0,1.0,&k));
      80          61 :     if (!symmlost) PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx));
      81          61 :     nconv = k;
      82             : 
      83             :     /* Update l */
      84          61 :     if (eps->reason != EPS_CONVERGED_ITERATING || breakdown) l = 0;
      85             :     else {
      86          48 :       l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
      87          48 :       l = PetscMin(l,t);
      88          48 :       PetscCall(DSGetTruncateSize(eps->ds,k,t,&l));
      89             :     }
      90          61 :     if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
      91          61 :     if (l) PetscCall(PetscInfo(eps,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
      92             : 
      93          61 :     if (eps->reason == EPS_CONVERGED_ITERATING) {
      94          48 :       PetscCheck(!breakdown,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Breakdown in Indefinite Krylov-Schur (beta=%g)",(double)beta);
      95             :       /* Prepare the Rayleigh quotient for restart */
      96          48 :       PetscCall(DSTruncate(eps->ds,k+l,PETSC_FALSE));
      97             :     }
      98             :     /* Update the corresponding vectors V(:,idx) = V*Q(:,idx) */
      99          61 :     PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&U));
     100          61 :     PetscCall(BVMultInPlace(eps->V,U,eps->nconv,k+l));
     101          61 :     PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&U));
     102             : 
     103             :     /* Move restart vector and update signature */
     104          61 :     if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) {
     105          48 :       PetscCall(BVCopyColumn(eps->V,nv,k+l));
     106          48 :       PetscCall(BVSetActiveColumns(eps->V,0,k+l));
     107          48 :       PetscCall(DSGetMatAndColumn(eps->ds,DS_MAT_D,0,&D,&vomega));
     108          48 :       PetscCall(BVSetSignature(eps->V,vomega));
     109          48 :       PetscCall(DSRestoreMatAndColumn(eps->ds,DS_MAT_D,0,&D,&vomega));
     110             :     }
     111             : 
     112          61 :     eps->nconv = k;
     113          74 :     PetscCall(EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv));
     114             :   }
     115          13 :   PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
     116          13 :   PetscFunctionReturn(PETSC_SUCCESS);
     117             : }