--- a/MatLib.mod.m2pp
+++ b/MatLib.mod.m2pp
@@ -64,6 +64,8 @@
(* 25.10.17, MRi: Einfuegen von TransposeMN aus Testumgebung *)
(* 07.12.17, MRi: Einfuegen von StuerzMN,NeumaierSum & NeumaierProdSum *)
(* aus Testumgebung *)
+ (* 30.05.18, MRi: Einfuegen von TriInfNorm *)
+ (* 10.06.18, MRi: ZufallMat um den nicht-quadratischen Fall ergaenzt *)
(*------------------------------------------------------------------------*)
(* Offene Punkte *)
(* *)
@@ -89,7 +91,7 @@
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
- (* $Id: MatLib.mod.m2pp,v 1.15 2018/01/16 09:19:51 mriedl Exp mriedl $ *)
+ (* $Id: MatLib.mod.m2pp,v 1.16 2018/06/09 13:29:13 mriedl Exp mriedl $ *)
FROM SYSTEM IMPORT TSIZE,ADR;
FROM Storage IMPORT ALLOCATE,DEALLOCATE;
@@ -129,21 +131,21 @@
<* ELSE *>
EXPORT Itab,IJtab,InvIJtab,Itab0,IJtab0;
<* END *>
-
+
VAR Index : ARRAY [0..MaxDim] OF CARDINAL;
VAR Index0 : ARRAY [0..MaxDim-1] OF CARDINAL;
-
+
PROCEDURE Itab(i : CARDINAL) : CARDINAL; BEGIN RETURN Index[i]; END Itab;
-
+
PROCEDURE IJtab(i,j : CARDINAL) : CARDINAL;
-
+
BEGIN
IF (i > j) THEN RETURN Index[i] + j; ELSE RETURN Index[j] + i; END;
END IJtab;
-
+
PROCEDURE InvIJtab( ij : CARDINAL;
VAR i,j : CARDINAL);
-
+
VAR k,l,m : CARDINAL;
BEGIN
k:=1; l:=MaxDim+1;
@@ -164,13 +166,13 @@
RETURN (i*(i+1)) DIV 2;
END;
END Itab0;
-
+
PROCEDURE IJtab0(i,j : CARDINAL) : CARDINAL;
-
+
BEGIN
IF (i > j) THEN RETURN Index0[i] + j; ELSE RETURN Index0[j] + i; END;
END IJtab0;
-
+
VAR i,ii : CARDINAL;
BEGIN
ii:=0; Index[0]:=0;
@@ -259,8 +261,8 @@
InitLRArray(M[i],x,n);
END;
ELSIF (CAP(Form) = 'D') THEN
- IF (n = 1) THEN
- M[0,0]:=x;
+ IF (n = 1) THEN
+ M[0,0]:=x;
ELSE
M[0,0]:=x; FOR j:=1 TO n-1 DO M[0,j]:=0.0; END;
FOR i:=1 TO n-1 DO
@@ -299,8 +301,8 @@
InitLCArray(M[i],x,n);
END;
ELSIF (CAP(Form) = 'D') THEN
- IF (n = 1) THEN
- M[0,0]:=x;
+ IF (n = 1) THEN
+ M[0,0]:=x;
ELSE
M[0,0]:=x; FOR j:=1 TO n-1 DO M[0,j]:=zero; END;
FOR i:=1 TO n-1 DO
@@ -547,7 +549,7 @@
DeAllocMat(T,M,N);
END StuerzMN;
-PROCEDURE Transpose(VAR A : ARRAY OF ARRAY OF LONGREAL;
+PROCEDURE Transpose(VAR A : ARRAY OF ARRAY OF LONGREAL;
n : CARDINAL);
(*----------------------------------------------------------------*)
@@ -560,23 +562,23 @@
VAR i,j,nb : CARDINAL;
t : LONGREAL;
BEGIN
- nb := 0;
+ nb := 0;
WHILE nb+block-1 < n DO
FOR i:=nb TO nb+block-1 DO
FOR j:=i+1 TO nb+block-1 DO
- t:=A[i,j]; A[i,j]:=A[j,i]; A[j,i]:=t;
+ t:=A[i,j]; A[i,j]:=A[j,i]; A[j,i]:=t;
END;
END;
FOR i:=nb+block TO n-1 DO
FOR j:=nb TO nb+block-1 DO
- t:=A[i,j]; A[i,j]:=A[j,i]; A[j,i]:=t;
+ t:=A[i,j]; A[i,j]:=A[j,i]; A[j,i]:=t;
END;
END;
nb:=nb+block;
END; (* WHILE *)
FOR i:=nb TO n-1 DO (* clean up code *)
FOR j:=i+1 TO n-1 DO
- t:=A[i,j]; A[i,j]:=A[j,i]; A[j,i]:=t;
+ t:=A[i,j]; A[i,j]:=A[j,i]; A[j,i]:=t;
END;
END;
END Transpose;
@@ -596,7 +598,7 @@
VAR r,c,i,j : CARDINAL;
BEGIN
- r := re - rb;
+ r := re - rb;
c := ce - cb;
IF (r <= BlockSize) AND (c <= BlockSize) THEN
FOR i:=rb TO re-1 DO
@@ -628,7 +630,7 @@
FOR i:=0 TO dim-1 DO
FOR j:=i+1 TO dim-1 DO
Zw:=Mat[i,j]; Mat[i,j]:=Mat[j,i];
- IF Conjug THEN
+ IF Conjug THEN
Mat[i,j] := LongComplexMath.conj(Mat[i,j]);
Zw := LongComplexMath.conj(Zw);
END;
@@ -978,7 +980,7 @@
INC(ij,i);
FOR j:=i+1 TO Nm1 DO (* obere Haelfte *)
Sum:=Sum + A[ij]*X[j];
- INC(jj);
+ INC(jj);
INC(ij,jj);
END;
Y[i]:=Sum;
@@ -1037,11 +1039,11 @@
Aik0 := A[i,0];
FOR j:=0 TO N-1 DO C[i,j] := Aik0*B[0,j]; END;
IF (K2 > 1) THEN
- Aik1 := A[i,1];
+ Aik1 := A[i,1];
FOR j:=0 TO N-1 DO C[i,j]:=C[i,j] + Aik1*B[1,j]; END;
END;
IF (K2 > 2) THEN
- Aik2 := A[i,2];
+ Aik2 := A[i,2];
FOR j:=0 TO N-1 DO C[i,j]:=C[i,j] + Aik2*B[2,j]; END;
END;
END;
@@ -1051,7 +1053,7 @@
END;
END;
FOR i:=0 TO M-1 DO
- FOR k:=Kmin TO K-1 BY level DO
+ FOR k:=Kmin TO K-1 BY level DO
Aik3 := A[i,k-3]; Aik2 := A[i,k-2];
Aik1 := A[i,k-1]; Aik0 := A[i,k-0];
FOR j:=0 TO N-1 DO
@@ -1085,7 +1087,7 @@
FOR i:=0 TO Mm1 DO (* kij *)
A0i:=A[0,i]; (* kji C *)
FOR j:=0 TO Nm1 DO (*----------------------*)
- C[i,j]:=A0i*B[0,j];
+ C[i,j]:=A0i*B[0,j];
END;
END;
ELSE
@@ -1093,7 +1095,7 @@
FOR j:=0 TO Nm1 DO C[i,j]:=0.0; END;
END;
END;
- FOR k:=Kmin TO Km1 BY 2 DO
+ FOR k:=Kmin TO Km1 BY 2 DO
FOR i:=0 TO Mm1 DO
Ak1i:=A[k-1,i];
Ak0i:=A[k-0,i];
@@ -1127,7 +1129,7 @@
FOR i:=0 TO Mm1 DO (* kij 2A 2B *)
Ai0:=A[i,0]; (* kji C 2A 2B *)
FOR j:=0 TO Nm1 DO (*----------------------*)
- C[i,j]:=Ai0*B[j,0];
+ C[i,j]:=Ai0*B[j,0];
END;
END;
ELSE
@@ -1135,7 +1137,7 @@
FOR j:=0 TO Nm1 DO C[i,j]:=0.0; END;
END;
END;
- FOR k:=Kmin TO Km1 BY 2 DO
+ FOR k:=Kmin TO Km1 BY 2 DO
FOR i:=0 TO Mm1 DO
Aik1 := A[i,k-1];
Aik0 := A[i,k-0];
@@ -1243,7 +1245,7 @@
(*======================= Matrix x SUPERVEKTOR =============================*)
PROCEDURE MatSvProdNN( dim : CARDINAL;
- VAR C : ARRAY OF ARRAY OF LONGREAL;
+ VAR C : ARRAY OF ARRAY OF LONGREAL;
VAR A : ARRAY OF ARRAY OF LONGREAL;
VAR Bp : ARRAY OF LONGREAL); (* SUPERVEKTOR *)
@@ -1283,7 +1285,7 @@
END MatSvProdNN;
PROCEDURE MatSvProdTN( dim : CARDINAL;
- VAR C : ARRAY OF ARRAY OF LONGREAL;
+ VAR C : ARRAY OF ARRAY OF LONGREAL;
VAR A : ARRAY OF ARRAY OF LONGREAL;
VAR Bp : ARRAY OF LONGREAL); (* SUPERVEKTOR *)
@@ -1299,7 +1301,7 @@
FOR j:=0 TO dim DO C[i,j]:=0.0; END;
END;
kk:=0;
- FOR k:=0 TO dim-1 DO
+ FOR k:=0 TO dim-1 DO
(* store colum B[k] in Bk to enable fast access in inner j loop *)
kj:=kk;
FOR j:=0 TO k DO Bk^[j]:=Bp[kj]; INC(kj); END;
@@ -1472,8 +1474,8 @@
Sum : LONGCOMPLEX;
BEGIN
Sum:=zero;
- FOR i:=0 TO dim-1 DO
- Sum:=Sum + A[i]*B[i];
+ FOR i:=0 TO dim-1 DO
+ Sum:=Sum + A[i]*B[i];
END;
RETURN Sum;
END CSkalarProd;
@@ -1530,7 +1532,7 @@
VAR ifehl : CARDINAL);
VAR i,j,ii : CARDINAL;
-BEGIN
+BEGIN
ifehl:=0;
IF (in > adimj) OR (jn > adimj) THEN
ifehl:=1;
@@ -1554,6 +1556,26 @@
END;
END;
END InsertMat;
+
+PROCEDURE TriInfNorm( n : CARDINAL;
+ VAR D,E : ARRAY OF LONGREAL) : LONGREAL;
+
+ VAR norm,t1,t2 : LONGREAL;
+ i : CARDINAL;
+BEGIN
+ IF (n = 1) THEN
+ norm := ABS(D[0]);
+ ELSE
+ t1 := ABS(D[0]) + ABS(E[0]);
+ t2 := ABS(D[n-1]) + ABS(E[n-2]);
+ IF (t1 > t2) THEN norm:=t1; ELSE norm:=t2; END;
+ FOR i:=1 TO n-2 DO
+ t1 := ABS(D[i]) + ABS(E[i]) + ABS(E[i-1]);
+ IF (t1 > norm) THEN norm:=t1; END;
+ END;
+ END;
+ RETURN norm;
+END TriInfNorm;
PROCEDURE MatNormL1(VAR A : ARRAY OF ARRAY OF LONGREAL;
N : INTEGER;
@@ -1582,7 +1604,7 @@
CONST MinLR = 3.834E-20;
MaxLR = 1.304E+19;
-
+
VAR i : CARDINAL;
S1,S2,S3,zw : LONGREAL;
Giant,AbsX,Max1,Max3 : LONGREAL;
@@ -1623,8 +1645,8 @@
IF (S2 >= Max3) THEN
Norm := sqrt(S2*(1.0 + (Max3 / S2)*(Max3*S3)));
END;
- IF (S2 < Max3) THEN
- Norm := sqrt(Max3*((S2 / Max3) + (Max3*S3)));
+ IF (S2 < Max3) THEN
+ Norm := sqrt(Max3*((S2 / Max3) + (Max3*S3)));
END;
ELSE
Norm := Max3*sqrt(S3);
@@ -1655,8 +1677,8 @@
END;
Norm := Max;
END;
- IF (Norm > (100.0/MAX(LONGREAL))) THEN
- Norm := 1.0 / Norm;
+ IF (Norm > (100.0/MAX(LONGREAL))) THEN
+ Norm := 1.0 / Norm;
IF (dim > klein) THEN
dscal(dim,Norm,X[0],1);
ELSE
@@ -1874,50 +1896,58 @@
FOR i:=0 TO dim-1 DO V[i] := 1.0 - 2.0*Zufall(); END;
END ZufallVec;
-PROCEDURE ZufallMat(VAR M : ARRAY OF ARRAY OF LONGREAL;
- dim : CARDINAL;
+PROCEDURE ZufallMat(VAR A : ARRAY OF ARRAY OF LONGREAL;
+ m,n : CARDINAL;
sym : BOOLEAN);
VAR i,j : CARDINAL;
x : LONGREAL;
BEGIN
- IF (dim = 0) THEN RETURN; END;
- IF (dim-1 > HIGH(M[0])) THEN
- Errors.ErrOut('dim > HIGH(MAT) (ZufallMat)');
+ Errors.Fehler:=FALSE;
+ IF (m = 0) AND (n = 0) THEN RETURN; END;
+ IF (m-1 > HIGH(A)) THEN
+ Errors.ErrOut('m > HIGH(A) (ZufallMat)');
+ Errors.Fehler:=TRUE;
RETURN;
END;
+ IF (n-1 > HIGH(A[0])) THEN
+ Errors.ErrOut('n > HIGH(A[0]) (ZufallMat)');
+ Errors.Fehler:=TRUE;
+ RETURN;
+ END;
+ IF sym AND (m # n) THEN
+ Errors.ErrOut('sym = wahr und m # n (ZufallMat)');
+ Errors.Fehler:=TRUE;
+ RETURN;
+ END;
(* InitZufall; *)
- IF (dim = 1) THEN M[0,0]:=1.0 - 2.0*Zufall(); RETURN; END;
+ IF (m = 1) AND (n = 1) THEN A[0,0]:=1.0 - 2.0*Zufall(); RETURN; END;
IF sym THEN
- FOR i:=0 TO dim-1 DO
- FOR j:=i+1 TO dim-1 DO
+ FOR i:=0 TO n-1 DO
+ FOR j:=i+1 TO n-1 DO
x:=1.0 - 2.0*Zufall();
- M[i,j]:=x; M[j,i]:=x;
- END;
- M[i,i]:= 1.0 - 2.0*Zufall();
+ A[i,j]:=x; A[j,i]:=x;
+ END;
+ A[i,i]:= 1.0 - 2.0*Zufall();
END;
ELSE
- FOR i:=0 TO dim-1 DO
- FOR j:=0 TO dim-1 DO M[i,j] := 1.0 - 2.0*Zufall(); END;
+ FOR i:=0 TO m-1 DO
+ FOR j:=0 TO n-1 DO A[i,j] := 1.0 - 2.0*Zufall(); END;
END;
END;
END ZufallMat;
-PROCEDURE ZufallSv(VAR SV : SUPERVEKTOR;
+PROCEDURE ZufallSv(VAR SV : ARRAY OF LONGREAL; (* SUPERVEKTOR *)
dim : CARDINAL);
VAR i,j,ij : CARDINAL;
BEGIN
- IF (dim > MaxDim) THEN
- Errors.ErrOut('dim > MaxDim (ZufallSv)');
- dim:=MaxDim;
- END;
(* InitZufall; *)
ij:=0;
FOR i:=1 TO dim DO
- FOR j:=1 TO i DO INC(ij); SV[ij]:=1.0 - 2.0*Zufall(); END;
+ FOR j:=1 TO i DO SV[ij]:=1.0 - 2.0*Zufall(); INC(ij); END;
END;
END ZufallSv;