DEFINITION MODULE LibDBlasM2;
(*------------------------------------------------------------------------*)
(* Basic Linear Algebra Prozeduren in Modula-2 mit M2-konformen *)
(* Aufrufen der Parameter. *)
(* Some level 1,2 & 3 BLAS routines implemented in Modula-2 pathing the *)
(* arguments in a "M2 manner". *)
(* *)
(* These routine use "normal" M2 mechanism for handing over open array *)
(* parameters. In Fortran and "C" sometime a array element is passed to *)
(* an open array field which is not how it works Modula-2 and Oberon. *)
(* For that purpose there is a Module "LibDBlas" (without the M2 at the *)
(* end) which is doing some address translation to be used when trans- *)
(* lating Fortran code to Modula-2. *)
(* Furthermore there is a "foreign" definition module LibDBlasF77 for *)
(* directly call the Fortran versions (only level 1 BLAS) - but here all *)
(* arguments had to be declared as "VAR" to confirm with Fortran calling *)
(* conventions. *)
(* *)
(* Routine provided so far are: *)
(* *)
(* dnrm2 : Euklidian Norm of a vector *)
(* dswap : swap two vectors *)
(* dcopy : copy a vector to a vector *)
(* drot : plane rotation *)
(* drotg : construct givens rotation *)
(* dscal : scale vector by a constant *)
(* daxpy : constant times a vector plus a vector *)
(* ddot : dot product of two vectors *)
(* idamax : index of vector element with largest absolute value *)
(* idamin : index of vector element with smallest absolute value *)
(* dasum : sum of vector elements *)
(* dgemv : matrix vector operations *)
(* dgemm : matrix matrix operations *)
(* dger : dyadic product of two vectors *)
(* dgemv : matrix vector operations (complex) *)
(* zgemm : matrix matrix operations (complex) *)
(* *)
(* Additional Routines *)
(* *)
(* SumVek : Sum vector elements *)
(* AbsSumVek : Sum of absolute values of vector element *)
(* ENorm : Euklidian norm of a vector *)
(*------------------------------------------------------------------------*)
(* $Id: LibDBlasM2.def,v 1.12 2018/09/12 13:20:49 mriedl Exp mriedl $ *)
FROM Deklera IMPORT FLOAT,CFLOAT; (* REAL/COMPLEX Type *)
PROCEDURE SumVek(VAR X : ARRAY OF FLOAT;
s,e : CARDINAL) : FLOAT;
(*----------------------------------------------------------------*)
(* Berechnet \sum_{i=s}^e X_i, wobei der erste (m"ogliche) *)
(* Index 1 ist. *)
(* *)
(* Summation des Vektors X nach Kahans Summationsformel. *)
(* Sehr pr"aziser, aber aufwendiger Summationsalgorithmus. *)
(* Vorsicht bei einigen Optimierer, es kann sein, da3 das berech- *)
(* nete "gard-digit" (c) wegoptimiert wird !!! *)
(*----------------------------------------------------------------*)
PROCEDURE AbsSumVek(VAR X : ARRAY OF FLOAT;
s,e : CARDINAL) : FLOAT;
(*----------------------------------------------------------------*)
(* Berechnet \sum_{i=s}^e |X_i|, wobei der erste (m"ogliche) *)
(* Index 1 ist. *)
(* Bis auf die Absolutwerte identisch mit SumVek. *)
(*----------------------------------------------------------------*)
PROCEDURE ENorm(VAR X : ARRAY OF FLOAT;
s,e : CARDINAL) : FLOAT;
(*----------------------------------------------------------------*)
(* Berechnet \sqrt{ \sum_{i=s}^e X_i^2} , wobei der erste *)
(* (m"ogliche) Index s=1 ist. *)
(* *)
(* The euclidian morm is computed by accumulating the sum of *)
(* squares in three different sums. The sums of squares for the *)
(* small and large components are scaled so that no overflows *)
(* occur. Non-destructive underflows are permitted. Underflows *)
(* and overflows do not occur in the computation of the unscaled *)
(* sum of squares for the intermediate components. *)
(* The definitions of small, intermediate and large components *)
(* depend on two constants, MinLR and MaxLR. The main *)
(* restrictions on these constants are that MinLR**2 not *)
(* underflow ans MaxLR**2 not overflow. The constants *)
(* given here are suitable for every known computer. *)
(* *)
(* Argonne national laboratory. MINPACK project. March 1980. *)
(* Burton S. Garbow, Kenneth E. Hillstrom, Jorge J. More *)
(* *)
(* Anpassung an Modula-2 und "Uberarbeitung : M.Riedl, 23.11.93. *)
(* Berechnet \sqrt{ \sum_{i=s}^e X_i^2} *)
(*----------------------------------------------------------------*)
PROCEDURE dnrm2( n : INTEGER;
VAR X : ARRAY OF FLOAT;
IncX : INTEGER): FLOAT;
(*----------------------------------------------------------------*)
(* dnrm2 returns the euclidean norm of a vector via the function *)
(* name, so that dnrm2 := sqrt( x'*x ) *)
(*----------------------------------------------------------------*)
PROCEDURE dswap( N : CARDINAL;
VAR X : ARRAY OF FLOAT;
incX : INTEGER;
VAR Y : ARRAY OF FLOAT;
incY : INTEGER);
(*----------------------------------------------------------------*)
(* Interchanges two vectors using unrolled loops for increments *)
(* equal to 1. *)
(*----------------------------------------------------------------*)
PROCEDURE dcopy( N : INTEGER;
VAR X : ARRAY OF FLOAT;
IncX : INTEGER;
VAR Y : ARRAY OF FLOAT;
IncY : INTEGER);
(*----------------------------------------------------------------*)
(* copies a vector, X, to a vector, Y. *)
(* uses unrolled loops for increments equal to one. *)
(* jack dongarra, linpack, 3/11/78. *)
(* MRi, Modula-2 10.04.16 *)
(*----------------------------------------------------------------*)
PROCEDURE drot( N : INTEGER;
VAR X : ARRAY OF FLOAT;
incX : INTEGER;
VAR Y : ARRAY OF FLOAT;
incY : INTEGER;
c,s : FLOAT);
(*---------------------------------------------------------------*)
(* Applies a plane rotation. *)
(* Jack Dongarra, linpack, 3/11/78. *)
(*---------------------------------------------------------------*)
PROCEDURE drotg(VAR da : FLOAT;
VAR db : FLOAT;
VAR c : FLOAT;
VAR s : FLOAT);
(*---------------------------------------------------------------*)
(* Construct a Givens plane rotation *)
(* Jack Dongarra, linpack, 3/11/78. *)
(*---------------------------------------------------------------*)
PROCEDURE dscal( n : INTEGER;
da : FLOAT;
VAR dx : ARRAY OF FLOAT;
incx : INTEGER);
(*----------------------------------------------------------------*)
(* Scales a vector by a constant (UNROLLED version) *)
(* *)
(* Jack Dongarra, linpack, 3/11/78. *)
(*----------------------------------------------------------------*)
PROCEDURE daxpy( n : INTEGER;
da : FLOAT;
VAR dx : ARRAY OF FLOAT;
incx : INTEGER;
VAR dy : ARRAY OF FLOAT;
incy : INTEGER);
(*----------------------------------------------------------------*)
(* constant times a vector plus a vector (UNROLLED version). *)
(* *)
(* Jack Dongarra, linpack, 3/11/78. *)
(*----------------------------------------------------------------*)
PROCEDURE ddot( N : INTEGER;
VAR X : ARRAY OF FLOAT;
IncX : INTEGER;
VAR Y : ARRAY OF FLOAT;
IncY : INTEGER) : FLOAT;
(*----------------------------------------------------------------*)
(* Forms the dot product of two vectors. Uses unrolled loops for *)
(* increments equal to one. *)
(* Implementation : Jack Dongarra, linpack, 3/11/78. *)
(* Adopted to Modula-2, MRi, 06.09.2015 *)
(*----------------------------------------------------------------*)
PROCEDURE idamax( n : INTEGER;
VAR dx : ARRAY OF FLOAT;
incx : INTEGER) : INTEGER;
(*----------------------------------------------------------------*)
(* Finds the index of element having max. absolute value. *)
(* Jack Dongarra, linpack, 3/11/78. *)
(* Please note that indexing is starting with 0 on M2 open arrays *)
(* so if X is declared [1..n] you need to add 1 to the result. *)
(*----------------------------------------------------------------*)
PROCEDURE idamin( n : INTEGER;
VAR X : ARRAY OF FLOAT;
IncX : INTEGER) : INTEGER;
(*----------------------------------------------------------------*)
(* Finds the index of element having min. absolute value. *)
(* Jack Dongarra, linpack, 3/11/78. *)
(* Please note that indexing is starting with 0 on M2 open arrays *)
(* so if X is declared [0..n-1] you need to substract 1 from the *)
(* result. This is done to be compatible with the calls to the *)
(* Fortran equivalent routines. *)
(*----------------------------------------------------------------*)
PROCEDURE dasum( dim : CARDINAL;
VAR X : ARRAY OF FLOAT;
Inc : CARDINAL; (* Inkrementwert >= 1*)
Unrolled : BOOLEAN) : FLOAT;
(*----------------------------------------------------------------*)
(* Berechnet die Summe der Absolutwerte der im Feld X gespeich- *)
(* erten Zahlen. *)
(*----------------------------------------------------------------*)
PROCEDURE dgemv( Trans : CHAR;
M,N : INTEGER;
Alpha : FLOAT;
VAR A : ARRAY OF ARRAY OF FLOAT;
LDA : INTEGER;
VAR X : ARRAY OF FLOAT;
IncX : INTEGER;
Beta : FLOAT;
VAR Y : ARRAY OF FLOAT;
IncY : INTEGER);
(*----------------------------------------------------------------*)
(* Purpose *)
(* ======= *)
(* *)
(* DGEMV performs one of the matrix-vector operations *)
(* *)
(* y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, *)
(* *)
(* where alpha and beta are scalars, x and y are vectors and A *)
(* is an m by n matrix. *)
(* *)
(* Parameters *)
(* ========== *)
(* *)
(* Trans : On entry, Trans specifies the operation to be *)
(* performed as follows: *)
(* *)
(* 'N' or 'n' y := alpha*A *x + beta*y. *)
(* 'T' or 't' y := alpha*A'*x + beta*y. *)
(* 'C' or 'c' y := alpha*A'*x + beta*y. *)
(* *)
(* M : On entry, M specifies the number of rows of the matrix *)
(* A. M must be at least zero. *)
(* N : On entry, N specifies the number of columns of the *)
(* matrix A. *)
(* N must be at least zero. *)
(* Alpha : On entry, Alpha specifies the scalar alpha. *)
(* A : DOUBLE PRECISION array of DIMENSION ( LDA, n ). *)
(* Before entry, the leading m by n part of the array A *)
(* must contain the matrix of coefficients. *)
(* Unchanged on exit. *)
(* LDA : On entry, LDA specifies the first dimension of A as *)
(* declared in the calling (sub) program. LDA must be at *)
(* least max( 1, m ). *)
(* X : DOUBLE PRECISION array of DIMENSION at least *)
(* ( 1 + ( n - 1 )*abs( IncX ) ) when Trans = 'N' or 'n' *)
(* and at least *)
(* ( 1 + ( m - 1 )*abs( IncX ) ) otherwise. *)
(* Before entry, the incremented array X must contain the *)
(* vector x. *)
(* Unchanged on exit. *)
(* IncX : On entry, IncX specifies the increment for the *)
(* elements of X. IncX must not be zero. *)
(* Beta : On entry, Beta specifies the scalar beta. When Beta is *)
(* supplied as zero then Y need not be set on input. *)
(* Y : DOUBLE PRECISION array of DIMENSION at least *)
(* ( 1 + ( m - 1 )*abs( IncY ) ) when Trans = 'N' or 'n' *)
(* and at least *)
(* ( 1 + ( n - 1 )*abs( IncY ) ) otherwise. *)
(* Before entry with Beta non-zero, the incremented array *)
(* Y must contain the vector y. On exit, Y is overwritten *)
(* by the updated vector y. *)
(* IncY : On entry, IncY specifies the increment for the *)
(* elements of Y. IncY must not be zero. *)
(* *)
(* Level 2 Blas routine. *)
(* *)
(* -- Written on 22-October-1986. *)
(* Jack Dongarra, Argonne National Lab. *)
(* Jeremy Du Croz, Nag Central Office. *)
(* Sven Hammarling, Nag Central Office. *)
(* Richard Hanson, Sandia National Labs. *)
(* *)
(* -- Portiert nach M2 von M. Riedl, August 1995. *)
(*----------------------------------------------------------------*)
PROCEDURE dgemm( TransA,TransB : CHAR;
M,N,K : INTEGER;
Alpha : FLOAT;
VAR A : ARRAY OF ARRAY OF FLOAT;
LDA : INTEGER;
VAR B : ARRAY OF ARRAY OF FLOAT;
LDB : INTEGER;
Beta : FLOAT;
VAR C : ARRAY OF ARRAY OF FLOAT;
LDC : INTEGER);
(*----------------------------------------------------------------*)
(* Purpose *)
(* ======= *)
(* *)
(* DGEMM performs one of the matrix-matrix operations *)
(* *)
(* C := alpha*op( A )*op( B ) + beta*C, *)
(* *)
(* where op( X ) is one of *)
(* *)
(* op( X ) = X or op( X ) = X', *)
(* *)
(* Alpha and Beta are scalars, and A, B and C are matrices, with *)
(* op( A ) an m by k matrix, op( B ) a k by n matrix and C an *)
(* m by n matrix. *)
(* *)
(* Parameters *)
(* ========== *)
(* *)
(* TransA : On entry, TransA specifies the form of op( A ) to be *)
(* used in the matrix multiplication as follows: *)
(* *)
(* TransA = 'N' or 'n', op( A ) = A. *)
(* TransA = 'T' or 't', op( A ) = A'. *)
(* TransA = 'C' or 'c', op( A ) = A'. *)
(* *)
(* TransB : On entry, TransB specifies the form of op( B ) to be *)
(* used in the matrix multiplication as follows: *)
(* *)
(* TransB = 'N' or 'n', op( B ) = B. *)
(* TransB = 'T' or 't', op( B ) = B'. *)
(* TransB = 'C' or 'c', op( B ) = B'. *)
(* *)
(* M : On entry, M specifies the number of rows of the *)
(* matrix op(A) and of the matrix C. M must be at least *)
(* zero. *)
(* N : On entry, N specifies the number of columns of the *)
(* matrix op( B ) and the number of columns of the *)
(* matrix C. N must be at least zero. *)
(* K : On entry, K specifies the number of columns of the *)
(* matrix op( A ) and the number of rows of the matrix *)
(* op( B ). K must be at least zero. *)
(* Alpha : On entry, Alpha specifies the scalar alpha. *)
(* A : FLOAT array of DIMENSION ( LDA, ka ), where ka is *)
(* k when TransA = 'N' or 'n', and is m otherwise. *)
(* Before entry with TransA = 'N' or 'n', the leading *)
(* m by k part of the array A must contain the matrix *)
(* A, otherwise the leading k by m part of the array *)
(* A must contain the matrix A. *)
(* Unchanged on exit. *)
(* LDA : On entry, LDA specifies the first dimension of A as *)
(* declared in the calling (sub) program. When TransA = *)
(* 'N' or 'n' then LDA must be at least max( 1, m ), *)
(* otherwise LDA must be at least max( 1, k ). *)
(* B : FLOAT array of DIMENSION ( LDB, kb ), *)
(* where kb is n when TransB = 'N' or 'n', and is k *)
(* otherwise. Before entry with TransB = 'N' or 'n', *)
(* the leading k by n part of the array B must *)
(* contain the matrix B, otherwise the leading n by k *)
(* part of the array B must contain the matrix B. *)
(* Unchanged on exit. *)
(* LDB : On entry, LDB specifies the first dimension of B as *)
(* declared in the calling (sub) program. When TransB = *)
(* 'N' or 'n' then LDB must be at least max( 1, k ), *)
(* otherwise LDB must be at least max( 1, n ). *)
(* Beta : On entry, Beta specifies the scalar beta. When Beta *)
(* is supplied as zero then C need not be set on input. *)
(* C : FLOAT array of DIMENSION ( LDC, n ). *)
(* Before entry, the leading m by n part of the array *)
(* C must contain the matrix C, except when beta is *)
(* zero, in which case C need not be set on entry. *)
(* On exit, the array C is overwritten by the m by n *)
(* matrix ( alpha*op( A )*op( B ) + beta*C ). *)
(* LDC : On entry, LDC specifies the first dimension of C as *)
(* declared in the calling (sub) program. LDC must be *)
(* at least max( 1, m ). *)
(* *)
(* Level 3 Blas routine. *)
(* *)
(* -- Written on 8-February-1989. *)
(* Jack Dongarra, Argonne National Laboratory. *)
(* Iain Duff, AERE Harwell. *)
(* Jeremy Du Croz, Numerical Algorithms Group Ltd. *)
(* Sven Hammarling, Numerical Algorithms Group Ltd. *)
(* *)
(* -- An M2 angepasst von M. Riedl, August 1995. *)
(* *)
(* Alternative, schnellere Routinen sind in MatLib zu finden *)
(* (MatMatProd{NN|NT|TN,TT}) *)
(*----------------------------------------------------------------*)
PROCEDURE dger( m,n : CARDINAL;
Alpha : FLOAT;
VAR X : ARRAY OF FLOAT;
IncX : CARDINAL;
VAR Y : ARRAY OF FLOAT;
IncY : CARDINAL;
VAR A : ARRAY OF ARRAY OF FLOAT;
lda : CARDINAL);
(*--------------------------------------------------------------*)
(* DGER performs the rank 1 operation *)
(* *)
(* A := Alpha*(x;y) + A, *)
(* *)
(* where Alpha is a scalar, x is an m element vector, y is *)
(* an n element vector and A is an m by n matrix. *)
(* (x;y) denotes a dyadic product of the vectors x and y. *)
(* *)
(* -- Written on 22-October-1986. *)
(* Jack Dongarra, Argonne National Lab. *)
(* Jeremy Du Croz, Nag Central Office. *)
(* Sven Hammarling, Nag Central Office. *)
(* Richard Hanson, Sandia National Labs. *)
(* *)
(* Anpassung an Modula-2 : M. Riedl, Aug. 1996. *)
(* *)
(* PARAMETER *)
(* ========= *)
(* *)
(* X : Array of dimension at least (1 + ( m - 1 )* IncX) . *)
(* Before entry, the incremented array X must contain *)
(* the m element vector x. *)
(* Unchanged on exit. *)
(* IncX : On entry, IncX specifies the increment for the *)
(* elements of X. *)
(* Y : Array of dimension at least ( 1 + ( n - 1 )* IncY ). *)
(* Before entry, the incremented array Y must contain *)
(* the n element vector y. *)
(* Unchanged on exit. *)
(* IncY : On entry, IncY specifies the increment for the *)
(* elements of Y. *)
(* A : Array of DIMENSION ( m, n ). *)
(* Before entry, the leading m by n part of the array A *)
(* must contain the matrix of coefficients. On exit, A *)
(* is overwritten by the updated matrix. *)
(* *)
(*--------------------------------------------------------------*)
PROCEDURE zswap( N : CARDINAL;
VAR X : ARRAY OF CFLOAT;
IncX : INTEGER;
VAR Y : ARRAY OF CFLOAT;
IncY : INTEGER);
(*----------------------------------------------------------------*)
(* Swap complex vectors X and Y *)
(*----------------------------------------------------------------*)
PROCEDURE zcopy( N : INTEGER;
VAR X : ARRAY OF CFLOAT;
IncX : INTEGER;
VAR Y : ARRAY OF CFLOAT;
IncY : INTEGER);
(*----------------------------------------------------------------*)
(* copies a vector, x, to a vector, y. *)
(* uses unrolled loops for increments equal to one. *)
(* jack dongarra, linpack, 3/11/78. *)
(* MRi, Modula-2 10.04.16 | 09.09.18 (complex version) *)
(*----------------------------------------------------------------*)
PROCEDURE zdotc( N : INTEGER;
VAR X : ARRAY OF CFLOAT;
IncX : INTEGER;
VAR Y : ARRAY OF CFLOAT;
IncY : INTEGER) : CFLOAT;
(*----------------------------------------------------------------*)
(* Forms the dot product of two vectors. Uses unrolled loops for *)
(* increments equal to one. *)
(*----------------------------------------------------------------*)
PROCEDURE dznrm2( N : INTEGER;
VAR X : ARRAY OF CFLOAT;
IncX : INTEGER) : FLOAT;
(*----------------------------------------------------------------*)
(* dznrm2 returns the euclidean norm of a vector so that *)
(* dznrm2 := sqrt( X**H*X ) *)
(*----------------------------------------------------------------*)
PROCEDURE zscal( n : INTEGER;
da : CFLOAT;
VAR dx : ARRAY OF CFLOAT;
IncX : INTEGER);
(*----------------------------------------------------------------*)
(* Scales a vector by a constant (UNROLLED version) *)
(*----------------------------------------------------------------*)
PROCEDURE zaxpy( n : INTEGER;
da : CFLOAT;
VAR X : ARRAY OF CFLOAT;
IncX : INTEGER;
VAR Y : ARRAY OF CFLOAT;
IncY : INTEGER);
(*----------------------------------------------------------------*)
(* constant times a vector plus a vector *)
(*----------------------------------------------------------------*)
PROCEDURE zdrot( N : INTEGER;
VAR X : ARRAY OF CFLOAT;
IncX : INTEGER;
VAR Y : ARRAY OF CFLOAT;
IncY : INTEGER;
c,s : FLOAT);
(*----------------------------------------------------------------*)
(* Applies a plane rotation, where the cos and sin (c and s) are *)
(* real and the vectors cx and cy are complex. *)
(*----------------------------------------------------------------*)
PROCEDURE zgemv( Trans : CHAR;
M,N : INTEGER;
Alpha : CFLOAT;
VAR A : ARRAY OF ARRAY OF CFLOAT;
lda : INTEGER;
VAR X : ARRAY OF CFLOAT;
IncX : INTEGER;
Beta : CFLOAT;
VAR Y : ARRAY OF CFLOAT;
IncY : INTEGER);
(*----------------------------------------------------------------*)
(* Performs one of the matrix-vector operations *)
(* *)
(* Y = Alpha* A *x + Beta*Y, or *)
(* Y = Alpha* A' *x + Beta*Y, or *)
(* Y = Alpha*conjg(A')*x + Beta*Y,; *)
(* *)
(* where Alpha and Beta are scalars, X and Y are vectors *)
(* and A is an M by N matrix. *)
(* *)
(* parameters *)
(* *)
(* Trans : trans specifies the operation to be performed as *)
(* follows: *)
(* *)
(* trans = 'N' or 'n' Y = Alpha*A*X + Beta*Y *)
(* trans = 'T' or 't' Y = Alpha*A'*X + Beta*Y *)
(* trans = 'C' or 'c' Y = Alpha*conjg(A')*X + Beta*Y *)
(* *)
(* *)
(* M : M specifies the number of rows of the matrix A. *)
(* M must be at least zero. *)
(* N : N specifies the number of columns of the matrix A. *)
(* N must be at least zero. *)
(* Alpha : Alpha specifies the scalar Alpha *)
(* A : array of dimension (lda,N) *)
(* Before entry, the leading M by N part of the *)
(* array A must contain the matrix of coefficients. *)
(* Unchanged on exit. *)
(* lda : on entry, lda specifies the first dimension *)
(* of A as declared in the calling (sub) program. *)
(* lda must be at least max(1,M). *)
(* X : array of dimension at least (1+(n-1 )*abs(incx)) *)
(* when trans = 'n' or 'N' and at least *)
(* (1+(m-1 )*abs(incx)) otherwise. *)
(* Before entry, the incremented array x must contain *)
(* the vector X. Unchanged on exit. *)
(* IncX : IncX specifies the increment for the elements of *)
(* X. IncX must not be zero. *)
(* Beta : Beta specifies the scalar Beta. When Beta is *)
(* supplied as zero then y need not be set on input. *)
(* Unchanged on exit. *)
(* Y : array of dimension at least (1+(m-1)*abs(incy)) *)
(* when trans = 'N' or 'n' and at least *)
(* (1+(n-1 *abs(incy)) otherwise. *)
(* Before entry with beta non-zero, the incremented *)
(* array Y must contain the vector Y. On exit, Y is *)
(* overwritten by the updated vector Y. *)
(* IncY : IncY specifies the increment for the elements of *)
(* Y. IncY must not be zero. *)
(* *)
(* level 2 blas routine. *)
(* *)
(* -- written on 22-october-1986. *)
(* Jack Dongarra, Argonne National Lab. *)
(* Jeremy du Croz, NAG central office. *)
(* Sven Hammarling, NAG central office. *)
(* Richard Hanson, Sandia National Labs. *)
(*----------------------------------------------------------------*)
PROCEDURE zgemm( TransA,TransB : CHAR;
M,N,K : INTEGER;
Alpha : CFLOAT;
VAR A : ARRAY OF ARRAY OF CFLOAT;
LDA : INTEGER;
VAR B : ARRAY OF ARRAY OF CFLOAT;
LDB : INTEGER;
Beta : CFLOAT;
VAR C : ARRAY OF ARRAY OF CFLOAT;
LDC : INTEGER);
(*----------------------------------------------------------------*)
(* Purpose *)
(* ======= *)
(* *)
(* ZGEMM performs one of the matrix-matrix operations *)
(* *)
(* C := alpha*op( A )*op( B ) + beta*C, *)
(* *)
(* where op( X ) is one of *)
(* *)
(* op( X ) = X or op( X ) = X', *)
(* *)
(* Alpha and Beta are scalars, and A, B and C are matrices, with *)
(* op( A ) an m by k matrix, op( B ) a k by n matrix and C an *)
(* m by n matrix. *)
(* *)
(* Parameters *)
(* ========== *)
(* *)
(* TransA : On entry, TransA specifies the form of op( A ) to be *)
(* used in the matrix multiplication as follows: *)
(* *)
(* TransA = 'N' or 'n', op( A ) = A. *)
(* TransA = 'T' or 't', op( A ) = A'. *)
(* TransA = 'C' or 'c', op( A ) = conj(A'). *)
(* *)
(* TransB : On entry, TransB specifies the form of op( B ) to be *)
(* used in the matrix multiplication as follows: *)
(* *)
(* TransB = 'N' or 'n', op( B ) = B. *)
(* TransB = 'T' or 't', op( B ) = B'. *)
(* TransB = 'C' or 'c', op( B ) = conj(B)'. *)
(* *)
(* M : On entry, M specifies the number of rows of the *)
(* matrix op(A) and of the matrix C. M must be at least *)
(* zero. *)
(* N : On entry, N specifies the number of columns of the *)
(* matrix op( B ) and the number of columns of the *)
(* matrix C. N must be at least zero. *)
(* K : On entry, K specifies the number of columns of the *)
(* matrix op( A ) and the number of rows of the matrix *)
(* op( B ). K must be at least zero. *)
(* Alpha : On entry, Alpha specifies the scalar alpha. *)
(* A : COMPLEX array of DIMENSION ( LDA, ka ), where ka is *)
(* k when TransA = 'N' or 'n', and is m otherwise. *)
(* Before entry with TransA = 'N' or 'n', the leading *)
(* m by k part of the array A must contain the matrix *)
(* A, otherwise the leading k by m part of the array *)
(* A must contain the matrix A. *)
(* Unchanged on exit. *)
(* LDA : On entry, LDA specifies the first dimension of A as *)
(* declared in the calling (sub) program. When TransA = *)
(* 'N' or 'n' then LDA must be at least max( 1, m ), *)
(* otherwise LDA must be at least max( 1, k ). *)
(* B : COMPLEX array of DIMENSION ( LDB, kb ), *)
(* where kb is n when TransB = 'N' or 'n', and is k *)
(* otherwise. Before entry with TransB = 'N' or 'n', *)
(* the leading k by n part of the array B must *)
(* contain the matrix B, otherwise the leading n by k *)
(* part of the array B must contain the matrix B. *)
(* Unchanged on exit. *)
(* LDB : On entry, LDB specifies the first dimension of B as *)
(* declared in the calling (sub) program. When TransB = *)
(* 'N' or 'n' then LDB must be at least max( 1, k ), *)
(* otherwise LDB must be at least max( 1, n ). *)
(* Beta : On entry, Beta specifies the scalar beta. When Beta *)
(* is supplied as zero then C need not be set on input. *)
(* C : COMPLEX array of DIMENSION ( LDC, n ). *)
(* Before entry, the leading m by n part of the array *)
(* C must contain the matrix C, except when beta is *)
(* zero, in which case C need not be set on entry. *)
(* On exit, the array C is overwritten by the m by n *)
(* matrix ( alpha*op( A )*op( B ) + beta*C ). *)
(* LDC : On entry, LDC specifies the first dimension of C as *)
(* declared in the calling (sub) program. LDC must be *)
(* at least max( 1, m ). *)
(* *)
(* Level 3 Blas routine. *)
(* *)
(* -- Written on 8-February-1989. *)
(* Jack Dongarra, Argonne National Laboratory. *)
(* Iain Duff, AERE Harwell. *)
(* Jeremy Du Croz, Numerical Algorithms Group Ltd. *)
(* Sven Hammarling, Numerical Algorithms Group Ltd. *)
(*----------------------------------------------------------------*)
END LibDBlasM2.