IMPLEMENTATION MODULE PMatLib;
(*------------------------------------------------------------------------*)
(* Einige Matrix-Operationen f"ur den Datentyp PMATRIX *)
(* Some matrix operations for data typ Deklera.PMATRIX *)
(*------------------------------------------------------------------------*)
(* Letzte Bearbeitung: *)
(* *)
(* 30.09.96, MRi: Durchsicht *)
(* 15.09.15, MRi: Anpassung an ISO M2, Bereinigung von nicht ben"otigten *)
(* Importen. *)
(*------------------------------------------------------------------------*)
(* Offene Punkte: *)
(* *)
(* - PMatVek und PMatMat nochmals gegen LibDBlas dgemv & dgemm abgelichen *)
(*------------------------------------------------------------------------*)
(* Implementation : Michael Riedl *)
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
(* $Id: PMatLib.mod,v 1.3 2016/04/15 15:42:59 mriedl Exp mriedl $ *)
FROM SYSTEM IMPORT TSIZE;
FROM Deklera IMPORT PMATRIX;
IMPORT Errors;
FROM DynMat IMPORT AllocMat,DeAllocMat;
FROM MatLib IMPORT ENorm;
PROCEDURE NormVek2(VAR X : ARRAY OF LONGREAL;
IStart,IEnd : CARDINAL;
MaxNorm : BOOLEAN);
(*-------------------------------------------------------------*)
(* Normiert den Vektor X im Bereich X[IStart] bis X[IEnd]. *)
(* Wenn MaxNorm = TRUE so, da3 das Maximale Element die Gr"o3e *)
(* 1 hat, ansonsten nach der Euklidschen Norm auf 1 *)
(*-------------------------------------------------------------*)
VAR i : CARDINAL;
Norm,Max : LONGREAL;
BEGIN
(* TIO.WrCard(IStart,7); TIO.WrCard(IEnd,7); TIO.WrLn; *)
IF NOT MaxNorm THEN
Norm := 1.0 / ENorm(X,IStart,IEnd);
FOR i:=IStart-1 TO IEnd-1 DO X[i]:=Norm*X[i]; END;
ELSE
Max:=X[IStart-1];
FOR i:=IStart TO IEnd-1 DO
IF (X[i] > Max) THEN Max:=X[i]; END;
END;
IF (Max # 0.0) THEN
Norm := 1.0 / Max;
FOR i:=IStart-1 TO IEnd-1 DO X[i]:=X[i] * Norm; END;
END;
END;
END NormVek2;
PROCEDURE PMatVek(VAR A : PMATRIX;
m,n : CARDINAL;
Alpha : LONGREAL;
VAR X : ARRAY OF LONGREAL;
Beta : LONGREAL;
VAR Y : ARRAY OF LONGREAL;
TransA : CHAR);
VAR Temp : LONGREAL;
i,j,max : CARDINAL;
BEGIN
IF (CAP(TransA) = "N") THEN (* Y = A x *)
max := m-1;
ELSE
max := n-1;
END;
IF (Beta = 0.0) THEN
FOR i:=0 TO max DO Y[i]:=0.0; END;
ELSE
FOR i:=0 TO max DO Y[i]:=Beta*Y[i]; END;
END;
IF (CAP(TransA) = "N") THEN (* Y = A x *)
FOR i:=1 TO m DO
Temp := 0.0;
FOR j:=1 TO n DO
Temp:=Temp + A^[i]^[j]*X[j-1];
END;
Y[i-1]:=Y[i-1] + Alpha*Temp;
END;
ELSE (* Y = A' x *)
FOR i:=1 TO n DO
Temp := 0.0;
FOR j:=1 TO m DO
Temp:=Temp + A^[j]^[i]*X[j-1];
END;
Y[i-1]:=Y[i-1] + Alpha*Temp;
END;
END;
END PMatVek;
PROCEDURE PMatMult(VAR C : PMATRIX;
VAR A : PMATRIX;
VAR B : PMATRIX;
ma,na : CARDINAL;
mb,nb : CARDINAL;
Alpha : LONGREAL;
Beta : LONGREAL;
TransA : CHAR;
TransB : CHAR);
CONST FFlag1 = "Matrixdimensionen nicht kompatibel (PMatMult) !)";
FFlag2 = "Kein Freispeicher vorhanden (PMatMult) !";
VAR i,j,k,N,M : CARDINAL;
S,Temp : LONGREAL;
Ctemp : PMATRIX;
TrA,TrB : BOOLEAN;
BEGIN
TrA := CAP(TransA) = "T";
TrB := CAP(TransB) = "T";
M:=0; N:=0; (* Wg. Compilerwarnung *)
IF NOT TrA AND NOT TrB THEN (* C = Alpha (A B) + Beta C *)
IF (na # mb) THEN Errors.FatalError(FFlag1); END;
M:=ma; N:=nb;
ELSIF TrA AND NOT TrB THEN (* C = Alpha (A'B) + Beta C *)
IF (ma # mb) THEN Errors.FatalError(FFlag1); END;
M:=na; N:=nb;
ELSIF NOT TrA AND TrB THEN (* C = Alpha (A B') + Beta C *)
IF (na # nb) THEN Errors.FatalError(FFlag1); END;
M:=ma; N:=mb;
ELSIF TrA AND TrB THEN (* C = Alpha (A'B') + Beta C *)
IF (ma # nb) THEN Errors.FatalError(FFlag1); END;
M:=na; N:=mb;
END;
IF (Alpha = 0.0) THEN
IF (Beta = 0.0) THEN
FOR i:=1 TO M DO
FOR j:=1 TO N DO C^[i]^[j]:=0.0; END;
END;
ELSE
FOR i:=1 TO M DO
FOR j:=1 TO N DO C^[i]^[j] := Beta*C^[i]^[j]; END;
END;
END;
ELSE (* Alpha # 0.0 *)
AllocMat(Ctemp,M,N);
IF (Ctemp = NIL) THEN
Errors.FatalError(FFlag2);
END;
(**** F"uhre die Berechnungen durch ! ****)
IF NOT TrA AND NOT TrB THEN (* C = Alpha (A B) + Beta C *)
FOR i:=1 TO ma DO (* Schleifenstruktur OK *)
IF (Beta = 0.0) THEN
FOR j:=1 TO nb DO Ctemp^[i]^[j] := 0.0; END;
ELSE
FOR j:=1 TO nb DO Ctemp^[i]^[j] := Beta*C^[i]^[j]; END;
END;
FOR j:=1 TO na DO
IF (A^[i]^[j] # 0.0) THEN
Temp := Alpha*A^[i]^[j];
FOR k:=1 TO nb DO
Ctemp^[i]^[k]:=Ctemp^[i]^[k] + Temp*B^[j]^[k];
END;
END;
END;
END; (* FOR i *)
ELSIF TrA AND NOT TrB THEN (* C = Alpha (A'B) + Beta C *)
FOR i:=1 TO na DO (* C = Alpha (A' B) + Beta C *)
FOR j:=1 TO nb DO
S:=0.0;
FOR k:=1 TO mb DO
S:=S + A^[k]^[i]*B^[k]^[j];
END;
IF (Beta = 0.0) THEN
Ctemp^[i]^[j] := Alpha*S;
ELSE
Ctemp^[i]^[j] := Alpha*S + Beta*C^[i]^[j];
END;
END;
END; (* FOR i *)
ELSIF NOT TrA AND TrB THEN (* C = Alpha (A B') + Beta C *)
FOR i:=1 TO ma DO (* Schleifenstruktur OK *)
FOR j:=1 TO mb DO
S:=0.0;
FOR k:=1 TO nb (*na*) DO
S:=S + A^[i]^[k]*B^[j]^[k];
END;
IF (Beta = 0.0) THEN
Ctemp^[i]^[j] := Alpha*S;
ELSE
Ctemp^[i]^[j] := Alpha*S + Beta*C^[i]^[j];
END;
END;
END; (* FOR i *)
ELSIF TrA AND TrB THEN (* C = Alpha (A'B') + Beta C *)
FOR i:=1 TO na DO
FOR j:=1 TO mb DO
S:=0.0;
FOR k:=1 TO nb DO
S:=S + A^[k]^[i]*B^[j]^[k];
END;
IF (Beta = 0.0) THEN
Ctemp^[i]^[j] := Alpha*S;
ELSE
Ctemp^[i]^[j] := Alpha*S + Beta*C^[i]^[j]
END;
END;
END;
ELSE
Errors.FatalError("Unzulaessige Option (PMatMul) !");
END; (* IF *)
FOR i:=1 TO M DO
FOR j:=1 TO N DO
C^[i]^[j]:=Ctemp^[i]^[j];
END;
END;
DeAllocMat(Ctemp,M,N);
END; (* IF Alpha *)
END PMatMult;
PROCEDURE MinMaxNorm(VAR X : ARRAY OF LONGREAL;
n : CARDINAL);
VAR i : CARDINAL;
Min,Max,Norm : LONGREAL;
BEGIN
Min := X[0]; Max := X[0];
FOR i:=1 TO n-1 DO
IF (X[i] < Min) THEN Min:=X[i]; END;
IF (X[i] > Max) THEN Max:=X[i]; END;
END;
Norm := 1.0 / (Max - Min);
FOR i:=0 TO n-1 DO
X[i]:=(X[i] - Min) * Norm;
END;
(*
Min := X[0];
FOR i:=1 TO n-1 DO
IF (X[i] < Min) THEN Min:=X[i]; END;
END;
FOR i:=0 TO n-1 DO
X[i]:=X[i] - Min;
END;
Max := X[0];
FOR i:=1 TO n-1 DO
IF (X[i] > Max) THEN Max:=X[i]; END;
END;
Norm := 1.0 / Max;
FOR i:=0 TO n-1 DO
X[i]:=X[i] * Norm;
END;
*)
END MinMaxNorm;
PROCEDURE dger( m,n : CARDINAL;
Alpha : LONGREAL;
VAR X : ARRAY OF LONGREAL;
IncX : CARDINAL;
VAR Y : ARRAY OF LONGREAL;
IncY : CARDINAL;
VAR A : PMATRIX);
VAR i,j,ix,jy : CARDINAL;
AlphaY : LONGREAL;
BEGIN
IF (IncX = 1) AND (IncY = 1) THEN
FOR i:=1 TO m DO
IF (X[i-1] # 0.0) THEN
AlphaY := Alpha*X[i-1];
FOR j:=1 TO n DO
A^[i]^[j]:=A^[i]^[j] + AlphaY*Y[j-1];
END;
END;
END;
ELSE
Errors.ErrOut("Inkremente ungleich 1 nicht getestet (dger) !");
ix := 0;
FOR i:=1 TO m DO
IF (X[i-1] # 0.0) THEN
AlphaY := Alpha*X[ix];
jy := 0;
FOR j:=1 TO n DO
A^[i]^[j]:=A^[i]^[j] + AlphaY*Y[jy];
INC(jy,IncY);
END;
END;
INC(ix,IncX);
END;
END;
END dger;
END PMatLib.