Skip to content

Instantly share code, notes, and snippets.

@stiv-yakovenko

stiv-yakovenko/rust_native_real_eigen Secret

Created September 2, 2017 17:56
Show Gist options
  • Save stiv-yakovenko/2c59c9588b615556f20d4a2537f6cd8d to your computer and use it in GitHub Desktop.
Save stiv-yakovenko/2c59c9588b615556f20d4a2537f6cd8d to your computer and use it in GitHub Desktop.
Download ZIP
rust compute eigen values with no c/c++ code dependency
Raw
rust_native_real_eigen
/*
This piece of code is transpiled EigenvalueDecomposition.java from Jama framework.
I had to do this, because I haven't found any buildable rust opensource project
to calculate real eigen values. There are many packages which are usually bound to
openBLAS, but I can't build them with gnu toolchain, and that's the only toolchain, which
allows gdb (and thus IDE) debug nowadays.
Quality code is far from perfect, hopefully someone will appreciate my one day of
manual code conversion nightmare and mention me in the source code.
Stepan Yakovenko,
https://github.com/stiv-yakovenko
*/
use rulinalg::matrix::BaseMatrix;
use num_traits::float::Float;
extern crate num_traits;
extern crate rulinalg;
use rulinalg::matrix::Matrix;
fn cdiv(xr: f64, xi: f64, yr: f64, yi: f64) -> (f64, f64) {
let r: f64;
let d: f64;
if yr.abs() > yi.abs() {
r = yi / yr;
d = yr + r * yi;
((xr + r * xi) / d, (xi - r * xr) / d)
} else {
r = yr / yi;
d = yi + r * yr;
((r * xr + xi) / d, (r * xi - xr) / d)
}
}
pub fn max(a: usize, b: usize) -> usize {
if a > b { a } else { b }
}
pub fn maxi16(a: i16, b: i16) -> i16 {
if a > b { a } else { b }
}
pub fn min(a: usize, b: usize) -> usize {
if a < b { a } else { b }
}
pub fn hqr2(n_in: usize, H: &mut Matrix<f64>, V: &mut Matrix<f64>, d: &mut Vec<f64>, e: &mut Vec<f64>) {
// 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.
// Initialize
let nn = n_in;
let mut n = nn as i16 - 1 ;
let low = 0;
let high = nn - 1;
let eps = (2.0).powf(-52.0);
let mut exshift = 0.0;
let mut p = 0.;
let mut q = 0.;
let mut r = 0.;
let mut s = 0.;
let mut z = 0.;
let mut t;
let mut w;
let mut x;
let mut y;
// Store roots isolated by balanc and compute matrix norm
let mut norm = 0.0;
let mut i = 0 as usize;
while i < nn {
if i < low || i > high {
d[i] = H[[i , i]];
e[i] = 0.0;
}
let mut j = maxi16(i as i16 - 1, 0) as usize;
while j < nn {
norm = norm + (H[[i, j]]).abs();
j = j + 1;
}
i = i + 1;
}
// Outer loop over eigenvalue index
let mut iter = 0;
while n >= low as i16 {
// Look for single small sub-diagonal element
let mut l = n;
while l > low as i16 {
s = (H[[l as usize - 1, l as usize - 1]]).abs() + (H[[l as usize , l as usize ]]).abs();
if s == 0.0 {
s = norm;
}
if (H[[l as usize , l as usize - 1]]).abs() < eps * s {
break;
}
l = l - 1;
}
// Check for convergence
// One root found
if l == n {
H[[n as usize , n as usize ]] = H[[n as usize , n as usize ]] + exshift;
d[n as usize ] = H[[n as usize , n as usize ]];
e[n as usize ] = 0.0;
n = n - 1;
iter = 0;
// Two roots found
} else if l == n - 1 {
w = H[[n as usize , n as usize - 1]] * H[[n as usize- 1, n as usize]];
p = (H[[n as usize- 1, n as usize- 1]] - H[[n as usize, n as usize]]) / 2.0;
q = p * p + w;
z = (q).abs().sqrt();
H[[n as usize, n as usize]] = H[[n as usize, n as usize]] + exshift;
H[[n as usize- 1, n as usize- 1]] = H[[n as usize- 1, n as usize- 1]] + exshift;
x = H[[n as usize, n as usize]];
// Real pair
if q >= 0. {
if p >= 0. {
z = p + z;
} else {
z = p - z;
}
d[n as usize - 1] = x + z;
d[n as usize] = d[n as usize - 1];
if z != 0.0 {
d[n as usize] = x - w / z;
}
e[n as usize - 1] = 0.0;
e[n as usize] = 0.0;
x = H[[n as usize, n as usize- 1]];
s = (x).abs() + (z).abs();
p = x / s;
q = z / s;
r = (p * p + q * q).sqrt();
p = p / r;
q = q / r;
// Row modification
let mut j = n - 1;
while j < nn as i16 {
z = H[[n as usize - 1, j as usize]];
H[[n as usize - 1, j as usize]] = q * z + p * H[[n as usize, j as usize]];
H[[n as usize, j as usize]] = q * H[[n as usize, j as usize]] - p * z;
j = j + 1;
}
// Column modification
let mut i = 0;
while i <= n {
z = H[[i as usize, n as usize - 1]];
H[[i as usize, n as usize - 1]] = q * z + p * H[[i as usize, n as usize]];
H[[i as usize, n as usize]] = q * H[[i as usize, n as usize]] - p * z;
i = i + 1;
}
// Accumulate transformations
let mut i = low;
while i <= high {
z = V[[i as usize, n as usize - 1]];
V[[i as usize, n as usize - 1]] = q * z + p * V[[i as usize, n as usize]];
V[[i as usize, n as usize]] = q * V[[i as usize, n as usize]] - p * z;
i = i + 1;
}
// Complex pair
} else {
d[n as usize - 1] = x + p;
d[n as usize] = x + p;
e[n as usize - 1] = z;
e[n as usize] = -z;
}
n = n - 2;
iter = 0;
// No convergence yet
} else {
// Form shift
x = H[[n as usize, n as usize]];
y = 0.0;
w = 0.0;
if l < n {
y = H[[n as usize - 1, n as usize- 1]];
w = H[[n as usize, n as usize- 1]] * H[[n as usize - 1, n as usize]];
}
// Wilkinson's original ad hoc shift
if iter == 10 {
exshift += x;
let mut i = low;
while i <= n as usize{
H[[i, i]] -= x;
i = i + 1;
}
s = (H[[n as usize, n as usize- 1]]).abs() + (H[[n as usize - 1, n as usize- 2]]).abs();
y = 0.75 * s;
x = y;
w = -0.4375 * s * s;
}
// MATLAB's new ad hoc shift
if iter == 30 {
s = (y - x) / 2.0;
s = s * s + w;
if s > 0. {
s = s.sqrt();
if y < x {
s = -s;
}
s = x - w / ((y - x) / 2.0 + s);
let mut i = low;
while i <= n as usize{
H[[i, i]] -= s;
i = i + 1;
}
exshift += s;
x = 0.964;
y = x;
w = y;
}
}
iter = iter + 1; // (Could check iteration count here.)
// Look for two consecutive small sub-diagonal elements
let mut m = n - 2;
while m >= l {
z = H[[m as usize, m as usize]];
r = x - z;
s = y - z;
p = (r * s - w) / H[[m as usize+ 1, m as usize]] + H[[m as usize, m as usize + 1]];
q = H[[m as usize+ 1, m as usize + 1]] - z - r - s;
r = H[[m as usize+ 2, m as usize+ 1]];
s = (p).abs() + (q).abs() + (r).abs();
p = p / s;
q = q / s;
r = r / s;
if m == l {
break;
}
if H[[m as usize, m as usize- 1]].abs() * (q).abs() + (r).abs() <
eps * ((p).abs() * ((H[[m as usize - 1, m as usize- 1]]).abs() + (z).abs() +
(H[[m as usize+ 1, m as usize+ 1]]).abs()
)) {
break;
}
m = m - 1;
}
let mut i = m + 2;
while i <= n {
H[[i as usize, i as usize - 2]] = 0.0;
if i > m + 2 {
H[[i as usize, i as usize- 3]] = 0.0;
}
i = i + 1;
}
// Double QR step involving rows l:n and columns m:n
let mut k = m;
while k <= n - 1 {
let notlast = if k != n - 1 { true } else { false };
if k != m {
p = H[[k as usize, k as usize - 1]];
q = H[[k as usize+ 1, k as usize - 1]];
r = if notlast { H[[k as usize + 2, k as usize- 1]] } else { 0.0 };
x = (p).abs() + (q).abs() + (r).abs();
if x == 0.0 {
continue;
}
p = p / x;
q = q / x;
r = r / x;
}
s = (p * p + q * q + r * r).sqrt();
if p < 0. {
s = -s;
}
if s != 0. {
if k != m {
H[[k as usize, k as usize - 1]] = -s * x;
} else if l != m {
H[[k as usize, k as usize - 1]] = -H[[k as usize, k as usize - 1]];
}
p = p + s;
x = p / s;
y = q / s;
z = r / s;
q = q / p;
r = r / p;
// Row modification
let mut j = k;
while j < nn as i16 {
p = H[[k as usize, j as usize]] + q * H[[k as usize+ 1, j as usize]];
if notlast {
p = p + r * H[[k as usize+ 2, j as usize]];
H[[k as usize+ 2, j as usize]] = H[[k as usize+ 2, j as usize]] - p * z;
}
H[[k as usize, j as usize]] = H[[k as usize, j as usize]] - p * x;
H[[k as usize+ 1, j as usize]] = H[[k as usize+ 1, j as usize]] - p * y;
j = j + 1;
}
// Column modification
let mut i = 0;
while i <= min(n as usize, k as usize + 3) {
p = x * H[[i, k as usize]] + y * H[[i as usize, k as usize+ 1]];
if notlast {
p = p + z * H[[i, k as usize+ 2]];
H[[i, k as usize+ 2]] = H[[i, k as usize+ 2]] - p * r;
}
H[[i, k as usize]] = H[[i, k as usize]] - p;
H[[i, k as usize+ 1]] = H[[i, k as usize+ 1]] - p * q;
i = i + 1;
}
// Accumulate transformations
let mut i = low;
while i <= high {
p = x * V[[i, k as usize]] + y * V[[i, k as usize + 1]];
if notlast {
p = p + z * V[[i as usize, k as usize+ 2]];
V[[i as usize, k as usize + 2]] = V[[i as usize, k as usize + 2]] - p * r;
}
V[[i, k as usize]] = V[[i, k as usize]] - p;
V[[i, k as usize + 1]] = V[[i, k as usize + 1]] - p * q;
i = i + 1;
}
} // (s != 0)
k = k + 1;
} // k loop
} // check convergence
} // while n >= low
// Backsubstitute to find vectors of upper triangular form
if norm == 0.0 {
return;
}
n = nn as i16 - 1;
while n >= 0 {
p = d[n as usize];
q = e[n as usize];
// Real vector
if q == 0. {
let mut l = n;
H[[n as usize, n as usize]] = 1.0;
let mut i = n as i16 - 1;
while i >= 0 {
w = H[[i as usize, i as usize]] - p;
r = 0.0;
let mut j = l;
while j <= n {
r = r + H[[i as usize, j as usize]] * H[[j as usize, n as usize]];
j = j + 1;
}
if e[i as usize] < 0.0 {
z = w;
s = r;
} else {
l = i;
if e[i as usize] == 0.0 {
if w != 0.0 {
H[[i as usize, n as usize]] = -r / w;
} else {
H[[i as usize, n as usize]] = -r / (eps * norm);
}
// Solve real equations
} else {
x = H[[i as usize, i as usize+ 1]];
y = H[[i as usize+ 1, i as usize]];
q = (d[i as usize] - p) * (d[i as usize] - p) + e[i as usize] * e[i as usize];
t = (x * s - z * r) / q;
H[[i as usize, n as usize]] = t;
if (x).abs() > (z).abs() {
H[[i as usize + 1, n as usize]] = (-r - w * t) / x;
} else {
H[[i as usize + 1, n as usize]] = (-s - y * t) / z;
}
}
// Overflow control
t = H[[i as usize, n as usize]];
if (eps * t).abs() * t > 1. {
let mut j = i;
while j <= n as i16 {
H[[j as usize, n as usize]] = H[[j as usize, n as usize]] / t;
j = j + 1;
}
}
}
i = i - 1;
}
// Complex vector
} else if q < 0. {
let mut l = n - 1;
// Last vector component imaginary so matrix is triangular
if (H[[n as usize, n as usize- 1]]).abs() > (H[[n as usize - 1, n as usize]]).abs() {
H[[n as usize- 1, n as usize- 1]] = q / H[[n as usize, n as usize- 1]];
H[[n as usize- 1, n as usize]] = -(H[[n as usize, n as usize]] - p) / H[[n as usize, n as usize - 1]];
} else {
let (cdivr, cdivi) = cdiv(0.0, -H[[n as usize - 1, n as usize]], H[[n as usize - 1, n as usize - 1]] - p, q);
H[[n as usize- 1, n as usize- 1]] = cdivr;
H[[n as usize - 1, n as usize]] = cdivi;
}
H[[n as usize, n as usize- 1]] = 0.0;
H[[n as usize, n as usize]] = 1.0;
let mut i = n - 2;
while i >= 0 {
let mut ra=0.;
let mut sa=0.;
let mut vr;
let vi;
let mut j = l;
while j <= n {
ra = ra + H[[i as usize, j as usize]] * H[[j as usize, n as usize - 1]];
sa = sa + H[[i as usize, j as usize]] * H[[j as usize, n as usize]];
j = j + 1;
}
w = H[[i as usize, i as usize]] - p;
if e[i as usize] < 0.0 {
z = w;
r = ra;
s = sa;
} else {
l = i;
if e[i as usize] == 0. {
let (cdivr, cdivi) = cdiv(-ra, -sa, w, q);
H[[i as usize, n as usize- 1]] = cdivr;
H[[i as usize, n as usize]] = cdivi;
} else {
// Solve complex equations
x = H[[i as usize, i as usize+ 1]];
y = H[[i as usize + 1, i as usize]];
vr = (d[i as usize] - p) * (d[i as usize] - p) + e[i as usize]* e[i as usize] - q * q;
vi = (d[i as usize] - p) * 2.0 * q;
if vr == 0.0 && vi == 0.0 {
vr = eps * norm * ((w).abs() + (q).abs() +
(x).abs() + (y).abs() + (z)).abs();
}
let (cdivr, cdivi) = cdiv(x * r - z * ra + q * sa, x * s - z * sa - q * ra, vr, vi);
H[[i as usize, n as usize- 1]] = cdivr;
H[[i as usize, n as usize]] = cdivi;
if (x).abs() > ((z).abs() + (q).abs()) {
H[[i as usize + 1, n as usize- 1]] = (-ra - w * H[[i as usize, n as usize - 1]] + q * H[[i as usize, n as usize]]) / x;
H[[i as usize+ 1, n as usize]] = (-sa - w * H[[i as usize, n as usize]] - q * H[[i as usize, n as usize - 1]]) / x;
} else {
let (cdivr, cdivi) = cdiv(-r - y * H[[i as usize, n as usize - 1]], -s - y * H[[i as usize, n as usize]], z, q);
H[[i as usize + 1, n as usize - 1]] = cdivr;
H[[i as usize + 1, n as usize]] = cdivi;
}
}
// Overflow control
t = (H[[i as usize, n as usize - 1]]).abs().max(H[[i as usize, n as usize]].abs());
if (eps * t) * t > 1. {
let mut j = i;
while j <= n {
j = j + 1;
H[[j as usize, n as usize - 1]] = H[[j as usize, n as usize- 1]] / t;
H[[j as usize, n as usize]] = H[[j as usize, n as usize]] / t;
}
}
}
i = i - 1;
}
}
n = n - 1;
}
// Vectors of isolated roots
let mut i = 0;
while i < nn {
if i < low || i > high {
let mut j = i;
while j < nn {
V[[i, j]] = H[[i, j]];
j = j + 1;
}
}
i = i + 1;
}
// Back transformation to get eigenvectors of original matrix
let mut j = nn as i16 - 1;
println!("low {:?}",low);
while j >= low as i16{
let mut i = low;
while i <= high {
z = 0.0;
let mut k = low;
while k <= min(j as usize, high) {
z = z + V[[i, k]] * H[[k, j as usize]];
k = k + 1;
}
V[[i, j as usize]] = z;
i = i + 1;
}
j = j - 1;
}
}
pub fn orthes(m: &mut Matrix<f64>, H: &mut Matrix<f64>, V: &mut Matrix<f64>) {
// This is derived from the Algol procedures orthes and ortran,
// by Martin and Wilkinson, Handbook for Auto. Comp.,
// Vol.ii-Linear Algebra, and the corresponding
// Fortran subroutines in EISPACK.
let low = 0;
let n = m.cols();
let high = n - 1;
let mut m = low + 1;
let mut ort = vec!(0.; n);
// for (int
// m = low + 1;
// m < = high - 1;
// m + + )
while m < high - 1 {
// Scale column.
let mut scale = 0.0;
let mut i = m;
//for (int i = m; i < = high; i + +)
while i <= high {
scale = scale + (H[[i, m - 1]]).abs();
i = i + 1;
}
if scale != 0.0 {
// Compute Householder transformation.
let mut h = 0.0;
let mut i = high;
while i >= m {
ort[i] = H[[i, m - 1]] / scale;
h += ort[i] * ort[i];
i = i - 1;
}
let mut g = h.sqrt();
if ort[m] > 0. {
g = -g;
}
h = h - ort[m] * g;
ort[m] = ort[m] - g;
// Apply Householder similarity transformation
// H = (I-u*u'/h)*H*(I-u*u')/h)
let mut j = m;
while j < n {
let mut f = 0.0;
let mut i = high;
while i >= m {
f += ort[i] * H[[i, j]];
i = i - 1;
}
f = f / h;
let mut i = m;
while
i <= high {
H[[i, j]] -= f * ort[i];
i = i + 1;
}
j = j + 1;
}
let mut i = 0;
while i <= high {
let mut f = 0.0;
let mut j = high;
while j >= m {
f += ort[j] * H[[i, j]];
j = j - 1;
}
f = f / h;
let mut j = m;
while j <= high {
H[[i, j]] -= f * ort[j];
j = j + 1;
}
i = i + 1;
}
ort[m] = scale * ort[m];
H[[m, m - 1]] = scale * g;
}
m = m + 1;
}
// Accumulate transformations (Algol's ortran).
for i in 0..n {
for j in 0..n {
V[[i, j]] = if i == j { 1.0 } else { 0.0 };
}
}
let mut m = high - 1;
while m >= low + 1 {
if H[[m, m - 1]] != 0.0 {
let mut i = m + 1;
while i <= high {
ort[i] = H[[i, m - 1]];
i = i + 1;
}
let mut j = m;
while j <= high {
let mut g = 0.0;
let mut i = m;
while i <= high {
g += ort[i] * V[[i, j]];
i = i + 1;
}
// Double division avoids possible underflow
g = (g / ort[m]) / H[[m, m - 1]];
let mut i = m;
while i <= high {
V[[i, j]] += g * ort[i];
i = i + 1;
}
j = j + 1;
}
}
m = m - 1;
}
}
fn calc_real_eigen(m:&mut Matrix<f64>) {
let n = m.cols();
let mut H = Matrix::new(m.cols(),m.cols(),vec!(0.;m.cols()*m.cols()));
let mut V = Matrix::new(m.cols(),m.cols(),vec!(0.;m.cols()*m.cols()));
//let mut ort = vec!(0.;n);
let mut d = vec!(0.;n);
let mut e = vec!(0.;n);
for i in 0..n {
for j in 0..n {
H[[i, j]] = m[[i,j]];
}
}
orthes(m,&mut H,&mut V);
hqr2(n,&mut H,&mut V,&mut d, &mut e);
println!("e={:?}",d);
}
fn main() {
println!("Hello, world!");
let mut m = Matrix::new(4, 4, vec![
0.0;16
]);
m[[1,0]]=1.;
m[[2,1]]=1.;
m[[3,2]]=1.;
m[[0,3]]=-9.;
m[[2,3]]=10.;
calc_real_eigen(&mut m);
// let m = Matrix::new(10);
// *m.get(5, 5) = 7.;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment