EXPLANATION  FILE FOR PROGRAM SYSMAT
			  ====================================

	
	Program Sysmat solves a linear system AX = B or can only calculate the inverse
	of the square matrix A(N,N) or its determinant. This program (from BIBLI [13])
	is organized to minimize the used memory space.

	1. Numerical Method
           ----------------

	Let be A, square matrix of size N x N, and B, another matrix of size N x M.
	We want to solve equation A X = B, where X is a matrix of unknowns of size N x M
	(this allows simultaneously solving several right side data).
	
	To minimize the used memory space, when calculations terminate, the inverse of A
	will be stored in place of A, and the solution matrix X = A^-1 B in place of B.
	
	We use here the pivotal Gauss-Jordan method with choice of greatest available
	pivot at each step.
	
	The principle of that method is to make a series of transformations on A matrix,
	lambda1, lambda2,...,lambdan, so that A becomes the unity matrix.
	
	So after k-1 transformations, the A matrix becomes:
	
		| 1 0 ... 0 | a'1k ... a'1n   |
		| 0 1 ... 0 | a'2k ... a'2n   |
		| --------- | --------------- |
		| 0 0 ... 1 | a'k-1k.. a'k-1n |
		| --------------------------- | = lambda    lambda    ... lambda  A
		| 0 0 ... 0 | a'kk ... a'kn   |         k-1       k-2           1
		| 0 0 ... 0 | --------------- |
		| --------- | --------------- |
		| 0 0 ... 0 | a'nk ... a'nn   |  

	The inverse of A is then easily obtained by:

		Inv(A) = lambda  lambda    ... lambda  I    (I = unity matrix)
		               n       n-1           1

	and the solution matrix X is obtained by:

		X = lambda  lambda    ... lambda  B
		          n       n-1           1

                     (k)      (k) 
	Let us note a    and b    the elements of matrices A and B after the k th 
	             ij       ij
	
	transformation k=0 means initial A and B matrices).
	
	This lambda  transformation makes following operations:
                   k

                 (k)     (k-1)     (k-1)  (k-1)    (k-1)
	1)	a    =  a      -  a      a      / a
	         ij      ij        ik     kj       kk

		for  | i=1,...,k-1, k+1,...,n
		     | j=k,...,n
			 
                 (k)     (k-1)     (k-1)  (k-1)    (k-1)
	2)	b    =  b      -  a      b      / a
	         ij      ij        ik     kj       kk

		for  | i=1,...,k-1, k+1,...,n
		     | j=1,...,m

                 (k)     (k-1)
	3)	a    =  a
	         ij      ij

		for  | i=1,...,n
		     | j=1,...,k-1

                 (k)     (k-1)    (k-1)
	4)	a    =  a      / a
	         kj      kj       kk

		for  j=1,...,n

                 (k)     (k-1)    (k-1)
	5)	b    =  b      / a
	         kj      kj       kk

		for  j=1,...,m

        For each step we must divide by a pivot, a  ; To minimize round up errors,
                                                  kk
        at each step we choose as pivot the greatest element in sub-matrix:


                | a'   ... a'   |
                |   kk       kn |
                | ------------- |
                | a'   ... a'   |
                |  nk       nn  |

        and we make the necessary line and column permutations to bring this element
        in k th position on the main diagonal. Then we must "rearrange" the obtained
        matrix to take into account these permutations.


        2. Programming technique
           ---------------------

        So the calculation is made of n steps (k=1,...,n).

        As we want to store Inv(I) in place of A, and X in place of B, the used
        algorithm is rather complex and we advise the reader to check it with a
        simple example, for instance of size n=3.

        Let us note li and ci the lines and columns of A and B. At step k, we do the
        following operations:

        1) - choose a pivot: we take the greatest elemnt (in absolute value) a ij,
             where i,j = k,...,n.  Let be a      this pivot stored in variable PV.
                                           ik jk

           - update the determinant DET = DET * S * PV, where S depends on permutations
             to make:
                               S = Sl * Sc with  | Sl =  1 if ik = k
                                                 | Sl = -1 if ik <> k
                                                 |
                                                 | Sc =  1 if jk = k
                                                 | Sc = -1 if jk <> k

           - store ik and jk in two vectors Pl and Pc:   Pl(k)=ik, Pc(k)=jk

        2) Locate pivot to position (k,k) by:

                  a) exchange lines in A and B:  lk <--> l ik

                  b) exchange columns in A:  ck <--> c jk

        3) store column k of matrix A in vector Cs:

                 Cs(i) = a ik (i=1,...,k-1,k+1,...,n)
                 Cs(k) = 0

        4) replace column k of A by a ik = delta ik

                 delta ik = 1 if i=k, delta ik = 0 if i <> k

        5) replace lines k of A and B by l / PV
                                          k

        6) transform the other lines of A and B as follows:

                 l  --> l  - Cs(j) l    (j=1,...,k-1,k+1,...,n)
                  j      j          k

        7) rearrange A and B matrices to take into account the made permutations:

           - Exchange lines of A and B in order:

                      l    <-->  l
                       n          Pc(n)

                      l    <-->  l
                       n-1        Pc(n-1)

                      ---         ---

                      l    <-->  l
                       1          Pc(1)

           - Exchange columns of A in order:

                      c    <-->  c
                       n          Pl(n)

                      c    <-->  c
                       n-1        Pl(n-1)

                      ---         ---

                      c    <-->  c
                       1          Pl(1)


---------------------------------------------------------------------------------------
End of file sysmat.txt