📄 mes_eig.h
字号:
}
///////////////////////////////////////////////////////////////////////////////
// schur -- computes the Schur decomposition of the matrix A in situ
// -- optionally, gives Q matrix such that Q^T.A.Q is upper triangular
// -- returns upper triangular Schur matrix
// Note: The precision of the eigenvector and temporary vectors is
// critical in stability of the algorithm.
// Must be 'double'.
// Must be more careful in handling mathematical expressions
// like x/sqrt(x^2+y^2).
// Note: This implementation is very sensitive to the condition number of
// the input matrix.
///////////////////////////////////////////////////////////////////////////////
template <typename T> int Matrix<T>::mes_schur(Matrix<T>& AA,Matrix<T>& QQ,int maxIter){
int compQ=1;
int i, j, k, k_min, k_max, k_tmp, n, split, iter, isOK=1;
double beta2, c, discrim, dummy, nu1, s, tmp, x, y, z;
double sqrt_macheps;
static Vector<double> diag, beta;
if(AA.IsComplex())
mes_error(E_NEED_REAL,"mes_schur");
if(AA.IsDouble()==false)
mes_warn(WARN_SINGLE_PRECISION,"mes_schur");
if ( NumRows(AA) != NumCols(AA) || ( NumRows(QQ) != NumCols(QQ) ) )
mes_error(E_SQUARE,"schur");
if ( NumRows(QQ) != NumRows(AA) )
mes_error(E_SIZES,"schur");
n = NumCols(AA);
diag.Resize(NumCols(AA));
beta.Resize(NumCols(AA));
/* compute Hessenberg form */
isOK=Matrix<double>::mes_Hfactor(AA,diag,beta);
if(isOK==0) return isOK;
/* save QQ if necessary */
if ( compQ ){
isOK=Matrix<double>::mes_makeHQ(AA,diag,beta,QQ);
if(isOK==0) return isOK;
}
isOK=Matrix<double>::mes_makeH(AA,AA);
if(isOK==0) return isOK;
sqrt_macheps = sqrt(MACH_EPS);
k_min = 1;
while ( k_min <= n ){
double a00, a01, a10, a11;
double scale, t, numer, denom;
/* find k_max to suit: submatrix k_min..k_max should be irreducible */
k_max = n;
for ( k = k_min; k < k_max; k++ ){
if ( AA(k+1,k) == 0.0 ){
k_max = k; break;
}
}
if ( k_max <= k_min ){
k_min = k_max + 1;
continue; /* outer loop */
}
/* check to see if we have a 2 x 2 block with complex eigenvalues */
if ( k_max == k_min + 1 ){
a00 = AA(k_min,k_min);
a01 = AA(k_min,k_max);
a10 = AA(k_max,k_min);
a11 = AA(k_max,k_max);
tmp = a00 - a11;
discrim = tmp*tmp + 4*a01*a10;
if ( discrim < 0.0 ){ // yes -- e-vals are complex
// -- put 2 x 2 block in form [a b; c a];
// then eigenvalues have real part a & imag part sqrt(|bc|)
numer = - tmp;
double atmp = a01+a10;
denom = ( atmp >= 0.0 ) ? (atmp + pythagoras(atmp,tmp)) : (atmp - pythagoras(atmp,tmp));
if ( denom != 0.0 ){ // t = s/c = numer/denom
t = numer/denom;
scale = c = 1.0/pythagoras(1,t);
s = -c*t;
}
else{
c = 1.0;
s = 0.0;
}
mes_rot_cols(AA,k_min,k_max,c,s,AA);
mes_rot_rows(AA,k_min,k_max,c,s,AA);
if ( compQ )
mes_rot_cols(QQ,k_min,k_max,c,s,QQ);
k_min = k_max + 1;
continue;
}
else {// discrim >= 0; i.e. block has two real eigenvalues
// no -- e-vals are not complex;
// split 2 x 2 block and continue
// s/c = numer/denom
numer = ( tmp >= 0.0 ) ? - tmp - sqrt(discrim) : - tmp + sqrt(discrim);
denom = 2*a01;
if ( Abs(numer) < Abs(denom) ){ // t = s/c = numer/denom
t = numer/denom;
scale = c = 1.0/pythagoras(1,t);
s = -c*t;
}
else if ( numer != 0.0 ){ // t = c/s = denom/numer
t = denom/numer;
scale = 1.0/pythagoras(1,t);
c = Abs(t)*scale;
s = ( t >= 0.0 ) ? -scale : scale;
}
else{ // numer == denom == 0
c = 0.0;
s = -1.0;
}
mes_rot_cols(AA,k_min,k_max,c,s,AA);
mes_rot_rows(AA,k_min,k_max,c,s,AA);
if ( compQ )
mes_rot_cols(QQ,k_min,k_max,c,s,QQ);
k_min = k_max + 1; // go to next block
continue;
}
}
/* now have r x r block with r >= 2: apply Francis QR step until block splits */
split = 0;
iter = 0;
while ( ! split ){
iter++;
/* set up Wilkinson/Francis complex shift */
k_tmp = k_max - 1;
a00 = AA(k_tmp,k_tmp);
a01 = AA(k_tmp,k_max);
a10 = AA(k_max,k_tmp);
a11 = AA(k_max,k_max);
/* treat degenerate cases differently
-- if there are still no splits after five iterations
and the bottom 2 x 2 looks degenerate, force it to split */
if ( iter >= 5 &&
Abs(a00-a11) < sqrt_macheps*(Abs(a00)+Abs(a11)) &&
(Abs(a01) < sqrt_macheps*(Abs(a00)+Abs(a11)) ||
Abs(a10) < sqrt_macheps*(Abs(a00)+Abs(a11))) ) {
if ( Abs(a01) < sqrt_macheps*(Abs(a00)+Abs(a11)) )
AA[k_tmp][k_max]=0.0;
if ( Abs(a10) < sqrt_macheps*(Abs(a00)+Abs(a11)) ){
AA[k_max][k_tmp]=0.0;
split = 1;
continue;
}
}
s = a00 + a11;
t = a00*a11 - a01*a10;
/* break loop if a 2 x 2 complex block */
if ( k_max == k_min + 1 && s*s < 4.0*t ){
split = 1;
continue;
}
/* perturb shift if convergence is slow */
if ( (iter % 10) == 0 ) {
s += iter*0.02;
t += iter*0.02;
}
/* set up Householder transformations */
k_tmp = k_min + 1;
a00 = AA(k_min,k_min);
a01 = AA(k_min,k_tmp);
a10 = AA(k_tmp,k_min);
a11 = AA(k_tmp,k_tmp);
x = a00*a00 + a01*a10 - s*a00 + t;
y = a10*(a00+a11-s);
if ( k_min + 2 <= k_max )
z = a10*AA[k_min+2][k_tmp];
else
z = 0.0;
for ( k = k_min; k <= k_max-1; k++ ){
if ( k < k_max - 1 ){
mes_hhldr3(x,y,z,&nu1,&beta2,&dummy);
mes_hhldr3cols(AA,k,local_max(k-1,1),beta2,nu1,y,z);
mes_hhldr3rows(AA,k,local_min(n,k+3),beta2,nu1,y,z);
if ( compQ )
mes_hhldr3rows(QQ,k,n,beta2,nu1,y,z);
}
else{
mes_givens(x,y,&c,&s);
mes_rot_cols(AA,k,k+1,c,s,AA);
mes_rot_rows(AA,k,k+1,c,s,AA);
if ( compQ )
mes_rot_cols(QQ,k,k+1,c,s,QQ);
}
x = AA(k+1,k);
if ( k <= k_max - 2 )
y = AA(k+2,k);
else
y = 0.0;
if ( k <= k_max - 3 )
z = AA(k+3,k);
else
z = 0.0;
}
for ( k = k_min; k <= k_max-2; k++ ){
/* zero appropriate sub-diagonals */
AA[k+2][k]=0.0;
if ( k < k_max-2 )
AA[k+3][k]=0.0;
}
/* test to see if matrix should split */
for ( k = k_min; k < k_max; k++ ){
if ( (Abs(AA[k+1][k])) < (MACH_EPS*(Abs(AA[k][k])+Abs(AA[k+1][k+1]))) ){
AA[k+1][k] = 0.0;
split = 1;
}
}
if(iter==maxIter){
cerr << "ERROR: Too many iterations in mes_schur().\n";
isOK=0;
assert(0);
return isOK;
}
}
}
/* polish up AA by zeroing strictly lower triangular elements
and small sub-diagonal elements */
for ( i = 1; i <= NumRows(AA); i++ ){
for ( j = 1; j < i-1; j++ ){
AA[i][j] = 0.0;
}
}
for ( i = 1; i <= NumRows(AA) - 1; i++ ){
if ( (Abs(AA[i+1][i])) < (MACH_EPS*(Abs(AA[i][i])+Abs(AA[i+1][i+1]))) ){
AA[i+1][i] = 0.0;
}
}
return isOK;
}
///////////////////////////////////////////////////////////////////////////////
// schur_vals -- compute real & imaginary parts of eigenvalues
// -- assumes T contains a block upper triangular matrix
// as produced by schur()
// -- real parts stored in real_pt, imaginary parts in imag_pt
///////////////////////////////////////////////////////////////////////////////
template <typename T> int Matrix<T>::mes_schur_evals(const Matrix<T>& TT,Vector<T>& real_pt,Vector<T>& imag_pt){
int i, n;
double discrim;
double diff, sum, tmp;
if ( NumRows(TT) != NumCols(TT) )
mes_error(E_SQUARE,"schur_evals");
n = NumCols(TT);
real_pt.Resize((int)n);
imag_pt.Resize((int)n);
i = 1;
while ( i <= n ){
if ( i <= n-1 && TT[i+1][i] != 0.0 ){ /* should be a complex eigenvalue */
sum = 0.5*(TT[i][i]+TT[i+1][i+1]);
diff = 0.5*(TT[i][i]-TT[i+1][i+1]);
discrim = diff*diff + TT[i][i+1]*TT[i+1][i];
if ( discrim < 0.0 ){ /* yes -- complex e-vals */
real_pt[i] = real_pt[i+1] = sum;
imag_pt[i] = sqrt(-discrim);
imag_pt[i+1] = - imag_pt[i];
}
else{ /* no -- actually both real */
tmp = sqrt(discrim);
real_pt[i] = sum + tmp;
real_pt[i+1] = sum - tmp;
imag_pt[i] = imag_pt[i+1] = 0.0;
}
i += 2;
}
else{ /* real eigenvalue */
real_pt[i] = TT[i][i];
imag_pt[i] = 0.0;
i++;
}
}
return 1;
}
///////////////////////////////////////////////////////////////////////////////
// schur_vecs -- returns eigenvectors computed from the real Schur
// decomposition of a matrix
// -- T is the block upper triangular Schur matrix
// -- Q is the orthognal matrix where A = Q.T.Q^T
// -- if Q is null, the eigenvectors of T are returned
// -- X_re is the real part of the matrix of eigenvectors,
// and X_im is the imaginary part of the matrix.
// -- X_re is returned
// Note: The precision of the eigenvector and temporary vectors is
// critical in stability of the algorithm.
// Must be 'double'.
///////////////////////////////////////////////////////////////////////////////
template <typename T> int Matrix<T>::mes_schur_vecs(const Matrix<T>& TT,const Matrix<T>& QQ,Matrix<T>& XX_re,Matrix<T>& XX_im){
int i, j, limit;
double t11_re, t11_im, t12, t21, t22_re, t22_im;
double l_re, l_im, det_re, det_im, invdet_re, invdet_im,
val1_re, val1_im, val2_re, val2_im,
tmp_val1_re, tmp_val1_im, tmp_val2_re, tmp_val2_im;
double sum, diff, discrim, magdet, norm, scale;
static Vector<double> tmp1_re, tmp1_im, tmp2_re, tmp2_im;
if ( NumRows(TT) != NumCols(TT) || NumRows(XX_re) != NumCols(XX_re) || ( NumRows(QQ) != NumCols(QQ) ) || ( NumRows(XX_im) != NumCols(XX_im) ) )
mes_error(E_SQUARE,"schur_vecs");
if ( NumRows(TT) != NumRows(XX_re) || ( NumRows(TT) != NumRows(QQ) ) || ( NumRows(TT) != NumRows(XX_im) ) )
mes_error(E_SIZES,"schur_vecs");
tmp1_re.Resize(NumRows(TT));
tmp1_im.Resize(NumRows(TT));
tmp2_re.Resize(NumRows(TT));
tmp2_im.Resize(NumRows(TT));
i = 1;
while ( i <= NumRows(TT) ){
if ( i <= NumRows(TT)-1 && TT[i+1][i] != 0.0 ){ /* complex eigenvalue */
sum = 0.5*(TT[i][i]+TT[i+1][i+1]);
diff = 0.5*(TT[i][i]-TT[i+1][i+1]);
discrim = diff*diff + TT[i][i+1]*TT[i+1][i];
l_re = l_im = 0.0;
if ( discrim < 0.0 ){ /* yes -- complex e-vals */
l_re = sum;
l_im = sqrt(-discrim);
}
else /* not correct Real Schur form */
mes_error(E_RANGE,"schur_vecs");
}
else{
l_re = TT[i][i];
l_im = 0.0;
}
tmp1_im=0;
Matrix<double>::mes_v_rand(tmp1_re);
tmp1_re=MACH_EPS*tmp1_re;
/* solve (T-l.I)x = tmp1 */
limit = ( l_im != 0.0 ) ? i+1 : i;
for ( j = limit+1; j <= NumRows(TT); j++ )
tmp1_re[j] = 0.0;
j = limit;
while ( j >= 1 ){
double sum;
int k;
if ( j > 1 && TT[j][j-1] != 0.0 ){ /* 2 x 2 diagonal block */
for(sum=0,k=j+1;k<=limit;k++) sum += tmp1_re[k]*TT[j-1][k];
val1_re = tmp1_re[j-1] - sum;
for(sum=0,k=j+1;k<=limit;k++) sum += tmp1_im[k]*TT[j-1][k];
val1_im = tmp1_im[j-1] - sum;
for(sum=0,k=j+1;k<=limit;k++) sum += tmp1_re[k]*TT[j][k];
val2_re = tmp1_re[j] - sum;
for(sum=0,k=j+1;k<=limit;k++) sum += tmp1_im[k]*TT[j][k];
val2_im = tmp1_im[j] - sum;
t11_re = TT[j-1][j-1] - l_re;
t11_im = - l_im;
t22_re = TT[j][j] - l_re;
t22_im = - l_im;
t12 = TT[j-1][j];
t21 = TT[j][j-1];
scale = Abs(TT[j-1][j-1]) + Abs(TT[j][j]) + Abs(t12) + Abs(t21) + Abs(l_re) + Abs(l_im);
det_re = t11_re*t22_re - t11_im*t22_im - t12*t21;
det_im = t11_re*t22_im + t11_im*t22_re;
magdet = det_re*det_re+det_im*det_im;
if ( sqrt(magdet) < MACH_EPS*scale ){
det_re = MACH_EPS*scale;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -