Subversion Repositories slepc-dev

/*                      

   SLEPc eigensolver: "krylovschur"

   Method: Krylov-Schur

   Algorithm:

       Single-vector Krylov-Schur method for both symmetric and non-symmetric

       problems.

   References:

       [1] "Krylov-Schur Methods in SLEPc", SLEPc Technical Report STR-7,

           available at http://www.grycap.upv.es/slepc.

       [2] G.W. Stewart, "A Krylov-Schur Algorithm for Large Eigenproblems",

           SIAM J. Matrix Analysis and App., 23(3), pp. 601-614, 2001.

   Last update: Oct 2006

   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

      SLEPc - Scalable Library for Eigenvalue Problem Computations

      Copyright (c) 2002-2007, Universidad Politecnica de Valencia, Spain

      This file is part of SLEPc. See the README file for conditions of use

      and additional information.

   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

*/

#include "src/eps/epsimpl.h"                /*I "slepceps.h" I*/

#include "slepcblaslapack.h"

#undef __FUNCT__  

#define __FUNCT__ "EPSSetUp_KRYLOVSCHUR"

PetscErrorCode EPSSetUp_KRYLOVSCHUR(EPS eps)

{

  PetscErrorCode ierr;

  PetscInt       N;

  PetscFunctionBegin;

  ierr = VecGetSize(eps->vec_initial,&N);CHKERRQ(ierr);

  if (eps->ncv) {

    if (eps->ncv<eps->nev+1) SETERRQ(1,"The value of ncv must be at least nev+1");

  }

  else eps->ncv = PetscMin(N,PetscMax(2*eps->nev,eps->nev+15));

  if (!eps->max_it) eps->max_it = PetscMax(100,2*N/eps->ncv);

  if (eps->ishermitian && (eps->which==EPS_LARGEST_IMAGINARY || eps->which==EPS_SMALLEST_IMAGINARY))

    SETERRQ(1,"Wrong value of eps->which");

  if (!eps->projection) {

    ierr = EPSSetProjection(eps,EPS_RITZ);CHKERRQ(ierr);

  } else if (eps->projection!=EPS_RITZ && eps->projection!=EPS_HARMONIC) {

    SETERRQ(PETSC_ERR_SUP,"Unsupported projection type\n");

  }

  ierr = EPSAllocateSolution(eps);CHKERRQ(ierr);

  ierr = PetscFree(eps->T);CHKERRQ(ierr);

  ierr = PetscMalloc(eps->ncv*eps->ncv*sizeof(PetscScalar),&eps->T);CHKERRQ(ierr);

  ierr = EPSDefaultGetWork(eps,2);CHKERRQ(ierr);

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "EPSTranslateHarmonic"

/*

   EPSTranslateHarmonic - Computes a translation of the Krylov decomposition

   in order to perform a harmonic projection.

   On input:

     A Krylov decomposition

                    OP * U = U * S + u * b^T

     S is the projected matrix (leading dimension is m)

     [U, u] is the basis of the Krylov subspace

     b is assumed to be beta*e_m^T

   Workspace:

     B is workspace to store a working copy of S and a working vector (last column)

     ipiv is workspace for pivots (int of length m)

   On output:

     S is updated as S + g*b', with g = (B-sigma*eye(m))'\b

     u is updated as u - U*g

     u is renormalized so beta is updated accordingly

*/

PetscErrorCode EPSTranslateHarmonic(EPS eps,PetscScalar *S,int m,PetscScalar *B,Vec *U,Vec u,PetscReal *beta,int *ipiv)

{

#if defined(PETSC_MISSING_LAPACK_GETRF) || defined(PETSC_MISSING_LAPACK_GETRS) 

  PetscFunctionBegin;

  SETERRQ(PETSC_ERR_SUP,"GETRF,GETRS - Lapack routines are unavailable.");

#else

  PetscErrorCode ierr;

  PetscBLASInt   info,one = 1;

  PetscReal      gamma;

  int            i;

  PetscFunctionBegin;

  /* Copy S to workspace B */

  ierr = PetscMemcpy(B,S,m*m*sizeof(PetscScalar));CHKERRQ(ierr);

  /* Last columns of B stores vectors b and g */

  ierr = PetscMemzero(B+m*m,m*sizeof(PetscScalar));CHKERRQ(ierr);

  B[(m-1)+m*m] = (PetscScalar)(*beta);

  /* g = (B-sigma*eye(m))'\b */

  for (i=0;i<m;i++)

    B[i+i*m] -= eps->target;

  LAPACKgetrf_(&m,&m,B,&m,ipiv,&info);

  if (info<0) SETERRQ(PETSC_ERR_LIB,"Bad argument to LU factorization");

  if (info>0) SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization");

  ierr = PetscLogFlops(2*m*m*m/3);CHKERRQ(ierr);

  LAPACKgetrs_("C",&m,&one,B,&m,ipiv,B+m*m,&m,&info);

  if (info) SETERRQ(PETSC_ERR_LIB,"GETRS - Bad solve");

  ierr = PetscLogFlops(2*m*m-m);CHKERRQ(ierr);

  /* S = S + g*b' */

  for (i=0;i<m;i++)

    S[i+(m-1)*m] = S[i+(m-1)*m] + B[i+m*m]*(*beta);

  /* u = u - U*g */

  for (i=0;i<m;i++)

    B[i+m*m] = -B[i+m*m];

  ierr = VecMAXPY(u,m,B+m*m,U);CHKERRQ(ierr);        

  /* Renormalize u */

  ierr = IPNorm(eps->ip,u,&gamma);CHKERRQ(ierr);

  ierr = VecScale(u,1.0/gamma);CHKERRQ(ierr);

  *beta = (*beta)*gamma;

  PetscFunctionReturn(0);

#endif

}

#undef __FUNCT__  

#define __FUNCT__ "EPSProjectedKSSymm"

/*

   EPSProjectedKSSym - Solves the projected eigenproblem in the Krylov-Schur

   method (symmetric case).

   On input:

     l is the number of vectors kept in previous restart (0 means first restart)

     S is the projected matrix (leading dimension is lds)

     Q is an orthogonal transformation matrix if l=0 (leading dimension is n)

   Workspace:

     ritz temporarily stores computed Ritz values

     perm is used for representing the permutation used for sorting values

   On output:

     S has (real) Schur form with diagonal blocks sorted appropriately

     Q contains the accumulated orthogonal transformations used in the process

*/

PetscErrorCode EPSProjectedKSSym(EPS eps,int l,PetscScalar *S,int lds,PetscScalar *Q,int n,PetscReal *ritz,int *perm)

{

  PetscErrorCode ierr;

  int            i,j;

  PetscFunctionBegin;

  /* Reduce S to diagonal form, S <- Q S Q' */

  if (l==0) {

    ierr = EPSDenseTridiagonal(n,S+eps->nconv*(lds+1),lds,ritz,Q+eps->nconv*n);CHKERRQ(ierr);

  } else {

    ierr = EPSDenseHEP(n,S+eps->nconv*(lds+1),lds,ritz,Q+eps->nconv*n);CHKERRQ(ierr);

  }

  /* Sort the remaining columns of the Schur form */

  if (eps->which == EPS_SMALLEST_REAL) {

    for (i=0;i<n;i++)

      eps->eigr[i+eps->nconv] = ritz[i];

  } else {

#ifdef PETSC_USE_COMPLEX

    for (i=0;i<n;i++)

      eps->eigr[i+eps->nconv] = ritz[i];

    ierr = EPSSortEigenvalues(n,eps->eigr+eps->nconv,eps->eigi,eps->which,n,perm);CHKERRQ(ierr);

#else

    ierr = EPSSortEigenvalues(n,ritz,eps->eigi+eps->nconv,eps->which,n,perm);CHKERRQ(ierr);

#endif

    for (i=0;i<n;i++)

      eps->eigr[i+eps->nconv] = ritz[perm[i]];

    ierr = PetscMemcpy(S,Q+eps->nconv*n,n*n*sizeof(PetscScalar));CHKERRQ(ierr);

    for (j=0;j<n;j++)

      for (i=0;i<n;i++)

        Q[(j+eps->nconv)*n+i] = S[perm[j]*n+i];

  }

  /* rebuild S from eigr */

  for (i=eps->nconv;i<eps->nv;i++) {

    S[i*(eps->ncv+1)] = eps->eigr[i];

    for (j=i+1;j<eps->ncv;j++)

      S[i*eps->ncv+j] = 0.0;

  }

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "EPSProjectedKSNonsymm"

/*

   EPSProjectedKSNonsym - Solves the projected eigenproblem in the Krylov-Schur

   method (non-symmetric case).

   On input:

     l is the number of vectors kept in previous restart (0 means first restart)

     S is the projected matrix (leading dimension is lds)

     Q is an orthogonal transformation matrix if l=0 (leading dimension is n)

   On output:

     S has (real) Schur form with diagonal blocks sorted appropriately

     Q contains the accumulated orthogonal transformations used in the process

*/

PetscErrorCode EPSProjectedKSNonsym(EPS eps,int l,PetscScalar *S,int lds,PetscScalar *Q,int n)

{

  PetscErrorCode ierr;

  int            i;

  PetscFunctionBegin;

  if (l==0) {

    ierr = PetscMemzero(Q,n*n*sizeof(PetscScalar));CHKERRQ(ierr);

    for (i=0;i<n;i++)

      Q[i*(n+1)] = 1.0;

  } else {

    /* Reduce S to Hessenberg form, S <- Q S Q' */

    ierr = EPSDenseHessenberg(n,eps->nconv,S,lds,Q);CHKERRQ(ierr);

  }

  /* Reduce S to (quasi-)triangular form, S <- Q S Q' */

  ierr = EPSDenseSchur(n,eps->nconv,S,lds,Q,eps->eigr,eps->eigi);CHKERRQ(ierr);

  /* Sort the remaining columns of the Schur form */

  if (eps->projection==EPS_HARMONIC) {

    ierr = EPSSortDenseSchurTarget(n,eps->nconv,S,lds,Q,eps->eigr,eps->eigi,eps->target);CHKERRQ(ierr);    

  } else {

    ierr = EPSSortDenseSchur(n,eps->nconv,S,lds,Q,eps->eigr,eps->eigi,eps->which);CHKERRQ(ierr);    

  }

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "EPSSolve_KRYLOVSCHUR"

PetscErrorCode EPSSolve_KRYLOVSCHUR(EPS eps)

{

  PetscErrorCode ierr;

  int            i,j,k,l,n,lwork,*perm;

  Vec            u=eps->work[1];

  PetscScalar    *S=eps->T,*B,*Q,*work;

  PetscReal      beta,*ritz;

  PetscTruth     breakdown;

  PetscFunctionBegin;

  ierr = PetscMemzero(S,eps->ncv*eps->ncv*sizeof(PetscScalar));CHKERRQ(ierr);

  ierr = PetscMalloc(eps->ncv*(eps->ncv+1)*sizeof(PetscScalar),&B);CHKERRQ(ierr);

  ierr = PetscMalloc(eps->ncv*eps->ncv*sizeof(PetscScalar),&Q);CHKERRQ(ierr);

  lwork = (eps->ncv+4)*eps->ncv;

  if (!eps->ishermitian) {

    ierr = PetscMalloc(lwork*sizeof(PetscScalar),&work);CHKERRQ(ierr);

  } else {

    ierr = PetscMalloc(eps->ncv*sizeof(PetscReal),&ritz);CHKERRQ(ierr);

  }

  ierr = PetscMalloc(eps->ncv*sizeof(int),&perm);CHKERRQ(ierr);

  /* Get the starting Arnoldi vector */

  ierr = EPSGetStartVector(eps,0,eps->V[0],PETSC_NULL);CHKERRQ(ierr);

  l = 0;

  /* Restart loop */

  while (eps->reason == EPS_CONVERGED_ITERATING) {

    eps->its++;

    /* Compute an nv-step Arnoldi factorization */

    eps->nv = eps->ncv;

    ierr = EPSBasicArnoldi(eps,PETSC_FALSE,S,eps->V,eps->nconv+l,&eps->nv,u,&beta,&breakdown);CHKERRQ(ierr);

    ierr = VecScale(u,1.0/beta);CHKERRQ(ierr);

    /* Compute translation of Krylov decomposition if harmonic projection used */

    if (eps->projection==EPS_HARMONIC) {

      ierr = EPSTranslateHarmonic(eps,S,eps->ncv,B,eps->V,u,&beta,perm);CHKERRQ(ierr);

    }

    /* Solve projected problem and compute residual norm estimates */

    if (eps->ishermitian) {

      n = eps->nv-eps->nconv;

      ierr = EPSProjectedKSSym(eps,l,S,eps->ncv,Q,n,ritz,perm);CHKERRQ(ierr);

      for (i=eps->nconv;i<eps->nv;i++)

        eps->errest[i] = beta*PetscAbsScalar(Q[(i+1)*n-1]) / PetscAbsScalar(eps->eigr[i]);

    } else { /* non-hermitian */

      n = eps->nv;

      ierr = EPSProjectedKSNonsym(eps,l,S,eps->ncv,Q,n);CHKERRQ(ierr);

      ierr = ArnoldiResiduals(S,eps->ncv,Q,beta,eps->nconv,n,eps->eigr,eps->eigi,eps->errest,work);CHKERRQ(ierr);

    }

    /* Check convergence */

    k = eps->nconv;

    while (k<eps->nv && eps->errest[k]<eps->tol) k++;    

    if (eps->its >= eps->max_it) eps->reason = EPS_DIVERGED_ITS;

    if (k >= eps->nev) eps->reason = EPS_CONVERGED_TOL;

    /* Update l */

    if (eps->reason != EPS_CONVERGED_ITERATING || breakdown) l = 0;

    else {

      l = (eps->nv-k)/2;

#if !defined(PETSC_USE_COMPLEX)

      if (S[(k+l-1)*(eps->ncv+1)+1] != 0.0) {

        if (k+l<eps->nv-1) l = l+1;

        else l = l-1;

      }

#endif

    }

    /* Update the corresponding vectors V(:,idx) = V*Q(:,idx) */

    for (i=eps->nconv;i<k+l;i++) {

      ierr = VecSet(eps->AV[i],0.0);CHKERRQ(ierr);

      if (!eps->ishermitian) {

        ierr = VecMAXPY(eps->AV[i],n,Q+i*n,eps->V);CHKERRQ(ierr);        

      } else {

        ierr = VecMAXPY(eps->AV[i],n,Q+i*n,eps->V+eps->nconv);CHKERRQ(ierr);

      }

    }

    for (i=eps->nconv;i<k+l;i++) {

      ierr = VecCopy(eps->AV[i],eps->V[i]);CHKERRQ(ierr);

    }

    eps->nconv = k;

    EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,eps->nv);

    if (eps->reason == EPS_CONVERGED_ITERATING) {

      if (breakdown) {

        /* start a new Arnoldi factorization */

        PetscInfo2(eps,"Breakdown in Krylov-Schur method (it=%i norm=%g)\n",eps->its,beta);

        ierr = EPSGetStartVector(eps,k,eps->V[k],&breakdown);CHKERRQ(ierr);

        if (breakdown) {

          eps->reason = EPS_DIVERGED_BREAKDOWN;

          PetscInfo(eps,"Unable to generate more start vectors\n");

        }

      } else {

        /* update the Arnoldi-Schur decomposition */

        for (i=k;i<k+l;i++) {

          S[i*eps->ncv+k+l] = Q[(i+1)*n-1]*beta;

        }

        if (eps->projection==EPS_HARMONIC) {

          /* force orthogonality of u */

          ierr = IPOrthogonalize(eps->ip,l,PETSC_NULL,eps->V,u,B+eps->ncv*eps->ncv,&beta,PETSC_NULL,eps->work[0]);CHKERRQ(ierr);

          /* S = S + g*b' */

          for (i=0;i<l;i++) {

            for (j=k;j<k+l;j++) {

              S[i+j*eps->ncv] += B[i+eps->ncv*eps->ncv]*S[k+l+j*eps->ncv];

            }

          }

          ierr = VecScale(u,1.0/beta);CHKERRQ(ierr);

          for (i=k;i<k+l;i++) {

            S[i*eps->ncv+k+l] *= beta;

          }

        }

        ierr = VecCopy(u,eps->V[k+l]);CHKERRQ(ierr);

      }

    }

  }

  ierr = PetscFree(B);CHKERRQ(ierr);

  ierr = PetscFree(Q);CHKERRQ(ierr);

  if (!eps->ishermitian) {

    ierr = PetscFree(work);CHKERRQ(ierr);

  } else {

    ierr = PetscFree(ritz);CHKERRQ(ierr);

  }

  ierr = PetscFree(perm);CHKERRQ(ierr);

  PetscFunctionReturn(0);

}

EXTERN_C_BEGIN

#undef __FUNCT__  

#define __FUNCT__ "EPSCreate_KRYLOVSCHUR"

PetscErrorCode EPSCreate_KRYLOVSCHUR(EPS eps)

{

  PetscFunctionBegin;

  eps->data                      = PETSC_NULL;

  eps->ops->solve                = EPSSolve_KRYLOVSCHUR;

  eps->ops->solvets              = PETSC_NULL;

  eps->ops->setup                = EPSSetUp_KRYLOVSCHUR;

  eps->ops->setfromoptions       = PETSC_NULL;

  eps->ops->destroy              = EPSDestroy_Default;

  eps->ops->view                 = PETSC_NULL;

  eps->ops->backtransform        = EPSBackTransform_Default;

  eps->ops->computevectors       = EPSComputeVectors_Schur;

  PetscFunctionReturn(0);

}

EXTERN_C_END