define(s_type,`real(skind)')
define(d_type,`real(dkind)')
define(c_type,`complex(skind)')
define(z_type,`complex(dkind)')

define(II,`intent(in)')
define(IIO,`intent(inout)')
define(DM1,`dimension(:)')
define(DM1a,`dimension(*)')
define(DM2,`dimension(:,:)')
define(OPT,`optional')

define(MAKEPUBLIC,`
ifdef(``ISPUBLIC'$1`'',`! $1 already public',`
public $1
define(``ISPUBLIC'$1`'')')')

define(BIND,`interface $1
  module procedure s$2_95
  module procedure d$2_95
  module procedure c$2_95
  module procedure z$2_95
end interface
MAKEPUBLIC($1)')
define(BINDR,`interface $1
  module procedure s$2_95
  module procedure d$2_95
end interface
MAKEPUBLIC($1)')
define(BINDC,`interface $1
  module procedure c$2_95
  module procedure z$2_95
end interface
MAKEPUBLIC($1)')

define(BTOPARG,`type(blas_$1_type),intent(in),optional:: $2
character:: v$2')

define(XOPARG,`$1_type,intent(in),optional:: $2
$1_type:: v$2')

define(CHARDEF,`
v$1 = "$2"
if (present($1)) v$1 = $1%c')
define(ONEDEF,`
v$1 = $2_one
if (present($1)) v$1 = $1')
define(ZERODEF,`
v$1 = $2_zero
if (present($1)) v$1 = $1')

define(LDCALC,`ld$1 = 1
if (size($1,2) > 1 .and. size($1,1) > 0) ld$1 = (loc($1(1,2))-loc($1(1,1)))/$2size')
define(INCCALC,`inc$1 = 1
if (size($1) > 1) inc$1 = (loc($1(2))-loc($1(1)))/$2size')


module blas_common
implicit none
integer,parameter:: skind = kind(1.e0), dkind = kind(1.d0)

! ones
s_type,parameter:: s_one = real(1,kind=skind)
d_type,parameter:: d_one = real(1,kind=dkind)
c_type,parameter:: c_one = cmplx(1,kind=skind)
z_type,parameter:: z_one = cmplx(1,kind=dkind)
! zeros
s_type,parameter:: s_zero = real(0,kind=skind)
d_type,parameter:: d_zero = real(0,kind=dkind)
c_type,parameter:: c_zero = cmplx(0,kind=skind)
z_type,parameter:: z_zero = cmplx(0,kind=dkind)
! enable these if your compiler supports sizeof
!integer,parameter:: ssize = sizeof(s_one)
!integer,parameter:: dsize = sizeof(d_one)
!integer,parameter:: csize = sizeof(c_one)
!integer,parameter:: zsize = sizeof(z_one)
! sizeof is less common extension, so provide sizes manually
integer,parameter:: ssize = 4 
integer,parameter:: dsize = 8 
integer,parameter:: csize = 2*ssize
integer,parameter:: zsize = 2*dsize

type blas_trans_type
  character:: c
end type
type(blas_trans_type),parameter:: blas_trans = blas_trans_type("T")
type(blas_trans_type),parameter:: blas_no_trans = blas_trans_type("N")
type(blas_trans_type),parameter:: blas_conj_trans = blas_trans_type("H")

type blas_uplo_type
  character:: c
end type
type(blas_uplo_type),parameter:: blas_lower = blas_uplo_type("L")
type(blas_uplo_type),parameter:: blas_upper = blas_uplo_type("U")

type blas_diag_type
  character:: c
end type
type(blas_diag_type),parameter:: blas_unit_diag = blas_diag_type("U")
type(blas_diag_type),parameter:: blas_non_unit_diag = blas_diag_type("N")

type blas_side_type
  character:: c
end type
type(blas_side_type),parameter:: blas_left_side = blas_side_type("L")
type(blas_side_type),parameter:: blas_right_side = blas_side_type("R")

end module

module blas95
use blas_common
!================================================================================
! This module provides Fortran 95 wrappers to the legacy BLAS. It allows a
! Fortran programmer to use BLAS conveniently using array sections, whole
! arrays and optional arguments. The Level 2 and 3 interfaces here are BLAST
! compliant.
!================================================================================

implicit none
private

public blas_trans, blas_no_trans, blas_conj_trans
public blas_lower, blas_upper
public blas_unit_diag, blas_non_unit_diag
public blas_left_side, blas_right_side 

BIND(GEMM,GEMV)
BIND(GBMV,GBMV)
BINDR(SYMM,SYMV)
BINDR(SBMV,SBMV)
BINDR(SPMV,SPMV)
BINDC(HEMM,HEMV)
BINDC(HBMV,HBMV)
BINDC(HPMV,HPMV)
BIND(TRMM,TRMV)
BIND(TBMV,TBMV)
BIND(TPMV,TPMV)
BIND(TRSM,TRSV)
BIND(TBSV,TBSV)
BIND(TPSV,TPSV)
BINDR(GEMM,GER)
BINDC(GEMM,GERUC)
BINDR(SYRK,SYR)
BINDR(SPR,SPR)
BINDC(HERK,HER)
BINDC(HPR,HPR)
BINDR(SYR2K,SYR2)
BINDR(SPR2,SPR)
BINDC(HER2K,HER2)
BINDC(HPR2,HPR)
BIND(GEMM,GEMM)
BIND(SYMM,SYMM)
BINDC(HEMM,HEMM)
BIND(TRMM,TRMM)
BIND(TRSM,TRSM)
BINDR(SYRK,SYRK)
BINDC(HERK,HERK)
BINDR(SYR2K,SYR2K)
BINDC(HER2K,HER2K)
  
public blas_OK_vector
interface blas_OK_vector
  module procedure sblas_OK_vector
  module procedure dblas_OK_vector
  module procedure cblas_OK_vector
  module procedure zblas_OK_vector
end interface
public blas_OK_matrix
interface blas_OK_matrix
  module procedure sblas_cont_vector
  module procedure dblas_cont_vector
  module procedure cblas_cont_vector
  module procedure zblas_cont_vector
  module procedure sblas_OK_matrix
  module procedure dblas_OK_matrix
  module procedure cblas_OK_matrix
  module procedure zblas_OK_matrix
end interface

contains

define(blas_OK_vec,`
function $1blas_OK_vector(vec) result(OK)
$1_type,dimension(:):: vec
logical:: OK
OK = .true.
if (size(vec) > 1) then
  OK = mod(loc(vec(2))-loc(vec(1)),$1size) == 0
end if
end function')
blas_OK_vec(s)
blas_OK_vec(d)
blas_OK_vec(c)
blas_OK_vec(z)

define(blas_OK_mat,`
function $1blas_OK_matrix(mat) result(OK)
$1_type,dimension(:,:):: mat
logical:: OK
integer:: inc2
OK = .true.
if (size(mat,1) > 1 .and. size(mat,2) > 0) then
  OK = loc(mat(2,1))-loc(mat(1,1)) == $1size
end if
if (size(mat,2) > 1 .and. size(mat,1) > 0) then
  inc2 = loc(mat(1,2))-loc(mat(1,1))
  OK = OK .and. mod(inc2,$1size) == 0 .and. inc2 > 0
end if
end function')
blas_OK_mat(s)
blas_OK_mat(d)
blas_OK_mat(c)
blas_OK_mat(z)

define(blas_cont_vec,`
function $1blas_cont_vector(mat) result(OK)
$1_type,dimension(:):: mat
logical:: OK
OK = .true.
if (size(mat) > 1) OK = loc(mat(2))-loc(mat(1)) == $1size
end function')
blas_cont_vec(s)
blas_cont_vec(d)
blas_cont_vec(c)
blas_cont_vec(z)

!================================================================================
! LEVEL 2 BLAS routines
!================================================================================

!========================================
! GEMV implementation
define(GEMV_95,`
subroutine $1GEMV_95(a,b,c,transa,alpha,beta)
$1_type,DM2,II:: a
$1_type,DM1,II:: b
$1_type,DM1,IIO:: c
BTOPARG(trans,transa)
XOPARG($1,alpha)
XOPARG($1,beta)
integer:: lda,incb,incc
integer:: m,n
CHARDEF(transa,N)
ONEDEF(alpha,$1)
ZERODEF(beta,$1)
LDCALC(a,$1)
INCCALC(b,$1)
INCCALC(c,$1)
m = size(a,1)
n = size(a,2)
call $1GEMV(vtransa,m,n,valpha,a(1,1),lda,b(1),incb,vbeta,c(1),incc)
end subroutine')

GEMV_95(s)
GEMV_95(d)
GEMV_95(c)
GEMV_95(z)

!========================================
! GBMV implementation
define(GBMV_95,`
subroutine $1GBMV_95(a,m,kl,x,y,trans,alpha,beta)
$1_type,DM2,II:: a
$1_type,DM1,II:: x
$1_type,DM1,IIO:: y
integer,II:: m,kl
BTOPARG(trans,trans)
XOPARG($1,alpha)
XOPARG($1,beta)
integer:: n,ku,lda,incx,incy
CHARDEF(trans,N)
ONEDEF(alpha,$1)
ZERODEF(beta,$1)
n = size(x)
ku = size(a,1)-1-kl
LDCALC(a,$1)
INCCALC(x,$1)
INCCALC(y,$1)
call $1GBMV(vtrans,m,n,kl,ku,valpha,a(1,1),lda,x(1),incx,vbeta,y(1),incy)
end subroutine')
GBMV_95(s)
GBMV_95(d)
GBMV_95(c)
GBMV_95(z)

!========================================
! SYMV implementation
define(SYMV_95,`
subroutine $1SYMV_95(a,b,c,uplo,alpha,beta)
$1_type,DM2,II:: a
$1_type,DM1,II:: b
$1_type,DM1,IIO:: c
BTOPARG(uplo,uplo)
XOPARG($1,alpha)
XOPARG($1,beta)
integer:: lda,incb,incc
integer:: n
CHARDEF(uplo,U)
ONEDEF(alpha,$1)
ZERODEF(beta,$1)
LDCALC(a,$1)
INCCALC(b,$1)
INCCALC(c,$1)
n = size(a,2)
call $1SYMV(vuplo,n,valpha,a(1,1),lda,b(1),incb,vbeta,c(1),incc)
end subroutine')

SYMV_95(s)
SYMV_95(d)

!========================================
! SBMV implementation
define(SBMV_95,`
subroutine $1SBMV_95(a,x,y,uplo,alpha,beta)
$1_type,DM2,II:: a
$1_type,DM1,II:: x
$1_type,DM1,IIO:: y
BTOPARG(uplo,uplo)
XOPARG($1,alpha)
XOPARG($1,beta)
integer:: n,k,lda,incx,incy
CHARDEF(uplo,U)
ONEDEF(alpha,$1)
ZERODEF(beta,$1)
n = size(x)
k = (size(a,1)-1)/2
LDCALC(a,$1)
INCCALC(x,$1)
INCCALC(y,$1)
call $1SBMV(vuplo,n,k,valpha,a(1,1),lda,x(1),incx,vbeta,y(1),incy)
end subroutine')
SBMV_95(s)
SBMV_95(d)

!========================================
! SPMV implementation
define(SPMV_95,`
subroutine $1SPMV_95(ap,x,y,uplo,alpha,beta)
$1_type,DM1a,II:: ap
$1_type,DM1,II:: x
$1_type,DM1,IIO:: y
BTOPARG(uplo,uplo)
XOPARG($1,alpha)
XOPARG($1,beta)
integer:: n,incx,incy
CHARDEF(uplo,U)
ONEDEF(alpha,$1)
ZERODEF(beta,$1)
n = size(x)
INCCALC(x,$1)
INCCALC(y,$1)
call $1SPMV(vuplo,n,valpha,ap(1),x(1),incx,vbeta,y(1),incy)
end subroutine')
SPMV_95(s)
SPMV_95(d)

!========================================
! HEMV implementation
define(HEMV_95,`
subroutine $1HEMV_95(a,b,c,uplo,alpha,beta)
$1_type,DM2,II:: a
$1_type,DM1,II:: b
$1_type,DM1,IIO:: c
BTOPARG(uplo,uplo)
XOPARG($1,alpha)
XOPARG($1,beta)
integer:: lda,incb,incc
integer:: n
CHARDEF(uplo,U)
ONEDEF(alpha,$1)
ZERODEF(beta,$1)

LDCALC(a,$1)
INCCALC(b,$1)
INCCALC(c,$1)
n = size(a,2)
call $1HEMV(vuplo,n,valpha,a(1,1),lda,b(1),incb,vbeta,c(1),incc)
end subroutine')

HEMV_95(c)
HEMV_95(z)

!========================================
! HBMV implementation
define(HBMV_95,`
subroutine $1HBMV_95(a,x,y,uplo,alpha,beta)
$1_type,DM2,II:: a
$1_type,DM1,II:: x
$1_type,DM1,IIO:: y
BTOPARG(uplo,uplo)
XOPARG($1,alpha)
XOPARG($1,beta)
integer:: n,k,lda,incx,incy
CHARDEF(uplo,U)
ONEDEF(alpha,$1)
ZERODEF(beta,$1)
n = size(x)
k = (size(a,1)-1)/2
LDCALC(a,$1)
INCCALC(x,$1)
INCCALC(y,$1)
call $1HBMV(vuplo,n,k,valpha,a(1,1),lda,x(1),incx,vbeta,y(1),incy)
end subroutine')
HBMV_95(c)
HBMV_95(z)

!========================================
! HPMV implementation
define(HPMV_95,`
subroutine $1HPMV_95(ap,x,y,uplo,alpha,beta)
$1_type,DM1a,II:: ap
$1_type,DM1,II:: x
$1_type,DM1,IIO:: y
BTOPARG(uplo,uplo)
XOPARG($1,alpha)
XOPARG($1,beta)
integer:: n,incx,incy
CHARDEF(uplo,U)
ONEDEF(alpha,$1)
ZERODEF(beta,$1)
n = size(x)
INCCALC(x,$1)
INCCALC(y,$1)
call $1HPMV(vuplo,n,valpha,ap(1),x(1),incx,vbeta,y(1),incy)
end subroutine')
HPMV_95(c)
HPMV_95(z)

!========================================
! TRMV implementation
define(TRMV_95,`
subroutine $1TRMV_95(t,b,side,uplo,transt,diag)
$1_type,DM2,II:: t
$1_type,DM1,IIO:: b
BTOPARG(side,side)
BTOPARG(uplo,uplo)
BTOPARG(trans,transt)
BTOPARG(diag,diag)
integer:: ldt,incb
integer:: n
CHARDEF(side,L)
CHARDEF(transt,N)
CHARDEF(uplo,U)
CHARDEF(diag,N)
LDCALC(t,$1)
INCCALC(b,$1)
n = size(t,1)
call $1TRMV(vuplo,vtranst,vdiag,n,t(1,1),ldt,b(1),incb)
end subroutine')

TRMV_95(s)
TRMV_95(d)
TRMV_95(c)
TRMV_95(z)

!========================================
! TBMV implementation
define(TBMV_95,`
subroutine $1TBMV_95(a,x,y,uplo,trans,diag)
$1_type,DM2,II:: a
$1_type,DM1,II:: x
$1_type,DM1,IIO:: y
BTOPARG(uplo,uplo)
BTOPARG(trans,trans)
BTOPARG(diag,diag)
integer:: n,k,lda,incx
CHARDEF(uplo,U)
CHARDEF(trans,N)
CHARDEF(diag,N)
n = size(x)
k = size(a,1)-1
LDCALC(a,$1)
INCCALC(x,$1)
call $1TBMV(vuplo,vtrans,vdiag,n,k,a(1,1),lda,x(1),incx)
end subroutine')
TBMV_95(s)
TBMV_95(d)
TBMV_95(c)
TBMV_95(z)

!========================================
! TPMV implementation
define(TPMV_95,`
subroutine $1TPMV_95(ap,x,y,uplo,trans,diag)
$1_type,DM1a,II:: ap
$1_type,DM1,II:: x
$1_type,DM1,IIO:: y
BTOPARG(uplo,uplo)
BTOPARG(trans,trans)
BTOPARG(diag,diag)
integer:: n,k,incx,incy
CHARDEF(uplo,U)
CHARDEF(trans,N)
CHARDEF(diag,N)
n = size(x)
INCCALC(x,$1)
INCCALC(y,$1)
call $1TPMV(vuplo,vtrans,vdiag,n,ap(1),x(1),incx,y(1),incy)
end subroutine')
TPMV_95(s)
TPMV_95(d)
TPMV_95(c)
TPMV_95(z)

!========================================
! TRSV implementation
define(TRSV_95,`
subroutine $1TRSV_95(t,b,side,uplo,transt,diag)
$1_type,DM2,II:: t
$1_type,DM1,IIO:: b
BTOPARG(side,side)
BTOPARG(uplo,uplo)
BTOPARG(trans,transt)
BTOPARG(diag,diag)
integer:: ldt,incb
integer:: n
CHARDEF(side,L)
CHARDEF(transt,N)
CHARDEF(uplo,U)
CHARDEF(diag,N)
LDCALC(t,$1)
INCCALC(b,$1)
n = size(t,1)
call $1TRSV(vuplo,vtranst,vdiag,n,t(1,1),ldt,b(1),incb)
end subroutine')

TRSV_95(s)
TRSV_95(d)
TRSV_95(c)
TRSV_95(z)

!========================================
! TBSV implementation
define(TBSV_95,`
subroutine $1TBSV_95(a,x,y,uplo,trans,diag)
$1_type,DM2,II:: a
$1_type,DM1,II:: x
$1_type,DM1,IIO:: y
BTOPARG(uplo,uplo)
BTOPARG(trans,trans)
BTOPARG(diag,diag)
integer:: n,k,lda,incx
CHARDEF(uplo,U)
CHARDEF(trans,N)
CHARDEF(diag,N)
n = size(x)
k = size(a,1)-1
LDCALC(a,$1)
INCCALC(x,$1)
call $1TBSV(vuplo,vtrans,vdiag,n,k,a(1,1),lda,x(1),incx)
end subroutine')
TBSV_95(s)
TBSV_95(d)
TBSV_95(c)
TBSV_95(z)

!========================================
! TPSV implementation
define(TPSV_95,`
subroutine $1TPSV_95(ap,x,y,uplo,trans,diag)
$1_type,DM1a,II:: ap
$1_type,DM1,II:: x
$1_type,DM1,IIO:: y
BTOPARG(uplo,uplo)
BTOPARG(trans,trans)
BTOPARG(diag,diag)
integer:: n,k,incx,incy
CHARDEF(uplo,U)
CHARDEF(trans,N)
CHARDEF(diag,N)
n = size(x)
INCCALC(x,$1)
INCCALC(y,$1)
call $1TPSV(vuplo,vtrans,vdiag,n,ap(1),x(1),incx,y(1),incy)
end subroutine')
TPSV_95(s)
TPSV_95(d)
TPSV_95(c)
TPSV_95(z)

!========================================
! GER implementation
define(GER_95,`
subroutine $1GER_95(a,b,c,alpha)
$1_type,DM1,II:: a
$1_type,DM1,II:: b
$1_type,DM2,IIO:: c
XOPARG($1,alpha)
integer:: ldc,inca,incb
integer:: m,n
ONEDEF(alpha,$1)
LDCALC(c,$1)
INCCALC(a,$1)
INCCALC(b,$1)
m = size(c,1)
n = size(c,2)
call $1GER(m,n,valpha,a(1),inca,b(1),incb,c(1,1),ldc)
end subroutine')

GER_95(s)
GER_95(d)

define(GERUC_95,`
subroutine $1GERUC_95(a,b,c,transb,alpha)
$1_type,DM1,II:: a
$1_type,DM1,II:: b
$1_type,DM2,IIO:: c
BTOPARG(trans,transb)
XOPARG($1,alpha)
integer:: ldc,inca,incb
integer:: m,n
CHARDEF(transb,T)
ONEDEF(alpha,$1)
LDCALC(c,$1)
INCCALC(a,$1)
INCCALC(b,$1)
m = size(c,1)
n = size(c,2)
if (vtransb == "H") then
  call $1GERC(m,n,valpha,a(1),inca,b(1),incb,c(1,1),ldc)
else 
  call $1GERU(m,n,valpha,a(1),inca,b(1),incb,c(1,1),ldc)
end if
end subroutine')

GERUC_95(c)
GERUC_95(z)

!========================================
! SYR implementation
define(SYR_95,`
subroutine $1SYR_95(a,c,uplo,alpha)
$1_type,DM1,II:: a
$1_type,DM2,IIO:: c
BTOPARG(uplo,uplo)
XOPARG($1,alpha)
integer:: n,ldc,inca
CHARDEF(uplo,U)
ONEDEF(alpha,$1)
LDCALC(c,$1)
INCCALC(a,$1)
n = size(a)
call $1SYR(vuplo,n,valpha,a(1),inca,c(1,1),ldc)
end subroutine')
SYR_95(s)
SYR_95(d)

!========================================
! SPR implementation
define(SPR_95,`
subroutine $1SPR_95(x,ap,uplo,alpha)
$1_type,DM1,II:: x
$1_type,DM1,IIO:: ap
BTOPARG(uplo,uplo)
XOPARG($1,alpha)
integer:: n,incx
CHARDEF(uplo,U)
ONEDEF(alpha,$1)
INCCALC(x,$1)
n = size(x)
call $1SPR(vuplo,n,valpha,x(1),incx,ap(1))
end subroutine')
SPR_95(s)
SPR_95(d)

!========================================
! HER implementation
define(HER_95,`
subroutine $1HER_95(a,c,uplo,alpha)
$1_type,DM1,II:: a
$1_type,DM2,IIO:: c
BTOPARG(uplo,uplo)
XOPARG($1,alpha)
integer:: n,ldc,inca
CHARDEF(uplo,U)
ONEDEF(alpha,$1)
LDCALC(c,$1)
INCCALC(a,$1)
n = size(a)
call $1HER(vuplo,n,valpha,a(1),inca,c(1,1),ldc)
end subroutine')
HER_95(c)
HER_95(z)

!========================================
! HPR implementation
define(HPR_95,`
subroutine $1HPR_95(x,ap,uplo,alpha)
$1_type,DM1,II:: x
$1_type,DM1,IIO:: ap
BTOPARG(uplo,uplo)
XOPARG($1,alpha)
integer:: n,incx
CHARDEF(uplo,U)
ONEDEF(alpha,$1)
INCCALC(x,$1)
n = size(x)
call $1HPR(vuplo,n,valpha,x(1),incx,ap(1))
end subroutine')
HPR_95(c)
HPR_95(z)

!========================================
! SYR2 implementation
define(SYR2_95,`
subroutine $1SYR2_95(a,b,c,uplo,alpha)
$1_type,DM1,II:: a,b
$1_type,DM2,IIO:: c
BTOPARG(uplo,uplo)
XOPARG($1,alpha)
integer:: n,inca,incb,ldc
CHARDEF(uplo,U)
ONEDEF(alpha,$1)
INCCALC(a,$1)
INCCALC(b,$1)
LDCALC(c,$1)
n = size(c,1)
call $1SYR2(vuplo,n,valpha,a(1),inca,b(1),incb,c(1,1),ldc)
end subroutine')
SYR2_95(s)
SYR2_95(d)

!========================================
! SPR2 implementation
define(SPR2_95,`
subroutine $1SPR2_95(x,y,ap,uplo,alpha)
$1_type,DM1,II:: x,y
$1_type,DM1,IIO:: ap
BTOPARG(uplo,uplo)
XOPARG($1,alpha)
integer:: n,incx,incy
CHARDEF(uplo,U)
ONEDEF(alpha,$1)
INCCALC(x,$1)
INCCALC(y,$1)
n = size(x)
call $1SPR2(vuplo,n,valpha,x(1),incx,y(1),incy,ap(1))
end subroutine')
SPR2_95(s)
SPR2_95(d)

!========================================
! HER2 implementation
define(HER2_95,`
subroutine $1HER2_95(a,b,c,uplo,alpha)
$1_type,DM1,II:: a,b
$1_type,DM2,IIO:: c
BTOPARG(uplo,uplo)
XOPARG($1,alpha)
integer:: n,inca,incb,ldc
CHARDEF(uplo,U)
ONEDEF(alpha,$1)
INCCALC(a,$1)
INCCALC(b,$1)
LDCALC(c,$1)
n = size(c,1)
call $1HER2(vuplo,n,valpha,a(1),inca,b(1),incb,c(1,1),ldc)
end subroutine')
HER2_95(c)
HER2_95(z)

!========================================
! HPR2 implementation
define(HPR2_95,`
subroutine $1HPR2_95(x,y,ap,uplo,alpha)
$1_type,DM1,II:: x,y
$1_type,DM1,IIO:: ap
BTOPARG(uplo,uplo)
XOPARG($1,alpha)
integer:: n,incx,incy
CHARDEF(uplo,U)
ONEDEF(alpha,$1)
INCCALC(x,$1)
INCCALC(y,$1)
n = size(x)
call $1HPR2(vuplo,n,valpha,x(1),incx,y(1),incy,ap(1))
end subroutine')
HPR2_95(c)
HPR2_95(z)


!================================================================================
! LEVEL 3 BLAS routines
!================================================================================

!========================================
! GEMM implementation
define(GEMM_95,`
subroutine $1GEMM_95(a,b,c,transa,transb,alpha,beta)
$1_type,DM2,II:: a
$1_type,DM2,II:: b
$1_type,DM2,IIO:: c
BTOPARG(trans,transa)
BTOPARG(trans,transb)
XOPARG($1,alpha)
XOPARG($1,beta)
integer:: lda,ldb,ldc
integer:: m,n,k
CHARDEF(transa,N)
CHARDEF(transb,N)
ONEDEF(alpha,$1)
ZERODEF(beta,$1)
LDCALC(a,$1)
LDCALC(b,$1)
LDCALC(c,$1)
m = size(c,1)
n = size(c,2)
k = merge(size(a,1),size(a,2),vtransa == "T")
call $1GEMM(vtransa,vtransb,m,n,k,valpha,a(1,1),lda,b(1,1),ldb,vbeta,c(1,1),ldc)
end subroutine')

GEMM_95(s)
GEMM_95(d)
GEMM_95(c)
GEMM_95(z)

!========================================
! SYMM implementation
define(SYMM_95,`
subroutine $1SYMM_95(a,b,c,side,uplo,alpha,beta)
$1_type,DM2,II:: a
$1_type,DM2,II:: b
$1_type,DM2,IIO:: c
BTOPARG(side,side)
BTOPARG(uplo,uplo)
XOPARG($1,alpha)
XOPARG($1,beta)
integer:: lda,ldb,ldc
integer:: m,n,k
CHARDEF(side,L)
CHARDEF(uplo,U)
ONEDEF(alpha,$1)
ZERODEF(beta,$1)
LDCALC(a,$1)
LDCALC(b,$1)
LDCALC(c,$1)
m = size(c,1)
n = size(c,2)
call $1SYMM(vside,vuplo,m,n,valpha,a(1,1),lda,b(1,1),ldb,vbeta,c(1,1),ldc)
end subroutine')
SYMM_95(s)
SYMM_95(d)
SYMM_95(c)
SYMM_95(z)

!========================================
! HEMM implementation
define(HEMM_95,`
subroutine $1HEMM_95(a,b,c,side,uplo,alpha,beta)
$1_type,DM2,II:: a
$1_type,DM2,II:: b
$1_type,DM2,IIO:: c
BTOPARG(side,side)
BTOPARG(uplo,uplo)
XOPARG($1,alpha)
XOPARG($1,beta)
integer:: lda,ldb,ldc
integer:: m,n,k
CHARDEF(side,L)
CHARDEF(uplo,U)
ONEDEF(alpha,$1)
ZERODEF(beta,$1)
LDCALC(a,$1)
LDCALC(b,$1)
LDCALC(c,$1)
m = size(c,1)
n = size(c,2)
call $1HEMM(vside,vuplo,m,n,valpha,a(1,1),lda,b(1,1),ldb,vbeta,c(1,1),ldc)
end subroutine')
HEMM_95(c)
HEMM_95(z)

!========================================
! TRMM implementation
define(TRMM_95,`
subroutine $1TRMM_95(t,b,side,uplo,transt,diag,alpha)
$1_type,DM2,II:: t
$1_type,DM2,IIO:: b
BTOPARG(side,side)
BTOPARG(uplo,uplo)
BTOPARG(trans,transt)
BTOPARG(diag,diag)
XOPARG($1,alpha)
integer:: ldt,ldb
integer:: m,n
CHARDEF(side,L)
CHARDEF(transt,N)
CHARDEF(uplo,U)
CHARDEF(diag,N)
ONEDEF(alpha,$1)
LDCALC(t,$1)
LDCALC(b,$1)
m = size(b,1)
n = size(b,2)
call $1TRMM(vside,vuplo,vtranst,vdiag,m,n,valpha,t(1,1),ldt,b(1,1),ldb)
end subroutine')

TRMM_95(s)
TRMM_95(d)
TRMM_95(c)
TRMM_95(z)

!========================================
! TRSM implementation
define(TRSM_95,`
subroutine $1TRSM_95(t,b,side,uplo,transt,diag,alpha)
$1_type,DM2,II:: t
$1_type,DM2,IIO:: b
BTOPARG(side,side)
BTOPARG(uplo,uplo)
BTOPARG(trans,transt)
BTOPARG(diag,diag)
XOPARG($1,alpha)
integer:: ldt,ldb
integer:: m,n
CHARDEF(side,L)
CHARDEF(transt,N)
CHARDEF(uplo,U)
CHARDEF(diag,N)
ONEDEF(alpha,$1)
LDCALC(t,$1)
LDCALC(b,$1)
m = size(b,1)
n = size(b,2)
call $1TRSM(vside,vuplo,vtranst,vdiag,m,n,valpha,t(1,1),ldt,b(1,1),ldb)
end subroutine')

TRSM_95(s)
TRSM_95(d)
TRSM_95(c)
TRSM_95(z)

!========================================
! SYRK implementation
define(SYRK_95,`
subroutine $1SYRK_95(a,c,uplo,trans,alpha,beta)
$1_type,DM2,II:: a
$1_type,DM2,IIO:: c
BTOPARG(uplo,uplo)
BTOPARG(trans,trans)
XOPARG($1,alpha)
XOPARG($1,beta)
integer:: lda,ldc
integer:: n,k
CHARDEF(uplo,U)
CHARDEF(trans,N)
ONEDEF(alpha,$1)
ZERODEF(beta,$1)
LDCALC(a,$1)
LDCALC(c,$1)
n = size(c,1)
if (vtrans == "N") then
  k = size(a,2)
else
  k = size(a,1)
end if
call $1SYRK(vuplo,vtrans,n,k,valpha,a(1,1),lda,vbeta,c(1,1),ldc)
end subroutine')
SYRK_95(s)
SYRK_95(d)
SYRK_95(c)
SYRK_95(z)

!========================================
! HERK implementation
define(HERK_95,`
subroutine $1HERK_95(a,c,uplo,trans,alpha,beta)
$1_type,DM2,II:: a
$1_type,DM2,IIO:: c
BTOPARG(uplo,uplo)
BTOPARG(trans,trans)
XOPARG($1,alpha)
XOPARG($1,beta)
integer:: lda,ldc
integer:: n,k
CHARDEF(uplo,U)
CHARDEF(trans,N)
ONEDEF(alpha,$1)
ZERODEF(beta,$1)
LDCALC(a,$1)
LDCALC(c,$1)
n = size(c,1)
if (vtrans == "N") then
  k = size(a,2)
else
  k = size(a,1)
end if
call $1HERK(vuplo,vtrans,n,k,valpha,a(1,1),lda,vbeta,c(1,1),ldc)
end subroutine')
HERK_95(c)
HERK_95(z)

!========================================
! SYR2K implementation
define(SYR2K_95,`
subroutine $1SYR2K_95(a,b,c,uplo,trans,alpha,beta)
$1_type,DM2,II:: a,b
$1_type,DM2,IIO:: c
BTOPARG(uplo,uplo)
BTOPARG(trans,trans)
XOPARG($1,alpha)
XOPARG($1,beta)
integer:: lda,ldb,ldc
integer:: n,k
CHARDEF(uplo,U)
CHARDEF(trans,N)
ONEDEF(alpha,$1)
ZERODEF(beta,$1)
LDCALC(a,$1)
LDCALC(b,$1)
LDCALC(c,$1)
n = size(c,1)
if (vtrans == "N") then
  k = size(a,2)
else
  k = size(a,1)
end if
call $1SYR2K(vuplo,vtrans,n,k,valpha,a(1,1),lda,b(1,1),ldb,vbeta,c(1,1),ldc)
end subroutine')
SYR2K_95(s)
SYR2K_95(d)
SYR2K_95(c)
SYR2K_95(z)

!========================================
! HER2K implementation
define(HER2K_95,`
subroutine $1HER2K_95(a,b,c,uplo,trans,alpha,beta)
$1_type,DM2,II:: a,b
$1_type,DM2,IIO:: c
BTOPARG(uplo,uplo)
BTOPARG(trans,trans)
XOPARG($1,alpha)
XOPARG($1,beta)
integer:: lda,ldb,ldc
integer:: n,k
CHARDEF(uplo,U)
CHARDEF(trans,N)
ONEDEF(alpha,$1)
ZERODEF(beta,$1)
LDCALC(a,$1)
LDCALC(b,$1)
LDCALC(c,$1)
n = size(c,1)
if (vtrans == "N") then
  k = size(a,2)
else
  k = size(a,1)
end if
call $1HER2K(vuplo,vtrans,n,k,valpha,a(1,1),lda,b(1,1),ldb,vbeta,c(1,1),ldc)
end subroutine')
HER2K_95(c)
HER2K_95(z)

end module blas95