DEFINITION MODULE LinLib;
(*------------------------------------------------------------------------*)
(* Modul zur L"osung linearer Gleichungssystem *)
(* Modul for solving lineare systems of equations. *)
(*------------------------------------------------------------------------*)
(* Implementation : Michael Riedl *)
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
(* $Id: LinLib.def,v 1.2 2017/10/04 07:02:51 mriedl Exp mriedl $ *)
FROM Deklera IMPORT VEKTOR,MATRIX,CARDVEKTOR,CVEKTOR,CMATRIX,PMATRIX;
PROCEDURE CondHager(VAR A : ARRAY OF ARRAY OF LONGREAL;
N : INTEGER;
VAR Cond : LONGREAL;
VAR iFehl : INTEGER);
(*-----------------------------------------------------------------*)
(* Estimates the L1 condition number of a matrix *)
(* *)
(* Hager, William; Condition estimates, SIAM journal on scientific *)
(* and statistical computing, 5,2, 311-316 (1984) *)
(* *)
(* This procedure is based on a Fortran version by John Burkardt *)
(* *)
(* iFehl = 1 : Matrix A cannot be factorizised *)
(* iFehl = 2 : not enought free memory *)
(* *)
(* In case iFehl # 0 Cond is MAX(LONGREAL) *)
(*-----------------------------------------------------------------*)
PROCEDURE Zerlege(VAR A : MATRIX; (* Koeffizientenmatrix *)
VAR IndexI : CARDVEKTOR;
VAR IndexJ : CARDVEKTOR;
VAR RPivot : VEKTOR;
VAR Det : LONGREAL; (* Determinante der Matrix *)
dim : CARDINAL; (* Dimension der qadr. Matrix *)
MaxPivot : BOOLEAN); (* w=Totalpivot, f=Zeilenpivot *)
(*----------------------------------------------------------------*)
(* Zerlegt die Matrix A nach dem Gaussschen Eliminationsverfahren *)
(* In den beiden Feldern IndexI und IndexJ werden dabei die *)
(* Informationen \"uber das Zeilen- bzw. Spaltenvertauschen f"ur *)
(* die Routine BerechneLoesung gespeichert. *)
(* Wenn A singul\"ar oder unl\"osbar wird die globale Variable *)
(* Fehler auf TRUE gesetzt. *)
(*----------------------------------------------------------------*)
PROCEDURE BerechneLoesung(VAR A : MATRIX;
VAR X : VEKTOR;
C : VEKTOR;
VAR RPivot : VEKTOR;
VAR IndexI : CARDVEKTOR;
VAR IndexJ : CARDVEKTOR;
dim : CARDINAL);
(*----------------------------------------------------------------*)
(* Berechnet den L\"osungsvektor X der Gleichung A X = C mit der *)
(* in der Zerlege nach Gauss zerlegten Matrix A *)
(* IndexI : Zeilenvertauschungsinformationen. *)
(* IndexJ : Spaltenvertauschungsinformationen. IndexJ_k = k wenn *)
(* nur Zeilenpivot. *)
(*----------------------------------------------------------------*)
PROCEDURE Det(VAR A : MATRIX; dim : CARDINAL) : LONGREAL;
(*----------------------------------------------------------------*)
(* Berechnet die Determinante von A, wobei die Koeffizienten- *)
(* matrix A nicht zerst\"ort wird. *)
(* Wenn Fehler auftreten, wird MAX(LONGREAL) zur\"uckgegeben und *)
(* die globale Variable Fehler auf TRUE gesetzt. *)
(*----------------------------------------------------------------*)
PROCEDURE Gauss(VAR A : MATRIX; (* Koeffizientenmatrix *)
VAR X : VEKTOR; (* L\"osungsvektor des LG *)
C : VEKTOR; (* Konstantenvektor des LG *)
dim : CARDINAL;
VAR Det : LONGREAL; (* Determinante von A *)
MaxIter : CARDINAL;
PivStrat : CARDINAL; (* Pivotstrategie *)
VAR iFehl : INTEGER); (* Fehlernummer *)
(*----------------------------------------------------------------*)
(* Berechnet die L\"osung des Gleichngssystems Ax = c durch *)
(* Gausselemination. Die Pivotstrategie ist dabei entweder eine *)
(* Totalpivotierung f\"ur PivStrat = 1, ansonsten wird (zeit - *)
(* sparend) mit Zeilenpivotierung gearbeitet. *)
(* Die Matrix A wird (au\3er bei iFehl = 1) nicht ver\"andert. *)
(* Wenn MaxIter > 1 wird der L\"osungsvektor mit dem Residuen- *)
(* vektor r = (Ax - c) nach A x2 = r und x = x-x2 MaxIter-mal *)
(* nachiteriert. *)
(* Wenn die Matrix singul\"ar oder sonstwie nicht l\"osbar ist, *)
(* wird die globale Variable Fehler auf TRUE gesetzt. *)
(* iFehl = 0 : kein Fehler aufgetreten. *)
(* 1 : Koeffizentenmatrix zerst\"ort (kein Nachiterieren) *)
(* 2 : Iteration nicht konvergent. *)
(* 3 : Koeffizentenmatix singul\"ar. *)
(* 4 : Dimension der Matrix < 1 *)
(*----------------------------------------------------------------*)
PROCEDURE Householder(VAR A : MATRIX;
VAR X : ARRAY OF LONGREAL;
C : ARRAY OF LONGREAL;
n : CARDINAL; (* Anzahl der Gleichungen, n >= m *)
m : CARDINAL); (* Anzahl der Unbekannten *)
(*----------------------------------------------------------------*)
(* Lit.: J.Stoer, 'Einf\"uhrung in die numerische Mathematik 1', *)
(* 4.Auflage, Berlin, Springer (1983). *)
(* L\"o\3t lineares Gleichungssystem A*X = C. *)
(* Wenn die Matrix singul\"ar oder sonstwie nicht l\"osbar ist, *)
(* wird die globale Variable Fehler auf TRUE gesetzt. Das Ver- *)
(* fahren kann insbesondere dann angewandt werden, wenn mehr *)
(* Gleichungen als Unbekannte vorhanden sind und nur diese zu- *)
(* s\"atzlichen Gleichungen eine eindeutige Bestimmung des *)
(* L\"osungsvektors zulassen. *)
(*----------------------------------------------------------------*)
PROCEDURE PHouseholder(VAR A : PMATRIX;
VAR X : ARRAY OF LONGREAL;
VAR C : ARRAY OF LONGREAL;
m : CARDINAL; (* Anzahl d. Gleichungen, m >= n *)
n : CARDINAL); (* Anzahl d. Unbekannten *)
(*----------------------------------------------------------------*)
(* Diese Routine ist mit der obenstehenden Householderroutine bis *)
(* auf die Form der Darstellung der Matrix A identisch. *)
(*----------------------------------------------------------------*)
PROCEDURE HhLin(VAR A : ARRAY OF LONGREAL;
VAR X : ARRAY OF LONGREAL;
VAR C : ARRAY OF LONGREAL;
n : CARDINAL; (* Zahl d. Gleichungen, n >= m *)
m : CARDINAL); (* Zahl d. Unbekannten *)
(*----------------------------------------------------------------*)
(* Die PROCEDURE HhLin ist, bis auf die Form der \"Ubergabe *)
(* der Parameter, identisch mit der PROCEDURE Householder. *)
(* *)
(* Die Koeffizientenmatrix A mu\3 als eindimensionales, ge - *)
(* packtes Feld der physikalischen Dimension n*m \"ubergeben *)
(* werden. *)
(*----------------------------------------------------------------*)
PROCEDURE Invert(VAR Ainv : ARRAY OF ARRAY OF LONGREAL; (* Invertiertes A *)
VAR A : ARRAY OF ARRAY OF LONGREAL; (* Zu invertierend *)
VAR Det : LONGREAL; (* Determinante der Matrix *)
dim : CARDINAL; (* Dimension der Matrix *)
VAR iFehl : INTEGER); (* Fehlercode *)
(*----------------------------------------------------------------*)
(* Berechnet die Inverse zu der Matrix A, wobei A durch die *)
(* Inverse \"uberschrieben wird, wenn der Parameter Ainv mit A *)
(* identisch ist. *)
(* Fehlermeldung, wenn A singul\"ar, die globale Variable Fehler *)
(* wird dann auf TRUE gesetzt. *)
(* *)
(* iFehl : 0 - kein Fehler *)
(* -1 - Matrix A ist (numerisch) singulart *)
(* 1 - Matrix A ist (numerisch) schlekt konditioniert *)
(* 2 - Speicherproblem - Abbruch *)
(* *)
(* Wenn iFehl eingangsseitig auf den Wert 99 gesetzt wird werden *)
(* eventuell Fehlermeldungen auf dem Fehlerkanal ausgegeben. *)
(* *)
(* Alternative Routine: InvertMat *)
(*----------------------------------------------------------------*)
PROCEDURE InvertMat(VAR Ainv : ARRAY OF ARRAY OF LONGREAL;
VAR A : ARRAY OF ARRAY OF LONGREAL;
VAR det : LONGREAL;
dim : CARDINAL);
(*----------------------------------------------------------------*)
(* Berechnet das Inverse der Matrix A und gibt dieses in Ainv *)
(* zurueck. Routine stuetzt sich auf LibDBlas und LinPack ab *)
(*----------------------------------------------------------------*)
PROCEDURE InvertSV( N : CARDINAL;
VAR A : ARRAY OF LONGREAL; (* SUPERVEKTOR *)
VAR iFehl : INTEGER);
(*----------------------------------------------------------------*)
(* Gauss-Jordan algorithm for the in situ inversion of a *)
(* positive definite matrix. Only the lower triangle is *)
(* transformed which is stored row by row in the vector *)
(* A[0:n*(n+1)/2-1]. *)
(* *)
(* iFehl : -1 wenn A singulaer oder nicht positiv definit ist *)
(* *)
(* Bauer/Reinsch, pp. 47 im "Handbuch", Alogol60 routine gidef2 *)
(*----------------------------------------------------------------*)
PROCEDURE Jacobi(VAR A : MATRIX; (* <== *)
VAR X : VEKTOR; (* <==> *)
VAR C : ARRAY OF LONGREAL; (* <== *)
dim : CARDINAL;
Tol : LONGREAL;
Neu : CARDINAL;
MaxIter : CARDINAL;
VAR iFehl : INTEGER);
(*----------------------------------------------------------------*)
(* Loest den Satz linearer Gleichungssystem A X = C nach der *)
(* Jacobi-Methode. *)
(* *)
(* Tol : Genauigkeitsangabe f��r *)
(* $\delta = \sum_{i=1}^dim \| X_i^k - X_i^{k-1} \|$ *)
(* der gefundenen L"osungsvektoren zweier aufeinander *)
(* folgenden Iterationen k,k-1. Wenn $\delta < Tol$ *)
(* gilt die L"osung als gefunden *)
(* Neu: *)
(* 1 : Der L��sungsvektor wird neu mit Nullen initialisiert *)
(* 0 : Der L��sungsvektor wird von der rufenden Prouzedur *)
(* "ubernommen. *)
(* MaxIter: *)
(* 0 : Die Routine setzt MaxIter auf 8*(dim**3) *)
(* >0 : Die Maximalzahl der Iterationen wird durch den Auf- *)
(* rufenden Prozess bestimmt. *)
(* iFehl: *)
(* 0 : Es ist kein Fehler aufgetreten. *)
(* 1 : Matrix ist nicht diagonal dominant, aber es wurde *)
(* eine L��sung gefunden (h\"oherer Zeitaufwand). *)
(* 2 : Die Matrix ist diagonal dominat aber MaxIter *)
(* wurde "uberschritten. Die L"osung kann ungenau sein. *)
(* 3 : Matrix ist nicht diagonal dominant und MaxIter *)
(* wurde "uberschritten. *)
(* 4 : Der L"oesung divergierte ueber die letzten zwei *)
(* Iterationen. *)
(* 5 : Bei der Berechnung des L"oesungsvektor droht ein *)
(* Ueberlauf - die L"oeung ist nicht brauchbar. *)
(* 6 : Ein Diagonalelement von A ist Null - es kann keine *)
(* L"osung berechnet werden. *)
(* F"ur iFehl > 0 wird die globale Variable Errors.Fehlerflag *)
(* mit einer entsprechenden Fehlermeldung belegt. *)
(* F"ur alle Fehlercodes gr"o\3er 3 wurde kein L"osung gefunden. *)
(*----------------------------------------------------------------*)
PROCEDURE GaussSeidel(VAR A : MATRIX; (* Koeffizinten - Matrix *)
VAR X : VEKTOR; (* L\"osungsvektor *)
VAR C : VEKTOR;
dim : CARDINAL;
genau : LONGREAL; (* Abbruchkriterium *)
MaxIter : CARDINAL;
Neu : CARDINAL;
FehlMeld : BOOLEAN; (* Fehlermeldungen j/n 1/0 *)
VAR Iter : CARDINAL;
VAR iFehl : INTEGER);
(*----------------------------------------------------------------*)
(* L\"o\3t das (diagonal dominante) lineare Gleichungssystem *)
(* A X = C nach dem (iterativen) Gauss - Seidel Algorithmus. *)
(* (E - L)x^{k+1} = C + Ux^k mit A = E-U-L. *)
(* Neu = *)
(* 0 : Der Startvektor f\"ur X wird vorgegeben. *)
(* 1 : Der Startvektor f\"ur X wird neu berechnet. *)
(* iFehl = *)
(* 0 : Es ist kein Fehler aufgetreten. *)
(* 1 : Es ist kein Fehler aufgetreten, die Matrix ist aber *)
(* nicht diagonal dominant. (h\"oherer Zeitaufwand) *)
(* 2 : MaxIter ist \"uberschritten *)
(* 3 : Problem ist mit dieser Routine nicht l\"o\3bar. *)
(* 6 : Ein Diagonalelement von A ist Null - es kann keine *)
(* L"osung berechnet werden. *)
(* In allen F\"allen, in denen FehlZahl > 0 werden auch die *)
(* globalen Variablen Fehler und Fehlerflag entsprechend *)
(* gesetzt. Wenn FehlZahl = 1 kann die Routine mit (Neu = 0) *)
(* wieder gerufen werden. *)
(* FehlMeld = *)
(* FALSE : Fehler werden nicht auf dem Bildschirm angezeigt. *)
(* TRUE : Alle Fehler werden auf dem Bildschirm angezeigt. *)
(* *)
(* Hinweis: Fuer ein ueberdeterminiertes Gleichungssystem *)
(* konvergiert der Algorithmus zu eine nicht eindeutigen *)
(* Loesung. Diese Situation kann nicht erkannt werden. *)
(* *)
(* Lit.: A.Jenninger, 'Matrix Computations for Engineers and *)
(* Sientists', J.Wiley & Sons, London 1977. *)
(* S.J.Leon, 'Linear Algebra', Maxmillian Pub. Co., *)
(* New York 1980. *)
(*----------------------------------------------------------------*)
PROCEDURE GaussJordan( N,M : CARDINAL;
VAR A : ARRAY OF ARRAY OF LONGREAL;
VAR B : ARRAY OF ARRAY OF LONGREAL;
inv : BOOLEAN;
VAR Det : LONGREAL;
VAR iFehl : INTEGER);
(*----------------------------------------------------------------*)
(* Solving a linear matrix system Ax = B with Gauss-Jordan method *)
(* using full pivoting at each step. During the process, original *)
(* matrices A and B are destroyed *)
(* *)
(* On Input: *)
(* *)
(* N : size of matrix A *)
(* M : number of simultaneous equations to be solved *)
(* A : matrix of size NxN *)
(* B : matrix of size NxM *)
(* inv : if true only the inverse of A is calculated *)
(* *)
(* On Output: *)
(* *)
(* A : inverse of A *)
(* Det : determinant of A *)
(* B : solution matrix of size NxM *)
(* *)
(* Based on a Pascal version provided in tpmath library *)
(*----------------------------------------------------------------*)
PROCEDURE Cholesky(VAR A : ARRAY OF LONGREAL; (* SUPERVEKTOR *)
VAR X : ARRAY OF LONGREAL; (* L\"osungsvektor *)
C : ARRAY OF LONGREAL;
dim : CARDINAL;
neu : BOOLEAN;
VAR Det : LONGREAL; (* Deteriminate von A *)
VAR iFehl : INTEGER);
(*----------------------------------------------------------------*)
(* Berechnet die L\"osung X eines linearen Gleichungsystem *)
(* A X = C mit symmetrischer, positiv definiten Koeffizienten - *)
(* matrix A in Superlinearform. *)
(* Ist A nicht positiv definit, so werden die globale Variabel *)
(* Fehler und Fehlerflag entsprechend gesetzt und auf Det wird *)
(* MAX(LONGREAL) zur\"uckgegeben. *)
(* *)
(* A : Koeffizientenmatrix, ausgangsseitig die Cholesky- *)
(* zerlegung von A *)
(* C : Koeffizientenvektor (rechte Seite des Gl.Sy.) *)
(* X : Loesungsvektor *)
(* dim : Dimension des Problems *)
(* neu : Wenn wahr wird die Zerlegung durchgefuehrt und der *)
(* der Loesungsvektor X berechnet *)
(* : Wenn falsch wird die Zerlegung nicht durchgefuehrt *)
(* und nur ein neuer Loesungsvektor berechnet. Dazu *)
(* muss zumindest ein Aufruf mit neu=TRUE vorausge- *)
(* gangen sein und A nicht mehr veraendert worden sein *)
(* Det : Determinate der Matrix A *)
(* iFehl : iFehl = 0 - keine Fehler *)
(* iFehl > 0 - A ist (eventuell aus numerischen *)
(* Gruenden, nicht positiv definit. iFehl *)
(* ist dann gleich dem Index der Gleich- *)
(* Gleichung bei der der Fehler aufge- *)
(* treten ist *)
(* iFehl = -1 - Dimensionierungsfehler *)
(*----------------------------------------------------------------*)
PROCEDURE CholDet1( N : CARDINAL;
VAR A : ARRAY OF ARRAY OF LONGREAL;
VAR p : ARRAY OF LONGREAL;
VAR d1 : LONGREAL;
VAR d2 : INTEGER;
VAR iFehl : INTEGER);
(*----------------------------------------------------------------*)
(* The upper triangle of a positive definite symmetrie matrix, *)
(* A, is stored in the upper triangle of an NxN array *)
(* A[i,j], i= 1 (1) n,j= 1 (1) n. *)
(* The Cholesky deeomposition A =LU, where U is the transpose of *)
(* L, is performed and stored in the remainder of the array A *)
(* except for the reciproeals of the diagonal elements which are *)
(* stored in p[i], i = 1 (1) n, instead of the elements *)
(* themselves. A is retained so that the solution obtained can *)
(* be subsequently improved. The determinant, *)
(* d1 * 2 ^ d2, of A is also computed. The proeedure will fail if *)
(* A, modified by the rounding errors, is not positive definite. *)
(* In that case iFehl will be equal to the index of the equation *)
(* where the error occured, 0 otherwise *)
(* *)
(* "Handbuch" routine choldet1 *)
(*----------------------------------------------------------------*)
PROCEDURE CholSol1( N : CARDINAL;
r : CARDINAL;
VAR A : ARRAY OF ARRAY OF LONGREAL;
VAR p : ARRAY OF LONGREAL;
VAR B : ARRAY OF ARRAY OF LONGREAL;
VAR X : ARRAY OF ARRAY OF LONGREAL);
(*----------------------------------------------------------------*)
(* Solves A X = B, where A is a positive definite symmetrie *)
(* matrix and B is an n x r matrix of r right-hand sides. The *)
(* proeedure cholsol1 must be preeeded by choldetl in which L is *)
(* produced in A[i,j] and p[i], from A. AX=B is solved in two *)
(* steps, Ly=B and UX=y. The matrix B is retained in order to *)
(* facilitate the refinement of X, but X is overwritten on y. *)
(* However, X and B can be identified in the call of the *)
(* procedure *)
(* *)
(* "Handbuch" routine cholsol1 *)
(*----------------------------------------------------------------*)
PROCEDURE CZerlege(VAR A : CMATRIX; (* Koeffizientenmatrix *)
VAR Index : CARDVEKTOR;
VAR Det : LONGCOMPLEX;
VAR EDet : INTEGER; (* Det(A) = Det*2^EDet *)
dim : CARDINAL); (* Dimension der qadr. Matrix *)
(*----------------------------------------------------------------*)
(* Zerlegt die komplexe Matrix A. (LU = A) mit Zeilenpivotierung *)
(* Vertauschungen werden in Index angezeigt *)
(* Fehlermeldungen wenn A singul\"ar oder unl\"osbar. *)
(*----------------------------------------------------------------*)
PROCEDURE CBerechneLoesung(VAR A : CMATRIX;
VAR X : CVEKTOR;
C : CVEKTOR;
Index : CARDVEKTOR;
dim : CARDINAL);
(*----------------------------------------------------------------*)
(* Berechnet den L\"osungsvektor X der Gleichung A x = c mit *)
(* der in der CZerlege LU-zerlegten Matrix A. *)
(* Die Routine ist das komplexe Analogon zu BerechneLoesung. *)
(*----------------------------------------------------------------*)
PROCEDURE CDet(dim : CARDINAL; VAR A : CMATRIX) : LONGCOMPLEX;
(*----------------------------------------------------------------*)
(* Berechnet die komplexe Determinante von A. *)
(* Die Koeffizientenmatrix A wird nach M\"oglichkeit nicht zer- *)
(* st\"ort, andernfalls wird die globale Variable Fehler auf TRUE *)
(* gesetzt und eine Fehlermeldung ausgeschrieben. *)
(*----------------------------------------------------------------*)
PROCEDURE CLoeseGlSy(VAR A : CMATRIX;
VAR X : CVEKTOR; (* L\"osungsvektor. *)
C : CVEKTOR;
VAR Det : LONGCOMPLEX; (* Determinte von A *)
dim : CARDINAL;
VAR iFehl : INTEGER);
(*----------------------------------------------------------------*)
(* Berechnet die L\"osung des komplexen Gleichngssystems Ax = c *)
(* Wenn die Matrix singul\"ar oder sonstwie nicht l\"osbar *)
(* ist, wird dir globale Variable Fehler auf TRUE gesetzt. *)
(* *)
(* iFehl : 0 wenn alles in Ordnung *)
(* : 1 wenn dim < 1 *)
(* : 2 wenn die Koeffizientenmatrix (ev. aus numerischen *)
(* Gruenden) singulaer ist *)
(*----------------------------------------------------------------*)
END LinLib.