DEFINITION MODULE MatLib;

  (*------------------------------------------------------------------------*)
  (* Bibiliotek einiger Matrix- und Vektoroperationen.                      *)
  (* Provides some matrix and vector operations.                            *)
  (*------------------------------------------------------------------------*)
  (* Implementation : Michael Riedl                                         *)
  (* Licence        : GNU Lesser General Public License (LGPL)              *)
  (*------------------------------------------------------------------------*)

  (* $Id: MatLib.def,v 1.12 2017/10/29 09:54:35 mriedl Exp mriedl $ *)

FROM Deklera IMPORT VEKTOR,MATRIX,SUPERVEKTOR,CMATRIX;

PROCEDURE InitIndex(VAR Index : ARRAY OF CARDINAL;
                        high  : CARDINAL);

          (*----------------------------------------------------------------*)
          (* Vorbelegt II[i] mit i*(i-1) DIV 2 , i=1..HIGH(Index) + 1       *)
          (*----------------------------------------------------------------*)

PROCEDURE Itab(i : CARDINAL) : CARDINAL;

          (*----------------------------------------------------------------*)
          (* Berechnet i*(i-1) DIV 2                                        *)
          (*----------------------------------------------------------------*)

PROCEDURE IJtab(i,j : CARDINAL) : CARDINAL;

          (*----------------------------------------------------------------*)
          (* Berechnet ITab(i) + j f\"ur i > j, sonst ITab(j) + i           *)
          (*----------------------------------------------------------------*)

PROCEDURE InvIJtab(    ij  : CARDINAL;
                   VAR i,j : CARDINAL);

          (*----------------------------------------------------------------*)
          (* Berechnet aus dem Paarindex ij die Indizes i,j mit i >= j      *)
          (*----------------------------------------------------------------*)

PROCEDURE Itab0(i : CARDINAL) : CARDINAL;

          (*----------------------------------------------------------------*)
          (* Analogon fuer Itab fuer Felder mit Indexstart 0                *)
          (*----------------------------------------------------------------*)

PROCEDURE IJtab0(i,j : CARDINAL) : CARDINAL;

          (*----------------------------------------------------------------*)
          (* Analogon fuer IJtab fuer Felder mit Indexstart 0               *)
          (*----------------------------------------------------------------*)

PROCEDURE InitLRArray(VAR X : ARRAY OF LONGREAL;
                          y : LONGREAL;
                          n : CARDINAL);

          (*----------------------------------------------------------------*)
          (* Initialisierrt das Feld X mit dem Wert y im Bereich X[0..n-1]  *)
          (* Benutzt ausgerollte Schleifen fuer n > 32                      *)
          (*                                                                *)
          (* Initialize Array X by zero, N beeing N the number of elements. *)
          (* Used unrolled loops for n > 32                                 *)
          (*----------------------------------------------------------------*)

PROCEDURE InitLCArray(VAR X : ARRAY OF LONGCOMPLEX;
                          y : LONGCOMPLEX;
                          n : CARDINAL);

          (*----------------------------------------------------------------*)
          (* Komplexe Version von InitLRArray - siehe dort                  *)
          (*----------------------------------------------------------------*)

PROCEDURE InitCMat(VAR A    : ARRAY OF ARRAY OF LONGCOMPLEX;
                       m,n  : CARDINAL;
                       Form : CHAR;
                       x    : LONGCOMPLEX);

          (*----------------------------------------------------------------*)
          (* Komplexe Version von InitMat - siehe dort                      *)
          (*----------------------------------------------------------------*)

PROCEDURE InitMat(VAR A    : ARRAY OF ARRAY OF LONGREAL;
                      m,n  : CARDINAL;
                      Form : CHAR;
                      x    : LONGREAL);

          (*----------------------------------------------------------------*)
          (* Vorbelegt die Matrix A[0..m-1,0..n-1] mit dem Wert x.          *)
          (* Fehlermeldung, wenn                                            *)
          (*   m-1 > HIGH(M) und/oder n-1 > HIGH(A[0])                      *) 
          (*   Form = "D" oder "d" und m # n                                *) 
          (* in Errors.Fehler und Errors.Fehlerflag.                        *)
          (* Form =                                                         *)
          (*   "a","A" : Alle Elemente werden mit x vorbelegt.              *)
          (*   "d","D" : Nur die Diagonale wird mit x vorbelegt, der Rest   *)
          (*             mit Null                                           *)
          (* Ist Form nicht CAP("a") oder CAP("d"), dann werden Vorgabe-    *)
          (* gem"a3 alle Elemente mit Null besetzt.                         *)
          (*----------------------------------------------------------------*)


PROCEDURE InitSV(VAR Sv   : SUPERVEKTOR;
                     dim  : CARDINAL;
                     Form : CHAR;
                     x    : LONGREAL);

          (*----------------------------------------------------------------*)
          (* Vorbelegt den Supervektor Sv mit dem Wert x.                   *)
          (* Erkl"ahrung : Siehe Routine VorbelegMat.                       *)
          (*----------------------------------------------------------------*)

PROCEDURE CopyMat(VAR B   : ARRAY OF ARRAY OF LONGREAL;  (*  ==> *)
                  VAR A   : ARRAY OF ARRAY OF LONGREAL;  (* <==  *)
                      m,n : CARDINAL);

          (*----------------------------------------------------------------*)
          (* Kopiert die Matrix A in die Matrix B (B:=A)                    *)
          (*----------------------------------------------------------------*)

PROCEDURE VektorinSV(VAR SV   : SUPERVEKTOR;
                         Vek  : VEKTOR;
                         n    : CARDINAL);

          (*----------------------------------------------------------------*)
          (* Addiert den Vektor Vek auf die Diagonale des Supervektors SV   *)
          (*----------------------------------------------------------------*)

PROCEDURE SvToMat(VAR SV     : ARRAY OF LONGREAL; (* SUPERVEKTOR *) (*  ==> *)
                  VAR Mat    : ARRAY OF ARRAY OF LONGREAL;          (* <==  *)
                      dim    : CARDINAL;
                      steuer : INTEGER);

          (*----------------------------------------------------------------*)
          (* Prozedur zur Transformation eines SUPERVEKTORs in eine normale *)
          (* quadratische MATRIX.                                           *)
          (*                                                                *)
          (* steuer =                                                       *)
          (*   -2 : Der Supervektor SV wird als untere Dreiecksmatrix in    *)
          (*        MAT eingesetzt, die anderen Matrixelemente werden mit   *)
          (*        Null vorbelegt                                          *)
          (*    1 : Der Supervektor SV wird als untere Dreiecksmatrix in    *)
          (*        MAT eingesetzt, die anderen Matrixelemente werden nicht *)
          (*        ver"andert                                              *)
          (*    0 : Der Supervektor SV wird so in MAT eingesetzt, da3 eine  *)
          (*        symmetrische, vollbesetzte Matrix entsteht.             *)
          (*    1 : Der Supervektor SV wird als obere Dreiecksmatrix in MAT *)
          (*        eingesetzt, die anderen Matrixelemente werden nicht     *)
          (*        veraendert.                                             *)
          (*    2 : Der Supervektor SV wird als obere Dreiecksmatrix in MAT *)
          (*        MAT eingesetzt, die anderen Matrixelemente werden mit   *)
          (*        Null vorbelegt                                          *)
          (*----------------------------------------------------------------*)

PROCEDURE MatToSv(VAR Mat    : MATRIX;       (*  ==> *)
                  VAR Sv     : SUPERVEKTOR;  (* <==  *)
                      dim    : CARDINAL;
                      steuer : INTEGER);

          (*----------------------------------------------------------------*)
          (* Prozedur zur Transformation einer auf einem zweidimensionalen  *)
          (* Feld gespeicherten (symmetrischen) Matrix in einen Supervektor *)
          (*                                                                *)
          (* steuer =                                                       *)
          (*   -1 : Die obere  Dreiecksmatrix wird in SV eingesetzt.        *)
          (*    0 : Die untere Dreiecksmatrix wird in SV eingesetzt.        *)
          (*    1 : Die untere Dreiecksmatrix wird in SV eingesetzt, aber   *)
          (*        es wird die (globale) Variable Errors.Fehler auf TRUE   *)
          (*        gesetzt und eine entsprechnde Fehlermeldung in          *)
          (*        Errors.Fehlerflag generiert, wenn die Matrix nicht sym- *)
          (*        metrisch ist. Im Fehlerfall ist das Ergebniss (Sv)      *)
          (*        unbestimmt.                                             *)
          (*----------------------------------------------------------------*)

PROCEDURE TriMat(VAR Mat    : MATRIX;    (* resultierende Matrix   *)
                 VAR HD,ND  : VEKTOR;    (* Haupt - Nebendiagonale *)
                     dim    : CARDINAL;
                     Form   : CARDINAL);

          (*----------------------------------------------------------------*)
          (* Wenn Form = 0 Transformation der in HD,ND gespeicherte         *)
          (*               tridiagonale Matrix in normale Matrixform.       *)
          (* Wenn Form = 1 wird aus der Matrix HD und ND abgespalten.       *)
          (*----------------------------------------------------------------*)

PROCEDURE TriSv(VAR SV    : SUPERVEKTOR; (* resultierende Matrix   *)
                VAR HD,ND : VEKTOR;      (* Haupt - Nebendiagonale *)
                    dim   : CARDINAL;
                    Form   : CARDINAL);

          (*----------------------------------------------------------------*)
          (* Wenn Form = 0 : Transformation der in HD,ND gespeicherte       *)
          (*                 tridiagonalen Matrix in Superlinearform.       *)
          (* Wenn Form = 1 : Abspalten der Haupt- und ersten Nebendiagonale *)
          (*                 von SV in HD,ND.                               *)
          (*----------------------------------------------------------------*)

PROCEDURE Stuerz(VAR Mat  : ARRAY OF ARRAY OF LONGREAL;
                     dim  : CARDINAL);

          (*----------------------------------------------------------------*)
          (* Spiegelt die Matrix Mat an der Hauptdiagonalen, d.h. die       *)
          (* Matrix wird transponiert.                                      *)
          (*----------------------------------------------------------------*)

PROCEDURE StuerzMN(VAR A   : ARRAY OF ARRAY OF LONGREAL;
                       M,N : CARDINAL);

          (*----------------------------------------------------------------*)
          (* Spiegelt die unsymmetrische Matrix A an der Hauptdiagonalen.   *)
          (* Die transponiert Matrix hat die Dimension N,M, A muss diesen   *)
          (* Speicherplatz bereithalten.                                    *)
          (*----------------------------------------------------------------*)

PROCEDURE Transpose(VAR A : ARRAY OF ARRAY OF LONGREAL; 
                        n : CARDINAL);

          (*----------------------------------------------------------------*)
          (* Spiegelt die quadratische Matrix Mat an der Hauptdiagonalen.   *)
          (* Diese Routine sollte fuer grosse n der Routien Stuerz vorge-   *)
          (* zogen werden                                                   *)
          (*----------------------------------------------------------------*)

PROCEDURE TransposeMN(    M,N : CARDINAL;
                      VAR A,T : ARRAY OF ARRAY OF LONGREAL);

          (*----------------------------------------------------------------*)
          (* Transpose of non-square matrix A of size M x N with additional *)
          (* storage for transposed matrix T                                *)
          (*----------------------------------------------------------------*)

PROCEDURE CStuerz(VAR Mat    : ARRAY OF ARRAY OF LONGCOMPLEX;
                      dim    : CARDINAL;
                      Conjug : BOOLEAN);

          (*----------------------------------------------------------------*)
          (* Spiegelt die komplexe Matrix Mat an der Hauptdiagonalen,       *)
          (* d.h. die Matrix wird transponiert.                             *)
          (* Wenn Conjug = TRUE wird beim Spiegeln das konjugiert           *)
          (* komplexe gesetzt (d.h. es wird die adjungierte Matrix erzeugt) *)
          (* Diese wird auch als hermitesch transponierte Matrix oder       *)
          (* transponiert-konjugierte Matrix bezeichnet.                    *)
          (*----------------------------------------------------------------*)

PROCEDURE StuerzRev(VAR A   : MATRIX;
                        dim : CARDINAL);

          (*----------------------------------------------------------------*)
          (* 'St\"urzt' die Matrix A derart, da\3 (f\"ur dim = 3)           *)
          (*   A[1,1] A[1,2] A[1,3]       A'[3,3] A'[2,3] A'[1,3]           *)
          (*   A[2,1] A[2,2] A[2,3]  ==>  A'[3,2] A'[1,2] A'[1,2]           *)
          (*   A[3,1] A[3,2] A[3,3]       A'[3,1] A'[2,1] A'[1,1]           *)
          (* d.h. es wird an der Diagonalen A[1,3], A[2,2] A[3,1]           *)
          (* gest"urzt.                                                     *)
          (*----------------------------------------------------------------*)

PROCEDURE IntSort(VAR IVek : ARRAY OF INTEGER;  (* Zu Sortierender Vektor   *)
                  VAR I    : ARRAY OF CARDINAL; (* Urspr"ungl. Positionen   *)
                      dim  : CARDINAL;          (* Dimension von IVek       *)
                      ord  : INTEGER);          (* Steuerparameter          *)

           (*---------------------------------------------------------------*)
           (* Sortiert den Vektor X nach der Gr"o3e der Elemente, wobei in  *)
           (* I die urspr"unglichen Positionen abgespeichert werden.        *)
           (*                                                               *)
           (* ord : -2 : Ansteigend nach Absolutwerten.                     *)
           (* ord : -1 : Ansteigend nach Werten.                            *)
           (* ord :  1 : Absteigend nach Werten.                            *)
           (* ord :  2 : Absteigend nach Absolutwerten.                     *)
           (*---------------------------------------------------------------*)

PROCEDURE QuickSort(VAR VEK  : ARRAY OF LONGREAL;
                    VAR II   : ARRAY OF CARDINAL;
                        dim  : CARDINAL;
                        ord  : INTEGER);

          (*----------------------------------------------------------------*)
          (* Sortiert Real-Feld und merkt sich auf II die  vorherigen       *)
          (* Positionen.                                                    *)
          (* ord      = -1 : Kleinstes Element zuerst.                      *)
          (* ord      =  1 : Gr\"o\3tes Element zuerst.                     *)
          (* ABS(ord) #  1 : Keine Ordnung                                  *)
          (*----------------------------------------------------------------*)

PROCEDURE Sort(VAR VEK   : VEKTOR;
                   dim   : CARDINAL;
                   ord   : INTEGER);

          (*----------------------------------------------------------------*)
          (* Sortiert VEKTOR VEK nach Gr"o3e der Elemente.                  *)
          (* Wenn ord =  1 in absteigender Reihenfolge                      *)
          (* Wenn ord #  1 in ansteigender Reihenfolge                      *)
          (*----------------------------------------------------------------*)

PROCEDURE SortVek(VAR X   : ARRAY OF LONGREAL; (* Zu Sortierender Vektor *)
                  VAR I   : ARRAY OF CARDINAL; (* Urspr"ungl. Positionen *)
                      dim : CARDINAL;          (* Dimension von X        *)
                      ord : INTEGER);          (* Steuerparameter        *)

          (*----------------------------------------------------------------*)
          (* Sortiert den Vektor X nach der Gr"o3e der Elemente, wobei in   *)
          (* I die urspr"unglichen Positionen abgespeichert werden.         *)
          (*                                                                *)
          (* ord : -2 : Ansteigend nach Absolutwerten.                      *)
          (* ord : -1 : Ansteigend nach Werten.                             *)
          (* ord :  1 : Absteigend nach Werten.                             *)
          (* ord :  2 : Absteigend nach Absolutwerten.                      *)
          (*                                                                *)
          (* Anzahl Vergleiche O(n*n), Anzahl Vertauschungen O(n)           *)
          (* daher nicht fuer sehr grosse Werte von dim benutzen            *)
          (*----------------------------------------------------------------*)

PROCEDURE CSortVek(VAR X   : ARRAY OF LONGCOMPLEX; (* Zu Sortierend         *)
                   VAR I   : ARRAY OF CARDINAL;  (* Urspr\"ungl. Positionen *)
                       dim : CARDINAL;           (* Dimension von X         *)
                       ord : INTEGER);           (* Steuerparameter         *)

          (*----------------------------------------------------------------*)
          (* Complexwertige Verions von SortVek                             *)
          (*----------------------------------------------------------------*)

PROCEDURE MatVekProd(VAR Y     : ARRAY OF LONGREAL;
                     VAR A     : ARRAY OF ARRAY OF LONGREAL;
                         X     : ARRAY OF LONGREAL;
                         M,N   : CARDINAL;
                         trans : BOOLEAN);

          (*----------------------------------------------------------------*)
          (* Berechent fuer ein quadratische Matrix A den Vektor Y = A*X    *)
          (* oder, fuer trans = TRUE, Y = A'*X wobei A' die transponierte   *)
          (* Matrix A ist.                                                  *)
          (*                                                                *)
          (*  A : M x N matrix, A[0..M-1][0..N-1]                           *)
          (*  X : N x 1 vector, X[0..N-1] fuer trans = falsch               *)
          (*      M x 1 vector, X[0..M-1] fuer trans = wahr                 *)
          (*  Y : M x 1 vector, X[0..N-1] fuer trans = falsch               *)
          (*      N x 1 vector, X[0..M-1] fuer trans = wahr                 *)
          (*----------------------------------------------------------------*)

PROCEDURE CMatVekProd(VAR X   : ARRAY OF LONGCOMPLEX;
                      VAR Mat : ARRAY OF ARRAY OF LONGCOMPLEX;
                          Vek : ARRAY OF LONGCOMPLEX;
                          dim : CARDINAL); (* Dim d. Matrix *)

          (*----------------------------------------------------------------*)
          (* Mulipliziet den komplexen Vektor Vek mit der komplexen Matrix  *)
          (* Mat. (X = Mat*Vek)                                             *)
          (*----------------------------------------------------------------*)

PROCEDURE SvVekProd(VAR Y   : ARRAY OF LONGREAL; (*  ==> *)
                    VAR A   : ARRAY OF LONGREAL; (* <==  SUPERVEKTOR *)
                    VAR X   : ARRAY OF LONGREAL; (* <==  *)
                        dim : CARDINAL);

          (*----------------------------------------------------------------*)
          (* Berechent fuer ein linear gespeicherte quadratische Matrix A   *)
          (* den Vektor Y = A*X                                             *)
          (*----------------------------------------------------------------*)

PROCEDURE MatMatProd(VAR C   : MATRIX;        (* Resultierende MATRIX *)
                     VAR A   : MATRIX;
                     VAR B   : MATRIX;
                         dim : CARDINAL);

          (*----------------------------------------------------------------*)
          (* Berechnet die Matrixmultiplikation C = A*B                     *)
          (* Die Matrix C kann mit A oder B identisch sein, wenn dies       *)
          (* keine Anforderung ist kann auch die schnellere Routine         *)
          (* MatMatProdNN genutzt werden. Eine andrere Alternative ist      *)
          (* LibDBlas.dgemm                                                 *)
          (*----------------------------------------------------------------*)

PROCEDURE MatMatProdNN(    M,N,K : CARDINAL;
                       VAR C     : ARRAY OF ARRAY OF LONGREAL;
                       VAR A,B   : ARRAY OF ARRAY OF LONGREAL);

          (*----------------------------------------------------------------*)
          (* MatMatProdNN calculates the product of matrix A times matrix   *)
          (* B, the resulting matrix is stored in C. The procedure uses     *)
          (* unrolled loops.                                                *)
          (*                                                                *)
          (*  ==> M : row dimension of matix A and resulting matrix C       *)
          (*  ==> K : row dimension of matix B and length of a single       *)
          (*          vector in matrixs A                                   *)
          (*  ==> N : length of a single vector in matrixs B                *)
          (*  ==> A : is a M by K matrix A[0..M-1][0..K-1]                  *)
          (*  ==> B : is a K by N matrix B[0..K-1][0..N-1]                  *)
          (*                                                                *)
          (* <==  C : is a M by N matrix C[0..M-1][0..N-1] (result)         *)
          (*                                                                *)
          (*          (A11 A12 ... A1K)   (B11 B12 B13 B14 ... B1N)         *)
          (*  (C) =   (A21 A22 ... A2K) x (B21 B22 B23 B24 ... B2N)         *)
          (*          (A31 A32 ... A3K)   (BK1 BK2 BK3 BK4 ... BKN)         *)
          (*          (AM1 AM2 ... AMK)                                     *)
          (*----------------------------------------------------------------*)

PROCEDURE MatMatProdTN(    M,N,K : CARDINAL;
                       VAR C     : ARRAY OF ARRAY OF LONGREAL;
                       VAR A,B   : ARRAY OF ARRAY OF LONGREAL);

          (*----------------------------------------------------------------*)
          (* This procedure computes C = A(transpose)*B using unrolled      *)
          (* loops                                                          *)
          (*                                                                *)
          (* Input Paramaters                                               *)
          (*                                                                *)
          (*   m  : number of rows in C, no. of columes in A                *)
          (*   n  : no. of columns in C and B                               *)
          (*   k  : no. of rows in A and B                                  *)
          (*   A  : first matrix,  K x M                                    *)
          (*   B  : second matrix, K x N                                    *)
          (*                                                                *)
          (* Output Parameter                                               *)
          (*                                                                *)
          (*   C  : A(transpose)*B (M x N)                                  *)
          (*----------------------------------------------------------------*)

PROCEDURE MatMatProdNT(    M,N,K : CARDINAL;
                       VAR C     : ARRAY OF ARRAY OF LONGREAL;
                       VAR A,B   : ARRAY OF ARRAY OF LONGREAL);

          (*----------------------------------------------------------------*)
          (* This procedure computes C = A * B(transpose) using unrolled    *)
          (* loops                                                          *)
          (*                                                                *)
          (* Input Paramaters                                               *)
          (*                                                                *)
          (*   m  : number of rows in C and A                               *)
          (*   n  : no. of columns in C and rows of B                       *)
          (*   k  : no. of columns in A and B                               *)
          (*   A  : first matrix,  M x K                                    *)
          (*   B  : second matrix, N x K                                    *)
          (*                                                                *)
          (* Output Parameter                                               *)
          (*                                                                *)
          (*   C  : A*B(transpose) (M x N)                                  *)
          (*----------------------------------------------------------------*)

PROCEDURE MatMatProdTT(    M,N,K : CARDINAL;
                       VAR C     : ARRAY OF ARRAY OF LONGREAL;
                       VAR A,B   : ARRAY OF ARRAY OF LONGREAL);

          (*----------------------------------------------------------------*)
          (* This procedure computes c=a(transpose)*b(transpose) using      *)
          (* unrolled loops. It transposes B in case N = K to gain speed    *)
          (* as one cannot find an chache optimal loop orderung for         *)
          (* A^t x B^t                                                      *)
          (*                                                                *)
          (* Input Paramaters                                               *)
          (*                                                                *)
          (*   m  : number of rows in C and columes in A                    *)
          (*   n  : no. of columns in C and B                               *)
          (*   k  : no. of rows in A and columes in B                       *)
          (*   A  : first matrix,  K x M                                    *)
          (*   B  : second matrix, N x K (unchanges even if transposed !)   *)
          (*                                                                *)
          (* Output Parameter                                               *)
          (*                                                                *)
          (*   C  : A(transpose)*B(transpose) (M x N)                       *)
          (*----------------------------------------------------------------*)

PROCEDURE CMatMatProd(VAR C   : CMATRIX;  (* Resultierende MATRIX *)
                      VAR A   : CMATRIX;
                      VAR B   : CMATRIX;
                          dim : CARDINAL);

          (*----------------------------------------------------------------*)
          (* Berechnet die komplexe Matrixmultiplikation C = A*B            *)
          (* Die Matrix C kann mit A oder B identisch sein, wenn dies       *)
          (* keine Anforderung ist kann auch universellere Routine          *)
          (* LibDBlas.zgemm genutzt werden.                                 *)
          (*----------------------------------------------------------------*)

PROCEDURE MatSvProdNN(    dim : CARDINAL;
                      VAR C   : ARRAY OF ARRAY OF LONGREAL; 
                      VAR A   : ARRAY OF ARRAY OF LONGREAL;
                      VAR Bp  : ARRAY OF LONGREAL); (* SUPERVEKTOR *)

          (*----------------------------------------------------------------*)
          (* Matrixmultiplikation C = A x B mit quadratischer, eventuell    *)
          (* unsymmetrischer Matrix A und symmetrischer Matrix B in         *)
          (* gepackter Form jeweils der Dimension "dim".                    *)
          (* Die Matrizen C und A koennen auf denselben Speicherplatz       *)
          (* verweisen.                                                     *)
          (*----------------------------------------------------------------*)

PROCEDURE MatSvProdTN(    dim : CARDINAL;
                      VAR C   : ARRAY OF ARRAY OF LONGREAL;  
                      VAR A   : ARRAY OF ARRAY OF LONGREAL;
                      VAR Bp   : ARRAY OF LONGREAL); (* SUPERVEKTOR *)

          (*----------------------------------------------------------------*)
          (* Matrixmultiplikation C = A^t x B mit quadratischer,            *)
          (* unsymmetrischer transponierter Matrix A und symmetrischer      *)
          (* Matrix B in gepackter Form jeweils der Dimension "dim".        *)
          (* Die Matrizen C und A koennen NICHT auf denselben Speicherplatz *)
          (* verweisen.                                                     *)
          (*----------------------------------------------------------------*)

PROCEDURE SvMatProdNN(    dim : CARDINAL;
                      VAR C   : ARRAY OF ARRAY OF LONGREAL;
                      VAR Ap  : ARRAY OF LONGREAL; (* SUPERVEKTOR *)
                      VAR B   : ARRAY OF ARRAY OF LONGREAL);

          (*----------------------------------------------------------------*)
          (* Matrixmultiplikation C = A x B mit symmetrische, quadratischer *)
          (* Matrix A in gepackter Form und quadratischer unsymmetrischer   *)
          (* Matrix B, jeweils der Dimension "dim".                         *)
          (* Die Matrizen C und B koennen NICHT auf denselben Speicherplatz *)
          (* verweisen.                                                     *)
          (*----------------------------------------------------------------*)

PROCEDURE SvMatProdNT(    dim : CARDINAL;
                      VAR C   : ARRAY OF ARRAY OF LONGREAL;
                      VAR Ap  : ARRAY OF LONGREAL; (* SUPERVEKTOR *)
                      VAR B   : ARRAY OF ARRAY OF LONGREAL);

          (*----------------------------------------------------------------*)
          (* Matrixmultiplikation C = A x B^t mit symmetrische, quadrat-    *)
          (* ischer Matrix A in gepackter Form und quadratischer unsym-     *)
          (* metrischer gestuerzter Matrix B, jeweils der Dimension "dim".  *)
          (* Die Matrizen C und B koennen NICHT auf denselben Speicherplatz *)
          (* verweisen.                                                     *)
          (*----------------------------------------------------------------*)

PROCEDURE SvSvProdNN(    dim    : CARDINAL;
                     VAR C      : ARRAY OF ARRAY OF LONGREAL;
                     VAR Ap,Bp  : ARRAY OF LONGREAL); (* SUPERVEKTOR *)

          (*----------------------------------------------------------------*)
          (* Matrixmultiplikation C = A x B mit symmetrischen quadratischen *)
          (* Matrizen A und B in gepackter Form jeweils der Dimension "dim" *)
          (*----------------------------------------------------------------*)

PROCEDURE SkalarProd(VAR A,B : ARRAY OF LONGREAL;
                         dim : CARDINAL) : LONGREAL;

          (*----------------------------------------------------------------*)
          (* Berechent das Skalarprodukt zweier Vektoren                    *)
          (*----------------------------------------------------------------*)

PROCEDURE CSkalarProd(VAR A,B : ARRAY OF LONGCOMPLEX;
                          dim : CARDINAL) : LONGCOMPLEX;

          (*----------------------------------------------------------------*)
          (* Berechent das Skalarprodukt zweier komplexer Vektoren          *)
          (*----------------------------------------------------------------*)

PROCEDURE MatSkalProd(VAR MAT : MATRIX;
                          x   : LONGREAL;
                          dim : CARDINAL);

          (*----------------------------------------------------------------*)
          (* Multipliziet Matrix MAT mit x                                  *)
          (*----------------------------------------------------------------*)

PROCEDURE KreuzProd(VAR Mat   : MATRIX;    (* Resultierende Matrix *)
                    VAR Vec1  : VEKTOR;
                    VAR Vec2  : VEKTOR;
                        dim   : CARDINAL);  (* L\"ange der Vektoren *)

          (*----------------------------------------------------------------*)
          (* Ber. das Kreuzprodukt zweier Vektoren                          *)
          (* Dabei Entsteht eine Matrix !                                   *)
          (*----------------------------------------------------------------*)

PROCEDURE SpaltVek(VAR SpVek : VEKTOR;
                   VAR MAT   : MATRIX;
                       dim   : CARDINAL;
                       j     : CARDINAL); (* j.-Spalte wird SpVek *)

          (*----------------------------------------------------------------*)
          (* Der Spaltenvektor j der Matrix MAT wird auf                    *)
          (* SpVek gelegt.                                                  *)
          (*----------------------------------------------------------------*)

PROCEDURE VekInMat(VAR VEK : VEKTOR;
                   VAR MAT : MATRIX;
                       j   : CARDINAL;
                       dim : CARDINAL);

          (*----------------------------------------------------------------*)
          (* Setzt den Vektor VEK als j-te Spalte in die                    *)
          (* Matrix MAT ein.                                                *)
          (*----------------------------------------------------------------*)

PROCEDURE InsertMat(VAR A     : MATRIX;
                        adimi : CARDINAL; (* first dimension of A *)
                        adimj : CARDINAL; (* second dimension of A *)
                    VAR B     : MATRIX;   (* Matrix to be inserted *)
                        in,jn : CARDINAL; (* 1. & 2. dimension ob B *)
                        is,js : CARDINAL; (* Startpoints where to insert *)
                    VAR ifehl : CARDINAL);

          (*----------------------------------------------------------------*)
          (* Einsetzen einer Matrix B in die Matix A an den Startpunkten    *)
          (* A[is,js].                                                      *)
          (*                                                                *)
          (* Fehlermeldungen:                                               *)
          (*                                                                *)
          (*   ifehl =   0 : alles in Ordnung                               *)
          (*   ifehl =   1 : Dimensionen von B groesser als die von A       *)
          (*   ifehl =  10 : Startpunkt is falsch gewaehlt                  *)
          (*   ifehl = 100 : Startpunkt js falsch gewaehlt                  *)
          (*----------------------------------------------------------------*)

PROCEDURE TriInfNorm(    n   : CARDINAL;
                     VAR D,E : ARRAY OF LONGREAL) : LONGREAL;

          (*----------------------------------------------------------------*)
          (* Berechne die Inf-Norm (maximale Zeilensumme) der symmetrischen *)
          (* Tridiagonalmatrix, in D,E uebergeben wurde. Diese ist fuer     *)
          (* diesen speziellen Fall identisch mit der 1-Norm (maximale      *)
          (* Spaltensumme) der Tridiagonalmatrix.                           *)
          (*                                                                *)
          (* Calculates the Infiniti-Norm (maximal row sum) of a symmetric  *)
          (* trigiagonal matrix provided in vectors D,E. This norm is       *)
          (* identical to the 1-Norm (maximal colume sum) in the special    *)
          (* case of a symmetic matrix                                      *)
          (*----------------------------------------------------------------*)

PROCEDURE MatNormL1(VAR A    : ARRAY OF ARRAY OF LONGREAL;
                        N    : INTEGER;
                    VAR Norm : LONGREAL);

         (*-----------------------------------------------------------------*)
         (* Berechnet die L1 Norm der Matrix A als maximale Zeilensumme der *)
         (* Absolutwerte.                                                   *)
         (*-----------------------------------------------------------------*)

PROCEDURE ENorm(VAR X   : ARRAY OF LONGREAL;
                    s,e : CARDINAL) : LONGREAL;

          (*----------------------------------------------------------------*)
          (* Berechnet \sqrt{ \sum_{i=s}^e X_i^2} , wobei der erste         *)
          (* (m"ogliche) Index s=1 ist.                                     *)
          (*                                                                *)
          (* The euclidian morm is computed by accumulating the sum of      *)
          (* squares in three different sums. The sums of squares for the   *)
          (* small and large components are scaled so that no overflows     *)
          (* occur. Non-destructive underflows are permitted. Underflows    *)
          (* and overflows do not occur in the computation of the unscaled  *)
          (* sum of squares for the intermediate components.                *)
          (* The definitions of small, intermediate and large components    *)
          (* depend on two constants, MinLR and MaxLR. The main             *)
          (* restrictions on these constants are that MinLR**2 not          *)
          (* underflow ans MaxLR**2 not overflow. The constants             *)
          (* given here are suitable for every known computer.              *)
          (*                                                                *)
          (* Argonne national laboratory. MINPACK project. March 1980.      *)
          (* Burton S. Garbow, Kenneth E. Hillstrom, Jorge J. More          *)
          (*                                                                *)
          (* Anpassung an Modula-2 und "Uberarbeitung : M.Riedl, 23.11.93.  *)
          (* Berechnet \sqrt{ \sum_{i=s}^e X_i^2}                           *)
          (*----------------------------------------------------------------*)

PROCEDURE NormVek(VAR X       : ARRAY OF LONGREAL;
                      dim     : CARDINAL;
                      MaxNorm : BOOLEAN);

          (*----------------------------------------------------------------*)
          (* Normiert den Vektor X im Bereich X[1] bis X[dim].              *)
          (* Wenn MaxNorm = TRUE so, da3 das Maximale Element die Gr"o3e    *)
          (* 1 hat, ansonsten nach der Euklidschen Norm auf 1               *)
          (*----------------------------------------------------------------*)

PROCEDURE NormMat(VAR A       : ARRAY OF ARRAY OF LONGREAL;
                      dim     : CARDINAL;
                      VekAnz  : CARDINAL;
                      MaxNorm : BOOLEAN);

          (*----------------------------------------------------------------*)
          (* MaxNorm =                                                      *)
          (*   FALSE : Die Spalten von MAT werden auf 1.0 normiert          *)
          (*           (euklidische Norm).                                  *)
          (*   TRUE  : Die Spalten von MAT werden so normiert,              *)
          (*           da\3 das betragsgr\"o\3te Element 1.0 wird           *)
          (*           (Maximumnorm) .                                      *)
          (*----------------------------------------------------------------*)

PROCEDURE CNormMat(VAR Mat    : ARRAY OF ARRAY OF LONGCOMPLEX;
                       dim    : CARDINAL;
                       VekAnz : CARDINAL;
                       MaxNorm: BOOLEAN);
          (*----------------------------------------------------------------*)
          (* Komplexes Analogon von NormMat, siehe dort.                    *)
          (*----------------------------------------------------------------*)

PROCEDURE MatSpur(VAR Spur : LONGREAL;
                  VAR MAT  : MATRIX;
                      dim  : CARDINAL);

          (*----------------------------------------------------------------*)
          (* Berechnet die Spur einer quadrarischen Matrix MAT.             *)
          (*----------------------------------------------------------------*)

PROCEDURE SvSpur(VAR Spur : LONGREAL;
                 VAR SV   : SUPERVEKTOR;
                     dim  : CARDINAL);

          (*----------------------------------------------------------------*)
          (* Berechnet die Spur einer symmetrischen                         *)
          (* Matrix SV in Superlinearform.                                  *)
          (*----------------------------------------------------------------*)

PROCEDURE MinMaxVek(VAR Min,Max : LONGREAL;
                    VAR Vek     : ARRAY OF LONGREAL;
                        dim     : CARDINAL;
                        Form    : CARDINAL);

          (*----------------------------------------------------------------*)
          (* Liefert das Minimum und das Maximum eines Vektor.              *)
          (* Wenn Form = 1 werden nur die Absolutwerte betrachtet,          *)
          (* wobei Min,Max VORZEICHENBEHAFTET zur"uckgegeben werden.        *)
          (*----------------------------------------------------------------*)

PROCEDURE MinMaxMat(VAR  Min,Max : LONGREAL;
                    VAR  Mat     : MATRIX;
                         dim     : CARDINAL;  (* Dimension der Matrix *)
                         Form    : CARDINAL);

          (*----------------------------------------------------------------*)
          (* Ermittelt das kleinste und gr"o3te Element einer Matrix.       *)
          (* Wenn nur die Absolutwerte gew"unscht werden, werden diese      *)
          (* jedoch VORZEICHENBEHAFTET zur"uckgegeben !                     *)
          (* Wenn Form =                                                    *)
          (*   0 : alle Elemente.                                           *)
          (*   1 : nur die Absolutwerte aller Elemente.                     *)
          (*   2 : nur Nebendiagonalelemente.                               *)
          (*   3 : nur Absolutwerte der Nebendiagonale.                     *)
          (*   4 : nur die Diagonalelemente.                                *)
          (*   5 : nur Absolutwerte der Diagonale.                          *)
          (* werden ber"ucksichtigt.                                        *)
          (*----------------------------------------------------------------*)

PROCEDURE MinMaxSv(VAR  Min,Max : LONGREAL;
                   VAR  SV      : ARRAY OF LONGREAL; (* SUPERVEKTOR *)
                        dim     : CARDINAL;  (* Dimension der Matrix *)
                        Form    : CARDINAL);

          (*----------------------------------------------------------------*)
          (* Ermittelt das kleinste und gr"o3te Element einer Matrix in     *)
          (* Superlinearform. Wenn nur die Absolutwerte gew"unscht werden,  *)
          (* so werden diese jedoch VORZEICHENBEHAFTET zur"uckgegeben !     *)
          (* Wenn Form =                                                    *)
          (*   0 : alle Elemente                                            *)
          (*   1 : nur die Absolutwerte                                     *)
          (*   2 : nur Nichtdiagonalelemente                                *)
          (*   3 : nur Absolutwerte der Nichtdiagonalelemente               *)
          (*   4 : nur die Diagonalelemente                                 *)
          (*   5 : nur Absolutwerte der Diagonale                           *)
          (* werden ber"ucksichtigt.                                        *)
          (*----------------------------------------------------------------*)

PROCEDURE ZufallVec(VAR V   : ARRAY OF LONGREAL;
                        dim : CARDINAL);

          (*----------------------------------------------------------------*)
          (* Erzeugt einen Vektor mit Zufallszahlen im Bereich [-1,1]       *)
          (*                                                                *)
          (* Generate a vector with random numbers in the interval [-1,1]   *)
          (*----------------------------------------------------------------*)

PROCEDURE ZufallMat(VAR A   : ARRAY OF ARRAY OF LONGREAL;
                        m,n : CARDINAL;
                        sym : BOOLEAN);

          (*----------------------------------------------------------------*)
          (* Erzeugt eine m x n Matrix mit Zufallszahlen im Bereich [-1,1]. *)
          (* Ist der Parameter sym auf wahr gesetzt wird, unter der Voraus- *)
          (* setzung dass m = n ist, eine symmetrische Matrix erzeugt.      *)
          (* Die Variable Errors.Fehler ist im Fehlerfall wahr.             *)
          (*                                                                *)
          (* Generate a m x n matrix A o with random numbers in the         *)
          (* interval [-1,1]. If the parameter sym is true a symmetric      *)
          (* matrix is generated under the condition that m is equal to n.  *)
          (* The variable Errors.Fehler will be set to true in case of an   *)
          (* error.                                                         *)
          (*----------------------------------------------------------------*)

PROCEDURE ZufallSv(VAR SV  : ARRAY OF LONGREAL; (* SUPERVEKTOR *)
                       dim : CARDINAL);

          (*----------------------------------------------------------------*)
          (* Erzeugt einen Supervektor mit Zufallszahlen in [-1,1]          *)
          (*                                                                *)
          (* Generate a supervector with random numbers (interval [-1,1])   *)
          (*----------------------------------------------------------------*)

PROCEDURE SumVek(VAR X   : ARRAY OF LONGREAL;
                     s,e : CARDINAL) : LONGREAL;

          (*----------------------------------------------------------------*)
          (* Berechnet \sum_{i=s}^e X_i, wobei der erste (m"ogliche)        *)
          (* Index 1 ist.                                                   *)
          (*                                                                *)
          (* Summation des Vektors X nach Kahans Summationsformel.          *)
          (* Sehr pr"aziser, aber aufwendiger Summationsalgorithmus.        *)
          (* Vorsicht bei einigen Optimierer, es kann sein, da3 das berech- *)
          (* nete "gard-digit" (c) wegoptimiert wird !!!                    *)
          (*----------------------------------------------------------------*)

PROCEDURE AbsSumVek(VAR X   : ARRAY OF LONGREAL;
                        s,e : CARDINAL) : LONGREAL;

          (*----------------------------------------------------------------*)
          (* Berechnet \sum_{i=s}^e |X_i|, wobei der erste (m"ogliche)      *)
          (* Index 1 ist.                                                   *)
          (* Bis auf die Absolutwerte identisch mit SumVek.                 *)
          (*----------------------------------------------------------------*)

PROCEDURE NeumaierSum(VAR X : ARRAY OF LONGREAL;
                          N : CARDINAL) : LONGREAL;

          (*----------------------------------------------------------------*)
          (* Summation des Vektors X nach der Neumaier Summationsformel.    *)
          (* Sehr pr"aziser, aber aufwendiger Summationsalgorithmus.        *)
          (* Vorsicht bei einigen Optimierer, es kann sein, da3 das berech- *)
          (* nete "gard-digit" (c) wegoptimiert wird !!!                    *)
          (*----------------------------------------------------------------*)

PROCEDURE NeumaierProdSum(VAR X,Y : ARRAY OF LONGREAL;
                              N   : CARDINAL) : LONGREAL;

          (*---------------------------------------------------------------*)
          (* Berechne die Summe des paarweisen Produkts der Vektoren       *)
          (* X und Y nach dem Neumaierschen Summationsalgorithmus          *)
          (*---------------------------------------------------------------*)

PROCEDURE SumVekUR(VAR X   : ARRAY OF LONGREAL;
                       a,e : CARDINAL) : LONGREAL;

          (*----------------------------------------------------------------*)
          (* Schnellere aber weniger praezise Version von SumVek, benutzt   *)
          (* ausgerollte Schleifen aber nicht Kahans Summationsformel       *)
          (*----------------------------------------------------------------*)

PROCEDURE AbsSumVekUR(VAR X   : ARRAY OF LONGREAL;
                          a,e : CARDINAL) : LONGREAL;

          (*----------------------------------------------------------------*)
          (* Bis auf die Absolutwerte identisch mit SumVekUR                *)
          (*----------------------------------------------------------------*)

END MatLib.
