(* ----------------------------------------------------------------------------
 * $Id: LINALGRN.mi,v 1.5 1995/11/05 09:24:14 kredel Exp $
 * ----------------------------------------------------------------------------
 * This file is part of MAS.
 * ----------------------------------------------------------------------------
 * Copyright (c) 1989 - 1992 Universitaet Passau
 * ----------------------------------------------------------------------------
 * $Log: LINALGRN.mi,v $
 * Revision 1.5  1995/11/05 09:24:14  kredel
 * Fixed comments and code.
 *
 * Revision 1.4  1992/10/15  16:29:16  kredel
 * Changed rcsid variable
 *
 * Revision 1.3  1992/08/21  13:52:13  kredel
 * Corrected error in RNMDET: wrong check for square matrix.
 *
 * Revision 1.2  1992/02/12  17:33:11  pesch
 * Moved CONST definition to the right place
 *
 * Revision 1.1  1992/01/22  15:12:51  kredel
 * Initial revision
 *
 * ----------------------------------------------------------------------------
 *)

IMPLEMENTATION MODULE LINALGRN;

(* MAS Linear Algebra Rational Number Implementation Module. *)


(*-------------------------------------------------------------------------

Implementation Module: Linear Algebra Rational Number 

Programmierpraktikum SS 1990: JUERGEN MUELLER, Heinz Kredel 

Contents: Standartoperations for matrix and vector manipulations,  
          Gaussian elimination, LU-Decomposition, Solve,    
          Inversion, Nullspace, Determinant. 

Stand : 27.09.90, Juergen Mueller                         
        from then date of the file, Heinz Kredel             

Remarks: A vector is represented as a list of elements.
         The elements may be integers or rational numbers.
         A matix is represented as a list of vectors.
         In most circumstances these vectors are interpreted row 
         vectors, but they can also be interpreted as column vectors.

--------------------------------------------------------------------------*)

  FROM MASELEM IMPORT GAMMAINT, MASMIN, MASEVEN;

  FROM MASBIOS IMPORT SWRITE, BLINES, CREADB, BKSP, MASORD, DIGIT;

  FROM MASSTOR IMPORT LIST, SIL, BETA, 
                      ADV, FIRST, RED, INV, LIST1, COMP, LENGTH;

  FROM MASERR IMPORT ERROR, severe;

  FROM SACLIST IMPORT LIST2, CCONC, CINV, CONC;

  FROM SACRN IMPORT RNINV, RNABS, RNCOMP, RNWRIT, RNNEG, RNDIF, RNQ, 
                    RNDEN, RNNUM, RNRED, RNINT, RNPROD, RNSUM;

  FROM MASRN IMPORT RNMAX, RNDWR;

  FROM MASAPF IMPORT RNDRD, APFRN, APWRIT;

CONST rcsidi = "$Id: LINALGRN.mi,v 1.5 1995/11/05 09:24:14 kredel Exp $";
CONST copyrighti = "Copyright (c) 1989 - 1992 Universitaet Passau";

 

TYPE F2 = PROCEDURE(LIST,LIST): LIST;

(*The TDI compiler accepts only global functions in expressions. *)   

PROCEDURE rnhi(i,j: LIST): LIST; (*Hilbert*)
BEGIN RETURN(RNRED(1,i+j-1)) END rnhi;

PROCEDURE rnum(i,j: LIST): LIST; (*unit matrix*)
BEGIN IF i = j THEN RETURN(RNINT(1)) ELSE RETURN(0) END;
(*9*) END rnum;

(*---------------- arbitrary matices --------------------*)

PROCEDURE MDIM(M : LIST): LIST;
(*Matrix dimension. M is a matrix. MDIM returns
max( row, column) of M. *)
VAR   zeile, spalte : LIST;
BEGIN
(*1*) zeile := LENGTH(FIRST(M));
      spalte := LENGTH(M);
      IF spalte < zeile THEN RETURN(zeile) 
                        ELSE RETURN(spalte) END;
(*9*) END MDIM;


PROCEDURE MGET(M, k, l : LIST): LIST;
(*Matrix get. M is a matrix. k, l are integers, 0 le k le  rows(M),
0 le l le  columns(M). MGET returns the element M(k,l) of matrix M. *)
VAR   i, j, h1, rest : LIST;
BEGIN
(*1*) rest := M; h1 := SIL;
      i := 0; j := 0;
(*2*) (*Suche nach Eintrag in der gewuenschter Zeile*)
      WHILE i < k DO
            i := i+1;
            ADV(rest, h1, rest);
            END;
      rest := h1;
(*3*) (*Suche nach Eintrag in der gewuenschten Spalte*)
      WHILE j < l DO
            j := j+1;
            ADV(rest, h1, rest);
            END;
      RETURN(h1);
(*9*) END MGET;


PROCEDURE MSET(M, k, l, x : LIST): LIST;
(*Matrix set. M is a matrix. k, l are integers, 0 le k le  rows(M),
0 le l le  columns(M). MSET sets the element M(k,l) to x. 
The new matrix is returned. *)
VAR    i, j, h1, rest1, rest2, neuk, neul, neuhilf, 
       egal, neu : LIST;
BEGIN
(*1*) rest1 := M; neuk := SIL;
      i := 0;
      j := 0;
(*2*) (*Suche k-te Zeile*)
      WHILE i < k DO
            i := i+1;
            ADV(rest1, h1, rest1);
            neuk := COMP(h1,neuk);
            END;
      rest2 := h1;
      neul := SIL;
(*3*) (*Suche l-te Spalte*)
      WHILE j < l DO
            j := j+1;
            ADV(rest2, h1, rest2);
            neul := COMP(h1,neul);
            END;
      neu := COMP(x,rest2); (*tausche x aus*)
(*4*) (*und fuege neu zusammen*)
      ADV(neul, egal, neul);
      WHILE neul <> SIL DO
            ADV(neul, neuhilf, neul);
            neu := COMP(neuhilf,neu);
            END;
      neu := LIST1(neu);
(*5*) (*restauriere Matrix*) 
      ADV(neuk, egal, neuk);
      WHILE neuk<>SIL DO
            ADV(neuk, neuhilf, neuk);
            neu := COMP(neuhilf,neu);
            END;
      neu := INV(neu);
(*6*) WHILE rest1 <> SIL DO
            ADV(rest1, neuhilf, rest1);
            neu := COMP(neuhilf,neu);
            END;
      neu := INV(neu);
      RETURN(neu);
(*9*) END MSET;


PROCEDURE VDELEL(V, i : LIST): LIST;
(*Vector delete element. V is a vector. The i-th element of V 
is deleted. 0 le i le length(V). *)
VAR   U, VP, v, j : LIST;
BEGIN
(*1*) IF i <= 0 THEN RETURN(V) END;
      IF V = SIL THEN RETURN(V) END;
      VP := V; j := 0; U := SIL;
(*2*) REPEAT
             j := j+1;
             IF VP = SIL THEN RETURN(V); END;
             ADV(VP, v, VP);
             U := COMP(v,U);
             UNTIL j=i;
(*3*) U := RED(U); U := INV(U); U := CONC(U,VP);
      RETURN(U);
(*9*) END VDELEL;


PROCEDURE MDELCOL(M, i : LIST): LIST;
(*Matrix delete column. M is a vector of row vectors. In each 
row the i-th element is deleted, 0 le i le  columns(M). The new 
matrix is returned. *)
VAR   N, MP, V : LIST;
BEGIN
(*1*) IF i <= 0 THEN RETURN(M) END;
      IF M = SIL THEN RETURN(M) END;
      MP := M; N := SIL;
(*2*) WHILE MP # SIL DO
            ADV(MP, V, MP);
            V := VDELEL(V,i);
            N := COMP(V,N);
            END;
(*3*) N := INV(N);
      RETURN(N);
(*9*) END MDELCOL;


PROCEDURE MMINOR(M, i, j : LIST): LIST;
(*Matrix minor. M is a vector of row vectors. The i-th column, 
0 le i le  rows(M), and in each remaining row the j-th element,
0 le j le  columns(M), is deleted. *)
VAR   MP: LIST;
BEGIN
(*1*) MP:=VDELEL(M,j); MP:=MDELCOL(MP,i); 
      RETURN(MP);
(*9*) END MMINOR;


PROCEDURE MTRANS(M : LIST): LIST;
(*Matrix transpose. M is a matrix. The transposed matrix is returned. *)
VAR   A, B, C, MT, a, b : LIST;
BEGIN
(*1*) B:=M; MT:=SIL;
      WHILE B # SIL DO 
            A:=B; B:=SIL; C:=SIL; 
            WHILE A # SIL DO ADV(A,a,A); 
                  IF a # SIL THEN ADV(a,b,a); 
                     B:=COMP(a,B); C:=COMP(b,C) END;
                  END;
            C:=INV(C); 
            IF C # SIL THEN MT:=COMP(C,MT) END;
            B:=INV(B);
            END;
      MT:=INV(MT); RETURN(MT);
(*9*) END MTRANS;


PROCEDURE VEL(a, n : LIST): LIST;
(*Vector elements. A vector of length n with elements a is returned. *)
VAR   A, k : LIST;
BEGIN
(*1*) A := SIL; k := 1;
(*2*) WHILE k<=n DO k := k+1;
            A := COMP(a,A); END;
(*3*) RETURN(A);
(*9*) END VEL;


PROCEDURE MFILL(M, m, n: LIST): LIST;
(*Matrix fill. M is an upper triangular matrix. A (m x n) matrix 
with zeros in the lower triangular part is returned. *)
VAR   i, j, V, MP, np: LIST;
BEGIN
(*1*) MP:=SIL; i:=0;
(*2*) WHILE i < m DO i:=i+1; j:=0; 
            IF M # SIL THEN ADV(M,V,M); np:=n-LENGTH(V);  
                       ELSE V:=SIL; np:=n END;
            WHILE j < np DO j:=j+1; V:=COMP(0,V) END;
            MP:=COMP(V,MP);
            END;
(*4*) MP:=INV(MP); RETURN(MP);
(*9*) END MFILL;


PROCEDURE MRANG(U: LIST): LIST;
(*Matrix rang. U is an upper triangular matrix from a
LU-decomposition. The rang of U is returned. *)
VAR   r: LIST; 
BEGIN
(*1*) r:=LENGTH(U); 
      RETURN(r);
(*9*) END MRANG;


PROCEDURE MAT(f: F2; m, n: LIST): LIST;
(*Matrix. An (m x n) matrix with elements given by f(i,j) is returned. *)
VAR   i, j, V, M, v: LIST;
BEGIN
(*1*) M:=SIL;
      IF (m <= 0) OR (n <= 0) THEN RETURN(M) END;
(*2*) i:=0;
      WHILE i < m DO i:=i+1; j:=0; V:=SIL; 
            WHILE j < n DO j:=j+1;
                  v:=f(i,j); V:=COMP(v,V) END;
            V:=INV(V); M:=COMP(V,M);
            END;
(*4*) M:=INV(M); RETURN(M);
(*9*) END MAT;


(*---------------- rational number matices --------------------*)

PROCEDURE RNMHILBERT(m, n : LIST): LIST;
(*Rational number matrix Hilbert. m, n integer. A (m x n) rational 
number Hilbert matrix is returned. *)
BEGIN
(*1*) RETURN(MAT(rnhi,m,n));
(*9*) END RNMHILBERT;


PROCEDURE RNUM(m, n : LIST): LIST;
(*Rational number unit matrix. m, n integer. A (m x n) rational 
number unit matrix is returned. *)
BEGIN
(*1*) RETURN(MAT(rnum,m,n));
(*9*) END RNUM;


PROCEDURE RNVWRITE(A, s : LIST);
(*Rational number vector write. A is a rational number vector. 
A is written to the output stream. The rational numbers are written 
as rational numbers if s = -1, as decimal approximations, with s decimal 
digits if s >= 0 or in floating point format if s < -1. *)
VAR   a : LIST;
BEGIN
(*1*) SWRITE("(");
(*2*) WHILE A # SIL DO
            ADV(A, a, A);
            IF s < 0 THEN 
               IF s = -1 THEN RNWRIT(a) 
                  ELSE a:=APFRN(a); APWRIT(a); END; 
               ELSE RNDWR(a,s); END;
            IF A # SIL THEN SWRITE(", ") END;
            END;
(*3*) SWRITE(")");
(*9*) END RNVWRITE;


PROCEDURE RNVREAD(): LIST;
(*Rational number vector read. A rational number vector is 
read from the input stream, and returned. *)
VAR   c, a, A: LIST;
BEGIN
(*1*) c:=CREADB(); A:=SIL;
      IF c # MASORD("(") THEN ERROR(severe,"RNVREAD ( expected."); 
         BKSP; RETURN(A) END;
(*2*) LOOP c:=CREADB(); 
           IF c = MASORD(")") THEN EXIT END; BKSP; 
           a:=RNDRD(); A:=COMP(a,A);
           c:=CREADB();               
           IF c # MASORD(",") THEN BKSP END;
           END;
(*3*) A:=INV(A); RETURN(A); 
(*9*) END RNVREAD;


PROCEDURE RNMWRITE(A, s : LIST);
(*Rational number matrix write. A is a rational number matrix. 
A is written to the output stream. The rational numbers are written 
as rational numbers if s = -1, as decimal approximations, with s decimal 
digits if s >= 0 or in floating point format if s < -1. *)
VAR   a : LIST;
BEGIN
(*1*) BLINES(0); SWRITE("(");
(*2*) WHILE A # SIL DO
            ADV(A, a, A);
            RNVWRITE(a,s);
            IF A # SIL THEN SWRITE(", "); BLINES(0) END;
            END;
(*3*) SWRITE(")"); BLINES(0);
(*9*) END RNMWRITE;


PROCEDURE RNMREAD(): LIST;
(*Rational number matrix read. A rational number matrix is 
read from the input stream, and returned. *)
VAR   c, a, A: LIST;
BEGIN
(*1*) c:=CREADB(); A:=SIL;
      IF c # MASORD("(") THEN ERROR(severe,"RNMREAD ( expected."); 
         BKSP; RETURN(A) END;
(*2*) LOOP c:=CREADB(); 
           IF c = MASORD(")") THEN EXIT END;
           IF c = MASORD("(") THEN BKSP; a:=RNVREAD(); 
              A:=COMP(a,A);
              c:=CREADB();               
              IF c # MASORD(",") THEN BKSP END;
              END;
           END;
(*3*) A:=INV(A); RETURN(A); 
(*9*) END RNMREAD;


PROCEDURE RNVFIV(A : LIST): LIST;
(*Rational number vector from integer vector. A is an integer 
vector. A rational number vector with denominators 1 
and nominators equal to the elements of A is returned. *)
VAR   a, c, C : LIST;
BEGIN
(*1*) C := SIL; a := SIL; c := SIL;
(*2*) WHILE A # SIL DO
            ADV(A, a, A);
            c := RNINT(a);
            C := COMP(c,C);
            END;
(*2*) C := INV(C);
      RETURN(C);
(*2*) END RNVFIV;


PROCEDURE RNMFIM(M : LIST): LIST;
(*Rational number matrix from integer matrix. A is an integer 
matrix. A rational number matrix with denominators 1 
and nominators equal to the elements of A is returned. *)
VAR   hilf, neu : LIST;
BEGIN
(*1*) hilf := SIL;
(*2*) WHILE M # SIL DO
            ADV(M, neu, M);
            neu := RNVFIV(neu);
            hilf := COMP(neu,hilf);
            END;
(*3*) hilf := INV(hilf);
      RETURN(hilf);
(*9*) END RNMFIM;


PROCEDURE RNVDIF(A, B : LIST): LIST;
(*Rational number vector difference. A and B are rational number 
vectors. The rational number vector C = A - B is returned. *)
VAR   a, b, c, C, test : LIST;
BEGIN
(*1*) C := SIL;
      WHILE A # SIL DO
            ADV(A, a, A);
            ADV(B, b, B);
            c := RNDIF(a,b);
            C := COMP(c,C);
            END;
      C := INV(C); RETURN(C);
(*9*) END RNVDIF;


PROCEDURE RNVQ(A, B : LIST): LIST;
(*Rational number vector quotient. A and B are rational number vectors.
The rational number vector C = A / FIRST(B) is returned. *)
    VAR   a, b, c, C : LIST;
    BEGIN
      C := SIL;
      IF B#SIL THEN
        WHILE A#SIL DO
          ADV(A, a, A);
          c := RNQ(a,FIRST(B));
          C := COMP(c,C);
        END;
        C := INV(C);
        RETURN(C);
      ELSE
        RETURN(A);
      END;
    END RNVQ;


PROCEDURE RNVQF(A : LIST): LIST;
(*Rational number vector quotient. A is a rational number vector.  
The rational number vector C = A / FIRST(A) is returned. *)
    VAR   b : LIST;
    BEGIN
      IF b#LIST(0) THEN
        A := RNVQ(A,A);
      END;
      RETURN(A);
    END RNVQF;


PROCEDURE RNVVSUM(A, B : LIST): LIST;
(*Rational number vector vector sum. A and B are rational number vectors. 
A rational number vector C = A + B is returned. *)
VAR   a, b, c, C, test : LIST;
BEGIN
(*1*) C := SIL;
(*2*) WHILE A # SIL DO
            ADV(A, a, A);
            ADV(B, b, B);
            c := RNSUM(a,b);
            C := COMP(c,C);
            END;
(*3*) C := INV(C); RETURN(C);
(*9*) END RNVVSUM;


PROCEDURE RNVSVSUM(A, B : LIST): LIST;
(*Rational number vector scalar sum. A and B are rational number vectors.
A rational number vector C = A + FIRST(B) is returned. *)
VAR   a, b, c, C : LIST;
BEGIN
(*1*) C := SIL; b:=FIRST(B); 
(*2*) WHILE A # SIL DO
            ADV(A, a, A);
            c := RNSUM(a,b);
            C := COMP(c,C);
            END;
(*3*) C := INV(C); RETURN(C);
(*9*) END RNVSVSUM;


PROCEDURE RNVSSUM(A : LIST): LIST;
(*Rational number vector scalar sum. A is a rational number vector.
A rational number C = a1 + a2 + ... + an is returned. *)
VAR   a, c : LIST;
BEGIN
(*1*) a := SIL; c := 0;
(*2*) WHILE A # SIL DO
            ADV(A, a, A);
            c := RNSUM(a,c);
            END;
      RETURN(c);
(*2*) END RNVSSUM;


PROCEDURE RNVSVPROD(A, B : LIST): LIST;
(*Rational number vector scalar vector product. A and B are 
rational number vectors. A rational number vector 
C = (a1*FIRST(B), ..., an*FIRST(B)) is returned. *)
VAR   a, b, c, C : LIST;
BEGIN
(*1*) C := SIL; b:=FIRST(B);
(*2*) WHILE A # SIL DO
            ADV(A, a, A);
            c := RNPROD(a,b);
            C := COMP(c,C);
            END;
(*3*) C := INV(C); RETURN(C);
(*9*) END RNVSVPROD;


PROCEDURE RNVVPROD(A, B : LIST): LIST;
(*Rational number vector vector product. A and B are 
rational number vectors. A rational number vector 
C = (a1*b1, ..., an*bn) is returned. *)
VAR   a, b, c, C : LIST;
BEGIN
(*1*) C := SIL;
(*2*) WHILE A # SIL DO
            ADV(A, a, A);
            ADV(B, b, B);
            c := RNPROD(a,b);
            C := COMP(c,C);
            END;
(*3*) C := INV(C); RETURN(C);
(*9*) END RNVVPROD;


PROCEDURE RNVSPROD(A, B : LIST): LIST;
(*Rational number vector scalar product. A and B are rational 
number vectors. A rational number C = a1*b1 + ... + an*bn is
returned. *)
VAR   a, b, c, C : LIST;
BEGIN
(*1*) C := 0;
(*2*) WHILE A # SIL DO
            ADV(A, a, A);
            ADV(B, b, B);
            c := RNPROD(a,b);
            C := RNSUM(c,C);
            END;
(*3*) RETURN(C);
(*9*) END RNVSPROD;


PROCEDURE RNVMAX(M : LIST): LIST;
(*Rational number vector maximum norm. M is a rational number 
vector. A rational number a = maximum absolute value M(i) 
is returned. *)
VAR   g, zeile, element, max : LIST;
BEGIN
(*1*) max := 0;
(*2*) WHILE M # SIL DO
            ADV(M, element, M);
            element := RNABS(element);
            max:=RNMAX(max,element);  
            END;
(*3*) RETURN(max);
(*9*) END RNVMAX;


PROCEDURE RNMSUM(A, B : LIST): LIST;
(*Rational number matrix sum. A and B are rational number
matrices. A rational number matrix C = A + B is returned. *)
VAR   a, b, c, C : LIST;
BEGIN
(*1*) IF A=SIL THEN RETURN(B) END;
      IF B=SIL THEN RETURN(A) END;
      C := SIL;
(*2*) WHILE A#SIL DO
            ADV(A, a, A);
            ADV(B, b, B);
            c := RNVVSUM(a,b);
            C := COMP(c,C);
            END;
(*3*) C := INV(C); RETURN(C);
(*9*) END RNMSUM;


PROCEDURE RNMDIF(A, B : LIST): LIST;
(*Rational number matrix difference. A and B are rational number
matrices. A rational number matrix C = A - B is returned. *)
VAR   a, b, c, C : LIST;
BEGIN
(*1*) IF A=SIL THEN
         B := RNSMPROD(B,RNINT(-1));
         RETURN(B);
         END;
      IF B=SIL THEN RETURN(A); END;
      C := SIL;
(*2*) WHILE A # SIL DO
            ADV(A, a, A);
            ADV(B, b, B);
            c := RNVDIF(a,b);
            C := COMP(c,C);
            END;
(*3*) C := INV(C); RETURN(C);
(*9*) END RNMDIF;


PROCEDURE RNSMPROD(A, B : LIST): LIST;
(*Rational number scalar and matrix product. A is a rational 
number matrix. B is a rational number. A rational number matrix 
C = A * B is returned. *)
VAR   a, b, c, C : LIST;
BEGIN
(*1*) IF A=SIL THEN RETURN(A) END;
      C := SIL; b:=LIST1(B);
(*2*) WHILE A # SIL DO
            ADV(A, a, A);
            c := RNVSVPROD(a,b); (* geaendert *)
            C := COMP(c,C);
            END;
(*3*) C := INV(C); RETURN(C);
(*9*) END RNSMPROD;


PROCEDURE RNMMAX(M : LIST): LIST;
(*Rational number matrix maximum norm. M is a rational number 
matrix. A rational number a = maximum absolute value M(i,j) 
is returned. *)
VAR   g, zeile, element, max : LIST;
BEGIN
(*1*) max := 0;
(*2*) WHILE M # SIL DO
            ADV(M, zeile, M);
            element := RNVMAX(zeile);
            max:=RNMAX(max,element);  
            END;
(*3*) RETURN(max);
(*9*) END RNMMAX;


PROCEDURE RNMGE(M : LIST): LIST;
(*Rational number matrix Gaussian elimination. M is a (n x m) 
rational number matrix. A (n x m) rational number matrix
resulting from Gaussian elimination is returned.  
RNMGELUD is called. *)
VAR   l, u, i, j, MP : LIST;
BEGIN
(*1*) i:=LENGTH(M);
      IF M # SIL THEN j:=LENGTH(FIRST(M))
                 ELSE j:=0 END;
(*2*) RNMGELUD(M,l,u);
(*3*) MP:=MFILL(u,i,j);
      RETURN(MP);
(*9*) END RNMGE;


PROCEDURE RNMDETL(M : LIST): LIST;
(*Rational number matrix determinant, using Laplace expansion. 
M is a rational number matrix. The determinant of M is returned. *)
VAR   i, d, dp, s, N, MP, V, VP, v : LIST;
BEGIN
(*1*) d := 0;
      IF M = SIL THEN RETURN(d); END;
      ADV(M, V, MP);
      IF MP = SIL THEN (*not square matrix ? *)
         IF RED(V) # SIL THEN d:=0 ELSE d:=FIRST(V) END;
         RETURN(d); END;
      s := 1;
      i := 0;
(*2*) WHILE V # SIL DO
            ADV(V, v, V);
            i := i+1;
            IF v # 0 THEN
               N := MDELCOL(MP,i);
               dp := RNMDETL(N);
               dp := RNPROD(v,dp);
               IF s < 0 THEN dp := RNNEG(dp) END;
               d := RNSUM(d,dp);
               END;
            s := -s; END;
      RETURN(d);
(*9*) END RNMDETL;


(*-------------------*)


PROCEDURE RNVLC(a, A, b, B : LIST): LIST;
(*Rational number vector linear combination. A and B are rational 
number vectors. a and b are rational numbers. A rational number vector 
C = a*A + b*B is returned. *)
VAR    c, cp, C, ap, bp : LIST;
BEGIN
(*1*) C := SIL;
(*2*) WHILE A#SIL DO
            ADV(A, ap, A);  ADV(B, bp, B);
            c := RNPROD(a,ap); cp := RNPROD(b,bp);
            c := RNSUM(c,cp); C := COMP(c,C);
            END;
(*3*) C := INV(C); RETURN(C);
(*9*) END RNVLC;


PROCEDURE RNSVPROD(a, A : LIST): LIST;
(*Rational number vector product with scalar. A is a rational 
number vector. a is a rational number. A rational number vector
C = a*A is returned. *)
VAR   c, cp, C, ap : LIST;
BEGIN
(*1*) C := SIL;
(*2*) WHILE A#SIL DO
            ADV(A, ap, A);
            c := RNPROD(a,ap);
            C := COMP(c,C);
            END;
(*3*) C := INV(C); RETURN(C);
(*9*) END RNSVPROD;


PROCEDURE RNMPROD(A, B : LIST): LIST;
(*Rational number matrix product. A and B are rational number
matrices. A rational number matrix C = A * B is returned, if 
the number of coloums of A is equal to the number of rows of B, 
otherwise the empty matrix is returned. *)
VAR   H1, H2, C, C1, a, c, BP, BT, b : LIST;
BEGIN
(*1*) C := SIL;
      IF A = SIL THEN RETURN(C) END;
      IF B = SIL THEN RETURN(C) END;
      H1 := LENGTH(FIRST(A));
      H2 := LENGTH(B);
      IF H1 # H2 THEN RETURN(C) END;
(*2*) BT:=MTRANS(B); (*transpose B. *)
      WHILE A # SIL DO ADV(A, a, A);
            C1:=SIL; BP:=BT;
            WHILE BP # SIL DO ADV(BP,b,BP);
                  c:=RNVSPROD(a,b);   
                  C1:=COMP(c,C1);
                  END;
            C1 := INV(C1);
            C := COMP(C1,C);
            END;
(*4*) C := INV(C); RETURN(C);
(*9*) END RNMPROD;


(*--------------------------*)


PROCEDURE RNMDET(M : LIST): LIST;
(*Rational number matrix determinant, using Gaussian elimination. 
M is a rational number matrix. The determinant of M is returned. *)
VAR   d, s, i, C, M1, MP, a, b, c, MP1, V, e: LIST;
BEGIN
(*1*) M1 := M; d := 0; s := 1;
      IF M1 = SIL THEN RETURN(d) END;
      V:=FIRST(M1); 
      IF V = SIL THEN RETURN(d) END;
      IF LENGTH(M1) # LENGTH(V) THEN RETURN(d) END; (* not square *) 
      e:=RNINT(1); d := e;
(*2*) WHILE M1 # SIL DO
            MP1 := M1; M1 := SIL;
            a := 0; i:=0;
            (*search pivot row and count sign changes. *)
            WHILE (MP1 # SIL) AND (a = 0) DO
                  ADV(MP1, C, MP1); i:=i+1;
                  ADV(C, a, C);
                  IF a = 0 THEN M1 := COMP(C,M1) END;
                  END;
            IF MASEVEN(i) THEN s:=-s END;
            IF a = 0 THEN d:=0; RETURN(d); END;
            (*rest matrix. *)
            d:=RNPROD(d,a);
            WHILE MP1 # SIL DO
                  ADV(MP1, V, MP1);
                  ADV(V, b, V);
                  IF b # 0 THEN b:=RNNEG(RNQ(b,a));
                     V := RNVLC(e,V,b,C) END;
                  M1 := COMP(V,M1);
                  END;
            M1 := INV(M1); 
            IF M1 # SIL THEN (*not square matrix ? *)
               IF FIRST(M1) = SIL THEN d := 0; RETURN(d) END;
               END;
            END;
(*3*) IF s < 0 THEN d:=RNNEG(d) END;
      RETURN(d);
(*9*) END RNMDET;


PROCEDURE RNMGELUD(M : LIST; VAR L, U: LIST);
(*Rational number matrix Gaussian elimination LU-decomposition. 
M is a rational number matrix represented rowwise. L is a lower 
triangular rational number matrix represented columnwise.
U is an upper triangular rational number matrix represented rowwise.
M = L * U for appropriate modifications of L and U. 
The pivot operations are also recorded in L. *)
VAR   C, M1, MP, L1, a, b, c, MP1, V, e, n, Z : LIST;
BEGIN
(*1*) M1 := M; MP := SIL; L := SIL; U := SIL;
      IF M1 = SIL THEN RETURN END;
      IF FIRST(M1) = SIL THEN RETURN END;
      e := RNINT(1);
      n := RNINT(0);
(*2*) WHILE M1 # SIL DO
            L1 := SIL; MP1 := M1; M1 := SIL;
            a := 0;
            (*search pivot row. *)
            WHILE (MP1 # SIL) AND (a = 0) DO
                  ADV(MP1, C, MP1);
                  ADV(C, a, C);
                  IF a = 0 THEN M1 := COMP(C,M1); L1 := COMP(a,L1); END;
                  END;
            IF a = 0 THEN (*C:=COMP(a,C); M1:=RED(M1); *)
               ELSE  
               a := RNQ(e,a);
               L1 := COMP(a,L1);
               C := RNSVPROD(a,C);
               WHILE MP1 # SIL DO
                     ADV(MP1, V, MP1);
                     ADV(V, b, V);
                     IF b # 0 THEN V := RNVLC(e,V,RNNEG(b),C) END;
                     M1 := COMP(V,M1);
                     L1 := COMP(b,L1);
                     END;
               C := COMP(e,C);
               MP := COMP(C,MP); 
               END;
            M1 := INV(M1); 
            L1 := INV(L1); L := COMP(L1,L);
            IF M1 # SIL THEN (*if not square matrix *)
               IF FIRST(M1) = SIL THEN M1 := SIL END;
               END;
            END;
(*3*) U := INV(MP); L := INV(L);
(*9*) END RNMGELUD;


PROCEDURE RNMLT(L, b : LIST): LIST;
(*Rational matrix lower triangular matrix transformation. 
L is a lower triangular rational number matrix represented 
columnwise as generated by RNMGELUD. b is a rational number 
vector. A rational number vector u = L * b is returned, 
such that if M * x = b and M = L * U, then U * x = u. *)
VAR   u, L1, b1, c, l, d, a : LIST;
BEGIN
(*1*) u := SIL;
(*2*) WHILE L # SIL DO
            ADV(L, L1, L);
            b1 := b; b := SIL;
            l := 0;
            WHILE (L1 # SIL) AND (l = 0) DO
                  ADV(L1, l, L1);
                  ADV(b1, c, b1);
                  IF l = 0 THEN b := COMP(c,b) END;
                  END;
            IF l = 0 
               THEN (*u := COMP(c,u);*) 
               ELSE (* l # 0 *)
                    c := RNPROD(c,l);
                    u := COMP(c,u);
                    WHILE L1 # SIL DO
                          ADV(L1, l, L1);
                          ADV(b1, d, b1);
                          a := RNPROD(l,c);
                          a := RNDIF(d,a);
                          b := COMP(a,b);
                          END;
                    END;
            b := INV(b);
            END;
(*3*) WHILE b <> SIL DO ADV(b,c,b); (*unsolvable ? *)
            IF c # 0 THEN RETURN(SIL) END;
            END;
(*4*) u := INV(u); RETURN(u);
(*9*) END RNMLT;


PROCEDURE RNMUT(U, b : LIST): LIST;
(*Rational matrix upper triangular matrix transformation. 
U is an upper triangular rational number matrix represented rowwise
as generated by RNMGELUD. b is a rational number vector 
b = L * b' as generated by RNMLT. A rational number vector x, 
such that U * x = b is returned. If no such x exists, then an 
empty vector is returned. If more than one such x exists, then 
for free x(i), x(i) = 0 is taken. *)
VAR   UP, x, U1, a, c, y, u, xp, b1: LIST;
BEGIN
(*1*) x:=SIL; 
      IF U = SIL THEN RETURN(x) END;
      UP := CINV(U); b1 := CINV(b);
(*2*) WHILE UP # SIL DO
            ADV(UP, U1, UP);
            U1 := CINV(U1);
            y := 0;
            xp := x;
            WHILE xp # SIL DO
                  ADV(xp, c, xp);
                  ADV(U1, u, U1);
                  u := RNPROD(u,c);
                  y := RNSUM(y,u);
                  END;
            xp:=SIL; 
            WHILE RED(U1) # SIL DO
                  U1:=RED(U1);
                  xp:=COMP(0,xp); (* arbitrary solution. *)
                  END;
            ADV(U1, u, U1); 
            ADV(b1, c, b1);
            c := RNDIF(c,y);
            (* u = 0, should not occur. *) 
            a := RNQ(c,u); 
            xp := COMP(a,xp); xp:=INV(xp);
            x := CONC(x,xp);
            END;
(*3*) x := INV(x); RETURN(x);
(*9*) END RNMUT;


PROCEDURE RNMSDS(L, U, b : LIST): LIST;
(*Rational number matrix solve decomposed system. 
L is a lower triangular rational number matrix represented 
columnwise, U is an upper triangular rational number matrix 
represented rowwise. L and U as generated by RNMGELUD.
If M = L * U, then a rational number vector x, such that 
M * x = b is returned. If no such x exists, then an empty 
vector is returned. If more than one such x exists, then 
for free x(i), x(i) = 0 is taken. *)
VAR   u, x : LIST;
BEGIN
(*1*) u := RNMLT(L,b);
      IF u # SIL THEN x := RNMUT(U,u);
                 ELSE x := SIL; (*unsolvable*) END;
(*3*) RETURN(x);
(*9*) END RNMSDS;


PROCEDURE RNMINV(A : LIST): LIST;
(*Rational number matrix inversion. A is a rational number matrix 
represented rowwise. If it exists, the inverse matrix of A is 
returned, otherwise an empty matrix is returned. *)
VAR   B, l, u, n, E, e, x : LIST;
BEGIN
(*1*) n := LENGTH(A); B := SIL;
      IF n <= 0 THEN RETURN(B) END;
      IF n # LENGTH(FIRST(A)) THEN RETURN(B) END;
(*2*) (*LU-decomposition. *)
      RNMGELUD(A,l,u); 
(*3*) (*solve A Ainv = E. *)
      E := RNUM(n,n);
      WHILE E # SIL DO
            ADV(E, e, E);
            x := RNMSDS(l,u,e);
            IF x = SIL THEN RETURN(SIL) END; (*singular *)
            B := COMP(x,B);
            END;
(*4*) B := INV(B); RETURN(B);
(*9*) END RNMINV;


PROCEDURE RNMUNS(U : LIST): LIST;
(*Rational number matrix upper triangular matrix solution null space. 
U is an upper triangular rational number matrix represented rowwise
as generated by RNMGELUD. A matrix X of linear independent rational 
number vectors x is returned, such that for each x in X, U * x = 0 holds. 
If only x = 0 satisfies the condition U * x = 0, then the 
matrix X is empty. *)
VAR   UP, V, N, Z, ZP, e, a, c, u, up, b, i, n, m, j: LIST; 
BEGIN
(*1*) UP:=CINV(U); N:=SIL; Z:=SIL; i:=0; e:=RNINT(1);
(*2*) WHILE UP # SIL DO ADV(UP,u,UP); 
            n:=LENGTH(u)-1;
            IF i < n THEN ADV(u,b,u); 
               WHILE i < n DO i:=i+1; j:=i; 
                     ZP:=COMP(e,Z); Z:=COMP(0,Z); 
                     WHILE j < n DO j:=j+1; ZP:=COMP(0,ZP) END;
                     a:=RNVSPROD(u,ZP); 
                     c:=RNNEG(a); a:=RNQ(c,b); 
                     j:=j+1; ZP:=COMP(a,ZP);
                     V:=UP; 
                     WHILE V # SIL DO ADV(V,up,V);
                           m:=LENGTH(up)-1; 
                           WHILE j < m DO j:=j+1; ZP:=COMP(0,ZP) END;
                           ADV(up,c,up); 
                           a:=RNVSPROD(up,ZP); 
                           a:=RNNEG(a); a:=RNQ(a,c); 
                           j:=j+1; ZP:=COMP(a,ZP);
                           END;
                     N:=COMP(ZP,N); 
                     END;
               END; 
            i:=n+1;    
            END;
      N:=INV(N); RETURN(N);
(*9*) END RNMUNS;


END LINALGRN.

(* -EOF- *)