--- a/LibDBlasM2.mod
+++ b/LibDBlasM2.mod
@@ -11,7 +11,7 @@
(* 23.11.93, MRi: Erstellen der 1. Version *)
(* Aug. 95, MRi: Erweiterung *)
(* 03.11.95, MRi: Durchsicht *)
- (* 13.10.15, MRi: Ersetzen von FLOAT(#) durch VAL(FLOAT,#) *)
+ (* 13.10.15, MRi: Ersetzen von REAL8(#) durch VAL(REAL8,#) *)
(* 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 *)
@@ -28,6 +28,7 @@
(* 29.10.17, MRi: Rolle von M,N in dgemm und dgemv an zgemm angepasst *)
(* 11.09.18, MRi: zgemv eingefuegt *)
(* 12.09.18, MRi: zswap,zcopy,zdotc,dznrm2,zscal,zaxpy,zdrot eingefuegt *)
+ (* 01.02.19, MRi: Abfragen in idamx,idamin korrigiert *)
(*------------------------------------------------------------------------*)
(* Testroutinen *)
(* *)
@@ -47,11 +48,10 @@
(* $Id: LibDBlasM2.mod,v 1.11 2018/09/12 13:20:49 mriedl Exp mriedl $ *)
IMPORT SYSTEM;
-FROM Deklera IMPORT FLOAT,CFLOAT; (* REAL/COMPLEX Type *)
FROM LongMath IMPORT sqrt;
FROM LongComplexMath IMPORT zero,one,conj,scalarMult;
IMPORT Errors;
-FROM F77func IMPORT MAX0;
+FROM F77func IMPORT REAL8,COMPLEX16,MAX0;
PROCEDURE IMax(a,b : INTEGER) : INTEGER;
@@ -60,11 +60,11 @@
IF (a > b) THEN RETURN a; ELSE RETURN b; END;
END IMax;
-PROCEDURE SumVek(VAR X : ARRAY OF FLOAT;
- a,e : CARDINAL) : FLOAT;
+PROCEDURE SumVek(VAR X : ARRAY OF REAL8;
+ a,e : CARDINAL) : REAL8;
VAR i : CARDINAL;
- s,c,y,t : FLOAT;
+ s,c,y,t : REAL8;
BEGIN
s := X[a-1];
c := 0.0;
@@ -77,11 +77,11 @@
RETURN s;
END SumVek;
-PROCEDURE AbsSumVek(VAR X : ARRAY OF FLOAT;
- a,e : CARDINAL) : FLOAT;
+PROCEDURE AbsSumVek(VAR X : ARRAY OF REAL8;
+ a,e : CARDINAL) : REAL8;
VAR i : CARDINAL;
- s,c,y,t : FLOAT;
+ s,c,y,t : REAL8;
BEGIN
s := ABS(X[a-1]);
@@ -95,8 +95,8 @@
RETURN s;
END AbsSumVek;
-PROCEDURE ENorm(VAR X : ARRAY OF FLOAT;
- a,e : CARDINAL) : FLOAT;
+PROCEDURE ENorm(VAR X : ARRAY OF REAL8;
+ a,e : CARDINAL) : REAL8;
(*---------------------------------------------------------------*)
(* Argonne national laboratory. MINPACK project. March 1980. *)
@@ -109,13 +109,13 @@
MaxLR = 1.304E+19;
VAR i : CARDINAL;
- S1,S2,S3,zw : FLOAT;
- Giant,AbsX,Max1,Max3 : FLOAT;
- Norm : FLOAT;
+ S1,S2,S3,zw : REAL8;
+ Giant,AbsX,Max1,Max3 : REAL8;
+ Norm : REAL8;
BEGIN
S1:=0.0; S2:=0.0; S3:=0.0;
Max1:=0.0; Max3:=0.0;
- Giant := MaxLR / VAL(FLOAT,e-a+1);
+ Giant := MaxLR / VAL(REAL8,e-a+1);
FOR i:=a-1 TO e-1 DO
AbsX := ABS(X[i]);
IF (AbsX < MinLR) THEN (* Summiere kleine Elemente. *)
@@ -163,8 +163,8 @@
(*========================================================================*)
PROCEDURE dnrm2( n : INTEGER;
- VAR X : ARRAY OF FLOAT;
- IncX : INTEGER): FLOAT;
+ VAR X : ARRAY OF REAL8;
+ IncX : INTEGER): REAL8;
(*-----------------------------------------------------------------*)
(* This version written on 25-October-1982. *)
@@ -174,8 +174,8 @@
(*-----------------------------------------------------------------*)
VAR i,ix : INTEGER;
- absxi,norm,scale,ssq : FLOAT;
- zw : FLOAT;
+ absxi,norm,scale,ssq : REAL8;
+ zw : REAL8;
BEGIN
IF (n < 1) OR (IncX < 1) THEN
norm := 0.0
@@ -205,9 +205,9 @@
END dnrm2;
PROCEDURE dswap( N : CARDINAL;
- VAR X : ARRAY OF FLOAT;
+ VAR X : ARRAY OF REAL8;
incX : INTEGER;
- VAR Y : ARRAY OF FLOAT;
+ VAR Y : ARRAY OF REAL8;
incY : INTEGER);
(*----------------------------------------------------------------*)
@@ -216,8 +216,8 @@
CONST level = 4;
- VAR Xi : FLOAT;
- Xtmp : ARRAY [0..level-1] OF FLOAT;
+ VAR Xi : REAL8;
+ Xtmp : ARRAY [0..level-1] OF REAL8;
ix,iy : INTEGER;
i,m : CARDINAL;
BEGIN
@@ -248,9 +248,9 @@
END dswap;
PROCEDURE dcopy( N : INTEGER;
- VAR X : ARRAY OF FLOAT;
+ VAR X : ARRAY OF REAL8;
IncX : INTEGER;
- VAR Y : ARRAY OF FLOAT;
+ VAR Y : ARRAY OF REAL8;
IncY : INTEGER);
(*----------------------------------------------------------------*)
@@ -286,18 +286,18 @@
END dcopy;
PROCEDURE drot( N : INTEGER;
- VAR X : ARRAY OF FLOAT;
+ VAR X : ARRAY OF REAL8;
incX : INTEGER;
- VAR Y : ARRAY OF FLOAT;
+ VAR Y : ARRAY OF REAL8;
incY : INTEGER;
- c,s : FLOAT);
+ c,s : REAL8);
(*----------------------------------------------------------------*)
(* Adopted to Modula-2, MRi, 31.03.2016 *)
(*----------------------------------------------------------------*)
VAR i,ix,iy : INTEGER;
- x1,y1 : FLOAT;
+ x1,y1 : REAL8;
BEGIN
IF (incX = 1) AND (incY = 1) THEN
FOR i:=0 TO N-1 DO
@@ -320,16 +320,16 @@
END;
END drot;
-PROCEDURE drotg(VAR da : FLOAT;
- VAR db : FLOAT;
- VAR c : FLOAT;
- VAR s : FLOAT);
+PROCEDURE drotg(VAR da : REAL8;
+ VAR db : REAL8;
+ VAR c : REAL8;
+ VAR s : REAL8);
(*----------------------------------------------------------------*)
(* Adopted to Modula-2, MRi, 31.03.2016 *)
(*----------------------------------------------------------------*)
- VAR r,z,roe,daos,dbos,scale : FLOAT;
+ VAR r,z,roe,daos,dbos,scale : REAL8;
BEGIN
IF (ABS(da) > ABS(db)) THEN roe := da; ELSE roe := db; END;
scale := ABS(da) + ABS(db);
@@ -358,8 +358,8 @@
END drotg;
PROCEDURE dscal( n : INTEGER;
- da : FLOAT;
- VAR dx : ARRAY OF FLOAT;
+ da : REAL8;
+ VAR dx : ARRAY OF REAL8;
incx : INTEGER);
(*----------------------------------------------------------------*)
@@ -394,10 +394,10 @@
END dscal;
PROCEDURE daxpy( n : INTEGER;
- da : FLOAT;
- VAR dx : ARRAY OF FLOAT;
+ da : REAL8;
+ VAR dx : ARRAY OF REAL8;
incx : INTEGER;
- VAR dy : ARRAY OF FLOAT;
+ VAR dy : ARRAY OF REAL8;
incy : INTEGER);
(*----------------------------------------------------------------*)
@@ -438,10 +438,10 @@
END daxpy;
PROCEDURE ddot( N : INTEGER;
- VAR X : ARRAY OF FLOAT;
+ VAR X : ARRAY OF REAL8;
IncX : INTEGER;
- VAR Y : ARRAY OF FLOAT;
- IncY : INTEGER) : FLOAT;
+ VAR Y : ARRAY OF REAL8;
+ IncY : INTEGER) : REAL8;
(*----------------------------------------------------------------*)
(* Forms the dot product of two vectors. Uses unrolled loops for *)
@@ -451,7 +451,7 @@
CONST veclen = 8;
- VAR dtemp : FLOAT;
+ VAR dtemp : REAL8;
i,ix,iy,m : INTEGER;
BEGIN
IF (N <= 0) THEN RETURN 0.0; END;
@@ -485,19 +485,19 @@
END ddot;
PROCEDURE idamax( n : INTEGER;
- VAR dx : ARRAY OF FLOAT;
+ VAR dx : ARRAY OF REAL8;
incx : INTEGER) : INTEGER;
(*----------------------------------------------------------------*)
(* Finds the index of element having max. absolute value. *)
(*----------------------------------------------------------------*)
- VAR dmax : FLOAT;
+ VAR dmax : REAL8;
i,ix,itemp : INTEGER;
BEGIN
- IF (n < 1) THEN RETURN -1; END;
- IF (n = 1) THEN RETURN 0; END;
+ IF (n < 1) THEN RETURN 0; END;
+ IF (n = 1) THEN RETURN 1; END;
itemp:=0;
IF (incx <> 1) THEN (* code for increment not equal to 1 *)
dmax := ABS(dx[0]);
@@ -522,18 +522,18 @@
END idamax;
PROCEDURE idamin( n : INTEGER;
- VAR X : ARRAY OF FLOAT;
+ VAR X : ARRAY OF REAL8;
IncX : INTEGER) : INTEGER;
(*----------------------------------------------------------------*)
(* Finds the index of element having min. absolute value. *)
(*----------------------------------------------------------------*)
- VAR dmin : FLOAT;
+ VAR dmin : REAL8;
i,ix,itemp : INTEGER;
BEGIN
- IF (n < 1) THEN RETURN -1; END;
- IF (n = 1) THEN RETURN 0; END;
+ IF (n < 1) THEN RETURN 0; END;
+ IF (n = 1) THEN RETURN 1; END;
itemp:=0;
IF (IncX <> 1) THEN (* code for Increment not equal to 1 *)
@@ -559,14 +559,14 @@
END idamin;
PROCEDURE dasum( dim : CARDINAL;
- VAR X : ARRAY OF FLOAT;
+ VAR X : ARRAY OF REAL8;
Inc : CARDINAL; (* Inkrementwert >= 1*)
- Unrolled : BOOLEAN) : FLOAT;
+ Unrolled : BOOLEAN) : REAL8;
CONST max = 6; (* 8 *)
VAR i,m : CARDINAL;
- Sum : FLOAT;
+ Sum : REAL8;
BEGIN
IF (dim = 0) THEN RETURN 0.0; END;
IF (Inc < 1) THEN
@@ -602,16 +602,16 @@
PROCEDURE dgemv ( Trans : CHAR;
M,N : INTEGER;
- Alpha : FLOAT;
- VAR A : ARRAY OF ARRAY OF FLOAT; (* [1..LDA][1..*] *)
+ Alpha : REAL8;
+ VAR A : ARRAY OF ARRAY OF REAL8; (* [1..LDA][1..*] *)
LDA : INTEGER;
- VAR X : ARRAY OF FLOAT; (* [1..*] *)
+ VAR X : ARRAY OF REAL8; (* [1..*] *)
IncX : INTEGER;
- Beta : FLOAT;
- VAR Y : ARRAY OF FLOAT; (* [1..*] *)
+ Beta : REAL8;
+ VAR Y : ARRAY OF REAL8; (* [1..*] *)
IncY : INTEGER);
- VAR Temp : FLOAT;
+ VAR Temp : REAL8;
i,j,iy,jx,jy,kx,ky,Info,LenX,LenY : INTEGER;
BEGIN (* Testen *)
(* Test the input parameters. *)
@@ -742,18 +742,18 @@
PROCEDURE dgemm( TransA,TransB : CHAR;
M,N,K : INTEGER;
- Alpha : FLOAT;
- VAR A : ARRAY OF ARRAY OF FLOAT; (* [1..LDA][1..*] *)
+ Alpha : REAL8;
+ VAR A : ARRAY OF ARRAY OF REAL8; (* [1..LDA][1..*] *)
LDA : INTEGER;
- VAR B : ARRAY OF ARRAY OF FLOAT; (* [1..LDB][1..*] *)
+ VAR B : ARRAY OF ARRAY OF REAL8; (* [1..LDB][1..*] *)
LDB : INTEGER;
- Beta : FLOAT;
- VAR C : ARRAY OF ARRAY OF FLOAT; (* [1..LDC][1..*] *)
+ Beta : REAL8;
+ VAR C : ARRAY OF ARRAY OF REAL8; (* [1..LDC][1..*] *)
LDC : INTEGER);
VAR NotA,NotB : BOOLEAN;
i,j,k,Info,NRowA,NRowB : INTEGER;
- Temp,Aki,Aik : FLOAT;
+ Temp,Aki,Aik : REAL8;
Nm1,Mm1,Km1 : INTEGER;
BEGIN
(* Set NotA and NotB as true if A and B respectively are not *)
@@ -899,16 +899,16 @@
END dgemm;
PROCEDURE dger( m,n : CARDINAL;
- Alpha : FLOAT;
- VAR X : ARRAY OF FLOAT;
+ Alpha : REAL8;
+ VAR X : ARRAY OF REAL8;
IncX : CARDINAL;
- VAR Y : ARRAY OF FLOAT;
+ VAR Y : ARRAY OF REAL8;
IncY : CARDINAL;
- VAR A : ARRAY OF ARRAY OF FLOAT;
+ VAR A : ARRAY OF ARRAY OF REAL8;
lda : CARDINAL);
VAR i,j,ix,jy : CARDINAL;
- AlphaX : FLOAT;
+ AlphaX : REAL8;
BEGIN
IF (lda-1 > HIGH(A)) THEN
Errors.Fehlerflag:="lda-1 > HIGH(A) (LibDBlas.dger) !";
@@ -945,15 +945,15 @@
(*=========================== Complex valued routines =====================*)
PROCEDURE zswap( N : CARDINAL;
- VAR X : ARRAY OF CFLOAT;
+ VAR X : ARRAY OF COMPLEX16;
IncX : INTEGER;
- VAR Y : ARRAY OF CFLOAT;
+ VAR Y : ARRAY OF COMPLEX16;
IncY : INTEGER);
CONST veclen = 4;
- VAR Xi : CFLOAT;
- Xtmp : ARRAY [0..veclen-1] OF CFLOAT;
+ VAR Xi : COMPLEX16;
+ Xtmp : ARRAY [0..veclen-1] OF COMPLEX16;
ix,iy : INTEGER;
i,m : CARDINAL;
BEGIN
@@ -986,9 +986,9 @@
END zswap;
PROCEDURE zcopy( N : INTEGER;
- VAR X : ARRAY OF CFLOAT;
+ VAR X : ARRAY OF COMPLEX16;
IncX : INTEGER;
- VAR Y : ARRAY OF CFLOAT;
+ VAR Y : ARRAY OF COMPLEX16;
IncY : INTEGER);
(*----------------------------------------------------------------*)
@@ -1024,14 +1024,14 @@
END zcopy;
PROCEDURE zdotc( N : INTEGER;
- VAR X : ARRAY OF CFLOAT;
+ VAR X : ARRAY OF COMPLEX16;
IncX : INTEGER;
- VAR Y : ARRAY OF CFLOAT;
- IncY : INTEGER) : CFLOAT;
+ VAR Y : ARRAY OF COMPLEX16;
+ IncY : INTEGER) : COMPLEX16;
CONST veclen = 4;
- VAR dtemp : CFLOAT;
+ VAR dtemp : COMPLEX16;
i,ix,iy,m : INTEGER;
BEGIN
IF (N <= 0) THEN RETURN zero; END;
@@ -1044,7 +1044,10 @@
FOR i:=0 TO m-1 DO (* clean-up loop *)
dtemp:=dtemp + conj(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
@@ -1064,12 +1067,12 @@
END zdotc;
PROCEDURE dznrm2( N : INTEGER;
- VAR X : ARRAY OF CFLOAT;
- IncX : INTEGER) : FLOAT;
-
- VAR norm,scale,ssq,tmp : FLOAT;
+ VAR X : ARRAY OF COMPLEX16;
+ IncX : INTEGER) : REAL8;
+
+ VAR norm,scale,ssq,tmp : REAL8;
i,ix : INTEGER;
- zw : FLOAT;
+ zw : REAL8;
BEGIN
IF (N < 1) OR (IncX < 1) THEN
norm := 0.0;
@@ -1108,8 +1111,8 @@
END dznrm2;
PROCEDURE zscal( n : INTEGER;
- da : CFLOAT;
- VAR X : ARRAY OF CFLOAT;
+ da : COMPLEX16;
+ VAR X : ARRAY OF COMPLEX16;
IncX : INTEGER);
CONST veclen = 4;
@@ -1140,10 +1143,10 @@
END zscal;
PROCEDURE zaxpy( n : INTEGER;
- da : CFLOAT;
- VAR X : ARRAY OF CFLOAT;
+ da : COMPLEX16;
+ VAR X : ARRAY OF COMPLEX16;
IncX : INTEGER;
- VAR Y : ARRAY OF CFLOAT;
+ VAR Y : ARRAY OF COMPLEX16;
IncY : INTEGER);
CONST veclen = 4;
@@ -1182,14 +1185,14 @@
END zaxpy;
PROCEDURE zdrot( N : INTEGER;
- VAR X : ARRAY OF CFLOAT;
+ VAR X : ARRAY OF COMPLEX16;
IncX : INTEGER;
- VAR Y : ARRAY OF CFLOAT;
+ VAR Y : ARRAY OF COMPLEX16;
IncY : INTEGER;
- c,s : FLOAT);
+ c,s : REAL8);
VAR i,ix,iy : INTEGER;
- tmp : CFLOAT;
+ tmp : COMPLEX16;
BEGIN
IF (IncX = 1) AND (IncY = 1) THEN
FOR i:=0 TO N-1 DO
@@ -1212,18 +1215,18 @@
PROCEDURE zgemv( Trans : CHAR;
M,N : INTEGER;
- Alpha : CFLOAT;
- VAR A : ARRAY OF ARRAY OF CFLOAT;
+ Alpha : COMPLEX16;
+ VAR A : ARRAY OF ARRAY OF COMPLEX16;
lda : INTEGER;
- VAR X : ARRAY OF CFLOAT;
+ VAR X : ARRAY OF COMPLEX16;
IncX : INTEGER;
- Beta : CFLOAT;
- VAR Y : ARRAY OF CFLOAT;
+ Beta : COMPLEX16;
+ VAR Y : ARRAY OF COMPLEX16;
IncY : INTEGER);
VAR i,j,iy,jx,kx,ky : INTEGER;
lenx,leny,info : INTEGER;
- s : CFLOAT;
+ s : COMPLEX16;
noconj : BOOLEAN;
BEGIN
(* test the input parameters *)
@@ -1377,13 +1380,13 @@
PROCEDURE zgemm( TransA,TransB : CHAR;
M,N,K : INTEGER;
- Alpha : CFLOAT;
- VAR A : ARRAY OF ARRAY OF CFLOAT;
+ Alpha : COMPLEX16;
+ VAR A : ARRAY OF ARRAY OF COMPLEX16;
LDA : INTEGER;
- VAR B : ARRAY OF ARRAY OF CFLOAT;
+ VAR B : ARRAY OF ARRAY OF COMPLEX16;
LDB : INTEGER;
- Beta : CFLOAT;
- VAR C : ARRAY OF ARRAY OF CFLOAT;
+ Beta : COMPLEX16;
+ VAR C : ARRAY OF ARRAY OF COMPLEX16;
LDC : INTEGER);
CONST zero = CMPLX(0.0,0.0);
@@ -1391,7 +1394,7 @@
VAR NotA,NotB,ConjA,ConjB : BOOLEAN;
i,j,k,Info,NRowA,NRowB : INTEGER;
- Temp,Aki,Aik : CFLOAT;
+ Temp,Aki,Aik : COMPLEX16;
Nm1,Mm1,Km1 : INTEGER;
BEGIN
(* Set NotA and NotB as true if A and B respectively are not *)