Numerik / sources / Diff of /LibDBlasM2.mod

Diff of /LibDBlasM2.mod [bfc22e] .. [28b809]

Switch to unified view


...
  (* 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 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      *)
  (*                dnrm2 (aus dnrm2.f)                                     *)
  (* 28.01.16, MRi: Umstellen von dgemv und dgemm auf "Open Array".         *)
...
  (* 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     *)
  (* 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                                                           *)
  (*                                                                        *)
  (* - dgemv in TstDGEMV                                                    *)
  (* - zgemv in TstZGEMV                                                    *)
...
  (*------------------------------------------------------------------------*)

  (* $Id: LibDBlasM2.mod,v 1.11 2018/09/12 13:20:49 mriedl Exp mriedl $ *)

                     IMPORT SYSTEM;

FROM LongMath        IMPORT sqrt;
FROM LongComplexMath IMPORT zero,one,conj,scalarMult;
                     IMPORT Errors;
FROM F77func         IMPORT REAL8,COMPLEX16,MAX0;


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 REAL8;
                     a,e : CARDINAL) : REAL8;

          VAR i       : CARDINAL;
              s,c,y,t : REAL8;
BEGIN
      s := X[a-1];
      c := 0.0;
      FOR i:=a TO e-1 DO
        y := X[i] - c;
...
        s := t;
      END;
      RETURN s;
END SumVek;

PROCEDURE AbsSumVek(VAR X   : ARRAY OF REAL8;
                        a,e : CARDINAL) : REAL8;

          VAR i       : CARDINAL;
              s,c,y,t : REAL8;

BEGIN
      s := ABS(X[a-1]);
      c := 0.0;
      FOR i:=a TO e-1 DO
...
        s := t;
      END;
      RETURN s;
END AbsSumVek;

PROCEDURE ENorm(VAR X   : ARRAY OF REAL8;
                    a,e : CARDINAL)          : REAL8;

          (*---------------------------------------------------------------*)
          (* Argonne national laboratory. MINPACK project. March 1980.     *)
          (* Burton S. Garbow, Kenneth E. Hillstrom, Jorge J. More         *)
          (*                                                               *)
...

          CONST MinLR = 3.834E-20;
                MaxLR = 1.304E+19;
 
          VAR   i                    : CARDINAL;
                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(REAL8,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;
...
  (*========================================================================*)
  (* Pure BLAS based routines                                               *)
  (*========================================================================*)

PROCEDURE dnrm2(    n    : INTEGER;
                VAR X    : ARRAY OF REAL8;
                    IncX : INTEGER): REAL8;

         (*-----------------------------------------------------------------*)
         (* 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 : REAL8;
             zw                   : REAL8;
BEGIN
      IF (n < 1) OR (IncX < 1) THEN
        norm := 0.0 
      ELSIF (n = 1) THEN
        norm := ABS(X[0]);
...
      END;
      RETURN norm;
END dnrm2;

PROCEDURE dswap(    N    : CARDINAL;
                VAR X    : ARRAY OF REAL8;  
                    incX : INTEGER;  
                VAR Y    : ARRAY OF REAL8;
                    incY : INTEGER);

          (*----------------------------------------------------------------*)
          (* Adopted to Modula-2, MRi, 30.03.2016                           *)
          (*----------------------------------------------------------------*)

          CONST level = 4;

          VAR   Xi       : REAL8; 
                Xtmp     : ARRAY [0..level-1] OF REAL8;
                ix,iy    : INTEGER;
                i,m      : CARDINAL;
BEGIN
      IF (N = 0) THEN RETURN; END;

...
        END;
      END;
END dswap;

PROCEDURE dcopy(    N    : INTEGER;
                VAR X    : ARRAY OF REAL8;
                    IncX : INTEGER;
                VAR Y    : ARRAY OF REAL8;
                    IncY : INTEGER);

          (*----------------------------------------------------------------*)
          (* Adopted to Modula-2, MRi, 04.04.2016                           *)
          (*----------------------------------------------------------------*)
...
        END;
     END;
END dcopy;

PROCEDURE drot(    N    : INTEGER;  
               VAR X    : ARRAY OF REAL8;  
                   incX : INTEGER;
               VAR Y    : ARRAY OF REAL8;
                   incY : INTEGER; 
                   c,s  : REAL8);

          (*----------------------------------------------------------------*)
          (* Adopted to Modula-2, MRi, 31.03.2016                           *)
          (*----------------------------------------------------------------*)

          VAR i,ix,iy : INTEGER;
              x1,y1   : REAL8; 
BEGIN
      IF (incX = 1) AND (incY = 1) THEN
        FOR i:=0 TO N-1 DO
          x1 := X[i]; 
          y1 := Y[i]; 
...
          iy := iy + incY; 
        END;
      END;
END drot;
      
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 : REAL8; 
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;
...
      END;
      da := r; db := z; 
END drotg;

PROCEDURE dscal(    n    : INTEGER; 
                    da   : REAL8;
                VAR dx   : ARRAY OF REAL8; 
                    incx : INTEGER);

          (*----------------------------------------------------------------*)
          (* Scales a vector by a constant (UNROLLED version)               *)
          (*----------------------------------------------------------------*)
...
        END;
      END;
END dscal;

PROCEDURE daxpy(    n    : INTEGER;
                    da   : REAL8;
                VAR dx   : ARRAY OF REAL8; 
                    incx : INTEGER;
                VAR dy   : ARRAY OF REAL8; 
                    incy : INTEGER);

          (*----------------------------------------------------------------*)
          (* constant times a vector plus a vector                          *)
          (*----------------------------------------------------------------*)
...
        END;
      END;
END daxpy;

PROCEDURE ddot(    N    : INTEGER;
               VAR X    : ARRAY OF REAL8;
                   IncX : INTEGER;
               VAR Y    : ARRAY OF REAL8;
                   IncY : INTEGER) : REAL8;

          (*----------------------------------------------------------------*)
          (* 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     : REAL8;
                i,ix,iy,m : INTEGER;
BEGIN
      IF (N <= 0) THEN RETURN 0.0; END;
      dtemp := 0.0;
      IF (IncX = 1) AND (IncY = 1) THEN
...
      END;
      RETURN dtemp;
END ddot;

PROCEDURE idamax(    n    : INTEGER;
                 VAR dx   : ARRAY OF REAL8; 
                     incx : INTEGER) : INTEGER;

          (*----------------------------------------------------------------*)
          (* Finds the index of element having max. absolute value.         *)
          (*----------------------------------------------------------------*)

          VAR dmax       : REAL8;
              i,ix,itemp : INTEGER;

BEGIN
      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]);
        ix := 1 + incx;
        FOR i:=1 TO n-1 DO
...
      END;
      RETURN itemp;
END idamax;

PROCEDURE idamin(    n    : INTEGER;
                 VAR X    : ARRAY OF REAL8;
                     IncX : INTEGER) : INTEGER;

          (*----------------------------------------------------------------*)
          (* Finds the index of element having min. absolute value.         *)
          (*----------------------------------------------------------------*)

          VAR dmin       : REAL8;
              i,ix,itemp : INTEGER;
BEGIN
      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 *)
        dmin := ABS(X[0]);
        ix := 1 + IncX;
...
      END;
      RETURN itemp + 1;
END idamin;

PROCEDURE dasum(    dim      : CARDINAL;
                VAR X        : ARRAY OF REAL8;
                    Inc      : CARDINAL; (* Inkrementwert >= 1*)
                    Unrolled : BOOLEAN) : REAL8;

          CONST max      = 6; (* 8 *)

          VAR   i,m : CARDINAL;
                Sum : REAL8;
BEGIN
      IF (dim = 0) THEN RETURN 0.0; END;
      IF (Inc < 1) THEN 
        Errors.FatalError("Inkrement kleiner 1 (LibDBlas.dasum) !"); 
      END;
...
      RETURN Sum;
END dasum;

PROCEDURE dgemv (   Trans  : CHAR;
                    M,N    : INTEGER;
                    Alpha  : REAL8;
                VAR A      : ARRAY OF ARRAY OF REAL8;   (* [1..LDA][1..*] *)
                    LDA    : INTEGER;
                VAR X      : ARRAY OF REAL8;   (* [1..*] *)
                    IncX   : INTEGER;
                    Beta   : REAL8;
                VAR Y      : ARRAY OF REAL8;   (* [1..*] *)
                    IncY   : INTEGER);

          VAR Temp                              : REAL8;
              i,j,iy,jx,jy,kx,ky,Info,LenX,LenY : INTEGER;
BEGIN (* Testen *)
      (* Test the input parameters. *)

      IF (NOT (CAP(Trans) = 'N')) AND (NOT (CAP(Trans) = 'T')) THEN
...
      END; (* IF Trans *)
END dgemv;

PROCEDURE dgemm(    TransA,TransB : CHAR;
                    M,N,K         : INTEGER;
                    Alpha         : REAL8;
                VAR A             : ARRAY OF ARRAY OF REAL8; (* [1..LDA][1..*] *)
                    LDA           : INTEGER;
                VAR B             : ARRAY OF ARRAY OF REAL8; (* [1..LDB][1..*] *)
                    LDB           : INTEGER;
                    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           : REAL8;
               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.      *)
...

      END; (* IF NotB *)
END dgemm;

PROCEDURE dger(    m,n   : CARDINAL;
                   Alpha : REAL8;
               VAR X     : ARRAY OF REAL8;
                   IncX  : CARDINAL;
               VAR Y     : ARRAY OF REAL8;
                   IncY  : CARDINAL;
               VAR A     : ARRAY OF ARRAY OF REAL8;
                   lda   : CARDINAL);

           VAR i,j,ix,jy : CARDINAL;
               AlphaX    : REAL8;
BEGIN
      IF (lda-1 > HIGH(A)) THEN
        Errors.Fehlerflag:="lda-1 > HIGH(A) (LibDBlas.dger) !";
        Errors.Fehler:=TRUE;
        RETURN;
...
END dger;

(*=========================== Complex valued routines =====================*)

PROCEDURE zswap(    N    : CARDINAL;
                VAR X    : ARRAY OF COMPLEX16;
                    IncX : INTEGER;
                VAR Y    : ARRAY OF COMPLEX16;
                    IncY : INTEGER);

          CONST veclen = 4;

          VAR   Xi    : COMPLEX16;
                Xtmp  : ARRAY [0..veclen-1] OF COMPLEX16;
                ix,iy : INTEGER;
                i,m   : CARDINAL;
BEGIN
      IF (N = 0) THEN RETURN; END;

...
        END;
      END;
END zswap;

PROCEDURE zcopy(    N    : INTEGER;
                VAR X    : ARRAY OF COMPLEX16;
                    IncX : INTEGER;
                VAR Y    : ARRAY OF COMPLEX16;
                    IncY : INTEGER);

          (*----------------------------------------------------------------*)
          (* Adopted to Modula-2, MRi, 04.04.2016, complex version 09.08.18 *)
          (*----------------------------------------------------------------*)
...
        END;
      END;
END zcopy;

PROCEDURE zdotc(    N    : INTEGER;
                VAR X    : ARRAY OF COMPLEX16;
                    IncX : INTEGER;
                VAR Y    : ARRAY OF COMPLEX16;
                    IncY : INTEGER) : COMPLEX16;

          CONST veclen = 4;

          VAR   dtemp     : COMPLEX16;
                i,ix,iy,m : INTEGER;
BEGIN
      IF (N <= 0) THEN RETURN zero; END;

      dtemp := zero;
...
        m := N MOD veclen;
        IF (m # 0) THEN
          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
          dtemp:=dtemp + conj(X[i+0])*Y[i+0] + conj(X[i+1])*Y[i+1] +
                         conj(X[i+2])*Y[i+2] + conj(X[i+3])*Y[i+3];
...
      END;
      RETURN dtemp;
END zdotc;

PROCEDURE dznrm2(    N    : INTEGER;
                 VAR X    : ARRAY OF COMPLEX16;
                     IncX : INTEGER) : REAL8;

          VAR norm,scale,ssq,tmp : REAL8;
              i,ix               : INTEGER;
              zw                 : REAL8;
BEGIN
      IF (N < 1) OR (IncX < 1) THEN
        norm := 0.0;
      ELSE
        scale := 0.0;
...
      END; (* IF *)
      RETURN norm;
END dznrm2;

PROCEDURE zscal(    n    : INTEGER;
                    da   : COMPLEX16;
                VAR X    : ARRAY OF COMPLEX16;
                    IncX : INTEGER);

          CONST veclen = 4;

          VAR   i,m    : INTEGER;
...
        END;
      END;
END zscal;

PROCEDURE zaxpy(    n    : INTEGER;
                    da   : COMPLEX16;
                VAR X    : ARRAY OF COMPLEX16;
                    IncX : INTEGER;
                VAR Y    : ARRAY OF COMPLEX16;
                    IncY : INTEGER);

          CONST veclen    = 4;

          VAR   i,ix,iy,m : INTEGER;
...
        END;
      END;
END zaxpy;

PROCEDURE zdrot(    N    : INTEGER;
                VAR X    : ARRAY OF COMPLEX16;
                    IncX : INTEGER;
                VAR Y    : ARRAY OF COMPLEX16;
                    IncY : INTEGER;
                    c,s  : REAL8);

          VAR i,ix,iy : INTEGER;
              tmp     : COMPLEX16;
BEGIN
      IF (IncX = 1) AND (IncY = 1) THEN
        FOR i:=0 TO N-1 DO
          tmp    := scalarMult( c,X[i]) + scalarMult(s,Y[i]);
          Y[i] := scalarMult(-s,X[i]) + scalarMult(c,Y[i]);
...
      END;
END zdrot;

PROCEDURE zgemv(    Trans : CHAR; 
                    M,N   : INTEGER;
                    Alpha : COMPLEX16;
                VAR A     : ARRAY OF ARRAY OF COMPLEX16;
                    lda   : INTEGER;
                VAR X     : ARRAY OF COMPLEX16;
                    IncX  : INTEGER;
                    Beta  : COMPLEX16;
                VAR Y     : ARRAY OF COMPLEX16;
                    IncY  : INTEGER);

          VAR i,j,iy,jx,kx,ky : INTEGER;
              lenx,leny,info  : INTEGER;
              s               : COMPLEX16;
              noconj          : BOOLEAN;
BEGIN
      (* test the input parameters *)

      info := 0;
...
      END; (* IF Trans *)
END zgemv;

PROCEDURE zgemm(    TransA,TransB : CHAR;
                    M,N,K         : INTEGER;
                    Alpha         : COMPLEX16;
                VAR A             : ARRAY OF ARRAY OF COMPLEX16;
                    LDA           : INTEGER;
                VAR B             : ARRAY OF ARRAY OF COMPLEX16;
                    LDB           : INTEGER;
                    Beta          : COMPLEX16;
                VAR C             : ARRAY OF ARRAY OF COMPLEX16;
                    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           : COMPLEX16;
                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. *)

	a/LibDBlasM2.mod		b/LibDBlasM2.mod
	...		...
9	(* Letzte Bearbeitung *)	9	(* Letzte Bearbeitung *)
10	(* *)	10	(* *)
11	(* 23.11.93, MRi: Erstellen der 1. Version *)	11	(* 23.11.93, MRi: Erstellen der 1. Version *)
12	(* Aug. 95, MRi: Erweiterung *)	12	(* Aug. 95, MRi: Erweiterung *)
13	(* 03.11.95, MRi: Durchsicht *)	13	(* 03.11.95, MRi: Durchsicht *)
14	(* 13.10.15, MRi: Ersetzen von FLOAT(#) durch VAL(FLOAT,#) *)	14	(* 13.10.15, MRi: Ersetzen von REAL8(#) durch VAL(REAL8,#) *)
15	(* 30.11.15, MRi: Erstelle der Routinen dscal,daxpy,idmax aus Linpack *)	15	(* 30.11.15, MRi: Erstelle der Routinen dscal,daxpy,idmax aus Linpack *)
16	(* Benchmark aus Beispielen von Stony Brook M2 *)	16	(* Benchmark aus Beispielen von Stony Brook M2 *)
17	(* 01.12.15, MRi: Umbenennen von BlasLib auf LibDBlas, einfuegen von *)	17	(* 01.12.15, MRi: Umbenennen von BlasLib auf LibDBlas, einfuegen von *)
18	(* dnrm2 (aus dnrm2.f) *)	18	(* dnrm2 (aus dnrm2.f) *)
19	(* 28.01.16, MRi: Umstellen von dgemv und dgemm auf "Open Array". *)	19	(* 28.01.16, MRi: Umstellen von dgemv und dgemm auf "Open Array". *)
	...		...
26	(* 21.10.16, MRi: idamin eingefuegt *)	26	(* 21.10.16, MRi: idamin eingefuegt *)
27	(* 28.10.17, MRi: zgemm eingefuegt *)	27	(* 28.10.17, MRi: zgemm eingefuegt *)
28	(* 29.10.17, MRi: Rolle von M,N in dgemm und dgemv an zgemm angepasst *)	28	(* 29.10.17, MRi: Rolle von M,N in dgemm und dgemv an zgemm angepasst *)
29	(* 11.09.18, MRi: zgemv eingefuegt *)	29	(* 11.09.18, MRi: zgemv eingefuegt *)
30	(* 12.09.18, MRi: zswap,zcopy,zdotc,dznrm2,zscal,zaxpy,zdrot eingefuegt *)	30	(* 12.09.18, MRi: zswap,zcopy,zdotc,dznrm2,zscal,zaxpy,zdrot eingefuegt *)
		31	(* 01.02.19, MRi: Abfragen in idamx,idamin korrigiert *)
31	(------------------------------------------------------------------------)	32	(------------------------------------------------------------------------)
32	(* Testroutinen *)	33	(* Testroutinen *)
33	(* *)	34	(* *)
34	(* - dgemv in TstDGEMV *)	35	(* - dgemv in TstDGEMV *)
35	(* - zgemv in TstZGEMV *)	36	(* - zgemv in TstZGEMV *)
	...		...
45	(------------------------------------------------------------------------)	46	(------------------------------------------------------------------------)
46		47
47	(* $Id: LibDBlasM2.mod,v 1.11 2018/09/12 13:20:49 mriedl Exp mriedl $ *)	48	(* $Id: LibDBlasM2.mod,v 1.11 2018/09/12 13:20:49 mriedl Exp mriedl $ *)
48		49
49	IMPORT SYSTEM;	50	IMPORT SYSTEM;
50	FROM Deklera IMPORT FLOAT,CFLOAT; (* REAL/COMPLEX Type *)
51	FROM LongMath IMPORT sqrt;	51	FROM LongMath IMPORT sqrt;
52	FROM LongComplexMath IMPORT zero,one,conj,scalarMult;	52	FROM LongComplexMath IMPORT zero,one,conj,scalarMult;
53	IMPORT Errors;	53	IMPORT Errors;
54	FROM F77func IMPORT MAX0;	54	FROM F77func IMPORT REAL8,COMPLEX16,MAX0;
55		55
56		56
57	PROCEDURE IMax(a,b : INTEGER) : INTEGER;	57	PROCEDURE IMax(a,b : INTEGER) : INTEGER;
58		58
59	BEGIN	59	BEGIN
60	IF (a > b) THEN RETURN a; ELSE RETURN b; END;	60	IF (a > b) THEN RETURN a; ELSE RETURN b; END;
61	END IMax;	61	END IMax;
62		62
63	PROCEDURE SumVek(VAR X : ARRAY OF FLOAT;	63	PROCEDURE SumVek(VAR X : ARRAY OF REAL8;
64	a,e : CARDINAL) : FLOAT;	64	a,e : CARDINAL) : REAL8;
65		65
66	VAR i : CARDINAL;	66	VAR i : CARDINAL;
67	s,c,y,t : FLOAT;	67	s,c,y,t : REAL8;
68	BEGIN	68	BEGIN
69	s := X[a-1];	69	s := X[a-1];
70	c := 0.0;	70	c := 0.0;
71	FOR i:=a TO e-1 DO	71	FOR i:=a TO e-1 DO
72	y := X[i] - c;	72	y := X[i] - c;
	...		...
75	s := t;	75	s := t;
76	END;	76	END;
77	RETURN s;	77	RETURN s;
78	END SumVek;	78	END SumVek;
79		79
80	PROCEDURE AbsSumVek(VAR X : ARRAY OF FLOAT;	80	PROCEDURE AbsSumVek(VAR X : ARRAY OF REAL8;
81	a,e : CARDINAL) : FLOAT;	81	a,e : CARDINAL) : REAL8;
82		82
83	VAR i : CARDINAL;	83	VAR i : CARDINAL;
84	s,c,y,t : FLOAT;	84	s,c,y,t : REAL8;
85		85
86	BEGIN	86	BEGIN
87	s := ABS(X[a-1]);	87	s := ABS(X[a-1]);
88	c := 0.0;	88	c := 0.0;
89	FOR i:=a TO e-1 DO	89	FOR i:=a TO e-1 DO
	...		...
93	s := t;	93	s := t;
94	END;	94	END;
95	RETURN s;	95	RETURN s;
96	END AbsSumVek;	96	END AbsSumVek;
97		97
98	PROCEDURE ENorm(VAR X : ARRAY OF FLOAT;	98	PROCEDURE ENorm(VAR X : ARRAY OF REAL8;
99	a,e : CARDINAL) : FLOAT;	99	a,e : CARDINAL) : REAL8;
100		100
101	(---------------------------------------------------------------)	101	(---------------------------------------------------------------)
102	(* Argonne national laboratory. MINPACK project. March 1980. *)	102	(* Argonne national laboratory. MINPACK project. March 1980. *)
103	(* Burton S. Garbow, Kenneth E. Hillstrom, Jorge J. More *)	103	(* Burton S. Garbow, Kenneth E. Hillstrom, Jorge J. More *)
104	(* *)	104	(* *)
	...		...
107		107
108	CONST MinLR = 3.834E-20;	108	CONST MinLR = 3.834E-20;
109	MaxLR = 1.304E+19;	109	MaxLR = 1.304E+19;
110		110
111	VAR i : CARDINAL;	111	VAR i : CARDINAL;
112	S1,S2,S3,zw : FLOAT;	112	S1,S2,S3,zw : REAL8;
113	Giant,AbsX,Max1,Max3 : FLOAT;	113	Giant,AbsX,Max1,Max3 : REAL8;
114	Norm : FLOAT;	114	Norm : REAL8;
115	BEGIN	115	BEGIN
116	S1:=0.0; S2:=0.0; S3:=0.0;	116	S1:=0.0; S2:=0.0; S3:=0.0;
117	Max1:=0.0; Max3:=0.0;	117	Max1:=0.0; Max3:=0.0;
118	Giant := MaxLR / VAL(FLOAT,e-a+1);	118	Giant := MaxLR / VAL(REAL8,e-a+1);
119	FOR i:=a-1 TO e-1 DO	119	FOR i:=a-1 TO e-1 DO
120	AbsX := ABS(X[i]);	120	AbsX := ABS(X[i]);
121	IF (AbsX < MinLR) THEN (* Summiere kleine Elemente. *)	121	IF (AbsX < MinLR) THEN (* Summiere kleine Elemente. *)
122	IF (AbsX > Max3) THEN	122	IF (AbsX > Max3) THEN
123	zw := Max3 / AbsX;	123	zw := Max3 / AbsX;
	...		...
161	(========================================================================)	161	(========================================================================)
162	(* Pure BLAS based routines *)	162	(* Pure BLAS based routines *)
163	(========================================================================)	163	(========================================================================)
164		164
165	PROCEDURE dnrm2( n : INTEGER;	165	PROCEDURE dnrm2( n : INTEGER;
166	VAR X : ARRAY OF FLOAT;	166	VAR X : ARRAY OF REAL8;
167	IncX : INTEGER): FLOAT;	167	IncX : INTEGER): REAL8;
168		168
169	(-----------------------------------------------------------------)	169	(-----------------------------------------------------------------)
170	(* This version written on 25-October-1982. *)	170	(* This version written on 25-October-1982. *)
171	(* Modified on 14-October-1993 to inline the call to DLASSQ. *)	171	(* Modified on 14-October-1993 to inline the call to DLASSQ. *)
172	(* Sven Hammarling, Nag Ltd. *)	172	(* Sven Hammarling, Nag Ltd. *)
173	(* Adapted to Modula-2 01.12.2015, M.Riedl *)	173	(* Adapted to Modula-2 01.12.2015, M.Riedl *)
174	(-----------------------------------------------------------------)	174	(-----------------------------------------------------------------)
175		175
176	VAR i,ix : INTEGER;	176	VAR i,ix : INTEGER;
177	absxi,norm,scale,ssq : FLOAT;	177	absxi,norm,scale,ssq : REAL8;
178	zw : FLOAT;	178	zw : REAL8;
179	BEGIN	179	BEGIN
180	IF (n < 1) OR (IncX < 1) THEN	180	IF (n < 1) OR (IncX < 1) THEN
181	norm := 0.0	181	norm := 0.0
182	ELSIF (n = 1) THEN	182	ELSIF (n = 1) THEN
183	norm := ABS(X[0]);	183	norm := ABS(X[0]);
	...		...
203	END;	203	END;
204	RETURN norm;	204	RETURN norm;
205	END dnrm2;	205	END dnrm2;
206		206
207	PROCEDURE dswap( N : CARDINAL;	207	PROCEDURE dswap( N : CARDINAL;
208	VAR X : ARRAY OF FLOAT;	208	VAR X : ARRAY OF REAL8;
209	incX : INTEGER;	209	incX : INTEGER;
210	VAR Y : ARRAY OF FLOAT;	210	VAR Y : ARRAY OF REAL8;
211	incY : INTEGER);	211	incY : INTEGER);
212		212
213	(----------------------------------------------------------------)	213	(----------------------------------------------------------------)
214	(* Adopted to Modula-2, MRi, 30.03.2016 *)	214	(* Adopted to Modula-2, MRi, 30.03.2016 *)
215	(----------------------------------------------------------------)	215	(----------------------------------------------------------------)
216		216
217	CONST level = 4;	217	CONST level = 4;
218		218
219	VAR Xi : FLOAT;	219	VAR Xi : REAL8;
220	Xtmp : ARRAY [0..level-1] OF FLOAT;	220	Xtmp : ARRAY [0..level-1] OF REAL8;
221	ix,iy : INTEGER;	221	ix,iy : INTEGER;
222	i,m : CARDINAL;	222	i,m : CARDINAL;
223	BEGIN	223	BEGIN
224	IF (N = 0) THEN RETURN; END;	224	IF (N = 0) THEN RETURN; END;
225		225
	...		...
246	END;	246	END;
247	END;	247	END;
248	END dswap;	248	END dswap;
249		249
250	PROCEDURE dcopy( N : INTEGER;	250	PROCEDURE dcopy( N : INTEGER;
251	VAR X : ARRAY OF FLOAT;	251	VAR X : ARRAY OF REAL8;
252	IncX : INTEGER;	252	IncX : INTEGER;
253	VAR Y : ARRAY OF FLOAT;	253	VAR Y : ARRAY OF REAL8;
254	IncY : INTEGER);	254	IncY : INTEGER);
255		255
256	(----------------------------------------------------------------)	256	(----------------------------------------------------------------)
257	(* Adopted to Modula-2, MRi, 04.04.2016 *)	257	(* Adopted to Modula-2, MRi, 04.04.2016 *)
258	(----------------------------------------------------------------)	258	(----------------------------------------------------------------)
	...		...
284	END;	284	END;
285	END;	285	END;
286	END dcopy;	286	END dcopy;
287		287
288	PROCEDURE drot( N : INTEGER;	288	PROCEDURE drot( N : INTEGER;
289	VAR X : ARRAY OF FLOAT;	289	VAR X : ARRAY OF REAL8;
290	incX : INTEGER;	290	incX : INTEGER;
291	VAR Y : ARRAY OF FLOAT;	291	VAR Y : ARRAY OF REAL8;
292	incY : INTEGER;	292	incY : INTEGER;
293	c,s : FLOAT);	293	c,s : REAL8);
294		294
295	(----------------------------------------------------------------)	295	(----------------------------------------------------------------)
296	(* Adopted to Modula-2, MRi, 31.03.2016 *)	296	(* Adopted to Modula-2, MRi, 31.03.2016 *)
297	(----------------------------------------------------------------)	297	(----------------------------------------------------------------)
298		298
299	VAR i,ix,iy : INTEGER;	299	VAR i,ix,iy : INTEGER;
300	x1,y1 : FLOAT;	300	x1,y1 : REAL8;
301	BEGIN	301	BEGIN
302	IF (incX = 1) AND (incY = 1) THEN	302	IF (incX = 1) AND (incY = 1) THEN
303	FOR i:=0 TO N-1 DO	303	FOR i:=0 TO N-1 DO
304	x1 := X[i];	304	x1 := X[i];
305	y1 := Y[i];	305	y1 := Y[i];
	...		...
318	iy := iy + incY;	318	iy := iy + incY;
319	END;	319	END;
320	END;	320	END;
321	END drot;	321	END drot;
322		322
323	PROCEDURE drotg(VAR da : FLOAT;	323	PROCEDURE drotg(VAR da : REAL8;
324	VAR db : FLOAT;	324	VAR db : REAL8;
325	VAR c : FLOAT;	325	VAR c : REAL8;
326	VAR s : FLOAT);	326	VAR s : REAL8);
327		327
328	(----------------------------------------------------------------)	328	(----------------------------------------------------------------)
329	(* Adopted to Modula-2, MRi, 31.03.2016 *)	329	(* Adopted to Modula-2, MRi, 31.03.2016 *)
330	(----------------------------------------------------------------)	330	(----------------------------------------------------------------)
331		331
332	VAR r,z,roe,daos,dbos,scale : FLOAT;	332	VAR r,z,roe,daos,dbos,scale : REAL8;
333	BEGIN	333	BEGIN
334	IF (ABS(da) > ABS(db)) THEN roe := da; ELSE roe := db; END;	334	IF (ABS(da) > ABS(db)) THEN roe := da; ELSE roe := db; END;
335	scale := ABS(da) + ABS(db);	335	scale := ABS(da) + ABS(db);
336	IF (scale # 0.0) THEN	336	IF (scale # 0.0) THEN
337	daos := da / scale;	337	daos := da / scale;
	...		...
356	END;	356	END;
357	da := r; db := z;	357	da := r; db := z;
358	END drotg;	358	END drotg;
359		359
360	PROCEDURE dscal( n : INTEGER;	360	PROCEDURE dscal( n : INTEGER;
361	da : FLOAT;	361	da : REAL8;
362	VAR dx : ARRAY OF FLOAT;	362	VAR dx : ARRAY OF REAL8;
363	incx : INTEGER);	363	incx : INTEGER);
364		364
365	(----------------------------------------------------------------)	365	(----------------------------------------------------------------)
366	(* Scales a vector by a constant (UNROLLED version) *)	366	(* Scales a vector by a constant (UNROLLED version) *)
367	(----------------------------------------------------------------)	367	(----------------------------------------------------------------)
	...		...
392	END;	392	END;
393	END;	393	END;
394	END dscal;	394	END dscal;
395		395
396	PROCEDURE daxpy( n : INTEGER;	396	PROCEDURE daxpy( n : INTEGER;
397	da : FLOAT;	397	da : REAL8;
398	VAR dx : ARRAY OF FLOAT;	398	VAR dx : ARRAY OF REAL8;
399	incx : INTEGER;	399	incx : INTEGER;
400	VAR dy : ARRAY OF FLOAT;	400	VAR dy : ARRAY OF REAL8;
401	incy : INTEGER);	401	incy : INTEGER);
402		402
403	(----------------------------------------------------------------)	403	(----------------------------------------------------------------)
404	(* constant times a vector plus a vector *)	404	(* constant times a vector plus a vector *)
405	(----------------------------------------------------------------)	405	(----------------------------------------------------------------)
	...		...
436	END;	436	END;
437	END;	437	END;
438	END daxpy;	438	END daxpy;
439		439
440	PROCEDURE ddot( N : INTEGER;	440	PROCEDURE ddot( N : INTEGER;
441	VAR X : ARRAY OF FLOAT;	441	VAR X : ARRAY OF REAL8;
442	IncX : INTEGER;	442	IncX : INTEGER;
443	VAR Y : ARRAY OF FLOAT;	443	VAR Y : ARRAY OF REAL8;
444	IncY : INTEGER) : FLOAT;	444	IncY : INTEGER) : REAL8;
445		445
446	(----------------------------------------------------------------)	446	(----------------------------------------------------------------)
447	(* Forms the dot product of two vectors. Uses unrolled loops for *)	447	(* Forms the dot product of two vectors. Uses unrolled loops for *)
448	(* increments equal to one. *)	448	(* increments equal to one. *)
449	(* Adopted to Modula-2, MRi, 06.09.2015 *)	449	(* Adopted to Modula-2, MRi, 06.09.2015 *)
450	(----------------------------------------------------------------)	450	(----------------------------------------------------------------)
451		451
452	CONST veclen = 8;	452	CONST veclen = 8;
453		453
454	VAR dtemp : FLOAT;	454	VAR dtemp : REAL8;
455	i,ix,iy,m : INTEGER;	455	i,ix,iy,m : INTEGER;
456	BEGIN	456	BEGIN
457	IF (N <= 0) THEN RETURN 0.0; END;	457	IF (N <= 0) THEN RETURN 0.0; END;
458	dtemp := 0.0;	458	dtemp := 0.0;
459	IF (IncX = 1) AND (IncY = 1) THEN	459	IF (IncX = 1) AND (IncY = 1) THEN
	...		...
483	END;	483	END;
484	RETURN dtemp;	484	RETURN dtemp;
485	END ddot;	485	END ddot;
486		486
487	PROCEDURE idamax( n : INTEGER;	487	PROCEDURE idamax( n : INTEGER;
488	VAR dx : ARRAY OF FLOAT;	488	VAR dx : ARRAY OF REAL8;
489	incx : INTEGER) : INTEGER;	489	incx : INTEGER) : INTEGER;
490		490
491	(----------------------------------------------------------------)	491	(----------------------------------------------------------------)
492	(* Finds the index of element having max. absolute value. *)	492	(* Finds the index of element having max. absolute value. *)
493	(----------------------------------------------------------------)	493	(----------------------------------------------------------------)
494		494
495	VAR dmax : FLOAT;	495	VAR dmax : REAL8;
496	i,ix,itemp : INTEGER;	496	i,ix,itemp : INTEGER;
497		497
498	BEGIN	498	BEGIN
499	IF (n < 1) THEN RETURN -1; END;	499	IF (n < 1) THEN RETURN 0; END;
500	IF (n = 1) THEN RETURN 0; END;	500	IF (n = 1) THEN RETURN 1; END;
501	itemp:=0;	501	itemp:=0;
502	IF (incx <> 1) THEN (* code for increment not equal to 1 *)	502	IF (incx <> 1) THEN (* code for increment not equal to 1 *)
503	dmax := ABS(dx[0]);	503	dmax := ABS(dx[0]);
504	ix := 1 + incx;	504	ix := 1 + incx;
505	FOR i:=1 TO n-1 DO	505	FOR i:=1 TO n-1 DO
	...		...
520	END;	520	END;
521	RETURN itemp;	521	RETURN itemp;
522	END idamax;	522	END idamax;
523		523
524	PROCEDURE idamin( n : INTEGER;	524	PROCEDURE idamin( n : INTEGER;
525	VAR X : ARRAY OF FLOAT;	525	VAR X : ARRAY OF REAL8;
526	IncX : INTEGER) : INTEGER;	526	IncX : INTEGER) : INTEGER;
527		527
528	(----------------------------------------------------------------)	528	(----------------------------------------------------------------)
529	(* Finds the index of element having min. absolute value. *)	529	(* Finds the index of element having min. absolute value. *)
530	(----------------------------------------------------------------)	530	(----------------------------------------------------------------)
531		531
532	VAR dmin : FLOAT;	532	VAR dmin : REAL8;
533	i,ix,itemp : INTEGER;	533	i,ix,itemp : INTEGER;
534	BEGIN	534	BEGIN
535	IF (n < 1) THEN RETURN -1; END;	535	IF (n < 1) THEN RETURN 0; END;
536	IF (n = 1) THEN RETURN 0; END;	536	IF (n = 1) THEN RETURN 1; END;
537		537
538	itemp:=0;	538	itemp:=0;
539	IF (IncX <> 1) THEN (* code for Increment not equal to 1 *)	539	IF (IncX <> 1) THEN (* code for Increment not equal to 1 *)
540	dmin := ABS(X[0]);	540	dmin := ABS(X[0]);
541	ix := 1 + IncX;	541	ix := 1 + IncX;
	...		...
557	END;	557	END;
558	RETURN itemp + 1;	558	RETURN itemp + 1;
559	END idamin;	559	END idamin;
560		560
561	PROCEDURE dasum( dim : CARDINAL;	561	PROCEDURE dasum( dim : CARDINAL;
562	VAR X : ARRAY OF FLOAT;	562	VAR X : ARRAY OF REAL8;
563	Inc : CARDINAL; (* Inkrementwert >= 1*)	563	Inc : CARDINAL; (* Inkrementwert >= 1*)
564	Unrolled : BOOLEAN) : FLOAT;	564	Unrolled : BOOLEAN) : REAL8;
565		565
566	CONST max = 6; (* 8 *)	566	CONST max = 6; (* 8 *)
567		567
568	VAR i,m : CARDINAL;	568	VAR i,m : CARDINAL;
569	Sum : FLOAT;	569	Sum : REAL8;
570	BEGIN	570	BEGIN
571	IF (dim = 0) THEN RETURN 0.0; END;	571	IF (dim = 0) THEN RETURN 0.0; END;
572	IF (Inc < 1) THEN	572	IF (Inc < 1) THEN
573	Errors.FatalError("Inkrement kleiner 1 (LibDBlas.dasum) !");	573	Errors.FatalError("Inkrement kleiner 1 (LibDBlas.dasum) !");
574	END;	574	END;
	...		...
600	RETURN Sum;	600	RETURN Sum;
601	END dasum;	601	END dasum;
602		602
603	PROCEDURE dgemv ( Trans : CHAR;	603	PROCEDURE dgemv ( Trans : CHAR;
604	M,N : INTEGER;	604	M,N : INTEGER;
605	Alpha : FLOAT;	605	Alpha : REAL8;
606	VAR A : ARRAY OF ARRAY OF FLOAT; (* [1..LDA][1..] )	606	VAR A : ARRAY OF ARRAY OF REAL8; (* [1..LDA][1..] )
607	LDA : INTEGER;	607	LDA : INTEGER;
608	VAR X : ARRAY OF FLOAT; (* [1..] )	608	VAR X : ARRAY OF REAL8; (* [1..] )
609	IncX : INTEGER;	609	IncX : INTEGER;
610	Beta : FLOAT;	610	Beta : REAL8;
611	VAR Y : ARRAY OF FLOAT; (* [1..] )	611	VAR Y : ARRAY OF REAL8; (* [1..] )
612	IncY : INTEGER);	612	IncY : INTEGER);
613		613
614	VAR Temp : FLOAT;	614	VAR Temp : REAL8;
615	i,j,iy,jx,jy,kx,ky,Info,LenX,LenY : INTEGER;	615	i,j,iy,jx,jy,kx,ky,Info,LenX,LenY : INTEGER;
616	BEGIN (* Testen *)	616	BEGIN (* Testen *)
617	(* Test the input parameters. *)	617	(* Test the input parameters. *)
618		618
619	IF (NOT (CAP(Trans) = 'N')) AND (NOT (CAP(Trans) = 'T')) THEN	619	IF (NOT (CAP(Trans) = 'N')) AND (NOT (CAP(Trans) = 'T')) THEN
	...		...
740	END; (* IF Trans *)	740	END; (* IF Trans *)
741	END dgemv;	741	END dgemv;
742		742
743	PROCEDURE dgemm( TransA,TransB : CHAR;	743	PROCEDURE dgemm( TransA,TransB : CHAR;
744	M,N,K : INTEGER;	744	M,N,K : INTEGER;
745	Alpha : FLOAT;	745	Alpha : REAL8;
746	VAR A : ARRAY OF ARRAY OF FLOAT; (* [1..LDA][1..] )	746	VAR A : ARRAY OF ARRAY OF REAL8; (* [1..LDA][1..] )
747	LDA : INTEGER;	747	LDA : INTEGER;
748	VAR B : ARRAY OF ARRAY OF FLOAT; (* [1..LDB][1..] )	748	VAR B : ARRAY OF ARRAY OF REAL8; (* [1..LDB][1..] )
749	LDB : INTEGER;	749	LDB : INTEGER;
750	Beta : FLOAT;	750	Beta : REAL8;
751	VAR C : ARRAY OF ARRAY OF FLOAT; (* [1..LDC][1..] )	751	VAR C : ARRAY OF ARRAY OF REAL8; (* [1..LDC][1..] )
752	LDC : INTEGER);	752	LDC : INTEGER);
753		753
754	VAR NotA,NotB : BOOLEAN;	754	VAR NotA,NotB : BOOLEAN;
755	i,j,k,Info,NRowA,NRowB : INTEGER;	755	i,j,k,Info,NRowA,NRowB : INTEGER;
756	Temp,Aki,Aik : FLOAT;	756	Temp,Aki,Aik : REAL8;
757	Nm1,Mm1,Km1 : INTEGER;	757	Nm1,Mm1,Km1 : INTEGER;
758	BEGIN	758	BEGIN
759	(* Set NotA and NotB as true if A and B respectively are not *)	759	(* Set NotA and NotB as true if A and B respectively are not *)
760	(* transposed and set NRowA, NColA and NRowB as the number of rows *)	760	(* transposed and set NRowA, NColA and NRowB as the number of rows *)
761	(* and columns of A and the number of rows of B respectively. *)	761	(* and columns of A and the number of rows of B respectively. *)
	...		...
897		897
898	END; (* IF NotB *)	898	END; (* IF NotB *)
899	END dgemm;	899	END dgemm;
900		900
901	PROCEDURE dger( m,n : CARDINAL;	901	PROCEDURE dger( m,n : CARDINAL;
902	Alpha : FLOAT;	902	Alpha : REAL8;
903	VAR X : ARRAY OF FLOAT;	903	VAR X : ARRAY OF REAL8;
904	IncX : CARDINAL;	904	IncX : CARDINAL;
905	VAR Y : ARRAY OF FLOAT;	905	VAR Y : ARRAY OF REAL8;
906	IncY : CARDINAL;	906	IncY : CARDINAL;
907	VAR A : ARRAY OF ARRAY OF FLOAT;	907	VAR A : ARRAY OF ARRAY OF REAL8;
908	lda : CARDINAL);	908	lda : CARDINAL);
909		909
910	VAR i,j,ix,jy : CARDINAL;	910	VAR i,j,ix,jy : CARDINAL;
911	AlphaX : FLOAT;	911	AlphaX : REAL8;
912	BEGIN	912	BEGIN
913	IF (lda-1 > HIGH(A)) THEN	913	IF (lda-1 > HIGH(A)) THEN
914	Errors.Fehlerflag:="lda-1 > HIGH(A) (LibDBlas.dger) !";	914	Errors.Fehlerflag:="lda-1 > HIGH(A) (LibDBlas.dger) !";
915	Errors.Fehler:=TRUE;	915	Errors.Fehler:=TRUE;
916	RETURN;	916	RETURN;
	...		...
943	END dger;	943	END dger;
944		944
945	(=========================== Complex valued routines =====================)	945	(=========================== Complex valued routines =====================)
946		946
947	PROCEDURE zswap( N : CARDINAL;	947	PROCEDURE zswap( N : CARDINAL;
948	VAR X : ARRAY OF CFLOAT;	948	VAR X : ARRAY OF COMPLEX16;
949	IncX : INTEGER;	949	IncX : INTEGER;
950	VAR Y : ARRAY OF CFLOAT;	950	VAR Y : ARRAY OF COMPLEX16;
951	IncY : INTEGER);	951	IncY : INTEGER);
952		952
953	CONST veclen = 4;	953	CONST veclen = 4;
954		954
955	VAR Xi : CFLOAT;	955	VAR Xi : COMPLEX16;
956	Xtmp : ARRAY [0..veclen-1] OF CFLOAT;	956	Xtmp : ARRAY [0..veclen-1] OF COMPLEX16;
957	ix,iy : INTEGER;	957	ix,iy : INTEGER;
958	i,m : CARDINAL;	958	i,m : CARDINAL;
959	BEGIN	959	BEGIN
960	IF (N = 0) THEN RETURN; END;	960	IF (N = 0) THEN RETURN; END;
961		961
	...		...
984	END;	984	END;
985	END;	985	END;
986	END zswap;	986	END zswap;
987		987
988	PROCEDURE zcopy( N : INTEGER;	988	PROCEDURE zcopy( N : INTEGER;
989	VAR X : ARRAY OF CFLOAT;	989	VAR X : ARRAY OF COMPLEX16;
990	IncX : INTEGER;	990	IncX : INTEGER;
991	VAR Y : ARRAY OF CFLOAT;	991	VAR Y : ARRAY OF COMPLEX16;
992	IncY : INTEGER);	992	IncY : INTEGER);
993		993
994	(----------------------------------------------------------------)	994	(----------------------------------------------------------------)
995	(* Adopted to Modula-2, MRi, 04.04.2016, complex version 09.08.18 *)	995	(* Adopted to Modula-2, MRi, 04.04.2016, complex version 09.08.18 *)
996	(----------------------------------------------------------------)	996	(----------------------------------------------------------------)
	...		...
1022	END;	1022	END;
1023	END;	1023	END;
1024	END zcopy;	1024	END zcopy;
1025		1025
1026	PROCEDURE zdotc( N : INTEGER;	1026	PROCEDURE zdotc( N : INTEGER;
1027	VAR X : ARRAY OF CFLOAT;	1027	VAR X : ARRAY OF COMPLEX16;
1028	IncX : INTEGER;	1028	IncX : INTEGER;
1029	VAR Y : ARRAY OF CFLOAT;	1029	VAR Y : ARRAY OF COMPLEX16;
1030	IncY : INTEGER) : CFLOAT;	1030	IncY : INTEGER) : COMPLEX16;
1031		1031
1032	CONST veclen = 4;	1032	CONST veclen = 4;
1033		1033
1034	VAR dtemp : CFLOAT;	1034	VAR dtemp : COMPLEX16;
1035	i,ix,iy,m : INTEGER;	1035	i,ix,iy,m : INTEGER;
1036	BEGIN	1036	BEGIN
1037	IF (N <= 0) THEN RETURN zero; END;	1037	IF (N <= 0) THEN RETURN zero; END;
1038		1038
1039	dtemp := zero;	1039	dtemp := zero;
	...		...
1042	m := N MOD veclen;	1042	m := N MOD veclen;
1043	IF (m # 0) THEN	1043	IF (m # 0) THEN
1044	FOR i:=0 TO m-1 DO (* clean-up loop *)	1044	FOR i:=0 TO m-1 DO (* clean-up loop *)
1045	dtemp:=dtemp + conj(X[i])*Y[i];	1045	dtemp:=dtemp + conj(X[i])*Y[i];
1046	END;	1046	END;
		1047	(****
1047	IF (N < veclen) THEN RETURN dtemp; END;	1048	IF (N < veclen) THEN RETURN dtemp; END;
		1049	***)
		1050
1048	END;	1051	END;
1049	(* i := m - veclen; *)	1052	(* i := m - veclen; *)
1050	FOR i:=m TO N-1 BY veclen DO	1053	FOR i:=m TO N-1 BY veclen DO
1051	dtemp:=dtemp + conj(X[i+0])Y[i+0] + conj(X[i+1])Y[i+1] +	1054	dtemp:=dtemp + conj(X[i+0])Y[i+0] + conj(X[i+1])Y[i+1] +
1052	conj(X[i+2])Y[i+2] + conj(X[i+3])Y[i+3];	1055	conj(X[i+2])Y[i+2] + conj(X[i+3])Y[i+3];
	...		...
1062	END;	1065	END;
1063	RETURN dtemp;	1066	RETURN dtemp;
1064	END zdotc;	1067	END zdotc;
1065		1068
1066	PROCEDURE dznrm2( N : INTEGER;	1069	PROCEDURE dznrm2( N : INTEGER;
1067	VAR X : ARRAY OF CFLOAT;	1070	VAR X : ARRAY OF COMPLEX16;
1068	IncX : INTEGER) : FLOAT;	1071	IncX : INTEGER) : REAL8;
1069		1072
1070	VAR norm,scale,ssq,tmp : FLOAT;	1073	VAR norm,scale,ssq,tmp : REAL8;
1071	i,ix : INTEGER;	1074	i,ix : INTEGER;
1072	zw : FLOAT;	1075	zw : REAL8;
1073	BEGIN	1076	BEGIN
1074	IF (N < 1) OR (IncX < 1) THEN	1077	IF (N < 1) OR (IncX < 1) THEN
1075	norm := 0.0;	1078	norm := 0.0;
1076	ELSE	1079	ELSE
1077	scale := 0.0;	1080	scale := 0.0;
	...		...
1106	END; (* IF *)	1109	END; (* IF *)
1107	RETURN norm;	1110	RETURN norm;
1108	END dznrm2;	1111	END dznrm2;
1109		1112
1110	PROCEDURE zscal( n : INTEGER;	1113	PROCEDURE zscal( n : INTEGER;
1111	da : CFLOAT;	1114	da : COMPLEX16;
1112	VAR X : ARRAY OF CFLOAT;	1115	VAR X : ARRAY OF COMPLEX16;
1113	IncX : INTEGER);	1116	IncX : INTEGER);
1114		1117
1115	CONST veclen = 4;	1118	CONST veclen = 4;
1116		1119
1117	VAR i,m : INTEGER;	1120	VAR i,m : INTEGER;
	...		...
1138	END;	1141	END;
1139	END;	1142	END;
1140	END zscal;	1143	END zscal;
1141		1144
1142	PROCEDURE zaxpy( n : INTEGER;	1145	PROCEDURE zaxpy( n : INTEGER;
1143	da : CFLOAT;	1146	da : COMPLEX16;
1144	VAR X : ARRAY OF CFLOAT;	1147	VAR X : ARRAY OF COMPLEX16;
1145	IncX : INTEGER;	1148	IncX : INTEGER;
1146	VAR Y : ARRAY OF CFLOAT;	1149	VAR Y : ARRAY OF COMPLEX16;
1147	IncY : INTEGER);	1150	IncY : INTEGER);
1148		1151