MatLib.mod.m2pp
2103 lines (1908 with data), 71.2 kB
IMPLEMENTATION MODULE MatLib;
(*========================================================================*)
(* HINWEIS : Bitte nur die Datei MatLib.mod.m2pp editieren *)
(*========================================================================*)
(* Es sind 2 Versionen enthalten die mit *)
(* *)
(* m2pp -D __{Parameter}__ < MatLib.mod.m2pp > MatLib.mod *)
(* *)
(* mit Parameter = {XDS|GM2} erzeugt werden koennen. *)
(* *)
(* XDS : Der Export von Objekten aus lokalen Modulen wird mit *)
(* "IMPORT" angezeigt - das scheint nicht ganz "regelkonform" *)
(* GM2 : Normales "EXPORT" wird genutzt (ELSE-Fall). *)
(* *)
(* ansonsten gibt es keine Aenderungen am Quellcode *)
(*------------------------------------------------------------------------*)
(* Bibiliotek einiger Matrix- und Vektoroperationen. *)
(* Module provides some matrix and vector operations. *)
(*------------------------------------------------------------------------*)
(* Letzte Bearbeitung *)
(* *)
(* 09.11.95, MRi: Diverse Korrekturen *)
(* 12.02.15, MRi: Prozedur InsertMat eingefuegt *)
(* 05.05.15, MRi: Prozedur InitLRArray eingefuegt *)
(* 13.10.15, MRi: Ersetzen von LFLOAT(#) durch VAL(LONGREAL,#) *)
(* 23.01.16, MRi: Ersetzen von IMPORT Random durch IMPORT RandomLib *)
(* 23.01.16, MRi: Umstellen von InitMat und ZufallMat auf "OPEN ARRAY" *)
(* 03.02.16, MRi: Umstellen von Stuerz "OPEN ARRAY" *)
(* Umstellen von MatVekProd auf "OPEN ARRY", Erweiterung *)
(* auf transponierten Fall. *)
(* 15.04.16, MRi: Kommentarfelder im Definitionsteil angepasst. Prozedur *)
(* CSkalProd reaktiviert auf Basis von Aufrufen von *)
(* CMathLib Funktionen, SkalProd auf offene Felder umgest. *)
(* In der Prozedur InitLRArray den Parameter y eingefuegt *)
(* in InitMat dim durch m,n ersetzt, groessenabfragen *)
(* korrigiert und Aufrufe von InitLRArray eingefuegt *)
(* NormVek ueberarbeitet *)
(* 17.04.16, MRi: Prozedur erneut aus LibDBlas uebernommen, ebenso *)
(* SumVek & AbsSumVek (passen hier doch besser hin). *)
(* 19.04.16, MRi: Prozedur ZufallVec eingefuegt. *)
(* 21.04.16, MRi: Prozedur CopyMat eingefuegt. *)
(* 25.04.16, MRi: Prozedur SvVekProd eingefuegt *)
(* 02.05.16, MRi: Prozedur NormMat auf offene Felder umgestellt, ruft nun *)
(* nur noch NormVek auf. In NormVek idamax eingefuegt *)
(* 18.07.16, MRi: Prozeduren SumVekUR und AbsSumVekUR eingefuegt *)
(* 27.12.16, MRi: Prozeduren MatNormL1 eingefuegt *)
(* 11.01.17, MRi: Prozeduren Itab0 und IJtab0 eingefuegt, SvToMat auf *)
(* offene Felder umgestellt und Steuerparameter sauber *)
(* definiert (war etwas "chaotisch" ;-) *)
(* 26.05.17, MRi: Ersetzen der Routinen SvSvProd, MatSvProd, SvMatProd *)
(* durch SvSvProdNN, MatSvProdNN,MatSvProdTN,SvMatProdNN *)
(* und SvMatProdNT. Prozedure InitIndex veraendert, *)
(* SvToMat auf offene Felder umgestellt und Parameter *)
(* "steuer" logischer Aufgebaut. *)
(* 26.05.17, MRi: Erfuegen der Prozeduren MatMatProd{NN|NT|TN|TT}, *)
(* Transpose, MatMatProdNNl4 (erstellt 04.01.16) *)
(* 23.07.17, MRi: Erstellen der ersten Version von TransposeMN *)
(* 24.09.17, MRi: Umbenennen von LONGCOMPLEX in LONGCMPLX sowie *)
(* CVEKTOR in CRVEKTOR und CMATRIX in CRMATRIX *)
(* 26.09.17, MRi: CStuerz auf ISO Komplex umgestelle, InitLCArray und *)
(* InitCMat hinzugefuegt *)
(* 29.09.17, MRi: Umstellung auf ISO COMPLEX *)
(* 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 *)
(* 12.09.18, MRi: CMatVekProd erweitert und an MatVekProd angepasst *)
(*------------------------------------------------------------------------*)
(* Offene Punkte *)
(* *)
(* - Nachsehen ob InitMat nicht "gelitten" hat *)
(* - Komplexe Typen auf ISO M2 umstellen *)
(* - SvToMat auf offene Felder umstellen, Code für ABS(steuer) = 1 *)
(* nutzen (ist auskommentiert) *)
(* - In SvVekProd erste Schleife eventuell ausrollen *)
(* - SvToMat pruefen und ggf. korrigieren *)
(* - CMatMatProd nochmal kurz pruefen *)
(*------------------------------------------------------------------------*)
(* Testroutinen *)
(* *)
(* - TestMatLib : (NormMat,NormVek,InitMat) *)
(* - TstMatSV : (SvVekProdNN, MatSvProd{NN|TN}, SvMatProd{NN|NT}) *)
(* - TstMaTran2 : TransposeMN *)
(* - TstRealMaMult : Matrix-Matrix Muliplitation (dgemm, Fortran dgemm) & *)
(* - TstCmplxMaMul : Matrix-Matrix Muliplitation (zgemm, CMatMatProd) *)
(* - TstMVP : Matrix-Vektor Muliplitation (dgemv, MatVekProd) *)
(* - TstMaTran2 : TransposeMN *)
(*------------------------------------------------------------------------*)
(* Implementation : Michael Riedl *)
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
(* $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;
FROM Deklera IMPORT VEKTOR,MATRIX,SUPERVEKTOR,PMATRIX,MaxDim,
CVEKTOR,CMATRIX;
FROM Errors IMPORT Fehlerflag,Fehler;
IMPORT Errors;
FROM LongMath IMPORT sqrt;
IMPORT LongComplexMath;
FROM RandomLib IMPORT Zufall;
FROM LibDBlas IMPORT idamax,dscal;
IMPORT LibDBlasM2;
FROM DynMat IMPORT AllocMat,DeAllocMat;
CONST zero = CMPLX(0.0,0.0);
(* IF NOT DEFINED(__BAxx__) THEN *) (* NEW __BAxx__+ *) (* END *)
<* IF DEFINED(__BLAS__) THEN *>
FROM LibDBlas IMPORT dcopy,ddot,daxpy;
<* END *>
PROCEDURE InitIndex(VAR Index : ARRAY OF CARDINAL;
high : CARDINAL);
VAR i,ii : CARDINAL;
BEGIN
ii:=0; FOR i:=0 TO high-1 DO Index[i]:=ii; INC(ii,i+1); END;
END InitIndex;
MODULE Tabelle; (* Kapselt Index *)
IMPORT MaxDim;
<* IF (__XDS__) THEN *>
IMPORT Itab,IJtab,InvIJtab,Itab0,IJtab0; (* EXPORT *)
<* 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;
REPEAT
m := (k + l) DIV 2;
IF (ij < Index[m]) THEN l:=m ELSE k:=m+1 END;
UNTIL (l <= k);
i:=l-1; IF (ij = Index[i]) THEN DEC(i); END;
j := ij - Index[i];
END InvIJtab;
PROCEDURE Itab0(i : CARDINAL) : CARDINAL;
BEGIN
IF (i < MaxDim) THEN
RETURN Index0[i];
ELSE
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;
FOR i:=0 TO MaxDim-1 DO Index[i+1]:=ii; INC(ii,i+1) END;
ii:=1; Index0[0]:=0;
FOR i:=1 TO MaxDim-1 DO Index0[i]:=ii; INC(ii,i+1) END;
END Tabelle;
PROCEDURE InitLRArray(VAR X : ARRAY OF LONGREAL;
y : LONGREAL;
n : CARDINAL);
VAR i,mod : CARDINAL;
unrolled : BOOLEAN;
BEGIN
unrolled := n > 32;
IF NOT unrolled THEN
FOR i:=0 TO n-1 DO X[i]:=y; END;
ELSE
IF (n > 0) THEN
mod := n MOD 4;
IF (mod # 0) THEN
FOR i:=0 TO mod-1 DO X[i]:=y; END;
END;
IF (n > 3) THEN (* Unrolled loop *)
FOR i:=mod TO n-1 BY 4 DO
X[i ]:=y; X[i+1]:=y;
X[i+2]:=y; X[i+3]:=y;
END;
END;
END;
END; (* IF unrolled *)
END InitLRArray;
PROCEDURE InitLCArray(VAR X : ARRAY OF LONGCOMPLEX;
y : LONGCOMPLEX;
n : CARDINAL);
VAR i,mod : CARDINAL;
unrolled : BOOLEAN;
BEGIN
unrolled := n > 32;
IF NOT unrolled THEN
FOR i:=0 TO n-1 DO X[i]:=y; END;
ELSE
IF (n > 0) THEN
mod := n MOD 4;
IF (mod # 0) THEN
FOR i:=0 TO mod-1 DO X[i]:=y; END;
END;
IF (n > 3) THEN (* Unrolled loop *)
FOR i:=mod TO n-1 BY 4 DO
X[i ]:=y; X[i+1]:=y;
X[i+2]:=y; X[i+3]:=y;
END;
END;
END;
END; (* IF unrolled *)
END InitLCArray;
PROCEDURE InitMat(VAR M : ARRAY OF ARRAY OF LONGREAL;
m,n : CARDINAL;
Form : CHAR;
x : LONGREAL);
VAR i,j : CARDINAL;
BEGIN
Errors.Fehler:=FALSE;
IF (n = 0) THEN RETURN; END;
IF (m-1 > HIGH(M)) OR (n-1 > HIGH(M[0])) THEN
Errors.Fehler:=TRUE;
Errors.Fehlerflag:='Dimensionierungsfehler (MatLib.InitMat)';
Errors.ErrOut(Errors.Fehlerflag);
IF (m-1 > HIGH(M)) THEN m:=HIGH(M)-1; END;
IF (n-1 > HIGH(M[0])) THEN n:=HIGH(M[0])-1; END;
ELSIF (CAP(Form) = 'D') AND (m # n) THEN
Errors.Fehler:=TRUE;
Errors.Fehlerflag:='Matrix nicht symmetrisch (MatLib.InitMat)';
Errors.ErrOut(Errors.Fehlerflag);
END;
Form:=CAP(Form); (* Form immer gro3 *)
IF (Form # 'A') AND (Form # 'D') THEN Form:='A'; x:=0.0; END;
IF (CAP(Form) = 'A') THEN
FOR i:=0 TO m-1 DO
InitLRArray(M[i],x,n);
END;
ELSIF (CAP(Form) = 'D') THEN
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
FOR j:=0 TO i-1 DO M[i,j]:=0.0; END;
M[i,i]:=x;
FOR j:=i+1 TO n-1 DO M[i,j]:=0.0; END;
END;
END;
END; (* IF form *)
END InitMat;
PROCEDURE InitCMat(VAR M : ARRAY OF ARRAY OF LONGCOMPLEX;
m,n : CARDINAL;
Form : CHAR;
x : LONGCOMPLEX);
VAR i,j : CARDINAL;
BEGIN
Errors.Fehler:=FALSE;
IF (n = 0) THEN RETURN; END;
IF (m-1 > HIGH(M)) OR (n-1 > HIGH(M[0])) THEN
Errors.Fehler:=TRUE;
Errors.Fehlerflag:='Dimensionierungsfehler (MatLib.InitMat)';
Errors.ErrOut(Errors.Fehlerflag);
IF (m-1 > HIGH(M)) THEN m:=HIGH(M)-1; END;
IF (n-1 > HIGH(M[0])) THEN n:=HIGH(M[0])-1; END;
ELSIF (CAP(Form) = 'D') AND (m # n) THEN
Errors.Fehler:=TRUE;
Errors.Fehlerflag:='Matrix nicht symmetrisch (MatLib.InitMat)';
Errors.ErrOut(Errors.Fehlerflag);
END;
Form:=CAP(Form); (* Form immer gro3 *)
IF (Form # 'A') AND (Form # 'D') THEN Form:='A'; x:=zero; END;
IF (CAP(Form) = 'A') THEN
FOR i:=0 TO m-1 DO
InitLCArray(M[i],x,n);
END;
ELSIF (CAP(Form) = 'D') THEN
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
FOR j:=0 TO i-1 DO M[i,j]:=zero; END;
M[i,i]:=x;
FOR j:=i+1 TO n-1 DO M[i,j]:=zero; END;
END;
END;
END; (* IF form *)
END InitCMat;
PROCEDURE InitSV(VAR Sv : SUPERVEKTOR;
dim : CARDINAL;
Form : CHAR;
x : LONGREAL);
VAR i,j,ij : CARDINAL;
BEGIN
Fehler:=FALSE;
IF (dim > MaxDim) THEN
Fehler:=TRUE; Fehlerflag:='dim > MaxDim (InitSV)';
Errors.ErrOut(Fehlerflag);
RETURN; (* !!! *)
END;
Form:=CAP(Form); (* Form immer gro\3 *)
IF (Form # 'A') AND (Form # 'D') THEN Form:='A'; x:=0.0; END;
IF (Form = 'A') THEN
ij:=0;
FOR i:=1 TO dim DO FOR j:=1 TO i DO INC(ij); Sv[ij]:=x; END; END;
ELSIF (Form = 'D') THEN
ij:=0;
FOR i:=1 TO dim DO
FOR j:=1 TO i-1 DO INC(ij); Sv[ij]:=0.0; END;
INC(ij); Sv[ij]:=x;
END;
END; (* IF Form *)
END InitSV;
PROCEDURE CopyMat(VAR B : ARRAY OF ARRAY OF LONGREAL; (* ==> *)
VAR A : ARRAY OF ARRAY OF LONGREAL; (* <== *)
m,n : CARDINAL);
VAR i : CARDINAL;
BEGIN
FOR i:=0 TO m-1 DO
LibDBlasM2.dcopy(n,A[i],1,B[i],1); (* copy A to B *)
END;
END CopyMat;
PROCEDURE VektorinSV(VAR SV : SUPERVEKTOR;
Vek : VEKTOR;
n : CARDINAL);
VAR i,ii : CARDINAL;
BEGIN
ii:=0; FOR i:=1 TO n DO INC(ii,i); SV[ii]:=SV[ii]+Vek[i]; END;
END VektorinSV;
PROCEDURE SvToMat(VAR SV : ARRAY OF LONGREAL; (* SUPERVEKTOR *)
VAR Mat : ARRAY OF ARRAY OF LONGREAL;
dim : CARDINAL;
steuer : INTEGER);
VAR i,j,ij : CARDINAL;
BEGIN
IF (ABS(steuer) > 2) THEN steuer:=0; END;
ij:=0;
IF (steuer = 0) THEN
FOR i:=0 TO dim-1 DO (* Belege untere Dreiecksmatrix mit SV. *)
FOR j:=0 TO i DO Mat[i,j]:=SV[ij]; INC(ij); END;
END;
FOR i:=0 TO dim-1 DO (* Belege obere Dreiecksmatrix mit SV. *)
FOR j:=i+1 TO dim-1 DO Mat[i,j]:=Mat[j,i]; END;
END;
ELSE
IF (steuer < 0) THEN
IF (steuer = -2) THEN
FOR i:=0 TO dim-1 DO (* Belege obere Dreiecksmatrix mit Null. *)
FOR j:=i+1 TO dim-1 DO Mat[i,j]:=0.0; END;
END;
END;
FOR i:=0 TO dim-1 DO (* Belege untere Dreiecksmatrix mit SV. *)
FOR j:=0 TO i DO Mat[i,j]:=SV[ij]; INC(ij); END;
END;
END;
IF (steuer > 0) THEN
IF (steuer = 2) THEN
FOR i:=1 TO dim-1 DO (* Belege untere Dreiecksmatrix mit Null. *)
FOR j:=0 TO i-1 DO Mat[i,j]:=0.0; END;
END;
END;
FOR i:=0 TO dim-1 DO (* Belege obere Dreiecksmatrix mit SV. *)
FOR j:=i TO dim-1 DO Mat[i,j]:=SV[IJtab0(i,j)]; END;
END;
END;
END; (* IF steuer = 0 *)
END SvToMat;
PROCEDURE MatToSv(VAR Mat : MATRIX;
VAR Sv : SUPERVEKTOR;
dim : CARDINAL;
steuer : INTEGER);
VAR i,j,ij : CARDINAL;
BEGIN
IF (ABS(steuer) > 1) THEN steuer:=0; END;
ij:=0;
Fehler:=FALSE;
IF (steuer < 0) THEN (* Testen !!! *)
FOR i:=1 TO dim DO
FOR j:=i TO dim DO INC(ij); Sv[ij]:=Mat[i,j]; END;
END;
ELSE
IF (steuer = 0) THEN
FOR i:=1 TO dim DO
FOR j:=1 TO i DO INC(ij); Sv[ij]:=Mat[i,j]; END;
END;
ELSE (* steuer = 1 *)
FOR i:=1 TO dim DO
FOR j:=1 TO i DO
IF (Mat[i,j] # Mat[j,i]) THEN
Fehler:=TRUE;
Fehlerflag := 'Matrix nicht symmetrisch (MatLib.MatToSv) !';
RETURN;
ELSE
INC(ij);
Sv[ij] := Mat[i,j];
END;
END; (* FOR j *)
END; (* FOR i *)
END; (* IF *)
END; (* IF *)
END MatToSv;
PROCEDURE TriMat(VAR Mat : MATRIX; (* resultierende Matrix *)
VAR HD,ND : VEKTOR; (* Haupt - Nebendiagonale *)
dim : CARDINAL;
Form : CARDINAL);
VAR i : CARDINAL;
BEGIN
IF (Form > 1) THEN Errors.FatalError(' Form > 1 (TriSv)'); END;
IF (Form = 0) THEN
InitMat(Mat,dim,dim,'a',0.0);
Mat[1,1]:=HD[1];
Mat[1,2]:=ND[1];
FOR i:=2 TO dim-1 DO
Mat[i,i]:=HD[i];
Mat[i,i-1]:=ND[i-1];
Mat[i,i+1]:=ND[i];
END;
Mat[dim,dim]:=HD[dim];
Mat[dim,dim-1]:=ND[dim-1];
ELSE
Mat[1,1]:=HD[1];
Mat[1,2]:=ND[1];
FOR i:=2 TO dim-1 DO
HD[i]:=Mat[i,i];
ND[i-1]:=Mat[i,i-1];
ND[i]:=Mat[i,i+1];
END;
HD[dim]:=Mat[dim,dim];
ND[dim-1]:=Mat[dim,dim-1];
ND[dim]:=0.0;
END;
END TriMat;
PROCEDURE TriSv(VAR SV : SUPERVEKTOR; (* resultierende Matrix *)
VAR HD,ND : VEKTOR; (* Haupt - Nebendiagonale *)
dim : CARDINAL;
Form : CARDINAL);
VAR i,ij,ii,SVdim : CARDINAL;
BEGIN
SVdim:=(dim*(dim+1)) DIV 2;
IF (Form > 1) THEN Errors.FatalError('Form > 1 (TriSv)'); END;
IF (Form = 0) THEN
FOR ij:=1 TO SVdim DO SV[ij]:=0.0; END;
SV[1]:=HD[1];
SV[2]:=ND[1];
ii:=1;
FOR i:=2 TO dim-1 DO
INC(ii,i);
SV[ii]:=HD[i];
SV[ii-1]:=ND[i-1];
END;
SV[SVdim]:=HD[dim];
SV[SVdim-1]:=ND[dim-1];
ELSE
HD[1]:=SV[1];
ND[2]:=SV[1];
ii:=1;
FOR i:=2 TO dim-1 DO
INC(ii,i);
HD[i]:=SV[ii];
ND[i-1]:=SV[ii-1];
END;
HD[dim]:=SV[SVdim];
ND[dim-1]:=SV[SVdim-1];
ND[dim]:=0.0;
END;
END TriSv;
PROCEDURE Stuerz(VAR Mat : ARRAY OF ARRAY OF LONGREAL;
dim : CARDINAL);
VAR i,j : CARDINAL;
Zw : LONGREAL;
BEGIN
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]; Mat[j,i]:=Zw;
END;
END;
END Stuerz;
PROCEDURE StuerzMN(VAR A : ARRAY OF ARRAY OF LONGREAL;
M,N : CARDINAL);
VAR i,j : CARDINAL;
T : PMATRIX;
BEGIN
IF (HIGH(A) < N) OR (HIGH(A[0]) < M) THEN
Errors.FatalError("Matrix A falsch dimensioniert (MatLib.StuerzMN)");
END;
AllocMat(T,M,N);
IF (T = NIL) THEN
Errors.FatalError("Kein Freispeicher vorhanden (MatLib.StuerzMN)");
END;
FOR i:=0 TO M-1 DO
FOR j:=0 TO N-1 DO T^[i+1]^[j+1]:=A[i,j]; END;
END;
FOR i:=N TO HIGH(A) DO
FOR j:=M TO HIGH(A[0]) DO
A[i,j]:=0.0;
END;
END;
FOR i:=0 TO M-1 DO
FOR j:=0 TO N-1 DO A[j,i]:=T^[i+1]^[j+1]; END;
END;
DeAllocMat(T,M,N);
END StuerzMN;
PROCEDURE Transpose(VAR A : ARRAY OF ARRAY OF LONGREAL;
n : CARDINAL);
(*----------------------------------------------------------------*)
(* The idea is comming form stackoverflow.com/questions/5200338/ *)
(* a-cache-efficient-matrix-transpose-program *)
(*----------------------------------------------------------------*)
CONST block = 64;
VAR i,j,nb : CARDINAL;
t : LONGREAL;
BEGIN
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;
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;
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;
END;
END;
END Transpose;
PROCEDURE TransposeMN( M,N : CARDINAL;
VAR A,T : ARRAY OF ARRAY OF LONGREAL);
(*----------------------------------------------------------------*)
(* Based on stackoverflow.com/questions/5200338 *)
(* "a-cache-efficient-matrix-transpose-program" *)
(*----------------------------------------------------------------*)
PROCEDURE CacheTranspose(rb,re : CARDINAL;
cb,ce : CARDINAL);
CONST BlockSize = 8;
VAR r,c,i,j : CARDINAL;
BEGIN
r := re - rb;
c := ce - cb;
IF (r <= BlockSize) AND (c <= BlockSize) THEN
FOR i:=rb TO re-1 DO
FOR j:=cb TO ce-1 DO
T[j,i] := A[i,j];
END;
END;
ELSIF (r >= c) THEN
CacheTranspose(rb, rb + (r DIV 2), cb, ce);
CacheTranspose(rb + (r DIV 2), re, cb, ce);
ELSE
CacheTranspose(rb, re, cb, cb + (c DIV 2));
CacheTranspose(rb, re, cb + (c DIV 2), ce);
END;
END CacheTranspose;
BEGIN
CacheTranspose(0,M+1, 0,N+1);
END TransposeMN;
PROCEDURE CStuerz(VAR Mat : ARRAY OF ARRAY OF LONGCOMPLEX;
dim : CARDINAL;
Conjug : BOOLEAN);
VAR i,j : CARDINAL;
Zw : LONGCOMPLEX;
BEGIN
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
Mat[i,j] := LongComplexMath.conj(Mat[i,j]);
Zw := LongComplexMath.conj(Zw);
END;
Mat[j,i]:=Zw;
END;
END;
END CStuerz;
PROCEDURE StuerzRev(VAR A : MATRIX;
dim : CARDINAL);
VAR i,j,p,q : CARDINAL;
x : LONGREAL;
BEGIN
q:=dim; (* q = dim-i+1 *)
FOR i:=1 TO dim DO
p:=dim; (* p = dim-j+1 *)
FOR j:=1 TO q DO x:=A[i,j]; A[i,j]:=A[p,q]; A[p,q]:=x; DEC(p); END;
DEC(q);
END;
END StuerzRev;
PROCEDURE IntSort(VAR IVek : ARRAY OF INTEGER; (* Zu Sortierender Vektor *)
VAR I : ARRAY OF CARDINAL; (* Urspr"ungliche Positionen *)
dim : CARDINAL; (* Dimension von IVek *)
ord : INTEGER); (* Steuerparameter *)
VAR i,u,o,imin,imax : CARDINAL;
Min,Max,inter : INTEGER;
Abs,Aufsteigend : BOOLEAN;
BEGIN
FOR i:=0 TO dim-1 DO I[i]:=i+1; END;
IF (ord = 0) THEN RETURN; END;
IF (dim <= 1) THEN RETURN; END;
IF (ABS(ord) > 2) THEN
Errors.ErrOut('ABS(ord) > 2 (IntSort)');
ord:=1;
END;
IF (ABS(ord) = 1) THEN Abs:=FALSE ELSE Abs:=TRUE END;
IF (ord > 0) THEN Aufsteigend:=FALSE ELSE Aufsteigend:=TRUE END;
u:=0; o:=dim-1;
REPEAT
imax:=u; Max:=IVek[u];
imin:=u; Min:=IVek[u];
IF Abs THEN Max:=ABS(Max); Min:=ABS(Min); END;
FOR i:=u+1 TO o DO
IF NOT Abs THEN
IF (IVek[i] > Max) THEN Max:=IVek[i]; imax:=i;
ELSIF (IVek[i] < Min) THEN Min:=IVek[i]; imin:=i; END;
ELSE
IF (ABS(IVek[i]) > Max) THEN Max:=ABS(IVek[i]); imax:=i;
ELSIF (ABS(IVek[i]) < Min) THEN Min:=ABS(IVek[i]); imin:=i; END;
END;
END; (* FOR i *)
IF Aufsteigend THEN
IF (imax <> o) THEN
inter:=IVek[o]; IVek[o]:=IVek[imax]; IVek[imax]:=inter;
i:=I[o]; I[o]:=I[imax]; I[imax]:=i;
END;
IF (imin = o) THEN
IF (imax = u) THEN imin:=u;
ELSE imin:=imax; END;
END;
IF (imin <> u) THEN
inter:=IVek[u]; IVek[u]:=IVek[imin]; IVek[imin]:=inter;
i:=I[u]; I[u]:=I[imin]; I[imin]:=i;
END;
ELSE (* Absteigend sortieren *)
IF (imax <> u) THEN
inter:=IVek[u]; IVek[u]:=IVek[imax]; IVek[imax]:=inter;
i:=I[u]; I[u]:=I[imax]; I[imax]:=i;
END;
IF (imin = u) THEN
IF (imax = o) THEN imin:=o;
ELSE imin:=imax; END;
END;
IF (imin <> o) THEN
inter:=IVek[o]; IVek[o]:=IVek[imin]; IVek[imin]:=inter;
i:=I[o]; I[o]:=I[imin]; I[imin]:=i;
END;
END; (* IF Absteigend *)
INC(u); DEC(o);
UNTIL (u >= o);
END IntSort;
PROCEDURE QuickSort(VAR VEK : ARRAY OF LONGREAL;
VAR II : ARRAY OF CARDINAL;
dim : CARDINAL;
ord : INTEGER);
PROCEDURE Sort(l,r : CARDINAL);
(*-----------------------------------------------------------*)
(* F\"uhrt den eigentlichen Qicksortalgorithmus durch, damit *)
(* Rekursion in QickSort m\"oglich ist. *)
(*-----------------------------------------------------------*)
VAR i,j,m : CARDINAL;
x,y : LONGREAL;
BEGIN
i:=l; j:=r;
x:=VEK[((l+r) DIV 2) - 1];
REPEAT
IF (ord = -1) THEN (* Absteigende Reihenfolge *)
WHILE (VEK[i-1] < x) DO INC(i); END;
WHILE (x < VEK[j-1]) DO DEC(j); END;
ELSE (* Aufsteigende Reihenfolge *)
WHILE (VEK[i-1] > x) DO INC(i); END;
WHILE (x > VEK[j-1]) DO DEC(j); END;
END;
IF (i <= j) THEN
y:=VEK[i-1]; VEK[i-1]:=VEK[j-1]; VEK[j-1]:=y;
m:=II[i-1]; II[i-1]:=II[j-1]; II[j-1]:=m;
INC(i); DEC(j);
END;
UNTIL (i > j);
IF (l < j) THEN Sort(l,j); END;
IF (i < r) THEN Sort(i,r); END;
END Sort;
VAR i : CARDINAL;
BEGIN
FOR i:=1 TO dim DO II[i-1]:=i; END;
IF (ord = -1) OR (ord = 1) THEN Sort(1,dim); END;
END QuickSort;
PROCEDURE Sort(VAR VEK : VEKTOR;
dim : CARDINAL;
ord : INTEGER);
VAR i,l,iminmax : CARDINAL;
minmax,Interim : LONGREAL;
BEGIN
FOR l:=1 TO dim DO
minmax:=VEK[l]; iminmax:=l;
FOR i:=l+1 TO dim DO
IF (ord = 1) THEN (* Berechne Maximum von VEK[n] *)
IF (VEK[i] > minmax) THEN
minmax:=VEK[i];
iminmax:=i;
END;
ELSE (* Berechne Minimum von VEK[n] *)
IF (VEK[i] < minmax) THEN
minmax:=VEK[i];
iminmax:=i;
END;
END; (* IF *)
END; (* FOR i *)
IF (iminmax <> l) THEN (* Vertausche VEKi mit VEKimax *)
Interim:=VEK[l];
VEK[l]:=VEK[iminmax];
VEK[iminmax]:=Interim;
END; (* IF *)
END; (* FOR l *)
END Sort;
PROCEDURE SortVek(VAR X : ARRAY OF LONGREAL; (* Zu Sortierender Vektor *)
VAR I : ARRAY OF CARDINAL; (* Urspr\"ungliche Positionen *)
dim : CARDINAL; (* Dimension von X *)
ord : INTEGER); (* Steuerparameter *)
VAR i,u,o,imin,imax : CARDINAL;
Min,Max,x : LONGREAL;
Abs,Aufsteigend : BOOLEAN;
BEGIN
FOR i:=0 TO dim-1 DO I[i]:=i+1; END;
IF (ord = 0) THEN RETURN; END;
IF (dim <= 1) THEN RETURN; END;
IF (ABS(ord) > 2) THEN
Errors.ErrOut('ABS(ord) > 2 (SortVek)');
ord:=1;
END;
IF (ABS(ord) = 1) THEN Abs:=FALSE ELSE Abs:=TRUE END;
IF (ord > 0) THEN Aufsteigend:=FALSE ELSE Aufsteigend:=TRUE END;
u:=0; o:=dim-1;
REPEAT
imax:=u; Max:=X[u];
imin:=u; Min:=X[u];
IF Abs THEN Max:=ABS(Max); Min:=ABS(Min); END;
FOR i:=u+1 TO o DO
IF NOT Abs THEN
IF (X[i] > Max) THEN Max:=X[i]; imax:=i;
ELSIF (X[i] < Min) THEN Min:=X[i]; imin:=i; END;
ELSE
IF (ABS(X[i]) > Max) THEN Max:=ABS(X[i]); imax:=i;
ELSIF (ABS(X[i]) < Min) THEN Min:=ABS(X[i]); imin:=i; END;
END;
END; (* FOR i *)
IF Aufsteigend THEN
IF (imax <> o) THEN
x:=X[o]; X[o]:=X[imax]; X[imax]:=x;
i:=I[o]; I[o]:=I[imax]; I[imax]:=i;
END;
IF (imin = o) THEN
IF (imax = u) THEN imin:=u;
ELSE imin:=imax; END;
END;
IF (imin <> u) THEN
x:=X[u]; X[u]:=X[imin]; X[imin]:=x;
i:=I[u]; I[u]:=I[imin]; I[imin]:=i;
END;
ELSE (* Absteigend sortieren *)
IF (imax <> u) THEN
x:=X[u]; X[u]:=X[imax]; X[imax]:=x;
i:=I[u]; I[u]:=I[imax]; I[imax]:=i;
END;
IF (imin = u) THEN
IF (imax = o) THEN imin:=o;
ELSE imin:=imax; END;
END;
IF (imin <> o) THEN
x:=X[o]; X[o]:=X[imin]; X[imin]:=x;
i:=I[o]; I[o]:=I[imin]; I[imin]:=i;
END;
END; (* IF Absteigend *)
INC(u); DEC(o);
UNTIL (u >= o);
END SortVek;
PROCEDURE CSortVek(VAR X : ARRAY OF LONGCOMPLEX; (* Zu Sortierend *)
VAR I : ARRAY OF CARDINAL; (* Urspr\"ungl. Positionen *)
dim : CARDINAL; (* Dimension von X *)
ord : INTEGER); (* Steuerparameter *)
VAR i,u,o,imin,imax : CARDINAL;
Min,Max : LONGREAL;
x : LONGCOMPLEX;
Aufsteigend : BOOLEAN;
AbsX : ARRAY [0..MaxDim-1] OF LONGREAL;
BEGIN
FOR i:=0 TO dim-1 DO I[i]:=i; END;
IF (ord = 0) THEN RETURN; END;
FOR i:=0 TO dim-1 DO AbsX[i]:=LongComplexMath.abs(X[i]); END;
IF (dim <= 1) THEN RETURN; END;
IF (ABS(ord) > 1) THEN
Errors.ErrOut('ABS(ord) > 1 (CSortVek)');
ord:=1;
END;
IF (ord > 0) THEN Aufsteigend:=FALSE; ELSE Aufsteigend:=TRUE; END;
u:=0; o:=dim-1;
REPEAT
imax:=u; Max:=AbsX[u];
imin:=u; Min:=AbsX[u];
FOR i:=u+1 TO o DO
IF (AbsX[i] > Max) THEN
Max:=AbsX[i]; imax:=i;
ELSIF (AbsX[i] < Min) THEN
Min:=AbsX[i]; imin:=i;
END;
END; (* FOR i *)
IF Aufsteigend THEN
IF (imax <> o) THEN
x:=X[o]; X[o]:=X[imax]; X[imax]:=x; AbsX[imax]:=AbsX[o];
i:=I[o]; I[o]:=I[imax]; I[imax]:=i;
END;
IF (imin = o) THEN
IF (imax = u) THEN imin:=u; ELSE imin:=imax; END;
END;
IF (imin <> u) THEN
x:=X[u]; X[u]:=X[imin]; X[imin]:=x; AbsX[imin]:=AbsX[u];
i:=I[u]; I[u]:=I[imin]; I[imin]:=i;
END;
ELSE (* Absteigend sortieren *)
IF (imax <> u) THEN
x:=X[u]; X[u]:=X[imax]; X[imax]:=x; AbsX[imax]:=AbsX[u];
i:=I[u]; I[u]:=I[imax]; I[imax]:=i;
END;
IF (imin = u) THEN
IF (imax = o) THEN imin:=o; ELSE imin:=imax; END;
END;
IF (imin <> o) THEN
x:=X[o]; X[o]:=X[imin]; X[imin]:=x; AbsX[imin]:=AbsX[o];
i:=I[o]; I[o]:=I[imin]; I[imin]:=i;
END;
END; (* IF Absteigend *)
INC(u); DEC(o);
UNTIL (u >= o);
END CSortVek;
PROCEDURE MatVekProd(VAR Y : ARRAY OF LONGREAL;
VAR A : ARRAY OF ARRAY OF LONGREAL;
X : ARRAY OF LONGREAL;
M,N : CARDINAL;
trans : BOOLEAN);
VAR i,j : CARDINAL;
Xi,Yi : LONGREAL;
BEGIN
IF (NOT trans) THEN (* y = A*x *)
FOR i:=0 TO M-1 DO
Yi := 0.0;
FOR j:=0 TO N-1 DO Yi:=Yi + A[i,j]*X[j]; END;
Y[i]:=Yi;
END;
ELSE (* y = A^{tr}*x *)
FOR j:=0 TO N-1 DO Y[j]:=0.0; END;
FOR i:=0 TO M-1 DO
IF (X[i] # 0.0) THEN
Xi := X[i];
FOR j:=0 TO N-1 DO Y[j]:=Y[j] + A[i,j]*Xi; END;
END;
END;
END;
END MatVekProd;
PROCEDURE CMatVekProd(VAR Y : ARRAY OF LONGCOMPLEX;
VAR A : ARRAY OF ARRAY OF LONGCOMPLEX;
X : ARRAY OF LONGCOMPLEX;
M,N : CARDINAL;
trans : CHAR);
VAR i,j : CARDINAL;
Xi,Yi : LONGCOMPLEX;
BEGIN
IF (CAP(trans) = "N") THEN (* y = A*x *)
FOR i:=0 TO M-1 DO
Yi := zero;
FOR j:=0 TO N-1 DO Yi:=Yi + A[i,j]*X[j]; END;
Y[i]:=Yi;
END;
ELSIF (CAP(trans) = "T") THEN (* y = A^{tr}*x *)
FOR j:=0 TO N-1 DO Y[j]:=zero; END;
FOR i:=0 TO M-1 DO
IF (X[i] # zero) THEN
Xi := X[i];
FOR j:=0 TO N-1 DO Y[j]:=Y[j] + A[i,j]*Xi; END;
END;
END;
ELSIF (CAP(trans) = "C") THEN (* y = conj(A^{tr})*x *)
FOR j:=0 TO N-1 DO Y[j]:=zero; END;
FOR i:=0 TO M-1 DO
IF (X[i] # zero) THEN
Xi := X[i];
FOR j:=0 TO N-1 DO
Y[j]:=Y[j] + LongComplexMath.conj(A[i,j])*Xi;
END;
END;
END;
ELSE
Errors.FatalError("trans # CAP({N|T|C}) (CMatVekProd)");
END;
END CMatVekProd;
PROCEDURE SvVekProd(VAR Y : ARRAY OF LONGREAL; (* ==> *)
VAR A : ARRAY OF LONGREAL; (* <== SUPERVEKTOR *)
VAR X : ARRAY OF LONGREAL; (* <== *)
dim : CARDINAL);
VAR Nm1,i,j,ii,jj,ij : CARDINAL;
Sum : LONGREAL;
BEGIN
Nm1 := dim-1;
ii:=0;
FOR i:=0 TO Nm1 DO
Sum:=0.0;
ij:=ii;
(* OFFEN : Naechste Schleife eventuell ausrollen *)
FOR j:=0 TO i DO (* untere Haelfte + Diagonale *)
Sum:=Sum + A[ij]*X[j]; INC(ij);
END;
jj:=i+1;
INC(ij,i);
FOR j:=i+1 TO Nm1 DO (* obere Haelfte *)
Sum:=Sum + A[ij]*X[j];
INC(jj);
INC(ij,jj);
END;
Y[i]:=Sum;
INC(ii,i+1);
END;
END SvVekProd;
PROCEDURE MatMatProd(VAR C : MATRIX; (* Resultierende MATRIX *)
VAR A : MATRIX;
VAR B : MATRIX;
dim : CARDINAL);
VAR i,j,k : CARDINAL;
x : LONGREAL;
Vek : VEKTOR;
B0 : POINTER TO MATRIX;
BEGIN
NEW(B0);
IF (B0 = NIL) THEN
Errors.FatalError("kein Freispeicher vorhanden(MatLib.MatMatProd) !");
END;
FOR i:=1 TO dim DO FOR j:=1 TO dim DO B0^[i,j]:=B[i,j]; END; END;
FOR i:=1 TO dim DO
FOR j:=1 TO dim DO Vek[j]:=0.0; END;
FOR j:=1 TO dim DO
x:=A[i,j];
IF (x # 0.0) THEN
FOR k:=1 TO dim DO Vek[k]:=Vek[k] + x*B0^[j,k]; END;
END;
END;
FOR j:=1 TO dim DO C[i,j]:=Vek[j]; END;
END;
DISPOSE(B0);
END MatMatProd;
PROCEDURE MatMatProdNN( M,N,K : CARDINAL;
VAR C : ARRAY OF ARRAY OF LONGREAL;
VAR A,B : ARRAY OF ARRAY OF LONGREAL);
(*----------------------------------------------------------------*)
(* A[0..M-1][0..K-1] x B[0..K-1][0..N-1] = C[0..M-1][0..N-1] *)
(* Schleifen werden nur vierfach ausgerollt da nur eine Zeilen *)
(* der Matrix B in den Zwischenspeicher passen muss *)
(*----------------------------------------------------------------*)
CONST level = 4; (* Schleifen werden bis "level" ausgerollt *)
VAR Aik0,Aik1,Aik2,Aik3 : LONGREAL;
i,j,k,K2,Kmin : CARDINAL;
BEGIN
K2 := (K MOD level);
Kmin := K2 + level - 1;
IF (K2 # 0) THEN
FOR i:=0 TO M-1 DO
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];
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];
FOR j:=0 TO N-1 DO C[i,j]:=C[i,j] + Aik2*B[2,j]; END;
END;
END;
ELSE (* initialize C with zero *)
FOR i:=0 TO M-1 DO
FOR j:=0 TO N-1 DO C[i,j]:=0.0; END;
END;
END;
FOR i:=0 TO M-1 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
C[i,j]:=C[i,j] + Aik3*B[k-3,j] + Aik2*B[k-2,j]
+ Aik1*B[k-1,j] + Aik0*B[k-0,j];
END;
END;
END;
END MatMatProdNN;
PROCEDURE MatMatProdTN( M,N,K : CARDINAL;
VAR C : ARRAY OF ARRAY OF LONGREAL;
VAR A,B : ARRAY OF ARRAY OF LONGREAL);
(*----------------------------------------------------------------*)
(* A[0..K-1][0..M-1] x B[0..K-1][0..N-1] = C[0..M-1][0..N-1] *)
(* Schleifen werden nur doppelt ausgerollt da zwei Zeilen der *)
(* Matrix B in den Zwischenspeicher passen mussen *)
(*----------------------------------------------------------------*)
VAR A0i,Ak1i,Ak0i : LONGREAL;
i,j,k,K2,Kmin : CARDINAL;
Mm1,Nm1,Km1 : CARDINAL;
BEGIN
Mm1:=M-1; (*----------------------*)
Nm1:=N-1; (* Schleifenanalyse *)
Km1:=K-1; (* ijk 2A 2B *)
K2 := K MOD 2; (* ikj 2A *)
Kmin := K2 + 1; (* jik C^ 2A 2B *)
IF (K2 # 0) THEN (* jki C 2B *)
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];
END;
END;
ELSE
FOR i:=0 TO Mm1 DO
FOR j:=0 TO Nm1 DO C[i,j]:=0.0; END;
END;
END;
FOR k:=Kmin TO Km1 BY 2 DO
FOR i:=0 TO Mm1 DO
Ak1i:=A[k-1,i];
Ak0i:=A[k-0,i];
FOR j:=0 TO Nm1 DO
C[i,j]:=C[i,j] + Ak1i*B[k-1,j] + Ak0i*B[k+0,j];
END;
END;
END;
END MatMatProdTN;
PROCEDURE MatMatProdNT( M,N,K : CARDINAL;
VAR C : ARRAY OF ARRAY OF LONGREAL;
VAR A,B : ARRAY OF ARRAY OF LONGREAL);
(*----------------------------------------------------------------*)
(* A[0..M-1][0..K-1] x B[0..N-1][0..K-1] = C[0..M-1][0..N-1] *)
(* Schleifen werden nur doppelt ausgerollt da zwei Zeilen der *)
(* Matrix B in den Zwischenspeicher passen mussen *)
(*----------------------------------------------------------------*)
VAR Ai0,Aik1,Aik0 : LONGREAL;
i,j,k,K2,Kmin : CARDINAL;
Mm1,Nm1,Km1 : CARDINAL;
BEGIN
Mm1:=M-1; (*----------------------*)
Nm1:=N-1; (* Schleifenanalyse *)
Km1:=K-1; (* ijk *)
K2 := K MOD 2; (* ikj 2B *)
Kmin := K2 + 1; (* jik C^ *)
IF (K2 # 0) THEN (* jki C 2A *)
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];
END;
END;
ELSE
FOR i:=0 TO Mm1 DO
FOR j:=0 TO Nm1 DO C[i,j]:=0.0; END;
END;
END;
FOR k:=Kmin TO Km1 BY 2 DO
FOR i:=0 TO Mm1 DO
Aik1 := A[i,k-1];
Aik0 := A[i,k-0];
FOR j:=0 TO Nm1 DO
C[i,j]:=C[i,j] + Aik1*B[j,k-1] + Aik0*B[j,k+0];
END;
END;
END;
END MatMatProdNT;
PROCEDURE MatMatProdTT( M,N,K : CARDINAL;
VAR C : ARRAY OF ARRAY OF LONGREAL;
VAR A,B : ARRAY OF ARRAY OF LONGREAL);
(*----------------------------------------------------------------*)
(* A[0..K-1][0..M-1] x B[0..N-1][0..K-1] = C[0..M-1][0..N-1] *)
(* Schleifen werden nur doppelt ausgerollt da zwei Zeilen der *)
(* Matrix A in den Zwischenspeicher passen mussen *)
(*----------------------------------------------------------------*)
VAR Bj0,Bjk1,Bjk0 : LONGREAL;
i,j,k,K2,Kmin : CARDINAL;
Mm1,Nm1,Km1 : CARDINAL;
BEGIN
IF (N = K) THEN (* Erheblich schneller - much faster !!! *)
IF (N < 32) THEN Stuerz(B,N); ELSE Transpose(B,N); END;
MatMatProdTN(M,N,K,C,A,B);
IF (N < 32) THEN Stuerz(B,N); ELSE Transpose(B,N); END;
RETURN;
END;
Mm1:=M-1; (*----------------------*)
Nm1:=N-1; (* Schleifenanalyse *)
Km1:=K-1; (* ijk 2A *)
K2 := K MOD 2; (* ikj 2A + 2B *)
Kmin := K2 + 1; (* jik 1C^ + 2A *)
IF (K2 # 0) THEN (* jki 1C + 2A *)
FOR j:=0 TO Nm1 DO (* kij 2A^ + 2B *)
Bj0:=B[j,0]; (* kji 1C + 2A *)
FOR i:=0 TO Mm1 DO (*----------------------*)
C[i,j]:=A[0,i]*Bj0;
END;
END;
ELSE
FOR i:=0 TO Mm1 DO
FOR j:=0 TO Nm1 DO C[i,j]:=0.0; END;
END;
END;
FOR k:=Kmin TO Km1 BY 2 DO
FOR j:=0 TO Nm1 DO
Bjk1 := B[j,k-1];
Bjk0 := B[j,k+0];
FOR i:=0 TO Mm1 DO
C[i,j]:=C[i,j] + A[k-1,i]*Bjk1 + A[k+0,i]*Bjk0;
END;
END;
END;
END MatMatProdTT;
PROCEDURE CMatMatProd(VAR C : CMATRIX; (* Resultierende MATRIX *)
VAR A : CMATRIX;
VAR B : CMATRIX;
dim : CARDINAL);
VAR i,j,k : CARDINAL;
x : LONGCOMPLEX;
Vek : CVEKTOR;
B0 : POINTER TO CMATRIX;
BEGIN
NEW(B0);
IF (B0 = NIL) THEN
Errors.FatalError("kein Freispeicher vorhanden(MatLib.CMatMatProd) !");
END;
FOR i:=1 TO dim DO FOR j:=1 TO dim DO B0^[i,j]:=B[i,j]; END; END;
FOR i:=1 TO dim DO
FOR j:=1 TO dim DO Vek[j]:=zero; END;
FOR j:=1 TO dim DO
x:=A[i,j];
IF (x # zero) THEN (* Vek[k]:=Vek[k] + x*B[j,k]; *)
IF (RE(x) = 0.0) THEN
FOR k:=1 TO dim DO
Vek[k]:=Vek[k] + CMPLX(-IM(x)*IM(B0^[j,k]),IM(x)*RE(B0^[j,k]));
(* Vek[k].re:=Vek[k].re - x.im*B0^[j,k].im; *)
(* Vek[k].im:=Vek[k].im + x.im*B0^[j,k].re; *)
END;
ELSIF (IM(x) = 0.0) THEN
FOR k:=1 TO dim DO
Vek[k]:=Vek[k] + CMPLX(RE(x)*RE(B0^[j,k]),RE(x)*IM(B0^[j,k]));
(* Vek[k].re:=Vek[k].re + x.re*B0^[j,k].re; *)
(* Vek[k].im:=Vek[k].im + x.re*B0^[j,k].im; *)
END;
ELSE
FOR k:=1 TO dim DO
Vek[k]:=Vek[k] + x*B0^[j,k];
END;
END;
END; (* IF *)
END; (* FOR j *)
FOR j:=1 TO dim DO C[i,j]:=Vek[j]; END;
END; (* FOR i *)
DISPOSE(B0);
END CMatMatProd;
(*======================= Matrix x SUPERVEKTOR =============================*)
PROCEDURE MatSvProdNN( dim : CARDINAL;
VAR C : ARRAY OF ARRAY OF LONGREAL;
VAR A : ARRAY OF ARRAY OF LONGREAL;
VAR Bp : ARRAY OF LONGREAL); (* SUPERVEKTOR *)
VAR i,j,k,jj,kj : CARDINAL;
Sum : LONGREAL;
Ci : POINTER TO ARRAY [0..MAX(INTEGER)-1] OF LONGREAL;
BEGIN
ALLOCATE(Ci,dim*TSIZE(LONGREAL)); (* Speicher fuer Zwischenergebniss *)
IF (Ci = NIL) THEN
Errors.FatalError("Kein Freispeicher vorhanden (MatSvProd)");
END;
FOR i:=0 TO dim-1 DO
jj:=0;
FOR j:=0 TO dim-1 DO
Sum:=0.0; kj:=jj;
(* optimal *)
IF (j > 0) THEN
<* IF DEFINED(__BLAS__) THEN *>
Sum:=ddot(j,A[i][0],1,Bp[jj],1);
INC(kj,j);
<* ELSE *>
FOR k:=0 TO j-1 DO Sum:=Sum + A[i,k]*Bp[kj]; INC(kj); END;
<* END *>
END;
(* nicht optimal ... *)
FOR k:=j TO dim-1 DO Sum:=Sum + A[i,k]*Bp[kj]; INC(kj,k+1); END;
Ci^[j]:=Sum;
INC(jj,j+1);
END;
<* IF DEFINED(__BLAS__) THEN *>
dcopy(VAL(INTEGER,dim),Ci^[0],1,C[i][0],1);
<* ELSE *>
FOR j:=0 TO dim-1 DO C[i,j]:=Ci^[j]; END;
<* END *>
END;
DEALLOCATE(Ci,dim*TSIZE(LONGREAL));
END MatSvProdNN;
PROCEDURE MatSvProdTN( dim : CARDINAL;
VAR C : ARRAY OF ARRAY OF LONGREAL;
VAR A : ARRAY OF ARRAY OF LONGREAL;
VAR Bp : ARRAY OF LONGREAL); (* SUPERVEKTOR *)
VAR i,j,k,kk,kj : CARDINAL;
Aki : LONGREAL;
Bk : POINTER TO ARRAY [0..MaxDim-1] OF LONGREAL;
BEGIN
ALLOCATE(Bk,dim*TSIZE(LONGREAL));
IF (Bk = NIL) THEN
Errors.FatalError("kein Freispeicher vorhanden (MatLib.SvSvProd) !");
END;
FOR i:=0 TO dim-1 DO
FOR j:=0 TO dim DO C[i,j]:=0.0; END;
END;
kk:=0;
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;
INC(kj,k);
FOR j:=k+1 TO dim-1 DO Bk^[j]:=Bp[kj]; INC(kj,j+1); END;
(* ideale Struktur *)
FOR i:=0 TO dim-1 DO
Aki:=A[k,i];
<* IF DEFINED(__BLAS__) THEN *>
daxpy(dim,Aki,Bk^[0],1,C[i][0],1);
<* ELSE *>
FOR j:=0 TO dim-1 DO C[i,j]:=C[i,j] + Aki*Bk^[j]; END;
<* END *>
END;
INC(kk,k+1);
END;
DEALLOCATE(Bk,dim*TSIZE(LONGREAL));
END MatSvProdTN;
(*======================= SUPERVEKTOR x Matrix =============================*)
PROCEDURE SvMatProdNN( dim : CARDINAL;
VAR C : ARRAY OF ARRAY OF LONGREAL;
VAR Ap : ARRAY OF LONGREAL; (* SUPERVEKTOR *)
VAR B : ARRAY OF ARRAY OF LONGREAL);
VAR i,j,k,ii,ik : CARDINAL;
Aik : LONGREAL;
BEGIN
ii:=0;
FOR i:=0 TO dim-1 DO (* Loops reorderd i,k,j *)
FOR j:=0 TO dim-1 DO C[i,j]:=0.0; END;
ik:=ii;
(* ideale Struktur *)
IF (i > 0) THEN
FOR k:=0 TO i-1 DO
Aik := Ap[ik];
<* IF DEFINED(__BLAS__) THEN *>
daxpy(dim,Aik,B[k][0],1,C[i][0],1);
<* ELSE *>
FOR j:=0 TO dim-1 DO C[i,j]:=C[i,j] + Aik*B[k,j]; END;
<* END *>
INC(ik);
END;
END;
FOR k:=i TO dim-1 DO
Aik := Ap[ik];
<* IF DEFINED(__BLAS__) THEN *>
daxpy(dim,Aik,B[k][0],1,C[i][0],1);
<* ELSE *>
FOR j:=0 TO dim-1 DO C[i,j]:=C[i,j] + Aik*B[k,j]; END;
<* END *>
INC(ik,k+1);
END;
INC(ii,i+1);
END;
END SvMatProdNN;
PROCEDURE SvMatProdNT( dim : CARDINAL;
VAR C : ARRAY OF ARRAY OF LONGREAL;
VAR Ap : ARRAY OF LONGREAL; (* SUPERVEKTOR *)
VAR B : ARRAY OF ARRAY OF LONGREAL);
VAR i,j,ii,ij : CARDINAL;
<* IF NOT DEFINED(__BLAS__) THEN *>
k : CARDINAL;
<* END *>
Sum : LONGREAL;
Ai,Ci : POINTER TO ARRAY [0..MAX(INTEGER)-1] OF LONGREAL;
BEGIN
ALLOCATE(Ci,dim*TSIZE(LONGREAL));
ALLOCATE(Ai,dim*TSIZE(LONGREAL));
IF (Ci = NIL) OR (Ai = NIL) THEN
Errors.FatalError("Kein Freispeicher vorhanden (SvMatProdNT)");
END;
ii:=0;
FOR i:=0 TO dim-1 DO
(* store colum A[i] in Ai to enable fast access in inner k loop *)
ij:=ii;
FOR j:=0 TO i DO Ai^[j]:=Ap[ij]; INC(ij); END;
INC(ij,i);
FOR j:=i+1 TO dim-1 DO Ai^[j]:=Ap[ij]; INC(ij,j+1); END;
(* ideale Struktur *)
FOR j:=0 TO dim-1 DO
<* IF DEFINED(__BLAS__) THEN *>
Sum:=ddot(dim,Ai^[0],1,B[j][0],1);
<* ELSE *>
Sum:=0.0;
FOR k:=0 TO dim-1 DO Sum:=Sum + Ai^[k]*B[j,k]; END;
<* END *>
Ci^[j]:=Sum;
END; (* FOR j *)
INC(ii,i+1);
<* IF DEFINED(__BLAS__) THEN *>
dcopy(VAL(INTEGER,dim),Ci^[0],1,C[i][0],1);
<* ELSE *>
FOR j:=0 TO dim-1 DO C[i,j] := Ci^[j]; END;
<* END *>
END; (* FOR i *)
DEALLOCATE(Ai,dim*TSIZE(LONGREAL));
DEALLOCATE(Ci,dim*TSIZE(LONGREAL));
END SvMatProdNT;
(*======================= SUPERVEKTOR x SUPERVEKTOR ========================*)
PROCEDURE SvSvProdNN( dim : CARDINAL;
VAR C : ARRAY OF ARRAY OF LONGREAL;
VAR Ap,Bp : ARRAY OF LONGREAL); (* SUPERVEKTOR *)
VAR i,j,k : CARDINAL;
ii,jj,ij,jk : CARDINAL;
Sum : LONGREAL;
Ai : POINTER TO ARRAY [0..MaxDim-1] OF LONGREAL;
BEGIN
ALLOCATE(Ai,dim*TSIZE(LONGREAL));
IF (Ai = NIL) THEN
Errors.FatalError("kein Freispeicher vorhanden (MatLib.SvSvProd) !");
END;
ii:=0;
FOR i:=0 TO dim-1 DO
(* store colum A[i] in Ai to enable fast access in inner k loops *)
ij:=ii;
FOR j:=0 TO i DO Ai^[j]:=Ap[ij]; INC(ij); END;
INC(ij,i);
FOR j:=i+1 TO dim-1 DO Ai^[j]:=Ap[ij]; INC(ij,j+1); END;
jj:=0;
FOR j:=0 TO dim-1 DO
jk:=jj;
(* data for A & B the following loop should fit in cache ... *)
<* IF DEFINED(__BLAS__) THEN *>
Sum:=ddot(j+1,Ai^[0],1,Bp[jj],1);
INC(jk,j+1);
<* ELSE *>
Sum:=0.0; (* Inline code ... *)
FOR k:=0 TO j DO Sum:=Sum + Ai^[k]*Bp[jk]; INC(jk); END;
<* END *>
INC(jk,j);
(* Only data for A in the following loop fit in cache ... *)
FOR k:=j+1 TO dim-1 DO
Sum:=Sum + Ai^[k]*Bp[jk];
INC(jk,k+1); (* we cannot avoide this "jumping ahead" ... *)
END;
C[i,j]:=Sum;
INC(jj,j+1);
END; (* FOR j *)
INC(ii,i+1);
END; (* FOR i *)
DEALLOCATE(Ai,dim*TSIZE(LONGREAL));
END SvSvProdNN;
(*================= Ende Matrixmultiplikation SUPERVEKTOR ==================*)
PROCEDURE SkalarProd(VAR A,B : ARRAY OF LONGREAL;
dim : CARDINAL) : LONGREAL;
VAR i : CARDINAL;
Sum : LONGREAL;
BEGIN
Sum:=0.0; FOR i:=0 TO dim-1 DO Sum:=Sum + A[i]*B[i]; END;
RETURN Sum;
END SkalarProd;
PROCEDURE CSkalarProd(VAR A,B : ARRAY OF LONGCOMPLEX;
dim : CARDINAL) : LONGCOMPLEX;
VAR i : CARDINAL;
Sum : LONGCOMPLEX;
BEGIN
Sum:=zero;
FOR i:=0 TO dim-1 DO
Sum:=Sum + A[i]*B[i];
END;
RETURN Sum;
END CSkalarProd;
PROCEDURE MatSkalProd(VAR MAT : MATRIX; (* Resultierende MATRIX *)
x : LONGREAL;
dim : CARDINAL);
VAR i,j : CARDINAL;
BEGIN
FOR i:=1 TO dim DO
FOR j:=1 TO dim DO MAT[i,j]:=x*MAT[i,j]; END;
END;
END MatSkalProd;
PROCEDURE KreuzProd(VAR Mat : MATRIX; (* Resultierende Matrix *)
VAR Vec1 : VEKTOR;
VAR Vec2 : VEKTOR;
dim : CARDINAL); (* L\"ange der Vektoren *)
VAR i,j : CARDINAL;
BEGIN
FOR i:=1 TO dim DO
FOR j:=1 TO dim DO Mat[i,j] := Vec1[i] * Vec2[j]; END;
END;
END KreuzProd;
PROCEDURE SpaltVek(VAR SpVek : VEKTOR;
VAR MAT : MATRIX;
dim : CARDINAL;
j : CARDINAL); (* j.-Spalte wird SpVek *)
VAR i : CARDINAL;
BEGIN
FOR i:=1 TO dim DO SpVek[i]:=MAT[i,j]; END;
END SpaltVek;
PROCEDURE VekInMat(VAR VEK : VEKTOR;
VAR MAT : MATRIX;
j : CARDINAL;
dim : CARDINAL);
VAR i : CARDINAL;
BEGIN
FOR i:=1 TO dim DO MAT[i,j]:=VEK[i]; END;
END VekInMat;
PROCEDURE InsertMat(VAR A : MATRIX;
adimi : CARDINAL; (* first dimension of A *)
adimj : CARDINAL; (* second dimension of A *)
VAR B : MATRIX; (* Matrix to be inserted *)
in,jn : CARDINAL; (* 1. & 2. dimension ob B *)
is,js : CARDINAL; (* Startpoints where to insert *)
VAR ifehl : CARDINAL);
VAR i,j,ii : CARDINAL;
BEGIN
ifehl:=0;
IF (in > adimj) OR (jn > adimj) THEN
ifehl:=1;
END;
IF (is + in - 1 > adimi) THEN
ifehl:=ifehl+10;
END;
IF (js + jn - 1 > adimj) THEN
ifehl:=ifehl+100;
END;
IF (ifehl # 0) THEN
Errors.ErrOut(
"InsertMat : Matrix B kann nicht in Matrix A eingesetzt werden");
RETURN;
END;
FOR i:=1 TO in DO
ii:=is+i-1;
FOR j:=1 TO jn DO
A[ii,js+j-1]:=B[i,j];
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;
VAR Norm : LONGREAL);
VAR i,j : INTEGER;
Sum : LONGREAL;
BEGIN
Norm:=0.0;
FOR i:=0 TO N-1 DO
Sum:=0.0;
FOR j:=0 TO N-1 DO Sum:=Sum + ABS(A[i,j]); END;
IF (Sum > Norm) THEN Norm:=Sum; END;
END;
END MatNormL1;
PROCEDURE ENorm(VAR X : ARRAY OF LONGREAL;
a,e : CARDINAL) : LONGREAL;
(*---------------------------------------------------------------*)
(* Argonne national laboratory. MINPACK project. March 1980. *)
(* Burton S. Garbow, Kenneth E. Hillstrom, Jorge J. More *)
(* *)
(* Anpassung an Modula-2 und "Uberarbeitung : M.Riedl, 23.11.93. *)
(*---------------------------------------------------------------*)
CONST MinLR = 3.834E-20;
MaxLR = 1.304E+19;
VAR i : CARDINAL;
S1,S2,S3,zw : LONGREAL;
Giant,AbsX,Max1,Max3 : LONGREAL;
Norm : LONGREAL;
BEGIN
S1:=0.0; S2:=0.0; S3:=0.0;
Max1:=0.0; Max3:=0.0;
Giant := MaxLR / VAL(LONGREAL,e-a+1);
FOR i:=a-1 TO e-1 DO
AbsX := ABS(X[i]);
IF (AbsX < MinLR) THEN (* Summiere kleine Elemente. *)
IF (AbsX > Max3) THEN
zw := Max3 / AbsX;
S3 := 1.0 + S3*(zw*zw);
Max3 := AbsX;
ELSE
IF (AbsX # 0.0) THEN zw := AbsX / Max3; S3:=S3 + zw*zw; END;
END;
ELSIF (AbsX > Giant) THEN (* Summiere gro\3e Elemente. *)
IF (AbsX > Max1) THEN
zw := Max1 / AbsX;
S1 := 1.0 + S1*(zw*zw);
Max1 := AbsX;
ELSE
zw := AbsX / Max1; S1:=S1 + zw*zw;
END;
ELSE (* Summiere mittlere Elemente. *)
S2 := S2 + AbsX*AbsX;
END;
END; (* FOR i *)
Norm:=0.0; (* Wg. XDS Compilerwarung, nicht noetig ! *)
IF (S1 # 0.0) THEN (* Berechne die Norm. *)
Norm := Max1*sqrt(S1 + (S2 / Max1) / Max1);
ELSE
IF (S2 # 0.0) THEN
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)));
END;
ELSE
Norm := Max3*sqrt(S3);
END;
END;
RETURN Norm;
END ENorm;
PROCEDURE NormVek(VAR X : ARRAY OF LONGREAL;
dim : CARDINAL;
MaxNorm : BOOLEAN);
CONST klein = 32;
VAR i : CARDINAL;
Norm,Max : LONGREAL;
BEGIN
IF NOT MaxNorm THEN
Norm := ENorm(X,1,dim);
ELSE
IF (dim > klein) THEN
Max:=ABS(X[idamax(dim,X[0],1)-1]);
ELSE
Max:=ABS(X[0]);
FOR i:=1 TO dim-1 DO
IF (ABS(X[i]) > Max) THEN Max:=ABS(X[i]); END;
END;
END;
Norm := Max;
END;
IF (Norm > (100.0/MAX(LONGREAL))) THEN
Norm := 1.0 / Norm;
IF (dim > klein) THEN
dscal(dim,Norm,X[0],1);
ELSE
FOR i:=0 TO dim-1 DO X[i]:=Norm*X[i]; END;
END;
END;
END NormVek;
PROCEDURE NormMat(VAR A : ARRAY OF ARRAY OF LONGREAL;
dim : CARDINAL;
VekAnz : CARDINAL;
MaxNorm : BOOLEAN);
VAR i : CARDINAL;
BEGIN
FOR i:=0 TO VekAnz-1 DO
NormVek(A[i],dim,MaxNorm);
END;
END NormMat;
PROCEDURE CNormMat(VAR Mat : ARRAY OF ARRAY OF LONGCOMPLEX;
dim : CARDINAL;
VekAnz : CARDINAL;
MaxNorm: BOOLEAN);
VAR i,j : CARDINAL;
Norm,Sum,Max,abs : LONGREAL;
BEGIN
Norm:=0.0;
FOR i:=0 TO VekAnz-1 DO
IF MaxNorm THEN
Max:=0.0;
FOR j:=0 TO dim-1 DO
IF ABS(RE(Mat[i,j])) > Max THEN Max:=ABS(RE(Mat[i,j])); END;
IF ABS(RE(Mat[i,j])) > Max THEN Max:=ABS(RE(Mat[i,j])); END;
END;
IF (Max # 0.0) THEN Norm := 1.0 / Max; END;
ELSE (* Euklidische Norm. *)
Sum:=0.0;
FOR j:=0 TO dim-1 DO
abs := LongComplexMath.abs(Mat[i,j]); Sum:=Sum + abs*abs;
END;
IF (Sum # 0.0) THEN Norm := 1.0 / sqrt(Sum); END;
END; (* IF *)
IF (Norm # 0.0) THEN
FOR j:=0 TO dim-1 DO
Mat[i,j]:=LongComplexMath.scalarMult(Norm,Mat[i,j]);
END;
END;
END; (* FOR i *)
END CNormMat;
PROCEDURE MatSpur(VAR Spur : LONGREAL;
VAR MAT : MATRIX;
dim : CARDINAL);
VAR i :CARDINAL;
BEGIN
Spur:=0.0; FOR i:=1 TO dim DO Spur:=Spur + MAT[i,i]; END;
END MatSpur;
PROCEDURE SvSpur(VAR Spur : LONGREAL;
VAR SV : SUPERVEKTOR;
dim : CARDINAL);
VAR i,ii : CARDINAL;
BEGIN
Spur:=0.0; ii:=0;
FOR i:=1 TO dim DO INC(ii,i); Spur:=Spur+SV[ii]; END;
END SvSpur;
PROCEDURE MinMaxVek(VAR Min,Max : LONGREAL;
VAR Vek : ARRAY OF LONGREAL;
dim : CARDINAL;
Form : CARDINAL);
VAR i : CARDINAL;
BEGIN
Min:=Vek[0]; Max:=Vek[0];
IF (Form = 1) THEN
FOR i:=1 TO dim-1 DO
IF (ABS(Vek[i]) < ABS(Min)) THEN Min:=Vek[i]; END;
IF (ABS(Vek[i]) > ABS(Max)) THEN Max:=Vek[i]; END;
END;
ELSE
FOR i:=1 TO dim-1 DO
IF (Vek[i] < Min) THEN Min:=Vek[i]; END;
IF (Vek[i] > Max) THEN Max:=Vek[i]; END;
END;
END; (* IF Form *)
END MinMaxVek;
PROCEDURE MinMaxMat(VAR Min,Max : LONGREAL;
VAR Mat : MATRIX;
dim : CARDINAL; (* Dimension der Matrix *)
Form : CARDINAL);
VAR i,j : CARDINAL;
BEGIN
IF (Form > 5) THEN Form:=0; END;
CASE Form OF
0 : Min:=Mat[1,1]; Max:=Mat[1,1];
FOR i:=1 TO dim DO
FOR j:=1 TO dim DO
IF (Mat[i,j] > Max) THEN Max:=Mat[i,j]; END;
IF (Mat[i,j] < Min) THEN Min:=Mat[i,j]; END;
END;
END; |
1 : Min:=Mat[1,1]; Max:=Mat[1,1];
FOR i:=1 TO dim DO
FOR j:=1 TO dim DO
IF (ABS(Mat[i,j]) > ABS(Max)) THEN Max:=Mat[i,j]; END;
IF (ABS(Mat[i,j]) < ABS(Min)) THEN Min:=Mat[i,j]; END;
END;
END; |
2 : Min:=Mat[2,1]; Max:=Mat[2,1];
FOR i:=1 TO dim DO
FOR j:=1 TO i-1 DO
IF (Mat[i,j] > Max) THEN Max:=Mat[i,j]; END;
IF (Mat[i,j] < Min) THEN Min:=Mat[i,j]; END;
END;
FOR j:=i+1 TO dim DO
IF (Mat[i,j] > Max) THEN Max:=Mat[i,j]; END;
IF (Mat[i,j] < Min) THEN Min:=Mat[i,j]; END;
END;
END; |
3 : Min:=Mat[2,1]; Max:=Mat[2,1];
FOR i:=1 TO dim DO
FOR j:=1 TO i-1 DO
IF (ABS(Mat[i,j]) > ABS(Max)) THEN Max:=Mat[i,j]; END;
IF (ABS(Mat[i,j]) < ABS(Min)) THEN Min:=Mat[i,j]; END;
END;
FOR j:=i+1 TO dim DO
IF (ABS(Mat[i,j]) > ABS(Max)) THEN Max:=Mat[i,j]; END;
IF (ABS(Mat[i,j]) < ABS(Min)) THEN Min:=Mat[i,j]; END;
END;
END; |
4 : Min:=Mat[1,1]; Max:=Mat[1,1];
FOR i:=2 TO dim DO
IF (Mat[i,i] > Max) THEN Max:=Mat[i,i]; END;
IF (Mat[i,i] < Min) THEN Min:=Mat[i,i]; END;
END; |
5 : Min:=Mat[1,1]; Max:=Mat[1,1];
FOR i:=2 TO dim DO
IF (ABS(Mat[i,i]) > ABS(Max)) THEN Max:=Mat[i,i]; END;
IF (ABS(Mat[i,i]) < ABS(Min)) THEN Min:=Mat[i,i]; END;
END;
END; (* CASE *)
END MinMaxMat;
PROCEDURE MinMaxSv(VAR Min,Max : LONGREAL;
VAR SV : ARRAY OF LONGREAL; (* SUPERVEKTOR *)
dim : CARDINAL; (* Dimension der Matrix *)
Form : CARDINAL);
VAR i,j,ij,ii,SVdim : CARDINAL;
BEGIN
SVdim:=(dim*(dim+1) DIV 2) - 1;
IF (Form > 5) THEN Form:=0; END;
CASE Form OF
0 : Min:=SV[0]; Max:=SV[0];
FOR ij:=1 TO SVdim DO
IF (SV[ij] > Max) THEN Max:=SV[ij]; END;
IF (SV[ij] < Min) THEN Min:=SV[ij]; END;
END; |
1 : Min:=SV[0]; Max:=SV[0];
FOR ij:=1 TO SVdim DO
IF (ABS(SV[ij]) > ABS(Max)) THEN Max:=SV[ij]; END;
IF (ABS(SV[ij]) < ABS(Min)) THEN Min:=SV[ij]; END;
END; |
2 : Min:=SV[1]; Max:=SV[1]; ij:=0;
FOR i:=2 TO dim DO
FOR j:=1 TO i-1 DO
INC(ij);
IF (SV[ij] > Max) THEN Max:=SV[ij]; END;
IF (SV[ij] < Min) THEN Min:=SV[ij]; END;
END;
INC(ij);
END; |
3 : Min:=SV[1]; Max:=SV[1]; ij:=0;
FOR i:=2 TO dim DO
FOR j:=1 TO i-1 DO
INC(ij);
IF (ABS(SV[ij]) > ABS(Max)) THEN Max:=SV[ij]; END;
IF (ABS(SV[ij]) < ABS(Min)) THEN Min:=SV[ij]; END;
END;
INC(ij);
END; |
4 : Min:=SV[0]; Max:=SV[0]; ii:=0;
FOR i:=2 TO dim DO
INC(ii,i);
IF (SV[ii] > Max) THEN Max:=SV[ii]; END;
IF (SV[ii] < Min) THEN Min:=SV[ii]; END;
END; |
5 : Min:=SV[0]; Max:=SV[0]; ii:=0;
FOR i:=2 TO dim DO
INC(ii,i);
IF (ABS(SV[ii]) > ABS(Max)) THEN Max:=SV[ii]; END;
IF (ABS(SV[ii]) < ABS(Min)) THEN Min:=SV[ii]; END;
END;
END; (* CASE *)
END MinMaxSv;
PROCEDURE ZufallVec(VAR V : ARRAY OF LONGREAL;
dim : CARDINAL);
VAR i : CARDINAL;
BEGIN
IF (dim = 0) THEN RETURN; END;
IF (dim-1 > HIGH(V)) THEN
Errors.ErrOut('dim > HIGH(MAT) (MatLib.ZufallVec)');
RETURN;
END;
FOR i:=0 TO dim-1 DO V[i] := 1.0 - 2.0*Zufall(); END;
END ZufallVec;
PROCEDURE ZufallMat(VAR A : ARRAY OF ARRAY OF LONGREAL;
m,n : CARDINAL;
sym : BOOLEAN);
VAR i,j : CARDINAL;
x : LONGREAL;
BEGIN
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 (m = 1) AND (n = 1) THEN A[0,0]:=1.0 - 2.0*Zufall(); RETURN; END;
IF sym THEN
FOR i:=0 TO n-1 DO
FOR j:=i+1 TO n-1 DO
x:=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 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 : ARRAY OF LONGREAL; (* SUPERVEKTOR *)
dim : CARDINAL);
VAR i,j,ij : CARDINAL;
BEGIN
(* InitZufall; *)
ij:=0;
FOR i:=1 TO dim DO
FOR j:=1 TO i DO SV[ij]:=1.0 - 2.0*Zufall(); INC(ij); END;
END;
END ZufallSv;
PROCEDURE SumVek(VAR X : ARRAY OF LONGREAL;
a,e : CARDINAL) : LONGREAL;
VAR i : CARDINAL;
s,c,y,t : LONGREAL;
BEGIN
s := X[a-1];
c := 0.0;
FOR i:=a TO e-1 DO
y := X[i] - c;
t := s + y;
c := (t - s) - y;
s := t;
END;
RETURN s;
END SumVek;
PROCEDURE AbsSumVek(VAR X : ARRAY OF LONGREAL;
a,e : CARDINAL) : LONGREAL;
VAR i : CARDINAL;
s,c,y,t : LONGREAL;
BEGIN
s := ABS(X[a-1]);
c := 0.0;
FOR i:=a TO e-1 DO
y := ABS(X[i]) - c;
t := s + y;
c := (t - s) - y;
s := t;
END;
RETURN s;
END AbsSumVek;
PROCEDURE NeumaierSum(VAR X : ARRAY OF LONGREAL;
N : CARDINAL) : LONGREAL;
VAR sum,c,t : LONGREAL;
i : CARDINAL;
BEGIN
sum := X[0];
c := 0.0;
FOR i:=1 TO N-1 DO
t := sum + X[i];
IF (ABS(sum) >= ABS(X[i])) THEN
c:=c + (sum - t) + X[i];
ELSE
c:=c + (X[i] - t) + sum;
END;
sum := t;
END;
RETURN sum + c;
END NeumaierSum;
PROCEDURE NeumaierProdSum(VAR X,Y : ARRAY OF LONGREAL;
N : CARDINAL) : LONGREAL;
VAR sum,c,t,xy : LONGREAL;
i : CARDINAL;
BEGIN
sum := X[0]*Y[0];
c := 0.0;
FOR i:=1 TO N-1 DO
xy := X[i]*Y[i];
t := sum + xy;
IF (ABS(sum) >= ABS(xy)) THEN
c:=c + (sum - t) + xy;
ELSE
c:=c + (xy - t) + sum;
END;
sum := t;
END;
RETURN sum + c;
END NeumaierProdSum;
PROCEDURE SumVekUR(VAR X : ARRAY OF LONGREAL;
a,e : CARDINAL) : LONGREAL;
VAR i,mod,n : CARDINAL;
unrolled : BOOLEAN;
s : LONGREAL;
BEGIN
n:=e-a+1;
unrolled := n > 32;
IF NOT unrolled THEN
s:=X[a-1];
FOR i:=a TO e-1 DO s:=s + X[i]; END;
ELSE
s:=0.0;
mod := (n MOD 4);
IF (mod # 0) THEN
FOR i:=a-1 TO a+mod-2 DO s:=s + X[i]; END;
END;
FOR i:=a+mod-1 TO e-1 BY 4 DO
s:=s+ X[i] + X[i+1] + X[i+2] + X[i+3];
END;
END; (* IF unrolled *)
RETURN s;
END SumVekUR;
PROCEDURE AbsSumVekUR(VAR X : ARRAY OF LONGREAL;
a,e : CARDINAL) : LONGREAL;
VAR i,mod,n : CARDINAL;
unrolled : BOOLEAN;
s : LONGREAL;
BEGIN
n:=e-a+1;
unrolled := n > 32;
IF NOT unrolled THEN
s:=ABS(X[a-1]);
FOR i:=a TO e-1 DO s:=s + ABS(X[i]); END;
ELSE
s:=0.0;
mod := (n MOD 4);
IF (mod # 0) THEN
FOR i:=a-1 TO a+mod-2 DO s:=s + ABS(X[i]); END;
END;
FOR i:=a+mod-1 TO e-1 BY 4 DO
s:=s+ ABS(X[i]) + ABS(X[i+1]) + ABS(X[i+2]) + ABS(X[i+3]);
END;
END; (* IF unrolled *)
RETURN s;
END AbsSumVekUR;
END MatLib.