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.

:oops: :x :twisted:
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

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

[danbaron]

[font=courier new][size=8pt]I found something that the compiler missed. If you look at the top of the program, you'll see that I used the compiler directive,
"EnableExplicit". That means that you are required to declare every variable. Now look at this little part of the program.


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

[font=courier new][size=8pt](The suffix, ".I", means the variable is declared as an integer, and, the suffix, ".D", means the variable is declared as a double.)

Anyway, one of the variables is declared twice. And the compiler didn't complain. I guess it doesn't really matter. My assumption is that there is only one
variable, "TEMP", since the program seems to work OK.

It seems hard for a compiler or interpreter to catch every error. I think the reason is, that so far, there is no way to automate the process. There is no
mathematical algorithm I know of, that can analyze a program and find every error. I think it still requires a human to define what constitutes an error, no
mathematical formula will do it. And, apparently, it is very difficult for humans to think of, and so check for, every possible error in advance, so that upon
initial release, a program is, "perfect". It seems that as the length of a program increases, the probability of ever finding every potential bug, goes to zero.
Some bugs might appear only once every 10 years, because the particular configuration of the input which causes them, is so unlikely. I think that may be where
the idea of "functional" languages, like Haskell, came from. If the input which a program uses cannot be modified by the program, then, most likely, finding all
of the potential bugs, would become much easier. If I remember correctly, one of the goals of the development of functional languages, is to be able to
mathematically prove that a program is correct, that it has no bugs. But, I think that in order for a proof to be theoretically possible - simply put, a program
must have no variables, only constants. Many programmers who do not work for universities, seem to find that requirement to be, quite limiting. On the other
hand, imagine a program with millions or billions of line of code, which, for instance, flies a pilotless passenger jet in the future. Maybe in that case, it
would be worth the extra effort that using a functional language requires, in order to be able to mathematically prove the program is bug-free, before it is
ever used in real life.

The other way to check a program that I can think of, would be (or is) to run it for a long time, using random input. In other words, the input for the program
you are testing, would be randomly generated by a second program - the two programs would interact. Suppose you run a program continuously for a year using random
input, and it never fails. Then, you could mathematically estimate the probability that it will fail in, say, 10 years of actual use. It would still be up to
humans to decide what the acceptable probability of failure is.

And, of course, errors and bugs are not the same thing, yes? To me, "errors", imply syntax and logic. I think, "bugs", are harder. To me, they usually mean that
program variables are related in such a complex way, that it becomes humanly impossible to find every possibility for disaster. (I agree that "disasters" are
relative. A disaster in a word processing program, might cause the program to crash. A disaster in a jet flying program, might cause the jet to crash.)

(I've gone on too long. I should have quit sooner. I'm babbling.)

:oops: :x :twisted:
Dan

[ErosOlmi]

Send them a bug report.

In any case I think it is a wanted behave. If you define the same variable name but with a different type, compiler will complain.

[danbaron]

[font=courier new][size=8pt]I filed a bug report, but I also do not think that they are unaware of the behavior.

There is an Air Force Base not far from my house. I remember an incident approximately 10 years ago. The Air Force had a new fighter jet, I forget which one.
The pilot landed the jet at the base. Apparently, everything on the jet was computerized. Nothing could be performed manually. After he landed the jet, I guess
he clicked the mouse to open the canopy (the plastic bubble) so he could get out. Nothing happened. There was a bug in the software. He sat on the runway inside
the jet for 5 hours. The fire department came and cut through the 1.25 inch thick (if I remember correctly) polycarbonate, to extricate him.

I can imagine that in the future, automobiles might be the same way. If you press the button to open the door and nothing happens, then you will have no choice
but to sit there, or to call for help. I don't like that prospect.

Imagine if passenger jetliners become totally computerized. What could potentially happen if one was in flight and it was struck by lightning, or by a weapon
which fired a directed electric pulse?

[Charles Pegge]

That pilot was lucky. There are quite a few aviation software horror stories in circulation. The software has to be tested with the same rigour as all the other components of the aircraft. I think this means you cannot rely on the theoretical correctness of your software and its compiler. Computers are finite-state machines and vulnerable to various limitations. (Not to mention getting zapped by high energy particles)

I think the best way to develop a reliable system is to co-develop test software alongside the main program. The test software may well end up much larger than the system being tested.

Charles

[danbaron]

[font=courier new][size=8pt]Imagine in the future, a giant passenger jet that is completely computerized. Everything on it, is controlled by its computer program.

Imagine that this is a part of the program's code.


IF V79Q14A = 1 THEN OPEN_ALL_OF_THE_DOORS()

[font=courier new][size=8pt]The programmers might be absolutely convinced that it is impossible for variable V79Q14A to equal 1, while the jet is in flight.

But, how many times does what is seemingly impossible occur during the execution of code?

When a change in one variable can affect the values of other variables, and there are many thousands of variables, I think it can become extremely difficult to
determine, what is or is not possible.

You can have long long chains of effect, --> A changes B, B changes C, .., Z47 changes QR48, .., NTW387 changes HD52, ..

Giant programs are made by maybe hundreds of programmers, supposedly collaborating. Maybe in the future it will be possible for computers to check millions of
lines of code, but one person could never do it, and, as soon as there are more than one person, miscommunication begins. And, if you think about it, how
reliable does it seem to divide up the checking of a program? Would you feel safe if, say, a program was divided into four parts, and one person checked each of
the parts? In that case, no one checks the program in its entirety. I agree that at least for now, "stress" testing is by far the best way to maximize
reliability, but, who will claim that it is infallible?

For a big computer program, a particular person may understand a particular part of it in detail. But, who understands all of it? I think the answer is, no one.

Going back to passenger jets being entirely computerized, imagine when the latest software is installed in the entire fleet, and, along with it comes a virus.
If it isn't detected, a thousand passenger jets could simultaneously fall out of the sky.

(The rumors are that even now, the newer big passenger jets can be remotely controlled, that what the pilot does can be overridden. The new big jets are, "fly
by wire", right? That is, all required mechanical force is provided by machinery, which is computer controlled. Who knows? Maybe the real reason that there are
still commercial pilots, is only to prevent the passengers from becoming terrified, and refusing to fly.)

:twisted: