'DOAN CHUONG TRINH TINH TOAN CAC GIA TRI RIENG CUA MTRAN
'HOAN THANH NGAY 05/03/2003. Trong ham nay khong xet den phan ao cua cac gia tri rieng
Private Sub Balanc(A() As Double, ByVal n As Integer)
Const RADIX = 2
Dim last As Integer, i As Integer, j As Integer
Dim s As Double, r As Double, g As Double
Dim f As Double, c As Double, sqrdx As Double
DoEvents: If mdl01.tiep_tuc = False Then Exit Sub
frmSTATUS.lblstatus01(1).Caption = mdl01.StrTbao & Chr(10) & "(Transform to Balanc matric)"
frmSTATUS.ProgressBar1.max = n
frmSTATUS.ProgressBar1.Value = 0
frmSTATUS.lblstatus01(1).Refresh
sqrdx = RADIX * RADIX
last = 0
While last = 0
frmSTATUS.ProgressBar1.Value = frmSTATUS.ProgressBar1.Value + 1
last = 1
For i = 1 To n
r = 0
c = 0
For j = 1 To n
If Not (i = j) Then
c = c + Abs(A(j, i))
r = r + Abs(A(i, j))
End If
Next
If (c And r) Then
g = r / RADIX
f = 1
s = c + r
Do While c < g
f = f * RADIX
c = c * sqrdx
Loop
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
For j = 1 To n
A(j, i) = A(j, i) * f
Next
End If
End If
Next
Wend
frmSTATUS.ProgressBar1.Value = frmSTATUS.ProgressBar1.max
End Sub
Private Sub Hessenberg(A() As Double, ByVal n As Integer)
Dim m As Integer, i As Integer, j As Integer
Dim x As Double, y As Double, tg As Double
frmSTATUS.lblstatus01(1).Caption = mdl01.StrTbao & Chr(10) & "(Transforming to Hessenberg matric form)"
frmSTATUS.ProgressBar1.max = n - 3
frmSTATUS.ProgressBar1.Value = 0
frmSTATUS.lblstatus01(1).Refresh
For m = 2 To n - 1
DoEvents
If mdl01.tiep_tuc = False Then Exit Sub
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
If Not (i = m) Then
For j = m - 1 To n
tg = A(i, j)
A(i, j) = A(m, j)
A(m, j) = tg
Next
For j = 1 To n
tg = A(j, i)
A(j, i) = A(j, m)
A(j, m) = tg
Next
End If
If Not (x = 0) Then
For i = m + 1 To n
y = A(i, m - 1) / x
For j = m To n
A(i, j) = A(i, j) - y * A(m, j)
Next
For j = 1 To n
A(j, m) = A(j, m) + y * A(j, i)
Next
Next
End If
frmSTATUS.ProgressBar1.Value = m - 2
Next
End Sub
Private Function Sign(ByVal A As Double, ByVal B As Double) As Double
If B > 0 Then
Sign = Abs(A)
Else
Sign = -Abs(A)
End If
End Function
Public Function Eig(A() As Double, ByVal n As Integer) As Double()
Dim nn As Integer, m As Integer, l As Integer
Dim k As Integer, j As Integer, its As Integer
Dim i As Integer, mmin As Integer, status As Integer
Dim z As Double, y As Double, x As Double
Dim w As Double, v As Double, u As Double
Dim t As Double, s As Double, r As Double
Dim p As Double, q As Double, anorm As Double
Dim wr() As Double
'========================================
Balanc A, n
Hessenberg A, n
'========================================
'BAT DAU TINH CAC GIA TRI RIENG:
'========================================
ReDim wr(1 To n) As Double
frmSTATUS.lblstatus01(1).Caption = mdl01.StrTbao & Chr(10) & "(Computing EigValues)"
frmSTATUS.ProgressBar1.max = n
frmSTATUS.ProgressBar1.Value = 0
frmSTATUS.lblstatus01(1).Refresh
anorm = Abs(A(1, 1))
For i = 2 To n
For j = (i - 1) To n
anorm = anorm + Abs(A(i, j))
Next
Next
status = 0
nn = n
t = 0
While nn >= 1
DoEvents
If mdl01.tiep_tuc = False Then Exit Function
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
x = A(nn, nn)
If l = nn Then
wr(nn) = x + t
nn = nn - 1
status = status + 1
frmSTATUS.ProgressBar1.Value = status
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 - 1) = x + z
wr(nn) = x + z
If Not (z = 0) Then wr(nn) = x - w / z
Else
wr(nn - 1) = x + p
wr(nn) = x + p
End If
nn = nn - 2
status = status + 2
frmSTATUS.ProgressBar1.Value = status
Else
If (its = 10) Or (its = 20) Then
t = t + x
For i = 1 To nn
A(i, i) = A(i, i) - x
Next
s = Abs(A(nn, nn - 1)) + Abs(A(nn - 1, nn - 2))
y = 0.75 * s
x = 0.75 * s
w = -0.4375 * s * s
End If
its = its + 1
For m = (nn - 2) To l 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 = l 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
For i = m + 2 To nn
A(i, i - 2) = 0
If Not (i = (m + 2)) Then A(i, i - 3) = 0
Next
For k = m To nn - 1
If Not (k = m) Then
p = A(k, k - 1)
q = A(k + 1, k - 1)
r = 0
If Not (k = (nn - 1)) Then r = A(k + 2, k - 1)
x = Abs(p) + Abs(q) + Abs(r)
If Not (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 Not (s = 0) Then
If k = m Then
If Not (l = m) Then
A(k, k - 1) = -A(k, k - 1)
End If
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 Not (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
If nn < k + 3 Then
mmin = nn
Else
mmin = k + 3
End If
For i = l To mmin
p = x * A(i, k) + y * A(i, k + 1)
If Not (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
End If
Next
End If
End If
Loop While l < (nn - 1)
Wend
Eig = wr
End Function