|
|
|
|
|
AN #176 - Mini Matrix Algebra |
|
by Natalius Kiedro
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
'@@@@@@@@@@@@@@@@@ @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ @@@@@@@@@@@@@@@@
'@@@@@@@@@@@@@@@@ @@@@@@@@@@@@@@@
'@@@@@@@@@@@@@@@ Mini Matrix Algebra Code for Bascom @@@@@@@@@@@@@@
'@@@@@@@@@@@@@@ Addition, Subtraction, Multiplication, Inversion @@@@@@@@@@@@@
'@@@@@@@@@@@@@@@ (Natalius Kiedro, March 2010) @@@@@@@@@@@@@@
'@@@@@@@@@@@@@@@@ @@@@@@@@@@@@@@@
'@@@@@@@@@@@@@@@@@ @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ @@@@@@@@@@@@@@@@
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@NK
'I missed matrix algebra because I recently went into a problem where I needed
'matrix multiplication in the AVR series (See my previous post on this issue).
'20 years ago I had matrix fun by writing a program which allowed to draw and to
'rotate 3D-pictures of biomolecules on a PC-screen. Another one from these days
'was on dynamic simulation and nonlinear data fitting on kinetic data. Things like
'3D-rotation and parameter optimization heavily depend on matrix algebra. Here is
'a little primer what it is, more information on thousands of pages in the inter-
'net.
'
'A matrix primer
'===============
'A matrix is a table of numbers. In most practical cases the numbers are organized
'in a square-like fashion such as this:
'
' | A(1,1) A(1,2) A(1,3) | | 1 2 3 |
' Matrix A = | A(2,1) A(2,2) A(2,3) | = | 4 5 6 | is a 3 x 3 matrix
' | A(3,1) A(3,2) A(3,3) | | 7 8 9 |
'
'The numbers A(i,j) in the matrix are called elements of a matrix. Element A(1,1)
'has the value 1, element A(2,3) the value 6 and so on. The variables i and j are
'indices of a matrix. i points to the row, and j points to the column where the
'element is to be found.
'
'Now, let's consider normal algebra. We have four basic operations everybody likes
'to work with:
'
'Addition: C = A + B
'Subtraction: C = A - B
'Multiplication: C = A * B
'Division: C = A / B
'
'Their equivalents in matrix algebra are called matrix addition, matrix subtraction
'matrix multiplication, and -OOPS- there is no term matrix division. In matrix
'algebra division is done stepwise. Let us do it in normal algebra to show how it
'is done in matrix algebra:
'
'C = 1 / B: In matrix algebra we call this matrix inversion.
'C = C * A: and this is matrix multiplication.
' Thus C = 1 / B * A = A / B
'
'So, matrix division is not an elementary operation here, because it can be broken
'down into matrix inversion and matrix multiplication.
'
'In normal algebra, 1, number one is special because:
'
'C = C * 1: (of course, Natalius, we know this!)
'C = C / 1: (come on, Natalius, do you want to bore us Bascomer?)
'
'No, I don't want to do bore anybody. I just want to share matrix fun with you!
'Like Neo did it some years ago. Let's reload the matrix which is the equivalent
'of number one:
'
' | 1 0 0 |
' | 0 1 0 | It is called identity matrix, others call it
' | 0 0 1 | unit matrix - and some like myself unity matrix.
'
'The path 1 1 1 through the matrix is coined matrix diagonal. The elements 0 are
'called off-diagonal elements. Special about a unity matrix is that all diagonal
'elements are one while all off-diagonal elements are zero.
'
'Why is number one and the unity matrix so special?
'
'Imagine, that multiplication and division were not available in Bascom and that
'you had to code these operations in assembly language yourself - starting from
'scratch. What would you do to proof that your algorithm is workin fine?
'Well, the easiest way were to multiply your input number by one or to devide it
'by one and to see that your number doesn't change - as expected. So, if we apply
'this recipy to matrix algebra, we expect that:
'
' Matrix A * Unity Matrix = Matrix A
'
'if the algorithm for matrix multiplication is working fine. If:
'
' Matrix A * Inverted Matrix A = Unity Matrix
'
'comes out, we know that our matrix inversion algorithm works fine. Simple, right?
'
'Matrix algebra in Bascom
'========================
'Matrix algebra is usually carried out using 2D-arrays which Bascom doesn't know.
'And it doesn't need to know it because assembly language sees computer memory
'always as a linear (1D) arrangement. Let us remember the arrangement of a 2D
'matrix by looking to the example in the beginning of the primer. In Bascom we
'arange the 2D-Matrix as a 1D-Matrix. So,Matrix A becomes:
'
'Matrix 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)|
'
'To convert a matrix algebra operation written for 2D matrices into "bascomized"
'1D matrices, one just needs to convert the indices i,j of elements A(i,j) into
'a single index ij pointing to the ij'th element A(ij) {A(1)..A(9)}.
'
'This is very simple: Here is an "empty" Matrix loop that does just this:
'
' Dim I As Byte, I0 As Byte, I As Byte, Ij As Byte, N As Byte
' ' N is the size of the matrix. 3 for a 3x3-, N for a NxN-matrix
' 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 'e.g for I = 2, I0 = 3:
' Ij = I0 + J '4, 5, 6
' Next J
' Next I
'
'Keep this "intuitive" variable notation for Ij or Ji or Ik or whatsover, because
'if N were 10, I were 6, and J were 3 then I0 were 50. It is easy to see that
'Ij = I0 + j because 53 = 50 + 3.
'
'The example and SUBs following are to play around in the simulator. I didn't
'config Serialout to let you see the output directly.
'
' Have matrix fun!
' Natalius
$regfile = "m128def.dat"
$crystal = 16000000
$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
'Filling matrix A with (pseudo)random numbers
'
For Iby = 1 To Mbyc
Iin = Rnd(1000)
Isi = Iin
A(iby) = Isi / 100
Next Iby
Print "Input matrix A"
Matprint A(1)
Matinv
Print "Inverted matrix B"
Matprint B(1)
Matmul
Print "Multiply matrix A with its inverse"
Matprint C(1)
Print "is the unity matrix if algo is correct"
'Note that limited precision of floating point operations will never lead to
'real zero when reestablishing the identity matrix from a matrix and its inverse.
'Zero shows up as either tiny positive or negative numbers. This is why -0.00
'and 0.00 are printed.
End 'end program
'*******************************************************************************
'******************************* Matrix Addition *******************************
'*******************************************************************************
Sub Matadd 'C = A + B
For Iby = 1 To Mbyc
C(iby) = A(iby) + B(iby)
Next Iby
End Sub
'*******************************************************************************
'******************************* Matrix Subtraction ****************************
'*******************************************************************************
Sub Matsub 'C = A - B
For Iby = 1 To Mbyc
C(iby) = A(iby) - B(iby)
Next Iby
End Sub
'*******************************************************************************
'*************************** Matrix Multiplication *****************************
'*******************************************************************************
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
'*******************************************************************************
'******************************* Matrix Inversion ******************************
'*******************************************************************************
Sub Matinv 'B = 1/A
'
'Copying the original matrix A into M1
'
For Iby = 1 To Mbyc
M1(iby) = A(iby)
Next Iby
'
'Filling two identity matrices: for i = j: M2(i,j) = 1 , M4(i,j) = 1
' for 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 'Rows from diagonal
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 'find largest element
Isi = Jsi '
Sby = Iby 'Memorize for swapping rows
End If
Next Iby
S0by = Sby - 1
S0by = S0by * Nbyc
For Iby = 1 To Nbyc
Riby = R0by + Iby 'Row R , column 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 "Singular matrix: No inversion possible"
End Sub
End If
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
'-----------------------------' Matrix multiplication 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 '---------------------------' Overwrite matrices M1,M2 with M3,M5
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 ' Reestablishing Identity Matrix M4
M4(ijby) = 1
Else
M4(ijby) = 0
End If
Next Jby
Next Iby
Next Rby
'-------------------------------' Matrix multiplication 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
'*******************************************************************************
'******************************* MatPrint **************************************
'*******************************************************************************
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
|
|
|
|
|
|
|