Advertisement  

Thursday, 27 July 2017
     
 
Main Menu
Home Home
Shop Shop
News News
BASCOM-AVR BASCOM-AVR
BASCOM-8051 BASCOM-8051
Products Products
Application Notes Application Notes
Publications Publications
Links Links
Support Center Support Center
Downloads Downloads
Forum Forum
Resellers Resellers
Contact Us Contact Us
Updates Updates
MCS Wiki MCS Wiki
Online Help
BASCOM-AVR Help BASCOM-AVR Help
BASCOM-8051 Help BASCOM-8051 Help
Contents in Cart
Show Cart
Your Cart is currently empty.
Search the Shop

Products Search

User Login
Username

Password

If you have problem after log in with disappeared login data, please press F5 in your browser

RSS News
 
     
 

 
   
     
 
AN #176 - Mini Matrix Algebra Print
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