Subversion Repositories slepc-dev

Compare Revisions

Ignore whitespace Rev 2828 → Rev 2829

/trunk/include/slepcblaslapack.h
110,6 → 110,7
#define LAPACKtgexc_ SLEPC_BLASLAPACK(tgexc,TGEXC)
#define LAPACKlag2_ SLEPC_BLASLAPACKREAL(lag2,LAG2)
#define LAPACKlasv2_ SLEPC_BLASLAPACKREAL(lasv2,LASV2)
#define LAPACKlartg_ SLEPC_BLASLAPACKREAL(lartg,LARTG)
#if !defined(PETSC_USE_COMPLEX)
#define LAPACKorghr_ SLEPC_BLASLAPACK(orghr,ORGHR)
#define LAPACKorgqr_ SLEPC_BLASLAPACK(orgqr,ORGQR)
232,6 → 233,7
extern void SLEPC_BLASLAPACK(gelqf,GELQF) (PetscBLASInt*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscScalar*,PetscBLASInt*,PetscBLASInt*);
extern void SLEPC_BLASLAPACKREAL(lag2,LAG2) (PetscReal*,PetscBLASInt*,PetscReal*,PetscBLASInt*,PetscReal*,PetscReal*,PetscReal*,PetscReal*,PetscReal*,PetscReal*);
extern void SLEPC_BLASLAPACKREAL(lasv2,LASV2) (PetscReal*,PetscReal*,PetscReal*,PetscReal*,PetscReal*,PetscReal*,PetscReal*,PetscReal*,PetscReal*);
extern void SLEPC_BLASLAPACKREAL(lartg,LARTG) (PetscReal*,PetscReal*,PetscReal*,PetscReal*,PetscReal*);
extern void BLASrot_(PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscReal*,PetscReal*);
#if !defined(PETSC_USE_COMPLEX)
extern void SLEPC_BLASLAPACK(tgexc,TGEXC) (PetscBLASInt*,PetscBLASInt*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscBLASInt*);
293,6 → 295,7
extern void PETSC_STDCALL SLEPC_BLASLAPACK(gelqf,GELQF) (PetscBLASInt*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscScalar*,PetscBLASInt*,PetscBLASInt*);
extern void PETSC_STDCALL SLEPC_BLASLAPACKREAL(lag2,LAG2) (PetscReal*,PetscBLASInt*,PetscReal*,PetscBLASInt*,PetscReal*,PetscReal*,PetscReal*,PetscReal*,PetscReal*,PetscReal*);
extern void PETSC_STDCALL SLEPC_BLASLAPACKREAL(lasv2,LASV2) (PetscReal*,PetscReal*,PetscReal*,PetscReal*,PetscReal*,PetscReal*,PetscReal*,PetscReal*,PetscReal*);
extern void PETSC_STDCALL SLEPC_BLASLAPACKREAL(lartg,LARTG) (PetscReal*,PetscReal*,PetscReal*,PetscReal*,PetscReal*);
extern void PETSC_STDCALL BLASrot_(PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscReal*,PetscReal*);
#if !defined(PETSC_USE_COMPLEX)
extern void PETSC_STDCALL SLEPC_BLASLAPACK(tgexc,TGEXC) (PetscBLASInt*,PetscBLASInt*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscBLASInt*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscBLASInt*);
/trunk/config/lapack.py
61,11 → 61,11
 
# LAPACK functions which are always used in real version
if petscconf.PRECISION == 'single':
functions += ['sstevr','sbdsdc','slamch','slag2','slasv2']
functions += ['sstevr','sbdsdc','slamch','slag2','slasv2','slartg']
elif petscconf.PRECISION == '__float128':
functions += ['qstevr','qbdsdc','qlamch','qlag2','qlasv2']
functions += ['qstevr','qbdsdc','qlamch','qlag2','qlasv2','qlartg']
else:
functions += ['dstevr','dbdsdc','dlamch','dlag2','dlasv2']
functions += ['dstevr','dbdsdc','dlamch','dlag2','dlasv2','dlartg']
# check for all functions at once
all = []
/trunk/src/ps/pshep.c
193,12 → 193,106
}
 
#undef __FUNCT__
#define __FUNCT__ "PSIntermediate_HEP"
#define __FUNCT__ "ArrowTridiag"
/*
Reduce to tridiagonal form
ARROWTRIDIAG reduces a symmetric arrowhead matrix of the form
 
[ d 0 0 0 e ]
[ 0 d 0 0 e ]
A = [ 0 0 d 0 e ]
[ 0 0 0 d e ]
[ e e e e d ]
 
to tridiagonal form
 
[ d e 0 0 0 ]
[ e d e 0 0 ]
T = Q'*A*Q = [ 0 e d e 0 ],
[ 0 0 e d e ]
[ 0 0 0 e d ]
 
where Q is an orthogonal matrix. Rutishauser's algorithm is used to
perform the reduction, which requires O(n**2) flops. The accumulation
of the orthogonal factor Q, however, requires O(n**3) flops.
 
Arguments
=========
 
N (input) INTEGER
The order of the matrix A. N >= 0.
 
D (input/output) DOUBLE PRECISION array, dimension (N)
On entry, the diagonal entries of the matrix A to be
reduced.
On exit, the diagonal entries of the reduced matrix T.
 
E (input/output) DOUBLE PRECISION array, dimension (N-1)
On entry, the off-diagonal entries of the matrix A to be
reduced.
On exit, the subdiagonal entries of the reduced matrix T.
 
Q (input/output) DOUBLE PRECISION array, dimension (LDQ, N)
On exit, the orthogonal matrix Q.
 
LDQ (input) INTEGER
The leading dimension of the array Q.
 
Note
====
Based on Fortran code contributed by Daniel Kressner
*/
static PetscErrorCode PSIntermediate_HEP(PS ps)
static PetscErrorCode ArrowTridiag(PetscBLASInt *n,PetscReal *d,PetscReal *e,PetscScalar *Q,PetscBLASInt *ldq)
{
PetscBLASInt i,j,j2,ld=*ldq,one=1;
PetscReal c,s,p,off,temp;
 
PetscFunctionBegin;
if (*n<=2) PetscFunctionReturn(0);
 
for (j=0;j<*n-2;j++) {
 
/* Eliminate entry e(j) by a rotation in the planes (j,j+1) */
temp = e[j+1];
LAPACKlartg_(&temp,&e[j],&c,&s,&e[j+1]);
s = -s;
 
/* Apply rotation to diagonal elements */
temp = d[j+1];
e[j] = c*s*(temp-d[j]);
d[j+1] = s*s*d[j] + c*c*temp;
d[j] = c*c*d[j] + s*s*temp;
 
/* Apply rotation to Q */
j2 = j+2;
BLASrot_(&j2,Q+j*ld,&one,Q+(j+1)*ld,&one,&c,&s);
 
/* Chase newly introduced off-diagonal entry to the top left corner */
for (i=j-1;i>=0;i--) {
off = -s*e[i];
e[i] = c*e[i];
temp = e[i+1];
LAPACKlartg_(&temp,&off,&c,&s,&e[i+1]);
s = -s;
temp = (d[i]-d[i+1])*s - 2.0*c*e[i];
p = s*temp;
d[i+1] += p;
d[i] -= p;
e[i] = -e[i] - c*temp;
j2 = j+2;
BLASrot_(&j2,Q+i*ld,&one,Q+(i+1)*ld,&one,&c,&s);
}
}
PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "PSIntermediate_HEP_Flip"
/*
Reduce to tridiagonal form by flipping the matrix and using standard
LAPACK tridiagonal reduction
*/
static PetscErrorCode PSIntermediate_HEP_Flip(PS ps)
{
#if defined(SLEPC_MISSING_LAPACK_SYTRD) || defined(SLEPC_MISSING_LAPACK_ORGTR)
PetscFunctionBegin;
SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"SYTRD/ORGTR - Lapack routine is unavailable.");
289,6 → 383,66
}
 
#undef __FUNCT__
#define __FUNCT__ "PSIntermediate_HEP"
/*
Reduce to tridiagonal form by means of ArrowTridiag.
*/
static PetscErrorCode PSIntermediate_HEP(PS ps)
{
#if defined(SLEPC_MISSING_LAPACK_SYTRD) || defined(SLEPC_MISSING_LAPACK_ORGTR)
PetscFunctionBegin;
SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"SYTRD/ORGTR - Lapack routine is unavailable.");
#else
PetscErrorCode ierr;
PetscInt i,j;
PetscBLASInt n1,n2,n3,lwork,info,l,n,ld,off;
PetscScalar *A,*S,*Q,*work,*tau;
PetscReal *d,*e;
 
PetscFunctionBegin;
n = PetscBLASIntCast(ps->n);
l = PetscBLASIntCast(ps->l);
ld = PetscBLASIntCast(ps->ld);
n1 = PetscBLASIntCast(ps->k-l+1); /* size of leading block, excluding locked */
n2 = PetscBLASIntCast(n-ps->k-1); /* size of trailing block */
n3 = n1+n2;
off = l+l*ld;
A = ps->mat[PS_MAT_A];
Q = ps->mat[PS_MAT_Q];
d = ps->rmat[PS_MAT_T];
e = ps->rmat[PS_MAT_T]+ld;
ierr = PetscMemzero(Q,ld*ld*sizeof(PetscScalar));CHKERRQ(ierr);
for (i=0;i<n;i++) Q[i+i*ld] = 1.0;
 
if (ps->compact) {
 
if (ps->state<PS_STATE_INTERMEDIATE) ArrowTridiag(&n1,d+l,e+l,Q+off,&ld);
 
} else {
 
for (i=0;i<l;i++) { d[i] = PetscRealPart(A[i+i*ld]); e[i] = 0.0; }
 
if (ps->state<PS_STATE_INTERMEDIATE) {
ierr = PSCopyMatrix_Private(ps,PS_MAT_Q,PS_MAT_A);CHKERRQ(ierr);
ierr = PSAllocateWork_Private(ps,ld+ld*ld,0,0);CHKERRQ(ierr);
tau = ps->work;
work = ps->work+ld;
lwork = ld*ld;
LAPACKsytrd_("L",&n3,Q+off,&ld,d+l,e+l,tau,work,&lwork,&info);
if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in Lapack xSYTRD %d",info);
LAPACKorgtr_("L",&n3,Q+off,&ld,tau,work,&lwork,&info);
if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in Lapack xORGTR %d",info);
} else {
/* copy tridiagonal to d,e */
for (i=l;i<n;i++) d[i] = PetscRealPart(A[i+i*ld]);
for (i=l;i<n-1;i++) e[i] = PetscRealPart(A[(i+1)+i*ld]);
}
}
PetscFunctionReturn(0);
#endif
}
 
#undef __FUNCT__
#define __FUNCT__ "PSSolve_HEP_QR"
PetscErrorCode PSSolve_HEP_QR(PS ps,PetscScalar *wr,PetscScalar *wi)
{