IMPLEMENTATION MODULE LibDBlasM2;
(*------------------------------------------------------------------------*)
(* Basic Linear Algebra Prozeduren in Modula-2 mit M2-konformen *)
(* Aufrufen der Parameter. *)
(* Some level 1,2 & 3 BLAS routines implemented in Modula-2 pathing the *)
(* arguments in a "M2 manner". *)
(*------------------------------------------------------------------------*)
(* Letzte Bearbeitung *)
(* *)
(* 23.11.93, MRi: Erstellen der 1. Version *)
(* Aug. 95, MRi: Erweiterung *)
(* 03.11.95, MRi: Durchsicht *)
(* 13.10.15, MRi: Ersetzen von LFLOAT(#) durch VAL(LONGREAL,#) *)
(* 30.11.15, MRi: Erstelle der Routinen dscal,daxpy,idmax aus Linpack *)
(* Benchmark aus Beispielen von Stony Brook M2 *)
(* 01.12.15, MRi: Umbenennen von BlasLib auf LibDBlas, einfuegen von *)
(* dnrm2 (aus dnrm2.f) *)
(* 28.01.16, MRi: Umstellen von dgemv und dgemm auf "Open Array". *)
(* und anpassen wegen zeilenorientiertem Speichermodell *)
(* 01.02.16, MRi: dgemm partiell neu kodiert, alle Tests fehlerfrei *)
(* 02.02.16, MRi: dgemv partiell neu kodiert, alle Tests fehlerfrei *)
(* 30.03.16, MRi: dswap eingefuegt *)
(* 31.03.16, MRi: drot und drotg eingefuegt *)
(* 13.04.16, MRi: dcopy eingefuegt *)
(* 21.10.16, MRi: idamin eingefuegt *)
(* 28.10.17, MRi: zgemm eingefuegt *)
(* 29.10.17, MRi: Rolle von M,N in dgemm und dgemv an zgemm angepasst *)
(*------------------------------------------------------------------------*)
(* Testroutinen *)
(* *)
(* - zgemm in TstCmplxMaMul *)
(*------------------------------------------------------------------------*)
(* Offene Punkte: *)
(* *)
(* - Weiteres austesten von dgemm und besonders dgemv *)
(* - Austesten von drot und drotg *)
(*------------------------------------------------------------------------*)
(* Implementation : Michael Riedl *)
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
(* $Id: LibDBlasM2.mod,v 1.9 2017/10/29 09:55:11 mriedl Exp mriedl $ *)
IMPORT SYSTEM;
FROM Deklera IMPORT FLOAT; (* REAL type *)
FROM LongMath IMPORT sqrt;
FROM LongComplexMath IMPORT conj;
IMPORT Errors;
TYPE PVEKTOR = POINTER TO ARRAY [0..MAX(INTEGER)-1] OF FLOAT;
PROCEDURE IMax(a,b : INTEGER) : INTEGER;
BEGIN
IF (a > b) THEN RETURN a; ELSE RETURN b; END;
END IMax;
PROCEDURE SumVek(VAR X : ARRAY OF FLOAT;
a,e : CARDINAL) : FLOAT;
VAR i : CARDINAL;
s,c,y,t : FLOAT;
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 FLOAT;
a,e : CARDINAL) : FLOAT;
VAR i : CARDINAL;
s,c,y,t : FLOAT;
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 ENorm(VAR X : ARRAY OF FLOAT;
a,e : CARDINAL) : FLOAT;
(*---------------------------------------------------------------*)
(* 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 : FLOAT;
Giant,AbsX,Max1,Max3 : FLOAT;
Norm : FLOAT;
BEGIN
S1:=0.0; S2:=0.0; S3:=0.0;
Max1:=0.0; Max3:=0.0;
Giant := MaxLR / VAL(FLOAT,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;
(*========================================================================*)
(* Pure BLAS based routines *)
(*========================================================================*)
PROCEDURE dnrm2( n : INTEGER;
VAR X : ARRAY OF FLOAT;
IncX : INTEGER): FLOAT;
(*-----------------------------------------------------------------*)
(* This version written on 25-October-1982. *)
(* Modified on 14-October-1993 to inline the call to DLASSQ. *)
(* Sven Hammarling, Nag Ltd. *)
(* Adapted to Modula-2 01.12.2015, M.Riedl *)
(*-----------------------------------------------------------------*)
VAR i,ix : INTEGER;
absxi,norm,scale,ssq : FLOAT;
zw : FLOAT;
BEGIN
IF (n < 1) OR (IncX < 1) THEN
norm := 0.0
ELSIF (n = 1) THEN
norm := ABS(X[0]);
ELSE
scale := 0.0;
ssq := 1.0;
ix := 0;
FOR i:=0 TO (n-1) DO
IF (X[ix] # 0.0) THEN
absxi := ABS(X[ix]);
IF (scale < absxi) THEN
zw := scale / absxi;
ssq := 1.0 + ssq*zw*zw;
scale := absxi;
ELSE
zw := absxi / scale;
ssq := ssq + zw*zw;
END
END;
INC(ix,IncX);
END;
norm := scale*sqrt(ssq);
END;
RETURN norm;
END dnrm2;
PROCEDURE dswap( N : CARDINAL;
VAR X : ARRAY OF FLOAT;
incX : INTEGER;
VAR Y : ARRAY OF FLOAT;
incY : INTEGER);
(*----------------------------------------------------------------*)
(* Adopted to Modula-2, MRi, 30.03.2016 *)
(*----------------------------------------------------------------*)
CONST level = 4;
VAR Xi : FLOAT;
Xtmp : ARRAY [0..level-1] OF FLOAT;
ix,iy : INTEGER;
i,m : CARDINAL;
BEGIN
IF (N = 0) THEN RETURN; END;
IF (incX = 1) AND (incY = 1) THEN
m := N MOD level;
IF (m # 0) THEN
FOR i:=0 TO m-1 DO (* Clean up loop *)
Xi := X[i]; X[i] := Y[i]; Y[i] := Xi;
END;
END;
IF (N < level) THEN RETURN; END;
FOR i:=m TO N-1 BY level DO
Xtmp[0] := X[i+0]; X[i+0] := Y[i+0]; Y[i+0] := Xtmp[0];
Xtmp[1] := X[i+1]; X[i+1] := Y[i+1]; Y[i+1] := Xtmp[1];
Xtmp[2] := X[i+2]; X[i+2] := Y[i+2]; Y[i+2] := Xtmp[2];
Xtmp[3] := X[i+3]; X[i+3] := Y[i+3]; Y[i+3] := Xtmp[3];
END;
ELSE (* code for unequal increments or equal increments # 1 *)
IF (incX > 0) THEN ix:=0; ELSE ix:=(1 - VAL(INTEGER,N))*incX; END;
IF (incY > 0) THEN iy:=0; ELSE iy:=(1 - VAL(INTEGER,N))*incY; END;
FOR i:=0 TO N-1 DO
Xi := X[ix]; X[ix] := Y[iy]; Y[iy] := Xi;
INC(ix,incX); INC(iy,incY);
END;
END;
END dswap;
PROCEDURE dcopy( N : INTEGER;
VAR X : ARRAY OF FLOAT;
IncX : INTEGER;
VAR Y : ARRAY OF FLOAT;
IncY : INTEGER);
(*----------------------------------------------------------------*)
(* Adopted to Modula-2, MRi, 04.04.2016 *)
(*----------------------------------------------------------------*)
VAR i,ix,iy,m : INTEGER;
BEGIN
IF (N <= 0) THEN RETURN; END;
IF (IncX = 1) AND (IncY = 1) THEN
(* code for both increments equal to 1 *)
m := (N MOD 8);
IF (m # 0) THEN (* Clean-up loop *)
FOR i:=0 TO m-1 DO Y[i] := X[i]; END;
IF (N < 8) THEN RETURN; END
END;
FOR i:=m TO N-1 BY 8 DO
Y[i+0] := X[i+0]; Y[i+1] := X[i+1];
Y[i+2] := X[i+2]; Y[i+3] := X[i+3];
Y[i+4] := X[i+4]; Y[i+5] := X[i+5];
Y[i+6] := X[i+6]; Y[i+7] := X[i+7];
END;
ELSE
(* code for unequal increments or equal increments not equal to 1 *)
IF (IncX > 0) THEN ix:=0; ELSE ix:=(1 - VAL(INTEGER,N))*IncX; END;
IF (IncY > 0) THEN iy:=0; ELSE iy:=(1 - VAL(INTEGER,N))*IncY; END;
FOR i:=0 TO N-1 DO
Y[iy] := X[ix];
INC(ix,IncX); INC(iy,IncY);
END;
END;
END dcopy;
PROCEDURE drot( N : INTEGER;
VAR X : ARRAY OF FLOAT;
incX : INTEGER;
VAR Y : ARRAY OF FLOAT;
incY : INTEGER;
c,s : FLOAT);
(*----------------------------------------------------------------*)
(* Adopted to Modula-2, MRi, 31.03.2016 *)
(*----------------------------------------------------------------*)
VAR i,ix,iy : INTEGER;
x1,y1 : FLOAT;
BEGIN
IF (incX = 1) AND (incY = 1) THEN
FOR i:=0 TO N-1 DO
x1 := X[i];
y1 := Y[i];
X[i] := c * x1 + s * y1;
Y[i] := -s * x1 + c * y1;
END;
ELSE
IF (incX > 0) THEN ix:=0; ELSE ix:=(1 - VAL(INTEGER,N))*incX; END;
IF (incY > 0) THEN iy:=0; ELSE iy:=(1 - VAL(INTEGER,N))*incY; END;
FOR i:=0 TO N-1 DO
x1 := X[ix];
y1 := Y[iy];
X[ix] := c * x1 + s * y1;
Y[iy] := -s * x1 + c * y1;
ix := ix + incX;
iy := iy + incY;
END;
END;
END drot;
PROCEDURE drotg(VAR da : FLOAT;
VAR db : FLOAT;
VAR c : FLOAT;
VAR s : FLOAT);
(*----------------------------------------------------------------*)
(* Adopted to Modula-2, MRi, 31.03.2016 *)
(*----------------------------------------------------------------*)
VAR r,z,roe,daos,dbos,scale : FLOAT;
BEGIN
IF (ABS(da) > ABS(db)) THEN roe := da; ELSE roe := db; END;
scale := ABS(da) + ABS(db);
IF (scale # 0.0) THEN
daos := da / scale;
dbos := db / scale;
r := scale*sqrt(daos*daos + dbos*dbos);
IF (roe < 0.0) THEN r:=-r; END; (* r:=DSIGN(1.0,roe)*r; *)
c := da / r;
s := db / r;
z := 1.0;
(*
IF( DABS(DA) .GT. DABS(DB) ) Z = S
IF( DABS(DB) .GE. DABS(DA) .AND. C .NE. 0.0D0 ) Z = 1.0D0/C
*)
IF (ABS(da) > ABS(db)) THEN
z := s;
ELSIF (c # 0.0) THEN
z := 1.0 / c;
END;
ELSE
c := 1.0; s := 0.0;
r := 0.0; z := 0.0;
END;
da := r; db := z;
END drotg;
PROCEDURE dscal( n : INTEGER;
da : FLOAT;
VAR dx : ARRAY OF FLOAT;
incx : INTEGER);
(*----------------------------------------------------------------*)
(* Scales a vector by a constant (UNROLLED version) *)
(*----------------------------------------------------------------*)
VAR i,m : INTEGER;
BEGIN
IF (n <= 0) THEN
RETURN;
END;
IF (incx <> 1) THEN (* code for increment not equal to 1 *)
FOR i:=0 TO n-1 DO
dx[i*incx] := da*dx[i*incx];
END;
ELSE (* code for increment equal to 1 *)
m := n REM 5;
IF (m <> 0) THEN
FOR i:=0 TO m-1 DO dx[i] := da*dx[i]; END;
IF (n < 5) THEN
RETURN;
END;
END;
FOR i:=m TO n-1 BY 5 DO
dx[i] := da*dx[i];
dx[i+1] := da*dx[i+1];
dx[i+2] := da*dx[i+2];
dx[i+3] := da*dx[i+3];
dx[i+4] := da*dx[i+4];
END;
END;
END dscal;
PROCEDURE daxpy( n : INTEGER;
da : FLOAT;
VAR dx : ARRAY OF FLOAT;
incx : INTEGER;
VAR dy : ARRAY OF FLOAT;
incy : INTEGER);
(*----------------------------------------------------------------*)
(* constant times a vector plus a vector *)
(*----------------------------------------------------------------*)
VAR i,ix,iy,m : INTEGER;
BEGIN
IF (n <= 0 ) THEN RETURN; END;
IF (da = 0.0) THEN RETURN; END;
IF ((incx <> 1) OR (incy <> 1)) THEN
(* code for unequal increments or equal increments <> 1 *)
ix := 1;
iy := 1;
IF (incx < 0) THEN ix := (-n+1)*incx + 1; END;
IF (incy < 0) THEN iy := (-n+1)*incy + 1; END;
FOR i:=0 TO n-1 DO
dy[iy] := dy[iy] + da*dx[ix];
INC(ix, incx);
INC(iy, incy);
END;
ELSE (* code for both increments equal to 1 *)
m := n REM 4;
IF (m <> 0) THEN
FOR i:=0 TO m-1 DO
dy[i] := dy[i] + da*dx[i];
END;
IF (n < 4) THEN RETURN; END;
END;
FOR i:=m TO n-1 BY 4 DO
dy[i+0] := dy[i+0] + da*dx[i+0];
dy[i+1] := dy[i+1] + da*dx[i+1];
dy[i+2] := dy[i+2] + da*dx[i+2];
dy[i+3] := dy[i+3] + da*dx[i+3];
END;
END;
END daxpy;
PROCEDURE ddot( N : INTEGER;
VAR X : ARRAY OF FLOAT;
IncX : INTEGER;
VAR Y : ARRAY OF FLOAT;
IncY : INTEGER) : FLOAT;
(*----------------------------------------------------------------*)
(* Forms the dot product of two vectors. Uses unrolled loops for *)
(* increments equal to one. *)
(* Adopted to Modula-2, MRi, 06.09.2015 *)
(*----------------------------------------------------------------*)
CONST veclen = 8;
VAR dtemp : FLOAT;
i,ix,iy,m : INTEGER;
BEGIN
IF (N <= 0) THEN RETURN 0.0; END;
dtemp := 0.0;
IF (IncX = 1) AND (IncY = 1) THEN
(* code for both increments equal to 1 *)
m := N MOD veclen;
IF (m # 0) THEN
FOR i:=0 TO m-1 DO (* clean-up loop *)
dtemp:=dtemp + X[i]*Y[i];
END;
IF (N < veclen) THEN RETURN dtemp; END;
END;
(* i := m - veclen; *)
FOR i:=m TO N-1 BY veclen DO
dtemp:=dtemp + X[i+0]*Y[i+0] + X[i+1]*Y[i+1] +
X[i+2]*Y[i+2] + X[i+3]*Y[i+3] +
X[i+4]*Y[i+4] + X[i+5]*Y[i+5] +
X[i+6]*Y[i+6] + X[i+7]*Y[i+7];
END;
ELSE
(* code for unequal increments or equal increments not equal to 1 *)
ix := 0; IF (IncX < 0) THEN ix := (1-N)*IncX; END;
iy := 0; IF (IncY < 0) THEN iy := (1-N)*IncY; END;
FOR i:=1 TO N DO
dtemp:=dtemp + X[ix]*Y[iy];
INC(ix,IncX); INC(iy,IncY);
END;
END;
RETURN dtemp;
END ddot;
PROCEDURE idamax( n : INTEGER;
VAR dx : ARRAY OF FLOAT;
incx : INTEGER) : INTEGER;
(*----------------------------------------------------------------*)
(* Finds the index of element having max. absolute value. *)
(*----------------------------------------------------------------*)
VAR dmax : FLOAT;
i,ix,itemp : INTEGER;
BEGIN
IF (n < 1) THEN RETURN -1; END;
IF (n = 1) THEN RETURN 0; END;
itemp:=0;
IF (incx <> 1) THEN (* code for increment not equal to 1 *)
dmax := ABS(dx[0]);
ix := 1 + incx;
FOR i:=1 TO n-1 DO
IF (ABS(dx[ix]) > dmax) THEN
itemp := i;
dmax := ABS(dx[ix]);
END;
INC(ix,incx);
END;
ELSE (* code for increment equal to 1 *)
dmax := ABS(dx[0]);
FOR i:=1 TO n-1 DO
IF (ABS(dx[i]) > dmax) THEN
itemp := i;
dmax := ABS(dx[i]);
END;
END;
END;
RETURN itemp;
END idamax;
PROCEDURE idamin( n : INTEGER;
VAR X : ARRAY OF FLOAT;
IncX : INTEGER) : INTEGER;
(*----------------------------------------------------------------*)
(* Finds the index of element having min. absolute value. *)
(*----------------------------------------------------------------*)
VAR dmin : FLOAT;
i,ix,itemp : INTEGER;
XX : PVEKTOR;
BEGIN
IF (n < 1) THEN RETURN -1; END;
IF (n = 1) THEN RETURN 0; END;
XX:=SYSTEM.CAST(PVEKTOR,SYSTEM.ADR(X));
itemp:=0;
IF (IncX <> 1) THEN (* code for Increment not equal to 1 *)
dmin := ABS(XX^[0]);
ix := 1 + IncX;
FOR i:=1 TO n-1 DO
IF (ABS(XX^[ix]) < dmin) THEN
itemp := i;
dmin := ABS(XX^[ix]);
END;
INC(ix,IncX);
END;
ELSE (* code for increment equal to 1 *)
dmin := ABS(XX^[0]);
FOR i:=1 TO n-1 DO
IF (ABS(XX^[i]) < dmin) THEN
itemp := i;
dmin := ABS(XX^[i]);
END;
END;
END;
RETURN itemp + 1;
END idamin;
PROCEDURE dasum( dim : CARDINAL;
VAR X : ARRAY OF FLOAT;
Inc : CARDINAL; (* Inkrementwert >= 1*)
Unrolled : BOOLEAN) : FLOAT;
CONST max = 6; (* 8 *)
VAR i,m : CARDINAL;
Sum : FLOAT;
BEGIN
IF (dim = 0) THEN RETURN 0.0; END;
IF (Inc < 1) THEN
Errors.FatalError("Inkrement kleiner 1 (LibDBlas.dasum) !");
END;
Sum:=0.0;
IF (Inc = 1) THEN
IF Unrolled THEN
(* n:=dim DIV max; *)
m:=dim DIV max;
IF (m > 0) THEN
FOR i:=0 TO m-1 DO Sum:=Sum + ABS(X[i]); END;
END;
FOR i:=m TO dim-1 BY max DO
Sum:=Sum + ABS(X[i ]) + ABS(X[i+1])
+ ABS(X[i+2]) + ABS(X[i+3])
+ ABS(X[i+4]) + ABS(X[i+5])
(* + ABS(X[i+6]) + ABS(X[i+7])*);
END;
ELSE
FOR i:=0 TO dim-1 DO Sum:=Sum + ABS(X[i]); END;
END;
ELSE
i:=0;
WHILE (i < dim) DO
Sum:=Sum + ABS(X[i]);
INC(i,Inc);
END;
END;
RETURN Sum;
END dasum;
PROCEDURE dgemv ( Trans : CHAR;
M,N : INTEGER;
Alpha : FLOAT;
VAR A : ARRAY OF ARRAY OF FLOAT; (* [1..LDA][1..*] *)
LDA : INTEGER;
VAR X : ARRAY OF FLOAT; (* [1..*] *)
IncX : INTEGER;
Beta : FLOAT;
VAR Y : ARRAY OF FLOAT; (* [1..*] *)
IncY : INTEGER);
VAR Temp : FLOAT;
i,j,iy,jx,jy,kx,ky,Info,LenX,LenY : INTEGER;
BEGIN (* Testen *)
(* Test the input parameters. *)
IF (CAP(Trans) = 'C') THEN
Info := 1;
ELSIF (M < 0) THEN
Info := 2;
ELSIF (N < 0) THEN
Info := 3;
ELSIF (LDA < IMax(1,M)) THEN
Info := 6;
ELSIF (IncX = 0) THEN
Info := 8;
ELSIF (IncY = 0) THEN
Info := 11;
ELSE
Info := 0;
END;
IF (Info <> 0) THEN (* Fehlerbehandlung ! *)
Errors.WriteLn; Errors.WriteLn;
Errors.WriteString("Fehler in LibDBlas.dgemv, Info = ");
Errors.WriteInt(Info);
Errors.WriteLn; Errors.WriteLn;
RETURN;
END;
(* Quick return if possible. *)
IF (M = 0) OR (N = 0) OR ((Alpha = 0.0) AND (Beta = 1.0)) THEN
RETURN;
END;
(* Set LenX and LenY, the lengths of the vectors x and y, and set *)
(* up the start points in X and Y. *)
IF (CAP(Trans) = 'N') THEN
LenX := N;
LenY := M;
ELSE
LenX := M;
LenY := N;
END;
IF (IncX > 0) THEN
kx := 1;
ELSE
kx := 1 - (LenX - 1)*IncX;
END;
IF (IncY > 0) THEN
ky := 1;
ELSE
ky := 1 - (LenY - 1)*IncY;
END;
(* Start the operations. In this version the elements of A are *)
(* accessed sequentially with one pass through A. *)
IF (Beta # 1.0) THEN (* First form y := beta*y. *)
IF (IncY = 1) THEN
IF (Beta = 0.0) THEN
FOR i:=0 TO LenY-1 DO Y[i]:=0.0; END;
ELSE
FOR i:=0 TO LenY-1 DO Y[i]:=Y[i]*Beta; END;
END;
ELSE
iy := ky - 1;
IF (Beta = 0.0) THEN
FOR i:=0 TO LenY-1 DO Y[iy]:=0.0; INC(iy,IncY); END;
ELSE
FOR i:=0 TO LenY-1 DO Y[iy]:=Y[iy]*Beta; INC(iy,IncY); END;
END;
END;
END;
IF (Alpha = 0.0) THEN RETURN; END;
(* ****************************** *)
IF (CAP(Trans) = 'N') THEN (* Form y := alpha*A*x + y. *)
jy := ky-1;
IF (IncY = 1) THEN
FOR i:=0 TO M-1 DO
Temp := 0.0;
FOR j:=0 TO N-1 DO Temp:=Temp + A[i,j]*X[j]; END;
Y[jy]:=Y[jy] + Alpha*Temp;
INC(jy,IncY);
END;
ELSE
FOR i:=0 TO M-1 DO
Temp := 0.0;
jx := kx-1;
FOR j:=0 TO N-1 DO
Temp:=Temp + A[i,j]*X[jx];
INC(jx,IncX);
END;
Y[jy]:=Y[jy] + Alpha*Temp;
INC(jy,IncY);
END;
END;
ELSE (* Form y := alpha*A'*x + y. *)
jx := kx-1;
IF (IncX = 1) THEN
FOR j:=0 TO M-1 DO
IF (X[jx] # 0.0) THEN
Temp := Alpha*X[jx];
FOR i:=0 TO N-1 DO Y[i]:=Y[i] + Temp*A[j,i]; END;
END;
INC(jx,IncX);
END
ELSE
FOR j:=0 TO M-1 DO
IF (X[jx] # 0.0) THEN
Temp := Alpha*X[jx];
iy := ky-1;
FOR i:=0 TO N-1 DO
Y[iy]:=Y[iy] + Temp*A[j,i];
INC(iy,IncY);
END;
END;
INC(jx,IncX);
END;
END; (* ELSE *)
END; (* IF Trans *)
END dgemv;
PROCEDURE dgemm( TransA,TransB : CHAR;
M,N,K : INTEGER;
Alpha : FLOAT;
VAR A : ARRAY OF ARRAY OF FLOAT; (* [1..LDA][1..*] *)
LDA : INTEGER;
VAR B : ARRAY OF ARRAY OF FLOAT; (* [1..LDB][1..*] *)
LDB : INTEGER;
Beta : FLOAT;
VAR C : ARRAY OF ARRAY OF FLOAT; (* [1..LDC][1..*] *)
LDC : INTEGER);
VAR NotA,NotB : BOOLEAN;
i,j,k,Info,NRowA,NRowB : INTEGER;
Temp,Aki,Aik : FLOAT;
Nm1,Mm1,Km1 : INTEGER;
BEGIN
(* Set NotA and NotB as true if A and B respectively are not *)
(* transposed and set NRowA, NColA and NRowB as the number of rows *)
(* and columns of A and the number of rows of B respectively. *)
NotA := CAP(TransA) = 'N';
NotB := CAP(TransB) = 'N';
IF (NotA) THEN
NRowA := M; (* NColA := K; not used MR, 01.12.15 *)
ELSE
NRowA := K; (* NColA := M; not used *)
END;
IF (NotB) THEN
NRowB := K;
ELSE
NRowB := N;
END;
(* Test the input parameters. *)
Info := 0;
IF (NOT NotA) AND (CAP(TransA) # 'C') AND (CAP(TransA) # 'T') THEN
Info := 1;
ELSIF (NOT NotB) AND (CAP(TransB) # 'C') AND (CAP(TransB) # 'T') THEN
Info := 2;
ELSIF (M < 0) THEN
Info := 3;
ELSIF (N < 0) THEN
Info := 4;
ELSIF (K < 0) THEN
Info := 5;
ELSIF (LDA < IMax(1,NRowA)) THEN
Info := 8;
ELSIF (LDB < IMax(1,NRowB)) THEN
Info := 10;
ELSIF (LDC < IMax(1,M)) THEN
Info := 13
END;
IF (Info <> 0) THEN (* Fehlerbehandlung ! *)
Errors.WriteLn; Errors.WriteLn;
Errors.WriteString("Fehler in LibDBlas.dgemm, Info = ");
Errors.WriteInt(Info);
Errors.WriteLn; Errors.WriteLn;
RETURN;
END;
(* Quick return if possible. *)
IF (M = 0) OR (N = 0) OR
(((Alpha = 0.0) OR (K = 0)) AND (Beta = 1.0))
THEN
RETURN;
END;
IF (Beta = 0.0) THEN
FOR i:=0 TO M-1 DO
FOR j:=0 TO N-1 DO
C[i,j]:=0.0;
END;
END;
ELSE
FOR i:=0 TO M-1 DO
FOR j:=0 TO N-1 DO
C[i,j]:=C[i,j]*Beta;
END;
END;
END;
IF (Alpha = 0.0) THEN
RETURN; (* !!! *)
END;
(* Start the operations. *)
Mm1:=M-1; Nm1:=N-1; Km1:=K-1;
IF (NotB) THEN
IF (NotA) THEN (* Form C := alpha*A*B + beta*C. *)
FOR i:=0 TO Mm1 DO (* i,k,j *)
FOR k:=0 TO Km1 DO
IF (A[i,k] # 0.0) THEN
Aik := Alpha*A[i,k];
FOR j:=0 TO Nm1 DO
C[i,j]:=C[i,j] + Aik*B[k,j];
END;
END;
END; (* FOR k *)
END; (* FOR i *)
ELSE (* Form C := alpha*A'*B + beta*C *)
FOR k:=0 TO Km1 DO (* k,i,j *)
FOR i:=0 TO Mm1 DO
Aki:=Alpha*A[k,i];
FOR j:=0 TO Nm1 DO
C[i,j]:=C[i,j] + Aki*B[k,j];
END;
END; (* FOR i *)
END; (* FOR k *)
END; (* IF NotA *)
ELSE (* NOT NotB *)
IF (NotA) THEN (* Form C := alpha*A*B' + beta*C *)
FOR i:=0 TO Mm1 DO (* i,j,k *)
FOR j:=0 TO Nm1 DO
Temp:=0.0;
FOR k:=0 TO Km1 DO
Temp:=Temp + A[i,k]*B[j,k];
END;
IF (Beta = 0.0) THEN
C[i,j] := Alpha*Temp;
ELSE
C[i,j] := Alpha*Temp + C[i,j];
END;
END; (* FOR j *)
END; (* FOR i *)
ELSE (* Form C := alpha*A'*B' + beta*C *)
FOR i:=0 TO Mm1 DO (* i,j,k *)
FOR j:=0 TO Nm1 DO
Temp:=0.0;
FOR k:=0 TO Km1 DO
Temp:=Temp + A[k,i]*B[j,k];
END;
IF (Beta = 0.0) THEN
C[i,j] := Alpha*Temp;
ELSE
C[i,j] := Alpha*Temp + C[i,j];
END;
END; (* FOR j *)
END; (* FOR i *)
END; (* IF NotA *)
END; (* IF NotB *)
END dgemm;
PROCEDURE zgemm( TransA,TransB : CHAR;
M,N,K : INTEGER;
Alpha : LONGCOMPLEX;
VAR A : ARRAY OF ARRAY OF LONGCOMPLEX;
LDA : INTEGER;
VAR B : ARRAY OF ARRAY OF LONGCOMPLEX;
LDB : INTEGER;
Beta : LONGCOMPLEX;
VAR C : ARRAY OF ARRAY OF LONGCOMPLEX;
LDC : INTEGER);
CONST zero = CMPLX(0.0,0.0);
one = CMPLX(1.0,0.0);
VAR NotA,NotB,ConjA,ConjB : BOOLEAN;
i,j,k,Info,NRowA,NRowB : INTEGER;
Temp,Aki,Aik : LONGCOMPLEX;
Nm1,Mm1,Km1 : INTEGER;
BEGIN
(* Set NotA and NotB as true if A and B respectively are not *)
(* transposed and set NRowA, NColA and NRowB as the number of rows *)
(* and columns of A and the number of rows of B respectively. *)
NotA := (CAP(TransA) = 'N');
NotB := (CAP(TransB) = 'N');
ConjA := (CAP(TransA) = 'C');
ConjB := (CAP(TransB) = 'C');
IF (NotA) THEN
NRowA := M; (* NColA := K; not used MR, 01.12.15 *)
ELSE
NRowA := K; (* NColA := M; not used *)
END;
IF (NotB) THEN
NRowB := K;
ELSE
NRowB := N;
END;
(* Test the input parameters. *)
Info := 0;
IF ((NOT NotA) AND (NOT ConjA) AND (CAP(TransA) # 'T')) THEN
Info := 1;
ELSIF ((NOT NotB) AND (NOT ConjB) AND (CAP(TransB) # 'T')) THEN
Info := 2;
ELSIF (M < 0) THEN
Info := 3;
ELSIF ( N < 0) THEN
Info := 4;
ELSIF (K < 0) THEN
Info := 5;
ELSIF (LDA < IMax(1,NRowA)) THEN
Info := 8;
ELSIF (LDB < IMax(1,NRowB)) THEN
Info := 10;
ELSIF (LDC < IMax(1,M)) THEN
Info := 13;
END;
IF (Info <> 0) THEN (* Fehlerbehandlung ! *)
Errors.WriteLn; Errors.WriteLn;
Errors.WriteString("Fehler in LibDBlas.dgemm, Info = ");
Errors.WriteInt(Info);
Errors.WriteLn; Errors.WriteLn;
RETURN;
END;
(* Quick return if possible. *)
IF (M = 0) OR (N = 0) OR
(((Alpha = zero) OR (K = 0)) AND (Beta = one))
THEN
RETURN;
END;
IF (Beta = zero) THEN
FOR i:=0 TO M-1 DO
FOR j:=0 TO N-1 DO C[i,j]:=zero; END;
END;
ELSE
FOR i:=0 TO M-1 DO
FOR j:=0 TO N-1 DO C[i,j]:=C[i,j]*Beta; END;
END;
END; (* IF *)
IF (Alpha = zero) THEN
RETURN; (* !!! *)
END;
Nm1:=N-1; Mm1:=M-1; Km1:=K-1;
(* start the operations *)
IF (NotB) THEN
IF (NotA) THEN (* form C := alpha * A x B + beta*C *)
FOR i:=0 TO Mm1 DO (* i,k,j *)
FOR k:=0 TO Km1 DO
IF (A[i,k] # zero) THEN
Aik := Alpha*A[i,k];
FOR j:=0 TO Nm1 DO C[i,j]:=C[i,j] + Aik*B[k,j]; END;
END;
END; (* FOR k *)
END; (* FOR i *)
ELSIF (ConjA) THEN (* form C := alpha * Conj(A) x B + beta*C *)
FOR k:=0 TO Km1 DO (* k,i,j *)
FOR i:=0 TO Mm1 DO
Aki:=Alpha*conj(A[k,i]);
FOR j:=0 TO Nm1 DO C[i,j]:=C[i,j] + Aki*B[k,j]; END;
END; (* FOR i *)
END; (* FOR k *)
ELSE (* form C := alpha * A' x B + beta*C *)
FOR k:=0 TO Km1 DO (* k,i,j *)
FOR i:=0 TO Mm1 DO
Aki:=Alpha*A[k,i];
FOR j:=0 TO Nm1 DO C[i,j]:=C[i,j] + Aki*B[k,j]; END;
END; (* FOR i *)
END; (* FOR k *)
END; (* IF NotA | ConjA *)
ELSIF (ConjB) THEN
IF (NotA) THEN (* form C := alpha * A x Conj(B) + beta*C *)
FOR i:=0 TO Mm1 DO (* i,j,k *)
FOR j:=0 TO Nm1 DO
Temp:=zero;
FOR k:=0 TO Km1 DO Temp:=Temp + A[i,k]*conj(B[j,k]); END;
IF (Beta = zero) THEN
C[i,j] := Alpha*Temp;
ELSE
C[i,j] := Alpha*Temp + C[i,j];
END;
END;
END; (* FOR i *)
ELSIF (ConjA) THEN (* form C := alpha * Conj(A) x Conj(B) + beta*C *)
FOR i:=0 TO Mm1 DO (* i,j,k *)
FOR j:=0 TO Nm1 DO
Temp:=zero;
FOR k:=0 TO Km1 DO Temp:=Temp + conj(A[k,i])*conj(B[j,k]); END;
IF (Beta = zero) THEN
C[i,j] := Alpha*Temp;
ELSE
C[i,j] := Alpha*Temp + C[i,j];
END;
END; (* FOR j *)
END; (* FOR i *)
ELSE (* form C := alpha * A' x Conj(B) + beta*C *)
FOR i:=0 TO Mm1 DO (* i,j,k *)
FOR j:=0 TO Nm1 DO
Temp:=zero;
FOR k:=0 TO Km1 DO Temp:=Temp + A[k,i]*conj(B[j,k]); END;
IF (Beta = zero) THEN
C[i,j] := Alpha*Temp;
ELSE
C[i,j] := Alpha*Temp + C[i,j];
END;
END; (* FOR j *)
END; (* FOR i *)
END; (* IF NotA | ConjA *)
ELSE (* NOT NotB *)
IF (NotA) THEN (* form C := alpha * A x B' + beta*C *)
FOR i:=0 TO Mm1 DO (* i,j,k *)
FOR j:=0 TO Nm1 DO
Temp:=zero;
FOR k:=0 TO Km1 DO Temp:=Temp + A[i,k]*B[j,k]; END;
IF (Beta = zero) THEN
C[i,j] := Alpha*Temp;
ELSE
C[i,j] := Alpha*Temp + C[i,j];
END;
END; (* FOR j *)
END; (* FOR i *)
ELSIF (ConjA) THEN (* form C := alpha * Conj(A) x B' + beta*C *)
FOR i:=0 TO Mm1 DO (* i,j,k *)
FOR j:=0 TO Nm1 DO
Temp:=zero;
FOR k:=0 TO Km1 DO Temp:=Temp + conj(A[k,i])*B[j,k]; END;
IF (Beta = zero) THEN
C[i,j] := Alpha*Temp;
ELSE
C[i,j] := Alpha*Temp + C[i,j];
END;
END; (* FOR j *)
END; (* FOR i *)
ELSE (* Form C := alpha * A' x B' + beta*C *)
FOR i:=0 TO Mm1 DO (* i,j,k *)
FOR j:=0 TO Nm1 DO
Temp:=zero;
FOR k:=0 TO Km1 DO Temp:=Temp + A[k,i]*B[j,k]; END;
IF (Beta = zero) THEN
C[i,j] := Alpha*Temp;
ELSE
C[i,j] := Alpha*Temp + C[i,j];
END;
END; (* FOR j *)
END; (* FOR i *)
END; (* IF NotA | ConjA *)
END; (* IF NotB | ConjB *)
END zgemm;
PROCEDURE dger( m,n : CARDINAL;
Alpha : FLOAT;
VAR X : ARRAY OF FLOAT;
IncX : CARDINAL;
VAR Y : ARRAY OF FLOAT;
IncY : CARDINAL;
VAR A : ARRAY OF ARRAY OF FLOAT;
lda : CARDINAL);
VAR i,j,ix,jy : CARDINAL;
AlphaX : FLOAT;
BEGIN
IF (lda-1 > HIGH(A)) THEN
Errors.Fehlerflag:="lda-1 > HIGH(A) (LibDBlas.dger) !";
Errors.Fehler:=TRUE;
RETURN;
END;
IF (m = 0) OR (n = 0) OR (Alpha = 0.0) THEN RETURN; END;
IF (IncX = 1) AND (IncY = 1) THEN
FOR i:=0 TO m-1 DO
IF (X[i] # 0.0) THEN
AlphaX := Alpha*X[i];
FOR j:=0 TO n-1 DO
A[i,j]:=A[i,j] + AlphaX*Y[j];
END;
END;
END;
ELSE
ix := 0;
FOR i:=0 TO m-1 DO
IF (X[i] # 0.0) THEN
AlphaX := Alpha*X[ix];
jy := 0;
FOR j:=0 TO n-1 DO
A[i,j]:=A[i,j] + AlphaX*Y[jy];
INC(jy,IncY);
END;
END;
INC(ix,IncX);
END;
END;
END dger;
END LibDBlasM2.