/*
This file contains some simple default routines for common operations.
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
SLEPc - Scalable Library for Eigenvalue Problem Computations
Copyright (c) 2002-2010, Universidad 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 <private/epsimpl.h> /*I "slepceps.h" I*/
#include <slepcblaslapack.h>
#undef __FUNCT__
#define __FUNCT__ "EPSReset_Default"
PetscErrorCode EPSReset_Default(EPS eps)
{
PetscErrorCode ierr;
PetscFunctionBegin;
PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
ierr = EPSDefaultFreeWork(eps);CHKERRQ(ierr);
ierr = EPSFreeSolution(eps);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "EPSBackTransform_Default"
PetscErrorCode EPSBackTransform_Default(EPS eps)
{
PetscErrorCode ierr;
PetscFunctionBegin;
PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
ierr = STBackTransform(eps->OP,eps->nconv,eps->eigr,eps->eigi);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "EPSComputeVectors_Default"
/*
EPSComputeVectors_Default - Compute eigenvectors from the vectors
provided by the eigensolver. This version just copies the vectors
and is intended for solvers such as power that provide the eigenvector.
*/
PetscErrorCode EPSComputeVectors_Default(EPS eps)
{
PetscFunctionBegin;
eps->evecsavailable = PETSC_TRUE;
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "EPSComputeVectors_Hermitian"
/*
EPSComputeVectors_Hermitian - Copies the Lanczos vectors as eigenvectors
using purification for generalized eigenproblems.
*/
PetscErrorCode EPSComputeVectors_Hermitian(EPS eps)
{
PetscErrorCode ierr;
PetscInt i;
PetscReal norm;
Vec w;
PetscFunctionBegin;
if (eps->isgeneralized) {
/* Purify eigenvectors */
ierr = VecDuplicate(eps->V[0],&w);CHKERRQ(ierr);
for (i=0;i<eps->nconv;i++) {
ierr = VecCopy(eps->V[i],w);CHKERRQ(ierr);
ierr = STApply(eps->OP,w,eps->V[i]);CHKERRQ(ierr);
ierr = IPNorm(eps->ip,eps->V[i],&norm);CHKERRQ(ierr);
ierr = VecScale(eps->V[i],1.0/norm);CHKERRQ(ierr);
}
ierr = VecDestroy(&w);CHKERRQ(ierr);
}
eps->evecsavailable = PETSC_TRUE;
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "EPSComputeVectors_Schur"
/*
EPSComputeVectors_Schur - Compute eigenvectors from the vectors
provided by the eigensolver. This version is intended for solvers
that provide Schur vectors. Given the partial Schur decomposition
OP*V=V*T, the following steps are performed:
1) compute eigenvectors of T: T*Z=Z*D
2) compute eigenvectors of OP: X=V*Z
If left eigenvectors are required then also do Z'*Tl=D*Z', Y=W*Z
*/
PetscErrorCode EPSComputeVectors_Schur(EPS eps)
{
#if defined(SLEPC_MISSING_LAPACK_TREVC)
SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_SUP,"TREVC - Lapack routine is unavailable.");
#else
PetscErrorCode ierr;
PetscInt i;
PetscBLASInt ncv,nconv,mout,info,one = 1;
PetscScalar *Z,*work,tmp;
#if defined(PETSC_USE_COMPLEX)
PetscReal *rwork;
#else
PetscReal normi;
#endif
PetscReal norm;
Vec w;
PetscFunctionBegin;
ncv = PetscBLASIntCast(eps->ncv);
nconv = PetscBLASIntCast(eps->nconv);
if (eps->ishermitian) {
ierr = EPSComputeVectors_Hermitian(eps);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
ierr = PetscMalloc(nconv*nconv*sizeof(PetscScalar),&Z);CHKERRQ(ierr);
ierr = PetscMalloc(3*nconv*sizeof(PetscScalar),&work);CHKERRQ(ierr);
#if defined(PETSC_USE_COMPLEX)
ierr = PetscMalloc(nconv*sizeof(PetscReal),&rwork);CHKERRQ(ierr);
#endif
/* right eigenvectors */
#if !defined(PETSC_USE_COMPLEX)
LAPACKtrevc_("R","A",PETSC_NULL,&nconv,eps->T,&ncv,PETSC_NULL,&nconv,Z,&nconv,&nconv,&mout,work,&info);
#else
LAPACKtrevc_("R","A",PETSC_NULL,&nconv,eps->T,&ncv,PETSC_NULL,&nconv,Z,&nconv,&nconv,&mout,work,rwork,&info);
#endif
if (info) SETERRQ1(((PetscObject)eps)->comm,PETSC_ERR_LIB,"Error in Lapack xTREVC %i",info);
/* normalize eigenvectors (when not using purification nor balancing)*/
if (!(eps->ispositive || (eps->balance!=EPS_BALANCE_NONE && eps->D))) {
for (i=0;i<eps->nconv;i++) {
#if !defined(PETSC_USE_COMPLEX)
if (eps->eigi[i] != 0.0) {
norm = BLASnrm2_(&nconv,Z+i*nconv,&one);
normi = BLASnrm2_(&nconv,Z+(i+1)*nconv,&one);
tmp = 1.0 / SlepcAbsEigenvalue(norm,normi);
BLASscal_(&nconv,&tmp,Z+i*nconv,&one);
BLASscal_(&nconv,&tmp,Z+(i+1)*nconv,&one);
i++;
} else
#endif
{
norm = BLASnrm2_(&nconv,Z+i*nconv,&one);
tmp = 1.0 / norm;
BLASscal_(&nconv,&tmp,Z+i*nconv,&one);
}
}
}
/* AV = V * Z */
ierr = SlepcUpdateVectors(eps->nconv,eps->V,0,eps->nconv,Z,eps->nconv,PETSC_FALSE);CHKERRQ(ierr);
/* Purify eigenvectors */
if (eps->ispositive) {
ierr = VecDuplicate(eps->V[0],&w);CHKERRQ(ierr);
for (i=0;i<eps->nconv;i++) {
ierr = VecCopy(eps->V[i],w);CHKERRQ(ierr);
ierr = STApply(eps->OP,w,eps->V[i]);CHKERRQ(ierr);
}
ierr = VecDestroy(&w);CHKERRQ(ierr);
}
/* Fix eigenvectors if balancing was used */
if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
for (i=0;i<eps->nconv;i++) {
ierr = VecPointwiseDivide(eps->V[i],eps->V[i],eps->D);CHKERRQ(ierr);
}
}
/* normalize eigenvectors (when using purification or balancing) */
if (eps->ispositive || (eps->balance!=EPS_BALANCE_NONE && eps->D)) {
for (i=0;i<eps->nconv;i++) {
#if !defined(PETSC_USE_COMPLEX)
if (eps->eigi[i] != 0.0) {
ierr = VecNorm(eps->V[i],NORM_2,&norm);CHKERRQ(ierr);
ierr = VecNorm(eps->V[i+1],NORM_2,&normi);CHKERRQ(ierr);
tmp = 1.0 / SlepcAbsEigenvalue(norm,normi);
ierr = VecScale(eps->V[i],tmp);CHKERRQ(ierr);
ierr = VecScale(eps->V[i+1],tmp);CHKERRQ(ierr);
i++;
} else
#endif
{
ierr = VecNormalize(eps->V[i],PETSC_NULL);CHKERRQ(ierr);
}
}
}
/* left eigenvectors */
if (eps->leftvecs) {
#if !defined(PETSC_USE_COMPLEX)
LAPACKtrevc_("R","A",PETSC_NULL,&nconv,eps->Tl,&ncv,PETSC_NULL,&nconv,Z,&nconv,&nconv,&mout,work,&info);
#else
LAPACKtrevc_("R","A",PETSC_NULL,&nconv,eps->Tl,&ncv,PETSC_NULL,&nconv,Z,&nconv,&nconv,&mout,work,rwork,&info);
#endif
if (info) SETERRQ1(((PetscObject)eps)->comm,PETSC_ERR_LIB,"Error in Lapack xTREVC %i",info);
/* AW = W * Z */
ierr = SlepcUpdateVectors(eps->nconv,eps->W,0,eps->nconv,Z,eps->nconv,PETSC_FALSE);CHKERRQ(ierr);
}
ierr = PetscFree(Z);CHKERRQ(ierr);
ierr = PetscFree(work);CHKERRQ(ierr);
#if defined(PETSC_USE_COMPLEX)
ierr = PetscFree(rwork);CHKERRQ(ierr);
#endif
eps->evecsavailable = PETSC_TRUE;
PetscFunctionReturn(0);
#endif
}
#undef __FUNCT__
#define __FUNCT__ "EPSDefaultGetWork"
/*
EPSDefaultGetWork - Gets a number of work vectors.
*/
PetscErrorCode EPSDefaultGetWork(EPS eps,PetscInt nw)
{
PetscErrorCode ierr;
PetscFunctionBegin;
if (eps->nwork != nw) {
ierr = SlepcVecDestroyVecs(eps->nwork,&eps->work);CHKERRQ(ierr);
eps->nwork = nw;
ierr = SlepcVecDuplicateVecs(eps->V[0],nw,&eps->work);CHKERRQ(ierr);
ierr = PetscLogObjectParents(eps,nw,eps->work);
}
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "EPSDefaultFreeWork"
/*
EPSDefaultFreeWork - Free work vectors.
*/
PetscErrorCode EPSDefaultFreeWork(EPS eps)
{
PetscErrorCode ierr;
PetscFunctionBegin;
PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
ierr = SlepcVecDestroyVecs(eps->nwork,&eps->work);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "EPSEigRelativeConverged"
/*
EPSDefaultConverged - Checks convergence relative to the eigenvalue.
*/
PetscErrorCode EPSEigRelativeConverged(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
{
PetscReal w;
PetscFunctionBegin;
w = SlepcAbsEigenvalue(eigr,eigi);
*errest = res;
if (w > res) *errest = res / w;
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "EPSAbsoluteConverged"
/*
EPSAbsoluteConverged - Checks convergence absolutely.
*/
PetscErrorCode EPSAbsoluteConverged(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
{
PetscFunctionBegin;
*errest = res;
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "EPSNormRelativeConverged"
/*
EPSNormRelativeConverged - Checks convergence relative to the eigenvalue and
the matrix norms.
*/
PetscErrorCode EPSNormRelativeConverged(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
{
PetscReal w;
PetscFunctionBegin;
w = SlepcAbsEigenvalue(eigr,eigi);
*errest = res / (eps->nrma + w*eps->nrmb);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "EPSComputeTrueResidual"
/*
EPSComputeTrueResidual - Computes the true residual norm of a given Ritz pair:
||r|| = ||A*x - lambda*B*x||
where lambda is the Ritz value and x is the Ritz vector.
Real lambda:
lambda = eigr
x = V*Z (V is an array of nv vectors, Z has length nv)
Complex lambda:
lambda = eigr+i*eigi
x = V*Z[0*nv]+i*V*Z[1*nv] (Z has length 2*nv)
*/
PetscErrorCode EPSComputeTrueResidual(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscScalar *Z,Vec *V,PetscInt nv,PetscReal *resnorm)
{
PetscErrorCode ierr;
Vec x,y,z=0;
PetscReal norm;
PetscFunctionBegin;
/* allocate workspace */
ierr = VecDuplicate(V[0],&x);CHKERRQ(ierr);
ierr = VecDuplicate(V[0],&y);CHKERRQ(ierr);
if (!eps->ishermitian && eps->ispositive) { ierr = VecDuplicate(V[0],&z);CHKERRQ(ierr); }
/* compute eigenvector */
ierr = SlepcVecMAXPBY(x,0.0,1.0,nv,Z,V);CHKERRQ(ierr);
/* purify eigenvector in positive generalized problems */
if (eps->ispositive) {
ierr = STApply(eps->OP,x,y);CHKERRQ(ierr);
if (eps->ishermitian) {
ierr = IPNorm(eps->ip,y,&norm);CHKERRQ(ierr);
} else {
ierr = VecNorm(y,NORM_2,&norm);CHKERRQ(ierr);
}
ierr = VecScale(y,1.0/norm);CHKERRQ(ierr);
ierr = VecCopy(y,x);CHKERRQ(ierr);
}
/* fix eigenvector if balancing is used */
if (!eps->ishermitian && eps->balance!=EPS_BALANCE_NONE && eps->D) {
ierr = VecPointwiseDivide(x,x,eps->D);CHKERRQ(ierr);
ierr = VecNormalize(x,&norm);CHKERRQ(ierr);
}
#if !defined(PETSC_USE_COMPLEX)
/* compute imaginary part of eigenvector */
if (!eps->ishermitian && eigi != 0.0) {
ierr = SlepcVecMAXPBY(y,0.0,1.0,nv,Z+nv,V);CHKERRQ(ierr);
if (eps->ispositive) {
ierr = STApply(eps->OP,y,z);CHKERRQ(ierr);
ierr = VecNorm(z,NORM_2,&norm);CHKERRQ(ierr);
ierr = VecScale(z,1.0/norm);CHKERRQ(ierr);
ierr = VecCopy(z,y);CHKERRQ(ierr);
}
if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
ierr = VecPointwiseDivide(y,y,eps->D);CHKERRQ(ierr);
ierr = VecNormalize(y,&norm);CHKERRQ(ierr);
}
}
#endif
/* compute relative error and update convergence flag */
ierr = EPSComputeResidualNorm_Private(eps,eigr,eigi,x,y,resnorm);
/* free workspace */
ierr = VecDestroy(&x);CHKERRQ(ierr);
ierr = VecDestroy(&y);CHKERRQ(ierr);
ierr = VecDestroy(&z);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "EPSBuildBalance_Krylov"
/*
EPSBuildBalance_Krylov - uses a Krylov subspace method to compute the
diagonal matrix to be applied for balancing in non-Hermitian problems.
*/
PetscErrorCode EPSBuildBalance_Krylov(EPS eps)
{
Vec z,p,r;
PetscInt i,j;
PetscReal norma;
PetscScalar *pz,*pD;
const PetscScalar *pr,*pp;
PetscErrorCode ierr;
PetscFunctionBegin;
ierr = VecDuplicate(eps->V[0],&r);CHKERRQ(ierr);
ierr = VecDuplicate(eps->V[0],&p);CHKERRQ(ierr);
ierr = VecDuplicate(eps->V[0],&z);CHKERRQ(ierr);
ierr = VecSet(eps->D,1.0);CHKERRQ(ierr);
for (j=0;j<eps->balance_its;j++) {
/* Build a random vector of +-1's */
ierr = SlepcVecSetRandom(z,eps->rand);CHKERRQ(ierr);
ierr = VecGetArray(z,&pz);CHKERRQ(ierr);
for (i=0;i<eps->nloc;i++) {
if (PetscRealPart(pz[i])<0.5) pz[i]=-1.0;
else pz[i]=1.0;
}
ierr = VecRestoreArray(z,&pz);CHKERRQ(ierr);
/* Compute p=DA(D\z) */
ierr = VecPointwiseDivide(r,z,eps->D);CHKERRQ(ierr);
ierr = STApply(eps->OP,r,p);CHKERRQ(ierr);
ierr = VecPointwiseMult(p,p,eps->D);CHKERRQ(ierr);
if (j==0) {
/* Estimate the matrix inf-norm */
ierr = VecAbs(p);CHKERRQ(ierr);
ierr = VecMax(p,PETSC_NULL,&norma);CHKERRQ(ierr);
}
if (eps->balance == EPS_BALANCE_TWOSIDE) {
/* Compute r=D\(A'Dz) */
ierr = VecPointwiseMult(z,z,eps->D);CHKERRQ(ierr);
ierr = STApplyTranspose(eps->OP,z,r);CHKERRQ(ierr);
ierr = VecPointwiseDivide(r,r,eps->D);CHKERRQ(ierr);
}
/* Adjust values of D */
ierr = VecGetArrayRead(r,&pr);CHKERRQ(ierr);
ierr = VecGetArrayRead(p,&pp);CHKERRQ(ierr);
ierr = VecGetArray(eps->D,&pD);CHKERRQ(ierr);
for (i=0;i<eps->nloc;i++) {
if (eps->balance == EPS_BALANCE_TWOSIDE) {
if (PetscAbsScalar(pp[i])>eps->balance_cutoff*norma && pr[i]!=0.0)
pD[i] *= sqrt(PetscAbsScalar(pr[i]/pp[i]));
} else {
if (pp[i]!=0.0) pD[i] *= 1.0/PetscAbsScalar(pp[i]);
}
}
ierr = VecRestoreArrayRead(r,&pr);CHKERRQ(ierr);
ierr = VecRestoreArrayRead(p,&pp);CHKERRQ(ierr);
ierr = VecRestoreArray(eps->D,&pD);CHKERRQ(ierr);
}
ierr = VecDestroy(&r);CHKERRQ(ierr);
ierr = VecDestroy(&p);CHKERRQ(ierr);
ierr = VecDestroy(&z);CHKERRQ(ierr);
PetscFunctionReturn(0);
}