Author |
Topic: Matrix multiplication (dot product) (Read 1079 times) |
|
Richard Russell
Administrator
member is offline
 
Posts: 1348
|
 |
Matrix multiplication (dot product)
« Thread started on: Nov 11th, 2016, 11:47am » |
|
Neither LB 4 nor LBB supports matrix multiplication (dot product) as a native operator, but BBC BASIC does so that makes it easy to incorporate it in an LBB program. This is exactly what I did as one of my speed-up modifications to Bluatigro's CG animation demo.
His original code used a multiply routine similar to this:
Code:' matrix( c ) = matrix( a ) . matrix( b )
for i = 0 to 3
for j = 0 to 3
c(i,j) = 0
for k = 0 to 3
c(i,j) = c(i,j) + a(i,k) * b(k,j)
next k
next j
next i
In my LBB version I used this as a (very) much faster replacement:
Code:
You might be surprised to discover how often matrix multiplication is useful. 3D manipulations (rotation, translation, scaling) as used by Bluatigro are probably some of the most common applications, but there are many others.
Richard.
|
|
 Logged
|
|
|
|
CryptoMan
New Member
member is offline
Gender: 
Posts: 46
|
|
Re: Matrix multiplication (dot product)
« Reply #1 on: Nov 18th, 2016, 10:43pm » |
|
How about Inverse of a square matrıx?
|
|
Logged
|
|
|
|
Richard Russell
Administrator
member is offline
Posts: 1348
|
|
Re: Matrix multiplication (dot product)
« Reply #2 on: Nov 19th, 2016, 09:59am » |
|
on Nov 18th, 2016, 10:43pm, CryptoMan wrote:| How about Inverse of a square matrıx? |
|
There's no built-in matrix inversion, but it's not that difficult to code in BASIC (the complexity depends on whether you need to choose the pivots carefully to maximize accuracy and avoid division-by-zero errors). Of course it won't be as fast as a native function would be.
I've got a matrix inversion routine in BBC BASIC but not in Liberty BASIC. Rather surprisingly, Rosetta Code doesn't seem to have a matrix inversion task either.
If I were to create a Programming Challenge to write a fast, accurate and stable matrix inversion in LB, do you think anybody other than me would submit an entry? 
Richard.
|
|
Logged
|
|
|
|
CryptoMan
New Member
member is offline
Gender:
Posts: 46
|
|
Re: Matrix multiplication (dot product)
« Reply #3 on: Nov 19th, 2016, 10:41pm » |
|
I would. In my past iife I used to work a lot on that. In those days memory was limited and we were spending more effort in storşng the matrices in smaller spacesdue to the nature of our matrices which were symetric around the diagonal. Today, there is obscene amount of memory is available.
Ofcourse, I know it can be done in Basic but wondering if there is a native version maybe even by accessing libraries on graphics processor cards which can do it lightning fast.
Probably, there are some 3rd party DLLs available for full set of matrix arithmetic. Just asked as a matter of curiosity if more native matrix functions
Available. I don't even remember the theory of Finite Elements now to make any use of it. Need to open the texr books and remember from 35 years ago.
I used to write it on TI59 and HP41CV programmable calculators and solving some of the problems in exams very quickly during the days when they did not expect such things from calculators and not forbidding them from exams. Later when IBM PC came in 1981 it was so simple and huge so we translated lots of code from IBM mainframe FORTRAN to BASIC, Pascal and C. So, it is always interesting for me to see what can be now done with all this additional power we have today. Some amazing things should be possible.
|
|
Logged
|
|
|
|
Richard Russell
Administrator
member is offline
Posts: 1348
|
|
Re: Matrix multiplication (dot product)
« Reply #4 on: Nov 20th, 2016, 11:48am » |
|
on Nov 19th, 2016, 10:41pm, CryptoMan wrote:| the nature of our matrices which were symetric around the diagonal. |
|
Mine usually are too, I suspect it follows when they are used for least-squares regression. Inverting such a symmetrical matrix is less demanding, I think, and that supports my feeling that there isn't a 'one size fits all' inversion routine. That may be one reason why inversion wasn't incorporated as a native operation in BBC BASIC 30 years ago (I was there at the time, but I don't think the BBC expressed a view either way).
Nevertheless embedded BBC BASIC code for array arithmetic can be very useful; one neat thing you can do is to normalise a vector to unit length in one fast operation:
Code:
That can be valuable in some 3D graphics applications.
Quote:| I know it can be done in Basic but wondering if there is a native version maybe even by accessing libraries |
|
In my experience matrix inversion hasn't been time-critical. Even in a real-time application like free-d2 (see the BBC BASIC Birthday page) it was acceptable to do the inversion in interpreted BASIC.
You may be interested in this example of using array arithmetic to fit a polynomial function to a large set of data points.
Richard.
|
|
Logged
|
|
|
|
Optimax
New Member
member is offline
Posts: 11
|
|
Re: Matrix multiplication (dot product)
« Reply #5 on: Nov 21st, 2016, 07:07am » |
|
A long time ago, I wrote an "inverse matrix" in Rapidq, that I later ported to LB. It uses the method of "minors" elaborated by Gauss and Jordan. There is an other method called "the simplex" that I never could find back.
Strangely there appears at the end a "-0.00" which is anyway equal to zero; could maybe come from rounding ?
Code:'GAUSS-JORDAN MATRICE INVERSE
'=============================
'A original matrix, InvA inverse of A, B extended matrix, P = A * InvA
'let's have a square matrix A (3x3) '2, 1, -2
'3 2 2
'5 4 3
'--- creating the extended matrix B = A + unit matrix
DATA 2, 1, -2, 1, 0, 0
DATA 3, 2, 2, 0, 1, 0
DATA 5, 4, 3, 0, 0, 1
n = 3
DIM A(n+1, n+1)
DIM B(n+1, n*2)
PRINT "EXTENDED MATRIX B :": PRINT
FOR i = 1 TO n
FOR j = 1 TO 2*n
READ m
B(i,j) = m
PRINT "", B(i,j); 'and displaying it
NEXT j
PRINT: PRINT
NEXT i
PRINT
'--- A is the left part of B, extraction
PRINT "ORIGINAL MATRIX A, extracted from B :"
PRINT
FOR i = 1 TO n
FOR j = 1 TO n
A(i,j) = B(i,j)
PRINT "", A(i,j); 'and displaying it
NEXT j
PRINT: PRINT
NEXT i
PRINT
'--- calculation begins here
'i = row index 'k = other row index (not-i) 'j = column index
'intermediary variables p, q, r to avoid divisions like A(i,j)/A(m,n) e.g., giving unpredictable results
FOR i = 1 TO n
'reducing each diagonal term B(i,i) to 1 by dividing the hole row by B(i,i)
p = B(i, i)
FOR j = 1 TO 2*n
q = B(i, j)
B(i, j) = q / p
NEXT j
'multiply row i by B(k,i)
'substract row i from row k
FOR k = 1 TO n
p = B(k, i)
IF k <> i THEN
FOR j = 1 TO 2*n
q = B(k, j)
r = B(i, j)
B(k, j) = q - p * r
NEXT j
END IF
NEXT k
NEXT i
'--- end of calculation, let's see matrix B now
PRINT "INTERMEDIAIRY MATRIX B :": PRINT
FOR i = 1 TO n
FOR j= 1 TO 2*n
m = B(i, j)
PRINT "", USING("###.###", m);
NEXT j
PRINT: PRINT
NEXT i
PRINT
'---INVERSE of a, i.e. InvA, is the *right* part of B
PRINT "INVERSE of A (rounded) :": PRINT
DIM InvA(n+1,n+1)
FOR i = 1 TO n
FOR j = n+1 TO 2*n 'copying the right part of B
InvA(i, j-n) = B(i,j)
m = InvA(i, j-n)
PRINT "", USING("###.###", m); 'rounding
NEXT j
PRINT: PRINT
NEXT i
PRINT
'--- Checking the result: A * InvA must guve P = the unit matrix
PRINT "CHECKING for InvA: A * Inv must give the unit matrix P"
PRINT
DIM P(n+1, n+1) 'empty
FOR i = 1 TO n 'for each row i o f A and P
FOR j = 1 TO n 'for each horizontal term j of row i of P
FOR k = 1 TO n 'for each vertical term k of column j of InvA
x = A(k,i): y = InvA(j,k)
P(i, j) = P(i, j) + x*y
NEXT k
m = P(i, j)
PRINT "", USING("###.###", m);
NEXT j
PRINT: PRINT
NEXT i
END
|
|
Logged
|
|
|
|
Richard Russell
Administrator
member is offline
Posts: 1348
|
|
Re: Matrix multiplication (dot product)
« Reply #6 on: Nov 21st, 2016, 09:43am » |
|
on Nov 21st, 2016, 07:07am, Optimax wrote:| A long time ago, I wrote an "inverse matrix" in Rapidq, that I later ported to LB. It uses the method of "minors" elaborated by Gauss and Jordan. |
|
The problem with that method is that there's no 'pivot selection': you simply work down the diagonal from top-left to bottom-right. This can result in inaccuracy, or worse, if you are unlucky with your input matrix.
That can easily be demonstrated by changing your test matrix so that the top-left element is 0 rather than 2:
Code:DATA 0, 1, -2, 1, 0, 0
DATA 3, 2, 2, 0, 1, 0
DATA 5, 4, 3, 0, 0, 1
Now when I run your program it reports a Division by zero error! It's not that the matrix is non-invertable, in fact the inverse matrix is this:
Code: -0.4 1 -0.4
0.2 -2 1.2
0.4 1 -0.6
The BBC BASIC code I referred to does attempt to be more 'intelligent' in its pivot selection, and it is successful in this case, but I make no claims for it in other respects.
Quote:| Strangely there appears at the end a "-0.00" |
|
I don't see that here, but it's a perfectly reasonable thing to output (representing a negative number which is zero when rounded to two decimal places)
Richard.
|
|
Logged
|
|
|
|
Optimax
New Member
member is offline
Posts: 11
|
|
Re: Matrix multiplication (dot product)
« Reply #7 on: Nov 25th, 2016, 12:30pm » |
|
You are right.
Thanks for advice.
I should have a look at that pivot thing.
|
|
Logged
|
|
|
|
|
