/*
|
/*
|
This file contains routines for handling small-size dense problems.
|
This file contains routines for handling small-size dense problems.
|
All routines are simply wrappers to LAPACK routines. Matrices passed in
|
All routines are simply wrappers to LAPACK routines. Matrices passed in
|
as arguments are assumed to be square matrices stored in column-major
|
as arguments are assumed to be square matrices stored in column-major
|
format with a leading dimension equal to the number of rows.
|
format with a leading dimension equal to the number of rows.
|
|
|
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
SLEPc - Scalable Library for Eigenvalue Problem Computations
|
SLEPc - Scalable Library for Eigenvalue Problem Computations
|
Copyright (c) 2002-2007, Universidad Politecnica de Valencia, Spain
|
Copyright (c) 2002-2009, Universidad Politecnica de Valencia, Spain
|
|
|
This file is part of SLEPc. See the README file for conditions of use
|
This file is part of SLEPc.
|
and additional information.
|
|
|
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 "private/epsimpl.h" /*I "slepceps.h" I*/
|
#include "private/epsimpl.h" /*I "slepceps.h" I*/
|
#include "slepcblaslapack.h"
|
#include "slepcblaslapack.h"
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "EPSDenseNHEP"
|
#define __FUNCT__ "EPSDenseNHEP"
|
/*@
|
/*@
|
EPSDenseNHEP - Solves a dense standard non-Hermitian Eigenvalue Problem.
|
EPSDenseNHEP - Solves a dense standard non-Hermitian Eigenvalue Problem.
|
|
|
Not Collective
|
Not Collective
|
|
|
Input Parameters:
|
Input Parameters:
|
+ n - dimension of the eigenproblem
|
+ n - dimension of the eigenproblem
|
- A - pointer to the array containing the matrix values
|
- A - pointer to the array containing the matrix values
|
|
|
Output Parameters:
|
Output Parameters:
|
+ w - pointer to the array to store the computed eigenvalues
|
+ w - pointer to the array to store the computed eigenvalues
|
. wi - imaginary part of the eigenvalues (only when using real numbers)
|
. wi - imaginary part of the eigenvalues (only when using real numbers)
|
. V - pointer to the array to store right eigenvectors
|
. V - pointer to the array to store right eigenvectors
|
- W - pointer to the array to store left eigenvectors
|
- W - pointer to the array to store left eigenvectors
|
|
|
Notes:
|
Notes:
|
If either V or W are PETSC_NULL then the corresponding eigenvectors are
|
If either V or W are PETSC_NULL then the corresponding eigenvectors are
|
not computed.
|
not computed.
|
|
|
Matrix A is overwritten.
|
Matrix A is overwritten.
|
|
|
This routine uses LAPACK routines xGEEVX.
|
This routine uses LAPACK routines xGEEVX.
|
|
|
Level: developer
|
Level: developer
|
|
|
.seealso: EPSDenseGNHEP(), EPSDenseHEP(), EPSDenseGHEP()
|
.seealso: EPSDenseGNHEP(), EPSDenseHEP(), EPSDenseGHEP()
|
@*/
|
@*/
|
PetscErrorCode EPSDenseNHEP(PetscInt n_,PetscScalar *A,PetscScalar *w,PetscScalar *wi,PetscScalar *V,PetscScalar *W)
|
PetscErrorCode EPSDenseNHEP(PetscInt n_,PetscScalar *A,PetscScalar *w,PetscScalar *wi,PetscScalar *V,PetscScalar *W)
|
{
|
{
|
#if defined(SLEPC_MISSING_LAPACK_GEEVX)
|
#if defined(SLEPC_MISSING_LAPACK_GEEVX)
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
SETERRQ(PETSC_ERR_SUP,"GEEVX - Lapack routine is unavailable.");
|
SETERRQ(PETSC_ERR_SUP,"GEEVX - Lapack routine is unavailable.");
|
#else
|
#else
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
PetscReal abnrm,*scale,dummy;
|
PetscReal abnrm,*scale,dummy;
|
PetscScalar *work;
|
PetscScalar *work;
|
PetscBLASInt ilo,ihi,n,lwork,info;
|
PetscBLASInt ilo,ihi,n,lwork,info;
|
const char *jobvr,*jobvl;
|
const char *jobvr,*jobvl;
|
#if defined(PETSC_USE_COMPLEX)
|
#if defined(PETSC_USE_COMPLEX)
|
PetscReal *rwork;
|
PetscReal *rwork;
|
#else
|
#else
|
PetscBLASInt idummy;
|
PetscBLASInt idummy;
|
#endif
|
#endif
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
ierr = PetscLogEventBegin(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
ierr = PetscLogEventBegin(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
n = PetscBLASIntCast(n_);
|
n = PetscBLASIntCast(n_);
|
lwork = PetscBLASIntCast(4*n_);
|
lwork = PetscBLASIntCast(4*n_);
|
if (V) jobvr = "V";
|
if (V) jobvr = "V";
|
else jobvr = "N";
|
else jobvr = "N";
|
if (W) jobvl = "V";
|
if (W) jobvl = "V";
|
else jobvl = "N";
|
else jobvl = "N";
|
ierr = PetscMalloc(lwork*sizeof(PetscScalar),&work);CHKERRQ(ierr);
|
ierr = PetscMalloc(lwork*sizeof(PetscScalar),&work);CHKERRQ(ierr);
|
ierr = PetscMalloc(n*sizeof(PetscReal),&scale);CHKERRQ(ierr);
|
ierr = PetscMalloc(n*sizeof(PetscReal),&scale);CHKERRQ(ierr);
|
#if defined(PETSC_USE_COMPLEX)
|
#if defined(PETSC_USE_COMPLEX)
|
ierr = PetscMalloc(2*n*sizeof(PetscReal),&rwork);CHKERRQ(ierr);
|
ierr = PetscMalloc(2*n*sizeof(PetscReal),&rwork);CHKERRQ(ierr);
|
LAPACKgeevx_("B",jobvl,jobvr,"N",&n,A,&n,w,W,&n,V,&n,&ilo,&ihi,scale,&abnrm,&dummy,&dummy,work,&lwork,rwork,&info);
|
LAPACKgeevx_("B",jobvl,jobvr,"N",&n,A,&n,w,W,&n,V,&n,&ilo,&ihi,scale,&abnrm,&dummy,&dummy,work,&lwork,rwork,&info);
|
if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack ZGEEVX %d",info);
|
if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack ZGEEVX %d",info);
|
ierr = PetscFree(rwork);CHKERRQ(ierr);
|
ierr = PetscFree(rwork);CHKERRQ(ierr);
|
#else
|
#else
|
LAPACKgeevx_("B",jobvl,jobvr,"N",&n,A,&n,w,wi,W,&n,V,&n,&ilo,&ihi,scale,&abnrm,&dummy,&dummy,work,&lwork,&idummy,&info);
|
LAPACKgeevx_("B",jobvl,jobvr,"N",&n,A,&n,w,wi,W,&n,V,&n,&ilo,&ihi,scale,&abnrm,&dummy,&dummy,work,&lwork,&idummy,&info);
|
if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack DGEEVX %d",info);
|
if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack DGEEVX %d",info);
|
#endif
|
#endif
|
ierr = PetscFree(work);CHKERRQ(ierr);
|
ierr = PetscFree(work);CHKERRQ(ierr);
|
ierr = PetscFree(scale);CHKERRQ(ierr);
|
ierr = PetscFree(scale);CHKERRQ(ierr);
|
ierr = PetscLogEventEnd(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
ierr = PetscLogEventEnd(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
#endif
|
#endif
|
}
|
}
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "EPSDenseGNHEP"
|
#define __FUNCT__ "EPSDenseGNHEP"
|
/*@
|
/*@
|
EPSDenseGNHEP - Solves a dense Generalized non-Hermitian Eigenvalue Problem.
|
EPSDenseGNHEP - Solves a dense Generalized non-Hermitian Eigenvalue Problem.
|
|
|
Not Collective
|
Not Collective
|
|
|
Input Parameters:
|
Input Parameters:
|
+ n - dimension of the eigenproblem
|
+ n - dimension of the eigenproblem
|
. A - pointer to the array containing the matrix values for A
|
. A - pointer to the array containing the matrix values for A
|
- B - pointer to the array containing the matrix values for B
|
- B - pointer to the array containing the matrix values for B
|
|
|
Output Parameters:
|
Output Parameters:
|
+ w - pointer to the array to store the computed eigenvalues
|
+ w - pointer to the array to store the computed eigenvalues
|
. wi - imaginary part of the eigenvalues (only when using real numbers)
|
. wi - imaginary part of the eigenvalues (only when using real numbers)
|
. V - pointer to the array to store right eigenvectors
|
. V - pointer to the array to store right eigenvectors
|
- W - pointer to the array to store left eigenvectors
|
- W - pointer to the array to store left eigenvectors
|
|
|
Notes:
|
Notes:
|
If either V or W are PETSC_NULL then the corresponding eigenvectors are
|
If either V or W are PETSC_NULL then the corresponding eigenvectors are
|
not computed.
|
not computed.
|
|
|
Matrices A and B are overwritten.
|
Matrices A and B are overwritten.
|
|
|
This routine uses LAPACK routines xGGEVX.
|
This routine uses LAPACK routines xGGEVX.
|
|
|
Level: developer
|
Level: developer
|
|
|
.seealso: EPSDenseNHEP(), EPSDenseHEP(), EPSDenseGHEP()
|
.seealso: EPSDenseNHEP(), EPSDenseHEP(), EPSDenseGHEP()
|
@*/
|
@*/
|
PetscErrorCode EPSDenseGNHEP(PetscInt n_,PetscScalar *A,PetscScalar *B,PetscScalar *w,PetscScalar *wi,PetscScalar *V,PetscScalar *W)
|
PetscErrorCode EPSDenseGNHEP(PetscInt n_,PetscScalar *A,PetscScalar *B,PetscScalar *w,PetscScalar *wi,PetscScalar *V,PetscScalar *W)
|
{
|
{
|
#if defined(SLEPC_MISSING_LAPACK_GGEVX)
|
#if defined(SLEPC_MISSING_LAPACK_GGEVX)
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
SETERRQ(PETSC_ERR_SUP,"GGEVX - Lapack routine is unavailable.");
|
SETERRQ(PETSC_ERR_SUP,"GGEVX - Lapack routine is unavailable.");
|
#else
|
#else
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
PetscReal *rscale,*lscale,abnrm,bbnrm,dummy;
|
PetscReal *rscale,*lscale,abnrm,bbnrm,dummy;
|
PetscScalar *alpha,*beta,*work;
|
PetscScalar *alpha,*beta,*work;
|
PetscInt i;
|
PetscInt i;
|
PetscBLASInt ilo,ihi,idummy,info,n;
|
PetscBLASInt ilo,ihi,idummy,info,n;
|
const char *jobvr,*jobvl;
|
const char *jobvr,*jobvl;
|
#if defined(PETSC_USE_COMPLEX)
|
#if defined(PETSC_USE_COMPLEX)
|
PetscReal *rwork;
|
PetscReal *rwork;
|
PetscBLASInt lwork;
|
PetscBLASInt lwork;
|
#else
|
#else
|
PetscReal *alphai;
|
PetscReal *alphai;
|
PetscBLASInt lwork;
|
PetscBLASInt lwork;
|
#endif
|
#endif
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
ierr = PetscLogEventBegin(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
ierr = PetscLogEventBegin(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
n = PetscBLASIntCast(n_);
|
n = PetscBLASIntCast(n_);
|
#if defined(PETSC_USE_COMPLEX)
|
#if defined(PETSC_USE_COMPLEX)
|
lwork = PetscBLASIntCast(2*n_);
|
lwork = PetscBLASIntCast(2*n_);
|
#else
|
#else
|
lwork = PetscBLASIntCast(6*n_);
|
lwork = PetscBLASIntCast(6*n_);
|
#endif
|
#endif
|
if (V) jobvr = "V";
|
if (V) jobvr = "V";
|
else jobvr = "N";
|
else jobvr = "N";
|
if (W) jobvl = "V";
|
if (W) jobvl = "V";
|
else jobvl = "N";
|
else jobvl = "N";
|
ierr = PetscMalloc(n*sizeof(PetscScalar),&alpha);CHKERRQ(ierr);
|
ierr = PetscMalloc(n*sizeof(PetscScalar),&alpha);CHKERRQ(ierr);
|
ierr = PetscMalloc(n*sizeof(PetscScalar),&beta);CHKERRQ(ierr);
|
ierr = PetscMalloc(n*sizeof(PetscScalar),&beta);CHKERRQ(ierr);
|
ierr = PetscMalloc(n*sizeof(PetscReal),&rscale);CHKERRQ(ierr);
|
ierr = PetscMalloc(n*sizeof(PetscReal),&rscale);CHKERRQ(ierr);
|
ierr = PetscMalloc(n*sizeof(PetscReal),&lscale);CHKERRQ(ierr);
|
ierr = PetscMalloc(n*sizeof(PetscReal),&lscale);CHKERRQ(ierr);
|
ierr = PetscMalloc(lwork*sizeof(PetscScalar),&work);CHKERRQ(ierr);
|
ierr = PetscMalloc(lwork*sizeof(PetscScalar),&work);CHKERRQ(ierr);
|
#if defined(PETSC_USE_COMPLEX)
|
#if defined(PETSC_USE_COMPLEX)
|
ierr = PetscMalloc(6*n*sizeof(PetscReal),&rwork);CHKERRQ(ierr);
|
ierr = PetscMalloc(6*n*sizeof(PetscReal),&rwork);CHKERRQ(ierr);
|
LAPACKggevx_("B",jobvl,jobvr,"N",&n,A,&n,B,&n,alpha,beta,W,&n,V,&n,&ilo,&ihi,lscale,rscale,&abnrm,&bbnrm,&dummy,&dummy,work,&lwork,rwork,&idummy,&idummy,&info);
|
LAPACKggevx_("B",jobvl,jobvr,"N",&n,A,&n,B,&n,alpha,beta,W,&n,V,&n,&ilo,&ihi,lscale,rscale,&abnrm,&bbnrm,&dummy,&dummy,work,&lwork,rwork,&idummy,&idummy,&info);
|
if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack ZGGEVX %d",info);
|
if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack ZGGEVX %d",info);
|
for (i=0;i<n;i++) {
|
for (i=0;i<n;i++) {
|
w[i] = alpha[i]/beta[i];
|
w[i] = alpha[i]/beta[i];
|
}
|
}
|
ierr = PetscFree(rwork);CHKERRQ(ierr);
|
ierr = PetscFree(rwork);CHKERRQ(ierr);
|
#else
|
#else
|
ierr = PetscMalloc(n*sizeof(PetscReal),&alphai);CHKERRQ(ierr);
|
ierr = PetscMalloc(n*sizeof(PetscReal),&alphai);CHKERRQ(ierr);
|
LAPACKggevx_("B",jobvl,jobvr,"N",&n,A,&n,B,&n,alpha,alphai,beta,W,&n,V,&n,&ilo,&ihi,lscale,rscale,&abnrm,&bbnrm,&dummy,&dummy,work,&lwork,&idummy,&idummy,&info);
|
LAPACKggevx_("B",jobvl,jobvr,"N",&n,A,&n,B,&n,alpha,alphai,beta,W,&n,V,&n,&ilo,&ihi,lscale,rscale,&abnrm,&bbnrm,&dummy,&dummy,work,&lwork,&idummy,&idummy,&info);
|
if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack DGGEVX %d",info);
|
if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack DGGEVX %d",info);
|
for (i=0;i<n;i++) {
|
for (i=0;i<n;i++) {
|
w[i] = alpha[i]/beta[i];
|
w[i] = alpha[i]/beta[i];
|
wi[i] = alphai[i]/beta[i];
|
wi[i] = alphai[i]/beta[i];
|
}
|
}
|
ierr = PetscFree(alphai);CHKERRQ(ierr);
|
ierr = PetscFree(alphai);CHKERRQ(ierr);
|
#endif
|
#endif
|
ierr = PetscFree(alpha);CHKERRQ(ierr);
|
ierr = PetscFree(alpha);CHKERRQ(ierr);
|
ierr = PetscFree(beta);CHKERRQ(ierr);
|
ierr = PetscFree(beta);CHKERRQ(ierr);
|
ierr = PetscFree(rscale);CHKERRQ(ierr);
|
ierr = PetscFree(rscale);CHKERRQ(ierr);
|
ierr = PetscFree(lscale);CHKERRQ(ierr);
|
ierr = PetscFree(lscale);CHKERRQ(ierr);
|
ierr = PetscFree(work);CHKERRQ(ierr);
|
ierr = PetscFree(work);CHKERRQ(ierr);
|
ierr = PetscLogEventEnd(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
ierr = PetscLogEventEnd(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
#endif
|
#endif
|
}
|
}
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "EPSDenseHEP"
|
#define __FUNCT__ "EPSDenseHEP"
|
/*@
|
/*@
|
EPSDenseHEP - Solves a dense standard Hermitian Eigenvalue Problem.
|
EPSDenseHEP - Solves a dense standard Hermitian Eigenvalue Problem.
|
|
|
Not Collective
|
Not Collective
|
|
|
Input Parameters:
|
Input Parameters:
|
+ n - dimension of the eigenproblem
|
+ n - dimension of the eigenproblem
|
. A - pointer to the array containing the matrix values
|
. A - pointer to the array containing the matrix values
|
- lda - leading dimension of A
|
- lda - leading dimension of A
|
|
|
Output Parameters:
|
Output Parameters:
|
+ w - pointer to the array to store the computed eigenvalues
|
+ w - pointer to the array to store the computed eigenvalues
|
- V - pointer to the array to store the eigenvectors
|
- V - pointer to the array to store the eigenvectors
|
|
|
Notes:
|
Notes:
|
If V is PETSC_NULL then the eigenvectors are not computed.
|
If V is PETSC_NULL then the eigenvectors are not computed.
|
|
|
Matrix A is overwritten.
|
Matrix A is overwritten.
|
|
|
This routine uses LAPACK routines DSYEVR or ZHEEVR.
|
This routine uses LAPACK routines DSYEVR or ZHEEVR.
|
|
|
Level: developer
|
Level: developer
|
|
|
.seealso: EPSDenseNHEP(), EPSDenseGNHEP(), EPSDenseGHEP()
|
.seealso: EPSDenseNHEP(), EPSDenseGNHEP(), EPSDenseGHEP()
|
@*/
|
@*/
|
PetscErrorCode EPSDenseHEP(PetscInt n_,PetscScalar *A,PetscInt lda_,PetscReal *w,PetscScalar *V)
|
PetscErrorCode EPSDenseHEP(PetscInt n_,PetscScalar *A,PetscInt lda_,PetscReal *w,PetscScalar *V)
|
{
|
{
|
#if defined(SLEPC_MISSING_LAPACK_SYEVR) || defined(SLEPC_MISSING_LAPACK_HEEVR)
|
#if defined(SLEPC_MISSING_LAPACK_SYEVR) || defined(SLEPC_MISSING_LAPACK_HEEVR)
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
SETERRQ(PETSC_ERR_SUP,"DSYEVR/ZHEEVR - Lapack routine is unavailable.");
|
SETERRQ(PETSC_ERR_SUP,"DSYEVR/ZHEEVR - Lapack routine is unavailable.");
|
#else
|
#else
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
PetscReal abstol = 0.0,vl,vu;
|
PetscReal abstol = 0.0,vl,vu;
|
PetscScalar *work;
|
PetscScalar *work;
|
PetscBLASInt il,iu,m,*isuppz,*iwork,n,lda,liwork,info;
|
PetscBLASInt il,iu,m,*isuppz,*iwork,n,lda,liwork,info;
|
const char *jobz;
|
const char *jobz;
|
#if defined(PETSC_USE_COMPLEX)
|
#if defined(PETSC_USE_COMPLEX)
|
PetscReal *rwork;
|
PetscReal *rwork;
|
PetscBLASInt lwork,lrwork;
|
PetscBLASInt lwork,lrwork;
|
#else
|
#else
|
PetscBLASInt lwork;
|
PetscBLASInt lwork;
|
#endif
|
#endif
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
ierr = PetscLogEventBegin(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
ierr = PetscLogEventBegin(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
n = PetscBLASIntCast(n_);
|
n = PetscBLASIntCast(n_);
|
lda = PetscBLASIntCast(lda_);
|
lda = PetscBLASIntCast(lda_);
|
liwork = PetscBLASIntCast(10*n_);
|
liwork = PetscBLASIntCast(10*n_);
|
#if defined(PETSC_USE_COMPLEX)
|
#if defined(PETSC_USE_COMPLEX)
|
lwork = PetscBLASIntCast(18*n_);
|
lwork = PetscBLASIntCast(18*n_);
|
lrwork = PetscBLASIntCast(24*n_);
|
lrwork = PetscBLASIntCast(24*n_);
|
#else
|
#else
|
lwork = PetscBLASIntCast(26*n_);
|
lwork = PetscBLASIntCast(26*n_);
|
#endif
|
#endif
|
if (V) jobz = "V";
|
if (V) jobz = "V";
|
else jobz = "N";
|
else jobz = "N";
|
ierr = PetscMalloc(2*n*sizeof(PetscBLASInt),&isuppz);CHKERRQ(ierr);
|
ierr = PetscMalloc(2*n*sizeof(PetscBLASInt),&isuppz);CHKERRQ(ierr);
|
ierr = PetscMalloc(lwork*sizeof(PetscScalar),&work);CHKERRQ(ierr);
|
ierr = PetscMalloc(lwork*sizeof(PetscScalar),&work);CHKERRQ(ierr);
|
ierr = PetscMalloc(liwork*sizeof(PetscBLASInt),&iwork);CHKERRQ(ierr);
|
ierr = PetscMalloc(liwork*sizeof(PetscBLASInt),&iwork);CHKERRQ(ierr);
|
#if defined(PETSC_USE_COMPLEX)
|
#if defined(PETSC_USE_COMPLEX)
|
ierr = PetscMalloc(lrwork*sizeof(PetscReal),&rwork);CHKERRQ(ierr);
|
ierr = PetscMalloc(lrwork*sizeof(PetscReal),&rwork);CHKERRQ(ierr);
|
LAPACKsyevr_(jobz,"A","L",&n,A,&lda,&vl,&vu,&il,&iu,&abstol,&m,w,V,&n,isuppz,work,&lwork,rwork,&lrwork,iwork,&liwork,&info);
|
LAPACKsyevr_(jobz,"A","L",&n,A,&lda,&vl,&vu,&il,&iu,&abstol,&m,w,V,&n,isuppz,work,&lwork,rwork,&lrwork,iwork,&liwork,&info);
|
if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack ZHEEVR %d",info);
|
if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack ZHEEVR %d",info);
|
ierr = PetscFree(rwork);CHKERRQ(ierr);
|
ierr = PetscFree(rwork);CHKERRQ(ierr);
|
#else
|
#else
|
LAPACKsyevr_(jobz,"A","L",&n,A,&lda,&vl,&vu,&il,&iu,&abstol,&m,w,V,&n,isuppz,work,&lwork,iwork,&liwork,&info);
|
LAPACKsyevr_(jobz,"A","L",&n,A,&lda,&vl,&vu,&il,&iu,&abstol,&m,w,V,&n,isuppz,work,&lwork,iwork,&liwork,&info);
|
if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack DSYEVR %d",info);
|
if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack DSYEVR %d",info);
|
#endif
|
#endif
|
ierr = PetscFree(isuppz);CHKERRQ(ierr);
|
ierr = PetscFree(isuppz);CHKERRQ(ierr);
|
ierr = PetscFree(work);CHKERRQ(ierr);
|
ierr = PetscFree(work);CHKERRQ(ierr);
|
ierr = PetscFree(iwork);CHKERRQ(ierr);
|
ierr = PetscFree(iwork);CHKERRQ(ierr);
|
ierr = PetscLogEventEnd(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
ierr = PetscLogEventEnd(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
#endif
|
#endif
|
}
|
}
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "EPSDenseGHEP"
|
#define __FUNCT__ "EPSDenseGHEP"
|
/*@
|
/*@
|
EPSDenseGHEP - Solves a dense Generalized Hermitian Eigenvalue Problem.
|
EPSDenseGHEP - Solves a dense Generalized Hermitian Eigenvalue Problem.
|
|
|
Not Collective
|
Not Collective
|
|
|
Input Parameters:
|
Input Parameters:
|
+ n - dimension of the eigenproblem
|
+ n - dimension of the eigenproblem
|
. A - pointer to the array containing the matrix values for A
|
. A - pointer to the array containing the matrix values for A
|
- B - pointer to the array containing the matrix values for B
|
- B - pointer to the array containing the matrix values for B
|
|
|
Output Parameters:
|
Output Parameters:
|
+ w - pointer to the array to store the computed eigenvalues
|
+ w - pointer to the array to store the computed eigenvalues
|
- V - pointer to the array to store the eigenvectors
|
- V - pointer to the array to store the eigenvectors
|
|
|
Notes:
|
Notes:
|
If V is PETSC_NULL then the eigenvectors are not computed.
|
If V is PETSC_NULL then the eigenvectors are not computed.
|
|
|
Matrices A and B are overwritten.
|
Matrices A and B are overwritten.
|
|
|
This routine uses LAPACK routines DSYGVD or ZHEGVD.
|
This routine uses LAPACK routines DSYGVD or ZHEGVD.
|
|
|
Level: developer
|
Level: developer
|
|
|
.seealso: EPSDenseNHEP(), EPSDenseGNHEP(), EPSDenseHEP()
|
.seealso: EPSDenseNHEP(), EPSDenseGNHEP(), EPSDenseHEP()
|
@*/
|
@*/
|
PetscErrorCode EPSDenseGHEP(PetscInt n_,PetscScalar *A,PetscScalar *B,PetscReal *w,PetscScalar *V)
|
PetscErrorCode EPSDenseGHEP(PetscInt n_,PetscScalar *A,PetscScalar *B,PetscReal *w,PetscScalar *V)
|
{
|
{
|
#if defined(SLEPC_MISSING_LAPACK_SYGVD) || defined(SLEPC_MISSING_LAPACK_HEGVD)
|
#if defined(SLEPC_MISSING_LAPACK_SYGVD) || defined(SLEPC_MISSING_LAPACK_HEGVD)
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
SETERRQ(PETSC_ERR_SUP,"DSYGVD/ZHEGVD - Lapack routine is unavailable.");
|
SETERRQ(PETSC_ERR_SUP,"DSYGVD/ZHEGVD - Lapack routine is unavailable.");
|
#else
|
#else
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
PetscScalar *work;
|
PetscScalar *work;
|
PetscBLASInt itype = 1,*iwork,info,n,
|
PetscBLASInt itype = 1,*iwork,info,n,
|
liwork;
|
liwork;
|
const char *jobz;
|
const char *jobz;
|
#if defined(PETSC_USE_COMPLEX)
|
#if defined(PETSC_USE_COMPLEX)
|
PetscReal *rwork;
|
PetscReal *rwork;
|
PetscBLASInt lwork,lrwork;
|
PetscBLASInt lwork,lrwork;
|
#else
|
#else
|
PetscBLASInt lwork;
|
PetscBLASInt lwork;
|
#endif
|
#endif
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
ierr = PetscLogEventBegin(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
ierr = PetscLogEventBegin(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
n = PetscBLASIntCast(n_);
|
n = PetscBLASIntCast(n_);
|
if (V) {
|
if (V) {
|
jobz = "V";
|
jobz = "V";
|
liwork = PetscBLASIntCast(5*n_+3);
|
liwork = PetscBLASIntCast(5*n_+3);
|
#if defined(PETSC_USE_COMPLEX)
|
#if defined(PETSC_USE_COMPLEX)
|
lwork = PetscBLASIntCast(n_*n_+2*n_);
|
lwork = PetscBLASIntCast(n_*n_+2*n_);
|
lrwork = PetscBLASIntCast(2*n_*n_+5*n_+1);
|
lrwork = PetscBLASIntCast(2*n_*n_+5*n_+1);
|
#else
|
#else
|
lwork = PetscBLASIntCast(2*n_*n_+6*n_+1);
|
lwork = PetscBLASIntCast(2*n_*n_+6*n_+1);
|
#endif
|
#endif
|
} else {
|
} else {
|
jobz = "N";
|
jobz = "N";
|
liwork = 1;
|
liwork = 1;
|
#if defined(PETSC_USE_COMPLEX)
|
#if defined(PETSC_USE_COMPLEX)
|
lwork = PetscBLASIntCast(n_+1);
|
lwork = PetscBLASIntCast(n_+1);
|
lrwork = PetscBLASIntCast(n_);
|
lrwork = PetscBLASIntCast(n_);
|
#else
|
#else
|
lwork = PetscBLASIntCast(2*n_+1);
|
lwork = PetscBLASIntCast(2*n_+1);
|
#endif
|
#endif
|
}
|
}
|
ierr = PetscMalloc(lwork*sizeof(PetscScalar),&work);CHKERRQ(ierr);
|
ierr = PetscMalloc(lwork*sizeof(PetscScalar),&work);CHKERRQ(ierr);
|
ierr = PetscMalloc(liwork*sizeof(PetscBLASInt),&iwork);CHKERRQ(ierr);
|
ierr = PetscMalloc(liwork*sizeof(PetscBLASInt),&iwork);CHKERRQ(ierr);
|
#if defined(PETSC_USE_COMPLEX)
|
#if defined(PETSC_USE_COMPLEX)
|
ierr = PetscMalloc(lrwork*sizeof(PetscReal),&rwork);CHKERRQ(ierr);
|
ierr = PetscMalloc(lrwork*sizeof(PetscReal),&rwork);CHKERRQ(ierr);
|
LAPACKsygvd_(&itype,jobz,"U",&n,A,&n,B,&n,w,work,&lwork,rwork,&lrwork,iwork,&liwork,&info);
|
LAPACKsygvd_(&itype,jobz,"U",&n,A,&n,B,&n,w,work,&lwork,rwork,&lrwork,iwork,&liwork,&info);
|
if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack ZHEGVD %d",info);
|
if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack ZHEGVD %d",info);
|
ierr = PetscFree(rwork);CHKERRQ(ierr);
|
ierr = PetscFree(rwork);CHKERRQ(ierr);
|
#else
|
#else
|
LAPACKsygvd_(&itype,jobz,"U",&n,A,&n,B,&n,w,work,&lwork,iwork,&liwork,&info);
|
LAPACKsygvd_(&itype,jobz,"U",&n,A,&n,B,&n,w,work,&lwork,iwork,&liwork,&info);
|
if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack DSYGVD %d",info);
|
if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack DSYGVD %d",info);
|
#endif
|
#endif
|
if (V) {
|
if (V) {
|
ierr = PetscMemcpy(V,A,n*n*sizeof(PetscScalar));CHKERRQ(ierr);
|
ierr = PetscMemcpy(V,A,n*n*sizeof(PetscScalar));CHKERRQ(ierr);
|
}
|
}
|
ierr = PetscFree(work);CHKERRQ(ierr);
|
ierr = PetscFree(work);CHKERRQ(ierr);
|
ierr = PetscFree(iwork);CHKERRQ(ierr);
|
ierr = PetscFree(iwork);CHKERRQ(ierr);
|
ierr = PetscLogEventEnd(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
ierr = PetscLogEventEnd(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
#endif
|
#endif
|
}
|
}
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "EPSDenseHessenberg"
|
#define __FUNCT__ "EPSDenseHessenberg"
|
/*@
|
/*@
|
EPSDenseHessenberg - Computes the Hessenberg form of a dense matrix.
|
EPSDenseHessenberg - Computes the Hessenberg form of a dense matrix.
|
|
|
Not Collective
|
Not Collective
|
|
|
Input Parameters:
|
Input Parameters:
|
+ n - dimension of the matrix
|
+ n - dimension of the matrix
|
. k - first active column
|
. k - first active column
|
- lda - leading dimension of A
|
- lda - leading dimension of A
|
|
|
Input/Output Parameters:
|
Input/Output Parameters:
|
+ A - on entry, the full matrix; on exit, the upper Hessenberg matrix (H)
|
+ A - on entry, the full matrix; on exit, the upper Hessenberg matrix (H)
|
- Q - on exit, orthogonal matrix of vectors A = Q*H*Q'
|
- Q - on exit, orthogonal matrix of vectors A = Q*H*Q'
|
|
|
Notes:
|
Notes:
|
Only active columns (from k to n) are computed.
|
Only active columns (from k to n) are computed.
|
|
|
Both A and Q are overwritten.
|
Both A and Q are overwritten.
|
|
|
This routine uses LAPACK routines xGEHRD and xORGHR/xUNGHR.
|
This routine uses LAPACK routines xGEHRD and xORGHR/xUNGHR.
|
|
|
Level: developer
|
Level: developer
|
|
|
.seealso: EPSDenseSchur(), EPSSortDenseSchur(), EPSDenseTridiagonal()
|
.seealso: EPSDenseSchur(), EPSSortDenseSchur(), EPSDenseTridiagonal()
|
@*/
|
@*/
|
PetscErrorCode EPSDenseHessenberg(PetscInt n_,PetscInt k,PetscScalar *A,PetscInt lda_,PetscScalar *Q)
|
PetscErrorCode EPSDenseHessenberg(PetscInt n_,PetscInt k,PetscScalar *A,PetscInt lda_,PetscScalar *Q)
|
{
|
{
|
#if defined(SLEPC_MISSING_LAPACK_GEHRD) || defined(SLEPC_MISSING_LAPACK_ORGHR) || defined(SLEPC_MISSING_LAPACK_UNGHR)
|
#if defined(SLEPC_MISSING_LAPACK_GEHRD) || defined(SLEPC_MISSING_LAPACK_ORGHR) || defined(SLEPC_MISSING_LAPACK_UNGHR)
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
SETERRQ(PETSC_ERR_SUP,"GEHRD,ORGHR/UNGHR - Lapack routines are unavailable.");
|
SETERRQ(PETSC_ERR_SUP,"GEHRD,ORGHR/UNGHR - Lapack routines are unavailable.");
|
#else
|
#else
|
PetscScalar *tau,*work;
|
PetscScalar *tau,*work;
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
PetscInt i,j;
|
PetscInt i,j;
|
PetscBLASInt ilo,lwork,info,n,lda;
|
PetscBLASInt ilo,lwork,info,n,lda;
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
ierr = PetscLogEventBegin(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
ierr = PetscLogEventBegin(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
n = PetscBLASIntCast(n_);
|
n = PetscBLASIntCast(n_);
|
lda = PetscBLASIntCast(lda_);
|
lda = PetscBLASIntCast(lda_);
|
ierr = PetscMalloc(n*sizeof(PetscScalar),&tau);CHKERRQ(ierr);
|
ierr = PetscMalloc(n*sizeof(PetscScalar),&tau);CHKERRQ(ierr);
|
lwork = n;
|
lwork = n;
|
ierr = PetscMalloc(lwork*sizeof(PetscScalar),&work);CHKERRQ(ierr);
|
ierr = PetscMalloc(lwork*sizeof(PetscScalar),&work);CHKERRQ(ierr);
|
ilo = PetscBLASIntCast(k+1);
|
ilo = PetscBLASIntCast(k+1);
|
LAPACKgehrd_(&n,&ilo,&n,A,&lda,tau,work,&lwork,&info);
|
LAPACKgehrd_(&n,&ilo,&n,A,&lda,tau,work,&lwork,&info);
|
if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack xGEHRD %d",info);
|
if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack xGEHRD %d",info);
|
for (j=0;j<n-1;j++) {
|
for (j=0;j<n-1;j++) {
|
for (i=j+2;i<n;i++) {
|
for (i=j+2;i<n;i++) {
|
Q[i+j*n] = A[i+j*lda];
|
Q[i+j*n] = A[i+j*lda];
|
A[i+j*lda] = 0.0;
|
A[i+j*lda] = 0.0;
|
}
|
}
|
}
|
}
|
LAPACKorghr_(&n,&ilo,&n,Q,&n,tau,work,&lwork,&info);
|
LAPACKorghr_(&n,&ilo,&n,Q,&n,tau,work,&lwork,&info);
|
if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack xORGHR %d",info);
|
if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack xORGHR %d",info);
|
ierr = PetscFree(tau);CHKERRQ(ierr);
|
ierr = PetscFree(tau);CHKERRQ(ierr);
|
ierr = PetscFree(work);CHKERRQ(ierr);
|
ierr = PetscFree(work);CHKERRQ(ierr);
|
ierr = PetscLogEventEnd(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
ierr = PetscLogEventEnd(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
#endif
|
#endif
|
}
|
}
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "EPSDenseSchur"
|
#define __FUNCT__ "EPSDenseSchur"
|
/*@
|
/*@
|
EPSDenseSchur - Computes the upper (quasi-)triangular form of a dense
|
EPSDenseSchur - Computes the upper (quasi-)triangular form of a dense
|
upper Hessenberg matrix.
|
upper Hessenberg matrix.
|
|
|
Not Collective
|
Not Collective
|
|
|
Input Parameters:
|
Input Parameters:
|
+ n - dimension of the matrix
|
+ n - dimension of the matrix
|
. k - first active column
|
. k - first active column
|
- ldh - leading dimension of H
|
- ldh - leading dimension of H
|
|
|
Input/Output Parameters:
|
Input/Output Parameters:
|
+ H - on entry, the upper Hessenber matrix; on exit, the upper
|
+ H - on entry, the upper Hessenber matrix; on exit, the upper
|
(quasi-)triangular matrix (T)
|
(quasi-)triangular matrix (T)
|
- Z - on entry, initial transformation matrix; on exit, orthogonal
|
- Z - on entry, initial transformation matrix; on exit, orthogonal
|
matrix of Schur vectors
|
matrix of Schur vectors
|
|
|
Output Parameters:
|
Output Parameters:
|
+ wr - pointer to the array to store the computed eigenvalues
|
+ wr - pointer to the array to store the computed eigenvalues
|
- wi - imaginary part of the eigenvalues (only when using real numbers)
|
- wi - imaginary part of the eigenvalues (only when using real numbers)
|
|
|
Notes:
|
Notes:
|
This function computes the (real) Schur decomposition of an upper
|
This function computes the (real) Schur decomposition of an upper
|
Hessenberg matrix H: H*Z = Z*T, where T is an upper (quasi-)triangular
|
Hessenberg matrix H: H*Z = Z*T, where T is an upper (quasi-)triangular
|
matrix (returned in H), and Z is the orthogonal matrix of Schur vectors.
|
matrix (returned in H), and Z is the orthogonal matrix of Schur vectors.
|
Eigenvalues are extracted from the diagonal blocks of T and returned in
|
Eigenvalues are extracted from the diagonal blocks of T and returned in
|
wr,wi. Transformations are accumulated in Z so that on entry it can
|
wr,wi. Transformations are accumulated in Z so that on entry it can
|
contain the transformation matrix associated to the Hessenberg reduction.
|
contain the transformation matrix associated to the Hessenberg reduction.
|
|
|
Only active columns (from k to n) are computed.
|
Only active columns (from k to n) are computed.
|
|
|
Both H and Z are overwritten.
|
Both H and Z are overwritten.
|
|
|
This routine uses LAPACK routines xHSEQR.
|
This routine uses LAPACK routines xHSEQR.
|
|
|
Level: developer
|
Level: developer
|
|
|
.seealso: EPSDenseHessenberg(), EPSSortDenseSchur(), EPSDenseTridiagonal()
|
.seealso: EPSDenseHessenberg(), EPSSortDenseSchur(), EPSDenseTridiagonal()
|
@*/
|
@*/
|
PetscErrorCode EPSDenseSchur(PetscInt n_,PetscInt k,PetscScalar *H,PetscInt ldh_,PetscScalar *Z,PetscScalar *wr,PetscScalar *wi)
|
PetscErrorCode EPSDenseSchur(PetscInt n_,PetscInt k,PetscScalar *H,PetscInt ldh_,PetscScalar *Z,PetscScalar *wr,PetscScalar *wi)
|
{
|
{
|
#if defined(SLEPC_MISSING_LAPACK_HSEQR)
|
#if defined(SLEPC_MISSING_LAPACK_HSEQR)
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
SETERRQ(PETSC_ERR_SUP,"HSEQR - Lapack routine is unavailable.");
|
SETERRQ(PETSC_ERR_SUP,"HSEQR - Lapack routine is unavailable.");
|
#else
|
#else
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
PetscBLASInt ilo,lwork,info,n,ldh;
|
PetscBLASInt ilo,lwork,info,n,ldh;
|
PetscScalar *work;
|
PetscScalar *work;
|
#if !defined(PETSC_USE_COMPLEX)
|
#if !defined(PETSC_USE_COMPLEX)
|
PetscInt j;
|
PetscInt j;
|
#endif
|
#endif
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
ierr = PetscLogEventBegin(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
ierr = PetscLogEventBegin(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
n = PetscBLASIntCast(n_);
|
n = PetscBLASIntCast(n_);
|
ldh = PetscBLASIntCast(ldh_);
|
ldh = PetscBLASIntCast(ldh_);
|
lwork = n;
|
lwork = n;
|
ierr = PetscMalloc(lwork*sizeof(PetscScalar),&work);CHKERRQ(ierr);
|
ierr = PetscMalloc(lwork*sizeof(PetscScalar),&work);CHKERRQ(ierr);
|
ilo = PetscBLASIntCast(k+1);
|
ilo = PetscBLASIntCast(k+1);
|
#if !defined(PETSC_USE_COMPLEX)
|
#if !defined(PETSC_USE_COMPLEX)
|
LAPACKhseqr_("S","V",&n,&ilo,&n,H,&ldh,wr,wi,Z,&n,work,&lwork,&info);
|
LAPACKhseqr_("S","V",&n,&ilo,&n,H,&ldh,wr,wi,Z,&n,work,&lwork,&info);
|
for (j=0;j<k;j++) {
|
for (j=0;j<k;j++) {
|
if (j==n-1 || H[j*ldh+j+1] == 0.0) {
|
if (j==n-1 || H[j*ldh+j+1] == 0.0) {
|
/* real eigenvalue */
|
/* real eigenvalue */
|
wr[j] = H[j*ldh+j];
|
wr[j] = H[j*ldh+j];
|
wi[j] = 0.0;
|
wi[j] = 0.0;
|
} else {
|
} else {
|
/* complex eigenvalue */
|
/* complex eigenvalue */
|
wr[j] = H[j*ldh+j];
|
wr[j] = H[j*ldh+j];
|
wr[j+1] = H[j*ldh+j];
|
wr[j+1] = H[j*ldh+j];
|
wi[j] = sqrt(PetscAbsReal(H[j*ldh+j+1])) *
|
wi[j] = sqrt(PetscAbsReal(H[j*ldh+j+1])) *
|
sqrt(PetscAbsReal(H[(j+1)*ldh+j]));
|
sqrt(PetscAbsReal(H[(j+1)*ldh+j]));
|
wi[j+1] = -wi[j];
|
wi[j+1] = -wi[j];
|
j++;
|
j++;
|
}
|
}
|
}
|
}
|
#else
|
#else
|
LAPACKhseqr_("S","V",&n,&ilo,&n,H,&ldh,wr,Z,&n,work,&lwork,&info);
|
LAPACKhseqr_("S","V",&n,&ilo,&n,H,&ldh,wr,Z,&n,work,&lwork,&info);
|
#endif
|
#endif
|
if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack xHSEQR %d",info);
|
if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack xHSEQR %d",info);
|
|
|
ierr = PetscFree(work);CHKERRQ(ierr);
|
ierr = PetscFree(work);CHKERRQ(ierr);
|
ierr = PetscLogEventEnd(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
ierr = PetscLogEventEnd(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
#endif
|
#endif
|
}
|
}
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "EPSSortDenseSchur"
|
#define __FUNCT__ "EPSSortDenseSchur"
|
/*@
|
/*@
|
EPSSortDenseSchur - Reorders the Schur decomposition computed by
|
EPSSortDenseSchur - Reorders the Schur decomposition computed by
|
EPSDenseSchur().
|
EPSDenseSchur().
|
|
|
Not Collective
|
Not Collective
|
|
|
Input Parameters:
|
Input Parameters:
|
+ n - dimension of the matrix
|
+ n - dimension of the matrix
|
. k - first active column
|
. k - first active column
|
. ldt - leading dimension of T
|
. ldt - leading dimension of T
|
- which - eigenvalue sort order
|
- which - eigenvalue sort order
|
|
|
Input/Output Parameters:
|
Input/Output Parameters:
|
+ T - the upper (quasi-)triangular matrix
|
+ T - the upper (quasi-)triangular matrix
|
. Z - the orthogonal matrix of Schur vectors
|
. Z - the orthogonal matrix of Schur vectors
|
. wr - pointer to the array to store the computed eigenvalues
|
. wr - pointer to the array to store the computed eigenvalues
|
- wi - imaginary part of the eigenvalues (only when using real numbers)
|
- wi - imaginary part of the eigenvalues (only when using real numbers)
|
|
|
Notes:
|
Notes:
|
This function reorders the eigenvalues in wr,wi located in positions k
|
This function reorders the eigenvalues in wr,wi located in positions k
|
to n according to the sort order specified in which. The Schur
|
to n according to the sort order specified in which. The Schur
|
decomposition Z*T*Z^T, is also reordered by means of rotations so that
|
decomposition Z*T*Z^T, is also reordered by means of rotations so that
|
eigenvalues in the diagonal blocks of T follow the same order.
|
eigenvalues in the diagonal blocks of T follow the same order.
|
|
|
Both T and Z are overwritten.
|
Both T and Z are overwritten.
|
|
|
This routine uses LAPACK routines xTREXC.
|
This routine uses LAPACK routines xTREXC.
|
|
|
Level: developer
|
Level: developer
|
|
|
.seealso: EPSDenseHessenberg(), EPSDenseSchur(), EPSDenseTridiagonal()
|
.seealso: EPSDenseHessenberg(), EPSDenseSchur(), EPSDenseTridiagonal()
|
@*/
|
@*/
|
PetscErrorCode EPSSortDenseSchur(PetscInt n_,PetscInt k,PetscScalar *T,PetscInt ldt_,PetscScalar *Z,PetscScalar *wr,PetscScalar *wi,EPSWhich which)
|
PetscErrorCode EPSSortDenseSchur(PetscInt n_,PetscInt k,PetscScalar *T,PetscInt ldt_,PetscScalar *Z,PetscScalar *wr,PetscScalar *wi,EPSWhich which)
|
{
|
{
|
#if defined(SLEPC_MISSING_LAPACK_TREXC)
|
#if defined(SLEPC_MISSING_LAPACK_TREXC)
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
SETERRQ(PETSC_ERR_SUP,"TREXC - Lapack routine is unavailable.");
|
SETERRQ(PETSC_ERR_SUP,"TREXC - Lapack routine is unavailable.");
|
#else
|
#else
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
PetscReal value,v;
|
PetscReal value,v;
|
PetscInt i,j;
|
PetscInt i,j;
|
PetscBLASInt ifst,ilst,info,pos,n,ldt;
|
PetscBLASInt ifst,ilst,info,pos,n,ldt;
|
#if !defined(PETSC_USE_COMPLEX)
|
#if !defined(PETSC_USE_COMPLEX)
|
PetscScalar *work;
|
PetscScalar *work;
|
#endif
|
#endif
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
ierr = PetscLogEventBegin(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
ierr = PetscLogEventBegin(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
n = PetscBLASIntCast(n_);
|
n = PetscBLASIntCast(n_);
|
ldt = PetscBLASIntCast(ldt_);
|
ldt = PetscBLASIntCast(ldt_);
|
#if !defined(PETSC_USE_COMPLEX)
|
#if !defined(PETSC_USE_COMPLEX)
|
ierr = PetscMalloc(n*sizeof(PetscScalar),&work);CHKERRQ(ierr);
|
ierr = PetscMalloc(n*sizeof(PetscScalar),&work);CHKERRQ(ierr);
|
#endif
|
#endif
|
|
|
for (i=k;i<n-1;i++) {
|
for (i=k;i<n-1;i++) {
|
switch(which) {
|
switch(which) {
|
case EPS_LARGEST_MAGNITUDE:
|
case EPS_LARGEST_MAGNITUDE:
|
case EPS_SMALLEST_MAGNITUDE:
|
case EPS_SMALLEST_MAGNITUDE:
|
value = SlepcAbsEigenvalue(wr[i],wi[i]);
|
value = SlepcAbsEigenvalue(wr[i],wi[i]);
|
break;
|
break;
|
case EPS_LARGEST_REAL:
|
case EPS_LARGEST_REAL:
|
case EPS_SMALLEST_REAL:
|
case EPS_SMALLEST_REAL:
|
value = PetscRealPart(wr[i]);
|
value = PetscRealPart(wr[i]);
|
break;
|
break;
|
case EPS_LARGEST_IMAGINARY:
|
case EPS_LARGEST_IMAGINARY:
|
case EPS_SMALLEST_IMAGINARY:
|
case EPS_SMALLEST_IMAGINARY:
|
#if !defined(PETSC_USE_COMPLEX)
|
#if !defined(PETSC_USE_COMPLEX)
|
value = PetscAbsReal(wi[i]);
|
value = PetscAbsReal(wi[i]);
|
#else
|
#else
|
value = PetscImaginaryPart(wr[i]);
|
value = PetscImaginaryPart(wr[i]);
|
#endif
|
#endif
|
break;
|
break;
|
default: SETERRQ(1,"Wrong value of which");
|
default: SETERRQ(1,"Wrong value of which");
|
}
|
}
|
pos = 0;
|
pos = 0;
|
for (j=i+1;j<n;j++) {
|
for (j=i+1;j<n;j++) {
|
switch(which) {
|
switch(which) {
|
case EPS_LARGEST_MAGNITUDE:
|
case EPS_LARGEST_MAGNITUDE:
|
case EPS_SMALLEST_MAGNITUDE:
|
case EPS_SMALLEST_MAGNITUDE:
|
v = SlepcAbsEigenvalue(wr[j],wi[j]);
|
v = SlepcAbsEigenvalue(wr[j],wi[j]);
|
break;
|
break;
|
case EPS_LARGEST_REAL:
|
case EPS_LARGEST_REAL:
|
case EPS_SMALLEST_REAL:
|
case EPS_SMALLEST_REAL:
|
v = PetscRealPart(wr[j]);
|
v = PetscRealPart(wr[j]);
|
break;
|
break;
|
case EPS_LARGEST_IMAGINARY:
|
case EPS_LARGEST_IMAGINARY:
|
case EPS_SMALLEST_IMAGINARY:
|
case EPS_SMALLEST_IMAGINARY:
|
#if !defined(PETSC_USE_COMPLEX)
|
#if !defined(PETSC_USE_COMPLEX)
|
v = PetscAbsReal(wi[j]);
|
v = PetscAbsReal(wi[j]);
|
#else
|
#else
|
v = PetscImaginaryPart(wr[j]);
|
v = PetscImaginaryPart(wr[j]);
|
#endif
|
#endif
|
break;
|
break;
|
default: SETERRQ(1,"Wrong value of which");
|
default: SETERRQ(1,"Wrong value of which");
|
}
|
}
|
switch(which) {
|
switch(which) {
|
case EPS_LARGEST_MAGNITUDE:
|
case EPS_LARGEST_MAGNITUDE:
|
case EPS_LARGEST_REAL:
|
case EPS_LARGEST_REAL:
|
case EPS_LARGEST_IMAGINARY:
|
case EPS_LARGEST_IMAGINARY:
|
if (v > value) {
|
if (v > value) {
|
value = v;
|
value = v;
|
pos = j;
|
pos = j;
|
}
|
}
|
break;
|
break;
|
case EPS_SMALLEST_MAGNITUDE:
|
case EPS_SMALLEST_MAGNITUDE:
|
case EPS_SMALLEST_REAL:
|
case EPS_SMALLEST_REAL:
|
case EPS_SMALLEST_IMAGINARY:
|
case EPS_SMALLEST_IMAGINARY:
|
if (v < value) {
|
if (v < value) {
|
value = v;
|
value = v;
|
pos = j;
|
pos = j;
|
}
|
}
|
break;
|
break;
|
default: SETERRQ(1,"Wrong value of which");
|
default: SETERRQ(1,"Wrong value of which");
|
}
|
}
|
#if !defined(PETSC_USE_COMPLEX)
|
#if !defined(PETSC_USE_COMPLEX)
|
if (wi[j] != 0) j++;
|
if (wi[j] != 0) j++;
|
#endif
|
#endif
|
}
|
}
|
if (pos) {
|
if (pos) {
|
ifst = PetscBLASIntCast(pos + 1);
|
ifst = PetscBLASIntCast(pos + 1);
|
ilst = PetscBLASIntCast(i + 1);
|
ilst = PetscBLASIntCast(i + 1);
|
#if !defined(PETSC_USE_COMPLEX)
|
#if !defined(PETSC_USE_COMPLEX)
|
LAPACKtrexc_("V",&n,T,&ldt,Z,&n,&ifst,&ilst,work,&info);
|
LAPACKtrexc_("V",&n,T,&ldt,Z,&n,&ifst,&ilst,work,&info);
|
#else
|
#else
|
LAPACKtrexc_("V",&n,T,&ldt,Z,&n,&ifst,&ilst,&info);
|
LAPACKtrexc_("V",&n,T,&ldt,Z,&n,&ifst,&ilst,&info);
|
#endif
|
#endif
|
if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack xTREXC %d",info);
|
if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack xTREXC %d",info);
|
|
|
for (j=k;j<n;j++) {
|
for (j=k;j<n;j++) {
|
#if !defined(PETSC_USE_COMPLEX)
|
#if !defined(PETSC_USE_COMPLEX)
|
if (j==n-1 || T[j*ldt+j+1] == 0.0) {
|
if (j==n-1 || T[j*ldt+j+1] == 0.0) {
|
/* real eigenvalue */
|
/* real eigenvalue */
|
wr[j] = T[j*ldt+j];
|
wr[j] = T[j*ldt+j];
|
wi[j] = 0.0;
|
wi[j] = 0.0;
|
} else {
|
} else {
|
/* complex eigenvalue */
|
/* complex eigenvalue */
|
wr[j] = T[j*ldt+j];
|
wr[j] = T[j*ldt+j];
|
wr[j+1] = T[j*ldt+j];
|
wr[j+1] = T[j*ldt+j];
|
wi[j] = sqrt(PetscAbsReal(T[j*ldt+j+1])) *
|
wi[j] = sqrt(PetscAbsReal(T[j*ldt+j+1])) *
|
sqrt(PetscAbsReal(T[(j+1)*ldt+j]));
|
sqrt(PetscAbsReal(T[(j+1)*ldt+j]));
|
wi[j+1] = -wi[j];
|
wi[j+1] = -wi[j];
|
j++;
|
j++;
|
}
|
}
|
#else
|
#else
|
wr[j] = T[j*(ldt+1)];
|
wr[j] = T[j*(ldt+1)];
|
#endif
|
#endif
|
}
|
}
|
}
|
}
|
#if !defined(PETSC_USE_COMPLEX)
|
#if !defined(PETSC_USE_COMPLEX)
|
if (wi[i] != 0) i++;
|
if (wi[i] != 0) i++;
|
#endif
|
#endif
|
}
|
}
|
|
|
#if !defined(PETSC_USE_COMPLEX)
|
#if !defined(PETSC_USE_COMPLEX)
|
ierr = PetscFree(work);CHKERRQ(ierr);
|
ierr = PetscFree(work);CHKERRQ(ierr);
|
#endif
|
#endif
|
ierr = PetscLogEventEnd(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
ierr = PetscLogEventEnd(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
|
|
#endif
|
#endif
|
}
|
}
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "EPSSortDenseSchurTarget"
|
#define __FUNCT__ "EPSSortDenseSchurTarget"
|
/*@
|
/*@
|
EPSSortDenseSchurTarget - Reorders the Schur decomposition computed by
|
EPSSortDenseSchurTarget - Reorders the Schur decomposition computed by
|
EPSDenseSchur().
|
EPSDenseSchur().
|
|
|
Not Collective
|
Not Collective
|
|
|
Input Parameters:
|
Input Parameters:
|
+ n - dimension of the matrix
|
+ n - dimension of the matrix
|
. k - first active column
|
. k - first active column
|
. ldt - leading dimension of T
|
. ldt - leading dimension of T
|
. target - the target value
|
. target - the target value
|
- which - eigenvalue sort order
|
- which - eigenvalue sort order
|
|
|
Input/Output Parameters:
|
Input/Output Parameters:
|
+ T - the upper (quasi-)triangular matrix
|
+ T - the upper (quasi-)triangular matrix
|
. Z - the orthogonal matrix of Schur vectors
|
. Z - the orthogonal matrix of Schur vectors
|
. wr - pointer to the array to store the computed eigenvalues
|
. wr - pointer to the array to store the computed eigenvalues
|
- wi - imaginary part of the eigenvalues (only when using real numbers)
|
- wi - imaginary part of the eigenvalues (only when using real numbers)
|
|
|
Notes:
|
Notes:
|
This function reorders the eigenvalues in wr,wi located in positions k
|
This function reorders the eigenvalues in wr,wi located in positions k
|
to n according to increasing distance to the target. The parameter which
|
to n according to increasing distance to the target. The parameter which
|
is used to determine if distance is relative to magnitude, real axis,
|
is used to determine if distance is relative to magnitude, real axis,
|
or imaginary axis. The Schur decomposition Z*T*Z^T, is also reordered
|
or imaginary axis. The Schur decomposition Z*T*Z^T, is also reordered
|
by means of rotations so that eigenvalues in the diagonal blocks of T
|
by means of rotations so that eigenvalues in the diagonal blocks of T
|
follow the same order.
|
follow the same order.
|
|
|
Both T and Z are overwritten.
|
Both T and Z are overwritten.
|
|
|
This routine uses LAPACK routines xTREXC.
|
This routine uses LAPACK routines xTREXC.
|
|
|
Level: developer
|
Level: developer
|
|
|
.seealso: EPSDenseHessenberg(), EPSDenseSchur(), EPSDenseTridiagonal()
|
.seealso: EPSDenseHessenberg(), EPSDenseSchur(), EPSDenseTridiagonal()
|
@*/
|
@*/
|
PetscErrorCode EPSSortDenseSchurTarget(PetscInt n_,PetscInt k,PetscScalar *T,PetscInt ldt_,PetscScalar *Z,PetscScalar *wr,PetscScalar *wi,PetscScalar target,EPSWhich which)
|
PetscErrorCode EPSSortDenseSchurTarget(PetscInt n_,PetscInt k,PetscScalar *T,PetscInt ldt_,PetscScalar *Z,PetscScalar *wr,PetscScalar *wi,PetscScalar target,EPSWhich which)
|
{
|
{
|
#if defined(SLEPC_MISSING_LAPACK_TREXC)
|
#if defined(SLEPC_MISSING_LAPACK_TREXC)
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
SETERRQ(PETSC_ERR_SUP,"TREXC - Lapack routine is unavailable.");
|
SETERRQ(PETSC_ERR_SUP,"TREXC - Lapack routine is unavailable.");
|
#else
|
#else
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
PetscReal value,v;
|
PetscReal value,v;
|
PetscInt i,j;
|
PetscInt i,j;
|
PetscBLASInt ifst,ilst,info,pos,n,ldt;
|
PetscBLASInt ifst,ilst,info,pos,n,ldt;
|
#if !defined(PETSC_USE_COMPLEX)
|
#if !defined(PETSC_USE_COMPLEX)
|
PetscScalar *work;
|
PetscScalar *work;
|
#endif
|
#endif
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
ierr = PetscLogEventBegin(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
ierr = PetscLogEventBegin(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
n = PetscBLASIntCast(n_);
|
n = PetscBLASIntCast(n_);
|
ldt = PetscBLASIntCast(ldt_);
|
ldt = PetscBLASIntCast(ldt_);
|
#if !defined(PETSC_USE_COMPLEX)
|
#if !defined(PETSC_USE_COMPLEX)
|
ierr = PetscMalloc(n*sizeof(PetscScalar),&work);CHKERRQ(ierr);
|
ierr = PetscMalloc(n*sizeof(PetscScalar),&work);CHKERRQ(ierr);
|
#endif
|
#endif
|
|
|
for (i=k;i<n-1;i++) {
|
for (i=k;i<n-1;i++) {
|
switch(which) {
|
switch(which) {
|
case EPS_LARGEST_MAGNITUDE:
|
case EPS_LARGEST_MAGNITUDE:
|
/* complex target only allowed if scalartype=complex */
|
/* complex target only allowed if scalartype=complex */
|
value = SlepcAbsEigenvalue(wr[i]-target,wi[i]);
|
value = SlepcAbsEigenvalue(wr[i]-target,wi[i]);
|
break;
|
break;
|
case EPS_LARGEST_REAL:
|
case EPS_LARGEST_REAL:
|
value = PetscAbsReal(PetscRealPart(wr[i]-target));
|
value = PetscAbsReal(PetscRealPart(wr[i]-target));
|
break;
|
break;
|
case EPS_LARGEST_IMAGINARY:
|
case EPS_LARGEST_IMAGINARY:
|
#if !defined(PETSC_USE_COMPLEX)
|
#if !defined(PETSC_USE_COMPLEX)
|
/* complex target only allowed if scalartype=complex */
|
/* complex target only allowed if scalartype=complex */
|
value = PetscAbsReal(wi[i]);
|
value = PetscAbsReal(wi[i]);
|
#else
|
#else
|
value = PetscAbsReal(PetscImaginaryPart(wr[i]-target));
|
value = PetscAbsReal(PetscImaginaryPart(wr[i]-target));
|
#endif
|
#endif
|
break;
|
break;
|
default: SETERRQ(1,"Wrong value of which");
|
default: SETERRQ(1,"Wrong value of which");
|
}
|
}
|
pos = 0;
|
pos = 0;
|
for (j=i+1;j<n;j++) {
|
for (j=i+1;j<n;j++) {
|
switch(which) {
|
switch(which) {
|
case EPS_LARGEST_MAGNITUDE:
|
case EPS_LARGEST_MAGNITUDE:
|
/* complex target only allowed if scalartype=complex */
|
/* complex target only allowed if scalartype=complex */
|
v = SlepcAbsEigenvalue(wr[j]-target,wi[j]);
|
v = SlepcAbsEigenvalue(wr[j]-target,wi[j]);
|
break;
|
break;
|
case EPS_LARGEST_REAL:
|
case EPS_LARGEST_REAL:
|
v = PetscAbsReal(PetscRealPart(wr[j]-target));
|
v = PetscAbsReal(PetscRealPart(wr[j]-target));
|
break;
|
break;
|
case EPS_LARGEST_IMAGINARY:
|
case EPS_LARGEST_IMAGINARY:
|
#if !defined(PETSC_USE_COMPLEX)
|
#if !defined(PETSC_USE_COMPLEX)
|
/* complex target only allowed if scalartype=complex */
|
/* complex target only allowed if scalartype=complex */
|
v = PetscAbsReal(wi[j]);
|
v = PetscAbsReal(wi[j]);
|
#else
|
#else
|
v = PetscAbsReal(PetscImaginaryPart(wr[j]-target));
|
v = PetscAbsReal(PetscImaginaryPart(wr[j]-target));
|
#endif
|
#endif
|
break;
|
break;
|
default: SETERRQ(1,"Wrong value of which");
|
default: SETERRQ(1,"Wrong value of which");
|
}
|
}
|
if (v < value) {
|
if (v < value) {
|
value = v;
|
value = v;
|
pos = j;
|
pos = j;
|
}
|
}
|
#if !defined(PETSC_USE_COMPLEX)
|
#if !defined(PETSC_USE_COMPLEX)
|
if (wi[j] != 0) j++;
|
if (wi[j] != 0) j++;
|
#endif
|
#endif
|
}
|
}
|
if (pos) {
|
if (pos) {
|
ifst = PetscBLASIntCast(pos + 1);
|
ifst = PetscBLASIntCast(pos + 1);
|
ilst = PetscBLASIntCast(i + 1);
|
ilst = PetscBLASIntCast(i + 1);
|
#if !defined(PETSC_USE_COMPLEX)
|
#if !defined(PETSC_USE_COMPLEX)
|
LAPACKtrexc_("V",&n,T,&ldt,Z,&n,&ifst,&ilst,work,&info);
|
LAPACKtrexc_("V",&n,T,&ldt,Z,&n,&ifst,&ilst,work,&info);
|
#else
|
#else
|
LAPACKtrexc_("V",&n,T,&ldt,Z,&n,&ifst,&ilst,&info);
|
LAPACKtrexc_("V",&n,T,&ldt,Z,&n,&ifst,&ilst,&info);
|
#endif
|
#endif
|
if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack xTREXC %d",info);
|
if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack xTREXC %d",info);
|
|
|
for (j=k;j<n;j++) {
|
for (j=k;j<n;j++) {
|
#if !defined(PETSC_USE_COMPLEX)
|
#if !defined(PETSC_USE_COMPLEX)
|
if (j==n-1 || T[j*ldt+j+1] == 0.0) {
|
if (j==n-1 || T[j*ldt+j+1] == 0.0) {
|
/* real eigenvalue */
|
/* real eigenvalue */
|
wr[j] = T[j*ldt+j];
|
wr[j] = T[j*ldt+j];
|
wi[j] = 0.0;
|
wi[j] = 0.0;
|
} else {
|
} else {
|
/* complex eigenvalue */
|
/* complex eigenvalue */
|
wr[j] = T[j*ldt+j];
|
wr[j] = T[j*ldt+j];
|
wr[j+1] = T[j*ldt+j];
|
wr[j+1] = T[j*ldt+j];
|
wi[j] = sqrt(PetscAbsReal(T[j*ldt+j+1])) *
|
wi[j] = sqrt(PetscAbsReal(T[j*ldt+j+1])) *
|
sqrt(PetscAbsReal(T[(j+1)*ldt+j]));
|
sqrt(PetscAbsReal(T[(j+1)*ldt+j]));
|
wi[j+1] = -wi[j];
|
wi[j+1] = -wi[j];
|
j++;
|
j++;
|
}
|
}
|
#else
|
#else
|
wr[j] = T[j*(ldt+1)];
|
wr[j] = T[j*(ldt+1)];
|
#endif
|
#endif
|
}
|
}
|
}
|
}
|
#if !defined(PETSC_USE_COMPLEX)
|
#if !defined(PETSC_USE_COMPLEX)
|
if (wi[i] != 0) i++;
|
if (wi[i] != 0) i++;
|
#endif
|
#endif
|
}
|
}
|
|
|
#if !defined(PETSC_USE_COMPLEX)
|
#if !defined(PETSC_USE_COMPLEX)
|
ierr = PetscFree(work);CHKERRQ(ierr);
|
ierr = PetscFree(work);CHKERRQ(ierr);
|
#endif
|
#endif
|
ierr = PetscLogEventEnd(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
ierr = PetscLogEventEnd(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
|
|
#endif
|
#endif
|
}
|
}
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "EPSDenseTridiagonal"
|
#define __FUNCT__ "EPSDenseTridiagonal"
|
/*@
|
/*@
|
EPSDenseTridiagonal - Solves a real tridiagonal Hermitian Eigenvalue Problem.
|
EPSDenseTridiagonal - Solves a real tridiagonal Hermitian Eigenvalue Problem.
|
|
|
Not Collective
|
Not Collective
|
|
|
Input Parameters:
|
Input Parameters:
|
+ n - dimension of the eigenproblem
|
+ n - dimension of the eigenproblem
|
. A - pointer to the array containing the matrix values
|
. A - pointer to the array containing the matrix values
|
- lda - leading dimension of A
|
- lda - leading dimension of A
|
|
|
Output Parameters:
|
Output Parameters:
|
+ w - pointer to the array to store the computed eigenvalues
|
+ w - pointer to the array to store the computed eigenvalues
|
- V - pointer to the array to store the eigenvectors
|
- V - pointer to the array to store the eigenvectors
|
|
|
Notes:
|
Notes:
|
If V is PETSC_NULL then the eigenvectors are not computed.
|
If V is PETSC_NULL then the eigenvectors are not computed.
|
|
|
This routine use LAPACK routines DSTEVR.
|
This routine use LAPACK routines DSTEVR.
|
|
|
Level: developer
|
Level: developer
|
|
|
.seealso: EPSDenseNHEP(), EPSDenseHEP(), EPSDenseGNHEP(), EPSDenseGHEP()
|
.seealso: EPSDenseNHEP(), EPSDenseHEP(), EPSDenseGNHEP(), EPSDenseGHEP()
|
@*/
|
@*/
|
PetscErrorCode EPSDenseTridiagonal(PetscInt n_,PetscReal *D,PetscReal *E,PetscReal *w,PetscScalar *V)
|
PetscErrorCode EPSDenseTridiagonal(PetscInt n_,PetscReal *D,PetscReal *E,PetscReal *w,PetscScalar *V)
|
{
|
{
|
#if defined(SLEPC_MISSING_LAPACK_STEVR)
|
#if defined(SLEPC_MISSING_LAPACK_STEVR)
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
SETERRQ(PETSC_ERR_SUP,"STEVR - Lapack routine is unavailable.");
|
SETERRQ(PETSC_ERR_SUP,"STEVR - Lapack routine is unavailable.");
|
#else
|
#else
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
PetscReal abstol = 0.0,vl,vu,*work;
|
PetscReal abstol = 0.0,vl,vu,*work;
|
PetscBLASInt il,iu,m,*isuppz,n,lwork,*iwork,liwork,info;
|
PetscBLASInt il,iu,m,*isuppz,n,lwork,*iwork,liwork,info;
|
const char *jobz;
|
const char *jobz;
|
#if defined(PETSC_USE_COMPLEX)
|
#if defined(PETSC_USE_COMPLEX)
|
PetscInt i,j;
|
PetscInt i,j;
|
PetscReal *VV;
|
PetscReal *VV;
|
#endif
|
#endif
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
ierr = PetscLogEventBegin(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
ierr = PetscLogEventBegin(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
n = PetscBLASIntCast(n_);
|
n = PetscBLASIntCast(n_);
|
lwork = PetscBLASIntCast(20*n_);
|
lwork = PetscBLASIntCast(20*n_);
|
liwork = PetscBLASIntCast(10*n_);
|
liwork = PetscBLASIntCast(10*n_);
|
if (V) {
|
if (V) {
|
jobz = "V";
|
jobz = "V";
|
#if defined(PETSC_USE_COMPLEX)
|
#if defined(PETSC_USE_COMPLEX)
|
ierr = PetscMalloc(n*n*sizeof(PetscReal),&VV);CHKERRQ(ierr);
|
ierr = PetscMalloc(n*n*sizeof(PetscReal),&VV);CHKERRQ(ierr);
|
#endif
|
#endif
|
} else jobz = "N";
|
} else jobz = "N";
|
ierr = PetscMalloc(2*n*sizeof(PetscBLASInt),&isuppz);CHKERRQ(ierr);
|
ierr = PetscMalloc(2*n*sizeof(PetscBLASInt),&isuppz);CHKERRQ(ierr);
|
ierr = PetscMalloc(lwork*sizeof(PetscReal),&work);CHKERRQ(ierr);
|
ierr = PetscMalloc(lwork*sizeof(PetscReal),&work);CHKERRQ(ierr);
|
ierr = PetscMalloc(liwork*sizeof(PetscBLASInt),&iwork);CHKERRQ(ierr);
|
ierr = PetscMalloc(liwork*sizeof(PetscBLASInt),&iwork);CHKERRQ(ierr);
|
#if defined(PETSC_USE_COMPLEX)
|
#if defined(PETSC_USE_COMPLEX)
|
LAPACKstevr_(jobz,"A",&n,D,E,&vl,&vu,&il,&iu,&abstol,&m,w,VV,&n,isuppz,work,&lwork,iwork,&liwork,&info);
|
LAPACKstevr_(jobz,"A",&n,D,E,&vl,&vu,&il,&iu,&abstol,&m,w,VV,&n,isuppz,work,&lwork,iwork,&liwork,&info);
|
#else
|
#else
|
LAPACKstevr_(jobz,"A",&n,D,E,&vl,&vu,&il,&iu,&abstol,&m,w,V,&n,isuppz,work,&lwork,iwork,&liwork,&info);
|
LAPACKstevr_(jobz,"A",&n,D,E,&vl,&vu,&il,&iu,&abstol,&m,w,V,&n,isuppz,work,&lwork,iwork,&liwork,&info);
|
#endif
|
#endif
|
if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack DSTEVR %d",info);
|
if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack DSTEVR %d",info);
|
#if defined(PETSC_USE_COMPLEX)
|
#if defined(PETSC_USE_COMPLEX)
|
if (V) {
|
if (V) {
|
for (i=0;i<n;i++)
|
for (i=0;i<n;i++)
|
for (j=0;j<n;j++)
|
for (j=0;j<n;j++)
|
V[i*n+j] = VV[i*n+j];
|
V[i*n+j] = VV[i*n+j];
|
ierr = PetscFree(VV);CHKERRQ(ierr);
|
ierr = PetscFree(VV);CHKERRQ(ierr);
|
}
|
}
|
#endif
|
#endif
|
ierr = PetscFree(isuppz);CHKERRQ(ierr);
|
ierr = PetscFree(isuppz);CHKERRQ(ierr);
|
ierr = PetscFree(work);CHKERRQ(ierr);
|
ierr = PetscFree(work);CHKERRQ(ierr);
|
ierr = PetscFree(iwork);CHKERRQ(ierr);
|
ierr = PetscFree(iwork);CHKERRQ(ierr);
|
ierr = PetscLogEventEnd(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
ierr = PetscLogEventEnd(EPS_Dense,0,0,0,0);CHKERRQ(ierr);
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
#endif
|
#endif
|
}
|
}
|
|
|
|
|