stdlib_linalg_blas_z.fypp Source File


Source Code

#:include "common.fypp" 
module stdlib_linalg_blas_z
     use stdlib_linalg_constants
     use stdlib_linalg_blas_aux
     use stdlib_linalg_blas_s
     use stdlib_linalg_blas_c
     use stdlib_linalg_blas_d
     implicit none(type,external)
     private


     public :: sp,dp,qp,lk,ilp
     public :: stdlib_zaxpy
     public :: stdlib_zcopy
     public :: stdlib_zdotc
     public :: stdlib_zdotu
     public :: stdlib_zdrot
     public :: stdlib_zdscal
     public :: stdlib_zgbmv
     public :: stdlib_zgemm
     public :: stdlib_zgemv
     public :: stdlib_zgerc
     public :: stdlib_zgeru
     public :: stdlib_zhbmv
     public :: stdlib_zhemm
     public :: stdlib_zhemv
     public :: stdlib_zher
     public :: stdlib_zher2
     public :: stdlib_zher2k
     public :: stdlib_zherk
     public :: stdlib_zhpmv
     public :: stdlib_zhpr
     public :: stdlib_zhpr2
     public :: stdlib_zrotg
     public :: stdlib_zscal
     public :: stdlib_zswap
     public :: stdlib_zsymm
     public :: stdlib_zsyr2k
     public :: stdlib_zsyrk
     public :: stdlib_ztbmv
     public :: stdlib_ztbsv
     public :: stdlib_ztpmv
     public :: stdlib_ztpsv
     public :: stdlib_ztrmm
     public :: stdlib_ztrmv
     public :: stdlib_ztrsm
     public :: stdlib_ztrsv

     ! 64-bit real constants 
     real(dp),    parameter, private ::     negone = -1.00_dp
     real(dp),    parameter, private ::       zero = 0.00_dp
     real(dp),    parameter, private ::       half = 0.50_dp
     real(dp),    parameter, private ::        one = 1.00_dp
     real(dp),    parameter, private ::        two = 2.00_dp
     real(dp),    parameter, private ::      three = 3.00_dp
     real(dp),    parameter, private ::       four = 4.00_dp
     real(dp),    parameter, private ::      eight = 8.00_dp
     real(dp),    parameter, private ::        ten = 10.00_dp

     ! 64-bit complex constants 
     complex(dp), parameter, private :: czero   = ( 0.0_dp,0.0_dp)
     complex(dp), parameter, private :: chalf   = ( 0.5_dp,0.0_dp)
     complex(dp), parameter, private :: cone    = ( 1.0_dp,0.0_dp)
     complex(dp), parameter, private :: cnegone = (-1.0_dp,0.0_dp)

     ! 64-bit scaling constants 
     integer,     parameter, private :: maxexp = maxexponent(zero) 
     integer,     parameter, private :: minexp = minexponent(zero) 
     real(dp),    parameter, private :: rradix = real(radix(zero),dp) 
     real(dp),    parameter, private :: ulp    = epsilon(zero) 
     real(dp),    parameter, private :: eps    = ulp*half 
     real(dp),    parameter, private :: safmin = rradix**max(minexp-1,1-maxexp) 
     real(dp),    parameter, private :: safmax = one/safmin 
     real(dp),    parameter, private :: smlnum = safmin/ulp 
     real(dp),    parameter, private :: bignum = safmax*ulp 
     real(dp),    parameter, private :: rtmin  = sqrt(smlnum) 
     real(dp),    parameter, private :: rtmax  = sqrt(bignum) 

     ! 64-bit Blue's scaling constants 
     ! ssml>=1/s and sbig==1/S with s,S as defined in https://doi.org/10.1145/355769.355771 
     real(dp),    parameter, private :: tsml   = rradix**ceiling((minexp-1)*half) 
     real(dp),    parameter, private :: tbig   = rradix**floor((maxexp-digits(zero)+1)*half) 
     real(dp),    parameter, private :: ssml   = rradix**(-floor((minexp-digits(zero))*half)) 
     real(dp),    parameter, private :: sbig   = rradix**(-ceiling((maxexp+digits(zero)-1)*half)) 


     contains


     pure subroutine stdlib_zaxpy(n,za,zx,incx,zy,incy)
     !! ZAXPY constant times a vector plus a vector.
        ! -- reference blas level1 routine --
        ! -- reference blas is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           ! Scalar Arguments 
           complex(dp), intent(in) :: za
           integer(ilp), intent(in) :: incx, incy, n
           ! Array Arguments 
           complex(dp), intent(in) :: zx(*)
           complex(dp), intent(inout) :: zy(*)
        ! =====================================================================
           ! Local Scalars 
           integer(ilp) :: i, ix, iy
           if (n<=0) return
           if (stdlib_cabs1(za)==0.0_dp) return
           if (incx==1 .and. incy==1) then
              ! code for both increments equal to 1
              do i = 1,n
                 zy(i) = zy(i) + za*zx(i)
              end do
           else
              ! code for unequal increments or equal increments
                ! not equal to 1
              ix = 1
              iy = 1
              if (incx<0) ix = (-n+1)*incx + 1
              if (incy<0) iy = (-n+1)*incy + 1
              do i = 1,n
                 zy(iy) = zy(iy) + za*zx(ix)
                 ix = ix + incx
                 iy = iy + incy
              end do
           end if
           return
     end subroutine stdlib_zaxpy


     pure subroutine stdlib_zcopy(n,zx,incx,zy,incy)
     !! ZCOPY copies a vector, x, to a vector, y.
        ! -- reference blas level1 routine --
        ! -- reference blas is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           ! Scalar Arguments 
           integer(ilp), intent(in) :: incx, incy, n
           ! Array Arguments 
           complex(dp), intent(in) :: zx(*)
           complex(dp), intent(out) :: zy(*)
        ! =====================================================================
           ! Local Scalars 
           integer(ilp) :: i, ix, iy
           if (n<=0) return
           if (incx==1 .and. incy==1) then
              ! code for both increments equal to 1
              do i = 1,n
               zy(i) = zx(i)
              end do
           else
              ! code for unequal increments or equal increments
                ! not equal to 1
              ix = 1
              iy = 1
              if (incx<0) ix = (-n+1)*incx + 1
              if (incy<0) iy = (-n+1)*incy + 1
              do i = 1,n
                 zy(iy) = zx(ix)
                 ix = ix + incx
                 iy = iy + incy
              end do
           end if
           return
     end subroutine stdlib_zcopy


     pure complex(dp) function stdlib_zdotc(n,zx,incx,zy,incy)
     !! ZDOTC forms the dot product of two complex vectors
     !! ZDOTC = X^H * Y
        ! -- reference blas level1 routine --
        ! -- reference blas is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           ! Scalar Arguments 
           integer(ilp), intent(in) :: incx, incy, n
           ! Array Arguments 
           complex(dp), intent(in) :: zx(*), zy(*)
        ! =====================================================================
           ! Local Scalars 
           complex(dp) :: ztemp
           integer(ilp) :: i, ix, iy
           ! Intrinsic Functions 
           intrinsic :: conjg
           ztemp = (0.0_dp,0.0_dp)
           stdlib_zdotc = (0.0_dp,0.0_dp)
           if (n<=0) return
           if (incx==1 .and. incy==1) then
              ! code for both increments equal to 1
              do i = 1,n
                 ztemp = ztemp + conjg(zx(i))*zy(i)
              end do
           else
              ! code for unequal increments or equal increments
                ! not equal to 1
              ix = 1
              iy = 1
              if (incx<0) ix = (-n+1)*incx + 1
              if (incy<0) iy = (-n+1)*incy + 1
              do i = 1,n
                 ztemp = ztemp + conjg(zx(ix))*zy(iy)
                 ix = ix + incx
                 iy = iy + incy
              end do
           end if
           stdlib_zdotc = ztemp
           return
     end function stdlib_zdotc


     pure complex(dp) function stdlib_zdotu(n,zx,incx,zy,incy)
     !! ZDOTU forms the dot product of two complex vectors
     !! ZDOTU = X^T * Y
        ! -- reference blas level1 routine --
        ! -- reference blas is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           ! Scalar Arguments 
           integer(ilp), intent(in) :: incx, incy, n
           ! Array Arguments 
           complex(dp), intent(in) :: zx(*), zy(*)
        ! =====================================================================
           ! Local Scalars 
           complex(dp) :: ztemp
           integer(ilp) :: i, ix, iy
           ztemp = (0.0_dp,0.0_dp)
           stdlib_zdotu = (0.0_dp,0.0_dp)
           if (n<=0) return
           if (incx==1 .and. incy==1) then
              ! code for both increments equal to 1
              do i = 1,n
                 ztemp = ztemp + zx(i)*zy(i)
              end do
           else
              ! code for unequal increments or equal increments
                ! not equal to 1
              ix = 1
              iy = 1
              if (incx<0) ix = (-n+1)*incx + 1
              if (incy<0) iy = (-n+1)*incy + 1
              do i = 1,n
                 ztemp = ztemp + zx(ix)*zy(iy)
                 ix = ix + incx
                 iy = iy + incy
              end do
           end if
           stdlib_zdotu = ztemp
           return
     end function stdlib_zdotu


     pure subroutine stdlib_zdrot( n, zx, incx, zy, incy, c, s )
     !! Applies a plane rotation, where the cos and sin (c and s) are real
     !! and the vectors cx and cy are complex.
     !! jack dongarra, linpack, 3/11/78.
        ! -- reference blas level1 routine --
        ! -- reference blas is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           ! Scalar Arguments 
           integer(ilp), intent(in) :: incx, incy, n
           real(dp), intent(in) :: c, s
           ! Array Arguments 
           complex(dp), intent(inout) :: zx(*), zy(*)
       ! =====================================================================
           ! Local Scalars 
           integer(ilp) :: i, ix, iy
           complex(dp) :: ctemp
           ! Executable Statements 
           if( n<=0 )return
           if( incx==1 .and. incy==1 ) then
              ! code for both increments equal to 1
              do i = 1, n
                 ctemp = c*zx( i ) + s*zy( i )
                 zy( i ) = c*zy( i ) - s*zx( i )
                 zx( i ) = ctemp
              end do
           else
              ! code for unequal increments or equal increments not equal
                ! to 1
              ix = 1
              iy = 1
              if( incx<0 )ix = ( -n+1 )*incx + 1
              if( incy<0 )iy = ( -n+1 )*incy + 1
              do i = 1, n
                 ctemp = c*zx( ix ) + s*zy( iy )
                 zy( iy ) = c*zy( iy ) - s*zx( ix )
                 zx( ix ) = ctemp
                 ix = ix + incx
                 iy = iy + incy
              end do
           end if
           return
     end subroutine stdlib_zdrot


     pure subroutine stdlib_zdscal(n,da,zx,incx)
     !! ZDSCAL scales a vector by a constant.
        ! -- reference blas level1 routine --
        ! -- reference blas is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           ! Scalar Arguments 
           real(dp), intent(in) :: da
           integer(ilp), intent(in) :: incx, n
           ! Array Arguments 
           complex(dp), intent(inout) :: zx(*)
        ! =====================================================================
           ! Local Scalars 
           integer(ilp) :: i, nincx
           ! Intrinsic Functions 
           intrinsic :: cmplx
           if (n<=0 .or. incx<=0) return
           if (incx==1) then
              ! code for increment equal to 1
              do i = 1,n
                 zx(i) = cmplx(da,0.0_dp,KIND=dp)*zx(i)
              end do
           else
              ! code for increment not equal to 1
              nincx = n*incx
              do i = 1,nincx,incx
                 zx(i) = cmplx(da,0.0_dp,KIND=dp)*zx(i)
              end do
           end if
           return
     end subroutine stdlib_zdscal


     pure subroutine stdlib_zgbmv(trans,m,n,kl,ku,alpha,a,lda,x,incx,beta,y,incy)
     !! ZGBMV performs one of the matrix-vector operations
     !! y := alpha*A*x + beta*y,   or   y := alpha*A**T*x + beta*y,   or
     !! y := alpha*A**H*x + beta*y,
     !! where alpha and beta are scalars, x and y are vectors and A is an
     !! m by n band matrix, with kl sub-diagonals and ku super-diagonals.
        ! -- reference blas level2 routine --
        ! -- reference blas is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           ! Scalar Arguments 
           complex(dp), intent(in) :: alpha, beta
           integer(ilp), intent(in) :: incx, incy, kl, ku, lda, m, n
           character, intent(in) :: trans
           ! Array Arguments 
           complex(dp), intent(in) :: a(lda,*), x(*)
           complex(dp), intent(inout) :: y(*)
        ! =====================================================================
           
           
           ! Local Scalars 
           complex(dp) :: temp
           integer(ilp) :: i, info, ix, iy, j, jx, jy, k, kup1, kx, ky, lenx, leny
           logical(lk) :: noconj
           ! Intrinsic Functions 
           intrinsic :: conjg,max,min
           ! test the input parameters.
           info = 0
           if (.not.stdlib_lsame(trans,'N') .and. .not.stdlib_lsame(trans,'T') &
                     .and..not.stdlib_lsame(trans,'C')) then
               info = 1
           else if (m<0) then
               info = 2
           else if (n<0) then
               info = 3
           else if (kl<0) then
               info = 4
           else if (ku<0) then
               info = 5
           else if (lda< (kl+ku+1)) then
               info = 8
           else if (incx==0) then
               info = 10
           else if (incy==0) then
               info = 13
           end if
           if (info/=0) then
               call stdlib_xerbla('ZGBMV ',info)
               return
           end if
           ! quick return if possible.
           if ((m==0) .or. (n==0) .or.((alpha==czero).and. (beta==cone))) return
           noconj = stdlib_lsame(trans,'T')
           ! set  lenx  and  leny, the lengths of the vectors x and y, and set
           ! up the start points in  x  and  y.
           if (stdlib_lsame(trans,'N')) then
               lenx = n
               leny = m
           else
               lenx = m
               leny = n
           end if
           if (incx>0) then
               kx = 1
           else
               kx = 1 - (lenx-1)*incx
           end if
           if (incy>0) then
               ky = 1
           else
               ky = 1 - (leny-1)*incy
           end if
           ! start the operations. in this version the elements of a are
           ! accessed sequentially with cone pass through the band part of a.
           ! first form  y := beta*y.
           if (beta/=cone) then
               if (incy==1) then
                   if (beta==czero) then
                       do i = 1,leny
                           y(i) = czero
                       end do
                   else
                       do i = 1,leny
                           y(i) = beta*y(i)
                       end do
                   end if
               else
                   iy = ky
                   if (beta==czero) then
                       do i = 1,leny
                           y(iy) = czero
                           iy = iy + incy
                       end do
                   else
                       do i = 1,leny
                           y(iy) = beta*y(iy)
                           iy = iy + incy
                       end do
                   end if
               end if
           end if
           if (alpha==czero) return
           kup1 = ku + 1
           if (stdlib_lsame(trans,'N')) then
              ! form  y := alpha*a*x + y.
               jx = kx
               if (incy==1) then
                   do j = 1,n
                       temp = alpha*x(jx)
                       k = kup1 - j
                       do i = max(1,j-ku),min(m,j+kl)
                           y(i) = y(i) + temp*a(k+i,j)
                       end do
                       jx = jx + incx
                   end do
               else
                   do j = 1,n
                       temp = alpha*x(jx)
                       iy = ky
                       k = kup1 - j
                       do i = max(1,j-ku),min(m,j+kl)
                           y(iy) = y(iy) + temp*a(k+i,j)
                           iy = iy + incy
                       end do
                       jx = jx + incx
                       if (j>ku) ky = ky + incy
                   end do
               end if
           else
              ! form  y := alpha*a**t*x + y  or  y := alpha*a**h*x + y.
               jy = ky
               if (incx==1) then
                   do j = 1,n
                       temp = czero
                       k = kup1 - j
                       if (noconj) then
                           do i = max(1,j-ku),min(m,j+kl)
                               temp = temp + a(k+i,j)*x(i)
                           end do
                       else
                           do i = max(1,j-ku),min(m,j+kl)
                               temp = temp + conjg(a(k+i,j))*x(i)
                           end do
                       end if
                       y(jy) = y(jy) + alpha*temp
                       jy = jy + incy
                   end do
               else
                   do j = 1,n
                       temp = czero
                       ix = kx
                       k = kup1 - j
                       if (noconj) then
                           do i = max(1,j-ku),min(m,j+kl)
                               temp = temp + a(k+i,j)*x(ix)
                               ix = ix + incx
                           end do
                       else
                           do i = max(1,j-ku),min(m,j+kl)
                               temp = temp + conjg(a(k+i,j))*x(ix)
                               ix = ix + incx
                           end do
                       end if
                       y(jy) = y(jy) + alpha*temp
                       jy = jy + incy
                       if (j>ku) kx = kx + incx
                   end do
               end if
           end if
           return
     end subroutine stdlib_zgbmv


     pure subroutine stdlib_zgemm(transa,transb,m,n,k,alpha,a,lda,b,ldb,beta,c,ldc)
     !! 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**T   or   op( X ) = X**H,
     !! 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.
        ! -- reference blas level3 routine --
        ! -- reference blas is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           ! Scalar Arguments 
           complex(dp), intent(in) :: alpha, beta
           integer(ilp), intent(in) :: k, lda, ldb, ldc, m, n
           character, intent(in) :: transa, transb
           ! Array Arguments 
           complex(dp), intent(in) :: a(lda,*), b(ldb,*)
           complex(dp), intent(inout) :: c(ldc,*)
        ! =====================================================================
           ! Intrinsic Functions 
           intrinsic :: conjg,max
           ! Local Scalars 
           complex(dp) :: temp
           integer(ilp) :: i, info, j, l, nrowa, nrowb
           logical(lk) :: conja, conjb, nota, notb
           
           
           ! set  nota  and  notb  as  true if  a  and  b  respectively are not
           ! conjugated or transposed, set  conja and conjb  as true if  a  and
           ! b  respectively are to be  transposed but  not conjugated  and set
           ! nrowa and nrowb  as the number of rows  of  a  and  b  respectively.
           nota = stdlib_lsame(transa,'N')
           notb = stdlib_lsame(transb,'N')
           conja = stdlib_lsame(transa,'C')
           conjb = stdlib_lsame(transb,'C')
           if (nota) then
               nrowa = m
           else
               nrowa = k
           end if
           if (notb) then
               nrowb = k
           else
               nrowb = n
           end if
           ! test the input parameters.
           info = 0
           if ((.not.nota) .and. (.not.conja) .and.(.not.stdlib_lsame(transa,'T'))) then
               info = 1
           else if ((.not.notb) .and. (.not.conjb) .and.(.not.stdlib_lsame(transb,'T'))) &
                     then
               info = 2
           else if (m<0) then
               info = 3
           else if (n<0) then
               info = 4
           else if (k<0) then
               info = 5
           else if (lda<max(1,nrowa)) then
               info = 8
           else if (ldb<max(1,nrowb)) then
               info = 10
           else if (ldc<max(1,m)) then
               info = 13
           end if
           if (info/=0) then
               call stdlib_xerbla('ZGEMM ',info)
               return
           end if
           ! quick return if possible.
           if ((m==0) .or. (n==0) .or.(((alpha==czero).or. (k==0)).and. (beta==cone))) &
                     return
           ! and when  alpha.eq.czero.
           if (alpha==czero) then
               if (beta==czero) then
                   do j = 1,n
                       do i = 1,m
                           c(i,j) = czero
                       end do
                   end do
               else
                   do j = 1,n
                       do i = 1,m
                           c(i,j) = beta*c(i,j)
                       end do
                   end do
               end if
               return
           end if
           ! start the operations.
           if (notb) then
               if (nota) then
                 ! form  c := alpha*a*b + beta*c.
                   do j = 1,n
                       if (beta==czero) then
                           do i = 1,m
                               c(i,j) = czero
                           end do
                       else if (beta/=cone) then
                           do i = 1,m
                               c(i,j) = beta*c(i,j)
                           end do
                       end if
                       do l = 1,k
                           temp = alpha*b(l,j)
                           do i = 1,m
                               c(i,j) = c(i,j) + temp*a(i,l)
                           end do
                       end do
                   end do
               else if (conja) then
                 ! form  c := alpha*a**h*b + beta*c.
                   do j = 1,n
                       do i = 1,m
                           temp = czero
                           do l = 1,k
                               temp = temp + conjg(a(l,i))*b(l,j)
                           end do
                           if (beta==czero) then
                               c(i,j) = alpha*temp
                           else
                               c(i,j) = alpha*temp + beta*c(i,j)
                           end if
                       end do
                   end do
               else
                 ! form  c := alpha*a**t*b + beta*c
                   do j = 1,n
                       do i = 1,m
                           temp = czero
                           do l = 1,k
                               temp = temp + a(l,i)*b(l,j)
                           end do
                           if (beta==czero) then
                               c(i,j) = alpha*temp
                           else
                               c(i,j) = alpha*temp + beta*c(i,j)
                           end if
                       end do
                   end do
               end if
           else if (nota) then
               if (conjb) then
                 ! form  c := alpha*a*b**h + beta*c.
                   do j = 1,n
                       if (beta==czero) then
                           do i = 1,m
                               c(i,j) = czero
                           end do
                       else if (beta/=cone) then
                           do i = 1,m
                               c(i,j) = beta*c(i,j)
                           end do
                       end if
                       do l = 1,k
                           temp = alpha*conjg(b(j,l))
                           do i = 1,m
                               c(i,j) = c(i,j) + temp*a(i,l)
                           end do
                       end do
                   end do
               else
                 ! form  c := alpha*a*b**t + beta*c
                   do j = 1,n
                       if (beta==czero) then
                           do i = 1,m
                               c(i,j) = czero
                           end do
                       else if (beta/=cone) then
                           do i = 1,m
                               c(i,j) = beta*c(i,j)
                           end do
                       end if
                       do l = 1,k
                           temp = alpha*b(j,l)
                           do i = 1,m
                               c(i,j) = c(i,j) + temp*a(i,l)
                           end do
                       end do
                   end do
               end if
           else if (conja) then
               if (conjb) then
                 ! form  c := alpha*a**h*b**h + beta*c.
                   do j = 1,n
                       do i = 1,m
                           temp = czero
                           do l = 1,k
                               temp = temp + conjg(a(l,i))*conjg(b(j,l))
                           end do
                           if (beta==czero) then
                               c(i,j) = alpha*temp
                           else
                               c(i,j) = alpha*temp + beta*c(i,j)
                           end if
                       end do
                   end do
               else
                 ! form  c := alpha*a**h*b**t + beta*c
                   do j = 1,n
                       do i = 1,m
                           temp = czero
                           do l = 1,k
                               temp = temp + conjg(a(l,i))*b(j,l)
                           end do
                           if (beta==czero) then
                               c(i,j) = alpha*temp
                           else
                               c(i,j) = alpha*temp + beta*c(i,j)
                           end if
                       end do
                   end do
               end if
           else
               if (conjb) then
                 ! form  c := alpha*a**t*b**h + beta*c
                   do j = 1,n
                       do i = 1,m
                           temp = czero
                           do l = 1,k
                               temp = temp + a(l,i)*conjg(b(j,l))
                           end do
                           if (beta==czero) then
                               c(i,j) = alpha*temp
                           else
                               c(i,j) = alpha*temp + beta*c(i,j)
                           end if
                       end do
                   end do
               else
                 ! form  c := alpha*a**t*b**t + beta*c
                   do j = 1,n
                       do i = 1,m
                           temp = czero
                           do l = 1,k
                               temp = temp + a(l,i)*b(j,l)
                           end do
                           if (beta==czero) then
                               c(i,j) = alpha*temp
                           else
                               c(i,j) = alpha*temp + beta*c(i,j)
                           end if
                       end do
                   end do
               end if
           end if
           return
     end subroutine stdlib_zgemm


     pure subroutine stdlib_zgemv(trans,m,n,alpha,a,lda,x,incx,beta,y,incy)
     !! ZGEMV performs one of the matrix-vector operations
     !! y := alpha*A*x + beta*y,   or   y := alpha*A**T*x + beta*y,   or
     !! y := alpha*A**H*x + beta*y,
     !! where alpha and beta are scalars, x and y are vectors and A is an
     !! m by n matrix.
        ! -- reference blas level2 routine --
        ! -- reference blas is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           ! Scalar Arguments 
           complex(dp), intent(in) :: alpha, beta
           integer(ilp), intent(in) :: incx, incy, lda, m, n
           character, intent(in) :: trans
           ! Array Arguments 
           complex(dp), intent(in) :: a(lda,*), x(*)
           complex(dp), intent(inout) :: y(*)
        ! =====================================================================
           
           
           ! Local Scalars 
           complex(dp) :: temp
           integer(ilp) :: i, info, ix, iy, j, jx, jy, kx, ky, lenx, leny
           logical(lk) :: noconj
           ! Intrinsic Functions 
           intrinsic :: conjg,max
           ! test the input parameters.
           info = 0
           if (.not.stdlib_lsame(trans,'N') .and. .not.stdlib_lsame(trans,'T') &
                     .and..not.stdlib_lsame(trans,'C')) then
               info = 1
           else if (m<0) then
               info = 2
           else if (n<0) then
               info = 3
           else if (lda<max(1,m)) then
               info = 6
           else if (incx==0) then
               info = 8
           else if (incy==0) then
               info = 11
           end if
           if (info/=0) then
               call stdlib_xerbla('ZGEMV ',info)
               return
           end if
           ! quick return if possible.
           if ((m==0) .or. (n==0) .or.((alpha==czero).and. (beta==cone))) return
           noconj = stdlib_lsame(trans,'T')
           ! set  lenx  and  leny, the lengths of the vectors x and y, and set
           ! up the start points in  x  and  y.
           if (stdlib_lsame(trans,'N')) then
               lenx = n
               leny = m
           else
               lenx = m
               leny = n
           end if
           if (incx>0) then
               kx = 1
           else
               kx = 1 - (lenx-1)*incx
           end if
           if (incy>0) then
               ky = 1
           else
               ky = 1 - (leny-1)*incy
           end if
           ! start the operations. in this version the elements of a are
           ! accessed sequentially with cone pass through a.
           ! first form  y := beta*y.
           if (beta/=cone) then
               if (incy==1) then
                   if (beta==czero) then
                       do i = 1,leny
                           y(i) = czero
                       end do
                   else
                       do i = 1,leny
                           y(i) = beta*y(i)
                       end do
                   end if
               else
                   iy = ky
                   if (beta==czero) then
                       do i = 1,leny
                           y(iy) = czero
                           iy = iy + incy
                       end do
                   else
                       do i = 1,leny
                           y(iy) = beta*y(iy)
                           iy = iy + incy
                       end do
                   end if
               end if
           end if
           if (alpha==czero) return
           if (stdlib_lsame(trans,'N')) then
              ! form  y := alpha*a*x + y.
               jx = kx
               if (incy==1) then
                   do j = 1,n
                       temp = alpha*x(jx)
                       do i = 1,m
                           y(i) = y(i) + temp*a(i,j)
                       end do
                       jx = jx + incx
                   end do
               else
                   do j = 1,n
                       temp = alpha*x(jx)
                       iy = ky
                       do i = 1,m
                           y(iy) = y(iy) + temp*a(i,j)
                           iy = iy + incy
                       end do
                       jx = jx + incx
                   end do
               end if
           else
              ! form  y := alpha*a**t*x + y  or  y := alpha*a**h*x + y.
               jy = ky
               if (incx==1) then
                   do j = 1,n
                       temp = czero
                       if (noconj) then
                           do i = 1,m
                               temp = temp + a(i,j)*x(i)
                           end do
                       else
                           do i = 1,m
                               temp = temp + conjg(a(i,j))*x(i)
                           end do
                       end if
                       y(jy) = y(jy) + alpha*temp
                       jy = jy + incy
                   end do
               else
                   do j = 1,n
                       temp = czero
                       ix = kx
                       if (noconj) then
                           do i = 1,m
                               temp = temp + a(i,j)*x(ix)
                               ix = ix + incx
                           end do
                       else
                           do i = 1,m
                               temp = temp + conjg(a(i,j))*x(ix)
                               ix = ix + incx
                           end do
                       end if
                       y(jy) = y(jy) + alpha*temp
                       jy = jy + incy
                   end do
               end if
           end if
           return
     end subroutine stdlib_zgemv


     pure subroutine stdlib_zgerc(m,n,alpha,x,incx,y,incy,a,lda)
     !! ZGERC performs the rank 1 operation
     !! A := alpha*x*y**H + 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.
        ! -- reference blas level2 routine --
        ! -- reference blas is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           ! Scalar Arguments 
           complex(dp), intent(in) :: alpha
           integer(ilp), intent(in) :: incx, incy, lda, m, n
           ! Array Arguments 
           complex(dp), intent(inout) :: a(lda,*)
           complex(dp), intent(in) :: x(*), y(*)
        ! =====================================================================
           
           ! Local Scalars 
           complex(dp) :: temp
           integer(ilp) :: i, info, ix, j, jy, kx
           ! Intrinsic Functions 
           intrinsic :: conjg,max
           ! test the input parameters.
           info = 0
           if (m<0) then
               info = 1
           else if (n<0) then
               info = 2
           else if (incx==0) then
               info = 5
           else if (incy==0) then
               info = 7
           else if (lda<max(1,m)) then
               info = 9
           end if
           if (info/=0) then
               call stdlib_xerbla('ZGERC ',info)
               return
           end if
           ! quick return if possible.
           if ((m==0) .or. (n==0) .or. (alpha==czero)) return
           ! start the operations. in this version the elements of a are
           ! accessed sequentially with cone pass through a.
           if (incy>0) then
               jy = 1
           else
               jy = 1 - (n-1)*incy
           end if
           if (incx==1) then
               do j = 1,n
                   if (y(jy)/=czero) then
                       temp = alpha*conjg(y(jy))
                       do i = 1,m
                           a(i,j) = a(i,j) + x(i)*temp
                       end do
                   end if
                   jy = jy + incy
               end do
           else
               if (incx>0) then
                   kx = 1
               else
                   kx = 1 - (m-1)*incx
               end if
               do j = 1,n
                   if (y(jy)/=czero) then
                       temp = alpha*conjg(y(jy))
                       ix = kx
                       do i = 1,m
                           a(i,j) = a(i,j) + x(ix)*temp
                           ix = ix + incx
                       end do
                   end if
                   jy = jy + incy
               end do
           end if
           return
     end subroutine stdlib_zgerc


     pure subroutine stdlib_zgeru(m,n,alpha,x,incx,y,incy,a,lda)
     !! ZGERU performs the rank 1 operation
     !! A := alpha*x*y**T + 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.
        ! -- reference blas level2 routine --
        ! -- reference blas is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           ! Scalar Arguments 
           complex(dp), intent(in) :: alpha
           integer(ilp), intent(in) :: incx, incy, lda, m, n
           ! Array Arguments 
           complex(dp), intent(inout) :: a(lda,*)
           complex(dp), intent(in) :: x(*), y(*)
        ! =====================================================================
           
           ! Local Scalars 
           complex(dp) :: temp
           integer(ilp) :: i, info, ix, j, jy, kx
           ! Intrinsic Functions 
           intrinsic :: max
           ! test the input parameters.
           info = 0
           if (m<0) then
               info = 1
           else if (n<0) then
               info = 2
           else if (incx==0) then
               info = 5
           else if (incy==0) then
               info = 7
           else if (lda<max(1,m)) then
               info = 9
           end if
           if (info/=0) then
               call stdlib_xerbla('ZGERU ',info)
               return
           end if
           ! quick return if possible.
           if ((m==0) .or. (n==0) .or. (alpha==czero)) return
           ! start the operations. in this version the elements of a are
           ! accessed sequentially with cone pass through a.
           if (incy>0) then
               jy = 1
           else
               jy = 1 - (n-1)*incy
           end if
           if (incx==1) then
               do j = 1,n
                   if (y(jy)/=czero) then
                       temp = alpha*y(jy)
                       do i = 1,m
                           a(i,j) = a(i,j) + x(i)*temp
                       end do
                   end if
                   jy = jy + incy
               end do
           else
               if (incx>0) then
                   kx = 1
               else
                   kx = 1 - (m-1)*incx
               end if
               do j = 1,n
                   if (y(jy)/=czero) then
                       temp = alpha*y(jy)
                       ix = kx
                       do i = 1,m
                           a(i,j) = a(i,j) + x(ix)*temp
                           ix = ix + incx
                       end do
                   end if
                   jy = jy + incy
               end do
           end if
           return
     end subroutine stdlib_zgeru


     pure subroutine stdlib_zhbmv(uplo,n,k,alpha,a,lda,x,incx,beta,y,incy)
     !! ZHBMV performs the matrix-vector  operation
     !! y := alpha*A*x + beta*y,
     !! where alpha and beta are scalars, x and y are n element vectors and
     !! A is an n by n hermitian band matrix, with k super-diagonals.
        ! -- reference blas level2 routine --
        ! -- reference blas is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           ! Scalar Arguments 
           complex(dp), intent(in) :: alpha, beta
           integer(ilp), intent(in) :: incx, incy, k, lda, n
           character, intent(in) :: uplo
           ! Array Arguments 
           complex(dp), intent(in) :: a(lda,*), x(*)
           complex(dp), intent(inout) :: y(*)
        ! =====================================================================
           
           
           ! Local Scalars 
           complex(dp) :: temp1, temp2
           integer(ilp) :: i, info, ix, iy, j, jx, jy, kplus1, kx, ky, l
           ! Intrinsic Functions 
           intrinsic :: real,conjg,max,min
           ! test the input parameters.
           info = 0
           if (.not.stdlib_lsame(uplo,'U') .and. .not.stdlib_lsame(uplo,'L')) then
               info = 1
           else if (n<0) then
               info = 2
           else if (k<0) then
               info = 3
           else if (lda< (k+1)) then
               info = 6
           else if (incx==0) then
               info = 8
           else if (incy==0) then
               info = 11
           end if
           if (info/=0) then
               call stdlib_xerbla('ZHBMV ',info)
               return
           end if
           ! quick return if possible.
           if ((n==0) .or. ((alpha==czero).and. (beta==cone))) return
           ! set up the start points in  x  and  y.
           if (incx>0) then
               kx = 1
           else
               kx = 1 - (n-1)*incx
           end if
           if (incy>0) then
               ky = 1
           else
               ky = 1 - (n-1)*incy
           end if
           ! start the operations. in this version the elements of the array a
           ! are accessed sequentially with cone pass through a.
           ! first form  y := beta*y.
           if (beta/=cone) then
               if (incy==1) then
                   if (beta==czero) then
                       do i = 1,n
                           y(i) = czero
                       end do
                   else
                       do i = 1,n
                           y(i) = beta*y(i)
                       end do
                   end if
               else
                   iy = ky
                   if (beta==czero) then
                       do i = 1,n
                           y(iy) = czero
                           iy = iy + incy
                       end do
                   else
                       do i = 1,n
                           y(iy) = beta*y(iy)
                           iy = iy + incy
                       end do
                   end if
               end if
           end if
           if (alpha==czero) return
           if (stdlib_lsame(uplo,'U')) then
              ! form  y  when upper triangle of a is stored.
               kplus1 = k + 1
               if ((incx==1) .and. (incy==1)) then
                   do j = 1,n
                       temp1 = alpha*x(j)
                       temp2 = czero
                       l = kplus1 - j
                       do i = max(1,j-k),j - 1
                           y(i) = y(i) + temp1*a(l+i,j)
                           temp2 = temp2 + conjg(a(l+i,j))*x(i)
                       end do
                       y(j) = y(j) + temp1*real(a(kplus1,j),KIND=dp) + alpha*temp2
                   end do
               else
                   jx = kx
                   jy = ky
                   do j = 1,n
                       temp1 = alpha*x(jx)
                       temp2 = czero
                       ix = kx
                       iy = ky
                       l = kplus1 - j
                       do i = max(1,j-k),j - 1
                           y(iy) = y(iy) + temp1*a(l+i,j)
                           temp2 = temp2 + conjg(a(l+i,j))*x(ix)
                           ix = ix + incx
                           iy = iy + incy
                       end do
                       y(jy) = y(jy) + temp1*real(a(kplus1,j),KIND=dp) + alpha*temp2
                       jx = jx + incx
                       jy = jy + incy
                       if (j>k) then
                           kx = kx + incx
                           ky = ky + incy
                       end if
                   end do
               end if
           else
              ! form  y  when lower triangle of a is stored.
               if ((incx==1) .and. (incy==1)) then
                   do j = 1,n
                       temp1 = alpha*x(j)
                       temp2 = czero
                       y(j) = y(j) + temp1*real(a(1,j),KIND=dp)
                       l = 1 - j
                       do i = j + 1,min(n,j+k)
                           y(i) = y(i) + temp1*a(l+i,j)
                           temp2 = temp2 + conjg(a(l+i,j))*x(i)
                       end do
                       y(j) = y(j) + alpha*temp2
                   end do
               else
                   jx = kx
                   jy = ky
                   do j = 1,n
                       temp1 = alpha*x(jx)
                       temp2 = czero
                       y(jy) = y(jy) + temp1*real(a(1,j),KIND=dp)
                       l = 1 - j
                       ix = jx
                       iy = jy
                       do i = j + 1,min(n,j+k)
                           ix = ix + incx
                           iy = iy + incy
                           y(iy) = y(iy) + temp1*a(l+i,j)
                           temp2 = temp2 + conjg(a(l+i,j))*x(ix)
                       end do
                       y(jy) = y(jy) + alpha*temp2
                       jx = jx + incx
                       jy = jy + incy
                   end do
               end if
           end if
           return
     end subroutine stdlib_zhbmv


     pure subroutine stdlib_zhemm(side,uplo,m,n,alpha,a,lda,b,ldb,beta,c,ldc)
     !! ZHEMM performs one of the matrix-matrix operations
     !! C := alpha*A*B + beta*C,
     !! or
     !! C := alpha*B*A + beta*C,
     !! where alpha and beta are scalars, A is an hermitian matrix and  B and
     !! C are m by n matrices.
        ! -- reference blas level3 routine --
        ! -- reference blas is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           ! Scalar Arguments 
           complex(dp), intent(in) :: alpha, beta
           integer(ilp), intent(in) :: lda, ldb, ldc, m, n
           character, intent(in) :: side, uplo
           ! Array Arguments 
           complex(dp), intent(in) :: a(lda,*), b(ldb,*)
           complex(dp), intent(inout) :: c(ldc,*)
        ! =====================================================================
           ! Intrinsic Functions 
           intrinsic :: real,conjg,max
           ! Local Scalars 
           complex(dp) :: temp1, temp2
           integer(ilp) :: i, info, j, k, nrowa
           logical(lk) :: upper
           
           
           ! set nrowa as the number of rows of a.
           if (stdlib_lsame(side,'L')) then
               nrowa = m
           else
               nrowa = n
           end if
           upper = stdlib_lsame(uplo,'U')
           ! test the input parameters.
           info = 0
           if ((.not.stdlib_lsame(side,'L')) .and. (.not.stdlib_lsame(side,'R'))) then
               info = 1
           else if ((.not.upper) .and. (.not.stdlib_lsame(uplo,'L'))) then
               info = 2
           else if (m<0) then
               info = 3
           else if (n<0) then
               info = 4
           else if (lda<max(1,nrowa)) then
               info = 7
           else if (ldb<max(1,m)) then
               info = 9
           else if (ldc<max(1,m)) then
               info = 12
           end if
           if (info/=0) then
               call stdlib_xerbla('ZHEMM ',info)
               return
           end if
           ! quick return if possible.
           if ((m==0) .or. (n==0) .or.((alpha==czero).and. (beta==cone))) return
           ! and when  alpha.eq.czero.
           if (alpha==czero) then
               if (beta==czero) then
                   do j = 1,n
                       do i = 1,m
                           c(i,j) = czero
                       end do
                   end do
               else
                   do j = 1,n
                       do i = 1,m
                           c(i,j) = beta*c(i,j)
                       end do
                   end do
               end if
               return
           end if
           ! start the operations.
           if (stdlib_lsame(side,'L')) then
              ! form  c := alpha*a*b + beta*c.
               if (upper) then
                   do j = 1,n
                       do i = 1,m
                           temp1 = alpha*b(i,j)
                           temp2 = czero
                           do k = 1,i - 1
                               c(k,j) = c(k,j) + temp1*a(k,i)
                               temp2 = temp2 + b(k,j)*conjg(a(k,i))
                           end do
                           if (beta==czero) then
                               c(i,j) = temp1*real(a(i,i),KIND=dp) + alpha*temp2
                           else
                               c(i,j) = beta*c(i,j) + temp1*real(a(i,i),KIND=dp) +&
                                         alpha*temp2
                           end if
                       end do
                   end do
               else
                   do j = 1,n
                       do i = m,1,-1
                           temp1 = alpha*b(i,j)
                           temp2 = czero
                           do k = i + 1,m
                               c(k,j) = c(k,j) + temp1*a(k,i)
                               temp2 = temp2 + b(k,j)*conjg(a(k,i))
                           end do
                           if (beta==czero) then
                               c(i,j) = temp1*real(a(i,i),KIND=dp) + alpha*temp2
                           else
                               c(i,j) = beta*c(i,j) + temp1*real(a(i,i),KIND=dp) +&
                                         alpha*temp2
                           end if
                       end do
                   end do
               end if
           else
              ! form  c := alpha*b*a + beta*c.
               loop_170: do j = 1,n
                   temp1 = alpha*real(a(j,j),KIND=dp)
                   if (beta==czero) then
                       do i = 1,m
                           c(i,j) = temp1*b(i,j)
                       end do
                   else
                       do i = 1,m
                           c(i,j) = beta*c(i,j) + temp1*b(i,j)
                       end do
                   end if
                   do k = 1,j - 1
                       if (upper) then
                           temp1 = alpha*a(k,j)
                       else
                           temp1 = alpha*conjg(a(j,k))
                       end if
                       do i = 1,m
                           c(i,j) = c(i,j) + temp1*b(i,k)
                       end do
                   end do
                   do k = j + 1,n
                       if (upper) then
                           temp1 = alpha*conjg(a(j,k))
                       else
                           temp1 = alpha*a(k,j)
                       end if
                       do i = 1,m
                           c(i,j) = c(i,j) + temp1*b(i,k)
                       end do
                   end do
               end do loop_170
           end if
           return
     end subroutine stdlib_zhemm


     pure subroutine stdlib_zhemv(uplo,n,alpha,a,lda,x,incx,beta,y,incy)
     !! ZHEMV performs the matrix-vector  operation
     !! y := alpha*A*x + beta*y,
     !! where alpha and beta are scalars, x and y are n element vectors and
     !! A is an n by n hermitian matrix.
        ! -- reference blas level2 routine --
        ! -- reference blas is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           ! Scalar Arguments 
           complex(dp), intent(in) :: alpha, beta
           integer(ilp), intent(in) :: incx, incy, lda, n
           character, intent(in) :: uplo
           ! Array Arguments 
           complex(dp), intent(in) :: a(lda,*), x(*)
           complex(dp), intent(inout) :: y(*)
        ! =====================================================================
           
           
           ! Local Scalars 
           complex(dp) :: temp1, temp2
           integer(ilp) :: i, info, ix, iy, j, jx, jy, kx, ky
           ! Intrinsic Functions 
           intrinsic :: real,conjg,max
           ! test the input parameters.
           info = 0
           if (.not.stdlib_lsame(uplo,'U') .and. .not.stdlib_lsame(uplo,'L')) then
               info = 1
           else if (n<0) then
               info = 2
           else if (lda<max(1,n)) then
               info = 5
           else if (incx==0) then
               info = 7
           else if (incy==0) then
               info = 10
           end if
           if (info/=0) then
               call stdlib_xerbla('ZHEMV ',info)
               return
           end if
           ! quick return if possible.
           if ((n==0) .or. ((alpha==czero).and. (beta==cone))) return
           ! set up the start points in  x  and  y.
           if (incx>0) then
               kx = 1
           else
               kx = 1 - (n-1)*incx
           end if
           if (incy>0) then
               ky = 1
           else
               ky = 1 - (n-1)*incy
           end if
           ! start the operations. in this version the elements of a are
           ! accessed sequentially with cone pass through the triangular part
           ! of a.
           ! first form  y := beta*y.
           if (beta/=cone) then
               if (incy==1) then
                   if (beta==czero) then
                       do i = 1,n
                           y(i) = czero
                       end do
                   else
                       do i = 1,n
                           y(i) = beta*y(i)
                       end do
                   end if
               else
                   iy = ky
                   if (beta==czero) then
                       do i = 1,n
                           y(iy) = czero
                           iy = iy + incy
                       end do
                   else
                       do i = 1,n
                           y(iy) = beta*y(iy)
                           iy = iy + incy
                       end do
                   end if
               end if
           end if
           if (alpha==czero) return
           if (stdlib_lsame(uplo,'U')) then
              ! form  y  when a is stored in upper triangle.
               if ((incx==1) .and. (incy==1)) then
                   do j = 1,n
                       temp1 = alpha*x(j)
                       temp2 = czero
                       do i = 1,j - 1
                           y(i) = y(i) + temp1*a(i,j)
                           temp2 = temp2 + conjg(a(i,j))*x(i)
                       end do
                       y(j) = y(j) + temp1*real(a(j,j),KIND=dp) + alpha*temp2
                   end do
               else
                   jx = kx
                   jy = ky
                   do j = 1,n
                       temp1 = alpha*x(jx)
                       temp2 = czero
                       ix = kx
                       iy = ky
                       do i = 1,j - 1
                           y(iy) = y(iy) + temp1*a(i,j)
                           temp2 = temp2 + conjg(a(i,j))*x(ix)
                           ix = ix + incx
                           iy = iy + incy
                       end do
                       y(jy) = y(jy) + temp1*real(a(j,j),KIND=dp) + alpha*temp2
                       jx = jx + incx
                       jy = jy + incy
                   end do
               end if
           else
              ! form  y  when a is stored in lower triangle.
               if ((incx==1) .and. (incy==1)) then
                   do j = 1,n
                       temp1 = alpha*x(j)
                       temp2 = czero
                       y(j) = y(j) + temp1*real(a(j,j),KIND=dp)
                       do i = j + 1,n
                           y(i) = y(i) + temp1*a(i,j)
                           temp2 = temp2 + conjg(a(i,j))*x(i)
                       end do
                       y(j) = y(j) + alpha*temp2
                   end do
               else
                   jx = kx
                   jy = ky
                   do j = 1,n
                       temp1 = alpha*x(jx)
                       temp2 = czero
                       y(jy) = y(jy) + temp1*real(a(j,j),KIND=dp)
                       ix = jx
                       iy = jy
                       do i = j + 1,n
                           ix = ix + incx
                           iy = iy + incy
                           y(iy) = y(iy) + temp1*a(i,j)
                           temp2 = temp2 + conjg(a(i,j))*x(ix)
                       end do
                       y(jy) = y(jy) + alpha*temp2
                       jx = jx + incx
                       jy = jy + incy
                   end do
               end if
           end if
           return
     end subroutine stdlib_zhemv


     pure subroutine stdlib_zher(uplo,n,alpha,x,incx,a,lda)
     !! ZHER performs the hermitian rank 1 operation
     !! A := alpha*x*x**H + A,
     !! where alpha is a real scalar, x is an n element vector and A is an
     !! n by n hermitian matrix.
        ! -- reference blas level2 routine --
        ! -- reference blas is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           ! Scalar Arguments 
           real(dp), intent(in) :: alpha
           integer(ilp), intent(in) :: incx, lda, n
           character, intent(in) :: uplo
           ! Array Arguments 
           complex(dp), intent(inout) :: a(lda,*)
           complex(dp), intent(in) :: x(*)
        ! =====================================================================
           
           ! Local Scalars 
           complex(dp) :: temp
           integer(ilp) :: i, info, ix, j, jx, kx
           ! Intrinsic Functions 
           intrinsic :: real,conjg,max
           ! test the input parameters.
           info = 0
           if (.not.stdlib_lsame(uplo,'U') .and. .not.stdlib_lsame(uplo,'L')) then
               info = 1
           else if (n<0) then
               info = 2
           else if (incx==0) then
               info = 5
           else if (lda<max(1,n)) then
               info = 7
           end if
           if (info/=0) then
               call stdlib_xerbla('ZHER  ',info)
               return
           end if
           ! quick return if possible.
           if ((n==0) .or. (alpha==real(czero,KIND=dp))) return
           ! set the start point in x if the increment is not unity.
           if (incx<=0) then
               kx = 1 - (n-1)*incx
           else if (incx/=1) then
               kx = 1
           end if
           ! start the operations. in this version the elements of a are
           ! accessed sequentially with cone pass through the triangular part
           ! of a.
           if (stdlib_lsame(uplo,'U')) then
              ! form  a  when a is stored in upper triangle.
               if (incx==1) then
                   do j = 1,n
                       if (x(j)/=czero) then
                           temp = alpha*conjg(x(j))
                           do i = 1,j - 1
                               a(i,j) = a(i,j) + x(i)*temp
                           end do
                           a(j,j) = real(a(j,j),KIND=dp) + real(x(j)*temp,KIND=dp)
                       else
                           a(j,j) = real(a(j,j),KIND=dp)
                       end if
                   end do
               else
                   jx = kx
                   do j = 1,n
                       if (x(jx)/=czero) then
                           temp = alpha*conjg(x(jx))
                           ix = kx
                           do i = 1,j - 1
                               a(i,j) = a(i,j) + x(ix)*temp
                               ix = ix + incx
                           end do
                           a(j,j) = real(a(j,j),KIND=dp) + real(x(jx)*temp,KIND=dp)
                       else
                           a(j,j) = real(a(j,j),KIND=dp)
                       end if
                       jx = jx + incx
                   end do
               end if
           else
              ! form  a  when a is stored in lower triangle.
               if (incx==1) then
                   do j = 1,n
                       if (x(j)/=czero) then
                           temp = alpha*conjg(x(j))
                           a(j,j) = real(a(j,j),KIND=dp) + real(temp*x(j),KIND=dp)
                           do i = j + 1,n
                               a(i,j) = a(i,j) + x(i)*temp
                           end do
                       else
                           a(j,j) = real(a(j,j),KIND=dp)
                       end if
                   end do
               else
                   jx = kx
                   do j = 1,n
                       if (x(jx)/=czero) then
                           temp = alpha*conjg(x(jx))
                           a(j,j) = real(a(j,j),KIND=dp) + real(temp*x(jx),KIND=dp)
                           ix = jx
                           do i = j + 1,n
                               ix = ix + incx
                               a(i,j) = a(i,j) + x(ix)*temp
                           end do
                       else
                           a(j,j) = real(a(j,j),KIND=dp)
                       end if
                       jx = jx + incx
                   end do
               end if
           end if
           return
     end subroutine stdlib_zher


     pure subroutine stdlib_zher2(uplo,n,alpha,x,incx,y,incy,a,lda)
     !! ZHER2 performs the hermitian rank 2 operation
     !! A := alpha*x*y**H + conjg( alpha )*y*x**H + A,
     !! where alpha is a scalar, x and y are n element vectors and A is an n
     !! by n hermitian matrix.
        ! -- reference blas level2 routine --
        ! -- reference blas is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           ! Scalar Arguments 
           complex(dp), intent(in) :: alpha
           integer(ilp), intent(in) :: incx, incy, lda, n
           character, intent(in) :: uplo
           ! Array Arguments 
           complex(dp), intent(inout) :: a(lda,*)
           complex(dp), intent(in) :: x(*), y(*)
        ! =====================================================================
           
           ! Local Scalars 
           complex(dp) :: temp1, temp2
           integer(ilp) :: i, info, ix, iy, j, jx, jy, kx, ky
           ! Intrinsic Functions 
           intrinsic :: real,conjg,max
           ! test the input parameters.
           info = 0
           if (.not.stdlib_lsame(uplo,'U') .and. .not.stdlib_lsame(uplo,'L')) then
               info = 1
           else if (n<0) then
               info = 2
           else if (incx==0) then
               info = 5
           else if (incy==0) then
               info = 7
           else if (lda<max(1,n)) then
               info = 9
           end if
           if (info/=0) then
               call stdlib_xerbla('ZHER2 ',info)
               return
           end if
           ! quick return if possible.
           if ((n==0) .or. (alpha==czero)) return
           ! set up the start points in x and y if the increments are not both
           ! unity.
           if ((incx/=1) .or. (incy/=1)) then
               if (incx>0) then
                   kx = 1
               else
                   kx = 1 - (n-1)*incx
               end if
               if (incy>0) then
                   ky = 1
               else
                   ky = 1 - (n-1)*incy
               end if
               jx = kx
               jy = ky
           end if
           ! start the operations. in this version the elements of a are
           ! accessed sequentially with cone pass through the triangular part
           ! of a.
           if (stdlib_lsame(uplo,'U')) then
              ! form  a  when a is stored in the upper triangle.
               if ((incx==1) .and. (incy==1)) then
                   do j = 1,n
                       if ((x(j)/=czero) .or. (y(j)/=czero)) then
                           temp1 = alpha*conjg(y(j))
                           temp2 = conjg(alpha*x(j))
                           do i = 1,j - 1
                               a(i,j) = a(i,j) + x(i)*temp1 + y(i)*temp2
                           end do
                           a(j,j) = real(a(j,j),KIND=dp) +real(x(j)*temp1+y(j)*temp2,KIND=dp)
                                     
                       else
                           a(j,j) = real(a(j,j),KIND=dp)
                       end if
                   end do
               else
                   do j = 1,n
                       if ((x(jx)/=czero) .or. (y(jy)/=czero)) then
                           temp1 = alpha*conjg(y(jy))
                           temp2 = conjg(alpha*x(jx))
                           ix = kx
                           iy = ky
                           do i = 1,j - 1
                               a(i,j) = a(i,j) + x(ix)*temp1 + y(iy)*temp2
                               ix = ix + incx
                               iy = iy + incy
                           end do
                           a(j,j) = real(a(j,j),KIND=dp) +real(x(jx)*temp1+y(jy)*temp2,KIND=dp)
                                     
                       else
                           a(j,j) = real(a(j,j),KIND=dp)
                       end if
                       jx = jx + incx
                       jy = jy + incy
                   end do
               end if
           else
              ! form  a  when a is stored in the lower triangle.
               if ((incx==1) .and. (incy==1)) then
                   do j = 1,n
                       if ((x(j)/=czero) .or. (y(j)/=czero)) then
                           temp1 = alpha*conjg(y(j))
                           temp2 = conjg(alpha*x(j))
                           a(j,j) = real(a(j,j),KIND=dp) +real(x(j)*temp1+y(j)*temp2,KIND=dp)
                                     
                           do i = j + 1,n
                               a(i,j) = a(i,j) + x(i)*temp1 + y(i)*temp2
                           end do
                       else
                           a(j,j) = real(a(j,j),KIND=dp)
                       end if
                   end do
               else
                   do j = 1,n
                       if ((x(jx)/=czero) .or. (y(jy)/=czero)) then
                           temp1 = alpha*conjg(y(jy))
                           temp2 = conjg(alpha*x(jx))
                           a(j,j) = real(a(j,j),KIND=dp) +real(x(jx)*temp1+y(jy)*temp2,KIND=dp)
                                     
                           ix = jx
                           iy = jy
                           do i = j + 1,n
                               ix = ix + incx
                               iy = iy + incy
                               a(i,j) = a(i,j) + x(ix)*temp1 + y(iy)*temp2
                           end do
                       else
                           a(j,j) = real(a(j,j),KIND=dp)
                       end if
                       jx = jx + incx
                       jy = jy + incy
                   end do
               end if
           end if
           return
     end subroutine stdlib_zher2


     pure subroutine stdlib_zher2k(uplo,trans,n,k,alpha,a,lda,b,ldb,beta,c,ldc)
     !! ZHER2K performs one of the hermitian rank 2k operations
     !! C := alpha*A*B**H + conjg( alpha )*B*A**H + beta*C,
     !! or
     !! C := alpha*A**H*B + conjg( alpha )*B**H*A + beta*C,
     !! where  alpha and beta  are scalars with  beta  real,  C is an  n by n
     !! hermitian matrix and  A and B  are  n by k matrices in the first case
     !! and  k by n  matrices in the second case.
        ! -- reference blas level3 routine --
        ! -- reference blas is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           ! Scalar Arguments 
           complex(dp), intent(in) :: alpha
           real(dp), intent(in) :: beta
           integer(ilp), intent(in) :: k, lda, ldb, ldc, n
           character, intent(in) :: trans, uplo
           ! Array Arguments 
           complex(dp), intent(in) :: a(lda,*), b(ldb,*)
           complex(dp), intent(inout) :: c(ldc,*)
        ! =====================================================================
           ! Intrinsic Functions 
           intrinsic :: real,conjg,max
           ! Local Scalars 
           complex(dp) :: temp1, temp2
           integer(ilp) :: i, info, j, l, nrowa
           logical(lk) :: upper
           
           
           ! test the input parameters.
           if (stdlib_lsame(trans,'N')) then
               nrowa = n
           else
               nrowa = k
           end if
           upper = stdlib_lsame(uplo,'U')
           info = 0
           if ((.not.upper) .and. (.not.stdlib_lsame(uplo,'L'))) then
               info = 1
           else if ((.not.stdlib_lsame(trans,'N')) .and.(.not.stdlib_lsame(trans,'C'))) &
                     then
               info = 2
           else if (n<0) then
               info = 3
           else if (k<0) then
               info = 4
           else if (lda<max(1,nrowa)) then
               info = 7
           else if (ldb<max(1,nrowa)) then
               info = 9
           else if (ldc<max(1,n)) then
               info = 12
           end if
           if (info/=0) then
               call stdlib_xerbla('ZHER2K',info)
               return
           end if
           ! quick return if possible.
           if ((n==0) .or. (((alpha==czero).or.(k==0)).and. (beta==one))) return
           ! and when  alpha.eq.czero.
           if (alpha==czero) then
               if (upper) then
                   if (beta==real(czero,KIND=dp)) then
                       do j = 1,n
                           do i = 1,j
                               c(i,j) = czero
                           end do
                       end do
                   else
                       do j = 1,n
                           do i = 1,j - 1
                               c(i,j) = beta*c(i,j)
                           end do
                           c(j,j) = beta*real(c(j,j),KIND=dp)
                       end do
                   end if
               else
                   if (beta==real(czero,KIND=dp)) then
                       do j = 1,n
                           do i = j,n
                               c(i,j) = czero
                           end do
                       end do
                   else
                       do j = 1,n
                           c(j,j) = beta*real(c(j,j),KIND=dp)
                           do i = j + 1,n
                               c(i,j) = beta*c(i,j)
                           end do
                       end do
                   end if
               end if
               return
           end if
           ! start the operations.
           if (stdlib_lsame(trans,'N')) then
              ! form  c := alpha*a*b**h + conjg( alpha )*b*a**h +
                         ! c.
               if (upper) then
                   do j = 1,n
                       if (beta==real(czero,KIND=dp)) then
                           do i = 1,j
                               c(i,j) = czero
                           end do
                       else if (beta/=one) then
                           do i = 1,j - 1
                               c(i,j) = beta*c(i,j)
                           end do
                           c(j,j) = beta*real(c(j,j),KIND=dp)
                       else
                           c(j,j) = real(c(j,j),KIND=dp)
                       end if
                       do l = 1,k
                           if ((a(j,l)/=czero) .or. (b(j,l)/=czero)) then
                               temp1 = alpha*conjg(b(j,l))
                               temp2 = conjg(alpha*a(j,l))
                               do i = 1,j - 1
                                   c(i,j) = c(i,j) + a(i,l)*temp1 +b(i,l)*temp2
                               end do
                               c(j,j) = real(c(j,j),KIND=dp) +real(a(j,l)*temp1+b(j,l)*temp2,&
                                         KIND=dp)
                           end if
                       end do
                   end do
               else
                   do j = 1,n
                       if (beta==real(czero,KIND=dp)) then
                           do i = j,n
                               c(i,j) = czero
                           end do
                       else if (beta/=one) then
                           do i = j + 1,n
                               c(i,j) = beta*c(i,j)
                           end do
                           c(j,j) = beta*real(c(j,j),KIND=dp)
                       else
                           c(j,j) = real(c(j,j),KIND=dp)
                       end if
                       do l = 1,k
                           if ((a(j,l)/=czero) .or. (b(j,l)/=czero)) then
                               temp1 = alpha*conjg(b(j,l))
                               temp2 = conjg(alpha*a(j,l))
                               do i = j + 1,n
                                   c(i,j) = c(i,j) + a(i,l)*temp1 +b(i,l)*temp2
                               end do
                               c(j,j) = real(c(j,j),KIND=dp) +real(a(j,l)*temp1+b(j,l)*temp2,&
                                         KIND=dp)
                           end if
                       end do
                   end do
               end if
           else
              ! form  c := alpha*a**h*b + conjg( alpha )*b**h*a +
                         ! c.
               if (upper) then
                   do j = 1,n
                       do i = 1,j
                           temp1 = czero
                           temp2 = czero
                           do l = 1,k
                               temp1 = temp1 + conjg(a(l,i))*b(l,j)
                               temp2 = temp2 + conjg(b(l,i))*a(l,j)
                           end do
                           if (i==j) then
                               if (beta==real(czero,KIND=dp)) then
                                   c(j,j) = real(alpha*temp1+conjg(alpha)*temp2,KIND=dp)
                               else
                                   c(j,j) = beta*real(c(j,j),KIND=dp) +real(alpha*temp1+conjg(&
                                             alpha)*temp2,KIND=dp)
                               end if
                           else
                               if (beta==real(czero,KIND=dp)) then
                                   c(i,j) = alpha*temp1 + conjg(alpha)*temp2
                               else
                                   c(i,j) = beta*c(i,j) + alpha*temp1 +conjg(alpha)*temp2
                               end if
                           end if
                       end do
                   end do
               else
                   do j = 1,n
                       do i = j,n
                           temp1 = czero
                           temp2 = czero
                           do l = 1,k
                               temp1 = temp1 + conjg(a(l,i))*b(l,j)
                               temp2 = temp2 + conjg(b(l,i))*a(l,j)
                           end do
                           if (i==j) then
                               if (beta==real(czero,KIND=dp)) then
                                   c(j,j) = real(alpha*temp1+conjg(alpha)*temp2,KIND=dp)
                               else
                                   c(j,j) = beta*real(c(j,j),KIND=dp) +real(alpha*temp1+conjg(&
                                             alpha)*temp2,KIND=dp)
                               end if
                           else
                               if (beta==real(czero,KIND=dp)) then
                                   c(i,j) = alpha*temp1 + conjg(alpha)*temp2
                               else
                                   c(i,j) = beta*c(i,j) + alpha*temp1 +conjg(alpha)*temp2
                               end if
                           end if
                       end do
                   end do
               end if
           end if
           return
     end subroutine stdlib_zher2k


     pure subroutine stdlib_zherk(uplo,trans,n,k,alpha,a,lda,beta,c,ldc)
     !! ZHERK performs one of the hermitian rank k operations
     !! C := alpha*A*A**H + beta*C,
     !! or
     !! C := alpha*A**H*A + beta*C,
     !! where  alpha and beta  are  real scalars,  C is an  n by n  hermitian
     !! matrix and  A  is an  n by k  matrix in the  first case and a  k by n
     !! matrix in the second case.
        ! -- reference blas level3 routine --
        ! -- reference blas is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           ! Scalar Arguments 
           real(dp), intent(in) :: alpha, beta
           integer(ilp), intent(in) :: k, lda, ldc, n
           character, intent(in) :: trans, uplo
           ! Array Arguments 
           complex(dp), intent(in) :: a(lda,*)
           complex(dp), intent(inout) :: c(ldc,*)
        ! =====================================================================
           ! Intrinsic Functions 
           intrinsic :: real,cmplx,conjg,max
           ! Local Scalars 
           complex(dp) :: temp
           real(dp) :: rtemp
           integer(ilp) :: i, info, j, l, nrowa
           logical(lk) :: upper
           
           ! test the input parameters.
           if (stdlib_lsame(trans,'N')) then
               nrowa = n
           else
               nrowa = k
           end if
           upper = stdlib_lsame(uplo,'U')
           info = 0
           if ((.not.upper) .and. (.not.stdlib_lsame(uplo,'L'))) then
               info = 1
           else if ((.not.stdlib_lsame(trans,'N')) .and.(.not.stdlib_lsame(trans,'C'))) &
                     then
               info = 2
           else if (n<0) then
               info = 3
           else if (k<0) then
               info = 4
           else if (lda<max(1,nrowa)) then
               info = 7
           else if (ldc<max(1,n)) then
               info = 10
           end if
           if (info/=0) then
               call stdlib_xerbla('ZHERK ',info)
               return
           end if
           ! quick return if possible.
           if ((n==0) .or. (((alpha==zero).or.(k==0)).and. (beta==one))) return
           ! and when  alpha.eq.zero.
           if (alpha==zero) then
               if (upper) then
                   if (beta==zero) then
                       do j = 1,n
                           do i = 1,j
                               c(i,j) = zero
                           end do
                       end do
                   else
                       do j = 1,n
                           do i = 1,j - 1
                               c(i,j) = beta*c(i,j)
                           end do
                           c(j,j) = beta*real(c(j,j),KIND=dp)
                       end do
                   end if
               else
                   if (beta==zero) then
                       do j = 1,n
                           do i = j,n
                               c(i,j) = zero
                           end do
                       end do
                   else
                       do j = 1,n
                           c(j,j) = beta*real(c(j,j),KIND=dp)
                           do i = j + 1,n
                               c(i,j) = beta*c(i,j)
                           end do
                       end do
                   end if
               end if
               return
           end if
           ! start the operations.
           if (stdlib_lsame(trans,'N')) then
              ! form  c := alpha*a*a**h + beta*c.
               if (upper) then
                   do j = 1,n
                       if (beta==zero) then
                           do i = 1,j
                               c(i,j) = zero
                           end do
                       else if (beta/=one) then
                           do i = 1,j - 1
                               c(i,j) = beta*c(i,j)
                           end do
                           c(j,j) = beta*real(c(j,j),KIND=dp)
                       else
                           c(j,j) = real(c(j,j),KIND=dp)
                       end if
                       do l = 1,k
                           if (a(j,l)/=cmplx(zero,KIND=dp)) then
                               temp = alpha*conjg(a(j,l))
                               do i = 1,j - 1
                                   c(i,j) = c(i,j) + temp*a(i,l)
                               end do
                               c(j,j) = real(c(j,j),KIND=dp) + real(temp*a(i,l),KIND=dp)
                           end if
                       end do
                   end do
               else
                   do j = 1,n
                       if (beta==zero) then
                           do i = j,n
                               c(i,j) = zero
                           end do
                       else if (beta/=one) then
                           c(j,j) = beta*real(c(j,j),KIND=dp)
                           do i = j + 1,n
                               c(i,j) = beta*c(i,j)
                           end do
                       else
                           c(j,j) = real(c(j,j),KIND=dp)
                       end if
                       do l = 1,k
                           if (a(j,l)/=cmplx(zero,KIND=dp)) then
                               temp = alpha*conjg(a(j,l))
                               c(j,j) = real(c(j,j),KIND=dp) + real(temp*a(j,l),KIND=dp)
                               do i = j + 1,n
                                   c(i,j) = c(i,j) + temp*a(i,l)
                               end do
                           end if
                       end do
                   end do
               end if
           else
              ! form  c := alpha*a**h*a + beta*c.
               if (upper) then
                   do j = 1,n
                       do i = 1,j - 1
                           temp = zero
                           do l = 1,k
                               temp = temp + conjg(a(l,i))*a(l,j)
                           end do
                           if (beta==zero) then
                               c(i,j) = alpha*temp
                           else
                               c(i,j) = alpha*temp + beta*c(i,j)
                           end if
                       end do
                       rtemp = zero
                       do l = 1,k
                           rtemp = rtemp + conjg(a(l,j))*a(l,j)
                       end do
                       if (beta==zero) then
                           c(j,j) = alpha*rtemp
                       else
                           c(j,j) = alpha*rtemp + beta*real(c(j,j),KIND=dp)
                       end if
                   end do
               else
                   do j = 1,n
                       rtemp = zero
                       do l = 1,k
                           rtemp = rtemp + conjg(a(l,j))*a(l,j)
                       end do
                       if (beta==zero) then
                           c(j,j) = alpha*rtemp
                       else
                           c(j,j) = alpha*rtemp + beta*real(c(j,j),KIND=dp)
                       end if
                       do i = j + 1,n
                           temp = zero
                           do l = 1,k
                               temp = temp + conjg(a(l,i))*a(l,j)
                           end do
                           if (beta==zero) then
                               c(i,j) = alpha*temp
                           else
                               c(i,j) = alpha*temp + beta*c(i,j)
                           end if
                       end do
                   end do
               end if
           end if
           return
     end subroutine stdlib_zherk


     pure subroutine stdlib_zhpmv(uplo,n,alpha,ap,x,incx,beta,y,incy)
     !! ZHPMV performs the matrix-vector operation
     !! y := alpha*A*x + beta*y,
     !! where alpha and beta are scalars, x and y are n element vectors and
     !! A is an n by n hermitian matrix, supplied in packed form.
        ! -- reference blas level2 routine --
        ! -- reference blas is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           ! Scalar Arguments 
           complex(dp), intent(in) :: alpha, beta
           integer(ilp), intent(in) :: incx, incy, n
           character, intent(in) :: uplo
           ! Array Arguments 
           complex(dp), intent(in) :: ap(*), x(*)
           complex(dp), intent(inout) :: y(*)
        ! =====================================================================
           
           
           ! Local Scalars 
           complex(dp) :: temp1, temp2
           integer(ilp) :: i, info, ix, iy, j, jx, jy, k, kk, kx, ky
           ! Intrinsic Functions 
           intrinsic :: real,conjg
           ! test the input parameters.
           info = 0
           if (.not.stdlib_lsame(uplo,'U') .and. .not.stdlib_lsame(uplo,'L')) then
               info = 1
           else if (n<0) then
               info = 2
           else if (incx==0) then
               info = 6
           else if (incy==0) then
               info = 9
           end if
           if (info/=0) then
               call stdlib_xerbla('ZHPMV ',info)
               return
           end if
           ! quick return if possible.
           if ((n==0) .or. ((alpha==czero).and. (beta==cone))) return
           ! set up the start points in  x  and  y.
           if (incx>0) then
               kx = 1
           else
               kx = 1 - (n-1)*incx
           end if
           if (incy>0) then
               ky = 1
           else
               ky = 1 - (n-1)*incy
           end if
           ! start the operations. in this version the elements of the array ap
           ! are accessed sequentially with cone pass through ap.
           ! first form  y := beta*y.
           if (beta/=cone) then
               if (incy==1) then
                   if (beta==czero) then
                       do i = 1,n
                           y(i) = czero
                       end do
                   else
                       do i = 1,n
                           y(i) = beta*y(i)
                       end do
                   end if
               else
                   iy = ky
                   if (beta==czero) then
                       do i = 1,n
                           y(iy) = czero
                           iy = iy + incy
                       end do
                   else
                       do i = 1,n
                           y(iy) = beta*y(iy)
                           iy = iy + incy
                       end do
                   end if
               end if
           end if
           if (alpha==czero) return
           kk = 1
           if (stdlib_lsame(uplo,'U')) then
              ! form  y  when ap contains the upper triangle.
               if ((incx==1) .and. (incy==1)) then
                   do j = 1,n
                       temp1 = alpha*x(j)
                       temp2 = czero
                       k = kk
                       do i = 1,j - 1
                           y(i) = y(i) + temp1*ap(k)
                           temp2 = temp2 + conjg(ap(k))*x(i)
                           k = k + 1
                       end do
                       y(j) = y(j) + temp1*real(ap(kk+j-1),KIND=dp) + alpha*temp2
                       kk = kk + j
                   end do
               else
                   jx = kx
                   jy = ky
                   do j = 1,n
                       temp1 = alpha*x(jx)
                       temp2 = czero
                       ix = kx
                       iy = ky
                       do k = kk,kk + j - 2
                           y(iy) = y(iy) + temp1*ap(k)
                           temp2 = temp2 + conjg(ap(k))*x(ix)
                           ix = ix + incx
                           iy = iy + incy
                       end do
                       y(jy) = y(jy) + temp1*real(ap(kk+j-1),KIND=dp) + alpha*temp2
                       jx = jx + incx
                       jy = jy + incy
                       kk = kk + j
                   end do
               end if
           else
              ! form  y  when ap contains the lower triangle.
               if ((incx==1) .and. (incy==1)) then
                   do j = 1,n
                       temp1 = alpha*x(j)
                       temp2 = czero
                       y(j) = y(j) + temp1*real(ap(kk),KIND=dp)
                       k = kk + 1
                       do i = j + 1,n
                           y(i) = y(i) + temp1*ap(k)
                           temp2 = temp2 + conjg(ap(k))*x(i)
                           k = k + 1
                       end do
                       y(j) = y(j) + alpha*temp2
                       kk = kk + (n-j+1)
                   end do
               else
                   jx = kx
                   jy = ky
                   do j = 1,n
                       temp1 = alpha*x(jx)
                       temp2 = czero
                       y(jy) = y(jy) + temp1*real(ap(kk),KIND=dp)
                       ix = jx
                       iy = jy
                       do k = kk + 1,kk + n - j
                           ix = ix + incx
                           iy = iy + incy
                           y(iy) = y(iy) + temp1*ap(k)
                           temp2 = temp2 + conjg(ap(k))*x(ix)
                       end do
                       y(jy) = y(jy) + alpha*temp2
                       jx = jx + incx
                       jy = jy + incy
                       kk = kk + (n-j+1)
                   end do
               end if
           end if
           return
     end subroutine stdlib_zhpmv


     pure subroutine stdlib_zhpr(uplo,n,alpha,x,incx,ap)
     !! ZHPR performs the hermitian rank 1 operation
     !! A := alpha*x*x**H + A,
     !! where alpha is a real scalar, x is an n element vector and A is an
     !! n by n hermitian matrix, supplied in packed form.
        ! -- reference blas level2 routine --
        ! -- reference blas is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           ! Scalar Arguments 
           real(dp), intent(in) :: alpha
           integer(ilp), intent(in) :: incx, n
           character, intent(in) :: uplo
           ! Array Arguments 
           complex(dp), intent(inout) :: ap(*)
           complex(dp), intent(in) :: x(*)
        ! =====================================================================
           
           ! Local Scalars 
           complex(dp) :: temp
           integer(ilp) :: i, info, ix, j, jx, k, kk, kx
           ! Intrinsic Functions 
           intrinsic :: real,conjg
           ! test the input parameters.
           info = 0
           if (.not.stdlib_lsame(uplo,'U') .and. .not.stdlib_lsame(uplo,'L')) then
               info = 1
           else if (n<0) then
               info = 2
           else if (incx==0) then
               info = 5
           end if
           if (info/=0) then
               call stdlib_xerbla('ZHPR  ',info)
               return
           end if
           ! quick return if possible.
           if ((n==0) .or. (alpha==real(czero,KIND=dp))) return
           ! set the start point in x if the increment is not unity.
           if (incx<=0) then
               kx = 1 - (n-1)*incx
           else if (incx/=1) then
               kx = 1
           end if
           ! start the operations. in this version the elements of the array ap
           ! are accessed sequentially with cone pass through ap.
           kk = 1
           if (stdlib_lsame(uplo,'U')) then
              ! form  a  when upper triangle is stored in ap.
               if (incx==1) then
                   do j = 1,n
                       if (x(j)/=czero) then
                           temp = alpha*conjg(x(j))
                           k = kk
                           do i = 1,j - 1
                               ap(k) = ap(k) + x(i)*temp
                               k = k + 1
                           end do
                           ap(kk+j-1) = real(ap(kk+j-1),KIND=dp) + real(x(j)*temp,KIND=dp)
                       else
                           ap(kk+j-1) = real(ap(kk+j-1),KIND=dp)
                       end if
                       kk = kk + j
                   end do
               else
                   jx = kx
                   do j = 1,n
                       if (x(jx)/=czero) then
                           temp = alpha*conjg(x(jx))
                           ix = kx
                           do k = kk,kk + j - 2
                               ap(k) = ap(k) + x(ix)*temp
                               ix = ix + incx
                           end do
                           ap(kk+j-1) = real(ap(kk+j-1),KIND=dp) + real(x(jx)*temp,KIND=dp)
                                     
                       else
                           ap(kk+j-1) = real(ap(kk+j-1),KIND=dp)
                       end if
                       jx = jx + incx
                       kk = kk + j
                   end do
               end if
           else
              ! form  a  when lower triangle is stored in ap.
               if (incx==1) then
                   do j = 1,n
                       if (x(j)/=czero) then
                           temp = alpha*conjg(x(j))
                           ap(kk) = real(ap(kk),KIND=dp) + real(temp*x(j),KIND=dp)
                           k = kk + 1
                           do i = j + 1,n
                               ap(k) = ap(k) + x(i)*temp
                               k = k + 1
                           end do
                       else
                           ap(kk) = real(ap(kk),KIND=dp)
                       end if
                       kk = kk + n - j + 1
                   end do
               else
                   jx = kx
                   do j = 1,n
                       if (x(jx)/=czero) then
                           temp = alpha*conjg(x(jx))
                           ap(kk) = real(ap(kk),KIND=dp) + real(temp*x(jx),KIND=dp)
                           ix = jx
                           do k = kk + 1,kk + n - j
                               ix = ix + incx
                               ap(k) = ap(k) + x(ix)*temp
                           end do
                       else
                           ap(kk) = real(ap(kk),KIND=dp)
                       end if
                       jx = jx + incx
                       kk = kk + n - j + 1
                   end do
               end if
           end if
           return
     end subroutine stdlib_zhpr


     pure subroutine stdlib_zhpr2(uplo,n,alpha,x,incx,y,incy,ap)
     !! ZHPR2 performs the hermitian rank 2 operation
     !! A := alpha*x*y**H + conjg( alpha )*y*x**H + A,
     !! where alpha is a scalar, x and y are n element vectors and A is an
     !! n by n hermitian matrix, supplied in packed form.
        ! -- reference blas level2 routine --
        ! -- reference blas is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           ! Scalar Arguments 
           complex(dp), intent(in) :: alpha
           integer(ilp), intent(in) :: incx, incy, n
           character, intent(in) :: uplo
           ! Array Arguments 
           complex(dp), intent(inout) :: ap(*)
           complex(dp), intent(in) :: x(*), y(*)
        ! =====================================================================
           
           ! Local Scalars 
           complex(dp) :: temp1, temp2
           integer(ilp) :: i, info, ix, iy, j, jx, jy, k, kk, kx, ky
           ! Intrinsic Functions 
           intrinsic :: real,conjg
           ! test the input parameters.
           info = 0
           if (.not.stdlib_lsame(uplo,'U') .and. .not.stdlib_lsame(uplo,'L')) then
               info = 1
           else if (n<0) then
               info = 2
           else if (incx==0) then
               info = 5
           else if (incy==0) then
               info = 7
           end if
           if (info/=0) then
               call stdlib_xerbla('ZHPR2 ',info)
               return
           end if
           ! quick return if possible.
           if ((n==0) .or. (alpha==czero)) return
           ! set up the start points in x and y if the increments are not both
           ! unity.
           if ((incx/=1) .or. (incy/=1)) then
               if (incx>0) then
                   kx = 1
               else
                   kx = 1 - (n-1)*incx
               end if
               if (incy>0) then
                   ky = 1
               else
                   ky = 1 - (n-1)*incy
               end if
               jx = kx
               jy = ky
           end if
           ! start the operations. in this version the elements of the array ap
           ! are accessed sequentially with cone pass through ap.
           kk = 1
           if (stdlib_lsame(uplo,'U')) then
              ! form  a  when upper triangle is stored in ap.
               if ((incx==1) .and. (incy==1)) then
                   do j = 1,n
                       if ((x(j)/=czero) .or. (y(j)/=czero)) then
                           temp1 = alpha*conjg(y(j))
                           temp2 = conjg(alpha*x(j))
                           k = kk
                           do i = 1,j - 1
                               ap(k) = ap(k) + x(i)*temp1 + y(i)*temp2
                               k = k + 1
                           end do
                           ap(kk+j-1) = real(ap(kk+j-1),KIND=dp) +real(x(j)*temp1+y(j)*temp2,&
                                     KIND=dp)
                       else
                           ap(kk+j-1) = real(ap(kk+j-1),KIND=dp)
                       end if
                       kk = kk + j
                   end do
               else
                   do j = 1,n
                       if ((x(jx)/=czero) .or. (y(jy)/=czero)) then
                           temp1 = alpha*conjg(y(jy))
                           temp2 = conjg(alpha*x(jx))
                           ix = kx
                           iy = ky
                           do k = kk,kk + j - 2
                               ap(k) = ap(k) + x(ix)*temp1 + y(iy)*temp2
                               ix = ix + incx
                               iy = iy + incy
                           end do
                           ap(kk+j-1) = real(ap(kk+j-1),KIND=dp) +real(x(jx)*temp1+y(jy)*temp2,&
                                     KIND=dp)
                       else
                           ap(kk+j-1) = real(ap(kk+j-1),KIND=dp)
                       end if
                       jx = jx + incx
                       jy = jy + incy
                       kk = kk + j
                   end do
               end if
           else
              ! form  a  when lower triangle is stored in ap.
               if ((incx==1) .and. (incy==1)) then
                   do j = 1,n
                       if ((x(j)/=czero) .or. (y(j)/=czero)) then
                           temp1 = alpha*conjg(y(j))
                           temp2 = conjg(alpha*x(j))
                           ap(kk) = real(ap(kk),KIND=dp) +real(x(j)*temp1+y(j)*temp2,KIND=dp)
                                     
                           k = kk + 1
                           do i = j + 1,n
                               ap(k) = ap(k) + x(i)*temp1 + y(i)*temp2
                               k = k + 1
                           end do
                       else
                           ap(kk) = real(ap(kk),KIND=dp)
                       end if
                       kk = kk + n - j + 1
                   end do
               else
                   do j = 1,n
                       if ((x(jx)/=czero) .or. (y(jy)/=czero)) then
                           temp1 = alpha*conjg(y(jy))
                           temp2 = conjg(alpha*x(jx))
                           ap(kk) = real(ap(kk),KIND=dp) +real(x(jx)*temp1+y(jy)*temp2,KIND=dp)
                                     
                           ix = jx
                           iy = jy
                           do k = kk + 1,kk + n - j
                               ix = ix + incx
                               iy = iy + incy
                               ap(k) = ap(k) + x(ix)*temp1 + y(iy)*temp2
                           end do
                       else
                           ap(kk) = real(ap(kk),KIND=dp)
                       end if
                       jx = jx + incx
                       jy = jy + incy
                       kk = kk + n - j + 1
                   end do
               end if
           end if
           return
     end subroutine stdlib_zhpr2


     pure subroutine stdlib_zrotg( a, b, c, s )
     !! The computation uses the formulas
     !! |x| = sqrt( Re(x)**2 + Im(x)**2 )
     !! sgn(x) = x / |x|  if x /= 0
     !! = 1        if x  = 0
     !! c = |a| / sqrt(|a|**2 + |b|**2)
     !! s = sgn(a) * conjg(b) / sqrt(|a|**2 + |b|**2)
     !! When a and b are real and r /= 0, the formulas simplify to
     !! r = sgn(a)*sqrt(|a|**2 + |b|**2)
     !! c = a / r
     !! s = b / r
     !! the same as in DROTG when |a| > |b|.  When |b| >= |a|, the
     !! sign of c and s will be different from those computed by DROTG
     !! if the signs of a and b are not the same.
        ! -- reference blas level1 routine --
        ! -- reference blas is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
        ! Constants 
        integer, parameter :: wp = kind(1._dp)
        ! Scaling Constants 
        ! Scalar Arguments 
        real(dp), intent(out) :: c
        complex(dp), intent(inout) :: a
        complex(dp), intent(in) :: b
        complex(dp), intent(out) :: s
        ! Local Scalars 
        real(dp) :: d, f1, f2, g1, g2, h2, p, u, uu, v, vv, w
        complex(dp) :: f, fs, g, gs, r, t
        ! Intrinsic Functions 
        intrinsic :: abs,aimag,conjg,max,min,real,sqrt
        ! Statement Functions 
        real(dp) :: abssq
        ! Statement Function Definitions 
        abssq( t ) = real( t,KIND=dp)**2 + aimag( t )**2
        ! Executable Statements 
        f = a
        g = b
        if( g == czero ) then
           c = one
           s = czero
           r = f
        else if( f == czero ) then
           c = zero
           g1 = max( abs(real(g,KIND=dp)), abs(aimag(g)) )
           if( g1 > rtmin .and. g1 < rtmax ) then
              ! use unscaled algorithm
              g2 = abssq( g )
              d = sqrt( g2 )
              s = conjg( g ) / d
              r = d
           else
              ! use scaled algorithm
              u = min( safmax, max( safmin, g1 ) )
              uu = one / u
              gs = g*uu
              g2 = abssq( gs )
              d = sqrt( g2 )
              s = conjg( gs ) / d
              r = d*u
           end if
        else
           f1 = max( abs(real(f,KIND=dp)), abs(aimag(f)) )
           g1 = max( abs(real(g,KIND=dp)), abs(aimag(g)) )
     if( f1 > rtmin .and. f1 < rtmax .and.          g1 > rtmin .and. g1 < rtmax ) then
              ! use unscaled algorithm
              f2 = abssq( f )
              g2 = abssq( g )
              h2 = f2 + g2
              if( f2 > rtmin .and. h2 < rtmax ) then
                 d = sqrt( f2*h2 )
              else
                 d = sqrt( f2 )*sqrt( h2 )
              end if
              p = 1 / d
              c = f2*p
              s = conjg( g )*( f*p )
              r = f*( h2*p )
           else
              ! use scaled algorithm
              u = min( safmax, max( safmin, f1, g1 ) )
              uu = one / u
              gs = g*uu
              g2 = abssq( gs )
              if( f1*uu < rtmin ) then
                 ! f is not well-scaled when scaled by g1.
                 ! use a different scaling for f.
                 v = min( safmax, max( safmin, f1 ) )
                 vv = one / v
                 w = v * uu
                 fs = f*vv
                 f2 = abssq( fs )
                 h2 = f2*w**2 + g2
              else
                 ! otherwise use the same scaling for f and g.
                 w = one
                 fs = f*uu
                 f2 = abssq( fs )
                 h2 = f2 + g2
              end if
              if( f2 > rtmin .and. h2 < rtmax ) then
                 d = sqrt( f2*h2 )
              else
                 d = sqrt( f2 )*sqrt( h2 )
              end if
              p = 1 / d
              c = ( f2*p )*w
              s = conjg( gs )*( fs*p )
              r = ( fs*( h2*p ) )*u
           end if
        end if
        a = r
        return
     end subroutine stdlib_zrotg


     pure subroutine stdlib_zscal(n,za,zx,incx)
     !! ZSCAL scales a vector by a constant.
        ! -- reference blas level1 routine --
        ! -- reference blas is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           ! Scalar Arguments 
           complex(dp), intent(in) :: za
           integer(ilp), intent(in) :: incx, n
           ! Array Arguments 
           complex(dp), intent(inout) :: zx(*)
        ! =====================================================================
           ! Local Scalars 
           integer(ilp) :: i, nincx
           if (n<=0 .or. incx<=0) return
           if (incx==1) then
              ! code for increment equal to 1
              do i = 1,n
                 zx(i) = za*zx(i)
              end do
           else
              ! code for increment not equal to 1
              nincx = n*incx
              do i = 1,nincx,incx
                 zx(i) = za*zx(i)
              end do
           end if
           return
     end subroutine stdlib_zscal


     pure subroutine stdlib_zswap(n,zx,incx,zy,incy)
     !! ZSWAP interchanges two vectors.
        ! -- reference blas level1 routine --
        ! -- reference blas is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           ! Scalar Arguments 
           integer(ilp), intent(in) :: incx, incy, n
           ! Array Arguments 
           complex(dp), intent(inout) :: zx(*), zy(*)
        ! =====================================================================
           ! Local Scalars 
           complex(dp) :: ztemp
           integer(ilp) :: i, ix, iy
           if (n<=0) return
           if (incx==1 .and. incy==1) then
             ! code for both increments equal to 1
              do i = 1,n
                 ztemp = zx(i)
                 zx(i) = zy(i)
                 zy(i) = ztemp
              end do
           else
             ! code for unequal increments or equal increments not equal
               ! to 1
              ix = 1
              iy = 1
              if (incx<0) ix = (-n+1)*incx + 1
              if (incy<0) iy = (-n+1)*incy + 1
              do i = 1,n
                 ztemp = zx(ix)
                 zx(ix) = zy(iy)
                 zy(iy) = ztemp
                 ix = ix + incx
                 iy = iy + incy
              end do
           end if
           return
     end subroutine stdlib_zswap


     pure subroutine stdlib_zsymm(side,uplo,m,n,alpha,a,lda,b,ldb,beta,c,ldc)
     !! ZSYMM performs one of the matrix-matrix operations
     !! C := alpha*A*B + beta*C,
     !! or
     !! C := alpha*B*A + beta*C,
     !! where  alpha and beta are scalars, A is a symmetric matrix and  B and
     !! C are m by n matrices.
        ! -- reference blas level3 routine --
        ! -- reference blas is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           ! Scalar Arguments 
           complex(dp), intent(in) :: alpha, beta
           integer(ilp), intent(in) :: lda, ldb, ldc, m, n
           character, intent(in) :: side, uplo
           ! Array Arguments 
           complex(dp), intent(in) :: a(lda,*), b(ldb,*)
           complex(dp), intent(inout) :: c(ldc,*)
        ! =====================================================================
           ! Intrinsic Functions 
           intrinsic :: max
           ! Local Scalars 
           complex(dp) :: temp1, temp2
           integer(ilp) :: i, info, j, k, nrowa
           logical(lk) :: upper
           
           
           ! set nrowa as the number of rows of a.
           if (stdlib_lsame(side,'L')) then
               nrowa = m
           else
               nrowa = n
           end if
           upper = stdlib_lsame(uplo,'U')
           ! test the input parameters.
           info = 0
           if ((.not.stdlib_lsame(side,'L')) .and. (.not.stdlib_lsame(side,'R'))) then
               info = 1
           else if ((.not.upper) .and. (.not.stdlib_lsame(uplo,'L'))) then
               info = 2
           else if (m<0) then
               info = 3
           else if (n<0) then
               info = 4
           else if (lda<max(1,nrowa)) then
               info = 7
           else if (ldb<max(1,m)) then
               info = 9
           else if (ldc<max(1,m)) then
               info = 12
           end if
           if (info/=0) then
               call stdlib_xerbla('ZSYMM ',info)
               return
           end if
           ! quick return if possible.
           if ((m==0) .or. (n==0) .or.((alpha==czero).and. (beta==cone))) return
           ! and when  alpha.eq.czero.
           if (alpha==czero) then
               if (beta==czero) then
                   do j = 1,n
                       do i = 1,m
                           c(i,j) = czero
                       end do
                   end do
               else
                   do j = 1,n
                       do i = 1,m
                           c(i,j) = beta*c(i,j)
                       end do
                   end do
               end if
               return
           end if
           ! start the operations.
           if (stdlib_lsame(side,'L')) then
              ! form  c := alpha*a*b + beta*c.
               if (upper) then
                   do j = 1,n
                       do i = 1,m
                           temp1 = alpha*b(i,j)
                           temp2 = czero
                           do k = 1,i - 1
                               c(k,j) = c(k,j) + temp1*a(k,i)
                               temp2 = temp2 + b(k,j)*a(k,i)
                           end do
                           if (beta==czero) then
                               c(i,j) = temp1*a(i,i) + alpha*temp2
                           else
                               c(i,j) = beta*c(i,j) + temp1*a(i,i) +alpha*temp2
                           end if
                       end do
                   end do
               else
                   do j = 1,n
                       do i = m,1,-1
                           temp1 = alpha*b(i,j)
                           temp2 = czero
                           do k = i + 1,m
                               c(k,j) = c(k,j) + temp1*a(k,i)
                               temp2 = temp2 + b(k,j)*a(k,i)
                           end do
                           if (beta==czero) then
                               c(i,j) = temp1*a(i,i) + alpha*temp2
                           else
                               c(i,j) = beta*c(i,j) + temp1*a(i,i) +alpha*temp2
                           end if
                       end do
                   end do
               end if
           else
              ! form  c := alpha*b*a + beta*c.
               loop_170: do j = 1,n
                   temp1 = alpha*a(j,j)
                   if (beta==czero) then
                       do i = 1,m
                           c(i,j) = temp1*b(i,j)
                       end do
                   else
                       do i = 1,m
                           c(i,j) = beta*c(i,j) + temp1*b(i,j)
                       end do
                   end if
                   do k = 1,j - 1
                       if (upper) then
                           temp1 = alpha*a(k,j)
                       else
                           temp1 = alpha*a(j,k)
                       end if
                       do i = 1,m
                           c(i,j) = c(i,j) + temp1*b(i,k)
                       end do
                   end do
                   do k = j + 1,n
                       if (upper) then
                           temp1 = alpha*a(j,k)
                       else
                           temp1 = alpha*a(k,j)
                       end if
                       do i = 1,m
                           c(i,j) = c(i,j) + temp1*b(i,k)
                       end do
                   end do
               end do loop_170
           end if
           return
     end subroutine stdlib_zsymm


     pure subroutine stdlib_zsyr2k(uplo,trans,n,k,alpha,a,lda,b,ldb,beta,c,ldc)
     !! ZSYR2K performs one of the symmetric rank 2k operations
     !! C := alpha*A*B**T + alpha*B*A**T + beta*C,
     !! or
     !! C := alpha*A**T*B + alpha*B**T*A + beta*C,
     !! where  alpha and beta  are scalars,  C is an  n by n symmetric matrix
     !! and  A and B  are  n by k  matrices  in the  first  case  and  k by n
     !! matrices in the second case.
        ! -- reference blas level3 routine --
        ! -- reference blas is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           ! Scalar Arguments 
           complex(dp), intent(in) :: alpha, beta
           integer(ilp), intent(in) :: k, lda, ldb, ldc, n
           character, intent(in) :: trans, uplo
           ! Array Arguments 
           complex(dp), intent(in) :: a(lda,*), b(ldb,*)
           complex(dp), intent(inout) :: c(ldc,*)
        ! =====================================================================
           ! Intrinsic Functions 
           intrinsic :: max
           ! Local Scalars 
           complex(dp) :: temp1, temp2
           integer(ilp) :: i, info, j, l, nrowa
           logical(lk) :: upper
           
           
           ! test the input parameters.
           if (stdlib_lsame(trans,'N')) then
               nrowa = n
           else
               nrowa = k
           end if
           upper = stdlib_lsame(uplo,'U')
           info = 0
           if ((.not.upper) .and. (.not.stdlib_lsame(uplo,'L'))) then
               info = 1
           else if ((.not.stdlib_lsame(trans,'N')) .and.(.not.stdlib_lsame(trans,'T'))) &
                     then
               info = 2
           else if (n<0) then
               info = 3
           else if (k<0) then
               info = 4
           else if (lda<max(1,nrowa)) then
               info = 7
           else if (ldb<max(1,nrowa)) then
               info = 9
           else if (ldc<max(1,n)) then
               info = 12
           end if
           if (info/=0) then
               call stdlib_xerbla('ZSYR2K',info)
               return
           end if
           ! quick return if possible.
           if ((n==0) .or. (((alpha==czero).or.(k==0)).and. (beta==cone))) return
           ! and when  alpha.eq.czero.
           if (alpha==czero) then
               if (upper) then
                   if (beta==czero) then
                       do j = 1,n
                           do i = 1,j
                               c(i,j) = czero
                           end do
                       end do
                   else
                       do j = 1,n
                           do i = 1,j
                               c(i,j) = beta*c(i,j)
                           end do
                       end do
                   end if
               else
                   if (beta==czero) then
                       do j = 1,n
                           do i = j,n
                               c(i,j) = czero
                           end do
                       end do
                   else
                       do j = 1,n
                           do i = j,n
                               c(i,j) = beta*c(i,j)
                           end do
                       end do
                   end if
               end if
               return
           end if
           ! start the operations.
           if (stdlib_lsame(trans,'N')) then
              ! form  c := alpha*a*b**t + alpha*b*a**t + c.
               if (upper) then
                   do j = 1,n
                       if (beta==czero) then
                           do i = 1,j
                               c(i,j) = czero
                           end do
                       else if (beta/=cone) then
                           do i = 1,j
                               c(i,j) = beta*c(i,j)
                           end do
                       end if
                       do l = 1,k
                           if ((a(j,l)/=czero) .or. (b(j,l)/=czero)) then
                               temp1 = alpha*b(j,l)
                               temp2 = alpha*a(j,l)
                               do i = 1,j
                                   c(i,j) = c(i,j) + a(i,l)*temp1 +b(i,l)*temp2
                               end do
                           end if
                       end do
                   end do
               else
                   do j = 1,n
                       if (beta==czero) then
                           do i = j,n
                               c(i,j) = czero
                           end do
                       else if (beta/=cone) then
                           do i = j,n
                               c(i,j) = beta*c(i,j)
                           end do
                       end if
                       do l = 1,k
                           if ((a(j,l)/=czero) .or. (b(j,l)/=czero)) then
                               temp1 = alpha*b(j,l)
                               temp2 = alpha*a(j,l)
                               do i = j,n
                                   c(i,j) = c(i,j) + a(i,l)*temp1 +b(i,l)*temp2
                               end do
                           end if
                       end do
                   end do
               end if
           else
              ! form  c := alpha*a**t*b + alpha*b**t*a + c.
               if (upper) then
                   do j = 1,n
                       do i = 1,j
                           temp1 = czero
                           temp2 = czero
                           do l = 1,k
                               temp1 = temp1 + a(l,i)*b(l,j)
                               temp2 = temp2 + b(l,i)*a(l,j)
                           end do
                           if (beta==czero) then
                               c(i,j) = alpha*temp1 + alpha*temp2
                           else
                               c(i,j) = beta*c(i,j) + alpha*temp1 +alpha*temp2
                           end if
                       end do
                   end do
               else
                   do j = 1,n
                       do i = j,n
                           temp1 = czero
                           temp2 = czero
                           do l = 1,k
                               temp1 = temp1 + a(l,i)*b(l,j)
                               temp2 = temp2 + b(l,i)*a(l,j)
                           end do
                           if (beta==czero) then
                               c(i,j) = alpha*temp1 + alpha*temp2
                           else
                               c(i,j) = beta*c(i,j) + alpha*temp1 +alpha*temp2
                           end if
                       end do
                   end do
               end if
           end if
           return
     end subroutine stdlib_zsyr2k


     pure subroutine stdlib_zsyrk(uplo,trans,n,k,alpha,a,lda,beta,c,ldc)
     !! ZSYRK performs one of the symmetric rank k operations
     !! C := alpha*A*A**T + beta*C,
     !! or
     !! C := alpha*A**T*A + beta*C,
     !! where  alpha and beta  are scalars,  C is an  n by n symmetric matrix
     !! and  A  is an  n by k  matrix in the first case and a  k by n  matrix
     !! in the second case.
        ! -- reference blas level3 routine --
        ! -- reference blas is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           ! Scalar Arguments 
           complex(dp), intent(in) :: alpha, beta
           integer(ilp), intent(in) :: k, lda, ldc, n
           character, intent(in) :: trans, uplo
           ! Array Arguments 
           complex(dp), intent(in) :: a(lda,*)
           complex(dp), intent(inout) :: c(ldc,*)
        ! =====================================================================
           ! Intrinsic Functions 
           intrinsic :: max
           ! Local Scalars 
           complex(dp) :: temp
           integer(ilp) :: i, info, j, l, nrowa
           logical(lk) :: upper
           
           
           ! test the input parameters.
           if (stdlib_lsame(trans,'N')) then
               nrowa = n
           else
               nrowa = k
           end if
           upper = stdlib_lsame(uplo,'U')
           info = 0
           if ((.not.upper) .and. (.not.stdlib_lsame(uplo,'L'))) then
               info = 1
           else if ((.not.stdlib_lsame(trans,'N')) .and.(.not.stdlib_lsame(trans,'T'))) &
                     then
               info = 2
           else if (n<0) then
               info = 3
           else if (k<0) then
               info = 4
           else if (lda<max(1,nrowa)) then
               info = 7
           else if (ldc<max(1,n)) then
               info = 10
           end if
           if (info/=0) then
               call stdlib_xerbla('ZSYRK ',info)
               return
           end if
           ! quick return if possible.
           if ((n==0) .or. (((alpha==czero).or.(k==0)).and. (beta==cone))) return
           ! and when  alpha.eq.czero.
           if (alpha==czero) then
               if (upper) then
                   if (beta==czero) then
                       do j = 1,n
                           do i = 1,j
                               c(i,j) = czero
                           end do
                       end do
                   else
                       do j = 1,n
                           do i = 1,j
                               c(i,j) = beta*c(i,j)
                           end do
                       end do
                   end if
               else
                   if (beta==czero) then
                       do j = 1,n
                           do i = j,n
                               c(i,j) = czero
                           end do
                       end do
                   else
                       do j = 1,n
                           do i = j,n
                               c(i,j) = beta*c(i,j)
                           end do
                       end do
                   end if
               end if
               return
           end if
           ! start the operations.
           if (stdlib_lsame(trans,'N')) then
              ! form  c := alpha*a*a**t + beta*c.
               if (upper) then
                   do j = 1,n
                       if (beta==czero) then
                           do i = 1,j
                               c(i,j) = czero
                           end do
                       else if (beta/=cone) then
                           do i = 1,j
                               c(i,j) = beta*c(i,j)
                           end do
                       end if
                       do l = 1,k
                           if (a(j,l)/=czero) then
                               temp = alpha*a(j,l)
                               do i = 1,j
                                   c(i,j) = c(i,j) + temp*a(i,l)
                               end do
                           end if
                       end do
                   end do
               else
                   do j = 1,n
                       if (beta==czero) then
                           do i = j,n
                               c(i,j) = czero
                           end do
                       else if (beta/=cone) then
                           do i = j,n
                               c(i,j) = beta*c(i,j)
                           end do
                       end if
                       do l = 1,k
                           if (a(j,l)/=czero) then
                               temp = alpha*a(j,l)
                               do i = j,n
                                   c(i,j) = c(i,j) + temp*a(i,l)
                               end do
                           end if
                       end do
                   end do
               end if
           else
              ! form  c := alpha*a**t*a + beta*c.
               if (upper) then
                   do j = 1,n
                       do i = 1,j
                           temp = czero
                           do l = 1,k
                               temp = temp + a(l,i)*a(l,j)
                           end do
                           if (beta==czero) then
                               c(i,j) = alpha*temp
                           else
                               c(i,j) = alpha*temp + beta*c(i,j)
                           end if
                       end do
                   end do
               else
                   do j = 1,n
                       do i = j,n
                           temp = czero
                           do l = 1,k
                               temp = temp + a(l,i)*a(l,j)
                           end do
                           if (beta==czero) then
                               c(i,j) = alpha*temp
                           else
                               c(i,j) = alpha*temp + beta*c(i,j)
                           end if
                       end do
                   end do
               end if
           end if
           return
     end subroutine stdlib_zsyrk


     pure subroutine stdlib_ztbmv(uplo,trans,diag,n,k,a,lda,x,incx)
     !! ZTBMV performs one of the matrix-vector operations
     !! x := A*x,   or   x := A**T*x,   or   x := A**H*x,
     !! where x is an n element vector and  A is an n by n unit, or non-unit,
     !! upper or lower triangular band matrix, with ( k + 1 ) diagonals.
        ! -- reference blas level2 routine --
        ! -- reference blas is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           ! Scalar Arguments 
           integer(ilp), intent(in) :: incx, k, lda, n
           character, intent(in) :: diag, trans, uplo
           ! Array Arguments 
           complex(dp), intent(in) :: a(lda,*)
           complex(dp), intent(inout) :: x(*)
        ! =====================================================================
           
           ! Local Scalars 
           complex(dp) :: temp
           integer(ilp) :: i, info, ix, j, jx, kplus1, kx, l
           logical(lk) :: noconj, nounit
           ! Intrinsic Functions 
           intrinsic :: conjg,max,min
           ! test the input parameters.
           info = 0
           if (.not.stdlib_lsame(uplo,'U') .and. .not.stdlib_lsame(uplo,'L')) then
               info = 1
           else if (.not.stdlib_lsame(trans,'N') .and. .not.stdlib_lsame(trans,'T') &
                     .and..not.stdlib_lsame(trans,'C')) then
               info = 2
           else if (.not.stdlib_lsame(diag,'U') .and. .not.stdlib_lsame(diag,'N')) then
               info = 3
           else if (n<0) then
               info = 4
           else if (k<0) then
               info = 5
           else if (lda< (k+1)) then
               info = 7
           else if (incx==0) then
               info = 9
           end if
           if (info/=0) then
               call stdlib_xerbla('ZTBMV ',info)
               return
           end if
           ! quick return if possible.
           if (n==0) return
           noconj = stdlib_lsame(trans,'T')
           nounit = stdlib_lsame(diag,'N')
           ! set up the start point in x if the increment is not unity. this
           ! will be  ( n - 1 )*incx   too small for descending loops.
           if (incx<=0) then
               kx = 1 - (n-1)*incx
           else if (incx/=1) then
               kx = 1
           end if
           ! start the operations. in this version the elements of a are
           ! accessed sequentially with cone pass through a.
           if (stdlib_lsame(trans,'N')) then
               ! form  x := a*x.
               if (stdlib_lsame(uplo,'U')) then
                   kplus1 = k + 1
                   if (incx==1) then
                       do j = 1,n
                           if (x(j)/=czero) then
                               temp = x(j)
                               l = kplus1 - j
                               do i = max(1,j-k),j - 1
                                   x(i) = x(i) + temp*a(l+i,j)
                               end do
                               if (nounit) x(j) = x(j)*a(kplus1,j)
                           end if
                       end do
                   else
                       jx = kx
                       do j = 1,n
                           if (x(jx)/=czero) then
                               temp = x(jx)
                               ix = kx
                               l = kplus1 - j
                               do i = max(1,j-k),j - 1
                                   x(ix) = x(ix) + temp*a(l+i,j)
                                   ix = ix + incx
                               end do
                               if (nounit) x(jx) = x(jx)*a(kplus1,j)
                           end if
                           jx = jx + incx
                           if (j>k) kx = kx + incx
                       end do
                   end if
               else
                   if (incx==1) then
                       do j = n,1,-1
                           if (x(j)/=czero) then
                               temp = x(j)
                               l = 1 - j
                               do i = min(n,j+k),j + 1,-1
                                   x(i) = x(i) + temp*a(l+i,j)
                               end do
                               if (nounit) x(j) = x(j)*a(1,j)
                           end if
                       end do
                   else
                       kx = kx + (n-1)*incx
                       jx = kx
                       do j = n,1,-1
                           if (x(jx)/=czero) then
                               temp = x(jx)
                               ix = kx
                               l = 1 - j
                               do i = min(n,j+k),j + 1,-1
                                   x(ix) = x(ix) + temp*a(l+i,j)
                                   ix = ix - incx
                               end do
                               if (nounit) x(jx) = x(jx)*a(1,j)
                           end if
                           jx = jx - incx
                           if ((n-j)>=k) kx = kx - incx
                       end do
                   end if
               end if
           else
              ! form  x := a**t*x  or  x := a**h*x.
               if (stdlib_lsame(uplo,'U')) then
                   kplus1 = k + 1
                   if (incx==1) then
                       do j = n,1,-1
                           temp = x(j)
                           l = kplus1 - j
                           if (noconj) then
                               if (nounit) temp = temp*a(kplus1,j)
                               do i = j - 1,max(1,j-k),-1
                                   temp = temp + a(l+i,j)*x(i)
                               end do
                           else
                               if (nounit) temp = temp*conjg(a(kplus1,j))
                               do i = j - 1,max(1,j-k),-1
                                   temp = temp + conjg(a(l+i,j))*x(i)
                               end do
                           end if
                           x(j) = temp
                       end do
                   else
                       kx = kx + (n-1)*incx
                       jx = kx
                       do j = n,1,-1
                           temp = x(jx)
                           kx = kx - incx
                           ix = kx
                           l = kplus1 - j
                           if (noconj) then
                               if (nounit) temp = temp*a(kplus1,j)
                               do i = j - 1,max(1,j-k),-1
                                   temp = temp + a(l+i,j)*x(ix)
                                   ix = ix - incx
                               end do
                           else
                               if (nounit) temp = temp*conjg(a(kplus1,j))
                               do i = j - 1,max(1,j-k),-1
                                   temp = temp + conjg(a(l+i,j))*x(ix)
                                   ix = ix - incx
                               end do
                           end if
                           x(jx) = temp
                           jx = jx - incx
                       end do
                   end if
               else
                   if (incx==1) then
                       do j = 1,n
                           temp = x(j)
                           l = 1 - j
                           if (noconj) then
                               if (nounit) temp = temp*a(1,j)
                               do i = j + 1,min(n,j+k)
                                   temp = temp + a(l+i,j)*x(i)
                               end do
                           else
                               if (nounit) temp = temp*conjg(a(1,j))
                               do i = j + 1,min(n,j+k)
                                   temp = temp + conjg(a(l+i,j))*x(i)
                               end do
                           end if
                           x(j) = temp
                       end do
                   else
                       jx = kx
                       do j = 1,n
                           temp = x(jx)
                           kx = kx + incx
                           ix = kx
                           l = 1 - j
                           if (noconj) then
                               if (nounit) temp = temp*a(1,j)
                               do i = j + 1,min(n,j+k)
                                   temp = temp + a(l+i,j)*x(ix)
                                   ix = ix + incx
                               end do
                           else
                               if (nounit) temp = temp*conjg(a(1,j))
                               do i = j + 1,min(n,j+k)
                                   temp = temp + conjg(a(l+i,j))*x(ix)
                                   ix = ix + incx
                               end do
                           end if
                           x(jx) = temp
                           jx = jx + incx
                       end do
                   end if
               end if
           end if
           return
     end subroutine stdlib_ztbmv


     pure subroutine stdlib_ztbsv(uplo,trans,diag,n,k,a,lda,x,incx)
     !! ZTBSV solves one of the systems of equations
     !! A*x = b,   or   A**T*x = b,   or   A**H*x = b,
     !! where b and x are n element vectors and A is an n by n unit, or
     !! non-unit, upper or lower triangular band matrix, with ( k + 1 )
     !! diagonals.
     !! No test for singularity or near-singularity is included in this
     !! routine. Such tests must be performed before calling this routine.
        ! -- reference blas level2 routine --
        ! -- reference blas is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           ! Scalar Arguments 
           integer(ilp), intent(in) :: incx, k, lda, n
           character, intent(in) :: diag, trans, uplo
           ! Array Arguments 
           complex(dp), intent(in) :: a(lda,*)
           complex(dp), intent(inout) :: x(*)
        ! =====================================================================
           
           ! Local Scalars 
           complex(dp) :: temp
           integer(ilp) :: i, info, ix, j, jx, kplus1, kx, l
           logical(lk) :: noconj, nounit
           ! Intrinsic Functions 
           intrinsic :: conjg,max,min
           ! test the input parameters.
           info = 0
           if (.not.stdlib_lsame(uplo,'U') .and. .not.stdlib_lsame(uplo,'L')) then
               info = 1
           else if (.not.stdlib_lsame(trans,'N') .and. .not.stdlib_lsame(trans,'T') &
                     .and..not.stdlib_lsame(trans,'C')) then
               info = 2
           else if (.not.stdlib_lsame(diag,'U') .and. .not.stdlib_lsame(diag,'N')) then
               info = 3
           else if (n<0) then
               info = 4
           else if (k<0) then
               info = 5
           else if (lda< (k+1)) then
               info = 7
           else if (incx==0) then
               info = 9
           end if
           if (info/=0) then
               call stdlib_xerbla('ZTBSV ',info)
               return
           end if
           ! quick return if possible.
           if (n==0) return
           noconj = stdlib_lsame(trans,'T')
           nounit = stdlib_lsame(diag,'N')
           ! set up the start point in x if the increment is not unity. this
           ! will be  ( n - 1 )*incx  too small for descending loops.
           if (incx<=0) then
               kx = 1 - (n-1)*incx
           else if (incx/=1) then
               kx = 1
           end if
           ! start the operations. in this version the elements of a are
           ! accessed by sequentially with cone pass through a.
           if (stdlib_lsame(trans,'N')) then
              ! form  x := inv( a )*x.
               if (stdlib_lsame(uplo,'U')) then
                   kplus1 = k + 1
                   if (incx==1) then
                       do j = n,1,-1
                           if (x(j)/=czero) then
                               l = kplus1 - j
                               if (nounit) x(j) = x(j)/a(kplus1,j)
                               temp = x(j)
                               do i = j - 1,max(1,j-k),-1
                                   x(i) = x(i) - temp*a(l+i,j)
                               end do
                           end if
                       end do
                   else
                       kx = kx + (n-1)*incx
                       jx = kx
                       do j = n,1,-1
                           kx = kx - incx
                           if (x(jx)/=czero) then
                               ix = kx
                               l = kplus1 - j
                               if (nounit) x(jx) = x(jx)/a(kplus1,j)
                               temp = x(jx)
                               do i = j - 1,max(1,j-k),-1
                                   x(ix) = x(ix) - temp*a(l+i,j)
                                   ix = ix - incx
                               end do
                           end if
                           jx = jx - incx
                       end do
                   end if
               else
                   if (incx==1) then
                       do j = 1,n
                           if (x(j)/=czero) then
                               l = 1 - j
                               if (nounit) x(j) = x(j)/a(1,j)
                               temp = x(j)
                               do i = j + 1,min(n,j+k)
                                   x(i) = x(i) - temp*a(l+i,j)
                               end do
                           end if
                       end do
                   else
                       jx = kx
                       do j = 1,n
                           kx = kx + incx
                           if (x(jx)/=czero) then
                               ix = kx
                               l = 1 - j
                               if (nounit) x(jx) = x(jx)/a(1,j)
                               temp = x(jx)
                               do i = j + 1,min(n,j+k)
                                   x(ix) = x(ix) - temp*a(l+i,j)
                                   ix = ix + incx
                               end do
                           end if
                           jx = jx + incx
                       end do
                   end if
               end if
           else
              ! form  x := inv( a**t )*x  or  x := inv( a**h )*x.
               if (stdlib_lsame(uplo,'U')) then
                   kplus1 = k + 1
                   if (incx==1) then
                       do j = 1,n
                           temp = x(j)
                           l = kplus1 - j
                           if (noconj) then
                               do i = max(1,j-k),j - 1
                                   temp = temp - a(l+i,j)*x(i)
                               end do
                               if (nounit) temp = temp/a(kplus1,j)
                           else
                               do i = max(1,j-k),j - 1
                                   temp = temp - conjg(a(l+i,j))*x(i)
                               end do
                               if (nounit) temp = temp/conjg(a(kplus1,j))
                           end if
                           x(j) = temp
                       end do
                   else
                       jx = kx
                       do j = 1,n
                           temp = x(jx)
                           ix = kx
                           l = kplus1 - j
                           if