--- a/LibDBlas.mod
+++ b/LibDBlas.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 LFLOAT(#) durch VAL(LONGREAL,#) *)
+ (* 13.10.15, MRi: Ersetzen von LREAL8(#) 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 *)
@@ -37,6 +37,8 @@
(* 21.06.18, MRi: Korrekturen in dznrm2 und zdotc *)
(* 11.09.18, MRi: Hinzufuegen von zcopy und zgemv (definition) *)
(* 30.09.18, MRi: Hinzufuegen von zdotu *)
+ (* 04.01.19, MRi: Korrektur in zaxpy bei Abfrage auf 0 von za *)
+ (* 01.02.19, MRi: Abfragen in idamx,idamin korrigiert *)
(*------------------------------------------------------------------------*)
(* Testroutinen *)
(* *)
@@ -53,18 +55,18 @@
(* $Id: LibDBlas.mod,v 1.6 2018/09/12 13:20:49 mriedl Exp mriedl $ *)
- IMPORT SYSTEM;
-FROM Deklera IMPORT FLOAT,CFLOAT; (* REAL,COMPLEX type *)
- IMPORT Errors;
-FROM LongMath IMPORT sqrt;
- IMPORT LongComplexMath;
+ IMPORT SYSTEM;
+ IMPORT Errors;
+FROM LongMath IMPORT sqrt;
+ IMPORT LongComplexMath;
+FROM LibDBlasL1F77 IMPORT REAL8,COMPLEX16;
CONST conj = LongComplexMath.conj;
scalarMult = LongComplexMath.scalarMult;
zero = LongComplexMath.zero;
-TYPE PVEKTOR = POINTER TO ARRAY [0..MAX(INTEGER)-1] OF FLOAT;
+TYPE PVEKTOR = POINTER TO ARRAY [0..MAX(INTEGER)-1] OF REAL8;
PZVEKTOR = POINTER TO ARRAY [0..MAX(INTEGER)-1] OF LONGCOMPLEX;
(*========================================================================*)
@@ -72,8 +74,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. *)
@@ -83,8 +85,8 @@
(*-----------------------------------------------------------------*)
VAR i,ix : INTEGER;
- absxi,norm,scale,ssq : FLOAT;
- zw : FLOAT;
+ absxi,norm,scale,ssq : REAL8;
+ zw : REAL8;
XX : PVEKTOR;
BEGIN
XX:=SYSTEM.CAST(PVEKTOR,SYSTEM.ADR(X));
@@ -116,9 +118,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);
(*----------------------------------------------------------------*)
@@ -127,8 +129,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;
XX,YY : PVEKTOR;
@@ -163,9 +165,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);
(*----------------------------------------------------------------*)
@@ -205,17 +207,17 @@
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;
XX,YY : PVEKTOR;
BEGIN
XX:=SYSTEM.CAST(PVEKTOR,SYSTEM.ADR(X));
@@ -242,17 +244,17 @@
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 : FLOAT;
- roe,daos,dbos,scale : FLOAT;
+ VAR r,z : REAL8;
+ roe,daos,dbos,scale : REAL8;
BEGIN
IF (ABS(da) > ABS(db)) THEN
roe := da;
@@ -285,8 +287,8 @@
END drotg;
PROCEDURE dscal( n : INTEGER;
- da : FLOAT;
- VAR dx : (* ARRAY OF *) FLOAT;
+ da : REAL8;
+ VAR dx : (* ARRAY OF *) REAL8;
IncX : INTEGER);
(*----------------------------------------------------------------*)
@@ -323,10 +325,10 @@
END dscal;
PROCEDURE daxpy( n : INTEGER;
- da : FLOAT;
- VAR X : (* ARRAY OF *) FLOAT;
+ da : REAL8;
+ VAR X : (* ARRAY OF *) REAL8;
IncX : INTEGER;
- VAR Y : (* ARRAY OF *) FLOAT;
+ VAR Y : (* ARRAY OF *) REAL8;
IncY : INTEGER);
(*----------------------------------------------------------------*)
@@ -371,10 +373,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 *)
@@ -384,7 +386,7 @@
CONST veclen = 8;
- VAR dtemp : FLOAT;
+ VAR dtemp : REAL8;
i,ix,iy,m : INTEGER;
XX,YY : PVEKTOR;
BEGIN
@@ -423,19 +425,19 @@
END ddot;
PROCEDURE idamax( n : INTEGER;
- VAR X : (* ARRAY OF *) FLOAT;
+ VAR X : (* 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;
XX : PVEKTOR;
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;
XX:=SYSTEM.CAST(PVEKTOR,SYSTEM.ADR(X));
@@ -463,19 +465,19 @@
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;
XX : PVEKTOR;
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;
XX:=SYSTEM.CAST(PVEKTOR,SYSTEM.ADR(X));
@@ -503,13 +505,13 @@
END idamin;
PROCEDURE dasum( dim : CARDINAL;
- VAR X : (* ARRAY OF *) FLOAT;
- IncX : CARDINAL) : FLOAT;
+ VAR X : (* ARRAY OF *) REAL8;
+ IncX : CARDINAL) : REAL8;
CONST max = 6; (* 8 *)
VAR i,ix,m : CARDINAL;
- Sum : FLOAT;
+ Sum : REAL8;
XX : PVEKTOR;
BEGIN
IF (dim = 0) THEN RETURN 0.0; END;
@@ -542,15 +544,15 @@
(*============================= complexe procedures ========================*)
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;
XX,YY : PZVEKTOR;
@@ -631,14 +633,14 @@
END zcopy;
PROCEDURE zdotu( 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;
XX,YY : PZVEKTOR;
BEGIN
@@ -675,14 +677,14 @@
END zdotu;
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;
XX,YY : PZVEKTOR;
BEGIN
@@ -720,11 +722,11 @@
PROCEDURE dznrm2( N : INTEGER;
VAR X : (* ARRAY OF *) LONGCOMPLEX;
- IncX : INTEGER) : LONGREAL;
+ IncX : INTEGER) : REAL8;
VAR norm,scale,ssq,tmp : LONGREAL;
i,ix : INTEGER;
- zw : FLOAT;
+ zw : REAL8;
XX : PZVEKTOR;
BEGIN
XX:=SYSTEM.CAST(PZVEKTOR,SYSTEM.ADR(X));
@@ -765,8 +767,8 @@
END dznrm2;
PROCEDURE zscal( n : INTEGER;
- da : CFLOAT;
- VAR dx : (* ARRAY OF *) CFLOAT;
+ da : COMPLEX16;
+ VAR dx : (* ARRAY OF *) COMPLEX16;
IncX : INTEGER);
CONST veclen = 4;
@@ -800,10 +802,10 @@
END zscal;
PROCEDURE zaxpy( n : INTEGER;
- da : CFLOAT;
- VAR X : (* ARRAY OF *) CFLOAT;
+ za : COMPLEX16;
+ VAR X : (* ARRAY OF *) COMPLEX16;
IncX : INTEGER;
- VAR Y : (* ARRAY OF *) CFLOAT;
+ VAR Y : (* ARRAY OF *) COMPLEX16;
IncY : INTEGER);
CONST veclen = 4;
@@ -811,49 +813,49 @@
VAR i,ix,iy,m : INTEGER;
XX,YY : PZVEKTOR;
BEGIN
- IF (n <= 0 ) THEN RETURN; END;
- IF (da = 0.0) THEN RETURN; END;
+ IF (n <= 0 ) THEN RETURN; END;
+ IF (RE(za) = 0.0) AND (IM(za) = 0.0) THEN RETURN; END;
XX:=SYSTEM.CAST(PZVEKTOR,SYSTEM.ADR(X));
YY:=SYSTEM.CAST(PZVEKTOR,SYSTEM.ADR(Y));
- IF ((IncX <> 1) OR (IncY <> 1)) THEN
+ 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];
+ YY^[iy] := YY^[iy] + za*XX^[ix];
INC(ix, IncX);
INC(iy, IncY);
END;
ELSE (* code for both increments equal to 1 *)
m := n REM veclen;
- IF (m <> 0) THEN
+ IF (m # 0) THEN
FOR i:=0 TO m-1 DO
- YY^[i] := YY^[i] + da*XX^[i];
+ YY^[i] := YY^[i] + za*XX^[i];
END;
IF (n < veclen) THEN RETURN; END;
END;
FOR i:=m TO n-1 BY veclen 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];
+ YY^[i+0] := YY^[i+0] + za*XX^[i+0];
+ YY^[i+1] := YY^[i+1] + za*XX^[i+1];
+ YY^[i+2] := YY^[i+2] + za*XX^[i+2];
+ YY^[i+3] := YY^[i+3] + za*XX^[i+3];
END;
END;
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;
XX,YY : PZVEKTOR;
BEGIN
XX:=SYSTEM.CAST(PZVEKTOR,SYSTEM.ADR(X));