DEFINITION MODULE LibDBlas;

  (*------------------------------------------------------------------------*)
  (* Basic Linear Algebra Prozeduren in Modula-2 wobei die Parameter        *)
  (* konform zur den Fortran-Routinen uebergeben werden.                    *)
  (* Some BLAS routines implemented in Modula-2 pathing the arguments in a  *)
  (* "Fortran manner".                                                      *)
  (*------------------------------------------------------------------------*)
  (* 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.                   *)
  (* The routines can be used as an alternative to those provided in Module *)
  (* "LibDBlasM2". They are doing some (ugly) address translation and can   *)
  (* be used when translating Fortran code to Modula-2.                     *)
  (* The routine in LibDBlasM2 use "normal" M2 mechanism for handing over   *)
  (* open array parameters and should be prefered for new code.             *)
  (* 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 vectors to 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                            *)
  (*                                                                        *)
  (*   complex procedures                                                   *)
  (*                                                                        *)
  (*   zswap     : swap two vectors                                         *)
  (*   zcopy     : copy vectors to vector                                   *)
  (*   zdotc     : dot product of two vectors                               *)
  (*   dznrm2    : euklidian norm of a vector                               *)
  (*   zscal     : scale vector by a constant                               *)
  (*   zaxpy     : constant times a vector plus a vector                    *)
  (*   zdrot     : applies a plane rotation                                 *)
  (*                                                                        *)
  (*   zgemv     : matrix matrix operations (complex)                       *)
  (*   zgemm     : matrix matrix operations (complex)                       *)
  (*                                                                        *)
  (*------------------------------------------------------------------------*)
  (* Hinweis:                                                               *)
  (* ========                                                               *)
  (* Bei den BLAS-Routinen werden Vektoren (ARRAY OF FLOAT) nur als "FLOAT" *)
  (* uebergeben. Dies entspricht der Benutzung innerhalb von Fortran-       *)
  (* Routinen. Die Adressen werden innehalb der Routinen dann als offene    *)
  (* Felder genutzt. Eine Ueberpruefung von Indizes kann so leider nicht    *)
  (* stattfinden.                                                           *)
  (*------------------------------------------------------------------------*)
  (* Implementation : Michael Riedl                                         *)
  (* Licence        : GNU Lesser General Public License (LGPL)              *)
  (*------------------------------------------------------------------------*)

  (* $Id: LibDBlas.def,v 1.7 2018/09/12 13:20:49 mriedl Exp mriedl $ *)

FROM Deklera IMPORT FLOAT,CFLOAT; (* REAL & COMPLEX type *)
             IMPORT LibDBlasM2; (* dgemv,dgemm,zgemm,dger *)

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 X    : (* 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 X    : (* ARRAY OF *) FLOAT; 
                    IncX : INTEGER;
                VAR Y    : (* 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 X    : (* 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 [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 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(    n        : CARDINAL;
                VAR X        : (* ARRAY OF *) FLOAT;
                    IncX     : CARDINAL) : 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);
 *)

CONST dgemv = LibDBlasM2.dgemv;

          (*----------------------------------------------------------------*)
          (* 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     : Real 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     : Real 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     : Real 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);
 *)

CONST dgemm = LibDBlasM2.dgemm;

          (*----------------------------------------------------------------*)
          (*  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      : Real 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      : Real 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      : Real 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}) wenn nicht eine "schnelle" BLAS      *)
          (* Bibiliotek (z.B. OpenBlas, ATLAS) eingebunden wird             *)
          (*                                                                *)
          (*  M : Anzahl der Zeilen der Matrix op(A) und der Matrix C       *)
          (*  N : Anzahl der Spalten der Matrix op(B) und der Matrix C      *)
          (*  K : Anzahl der Zeilen der Matrix op(A) und der Spalten        *)
          (*      der Matrix op(B)                                          *)
          (*----------------------------------------------------------------*)

(*
 * 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);
 *)

CONST dger = LibDBlasM2.dger;

           (*--------------------------------------------------------------*)
           (*  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.                *)
           (*                                                              *)
           (*--------------------------------------------------------------*)

(*============================= complexe procedures ========================*)

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);
 *)

CONST zgemv = LibDBlasM2.zgemv;

          (*----------------------------------------------------------------*)
          (* 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);
 *)

CONST zgemm = LibDBlasM2.zgemm;

          (*----------------------------------------------------------------*)
          (*  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 LibDBlas.
