| Line 55... |
Line 55... |
{
|
{
|
PetscInt i;
|
PetscInt i;
|
PetscReal norm_;
|
PetscReal norm_;
|
/* s-normalization */
|
/* s-normalization */
|
norm_ = 0.0;
|
norm_ = 0.0;
|
for(i=0;i<n;i++){norm_ += X[i]*s[i]*X[i];}
|
for(i=0;i<n;i++){norm_ += PetscRealPart(X[i]*s[i]*X[i]);}
|
if(norm_<0){norm_ = -PetscSqrtReal(-norm_);}
|
if(norm_<0){norm_ = -PetscSqrtReal(-norm_);}
|
else {norm_ = PetscSqrtReal(norm_);}
|
else {norm_ = PetscSqrtReal(norm_);}
|
for(i=0;i<n;i++)X[i] /= norm_;
|
for(i=0;i<n;i++)X[i] /= norm_;
|
if(norm) *norm = norm_;
|
if(norm) *norm = norm_;
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
| Line 219... |
Line 219... |
PetscFunctionBegin;
|
PetscFunctionBegin;
|
if (left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented for left eigenvectors");
|
if (left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented for left eigenvectors");
|
else mat = PS_MAT_Q;
|
else mat = PS_MAT_Q;
|
k_ = *k;
|
k_ = *k;
|
if(k_ < ps->n-1){
|
if(k_ < ps->n-1){
|
b = (ps->compact)?*(ps->rmat[PS_MAT_T]+ld+k_):*(ps->mat[PS_MAT_A]+(k_+1)+ld*k_);
|
b = (ps->compact)?*(ps->rmat[PS_MAT_T]+ld+k_):PetscRealPart(*(ps->mat[PS_MAT_A]+(k_+1)+ld*k_));
|
}else b = 0;
|
}else b = 0.0;
|
if(b==0){/* real */
|
if(b == 0.0){/* real */
|
ierr = PetscMemcpy(ps->mat[mat]+(k_)*ld,Q+(k_)*ld,ld*sizeof(PetscScalar));CHKERRQ(ierr);
|
ierr = PetscMemcpy(ps->mat[mat]+(k_)*ld,Q+(k_)*ld,ld*sizeof(PetscScalar));CHKERRQ(ierr);
|
}else{ /* complex block */
|
}else{ /* complex block */
|
if(ps->compact){
|
if(ps->compact){
|
s1 = *(ps->rmat[PS_MAT_T]+2*ld+k_);
|
s1 = *(ps->rmat[PS_MAT_T]+2*ld+k_);
|
d1 = *(ps->rmat[PS_MAT_T]+k_);
|
d1 = *(ps->rmat[PS_MAT_T]+k_);
|
| Line 314... |
Line 314... |
/* real eigenvalue */
|
/* real eigenvalue */
|
wr[j] = A[j+j*ld]/B[j+j*ld];
|
wr[j] = A[j+j*ld]/B[j+j*ld];
|
wi[j] = 0.0 ;
|
wi[j] = 0.0 ;
|
} else {
|
} else {
|
/* diagonal block */
|
/* diagonal block */
|
d1 = A[j+j*ld]/B[j+j*ld]; d2 = A[(j+1)+(j+1)*ld]/B[(j+1)+(j+1)*ld];
|
d1 = PetscRealPart(A[j+j*ld])/PetscRealPart(B[j+j*ld]);
|
e1 = A[j+(j+1)*ld]/B[j+j*ld]; e2 = A[(j+1)+j*ld]/B[(j+1)+(j+1)*ld];
|
d2 = PetscRealPart(A[(j+1)+(j+1)*ld])/PetscRealPart(B[(j+1)+(j+1)*ld]);
|
|
e1 = PetscRealPart(A[j+(j+1)*ld])/PetscRealPart(B[j+j*ld]);
|
|
e2 = PetscRealPart(A[(j+1)+j*ld])/PetscRealPart(B[(j+1)+(j+1)*ld]);
|
wr[j] = (d1+d2)/2; wr[j+1] = wr[j];
|
wr[j] = (d1+d2)/2; wr[j+1] = wr[j];
|
disc = (d1-d2)*(d1-d2) - (e1-e2)*(e1-e2);
|
disc = (d1-d2)*(d1-d2) - (e1-e2)*(e1-e2);
|
if (disc<0){ /* complex eigenvalues */
|
if (disc<0.0){ /* complex eigenvalues */
|
wi[j] = PetscSqrtReal(-disc)/2; wi[j+1] = -wi[j];
|
wi[j] = PetscSqrtReal(-disc)/2; wi[j+1] = -wi[j];
|
}else{ /* real eigenvalues */
|
}else{ /* real eigenvalues */
|
disc = PetscSqrtReal(disc)/2;
|
disc = PetscSqrtReal(disc)/2;
|
wr[j] = wr[j]+disc; wr[j+1]=wr[j+1]-disc; wi[j] = 0.0; wi[j+1] = 0.0;
|
wr[j] = wr[j]+disc; wr[j+1]=wr[j+1]-disc; wi[j] = 0.0; wi[j+1] = 0.0;
|
}
|
}
|
| Line 346... |
Line 348... |
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "PSlarth_private"
|
#define __FUNCT__ "PSlarth_private"
|
static PetscErrorCode PSlarth_private(PetscReal x1,PetscReal x2,PetscInt *type,PetscReal *c,PetscReal *s,PetscReal *r)
|
static PetscErrorCode PSlarth_private(PetscReal x1,PetscReal x2,PetscInt *type,PetscReal *c,PetscReal *s,PetscReal *r)
|
{
|
{
|
PetscReal t,n2,xa,xb;
|
PetscReal t,n2,xa,xb;
|
|
PetscInt type_;
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
if(x2==0) {
|
if(x2==0) {
|
*c = 1.0; *s = 0.0; *r = PetscAbsReal(x1); *type = 1;
|
*c = 1.0; *s = 0.0; *r = PetscAbsReal(x1); *type = 1;
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
| Line 358... |
Line 361... |
*c = 0; *s = 0; *r = 0; *type = 0;
|
*c = 0; *s = 0; *r = 0; *type = 0;
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
if(PetscAbsReal(x1)>PetscAbsReal(x2)){
|
if(PetscAbsReal(x1)>PetscAbsReal(x2)){
|
xa = x1; xb = x2; *type =1;
|
xa = x1; xb = x2; type_ =1;
|
}else{ xa = x2; xb = x1; *type =2; }
|
}else{ xa = x2; xb = x1; type_ =2;
|
|
}
|
t = xb/xa;
|
t = xb/xa;
|
n2 = PetscAbsReal(1 - t*t);
|
n2 = PetscAbsReal(1 - t*t);
|
*r = PetscSqrtReal(n2)*PetscAbsReal(xa);
|
*r = PetscSqrtReal(n2)*PetscAbsReal(xa);
|
*c = x1/(*r);
|
*c = x1/(*r);
|
*s = x2/(*r);
|
*s = x2/(*r);
|
|
if(type_ == 2) *r *= -1;
|
|
if(type) *type = type_;
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
/*
|
/*
|
| c -s|
|
| c -s|
|
| Line 376... |
Line 382... |
|c s|
|
|c s|
|
[X1 X2]|s c|
|
[X1 X2]|s c|
|
*/
|
*/
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "PSRoth_private"
|
#define __FUNCT__ "PSRoth_private"
|
static PetscErrorCode PSRoth_private(PetscInt n, PetscScalar *X1, PetscScalar *X2, PetscReal c, PetscReal s, PetscInt type)
|
static PetscErrorCode PSRoth_private(PetscInt n, PetscScalar *X1, PetscScalar *X2, PetscReal c, PetscReal s)
|
{
|
{
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
|
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
| Line 391... |
Line 397... |
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "ArrowTridiagDiag"
|
#define __FUNCT__ "ArrowTridiagDiag"
|
static PetscErrorCode ArrowTridiagDiag(PetscInt k,PetscInt n,PetscReal *d,PetscReal *e,PetscReal *Omega,PetscScalar *Q,PetscInt ldq)
|
static PetscErrorCode ArrowTridiagDiag(PetscInt k,PetscInt n,PetscReal *d,PetscReal *e,PetscReal *Omega,PetscScalar *Q,PetscInt ldq)
|
{
|
{
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
PetscBLASInt i,j,j2,ld=ldq,one=1;
|
PetscBLASInt j2,ld=ldq,one=1;
|
PetscReal c,s,p,off,temp,e1,e0,d1,d0;
|
PetscInt i,j,type;
|
|
PetscReal c,s,p,off,e1,e0,d1,d0,temp;
|
|
PetscErrorCode ierr;
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
if (n<=2) PetscFunctionReturn(0);
|
if (n<=2) PetscFunctionReturn(0);
|
|
|
for (j=0;j<n-2;j++) {
|
for (j=0;j<n-2;j++) {
|
|
|
/* Eliminate entry e(j) by a rotation in the planes (j,j+1) */
|
/* Eliminate entry e(j) by a rotation in the planes (j,j+1) */
|
|
type = (Omega[j]*Omega[j+1]>0.0)?1:-1;
|
|
if(PetscAbsReal(e[j+1]) < PetscAbsReal(e[j])) type = 2;
|
e0 = e[j]; e1 = e[j+1];
|
e0 = e[j]; e1 = e[j+1];
|
d0 = d[j]*Omega[j]; d1 = d[j+1]*Omega[j+1];
|
d0 = d[j]; d1 = d[j+1];
|
if( Omega[j+1]*Omega[j]>0.0){ /* unitary rotation */
|
temp = e[j+1];
|
LAPACKlartg_(&e1,&e0,&c,&s,&e[j+1]);
|
if( type > 0){ /* unitary rotator */
|
|
LAPACKlartg_(&temp,&e[j],&c,&s,&e[j+1]);
|
|
}else{/* hyperbolic rotator */
|
|
ierr = PSlarth_private(temp,e[j],PETSC_NULL,&c,&s,&e[j+1]);CHKERRQ(ierr);
|
|
}
|
|
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 diagonal elements */
|
/* Apply rotation to Q */
|
e[j] = Omega[j]*c*s*(d0-d1);
|
j2 = j+2;
|
d[j+1] = Omega[j+1]*(s*s*d0 + c*c*d1);
|
if(type == 0) {
|
d[j] = Omega[j]*(c*c*d0 + s*s*d1);
|
|
|
|
/* Apply rotation to Q */
|
|
j2 = j+2;
|
|
s = -s;
|
|
BLASrot_(&j2,Q+j*ld,&one,Q+(j+1)*ld,&one,&c,&s);
|
BLASrot_(&j2,Q+j*ld,&one,Q+(j+1)*ld,&one,&c,&s);
|
}else{
|
}else{
|
|
ierr = PSRoth_private(j2, Q+j*ld,Q+(j+1)*ld,c,s);CHKERRQ(ierr);
|
}
|
}
|
/* Chase newly introduced off-diagonal entry to the top left corner */
|
/* Chase newly introduced off-diagonal entry to the top left corner */
|
for (i=j-1;i>=0;i--) {
|
for (i=j-1;i>=0;i--) {
|
off = -s*e[i];
|
off = -type*s*e[i];
|
e[i] = c*e[i];
|
e[i] = c*e[i];
|
temp = e[i+1];
|
temp = e[i+1];
|
|
type = (Omega[i]*Omega[i+1]>0.0)?1:-1;
|
|
if(type > 0){
|
LAPACKlartg_(&temp,&off,&c,&s,&e[i+1]);
|
LAPACKlartg_(&temp,&off,&c,&s,&e[i+1]);
|
s = -s;
|
}else{
|
temp = (d[i]-d[i+1])*s - 2.0*c*e[i];
|
ierr = PSlarth_private(e[i+1],off,PETSC_NULL,&c,&s,&e[i+1]);CHKERRQ(ierr);
|
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);
|
|
}
|
}
|
|
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);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
|
|
| Line 556... |
Line 574... |
ierr = PSOrthog_private(s+ps->l, Q+j*ld+ps->l, ss[j],Q+i*ld+ps->l, PETSC_NULL,n1);CHKERRQ(ierr);
|
ierr = PSOrthog_private(s+ps->l, Q+j*ld+ps->l, ss[j],Q+i*ld+ps->l, PETSC_NULL,n1);CHKERRQ(ierr);
|
}
|
}
|
}
|
}
|
ierr = PSNormIndef_private(s+ps->l,Q+i*ld+ps->l,&d1,n1);CHKERRQ(ierr);
|
ierr = PSNormIndef_private(s+ps->l,Q+i*ld+ps->l,&d1,n1);CHKERRQ(ierr);
|
ss[i] = (d1<0.0)?-1:1;
|
ss[i] = (d1<0.0)?-1:1;
|
d[i] = wr[i]*ss[i]; e[i] = 0.0;
|
d[i] = PetscRealPart(wr[i])*ss[i]; e[i] = 0.0;
|
}else{
|
}else{
|
for(j=i-1;j>=ps->l;j--){
|
for(j=i-1;j>=ps->l;j--){
|
/* s-orthogonalization of Qi and Qi+1*/
|
/* s-orthogonalization of Qi and Qi+1*/
|
if(PetscAbsScalar(wr[j]-wr[i])<toldeg && PetscAbsScalar(PetscAbsScalar(wi[j])-PetscAbsScalar(wi[i]))<toldeg){
|
if(PetscAbsScalar(wr[j]-wr[i])<toldeg && PetscAbsScalar(PetscAbsScalar(wi[j])-PetscAbsScalar(wi[i]))<toldeg){
|
ierr = PSOrthog_private(s+ps->l, Q+j*ld+ps->l, ss[j],Q+i*ld+ps->l, PETSC_NULL,n1);CHKERRQ(ierr);
|
ierr = PSOrthog_private(s+ps->l, Q+j*ld+ps->l, ss[j],Q+i*ld+ps->l, PETSC_NULL,n1);CHKERRQ(ierr);
|
| Line 570... |
Line 588... |
ierr = PSNormIndef_private(s+ps->l,Q+i*ld+ps->l,&d1,n1);CHKERRQ(ierr);
|
ierr = PSNormIndef_private(s+ps->l,Q+i*ld+ps->l,&d1,n1);CHKERRQ(ierr);
|
ss[i] = (d1<0)?-1:1;
|
ss[i] = (d1<0)?-1:1;
|
ierr = PSOrthog_private(s+ps->l, Q+i*ld+ps->l, ss[i],Q+(i+1)*ld+ps->l, &h,n1);CHKERRQ(ierr);
|
ierr = PSOrthog_private(s+ps->l, Q+i*ld+ps->l, ss[i],Q+(i+1)*ld+ps->l, &h,n1);CHKERRQ(ierr);
|
ierr = PSNormIndef_private(s+ps->l,Q+(i+1)*ld+ps->l,&d2,n1);CHKERRQ(ierr);
|
ierr = PSNormIndef_private(s+ps->l,Q+(i+1)*ld+ps->l,&d2,n1);CHKERRQ(ierr);
|
ss[i+1] = (d2<0)?-1:1;
|
ss[i+1] = (d2<0)?-1:1;
|
d[i] = (wr[i]-wi[i]*h/d1)*ss[i]; d[i+1] = (wr[i]+wi[i]*h/d1)*ss[i+1];
|
d[i] = (PetscRealPart(wr[i])-PetscRealPart(wi[i])*h/d1)*ss[i];
|
e[i] = wi[i]*d2/d1*ss[i]; e[i+1] = 0.0;
|
d[i+1] = (PetscRealPart(wr[i])+PetscRealPart(wi[i])*h/d1)*ss[i+1];
|
|
e[i] = PetscRealPart(wi[i])*d2/d1*ss[i]; e[i+1] = 0.0;
|
i++;
|
i++;
|
}
|
}
|
}
|
}
|
for(i=ps->l;i<ps->n;i++) s[i] = ss[i];
|
for(i=ps->l;i<ps->n;i++) s[i] = ss[i];
|
ps->k = ps->l;
|
ps->k = ps->l;
|
| Line 708... |
Line 727... |
ierr = PSOrthog_private(s+ps->l, Q+j*ld+ps->l, ss[j],Q+i*ld+ps->l, PETSC_NULL,n1);CHKERRQ(ierr);
|
ierr = PSOrthog_private(s+ps->l, Q+j*ld+ps->l, ss[j],Q+i*ld+ps->l, PETSC_NULL,n1);CHKERRQ(ierr);
|
}
|
}
|
}
|
}
|
ierr = PSNormIndef_private(s+ps->l,Q+i*ld+ps->l,&h,n1);CHKERRQ(ierr);
|
ierr = PSNormIndef_private(s+ps->l,Q+i*ld+ps->l,&h,n1);CHKERRQ(ierr);
|
ss[i] = (h<0)?-1:1;
|
ss[i] = (h<0)?-1:1;
|
d[i] = wr[i]*ss[i]; e[i] = 0.0;
|
d[i] = PetscRealPart(wr[i])*ss[i]; e[i] = 0.0;
|
}else{
|
}else{
|
for(j=i-1;j>=ps->l;j--){
|
for(j=i-1;j>=ps->l;j--){
|
/* s-orthogonalization of Qi and Qi+1*/
|
/* s-orthogonalization of Qi and Qi+1*/
|
if(PetscAbsScalar(wr[j]-wr[i])<toldeg && PetscAbsReal(PetscAbsScalar(wi[j])-PetscAbsScalar(wi[i]))<toldeg){
|
if(PetscAbsScalar(wr[j]-wr[i])<toldeg && PetscAbsReal(PetscAbsScalar(wi[j])-PetscAbsScalar(wi[i]))<toldeg){
|
ierr = PSOrthog_private(s+ps->l, Q+j*ld+ps->l, ss[j],Q+i*ld+ps->l, PETSC_NULL,n1);CHKERRQ(ierr);
|
ierr = PSOrthog_private(s+ps->l, Q+j*ld+ps->l, ss[j],Q+i*ld+ps->l, PETSC_NULL,n1);CHKERRQ(ierr);
|
| Line 722... |
Line 741... |
ierr = PSNormIndef_private(s+ps->l,Q+i*ld+ps->l,&d1,n1);CHKERRQ(ierr);
|
ierr = PSNormIndef_private(s+ps->l,Q+i*ld+ps->l,&d1,n1);CHKERRQ(ierr);
|
ss[i] = (d1<0)?-1:1;
|
ss[i] = (d1<0)?-1:1;
|
ierr = PSOrthog_private(s+ps->l, Q+i*ld+ps->l, ss[i],Q+(i+1)*ld+ps->l, &h,n1);CHKERRQ(ierr);
|
ierr = PSOrthog_private(s+ps->l, Q+i*ld+ps->l, ss[i],Q+(i+1)*ld+ps->l, &h,n1);CHKERRQ(ierr);
|
ierr = PSNormIndef_private(s+ps->l,Q+(i+1)*ld+ps->l,&d2,n1);CHKERRQ(ierr);
|
ierr = PSNormIndef_private(s+ps->l,Q+(i+1)*ld+ps->l,&d2,n1);CHKERRQ(ierr);
|
ss[i+1] = (d2<0)?-1:1;
|
ss[i+1] = (d2<0)?-1:1;
|
d[i] = (wr[i]-wi[i]*h/d1)*ss[i]; d[i+1] = (wr[i]+wi[i]*h/d1)*ss[i+1];
|
d[i] = (PetscRealPart(wr[i])-PetscRealPart(wi[i])*h/d1)*ss[i];
|
e[i] = wi[i]*d2/d1*ss[i]; e[i+1] = 0.0;
|
d[i+1] = (PetscRealPart(wr[i])+PetscRealPart(wi[i])*h/d1)*ss[i+1];
|
|
e[i] = PetscRealPart(wi[i])*d2/d1*ss[i]; e[i+1] = 0.0;
|
i++;
|
i++;
|
}
|
}
|
}
|
}
|
for(i=ps->l;i<ps->n;i++) s[i] = ss[i];
|
for(i=ps->l;i<ps->n;i++) s[i] = ss[i];
|
ps->k = ps->l;
|
ps->k = ps->l;
|
| Line 807... |
Line 827... |
n = ps->n;
|
n = ps->n;
|
d = ps->rmat[PS_MAT_T];
|
d = ps->rmat[PS_MAT_T];
|
e = d + ps->ld;
|
e = d + ps->ld;
|
s = d + 2*ps->ld;
|
s = d + 2*ps->ld;
|
ierr = PSAllocateWork_Private(ps,ps->ld,ps->ld,ps->ld);CHKERRQ(ierr);
|
ierr = PSAllocateWork_Private(ps,ps->ld,ps->ld,ps->ld);CHKERRQ(ierr);
|
perm = ps->iwork;
|
perm = ps->perm;
|
ierr = PSSortEigenvalues_Private(ps,wr,wi,perm,comp_func,comp_ctx);CHKERRQ(ierr);
|
ierr = PSSortEigenvalues_Private(ps,wr,wi,perm,comp_func,comp_ctx);CHKERRQ(ierr);
|
ierr = PetscMemcpy(ps->work,wr,n*sizeof(PetscScalar));CHKERRQ(ierr);
|
ierr = PetscMemcpy(ps->work,wr,n*sizeof(PetscScalar));CHKERRQ(ierr);
|
for (i=ps->l;i<n;i++) {
|
for (i=ps->l;i<n;i++) {
|
wr[i] = *(ps->work + perm[i]);
|
wr[i] = *(ps->work + perm[i]);
|
}
|
}
|