--- a/LibDBlas.mod
+++ b/LibDBlas.mod
@@ -36,6 +36,7 @@
(* dznrm2,zdotc,zscal,zaxpy,zswap,zdrot *)
(* 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 *)
(*------------------------------------------------------------------------*)
(* Testroutinen *)
(* *)
@@ -629,6 +630,50 @@
END;
END zcopy;
+PROCEDURE zdotu( N : INTEGER;
+ VAR X : (* ARRAY OF *) CFLOAT;
+ IncX : INTEGER;
+ VAR Y : (* ARRAY OF *) CFLOAT;
+ IncY : INTEGER) : CFLOAT;
+
+ CONST veclen = 4;
+
+ VAR dtemp : CFLOAT;
+ i,ix,iy,m : INTEGER;
+ XX,YY : PZVEKTOR;
+BEGIN
+ IF (N <= 0) THEN RETURN zero; END;
+
+ XX:=SYSTEM.CAST(PZVEKTOR,SYSTEM.ADR(X));
+ YY:=SYSTEM.CAST(PZVEKTOR,SYSTEM.ADR(Y));
+
+ dtemp := zero;
+ 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];
+ 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 zdotu;
+
PROCEDURE zdotc( N : INTEGER;
VAR X : (* ARRAY OF *) CFLOAT;
IncX : INTEGER;