|
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 LFLOAT(#) durch VAL(LONGREAL,#) *)
|
14 |
(* 13.10.15, MRi: Ersetzen von FLOAT(#) durch VAL(FLOAT,#) *)
|
| 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". *)
|
|
... |
|
... |
| 24 |
(* 31.03.16, MRi: drot und drotg eingefuegt *)
|
24 |
(* 31.03.16, MRi: drot und drotg eingefuegt *)
|
| 25 |
(* 13.04.16, MRi: dcopy eingefuegt *)
|
25 |
(* 13.04.16, MRi: dcopy eingefuegt *)
|
| 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 *)
|
|
|
30 |
(* 12.09.18, MRi: zswap,zcopy,zdotc,dznrm2,zscal,zaxpy,zdrot eingefuegt *)
|
| 29 |
(*------------------------------------------------------------------------*)
|
31 |
(*------------------------------------------------------------------------*)
|
| 30 |
(* Testroutinen *)
|
32 |
(* Testroutinen *)
|
| 31 |
(* *)
|
33 |
(* *)
|
|
|
34 |
(* - dgemv in TstDGEMV *)
|
|
|
35 |
(* - zgemv in TstZGEMV *)
|
| 32 |
(* - zgemm in TstCmplxMaMul *)
|
36 |
(* - zgemm in TstCmplxMaMul *)
|
| 33 |
(*------------------------------------------------------------------------*)
|
37 |
(*------------------------------------------------------------------------*)
|
| 34 |
(* Offene Punkte: *)
|
38 |
(* Offene Punkte: *)
|
| 35 |
(* *)
|
39 |
(* *)
|
| 36 |
(* - Weiteres austesten von dgemm und besonders dgemv *)
|
40 |
(* - Weiteres austesten von dgemm und besonders dgemv *)
|
|
... |
|
... |
| 38 |
(*------------------------------------------------------------------------*)
|
42 |
(*------------------------------------------------------------------------*)
|
| 39 |
(* Implementation : Michael Riedl *)
|
43 |
(* Implementation : Michael Riedl *)
|
| 40 |
(* Licence : GNU Lesser General Public License (LGPL) *)
|
44 |
(* Licence : GNU Lesser General Public License (LGPL) *)
|
| 41 |
(*------------------------------------------------------------------------*)
|
45 |
(*------------------------------------------------------------------------*)
|
| 42 |
|
46 |
|
| 43 |
(* $Id: LibDBlasM2.mod,v 1.9 2017/10/29 09:55:11 mriedl Exp mriedl $ *)
|
47 |
(* $Id: LibDBlasM2.mod,v 1.11 2018/09/12 13:20:49 mriedl Exp mriedl $ *)
|
| 44 |
|
48 |
|
| 45 |
IMPORT SYSTEM;
|
49 |
IMPORT SYSTEM;
|
| 46 |
FROM Deklera IMPORT FLOAT; (* REAL type *)
|
50 |
FROM Deklera IMPORT FLOAT,CFLOAT; (* REAL/COMPLEX Type *)
|
| 47 |
FROM LongMath IMPORT sqrt;
|
51 |
FROM LongMath IMPORT sqrt;
|
| 48 |
FROM LongComplexMath IMPORT conj;
|
52 |
FROM LongComplexMath IMPORT zero,one,conj,scalarMult;
|
| 49 |
IMPORT Errors;
|
53 |
IMPORT Errors;
|
|
|
54 |
FROM F77func IMPORT MAX0;
|
| 50 |
|
55 |
|
| 51 |
IMPORT TIO;
|
|
|
| 52 |
|
|
|
| 53 |
TYPE PVEKTOR = POINTER TO ARRAY [0..MAX(INTEGER)-1] OF FLOAT;
|
|
|
| 54 |
|
56 |
|
| 55 |
PROCEDURE IMax(a,b : INTEGER) : INTEGER;
|
57 |
PROCEDURE IMax(a,b : INTEGER) : INTEGER;
|
| 56 |
|
58 |
|
| 57 |
BEGIN
|
59 |
BEGIN
|
| 58 |
IF (a > b) THEN RETURN a; ELSE RETURN b; END;
|
60 |
IF (a > b) THEN RETURN a; ELSE RETURN b; END;
|
|
... |
|
... |
| 182 |
ELSE
|
184 |
ELSE
|
| 183 |
scale := 0.0;
|
185 |
scale := 0.0;
|
| 184 |
ssq := 1.0;
|
186 |
ssq := 1.0;
|
| 185 |
ix := 0;
|
187 |
ix := 0;
|
| 186 |
FOR i:=0 TO n-1 DO
|
188 |
FOR i:=0 TO n-1 DO
|
| 187 |
TIO.WrStr("ix = "); TIO.WrInt(ix,1); TIO.WrLn;
|
|
|
| 188 |
IF (X[ix] # 0.0) THEN
|
189 |
IF (X[ix] # 0.0) THEN
|
| 189 |
absxi := ABS(X[ix]);
|
190 |
absxi := ABS(X[ix]);
|
| 190 |
IF (scale < absxi) THEN
|
191 |
IF (scale < absxi) THEN
|
| 191 |
zw := scale / absxi;
|
192 |
zw := scale / absxi;
|
| 192 |
ssq := 1.0 + ssq*zw*zw;
|
193 |
ssq := 1.0 + ssq*zw*zw;
|
|
... |
|
... |
| 528 |
(* Finds the index of element having min. absolute value. *)
|
529 |
(* Finds the index of element having min. absolute value. *)
|
| 529 |
(*----------------------------------------------------------------*)
|
530 |
(*----------------------------------------------------------------*)
|
| 530 |
|
531 |
|
| 531 |
VAR dmin : FLOAT;
|
532 |
VAR dmin : FLOAT;
|
| 532 |
i,ix,itemp : INTEGER;
|
533 |
i,ix,itemp : INTEGER;
|
| 533 |
XX : PVEKTOR;
|
|
|
| 534 |
BEGIN
|
534 |
BEGIN
|
| 535 |
IF (n < 1) THEN RETURN -1; END;
|
535 |
IF (n < 1) THEN RETURN -1; END;
|
| 536 |
IF (n = 1) THEN RETURN 0; END;
|
536 |
IF (n = 1) THEN RETURN 0; END;
|
| 537 |
|
537 |
|
| 538 |
XX:=SYSTEM.CAST(PVEKTOR,SYSTEM.ADR(X));
|
|
|
| 539 |
|
|
|
| 540 |
itemp:=0;
|
538 |
itemp:=0;
|
| 541 |
IF (IncX <> 1) THEN (* code for Increment not equal to 1 *)
|
539 |
IF (IncX <> 1) THEN (* code for Increment not equal to 1 *)
|
| 542 |
dmin := ABS(XX^[0]);
|
540 |
dmin := ABS(X[0]);
|
| 543 |
ix := 1 + IncX;
|
541 |
ix := 1 + IncX;
|
| 544 |
FOR i:=1 TO n-1 DO
|
542 |
FOR i:=1 TO n-1 DO
|
| 545 |
IF (ABS(XX^[ix]) < dmin) THEN
|
543 |
IF (ABS(X[ix]) < dmin) THEN
|
| 546 |
itemp := i;
|
544 |
itemp := i;
|
| 547 |
dmin := ABS(XX^[ix]);
|
545 |
dmin := ABS(X[ix]);
|
| 548 |
END;
|
546 |
END;
|
| 549 |
INC(ix,IncX);
|
547 |
INC(ix,IncX);
|
| 550 |
END;
|
548 |
END;
|
| 551 |
ELSE (* code for increment equal to 1 *)
|
549 |
ELSE (* code for increment equal to 1 *)
|
| 552 |
dmin := ABS(XX^[0]);
|
550 |
dmin := ABS(X[0]);
|
| 553 |
FOR i:=1 TO n-1 DO
|
551 |
FOR i:=1 TO n-1 DO
|
| 554 |
IF (ABS(XX^[i]) < dmin) THEN
|
552 |
IF (ABS(X[i]) < dmin) THEN
|
| 555 |
itemp := i;
|
553 |
itemp := i;
|
| 556 |
dmin := ABS(XX^[i]);
|
554 |
dmin := ABS(X[i]);
|
| 557 |
END;
|
555 |
END;
|
| 558 |
END;
|
556 |
END;
|
| 559 |
END;
|
557 |
END;
|
| 560 |
RETURN itemp + 1;
|
558 |
RETURN itemp + 1;
|
| 561 |
END idamin;
|
559 |
END idamin;
|
|
... |
|
... |
| 616 |
VAR Temp : FLOAT;
|
614 |
VAR Temp : FLOAT;
|
| 617 |
i,j,iy,jx,jy,kx,ky,Info,LenX,LenY : INTEGER;
|
615 |
i,j,iy,jx,jy,kx,ky,Info,LenX,LenY : INTEGER;
|
| 618 |
BEGIN (* Testen *)
|
616 |
BEGIN (* Testen *)
|
| 619 |
(* Test the input parameters. *)
|
617 |
(* Test the input parameters. *)
|
| 620 |
|
618 |
|
| 621 |
IF (CAP(Trans) = 'C') THEN
|
619 |
IF (NOT (CAP(Trans) = 'N')) AND (NOT (CAP(Trans) = 'T')) THEN
|
| 622 |
Info := 1;
|
620 |
Info := 1;
|
| 623 |
ELSIF (M < 0) THEN
|
621 |
ELSIF (M < 0) THEN
|
| 624 |
Info := 2;
|
622 |
Info := 2;
|
| 625 |
ELSIF (N < 0) THEN
|
623 |
ELSIF (N < 0) THEN
|
| 626 |
Info := 3;
|
624 |
Info := 3;
|
|
... |
|
... |
| 669 |
ky := 1 - (LenY - 1)*IncY;
|
667 |
ky := 1 - (LenY - 1)*IncY;
|
| 670 |
END;
|
668 |
END;
|
| 671 |
|
669 |
|
| 672 |
(* Start the operations. In this version the elements of A are *)
|
670 |
(* Start the operations. In this version the elements of A are *)
|
| 673 |
(* accessed sequentially with one pass through A. *)
|
671 |
(* accessed sequentially with one pass through A. *)
|
| 674 |
|
|
|
| 675 |
|
672 |
|
| 676 |
IF (Beta # 1.0) THEN (* First form y := beta*y. *)
|
673 |
IF (Beta # 1.0) THEN (* First form y := beta*y. *)
|
| 677 |
IF (IncY = 1) THEN
|
674 |
IF (IncY = 1) THEN
|
| 678 |
IF (Beta = 0.0) THEN
|
675 |
IF (Beta = 0.0) THEN
|
| 679 |
FOR i:=0 TO LenY-1 DO Y[i]:=0.0; END;
|
676 |
FOR i:=0 TO LenY-1 DO Y[i]:=0.0; END;
|
|
... |
|
... |
| 943 |
INC(ix,IncX);
|
940 |
INC(ix,IncX);
|
| 944 |
END;
|
941 |
END;
|
| 945 |
END;
|
942 |
END;
|
| 946 |
END dger;
|
943 |
END dger;
|
| 947 |
|
944 |
|
|
|
945 |
(*=========================== Complex valued routines =====================*)
|
|
|
946 |
|
|
|
947 |
PROCEDURE zswap( N : CARDINAL;
|
|
|
948 |
VAR X : ARRAY OF CFLOAT;
|
|
|
949 |
IncX : INTEGER;
|
|
|
950 |
VAR Y : ARRAY OF CFLOAT;
|
|
|
951 |
IncY : INTEGER);
|
|
|
952 |
|
|
|
953 |
CONST veclen = 4;
|
|
|
954 |
|
|
|
955 |
VAR Xi : CFLOAT;
|
|
|
956 |
Xtmp : ARRAY [0..veclen-1] OF CFLOAT;
|
|
|
957 |
ix,iy : INTEGER;
|
|
|
958 |
i,m : CARDINAL;
|
|
|
959 |
BEGIN
|
|
|
960 |
IF (N = 0) THEN RETURN; END;
|
|
|
961 |
|
|
|
962 |
IF (IncX = 1) AND (IncY = 1) THEN
|
|
|
963 |
m := N MOD veclen;
|
|
|
964 |
IF (m # 0) THEN
|
|
|
965 |
FOR i:=0 TO m-1 DO (* Clean up loop *)
|
|
|
966 |
Xi := X[i]; X[i] := Y[i]; Y[i] := Xi;
|
|
|
967 |
END;
|
|
|
968 |
END;
|
|
|
969 |
IF (N < veclen) THEN RETURN; END;
|
|
|
970 |
FOR i:=m TO N-1 BY veclen DO
|
|
|
971 |
Xtmp[0] := X[i+0]; Xtmp[1] := X[i+1];
|
|
|
972 |
Xtmp[2] := X[i+2]; Xtmp[3] := X[i+3];
|
|
|
973 |
X[i+0] := Y[i+0]; X[i+1] := Y[i+1];
|
|
|
974 |
X[i+2] := Y[i+2]; X[i+3] := Y[i+3];
|
|
|
975 |
Y[i+1] := Xtmp[1]; Y[i+0] := Xtmp[0];
|
|
|
976 |
Y[i+2] := Xtmp[2]; Y[i+3] := Xtmp[3];
|
|
|
977 |
END;
|
|
|
978 |
ELSE (* code for unequal increments or equal increments # 1 *)
|
|
|
979 |
IF (IncX > 0) THEN ix:=0; ELSE ix:=(1 - VAL(INTEGER,N))*IncX; END;
|
|
|
980 |
IF (IncY > 0) THEN iy:=0; ELSE iy:=(1 - VAL(INTEGER,N))*IncY; END;
|
|
|
981 |
FOR i:=0 TO N-1 DO
|
|
|
982 |
Xi := X[ix]; X[ix] := Y[iy]; Y[iy] := Xi;
|
|
|
983 |
INC(ix,IncX); INC(iy,IncY);
|
|
|
984 |
END;
|
|
|
985 |
END;
|
|
|
986 |
END zswap;
|
|
|
987 |
|
|
|
988 |
PROCEDURE zcopy( N : INTEGER;
|
|
|
989 |
VAR X : ARRAY OF CFLOAT;
|
|
|
990 |
IncX : INTEGER;
|
|
|
991 |
VAR Y : ARRAY OF CFLOAT;
|
|
|
992 |
IncY : INTEGER);
|
|
|
993 |
|
|
|
994 |
(*----------------------------------------------------------------*)
|
|
|
995 |
(* Adopted to Modula-2, MRi, 04.04.2016, complex version 09.08.18 *)
|
|
|
996 |
(*----------------------------------------------------------------*)
|
|
|
997 |
|
|
|
998 |
VAR i,ix,iy,m : INTEGER;
|
|
|
999 |
BEGIN
|
|
|
1000 |
IF (N <= 0) THEN RETURN; END;
|
|
|
1001 |
|
|
|
1002 |
IF (IncX = 1) AND (IncY = 1) THEN
|
|
|
1003 |
(* code for both increments equal to 1 *)
|
|
|
1004 |
m := (N MOD 8);
|
|
|
1005 |
IF (m # 0) THEN (* Clean-up loop *)
|
|
|
1006 |
FOR i:=0 TO m-1 DO Y[i] := X[i]; END;
|
|
|
1007 |
IF (N < 8) THEN RETURN; END
|
|
|
1008 |
END;
|
|
|
1009 |
FOR i:=m TO N-1 BY 8 DO
|
|
|
1010 |
Y[i+0] := X[i+0]; Y[i+1] := X[i+1];
|
|
|
1011 |
Y[i+2] := X[i+2]; Y[i+3] := X[i+3];
|
|
|
1012 |
Y[i+4] := X[i+4]; Y[i+5] := X[i+5];
|
|
|
1013 |
Y[i+6] := X[i+6]; Y[i+7] := X[i+7];
|
|
|
1014 |
END;
|
|
|
1015 |
ELSE
|
|
|
1016 |
(* code for unequal increments or equal increments not equal to 1 *)
|
|
|
1017 |
IF (IncX > 0) THEN ix:=0; ELSE ix:=(1 - VAL(INTEGER,N))*IncX; END;
|
|
|
1018 |
IF (IncY > 0) THEN iy:=0; ELSE iy:=(1 - VAL(INTEGER,N))*IncY; END;
|
|
|
1019 |
FOR i:=0 TO N-1 DO
|
|
|
1020 |
Y[iy] := X[ix];
|
|
|
1021 |
INC(ix,IncX); INC(iy,IncY);
|
|
|
1022 |
END;
|
|
|
1023 |
END;
|
|
|
1024 |
END zcopy;
|
|
|
1025 |
|
|
|
1026 |
PROCEDURE zdotc( N : INTEGER;
|
|
|
1027 |
VAR X : ARRAY OF CFLOAT;
|
|
|
1028 |
IncX : INTEGER;
|
|
|
1029 |
VAR Y : ARRAY OF CFLOAT;
|
|
|
1030 |
IncY : INTEGER) : CFLOAT;
|
|
|
1031 |
|
|
|
1032 |
CONST veclen = 4;
|
|
|
1033 |
|
|
|
1034 |
VAR dtemp : CFLOAT;
|
|
|
1035 |
i,ix,iy,m : INTEGER;
|
|
|
1036 |
BEGIN
|
|
|
1037 |
IF (N <= 0) THEN RETURN zero; END;
|
|
|
1038 |
|
|
|
1039 |
dtemp := zero;
|
|
|
1040 |
IF (IncX = 1) AND (IncY = 1) THEN
|
|
|
1041 |
(* code for both increments equal to 1 *)
|
|
|
1042 |
m := N MOD veclen;
|
|
|
1043 |
IF (m # 0) THEN
|
|
|
1044 |
FOR i:=0 TO m-1 DO (* clean-up loop *)
|
|
|
1045 |
dtemp:=dtemp + conj(X[i])*Y[i];
|
|
|
1046 |
END;
|
|
|
1047 |
IF (N < veclen) THEN RETURN dtemp; END;
|
|
|
1048 |
END;
|
|
|
1049 |
(* i := m - veclen; *)
|
|
|
1050 |
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] +
|
|
|
1052 |
conj(X[i+2])*Y[i+2] + conj(X[i+3])*Y[i+3];
|
|
|
1053 |
END;
|
|
|
1054 |
ELSE
|
|
|
1055 |
(* code for unequal increments or equal increments not equal to 1 *)
|
|
|
1056 |
ix := 0; IF (IncX < 0) THEN ix := (1-N)*IncX; END;
|
|
|
1057 |
iy := 0; IF (IncY < 0) THEN iy := (1-N)*IncY; END;
|
|
|
1058 |
FOR i:=1 TO N DO
|
|
|
1059 |
dtemp:=dtemp + conj(X[ix])*Y[iy];
|
|
|
1060 |
INC(ix,IncX); INC(iy,IncY);
|
|
|
1061 |
END;
|
|
|
1062 |
END;
|
|
|
1063 |
RETURN dtemp;
|
|
|
1064 |
END zdotc;
|
|
|
1065 |
|
|
|
1066 |
PROCEDURE dznrm2( N : INTEGER;
|
|
|
1067 |
VAR X : ARRAY OF CFLOAT;
|
|
|
1068 |
IncX : INTEGER) : FLOAT;
|
|
|
1069 |
|
|
|
1070 |
VAR norm,scale,ssq,tmp : FLOAT;
|
|
|
1071 |
i,ix : INTEGER;
|
|
|
1072 |
zw : FLOAT;
|
|
|
1073 |
BEGIN
|
|
|
1074 |
IF (N < 1) OR (IncX < 1) THEN
|
|
|
1075 |
norm := 0.0;
|
|
|
1076 |
ELSE
|
|
|
1077 |
scale := 0.0;
|
|
|
1078 |
ssq := 1.0;
|
|
|
1079 |
ix := 0;
|
|
|
1080 |
FOR i:=0 TO N-1 DO (* DO ix=1,1+(N-1)*IncX,IncX *)
|
|
|
1081 |
IF (RE(X[ix]) # 0.0) THEN
|
|
|
1082 |
tmp := ABS(RE(X[ix]));
|
|
|
1083 |
IF (scale < tmp) THEN
|
|
|
1084 |
zw := scale / tmp;
|
|
|
1085 |
ssq := 1.0 + ssq*(zw*zw);
|
|
|
1086 |
scale := tmp;
|
|
|
1087 |
ELSE
|
|
|
1088 |
zw := tmp / scale;
|
|
|
1089 |
ssq:=ssq + (zw*zw);
|
|
|
1090 |
END;
|
|
|
1091 |
END;
|
|
|
1092 |
IF (IM(X[ix]) # 0.0) THEN
|
|
|
1093 |
tmp := ABS(IM(X[ix]));
|
|
|
1094 |
IF (scale < tmp) THEN
|
|
|
1095 |
zw := scale / tmp;
|
|
|
1096 |
ssq := 1.0 + ssq*(zw*zw);
|
|
|
1097 |
scale := tmp;
|
|
|
1098 |
ELSE
|
|
|
1099 |
zw := tmp / scale;
|
|
|
1100 |
ssq:=ssq + (zw*zw);
|
|
|
1101 |
END;
|
|
|
1102 |
END;
|
|
|
1103 |
INC(ix,IncX);
|
|
|
1104 |
END; (* FOR *)
|
|
|
1105 |
norm := scale*sqrt(ssq);
|
|
|
1106 |
END; (* IF *)
|
|
|
1107 |
RETURN norm;
|
|
|
1108 |
END dznrm2;
|
|
|
1109 |
|
|
|
1110 |
PROCEDURE zscal( n : INTEGER;
|
|
|
1111 |
da : CFLOAT;
|
|
|
1112 |
VAR X : ARRAY OF CFLOAT;
|
|
|
1113 |
IncX : INTEGER);
|
|
|
1114 |
|
|
|
1115 |
CONST veclen = 4;
|
|
|
1116 |
|
|
|
1117 |
VAR i,m : INTEGER;
|
|
|
1118 |
BEGIN
|
|
|
1119 |
IF (n <= 0) THEN RETURN; END;
|
|
|
1120 |
|
|
|
1121 |
IF (IncX <> 1) THEN (* code for increment not equal to 1 *)
|
|
|
1122 |
FOR i:=0 TO n-1 DO
|
|
|
1123 |
X[i*IncX] := da*X[i*IncX];
|
|
|
1124 |
END;
|
|
|
1125 |
ELSE (* code for increment equal to 1 *)
|
|
|
1126 |
m := n REM veclen;
|
|
|
1127 |
IF (m <> 0) THEN
|
|
|
1128 |
FOR i:=0 TO m-1 DO X[i] := da*X[i]; END;
|
|
|
1129 |
IF (n < veclen) THEN
|
|
|
1130 |
RETURN;
|
|
|
1131 |
END;
|
|
|
1132 |
END;
|
|
|
1133 |
FOR i:=m TO n-1 BY veclen DO
|
|
|
1134 |
X[i+0] := da*X[i+0];
|
|
|
1135 |
X[i+1] := da*X[i+1];
|
|
|
1136 |
X[i+2] := da*X[i+2];
|
|
|
1137 |
X[i+3] := da*X[i+3];
|
|
|
1138 |
END;
|
|
|
1139 |
END;
|
|
|
1140 |
END zscal;
|
|
|
1141 |
|
|
|
1142 |
PROCEDURE zaxpy( n : INTEGER;
|
|
|
1143 |
da : CFLOAT;
|
|
|
1144 |
VAR X : ARRAY OF CFLOAT;
|
|
|
1145 |
IncX : INTEGER;
|
|
|
1146 |
VAR Y : ARRAY OF CFLOAT;
|
|
|
1147 |
IncY : INTEGER);
|
|
|
1148 |
|
|
|
1149 |
CONST veclen = 4;
|
|
|
1150 |
|
|
|
1151 |
VAR i,ix,iy,m : INTEGER;
|
|
|
1152 |
BEGIN
|
|
|
1153 |
IF (n <= 0 ) THEN RETURN; END;
|
|
|
1154 |
IF (da = 0.0) THEN RETURN; END;
|
|
|
1155 |
|
|
|
1156 |
IF ((IncX <> 1) OR (IncY <> 1)) THEN
|
|
|
1157 |
(* code for unequal increments or equal increments <> 1 *)
|
|
|
1158 |
ix := 1;
|
|
|
1159 |
iy := 1;
|
|
|
1160 |
IF (IncX < 0) THEN ix := (-n+1)*IncX + 1; END;
|
|
|
1161 |
IF (IncY < 0) THEN iy := (-n+1)*IncY + 1; END;
|
|
|
1162 |
FOR i:=0 TO n-1 DO
|
|
|
1163 |
Y[iy] := Y[iy] + da*X[ix];
|
|
|
1164 |
INC(ix, IncX);
|
|
|
1165 |
INC(iy, IncY);
|
|
|
1166 |
END;
|
|
|
1167 |
ELSE (* code for both increments equal to 1 *)
|
|
|
1168 |
m := n REM veclen;
|
|
|
1169 |
IF (m <> 0) THEN
|
|
|
1170 |
FOR i:=0 TO m-1 DO
|
|
|
1171 |
Y[i] := Y[i] + da*X[i];
|
|
|
1172 |
END;
|
|
|
1173 |
IF (n < veclen) THEN RETURN; END;
|
|
|
1174 |
END;
|
|
|
1175 |
FOR i:=m TO n-1 BY veclen DO
|
|
|
1176 |
Y[i] := Y[i] + da*X[i];
|
|
|
1177 |
Y[i+1] := Y[i+1] + da*X[i+1];
|
|
|
1178 |
Y[i+2] := Y[i+2] + da*X[i+2];
|
|
|
1179 |
Y[i+3] := Y[i+3] + da*X[i+3];
|
|
|
1180 |
END;
|
|
|
1181 |
END;
|
|
|
1182 |
END zaxpy;
|
|
|
1183 |
|
|
|
1184 |
PROCEDURE zdrot( N : INTEGER;
|
|
|
1185 |
VAR X : ARRAY OF CFLOAT;
|
|
|
1186 |
IncX : INTEGER;
|
|
|
1187 |
VAR Y : ARRAY OF CFLOAT;
|
|
|
1188 |
IncY : INTEGER;
|
|
|
1189 |
c,s : FLOAT);
|
|
|
1190 |
|
|
|
1191 |
VAR i,ix,iy : INTEGER;
|
|
|
1192 |
tmp : CFLOAT;
|
|
|
1193 |
BEGIN
|
|
|
1194 |
IF (IncX = 1) AND (IncY = 1) THEN
|
|
|
1195 |
FOR i:=0 TO N-1 DO
|
|
|
1196 |
tmp := scalarMult( c,X[i]) + scalarMult(s,Y[i]);
|
|
|
1197 |
Y[i] := scalarMult(-s,X[i]) + scalarMult(c,Y[i]);
|
|
|
1198 |
X[i] := tmp;
|
|
|
1199 |
END;
|
|
|
1200 |
ELSE
|
|
|
1201 |
IF (IncX > 0) THEN ix:=0; ELSE ix:=(1 - VAL(INTEGER,N))*IncX; END;
|
|
|
1202 |
IF (IncY > 0) THEN iy:=0; ELSE iy:=(1 - VAL(INTEGER,N))*IncY; END;
|
|
|
1203 |
FOR i:=0 TO N-1 DO
|
|
|
1204 |
tmp := scalarMult( c,X[ix]) + scalarMult(s,Y[iy]);
|
|
|
1205 |
Y[iy] := scalarMult(-s,X[ix]) + scalarMult(c,Y[iy]);
|
|
|
1206 |
X[iy] := tmp;
|
|
|
1207 |
ix := ix + IncX;
|
|
|
1208 |
iy := iy + IncY;
|
|
|
1209 |
END;
|
|
|
1210 |
END;
|
|
|
1211 |
END zdrot;
|
|
|
1212 |
|
|
|
1213 |
PROCEDURE zgemv( Trans : CHAR;
|
|
|
1214 |
M,N : INTEGER;
|
|
|
1215 |
Alpha : CFLOAT;
|
|
|
1216 |
VAR A : ARRAY OF ARRAY OF CFLOAT;
|
|
|
1217 |
lda : INTEGER;
|
|
|
1218 |
VAR X : ARRAY OF CFLOAT;
|
|
|
1219 |
IncX : INTEGER;
|
|
|
1220 |
Beta : CFLOAT;
|
|
|
1221 |
VAR Y : ARRAY OF CFLOAT;
|
|
|
1222 |
IncY : INTEGER);
|
|
|
1223 |
|
|
|
1224 |
VAR i,j,iy,jx,kx,ky : INTEGER;
|
|
|
1225 |
lenx,leny,info : INTEGER;
|
|
|
1226 |
s : CFLOAT;
|
|
|
1227 |
noconj : BOOLEAN;
|
|
|
1228 |
BEGIN
|
|
|
1229 |
(* test the input parameters *)
|
|
|
1230 |
|
|
|
1231 |
info := 0;
|
|
|
1232 |
IF (NOT (CAP(Trans) = 'N')) AND (NOT (CAP(Trans) = 'T')) AND
|
|
|
1233 |
(NOT (CAP(Trans) = 'C'))
|
|
|
1234 |
THEN
|
|
|
1235 |
info := 1;
|
|
|
1236 |
ELSIF (M < 0) THEN
|
|
|
1237 |
info := 2;
|
|
|
1238 |
ELSIF (N < 0) THEN
|
|
|
1239 |
info := 3;
|
|
|
1240 |
ELSIF (lda < MAX0(1,M)) THEN
|
|
|
1241 |
info := 6;
|
|
|
1242 |
ELSIF (IncX = 0) THEN
|
|
|
1243 |
info := 8;
|
|
|
1244 |
ELSIF (IncY = 0) THEN
|
|
|
1245 |
info := 11;
|
|
|
1246 |
END;
|
|
|
1247 |
IF (info # 0) THEN (* Fehlerbehandlung ! *)
|
|
|
1248 |
Errors.WriteLn; Errors.WriteLn;
|
|
|
1249 |
Errors.WriteString("Fehler in LibDBlas.zgemv, info = ");
|
|
|
1250 |
Errors.WriteInt(info);
|
|
|
1251 |
Errors.WriteLn; Errors.WriteLn;
|
|
|
1252 |
RETURN;
|
|
|
1253 |
END;
|
|
|
1254 |
|
|
|
1255 |
(* quick return if possible *)
|
|
|
1256 |
|
|
|
1257 |
IF ((M = 0) OR (N = 0) OR ((Alpha = zero) AND (Beta = one))) THEN
|
|
|
1258 |
RETURN
|
|
|
1259 |
END;
|
|
|
1260 |
|
|
|
1261 |
noconj := (CAP(Trans) = 'T');
|
|
|
1262 |
|
|
|
1263 |
(* set lenx and leny, the lengths of the vectors X and Y, *)
|
|
|
1264 |
(* and set up the start points in X and Y. *)
|
|
|
1265 |
|
|
|
1266 |
IF (CAP(Trans) = 'N') THEN
|
|
|
1267 |
lenx := N;
|
|
|
1268 |
leny := M;
|
|
|
1269 |
ELSE
|
|
|
1270 |
lenx := M;
|
|
|
1271 |
leny := N;
|
|
|
1272 |
END;
|
|
|
1273 |
|
|
|
1274 |
IF (IncX > 0) THEN
|
|
|
1275 |
kx := 1;
|
|
|
1276 |
ELSE
|
|
|
1277 |
kx := 1 - (lenx-1)*IncX;
|
|
|
1278 |
END;
|
|
|
1279 |
IF (IncY > 0) THEN
|
|
|
1280 |
ky := 1;
|
|
|
1281 |
ELSE
|
|
|
1282 |
ky := 1 - (leny-1)*IncY;
|
|
|
1283 |
END;
|
|
|
1284 |
|
|
|
1285 |
(* start the operations. In this version the elements of A are *)
|
|
|
1286 |
(* accessed sequentially with one pass through A *)
|
|
|
1287 |
|
|
|
1288 |
(* first form Y = Beta*Y *)
|
|
|
1289 |
|
|
|
1290 |
IF (Beta # one) THEN
|
|
|
1291 |
IF (IncY # 1) THEN
|
|
|
1292 |
iy := ky - 1;
|
|
|
1293 |
IF (Beta = zero) THEN
|
|
|
1294 |
FOR i:=0 TO leny-1 DO
|
|
|
1295 |
Y[iy] := zero;
|
|
|
1296 |
INC(iy,IncY);
|
|
|
1297 |
END; (* FOR *)
|
|
|
1298 |
ELSE
|
|
|
1299 |
FOR i:=0 TO leny-1 DO
|
|
|
1300 |
Y[iy] := Beta*Y[iy];
|
|
|
1301 |
INC(iy,IncY);
|
|
|
1302 |
END; (* FOR *)
|
|
|
1303 |
END; (* IF *)
|
|
|
1304 |
ELSE
|
|
|
1305 |
IF (Beta = zero) THEN
|
|
|
1306 |
FOR i:=0 TO leny-1 DO
|
|
|
1307 |
Y[i] := zero;
|
|
|
1308 |
END;
|
|
|
1309 |
ELSE
|
|
|
1310 |
FOR i:=0 TO leny-1 DO
|
|
|
1311 |
Y[i] := Beta*Y[i];
|
|
|
1312 |
END;
|
|
|
1313 |
END; (* IF *)
|
|
|
1314 |
END; (* IF *)
|
|
|
1315 |
END; (* IF *)
|
|
|
1316 |
|
|
|
1317 |
IF (Alpha = zero) THEN RETURN END;
|
|
|
1318 |
|
|
|
1319 |
IF (CAP(Trans) = 'N') THEN (* Form y := alpha*A*x + y. *)
|
|
|
1320 |
|
|
|
1321 |
|
|
|
1322 |
IF (IncX = 1) AND (IncY = 1) THEN
|
|
|
1323 |
FOR i:=0 TO M-1 DO
|
|
|
1324 |
s:=zero;
|
|
|
1325 |
FOR j:=0 TO N-1 DO
|
|
|
1326 |
s:=s + Alpha*A[i,j]*X[j];
|
|
|
1327 |
END;
|
|
|
1328 |
Y[i]:=Y[i] + s;
|
|
|
1329 |
END;
|
|
|
1330 |
ELSE
|
|
|
1331 |
iy:=ky-1;
|
|
|
1332 |
FOR i:=0 TO M-1 DO
|
|
|
1333 |
s:=zero;
|
|
|
1334 |
jx:=kx-1;
|
|
|
1335 |
FOR j:=0 TO N-1 DO
|
|
|
1336 |
s:=s + Alpha*A[i,j]*X[jx]; INC(jx,IncX);
|
|
|
1337 |
END;
|
|
|
1338 |
Y[iy]:=Y[iy] + s; INC(iy,IncY);
|
|
|
1339 |
END;
|
|
|
1340 |
END;
|
|
|
1341 |
|
|
|
1342 |
ELSE
|
|
|
1343 |
|
|
|
1344 |
(* Form Y = Alpha*A'*X + Y or Y = Alpha*conj(A')*x + y. *)
|
|
|
1345 |
|
|
|
1346 |
IF (IncX = 1) AND (IncY = 1) THEN
|
|
|
1347 |
FOR j:=0 TO M-1 DO
|
|
|
1348 |
IF noconj THEN
|
|
|
1349 |
FOR i:=0 TO N-1 DO
|
|
|
1350 |
Y[i]:=Y[i] + Alpha*A[j,i]*X[j];
|
|
|
1351 |
END;
|
|
|
1352 |
ELSE (* (CAP(Trans) = "C") *)
|
|
|
1353 |
FOR i:=0 TO N-1 DO
|
|
|
1354 |
Y[i]:=Y[i] + Alpha*conj(A[j,i])*X[j];
|
|
|
1355 |
END;
|
|
|
1356 |
END;
|
|
|
1357 |
END;
|
|
|
1358 |
ELSE
|
|
|
1359 |
jx:=kx-1;
|
|
|
1360 |
FOR j:=0 TO M-1 DO
|
|
|
1361 |
IF noconj THEN
|
|
|
1362 |
iy:=ky-1;
|
|
|
1363 |
FOR i:=0 TO N-1 DO
|
|
|
1364 |
Y[iy]:=Y[iy] + Alpha*A[j,i]*X[jx]; INC(iy,IncY);
|
|
|
1365 |
END;
|
|
|
1366 |
ELSE (* (CAP(Trans) = "C") *)
|
|
|
1367 |
iy:=ky-1;
|
|
|
1368 |
FOR i:=0 TO N-1 DO
|
|
|
1369 |
Y[iy]:=Y[iy] + Alpha*conj(A[j,i])*X[jx]; INC(iy,IncY);
|
|
|
1370 |
END;
|
|
|
1371 |
END;
|
|
|
1372 |
INC(jx,IncX);
|
|
|
1373 |
END; (* FOR j *)
|
|
|
1374 |
END;
|
|
|
1375 |
END; (* IF Trans *)
|
|
|
1376 |
END zgemv;
|
|
|
1377 |
|
| 948 |
PROCEDURE zgemm( TransA,TransB : CHAR;
|
1378 |
PROCEDURE zgemm( TransA,TransB : CHAR;
|
| 949 |
M,N,K : INTEGER;
|
1379 |
M,N,K : INTEGER;
|
| 950 |
Alpha : LONGCOMPLEX;
|
1380 |
Alpha : CFLOAT;
|
| 951 |
VAR A : ARRAY OF ARRAY OF LONGCOMPLEX;
|
1381 |
VAR A : ARRAY OF ARRAY OF CFLOAT;
|
| 952 |
LDA : INTEGER;
|
1382 |
LDA : INTEGER;
|
| 953 |
VAR B : ARRAY OF ARRAY OF LONGCOMPLEX;
|
1383 |
VAR B : ARRAY OF ARRAY OF CFLOAT;
|
| 954 |
LDB : INTEGER;
|
1384 |
LDB : INTEGER;
|
| 955 |
Beta : LONGCOMPLEX;
|
1385 |
Beta : CFLOAT;
|
| 956 |
VAR C : ARRAY OF ARRAY OF LONGCOMPLEX;
|
1386 |
VAR C : ARRAY OF ARRAY OF CFLOAT;
|
| 957 |
LDC : INTEGER);
|
1387 |
LDC : INTEGER);
|
| 958 |
|
1388 |
|
| 959 |
CONST zero = CMPLX(0.0,0.0);
|
1389 |
CONST zero = CMPLX(0.0,0.0);
|
| 960 |
one = CMPLX(1.0,0.0);
|
1390 |
one = CMPLX(1.0,0.0);
|
| 961 |
|
1391 |
|
| 962 |
VAR NotA,NotB,ConjA,ConjB : BOOLEAN;
|
1392 |
VAR NotA,NotB,ConjA,ConjB : BOOLEAN;
|
| 963 |
i,j,k,Info,NRowA,NRowB : INTEGER;
|
1393 |
i,j,k,Info,NRowA,NRowB : INTEGER;
|
| 964 |
Temp,Aki,Aik : LONGCOMPLEX;
|
1394 |
Temp,Aki,Aik : CFLOAT;
|
| 965 |
Nm1,Mm1,Km1 : INTEGER;
|
1395 |
Nm1,Mm1,Km1 : INTEGER;
|
| 966 |
BEGIN
|
1396 |
BEGIN
|
| 967 |
(* Set NotA and NotB as true if A and B respectively are not *)
|
1397 |
(* Set NotA and NotB as true if A and B respectively are not *)
|
| 968 |
(* transposed and set NRowA, NColA and NRowB as the number of rows *)
|
1398 |
(* transposed and set NRowA, NColA and NRowB as the number of rows *)
|
| 969 |
(* and columns of A and the number of rows of B respectively. *)
|
1399 |
(* and columns of A and the number of rows of B respectively. *)
|