Subversion Repositories slepc-dev

/*

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

   SLEPc - Scalable Library for Eigenvalue Problem Computations

   Copyright (c) 2002-2011, Universitat Politecnica de Valencia, Spain

   This file is part of SLEPc.

   SLEPc is free software: you can redistribute it and/or modify it under  the

   terms of version 3 of the GNU Lesser General Public License as published by

   the Free Software Foundation.

   SLEPc  is  distributed in the hope that it will be useful, but WITHOUT  ANY

   WARRANTY;  without even the implied warranty of MERCHANTABILITY or  FITNESS

   FOR  A  PARTICULAR PURPOSE. See the GNU Lesser General Public  License  for

   more details.

   You  should have received a copy of the GNU Lesser General  Public  License

   along with SLEPc. If not, see <http://www.gnu.org/licenses/>.

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

*/

#include <slepc-private/psimpl.h>      /*I "slepcps.h" I*/

#include <slepcblaslapack.h>

#undef __FUNCT__  

#define __FUNCT__ "PSAllocate_NHEP"

PetscErrorCode PSAllocate_NHEP(PS ps,PetscInt ld)

{

  PetscErrorCode ierr;

  PetscFunctionBegin;

  ierr = PSAllocateMat_Private(ps,PS_MAT_A);CHKERRQ(ierr);

  ierr = PSAllocateMat_Private(ps,PS_MAT_Q);CHKERRQ(ierr);

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "PSView_NHEP"

PetscErrorCode PSView_NHEP(PS ps,PetscViewer viewer)

{

  PetscErrorCode ierr;

  PetscFunctionBegin;

  ierr = PSViewMat_Private(ps,viewer,PS_MAT_A);CHKERRQ(ierr);

  if (ps->state>PS_STATE_INTERMEDIATE) {

    ierr = PSViewMat_Private(ps,viewer,PS_MAT_Q);CHKERRQ(ierr);

  }

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "PSVectors_NHEP_Eigen_Some"

PetscErrorCode PSVectors_NHEP_Eigen_Some(PS ps,PetscInt *k,PetscReal *rnorm,PetscBool left)

{

#if defined(SLEPC_MISSING_LAPACK_TREVC)

  PetscFunctionBegin;

  SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"TREVC - Lapack routine is unavailable.");

#else

  PetscErrorCode ierr;

  PetscInt       i;

  PetscBLASInt   mm,mout,info,ld,n,inc = 1;

  PetscScalar    tmp,done=1.0,zero=0.0;

  PetscReal      norm;

  PetscBool      iscomplex = PETSC_FALSE;

  PetscBLASInt   *select;

  PetscScalar    *A = ps->mat[PS_MAT_A];

  PetscScalar    *Q = ps->mat[PS_MAT_Q];

  PetscScalar    *X = ps->mat[PS_MAT_X];

  PetscScalar    *Y;

#if defined(PETSC_USE_COMPLEX)

  PetscReal      *rwork;

#endif

  PetscFunctionBegin;

  if (ps->state<PS_STATE_CONDENSED) SETERRQ(((PetscObject)ps)->comm,PETSC_ERR_ORDER,"Must call PSSolve() first");

  if (left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented for left eigenvectors");

  n  = PetscBLASIntCast(ps->n);

  ld = PetscBLASIntCast(ps->ld);

  if ((*k)<n-1 && A[(*k)+1+(*k)*ld]!=0.0) iscomplex = PETSC_TRUE;

  ierr = PSAllocateWork_Private(ps,0,0,ld);CHKERRQ(ierr);

  select = ps->iwork;

  for (i=0;i<n;i++) select[i] = (PetscBLASInt)PETSC_FALSE;

  /* Compute k-th eigenvector Y of A */

  Y = X+(*k)*ld;

  mm = iscomplex? 2: 1;

  select[*k] = (PetscBLASInt)PETSC_TRUE;

#if !defined(PETSC_USE_COMPLEX)

  if (iscomplex) select[(*k)+1] = (PetscBLASInt)PETSC_TRUE;

  ierr = PSAllocateWork_Private(ps,3*ld,0,0);CHKERRQ(ierr);

  LAPACKtrevc_("R","S",select,&n,A,&ld,PETSC_NULL,&ld,Y,&ld,&mm,&mout,ps->work,&info);

#else

  ierr = PSAllocateWork_Private(ps,2*ld,ld,0);CHKERRQ(ierr);

  LAPACKtrevc_("R","S",select,&n,A,&ld,PETSC_NULL,&ld,Y,&ld,&mm,&mout,ps->work,ps->rwork,&info);

#endif

  if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in Lapack xTREVC %d",info);

  if (mout != mm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Inconsistent arguments");

  ierr = PetscMemcpy(ps->work,Y,mout*ld*sizeof(PetscScalar));CHKERRQ(ierr);

  /* accumulate and normalize eigenvectors */

  BLASgemv_("N",&n,&n,&done,Q,&ld,ps->work,&inc,&zero,Y,&inc);

#if !defined(PETSC_USE_COMPLEX)

  if (iscomplex) BLASgemv_("N",&n,&n,&done,Q,&ld,ps->work+ld,&inc,&zero,Y+ld,&inc);

#endif

  norm = BLASnrm2_(&n,Y,&inc);

#if !defined(PETSC_USE_COMPLEX)

  tmp  = BLASnrm2_(&n,Y+ld,&inc);

  norm = SlepcAbsEigenvalue(norm,tmp);

#endif

  tmp = 1.0 / norm;

  BLASscal_(&n,&tmp,Y,&inc);

#if !defined(PETSC_USE_COMPLEX)

  BLASscal_(&n,&tmp,Y+ld,&inc);

#endif

  /* set output arguments */

  if (iscomplex) (*k)++;

  if (rnorm) {

    if (iscomplex) *rnorm = SlepcAbsEigenvalue(Y[n-1],Y[n-1+ld]);

    else *rnorm = PetscAbsScalar(Y[n-1]);

  }

  PetscFunctionReturn(0);

#endif

}

#undef __FUNCT__  

#define __FUNCT__ "PSVectors_NHEP_Eigen_All"

PetscErrorCode PSVectors_NHEP_Eigen_All(PS ps,PetscBool left)

{

#if defined(SLEPC_MISSING_LAPACK_TREVC)

  PetscFunctionBegin;

  SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"TREVC - Lapack routine is unavailable.");

#else

  PetscErrorCode ierr;

  PetscBLASInt   n,ld,mout,info;

  PetscScalar    *A = ps->mat[PS_MAT_A];

  PetscScalar    *Q = ps->mat[PS_MAT_Q];

  PetscScalar    *X,*Y;

  const char     *side;

  PetscFunctionBegin;

  if (ps->state<PS_STATE_CONDENSED) SETERRQ(((PetscObject)ps)->comm,PETSC_ERR_ORDER,"Must call PSSolve() first");

  n  = PetscBLASIntCast(ps->n);

  ld = PetscBLASIntCast(ps->ld);

  if (left) {

    X = PETSC_NULL;

    Y = ps->mat[PS_MAT_Y];

    side = "L";

    ierr = PetscMemcpy(Y,Q,ld*ld*sizeof(PetscScalar));CHKERRQ(ierr);

  } else {

    X = ps->mat[PS_MAT_X];

    Y = PETSC_NULL;

    side = "R";

    ierr = PetscMemcpy(X,Q,ld*ld*sizeof(PetscScalar));CHKERRQ(ierr);

  }

#if !defined(PETSC_USE_COMPLEX)

  ierr = PSAllocateWork_Private(ps,3*ld,0,0);CHKERRQ(ierr);

  LAPACKtrevc_(side,"B",PETSC_NULL,&n,A,&ld,Y,&ld,X,&ld,&n,&mout,ps->work,&info);

#else

  ierr = PSAllocateWork_Private(ps,2*ld,ld,0);CHKERRQ(ierr);

  LAPACKtrevc_(side,"B",PETSC_NULL,&n,A,&ld,Y,&ld,X,&ld,&n,&mout,ps->work,ps->rwork,&info);

#endif

  if (info) SETERRQ1(((PetscObject)ps)->comm,PETSC_ERR_LIB,"Error in Lapack xTREVC %i",info);

  PetscFunctionReturn(0);

#endif

}

#undef __FUNCT__  

#define __FUNCT__ "PSVectors_NHEP"

PetscErrorCode PSVectors_NHEP(PS ps,PSMatType mat,PetscInt *k,PetscReal *rnorm)

{

  PetscErrorCode ierr;

  PetscFunctionBegin;

  switch (mat) {

    case PS_MAT_X:

      if (k) {

        ierr = PSVectors_NHEP_Eigen_Some(ps,k,rnorm,PETSC_FALSE);CHKERRQ(ierr);

      } else {

        ierr = PSVectors_NHEP_Eigen_All(ps,PETSC_FALSE);CHKERRQ(ierr);

      }

      break;

    case PS_MAT_Y:

      if (k) {

        ierr = PSVectors_NHEP_Eigen_Some(ps,k,rnorm,PETSC_TRUE);CHKERRQ(ierr);

      } else {

        ierr = PSVectors_NHEP_Eigen_All(ps,PETSC_TRUE);CHKERRQ(ierr);

      }

      break;

    case PS_MAT_U:

    case PS_MAT_VT:

      SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented yet");

      break;

    default:

      SETERRQ(((PetscObject)ps)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");

  }

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "PSSolve_NHEP"

PetscErrorCode PSSolve_NHEP(PS ps,PetscScalar *wr,PetscScalar *wi)

{

#if defined(SLEPC_MISSING_LAPACK_GEHRD) || defined(SLEPC_MISSING_LAPACK_ORGHR) || defined(PETSC_MISSING_LAPACK_HSEQR)

  PetscFunctionBegin;

  SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GEHRD/ORGHR/HSEQR - Lapack routines are unavailable.");

#else

  PetscErrorCode ierr;

  PetscScalar    *work,*tau;

  PetscInt       i,j;

  PetscBLASInt   ilo,lwork,info,n,ld;

  PetscScalar    *A = ps->mat[PS_MAT_A];

  PetscScalar    *Q = ps->mat[PS_MAT_Q];

  PetscFunctionBegin;

  PetscValidPointer(wi,2);

  n   = PetscBLASIntCast(ps->n);

  ld  = PetscBLASIntCast(ps->ld);

  ilo = PetscBLASIntCast(ps->l+1);

  ierr = PSAllocateWork_Private(ps,ld+ld*ld,0,0);CHKERRQ(ierr);

  tau  = ps->work;

  work = ps->work+ld;

  lwork = ld*ld;

  /* initialize orthogonal matrix */

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

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

    Q[i+i*ld] = 1.0;

  if (n==1) PetscFunctionReturn(0);    

  /* reduce to upper Hessenberg form */

  if (ps->state<PS_STATE_INTERMEDIATE) {

    LAPACKgehrd_(&n,&ilo,&n,A,&ld,tau,work,&lwork,&info);

    if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in Lapack xGEHRD %d",info);

    for (j=0;j<n-1;j++) {

      for (i=j+2;i<n;i++) {

        Q[i+j*ld] = A[i+j*ld];

        A[i+j*ld] = 0.0;

      }

    }

    LAPACKorghr_(&n,&ilo,&n,Q,&ld,tau,work,&lwork,&info);

    if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in Lapack xORGHR %d",info);

  }

  /* compute the (real) Schur form */

#if !defined(PETSC_USE_COMPLEX)

  LAPACKhseqr_("S","V",&n,&ilo,&n,A,&ld,wr,wi,Q,&ld,work,&lwork,&info);

  for (j=0;j<ps->l;j++) {

    if (j==n-1 || A[j+1+j*ld] == 0.0) {

      /* real eigenvalue */

      wr[j] = A[j+j*ld];

      wi[j] = 0.0;

    } else {

      /* complex eigenvalue */

      wr[j] = A[j+j*ld];

      wr[j+1] = A[j+j*ld];

      wi[j] = PetscSqrtReal(PetscAbsReal(A[j+1+j*ld])) *

              PetscSqrtReal(PetscAbsReal(A[j+(j+1)*ld]));

      wi[j+1] = -wi[j];

      j++;

    }

  }

#else

  LAPACKhseqr_("S","V",&n,&ilo,&n,A,&ld,wr,Q,&ld,work,&lwork,&info);

#endif

  if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in Lapack xHSEQR %d",info);

  PetscFunctionReturn(0);

#endif

}

#undef __FUNCT__  

#define __FUNCT__ "PSSort_NHEP"

PetscErrorCode PSSort_NHEP(PS ps,PetscScalar *wr,PetscScalar *wi,PetscErrorCode (*comp_func)(PetscScalar,PetscScalar,PetscScalar,PetscScalar,PetscInt*,void*),void *comp_ctx)

{

#if defined(SLEPC_MISSING_LAPACK_TREXC)

  PetscFunctionBegin;

  SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"TREXC - Lapack routine is unavailable.");

#else

  PetscErrorCode ierr;

  PetscScalar    re,im;

  PetscInt       i,j,pos,result;

  PetscBLASInt   ifst,ilst,info,n,ld;

  PetscScalar    *T = ps->mat[PS_MAT_A];

  PetscScalar    *Q = ps->mat[PS_MAT_Q];

#if !defined(PETSC_USE_COMPLEX)

  PetscScalar    *work;

#endif

  PetscFunctionBegin;

  PetscValidPointer(wi,2);

  n  = PetscBLASIntCast(ps->n);

  ld = PetscBLASIntCast(ps->ld);

#if !defined(PETSC_USE_COMPLEX)

  ierr = PSAllocateWork_Private(ps,ld,0,0);CHKERRQ(ierr);

  work = ps->work;

#endif

  /* selection sort */

  for (i=ps->l;i<n-1;i++) {

    re = wr[i];

    im = wi[i];

    pos = 0;

    j=i+1; /* j points to the next eigenvalue */

#if !defined(PETSC_USE_COMPLEX)

    if (im != 0) j=i+2;

#endif

    /* find minimum eigenvalue */

    for (;j<n;j++) {

      ierr = (*comp_func)(re,im,wr[j],wi[j],&result,comp_ctx);CHKERRQ(ierr);

      if (result > 0) {

        re = wr[j];

        im = wi[j];

        pos = j;

      }

#if !defined(PETSC_USE_COMPLEX)

      if (wi[j] != 0) j++;

#endif

    }

    if (pos) {

      /* interchange blocks */

      ifst = PetscBLASIntCast(pos+1);

      ilst = PetscBLASIntCast(i+1);

#if !defined(PETSC_USE_COMPLEX)

      LAPACKtrexc_("V",&n,T,&ld,Q,&ld,&ifst,&ilst,work,&info);

#else

      LAPACKtrexc_("V",&n,T,&ld,Q,&ld,&ifst,&ilst,&info);

#endif

      if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in Lapack xTREXC %d",info);

      /* recover original eigenvalues from T matrix */

      for (j=i;j<n;j++) {

        wr[j] = T[j+j*ld];

#if !defined(PETSC_USE_COMPLEX)

        if (j<n-1 && T[j+1+j*ld] != 0.0) {

          /* complex conjugate eigenvalue */

          wi[j] = PetscSqrtReal(PetscAbsReal(T[j+1+j*ld])) *

                  PetscSqrtReal(PetscAbsReal(T[j+(j+1)*ld]));

          wr[j+1] = wr[j];

          wi[j+1] = -wi[j];

          j++;

        } else

#endif

        wi[j] = 0.0;

      }

    }

#if !defined(PETSC_USE_COMPLEX)

    if (wi[i] != 0) i++;

#endif

  }

  PetscFunctionReturn(0);

#endif 

}

#undef __FUNCT__  

#define __FUNCT__ "PSCond_NHEP"

PetscErrorCode PSCond_NHEP(PS ps,PetscReal *cond)

{

#if defined(PETSC_MISSING_LAPACK_GETRF) || defined(SLEPC_MISSING_LAPACK_GETRI) || defined(SLEPC_MISSING_LAPACK_LANGE) || defined(SLEPC_MISSING_LAPACK_LANHS)

  PetscFunctionBegin;

  SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRF/GETRI/LANGE/LANHS - Lapack routines are unavailable.");

#else

  PetscErrorCode ierr;

  PetscScalar    *work;

  PetscReal      *rwork;

  PetscBLASInt   *ipiv;

  PetscBLASInt   lwork,info,n,ld;

  PetscReal      hn,hin;

  PetscScalar    *A;

  PetscFunctionBegin;

  n  = PetscBLASIntCast(ps->n);

  ld = PetscBLASIntCast(ps->ld);

  lwork = 8*ld;

  ierr = PSAllocateWork_Private(ps,lwork,ld,ld);CHKERRQ(ierr);

  work  = ps->work;

  rwork = ps->rwork;

  ipiv  = ps->iwork;

  /* use workspace matrix W to avoid overwriting A */

  ierr = PSAllocateMat_Private(ps,PS_MAT_W);CHKERRQ(ierr);

  A = ps->mat[PS_MAT_W];

  ierr = PetscMemcpy(A,ps->mat[PS_MAT_A],sizeof(PetscScalar)*ps->ld*ps->ld);CHKERRQ(ierr);

  /* norm of A */

  if (ps->state<PS_STATE_INTERMEDIATE) hn = LAPACKlange_("I",&n,&n,A,&ld,rwork);

  else hn = LAPACKlanhs_("I",&n,A,&ld,rwork);

  /* norm of inv(A) */

  LAPACKgetrf_(&n,&n,A,&ld,ipiv,&info);

  if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in Lapack xGETRF %d",info);

  LAPACKgetri_(&n,A,&ld,ipiv,work,&lwork,&info);

  if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in Lapack xGETRI %d",info);

  hin = LAPACKlange_("I",&n,&n,A,&ld,rwork);

  *cond = hn*hin;

  PetscFunctionReturn(0);

#endif

}

#undef __FUNCT__  

#define __FUNCT__ "PSTranslateHarmonic_NHEP"

PetscErrorCode PSTranslateHarmonic_NHEP(PS ps,PetscScalar tau,PetscReal beta,PetscBool recover,PetscScalar *gin,PetscReal *gamma)

{

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

  PetscFunctionBegin;

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

#else

  PetscErrorCode ierr;

  PetscInt       i,j;

  PetscBLASInt   *ipiv,info,n,ld,one=1,ncol;

  PetscScalar    *A,*B,*Q,*g=gin,*ghat;

  PetscScalar    done=1.0,dmone=-1.0,dzero=0.0;

  PetscReal      gnorm;

  PetscFunctionBegin;

  n  = PetscBLASIntCast(ps->n);

  ld = PetscBLASIntCast(ps->ld);

  A  = ps->mat[PS_MAT_A];

  if (!recover) {

    ierr = PSAllocateWork_Private(ps,0,0,ld);CHKERRQ(ierr);

    ipiv = ps->iwork;

    if (!g) {

      ierr = PSAllocateWork_Private(ps,ld,0,0);CHKERRQ(ierr);

      g = ps->work;

    }

    /* use workspace matrix W to factor A-tau*eye(n) */

    ierr = PSAllocateMat_Private(ps,PS_MAT_W);CHKERRQ(ierr);

    B = ps->mat[PS_MAT_W];

    ierr = PetscMemcpy(B,A,sizeof(PetscScalar)*ld*ld);CHKERRQ(ierr);

    /* Vector g initialy stores b = beta*e_n^T */

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

    g[n-1] = beta;

    /* g = (A-tau*eye(n))'\b */

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

      B[i+i*ld] -= tau;

    LAPACKgetrf_(&n,&n,B,&ld,ipiv,&info);

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

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

    ierr = PetscLogFlops(2.0*n*n*n/3.0);CHKERRQ(ierr);

    LAPACKgetrs_("C",&n,&one,B,&ld,ipiv,g,&ld,&info);

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

    ierr = PetscLogFlops(2.0*n*n-n);CHKERRQ(ierr);

    /* A = A + g*b' */

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

      A[i+(n-1)*ld] += g[i]*beta;

  } else { /* recover */

    PetscValidPointer(g,6);

    ierr = PSAllocateWork_Private(ps,ld,0,0);CHKERRQ(ierr);

    ghat = ps->work;

    Q    = ps->mat[PS_MAT_Q];

    /* g^ = -Q(:,idx)'*g */

    ncol = PetscBLASIntCast(ps->l+ps->k);

    BLASgemv_("C",&n,&ncol,&dmone,Q,&ld,g,&one,&dzero,ghat,&one);

    /* A = A + g^*b' */

    for (i=0;i<ps->l+ps->k;i++)

      for (j=ps->l;j<ps->l+ps->k;j++)

        A[i+j*ld] += ghat[i]*Q[n-1+j*ld]*beta;

    /* g~ = (I-Q(:,idx)*Q(:,idx)')*g = g+Q(:,idx)*g^ */

    BLASgemv_("N",&n,&ncol,&done,Q,&ld,ghat,&one,&done,g,&one);

  }

  /* Compute gamma factor */

  if (gamma) {

    gnorm = 0.0;

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

      gnorm = gnorm + PetscRealPart(g[i]*PetscConj(g[i]));

    *gamma = PetscSqrtReal(1.0+gnorm);

  }

  PetscFunctionReturn(0);

#endif

}

EXTERN_C_BEGIN

#undef __FUNCT__  

#define __FUNCT__ "PSCreate_NHEP"

PetscErrorCode PSCreate_NHEP(PS ps)

{

  PetscFunctionBegin;

  ps->nmeth  = 1;

  ps->ops->allocate      = PSAllocate_NHEP;

  ps->ops->view          = PSView_NHEP;

  ps->ops->vectors       = PSVectors_NHEP;

  ps->ops->solve         = PSSolve_NHEP;

  ps->ops->sort          = PSSort_NHEP;

  ps->ops->cond          = PSCond_NHEP;

  ps->ops->translate     = PSTranslateHarmonic_NHEP;

  PetscFunctionReturn(0);

}

EXTERN_C_END