--- a/LinLib.mod
+++ b/LinLib.mod
@@ -42,9 +42,17 @@
(* 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 *)
(*------------------------------------------------------------------------*)
(* 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. *)
@@ -58,7 +66,7 @@
(* 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 *)
+ (* hier k��nnte die Divergenz-Konstante auch Verbesserungen *)
(* bringen, bitte testen. *)
(* *)
(* Test Jacobi GaussSeidel *)
@@ -103,24 +111,33 @@
(* 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;
+ (* $Id: LinLib.mod,v 1.13 2018/08/31 14:14:51 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;
@@ -133,11 +150,10 @@
(* 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;
+ 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
@@ -205,10 +221,10 @@
DEALLOCATE(b ,N*TSIZE(LONGREAL));
END CondHager;
-PROCEDURE Zerlege(VAR A : MATRIX; (* Koeffizientenmatrix *)
- VAR IndexI : CARDVEKTOR;
- VAR IndexJ : CARDVEKTOR;
- VAR RPivot : VEKTOR;
+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 *)
@@ -222,38 +238,39 @@
(* 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;
+ 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:=1 TO dim DO
+ 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 DO
- FOR j:=k TO dim DO
+ 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 DO
+ 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 DO
+ 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:=1 TO dim DO (* Vertausch A_{j,k},A_{j,jmax]} *)
+ 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;
@@ -262,50 +279,61 @@
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;
+ 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 DO
+ FOR i:=k+1 TO dim-1 DO
Aik:=A[i,k];
- FOR j:=k+1 TO dim DO A[i,j]:=A[i,j] - Aik*A[k,j]; END;
+ 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[1]; FOR k:=2 TO dim DO Det:=Det*Pivot[k]; END;
+ 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 : MATRIX; dim : CARDINAL) : LONGREAL;
+PROCEDURE Det(VAR A : ARRAY OF ARRAY OF LONGREAL;
+ dim : CARDINAL) : LONGREAL;
VAR i,j : CARDINAL;
- Ii,Ij : CARDVEKTOR;
+ Ii,Ij : POINTER TO ARRAY [0..MAXINT-1] OF CARDINAL;
Det : LONGREAL;
- RPivot : VEKTOR;
- Az : POINTER TO MATRIX;
-BEGIN
- IF (dim = 1) THEN RETURN A[1,1]; END;
+ 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;
- NEW(Az);
+ AllocMat(Az,dim,dim);
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);
+ 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 : MATRIX;
- VAR X : VEKTOR;
- C : VEKTOR;
- VAR RPivot : VEKTOR;
- VAR IndexI : CARDVEKTOR;
- VAR IndexJ : CARDVEKTOR;
+PROCEDURE BerechneLoesung(VAR A : ARRAY OF ARRAY OF LONGREAL;
+ VAR X : ARRAY OF LONGREAL;
+ C : ARRAY OF LONGREAL;
+ VAR RPivot : ARRAY OF LONGREAL;
+ VAR IndexI : ARRAY OF CARDINAL;
+ VAR IndexJ : ARRAY OF CARDINAL;
dim : CARDINAL);
(*----------------------------------------------------------------*)
@@ -318,40 +346,42 @@
VAR k,i,j : CARDINAL;
Zw : LONGREAL;
- Z : VEKTOR;
-BEGIN
- FOR k:=1 TO dim DO
+ Z : POINTER TO ARRAY [0..MAXINT-1] OF LONGREAL;
+BEGIN
+ ALLOCATE(Z,dim*TSIZE(LONGREAL));
+ FOR k:=0 TO dim-1 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;
+ FOR i:=k+1 TO dim-1 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;
+ Z^[dim-1]:=C[dim-1]; (* Berechne die L\"osungen X[i] *)
+ FOR k:=dim-2 TO 0 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;
+ FOR j:=0 TO dim-1 DO X[IndexJ[j]]:=Z^[j]; END;
+ DEALLOCATE(Z,dim*TSIZE(LONGREAL));
END BerechneLoesung;
-PROCEDURE Gauss(VAR A : MATRIX; (* Koeffizientenmatrix *)
- VAR X : VEKTOR; (* L\"osungsvektor des LG *)
- C : VEKTOR; (* Konstantenvektor des LG *)
+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 *)
+ VAR Det : LONGREAL; (* Determinante von A *)
MaxIter : CARDINAL;
- PivStrat : CARDINAL; (* Pivotstrategie *)
- VAR iFehl : INTEGER); (* Fehlernummer *)
+ 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;
+ B,RPivot : POINTER TO ARRAY [0..MAXINT-1] OF LONGREAL;
+ IndexI,IndexJ : POINTER TO ARRAY [0..MAXINT-1] OF CARDINAL;
MaxPivot : BOOLEAN;
- Az : POINTER TO MATRIX;
+ Az : PMATRIX;
BEGIN
IF (dim < 1) THEN
Fehler:=TRUE; Fehlerflag:='dim < 1 (LinLib.Gauss)';
@@ -360,73 +390,88 @@
END;
MaxPivot := (PivStrat = 1);
Fehler:=FALSE; iFehl:=0;
- NEW(Az);
+
+ 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,j]; END; END;
+ FOR i:=1 TO dim DO
+ FOR j:=1 TO dim DO Az^[i]^[j]:=A[i-1,j-1]; END;
+ END;
ELSE
- MaxIter:=1;
Fehler:=TRUE; Fehlerflag:='kein Freispeicher vorhanden (Gauss).';
iFehl:=1;
END;
- Zerlege(A,IndexI,IndexJ,RPivot,Det,dim,MaxPivot);
+ 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) !';
+ iFehl:=3; Fehlerflag:='Matrix singulaer (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 *)
+ 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,j]:=Az^[i,j]; END; END;
- DISPOSE(Az);
- END;
+ 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 : MATRIX;
+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,Alpha,Beta,Xi : LONGREAL;
- i,j,k : CARDINAL;
- D : VEKTOR;
+ 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:=1 TO m DO
+ FOR j:=0 TO m-1 DO
Sum:=0.0;
- FOR i:=j TO n DO Sum:=Sum + A[i,j]*A[i,j]; END;
+ 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;
- D[j]:=s;
+ (* 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)';
@@ -435,22 +480,23 @@
Beta:=1.0 / Alpha;
END;
A[j,j]:=A[j,j] - s;
- FOR k:=j+1 TO m DO
+ FOR k:=j+1 TO m-1 DO
Sum:=0.0;
- FOR i:=j TO n DO Sum:=Sum + A[i,j]*A[i,k]; END;
+ FOR i:=j TO n-1 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;
+ 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 DO Sum:=Sum + A[i,j]*C[i-1]; END;
+ FOR i:=j TO n-1 DO Sum:=Sum + A[i,j]*C[i]; END;
Sum:=Beta*Sum;
- FOR i:=j TO n DO C[i-1]:=C[i-1] + A[i,j]*Sum; END;
+ FOR i:=j TO n-1 DO C[i]:=C[i] + 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;
+ 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;
@@ -531,19 +577,17 @@
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;
+ 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;
- IF (m > 256) THEN
- FatalError('Anzahl der Gleichungen gr\"o\3er 256 (HhLin) !');
- END;
+ ALLOCATE(D ,n*TSIZE(LONGREAL));
ALLOCATE(Cloc,n*TSIZE(LONGREAL));
- IF (Cloc = NIL) THEN
+ 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;
@@ -558,7 +602,7 @@
RETURN; (* !!! *)
END;
IF (A[jj] < 0.0) THEN s:=sqrt(Sum) ELSE s:=-sqrt(Sum); END;
- D[j]:=s;
+ D^[j]:=s;
Alpha:=s*A[jj] - Sum;
IF (Alpha = 0.0) THEN
Fehler:=TRUE; Fehlerflag:='Matrix singul"ar (LinLib.HhLin)';
@@ -583,9 +627,10 @@
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];
+ 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 *)
@@ -595,7 +640,7 @@
VAR iFehl : INTEGER);
VAR i,j,k,pk : CARDINAL;
- P : POINTER TO ARRAY [0..MAX(INTEGER)-1] OF CARDINAL;
+ P : POINTER TO ARRAY [0..MAXINT-1] OF CARDINAL;
RecPivot,Max : LONGREAL;
Sum,q,Zw : LONGREAL;
Melden : BOOLEAN;
@@ -694,9 +739,9 @@
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;
+ 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));
@@ -731,7 +776,7 @@
VAR i,ii,ij,k,m : CARDINAL;
p,q : LONGREAL;
- h : POINTER TO ARRAY [1..MAX(INTEGER)] OF LONGREAL;
+ h : POINTER TO ARRAY [1..MAXINT] OF LONGREAL;
BEGIN
ALLOCATE(h,N*TSIZE(LONGREAL));
iFehl:=0;
@@ -763,16 +808,17 @@
DEALLOCATE(h,N*TSIZE(LONGREAL));
END InvertSV;
-PROCEDURE Jacobi(VAR A : MATRIX; (* <== *)
- VAR X : VEKTOR; (* <==> *)
- VAR C : ARRAY OF LONGREAL; (* <== *)
+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: VEKTOR; dim : CARDINAL) : BOOLEAN;
+ 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 *)
@@ -780,14 +826,14 @@
VAR i : CARDINAL; Max : LONGREAL;
BEGIN
- Max:=ABS(V[1]);
- FOR i:=2 TO dim DO
+ 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 : MATRIX;
+ PROCEDURE TestDiagDom(VAR A : ARRAY OF ARRAY OF LONGREAL;
dim : CARDINAL;
VAR DiagDom : BOOLEAN);
@@ -798,30 +844,33 @@
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 (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 : VEKTOR;
- Iter,i,j : CARDINAL;
- Delta,Sum : LONGREAL;
- DeltaAlt1,DeltaAlt2 : LONGREAL;
- DiagDom,ok,divergent : BOOLEAN;
+ 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:=1 TO dim DO
+ 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);
@@ -829,7 +878,7 @@
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 *)
+ FOR i:=0 TO dim-1 DO X[i]:=0.0; END; (* Initialize Starting Value *)
END;
DeltaAlt1:=Ueberlauf; DeltaAlt2:=MAX(LONGREAL);
@@ -837,21 +886,25 @@
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;
+ 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:=1 TO dim DO Delta:=Delta + ABS(X[i] - Xalt[i]); END;
+ 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
@@ -875,27 +928,28 @@
END;
END Jacobi;
-PROCEDURE GaussSeidel(VAR A : MATRIX; (* Koeffizinten - Matrix *)
- VAR X : VEKTOR; (* L\"osungsvektor *)
- VAR C : VEKTOR;
+PROCEDURE GaussSeidel(VAR A : ARRAY OF ARRAY OF LONGREAL;
+ VAR X : ARRAY OF LONGREAL;
+ VAR C : ARRAY OF LONGREAL;
dim : CARDINAL;
- genau : LONGREAL; (* Abbruchkriterium *)
+ genau : LONGREAL;
MaxIter : CARDINAL;
Neu : CARDINAL;
- FehlMeld : BOOLEAN; (* Fehlermeldungen j/n 1/0 *)
+ 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 : VEKTOR;
- Sum,Prod,RezDim : LONGREAL;
- DiagDom,Konvergenz : BOOLEAN;
- ok : BOOLEAN;
-BEGIN
- FOR i:=1 TO dim DO
+ 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;
@@ -903,49 +957,54 @@
RETURN;
END;
END;
- IF (dim = 1) THEN X[1]:=C[1]/A[1,1]; Iter:=0; iFehl:=0; RETURN 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
- 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 (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[1,1]); (* Test, ob Gleichungssystem l\"osbar. *)
- FOR i:=2 TO dim DO Sum:=Sum + ABS(A[i,i]); END;
+ 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:=1 TO dim DO
+ FOR i:=0 TO dim-1 DO
IF (ABS(A[i,i]) < genau*Sum) THEN
- Fehler:=TRUE; iFehl:=4;
- Fehlerflag:='Gleichungssystem nicht l"osbar (LinLib.GaussSeidel)';
+ 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:=1 TO dim DO
+ FOR i:=0 TO dim-1 DO
IF (ABS(A[i,i]) < Unterlauf) THEN (* 14.09.15 *)
- X[i]:=MachEps; Xalt[i]:=MachEps;
+ X[i]:=MachEps; Xalt^[i]:=MachEps;
ELSE
- X[i]:=A[i,i]; Xalt[i]:=X[i]; END;
+ 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];
+ 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;
+ X[i]:=MachEps; Xalt^[i]:=MachEps;
ELSE
- Z[i]:=X[i]; Xalt[i]:=X[i];
+ Z^[i]:=X[i]; Xalt^[i]:=X[i];
END;
END;
END;
@@ -953,57 +1012,61 @@
Prod:=1.0; Iter:=0;
LOOP (* Gauss - Seidel Iterationen. *)
IF (Iter < MaxIter) THEN INC(Iter); ELSE EXIT END; (* !!! *)
- FOR i:=1 TO dim DO
+ FOR i:=0 TO dim-1 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; (* ?? *)
+ 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:=1; (* Teste auf Konvergenz *)
+ i:=0; (* Teste auf Konvergenz *)
REPEAT
- Konvergenz:=(ABS(Xalt[i] - X[i]) < genau*ABS(X[i]));
+ Konvergenz:=(ABS(Xalt^[i] - X[i]) < genau*ABS(X[i]));
INC(i);
- UNTIL NOT Konvergenz OR (i > dim);
+ UNTIL NOT Konvergenz OR (i = dim);
IF Konvergenz THEN EXIT END; (* !!! *)
IF NOT DiagDom THEN (* Teste auf Divergenz ! *)
- i:=1;
+ i:=0;
REPEAT
- ok := (ABS(Xalt[i]) > Unterlauf); INC(i);
- UNTIL (i > dim) OR NOT ok;
+ 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;
+ 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:='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;
+ Fehlerflag:=
+ 'Matrix nicht diagonal dominant & Loesungfolge divergiert';
+ Append(Fehlerflag,' (LinLib.GaussSeidel)');
+ EXIT;
END;
END; (* IF *)
- FOR i:=1 TO dim DO Xalt[i]:=X[i]; END;
+ FOR i:=0 TO dim-1 DO Xalt^[i]:=X[i]; END;
END; (* LOOP *)
- IF (Iter <= MaxIter) THEN
+ 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 l"osbar (LinLib.GaussSeidel)';
- IF (Neu = 0) THEN FOR i:=1 TO dim DO X[i]:=Z[i]; END; END;
+ 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;
@@ -1313,43 +1376,45 @@
END;
END CholSol1;
-PROCEDURE CZerlege(VAR A : CMATRIX; (* Koeffizientenmatrix *)
- VAR Index : CARDVEKTOR;
+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 : VEKTOR;
-BEGIN
+ 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:=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;
+ 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:=1 TO dim DO
+ FOR k:=0 TO dim-1 DO
l:=k; z:=0.0;
- FOR m:=k TO dim DO
+ FOR m:=k TO dim-1 DO
CSum :=A[m,k];
- FOR j:=1 TO k-1 DO
- CSum:=CSum - A[m,j]*A[j,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];
+ 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
+ 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;
@@ -1367,6 +1432,7 @@
IF (Det = zero) THEN (* 03.10.17 *)
Fehler:=TRUE; Fehlerflag:='Matrix singulaer (LinLib.CZerlege)';
+ DEALLOCATE(SumVek,dim*TSIZE(LONGREAL));
RETURN;
END;
@@ -1379,20 +1445,21 @@
w:=16.0*w; DEC(EDet,4);
Det:=scalarMult(16.0,Det);
END;
- FOR l:=k+1 TO dim DO
+ FOR l:=k+1 TO dim-1 DO
CSum := A[k,l];
- FOR j:=1 TO k-1 DO
- CSum:=CSum - A[k,j]*A[j,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 : CMATRIX;
- VAR X : CVEKTOR;
- C : CVEKTOR;
- Index : CARDVEKTOR;
+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;
@@ -1400,58 +1467,75 @@
BEGIN
CVecPerm(Index,1,dim,C); (* Bringe C in die "richtige" Ordnung *)
- FOR k:=1 TO dim DO
+ FOR k:=0 TO dim-1 DO
C[k]:=C[k] / A[k,k]; (* Teile durch Pivot *)
- FOR j:=k+1 TO dim DO
+ FOR j:=k+1 TO dim-1 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
+ X[dim-1]:=C[dim-1];
+ FOR k:=dim-2 TO 0 BY -1 DO
CSum:=C[k];
- FOR j:=k+1 TO dim DO
+ 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 : CMATRIX) : LONGCOMPLEX;
-
- VAR Index : CARDVEKTOR;
- Det : LONGCOMPLEX;
- EDet : INTEGER;
- i,j : CARDINAL;
- Dete : LONGREAL;
- Az : POINTER TO CMATRIX;
+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;
- 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;
+ 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;
+ 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;
+ 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 : CMATRIX;
- VAR X : CVEKTOR; (* L\"osungsvektor. *)
- C : CVEKTOR;
- VAR Det : LONGCOMPLEX; (* Determinte von A *)
- dim : CARDINAL;
- VAR iFehl : INTEGER);
-
- VAR Index : CARDVEKTOR;
+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
@@ -1461,8 +1545,17 @@
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);
+ CZerlege(A,Index^,Det,EDet,dim);
IF Fehler THEN
iFehl:=2;
ErrOut(Fehlerflag);
@@ -1471,8 +1564,9 @@
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);
+ CBerechneLoesung(A,X,C,Index^,dim);
+
+ DEALLOCATE(Index,dim*TSIZE(CARDINAL));
END CLoeseGlSy;
END LinLib.
-