IMPLEMENTATION MODULE LinLib;
(*------------------------------------------------------------------------*)
(* Modul zur L"osung linearer Gleichungssystem. *)
(* Modul for solving lineare systems of equations. *)
(*------------------------------------------------------------------------*)
(* Testroutine TestLinSol mit 15 Testf"allen f端r Gauss,Householder, *)
(* Jacobi,GaussSeidel in trlgssym.ein vorhanden. Tests fuer Cholesky, *)
(* CholDet1 und CholSol1 in TstCholesky. Fuer GaussJordan in TstGaussJ. *)
(* Fuer CDet,CLoeseGlSys in Tests in TstCmplxLL *)
(*------------------------------------------------------------------------*)
(* Letzte Bearbeitung: *)
(* *)
(* 10.01.97, MRi: Durchsicht *)
(* 12.09.15, MRi: iFehl Parameter in GaussSeidel auf CARDINAL umgestellt *)
(* 11.09.15, MRi: Prozedure Jacobi (29.06.96, aus remez.mod) integtiert. *)
(* 13.09.15, MRi: Prozedure Jacobi hefig ueberarbeitet, einige Sicher- *)
(* heitsabfragen aktualisiert/hinzugefuegt. *)
(* 14.09.15, MRi: In Prozedur GaussSeidel einige Sicherheitsabfragen *)
(* hinzugefuegt (Diagonalelement oder Loesungsvektor- *)
(* element Null. Konstantenvektor C als VAR deklariert. *)
(* 13.10.15, MRi: Ersetzen von LFLOAT(#) durch VAL(LONGREAL,#) *)
(* 16.12.16, MRi: In Cholesky Vorzeichen des Loesungsvektors geaendert *)
(* 23.12.16, MRi: In Invert zusaetzlichen Test auf Singularitaet eingef. *)
(* und alten Test "heruntergestuft" - bekannte Testfaelle *)
(* werden jetzt richtig behandelt. *)
(* 21.12.16, MRi: Erstellen der ersten M2-Version von GaussJordan based *)
(* on a Pascal version provided in tpmath library, see *)
(* also Pascal version by Jean-Pierre Moreau *)
(* 27.12.16, MRi: Prozedur CondHager eingefuegt *)
(* 28.12.16, MRI: Zerlege & BerechneLoesung nach aussen sichtbar gemacht *)
(* 21.05.17, MRI: Invert auf offene Felder umgestellt, InvertMat erstellt *)
(* 26.07.17, MRi: Erstellen der ersten Version von InvertSV *)
(* 28.07.17, MRi: Cholesky auf offene Felder umgestellt *)
(* 30.07.17, MRi: Erstellen der ersten Version von CholDet1 & CholSol1 *)
(* 01.08.17, MRi: Umstellen von GaussJordan auf offene Felder *)
(* 02.08.17, MRi: Fehlerbehandlung in Invert und InverMat verbessert *)
(* 24.09.17, MRi: Umbenennen von LONGCOMPLEX in LONGCMPLX sowie *)
(* CVEKTOR in CRVEKTOR und CMATRIX in CRMATRIX *)
(* 03.10.17, MRi: CZerlege CBerechneLoesung,CLoeseGlSy,CDet auf *)
(* ISO COMPLEX umgestellt *)
(* 05.10.17, MRi: Pivotierungsproblem in CZerlege,CBerecneLoesung *)
(* durch Einfuehren der Routine CVecPerm und echter Ver- *)
(* tauschung des Indexvektors geloesst *)
(* 29.08.18, MRi: Zerlege,Det,Gauss,BerechneLoesung, Householder, Jacobi *)
(* und GaussSeidel auf offene Felder umgestellt *)
(* 30.08.18, MRi: CZerlege,CDet,CBerechneLoesung und CLoeseGlSy auf *)
(* offene Felder umgestellt *)
(* 09.09.18, MRi: In BerechneLoesung Parameter "C" auf VAR gesetzt *)
(*------------------------------------------------------------------------*)
(* Offene Punkte: *)
(* *)
(* - Die Routinen Housholder und PHouseholder sind von der Indizierung *)
(* her eine Katastrophe - bitte aendern so dass A gestuerzt vorgegeben *)
(* werden muss. *)
(*------------------------------------------------------------------------*)
(* *)
(* 13.09.15: Routine GaussSeidel fuehrt bei einigen nicht diagonal *)
(* dominanten Matrizen zu Abbruechen - testen. *)
(* 14.09.15: Probleme in GaussSeidel wenn Diagonalelement Null oder *)
(* Loesungsvektorelement Null werden nun abgefangen. Die Werte *)
(* von iFehl sollten noch durchgesehen werden - eventuell an *)
(* Jacobi-Routine anpassen. Es wird h"aufig auf iFehl=1 gesetzt *)
(* (keine diagonale Dominanz) obwohl keine vern"unftige L"osung *)
(* berechnet wurde. *)
(* - Einfuehren der Konstante "Divergent" bringt hier Besserung, *)
(* allerdings nicht f端r alle Testf"alle. *)
(* Tests die durchfallen sind 7 & 11 (singul"are Matrix) *)
(* ---> Jacobi meldet bei Test 7 nur "MaxIter "uberschritten" - *)
(* hier k旦nnte die Divergenz-Konstante auch Verbesserungen *)
(* bringen, bitte testen. *)
(* *)
(* Test Jacobi GaussSeidel *)
(* *)
(* 1 2 (> MaxIter) OK (#) 0 OK *)
(* 2 0 OK 0 OK *)
(* 3 1 (!DiagDom) OK 1 (!DiagDom) OK *)
(* 4 4 (divergent) OK 3 (divergent) OK *)
(* 5 4 (divergent) OK 3 (divergent) OK *)
(* 6 6 (A_ii = 0 ) OK 6 (A_ii = 0 ) OK *)
(* war 7 ?? 3 (> MaxIter) ==> 1 (!DiagDom ) ?? *)
(* ist 7 3 (divergent) OK ==> 1 (!DiagDom ) ?? *)
(* 8 4 (divergent) OK 3 (divergent) OK *)
(* 9 0 OK 0 OK *)
(* 10 4 (divergent) OK 3 (divergent) OK *)
(* singul"ar 11 4 (divergent) OK ==> 1 (!DiagDom) ?? *)
(* 12 0 OK 0 OK *)
(* 13 4 (divergent) OK 3 (divergent) OK *)
(* 14 4 (divergent) OK 3 (divergent) OK *)
(* 15 0 OK 0 OK *)
(* *)
(* (#) N"aherungsl"osung OK *)
(* *)
(* - Prozedur Invert liefert f端r Test 6 X=1,2,3 Det=2 *)
(* 0.0000 1.0000 1.0000 5.0000 *)
(* 3.0000 0.0000 1.0000 6.0000 *)
(* -1.0000 1.0000 0.0000 1.0000 *)
(* falsche Ergebnisse - wieso ??? Die Prozedur meldet die Matrix sei *)
(* singulaer - was nicht stimmt. *)
(* *)
(* WICHTIG !!! *)
(* *)
(* - Alle Routinen auf richtige Indizierung ueberpruefen !!! *)
(* Siehe Hinweise in ~/Modula-2/log.txt *)
(* *)
(* - Routine Invert berechnet einige Inverse nicht obwohl die Matrix *)
(* invertierbar ist und bricht bei einer singulaeren Maxtrix nicht ab *)
(* Testfaelle in jacobi.ein *)
(* *)
(*------------------------------------------------------------------------*)
(* Implementation : Michael Riedl *)
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
(* $Id: LinLib.mod,v 1.14 2018/09/09 21:24:42 mriedl Exp mriedl $ *)
FROM SYSTEM IMPORT TSIZE;
FROM Storage IMPORT ALLOCATE,DEALLOCATE;
FROM Deklera IMPORT PVEKTOR,PMATRIX;
FROM Errors IMPORT Fehler,Fehlerflag,ErrOut,FatalError;
FROM LongMath IMPORT sqrt;
FROM LMathLib IMPORT MaxCard,CardPot,MachEps,Small;
IMPORT LongComplexMath;
FROM MatLib IMPORT MatNormL1;
FROM LibDBlas IMPORT dcopy,idamax;
FROM LinPack IMPORT dgefa,dgedi,dgesl;
FROM DynMat IMPORT AllocMat,DeAllocMat;
FROM SortLib IMPORT CVecPerm;
FROM StringsLib IMPORT Append;
CONST zero = LongComplexMath.zero;
one = LongComplexMath.one;
scalarMult = LongComplexMath.scalarMult;
CABS = LongComplexMath.abs;
CONST Ueberlauf = VAL(LONGREAL,MAX(REAL))*VAL(LONGREAL,MAX(REAL));
Unterlauf = 1.0 / Ueberlauf;
CONST MAXINT = MAX(INTEGER);
PROCEDURE sqr(x : LONGREAL) : LONGREAL; BEGIN RETURN x*x; END sqr;
PROCEDURE CondHager(VAR A : ARRAY OF ARRAY OF LONGREAL;
N : INTEGER;
VAR Cond : LONGREAL;
VAR iFehl : INTEGER);
(*-----------------------------------------------------------------*)
(* This procedure is based on a Fortran version by John Burkardt *)
(*-----------------------------------------------------------------*)
VAR Atmp : PMATRIX;
b : POINTER TO ARRAY [0..MAXINT-1] OF LONGREAL;
pivot : POINTER TO ARRAY [0..MAXINT-1] OF INTEGER;
c1,c2,Norm : LONGREAL;
i,j,i1,i2,info : INTEGER;
BEGIN
iFehl := 0;
AllocMat(Atmp,N,N);
ALLOCATE(b ,N*TSIZE(LONGREAL));
ALLOCATE(pivot,N*TSIZE(INTEGER));
IF (Atmp = NIL) OR (b = NIL) OR (pivot = NIL) THEN
IF (Atmp # NIL) THEN DeAllocMat(Atmp,N,N); END;
IF (b # NIL) THEN DEALLOCATE(b,N*TSIZE(LONGREAL)); END;
Cond:=MAX(LONGREAL); iFehl:=2;
RETURN;
END;
i1 := -1; (* MAX(INTEGER); *)
c1 := 0.0;
c2 := 0.0;
FOR i:=0 TO N-1 DO (* Save matrix A on Atmp *)
FOR j:=0 TO N-1 DO Atmp^[i+1]^[j+1] := A[i,j]; END;
END;
(* factor the matrix *)
dgefa(A,N,N,pivot^,info);
IF (info = 0) THEN
FOR i:=0 TO N-1 DO b^[i] := 1.0 / FLOAT(N); END;
LOOP
dgesl(A,N,N,pivot^,b^,0);
c2 := 0.0;
FOR i:=0 TO N-1 DO c2:=c2 + ABS(b^[i]); END;
FOR i:=0 TO N-1 DO
IF (b^[i] >= 0.0) THEN b^[i]:=1.0; ELSE b^[i]:=-1.0; END;
END;
dgesl(A,N,N,pivot^,b^,1);
i2 := idamax(N,b^[0],1) - 1;
IF (i1 >= 0) AND((i1 = i2) OR (c2 <= c1)) THEN
EXIT;
END;
i1 := i2;
c1 := c2;
FOR i:=0 TO N-1 DO b^[i]:=0.0; END; b^[i1]:=1.0;
END; (* LOOP *)
ELSE
iFehl:=1;
Cond:=MAX(LONGREAL);
END;
FOR i:=0 TO N-1 DO
FOR j:=0 TO N-1 DO A[i,j]:=Atmp^[i+1]^[j+1]; END;
END;
IF (iFehl = 0) THEN
MatNormL1(A,N,N,Norm);
Cond := c2*Norm;
END;
DeAllocMat(Atmp,N,N);
DEALLOCATE(pivot,N*TSIZE(INTEGER));
DEALLOCATE(b ,N*TSIZE(LONGREAL));
END CondHager;
PROCEDURE Zerlege(VAR A : ARRAY OF ARRAY OF LONGREAL;
VAR IndexI : ARRAY OF CARDINAL;
VAR IndexJ : ARRAY OF CARDINAL;
VAR RPivot : ARRAY OF LONGREAL;
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. *)
(*----------------------------------------------------------------*)
VAR i,j,k,imax,jmax,ii : CARDINAL;
Vorzeichen : INTEGER;
Max,Aik,Zw : LONGREAL;
Pivot : POINTER TO ARRAY [0..MAXINT-1] OF LONGREAL;
BEGIN
ALLOCATE(Pivot,dim*TSIZE(LONGREAL));
FOR k:=0 TO dim-1 DO IndexI[k]:=k; IndexJ[k]:=k; END;
Vorzeichen:=1; Det:=1.0; Fehler:=FALSE;
FOR k:=0 TO dim-1 DO
Max:=ABS(A[k,k]); imax:=k; jmax:=k;
IF MaxPivot THEN (* Finde Pivot *)
FOR i:=k TO dim-1 DO
FOR j:=k TO dim-1 DO
IF (ABS(A[i,j]) > Max) THEN
Max:=ABS(A[i,j]); imax:=i; jmax:=j;
END;
END; (* FOR j *)
END; (* FOR i *)
ELSE (* Zeilenpivotierung *)
FOR i:=k+1 TO dim-1 DO
IF (ABS(A[i,k]) > Max) THEN Max:=ABS(A[i,k]); imax:=i; END;
END;
END; (* IF MaxPivot *)
IF (imax # k) THEN (* Vertausche Zeile k,imax *)
Vorzeichen:= - Vorzeichen;
FOR j:=k TO dim-1 DO
Zw:=A[k,j]; A[k,j]:=A[imax,j]; A[imax,j]:=Zw;
END; (* 'Merke' Zeilenvertauschen *)
IndexI[k]:=imax;
END;
IF (jmax # k) (* AND MaxPivot *) THEN (* Vertausche Spalte k,jmax *)
Vorzeichen:= - Vorzeichen;
FOR i:=0 TO dim-1 DO (* Vertausch A_{j,k},A_{j,jmax]} *)
Zw:=A[i,k]; A[i,k]:=A[i,jmax]; A[i,jmax]:=Zw;
END; (* 'Merke' Spaltenvertauschen *)
ii:=IndexJ[k]; IndexJ[k]:=IndexJ[jmax]; IndexJ[jmax]:=ii;
END;
IF (A[k,k] = 0.0) THEN (* Sicherung ! *)
Fehler:=TRUE; Fehlerflag:='Matrix singul"ar (LinLib.Zerlege).';
Det:=0.0; RETURN;
END;
Pivot^[k]:=A[k,k]; RPivot[k]:=1.0 / Pivot^[k];
FOR i:=k+1 TO dim-1 DO A[k,i]:=A[k,i]*RPivot[k]; END;
A[k,k]:=1.0;
FOR i:=k+1 TO dim-1 DO
Aik:=A[i,k];
FOR j:=k+1 TO dim-1 DO A[i,j]:=A[i,j] - Aik*A[k,j]; END;
(*A[i,k]:=0.0;*)
END; (* FOR i *)
END; (* FOR k *)
Det:=Pivot^[0]; FOR k:=1 TO dim-1 DO Det:=Det*Pivot^[k]; END;
Det:=Det*VAL(LONGREAL,Vorzeichen);
DEALLOCATE(Pivot,dim*TSIZE(LONGREAL));
END Zerlege;
PROCEDURE Det(VAR A : ARRAY OF ARRAY OF LONGREAL;
dim : CARDINAL) : LONGREAL;
VAR i,j : CARDINAL;
Ii,Ij : POINTER TO ARRAY [0..MAXINT-1] OF CARDINAL;
Det : LONGREAL;
RPivot : POINTER TO ARRAY [0..MAXINT-1] OF LONGREAL;
Az : PMATRIX;
BEGIN
IF (dim = 1) THEN RETURN A[0,0]; END;
IF (dim < 1) THEN
Fehler:=TRUE; Fehlerflag:='Dimension der Matrix < 1 (LinLib.Det)';
ErrOut(Fehlerflag); RETURN MAX(LONGREAL);
END;
AllocMat(Az,dim,dim);
IF (Az = NIL) THEN
Fehler:=TRUE; Fehlerflag:='Kein Freispeicher vorhanden (LinLib.Det).';
ErrOut(Fehlerflag); RETURN MAX(LONGREAL);
END;
ALLOCATE(Ii,dim*TSIZE(CARDINAL));
ALLOCATE(Ij,dim*TSIZE(CARDINAL));
ALLOCATE(RPivot,dim*TSIZE(LONGREAL));
FOR i:=0 TO dim-1 DO
FOR j:=0 TO dim-1 DO Az^[i+1]^[j+1]:=A[i,j]; END;
END;
Fehler:=FALSE;
Zerlege(A,Ii^,Ij^,RPivot^,Det,dim,TRUE);
FOR i:=1 TO dim DO FOR j:=1 TO dim DO A[i-1,j-1]:=Az^[i]^[j]; END; END;
DeAllocMat(Az,dim,dim);
DEALLOCATE(Ii,dim*TSIZE(CARDINAL));
DEALLOCATE(Ij,dim*TSIZE(CARDINAL));
DEALLOCATE(RPivot,dim*TSIZE(LONGREAL));
RETURN Det;
END Det;
PROCEDURE BerechneLoesung(VAR A : ARRAY OF ARRAY OF LONGREAL;
VAR X : ARRAY OF LONGREAL;
VAR C : ARRAY OF LONGREAL;
VAR RPivot : ARRAY OF LONGREAL;
VAR IndexI : ARRAY OF CARDINAL;
VAR IndexJ : ARRAY OF CARDINAL;
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. *)
(*----------------------------------------------------------------*)
VAR k,i,j : CARDINAL;
Zw : LONGREAL;
Z,Ctmp : POINTER TO ARRAY [0..MAXINT-1] OF LONGREAL;
BEGIN
ALLOCATE(Z, dim*TSIZE(LONGREAL));
ALLOCATE(Ctmp,dim*TSIZE(LONGREAL));
FOR k:=0 TO dim-1 DO Ctmp^[k]:=C[k]; END;
FOR k:=0 TO dim-1 DO
IF (IndexI[k] # k) THEN
Zw:=Ctmp^[IndexI[k]]; Ctmp^[IndexI[k]]:=Ctmp^[k]; Ctmp^[k]:=Zw;
END;
Ctmp^[k]:=Ctmp^[k]*RPivot[k]; (* Teile durch Pivot *)
FOR i:=k+1 TO dim-1 DO Ctmp^[i]:=Ctmp^[i] - A[i,k]*Ctmp^[k]; END;
END; (* FOR k *)
Z^[dim-1]:=Ctmp^[dim-1]; (* Berechne die L\"osungen X[i] *)
FOR k:=dim-2 TO 0 BY -1 DO
Z^[k]:=Ctmp^[k]; FOR j:=k+1 TO dim DO Z^[k]:=Z^[k] - A[k,j]*Z^[j]; END;
END; (* FOR k *)
FOR j:=0 TO dim-1 DO X[IndexJ[j]]:=Z^[j]; END;
DEALLOCATE(Z, dim*TSIZE(LONGREAL));
DEALLOCATE(Ctmp,dim*TSIZE(LONGREAL));
END BerechneLoesung;
PROCEDURE Gauss(VAR A : ARRAY OF ARRAY OF LONGREAL;
VAR X : ARRAY OF LONGREAL;
C : ARRAY OF LONGREAL;
dim : CARDINAL;
VAR Det : LONGREAL; (* Determinante von A *)
MaxIter : CARDINAL;
PivStrat : CARDINAL; (* Pivotstrategie *)
VAR iFehl : INTEGER); (* Fehlernummer *)
CONST eps = 100.0*MachEps;
VAR Iter,i,j : CARDINAL;
d0,d1,XMax,BMax : LONGREAL;
B,RPivot : POINTER TO ARRAY [0..MAXINT-1] OF LONGREAL;
IndexI,IndexJ : POINTER TO ARRAY [0..MAXINT-1] OF CARDINAL;
MaxPivot : BOOLEAN;
Az : PMATRIX;
BEGIN
IF (dim < 1) THEN
Fehler:=TRUE; Fehlerflag:='dim < 1 (LinLib.Gauss)';
ErrOut(Fehlerflag); iFehl:=4;
Det:=MAX(LONGREAL); RETURN;
END;
MaxPivot := (PivStrat = 1);
Fehler:=FALSE; iFehl:=0;
AllocMat(Az,dim,dim);
ALLOCATE(IndexI,dim*TSIZE(CARDINAL));
ALLOCATE(IndexJ,dim*TSIZE(CARDINAL));
ALLOCATE(B,dim*TSIZE(LONGREAL));
ALLOCATE(RPivot,dim*TSIZE(LONGREAL));
IF (Az # NIL) THEN
FOR i:=1 TO dim DO
FOR j:=1 TO dim DO Az^[i]^[j]:=A[i-1,j-1]; END;
END;
ELSE
Fehler:=TRUE; Fehlerflag:='kein Freispeicher vorhanden (Gauss).';
iFehl:=1;
END;
Zerlege(A,IndexI^,IndexJ^,RPivot^,Det,dim,MaxPivot);
IF Fehler THEN
iFehl:=3; Fehlerflag:='Matrix singulaer (LinLib.Gauss) !';
ErrOut(Fehlerflag);
ELSE
BerechneLoesung(A,X,C,RPivot^,IndexI^,IndexJ^,dim);
Iter:=1; d0:=0.0;
LOOP
INC(Iter); IF (Iter > MaxIter) THEN EXIT; END;
FOR i:=1 TO dim DO (* B = Az X *)
B^[i-1]:=0.0;
FOR j:=1 TO dim DO B^[i-1]:=B^[i-1] + Az^[i]^[j]*X[j-1]; END;
END;
(* Residuenvektor *)
FOR i:=1 TO dim DO B^[i-1] := B^[i-1] - C[i-1]; END;
BerechneLoesung(A,B^,B^,RPivot^,IndexI^,IndexJ^,dim);
FOR i:=1 TO dim DO X[i-1]:=X[i-1] - B^[i-1]; END;
XMax:=ABS(X[0]); BMax:=ABS(B^[0]);
FOR i:=2 TO dim DO
IF (ABS(X [i-1]) > XMax) THEN XMax:=ABS(X [i-1]); END;
IF (ABS(B^[i-1]) > BMax) THEN BMax:=ABS(B^[i-1]); END;
END;
IF (BMax < eps*XMax) THEN EXIT; END; (* !!! *)
d1:=BMax / XMax;
IF (Iter > 2) AND (d1 > 0.5*d0) THEN
Fehlerflag:='Loesung nicht konvergent (LinLib.Gauss).';
Fehler:=TRUE; iFehl:=2; EXIT; (* !!! *)
END;
d0:=d1;
END; (* LOOP *)
END; (* IF Fehler *)
IF (Az # NIL) THEN (* Speicher altes A zur\"uck *)
FOR i:=1 TO dim DO FOR j:=1 TO dim DO A[i-1,j-1]:=Az^[i]^[j]; END; END;
DeAllocMat(Az,dim,dim);
END;
DEALLOCATE(IndexI,dim*TSIZE(CARDINAL));
DEALLOCATE(IndexJ,dim*TSIZE(CARDINAL));
DEALLOCATE(B,dim*TSIZE(LONGREAL));
DEALLOCATE(RPivot,dim*TSIZE(LONGREAL));
END Gauss;
PROCEDURE Householder(VAR A : ARRAY OF ARRAY OF LONGREAL;
VAR X : ARRAY OF LONGREAL;
C : ARRAY OF LONGREAL;
n : CARDINAL; (* Anzahl der Gleichungen, n >= m *)
m : CARDINAL); (* Anzahl der Unbekannten *)
VAR Sum,s,Xi : LONGREAL;
Alpha,Beta : LONGREAL;
i,j,k : CARDINAL;
D : POINTER TO ARRAY [0..MAXINT-1] OF LONGREAL;
BEGIN
IF (n < m) THEN
FatalError('Zuwenige Gleichungen (n < m) (LinLib.Householder).');
END;
ALLOCATE(D,m*TSIZE(LONGREAL));
Fehler:=FALSE;
FOR j:=0 TO m-1 DO
Sum:=0.0;
FOR i:=j TO n-1 DO Sum:=Sum + A[i,j]*A[i,j]; END;
IF (Sum = 0.0) THEN
Fehler:=TRUE; Fehlerflag:='Matrix singul"ar (LinLib.Householder)';
RETURN; (* !!! *)
END;
IF (A[j,j] < 0.0) THEN s:=sqrt(Sum) ELSE s:=-sqrt(Sum); END;
(* Waehle s so dass s*Ajj negativ ist *)
D^[j]:=s;
Alpha:=s*A[j,j] - Sum;
IF (Alpha = 0.0) THEN
Fehler:=TRUE; Fehlerflag:='Matrix singul"ar (LinLib.Householder)';
RETURN; (* !!! *)
ELSE
Beta:=1.0 / Alpha;
END;
A[j,j]:=A[j,j] - s;
FOR k:=j+1 TO m-1 DO
Sum:=0.0;
FOR i:=j TO n-1 DO Sum:=Sum + A[i,j]*A[i,k]; END;
Sum:=Beta*Sum;
FOR i:=j TO n-1 DO A[i,k]:=A[i,k] + A[i,j]*Sum; END;
END;
Sum:=0.0;
FOR i:=j TO n-1 DO Sum:=Sum + A[i,j]*C[i]; END;
Sum:=Beta*Sum;
FOR i:=j TO n-1 DO C[i]:=C[i] + A[i,j]*Sum; END;
END; (* FOR j *)
FOR i:=m-1 TO 0 BY -1 DO
Xi:=C[i];
FOR j:=m-1 TO i+1 BY -1 DO Xi:=Xi - A[i,j]*X[j]; END;
X[i]:=Xi / D^[i];
END;
DEALLOCATE(D,m*TSIZE(LONGREAL));
END Householder;
PROCEDURE PHouseholder(VAR A : PMATRIX;
VAR X : ARRAY OF LONGREAL;
VAR C : ARRAY OF LONGREAL;
m : CARDINAL; (* Anzahl der Gleichungen, n >= m *)
n : CARDINAL); (* Anzahl der Unbekannten *)
(* C und X d"urfen nicht auf denselben Speicherplatz verweisen ! *)
VAR Sum,s,Alpha,Beta,Xi : LONGREAL;
i,j,k : CARDINAL;
D,Ctemp : PVEKTOR;
BEGIN
IF (m < n) THEN
ErrOut('Zuwenige Gleichungen (m < n) (PHouseholder).');
END;
ALLOCATE(D , n*TSIZE(LONGREAL));
ALLOCATE(Ctemp,MaxCard(m,n)*TSIZE(LONGREAL));
IF (D = NIL) OR (Ctemp = NIL) THEN
FatalError("Kein Freispeicher vorhanden (PHouseholder) !");
END;
FOR i:=1 TO m DO Ctemp^[i]:=C[i-1]; END;
Fehler:=FALSE;
FOR j:=1 TO n DO
Sum:=0.0;
FOR i:=j TO m DO Sum:=Sum + A^[i]^[j]*A^[i]^[j]; END;
IF (Sum = 0.0) THEN
Fehler:=TRUE;
Fehlerflag:='Matrix singul"ar (PHouseholder)';
RETURN; (* !!! *)
END;
IF (A^[j]^[j] < 0.0) THEN s:=sqrt(Sum) ELSE s:=-sqrt(Sum); END;
D^[j]:=s;
Alpha:=s*A^[j]^[j] - Sum;
IF (Alpha = 0.0) THEN
Fehler:=TRUE;
Fehlerflag:='Matrix singul"ar (PHouseholder)';
RETURN; (* !!! *)
ELSE
Beta:=1.0 / Alpha;
END;
A^[j]^[j]:=A^[j]^[j] - s;
FOR k:=j+1 TO n DO
Sum:=0.0;
FOR i:=j TO m DO Sum:=Sum + A^[i]^[j]*A^[i]^[k]; END;
Sum:=Beta*Sum;
FOR i:=j TO m DO A^[i]^[k]:=A^[i]^[k] + A^[i]^[j]*Sum; END;
END;
Sum:=0.0;
FOR i:=j TO m DO Sum:=Sum + A^[i]^[j]*Ctemp^[i]; END;
Sum:=Beta*Sum;
FOR i:=j TO m DO Ctemp^[i]:=Ctemp^[i] + A^[i]^[j]*Sum; END;
END; (* FOR j *)
FOR i:=n TO 1 BY -1 DO
Xi:=Ctemp^[i];
FOR j:=n TO i+1 BY -1 DO Xi:=Xi - A^[i]^[j]*X[j-1]; END;
X[i-1]:=Xi / D^[i];
END;
DEALLOCATE(D,n*TSIZE(LONGREAL));
DEALLOCATE(Ctemp,MaxCard(m,n)*TSIZE(LONGREAL));
END PHouseholder;
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 *)
PROCEDURE IndexIJ(i,j : CARDINAL) : CARDINAL;
(* FUSCH: IndexIJ ist sehr ineffizient !!! *)
BEGIN
RETURN (i-1)*m + j - 1;
END IndexIJ;
VAR Sum,s,Xi : LONGREAL;
Alpha,Beta : LONGREAL;
i,j,k,jj,ij,ik : CARDINAL;
D,Cloc : POINTER TO ARRAY [0..MAXINT-1] OF LONGREAL;
BEGIN
IF (n < m) THEN
FatalError('Zuwenige Gleichungen (n < m) (LinLib.HhLin).');
END;
ALLOCATE(D ,n*TSIZE(LONGREAL));
ALLOCATE(Cloc,n*TSIZE(LONGREAL));
IF (Cloc = NIL) OR (D = NIL) THEN
FatalError('Kein Freispeicher vorhanden (HhLin) !');
END;
FOR i:=0 TO n-1 DO Cloc^[i]:=C[i]; END;
Fehler:=FALSE;
FOR j:=1 TO m DO
jj:=IndexIJ(j,j);
Sum:=0.0;
FOR i:=j TO n DO ij:=IndexIJ(i,j); Sum:=Sum + A[ij]*A[ij]; END;
IF (Sum = 0.0) THEN
Fehler:=TRUE; Fehlerflag:='Matrix singul"ar (LinLib.HhLin)';
RETURN; (* !!! *)
END;
IF (A[jj] < 0.0) THEN s:=sqrt(Sum) ELSE s:=-sqrt(Sum); END;
D^[j]:=s;
Alpha:=s*A[jj] - Sum;
IF (Alpha = 0.0) THEN
Fehler:=TRUE; Fehlerflag:='Matrix singul"ar (LinLib.HhLin)';
RETURN; (* !!! *)
ELSE
Beta:=1.0 / Alpha;
END;
A[jj]:=A[jj] - s;
FOR k:=j+1 TO m DO
Sum:=0.0;
FOR i:=j TO n DO Sum:=Sum + A[IndexIJ(i,j)]*A[IndexIJ(i,k)]; END;
Sum:=Beta*Sum;
FOR i:=j TO n DO
ik:=IndexIJ(i,k); A[ik]:=A[ik] + A[IndexIJ(i,j)]*Sum;
END;
END;
Sum:=0.0;
FOR i:=j TO n DO Sum:=Sum + A[IndexIJ(i,j)]*Cloc^[i-1]; END;
Sum:=Beta*Sum;
FOR i:=j TO n DO Cloc^[i-1]:=Cloc^[i-1] + A[IndexIJ(i,j)]*Sum; END;
END; (* FOR j *)
FOR i:=m TO 1 BY -1 DO
Xi:=Cloc^[i-1];
FOR j:=m TO i+1 BY -1 DO Xi:=Xi - A[IndexIJ(i,j)]*X[j-1]; END;
X[i-1]:=Xi / D^[i];
END;
DEALLOCATE(Cloc,n*TSIZE(LONGREAL));
DEALLOCATE(D ,n*TSIZE(LONGREAL));
END HhLin;
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);
VAR i,j,k,pk : CARDINAL;
P : POINTER TO ARRAY [0..MAXINT-1] OF CARDINAL;
RecPivot,Max : LONGREAL;
Sum,q,Zw : LONGREAL;
Melden : BOOLEAN;
BEGIN
Melden := iFehl = 99;
Fehler:=FALSE; Fehlerflag:='OK'; iFehl := 0;
ALLOCATE(P,dim*TSIZE(CARDINAL));
IF (P = NIL) THEN
Fehler:=TRUE;
Fehlerflag:='Kein Freispeicher verfuegbar (LinLib.Invert)';
IF Melden THEN ErrOut(Fehlerflag); END;
iFehl := 2;
RETURN;
END;
Det:=1.0;
FOR i:=0 TO dim-1 DO FOR j:=0 TO dim-1 DO Ainv[i,j]:=A[i,j]; END; END;
FOR k:=0 TO dim-1 DO
Max:=0.0;
P^[k]:=0; (* ??? *)
FOR i:=k TO dim-1 DO
Sum:=0.0;
FOR j:=k TO dim-1 DO Sum:=Sum + ABS(Ainv[i,j]); END;
IF (ABS(Sum) < Unterlauf*Unterlauf) THEN (* 02.08.17 *)
Fehler:=TRUE;
Fehlerflag:='Matrix singulaer (LinLib.Invert).';
IF Melden THEN ErrOut(Fehlerflag); END;
DEALLOCATE(P,dim*TSIZE(CARDINAL));
iFehl := -1; Det:=0.0;
RETURN;
END;
q:=ABS(Ainv[i,k]) / Sum;
IF (q > Max) THEN Max:=q; P^[k]:=i; END;
IF (ABS(Max) < Unterlauf*Unterlauf) AND (iFehl = 0) THEN
(* Fehler:=TRUE; *)
Fehlerflag:='Matrix schlecht konditioniert (LinLib.Invert).';
IF Melden THEN ErrOut(Fehlerflag); END;
iFehl := 1;
END;
END; (* FOR i *)
pk:=P^[k];
IF (pk # k) THEN
FOR j:=0 TO dim-1 DO
Zw:=Ainv[k,j]; Ainv[k,j]:=Ainv[pk,j]; Ainv[pk,j]:=Zw;
END;
END; (* IF *)
Det:=Det*Ainv[k,k];
IF (ABS(Det) < Unterlauf) THEN (* 23.12.16 *)
Fehler:=TRUE; Fehlerflag:='Matrix singulaer (LinLib.Invert).';
IF Melden THEN ErrOut(Fehlerflag); END;
DEALLOCATE(P,dim*TSIZE(CARDINAL));
iFehl := -1; Det:=0.0;
RETURN;
END;
RecPivot:=1.0 / Ainv[k,k];
FOR j:=0 TO dim-1 DO
IF (j # k) THEN
Ainv[k,j]:=-Ainv[k,j]*RecPivot;
IF (k > 0) THEN
FOR i:=0 TO k-1 DO
Ainv[i,j]:=Ainv[i,j] + Ainv[i,k]*Ainv[k,j];
END;
END;
FOR i:=k+1 TO dim-1 DO
Ainv[i,j]:=Ainv[i,j] + Ainv[i,k]*Ainv[k,j];
END;
END; (* IF *)
END; (* FOR j *)
FOR i:=0 TO dim-1 DO Ainv[i,k]:=Ainv[i,k]*RecPivot; END;
Ainv[k,k]:=RecPivot;
END; (* FOR k *)
FOR k:=dim-2 TO 0 BY -1 DO
pk:=P^[k];
IF (pk # k) THEN
Det:=-Det;
FOR i:=0 TO dim-1 DO
Zw:=Ainv[i,k]; Ainv[i,k]:=Ainv[i,pk]; Ainv[i,pk]:=Zw;
END;
END;
END; (* FOR k *)
DEALLOCATE(P,dim*TSIZE(CARDINAL));
IF (ABS(Det) < Unterlauf*Unterlauf) THEN (* 23.12.16 *)
Fehler:=TRUE; Fehlerflag:='Matrix singulaer (LinLib.Invert).';
IF Melden THEN ErrOut(Fehlerflag); END;
iFehl := -1;
RETURN;
END;
END Invert;
PROCEDURE InvertMat(VAR Ainv : ARRAY OF ARRAY OF LONGREAL;
VAR A : ARRAY OF ARRAY OF LONGREAL;
VAR det : LONGREAL;
dim : CARDINAL);
VAR i : CARDINAL;
n,Ndim : INTEGER;
info : INTEGER;
ipvt : POINTER TO ARRAY [0..MAXINT-1] OF INTEGER;
work : POINTER TO ARRAY [0..MAXINT-1] OF LONGREAL;
Det : ARRAY [0..2] OF LONGREAL;
BEGIN
Fehler:=FALSE; Fehlerflag:='OK';
ALLOCATE(ipvt,dim*TSIZE(INTEGER));
IF (ipvt = NIL) THEN
FatalError("Kein Freispeicher vorhanden (LinLib.InvertMat)");
END;
n:=VAL(INTEGER,dim);
FOR i:=0 TO dim-1 DO (* Copy A to Ainv *)
dcopy(n,A[i][0],1,Ainv[i][0],1);
END;
Ndim:=VAL(INTEGER,HIGH(Ainv[0])+1);
dgefa(Ainv,Ndim,n,ipvt^,info);
IF (info = 0) THEN
ALLOCATE(work,dim*TSIZE(LONGREAL));
IF (work = NIL) THEN
DEALLOCATE(ipvt,dim*TSIZE(INTEGER));
FatalError("Kein Freispeicher vorhanden (LinLib.InvertMat)");
END;
dgedi(Ainv,Ndim,n,ipvt^,Det,work^,11);
det:=Det[2];
DEALLOCATE(work,dim*TSIZE(LONGREAL));
ELSE
Fehler:=TRUE; Fehlerflag:="Matrix singulaer (LinLib.InvertMat)";
det:=0.0;
END;
DEALLOCATE(ipvt,dim*TSIZE(INTEGER));
END InvertMat;
PROCEDURE InvertSV( N : CARDINAL;
VAR A : ARRAY OF LONGREAL; (* SUPERVEKTOR *)
VAR iFehl : INTEGER);
VAR i,ii,ij,k,m : CARDINAL;
p,q : LONGREAL;
h : POINTER TO ARRAY [1..MAXINT] OF LONGREAL;
BEGIN
ALLOCATE(h,N*TSIZE(LONGREAL));
iFehl:=0;
k:=N+1;
LOOP
DEC(k);
IF (k = 0) THEN EXIT; END; (* FOR k:=N TO 1 BY -1 DO *)
p := A[0];
IF (p <= 0.0) THEN
iFehl := -1;
EXIT;
ELSE
ii:=1; m:=1;
FOR i:=2 TO N DO
m := ii;
ii := ii + i; (* m =i*(i -1)/2, ii =i*(i +1)/2 *)
q := A[m];
h^[i] := -q / p;
IF (i > k) THEN h^[i] := -h^[i]; END;
FOR ij:=(m+1) TO ii-1 DO
A[ij-i] := A[ij] + q*h^[ij-m+1];
END;
END;
DEC(m); (* m=n*(n-1)/2-1, ii=n*(n+1)/2 *)
A[ii-1] := 1.0 / p;
FOR i:=2 TO N DO A[m+i-1] := h^[i]; END;
END;
END;
DEALLOCATE(h,N*TSIZE(LONGREAL));
END InvertSV;
PROCEDURE Jacobi(VAR A : ARRAY OF ARRAY OF LONGREAL; (* <== *)
VAR X : ARRAY OF LONGREAL; (* <==> *)
VAR C : ARRAY OF LONGREAL; (* <== *)
dim : CARDINAL;
Tol : LONGREAL;
Neu : CARDINAL;
MaxIter : CARDINAL;
VAR iFehl : INTEGER);
PROCEDURE IstUeberlauf(VAR V : ARRAY OF LONGREAL;
dim : CARDINAL) : BOOLEAN;
(* Ermittel das absolut gr"osste Element des L"osungs- *)
(* vektors. Ist dieses gr"oesser als Ueberlauf wird *)
(* TRUE zur"uckgegeben *)
VAR i : CARDINAL; Max : LONGREAL;
BEGIN
Max:=ABS(V[0]);
FOR i:=1 TO dim-1 DO
IF (ABS(V[i]) > Max) THEN Max:=ABS(V[i]); END;
END;
IF (Max > Ueberlauf) THEN RETURN TRUE; ELSE RETURN FALSE; END;
END IstUeberlauf;
PROCEDURE TestDiagDom(VAR A : ARRAY OF ARRAY OF LONGREAL;
dim : CARDINAL;
VAR DiagDom : BOOLEAN);
(* Teste, ob die Matrix A diagonal dominant ist *)
VAR i,j : CARDINAL;
Sum : LONGREAL;
BEGIN
DiagDom:=TRUE; i:=0;
REPEAT
Sum:=0.0;
IF (i > 0) THEN
FOR j:=0 TO i-1 DO Sum:=Sum + ABS(A[i,j]); END;
END;
FOR j:=i+1 TO dim-1 DO Sum:=Sum + ABS(A[i,j]); END;
IF (Sum > ABS(A[i,i])) THEN
DiagDom:=FALSE;
END;
INC(i);
UNTIL (i = dim) OR NOT DiagDom;
END TestDiagDom;
VAR Xalt : POINTER TO ARRAY [0..MAXINT-1] OF LONGREAL;
Iter,i,j : CARDINAL;
Delta,Sum : LONGREAL;
DeltaAlt1,DeltaAlt2 : LONGREAL;
DiagDom,ok,divergent : BOOLEAN;
BEGIN
iFehl:=0; Fehlerflag:="";
FOR i:=0 TO dim-1 DO
IF (ABS(A[i,i]) < Unterlauf) THEN
Fehlerflag:="Diagonalelement ist 0 (LinLib.Jacobi)";
iFehl:=6;
RETURN;
END;
END;
ALLOCATE(Xalt,dim*TSIZE(LONGREAL));
TestDiagDom(A,dim,DiagDom);
IF (MaxIter = 0) THEN MaxIter:=8*dim*dim*dim; END;
IF (Tol = 0.0) THEN Tol:=VAL(LONGREAL,dim)*10.0*MachEps; END;
IF (Neu = 1) THEN
FOR i:=0 TO dim-1 DO X[i]:=0.0; END; (* Initialize Starting Value *)
END;
DeltaAlt1:=Ueberlauf; DeltaAlt2:=MAX(LONGREAL);
divergent:=FALSE;
Iter:=0;
REPEAT (* Jacobiiterationen *)
INC(Iter);
FOR i:=0 TO dim-1 DO Xalt^[i]:=X[i]; END;
FOR i:=0 TO dim-1 DO
Sum := C[i];
IF (i > 0) THEN
FOR j:=0 TO i-1 DO Sum:=Sum - A[i,j]*Xalt^[j]; END;
END;
FOR j:=i+1 TO dim-1 DO Sum:=Sum - A[i,j]*Xalt^[j]; END;
X[i] := Sum / A[i,i];
END;
Delta:=0.0;
FOR i:=0 TO dim-1 DO Delta:=Delta + ABS(X[i] - Xalt^[i]); END;
IF (Iter > 8*dim) THEN (* Teste erst nach einigen Iterationen *)
divergent:=(Delta > DeltaAlt1) AND (Delta > DeltaAlt2);
END;
DeltaAlt2:=DeltaAlt1; DeltaAlt1:=Delta;
ok := Delta <= Tol;
UNTIL ok OR (Iter > MaxIter) OR divergent OR IstUeberlauf(X,dim);
DEALLOCATE(Xalt,dim*TSIZE(LONGREAL));
IF (Iter <= MaxIter) THEN (* Fehlerauswertung *)
IF NOT ok THEN
IF divergent THEN
Fehlerflag:="Loesungs divergiert (LinLib.Jacobi) !";
iFehl:=4;
ELSE (* IstUeberlauf = TRUE *)
Fehlerflag:="Ueberlauf (LinLib.Jacobi) !";
iFehl:=5;
END;
ELSIF NOT DiagDom THEN
Fehlerflag:='Matrix nicht diagonal dominant (LinLib.Jacobi)';
iFehl:=1;
END;
ELSE
Fehlerflag:="MaxIter 端berschritten (LinLib.Jacobi) !";
IF DiagDom THEN iFehl:=2; ELSE iFehl:=3; END;
END;
IF (iFehl # 0) THEN
ErrOut(Fehlerflag);
END;
END Jacobi;
PROCEDURE GaussSeidel(VAR A : ARRAY OF ARRAY OF LONGREAL;
VAR X : ARRAY OF LONGREAL;
VAR C : ARRAY OF LONGREAL;
dim : CARDINAL;
genau : LONGREAL;
MaxIter : CARDINAL;
Neu : CARDINAL;
FehlMeld : BOOLEAN;
VAR Iter : CARDINAL;
VAR iFehl : INTEGER);
CONST eps = 100.0*MachEps;
Divergent = 1.0E+08; (* Festlegung reine Willk"ur *)
VAR i,j : CARDINAL;
Xalt,Z : POINTER TO ARRAY [0..MAXINT-1] OF LONGREAL;
Sum,Prod,RezDim : LONGREAL;
DiagDom,Konvergenz : BOOLEAN;
ok : BOOLEAN;
BEGIN
iFehl:=0;
FOR i:=0 TO dim-1 DO
IF (ABS(A[i,i]) < Unterlauf) THEN
Fehlerflag:="Diagonalelement ist 0 (LinLib.GaussSeidel)";
IF FehlMeld THEN ErrOut(Fehlerflag); END;
iFehl:=6;
RETURN;
END;
END;
IF (dim = 1) THEN X[0]:=C[0]/A[0,0]; Iter:=0; RETURN END;
IF (genau <= 0.0) THEN genau:=eps; END;
IF (MaxIter = 0) THEN MaxIter:=MaxCard(dim*dim,40); END;
IF (Neu > 1) THEN Neu:=1; END;
ALLOCATE(Xalt,dim*TSIZE(LONGREAL));
ALLOCATE(Z ,dim*TSIZE(LONGREAL));
DiagDom:=TRUE; i:=0; (* Teste, ob Matrix diagonal dominant. *)
REPEAT
Sum:=0.0;
IF (i > 0) THEN
FOR j:=0 TO i-1 DO Sum:=Sum + ABS(A[i,j]); END;
END;
FOR j:=i+1 TO dim-1 DO Sum:=Sum + ABS(A[i,j]); END;
IF (Sum > ABS(A[i,i])) THEN
Fehlerflag:='Matrix nicht diagonal dominant (LinLib.GaussSeidel)';
DiagDom:=FALSE;
END;
INC(i);
UNTIL (i = dim) OR NOT DiagDom;
Sum:=ABS(A[0,0]); (* Test, ob Gleichungssystem l\"osbar. *)
FOR i:=1 TO dim-1 DO Sum:=Sum + ABS(A[i,i]); END;
RezDim:=1.0 / VAL(LONGREAL,dim); Sum:=Sum*RezDim;
FOR i:=0 TO dim-1 DO
IF (ABS(A[i,i]) < genau*Sum) THEN
Fehler:=TRUE; iFehl:=3;
Fehlerflag:='Gleichungssystem nicht loesbar (LinLib.GaussSeidel)';
IF FehlMeld THEN ErrOut(Fehlerflag); END;
RETURN;
END;
END;
IF (Neu = 1) THEN (* Neuer Startvektor *)
FOR i:=0 TO dim-1 DO
IF (ABS(A[i,i]) < Unterlauf) THEN (* 14.09.15 *)
X[i]:=MachEps; Xalt^[i]:=MachEps;
ELSE
X[i]:=A[i,i]; Xalt^[i]:=X[i]; END;
END;
ELSE (* Speicher den Startvektor in Z *)
FOR i:=0 TO dim-1 DO
Z^[i]:=X[i];
IF (ABS(X[i]) < Unterlauf) THEN (* 14.09.15 *)
X[i]:=MachEps; Xalt^[i]:=MachEps;
ELSE
Z^[i]:=X[i]; Xalt^[i]:=X[i];
END;
END;
END;
Prod:=1.0; Iter:=0;
LOOP (* Gauss - Seidel Iterationen. *)
IF (Iter < MaxIter) THEN INC(Iter); ELSE EXIT END; (* !!! *)
FOR i:=0 TO dim-1 DO
Sum:=0.0;
IF (i > 0) THEN
FOR j:=0 TO i-1 DO Sum:=Sum + A[i,j]*X[j]; END;
END;
FOR j:=i+1 TO dim-1 DO Sum:=Sum + A[i,j]*Xalt^[j]; END; (* ?? *)
X[i]:=(C[i] - Sum) / A[i,i];
END; (* FOR i *)
i:=0; (* Teste auf Konvergenz *)
REPEAT
Konvergenz:=(ABS(Xalt^[i] - X[i]) < genau*ABS(X[i]));
INC(i);
UNTIL NOT Konvergenz OR (i = dim);
IF Konvergenz THEN EXIT END; (* !!! *)
IF NOT DiagDom THEN (* Teste auf Divergenz ! *)
i:=0;
REPEAT
ok := (ABS(Xalt^[i]) > Unterlauf); INC(i);
UNTIL (i = dim) OR NOT ok;
IF ok THEN
Sum:=0.0;
FOR i:=0 TO dim-1 DO Sum:=Sum + ABS(X[i] / Xalt^[i]); END;
Prod:=Prod*Sum*RezDim;
END;
IF (Prod > Divergent) OR NOT ok THEN
Fehler:=TRUE; iFehl:=4;
Fehlerflag:=
'Matrix nicht diagonal dominant & Loesungfolge divergiert';
Append(Fehlerflag,' (LinLib.GaussSeidel)');
EXIT;
END;
END; (* IF *)
FOR i:=0 TO dim-1 DO Xalt^[i]:=X[i]; END;
END; (* LOOP *)
IF (iFehl = 4) THEN
IF (Neu = 0) THEN
FOR j:=0 TO dim-1 DO X[j]:=Z^[j]; END;
END;
ELSIF (Iter <= MaxIter) THEN
IF DiagDom THEN iFehl:=0; ELSE iFehl:=1; END;
ELSE
IF NOT DiagDom THEN
Fehler:=TRUE; iFehl:=3;
Fehlerflag:='Gleichungssystem nicht loesbar (LinLib.GaussSeidel)';
IF (Neu = 0) THEN FOR i:=1 TO dim DO X[i]:=Z^[i]; END; END;
ELSE
Fehler:=TRUE; iFehl:=2;
Fehlerflag:='MaxIter "uberschritten (LinLib.GaussSeidel)';
END;
END;
IF (iFehl # 0) AND FehlMeld THEN ErrOut(Fehlerflag); END;
DEALLOCATE(Xalt,dim*TSIZE(LONGREAL));
DEALLOCATE(Z ,dim*TSIZE(LONGREAL));
END GaussSeidel;
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);
CONST eps = Small*Small;
MAXINT = MAX(INTEGER)-1;
VAR Pi,Pj : POINTER TO ARRAY [0..MAXINT] OF CARDINAL;
C : POINTER TO ARRAY [0..MAXINT] OF LONGREAL;
Pivot,tmp : LONGREAL;
i,j,k,ik,jk : CARDINAL;
BEGIN
ALLOCATE(Pi,N*TSIZE(CARDINAL));
ALLOCATE(Pj,N*TSIZE(CARDINAL));
ALLOCATE(C ,N*TSIZE(LONGREAL));
IF (Pi = NIL) OR (Pj = NIL) OR (C = NIL) THEN
FatalError("Kein Freispeicher vorhanden (LinLib.GaussJordan)");
END;
IF (inv) THEN M:=0; END;
FOR i:=0 TO N-1 DO (* Initializations *)
Pi^[i] := 0;
Pj^[i] := 0;
C^[i] := 0.0;
END;
Det := 1.0;
iFehl := 0;
FOR k:=0 TO N-1 DO (* Main loop *)
(* Searching greatest pivot in submatrix A[k..n,k..n] *)
Pivot:=A[k,k];
ik:=k;
jk:=k;
FOR i:=k TO N-1 DO
FOR j:=k TO N-1 DO
IF (ABS(A[i,j]) > ABS(Pivot)) THEN
Pivot := A[i,j];
ik:=i;
jk:=j;
END;
END;
END;
Pi^[k]:=jk; Pj^[k]:=ik; (* Store pivot position *)
(* update determinant *)
IF (ik # k) THEN Det:= - Det;END;
IF (jk # k) THEN Det:= - Det;END;
Det:=Det*Pivot;
IF (ABS(Det) < eps) THEN
(* Too weak pivot ==> quasi-singular matrix *)
Fehlerflag:='System badly conditioned ? (GaussJordan)';
Fehler:=TRUE; iFehl:=2;
END;
IF (ik # k) THEN
FOR i:=0 TO N-1 DO (* exchange current row k with pivot row ik *)
tmp := A[ik,i];
A[ik,i] := A[k,i];
A[k,i] := tmp;
END;
END;
IF (M > 0) THEN
FOR i:=0 TO M-1 DO
tmp := B[ik,i];
B[ik,i] := B[k,i];
B[k,i] := tmp;
END;
END;
IF (jk # k) THEN
(* Exchange current column (K) with pivot column (Jk) *)
FOR i:=0 TO N-1 DO
tmp := A[i,jk];
A[i,jk] := A[i,k];
A[i,k] := tmp;
END;
END;
(* Store column K of matrix A into C *)
(* and set this column to zero *)
FOR i:=0 TO N-1 DO
C^[i] := A[i,k];
A[i,k] := 0.0;
END;
C^[k] := 0.0;
A[k,k]:= 1.0; (* Line k of matrix A is modified:*)
IF (ABS(Pivot) < eps) THEN
Fehlerflag:='Pivot too small (GaussJordan)';
Fehler:=TRUE; iFehl:=2;
RETURN;
END;
(* Transform pivot row *)
tmp := 1.0 / Pivot;
FOR i:=0 TO N-1 DO
A[k,i]:=A[k,i]*tmp;
END;
IF (M > 0) THEN
FOR i:=0 TO M-1 DO
B[k,i]:=B[k,i]*tmp;
END;
END;
FOR j:=0 TO N-1 DO (* Transform other rows *)
IF (j # k) THEN
FOR i:=0 TO N-1 DO
A[j,i]:=A[j,i] - C^[j]*A[k,i];
END;
IF (M > 0) THEN
FOR i:=0 TO M-1 DO
B[j,i]:=B[j,i] - C^[j]*B[k,i];
END;
END;
END;
END;
END; (* k *)
FOR i:=N-1 TO 0 BY -1 DO (* Exchange lines of inverse matrix *)
ik:=Pi^[i];
IF (ik # i) THEN
FOR j:=0 TO N-1 DO
tmp := A[i,j];
A[i,j] := A[ik,j];
A[ik,j] := tmp;
END;
IF (M > 0) THEN
FOR j:=0 TO M-1 DO
tmp := B[i,j];
B[i,j] := B[ik,j];
B[ik,j]:= tmp;
END;
END;
END;
END; (* i *)
FOR j:=N-1 TO 0 BY -1 DO (* exchange columns of inverse matrix *)
jk:=Pj^[j];
IF (jk # j) THEN
FOR i:=0 TO N-1 DO
tmp := A[i,j];
A[i,j] := A[i,jk];
A[i,jk] := tmp;
END;
END;
END;
IF (ABS(Det) < Small) THEN
Fehlerflag:='determinant equals zero (GaussJordan)';
Fehler:=TRUE; iFehl:=1;
END;
DEALLOCATE(Pi,N*TSIZE(CARDINAL));
DEALLOCATE(Pj,N*TSIZE(CARDINAL));
DEALLOCATE(C ,N*TSIZE(LONGREAL));
END GaussJordan;
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);
CONST tiny = Small*Small;
VAR i,j,k,ii,kk,ij,ik,ki,jk : CARDINAL;
Sum : LONGREAL;
BEGIN
Fehler:=FALSE; iFehl:=0;
IF (dim = 0) OR (dim > HIGH(X)-1) THEN
iFehl:= -1; Fehler:=TRUE;
Fehlerflag:='Dimensionierungsfehler (LinLib.Cholesky)';
Det:=MAX(LONGREAL);
RETURN;
END;
IF neu THEN
Det:=1.0; kk:=0;
FOR k:=1 TO dim DO
INC(kk,k);
IF (A[kk-1] <= tiny) THEN
iFehl:=k; Fehler:=TRUE;
Fehlerflag:='Matrix nicht positiv definit (LinLib.Cholesky)';
ErrOut(Fehlerflag);
Det:=MAX(LONGREAL);
RETURN;
END;
Det:=Det*A[kk-1];
A[kk-1]:=sqrt(A[kk-1]);
ik:=kk+k-1;
FOR i:=k+1 TO dim DO
A[ik]:=A[ik] / A[kk-1];
ij:=ik; jk:=kk+k-1;
FOR j:=k+1 TO i DO
INC(ij); A[ij]:=A[ij] - A[ik]*A[jk]; INC(jk,j);
END;
INC(ik,i);
END;
END; (* FOR k *)
END; (* IF neu *)
C[0]:=C[0] / A[0];
ij:=1;
FOR i:=1 TO dim-1 DO
Sum:=C[i];
FOR j:=0 TO i-1 DO Sum:=Sum - A[ij]*C[j]; INC(ij); END;
C[i]:=Sum / A[ij];
INC(ij);
END;
ii:=(dim*(dim+1) DIV 2) - 1;
FOR i:=dim-1 TO 0 BY -1 DO
Sum:=C[i]; ki:=ii+i+1;
FOR k:=i+1 TO dim-1 DO Sum:=Sum + A[ki]*X[k]; INC(ki,k+1); END;
X[i] := -Sum / A[ii]; DEC(ii,i+1);
END;
FOR i:=0 TO dim-1 DO X[i]:=-X[i]; END;
END Cholesky;
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);
VAR i,j,k : CARDINAL;
x : LONGREAL;
BEGIN
iFehl:=0;
d1 := 1.0;
d2 := 0;
FOR i:=0 TO N-1 DO
FOR j:=i TO N-1 DO
x:=A[i,j];
IF (i > 0) THEN
FOR k:=i-1 TO 0 BY -1 DO
x:=x - A[j,k]*A[i,k];
END;
END;
IF (j = i) THEN
d1:=d1*x;
IF (x = 0.0) THEN
d2:= 0;
iFehl:= -1; RETURN;
END;
WHILE (ABS(d1) >= 1.0) DO
d1:=d1*0.0625;
d2:=d2 + 4;
END;
WHILE (ABS(d1) < 0.0625) DO
d1:=d1 * 16.0;
d2:=d2 - 4;
END;
IF (x <= 0.0) THEN iFehl:= -1; RETURN; END; (* war x < 0 *)
p[i] := 1.0 / sqrt (x)
ELSE
A[j,i]:= x*p[i];
END;
END;
END;
END 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);
VAR i,j,k : CARDINAL;
z : LONGREAL;
BEGIN
FOR j:=0 TO r-1 DO
(* solution of Ly=b *)
X[0,j]:=B[0,j]*p[0];
FOR i:=1 TO N-1 DO
z:=B[i,j];
FOR k:=i-1 TO 0 BY -1 DO
z:=z - A[i,k]*X[k,j];
END;
X[i,j]:=z*p[i];
END;
(* solution of U x= y *)
FOR i:=N-1 TO 0 BY -1 DO
z:=X[i,j];
FOR k:=i+1 TO N-1 DO
z:=z- A[k,i]*X[k,j];
END;
X[i,j] := z*p[i];
END;
END;
END CholSol1;
PROCEDURE CZerlege(VAR A : ARRAY OF ARRAY OF LONGCOMPLEX;
VAR Index : ARRAY OF CARDINAL;
VAR Det : LONGCOMPLEX;
VAR EDet : INTEGER; (* Det(A) = Det*2^EDet *)
dim : CARDINAL); (* Dimension der qadr. Matrix *)
CONST eps = VAL(LONGREAL,1.0 / MAX(REAL)); (* *MAX(REAL)); *)
VAR j,k,l,m,i : CARDINAL;
CSum,Zw : LONGCOMPLEX;
w,x,y,z,Sum : LONGREAL;
SumVek : POINTER TO ARRAY [0..MAXINT-1] OF LONGREAL;
BEGIN
ALLOCATE(SumVek,dim*TSIZE(LONGREAL));
Fehler:=FALSE; EDet:=0;
FOR k:=0 TO dim-1 DO
Sum:=0.0; Index[k]:=k+1; (* k+1 ??? *)
FOR l:=0 TO dim-1 DO Sum:=Sum + sqr(RE(A[k,l])) + sqr(IM(A[k,l])); END;
SumVek^[k]:=Sum;
END;
Det := CMPLX(1.0,0.0);
FOR k:=0 TO dim-1 DO
l:=k; z:=0.0;
FOR m:=k TO dim-1 DO
CSum :=A[m,k];
IF (k > 0) THEN
FOR j:=0 TO k-1 DO CSum:=CSum - A[m,j]*A[j,k]; END;
END;
A[m,k] := CSum;
x:=(sqr(RE(CSum)) + sqr(IM(CSum))) / SumVek^[m];
IF (x > z) THEN z:=x; l:=m; END;
END; (* FOR m *)
IF (l # k) THEN (* Vertausche Zeilen *)
Det := - Det;
FOR j:=0 TO dim-1 DO
Zw:=A[k,j]; A[k,j]:=A[l,j]; A[l,j]:=Zw;
END;
i :=Index[k]; Index[k]:=Index[l]; Index[l]:=i;
END; (* IF l # k *)
x:=RE(A[k,k]); y:=IM(A[k,k]);
z:=x*x + y*y;
IF (ABS(z) < eps) THEN
(* Wegen Det := Det * CMPLX(x,y) *)
Fehler:=TRUE; Fehlerflag:='Matrix singulaer (LinLib.CZerlege)';
Det := zero;
RETURN;
END;
Det:=Det*CMPLX(x,y);
IF (Det = zero) THEN (* 03.10.17 *)
Fehler:=TRUE; Fehlerflag:='Matrix singulaer (LinLib.CZerlege)';
DEALLOCATE(SumVek,dim*TSIZE(LONGREAL));
RETURN;
END;
IF (ABS(RE(Det)) > ABS(IM(Det))) THEN w:=RE(Det); ELSE w:=IM(Det); END;
WHILE (ABS(w) >= 1.0) DO
w:=0.0625*w; INC(EDet,4);
Det:=scalarMult(0.0625,Det);
END;
WHILE (ABS(w) < 0.0625) DO
w:=16.0*w; DEC(EDet,4);
Det:=scalarMult(16.0,Det);
END;
FOR l:=k+1 TO dim-1 DO
CSum := A[k,l];
IF (k > 0) THEN
FOR j:=0 TO k-1 DO CSum:=CSum - A[k,j]*A[j,l]; END;
END;
A[k,l]:= scalarMult( (1.0/z) ,CSum*CMPLX(x,-y));
END; (* FOR i *)
END; (* FOR k *)
DEALLOCATE(SumVek,dim*TSIZE(LONGREAL));
END CZerlege;
PROCEDURE CBerechneLoesung(VAR A : ARRAY OF ARRAY OF LONGCOMPLEX;
VAR X : ARRAY OF LONGCOMPLEX;
C : ARRAY OF LONGCOMPLEX;
VAR Index : ARRAY OF CARDINAL;
dim : CARDINAL);
VAR k,j : CARDINAL;
CSum : LONGCOMPLEX;
BEGIN
CVecPerm(Index,1,dim,C); (* Bringe C in die "richtige" Ordnung *)
FOR k:=0 TO dim-1 DO
C[k]:=C[k] / A[k,k]; (* Teile durch Pivot *)
FOR j:=k+1 TO dim-1 DO
C[j]:=C[j] - A[j,k]*C[k];
END;
END; (* FOR k *)
X[dim-1]:=C[dim-1];
FOR k:=dim-2 TO 0 BY -1 DO
CSum:=C[k];
FOR j:=k+1 TO dim-1 DO
CSum:=CSum - A[k,j]*X[j];
END;
X[k]:=CSum;
END;
END CBerechneLoesung;
PROCEDURE CDet( dim : CARDINAL;
VAR A : ARRAY OF ARRAY OF LONGCOMPLEX) : LONGCOMPLEX;
VAR Index : POINTER TO ARRAY [0..MAXINT-1] OF CARDINAL;
Det : LONGCOMPLEX;
EDet : INTEGER;
i,j : CARDINAL;
Dete : LONGREAL;
AzRe,AzIm : PMATRIX;
BEGIN
Fehler:=FALSE;
ALLOCATE(Index,dim*TSIZE(CARDINAL));
AllocMat(AzRe,dim,dim);
AllocMat(AzIm,dim,dim);
IF (AzRe # NIL) AND (AzIm = NIL) THEN
DeAllocMat(AzRe,dim,dim);
END;
IF (AzRe # NIL) AND (AzIm # NIL) THEN
FOR i:=1 TO dim DO
FOR j:=1 TO dim DO
AzRe^[i]^[j]:=RE(A[i-1,j-1]);
AzIm^[i]^[j]:=IM(A[i-1,j-1]);
END;
END;
ELSE
ErrOut(' Koeffizenzenmatrix zerstoert (LinLib.Det).');
END;
CZerlege(A,Index^,Det,EDet,dim); Fehler:=FALSE; Fehlerflag[0]:=0C;
Dete:=CardPot(2.0,VAL(CARDINAL,ABS(EDet)));
IF (EDet < 0) THEN Dete:=1.0 / Dete; END;
Det := scalarMult(Dete,Det);
IF (AzRe # NIL) AND (AzIm # NIL) THEN
FOR i:=1 TO dim DO
FOR j:=1 TO dim DO
A[i-1,j-1]:=CMPLX(AzRe^[i]^[j],AzIm^[i]^[j]);
END;
END;
DeAllocMat(AzIm,dim,dim);
DeAllocMat(AzRe,dim,dim);
END;
DEALLOCATE(Index,dim*TSIZE(CARDINAL));
RETURN Det;
END CDet;
PROCEDURE CLoeseGlSy(VAR A : ARRAY OF ARRAY OF LONGCOMPLEX;
VAR X : ARRAY OF LONGCOMPLEX;
C : ARRAY OF LONGCOMPLEX;
VAR Det : LONGCOMPLEX; (* Determinte von A *)
dim : CARDINAL;
VAR iFehl : INTEGER);
VAR Index : POINTER TO ARRAY [0..MAXINT-1] OF CARDINAL;
x : LONGREAL;
EDet : INTEGER;
BEGIN
IF (dim < 1) THEN
Fehler:=TRUE; Fehlerflag:='dim < 1 (CLoeseGlSy)';
iFehl := 1;
Det := zero;
RETURN;
END;
ALLOCATE(Index,dim*TSIZE(CARDINAL));
IF (Index = NIL) THEN
Fehler:=TRUE; Fehlerflag:='Kein Freispeicher vorhanden (CLoeseGlSy)';
iFehl:=3;
Det := zero;
ErrOut(Fehlerflag);
END;
Fehler:=FALSE; iFehl:=0;
CZerlege(A,Index^,Det,EDet,dim);
IF Fehler THEN
iFehl:=2;
ErrOut(Fehlerflag);
RETURN;
END; (* !!! *)
x:=CardPot(2.0,VAL(CARDINAL,ABS(EDet)));
IF (EDet < 0) THEN x:=1.0 / x; END;
Det := scalarMult(x,Det);
CBerechneLoesung(A,X,C,Index^,dim);
DEALLOCATE(Index,dim*TSIZE(CARDINAL));
END CLoeseGlSy;
END LinLib.