'----------------------------------------------------------------------------------
'START PROGRAM HEADER
'----------------------------------------------------------------------------------

'FILE = MATINV.EBA

'EB CONSOLE PROGRAM

'THIS PROGRAM USES THE GAUSS-JORDAN METHOD TO INVERT MATRICES,
'AND TIMES THE INVERSIONS.

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

$INCLUDE "WINDOWS.INC"

'----------------------------------------------------------------------------------
'END PROGRAM HEADER
'----------------------------------------------------------------------------------
/*
-----------------------------------------------------------------------------------
START PROGRAM SUMMARY
-----------------------------------------------------------------------------------

This is an EB console program that uses the Gauss-Jordan method to
invert matrices.

For an N x N matrix, you set the value of the program constant, N.

All of the matrices used in the program are static. Each time you change
the value of N, you must re-compile the program.

When the program executes, the N x N matrix is randomly filled with
DOUBLEs, in the range of (-1000, 1000).

Then, the matrix is inverted.

The time to invert the matrix is calculated, and displayed in the
console window.

After the inversion, a check is made to determine if the inversion is
correct.

Let's call the original matrix, "A", and its inversion, "A'".

After the inversion, we have both A and A', i.e., the inversion does not
overwrite A.

In order to check the inversion, we multiply A and A'. Their product
should be the N x N identity matrix, "I".

We then multiply a vector of length N, with all values set equal to 1, by
the I we calculated.

According to the properties of the identity matrix, the multiplication
should not alter the vector, i.e., every value should still equal 1.

The check vector is written to a text file, so that you can verify that
the inversion worked. If all of the vector's values are very close to 1,
then the inversion is good. The elapsed time is also written to the text
file.

-----------------------------------------------------------------------------------
END PROGRAM PROGRAM SUMMARY
-----------------------------------------------------------------------------------
START PROGRAM OPERATION
-----------------------------------------------------------------------------------

The only thing you need to do, is set the constant, "MATSIZE", just
below. It corresponds to, "N", for an N x N matrix.

Then, compile the program as a "Console EXE", and run it.

After execution completes, you can open and look in the file,
"MATINV.TXT" (by default), to see if the inversion is correct.

-----------------------------------------------------------------------------------
END PROGRAM OPERATION
-----------------------------------------------------------------------------------
*/
'----------------------------------------------------------------------------------
'START USER SET CONSTANTS
'----------------------------------------------------------------------------------

'Set the size, "N", of the N x N matrix to be inverted.

CONST MATSIZE = 200

'Set the value to seed the random number generator, RANDSEED.
'RANDSEED cannot be 0, if it is, it will be changed to 1.
'If RANDSEED is negative, it will be set to the number of
'elapsed seconds in the current day.

CONST RANDSEED = -1

'After the matrix is inverted, a vector which checks the inversion,
'is calculated and written to a file, by default, "MATINV.TXT".

CONST FILESTRING = "MATINV.TXT"

'----------------------------------------------------------------------------------
'END USER SET CONSTANTS
'----------------------------------------------------------------------------------
'NOTHING SHOULD NEED TO BE CHANGED BELOW THIS POINT.
'----------------------------------------------------------------------------------
/*
-----------------------------------------------------------------------------------
START PROGRAM EXPLANATION
-----------------------------------------------------------------------------------
START EXPLANATION OF GAUSS-JORDAN METHOD OF MATRIX INVERSION
-----------------------------------------------------------------------------------

VECTORS AND DOT PRODUCTS

A "vector", is a one dimensional array.

The dot product, "DP", of two N dimensional vectors, is a
scalar (number), which is the sum of the products of the
values of corresponding indices.

Say, for instance, we have two 3 dimensional INT vectors, V1, and V2.

Let,

V1 = (4, -2, 6)
V2 = (-9, 7, 3)

The dot product of V1 and V2 can be represented in EB as

INT V1[4]
INT V2[4]
INT I, DP

V1 =  0, 4, -2, 6
V2 =  0, -9, 7, 3

DP = 0
FOR I = 1 TO 3
DP += V1[I] * V2[I]
NEXT I

Therefore, DP = 4 * -9 + -2 * 7 + 6 * 3 = -32.

(Here, for both vectors and matrices, we will ignore the 0 index.
So, for a 3 dimensional vector, "V", we will dimension it as V[4].)

In order for the dot product to be defined, the two vectors
must be of the same dimension (length).

-----------------------------------------------------------------------------------

MATRICES AND MATRIX MULTIPLICATION

A "matrix" is a two dimensional array.


Here is an example of a 2 x 2 matrix, call it "A".

A[1, 1] = 2
A[1, 2] = 3
A[2, 1] = 4
A[2, 2] = 5

We'll show it like this.
  
---------
| 2 | 3 |
|-------|
| 4 | 5 |
---------

Matrix multiplication of two 2 x 2 matrices, A, and B, giving the matrix, C,
can be shown in EB as

INT I, J, K
INT A[3, 3]
INT B[3, 3]
INT C[3, 3]

A[1, 1] = 2
A[1, 2] = 3
A[2, 1] = 4
A[2, 2] = 5

B[1, 1] = 9
B[1, 2] = 7
B[2, 1] = 6
B[2, 2] = 1

FOR I = 1 TO 2
FOR J = 1 TO 2
C[I, J] = 0
FOR K = 1 TO 2
C[I, J] += A[I, K] * B[K, J]
NEXT K
NEXT J
NEXT I

Matrix multiplication of two matrices, produces a third matrix.

The entry C[I, J], equals the dot product of the Ith row of A, with the
Jth column of B. From the definition of the dot product, this means that
the length of the Ith row of A, must equal the length of the Jth column
of B. The length of the Ith row of A, equals the number of columns of A.
The length of the Jth column of B, equals the number of rows of B.
Therefore, in order for the matrix product of two matrices to exist, the
number of columns in the first matrix, must equal the number of rows in
the second matrix. In that case, the matrices are called "conformable".

In order to be able to multiply two matrices, the dimensions of the
matrices must have the form

[M, N] * [N, P]

(For any matrix, the first dimension is the number of rows, and the
second dimension is the number of columns. Note also, that a N
dimensional vector, is also a matrix; with dimensions [N, 1] (column
vector), or [1, N] (row vector).)

The adjacent dimensions, "N", must be equal.
In that case, the dimensions of the product matrix will be [M, P].

Notice what happens if we reverse the multiplication.

[N, P] * [M, N]

Unless, P equals M, the matrices are "non-conformable", and cannot be
multiplied.

We can graphically show the multiplication, like this. Then, each entry
in the product matrix, is the dot product of the intersection of the row
of the first multiplying matrix, with the column of the second
multiplying matrix.

           ---------
           | 9 | 7 |
           |-------|
           | 6 | 1 |
           ---------
---------  ---------
| 2 | 3 |  | A | B |
|-------|  |-------|
| 4 | 5 |  | C | D |
---------  ---------

A = 2 * 9 + 3 * 6 = 18
B = 2 * 7 + 3 * 1 = 17
C = 4 * 9 + 5 * 6 = 66
D = 4 * 7 + 5 * 1 = 33

So,

---------   ---------   -----------
| 2 | 3 |   | 9 | 7 |   | 18 | 17 |
|-------| * |-------| = |---------|
| 4 | 5 |   | 6 | 1 |   | 66 | 33 |
---------   ---------   -----------

Let's show an example of multiplication of non-square matrices.

                   ---------
                   | 7 | 9 |
                   |-------|
                   |-3 |-6 |
                   |-------|
                   | 1 | 8 |
                   |-------|
                   | 5 |-4 | 
			        	   ---------
-----------------  ---------
| 8 | 3 | 2 | 6 |  | A | B |
|---------------|  |-------|
|-2 | 5 | 0 | 4 |  | C | D |
|---------------|  |-------|
| 9 | 1 |-7 | 11|  | E | F |
-----------------  ---------

The dimensions of the multiplying matrices are, [3, 4] and [4, 2].

So, the matrices are conformable, and the product matrix has the
dimensions, [3, 2].

A =   8  * 7 + 3 * (-3) +   2  * 1 +  6 *   5  =  79
B =   8  * 9 + 3 * (-6) +   2  * 8 +  6 * (-4) =  46
C = (-2) * 7 + 5 * (-3) +   0  * 1 +  4 *   5  =  -9
D = (-2) * 9 + 5 * (-6) +   0  * 8 +  4 * (-4) = -64
E =   9  * 7 + 1 * (-3) + (-7) * 1 + 11 *   5  = 108
F =   9  * 9 + 1 * (-6) + (-7) * 8 + 11 * (-4) = -25

But, if we reverse the multiplication, we get,

[4, 2] * [3, 4]

The matrices are non-conformable, and the matrix product does not exist.

Generally, if A and B are two matrices, then A * B, is not equal to
B * A. And one or both of the products may not exist. Except in special
cases, matrix multiplication is not commutative ( A * B is not equal to
B * A). This is even true for the multiplication of square matrices
having the same dimension.

-----------------------------------------------------------------------------------

THE IDENTITY MATRIX

The N dimensional identity matrix, is an N x N matrix with the following
properties.

I[i, j] = 0 if i <> j

I[i, j] = 1 if i = j

The two dimensional identity matrix is

---------
| 1 | 0 |
|-------|
| 0 | 1 |
---------

The diagonal extending from upper left to lower right, is called the
"main diagonal".

For any identity matrix, I, the entries on the main diagonal have the
value 1, and all other entries have the values 0.

The three dimensional identity matrix is

-------------  
| 1 | 0 | 0 |  
|-----------|  
| 0 | 1 | 0 |  
|-----------|  
| 0 | 0 | 1 |  
-------------  

The identity matrix has the property, that any matrix it multiplies, or
is multiplied by, is unchanged.

That is, for any matrix, A, if A and I are conformable for the
multiplication, A * I, then,

A * I = A

And, for any matrix, A, if A and I are conformable for the 
multiplication, I * A, then,

I * A = A

The identity matrix is the analog to the integer 1, in the
multiplication of numbers.

Here is an example

           ----------------- 
           | 1 | 2 | 3 | 4 |
           |---------------| 
           | 5 | 6 | 7 | 8 |
           ----------------- 
---------  ----------------- 
| 1 | 0 |  | A | B | C | D | 
|-------|  |---------------| 
| 0 | 1 |  | E | F | G | H |
---------  ----------------- 

You can check that

A = 1
B = 2
C = 3
D = 4
E = 5
F = 6
G = 7
H = 8

Here is an example of a row vector multiplying I

               ------------- 
               | 1 | 0 | 0 | 
               |-----------| 
               | 0 | 1 | 0 | 
               |-----------| 
               | 0 | 0 | 1 | 
               ------------- 
-------------  -------------		  
| 6 |-3 | 0 |  | A | B | C | 	  
-------------  -------------

You can check that

A =  6
B = -3
C =  0 	  

-----------------------------------------------------------------------------------

MATRICES WITHOUT INVERSES

We'll only consider the matrix inverse, for square matrices (N x N). Not
every square matrix has an inverse. A square matrix which does not have
an inverse, is called, "singular".

Here are two ways to show that a matrix is singular, and has no
inverse.


1)

If one row is a linear combination of one or more of the other rows,
then the matrix is singular.

For instance, call this matrix "A".

------------- 
|-3 | 5 | 6 | 
|-----------| 
|12 |17 | 9 | 
|-----------| 
| 2 | 9 | 7 | 
------------- 

This matrix is singular, because

row 2 = -2 * (row 1) + 3 * (row 3)


2)

If a non-zero vector, "V", can be found such that

A * V = 0

then A is singular. (Here "0" means the zero vector, all values equal
0.)

For the matrix A, above, it turns out that any vector V = (x, y, z) such
that

x =  (57/111) * z    (z can be any non-zero value)
y = (-99/111) * z

will give

A * V = 0

------------------------------------------------------------------------------

THE MATRIX INVERSE

For an N x N matrix, "A", its inverse is an N x N matrix that we'll call
"A'".

A' is the inverse of A, if

A * A' = A' * A = I   (the identity matrix)

The matrix inverse, is the analog of the reciprocal in the
multiplication of numbers. For instance

   7     *      1/7     =      1
(number)   (reciprocal)


Here is an example of a matrix and its inverse.

      A                     A'                    I
-------------   -------------------------   ------------- 
| 1 | 2 | 3 |   | -5/18 |  1/18 |  7/18 |   | 1 | 0 | 0 | 
|-----------|   |-----------------------|   |-----------| 
| 2 | 3 | 1 | * |  1/18 |  7/18 | -5/18 | = | 0 | 1 | 0 | 
|-----------|   |-----------------------|   |-----------| 
| 3 | 1 | 2 |   |  7/18 | -5/18 |  1/18 |   | 0 | 0 | 1 | 
-------------   -------------------------   ------------- 

A and A' are each the inverse of the other.

------------------------------------------------------------------------------

CALCULATING THE MATRIX INVERSE

We'll calculate the matrix inverse by the Gauss-Jordan Method.

The method relies on the three "Elementary Row Operations".

These operations are:


1)
Switching two rows of a matrix.

2)
Replacing a row by the same row multiplied by a non-zero constant.

3)
Replacing a row by the same row plus a constant multiple of another row.


The Gauss-Jordan Method uses the elementary row operations to transform a
matrix, "A", to the identity matrix, I. Simultaneously, it performs the
same operations that transform A to I, on another matrix which is
initially I. As A changes to I, I changes to A', the inverse of A.

We'll do an example, to make the method clear, then we'll show why the
method works.

Say, we want to find the inverse of the matrix,

     A
---------
| 2 | 3 |
|-------|
| 4 | 5 |
---------

To do that, we will use the elementary row operations to transform it
to
    I
---------
| 1 | 0 |
|-------|
| 0 | 1 |
---------

At the same time, every operation we perform on A, we will also perform
on a matrix which begins as I. The confusion should vanish as we perform
the procedure.

First, we'll write A and I next to each other.

    A           I
---------   ---------
| 2 | 3 |   | 1 | 0 |
|-------|   |-------|
| 4 | 5 |   | 0 | 1 |
---------   ---------

1)
For the 1st step, we'll interchange the first and second rows of A. This
step is not necessary in this case, but we'll do it anyway, just to show
that we can.  Remember, that whatever we do to A, we must also do to I,
so we'll switch its first and second row too. Now, we'll call A, "Am1",
and I, "Im1", to indicate that they have been modified.

   Am1          Im1
---------   ---------
| 4 | 5 |   | 0 | 1 |
|-------|   |-------|
| 2 | 3 |   | 1 | 0 |
---------   ---------

2)
Now, we'll scale the first row of Am1 by 1/4, and do the same to Im1.

   Am2             Im2
-----------   -----------
| 1 | 5/4 |   | 0 | 1/4 |
|---------|   |---------|
| 2 |  3  |   | 1 |  0  |
-----------   -----------

3)
Now, we'll replace row two by itself, minus 2 times row one.

   Am3             Im3
-----------   ------------
| 1 | 5/4 |   | 0 |  1/4 |
|---------|   |----------|
| 0 | 1/2 |   | 1 | -1/2 |
-----------   ------------

4)
Now, we'll scale row two by 2.

   Am4             Im4
-----------   ------------
| 1 | 5/4 |   | 0 |  1/4 |
|---------|   |----------|
| 0 |  1  |   | 2 |  -1  |
-----------   ------------

5)
Now, we'll replace row one by itself, minus 5/4 times row two.

   Am5              Im5
-----------   --------------
| 1 |  0  |   | -5/2 | 3/2 |
|---------|   |------------|
| 0 |  1  |   |   2  | -1  |
-----------   --------------

Now, we are done. Notice that A has become I. And, we claim that I has
become A'. We can multiply A and Im5 to find out. Remember, that their
product should be I.

                 Im5
           --------------
           | -5/2 | 3/2 |
           |------------|
           |   2  | -1  |
    A      --------------                             
---------  --------------
| 2 | 3 |  |   A  |  B  |
|-------|  |------------|
| 4 | 5 |  |   C  |  D  |
---------  --------------

A = 2 * (-5/2) + 3 *   2  = 1
B = 2 *  (3/2) + 3 * (-1) = 0
C = 4 * (-5/2) + 5 *   2  = 0
D = 4 *  (3/2) + 5 * (-1) = 1

So, we have succeeded.

    A              A'            I
---------   --------------   ---------
| 2 | 3 |   | -5/2 | 3/2 |   | 1 | 0 |
|-------| * |------------| = |-------|
| 4 | 5 |   |   2  | -1  |   | 0 | 1 |
---------   --------------   ---------

Or,

       A'            A           I     
--------------   ---------   --------- 
| -5/2 | 3/2 |   | 2 | 3 |   | 1 | 0 | 
|------------| * |-------| = |-------| 
|   2  | -1  |   | 4 | 5 |   | 0 | 1 | 
--------------   ---------   --------- 

So, we have shown the steps of the method. Now, we must show why it
works.

To show why the method works, first we'll show that each of the
elementary row operations we applied to A, was equivalent to multiplying
A by another matrix. We'll perform the 5 transformations we did to A,
again. Only, this time we'll do them by multiplying A by 5 matrices.


1)
   m1           A           Am1     
---------   ---------   --------- 
| 0 | 1 |   | 2 | 3 |   | 4 | 5 | 
|-------| * |-------| = |-------| 
| 1 | 0 |   | 4 | 5 |   | 2 | 3 | 
---------   ---------   --------- 


2)
     m2          Am1          Am2     
-----------   ---------   ----------- 
| 1/4 | 0 |   | 4 | 5 |   | 1 | 5/4 | 
|---------| * |-------| = |---------| 
|  0  | 1 |   | 2 | 3 |   | 2 |  3  | 
-----------   ---------   -----------

 
3)
    m3           Am2          Am3     
----------   -----------   ----------- 
|  1 | 0 |   | 1 | 5/4 |   | 1 | 5/4 | 
|--------| * |---------| = |---------| 
| -2 | 1 |   | 2 |  3  |   | 0 | 1/2 | 
----------   -----------   ----------- 


4)
    m4           Am3          Am4     
---------   -----------   ----------- 
| 1 | 0 |   | 1 | 5/4 |   | 1 | 5/4 | 
|-------| * |---------| = |---------| 
| 0 | 2 |   | 0 | 1/2 |   | 0 |  1  | 
---------   -----------   -----------

 
5)
     m5            Am4           Am5     
------------   -----------   ----------- 
| 1 | -5/4 |   | 1 | 5/4 |   | 1 |  0  | 
|----------| * |---------| = |---------| 
| 0 |   1  |   | 0 |  1  |   | 0 |  1  | 
------------   -----------   ----------- 

So, we see, that the transformations we did to A using the elementary row
operations, are equivalent to multiplying A by five matrices.

Let's see what we get, when we multiply the five matrices (m1 - m5)
together. Since m1 multiplied A first, and m5 multiplied A last, we'll
multiply the five matrices in reverse order.

     m5           m4          m3           m2            m1              A' 
------------   ---------   ----------   -----------   ---------   --------------
| 1 | -5/4 |   | 1 | 0 |   |  1 | 0 |   | 1/4 | 0 |   | 0 | 1 |   | -5/2 | 3/2 |
|----------| * |-------| * |--------| * |---------| * |-------| = |------------|
| 0 |   1  |   | 0 | 2 |   | -2 | 1 |   |  0  | 1 |   | 1 | 0 |   |   2  | -1  |
------------   ---------   ----------   -----------   ---------   --------------

We see that when we multiply m5 * m4 * m3 * m2 * m1, we get A', the
inverse of A. Therefore, the transformation of A to I using the
elementary row operations, is equivalent to multiplying A by A'.

But, remember, that we we applied the same elementary row operations to
a matrix which began as I. So, effectively, we multiplied I by A'. And
since, the product of multiplying a conforming matrix by I, is the
matrix, we end up with A'.

That is how the Gauss-Jordan method of matrix inversion works.

(Here is a third way to determine if a square matrix is singular, and
has no inverse. In that case, it will be impossible to transform it into
the identity matrix, by using the elementary row operations.)

-----------------------------------------------------------------------------------
END EXPLANATION OF GAUSS-JORDAN METHOD OF MATRIX INVERSION
-----------------------------------------------------------------------------------
START EXPLANATION OF PROGRAM'S USE OF 1D ARRAYS TO REPRESENT 2D MATRICES
-----------------------------------------------------------------------------------

My experience, is that in EB, if you dimension a (2D) matrix in the parameter
list of a subroutine, using a constant, the subroutine will not function
properly.

For instance, here I have defined a global constant, "D", and dimensioned a
global matrix, "A", using it. That works fine.

CONST D = 10
DEF A[D, D] AS DOUBLE

But, if I then write a subroutine like this, the program will link and
run, but the subroutine's behavior will be incorrect.
  
SUB PROCESSMATRIX(M[D, D] AS DOUBLE)
.
.
ENDSUB

To get correct behavior, I would have to do this.

SUB PROCESSMATRIX(M[10, 10] AS DOUBLE)
.
.
ENDSUB

However, for (1D) arrays, you do not have to put the size of the array
in a subroutine's parameter list. You can just do this.

SUB PROCESSMATRIX(M[] AS DOUBLE)
.
.
ENDSUB

Within the subroutine, you can still use the global constant, D.

Therefore, in the program, I represent 2D matrices, using 1D arrays.

I still use the 2D indexing, but it is transformed into 1D indexing.

So, in the program, instead of writing,

M[I, J] = X

I write,

SETARRAYVAL(M, D, D, I, J, X)

(Here, D, is the same global constant, as above.)

Instead of writing,

X = M[I, J]

I write,

X = GETARRAYVAL(M, D, D, I, J)

Instead of writing,

M[I, J] += X

I write,

ADDTOARRAYVAL(M, D, D, I, J, X)

instead of writing,

M[I, J] /= X

I write,

DIVIDEARRAYVAL(M, D, D, I, J, X)

-----------------------------------------------------------------------------------
END EXPLANATION OF PROGRAM'S USE OF 1D ARRAYS TO REPRESENT 2D MATRICES
-----------------------------------------------------------------------------------
START MISCELLANEOUS
-----------------------------------------------------------------------------------

1)
Since the matrix to be inverted is filled with random DOUBLEs, the
probability of it being singular, and therefore without an inverse, is
tiny.

2)
When using the Gauss-Jordan method to invert an N x N matrix, there are
N cubed (N ^ 3) calculations.

Therefore, the time to invert an N x N matrix is,

t(N x N) = C * (N ^ 3)

where, C is a constant.

For an M x M matrix, the time is,

t(M x M) = C * (M ^ 3)

Therefore,

t(M x M) / t(N x N) = (C * (M ^ 3)) / (C * (N ^ 3))

So,

t(M x M) = ((M / N) ^ 3) * t(N x N)

-----------------------------------------------------------------------------------
END MISCELLANEOUS
-----------------------------------------------------------------------------------
END PROGRAM EXPLANATION
-----------------------------------------------------------------------------------
*/
'----------------------------------------------------------------------------------
'START GLOBAL CONSTANTS AND VARIABLES
'----------------------------------------------------------------------------------

'utility constants

CONST MS = MATSIZE
CONST MSP1 = MS + 1
CONST NN = MS * MS
CONST NNP1 = NN + 1

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

'for function RAN2() *******************************************

CONST IM1       = 2147483563
CONST IM2       = 2147483399
CONST IMM1      = IM1 - 1
CONST IA1       = 40014
CONST IA2       = 40692
CONST IQ1       = 53668
CONST IQ2       = 52774
CONST IR1       = 12211
CONST IR2       = 3791
CONST NTAB      = 32
CONST NDIV      = 1 + IMM1 / NTAB
CONST DEPS      = 1.0E-16
INT   RAN2SEED  = 0
INT   IDUM2     = 0
INT   IY        = 0
CHAR  RAN2FLAG  = 0
INT   IV[NTAB]
INT   INITIALRAN2SEEDVAL

'for function RAN2() *******************************************

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

'Matrix to be inverted, "A".

DEF A[NNP1] AS DOUBLE

'Copy of A, "AM". During the inversion, AM will be transformed to the identity matrix, "I".

'AM = "A modified"

DEF AM[NNP1] AS DOUBLE

'Matrix to hold the inversion of A, "AI".

'AI = "A inversion"

DEF AI[NNP1] AS DOUBLE

'Identity matrix, which is calculated after the inversion, "ID", product of A times AI.

'ID = "identity matrix"

DEF ID[NNP1] AS DOUBLE

'utility vectors

DEF V1[MSP1] AS DOUBLE
DEF V2[MSP1] AS DOUBLE

'Note, that for the inversion of an N x N matrix, A, we will use four
'N x N matrices, A, AM, AI, and ID.

'A is the original matrix, which is unchanged.

'AM begins as a copy of A, and then is transformed into the N x N identity
'matrix, "I", using the Gauss-Jordan method.

'AI begins as the N x N identity matrix, I, and then is transformed into the
'inversion of A, using the Gauss-Jordan method.

'After the inversion, A and AI are multiplied together, and ID is
'used to hold their product. If the inversion succeeded, then ID should
'equal the N x N identity matrix, I.

'To check if the inversion is correct, ID is multiplied by an
'N-dimensional vector, with all values preset to 1. If, after the
'multiplication, all of the values still equal 1, then the
'inversion was successful.

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

FILE OUTFILE
SYSTEMTIME ST0, ST1
DOUBLE ET

'----------------------------------------------------------------------------------
'END GLOBAL CONSTANTS AND VARIABLES
'----------------------------------------------------------------------------------
'START PROGRAM
'----------------------------------------------------------------------------------

IF MATSIZE < 1 THEN END

OPENCONSOLE

SETRAN2SEED(RANDSEED)
INITIALIZEMATRICES()

PRINT "INVERTING MATRIX."
PRINT

STARTTIMER1()
INVERTMAT(AM, AI)
STOPTIMER1()

PRINT "DONE."
PRINT

CALCELAPSEDTIME()
SHOWELAPSEDTIME()

PRINT
PRINT "CALCULATING CHECK VECTOR."
PRINT

CALCCHECKVECTOR()
WRITECHECKVECTOR()

PRINT "DONE."
PRINT
PRINT "PRESS (ALMOST) ANY KEY TO QUIT."
PRINT

WAITFORKEYPRESS()

CLOSECONSOLE
END

'----------------------------------------------------------------------------------
'END PROGRAM
'----------------------------------------------------------------------------------
'START MATRIX INVERSION FUNCTIONS
'----------------------------------------------------------------------------------

SUB INVERTMAT(INVERT[] AS DOUBLE, INVERSION[] AS DOUBLE)
'Matrix inversion by the Gauss-Jordan method

INT ROW, COL
DOUBLE DIVISOR, MULTIPLIER

SETPIVOTS()

'Scale each row of AM, so that 1 appears on the pivot (AM[COL, COL]).

FOR COL = 1 TO MS
DIVISOR = GETARRAYVAL(INVERT, MS, MS, COL, COL)
IF DIVISOR = 0 THEN SINGULARMATRIX()
SCALEROWBYDIV(INVERT, COL, DIVISOR)
SCALEROWBYDIV(INVERSION, COL, DIVISOR)

'Add multiple of pivot row of AM, to rows below pivot (AM[COL, COL]), so
'the column below the pivot, is zeroed out.

FOR ROW = COL + 1 TO MS
MULTIPLIER = - GETARRAYVAL(INVERT, MS, MS, ROW, COL)
ADDROWTOROW(INVERT, ROW, COL, MULTIPLIER)
ADDROWTOROW(INVERSION, ROW, COL, MULTIPLIER)
NEXT ROW
NEXT COL

'Add multiple of pivot row of AM, to rows above pivot (AM[COL, COL]), so
'the column above the pivot, is zeroed out.

FOR COL = MS TO 2 STEP -1
FOR ROW = COL - 1 TO 1 STEP -1
MULTIPLIER = - GETARRAYVAL(INVERT, MS, MS, ROW, COL)
ADDROWTOROW(INVERT, ROW, COL, MULTIPLIER)
ADDROWTOROW(INVERSION, ROW, COL, MULTIPLIER)
NEXT ROW
NEXT COL

ENDSUB

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

SUB SETPIVOTS()
'This function is not strictly necessary. It switches the rows of AM,
'so that in each column, the value with the largest magnitude, appears
'in the pivot position (AM[I, I]). Doing this, makes the floating point
'calculations more accurate.

INT ROW, COL, SWITCHROW
DOUBLE TEST, MAX
FOR COL = 1 TO MS - 1
MAX = 0
FOR ROW = COL TO MS
TEST = ABS(GETARRAYVAL(AM, MS, MS, ROW, COL))
IF TEST > MAX
MAX = TEST
SWITCHROW = ROW
ENDIF
NEXT ROW
IF SWITCHROW > COL
SWITCHROWS(AM,COL,SWITCHROW)
SWITCHROWS(AI,COL,SWITCHROW)
ENDIF
NEXT COL
ENDSUB

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

SUB SINGULARMATRIX()
PRINT "MATRIX IS SINGULAR, AND HAS NO INVERSION."
PRINT
PRINT "PRESS (ALMOST) ANY KEY TO QUIT."
PRINT
WAITFORKEYPRESS()
CLOSECONSOLE
END
ENDSUB

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

SUB INITIALIZEMATRICES()
INT I, J

'Fill A with random DOUBLEs.

FOR I = 1 TO MS
FOR J = 1 TO MS
SETARRAYVAL(A, MS, MS, I, J, GETRANDOMVALUE())
NEXT J
NEXT I

'Set AM equal to A.

FOR I = 1 TO NN
AM[I] = A[I]
NEXT I

'Set AI equal to the identity matrix, "I".

FOR I = 1 TO MS
FOR J = 1 TO MS
SETARRAYVAL(AI, MS, MS, I, J, 0)
NEXT J
NEXT I

FOR I = 1 TO MS
SETARRAYVAL(AI, MS, MS, I, I, 1)
NEXT I

'Set all values of the check vector, equal to 1.

FOR I =1 TO MS
V1[I] = 1
NEXT I

ENDSUB

'----------------------------------------------------------------------------------
'END MATRIX INVERSION FUNCTIONS
'----------------------------------------------------------------------------------
'START CHECK VECTOR FUNCTIONS
'----------------------------------------------------------------------------------

SUB CALCCHECKVECTOR()

MATMATMULT(A, AI, ID)
MATVECMULT(ID, V1, V2)
ENDSUB

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

SUB WRITECHECKVECTOR()
INT I
STRING S1, S2, S3, S4
S2 = "    "
STRING T  = USING("#######.###",ET)

STRING S5 = APPEND$("'Inversion of", STR$(MATSIZE), " x", STR$(MATSIZE), " matrix,")
STRING S6 = APPEND$("'elapsed time = ", T, " seconds")


OPENFILE(OUTFILE,FILESTRING,"W")

WRITE OUTFILE, S5
WRITE OUTFILE, "'of uniformly distributed random DOUBLEs,"
WRITE OUTFILE, "'in the range of (-1000, 1000)."
WRITE OUTFILE, ""

WRITE OUTFILE, S6
WRITE OUTFILE, ""

WRITE OUTFILE, "'If the inversion was calculated correctly, then every"
WRITE OUTFILE, "'value shown below should be very close to 1."
WRITE OUTFILE, ""

FOR I = 1 TO MS
S1 = USING("0####",I)
S3 = USING("#.#################", V2[I])

S4 = APPEND$(S1, S2, S3)
WRITE OUTFILE, S4
NEXT I
CLOSEFILE(OUTFILE)
ENDSUB

'----------------------------------------------------------------------------------
'END CHECK VECTOR FUNCTIONS
'----------------------------------------------------------------------------------
'START MATRIX FUNCTIONS
'----------------------------------------------------------------------------------

SUB SWITCHROWS(M[] AS DOUBLE, INT R1, INT R2)
'interchanges two rows

INT I
DOUBLE TEMP
FOR I = 1 TO MS
TEMP = GETARRAYVAL(M, MS, MS, R1, I)
SETARRAYVAL(M, MS, MS, R1, I, GETARRAYVAL(M, MS, MS, R2, I))
SETARRAYVAL(M, MS, MS, R2, I, TEMP)
NEXT I
ENDSUB

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

SUB ADDROWTOROW(M[] AS DOUBLE, INT TOROW, INT ADDROW, DOUBLE ADDROWMUL)
'adds a multiple of one row to another

INT I
DOUBLE PROD
FOR I = 1 TO MS
PROD = GETARRAYVAL(M, MS, MS, ADDROW, I) * ADDROWMUL
ADDTOARRAYVAL(M, MS, MS, TOROW, I, PROD)
NEXT I
ENDSUB

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

SUB SCALEROWBYDIV(M[] AS DOUBLE, INT ROW, DOUBLE DIVISOR)
'scales a row by division

INT I
FOR I = 1 TO MS
DIVIDEARRAYVAL(M, MS, MS, ROW, I, DIVISOR)
NEXT I
ENDSUB

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

SUB MATMATMULT(M1[] AS DOUBLE, M2[] AS DOUBLE, M3[] AS DOUBLE)
'multiplication of two matrices

INT I, J, K
DOUBLE SUM
FOR I = 1 TO MS
FOR J = 1 TO MS
SUM = 0
FOR K = 1 TO MS
SUM += GETARRAYVAL(M1, MS, MS, I, K) * GETARRAYVAL(M2, MS, MS, K, J)
NEXT K
SETARRAYVAL(M3, MS, MS, I, J, SUM)
NEXT J
NEXT I
ENDSUB

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

SUB MATVECMULT(M[] AS DOUBLE, VEC1[] AS DOUBLE, VEC2[] AS DOUBLE)
'multiplication of a matrix and a vector

INT I, J
DOUBLE SUM
FOR I = 1 TO MS
SUM = 0
FOR J = 1 TO MS
SUM += GETARRAYVAL(M, MS, MS, I, J) * VEC1[J]
NEXT J
VEC2[I] = SUM
NEXT I
ENDSUB

'----------------------------------------------------------------------------------
'END MATRIX FUNCTIONS
'----------------------------------------------------------------------------------
'START 2D TO 1D ARRAY INDEXING CONVERSION FUNCTIONS
'----------------------------------------------------------------------------------
'Unlike most of the other functions in this program, these functions do
'not depend on any global constants, for instance, "MS". All data are
'passed into the functions as parameters. I have mixed feelings about
'using the global constants in the other functions. Since I made the
'global arrays static, I needed the global constants, in order to
'dimension the arrays. Therefore, it seemed redundant to pass them as
'parameters. However, for the functions in this section, it was easy to
'pass in all required data, and therefore make them generally usable. I
'admit that most of the other functions, will not work for the general
'cases.
'----------------------------------------------------------------------------------

SUB GETINDEX(INT NROWS, INT NCOLS, INT ROW, INT COL), INT
'assumes array indices begin at 1

'Check for out of range.
IF ROW > NROWS THEN RETURN 0
IF COL > NCOLS THEN RETURN 0

RETURN (ROW - 1) * NCOLS + COL
ENDSUB

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

SUB SETARRAYVAL(M[] AS DOUBLE, INT NROWS, INT NCOLS, INT ROW, INT COL, DOUBLE V)
INT INDEX
INDEX = GETINDEX(NROWS, NCOLS, ROW, COL)
M[INDEX] = V
ENDSUB

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

SUB ADDTOARRAYVAL(M[] AS DOUBLE, INT NROWS, INT NCOLS, INT ROW, INT COL, DOUBLE V)
INT INDEX
INDEX = GETINDEX(NROWS, NCOLS, ROW, COL)
M[INDEX] += V
ENDSUB

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

SUB DIVIDEARRAYVAL(M[] AS DOUBLE, INT NROWS, INT NCOLS, INT ROW, INT COL, DOUBLE V)
INT INDEX
INDEX = GETINDEX(NROWS, NCOLS, ROW, COL)
M[INDEX] /= V
ENDSUB

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

SUB GETARRAYVAL(M[] AS DOUBLE, INT NROWS, INT NCOLS, INT ROW, INT COL), DOUBLE
INT INDEX
INDEX = GETINDEX(NROWS, NCOLS, ROW, COL)
RETURN M[INDEX]
ENDSUB

'----------------------------------------------------------------------------------
'END 2D TO 1D ARRAY INDEXING CONVERSION FUNCTIONS
'----------------------------------------------------------------------------------
'START WINDOWS TIMER FUNCTIONS
'----------------------------------------------------------------------------------

SUB STARTTIMER1()
_GETLOCALTIME(ST0)
ENDSUB

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

SUB STOPTIMER1()
_GETLOCALTIME(ST1)
ENDSUB

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

SUB CALCELAPSEDTIME()
'returns elapsed time in seconds

CHAR BORROW = FALSE
WORD DAY0 = ST0.WDAY
WORD HOU0 = ST0.WHOUR
WORD MIN0 = ST0.WMINUTE
WORD SEC0 = ST0.WSECOND
WORD MIL0 = ST0.WMILLISECONDS
WORD DAY1 = ST1.WDAY
WORD HOU1 = ST1.WHOUR
WORD MIN1 = ST1.WMINUTE
WORD SEC1 = ST1.WSECOND
WORD MIL1 = ST1.WMILLISECONDS
WORD EHOU, EMIN, ESEC, EMIL

IF MIL0 > MIL1
BORROW = TRUE
MIL1 += 1000
ENDIF

EMIL = MIL1 - MIL0

IF BORROW = TRUE
SEC1 -= 1
BORROW = FALSE
ENDIF

IF SEC0 > SEC1
BORROW = TRUE
SEC1 += 60
ENDIF

ESEC = SEC1 - SEC0

IF BORROW = TRUE
MIN1 -= 1
BORROW = FALSE
ENDIF

IF MIN0 > MIN1
BORROW = TRUE
MIN1 += 60
ENDIF

EMIN = MIN1 - MIN0

IF BORROW = TRUE
HOU1 -= 1
BORROW = FALSE
ENDIF

IF HOU0 > HOU1
BORROW = TRUE
HOU1 += 24
ENDIF

EHOU = HOU1 - HOU0

ET = EHOU * 3600 + EMIN * 60 + ESEC + EMIL / 1000.0

ENDSUB

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

SUB SHOWELAPSEDTIME()

STRING T  = USING("#######.###",ET)
STRING TS = APPEND$("ELAPSED TIME = ",T," SECONDS")
PRINT TS
ENDSUB

'----------------------------------------------------------------------------------
'END WINDOWS TIMER FUNCTIONS
'----------------------------------------------------------------------------------
'START RANDOM NUMBER GENERATION FUNCTIONS
'----------------------------------------------------------------------------------

SUB SETRAN2SEED(INT I)
'for a random seed, set I negative.

IF I < 0
RAN2SEED = GETNUMSECONDS()
ELSE
RAN2SEED = I
ENDIF

INITIALRAN2SEEDVAL = RAN2SEED
IF INITIALRAN2SEEDVAL = 0 THEN INITIALRAN2SEEDVAL = 1

ENDSUB

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

SUB RAN2(), DOUBLE
'random double generator
'from "numerical Recipes in C", 2nd edition, p.282.
'0 < RAN2() < 1

DOUBLE AM   = 1.0 / IM1
DOUBLE RNMX = 1.0 - DEPS
INT J, K
DOUBLE TEMP

IF RAN2FLAG = 0
RAN2SEED = ABS(RAN2SEED)
IF RAN2SEED < 1 THEN RAN2SEED = 1
IDUM2 = RAN2SEED
FOR J = NTAB + 7 TO 0 STEP -1
K = RAN2SEED / IQ1
RAN2SEED = IA1 * (RAN2SEED - K * IQ1) - K * IR1
IF RAN2SEED < 0 THEN RAN2SEED += IM1
IF J < NTAB THEN IV[J] = RAN2SEED
NEXT J
IY = IV[0]
RAN2FLAG = 1
ENDIF

K = RAN2SEED / IQ1
RAN2SEED = IA1 * (RAN2SEED - K * IQ1) - K * IR1
IF RAN2SEED < 0 THEN RAN2SEED += IM1
K = IDUM2 / IQ2
IDUM2 = IA2 * (IDUM2 - K * IQ2) - K * IR2
IF IDUM2 < 0 THEN IDUM2 += IM2
J = IY / NDIV
IY = IV[J] - IDUM2
IV[J] = RAN2SEED
IF IY < 1 THEN IY += IMM1
TEMP = AM * IY
IF TEMP > RNMX
RETURN RNMX
ELSE
RETURN TEMP
ENDIF

'This statement never executes.
'Added, so compiler doesn't warn that a return value is expected.
RETURN 0

ENDSUB

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

SUB GETRANDOMVALUE(), DOUBLE
'returns uniformly distributed values, between -1000 and 1000

RETURN RAN2() * 2000 - 1000
ENDSUB

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

SUB GETNUMSECONDS(), INT
'returns total number of seconds, since midnight of current day

INT TOTSECS = 0
STRING TIMENOW = TIME$
STRING TEMP
TEMP = MID$(TIMENOW, 1, 2)
TOTSECS += VAL(TEMP) * 3600
TEMP = MID$(TIMENOW, 4, 2)
TOTSECS += VAL(TEMP) * 60
TEMP = MID$(TIMENOW, 7, 2)
TOTSECS += VAL(TEMP)
RETURN TOTSECS

ENDSUB

'----------------------------------------------------------------------------------
'END RANDOM NUMBER GENERATION FUNCTIONS
'----------------------------------------------------------------------------------
'START UTILITY FUNCTIONS
'----------------------------------------------------------------------------------

SUB WAITFORKEYPRESS()
WHILE TRUE
S=INKEY$
IF S > CHR$(0) THEN RETURN
ENDWHILE
ENDSUB

'----------------------------------------------------------------------------------
'END UTILITY FUNCTIONS
'----------------------------------------------------------------------------------