ð#Syntax10.Scn.FntÁÁMODULE MatrixL; (* JFS 05/01/99 Version 0.1 *)
This module provides a set of procedures to perform basic calculations on matrixes of longreals.
The result of any operation is returned in the R matrix, the other matrixes being operands.
BEWARE that as the matrixes are VAR they are altered by the procedure. This is particularly true
for GaussMod, Determinant, Invert and PseudoInvert PROCEDUREs where the inputs are changed
during the computations. The action of each PROCEDURE is described hereafter. Notice that the
BOOLEAN returned by the PROCEDUREs is an indication of the consistency of dimensions. If the
BOOLEAN is false it means that the dimensions of the matrixes doesn't match or that there is no
solution.
PROCEDURE ScalMult*(VAR M:ARRAY OF ARRAY OF LONGREAL; K:LONGREAL);
Return the matrix M time the scalar K: [M] <-K*[M]
PROCEDURE Add*(VAR M1,M2,R :ARRAY OF ARRAY OF LONGREAL):BOOLEAN;
Return the sum of M1 and M2 in R: [R] <-[M1]+[M2]
PROCEDURE Sub*(VAR M1,M2,R :ARRAY OF ARRAY OF LONGREAL):BOOLEAN;
Return the difference of M1 and M2 in R: [R] <-[M1]-[M2]
PROCEDURE Mult*(VAR M1,M2,R :ARRAY OF ARRAY OF LONGREAL):BOOLEAN;
Return the product of M1 and M2 in R: [R] <-[M1]*[M2]
PROCEDURE Transpose*(VAR M1,R : ARRAY OF ARRAY OF LONGREAL):BOOLEAN;
Return the transpose of M1 in R: [R] <-[M1]t
PROCEDURE Unity*(VAR R:ARRAY OF ARRAY OF LONGREAL):BOOLEAN;
Return an unity matrix in R: [R] <-[I]
PROCEDURE Display*(VAR M:ARRAY OF ARRAY OF LONGREAL;Acc:INTEGER);
Provides a basic matrix display capability. Each element of R is displayed with a field of Acc
figures as specified in Out.LongReal.
PROCEDURE Find(VAR M:ARRAY OF ARRAY OF LONGREAL; Rinit:LONGINT; VAR Lswap:LONGINT);
This is an auxiliary procedure to get the first line with a pivot # 0 after line Rinit. The result is
in Lswap. It is not intended to be exported from the module.
PROCEDURE Swap(VAR M:ARRAY OF ARRAY OF LONGREAL;Up,Down:LONGINT);
This is an auxiliary procedure to swap lines Up and Down of Matrix M. It is not intended to be
exported from the module. There is a more efficient way of swapping lines by using an inter
mediate array of indexes but that needs either to pass this auxiliary array in the procedure
to allow for open arrays to be used or to implement a more complex algorithm using a
dynamic list. For simplicity sake I have used the direct implementation.
PROCEDURE GaussMod*(VAR A,R: ARRAY OF ARRAY OF LONGREAL; VAR Rank:INTEGER):BOOLEAN;
This procedure solves the system [A]*[X]=[R] and return the solution in [R]: [R] <-(1/[A])*[R]
BEWARE that A is changed at the end of the computation.
In that particular case the Boolean indicator is not only significant for dimensions consistency
but also for a solvable system. In other words if the returned value is FALSE either the dimensi
ons doesn't match OR there is no solution to the system. The number of free vectors in [A] is
indicated by the value of Rank. At the end of the computation either A is changed in an unity
matrix or only a number of free vectors equal to Rank has been found.
The algorithm used is a modified version of the Gauss one. Usually the Gauss algorithm works
in two steps:
first: triangulate the system, then: solve it.
Here the algorithm has been modified in order to get in one pass a direct solution by turning
A in an unity matrix rather a triangular one.
PROCEDURE Determinant*(VAR A:ARRAY OF ARRAY OF LONGREAL;VAR Det:LONGREAL;VAR Rank:LONGINT):BOOLEAN;
This procedure returns in Det the determinant of A. The boolean indicator and Rank have the
same meanings than for GaussMod. The alogrithm used is based on a triangulation of A and on
the fact that the determinant of a triangular matrix is the product of the elements on the dia
gonal.
BEWARE that A is changed at the end of the computation.
PROCEDURE Invert*(VAR A,R:ARRAY OF ARRAY OF LONGREAL;Rank:INTEGER):BOOLEAN;
This procedure returns the inverse of A in R: [R] <-1/[A]
BEWARE that A is changed at the end of the computation.
That's mainly the Jordan algorithm but modified as the GaussMod above. In fact Invert uses
GaussMod with [R] set to a unity matrix of appropriate dimensions.
Notice that there is a alternative to this algorithm that make the inversion of A on itself.
This approach saves some memory (what is not critical today) but needs to store during the
computations a list of lines having been swapped in order to perform the inverse swapping at
the end before returning the result. The amount of computations is very similar to the actual
implementation but needs to pass an auxiliary array to store the swaps with the same com
ments than for GaussMod. That's the main reason for which this solution has not been used.
PROCEDURE PseudoInvert*(VAR A,At,P,I,R:ARRAY OF ARRAY OF LONGREAL):BOOLEAN;
This procedure returns the pseudoinvert of A in R. The pseudo invert is defined as the
matrix minimizing the quadratic error on the system: [A]*[X]=[B] when they are more
equations than unknowns in X. In other words if PsI[A] is the so called pseudo inverse of [A]
then [X]=PsI[A] minimizes ||[B]-[A][X]|| where || || stands for the euclidian norm.
At,P,I are intermedate matrixes necessary for the computations.
If [A] is of dimensions m lines n rows with m>=n then:
[At] is n,m
[P] is m,m
[I] is m,m
[R] is n,m
When calling the procedure only A has to be defined, the other matrixes are undefined as
long as they are of dimensions as defined above. At the end, At contents the transpose of A,
P contents the product At*A, I contents the inverse of P and R contents the pseudo inverse of
A. A is unchanged.
PROCEDURE Copy*(VAR A,B:ARRAY OF ARRAY OF LONGREAL):BOOLEAN;
This procedure simply duplicate A in B: [A] <-[B]
END MatrixL.