PureBasic --> Timing Matrix Inversions

[danbaron]

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

It was made with PureBasic, version 4.5.

The program fills an N x N matrix with random doubles in the range of [0, 2147483647], and inverts it using the same
algorithm that was used for the corresponding C, Java, and PowerBASIC 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.

From within the PureBasic IDE, there are three ways you can compile a program; "Compile with Debugger", "Compile without
Debugger", and, "Create Executable...". I used, "Create Executable...", for the times listed below. In other words, I
made the ".exe" file, and double clicked on it to time the inversions.

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

I'll compare it to my results for PowerBASIC.

A = PowerBASIC, using the internal function, "MAT INV()".
B = PowerBASIC, using my coding of the algorithm (a subroutine).
C = PureBasic, using my coding of the algorithm (a procedure).

The times on my machine were (seconds),

A B C
105.881 54.896 37.627
104.301 54.863 37.784
104.918 55.051 37.643

I attached the program file below. I called it, "MATINV.PB.TXT".

I also attached the ".exe" file, in case someone wants to try it who doesn't have PureBasic. I called it,
"MATINV.EXE.TXT". (The ".exe" file, is only 16 KB.)

If you download the files, then, remove the ".TXT" extensions.


Dan

;************************************************************************************************************

; FILE = "MATINV.PB"

; Made with the PureBasic, version 4.5.

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

; It fills an N x N matrix with random doubles in the range of [0, 2147483647],
; and inverts it using the procedure, "INV()", below.
; 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.

; 2010 - 08 - 09:
; The program uses the same inversion algorithm that was used in the corresponding, C, Java, and PowerBASIC programs.

;************************************************************************************************************

EnableExplicit

;Set the size of the N x N matrix to invert, on the next line.
#N = 1000

#NMO = #N - 1

#FN = "MATINV.TXT"

;************************************************************************************************************
;************************************************************************************************************
;************************************************************************************************************

Global DATE.S
Global TIME.S
Global T1.D
Global T2.D
Global TT.D
Global Dim A0.D(#NMO, #NMO)
Global Dim AI.D(#NMO, #NMO)
Global Dim ID.D(#NMO, #NMO)
Global Dim VA.D(#NMO)
Global Dim VB.D(#NMO)

;---------------------------------------------

Declare INITIALIZEARRAYS(Array A1.D(2), Array A2.D(2), Array V.D(1))
Declare CONSOLEOUTPUT(WHICHTIME.I, T.D)
Declare WRITEFILE(Array V.D(1), T.D)
Declare INV(Array A.D(2))
Declare MATMATMULT(Array M1.D(2), Array M2.D(2), Array M3.D(2))
Declare MATVECMULT(Array M1.D(2), Array V1.D(1), Array V2.D(1))

;---------------------------------------------

DATE = FormatDate("%yyyy-%mm-%dd", Date())
TIME = FormatDate("%hh:%ii:%ss", Date())

;---------------------------------------------

Define DUMB.S

;---------------------------------------------

OpenConsole()

INITIALIZEARRAYS(A0(), AI(), VA())

CONSOLEOUTPUT(1, TT)

T1 = ElapsedMilliseconds()
INV(AI())
T2 = ElapsedMilliseconds()
TT = (T2 - T1) / 1000

CONSOLEOUTPUT(2, TT)

MATMATMULT(A0(), AI(), ID())
MATVECMULT(ID(), VA(), VB())

WRITEFILE(VB(), TT)

CONSOLEOUTPUT(3, TT)
DUMB = Input()

;************************************************************************************************************
;************************************************************************************************************
;************************************************************************************************************

Procedure INITIALIZEARRAYS(Array A1.D(2), Array A2.D(2), Array V.D(1))
Define FIRST.I, LAST.I, I.I, J.I, D.D

FIRST = 0
LAST = ArraySize(V())

For I = FIRST To LAST
V(I) = 1
For J = FIRST To LAST
D = Random(2147483647)
A1(I, J) = D
A2(I, J) = D
Next
Next
EndProcedure

;************************************************************************************************************
;************************************************************************************************************
;************************************************************************************************************

Procedure CONSOLEOUTPUT(WHICHTIME.I, T.D)
Define SO.S, SN.S, ST.S

SN = Str(#N)
ST = StrD(T, 3)

Select WHICHTIME

Case 1
SO = "Inverting a " + SN + " by " + SN + " matrix.."
PrintN("")
PrintN(SO)
PrintN("")

Case 2
PrintN("Done.")
PrintN("")
SO = "elapsed time = " + ST + " seconds"
PrintN(SO)
PrintN("")
PrintN("Calculating and writing the check vector..")
PrintN("")

Case 3
PrintN("Done.")
PrintN("")
PrintN("Press 'Enter' to quit.")
PrintN("")

Default
PrintN("Something wrong.")
PrintN("")

EndSelect

EndProcedure

;************************************************************************************************************
;************************************************************************************************************
;************************************************************************************************************

Procedure WRITEFILE(Array V.D(1), T.D)
Define FIRST.I, LAST.I, I.I
Define S1.S, S2.S, SN.S, SO.S, ST.S

FIRST = 0
LAST = ArraySize(V())

SN = Str(#N)
ST = StrD(T, 3)

CreateFile(0, #FN)

S1 = "; FILE = " + "'" + #FN + "'"
WriteStringN(0, S1)
S1 = "; " + DATE + ", " + TIME + "."
WriteStringN(0, S1)
WriteStringN(0, "")
S1 = "; Inversion of a " + SN + " by " + SN + " matrix,"
WriteStringN(0, S1)
WriteStringN(0, "; filled with uniformly distributed random doubles,")
WriteStringN(0, "; in the range of [0, 2147483647].")
WriteStringN(0, "")
S1 = "; elapsed time = " + ST + " seconds"
WriteStringN(0, S1)
WriteStringN(0, "")
WriteStringN(0, "; If the inversion was calculated correctly, then every")
WriteStringN(0, "; value shown below should be very close to 1.")
WriteStringN(0, "")
WriteStringN(0, "Index      Value")

For I = FIRST To LAST
S1 = RSet(Str(I + 1), 4, "0")
S1 = S1 + "  "
S2 = StrD(V(I), 20)
S1 = S1 + S2
WriteStringN(0, S1)
Next

CloseFile(0)

EndProcedure

;************************************************************************************************************
;************************************************************************************************************
;************************************************************************************************************

Procedure INV(Array A.D(2))
;Uses the Gauss-Jordan Method.
;Inverts matrix, A(), "in-place".

Define PIVOT.I, ROW.I, COL.I, MAXROW.I, ROWSWITCH.I, SWITCHCOUNT.I, FIRST.I, LAST.I
Define TEST.D, TEMP.D, COLMAX.D, DIVISOR.D, FACTOR.D, TEMP.D

;---------------------------------------------

FIRST = 0
LAST = ArraySize(A(), 1)

Dim ROWSWITCHES.I(#NMO, 1)

SWITCHCOUNT = -1

For PIVOT = FIRST To LAST

COLMAX = 0
MAXROW = PIVOT

;---------------------------------------------

For ROW = PIVOT To LAST
TEST = Abs(A(ROW, PIVOT))
If( TEST > COLMAX)
COLMAX = TEST
MAXROW = ROW
EndIf
Next

;---------------------------------------------

If(MAXROW <> PIVOT)
For COL = FIRST To LAST
TEMP = A(PIVOT, COL)
A(PIVOT, COL) = A(MAXROW, COL)
A(MAXROW, COL) = TEMP
Next

SWITCHCOUNT = SWITCHCOUNT + 1
ROWSWITCHES(SWITCHCOUNT, 0) = PIVOT
ROWSWITCHES(SWITCHCOUNT, 1) = MAXROW
EndIf

;---------------------------------------------

DIVISOR = A(PIVOT, PIVOT)
A(PIVOT, PIVOT) = 1

For COL = FIRST To LAST
A(PIVOT, COL) = A(PIVOT, COL) / DIVISOR
Next

;---------------------------------------------

For ROW = FIRST To LAST

If(ROW <> PIVOT)

FACTOR = -A(ROW, PIVOT)
A(ROW, PIVOT) = 0

For COL = FIRST To LAST
A(ROW, COL) = A(ROW, COL) + A(PIVOT, COL) * FACTOR
Next

EndIf
Next

Next

;---------------------------------------------

For ROWSWITCH = SWITCHCOUNT To 1 Step -1
For ROW = FIRST To LAST
TEMP = A(ROW, ROWSWITCHES(ROWSWITCH, 0))
A(ROW, ROWSWITCHES(ROWSWITCH, 0)) = A(ROW, ROWSWITCHES(ROWSWITCH, 1))
A(ROW, ROWSWITCHES(ROWSWITCH, 1)) = TEMP
Next
Next

EndProcedure

;************************************************************************************************************
;************************************************************************************************************
;************************************************************************************************************
;************************************************************************************************

Procedure MATMATMULT(Array M1.D(2), Array M2.D(2), Array M3.D(2))
Define I.I, J.I, K.I, FIRST.I, LAST.I

FIRST = 0
LAST = ArraySize(M1(), 1)

For I = FIRST To LAST
For J = FIRST To LAST
M3(I, J) = 0; 
For K = FIRST To LAST
 M3(I, J) = M3(I, J) + M1(I, K) * M2(K, J);
Next
Next
Next
EndProcedure

;************************************************************************************************
;************************************************************************************************
;************************************************************************************************

Procedure MATVECMULT(Array M1.D(2), Array V1.D(1), Array V2.D(1))
Define I.I, J.I, K.I, FIRST.I, LAST.I

FIRST = 0;
LAST = ArraySize(M1(), 1)

For I = FIRST To LAST
V2(I) = 0;
For J = FIRST To LAST
V2(I) = V2(I) + M1(I, J) * V1(J) 
Next
Next
EndProcedure

;************************************************************************************************
;************************************************************************************************
;************************************************************************************************