Subversion Repositories slepc-dev

Compare Revisions

Ignore whitespace Rev 947 → Rev 948

/trunk/include/slepcblaslapack.h
62,6 → 62,7
 
#endif
 
#define BLAStrsm_ SLEPC_BLASLAPACK(trsm,TRSM)
#define LAPACKlaev2_ SLEPC_BLASLAPACK(laev2,LAEV2)
#define LAPACKgehrd_ SLEPC_BLASLAPACK(gehrd,GEHRD)
#define LAPACKlanhs_ SLEPC_BLASLAPACK(lanhs,LANHS)
72,15 → 73,18
#define LAPACKtrevc_ SLEPC_BLASLAPACK(trevc,TREVC)
#define LAPACKgeevx_ SLEPC_BLASLAPACK(geevx,GEEVX)
#define LAPACKggevx_ SLEPC_BLASLAPACK(ggevx,GEEVX)
#define LAPACKgelqf_ SLEPC_BLASLAPACK(gelqf,GELQF)
 
#if !defined(PETSC_USE_COMPLEX)
#define LAPACKorghr_ SLEPC_BLASLAPACK(orghr,ORGHR)
#define LAPACKsyevr_ SLEPC_BLASLAPACK(syevr,SYEVR)
#define LAPACKsygvd_ SLEPC_BLASLAPACK(sygvd,SYGVD)
#define LAPACKormlq_ SLEPC_BLASLAPACK(ormlq,ORMLQ)
#else
#define LAPACKorghr_ SLEPC_BLASLAPACK(unghr,UNGHR)
#define LAPACKsyevr_ SLEPC_BLASLAPACK(heevr,HEEVR)
#define LAPACKsygvd_ SLEPC_BLASLAPACK(hegvd,HEGVD)
#define LAPACKormlq_ SLEPC_BLASLAPACK(unmlq,UNMLQ)
#endif
 
#define LAPACKlamch_ SLEPC_BLASLAPACKREAL(lamch,LAMCH)
96,6 → 100,9
EXTERN void LAPACKorghr_(PetscBLASInt*,PetscBLASInt*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscScalar*,PetscBLASInt*,PetscBLASInt*);
EXTERN void LAPACKgetri_(PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscBLASInt*);
EXTERN void LAPACKstevr_(const char*,const char*,PetscBLASInt*,PetscReal*,PetscReal*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt*,PetscReal*,PetscBLASInt*,PetscReal*,PetscReal*,PetscBLASInt*,PetscBLASInt*,PetscReal*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt,PetscBLASInt);
EXTERN void LAPACKgelqf_(PetscBLASInt*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscScalar*,PetscBLASInt*,PetscBLASInt*);
EXTERN void BLAStrsm_(const char*,const char*,const char*,const char*,PetscBLASInt*,PetscBLASInt*,PetscScalar*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscBLASInt,PetscBLASInt,PetscBLASInt,PetscBLASInt);
EXTERN void LAPACKormlq_(const char*,const char*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt,PetscBLASInt);
 
#if !defined(PETSC_USE_COMPLEX)
EXTERN void LAPACKhseqr_(const char*,const char*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscScalar*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt,PetscBLASInt);
/trunk/config/lapack.py
7,16 → 7,16
 
def Check(conf):
 
functions = ['laev2','gehrd','lanhs','lange','getri','hseqr','trexc','trevc','geevx','ggevx']
functions = ['laev2','gehrd','lanhs','lange','getri','hseqr','trexc','trevc','geevx','ggevx','gelqf']
 
if petscconf.SCALAR == 'real':
functions += ['orghr','syevr','sygvd']
functions += ['orghr','syevr','sygvd','ormlq']
if petscconf.PRECISION == 'double':
prefix = 'd'
else:
prefix = 's'
else:
functions += ['unghr','heevr','hegvd']
functions += ['unghr','heevr','hegvd','unmlq']
if petscconf.PRECISION == 'double':
prefix = 'z'
else:
60,4 → 60,4
conf.write(' -DSLEPC_MISSING_LAPACK_' + i.upper())
conf.write('\n')
return missing
return missing
/trunk/src/eps/interface/orthog.c
3,6 → 3,7
See the SLEPc users manual for a detailed explanation.
*/
#include "src/eps/epsimpl.h" /*I "slepceps.h" I*/
#include "slepcblaslapack.h"
 
#undef __FUNCT__
#define __FUNCT__ "EPSQRDecomposition"
298,6 → 299,10
#define __FUNCT__ "EPSCGSBiOrthogonalization"
static PetscErrorCode EPSCGSBiOrthogonalization(EPS eps,int n,Vec *V,Vec *W,Vec v,PetscScalar *H,PetscReal *hnorm,PetscReal *norm)
{
#if defined(SLEPC_MISSING_LAPACK_GELQF) || defined(SLEPC_MISSING_LAPACK_ORMLQ)
PetscFunctionBegin;
SETERRQ(PETSC_ERR_SUP,"xGELQF - Lapack routine is unavailable.");
#else
PetscErrorCode ierr;
int j,ione=1,lwork,info;
PetscScalar shh[100],*lhh,*vw,*tau,one=1.0,*work;
317,58 → 322,21
lwork = n;
ierr = PetscMalloc(n*sizeof(PetscScalar),&tau);CHKERRQ(ierr);
ierr = PetscMalloc(lwork*sizeof(PetscScalar),&work);CHKERRQ(ierr);
dgelqf_(&n,&n,vw,&n,tau,work,&lwork,&info);
if (info<0) SETERRQ(1,"DGELQF");
LAPACKgelqf_(&n,&n,vw,&n,tau,work,&lwork,&info);
if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack xGELQF %i",info);
/*** First orthogonalization ***/
 
/* h = W^* v */
/* q = v - V h */
ierr = STMInnerProduct(eps->OP,n,v,W,H);CHKERRQ(ierr);
dtrsm_("L","L","N","N",&n,&ione,&one,vw,&n,H,&n);
dormlq_("L","N",&n,&ione,&n,vw,&n,tau,H,&n,work,&lwork,&info);
if (info<0) SETERRQ(1,"DORMLQ");
BLAStrsm_("L","L","N","N",&n,&ione,&one,vw,&n,H,&n,1,1,1,1);
LAPACKormlq_("L","N",&n,&ione,&n,vw,&n,tau,H,&n,work,&lwork,&info,1,1);
if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack xORMLQ %i",info);
ierr = VecSet(w,0.0);CHKERRQ(ierr);
ierr = VecMAXPY(w,n,H,V);CHKERRQ(ierr);
ierr = VecAXPY(v,-1.0,w);CHKERRQ(ierr);
 
#if 0
/* compute hnorm */
if (hnorm) {
*hnorm = 0.0;
for (j=0; j<n; j++) {
*hnorm += PetscRealPart(H[j] * PetscConj(H[j]));
}
*hnorm = sqrt(*hnorm);
}
/* compute norm of v for refinement or linear dependence checking */
if (eps->orthog_ref == EPS_ORTH_REFINE_IFNEEDED ||
(eps->orthog_ref == EPS_ORTH_REFINE_ALWAYS && hnorm) ) {
ierr = STNorm(eps->OP,v,norm);CHKERRQ(ierr);
}
 
/*** Second orthogonalization if necessary ***/
/* if ||q|| < eta ||h|| */
if ((eps->orthog_ref == EPS_ORTH_REFINE_IFNEEDED && *norm < eps->orthog_eta * *hnorm) ||
eps->orthog_ref == EPS_ORTH_REFINE_ALWAYS) {
PetscLogInfo((eps,"EPSClassicalGramSchmidtOrthogonalization:Performing iterative refinement wnorm %g hnorm %g\n",norm ? *norm : 0,hnorm ? *hnorm : 0));
 
/* s = W^* q */
/* q = q - V s ; h = h + s */
ierr = STMInnerProduct(eps->OP,n,v,W,lhh);CHKERRQ(ierr);
for (j=0;j<n;j++) {
H[j] += lhh[j];
}
ierr = VecSet(w,0.0);CHKERRQ(ierr);
ierr = VecMAXPY(w,n,lhh,V);CHKERRQ(ierr);
ierr = VecAXPY(v,-1.0,w);CHKERRQ(ierr);
 
if (hnorm) *hnorm = *norm;
}
#endif
 
/* compute norm of v */
if (norm) { ierr = STNorm(eps->OP,v,norm);CHKERRQ(ierr); }
378,6 → 346,7
ierr = PetscFree(work);CHKERRQ(ierr);
ierr = VecDestroy(w);CHKERRQ(ierr);
PetscFunctionReturn(0);
#endif
}
 
#undef __FUNCT__