--- a
+++ b/LinLib.mod
@@ -0,0 +1,1478 @@
+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 tlgs.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 *)
+ (*------------------------------------------------------------------------*)
+ (* Offene Punkte: *)
+ (* *)
+ (* *)
+ (* 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.2 2017/10/04 06:24:25 mriedl Exp mriedl $ *)
+
+FROM SYSTEM IMPORT TSIZE;
+FROM Storage IMPORT ALLOCATE,DEALLOCATE;
+FROM Deklera IMPORT VEKTOR,MATRIX,CARDVEKTOR,PVEKTOR,PMATRIX,
+ CVEKTOR,CMATRIX;
+FROM Errors IMPORT Fehler,Fehlerflag,ErrOut,FatalError;
+FROM LongMath IMPORT sqrt;
+FROM LMathLib IMPORT MaxCard,CardPot,MachEps,Small;
+FROM LongComplexMath IMPORT zero,scalarMult;
+FROM MatLib IMPORT MatNormL1;
+FROM LibDBlas IMPORT dcopy,idamax;
+FROM LinPack IMPORT dgefa,dgedi,dgesl;
+FROM DynMat IMPORT AllocMat,DeAllocMat;
+FROM SortLib IMPORT CVecPerm;
+
+CONST Ueberlauf = VAL(LONGREAL,MAX(REAL))*VAL(LONGREAL,MAX(REAL));
+ Unterlauf = 1.0 / Ueberlauf;
+
+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 *)
+ (*-----------------------------------------------------------------*)
+
+ CONST MaxI1 = MAX(INTEGER)-1;
+
+ VAR Atmp : PMATRIX;
+ b : POINTER TO ARRAY [0..MaxI1] OF LONGREAL;
+ pivot : POINTER TO ARRAY [0..MaxI1] 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,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 : 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. *)
+ (*----------------------------------------------------------------*)
+
+ VAR i,j,k,imax,jmax,ii : CARDINAL;
+ Vorzeichen : INTEGER;
+ Max,Aik,Zw : LONGREAL;
+ Pivot : VEKTOR;
+BEGIN
+ FOR k:=1 TO dim DO IndexI[k]:=k; IndexJ[k]:=k; END;
+ Vorzeichen:=1; Det:=1.0; Fehler:=FALSE;
+ FOR k:=1 TO dim DO
+ Max:=ABS(A[k,k]); imax:=k; jmax:=k;
+ IF MaxPivot THEN (* Finde Pivot *)
+ FOR i:=k TO dim DO
+ FOR j:=k TO dim 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 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 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:=1 TO dim 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 DO A[k,i]:=A[k,i]*RPivot[k]; END;
+ A[k,k]:=1.0;
+ FOR i:=k+1 TO dim DO
+ Aik:=A[i,k];
+ FOR j:=k+1 TO dim 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[1]; FOR k:=2 TO dim DO Det:=Det*Pivot[k]; END;
+ Det:=Det*VAL(LONGREAL,Vorzeichen);
+END Zerlege;
+
+PROCEDURE Det(VAR A : MATRIX; dim : CARDINAL) : LONGREAL;
+
+ VAR i,j : CARDINAL;
+ Ii,Ij : CARDVEKTOR;
+ Det : LONGREAL;
+ RPivot : VEKTOR;
+ Az : POINTER TO MATRIX;
+BEGIN
+ IF (dim = 1) THEN RETURN A[1,1]; END;
+ IF (dim < 1) THEN
+ Fehler:=TRUE; Fehlerflag:='Dimension der Matrix < 1 (LinLib.Det)';
+ ErrOut(Fehlerflag); RETURN MAX(LONGREAL);
+ END;
+ NEW(Az);
+ IF (Az = NIL) THEN
+ Fehler:=TRUE; Fehlerflag:='Kein Freispeicher vorhanden (LinLib.Det).';
+ ErrOut(Fehlerflag); RETURN MAX(LONGREAL);
+ END;
+ FOR i:=1 TO dim DO FOR j:=1 TO dim DO Az^[i,j]:=A[i,j]; END; END;
+ Zerlege(A,Ii,Ij,RPivot,Det,dim,TRUE); Fehler:=FALSE;
+ FOR i:=1 TO dim DO FOR j:=1 TO dim DO A[i,j]:=Az^[i,j]; END; END;
+ DISPOSE(Az);
+ RETURN Det;
+END Det;
+
+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. *)
+ (*----------------------------------------------------------------*)
+
+ VAR k,i,j : CARDINAL;
+ Zw : LONGREAL;
+ Z : VEKTOR;
+BEGIN
+ FOR k:=1 TO dim DO
+ IF (IndexI[k] # k) THEN
+ Zw:=C[IndexI[k]]; C[IndexI[k]]:=C[k]; C[k]:=Zw;
+ END;
+ C[k]:=C[k]*RPivot[k]; (* Teile durch Pivot *)
+ FOR i:=k+1 TO dim DO C[i]:=C[i] - A[i,k]*C[k]; END;
+ END; (* FOR k *)
+
+ Z[dim]:=C[dim]; (* Berechne die L\"osungen X[i] *)
+ FOR k:=dim-1 TO 1 BY -1 DO
+ Z[k]:=C[k]; FOR j:=k+1 TO dim DO Z[k]:=Z[k] - A[k,j]*Z[j]; END;
+ END; (* FOR k *)
+ FOR j:=1 TO dim DO X[IndexJ[j]]:=Z[j]; END;
+END BerechneLoesung;
+
+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 *)
+
+ CONST eps = 100.0*MachEps;
+
+ VAR Iter,i,j : CARDINAL;
+ d0,d1,XMax,BMax : LONGREAL;
+ B,RPivot : VEKTOR;
+ IndexI,IndexJ : CARDVEKTOR;
+ MaxPivot : BOOLEAN;
+ Az : POINTER TO MATRIX;
+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;
+ NEW(Az);
+ IF (Az # NIL) THEN
+ FOR i:=1 TO dim DO FOR j:=1 TO dim DO Az^[i,j]:=A[i,j]; END; END;
+ ELSE
+ MaxIter:=1;
+ Fehler:=TRUE; Fehlerflag:='kein Freispeicher vorhanden (Gauss).';
+ iFehl:=1;
+ END;
+ Zerlege(A,IndexI,IndexJ,RPivot,Det,dim,MaxPivot);
+ IF Fehler THEN
+ IF (Az # NIL) THEN A:=Az^; DISPOSE(Az); END;
+ iFehl:=3; Fehlerflag:='Matrix singul"ar (LinLib.Gauss) !';
+ ErrOut(Fehlerflag);
+ RETURN; (* !!! *)
+ END;
+ 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]:=0.0; FOR j:=1 TO dim DO B[i]:=B[i] + Az^[i,j]*X[j]; END;
+ END;
+ FOR i:=1 TO dim DO B[i] := B[i] - C[i]; END; (* Residuenvektor *)
+ BerechneLoesung(A,B,B,RPivot,IndexI,IndexJ,dim);
+ FOR i:=1 TO dim DO X[i]:=X[i] - B[i]; END;
+ XMax:=ABS(X[1]); BMax:=ABS(B[1]);
+ FOR i:=2 TO dim DO
+ IF (ABS(X[i]) > XMax) THEN XMax:=ABS(X[i]); END;
+ IF (ABS(B[i]) > BMax) THEN BMax:=ABS(B[i]); END;
+ END;
+ IF (BMax < eps*XMax) THEN EXIT; END; (* !!! *)
+ d1:=BMax / XMax;
+ IF (Iter > 2) AND (d1 > 0.5*d0) THEN
+ Fehlerflag:='L"osung nicht konvergent (LinLib.Gauss).';
+ Fehler:=TRUE; iFehl:=2; EXIT; (* !!! *)
+ END;
+ d0:=d1;
+ END; (* LOOP *)
+ IF (Az # NIL) THEN (* Speicher altes A zur\"uck *)
+ FOR i:=1 TO dim DO FOR j:=1 TO dim DO A[i,j]:=Az^[i,j]; END; END;
+ DISPOSE(Az);
+ END;
+END Gauss;
+
+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 *)
+
+ VAR Sum,s,Alpha,Beta,Xi : LONGREAL;
+ i,j,k : CARDINAL;
+ D : VEKTOR;
+BEGIN
+ IF (n < m) THEN
+ FatalError('Zuwenige Gleichungen (n < m) (LinLib.Householder).');
+ END;
+ Fehler:=FALSE;
+ FOR j:=1 TO m DO
+ Sum:=0.0;
+ FOR i:=j TO n 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;
+ 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 DO
+ Sum:=0.0;
+ FOR i:=j TO n DO Sum:=Sum + A[i,j]*A[i,k]; END;
+ Sum:=Beta*Sum;
+ FOR i:=j TO n DO A[i,k]:=A[i,k] + A[i,j]*Sum; END;
+ END;
+ Sum:=0.0;
+ FOR i:=j TO n DO Sum:=Sum + A[i,j]*C[i-1]; END;
+ Sum:=Beta*Sum;
+ FOR i:=j TO n DO C[i-1]:=C[i-1] + A[i,j]*Sum; END;
+ END; (* FOR j *)
+ FOR i:=m TO 1 BY -1 DO
+ Xi:=C[i-1];
+ FOR j:=m TO i+1 BY -1 DO Xi:=Xi - A[i,j]*X[j-1]; END;
+ X[i-1]:=Xi / D[i];
+ END;
+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,Alpha,Beta,Xi : LONGREAL;
+ i,j,k,jj,ij,ik : CARDINAL;
+ D : ARRAY [1..256] OF LONGREAL;
+ Cloc : POINTER TO ARRAY [0..255] OF LONGREAL;
+BEGIN
+ IF (n < m) THEN
+ FatalError('Zuwenige Gleichungen (n < m) (LinLib.HhLin).');
+ END;
+ IF (m > 256) THEN
+ FatalError('Anzahl der Gleichungen gr\"o\3er 256 (HhLin) !');
+ END;
+ ALLOCATE(Cloc,n*TSIZE(LONGREAL));
+ IF (Cloc = 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));
+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..MAX(INTEGER)-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..MAX(INTEGER)-1] OF INTEGER;
+ work : POINTER TO ARRAY [0..MAX(INTEGER)-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..MAX(INTEGER)] 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 : MATRIX; (* <== *)
+ VAR X : VEKTOR; (* <==> *)
+ VAR C : ARRAY OF LONGREAL; (* <== *)
+ dim : CARDINAL;
+ Tol : LONGREAL;
+ Neu : CARDINAL;
+ MaxIter : CARDINAL;
+ VAR iFehl : INTEGER);
+
+ PROCEDURE IstUeberlauf(VAR V: VEKTOR; 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[1]);
+ FOR i:=2 TO dim 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 : MATRIX;
+ 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
+ INC(i);
+ Sum:=0.0;
+ FOR j:=1 TO i-1 DO Sum:=Sum + ABS(A[i,j]); END;
+ FOR j:=i+1 TO dim DO Sum:=Sum + ABS(A[i,j]); END;
+ IF (Sum > ABS(A[i,i])) THEN
+ DiagDom:=FALSE;
+ END;
+ UNTIL (i = dim) OR NOT DiagDom;
+ END TestDiagDom;
+
+ VAR Xalt : VEKTOR;
+ Iter,i,j : CARDINAL;
+ Delta,Sum : LONGREAL;
+ DeltaAlt1,DeltaAlt2 : LONGREAL;
+ DiagDom,ok,divergent : BOOLEAN;
+BEGIN
+ iFehl:=0; Fehlerflag:="";
+ FOR i:=1 TO dim DO
+ IF (ABS(A[i,i]) < Unterlauf) THEN
+ Fehlerflag:="Diagonalelement ist 0 (LinLib.Jacobi)";
+ iFehl:=6;
+ RETURN;
+ END;
+ END;
+
+ 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:=1 TO dim 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:=1 TO dim DO Xalt[i]:=X[i]; END;
+ FOR i:=1 TO dim DO
+ Sum := C[i-1];
+ FOR j:=1 TO i-1 DO Sum:=Sum - A[i,j]*Xalt[j]; END;
+ FOR j:=i+1 TO dim DO Sum:=Sum - A[i,j]*Xalt[j]; END;
+ X[i] := Sum / A[i,i];
+ END;
+ Delta:=0.0;
+ FOR i:=1 TO dim 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);
+
+ 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 : 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);
+
+ CONST eps = 100.0*MachEps;
+ Divergent = 1.0E+08; (* Festlegung reine Willk"ur *)
+
+ VAR i,j : CARDINAL;
+ Xalt,Z : VEKTOR;
+ Sum,Prod,RezDim : LONGREAL;
+ DiagDom,Konvergenz : BOOLEAN;
+ ok : BOOLEAN;
+BEGIN
+ FOR i:=1 TO dim 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[1]:=C[1]/A[1,1]; Iter:=0; iFehl:=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;
+
+ DiagDom:=TRUE; i:=0; (* Teste, ob Matrix diagonal dominant. *)
+ REPEAT
+ INC(i);
+ Sum:=0.0;
+ FOR j:=1 TO i-1 DO Sum:=Sum + ABS(A[i,j]); END;
+ FOR j:=i+1 TO dim 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;
+ UNTIL (i = dim) OR NOT DiagDom;
+
+ Sum:=ABS(A[1,1]); (* Test, ob Gleichungssystem l\"osbar. *)
+ FOR i:=2 TO dim DO Sum:=Sum + ABS(A[i,i]); END;
+ RezDim:=1.0 / VAL(LONGREAL,dim); Sum:=Sum*RezDim;
+ FOR i:=1 TO dim DO
+ IF (ABS(A[i,i]) < genau*Sum) THEN
+ Fehler:=TRUE; iFehl:=4;
+ Fehlerflag:='Gleichungssystem nicht l"osbar (LinLib.GaussSeidel)';
+ IF FehlMeld THEN ErrOut(Fehlerflag); END;
+ RETURN;
+ END;
+ END;
+
+ IF (Neu = 1) THEN (* Neuer Startvektor *)
+ FOR i:=1 TO dim 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:=1 TO dim 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:=1 TO dim DO
+ Sum:=0.0;
+ FOR j:=1 TO i-1 DO Sum:=Sum + A[i,j]*X[j]; END;
+ FOR j:=i+1 TO dim DO Sum:=Sum + A[i,j]*Xalt[j]; END; (* ?? *)
+ X[i]:=(C[i] - Sum) / A[i,i];
+ END; (* FOR i *)
+ i:=1; (* 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:=1;
+ REPEAT
+ ok := (ABS(Xalt[i]) > Unterlauf); INC(i);
+ UNTIL (i > dim) OR NOT ok;
+ IF ok THEN
+ Sum:=0.0;
+ FOR i:=1 TO dim 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:='Loesungsfolge ist divergent (LinLib.GaussSeidel)';
+ IF FehlMeld THEN ErrOut(Fehlerflag); END;
+(*
+ Der folgene Code muss hinter das Ende des LOOP-Endes geschoben werden.
+ RETURN ist durch EXIT zu ersetzen. Die Auswertung des Fehlercodes sollte
+ ebenfalls verschoben werden.
+*)
+ IF (Neu = 0) THEN FOR j:=1 TO dim DO X[j]:=Z[j]; END; END;
+ RETURN;
+ END;
+ END; (* IF *)
+ FOR i:=1 TO dim DO Xalt[i]:=X[i]; END;
+ END; (* LOOP *)
+
+ IF (Iter <= MaxIter) THEN
+ IF DiagDom THEN iFehl:=0; ELSE iFehl:=1; END;
+ ELSE
+ IF NOT DiagDom THEN
+ Fehler:=TRUE; iFehl:=3;
+ Fehlerflag:='Gleichungssystem nicht l"osbar (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;
+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 : CMATRIX; (* Koeffizientenmatrix *)
+ VAR Index : CARDVEKTOR;
+ 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 : VEKTOR;
+BEGIN
+ Fehler:=FALSE; EDet:=0;
+
+ FOR k:=1 TO dim DO
+ Sum:=0.0; Index[k]:=k;
+ FOR l:=1 TO dim 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:=1 TO dim DO
+ l:=k; z:=0.0;
+ FOR m:=k TO dim DO
+ CSum :=A[m,k];
+ FOR j:=1 TO k-1 DO
+ CSum:=CSum - A[m,j]*A[j,k];
+ 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:=1 TO dim BY 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)';
+ 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 DO
+ CSum := A[k,l];
+ FOR j:=1 TO k-1 DO
+ CSum:=CSum - A[k,j]*A[j,l];
+ END;
+ A[k,l]:= scalarMult( (1.0/z) ,CSum*CMPLX(x,-y));
+ END; (* FOR i *)
+ END; (* FOR k *)
+END CZerlege;
+
+PROCEDURE CBerechneLoesung(VAR A : CMATRIX;
+ VAR X : CVEKTOR;
+ C : CVEKTOR;
+ Index : CARDVEKTOR;
+ dim : CARDINAL);
+
+ VAR k,j : CARDINAL;
+ CSum : LONGCOMPLEX;
+BEGIN
+ CVecPerm(Index,1,dim,C); (* Bringe C in die "richtige" Ordnung *)
+
+ FOR k:=1 TO dim DO
+ C[k]:=C[k] / A[k,k]; (* Teile durch Pivot *)
+ FOR j:=k+1 TO dim DO
+ C[j]:=C[j] - A[j,k]*C[k];
+ END;
+ END; (* FOR k *)
+
+ X[dim]:=C[dim];
+ FOR k:=dim-1 TO 1 BY -1 DO
+ CSum:=C[k];
+ FOR j:=k+1 TO dim DO
+ CSum:=CSum - A[k,j]*X[j];
+ END;
+ X[k]:=CSum;
+ END;
+END CBerechneLoesung;
+
+PROCEDURE CDet(dim : CARDINAL;VAR A : CMATRIX) : LONGCOMPLEX;
+
+ VAR Index : CARDVEKTOR;
+ Det : LONGCOMPLEX;
+ EDet : INTEGER;
+ i,j : CARDINAL;
+ Dete : LONGREAL;
+ Az : POINTER TO CMATRIX;
+BEGIN
+ Fehler:=FALSE;
+ NEW(Az);
+ IF (Az # NIL) THEN
+ FOR i:=1 TO dim DO FOR j:=1 TO dim DO Az^[i,j]:=A[i,j]; 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 (Az # NIL) THEN
+ FOR i:=1 TO dim DO FOR j:=1 TO dim DO A[i,j]:=Az^[i,j]; END; END;
+ DISPOSE(Az);
+ END;
+ RETURN Det;
+END CDet;
+
+PROCEDURE CLoeseGlSy(VAR A : CMATRIX;
+ VAR X : CVEKTOR; (* L\"osungsvektor. *)
+ C : CVEKTOR;
+ VAR Det : LONGCOMPLEX; (* Determinte von A *)
+ dim : CARDINAL;
+ VAR iFehl : INTEGER);
+
+ VAR Index : CARDVEKTOR;
+ x : LONGREAL;
+ EDet : INTEGER;
+BEGIN
+ IF (dim < 1) THEN
+ Fehler:=TRUE; Fehlerflag:='dim < 1 (CLoeseGlSy)';
+ iFehl := 1;
+ Det := zero;
+ RETURN;
+ 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);
+END CLoeseGlSy;
+
+END LinLib.
+