Option Explicit
Option Base 0

DefLng I-N
DefDbl A-H, O-Z

Const RADIX = 2#

' A library of eigenvalue subroutines.

' The main user routine is:

' getEigen(A(), n, RealEigenVal(), ImagEigenVal#(), EigenVectors(), ierr)

' where: A() is an n x n positive definite matrix,
' RealEigenVal & ImagEigenVal are the real and imaginary parts
' of the eigenvalues,
' Eigenvectors is the n x n matrix whose columns are the
' eigenvectors of A().

' ierr=0 --> results OK
' ierr=1 --> non-zero imaginary components, no eigenvectors computed
' ierr=2 --> results may be inaccurate.

Sub balanc(A(), n)
Dim sqrdx, last, i, r, C, j, f, s, g

sqrdx = RADIX * RADIX
last = 0

While (last = 0)
last = 1
For i = 1 To n
r = 0#
C = 0#
For j = 1 To n
If (j <> i) Then
C = C + Abs(A(j, i))
r = r + Abs(A(i, j))
End If
Next j

If (C > 0# And r > 0#) Then
g = r / RADIX
f = 1#
s = C + r
While (C < g)
f = f * RADIX
C = C * sqrdx
Wend

g = r * RADIX
While (C > g)
f = f / RADIX
C = C / sqrdx
Wend
If ((C + r) / f < (0.95 * s)) Then
last = 0
g = 1# / f
For j = 1 To n
A(i, j) = A(i, j) * g
Next j
For j = 1 To n
A(j, i) = A(j, i) * f
Next j
End If
End If
Next i
Wend

End Sub

Sub elmhes(A(), n)
Dim m, X, i, j, y

For m = 2 To n - 1
X = 0#
i = m
For j = m To n
If Abs(A(j, m - 1)) > Abs(X) Then
X = A(j, m - 1)
i = j
End If
Next j

If i <> m Then
For j = m - 1 To n
swap A(i, j), A(m, j)
Next j
For j = 1 To n
swap A(j, i), A(j, m)
Next j
End If

If (X <> 0#) Then
For i = m + 1 To n
y = A(i, m - 1)
If y <> 0# Then
y = y / X
A(i, m - 1) = y
For j = m To n
A(i, j) = A(i, j) - y * A(m, j)
Next j
For j = 1 To n
A(j, m) = A(j, m) + y * A(j, i)
Next j
End If
Next i
End If
Next m

End Sub

Function EigenValue(AA As Variant) As Variant
Dim i, j, k, n, ierr
Dim A As Variant, BB(), xx(), real(), not_real() As Double

A = AA
n = UBound(A, 1)

ReDim real(n), not_real(n), xx(1 To n, 1 To n), BB(n, n)
For i = 1 To n
For j = 1 To n
BB(i, j) = A(i, j)
Next j
Next i

getEigen BB(), n, real(), not_real(), xx(), ierr

If ierr = 1 Then
'MsgBox "Complex Eigenvalues"
EigenValue = 0#
Exit Function
Else
For i = 1 To n
For j = 1 To n
xx(i, j) = 0#
Next j
xx(i, i) = real(i)
Next i
If ierr = 2 Then
'MsgBox "Inaccurate results"
End If
End If

EigenValue = xx()

End Function

Function EigenVector(AA As Variant) As Variant
Dim i, j, k, n, ierr
Dim A As Variant, BB(), xx(), real(), not_real() As Double

A = AA
n = UBound(A, 1)

ReDim real(n), not_real(n), xx(1 To n, 1 To n), BB(n, n)
For i = 1 To n
For j = 1 To n
BB(i, j) = A(i, j)
Next j
Next i

getEigen BB(), n, real(), not_real(), xx(), ierr

If ierr = 1 Then
'MsgBox "Complex Eigenvalues"
EigenVector = 0#
Exit Function
ElseIf ierr = 2 Then
'MsgBox "Inaccurate results"
End If

EigenVector = xx()
End Function

Sub getEigen(AAkeep(), n, RealEigenVal(), ImagEigenVal() As Double, EigenVectors(), ierrflag)
Dim i, j, k, sum, oldsum, d, iflag
ReDim AA(n, n)
ReDim indx(n), temp(n)

ierrflag = 0
For i = 1 To n
ImagEigenVal#(i) = 0#
RealEigenVal(i) = 0#
For j = 1 To n
AA(i, j) = AAkeep(i, j)
EigenVectors(i, j) = 0#
Next j
Next i

' Generate eigen values for matrix AA. This destroys AA.
balanc AA(), n
elmhes AA(), n
hqr AA(), n, RealEigenVal(), ImagEigenVal#()

For i = 1 To n
If ImagEigenVal#(i) <> 0# Then
ierrflag = 1 'Imaginary eigenvalues
Exit Sub
End If
Next i

'Generate eigen vectors.
For k = 1 To n
'Restore AA.
For i = 1 To n
For j = 1 To n
AA(i, j) = AAkeep(i, j)
Next j
Next i

'Adjust AA for kth eigen value.
For i = 1 To n
AA(i, i) = AA(i, i) - RealEigenVal(k)
Next i

'Generate random estimate of eigen vector.
RANDOMIZE TIMER
For i = 1 To n
temp(i) = RND
Next i

'Normalize our guess...
sum = 0#
For i = 1 To n
sum = sum + temp(i) * temp(i)
Next i
sum = SQR(sum)
For i = 1 To n
temp(i) = temp(i) / sum
Next i
sum = 1#

ludcmp n, d, AA(), indx(), iflag
j = 0
Do
lubksb n, AA(), indx(), temp()
oldsum = sum

sum = 0#
For i = 1 To n
sum = sum + temp(i) * temp(i)
Next i
sum = SQR(sum)
For i = 1 To n
temp(i) = temp(i) / sum
Next i
j = j + 1
Loop While j < 3

'Save eigen vector as kth column of eigenVectors.
For i = 1 To n
EigenVectors(i, k) = temp(i)
Next i
Next k

mMult AAkeep(), EigenVectors(), AA(), n, n, n

'Find largest absolute error between Ax/lamda and x.
sum = 0#
For i = 1 To n
For j = 1 To n
If Abs((AA(i, j) / RealEigenVal(j)) - EigenVectors(i, j)) > sum Then
sum = Abs((AA(i, j) / RealEigenVal(j)) - EigenVectors(i, j))
End If
Next j
Next i

If sum > 0.0000000001 Then
ierrflag = 2 'Eigenvalues/eigenvectors are poor.
End If

Erase AA()
Erase indx(), temp()
End Sub

Sub hqr(A(), n, wr(), wi())

Dim anorm, i, j, k, nn, t, its, L, s, X
Dim y, p, q, z, r, m, u, v, w, mmin

anorm = Abs(A(1, 1))

For i = 2 To n
For j = i - 1 To n
anorm = anorm + Abs(A(i, j))
Next j
Next i

nn = n
t = 0#

While (nn >= 1)
its = 0
Do
For L = nn To 2 Step -1
s = Abs(A(L - 1, L - 1)) + Abs(A(L, L))
If s = 0# Then s = anorm
If (Abs(A(L, L - 1)) + s) = s Then Exit For
Next L
X = A(nn, nn)
If L = nn Then
wr(nn) = X + t
nn = nn - 1
wi(nn) = 0#
Else
y = A(nn - 1, nn - 1)
w = A(nn, nn - 1) * A(nn - 1, nn)
If L = nn - 1 Then
p = 0.5 * (y - X)
q = p * p + w
z = SQR(Abs(q))
X = X + t
If q >= 0# Then
z = p + SIGN(z, p)
wr(nn) = X + z
wr(nn - 1) = wr(nn)
If z <> 0# Then wr(nn) = X - w / z
wi(nn) = 0#
wi(nn - 1) = 0#
Else
wr(nn) = X + p
wr(nn - 1) = wr(nn)
wi(nn) = z
wi(nn - 1) = -z
End If
nn = nn - 2
Else
If its = 30 Then
'Print "Too many iterations in HQR"
Stop
End If
If its = 10 Or its = 20 Then
t = t + X
For i = 1 To nn
A(i, i) = A(i, i) - X
Next i
s = Abs(A(nn, nn - 1)) + Abs(A(nn - 1, nn - 2))
X = 0.75 * s
y = X
w = -0.4375 * s * s
End If
its = its + 1
For m = nn - 2 To 1 Step -1
z = A(m, m)
r = X - z
s = y - z
p = (r * s - w) / A(m + 1, m) + A(m, m + 1)
q = A(m + 1, m + 1) - z - r - s
r = A(m + 2, m + 1)
s = Abs(p) + Abs(q) + Abs(r)
p = p / s
q = q / s
r = r / s
If m = 1 Then Exit For
u = Abs(A(m, m - 1)) * (Abs(q) + Abs(r))
v = Abs(p) * (Abs(A(m - 1, m - 1)) + Abs(z) + Abs(A(m + 1, m + 1)))
If u + v = v Then Exit For
Next m
For i = m + 2 To nn
A(i, i - 2) = 0#
If i <> m + 2 Then A(i, i - 3) = 0!
Next i
For k = m To nn - 1
If k <> m Then
p = A(k, k - 1)
q = A(k + 1, k - 1)
r = 0#
If k <> nn - 1 Then r = A(k + 2, k - 1)
X = Abs(p) + Abs(q) + Abs(r)
If X <> 0# Then
p = p / X
q = q / X
r = r / X
End If
End If
s = SIGN(SQR(p * p + q * q + r * r), p)
If s <> 0# Then
If k = m Then
If L <> m Then A(k, k - 1) = -A(k, k - 1)
Else
A(k, k - 1) = -s * X
End If
p = p + s
X = p / s
y = q / s
z = r / s
q = q / p
r = r / p
For j = k To nn
p = A(k, j) + q * A(k + 1, j)
If k <> nn - 1 Then
p = p + r * A(k + 2, j)
A(k + 2, j) = A(k + 2, j) - p * z
End If
A(k + 1, j) = A(k + 1, j) - p * y
A(k, j) = A(k, j) - p * X
Next j
If nn < k + 3 Then
mmin = nn
Else
mmin = k + 3
End If
For i = 1 To mmin
p = X * A(i, k) + y * A(i, k + 1)
If k <> nn - 1 Then
p = p + z * A(i, k + 2)
A(i, k + 2) = A(i, k + 2) - p * r
End If
A(i, k + 1) = A(i, k + 1) - p * q
A(i, k) = A(i, k) - p
Next i
End If
Next k
End If
End If
Loop While (L < nn - 1)
Wend

End Sub

Function SIGN(A, B)
If B > 0# Then
SIGN = Abs(A)
Else
SIGN = -Abs(A)
End If
End Function

Sub swap(A, B)
Dim temp

temp = A
A = B
B = temp
End Sub