' Bascom的迷你矩陣代數代碼
' 加法、減法、乘法、倒置
' (Natalius Kiedro，2010年3月）
'我錯過了矩陣代數的學習，因為最近我在使用AVR系列（參見我之前關于這個問題的帖子）時遇到了一個需要用到‘矩陣乘法’的問題。
'20年前，我通過編寫一個程序來玩轉矩陣，這個程序可以在電腦屏幕上繪制并旋轉生物分子的3D圖像。那個
'時期的另一個項目是關于動力學數據的動態模擬和非線性數據擬合。像3D旋轉和參數優化這類事情，在很大
'程度上都依賴于矩陣代數。下面是一個關于矩陣的簡要介紹，互聯網上有成千上萬頁的資料可以提供更多信息。
'矩陣入門簡介’
'===============
'矩陣是一個數字表。在大多數情況下，數字以類似正方形的形式排列，如下所示：
' | A(1,1) A(1,2) A(1,3) | | 1 2 3 |
'矩陣A = | A(2,1) A(2,2) A(2,3) | = | 4 5 6 | 是一個3x3矩陣
' | A(3,1) A(3,2) A(3,3) | | 7 8 9 |
'矩陣中的數字A(i,j)被稱為矩陣的元素。元素A(1,1)的值為1，元素A(2,3)的值為6，以此類推。
'變量i和j是矩陣的索引，i指向行，j指向列，元素就位于該行和該列的交叉點上。
'現在，讓我們來看看普通代數。我們有四個每個人都喜歡使用的基本運算：
'加法：C = A + B
'減法：C = A - B
'乘法：C = A * B
'除法：C = A / B
'它們在矩陣代數中的等價運算分別被稱為矩陣加法、矩陣減法、矩陣乘法，以及——哎呀——沒有所謂的矩陣除法。
'在矩陣代數中，除法是分步進行的。讓我們用普通代數來演示一下，看看它在矩陣代數中是如何進行的：
'C = 1 / B：在矩陣代數中，我們稱之為矩陣求逆。
'C = C * A：這是矩陣乘法。
'因此，C = 1 / B * A = A / B（在矩陣代數中，這種表示并不嚴謹，僅用于說明概念）
'所以，在矩陣代數中，除法并不是一個基本運算，因為它可以分解為矩陣求逆和矩陣乘法。
'在普通代數中，數字1是特殊的，因為：
'C = C * 1：（當然，納塔利烏斯，我們都知道這一點！）
'C = C / 1：（拜托，納塔利烏斯，你想讓我們覺得無聊嗎，巴斯科默？）
'不，我不想讓任何人覺得無聊。我只是想和你分享矩陣的樂趣！就像尼奧幾年前做的那樣。
'讓我們重新加載這個矩陣，它是……的等價物（此處原文被截斷，但根據上下文，可以推測“
'它是……的等價物”可能是指矩陣運算在某種情況下的應用或等價表示）。
'數字1的特殊性在于：
'| 1 0 0 |
'| 0 1 0 | 它被稱為單位矩陣，也有人稱之為
'| 0 0 1 | 統一矩陣——當然，還有像我這樣的人喜歡稱之為單位矩陣。
'從矩陣的1,1,1路徑穿過的線被稱為矩陣的對角線。元素0被稱為非對角元素。單位矩陣的特別之處在于，
'其對角線上的元素全部為1，而非對角線上的元素全部為0。
'為什么數字1和單位矩陣如此特殊呢？
'想象一下，如果在Bascom中無法使用乘法和除法，而你需要從零開始，在匯編語言中自己編寫這些運算的代碼。
'你會如何證明你的算法是有效的呢？
'最簡單的方法就是將你的輸入數字乘以1或除以1，然后觀察數字是否保持不變——正如預期的那樣。
'因此，如果我們將這個原理應用到矩陣代數中，我們期望：
'矩陣A乘以單位矩陣 = 矩陣A
'如果矩陣乘法的算法是有效的。如果：
'矩陣A乘以矩陣A的逆矩陣 = 單位矩陣
'那么我們就知道我們的矩陣求逆算法是有效的。很簡單，對吧？
'Bascom中的矩陣代數
'========================
'矩陣代數通常使用二維數組進行計算，但Bascom并不了解二維數組。而且它也不需要了解，
'因為匯編語言總是將計算機內存視為線性（一維）排列的。讓我們回顧一下本教程開頭給出的二維矩陣排列的例子。
'在Bascom中，我們將二維矩陣當作一維矩陣來處理。因此，矩陣A變為：
'矩陣A = |A(1,1) A(1,2) A(1,3) A(2,1) A(2,2) A(2,3) A(3,1) A(3,2) A(3,3)|
'= | 1 2 3 4 5 6 7 8 9 |
'= | A(1) A(2) A(3) A(4) A(5) A(6) A(7) A(8) A(9)|
'為了將針對二維矩陣編寫的矩陣代數運算轉換為適用于“Bascom化”的一維矩陣的運算，
'我們只需將元素A(i,j)的索引i,j轉換為一個指向第ij個元素A(ij)（即A(1)至A(9)）的單索引ij。
'這非常簡單：以下是一個只執行此操作的“空”矩陣循環：
'Dim I As Byte, I0 As Byte, J As Byte, Ij As Byte, N As Byte
' N是矩陣的大小。對于3x3矩陣，N為3；對于NxN矩陣，N為相應的值
'N = 3
'For I = 1 TO N '1, 2, 3
'I0 = I - 1 '0, 1, 2
'I0 = I0 * N '0, 3, 6（計算當前行之前所有行的元素總數）
'For J = 1 TO N '例如，當I=2時，I0=3，則J=1,2,3
'Ij = I0 + J '計算一維數組中的索引位置，即4, 5, 6
'Next J
'Next I
'斷線內為AI給出的釋義：
'////////////////////////////////////////////////////////////////////////////////////////////////////////////
'在矩陣程序中，變量通常用于表示矩陣的維度、索引位置以及輔助計算。以下是對您提供的代碼中各個變量的具體含義的解釋：
' N As Byte：
'含義：N代表矩陣的大小，即矩陣的行數和列數。對于NxN矩陣，N的值是相同的。
'示例：對于3x3矩陣，N為3。
'I As Byte：
'含義：I用作矩陣行的索引，通常從1開始，到N結束。
'示例：對于3x3矩陣，I的取值為1, 2, 3。
'I0 As Byte：
'含義：I0是一個輔助變量，用于計算當前行之前所有行的元素總數。它首先通過I - 1得到當前行的相對位置（從0開始計數），
'然后乘以N（矩陣的列數/大小）來得到之前所有行的元素總數。
'示例：對于3x3矩陣，當I=1時，I0=0；當I=2時，I0=3（因為第一行有3個元素）；當I=3時，I0=6（因為前兩行共有6個元素）。
'J As Byte：
'含義：J用作矩陣列的索引，通常從1開始，到N結束。
'示例：對于3x3矩陣，當I=2（即處理第二行）時，J的取值為1, 2, 3。
'Ij As Byte：
'含義：Ij是一個輔助變量，用于計算當前元素在一維數組中的索引位置。它通過I0 + J得到，
'其中I0表示當前元素之前所有元素的數量，J表示當前元素在當前行中的位置。
'示例：對于3x3矩陣，當I=2（第二行）且J=1時，Ij=I0+J=3+1=4，表示第二行第一個元素在一維數組中的索引為4。
'這段代碼的目的是遍歷一個NxN矩陣的所有元素，并通過計算每個元素在一維數組中的索引位置，以便
'可以對這些元素進行操作（如讀取、修改等）。這種方法在處理需要將二維數組（如矩陣）轉換為一維數組進行處理的場景時非常有用。
'////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
'請保持這種“直觀”的變量表示法，如Ij、Ji、Ik等，因為當N=10、I=6、J=3時，I0=50。很容易看出Ij=I0+j，因為53=50+3。
'下面的示例和子程序可以在模擬器中運行。我沒有配置串行輸出，因此您無法直接看到輸出結果。
'祝您在矩陣運算中玩得開心！
'作者：納塔利烏斯
$regfile = "m328pdef.dat"
$crystal = 8000000
$baud = 19200
Const Nbyc = 4
Const Mbyc = 16 'Nbyc * Nbyc
Declare Sub Matadd()
Declare Sub Matsub()
Declare Sub Matmul()
Declare Sub Matinv()
Declare Sub Matprint(byval Ksi As Single)
Dim Iby As Byte , I0by As Byte , Ijby As Byte , Ikby As Byte , Irby As Byte
Dim Jby As Byte , J0by As Byte
Dim Kby As Byte , K0by As Byte , Kjby As Byte
Dim Rby As Byte , R0by As Byte , Riby As Byte , Rjby As Byte , Rrby As Byte , Rsby As Byte
Dim Sby As Byte , S0by As Byte , Siby As Byte , Srby As Byte , Ssby As Byte
Dim Iin As Integer
Dim Isi As Single , Jsi As Single , Ksi As Single , Msi As Single
Dim A(mbyc) As Single , B(mbyc) As Single , C(mbyc) As Single
'Temporary:
Dim M1(mbyc) As Single , M2(mbyc) As Single , M3(mbyc) As Single , M4(mbyc) As Single , M5(mbyc ) As Single
'用（偽）隨機數填充矩陣A
'
For Iby = 1 To Mbyc
Iin = Rnd(1000)
Isi = Iin
A(iby) = Isi / 100
Next Iby
Print "輸入矩陣A"
Matprint A(1)
Matinv
Print "倒置矩陣B"
Matprint B(1)
Matmul
Print "將矩陣A與其逆相乘"
Matprint C(1)
Print "如果算法是正確，那么它就是單位矩陣。"
'請注意，當從矩陣及其逆矩陣重新建立單位矩陣時，浮點運算的有限精度永遠不會導致真正的零出現。
'零會表現為極小的正數或負數。這就是為什么打印出的是-0.00和0.00的原因。
End 'end program
'*******************************************************************************
'******************************* 矩陣加法 *******************************
'*******************************************************************************
Sub Matadd 'C = A + B
For Iby = 1 To Mbyc
C(iby) = A(iby) + B(iby)
Next Iby
End Sub
'*******************************************************************************
'******************************* 矩陣減法 ****************************
'*******************************************************************************
Sub Matsub 'C = A - B
For Iby = 1 To Mbyc
C(iby) = A(iby) - B(iby)
Next Iby
End Sub
'*******************************************************************************
'*************************** 矩陣乘法*****************************
'*******************************************************************************
Sub Matmul 'C = A * B
For Iby = 1 To Nbyc
I0by = Iby - 1
I0by = I0by * Nbyc
For Jby = 1 To Nbyc
Ijby = I0by + Jby
C(ijby) = 0
For Kby = 1 To Nbyc
K0by = Kby - 1
K0by = K0by * Nbyc
Kjby = K0by + Jby
Ikby = I0by + Kby
Jsi = A(ikby) * B(kjby)
C(ijby) = C(ijby) + Jsi
Next Kby
Next Jby
Next Iby
End Sub
'*******************************************************************************
'******************************* 矩陣求逆 ******************************
'*******************************************************************************
Sub Matinv 'B = 1/A
'
'將原始矩陣A復制到M1中
'
For Iby = 1 To Mbyc
M1(iby) = A(iby)
Next Iby
'
'填充兩個單位矩陣：對于i等于j的情況，M2(i,j) = 1，M4(i,j) = 1；
'對于i不等于j的情況，M2(i,j) = 0，M4(i,j) = 0。
For Iby = 1 To Nbyc
I0by = Iby - 1
I0by = I0by * Nbyc
For Jby = 1 To Nbyc
Ijby = I0by + Jby
If Iby = Jby Then
M2(ijby) = 1
M4(ijby) = 1
Else
M2(ijby) = 0
M4(ijby) = 0
End If
Next Jby
Next Iby
For Rby = 1 To Nbyc 'Row byte RBy
R0by = Rby - 1 '0, 1, 2, 3..
R0by = R0by * Nbyc '0,Nbyc,2Nbyc,3Nbyc
Isi = 0
'------------------------------------'Pivotation中心點或軸
For Iby = Rby To Nbyc '矩陣或表格的對角線開始提取或構建行
I0by = Iby - 1
I0by = I0by * Nbyc
Irby = I0by + Rby 'index of element in 1D: Ir一維中元素的索引
Jsi = M1(irby)
Jsi = Abs(jsi)
If Jsi >= Isi Then '找到最大的元素
Isi = Jsi '
Sby = Iby '記住用于交換行的（操作或步驟）
End If
Next Iby
S0by = Sby - 1
S0by = S0by * Nbyc
For Iby = 1 To Nbyc
Riby = R0by + Iby 'Row R , column I “R”代表行號，而“I”代表列號。
Siby = S0by + Iby 'Row to swap, “要交換的行”
Swap M1(riby) , M1(siby)
Next Iby
Rrby = R0by + Rby
Jsi = M1(rrby)
Jsi = Abs(jsi)
If Jsi < 10e-30 Then
Print "奇異矩陣：無法進行求逆"
End If
'End Sub
Ssby = S0by + Sby
Rsby = R0by + Sby
Srby = S0by + Rby
M4(rrby) = 0
M4(ssby) = 0
M4(rsby) = 1
M4(srby) = 1
'-----------------------------------' 高斯-約旦消元法（Gauss-Jordan Elimination）
For Iby = 1 To Nbyc
I0by = Iby - 1
I0by = I0by * Nbyc
Irby = I0by + Rby
For Jby = 1 To Nbyc
J0by = Iby - 1
J0by = J0by * Nbyc
Ijby = I0by + Jby
Rjby = R0by + Jby
If Iby = Rby Then
If Iby = Jby Then
M3(rrby) = 1 / M1(rrby)
Else
M3(ijby) = M1(ijby) / M1(rrby)
M3(ijby) = 0 - M3(ijby)
End If
Elseif Jby = Rby Then
M3(irby) = M1(irby) / M1(rrby)
Else
Jsi = M1(rjby) * M1(irby)
Jsi = Jsi / M1(rrby)
M3(ijby) = M1(ijby) - Jsi
End If
Next Jby
Next Iby
'-----------------------------' 矩陣乘法 M5 = M4 × M2
For Iby = 1 To Nbyc
I0by = Iby - 1
I0by = I0by * Nbyc
For Jby = 1 To Nbyc
Ijby = I0by + Jby
M5(ijby) = 0
For Kby = 1 To Nbyc
K0by = Kby - 1
K0by = K0by * Nbyc
Kjby = K0by + Jby
Ikby = I0by + Kby
Jsi = M4(ikby) * M2(kjby)
M5(ijby) = M5(ijby) + Jsi
Next Kby
Next Jby
Next Iby
'---------------------------' 用M3,M5覆蓋矩陣M1,M2
For Iby = 1 To Nbyc
I0by = Iby - 1
I0by = I0by * Nbyc
For Jby = 1 To Nbyc
Ijby = I0by + Jby
M2(ijby) = M5(ijby)
M1(ijby) = M3(ijby)
If Iby = Jby Then ' 重建單位矩陣 M4
M4(ijby) = 1
Else
M4(ijby) = 0
End If
Next Jby
Next Iby
Next Rby
'-------------------------------' 矩陣乘法 B = M1×M2
For Iby = 1 To Nbyc
I0by = Iby - 1
I0by = I0by * Nbyc
For Jby = 1 To Nbyc
J0by = Jby - 1
J0by = J0by * Nbyc
Ijby = I0by + Jby
B(ijby) = 0
For Kby = 1 To Nbyc
K0by = Kby - 1
K0by = K0by * Nbyc
Ikby = I0by + Kby
Kjby = K0by + Jby
Jsi = M1(ikby) * M2(kjby)
B(ijby) = B(ijby) + Jsi
Next Kby
Next Jby
Next Iby
End Sub
'*******************************************************************************
'******************************* 打印矩陣 **************************************
'*******************************************************************************
Sub Matprint(byval Ksi As Single)
Jby = 0 : Kby = 0
If Ksi = A(1) Then
Kby = 1
Elseif Ksi = B(1) Then
Kby = 2
Elseif Ksi = C(1) Then
Kby = 3
End If
For Iby = 1 To Mbyc
Select Case Kby
Case 1
Msi = A(iby)
Case 2
Msi = B(iby)
Case 3
Msi = C(iby)
Case Else
Exit Sub
End Select
If Msi >= 0 Then
Print " ";
End If
Print Fusing(msi , "#.##") ; " ";
Incr Jby
If Jby = Nbyc Then
Print
Jby = 0
End If
Next Iby
End Sub
Proteus 仿真輸出結果：
|
