HotBasic --> Timing Matrix Inversions

[danbaron]

[font=courier new][size=8pt]Here is a HotBasic program which times matrix inversions.

(Note that HotBasic is compiled. I posted this topic under Scripting, because, I could not determine a more appropriate category.)

The program fills an N x N matrix with random doubles in the range of (0, 1), and inverts it using the same
algorithm that was used for the corresponding C, Java, PowerBASIC, and PureBasic programs.

The program runs in a console window, and indicates its progress in the window.

It times the inversion, and then writes the elapsed time, and a check vector to the file, "matinv.txt". If every value
of the N-dimensional check vector is very close to 1, then the inversion is correct.

I ran the program three times for a 1000 x 1000 matrix.

The times on my machine were (seconds),

1 147.187
2 147.328
3 147.312

I attached the program file below. I called it, "matinv.bas.txt".

You can download the trial version of HotBasic from the HotBasic website.
The trial version will compile console programs.

If you download the attached file, then, remove the ".txt" extension.

:oops: :x :twisted:
Dan


$optimize
' Optimize directive (previous line), must be the first line of the file.
'***************************************************************************************************

' FILE = "matinv.bas"

'***************************************************************************************************

$apptype console
$typecheck on

$define msize 1000
$define FN "matinv.txt"

'***************************************************************************************************

defdbl t1
defdbl t2
defdbl tt

dim a0(1 to msize, 1 to msize) as double
dim ai(1 to msize, 1 to msize) as double
dim id(1 to msize, 1 to msize) as double
dim va(1 to msize) as double
dim vb(1 to msize) as double

defint a0p = objptr(a0)
defint aip = objptr(ai)
defint idp = objptr(id)
defint vap = objptr(va)
defint vbp = objptr(vb)

'***************************************************************************************************

declare sub initialize(m1 as integer, m2 as integer, v1 as integer)
declare sub consoleoutput(whichtime as integer, t as double)
declare sub writefile(v as integer, t as double)
declare sub inv(a as integer)
declare sub matmatmult(m1 as integer, m2 as integer, m3 as integer)
declare sub matvecmult(m1 as integer, v1 as integer, v2 as integer)

'***************************************************************************************************

defstr wait

' You don't pass arrays, you pass pointers to arrays.
initialize(a0p, aip, vap)
consoleoutput(1, tt)

t1 = timer
inv(aip)
t2 = timer
tt = t2 - t1

consoleoutput(2, tt)

matmatmult(a0p, aip, idp)
matvecmult(idp, vap, vbp)

writefile(vbp, tt)

consoleoutput(3, tt)

input wait

end


'***************************************************************************************************

sub initialize(m1 as integer, m2 as integer, v1 as integer)
' The next two lines do pointer arithmetic.
' They get information about the array, found at specific offset locations.
' "iadd", means, "integer addition".
defint first = byref(iadd(v1, 24))
defint last = byref(iadd(v1, 44))

defint i, j
defdbl d

randomize timer

for i = first to last
arrayref(v1, i) = 1
for j = first to last
d = rnd
arrayref(m1, i, j) = d
arrayref(m2, i, j) = d
next
next
end sub

'***************************************************************************************************

sub consoleoutput(whichtime as integer, t as double)
defstr so, sn, st

sn = str$(msize)
st = str$(t)

select case whichtime

case 1
so = "Inverting a " + sn + " by " + sn + " matrix.."
print
print(so)
print

case 2
print("Done.")
print
so = "elapsed time = " + st + " seconds"
print(so)
print
print("Calculating and writing the check vector..")
print

case 3
print("Done.")
print
print("Press 'enter' to quit.")
print

case else
print("Something wrong.")
print

end select

end sub

'***************************************************************************************************

sub writefile(v as integer, t as double)
defint first, last, i, j, m, n
defstr s1, s2, sn, so, st, d0
dim outfile as file
dim d as date

d.update
d0 = "' " + str$(d.year) + "-" + str$(d.month) + "-" + str$(d.day) + ", " + str$(d.hour) + ":" + str$(d.minute) + "."

first = byref(iadd(v1, 24))
last = byref(iadd(v1, 44))

sn = str$(msize)
st = str$(t)

outfile.open(FN, 65535)

s1 = "' file = " + "'" + FN + "'"
outfile.writeline(s1)
outfile.writeline(d0)
outfile.writeline("")
s1 = "' Inversion of a " + sn + " by " + sn + " matrix,"
outfile.writeline(s1)
outfile.writeline("' filled with uniformly distributed random doubles,")
outfile.writeline("' in the range of (0, 1).")
outfile.writeline("")
s1 = "' elapsed time = " + st + " seconds"
outfile.writeline(s1)
outfile.writeline("")
outfile.writeline("' If the inversion was calculated correctly, then every")
outfile.writeline("' value shown below should be very close to 1.")
outfile.writeline("")
outfile.writeline("index value")

for i = first to last

s1 = str$(i)
m = len(s1)
n = 4 - m
for j = 1 to n
s1 = insert$(s1, "0", 1)
next

s1 = s1 + " "
s2 = str$(arrayref(v, i))
s1 = s1 + s2
outfile.writeline(s1)
next

outfile.close

end sub

'***************************************************************************************************

sub inv(a as integer)
' Uses the Gauss-Jordan Method.
' Inverts matrix, a(), "in-place".

defint pivot, row, col, maxrow, rowswitch, switchcount, first, last, i, j
defdbl test, temp, colmax, divisor, factor

'---------------------------------------------

first = byref(iadd(m1, 24))
last = byref(iadd(m1, 44))

dim rowswitches(1 to msize, 1 to 2) as integer

switchcount = 0

for pivot = first to last

colmax = 0
maxrow = pivot

'---------------------------------------------

for row = pivot to last
test = abs(arrayref(a, row, pivot))
if( test > colmax) then
colmax = test
maxrow = row
end if
next

'---------------------------------------------

if(maxrow <> pivot) then
for col = first to last
temp = arrayref(a, pivot, col)
arrayref(a, pivot, col) = arrayref(a, maxrow, col)
arrayref(a, maxrow, col) = temp
next

switchcount = switchcount + 1
rowswitches(switchcount, 1) = pivot
rowswitches(switchcount, 2) = maxrow
end if

'---------------------------------------------

divisor = arrayref(a, pivot, pivot)
arrayref(a, pivot, pivot) = 1

for col = first to last
arrayref(a, pivot, col) = arrayref(a, pivot, col) / divisor
next

'---------------------------------------------

for row = first to last

if(row <> pivot) then

factor = -arrayref(a, row, pivot)
arrayref(a, row, pivot) = 0

for col = first to last
arrayref(a, row, col) = arrayref(a, row, col) + arrayref(a, pivot, col) * factor
next

end if
next

next

'---------------------------------------------

for rowswitch = switchcount to 1 step -1
i = rowswitches(rowswitch, 1)
j = rowswitches(rowswitch, 2)
for row = first to last
temp = arrayref(a, row, i)
arrayref(a, row, i) = arrayref(a, row, j)
arrayref(a, row, j) = temp
next
next

end sub

'************************************************************************************************************

sub matmatmult(m1 as integer, m2 as integer, m3 as integer)
defint i, j, k, first, last

first = byref(iadd(m1, 24))
last = byref(iadd(m1, 44))

for i = first to last
for j = first to last
arrayref(m3, i, j) = 0
for k = first to last
arrayref(m3, i, j) = arrayref(m3, i, j) + arrayref(m1, i, k) * arrayref(m2, k, j)
next
next
next
end sub

'************************************************************************************************************

sub matvecmult(m1 as integer, v1 as integer, v2 as integer)
defint i, j, k, first, last

first = byref(iadd(m1, 24))
last = byref(iadd(m1, 44))

for i = first to last
arrayref(v2, i) = 0
for j = first to last
arrayref(v2, i) = arrayref(v2, i) + arrayref(m1, i, j) * arrayref(v1, j)
next
next
end sub

'************************************************************************************************************