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.