IMPLEMENTATION MODULE LibDBlas;
(*------------------------------------------------------------------------*)
(* Basic Linear Algebra Prozeduren in Modula-2 wobei die Parameter *)
(* konform zur den Fortran-Routinen uebergeben werden. *)
(* Some level 1 BLAS routines implemented in Modula-2 pathing the *)
(* arguments in a "Fortran 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 *)
(* 05.04.16, MRi: Umstellung der Uerbergabeparameter fuer Vektoren *)
(* auf die Startadresse im jeweiligen Feld um bessere *)
(* Portierbarkeit von Fortran Quellen zu erm��glichen. *)
(* 15.04.16, MRi: dcopy aus LibDBlasM2 uebernommen und angepasst *)
(* 10.05.16, MRi: dger aus PMatLib uebernommen und angepasst *)
(* 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 *)
(* 29.10.17, MRi: dgemv,dgemm,zgemm und dger sind nur noch Verweise auf *)
(* die entsprechenden Routinen in LibDBlasM2 *)
(*------------------------------------------------------------------------*)
(* Testroutinen *)
(* *)
(* - zgemm in TstCmplxMaMul *)
(*------------------------------------------------------------------------*)
(* Offene Punkte: *)
(* *)
(* - Weiteres austesten von dgemm und besonders dgemv *)
(* - Austesten von dger *)
(*------------------------------------------------------------------------*)
(* Implementation : Michael Riedl *)
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
(* $Id: LibDBlas.mod,v 1.4 2017/10/29 09:54:58 mriedl Exp mriedl $ *)
IMPORT SYSTEM;
FROM LongMath IMPORT sqrt;
IMPORT Errors;
FROM Deklera IMPORT FLOAT; (* REAL type *)
TYPE PVEKTOR = POINTER TO ARRAY [0..MAX(INTEGER)-1] OF FLOAT;
(*========================================================================*)
(* 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;
XX : PVEKTOR;
BEGIN
XX:=SYSTEM.CAST(PVEKTOR,SYSTEM.ADR(X));
IF (n < 1) OR (IncX < 1) THEN
norm := 0.0
ELSIF (n = 1) THEN
norm := ABS(XX^[0]);
ELSE
scale := 0.0;
ssq := 1.0;
ix := 0;
FOR i:=0 TO (n-1) DO
IF (XX^[ix] # 0.0) THEN
absxi := ABS(XX^[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;
XX,YY : PVEKTOR;
BEGIN
IF (N = 0) THEN RETURN; END;
XX:=SYSTEM.CAST(PVEKTOR,SYSTEM.ADR(X));
YY:=SYSTEM.CAST(PVEKTOR,SYSTEM.ADR(Y));
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 := XX^[i]; XX^[i] := YY^[i]; YY^[i] := Xi;
END;
END;
IF (N < level) THEN RETURN; END;
FOR i:=m TO N-1 BY level DO
Xtmp[0] := XX^[i+0]; XX^[i+0] := YY^[i+0]; YY^[i+0] := Xtmp[0];
Xtmp[1] := XX^[i+1]; XX^[i+1] := YY^[i+1]; YY^[i+1] := Xtmp[1];
Xtmp[2] := XX^[i+2]; XX^[i+2] := YY^[i+2]; YY^[i+2] := Xtmp[2];
Xtmp[3] := XX^[i+3]; XX^[i+3] := YY^[i+3]; YY^[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 := XX^[ix]; XX^[ix] := YY^[iy]; YY^[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;
XX,YY : PVEKTOR;
BEGIN
IF (N <= 0) THEN RETURN; END;
XX:=SYSTEM.CAST(PVEKTOR,SYSTEM.ADR(X));
YY:=SYSTEM.CAST(PVEKTOR,SYSTEM.ADR(Y));
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 YY^[i] := XX^[i]; END;
IF (N < 8) THEN RETURN; END
END;
FOR i:=m TO N-1 BY 8 DO
YY^[i+0] := XX^[i+0]; YY^[i+1] := XX^[i+1];
YY^[i+2] := XX^[i+2]; YY^[i+3] := XX^[i+3];
YY^[i+4] := XX^[i+4]; YY^[i+5] := XX^[i+5];
YY^[i+6] := XX^[i+6]; YY^[i+7] := XX^[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
YY^[iy] := XX^[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;
XX,YY : PVEKTOR;
BEGIN
XX:=SYSTEM.CAST(PVEKTOR,SYSTEM.ADR(X));
YY:=SYSTEM.CAST(PVEKTOR,SYSTEM.ADR(Y));
IF (IncX = 1) AND (IncY = 1) THEN
FOR i:=0 TO N-1 DO
x1 := XX^[i];
y1 := YY^[i];
XX^[i] := c * x1 + s * y1;
YY^[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 := XX^[ix];
y1 := YY^[iy];
XX^[ix] := c * x1 + s * y1;
YY^[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 : FLOAT;
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;
c := da / r;
s := db / r;
z := 1.0;
IF (ABS(da) > ABS(db)) THEN
z := s;
END;
IF (ABS(db) >= ABS(da)) AND (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;
XX : PVEKTOR;
BEGIN
IF (n <= 0) THEN RETURN; END;
XX:=SYSTEM.CAST(PVEKTOR,SYSTEM.ADR(dx));
IF (IncX <> 1) THEN (* code for increment not equal to 1 *)
FOR i:=0 TO n-1 DO
XX^[i*IncX] := da*XX^[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 XX^[i] := da*XX^[i]; END;
IF (n < 5) THEN
RETURN;
END;
END;
FOR i:=m TO n-1 BY 5 DO
XX^[i] := da*XX^[i];
XX^[i+1] := da*XX^[i+1];
XX^[i+2] := da*XX^[i+2];
XX^[i+3] := da*XX^[i+3];
XX^[i+4] := da*XX^[i+4];
END;
END;
END dscal;
PROCEDURE daxpy( n : INTEGER;
da : FLOAT;
VAR X : (* ARRAY OF *) FLOAT;
IncX : INTEGER;
VAR Y : (* ARRAY OF *) FLOAT;
IncY : INTEGER);
(*----------------------------------------------------------------*)
(* constant times a vector plus a vector *)
(*----------------------------------------------------------------*)
VAR i,ix,iy,m : INTEGER;
XX,YY : PVEKTOR;
BEGIN
IF (n <= 0 ) THEN RETURN; END;
IF (da = 0.0) THEN RETURN; END;
XX:=SYSTEM.CAST(PVEKTOR,SYSTEM.ADR(X));
YY:=SYSTEM.CAST(PVEKTOR,SYSTEM.ADR(Y));
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
YY^[iy] := YY^[iy] + da*XX^[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
YY^[i] := YY^[i] + da*XX^[i];
END;
IF (n < 4) THEN RETURN; END;
END;
FOR i:=m TO n-1 BY 4 DO
YY^[i] := YY^[i] + da*XX^[i];
YY^[i+1] := YY^[i+1] + da*XX^[i+1];
YY^[i+2] := YY^[i+2] + da*XX^[i+2];
YY^[i+3] := YY^[i+3] + da*XX^[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;
XX,YY : PVEKTOR;
BEGIN
IF (N <= 0) THEN RETURN 0.0; END;
XX:=SYSTEM.CAST(PVEKTOR,SYSTEM.ADR(X));
YY:=SYSTEM.CAST(PVEKTOR,SYSTEM.ADR(Y));
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 + XX^[i]*YY^[i];
END;
IF (N < veclen) THEN RETURN dtemp; END;
END;
(* i := m - veclen; *)
FOR i:=m TO N-1 BY veclen DO
dtemp:=dtemp + XX^[i+0]*YY^[i+0] + XX^[i+1]*YY^[i+1] +
XX^[i+2]*YY^[i+2] + XX^[i+3]*YY^[i+3] +
XX^[i+4]*YY^[i+4] + XX^[i+5]*YY^[i+5] +
XX^[i+6]*YY^[i+6] + XX^[i+7]*YY^[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 + XX^[ix]*YY^[iy];
INC(ix,IncX); INC(iy,IncY);
END;
END;
RETURN dtemp;
END ddot;
PROCEDURE idamax( n : INTEGER;
VAR X : (* ARRAY OF *) FLOAT;
IncX : INTEGER) : INTEGER;
(*----------------------------------------------------------------*)
(* Finds the index of element having max. absolute value. *)
(*----------------------------------------------------------------*)
VAR dmax : 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 *)
dmax := ABS(XX^[0]);
ix := 1 + IncX;
FOR i:=1 TO n-1 DO
IF (ABS(XX^[ix]) > dmax) THEN
itemp := i;
dmax := ABS(XX^[ix]);
END;
INC(ix,IncX);
END;
ELSE (* code for increment equal to 1 *)
dmax := ABS(XX^[0]);
FOR i:=1 TO n-1 DO
IF (ABS(XX^[i]) > dmax) THEN
itemp := i;
dmax := ABS(XX^[i]);
END;
END;
END;
RETURN itemp + 1;
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;
IncX : CARDINAL) : FLOAT;
CONST max = 6; (* 8 *)
VAR i,ix,m : CARDINAL;
Sum : FLOAT;
XX : PVEKTOR;
BEGIN
IF (dim = 0) THEN RETURN 0.0; END;
IF (IncX < 1) THEN
Errors.FatalError("Inkrement kleiner 1 (LibDBlas.dasum) !");
END;
XX:=SYSTEM.CAST(PVEKTOR,SYSTEM.ADR(X));
Sum:=0.0;
IF (IncX = 1) THEN
(* n:=dim DIV max; *)
m:=dim DIV max;
IF (m > 0) THEN
FOR i:=0 TO m-1 DO Sum:=Sum + ABS(XX^[i]); END;
END;
FOR i:=m TO dim-1 BY max DO
Sum:=Sum + ABS(XX^[i ]) + ABS(XX^[i+1])
+ ABS(XX^[i+2]) + ABS(XX^[i+3])
+ ABS(XX^[i+4]) + ABS(XX^[i+5])
(* + ABS(XX^[i+6]) + ABS(XX^[i+7])*);
END;
ELSE
ix:=0;
FOR i:=0 TO dim-1 DO Sum:=Sum + ABS(XX^[ix]); INC(ix,IncX); END;
END;
RETURN Sum;
END dasum;
END LibDBlas.