edu.mines.jtk.la.DMatrixEvd Maven / Gradle / Ivy
Show all versions of edu-mines-jtk Show documentation
/****************************************************************************
Copyright 2005, Colorado School of Mines and others.
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
****************************************************************************/
package edu.mines.jtk.la;
import static edu.mines.jtk.util.ArrayMath.*;
import edu.mines.jtk.util.Check;
/**
* Eigenvalue and eigenvector decomposition of a square matrix A.
*
* If A is symmetric, then A = V*D*V' where the matrix of eigenvalues D
* is diagonal and the matrix of eigenvectors V is orthogonal (V*V' = I).
*
* If A is not symmetric, then the eigenvalue matrix D is block diagonal
* with real eigenvalues in 1-by-1 blocks and any complex eigenvalues
* lambda + i*mu in 2-by-2 block [lambda, mu; -mu, lambda]. The columns
* of V represent the eigenvectors in the sense that A*V = V*D. The matrix
* V may be badly conditioned or even singular, so the validity of the
* equation A = V*D*inverse(V) depends on the condition number of V.
*
* This class was adapted from the package Jama, which was developed by
* Joe Hicklin, Cleve Moler, and Peter Webb of The MathWorks, Inc., and by
* Ronald Boisvert, Bruce Miller, Roldan Pozo, and Karin Remington of the
* National Institue of Standards and Technology.
* @author Dave Hale, Colorado School of Mines
* @version 2005.12.01
*/
public class DMatrixEvd {
/**
* Constructs an eigenvalue decomposition for the specified square matrix.
* @param a the square matrix
*/
public DMatrixEvd(DMatrix a) {
Check.argument(a.isSquare(),"matrix a is square");
double[][] aa = a.getArray();
int n = a.getN();
_n = n;
_v = new double[n][n];
_d = new double[n];
_e = new double[n];
if (a.isSymmetric()) {
for (int i=0; i0.0) {
da[i][i+1] = _e[i];
} else if (_e[i]<0.0) {
da[i][i-1] = _e[i];
}
}
return d;
}
/**
* Gets the real parts of the eigenvalues.
* @return array of real parts = real(diag(D)).
*/
public double[] getRealEigenvalues() {
return copy(_d);
}
/**
* Gets the imaginary parts of the eigenvalues
* @return array of imaginary parts = imag(diag(D))
*/
public double[] getImagEigenvalues() {
return copy(_e);
}
///////////////////////////////////////////////////////////////////////////
// private
private int _n; // row and column dimensions for square matrix V
private double[][] _v; // eigenvectors V
private double[] _d, _e; // eigenvalues
private double[][] _h; // nonsymmetric Hessenberg form
// Symmetric Householder reduction to tridiagonal form.
// This is derived from the Algol procedures tred2 by
// Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
// Fortran subroutine in EISPACK.
private void tred2() {
int n = _n;
for (int j=0; j0; --i) {
// Scale to avoid under/overflow.
double scale = 0.0;
double h = 0.0;
for (int k=0; k0.0)
g = -g;
_e[i] = scale*g;
h -= f * g;
_d[i-1] = f-g;
for (int j=0; jl) {
//int iter = 0;
do {
//++iter; // (Could check iteration count here.)
// Compute implicit shift
double g = _d[l];
double p = (_d[l+1] - g) / (2.0 * _e[l]);
double r = hypot(p,1.0);
if (p<0)
r = -r;
_d[l] = _e[l]/(p+r);
_d[l+1] = _e[l]*(p+r);
double dl1 = _d[l+1];
double h = g-_d[l];
for (int i=l+2; i=l; --i) {
c3 = c2;
c2 = c;
s2 = s;
g = c * _e[i];
h = c * p;
r = hypot(p,_e[i]);
_e[i+1] = s*r;
s = _e[i]/r;
c = p/r;
p = c*_d[i]-s*g;
_d[i+1] = h+s*(c*g+s*_d[i]);
// Accumulate transformation.
for (int k=0; keps*tst1);
}
_d[l] += f;
_e[l] = 0.0;
}
// Sort eigenvalues and corresponding vectors.
for (int i=0; i=m; --i) {
ort[i] = _h[i][m-1]/scale;
h += ort[i]*ort[i];
}
double g = sqrt(h);
if (ort[m]>0.0)
g = -g;
h -= ort[m]*g;
ort[m] -= g;
// Householder similarity transformation H = (I-u*u'/h)*H*(I-u*u')/h).
for (int j=m; j=m; --i)
f += ort[i]*_h[i][j];
f /= h;
for (int i=m; i<=high; ++i)
_h[i][j] -= f*ort[i];
}
for (int i=0; i<=high; ++i) {
double f = 0.0;
for (int j=high; j>=m; --j)
f += ort[j]*_h[i][j];
f /= h;
for (int j=m; j<=high; ++j)
_h[i][j] -= f*ort[j];
}
ort[m] = scale*ort[m];
_h[m][m-1] = scale*g;
}
}
// Accumulate transformations (Algol's ortran).
for (int i=0; i=low+1; --m) {
if (_h[m][m-1]!=0.0) {
for (int i=m+1; i<=high; ++i)
ort[i] = _h[i][m-1];
for (int j=m; j<=high; ++j) {
double g = 0.0;
for (int i=m; i<=high; ++i)
g += ort[i]*_v[i][j];
g = (g/ort[m])/_h[m][m-1]; // double division avoids underflow
for (int i=m; i<=high; ++i)
_v[i][j] += g*ort[i];
}
}
}
}
// Complex scalar division.
private double _cdivr, _cdivi;
private void cdiv(double xr, double xi, double yr, double yi) {
double r,d;
if (abs(yr) > abs(yi)) {
r = yi/yr;
d = yr+r*yi;
_cdivr = (xr+r*xi)/d;
_cdivi = (xi-r*xr)/d;
} else {
r = yr/yi;
d = yi+r*yr;
_cdivr = (r*xr+xi)/d;
_cdivi = (r*xi-xr)/d;
}
}
// Nonsymmetric reduction from Hessenberg to real Schur form.
// This is derived from the Algol procedure hqr2,
// by Martin and Wilkinson, Handbook for Auto. Comp.,
// Vol.ii-Linear Algebra, and the corresponding
// Fortran subroutine in EISPACK.
private void hqr2() {
// Initialize
int nn = _n;
int n = nn-1;
int low = 0;
int high = nn-1;
double eps = pow(2.0,-52.0);
double exshift = 0.0;
double p=0.0,q=0.0,r=0.0,s=0.0,z=0.0,t,w,x,y;
// Store roots isolated by balance and compute matrix norm
double norm = 0.0;
for (int i=0; ihigh) {
_d[i] = _h[i][i];
_e[i] = 0.0;
}
for (int j=max(i-1,0); j=low) {
// Look for single small sub-diagonal element
int l = n;
while (l>low) {
s = abs(_h[l-1][l-1])+abs(_h[l][l]);
if (s==0.0)
s = norm;
if (abs(_h[l][l-1])=0) {
if (p>=0) {
z = p+z;
} else {
z = p-z;
}
_d[n-1] = x+z;
_d[n] = _d[n-1];
if (z!=0.0)
_d[n] = x-w/z;
_e[n-1] = 0.0;
_e[n] = 0.0;
x = _h[n][n-1];
s = abs(x)+abs(z);
p = x/s;
q = z/s;
r = sqrt(p*p+q*q);
p /= r;
q /= r;
// Row modification
for (int j=n-1; j0.0) {
s = sqrt(s);
if (y=l) {
z = _h[m][m];
r = x-z;
s = y-z;
p = (r*s-w)/_h[m+1][m]+_h[m][m+1];
q = _h[m+1][m+1]-z-r-s;
r = _h[m+2][m+1];
s = abs(p)+abs(q)+abs(r);
p /= s;
q /= s;
r /= s;
if (m==l)
break;
if (abs(_h[m][m-1])*(abs(q)+abs(r)) <
eps*(abs(p)*(abs(_h[m-1][m-1])+abs(z)+abs(_h[m+1][m+1])))) {
break;
}
--m;
}
for (int i=m+2; i<=n; ++i) {
_h[i][i-2] = 0.0;
if (i>m+2)
_h[i][i-3] = 0.0;
}
// Double QR step involving rows l:n and columns m:n
for (int k=m; k<=n-1; ++k) {
boolean notlast = (k!=n-1);
if (k!=m) {
p = _h[k][k-1];
q = _h[k+1][k-1];
r = notlast?_h[k+2][k-1]:0.0;
x = abs(p)+abs(q)+abs(r);
if (x!=0.0) {
p /= x;
q /= x;
r /= x;
}
}
if (x==0.0)
break;
s = sqrt(p*p+q*q+r*r);
if (p<0.0)
s = -s;
if (s!=0.0) {
if (k!=m) {
_h[k][k-1] = -s*x;
} else if (l!=m) {
_h[k][k-1] = -_h[k][k-1];
}
p = p+s;
x = p/s;
y = q/s;
z = r/s;
q = q/p;
r = r/p;
// Row modification
for (int j=k; j=low)
// Backsubstitute to find vectors of upper triangular form
if (norm==0.0)
return;
for (n=nn-1; n>=0; --n) {
p = _d[n];
q = _e[n];
// If real vector, ...
if (q==0.0) {
int l = n;
_h[n][n] = 1.0;
for (int i=n-1; i>=0; --i) {
w = _h[i][i]-p;
r = 0.0;
for (int j=l; j<=n; ++j)
r = r+_h[i][j]*_h[j][n];
if (_e[i]<0.0) {
z = w;
s = r;
} else {
l = i;
if (_e[i]==0.0) {
if (w!=0.0) {
_h[i][n] = -r/w;
} else {
_h[i][n] = -r/(eps*norm);
}
}
// Solve real equations
else {
x = _h[i][i+1];
y = _h[i+1][i];
q = (_d[i]-p)*(_d[i]-p)+_e[i]*_e[i];
t = (x*s-z*r)/q;
_h[i][n] = t;
if (abs(x)>abs(z)) {
_h[i+1][n] = (-r-w*t)/x;
} else {
_h[i+1][n] = (-s-y*t)/z;
}
}
// Overflow control
t = abs(_h[i][n]);
if ((eps*t)*t>1.0) {
for (int j=i; j<=n; ++j)
_h[j][n] /= t;
}
}
}
}
// else if omplex vector, ...
else if (q<0.0) {
int l = n-1;
// Last vector component imaginary so matrix is triangular
if (abs(_h[n][n-1])>abs(_h[n-1][n])) {
_h[n-1][n-1] = q/_h[n][n-1];
_h[n-1][n] = -(_h[n][n]-p)/_h[n][n-1];
} else {
cdiv(0.0,-_h[n-1][n],_h[n-1][n-1]-p,q);
_h[n-1][n-1] = _cdivr;
_h[n-1][n] = _cdivi;
}
_h[n][n-1] = 0.0;
_h[n][n] = 1.0;
for (int i=n-2; i>=0; --i) {
double ra,sa,vr,vi;
ra = 0.0;
sa = 0.0;
for (int j=l; j<=n; ++j) {
ra += _h[i][j]*_h[j][n-1];
sa += _h[i][j]*_h[j][n];
}
w = _h[i][i]-p;
if (_e[i]<0.0) {
z = w;
r = ra;
s = sa;
} else {
l = i;
if (_e[i]==0.0) {
cdiv(-ra,-sa,w,q);
_h[i][n-1] = _cdivr;
_h[i][n] = _cdivi;
} else {
// Solve complex equations
x = _h[i][i+1];
y = _h[i+1][i];
vr = (_d[i]-p)*(_d[i]-p)+_e[i]*_e[i]-q*q;
vi = (_d[i]-p)*2.0*q;
if (vr==0.0 && vi==0.0)
vr = eps*norm*(abs(w)+abs(q)+abs(x)+abs(y)+abs(z));
cdiv(x*r-z*ra+q*sa,x*s-z*sa-q*ra,vr,vi);
_h[i][n-1] = _cdivr;
_h[i][n] = _cdivi;
if (abs(x)>(abs(z)+abs(q))) {
_h[i+1][n-1] = (-ra-w*_h[i][n-1]+q*_h[i][n])/x;
_h[i+1][n] = (-sa-w*_h[i][n]-q*_h[i][n-1])/x;
} else {
cdiv(-r-y*_h[i][n-1],-s-y*_h[i][n],z,q);
_h[i+1][n-1] = _cdivr;
_h[i+1][n] = _cdivi;
}
}
// Overflow control
t = max(abs(_h[i][n-1]),abs(_h[i][n]));
if ((eps*t)*t>1.0) {
for (int j=i; j< n; ++j) {
_h[j][n-1] /= t;
_h[j][n] /= t;
}
}
}
}
}
}
// Vectors of isolated roots
for (int i=0; ihigh) {
for (int j=i; j=low; --j) {
for (int i=low; i<=high; ++i) {
z = 0.0;
for (int k=low; k<=min(j,high); ++k)
z += _v[i][k]*_h[k][j];
_v[i][j] = z;
}
}
}
}