function cabs1 ( z )
!
!*******************************************************************************
!
!! CABS1 returns the L1 norm of a complex number.
!
!
!  Modified:
!
!    22 May 2002
!
!  Parameters:
!
!    Input, complex Z, the number whose norm is desired.
!
!    Output, real CABS1, the L1 norm of Z.
!
  implicit none
!
  real cabs1
  complex z
!
  cabs1 = abs ( real ( z ) ) + abs ( aimag ( z ) )

  return
end
subroutine caxpy ( n, ca, cx, incx, cy, incy )
!
!*******************************************************************************
!
!! CAXPY computes a constant times a vector plus a vector.
!
!
!     jack dongarra, linpack, 3/11/78.
!     modified 12/3/93, array(1) declarations changed to array(*)
!
!  Reference:
!
!    Lawson, Hanson, Kincaid, Krogh,
!    Basic Linear Algebra Subprograms for Fortran Usage,
!    Algorithm 539,
!    ACM Transactions on Mathematical Software,
!    Volume 5, Number 3, September 1979, pages 308-323.
!
  implicit none
!
  complex ca
  real cabs1
  complex cx(*)
  complex cy(*)
  integer i
  integer incx
  integer incy
  integer ix
  integer iy
  integer n
!
  if ( n <= 0 ) then
    return
  end if

  if ( cabs1 ( ca ) == 0.0e+00 ) then
    return
  end if

  if ( incx /= 1 .or. incy /= 1 ) then

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( -n + 1 ) * incx + 1
    end if

    if ( incy >= 0 ) then
      iy = 1
    else
      iy = ( -n + 1 ) * incy + 1
    end if

    do i = 1,n
      cy(iy) = cy(iy) + ca*cx(ix)
      ix = ix + incx
      iy = iy + incy
    end do

  else

    cy(1:n) = cy(1:n) + ca * cx(1:n)

  end if

  return
end
subroutine ccopy ( n, cx, incx, cy, incy )
!
!*******************************************************************************
!
!! CCOPY copies a vector, x, to a vector, y.
!
!
!     jack dongarra, linpack, 3/11/78.
!     modified 12/3/93, array(1) declarations changed to array(*)
!
!  Reference:
!
!    Lawson, Hanson, Kincaid, Krogh,
!    Basic Linear Algebra Subprograms for Fortran Usage,
!    Algorithm 539,
!    ACM Transactions on Mathematical Software,
!    Volume 5, Number 3, September 1979, pages 308-323.
!
  implicit none
!
  complex cx(*)
  complex cy(*)
  integer i
  integer incx
  integer incy
  integer ix
  integer iy
  integer n
!
  if ( n <= 0 ) then
    return
  end if

  if ( incx /= 1 .or. incy /= 1 ) then

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( -n + 1 ) * incx + 1
    end if

    if ( incy >= 0 ) then
      iy = 1
    else
      iy = ( -n + 1 ) * incy + 1
    end if

    do i = 1, n
      cy(iy) = cx(ix)
      ix = ix + incx
      iy = iy + incy
    end do

  else

    cy(1:n) = cx(1:n)

  end if

  return
end
function cdotc ( n, cx, incx, cy, incy )
!
!*******************************************************************************
!
!! CDOTC forms the dot product of two vectors, conjugating the first vector.
!
!
!     jack dongarra, linpack,  3/11/78.
!     modified 12/3/93, array(1) declarations changed to array(*)
!
!  Reference:
!
!    Lawson, Hanson, Kincaid, Krogh,
!    Basic Linear Algebra Subprograms for Fortran Usage,
!    Algorithm 539,
!    ACM Transactions on Mathematical Software,
!    Volume 5, Number 3, September 1979, pages 308-323.
!
  implicit none
!
  complex cdotc
  complex ctemp
  complex cx(*)
  complex cy(*)
  integer i
  integer incx
  integer incy
  integer ix
  integer iy
  integer n
!
  ctemp = ( 0.0, 0.0 )
  cdotc = ( 0.0, 0.0 )

  if ( n <= 0 ) then
    return
  end if

  if ( incx == 1 .and. incy == 1 ) then

    do i = 1,n
      ctemp = ctemp + conjg ( cx(i) ) * cy(i)
    end do

  else

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( -n + 1 ) * incx + 1
    end if

    if ( incy >= 0 ) then
      iy = 1
    else
      iy = ( -n + 1 ) * incy + 1
    end if

    do i = 1, n
      ctemp = ctemp + conjg ( cx(ix) ) * cy(iy)
      ix = ix + incx
      iy = iy + incy
    end do

  end if

  cdotc = ctemp

  return
end
function cdotu ( n, cx, incx, cy, incy )
!
!*******************************************************************************
!
!! CDOTU forms the (unconjugated) dot product of two vectors.
!
!
!  Discussion:
!
!    Using the FORTRAN 90 function DOT_PRODUCT with complex vectors
!    as arguments will give you the conjugated dot product!
!
!     jack dongarra, linpack, 3/11/78.
!     modified 12/3/93, array(1) declarations changed to array(*)
!
!  Reference:
!
!    Lawson, Hanson, Kincaid, Krogh,
!    Basic Linear Algebra Subprograms for Fortran Usage,
!    Algorithm 539,
!    ACM Transactions on Mathematical Software,
!    Volume 5, Number 3, September 1979, pages 308-323.
!
  implicit none
!
  complex cdotu
  complex ctemp
  complex cx(*)
  complex cy(*)
  integer i
  integer incx
  integer incy
  integer ix
  integer iy
  integer n
!
  ctemp = ( 0.0, 0.0 )
  cdotu = ( 0.0, 0.0 )

  if ( n <= 0 ) then
    return
  end if

  if ( incx == 1 .and. incy == 1 ) then

    do i = 1,n
      ctemp = ctemp + cx(i) * cy(i)
    end do

  else

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( -n + 1 ) * incx + 1
    end if

    if ( incy >= 0 ) then
      iy = 1
    else
      iy = ( -n + 1 ) * incy + 1
    end if

    do i = 1,n
      ctemp = ctemp + cx(ix) * cy(iy)
      ix = ix + incx
      iy = iy + incy
    end do

  end if

  cdotu = ctemp

  return
end
subroutine cgbmv ( trans, m, n, kl, ku, alpha, a, lda, x, incx, beta, y, incy )
!
!*******************************************************************************
!
!! CGBMV  performs one of the matrix-vector operations
!
!     y := alpha*A*x + beta*y,   or   y := alpha*A'*x + beta*y,   or
!
!     y := alpha*conjg( A' )*x + beta*y,
!
!  where alpha and beta are scalars, x and y are vectors and A is an
!  m by n band matrix, with kl sub-diagonals and ku super-diagonals.
!
!  Parameters:
!
!  TRANS  - character.
!           On entry, TRANS specifies the operation to be performed as
!           follows:
!
!              TRANS = 'N' or 'n'   y := alpha*A*x + beta*y.
!
!              TRANS = 'T' or 't'   y := alpha*A'*x + beta*y.
!
!              TRANS = 'C' or 'c'   y := alpha*conjg( A' )*x + beta*y.
!
!           Unchanged on exit.
!
!  M      - integer.
!           On entry, M specifies the number of rows of the matrix A.
!           M must be at least zero.
!           Unchanged on exit.
!
!  N      - integer.
!           On entry, N specifies the number of columns of the matrix A.
!           N must be at least zero.
!           Unchanged on exit.
!
!  KL     - integer.
!           On entry, KL specifies the number of sub-diagonals of the
!           matrix A. KL must satisfy  0  <=  KL.
!           Unchanged on exit.
!
!  KU     - integer.
!           On entry, KU specifies the number of super-diagonals of the
!           matrix A. KU must satisfy  0  <=  KU.
!           Unchanged on exit.
!
!  ALPHA  - complex         .
!           On entry, ALPHA specifies the scalar alpha.
!           Unchanged on exit.
!
!  A      - complex          array of DIMENSION ( LDA, n ).
!           Before entry, the leading ( kl + ku + 1 ) by n part of the
!           array A must contain the matrix of coefficients, supplied
!           column by column, with the leading diagonal of the matrix in
!           row ( ku + 1 ) of the array, the first super-diagonal
!           starting at position 2 in row ku, the first sub-diagonal
!           starting at position 1 in row ( ku + 2 ), and so on.
!           Elements in the array A that do not correspond to elements
!           in the band matrix (such as the top left ku by ku triangle)
!           are not referenced.
!           The following program segment will transfer a band matrix
!           from conventional full matrix storage to band storage:
!
!                 do j = 1, n
!                    k = kU + 1 - J
!                    do i = max ( 1, j - KU ), min ( M, J + KL )
!                       a( K + I,j) = matrix(i,j)
!                    end do
!                 end do
!
!           Unchanged on exit.
!
!  LDA    - integer.
!           On entry, LDA specifies the first dimension of A as declared
!           in the calling (sub) program. LDA must be at least
!           ( kl + ku + 1 ).
!           Unchanged on exit.
!
!  X      - complex          array of DIMENSION at least
!           ( 1 + ( n - 1 ) * abs( incx ) ) when TRANS = 'N' or 'n'
!           and at least
!           ( 1 + ( m - 1 ) * abs( incx ) ) otherwise.
!           Before entry, the incremented array X must contain the
!           vector x.
!           Unchanged on exit.
!
!  incx   - integer.
!           On entry, incx specifies the increment for the elements of
!           X. incx must not be zero.
!           Unchanged on exit.
!
!  BETA   - complex         .
!           On entry, BETA specifies the scalar beta. When BETA is
!           supplied as zero then Y need not be set on input.
!           Unchanged on exit.
!
!  Y      - complex          array of DIMENSION at least
!           ( 1 + ( m - 1 ) * abs( incy ) ) when TRANS = 'N' or 'n'
!           and at least
!           ( 1 + ( n - 1 ) * abs( incy ) ) otherwise.
!           Before entry, the incremented array Y must contain the
!           vector y. On exit, Y is overwritten by the updated vector y.
!
!
!  incy   - integer.
!           On entry, incy specifies the increment for the elements of
!           Y. incy must not be zero.
!           Unchanged on exit.
!
  implicit none
!
  complex alpha
  complex beta
      integer            incx, incy, kl, ku, lda, m, n
      character        trans
      complex            a( lda, * ), x( * ), y( * )
      complex            one
      parameter        ( one  = ( 1.0e+0, 0.0e+0 ) )
      complex            zero
      parameter        ( zero = ( 0.0e+0, 0.0e+0 ) )
      complex            temp
      integer            i, info, ix, iy, j, jx, jy, k, kup1, kx, ky, &
                         lenx, leny
      logical            noconj
      logical            lsame
!
!  Test the input.
!
      info = 0
      if   ( .not.lsame ( trans, 'N' ) .and.  &
               .not.lsame ( trans, 'T' ) .and.  &
               .not.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 xerbla ( 'cgbmv ', info )
    return
  end if
!
!  Quick return if possible.
!
  if ( ( m == 0 ).or.( n == 0 ) .or.  &
    ( ( alpha == zero ) .and. ( beta == one ) ) ) then
    return
  end if

      noconj = 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 ( 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 one pass through the band part of A.
!
!  First form  y := beta*y.
!
      if ( beta /= one ) then
         if ( incy == 1 ) then
            if ( beta == zero ) then
               do i = 1, leny
                  y(i) = zero
               end do
            else
               do i = 1, leny
                  y(i) = beta * y(i)
               end do
            end if
         else
            iy = ky
            if ( beta == zero ) then
               do i = 1, leny
                  y(iy) = zero
                  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 == zero ) then
         return
      end if

      kup1 = ku + 1
      if ( lsame ( trans, 'N' ) ) then
!
!  Form  y := alpha*A*x + y.
!
         jx = kx
         if ( incy == 1 ) then
            do j = 1, n
               if ( x(jx) /= zero ) then
                  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
               end if
               jx = jx + incx
            end do
         else
            do j = 1, n
               if ( x(jx) /= zero ) then
                  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
               end if
               jx = jx + incx
               if ( j > ku ) then
                  ky = ky + incy
               end if
            end do
         end if
      else
!
!  Form  y := alpha*A'*x + y  or  y := alpha*conjg( A' )*x + y.
!
         jy = ky
         if ( incx == 1 ) then
            do j = 1, n
               temp = zero
               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 = zero
               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 cgemm ( transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, &
  ldc )
!
!*******************************************************************************
!
!! CGEMM 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'   or   op( X ) = conjg( X' ),
!
!  alpha and beta are scalars, and A, B and C are matrices, with op( A )
!  an m by k matrix,  op( B )  a  k by n matrix and  C an m by n matrix.
!
!  Parameters:
!
!  TRANSA - character.
!           On entry, TRANSA specifies the form of op( A ) to be used in
!           the matrix multiplication as follows:
!
!              TRANSA = 'N' or 'n',  op( A ) = A.
!
!              TRANSA = 'T' or 't',  op( A ) = A'.
!
!              TRANSA = 'C' or 'c',  op( A ) = conjg( A' ).
!
!           Unchanged on exit.
!
!  TRANSB - character.
!           On entry, TRANSB specifies the form of op( B ) to be used in
!           the matrix multiplication as follows:
!
!              TRANSB = 'N' or 'n',  op( B ) = B.
!
!              TRANSB = 'T' or 't',  op( B ) = B'.
!
!              TRANSB = 'C' or 'c',  op( B ) = conjg( B' ).
!
!           Unchanged on exit.
!
!  M      - integer.
!           On entry,  M  specifies  the number  of rows  of the  matrix
!           op( A )  and of the  matrix C.  M  must  be at least  zero.
!           Unchanged on exit.
!
!  N      - integer.
!           On entry,  N  specifies the number  of columns of the matrix
!           op( B ) and the number of columns of the matrix C. N must be
!           at least zero.
!           Unchanged on exit.
!
!  K      - integer.
!           On entry,  K  specifies  the number of columns of the matrix
!           op( A ) and the number of rows of the matrix op( B ). K must
!           be at least  zero.
!           Unchanged on exit.
!
!  ALPHA  - complex         .
!           On entry, ALPHA specifies the scalar alpha.
!           Unchanged on exit.
!
!  A      - complex          array of DIMENSION ( LDA, ka ), where ka is
!           k  when  TRANSA = 'N' or 'n',  and is  m  otherwise.
!           Before entry with  TRANSA = 'N' or 'n',  the leading  m by k
!           part of the array  A  must contain the matrix A,  otherwise
!           the leading  k by m  part of the array  A  must contain  the
!           matrix A.
!           Unchanged on exit.
!
!  LDA    - integer.
!           On entry, LDA specifies the first dimension of A as declared
!           in the calling (sub) program. When  TRANSA = 'N' or 'n' then
!           LDA must be at least  max( 1, m ), otherwise  LDA must be at
!           least  max( 1, k ).
!           Unchanged on exit.
!
!  B      - complex          array of DIMENSION ( LDB, kb ), where kb is
!           n  when  TRANSB = 'N' or 'n',  and is  k  otherwise.
!           Before entry with  TRANSB = 'N' or 'n',  the leading  k by n
!           part of the array  B  must contain the matrix B,  otherwise
!           the leading  n by k  part of the array  B  must contain  the
!           matrix B.
!           Unchanged on exit.
!
!  LDB    - integer.
!           On entry, LDB specifies the first dimension of B as declared
!           in the calling (sub) program. When  TRANSB = 'N' or 'n' then
!           LDB must be at least  max( 1, k ), otherwise  LDB must be at
!           least  max( 1, n ).
!           Unchanged on exit.
!
!  BETA   - complex         .
!           On entry,  BETA  specifies the scalar  beta.  When  BETA  is
!           supplied as zero then C need not be set on input.
!           Unchanged on exit.
!
!  C      - complex          array of DIMENSION ( LDC, n ).
!           Before entry, the leading  m by n  part of the array  C must
!           contain the matrix C,  except when  beta  is zero, in which
!           case C need not be set on entry.
!           On exit, the array  C  is overwritten by the  m by n  matrix
!           ( alpha*op( A )*op( B ) + beta*C ).
!
!  LDC    - integer.
!           On entry, LDC specifies the first dimension of C as declared
!           in  the  calling  (sub)  program.   LDC  must  be  at  least
!           max( 1, m ).
!           Unchanged on exit.
!
  implicit none
!
  character transa
  character transb
      integer            m, n, k, lda, ldb, ldc
      complex            alpha, beta
      complex            a( lda, * ), b( ldb, * ), c( ldc, * )
      logical            lsame
      logical            conja, conjb, nota, notb
      integer            i, info, j, l, ncola, nrowa, nrowb
      complex            temp
      complex            one
      parameter        ( one  = ( 1.0e+0, 0.0e+0 ) )
      complex            zero
      parameter        ( zero = ( 0.0e+0, 0.0e+0 ) )
!
!     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, NCOLA and  NROWB  as the number of rows and  columns  of  A
!     and the number of rows of  B  respectively.
!
      nota  = lsame ( transa, 'N' )
      notb  = lsame ( transb, 'N' )
      conja = lsame ( transa, 'C' )
      conjb = lsame ( transb, 'C' )
      if ( nota ) then
         nrowa = m
         ncola = k
      else
         nrowa = k
         ncola = m
      end if

      if ( notb ) then
         nrowb = k
      else
         nrowb = n
      end if
!
!  Test the input.
!
      info = 0
      if (      ( .not.nota                 ) .and.  &
               ( .not.conja                ) .and.  &
               ( .not.lsame ( transa, 'T' ) )      ) then
         info = 1
      else if ( ( .not.notb                 ) .and.  &
               ( .not.conjb                ) .and.  &
               ( .not.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 xerbla ( 'cgemm ', info )
         return
      end if
!
!  Quick return if possible.
!
      if ( ( m == 0 ).or.( n == 0 ) .or.  &
          ( ( ( alpha == zero ).or.( k == 0 ) ) .and. ( beta == one ) ) ) then
         return
      end if
!
!  And when ALPHA == zero.
!
      if ( alpha == zero ) then
         if ( beta == zero ) then
            c(1:m,1:n) = zero
         else
            c(1:m,1:n) = beta * c(1:m,1:n)
         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 == zero ) then
                  do i = 1, m
                     c(i,j) = zero
                  end do
               else if ( beta/= one ) then
                  do i = 1, m
                     c(i,j) = beta * c(i,j)
                  end do
               end if
               do l = 1, k
                  if ( b(l,j) /= zero ) then
                     temp = alpha * b(l,j)
                     do i = 1, m
                        c(i,j) = c(i,j) + temp * a(i,l)
                     end do
                  end if
               end do
            end do
         else if ( conja ) then
!
!  Form  C := alpha*conjg( A' ) * b + beta*C.
!
            do j = 1, n
               do i = 1, m
                  temp = zero
                  do l = 1, k
                     temp = temp + conjg ( a(l,i) ) * b(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
         else
!
!  Form  C := alpha*A' * b + beta*C
!
            do j = 1, n
               do i = 1, m
                  temp = zero
                  do l = 1, k
                     temp = temp + a(l,i) * b(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
      else if ( nota ) then
         if ( conjb ) then
!
!  Form  C := alpha*A*conjg( B' ) + beta*C.
!
            do j = 1, n
               if ( beta == zero ) then
                  do i = 1, m
                     c(i,j) = zero
                  end do
               else if ( beta/= one ) then
                  do i = 1, m
                     c(i,j) = beta * c(i,j)
                  end do
               end if
               do l = 1, k
                  if ( b(j,l) /= zero ) then
                     temp = alpha * conjg ( b(j,l) )
                     do i = 1, m
                        c(i,j) = c(i,j) + temp * a(i,l)
                     end do
                  end if
               end do
            end do
         else
!
!  Form  C := alpha*A * b' + beta*C
!
            do j = 1, n
               if ( beta == zero ) then
                  do i = 1, m
                     c(i,j) = zero
                  end do
               else if ( beta /= one ) then
                  do i = 1, m
                     c(i,j) = beta * c(i,j)
                  end do
               end if
               do l = 1, k
                  if ( b(j,l) /= zero ) then
                     temp = alpha * b(j,l)
                     do i = 1, m
                        c(i,j) = c(i,j) + temp * a(i,l)
                     end do
                  end if
               end do
            end do
         end if
      else if ( conja ) then
         if ( conjb ) then
!
!  Form C := alpha*conjg( A' )*conjg( B' ) + beta*C.
!
            do j = 1, n
               do i = 1, m
                  temp = zero
                  do l = 1, k
                     temp = temp + conjg ( a(l,i) ) * conjg ( b(j,l) )
                  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
         else
!
!  Form C := alpha*conjg( A' ) * b' + beta*C
!
            do j = 1, n
               do i = 1, m
                  temp = zero
                  do l = 1, k
                     temp = temp + conjg ( a(l,i) ) * b(j,l)
                  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
      else
         if ( conjb ) then
!
!  Form  C := alpha*A'*conjg( B' ) + beta*C
!
            do j = 1, n
               do i = 1, m
                  temp = zero
                  do l = 1, k
                     temp = temp + a(l,i)*conjg ( b(j,l) )
                  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
         else
!
!  Form  C := alpha*A' * b' + beta*C
!
            do j = 1, n
               do i = 1, m
                  temp = zero
                  do l = 1, k
                     temp = temp + a(l,i) * b(j,l)
                  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 cgemv ( trans, m, n, alpha, a, lda, x, incx, beta, y, incy )
!
!*******************************************************************************
!
!! CGEMV performs one of the matrix-vector operations
!
!     y := alpha*A*x + beta*y,   or   y := alpha*A'*x + beta*y,   or
!
!     y := alpha*conjg( A' )*x + beta*y,
!
!  where alpha and beta are scalars, x and y are vectors and A is an
!  m by n matrix.
!
!  Parameters:
!
!  TRANS  - character.
!           On entry, TRANS specifies the operation to be performed as
!           follows:
!
!              TRANS = 'N' or 'n'   y := alpha*A*x + beta*y.
!
!              TRANS = 'T' or 't'   y := alpha*A'*x + beta*y.
!
!              TRANS = 'C' or 'c'   y := alpha*conjg( A' )*x + beta*y.
!
!           Unchanged on exit.
!
!  M      - integer.
!           On entry, M specifies the number of rows of the matrix A.
!           M must be at least zero.
!           Unchanged on exit.
!
!  N      - integer.
!           On entry, N specifies the number of columns of the matrix A.
!           N must be at least zero.
!           Unchanged on exit.
!
!  ALPHA  - complex         .
!           On entry, ALPHA specifies the scalar alpha.
!           Unchanged on exit.
!
!  A      - complex          array of DIMENSION ( LDA, n ).
!           Before entry, the leading m by n part of the array A must
!           contain the matrix of coefficients.
!           Unchanged on exit.
!
!  LDA    - integer.
!           On entry, LDA specifies the first dimension of A as declared
!           in the calling (sub) program. LDA must be at least
!           max( 1, m ).
!           Unchanged on exit.
!
!  X      - complex          array of DIMENSION at least
!           ( 1 + ( n - 1 ) * abs( incx ) ) when TRANS = 'N' or 'n'
!           and at least
!           ( 1 + ( m - 1 ) * abs( incx ) ) otherwise.
!           Before entry, the incremented array X must contain the
!           vector x.
!           Unchanged on exit.
!
!  incx   - integer.
!           On entry, incx specifies the increment for the elements of
!           X. incx must not be zero.
!           Unchanged on exit.
!
!  BETA   - complex         .
!           On entry, BETA specifies the scalar beta. When BETA is
!           supplied as zero then Y need not be set on input.
!           Unchanged on exit.
!
!  Y      - complex          array of DIMENSION at least
!           ( 1 + ( m - 1 ) * abs( incy ) ) when TRANS = 'N' or 'n'
!           and at least
!           ( 1 + ( n - 1 ) * abs( incy ) ) otherwise.
!           Before entry with BETA non-zero, the incremented array Y
!           must contain the vector y. On exit, Y is overwritten by the
!           updated vector y.
!
!  incy   - integer.
!           On entry, incy specifies the increment for the elements of
!           Y. incy must not be zero.
!           Unchanged on exit.
!
  implicit none
!
  complex            alpha, beta
      integer            incx, incy, lda, m, n
      character        trans
      complex            a( lda, * ), x( * ), y( * )
      complex            one
      parameter        ( one  = ( 1.0e+0, 0.0e+0 ) )
      complex            zero
      parameter        ( zero = ( 0.0e+0, 0.0e+0 ) )
      complex            temp
      integer            i, info, ix, iy, j, jx, jy, kx, ky, lenx, leny
      logical            noconj
      logical            lsame
!
!  Test the input.
!
      info = 0
      if   ( .not.lsame ( trans, 'N' ) .and.  &
               .not.lsame ( trans, 'T' ) .and.  &
               .not.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 xerbla ( 'cgemv ', info )
         return
      end if
!
!  Quick return if possible.
!
      if ( ( m == 0 ).or.( n == 0 ) .or.  &
          ( ( alpha == zero ) .and. ( beta == one ) ) ) then
         return
      end if

      noconj = 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 ( 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 one pass through A.
!
!  First form  y := beta*y.
!
      if ( beta /= one ) then
         if ( incy == 1 ) then
            if ( beta == zero ) then
               do i = 1, leny
                  y(i) = zero
               end do
            else
               do i = 1, leny
                  y(i) = beta * y(i)
               end do
            end if
         else
            iy = ky
            if ( beta == zero ) then
               do i = 1, leny
                  y(iy) = zero
                  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 == zero ) then
         return
      end if

      if ( lsame ( trans, 'N' ) ) then
!
!  Form  y := alpha*A*x + y.
!
         jx = kx
         if ( incy == 1 ) then
            do j = 1, n
               if ( x(jx) /= zero ) then
                  temp = alpha * x(jx)
                  do i = 1, m
                     y(i) = y(i) + temp * a(i,j)
                  end do
               end if
               jx = jx + incx
            end do
         else
            do j = 1, n
               if ( x(jx) /= zero ) then
                  temp = alpha * x(jx)
                  iy = ky
                  do i = 1, m
                     y(iy) = y(iy) + temp * a(i,j)
                     iy = iy + incy
                  end do
               end if
               jx = jx + incx
            end do
         end if
      else
!
!  Form  y := alpha*A'*x + y  or  y := alpha*conjg( A' )*x + y.
!
         jy = ky
         if ( incx == 1 ) then
            do j = 1, n
               temp = zero
               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 = zero
               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 cgerc ( m, n, alpha, x, incx, y, incy, a, lda )
!
!*******************************************************************************
!
!! CGERC performs the rank 1 operation A := alpha*x*conjg( y' ) + A.
!
!
!  Discussion:
!
!    alpha is a scalar, x is an m element vector, y is an n element
!    vector and A is an m by n matrix.
!
!  Parameters:
!
!  M      - integer.
!           On entry, M specifies the number of rows of the matrix A.
!           M must be at least zero.
!           Unchanged on exit.
!
!  N      - integer.
!           On entry, N specifies the number of columns of the matrix A.
!           N must be at least zero.
!           Unchanged on exit.
!
!  ALPHA  - complex         .
!           On entry, ALPHA specifies the scalar alpha.
!           Unchanged on exit.
!
!  X      - complex          array of dimension at least
!           ( 1 + ( m - 1 ) * abs( incx ) ).
!           Before entry, the incremented array X must contain the m
!           element vector x.
!           Unchanged on exit.
!
!  incx   - integer.
!           On entry, incx specifies the increment for the elements of
!           X. incx must not be zero.
!           Unchanged on exit.
!
!  Y      - complex          array of dimension at least
!           ( 1 + ( n - 1 ) * abs( incy ) ).
!           Before entry, the incremented array Y must contain the n
!           element vector y.
!           Unchanged on exit.
!
!  incy   - integer.
!           On entry, incy specifies the increment for the elements of
!           Y. incy must not be zero.
!           Unchanged on exit.
!
!  A      - complex          array of DIMENSION ( LDA, n ).
!           Before entry, the leading m by n part of the array A must
!           contain the matrix of coefficients. On exit, A is
!           overwritten by the updated matrix.
!
!  LDA    - integer.
!           On entry, LDA specifies the first dimension of A as declared
!           in the calling (sub) program. LDA must be at least
!           max( 1, m ).
!           Unchanged on exit.
!
  implicit none
!
  complex alpha
  integer incx
  integer incy
  integer lda
  integer m
  integer n
!     .. Array Arguments ..
      complex            a( lda, * ), x( * ), y( * )
!
      complex            zero
      parameter        ( zero = ( 0.0e+0, 0.0e+0 ) )
!     .. Local Scalars ..
      complex            temp
      integer            i, info, ix, j, jy, kx
!
!  Test the input.
!
      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 xerbla ( 'cgerc ', info )
         return
      end if
!
!  Quick return if possible.
!
  if ( ( m == 0 ).or.( n == 0 ) .or. ( alpha == zero ) ) then
    return
  end if
!
!  Start the operations. In this version the elements of A are
!  accessed sequentially with one 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) /= zero ) 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) /= zero ) 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 cgeru ( m, n, alpha, x, incx, y, incy, a, lda )
!
!*******************************************************************************
!
!! CGERU  performs the rank 1 operation
!
!     A := alpha*x*y' + A,
!
!  where alpha is a scalar, x is an m element vector, y is an n element
!  vector and A is an m by n matrix.
!
!  Parameters:
!
!  M      - integer.
!           On entry, M specifies the number of rows of the matrix A.
!           M must be at least zero.
!           Unchanged on exit.
!
!  N      - integer.
!           On entry, N specifies the number of columns of the matrix A.
!           N must be at least zero.
!           Unchanged on exit.
!
!  ALPHA  - complex         .
!           On entry, ALPHA specifies the scalar alpha.
!           Unchanged on exit.
!
!  X      - complex          array of dimension at least
!           ( 1 + ( m - 1 ) * abs( incx ) ).
!           Before entry, the incremented array X must contain the m
!           element vector x.
!           Unchanged on exit.
!
!  incx   - integer.
!           On entry, incx specifies the increment for the elements of
!           X. incx must not be zero.
!           Unchanged on exit.
!
!  Y      - complex          array of dimension at least
!           ( 1 + ( n - 1 ) * abs( incy ) ).
!           Before entry, the incremented array Y must contain the n
!           element vector y.
!           Unchanged on exit.
!
!  incy   - integer.
!           On entry, incy specifies the increment for the elements of
!           Y. incy must not be zero.
!           Unchanged on exit.
!
!  A      - complex          array of DIMENSION ( LDA, n ).
!           Before entry, the leading m by n part of the array A must
!           contain the matrix of coefficients. On exit, A is
!           overwritten by the updated matrix.
!
!  LDA    - integer.
!           On entry, LDA specifies the first dimension of A as declared
!           in the calling (sub) program. LDA must be at least
!           max( 1, m ).
!           Unchanged on exit.
!
!
!  Level 2 Blas routine.
!
!  -- Written on 22-October-1986.
!     Jack Dongarra, Argonne National Lab.
!     Jeremy Du Croz, Nag Central Office.
!     Sven Hammarling, Nag Central Office.
!     Richard Hanson, Sandia National Labs.
!
  implicit none
!
!     .. Scalar Arguments ..
      complex            alpha
      integer            incx, incy, lda, m, n
!     .. Array Arguments ..
      complex            a( lda, * ), x( * ), y( * )
!
      complex            zero
      parameter        ( zero = ( 0.0e+0, 0.0e+0 ) )
!     .. Local Scalars ..
      complex            temp
      integer            i, info, ix, j, jy, kx
!
!  Test the input.
!
      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 xerbla ( 'cgeru ', info )
         return
      end if
!
!  Quick return if possible.
!
      if ( ( m == 0 ).or.( n == 0 ) .or. ( alpha == zero ) ) then
         return
      end if
!
!  Start the operations. In this version the elements of A are
!  accessed sequentially with one 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) /= zero ) 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) /= zero ) 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 chbmv ( uplo, n, k, alpha, a, lda, x, incx, beta, y, incy )
!
!*******************************************************************************
!
!! CHBMV 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.
!
!  Parameters:
!
!  UPLO   - character.
!           On entry, UPLO specifies whether the upper or lower
!           triangular part of the band matrix A is being supplied as
!           follows:
!
!              UPLO = 'U' or 'u'   The upper triangular part of A is
!                                  being supplied.
!
!              UPLO = 'L' or 'l'   The lower triangular part of A is
!                                  being supplied.
!
!           Unchanged on exit.
!
!  N      - integer.
!           On entry, N specifies the order of the matrix A.
!           N must be at least zero.
!           Unchanged on exit.
!
!  K      - integer.
!           On entry, K specifies the number of super-diagonals of the
!           matrix A. K must satisfy  0  <=  K.
!           Unchanged on exit.
!
!  ALPHA  - complex         .
!           On entry, ALPHA specifies the scalar alpha.
!           Unchanged on exit.
!
!  A      - complex          array of DIMENSION ( LDA, n ).
!           Before entry with UPLO = 'U' or 'u', the leading ( k + 1 )
!           by n part of the array A must contain the upper triangular
!           band part of the hermitian matrix, supplied column by
!           column, with the leading diagonal of the matrix in row
!           ( k + 1 ) of the array, the first super-diagonal starting at
!           position 2 in row k, and so on. The top left k by k triangle
!           of the array A is not referenced.
!           The following program segment will transfer the upper
!           triangular part of a hermitian band matrix from conventional
!           full matrix storage to band storage:
!
!                 do j = 1, n
!                    m = k + 1 - J
!                    do i = max ( 1, j - K ), J
!                       a( M + I,j) = matrix(i,j)
!                    end do
!                 end do
!
!           Before entry with UPLO = 'L' or 'l', the leading ( k + 1 )
!           by n part of the array A must contain the lower triangular
!           band part of the hermitian matrix, supplied column by
!           column, with the leading diagonal of the matrix in row 1 of
!           the array, the first sub-diagonal starting at position 1 in
!           row 2, and so on. The bottom right k by k triangle of the
!           array A is not referenced.
!           The following program segment will transfer the lower
!           triangular part of a hermitian band matrix from conventional
!           full matrix storage to band storage:
!
!                 do j = 1, n
!                    m = 1 - J
!                    do i = J, min ( N, j + K )
!                       a( M + I,j) = matrix(i,j)
!                    end do
!                 end do
!
!           Note that the imaginary parts of the diagonal elements need
!           not be set and are assumed to be zero.
!           Unchanged on exit.
!
!  LDA    - integer.
!           On entry, LDA specifies the first dimension of A as declared
!           in the calling (sub) program. LDA must be at least
!           ( k + 1 ).
!           Unchanged on exit.
!
!  X      - complex          array of DIMENSION at least
!           ( 1 + ( n - 1 ) * abs( incx ) ).
!           Before entry, the incremented array X must contain the
!           vector x.
!           Unchanged on exit.
!
!  incx   - integer.
!           On entry, incx specifies the increment for the elements of
!           X. incx must not be zero.
!           Unchanged on exit.
!
!  BETA   - complex         .
!           On entry, BETA specifies the scalar beta.
!           Unchanged on exit.
!
!  Y      - complex          array of DIMENSION at least
!           ( 1 + ( n - 1 ) * abs( incy ) ).
!           Before entry, the incremented array Y must contain the
!           vector y. On exit, Y is overwritten by the updated vector y.
!
!  incy   - integer.
!           On entry, incy specifies the increment for the elements of
!           Y. incy must not be zero.
!           Unchanged on exit.
!
!
  implicit none
!
      complex            alpha, beta
      integer            incx, incy, k, lda, n
      character        uplo
      complex            a( lda, * ), x( * ), y( * )
      complex            one
      parameter        ( one  = ( 1.0e+0, 0.0e+0 ) )
      complex            zero
      parameter        ( zero = ( 0.0e+0, 0.0e+0 ) )
      complex            temp1, temp2
      integer            i, info, ix, iy, j, jx, jy, kplus1, kx, ky, l
      logical            lsame
!
!  Test the input.
!
      info = 0
      if   ( .not.lsame ( uplo, 'U' ) .and.  &
               .not.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 xerbla ( 'chbmv ', info )
         return
      end if
!
!  Quick return if possible.
!
      if ( ( n == 0 ).or.( ( alpha == zero ) .and. ( beta == one ) ) ) then
         return
      end if
!
!  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 one pass through A.
!
!  First form  y := beta*y.
!
      if ( beta/= one ) then
         if ( incy == 1 ) then
            if ( beta == zero ) then
               y(1:n) = zero
            else
               do i = 1, n
                  y(i) = beta * y(i)
               end do
            end if
         else
            iy = ky
            if ( beta == zero ) then
               do i = 1, n
                  y(iy) = zero
                  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 == zero ) then
         return
      end if

      if ( 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 = zero
               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) ) &
                               + alpha * temp2
            end do
         else
            jx = kx
            jy = ky
            do j = 1, n
               temp1 = alpha * x(jx)
               temp2 = zero
               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) ) &
                                 + 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 = zero
               y(j) = y(j) + temp1 * real( a( 1,j) )
               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   = zero
               y(jy) = y(jy) + temp1 * real( a( 1,j) )
               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 chemm ( side, uplo, m, n, alpha, a, lda, b, ldb, beta, c, ldc )
!
!*******************************************************************************
!
!! CHEMM 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.
!
!  Parameters:
!
!  SIDE   - character.
!           On entry,  SIDE  specifies whether  the  hermitian matrix A
!           appears on the  left or right  in the  operation as follows:
!
!              SIDE = 'L' or 'l'   C := alpha*A * b + beta*C,
!
!              SIDE = 'R' or 'r'   C := alpha * b*A + beta*C,
!
!           Unchanged on exit.
!
!  UPLO   - character.
!           On  entry,   UPLO  specifies  whether  the  upper  or  lower
!           triangular  part  of  the  hermitian  matrix A  is  to  be
!           referenced as follows:
!
!              UPLO = 'U' or 'u'   Only the upper triangular part of the
!                                  hermitian matrix is to be referenced.
!
!              UPLO = 'L' or 'l'   Only the lower triangular part of the
!                                  hermitian matrix is to be referenced.
!
!           Unchanged on exit.
!
!  M      - integer.
!           On entry,  M  specifies the number of rows of the matrix C.
!           M  must be at least zero.
!           Unchanged on exit.
!
!  N      - integer.
!           On entry, N specifies the number of columns of the matrix C.
!           N  must be at least zero.
!           Unchanged on exit.
!
!  ALPHA  - complex         .
!           On entry, ALPHA specifies the scalar alpha.
!           Unchanged on exit.
!
!  A      - complex          array of DIMENSION ( LDA, ka ), where ka is
!           m  when  SIDE = 'L' or 'l'  and is n  otherwise.
!           Before entry  with  SIDE = 'L' or 'l',  the  m by m  part of
!           the array  A  must contain the  hermitian matrix,  such that
!           when  UPLO = 'U' or 'u', the leading m by m upper triangular
!           part of the array  A  must contain the upper triangular part
!           of the  hermitian matrix and the  strictly  lower triangular
!           part of  A  is not referenced,  and when  UPLO = 'L' or 'l',
!           the leading  m by m  lower triangular part  of the  array  A
!           must  contain  the  lower triangular part  of the  hermitian
!           matrix and the  strictly upper triangular part of  A  is not
!           referenced.
!           Before entry  with  SIDE = 'R' or 'r',  the  n by n  part of
!           the array  A  must contain the  hermitian matrix,  such that
!           when  UPLO = 'U' or 'u', the leading n by n upper triangular
!           part of the array  A  must contain the upper triangular part
!           of the  hermitian matrix and the  strictly  lower triangular
!           part of  A  is not referenced,  and when  UPLO = 'L' or 'l',
!           the leading  n by n  lower triangular part  of the  array  A
!           must  contain  the  lower triangular part  of the  hermitian
!           matrix and the  strictly upper triangular part of  A  is not
!           referenced.
!           Note that the imaginary parts  of the diagonal elements need
!           not be set, they are assumed to be zero.
!           Unchanged on exit.
!
!  LDA    - integer.
!           On entry, LDA specifies the first dimension of A as declared
!           in the  calling (sub) program. When  SIDE = 'L' or 'l'  then
!           LDA must be at least  max( 1, m ), otherwise  LDA must be at
!           least max( 1, n ).
!           Unchanged on exit.
!
!  B      - complex          array of DIMENSION ( LDB, n ).
!           Before entry, the leading  m by n part of the array  B  must
!           contain the matrix B.
!           Unchanged on exit.
!
!  LDB    - integer.
!           On entry, LDB specifies the first dimension of B as declared
!           in  the  calling  (sub)  program.   LDB  must  be  at  least
!           max( 1, m ).
!           Unchanged on exit.
!
!  BETA   - complex         .
!           On entry,  BETA  specifies the scalar  beta.  When  BETA  is
!           supplied as zero then C need not be set on input.
!           Unchanged on exit.
!
!  C      - complex          array of DIMENSION ( LDC, n ).
!           Before entry, the leading  m by n  part of the array  C must
!           contain the matrix C,  except when  beta  is zero, in which
!           case C need not be set on entry.
!           On exit, the array  C  is overwritten by the  m by n updated
!           matrix.
!
!  LDC    - integer.
!           On entry, LDC specifies the first dimension of C as declared
!           in  the  calling  (sub)  program.   LDC  must  be  at  least
!           max( 1, m ).
!           Unchanged on exit.
!
  implicit none
!
      character        side, uplo
      integer            m, n, lda, ldb, ldc
      complex            alpha, beta
      complex            a( lda, * ), b( ldb, * ), c( ldc, * )
      logical            lsame
      logical            upper
      integer            i, info, j, k, nrowa
      complex            temp1, temp2
      complex            one
      parameter        ( one  = ( 1.0e+0, 0.0e+0 ) )
      complex            zero
      parameter        ( zero = ( 0.0e+0, 0.0e+0 ) )
!
!     Set NROWA as the number of rows of A.
!
      if ( lsame ( side, 'L' ) ) then
         nrowa = m
      else
         nrowa = n
      end if
      upper = lsame ( uplo, 'U' )
!
!  Test the input.
!
      info = 0

      if (      ( .not.lsame ( side, 'L' ) ) .and.  &
               ( .not.lsame ( side, 'R' ) )      ) then
         info = 1
      else if ( ( .not.upper              ) .and.  &
               ( .not.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 xerbla ( 'chemm ', info )
         return
      end if
!
!  Quick return if possible.
!
      if ( ( m == 0 ).or.( n == 0 ) .or.  &
          ( ( alpha == zero ) .and. ( beta == one ) ) ) then
         return
      end if
!
!  And when  alpha == zero.
!
      if ( alpha == zero ) then
         if ( beta == zero ) then
            c(1:m,1:n) = zero
         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 ( 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 = zero
                  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 == zero ) then
                     c(i,j) = temp1 * real( a(i,i) ) + &
                                 alpha * temp2
                  else
                     c(i,j) = beta * c(i,j)         + &
                                 temp1 * real( 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 = zero
                  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 == zero ) then
                     c(i,j) = temp1 * real( a(i,i) ) + &
                                 alpha * temp2
                  else
                     c(i,j) = beta * c(i,j)         + &
                                 temp1 * real( a(i,i) ) + &
                                 alpha * temp2
                  end if
               end do
            end do
         end if
      else
!
!  Form  C := alpha * b*A + beta*C.
!
         do j = 1, n
            temp1 = alpha * real( a(j,j) )
            if ( beta == zero ) 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
      end if

  return
end
subroutine chemv ( uplo, n, alpha, a, lda, x, incx, beta, y, incy )
!
!*******************************************************************************
!
!! CHEMV 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.
!
!  Parameters:
!
!  UPLO   - character.
!           On entry, UPLO specifies whether the upper or lower
!           triangular part of the array A is to be referenced as
!           follows:
!
!              UPLO = 'U' or 'u'   Only the upper triangular part of A
!                                  is to be referenced.
!
!              UPLO = 'L' or 'l'   Only the lower triangular part of A
!                                  is to be referenced.
!
!           Unchanged on exit.
!
!  N      - integer.
!           On entry, N specifies the order of the matrix A.
!           N must be at least zero.
!           Unchanged on exit.
!
!  ALPHA  - complex         .
!           On entry, ALPHA specifies the scalar alpha.
!           Unchanged on exit.
!
!  A      - complex          array of DIMENSION ( LDA, n ).
!           Before entry with  UPLO = 'U' or 'u', the leading n by n
!           upper triangular part of the array A must contain the upper
!           triangular part of the hermitian matrix and the strictly
!           lower triangular part of A is not referenced.
!           Before entry with UPLO = 'L' or 'l', the leading n by n
!           lower triangular part of the array A must contain the lower
!           triangular part of the hermitian matrix and the strictly
!           upper triangular part of A is not referenced.
!           Note that the imaginary parts of the diagonal elements need
!           not be set and are assumed to be zero.
!           Unchanged on exit.
!
!  LDA    - integer.
!           On entry, LDA specifies the first dimension of A as declared
!           in the calling (sub) program. LDA must be at least
!           max( 1, n ).
!           Unchanged on exit.
!
!  X      - complex          array of dimension at least
!           ( 1 + ( n - 1 ) * abs( incx ) ).
!           Before entry, the incremented array X must contain the n
!           element vector x.
!           Unchanged on exit.
!
!  incx   - integer.
!           On entry, incx specifies the increment for the elements of
!           X. incx must not be zero.
!           Unchanged on exit.
!
!  BETA   - complex         .
!           On entry, BETA specifies the scalar beta. When BETA is
!           supplied as zero then Y need not be set on input.
!           Unchanged on exit.
!
!  Y      - complex          array of dimension at least
!           ( 1 + ( n - 1 ) * abs( incy ) ).
!           Before entry, the incremented array Y must contain the n
!           element vector y. On exit, Y is overwritten by the updated
!           vector y.
!
!  incy   - integer.
!           On entry, incy specifies the increment for the elements of
!           Y. incy must not be zero.
!           Unchanged on exit.
!
  implicit none
!
  complex alpha
  complex beta
      integer            incx, incy, lda, n
      character        uplo
      complex            a( lda, * ), x( * ), y( * )
      complex            one
      parameter        ( one  = ( 1.0e+0, 0.0e+0 ) )
      complex            zero
      parameter        ( zero = ( 0.0e+0, 0.0e+0 ) )
      complex            temp1, temp2
      integer            i, info, ix, iy, j, jx, jy, kx, ky
      logical            lsame
!
!  Test the input.
!
      info = 0
      if   ( .not.lsame ( uplo, 'U' ) .and.  &
               .not.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 xerbla ( 'chemv ', info )
         return
      end if
!
!  Quick return if possible.
!
      if ( ( n == 0 ).or.( ( alpha == zero ) .and. ( beta == one ) ) ) then
         return
      end if
!
!  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 one pass through the triangular part
!  of A.
!
!     First form  y := beta*y.
!
      if ( beta/= one ) then
         if ( incy == 1 ) then
            if ( beta == zero ) then
               y(1:n) = zero
            else
               do i = 1, n
                  y(i) = beta * y(i)
               end do
            end if
         else
            iy = ky
            if ( beta == zero ) then
               do i = 1, n
                  y(iy) = zero
                  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 == zero ) &
         return
      if ( 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 = zero
               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) ) + alpha * temp2
            end do
         else
            jx = kx
            jy = ky
            do j = 1, n
               temp1 = alpha * x(jx)
               temp2 = zero
               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) ) + 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 = zero
               y(j) = y(j) + temp1 * real( a(j,j) )
               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   = zero
               y(jy) = y(jy) + temp1 * real( a(j,j) )
               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 cher ( uplo, n, alpha, x, incx, a, lda )
!
!*******************************************************************************
!
!! CHER performs the hermitian rank 1 operation
!
!     A := alpha*x*conjg( x' ) + A,
!
!  where alpha is a real scalar, x is an n element vector and A is an
!  n by n hermitian matrix.
!
!  Parameters:
!
!  UPLO   - character.
!           On entry, UPLO specifies whether the upper or lower
!           triangular part of the array A is to be referenced as
!           follows:
!
!              UPLO = 'U' or 'u'   Only the upper triangular part of A
!                                  is to be referenced.
!
!              UPLO = 'L' or 'l'   Only the lower triangular part of A
!                                  is to be referenced.
!
!           Unchanged on exit.
!
!  N      - integer.
!           On entry, N specifies the order of the matrix A.
!           N must be at least zero.
!           Unchanged on exit.
!
!  ALPHA  - real            .
!           On entry, ALPHA specifies the scalar alpha.
!           Unchanged on exit.
!
!  X      - complex          array of dimension at least
!           ( 1 + ( n - 1 ) * abs( incx ) ).
!           Before entry, the incremented array X must contain the n
!           element vector x.
!           Unchanged on exit.
!
!  incx   - integer.
!           On entry, incx specifies the increment for the elements of
!           X. incx must not be zero.
!           Unchanged on exit.
!
!  A      - complex          array of DIMENSION ( LDA, n ).
!           Before entry with  UPLO = 'U' or 'u', the leading n by n
!           upper triangular part of the array A must contain the upper
!           triangular part of the hermitian matrix and the strictly
!           lower triangular part of A is not referenced. On exit, the
!           upper triangular part of the array A is overwritten by the
!           upper triangular part of the updated matrix.
!           Before entry with UPLO = 'L' or 'l', the leading n by n
!           lower triangular part of the array A must contain the lower
!           triangular part of the hermitian matrix and the strictly
!           upper triangular part of A is not referenced. On exit, the
!           lower triangular part of the array A is overwritten by the
!           lower triangular part of the updated matrix.
!           Note that the imaginary parts of the diagonal elements need
!           not be set, they are assumed to be zero, and on exit they
!           are set to zero.
!
!  LDA    - integer.
!           On entry, LDA specifies the first dimension of A as declared
!           in the calling (sub) program. LDA must be at least
!           max( 1, n ).
!           Unchanged on exit.
!
  implicit none
!
  real  alpha
      integer            incx, lda, n
      character        uplo
      complex            a( lda, * ), x( * )
      complex            zero
      parameter        ( zero = ( 0.0e+0, 0.0e+0 ) )
      complex            temp
      integer            i, info, ix, j, jx, kx
      logical            lsame
!
!  Test the input.
!
      info = 0
      if   ( .not.lsame ( uplo, 'U' ) .and.  &
               .not.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 xerbla ( 'cher  ', info )
         return
      end if
!
!  Quick return if possible.
!
      if ( n == 0 .or. alpha == real ( zero ) ) then
         return
      end if
!
!  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 one pass through the triangular part
!  of A.
!
      if ( lsame ( uplo, 'U' ) ) then
!
!  Form  A  when A is stored in upper triangle.
!
         if ( incx == 1 ) then
            do j = 1, n
               if ( x(j) /= zero ) 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) ) + real( x(j) * temp )
               else
                  a(j,j) = real( a( j,j) )
               end if
            end do
         else
            jx = kx
            do j = 1, n
               if ( x(jx) /= zero ) 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) ) + real ( x(jx) * temp )
               else
                  a(j,j) = real( a( j,j) )
               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) /= zero ) then
                  temp = alpha * conjg ( x(j) )
                  a(j,j) = real( a( j,j) ) + real( temp * x(j) )
                  do i = j + 1, n
                     a(i,j) = a(i,j) + x(i) * temp
                  end do
               else
                  a(j,j) = real( a( j,j) )
               end if
            end do
         else
            jx = kx
            do j = 1, n
               if ( x(jx) /= zero ) then
                  temp = alpha * conjg ( x(jx) )
                  a(j,j) = real( a( j,j) ) + real( temp * x(jx) )
                  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) )
               end if
               jx = jx + incx
            end do
         end if
      end if

  return
end
subroutine cher2 ( uplo, n, alpha, x, incx, y, incy, a, lda )
!
!*******************************************************************************
!
!! CHER2 performs the hermitian rank 2 operation
!
!     A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + A,
!
!  where alpha is a scalar, x and y are n element vectors and A is an n
!  by n hermitian matrix.
!
!  Parameters:
!
!  UPLO   - character.
!           On entry, UPLO specifies whether the upper or lower
!           triangular part of the array A is to be referenced as
!           follows:
!
!              UPLO = 'U' or 'u'   Only the upper triangular part of A
!                                  is to be referenced.
!
!              UPLO = 'L' or 'l'   Only the lower triangular part of A
!                                  is to be referenced.
!
!           Unchanged on exit.
!
!  N      - integer.
!           On entry, N specifies the order of the matrix A.
!           N must be at least zero.
!           Unchanged on exit.
!
!  ALPHA  - complex         .
!           On entry, ALPHA specifies the scalar alpha.
!           Unchanged on exit.
!
!  X      - complex          array of dimension at least
!           ( 1 + ( n - 1 ) * abs( incx ) ).
!           Before entry, the incremented array X must contain the n
!           element vector x.
!           Unchanged on exit.
!
!  incx   - integer.
!           On entry, incx specifies the increment for the elements of
!           X. incx must not be zero.
!           Unchanged on exit.
!
!  Y      - complex          array of dimension at least
!           ( 1 + ( n - 1 ) * abs( incy ) ).
!           Before entry, the incremented array Y must contain the n
!           element vector y.
!           Unchanged on exit.
!
!  incy   - integer.
!           On entry, incy specifies the increment for the elements of
!           Y. incy must not be zero.
!           Unchanged on exit.
!
!  A      - complex          array of DIMENSION ( LDA, n ).
!           Before entry with  UPLO = 'U' or 'u', the leading n by n
!           upper triangular part of the array A must contain the upper
!           triangular part of the hermitian matrix and the strictly
!           lower triangular part of A is not referenced. On exit, the
!           upper triangular part of the array A is overwritten by the
!           upper triangular part of the updated matrix.
!           Before entry with UPLO = 'L' or 'l', the leading n by n
!           lower triangular part of the array A must contain the lower
!           triangular part of the hermitian matrix and the strictly
!           upper triangular part of A is not referenced. On exit, the
!           lower triangular part of the array A is overwritten by the
!           lower triangular part of the updated matrix.
!           Note that the imaginary parts of the diagonal elements need
!           not be set, they are assumed to be zero, and on exit they
!           are set to zero.
!
!  LDA    - integer.
!           On entry, LDA specifies the first dimension of A as declared
!           in the calling (sub) program. LDA must be at least
!           max( 1, n ).
!           Unchanged on exit.
!
  implicit none
!
      complex            alpha
      integer            incx, incy, lda, n
      character        uplo
      complex            a( lda, * ), x( * ), y( * )
      complex            zero
      parameter        ( zero = ( 0.0e+0, 0.0e+0 ) )
      complex            temp1, temp2
      integer            i, info, ix, iy, j, jx, jy, kx, ky
      logical            lsame
!
!  Test the input.
!
      info = 0
      if   ( .not.lsame ( uplo, 'U' ) .and.  &
               .not.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 xerbla ( 'cher2 ', info )
         return
      end if
!
!  Quick return if possible.
!
      if ( ( n == 0 ).or.( alpha == zero ) ) then
         return
      end if
!
!  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 one pass through the triangular part
!  of A.
!
      if ( 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) /= zero ).or.( y(j) /= zero ) ) 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) ) + &
                              real( x(j) * temp1 + y(j) * temp2 )
               else
                  a(j,j) = real( a( j,j) )
               end if
            end do
         else
            do j = 1, n
               if ( ( x(jx) /= zero ).or.( y(jy) /= zero ) ) 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) ) + &
                              real( x(jx) * temp1 + y(jy) * temp2 )
               else
                  a(j,j) = real( a( j,j) )
               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) /= zero ).or.( y(j) /= zero ) ) then
                  temp1     = alpha * conjg ( y(j) )
                  temp2     = conjg ( alpha * x(j) )
                  a(j,j) = real( a( j,j) ) + &
                              real( x(j) * temp1 + y(j) * temp2 )
                  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) )
               end if
            end do
         else
            do j = 1, n
               if ( ( x(jx) /= zero ).or.( y(jy) /= zero ) ) then
                  temp1     = alpha * conjg ( y(jy) )
                  temp2     = conjg ( alpha * x(jx) )
                  a(j,j) = real( a( j,j) ) + &
                              real( x(jx) * temp1 + y(jy) * temp2 )
                  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) )
               end if
               jx = jx + incx
               jy = jy + incy
            end do
         end if
      end if

  return
end
subroutine cher2k ( uplo, trans, n, k, alpha, a, lda, b, ldb, beta, c, ldc )
!
!*******************************************************************************
!
!! CHER2K  performs one of the hermitian rank 2k operations
!
!     C := alpha*A*conjg( B' ) + conjg( alpha ) * b*conjg( A' ) + beta*C,
!
!  or
!
!     C := alpha*conjg( A' ) * b + conjg( alpha )*conjg( B' )*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.
!
!  Parameters:
!
!  UPLO   - character.
!           On  entry,   UPLO  specifies  whether  the  upper  or  lower
!           triangular  part  of the  array  C  is to be  referenced  as
!           follows:
!
!              UPLO = 'U' or 'u'   Only the  upper triangular part of  C
!                                  is to be referenced.
!
!              UPLO = 'L' or 'l'   Only the  lower triangular part of  C
!                                  is to be referenced.
!
!           Unchanged on exit.
!
!  TRANS  - character.
!           On entry,  TRANS  specifies the operation to be performed as
!           follows:
!
!              TRANS = 'N' or 'n'    C := alpha*A*conjg( B' )          +
!                                         conjg( alpha ) * b*conjg( A' ) +
!                                         beta*C.
!
!              TRANS = 'C' or 'c'    C := alpha*conjg( A' ) * b          +
!                                         conjg( alpha )*conjg( B' )*A +
!                                         beta*C.
!
!           Unchanged on exit.
!
!  N      - integer.
!           On entry,  N specifies the order of the matrix C.  N must be
!           at least zero.
!           Unchanged on exit.
!
!  K      - integer.
!           On entry with  TRANS = 'N' or 'n',  K  specifies  the number
!           of  columns  of the  matrices  A and B,  and on  entry  with
!           TRANS = 'C' or 'c',  K  specifies  the number of rows of the
!           matrices  A and B.  K must be at least zero.
!           Unchanged on exit.
!
!  ALPHA  - complex         .
!           On entry, ALPHA specifies the scalar alpha.
!           Unchanged on exit.
!
!  A      - complex          array of DIMENSION ( LDA, ka ), where ka is
!           k  when  TRANS = 'N' or 'n',  and is  n  otherwise.
!           Before entry with  TRANS = 'N' or 'n',  the  leading  n by k
!           part of the array  A  must contain the matrix A,  otherwise
!           the leading  k by n  part of the array  A  must contain  the
!           matrix A.
!           Unchanged on exit.
!
!  LDA    - integer.
!           On entry, LDA specifies the first dimension of A as declared
!           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
!           then  LDA must be at least  max( 1, n ), otherwise  LDA must
!           be at least  max( 1, k ).
!           Unchanged on exit.
!
!  B      - complex          array of DIMENSION ( LDB, kb ), where kb is
!           k  when  TRANS = 'N' or 'n',  and is  n  otherwise.
!           Before entry with  TRANS = 'N' or 'n',  the  leading  n by k
!           part of the array  B  must contain the matrix B,  otherwise
!           the leading  k by n  part of the array  B  must contain  the
!           matrix B.
!           Unchanged on exit.
!
!  LDB    - integer.
!           On entry, LDB specifies the first dimension of B as declared
!           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
!           then  LDB must be at least  max( 1, n ), otherwise  LDB must
!           be at least  max( 1, k ).
!           Unchanged on exit.
!
!  BETA   - real            .
!           On entry, BETA specifies the scalar beta.
!           Unchanged on exit.
!
!  C      - complex          array of DIMENSION ( LDC, n ).
!           Before entry  with  UPLO = 'U' or 'u',  the leading  n by n
!           upper triangular part of the array C must contain the upper
!           triangular part  of the  hermitian matrix and the strictly
!           lower triangular part of C is not referenced.  On exit, the
!           upper triangular part of the array  C is overwritten by the
!           upper triangular part of the updated matrix.
!           Before entry  with  UPLO = 'L' or 'l',  the leading  n by n
!           lower triangular part of the array C must contain the lower
!           triangular part  of the  hermitian matrix and the strictly
!           upper triangular part of C is not referenced.  On exit, the
!           lower triangular part of the array  C is overwritten by the
!           lower triangular part of the updated matrix.
!           Note that the imaginary parts of the diagonal elements need
!           not be set,  they are assumed to be zero,  and on exit they
!           are set to zero.
!
!  LDC    - integer.
!           On entry, LDC specifies the first dimension of C as declared
!           in  the  calling  (sub)  program.   LDC  must  be  at  least
!           max( 1, n ).
!           Unchanged on exit.
!
  implicit none
!
  character        uplo, trans
  integer            n, k, lda, ldb, ldc
  real               beta
  complex            alpha
      complex            a( lda, * ), b( ldb, * ), c( ldc, * )
      logical            lsame
      logical            upper
      integer            i, info, j, l, nrowa
      complex            temp1, temp2
      real               one
      parameter        ( one  = 1.0e+0 )
      complex            zero
      parameter        ( zero = ( 0.0e+0, 0.0e+0 ) )
!
!  Test the input.
!
      if ( lsame ( trans, 'N' ) ) then
         nrowa = n
      else
         nrowa = k
      end if
      upper = lsame ( uplo, 'U' )

      info = 0
      if (      ( .not.upper               ) .and.  &
               ( .not.lsame ( uplo , 'L' ) )      ) then
         info = 1
      else if ( ( .not.lsame ( trans, 'N' ) ) .and.  &
               ( .not.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 xerbla ( 'cher2k', info )
         return
      end if
!
!  Quick return if possible.
!
      if ( ( n == 0 ).or. &
          ( ( ( alpha == zero ).or.( k == 0 ) ) .and. ( beta == one ) ) ) &
         return
!
!  And when  alpha == zero.
!
      if ( alpha == zero ) then
         if ( upper ) then
            if ( beta == real( zero ) ) then
               do j = 1, n
                  c(1:j,j) = zero
               end do
            else
               do j = 1, n
                  c(1:j-1,j) = beta * c(1:j-1,j)
                  c(j,j) = beta * real( c( j,j) )
               end do
            end if
         else
            if ( beta == real( zero ) ) then
               do j = 1, n
                  c(j:n,j) = zero
               end do
            else
               do j = 1, n
                  c(j,j) = beta * real( c( j,j) )
                  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 ( lsame ( trans, 'N' ) ) then
!
!  Form  C := alpha*A*conjg( B' ) + conjg( alpha ) * b*conjg( A' ) + C.
!
         if ( upper ) then
            do j = 1, n
               if ( beta == real( 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) )
               else
                  c(j,j) = real( c( j,j) )
               end if
               do l = 1, k
                  if ( ( a(j,l) /= zero ).or. ( b(j,l) /= zero ) ) 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) )         + &
                                 real( a(j,l) * temp1 + b(j,l) * temp2   )
                  end if
               end do
            end do
         else
            do j = 1, n
               if ( beta == real( zero ) ) then
                  do i = j, n
                     c(i,j) = zero
                  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) )
               else
                  c(j,j) = real( c( j,j) )
               end if
               do l = 1, k
                  if ( ( a(j,l) /= zero ).or. ( b(j,l) /= zero ) ) 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) )         + &
                                 real( a(j,l) * temp1 + b(j,l) * temp2   )
                  end if
               end do
            end do
         end if
      else
!
!  Form  C := alpha*conjg( A' ) * b + conjg( alpha )*conjg( B' )*A + C.
!
         if ( upper ) then
            do j = 1, n
               do i = 1, j
                  temp1 = zero
                  temp2 = zero
                  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( zero ) ) then
                        c(j,j) = real(        alpha   * temp1 + &
                                          conjg ( alpha ) * temp2   )
                     else
                        c(j,j) = beta * real( c( j,j) )         + &
                                    real(        alpha   * temp1 + &
                                          conjg ( alpha ) * temp2   )
                     end if
                  else
                     if ( beta == real( zero ) ) 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 = zero
                  temp2 = zero
                  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( zero ) ) then
                        c(j,j) = real(        alpha   * temp1 + &
                                          conjg ( alpha ) * temp2   )
                     else
                        c(j,j) = beta * real( c( j,j) )         + &
                                    real(        alpha   * temp1 + &
                                          conjg ( alpha ) * temp2   )
                     end if
                  else
                     if ( beta == real( zero ) ) 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 cherk ( uplo, trans, n, k, alpha, a, lda, beta, c, ldc )
!
!*******************************************************************************
!
!! CHERK performs one of the hermitian rank k operations
!
!     C := alpha*A*conjg( A' ) + beta*C,
!
!  or
!
!     C := alpha*conjg( A' )*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.
!
!  Parameters:
!
!  UPLO   - character.
!           On  entry,   UPLO  specifies  whether  the  upper  or  lower
!           triangular  part  of the  array  C  is to be  referenced  as
!           follows:
!
!              UPLO = 'U' or 'u'   Only the  upper triangular part of  C
!                                  is to be referenced.
!
!              UPLO = 'L' or 'l'   Only the  lower triangular part of  C
!                                  is to be referenced.
!
!           Unchanged on exit.
!
!  TRANS  - character.
!           On entry,  TRANS  specifies the operation to be performed as
!           follows:
!
!              TRANS = 'N' or 'n'   C := alpha*A*conjg( A' ) + beta*C.
!
!              TRANS = 'C' or 'c'   C := alpha*conjg( A' )*A + beta*C.
!
!           Unchanged on exit.
!
!  N      - integer.
!           On entry,  N specifies the order of the matrix C.  N must be
!           at least zero.
!           Unchanged on exit.
!
!  K      - integer.
!           On entry with  TRANS = 'N' or 'n',  K  specifies  the number
!           of  columns   of  the   matrix A,   and  on   entry   with
!           TRANS = 'C' or 'c',  K  specifies  the number of rows of the
!           matrix A.  K must be at least zero.
!           Unchanged on exit.
!
!  ALPHA  - real            .
!           On entry, ALPHA specifies the scalar alpha.
!           Unchanged on exit.
!
!  A      - complex          array of DIMENSION ( LDA, ka ), where ka is
!           k  when  TRANS = 'N' or 'n',  and is  n  otherwise.
!           Before entry with  TRANS = 'N' or 'n',  the  leading  n by k
!           part of the array  A  must contain the matrix A,  otherwise
!           the leading  k by n  part of the array  A  must contain  the
!           matrix A.
!           Unchanged on exit.
!
!  LDA    - integer.
!           On entry, LDA specifies the first dimension of A as declared
!           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
!           then  LDA must be at least  max( 1, n ), otherwise  LDA must
!           be at least  max( 1, k ).
!           Unchanged on exit.
!
!  BETA   - real            .
!           On entry, BETA specifies the scalar beta.
!           Unchanged on exit.
!
!  C      - complex          array of DIMENSION ( LDC, n ).
!           Before entry  with  UPLO = 'U' or 'u',  the leading  n by n
!           upper triangular part of the array C must contain the upper
!           triangular part  of the  hermitian matrix and the strictly
!           lower triangular part of C is not referenced.  On exit, the
!           upper triangular part of the array  C is overwritten by the
!           upper triangular part of the updated matrix.
!           Before entry  with  UPLO = 'L' or 'l',  the leading  n by n
!           lower triangular part of the array C must contain the lower
!           triangular part  of the  hermitian matrix and the strictly
!           upper triangular part of C is not referenced.  On exit, the
!           lower triangular part of the array  C is overwritten by the
!           lower triangular part of the updated matrix.
!           Note that the imaginary parts of the diagonal elements need
!           not be set,  they are assumed to be zero,  and on exit they
!           are set to zero.
!
!  LDC    - integer.
!           On entry, LDC specifies the first dimension of C as declared
!           in  the  calling  (sub)  program.   LDC  must  be  at  least
!           max( 1, n ).
!           Unchanged on exit.
!
!
  implicit none
!
      character        uplo, trans
      integer            n, k, lda, ldc
      real               alpha, beta
      complex            a( lda, * ), c( ldc, * )
      logical            lsame
      logical            upper
      integer            i, info, j, l, nrowa
      real               rtemp
      complex            temp
      real               one ,         zero
      parameter        ( one = 1.0e+0, zero = 0.0e+0 )
!
!  Test the input.
!
      if ( lsame ( trans, 'N' ) ) then
         nrowa = n
      else
         nrowa = k
      end if

      upper = lsame ( uplo, 'U' )

      info = 0
      if (      ( .not.upper               ) .and.  &
               ( .not.lsame ( uplo , 'L' ) )      ) then
         info = 1
      else if ( ( .not.lsame ( trans, 'N' ) ) .and.  &
               ( .not.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 xerbla ( 'cherk ', info )
         return
      end if
!
!  Quick return if possible.
!
      if ( ( n == 0 ).or. &
          ( ( ( alpha == zero ).or.( k == 0 ) ) .and. ( beta == one ) ) ) &
         return
!
!  And when  alpha == 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) )
               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) )
                  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 ( lsame ( trans, 'N' ) ) then
!
!  Form  C := alpha*A*conjg( A' ) + 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) )
               else
                  c(j,j) = real( c( j,j) )
               end if
               do l = 1, k
                  if ( a(j,l) /= cmplx ( zero ) ) 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)      ) + &
                                 real( temp * a(i,l) )
                  end if
               end do
            end do
         else
            do j = 1, n
               if ( beta == zero ) then
                  c(j:n,j) = zero
               else if ( beta/= one ) then
                  c(j,j) = beta * real( c( j,j) )
                  c(j+1:n,j) = beta * c(j+1:n,j)
               else
                  c(j,j) = real( c( j,j) )
               end if
               do l = 1, k
                  if ( a(j,l) /= cmplx ( zero ) ) then
                     temp = alpha * conjg ( a(j,l) )
                     c(j,j) = real( c( j,j)      )   + &
                                 real( temp * a(j,l) )
                     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*conjg( A' )*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) )
               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) )
               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 chpmv ( uplo, n, alpha, ap, x, incx, beta, y, incy )
!
!*******************************************************************************
!
!! CHPMV  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.
!
!  Parameters:
!
!  UPLO   - character.
!           On entry, UPLO specifies whether the upper or lower
!           triangular part of the matrix A is supplied in the packed
!           array AP as follows:
!
!              UPLO = 'U' or 'u'   The upper triangular part of A is
!                                  supplied in AP.
!
!              UPLO = 'L' or 'l'   The lower triangular part of A is
!                                  supplied in AP.
!
!           Unchanged on exit.
!
!  N      - integer.
!           On entry, N specifies the order of the matrix A.
!           N must be at least zero.
!           Unchanged on exit.
!
!  ALPHA  - complex         .
!           On entry, ALPHA specifies the scalar alpha.
!           Unchanged on exit.
!
!  AP     - complex          array of DIMENSION at least
!           ( ( n*( n + 1 ) )/2 ).
!           Before entry with UPLO = 'U' or 'u', the array AP must
!           contain the upper triangular part of the hermitian matrix
!           packed sequentially, column by column, so that ap( 1 )
!           contains a( 1, 1 ), ap( 2 ) and ap( 3 ) contain a( 1, 2 )
!           and a( 2, 2 ) respectively, and so on.
!           Before entry with UPLO = 'L' or 'l', the array AP must
!           contain the lower triangular part of the hermitian matrix
!           packed sequentially, column by column, so that ap( 1 )
!           contains a( 1, 1 ), ap( 2 ) and ap( 3 ) contain a( 2, 1 )
!           and a( 3, 1 ) respectively, and so on.
!           Note that the imaginary parts of the diagonal elements need
!           not be set and are assumed to be zero.
!           Unchanged on exit.
!
!  X      - complex          array of dimension at least
!           ( 1 + ( n - 1 ) * abs( incx ) ).
!           Before entry, the incremented array X must contain the n
!           element vector x.
!           Unchanged on exit.
!
!  incx   - integer.
!           On entry, incx specifies the increment for the elements of
!           X. incx must not be zero.
!           Unchanged on exit.
!
!  BETA   - complex         .
!           On entry, BETA specifies the scalar beta. When BETA is
!           supplied as zero then Y need not be set on input.
!           Unchanged on exit.
!
!  Y      - complex          array of dimension at least
!           ( 1 + ( n - 1 ) * abs( incy ) ).
!           Before entry, the incremented array Y must contain the n
!           element vector y. On exit, Y is overwritten by the updated
!           vector y.
!
!  incy   - integer.
!           On entry, incy specifies the increment for the elements of
!           Y. incy must not be zero.
!           Unchanged on exit.
!
!
  implicit none
!
      complex            alpha, beta
      integer            incx, incy, n
      character        uplo
      complex            ap( * ), x( * ), y( * )
      complex            one
      parameter        ( one  = ( 1.0e+0, 0.0e+0 ) )
      complex            zero
      parameter        ( zero = ( 0.0e+0, 0.0e+0 ) )
      complex            temp1, temp2
      integer            i, info, ix, iy, j, jx, jy, k, kk, kx, ky
      logical            lsame
!
!  Test the input.
!
      info = 0
      if   ( .not.lsame ( uplo, 'U' ) .and.  &
               .not.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 xerbla ( 'chpmv ', info )
         return
      end if
!
!     Quick return if possible.
!
      if ( ( n == 0 ).or.( ( alpha == zero ) .and. ( beta == one ) ) ) &
         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 one pass through AP.
!
!  First form  y := beta*y.
!
      if ( beta/= one ) then
         if ( incy == 1 ) then
            if ( beta == zero ) then
               y(1:n) = zero
            else
               y(1:n) = beta * y(1:n)
            end if
         else
            iy = ky
            if ( beta == zero ) then
               do i = 1, n
                  y(iy) = zero
                  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 == zero ) then
         return
      end if
      kk = 1
      if ( 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 = zero
               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 ) ) &
                               + alpha * temp2
               kk = kk + j
            end do
         else
            jx = kx
            jy = ky
            do j = 1, n
               temp1 = alpha * x(jx)
               temp2 = zero
               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 ) ) &
                                 + 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 = zero
               y(j) = y(j) + temp1 * real( ap( kk ) )
               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   = zero
               y(jy) = y(jy) + temp1 * real( ap( kk ) )
               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 chpr ( uplo, n, alpha, x, incx, ap )
!
!*******************************************************************************
!
!! CHPR    performs the hermitian rank 1 operation
!
!     A := alpha*x*conjg( x' ) + 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.
!
!  Parameters:
!
!  UPLO   - character.
!           On entry, UPLO specifies whether the upper or lower
!           triangular part of the matrix A is supplied in the packed
!           array AP as follows:
!
!              UPLO = 'U' or 'u'   The upper triangular part of A is
!                                  supplied in AP.
!
!              UPLO = 'L' or 'l'   The lower triangular part of A is
!                                  supplied in AP.
!
!           Unchanged on exit.
!
!  N      - integer.
!           On entry, N specifies the order of the matrix A.
!           N must be at least zero.
!           Unchanged on exit.
!
!  ALPHA  - real            .
!           On entry, ALPHA specifies the scalar alpha.
!           Unchanged on exit.
!
!  X      - complex          array of dimension at least
!           ( 1 + ( n - 1 ) * abs( incx ) ).
!           Before entry, the incremented array X must contain the n
!           element vector x.
!           Unchanged on exit.
!
!  incx   - integer.
!           On entry, incx specifies the increment for the elements of
!           X. incx must not be zero.
!           Unchanged on exit.
!
!  AP     - complex          array of DIMENSION at least
!           ( ( n*( n + 1 ) )/2 ).
!           Before entry with  UPLO = 'U' or 'u', the array AP must
!           contain the upper triangular part of the hermitian matrix
!           packed sequentially, column by column, so that ap( 1 )
!           contains a( 1, 1 ), ap( 2 ) and ap( 3 ) contain a( 1, 2 )
!           and a( 2, 2 ) respectively, and so on. On exit, the array
!           AP is overwritten by the upper triangular part of the
!           updated matrix.
!           Before entry with UPLO = 'L' or 'l', the array AP must
!           contain the lower triangular part of the hermitian matrix
!           packed sequentially, column by column, so that ap( 1 )
!           contains a( 1, 1 ), ap( 2 ) and ap( 3 ) contain a( 2, 1 )
!           and a( 3, 1 ) respectively, and so on. On exit, the array
!           AP is overwritten by the lower triangular part of the
!           updated matrix.
!           Note that the imaginary parts of the diagonal elements need
!           not be set, they are assumed to be zero, and on exit they
!           are set to zero.
!
  implicit none
!
      real               alpha
      integer            incx, n
      character        uplo
      complex            ap( * ), x( * )
      complex            zero
      parameter        ( zero = ( 0.0e+0, 0.0e+0 ) )
      complex            temp
      integer            i, info, ix, j, jx, k, kk, kx
      logical            lsame
!
!  Test the input.
!
      info = 0
      if   ( .not.lsame ( uplo, 'U' ) .and.  &
               .not.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 xerbla ( 'chpr  ', info )
         return
      end if
!
!  Quick return if possible.
!
      if ( ( n == 0 ).or.( alpha == real( zero ) ) ) &
         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 one pass through AP.
!
      kk = 1
      if ( lsame ( uplo, 'U' ) ) then
!
!  Form  A  when upper triangle is stored in AP.
!
         if ( incx == 1 ) then
            do j = 1, n
               if ( x(j) /= zero ) 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 ) ) &
                                     + real( x(j) * temp )
               else
                  ap( kk + j - 1 ) = real( ap( kk + j - 1 ) )
               end if
               kk = kk + j
            end do
         else
            jx = kx
            do j = 1, n
               if ( x(jx) /= zero ) 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 ) ) &
                                     + real( x(jx) * temp )
               else
                  ap( kk + j - 1 ) = real( ap( kk + j - 1 ) )
               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) /= zero ) then
                  temp = alpha * conjg ( x(j) )
                  ap( kk ) = real( ap( kk ) ) + real( temp * x(j) )
                  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 ) )
               end if
               kk = kk + n - j + 1
            end do
         else
            jx = kx
            do j = 1, n
               if ( x(jx) /= zero ) then
                  temp = alpha * conjg ( x(jx) )
                  ap( kk ) = real( ap( kk ) ) + real( temp * x(jx) )
                  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 ) )
               end if
               jx = jx + incx
               kk = kk + n - j + 1
            end do
         end if
      end if

  return
end
subroutine chpr2 ( uplo, n, alpha, x, incx, y, incy, ap )
!
!*******************************************************************************
!
!! CHPR2 performs the hermitian rank 2 operation
!
!     A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + 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.
!
!  Parameters:
!
!  UPLO   - character.
!           On entry, UPLO specifies whether the upper or lower
!           triangular part of the matrix A is supplied in the packed
!           array AP as follows:
!
!              UPLO = 'U' or 'u'   The upper triangular part of A is
!                                  supplied in AP.
!
!              UPLO = 'L' or 'l'   The lower triangular part of A is
!                                  supplied in AP.
!
!           Unchanged on exit.
!
!  N      - integer.
!           On entry, N specifies the order of the matrix A.
!           N must be at least zero.
!           Unchanged on exit.
!
!  ALPHA  - complex         .
!           On entry, ALPHA specifies the scalar alpha.
!           Unchanged on exit.
!
!  X      - complex          array of dimension at least
!           ( 1 + ( n - 1 ) * abs( incx ) ).
!           Before entry, the incremented array X must contain the n
!           element vector x.
!           Unchanged on exit.
!
!  incx   - integer.
!           On entry, incx specifies the increment for the elements of
!           X. incx must not be zero.
!           Unchanged on exit.
!
!  Y      - complex          array of dimension at least
!           ( 1 + ( n - 1 ) * abs( incy ) ).
!           Before entry, the incremented array Y must contain the n
!           element vector y.
!           Unchanged on exit.
!
!  incy   - integer.
!           On entry, incy specifies the increment for the elements of
!           Y. incy must not be zero.
!           Unchanged on exit.
!
!  AP     - complex          array of DIMENSION at least
!           ( ( n*( n + 1 ) )/2 ).
!           Before entry with  UPLO = 'U' or 'u', the array AP must
!           contain the upper triangular part of the hermitian matrix
!           packed sequentially, column by column, so that ap( 1 )
!           contains a( 1, 1 ), ap( 2 ) and ap( 3 ) contain a( 1, 2 )
!           and a( 2, 2 ) respectively, and so on. On exit, the array
!           AP is overwritten by the upper triangular part of the
!           updated matrix.
!           Before entry with UPLO = 'L' or 'l', the array AP must
!           contain the lower triangular part of the hermitian matrix
!           packed sequentially, column by column, so that ap( 1 )
!           contains a( 1, 1 ), ap( 2 ) and ap( 3 ) contain a( 2, 1 )
!           and a( 3, 1 ) respectively, and so on. On exit, the array
!           AP is overwritten by the lower triangular part of the
!           updated matrix.
!           Note that the imaginary parts of the diagonal elements need
!           not be set, they are assumed to be zero, and on exit they
!           are set to zero.
!
  implicit none
!
      complex            alpha
      integer            incx, incy, n
      character        uplo
      complex            ap( * ), x( * ), y( * )
      complex            zero
      parameter        ( zero = ( 0.0e+0, 0.0e+0 ) )
      complex            temp1, temp2
      integer            i, info, ix, iy, j, jx, jy, k, kk, kx, ky
      logical            lsame
!
!  Test the input.
!
      info = 0
      if   ( .not.lsame ( uplo, 'U' ) .and.  &
               .not.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 xerbla ( 'chpr2 ', info )
         return
      end if
!
!  Quick return if possible.
!
      if ( ( n == 0 ).or.( alpha == zero ) ) &
         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 one pass through AP.
!
      kk = 1
      if ( 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) /= zero ).or.( y(j) /= zero ) ) 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 ) ) + &
                                     real( x(j) * temp1 + y(j) * temp2 )
               else
                  ap( kk + j - 1 ) = real( ap( kk + j - 1 ) )
               end if
               kk = kk + j
            end do
         else
            do j = 1, n
               if ( ( x(jx) /= zero ).or.( y(jy) /= zero ) ) 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 ) ) + &
                                     real( x(jx) * temp1 + &
                                           y(jy) * temp2 )
               else
                  ap( kk + j - 1 ) = real( ap( kk + j - 1 ) )
               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) /= zero ).or.( y(j) /= zero ) ) then
                  temp1   = alpha * conjg ( y(j) )
                  temp2   = conjg ( alpha * x(j) )
                  ap( kk ) = real( ap( kk ) ) + &
                             real( x(j) * temp1 + y(j) * temp2 )
                  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 ) )
               end if
               kk = kk + n - j + 1
            end do
         else
            do j = 1, n
               if ( ( x(jx) /= zero ).or.( y(jy) /= zero ) ) then
                  temp1    = alpha * conjg ( y(jy) )
                  temp2    = conjg ( alpha * x(jx) )
                  ap( kk ) = real( ap( kk ) ) + &
                             real( x(jx) * temp1 + y(jy) * temp2 )
                  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 ) )
               end if
               jx = jx + incx
               jy = jy + incy
               kk = kk + n - j + 1
            end do
         end if
      end if

  return
end
subroutine crotg ( ca, cb, c, s )
!
!*******************************************************************************
!
!! CROTG ???
!
  implicit none
!
  complex alpha
  complex c
  complex ca
  complex cb
  real norm
  complex s
  real scale
!
  if ( cabs(ca) == 0.0 ) then
    c = 0.0
    s = ( 1.0, 0.0 )
    ca = cb
  else
    scale = cabs(ca) + cabs(cb)
    norm = scale * sqrt((cabs(ca/scale))**2 + (cabs(cb/scale))**2)
    alpha = ca /cabs(ca)
    c = cabs(ca) / norm
    s = alpha * conjg(cb) / norm
    ca = alpha * norm
  end if

  return
end
subroutine cscal ( n, ca, cx, incx )
!
!*******************************************************************************
!
!! CSCAL scales a vector by a constant.
!
!
!     jack dongarra, linpack,  3/11/78.
!     modified 3/93 to return if incx  <=  0.
!     modified 12/3/93, array(1) declarations changed to array(*)
!
!  Reference:
!
!    Lawson, Hanson, Kincaid, Krogh,
!    Basic Linear Algebra Subprograms for Fortran Usage,
!    Algorithm 539,
!    ACM Transactions on Mathematical Software,
!    Volume 5, Number 3, September 1979, pages 308-323.
!
  implicit none
!
  complex ca
  complex cx(*)
  integer i
  integer incx
  integer n
  integer nincx
!
  if ( n <= 0 .or. incx <= 0 ) then
    return
  end if

  if ( incx == 1 ) then

    cx(1:n) = ca * cx(1:n)

  else

    nincx = n*incx

    cx(1:nincx:incx) = ca * cx(1:nincx:incx)

  end if

  return
end
function csign ( z1, z2 )
!
!*******************************************************************************
!
!! CSIGN is a sort of complex transfer-of-sign function.
!
!
!  Modified:
!
!    30 May 2002
!
!  Parameters:
!
  implicit none
!
  complex csign
  complex z1
  complex z2
!
  csign = cabs ( z1 ) * ( z2 / cabs ( z2 ) )

  return
end
function csign1 ( z1, z2 )
!
!*******************************************************************************
!
!! CSIGN1 is a sort of complex transfer-of-sign function.
!
!
!  Modified:
!
!    23 May 2002
!
!  Parameters:
!
  implicit none
!
  real cabs1
  complex csign1
  complex z1
  complex z2
!
  csign1 = cabs1 ( z1 ) * ( z2 / cabs1 ( z2 ) )

  return
end
subroutine csrot ( n, cx, incx, cy, incy, c, s )
!
!*******************************************************************************
!
!! CSROT applies a plane rotation,
!
!
!  The cosine C and sine S are real
!     and the vectors cx and cy are complex.
!
!
!     jack dongarra, linpack, 3/11/78.
!
!  Reference:
!
!    Lawson, Hanson, Kincaid, Krogh,
!    Basic Linear Algebra Subprograms for Fortran Usage,
!    Algorithm 539,
!    ACM Transactions on Mathematical Software,
!    Volume 5, Number 3, September 1979, pages 308-323.
!
  implicit none
!
  real c
  complex ctemp
  complex cx(*)
  complex cy(*)
  integer i
  integer incx
  integer incy
  integer ix
  integer iy
  integer n
  real s
!
  if ( n <= 0 ) then
    return
  end if

  if (incx == 1 .and. incy == 1 ) then

    do i = 1,n
      ctemp = c*cx(i) + s*cy(i)
      cy(i) = c*cy(i) - s*cx(i)
      cx(i) = ctemp
    end do

  else

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( -n + 1 ) * incx + 1
    end if

    if ( incy >= 0 ) then
      iy = 1
    else
      iy = ( -n + 1 ) * incy + 1
    end if

    do i = 1, n
      ctemp = c * cx(ix) + s * cy(iy)
      cy(iy) = c * cy(iy) - s * cx(ix)
      cx(ix) = ctemp
      ix = ix + incx
      iy = iy + incy
    end do

  end if

  return
end
subroutine csscal ( n, sa, cx, incx )
!
!*******************************************************************************
!
!! CSSCAL scales a complex vector by a real constant.
!
!
!     jack dongarra, linpack, 3/11/78.
!     modified 3/93 to return if incx  <=  0.
!     modified 12/3/93, array(1) declarations changed to array(*)
!
!  Reference:
!
!    Lawson, Hanson, Kincaid, Krogh,
!    Basic Linear Algebra Subprograms for Fortran Usage,
!    Algorithm 539,
!    ACM Transactions on Mathematical Software,
!    Volume 5, Number 3, September 1979, pages 308-323.
!
  implicit none
!
  complex cx(*)
  integer i
  integer incx
  integer n
  integer nincx
  real sa
!
  if ( n <= 0 .or. incx <= 0 ) then
    return
  end if

  if ( incx == 1 ) then

    cx(1:n) = sa * cx(1:n)

  else

    nincx = n*incx

    cx(1:nincx:incx) = sa * cx(1:nincx:incx)

  end if

  return
end
subroutine cswap ( n, cx, incx, cy, incy )
!
!*******************************************************************************
!
!! CSWAP interchanges two vectors.
!
!
!     jack dongarra, linpack, 3/11/78.
!     modified 12/3/93, array(1) declarations changed to array(*)
!
!  Reference:
!
!    Lawson, Hanson, Kincaid, Krogh,
!    Basic Linear Algebra Subprograms for Fortran Usage,
!    Algorithm 539,
!    ACM Transactions on Mathematical Software,
!    Volume 5, Number 3, September 1979, pages 308-323.
!
  implicit none
!
  complex ctemp
  complex cx(*)
  complex cy(*)
  integer i
  integer incx
  integer incy
  integer ix
  integer iy
  integer n
!
  if ( n <= 0 ) then
    return
  end if

  if ( incx == 1 .and. incy == 1 ) then

    do i = 1,n
      ctemp = cx(i)
      cx(i) = cy(i)
      cy(i) = ctemp
    end do

  else

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( -n + 1 ) * incx + 1
    end if

    if ( incy >= 0 ) then
      iy = 1
    else
      iy = ( -n + 1 ) * incy + 1
    end if

    do i = 1, n
      ctemp = cx(ix)
      cx(ix) = cy(iy)
      cy(iy) = ctemp
      ix = ix + incx
      iy = iy + incy
    end do

  end if

  return
end
subroutine csymm ( side, uplo, m, n, alpha, a, lda, b, ldb, beta, c, ldc )
!
!*******************************************************************************
!
!! CSYMM 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.
!
!  Parameters:
!
!  SIDE   - character.
!           On entry,  SIDE  specifies whether  the  symmetric matrix A
!           appears on the  left or right  in the  operation as follows:
!
!              SIDE = 'L' or 'l'   C := alpha*A * b + beta*C,
!
!              SIDE = 'R' or 'r'   C := alpha * b*A + beta*C,
!
!           Unchanged on exit.
!
!  UPLO   - character.
!           On  entry,   UPLO  specifies  whether  the  upper  or  lower
!           triangular  part  of  the  symmetric  matrix A  is  to  be
!           referenced as follows:
!
!              UPLO = 'U' or 'u'   Only the upper triangular part of the
!                                  symmetric matrix is to be referenced.
!
!              UPLO = 'L' or 'l'   Only the lower triangular part of the
!                                  symmetric matrix is to be referenced.
!
!           Unchanged on exit.
!
!  M      - integer.
!           On entry,  M  specifies the number of rows of the matrix C.
!           M  must be at least zero.
!           Unchanged on exit.
!
!  N      - integer.
!           On entry, N specifies the number of columns of the matrix C.
!           N  must be at least zero.
!           Unchanged on exit.
!
!  ALPHA  - complex         .
!           On entry, ALPHA specifies the scalar alpha.
!           Unchanged on exit.
!
!  A      - complex          array of DIMENSION ( LDA, ka ), where ka is
!           m  when  SIDE = 'L' or 'l'  and is n  otherwise.
!           Before entry  with  SIDE = 'L' or 'l',  the  m by m  part of
!           the array  A  must contain the  symmetric matrix,  such that
!           when  UPLO = 'U' or 'u', the leading m by m upper triangular
!           part of the array  A  must contain the upper triangular part
!           of the  symmetric matrix and the  strictly  lower triangular
!           part of  A  is not referenced,  and when  UPLO = 'L' or 'l',
!           the leading  m by m  lower triangular part  of the  array  A
!           must  contain  the  lower triangular part  of the  symmetric
!           matrix and the  strictly upper triangular part of  A  is not
!           referenced.
!           Before entry  with  SIDE = 'R' or 'r',  the  n by n  part of
!           the array  A  must contain the  symmetric matrix,  such that
!           when  UPLO = 'U' or 'u', the leading n by n upper triangular
!           part of the array  A  must contain the upper triangular part
!           of the  symmetric matrix and the  strictly  lower triangular
!           part of  A  is not referenced,  and when  UPLO = 'L' or 'l',
!           the leading  n by n  lower triangular part  of the  array  A
!           must  contain  the  lower triangular part  of the  symmetric
!           matrix and the  strictly upper triangular part of  A  is not
!           referenced.
!           Unchanged on exit.
!
!  LDA    - integer.
!           On entry, LDA specifies the first dimension of A as declared
!           in the  calling (sub) program. When  SIDE = 'L' or 'l'  then
!           LDA must be at least  max( 1, m ), otherwise  LDA must be at
!           least max( 1, n ).
!           Unchanged on exit.
!
!  B      - complex          array of DIMENSION ( LDB, n ).
!           Before entry, the leading  m by n part of the array  B  must
!           contain the matrix B.
!           Unchanged on exit.
!
!  LDB    - integer.
!           On entry, LDB specifies the first dimension of B as declared
!           in  the  calling  (sub)  program.   LDB  must  be  at  least
!           max( 1, m ).
!           Unchanged on exit.
!
!  BETA   - complex         .
!           On entry,  BETA  specifies the scalar  beta.  When  BETA  is
!           supplied as zero then C need not be set on input.
!           Unchanged on exit.
!
!  C      - complex          array of DIMENSION ( LDC, n ).
!           Before entry, the leading  m by n  part of the array  C must
!           contain the matrix C,  except when  beta  is zero, in which
!           case C need not be set on entry.
!           On exit, the array  C  is overwritten by the  m by n updated
!           matrix.
!
!  LDC    - integer.
!           On entry, LDC specifies the first dimension of C as declared
!           in  the  calling  (sub)  program.   LDC  must  be  at  least
!           max( 1, m ).
!           Unchanged on exit.
!
  implicit none
!
!     .. Scalar Arguments ..
      character        side, uplo
      integer            m, n, lda, ldb, ldc
      complex            alpha, beta
!     .. Array Arguments ..
      complex            a( lda, * ), b( ldb, * ), c( ldc, * )
!     ..
      logical            lsame
      logical            upper
      integer            i, info, j, k, nrowa
      complex            temp1, temp2
      complex            one
      parameter        ( one  = ( 1.0e+0, 0.0e+0 ) )
      complex            zero
      parameter        ( zero = ( 0.0e+0, 0.0e+0 ) )
!
!  Set NROWA as the number of rows of A.
!
      if ( lsame ( side, 'L' ) ) then
         nrowa = m
      else
         nrowa = n
      end if

      upper = lsame ( uplo, 'U' )
!
!  Test the input.
!
      info = 0
      if (      ( .not.lsame ( side, 'L' ) ) .and.  &
               ( .not.lsame ( side, 'R' ) )      ) then
         info = 1
      else if ( ( .not.upper              ) .and.  &
               ( .not.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 xerbla ( 'csymm ', info )
         return
      end if
!
!  Quick return if possible.
!
      if ( ( m == 0 ).or.( n == 0 ) .or.  &
          ( ( alpha == zero ) .and. ( beta == one ) ) ) &
         return
!
!  And when  alpha == zero.
!
      if ( alpha == zero ) then
         if ( beta == zero ) then
            c(1:m,1:n) = zero
         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 ( 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 = zero
                  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 == zero ) 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 = zero
                  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 == zero ) 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.
!
         do j = 1, n

            temp1 = alpha * a(j,j)
            if ( beta == zero ) 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

      end if

  return
end
subroutine csyr2k ( uplo, trans, n, k, alpha, a, lda, b, ldb, beta, c, ldc )
!
!*******************************************************************************
!
!! CSYR2K  performs one of the symmetric rank 2k operations
!
!     C := alpha*A * b' + alpha*B*A' + beta*C,
!
!  or
!
!     C := alpha*A' * b + alpha*B'*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.
!
!  Parameters:
!
!  UPLO   - character.
!           On  entry,   UPLO  specifies  whether  the  upper  or  lower
!           triangular  part  of the  array  C  is to be  referenced  as
!           follows:
!
!              UPLO = 'U' or 'u'   Only the  upper triangular part of  C
!                                  is to be referenced.
!
!              UPLO = 'L' or 'l'   Only the  lower triangular part of  C
!                                  is to be referenced.
!
!           Unchanged on exit.
!
!  TRANS  - character.
!           On entry,  TRANS  specifies the operation to be performed as
!           follows:
!
!              TRANS = 'N' or 'n'    C := alpha*A * b' + alpha*B*A' +
!                                         beta*C.
!
!              TRANS = 'T' or 't'    C := alpha*A' * b + alpha*B'*A +
!                                         beta*C.
!
!           Unchanged on exit.
!
!  N      - integer.
!           On entry,  N specifies the order of the matrix C.  N must be
!           at least zero.
!           Unchanged on exit.
!
!  K      - integer.
!           On entry with  TRANS = 'N' or 'n',  K  specifies  the number
!           of  columns  of the  matrices  A and B,  and on  entry  with
!           TRANS = 'T' or 't',  K  specifies  the number of rows of the
!           matrices  A and B.  K must be at least zero.
!           Unchanged on exit.
!
!  ALPHA  - complex         .
!           On entry, ALPHA specifies the scalar alpha.
!           Unchanged on exit.
!
!  A      - complex          array of DIMENSION ( LDA, ka ), where ka is
!           k  when  TRANS = 'N' or 'n',  and is  n  otherwise.
!           Before entry with  TRANS = 'N' or 'n',  the  leading  n by k
!           part of the array  A  must contain the matrix A,  otherwise
!           the leading  k by n  part of the array  A  must contain  the
!           matrix A.
!           Unchanged on exit.
!
!  LDA    - integer.
!           On entry, LDA specifies the first dimension of A as declared
!           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
!           then  LDA must be at least  max( 1, n ), otherwise  LDA must
!           be at least  max( 1, k ).
!           Unchanged on exit.
!
!  B      - complex          array of DIMENSION ( LDB, kb ), where kb is
!           k  when  TRANS = 'N' or 'n',  and is  n  otherwise.
!           Before entry with  TRANS = 'N' or 'n',  the  leading  n by k
!           part of the array  B  must contain the matrix B,  otherwise
!           the leading  k by n  part of the array  B  must contain  the
!           matrix B.
!           Unchanged on exit.
!
!  LDB    - integer.
!           On entry, LDB specifies the first dimension of B as declared
!           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
!           then  LDB must be at least  max( 1, n ), otherwise  LDB must
!           be at least  max( 1, k ).
!           Unchanged on exit.
!
!  BETA   - complex         .
!           On entry, BETA specifies the scalar beta.
!           Unchanged on exit.
!
!  C      - complex          array of DIMENSION ( LDC, n ).
!           Before entry  with  UPLO = 'U' or 'u',  the leading  n by n
!           upper triangular part of the array C must contain the upper
!           triangular part  of the  symmetric matrix and the strictly
!           lower triangular part of C is not referenced.  On exit, the
!           upper triangular part of the array  C is overwritten by the
!           upper triangular part of the updated matrix.
!           Before entry  with  UPLO = 'L' or 'l',  the leading  n by n
!           lower triangular part of the array C must contain the lower
!           triangular part  of the  symmetric matrix and the strictly
!           upper triangular part of C is not referenced.  On exit, the
!           lower triangular part of the array  C is overwritten by the
!           lower triangular part of the updated matrix.
!
!  LDC    - integer.
!           On entry, LDC specifies the first dimension of C as declared
!           in  the  calling  (sub)  program.   LDC  must  be  at  least
!           max( 1, n ).
!           Unchanged on exit.
!
  implicit none
!
      character        uplo, trans
      integer            n, k, lda, ldb, ldc
      complex            alpha, beta
      complex            a( lda, * ), b( ldb, * ), c( ldc, * )
      logical            lsame
      logical            upper
      integer            i, info, j, l, nrowa
      complex            temp1, temp2
      complex            one
      parameter        ( one  = ( 1.0e+0, 0.0e+0 ) )
      complex            zero
      parameter        ( zero = ( 0.0e+0, 0.0e+0 ) )
!
!  Test the input.
!
      if ( lsame ( trans, 'N' ) ) then
         nrowa = n
      else
         nrowa = k
      end if
      upper = lsame ( uplo, 'U' )

      info = 0
      if (      ( .not.upper               ) .and.  &
               ( .not.lsame ( uplo , 'L' ) )      ) then
         info = 1
      else if ( ( .not.lsame ( trans, 'N' ) ) .and.  &
               ( .not.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 xerbla ( 'csyr2k', info )
         return
      end if
!
!  Quick return if possible.
!
      if ( ( n == 0 ).or. &
          ( ( ( alpha == zero ).or.( k == 0 ) ) .and. ( beta == one ) ) ) &
         return
!
!  And when  alpha == zero.
!
      if ( alpha == zero ) then
         if ( upper ) then
            if ( beta == zero ) then
               do j = 1, n
                  c(1:j,j) = zero
               end do
            else
               do j = 1, n
                  c(1:j,j) = beta * c(1:j,j)
               end do
            end if
         else
            if ( beta == zero ) then
               do j = 1, n
                  c(j:n,j) = zero
               end do
            else
               do j = 1, n
                  c(j:n,j) = beta * c(j:n,j)
               end do
            end if
         end if
         return
      end if
!
!  Start the operations.
!
      if ( lsame ( trans, 'N' ) ) then
!
!  Form  C := alpha*A * b' + alpha*B*A' + C.
!
         if ( upper ) then

            do j = 1, n

               if ( beta == zero ) then
                  c(1:j,j) = zero
               else if ( beta/= one ) then
                  c(1:j,j) = beta * c(1:j,j)
               end if

               do l = 1, k

                  if ( ( a(j,l) /= zero ).or. ( b(j,l) /= zero ) ) 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 == zero ) then
                  c(j:n,j) = zero
               else if ( beta/= one ) then
                  c(j:n,j) = beta * c(j:n,j)
               end if
               do l = 1, k
                  if ( ( a(j,l) /= zero ).or. ( b(j,l) /= zero ) ) 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' * b + alpha*B'*A + C.
!
         if ( upper ) then

            do j = 1, n
               do i = 1, j
                  temp1 = zero
                  temp2 = zero
                  do l = 1, k
                     temp1 = temp1 + a(l,i) * b(l,j)
                     temp2 = temp2 + b(l,i) * a(l,j)
                  end do
                  if ( beta == zero ) 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 = zero
                  temp2 = zero
                  do l = 1, k
                     temp1 = temp1 + a(l,i) * b(l,j)
                     temp2 = temp2 + b(l,i) * a(l,j)
                  end do
                  if ( beta == zero ) 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 csyrk ( uplo, trans, n, k, alpha, a, lda, beta, c, ldc )
!
!*******************************************************************************
!
!! CSYRK  performs one of the symmetric rank k operations
!
!     C := alpha*A*A' + beta*C,
!
!  or
!
!     C := alpha*A'*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.
!
!  Parameters:
!
!  UPLO   - character.
!           On  entry,   UPLO  specifies  whether  the  upper  or  lower
!           triangular  part  of the  array  C  is to be  referenced  as
!           follows:
!
!              UPLO = 'U' or 'u'   Only the  upper triangular part of  C
!                                  is to be referenced.
!
!              UPLO = 'L' or 'l'   Only the  lower triangular part of  C
!                                  is to be referenced.
!
!           Unchanged on exit.
!
!  TRANS  - character.
!           On entry,  TRANS  specifies the operation to be performed as
!           follows:
!
!              TRANS = 'N' or 'n'   C := alpha*A*A' + beta*C.
!
!              TRANS = 'T' or 't'   C := alpha*A'*A + beta*C.
!
!           Unchanged on exit.
!
!  N      - integer.
!           On entry,  N specifies the order of the matrix C.  N must be
!           at least zero.
!           Unchanged on exit.
!
!  K      - integer.
!           On entry with  TRANS = 'N' or 'n',  K  specifies  the number
!           of  columns   of  the   matrix A,   and  on   entry   with
!           TRANS = 'T' or 't',  K  specifies  the number of rows of the
!           matrix A.  K must be at least zero.
!           Unchanged on exit.
!
!  ALPHA  - complex         .
!           On entry, ALPHA specifies the scalar alpha.
!           Unchanged on exit.
!
!  A      - complex          array of DIMENSION ( LDA, ka ), where ka is
!           k  when  TRANS = 'N' or 'n',  and is  n  otherwise.
!           Before entry with  TRANS = 'N' or 'n',  the  leading  n by k
!           part of the array  A  must contain the matrix A,  otherwise
!           the leading  k by n  part of the array  A  must contain  the
!           matrix A.
!           Unchanged on exit.
!
!  LDA    - integer.
!           On entry, LDA specifies the first dimension of A as declared
!           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
!           then  LDA must be at least  max( 1, n ), otherwise  LDA must
!           be at least  max( 1, k ).
!           Unchanged on exit.
!
!  BETA   - complex         .
!           On entry, BETA specifies the scalar beta.
!           Unchanged on exit.
!
!  C      - complex          array of DIMENSION ( LDC, n ).
!           Before entry  with  UPLO = 'U' or 'u',  the leading  n by n
!           upper triangular part of the array C must contain the upper
!           triangular part  of the  symmetric matrix and the strictly
!           lower triangular part of C is not referenced.  On exit, the
!           upper triangular part of the array  C is overwritten by the
!           upper triangular part of the updated matrix.
!           Before entry  with  UPLO = 'L' or 'l',  the leading  n by n
!           lower triangular part of the array C must contain the lower
!           triangular part  of the  symmetric matrix and the strictly
!           upper triangular part of C is not referenced.  On exit, the
!           lower triangular part of the array  C is overwritten by the
!           lower triangular part of the updated matrix.
!
!  LDC    - integer.
!           On entry, LDC specifies the first dimension of C as declared
!           in  the  calling  (sub)  program.   LDC  must  be  at  least
!           max( 1, n ).
!           Unchanged on exit.
!
  implicit none
!
!     .. Scalar Arguments ..
      character        uplo, trans
      integer            n, k, lda, ldc
      complex            alpha, beta
!     .. Array Arguments ..
      complex            a( lda, * ), c( ldc, * )
!     ..
      logical            lsame
      logical            upper
      integer            i, info, j, l, nrowa
      complex            temp
      complex            one
      parameter        ( one  = ( 1.0e+0, 0.0e+0 ) )
      complex            zero
      parameter        ( zero = ( 0.0e+0, 0.0e+0 ) )
!
!  Test the input.
!
      if ( lsame ( trans, 'N' ) ) then
         nrowa = n
      else
         nrowa = k
      end if
      upper = lsame ( uplo, 'U' )

      info = 0
      if (      ( .not.upper               ) .and.  &
               ( .not.lsame ( uplo , 'L' ) )      ) then
         info = 1
      else if ( ( .not.lsame ( trans, 'N' ) ) .and.  &
               ( .not.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 xerbla ( 'csyrk ', info )
         return
      end if
!
!  Quick return if possible.
!
      if ( ( n == 0 ).or. &
          ( ( ( alpha == zero ).or.( k == 0 ) ) .and. ( beta == one ) ) ) &
         return
!
!  And when ALPHA == zero.
!
      if ( alpha == zero ) then
         if ( upper ) then
            if ( beta == zero ) then
               do j = 1, n
                  c(1:j,j) = zero
               end do
            else
               do j = 1, n
                  c(1:j,j) = beta * c(1:j,j)
               end do
            end if
         else
            if ( beta == zero ) then
               do j = 1, n
                  c(j:n,j) = zero
               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 ( lsame ( trans, 'N' ) ) then
!
!  Form C := alpha*A*A' + beta*C.
!
         if ( upper ) then
            do j = 1, n
               if ( beta == zero ) then
                  c(1:j,j) = zero
               else if ( beta/= one ) then
                  c(1:j,j) = beta * c(1:j,j)
               end if
               do l = 1, k
                  if ( a(j,l) /= zero ) 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 == zero ) then
                  c(j:n,j) = zero
               else if ( beta/= one ) then
                  do i = j, n
                     c(i,j) = beta * c(i,j)
                  end do
               end if
               do l = 1, k
                  if ( a(j,l) /= zero ) 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'*A + beta*C.
!
         if ( upper ) then
            do j = 1, n
               do i = 1, j
                  temp = zero
                  do l = 1, k
                     temp = temp + 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
         else
            do j = 1, n
               do i = j, n
                  temp = zero
                  do l = 1, k
                     temp = temp + 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 ctbmv ( uplo, trans, diag, n, k, a, lda, x, incx )
!
!*******************************************************************************
!
!! CTBMV  performs one of the matrix-vector operations
!
!     x := A*x,   or   x := A'*x,   or   x := conjg( A' )*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.
!
!  Parameters:
!
!  UPLO   - character.
!           On entry, UPLO specifies whether the matrix is an upper or
!           lower triangular matrix as follows:
!
!              UPLO = 'U' or 'u'   A is an upper triangular matrix.
!
!              UPLO = 'L' or 'l'   A is a lower triangular matrix.
!
!           Unchanged on exit.
!
!  TRANS  - character.
!           On entry, TRANS specifies the operation to be performed as
!           follows:
!
!              TRANS = 'N' or 'n'   x := A*x.
!
!              TRANS = 'T' or 't'   x := A'*x.
!
!              TRANS = 'C' or 'c'   x := conjg( A' )*x.
!
!           Unchanged on exit.
!
!  DIAG   - character.
!           On entry, DIAG specifies whether or not A is unit
!           triangular as follows:
!
!              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
!
!              DIAG = 'N' or 'n'   A is not assumed to be unit
!                                  triangular.
!
!           Unchanged on exit.
!
!  N      - integer.
!           On entry, N specifies the order of the matrix A.
!           N must be at least zero.
!           Unchanged on exit.
!
!  K      - integer.
!           On entry with UPLO = 'U' or 'u', K specifies the number of
!           super-diagonals of the matrix A.
!           On entry with UPLO = 'L' or 'l', K specifies the number of
!           sub-diagonals of the matrix A.
!           K must satisfy  0  <=  K.
!           Unchanged on exit.
!
!  A      - complex          array of DIMENSION ( LDA, n ).
!           Before entry with UPLO = 'U' or 'u', the leading ( k + 1 )
!           by n part of the array A must contain the upper triangular
!           band part of the matrix of coefficients, supplied column by
!           column, with the leading diagonal of the matrix in row
!           ( k + 1 ) of the array, the first super-diagonal starting at
!           position 2 in row k, and so on. The top left k by k triangle
!           of the array A is not referenced.
!           The following program segment will transfer an upper
!           triangular band matrix from conventional full matrix storage
!           to band storage:
!
!                 do j = 1, n
!                    m = k + 1 - J
!                    do i = max ( 1, j - K ), J
!                       a( M + I,j) = matrix(i,j)
!                    end do
!                 end do
!
!           Before entry with UPLO = 'L' or 'l', the leading ( k + 1 )
!           by n part of the array A must contain the lower triangular
!           band part of the matrix of coefficients, supplied column by
!           column, with the leading diagonal of the matrix in row 1 of
!           the array, the first sub-diagonal starting at position 1 in
!           row 2, and so on. The bottom right k by k triangle of the
!           array A is not referenced.
!           The following program segment will transfer a lower
!           triangular band matrix from conventional full matrix storage
!           to band storage:
!
!                 do j = 1, n
!                    m = 1 - J
!                    do i = J, min ( N, j + K )
!                       a( M + I,j) = matrix(i,j)
!                    end do
!                 end do
!
!           Note that when DIAG = 'U' or 'u' the elements of the array A
!           corresponding to the diagonal elements of the matrix are not
!           referenced, but are assumed to be unity.
!           Unchanged on exit.
!
!  LDA    - integer.
!           On entry, LDA specifies the first dimension of A as declared
!           in the calling (sub) program. LDA must be at least
!           ( k + 1 ).
!           Unchanged on exit.
!
!  X      - complex          array of dimension at least
!           ( 1 + ( n - 1 ) * abs( incx ) ).
!           Before entry, the incremented array X must contain the n
!           element vector x. On exit, X is overwritten with the
!           tranformed vector x.
!
!  incx   - integer.
!           On entry, incx specifies the increment for the elements of
!           X. incx must not be zero.
!           Unchanged on exit.
!
  implicit none
!
      integer            incx, k, lda, n
      character        diag, trans, uplo
      complex            a( lda, * ), x( * )
      complex            zero
      parameter        ( zero = ( 0.0e+0, 0.0e+0 ) )
      complex            temp
      integer            i, info, ix, j, jx, kplus1, kx, l
      logical            noconj, nounit
      logical            lsame
!
!  Test the input.
!
      info = 0
      if   ( .not.lsame ( uplo , 'U' ) .and.  &
               .not.lsame ( uplo , 'L' )      ) then
         info = 1
      else if ( .not.lsame ( trans, 'N' ) .and.  &
               .not.lsame ( trans, 'T' ) .and.  &
               .not.lsame ( trans, 'C' )      ) then
         info = 2
      else if ( .not.lsame ( diag , 'U' ) .and.  &
               .not.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 xerbla ( 'ctbmv ', info )
         return
      end if
!
!  Quick return if possible.
!
      if ( n == 0 ) then
         return
      end if

      noconj = lsame ( trans, 'T' )
      nounit = 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 one pass through A.
!
      if ( lsame ( trans, 'N' ) ) then
!
!  Form x := A*x.
!
         if ( lsame ( uplo, 'U' ) ) then
            kplus1 = k + 1
            if ( incx == 1 ) then
               do j = 1, n
                  if ( x(j) /= zero ) 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 ) then
                        x(j) = x(j) * a( kplus1,j)
                     end if
                  end if
               end do
            else
               jx = kx
               do j = 1, n
                  if ( x(jx) /= zero ) 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 ) then
                        x(jx) = x(jx) * a( kplus1,j)
                     end if
                  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) /= zero ) 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 ) then
                        x(j) = x(j) * a( 1,j)
                     end if
                  end if
               end do
            else
               kx = kx + ( n - 1 ) * incx
               jx = kx
               do j = n, 1, -1
                  if ( x(jx) /= zero ) 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 ) then
                     kx = kx - incx
                  end if
               end do
            end if
         end if
      else
!
!  Form  x := A'*x  or  x := conjg( A' )*x.
!
         if ( 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 ) then
                        temp = temp * a( kplus1,j)
                     end if
                     do i = j - 1, max( 1, j - k ), -1
                        temp = temp + a( l + i,j) * x(i)
                     end do
                  else
                     if ( nounit ) then
                        temp = temp * conjg ( a( kplus1,j) )
                     end if
                     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 ) then
                        temp = temp * a( kplus1,j)
                     end if
                     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 ) then
                        temp = temp*conjg ( a( kplus1,j) )
                     end if
                     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 ) then
                        temp = temp * a( 1,j)
                     end if
                     do i = j + 1, min( n, j + k )
                        temp = temp + a( l + i,j) * x(i)
                     end do
                  else
                     if ( nounit ) then
                        temp = temp * conjg ( a( 1,j) )
                     end if
                     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 ctbsv ( uplo, trans, diag, n, k, a, lda, x, incx )
!
!*******************************************************************************
!
!! CTBSV  solves one of the systems of equations
!
!     A*x = b,   or   A'*x = b,   or   conjg( A' )*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.
!
!  Parameters:
!
!  UPLO   - character.
!           On entry, UPLO specifies whether the matrix is an upper or
!           lower triangular matrix as follows:
!
!              UPLO = 'U' or 'u'   A is an upper triangular matrix.
!
!              UPLO = 'L' or 'l'   A is a lower triangular matrix.
!
!           Unchanged on exit.
!
!  TRANS  - character.
!           On entry, TRANS specifies the equations to be solved as
!           follows:
!
!              TRANS = 'N' or 'n'   A*x = b.
!
!              TRANS = 'T' or 't'   A'*x = b.
!
!              TRANS = 'C' or 'c'   conjg( A' )*x = b.
!
!           Unchanged on exit.
!
!  DIAG   - character.
!           On entry, DIAG specifies whether or not A is unit
!           triangular as follows:
!
!              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
!
!              DIAG = 'N' or 'n'   A is not assumed to be unit
!                                  triangular.
!
!           Unchanged on exit.
!
!  N      - integer.
!           On entry, N specifies the order of the matrix A.
!           N must be at least zero.
!           Unchanged on exit.
!
!  K      - integer.
!           On entry with UPLO = 'U' or 'u', K specifies the number of
!           super-diagonals of the matrix A.
!           On entry with UPLO = 'L' or 'l', K specifies the number of
!           sub-diagonals of the matrix A.
!           K must satisfy  0  <=  K.
!           Unchanged on exit.
!
!  A      - complex          array of DIMENSION ( LDA, n ).
!           Before entry with UPLO = 'U' or 'u', the leading ( k + 1 )
!           by n part of the array A must contain the upper triangular
!           band part of the matrix of coefficients, supplied column by
!           column, with the leading diagonal of the matrix in row
!           ( k + 1 ) of the array, the first super-diagonal starting at
!           position 2 in row k, and so on. The top left k by k triangle
!           of the array A is not referenced.
!           The following program segment will transfer an upper
!           triangular band matrix from conventional full matrix storage
!           to band storage:
!
!                 do j = 1, n
!                    m = k + 1 - J
!                    do i = max ( 1, j - K ), J
!                       a( M + I,j) = matrix(i,j)
!                    end do
!                 end do
!
!           Before entry with UPLO = 'L' or 'l', the leading ( k + 1 )
!           by n part of the array A must contain the lower triangular
!           band part of the matrix of coefficients, supplied column by
!           column, with the leading diagonal of the matrix in row 1 of
!           the array, the first sub-diagonal starting at position 1 in
!           row 2, and so on. The bottom right k by k triangle of the
!           array A is not referenced.
!           The following program segment will transfer a lower
!           triangular band matrix from conventional full matrix storage
!           to band storage:
!
!                 do j = 1, n
!                    m = 1 - J
!                    do i = J, min ( N, j + K )
!                       a( M + I,j) = matrix(i,j)
!                    end do
!                 end do
!
!           Note that when DIAG = 'U' or 'u' the elements of the array A
!           corresponding to the diagonal elements of the matrix are not
!           referenced, but are assumed to be unity.
!           Unchanged on exit.
!
!  LDA    - integer.
!           On entry, LDA specifies the first dimension of A as declared
!           in the calling (sub) program. LDA must be at least
!           ( k + 1 ).
!           Unchanged on exit.
!
!  X      - complex          array of dimension at least
!           ( 1 + ( n - 1 ) * abs( incx ) ).
!           Before entry, the incremented array X must contain the n
!           element right-hand side vector b. On exit, X is overwritten
!           with the solution vector x.
!
!  incx   - integer.
!           On entry, incx specifies the increment for the elements of
!           X. incx must not be zero.
!           Unchanged on exit.
!
  implicit none
!
      integer            incx, k, lda, n
      character        diag, trans, uplo
      complex            a( lda, * ), x( * )
      complex            zero
      parameter        ( zero = ( 0.0e+0, 0.0e+0 ) )
      complex            temp
      integer            i, info, ix, j, jx, kplus1, kx, l
      logical            noconj, nounit
      logical            lsame
!
!  Test the input.
!
      info = 0
      if   ( .not.lsame ( uplo , 'U' ) .and.  &
               .not.lsame ( uplo , 'L' )      ) then
         info = 1
      else if ( .not.lsame ( trans, 'N' ) .and.  &
               .not.lsame ( trans, 'T' ) .and.  &
               .not.lsame ( trans, 'C' )      ) then
         info = 2
      else if ( .not.lsame ( diag , 'U' ) .and.  &
               .not.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 xerbla ( 'ctbsv ', info )
         return
      end if
!
!  Quick return if possible.
!
      if ( n == 0 ) then
         return
      end if

      noconj = lsame ( trans, 'T' )
      nounit = 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 one pass through A.
!
      if ( lsame ( trans, 'N' ) ) then
!
!  Form  x := inv( A )*x.
!
         if ( lsame ( uplo, 'U' ) ) then
            kplus1 = k + 1
            if ( incx == 1 ) then
               do j = n, 1, -1
                  if ( x(j) /= zero ) then
                     l = kplus1 - j
                     if ( nounit ) then
                        x(j) = x(j)/a( kplus1,j)
                     end if
                     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) /= zero ) 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) /= zero ) then
                     l = 1 - j
                     if ( nounit ) then
                        x(j) = x(j) / a( 1,j)
                     end if
                     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) /= zero ) then
                     ix = kx
                     l  = 1  - j
                     if ( nounit ) then
                        x(jx) = x(jx)/a( 1,j)
                     end if
                     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' )*x  or  x := inv( conjg( A') )*x.
!
         if ( 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 ) then
                        temp = temp/a( kplus1,j)
                     end if
                  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 ( noconj ) then
                     do  i = max ( 1, j - k ), j - 1
                        temp = temp - a( l + i,j) * x(ix)
                        ix = ix + incx
                     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(ix)
                        ix = ix + incx
                     end do
                     if ( nounit ) then
                        temp = temp/conjg ( a( kplus1,j) )
                     end if
                  end if
                  x(jx) = temp
                  jx = jx + incx
                  if ( j > k ) then
                     kx = kx + incx
                  end if
               end do
            end if
         else
            if ( incx == 1 ) then
               do j = n, 1, -1
                  temp = x(j)
                  l    = 1      - j
                  if ( noconj ) then
                     do i = min ( n, j + k ), j + 1, -1
                        temp = temp - a( l + i,j) * x(i)
                     end do
                     if ( nounit ) &
                        temp = temp/a( 1,j)
                  else
                     do i = min ( n, j + k ), j + 1, -1
                        temp = temp - conjg ( a( l + i,j) ) * x(i)
                     end do
                     if ( nounit ) &
                        temp = temp/conjg ( a( 1,j) )
                  end if
                  x(j) = temp
               end do
            else
               kx = kx + ( n - 1 ) * incx
               jx = kx
               do j = n, 1, -1
                  temp = x(jx)
                  ix = kx
                  l = 1 - j
                  if ( noconj ) then
                     do i = min ( n, j + k ), j + 1, -1
                        temp = temp - a( l + i,j) * x(ix)
                        ix = ix - incx
                     end do
                     if ( nounit ) &
                        temp = temp/a( 1,j)
                  else
                     do i = min ( n, j + k ), j + 1, -1
                        temp = temp - conjg ( a( l + i,j) ) * x(ix)
                        ix = ix - incx
                     end do
                     if ( nounit ) then
                        temp = temp/conjg ( a( 1,j) )
                     end if
                  end if
                  x(jx) = temp
                  jx = jx - incx
                  if ( ( n - j )>=k ) &
                     kx = kx - incx
               end do
            end if
         end if
      end if

  return
end
subroutine ctpmv ( uplo, trans, diag, n, ap, x, incx )
!
!*******************************************************************************
!
!! CTPMV  performs one of the matrix-vector operations
!
!     x := A*x,   or   x := A'*x,   or   x := conjg( A' )*x,
!
!  where x is an n element vector and  A is an n by n unit, or non-unit,
!  upper or lower triangular matrix, supplied in packed form.
!
!  Parameters:
!
!  UPLO   - character.
!           On entry, UPLO specifies whether the matrix is an upper or
!           lower triangular matrix as follows:
!
!              UPLO = 'U' or 'u'   A is an upper triangular matrix.
!
!              UPLO = 'L' or 'l'   A is a lower triangular matrix.
!
!           Unchanged on exit.
!
!  TRANS  - character.
!           On entry, TRANS specifies the operation to be performed as
!           follows:
!
!              TRANS = 'N' or 'n'   x := A*x.
!
!              TRANS = 'T' or 't'   x := A'*x.
!
!              TRANS = 'C' or 'c'   x := conjg( A' )*x.
!
!           Unchanged on exit.
!
!  DIAG   - character.
!           On entry, DIAG specifies whether or not A is unit
!           triangular as follows:
!
!              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
!
!              DIAG = 'N' or 'n'   A is not assumed to be unit
!                                  triangular.
!
!           Unchanged on exit.
!
!  N      - integer.
!           On entry, N specifies the order of the matrix A.
!           N must be at least zero.
!           Unchanged on exit.
!
!  AP     - complex          array of DIMENSION at least
!           ( ( n*( n + 1 ) )/2 ).
!           Before entry with  UPLO = 'U' or 'u', the array AP must
!           contain the upper triangular matrix packed sequentially,
!           column by column, so that ap( 1 ) contains a( 1, 1 ),
!           ap( 2 ) and ap( 3 ) contain a( 1, 2 ) and a( 2, 2 )
!           respectively, and so on.
!           Before entry with UPLO = 'L' or 'l', the array AP must
!           contain the lower triangular matrix packed sequentially,
!           column by column, so that ap( 1 ) contains a( 1, 1 ),
!           ap( 2 ) and ap( 3 ) contain a( 2, 1 ) and a( 3, 1 )
!           respectively, and so on.
!           Note that when  DIAG = 'U' or 'u', the diagonal elements of
!           A are not referenced, but are assumed to be unity.
!           Unchanged on exit.
!
!  X      - complex          array of dimension at least
!           ( 1 + ( n - 1 ) * abs( incx ) ).
!           Before entry, the incremented array X must contain the n
!           element vector x. On exit, X is overwritten with the
!           tranformed vector x.
!
!  incx   - integer.
!           On entry, incx specifies the increment for the elements of
!           X. incx must not be zero.
!           Unchanged on exit.
!
  implicit none
!
      integer            incx, n
      character        diag, trans, uplo
      complex            ap( * ), x( * )
      complex            zero
      parameter        ( zero = ( 0.0e+0, 0.0e+0 ) )
      complex            temp
      integer            i, info, ix, j, jx, k, kk, kx
      logical            noconj, nounit
      logical            lsame
!
!  Test the input.
!
      info = 0
      if   ( .not.lsame ( uplo , 'U' ) .and.  &
               .not.lsame ( uplo , 'L' )      ) then
         info = 1
      else if ( .not.lsame ( trans, 'N' ) .and.  &
               .not.lsame ( trans, 'T' ) .and.  &
               .not.lsame ( trans, 'C' )      ) then
         info = 2
      else if ( .not.lsame ( diag , 'U' ) .and.  &
               .not.lsame ( diag , 'N' )      ) then
         info = 3
      else if ( n<0 ) then
         info = 4
      else if ( incx == 0 ) then
         info = 7
      end if
      if ( info /= 0 ) then
         call xerbla ( 'ctpmv ', info )
         return
      end if
!
!  Quick return if possible.
!
      if ( n == 0 ) then
         return
      end if

      noconj = lsame ( trans, 'T' )
      nounit = 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 AP are
!  accessed sequentially with one pass through AP.
!
      if ( lsame ( trans, 'N' ) ) then
!
!  Form  x:= A*x.
!
         if ( lsame ( uplo, 'U' ) ) then
            kk = 1
            if ( incx == 1 ) then
               do j = 1, n
                  if ( x(j) /= zero ) then
                     temp = x(j)
                     k    = kk
                     do i = 1, j - 1
                        x(i) = x(i) + temp * ap(k)
                        k      = k      + 1
                     end do
                     if ( nounit ) then
                        x(j) = x(j) * ap( kk + j - 1 )
                     end if
                  end if
                  kk = kk + j
               end do
            else
               jx = kx
               do j = 1, n
                  if ( x(jx) /= zero ) then
                     temp = x(jx)
                     ix = kx
                     do k = kk, kk + j - 2
                        x(ix) = x(ix) + temp * ap(k)
                        ix = ix + incx
                     end do
                     if ( nounit ) then
                        x(jx) = x(jx) * ap( kk + j - 1 )
                     end if
                  end if
                  jx = jx + incx
                  kk = kk + j
               end do
            end if
         else
            kk = ( n*( n + 1 ) )/2
            if ( incx == 1 ) then
               do j = n, 1, -1
                  if ( x(j) /= zero ) then
                     temp = x(j)
                     k    = kk
                     do i = n, j + 1, -1
                        x(i) = x(i) + temp * ap(k)
                        k      = k      - 1
                     end do
                     if ( nounit ) &
                        x(j) = x(j) * ap( kk - n + j )
                  end if
                  kk = kk - ( n - j + 1 )
               end do
            else
               kx = kx + ( n - 1 ) * incx
               jx = kx
               do j = n, 1, -1
                  if ( x(jx) /= zero ) then
                     temp = x(jx)
                     ix = kx
                     do k = kk, kk - ( n - ( j + 1 ) ), -1
                        x(ix) = x(ix) + temp * ap(k)
                        ix = ix - incx
                     end do
                     if ( nounit ) &
                        x(jx) = x(jx) * ap( kk - n + j )
                  end if
                  jx = jx - incx
                  kk = kk - ( n - j + 1 )
               end do
            end if
         end if
      else
!
!  Form  x := A'*x  or  x := conjg( A' )*x.
!
         if ( lsame ( uplo, 'U' ) ) then
            kk = ( n*( n + 1 ) )/2
            if ( incx == 1 ) then
               do j = n, 1, -1
                  temp = x(j)
                  k    = kk     - 1
                  if ( noconj ) then
                     if ( nounit ) then
                        temp = temp * ap( kk )
                     end if
                     do i = j - 1, 1, -1
                        temp = temp + ap(k) * x(i)
                        k    = k    - 1
                     end do
                  else
                     if ( nounit ) then
                        temp = temp*conjg ( ap( kk ) )
                     end if
                     do i = j - 1, 1, -1
                        temp = temp + conjg ( ap(k) ) * x(i)
                        k    = k    - 1
                     end do
                  end if
                  x(j) = temp
                  kk     = kk   - j
               end do
            else
               jx = kx + ( n - 1 ) * incx
               do j = n, 1, -1
                  temp = x(jx)
                  ix = jx
                  if ( noconj ) then
                     if ( nounit ) then
                        temp = temp * ap( kk )
                     end if
                     do k = kk - 1, kk - j + 1, -1
                        ix = ix - incx
                        temp = temp + ap(k) * x(ix)
                     end do
                  else
                     if ( nounit ) &
                        temp = temp*conjg ( ap( kk ) )
                     do k = kk - 1, kk - j + 1, -1
                        ix = ix - incx
                        temp = temp + conjg ( ap(k) ) * x(ix)
                     end do
                  end if
                  x(jx) = temp
                  jx = jx - incx
                  kk      = kk   - j
               end do
            end if
         else
            kk = 1
            if ( incx == 1 ) then
               do j = 1, n
                  temp = x(j)
                  k    = kk     + 1
                  if ( noconj ) then
                     if ( nounit ) &
                        temp = temp * ap( kk )
                     do i = j + 1, n
                        temp = temp + ap(k) * x(i)
                        k    = k    + 1
                     end do
                  else
                     if ( nounit ) &
                        temp = temp * conjg ( ap( kk ) )
                     do i = j + 1, n
                        temp = temp + conjg ( ap(k) ) * x(i)
                        k    = k    + 1
                     end do
                  end if
                  x(j) = temp
                  kk     = kk   + ( n - j + 1 )
               end do
            else
               jx = kx
               do j = 1, n
                  temp = x(jx)
                  ix = jx
                  if ( noconj ) then
                     if ( nounit ) &
                        temp = temp * ap( kk )
                     do k = kk + 1, kk + n - j
                        ix = ix + incx
                        temp = temp + ap(k) * x(ix)
                     end do
                  else
                     if ( nounit ) &
                        temp = temp*conjg ( ap( kk ) )
                     do k = kk + 1, kk + n - j
                        ix = ix + incx
                        temp = temp + conjg ( ap(k) ) * x(ix)
                     end do
                  end if
                  x(jx) = temp
                  jx = jx + incx
                  kk      = kk   + ( n - j + 1 )
               end do
            end if
         end if
      end if

  return
end
subroutine ctpsv ( uplo, trans, diag, n, ap, x, incx )
!
!*******************************************************************************
!
!! CTPSV  solves one of the systems of equations
!
!     A*x = b,   or   A'*x = b,   or   conjg( A' )*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 matrix, supplied in packed form.
!
!  No test for singularity or near-singularity is included in this
!  routine. Such tests must be performed before calling this routine.
!
!  Parameters:
!
!  UPLO   - character.
!           On entry, UPLO specifies whether the matrix is an upper or
!           lower triangular matrix as follows:
!
!              UPLO = 'U' or 'u'   A is an upper triangular matrix.
!
!              UPLO = 'L' or 'l'   A is a lower triangular matrix.
!
!           Unchanged on exit.
!
!  TRANS  - character.
!           On entry, TRANS specifies the equations to be solved as
!           follows:
!
!              TRANS = 'N' or 'n'   A*x = b.
!
!              TRANS = 'T' or 't'   A'*x = b.
!
!              TRANS = 'C' or 'c'   conjg( A' )*x = b.
!
!           Unchanged on exit.
!
!  DIAG   - character.
!           On entry, DIAG specifies whether or not A is unit
!           triangular as follows:
!
!              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
!
!              DIAG = 'N' or 'n'   A is not assumed to be unit
!                                  triangular.
!
!           Unchanged on exit.
!
!  N      - integer.
!           On entry, N specifies the order of the matrix A.
!           N must be at least zero.
!           Unchanged on exit.
!
!  AP     - complex          array of DIMENSION at least
!           ( ( n*( n + 1 ) )/2 ).
!           Before entry with  UPLO = 'U' or 'u', the array AP must
!           contain the upper triangular matrix packed sequentially,
!           column by column, so that ap( 1 ) contains a( 1, 1 ),
!           ap( 2 ) and ap( 3 ) contain a( 1, 2 ) and a( 2, 2 )
!           respectively, and so on.
!           Before entry with UPLO = 'L' or 'l', the array AP must
!           contain the lower triangular matrix packed sequentially,
!           column by column, so that ap( 1 ) contains a( 1, 1 ),
!           ap( 2 ) and ap( 3 ) contain a( 2, 1 ) and a( 3, 1 )
!           respectively, and so on.
!           Note that when  DIAG = 'U' or 'u', the diagonal elements of
!           A are not referenced, but are assumed to be unity.
!           Unchanged on exit.
!
!  X      - complex          array of dimension at least
!           ( 1 + ( n - 1 ) * abs( incx ) ).
!           Before entry, the incremented array X must contain the n
!           element right-hand side vector b. On exit, X is overwritten
!           with the solution vector x.
!
!  incx   - integer.
!           On entry, incx specifies the increment for the elements of
!           X. incx must not be zero.
!           Unchanged on exit.
!
  implicit none
!
      integer            incx, n
      character        diag, trans, uplo
      complex            ap( * ), x( * )
      complex            zero
      parameter        ( zero = ( 0.0e+0, 0.0e+0 ) )
      complex            temp
      integer            i, info, ix, j, jx, k, kk, kx
      logical            noconj, nounit
      logical            lsame
!
!  Test the input.
!
      info = 0
      if   ( .not.lsame ( uplo , 'U' ) .and.  &
               .not.lsame ( uplo , 'L' )      ) then
         info = 1
      else if ( .not.lsame ( trans, 'N' ) .and.  &
               .not.lsame ( trans, 'T' ) .and.  &
               .not.lsame ( trans, 'C' )      ) then
         info = 2
      else if ( .not.lsame ( diag , 'U' ) .and.  &
               .not.lsame ( diag , 'N' )      ) then
         info = 3
      else if ( n<0 ) then
         info = 4
      else if ( incx == 0 ) then
         info = 7
      end if
      if ( info/=0 ) then
         call xerbla ( 'ctpsv ', info )
         return
      end if
!
!  Quick return if possible.
!
      if ( n == 0 ) &
         return

      noconj = lsame ( trans, 'T' )
      nounit = 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 AP are
!  accessed sequentially with one pass through AP.
!
      if ( lsame ( trans, 'N' ) ) then
!
!  Form  x := inv( A )*x.
!
         if ( lsame ( uplo, 'U' ) ) then
            kk = ( n*( n + 1 ) )/2
            if ( incx == 1 ) then
               do j = n, 1, -1
                  if ( x(j) /= zero ) then
                     if ( nounit ) &
                        x(j) = x(j)/ap( kk )
                     temp = x(j)
                     k    = kk     - 1
                     do i = j - 1, 1, -1
                        x(i) = x(i) - temp * ap(k)
                        k      = k      - 1
                     end do
                  end if
                  kk = kk - j
               end do
            else
               jx = kx + ( n - 1 ) * incx
               do j = n, 1, -1
                  if ( x(jx) /= zero ) then
                     if ( nounit ) &
                        x(jx) = x(jx)/ap( kk )
                     temp = x(jx)
                     ix = jx
                     do k = kk - 1, kk - j + 1, -1
                        ix = ix - incx
                        x(ix) = x(ix) - temp * ap(k)
                     end do
                  end if
                  jx = jx - incx
                  kk = kk - j
               end do
            end if
         else
            kk = 1
            if ( incx == 1 ) then
               do j = 1, n
                  if ( x(j) /= zero ) then
                     if ( nounit ) &
                        x(j) = x(j)/ap( kk )
                     temp = x(j)
                     k    = kk     + 1
                     do i = j + 1, n
                        x(i) = x(i) - temp * ap(k)
                        k      = k      + 1
                     end do
                  end if
                  kk = kk + ( n - j + 1 )
               end do
            else
               jx = kx
               do j = 1, n
                  if ( x(jx) /= zero ) then
                     if ( nounit ) &
                        x(jx) = x(jx)/ap( kk )
                     temp = x(jx)
                     ix = jx
                     do k = kk + 1, kk + n - j
                        ix = ix + incx
                        x(ix) = x(ix) - temp * ap(k)
                     end do
                  end if
                  jx = jx + incx
                  kk = kk + ( n - j + 1 )
               end do
            end if
         end if
      else
!
!  Form  x := inv( A' )*x  or  x := inv( conjg( A' ) )*x.
!
         if ( lsame ( uplo, 'U' ) ) then
            kk = 1
            if ( incx == 1 ) then
               do j = 1, n
                  temp = x(j)
                  k    = kk
                  if ( noconj ) then
                     do i = 1, j - 1
                        temp = temp - ap(k) * x(i)
                        k    = k    + 1
                     end do
                     if ( nounit ) then
                        temp = temp/ap( kk + j - 1 )
                     end if
                  else
                     do i = 1, j - 1
                        temp = temp - conjg ( ap(k) ) * x(i)
                        k    = k    + 1
                     end do
                     if ( nounit ) &
                        temp = temp/conjg ( ap( kk + j - 1 ) )
                  end if
                  x(j) = temp
                  kk     = kk   + j
               end do
            else
               jx = kx
               do j = 1, n
                  temp = x(jx)
                  ix = kx
                  if ( noconj ) then
                     do k = kk, kk + j - 2
                        temp = temp - ap(k) * x(ix)
                        ix = ix + incx
                     end do
                     if ( nounit ) then
                        temp = temp/ap( kk + j - 1 )
                     end if
                  else
                     do k = kk, kk + j - 2
                        temp = temp - conjg ( ap(k) ) * x(ix)
                        ix = ix + incx
                     end do
                     if ( nounit ) then
                        temp = temp / conjg ( ap( kk + j - 1 ) )
                     end if
                  end if
                  x(jx) = temp
                  jx = jx + incx
                  kk      = kk   + j
               end do
            end if
         else
            kk = ( n * ( n + 1 ) )/2
            if ( incx == 1 ) then
               do j = n, 1, -1
                  temp = x(j)
                  k    = kk
                  if ( noconj ) then
                     do i = n, j + 1, -1
                        temp = temp - ap(k) * x(i)
                        k    = k    - 1
                     end do
                     if ( nounit ) &
                        temp = temp / ap( kk - n + j )
                  else
                     do i = n, j + 1, -1
                        temp = temp - conjg ( ap(k) ) * x(i)
                        k    = k    - 1
                     end do
                     if ( nounit ) &
                        temp = temp / conjg ( ap( kk - n + j ) )
                  end if
                  x(j) = temp
                  kk     = kk   - ( n - j + 1 )
               end do
            else
               kx = kx + ( n - 1 ) * incx
               jx = kx
               do j = n, 1, -1
                  temp = x(jx)
                  ix = kx
                  if ( noconj ) then
                     do k = kk, kk - ( n - ( j + 1 ) ), -1
                        temp = temp - ap(k) * x(ix)
                        ix = ix - incx
                     end do
                     if ( nounit ) then
                        temp = temp/ap( kk - n + j )
                     end if
                  else
                     do k = kk, kk - ( n - ( j + 1 ) ), -1
                        temp = temp - conjg ( ap(k) ) * x(ix)
                        ix = ix - incx
                     end do
                     if ( nounit ) then
                        temp = temp/conjg ( ap( kk - n + j ) )
                     end if
                  end if
                  x(jx) = temp
                  jx = jx - incx
                  kk      = kk   - ( n - j + 1 )
               end do
            end if
         end if
      end if

  return
end
subroutine ctrmm ( side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb )
!
!*******************************************************************************
!
!! CTRMM performs one of the matrix-matrix operations
!
!     B := alpha*op( A ) * b,   or   B := alpha*B*op( A )
!
!  where  alpha  is a scalar,  B  is an m by n matrix,  A  is a unit, or
!  non-unit,  upper or lower triangular matrix and  op( A )  is one  of
!
!     op( A ) = A   or   op( A ) = A'   or   op( A ) = conjg( A' ).
!
!  Parameters:
!
!  SIDE   - character.
!           On entry,  SIDE specifies whether  op( A ) multiplies B from
!           the left or right as follows:
!
!              SIDE = 'L' or 'l'   B := alpha*op( A ) * b.
!
!              SIDE = 'R' or 'r'   B := alpha * b*op( A ).
!
!           Unchanged on exit.
!
!  UPLO   - character.
!           On entry, UPLO specifies whether the matrix A is an upper or
!           lower triangular matrix as follows:
!
!              UPLO = 'U' or 'u'   A is an upper triangular matrix.
!
!              UPLO = 'L' or 'l'   A is a lower triangular matrix.
!
!           Unchanged on exit.
!
!  TRANSA - character.
!           On entry, TRANSA specifies the form of op( A ) to be used in
!           the matrix multiplication as follows:
!
!              TRANSA = 'N' or 'n'   op( A ) = A.
!
!              TRANSA = 'T' or 't'   op( A ) = A'.
!
!              TRANSA = 'C' or 'c'   op( A ) = conjg( A' ).
!
!           Unchanged on exit.
!
!  DIAG   - character.
!           On entry, DIAG specifies whether or not A is unit triangular
!           as follows:
!
!              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
!
!              DIAG = 'N' or 'n'   A is not assumed to be unit
!                                  triangular.
!
!           Unchanged on exit.
!
!  M      - integer.
!           On entry, M specifies the number of rows of B. M must be at
!           least zero.
!           Unchanged on exit.
!
!  N      - integer.
!           On entry, N specifies the number of columns of B.  N must be
!           at least zero.
!           Unchanged on exit.
!
!  ALPHA  - complex         .
!           On entry,  ALPHA specifies the scalar  alpha. When  alpha is
!           zero then  A is not referenced and  B need not be set before
!           entry.
!           Unchanged on exit.
!
!  A      - complex          array of DIMENSION ( LDA, k ), where k is m
!           when  SIDE = 'L' or 'l'  and is  n  when  SIDE = 'R' or 'r'.
!           Before entry  with  UPLO = 'U' or 'u',  the  leading  k by k
!           upper triangular part of the array  A must contain the upper
!           triangular matrix and the strictly lower triangular part of
!           A is not referenced.
!           Before entry  with  UPLO = 'L' or 'l',  the  leading  k by k
!           lower triangular part of the array  A must contain the lower
!           triangular matrix and the strictly upper triangular part of
!           A is not referenced.
!           Note that when  DIAG = 'U' or 'u',  the diagonal elements of
!           A  are not referenced either,  but are assumed to be  unity.
!           Unchanged on exit.
!
!  LDA    - integer.
!           On entry, LDA specifies the first dimension of A as declared
!           in the calling (sub) program.  When  SIDE = 'L' or 'l'  then
!           LDA  must be at least  max( 1, m ),  when  SIDE = 'R' or 'r'
!           then LDA must be at least max( 1, n ).
!           Unchanged on exit.
!
!  B      - complex          array of DIMENSION ( LDB, n ).
!           Before entry,  the leading  m by n part of the array  B must
!           contain the matrix B,  and  on exit  is overwritten  by the
!           transformed matrix.
!
!  LDB    - integer.
!           On entry, LDB specifies the first dimension of B as declared
!           in  the  calling  (sub)  program.   LDB  must  be  at  least
!           max( 1, m ).
!           Unchanged on exit.
!
  implicit none
!
!     .. Scalar Arguments ..
      character        side, uplo, transa, diag
      integer            m, n, lda, ldb
      complex            alpha
!     .. Array Arguments ..
      complex            a( lda, * ), b( ldb, * )
!     ..
      logical            lsame
      logical            lside, noconj, nounit, upper
      integer            i, info, j, k, nrowa
      complex            temp
      complex            one
      parameter        ( one  = ( 1.0e+0, 0.0e+0 ) )
      complex            zero
      parameter        ( zero = ( 0.0e+0, 0.0e+0 ) )
!
!  Test the input.
!
      lside  = lsame ( side  , 'L' )

      if ( lside ) then
         nrowa = m
      else
         nrowa = n
      end if

      noconj = lsame ( transa, 'T' )
      nounit = lsame ( diag  , 'N' )
      upper  = lsame ( uplo  , 'U' )

      info   = 0
      if (      ( .not.lside                ) .and.  &
               ( .not.lsame ( side  , 'R' ) )      ) then
         info = 1
      else if ( ( .not.upper                ) .and.  &
               ( .not.lsame ( uplo  , 'L' ) )      ) then
         info = 2
      else if ( ( .not.lsame ( transa, 'N' ) ) .and.  &
               ( .not.lsame ( transa, 'T' ) ) .and.  &
               ( .not.lsame ( transa, 'C' ) )      ) then
         info = 3
      else if ( ( .not.lsame ( diag  , 'U' ) ) .and.  &
               ( .not.lsame ( diag  , 'N' ) )      ) then
         info = 4
      else if ( m  <0               ) then
         info = 5
      else if ( n  <0               ) then
         info = 6
      else if ( lda<max ( 1, nrowa ) ) then
         info = 9
      else if ( ldb<max ( 1, m     ) ) then
         info = 11
      end if

      if ( info /= 0 ) then
         call xerbla ( 'ctrmm ', info )
         return
      end if
!
!  Quick return if possible.
!
      if ( n == 0 ) then
         return
      end if
!
!  And when  alpha == zero.
!
      if ( alpha == zero ) then
         b(1:m,1:n) = zero
         return
      end if
!
!  Start the operations.
!
      if ( lside ) then

         if ( lsame ( transa, 'N' ) ) then
!
!  Form  B := alpha*A * b.
!
            if ( upper ) then
               do j = 1, n
                  do k = 1, m
                     if ( b(k,j) /= zero ) then
                        temp = alpha * b(k,j)
                        do i = 1, k - 1
                           b(i,j) = b(i,j) + temp * a(i,k)
                        end do
                        if ( nounit ) then
                           temp = temp * a(k,k)
                        end if
                        b(k,j) = temp
                     end if
                  end do
               end do
            else
               do j = 1, n
                  do k = m, 1, -1
                     if ( b(k,j) /= zero ) then
                        temp = alpha * b(k,j)
                        b(k,j) = temp
                        if ( nounit ) then
                           b(k,j) = b(k,j) * a(k,k)
                        end if
                        do i = k + 1, m
                           b(i,j) = b(i,j) + temp * a(i,k)
                        end do
                     end if
                  end do
               end do
            end if
         else
!
!  Form  B := alpha*A' * b   or   B := alpha*conjg( A' )*B.
!
            if ( upper ) then
               do j = 1, n
                  do i = m, 1, -1
                     temp = b(i,j)
                     if ( noconj ) then
                        if ( nounit ) then
                           temp = temp * a(i,i)
                        end if
                        do k = 1, i - 1
                           temp = temp + a(k,i) * b(k,j)
                        end do
                     else
                        if ( nounit ) then
                           temp = temp*conjg ( a(i,i) )
                        end if
                        do k = 1, i - 1
                           temp = temp + conjg ( a(k,i) ) * b(k,j)
                        end do
                     end if
                     b(i,j) = alpha * temp
                  end do
               end do
            else
               do j = 1, n
                  do i = 1, m
                     temp = b(i,j)
                     if ( noconj ) then
                        if ( nounit ) then
                           temp = temp * a(i,i)
                        end if
                        do k = i + 1, m
                           temp = temp + a(k,i) * b(k,j)
                        end do
                     else
                        if ( nounit ) then
                           temp = temp * conjg ( a(i,i) )
                        end if
                        do k = i + 1, m
                           temp = temp + conjg ( a(k,i) ) * b(k,j)
                        end do
                     end if
                     b(i,j) = alpha * temp
                  end do
               end do
            end if
         end if
      else
         if ( lsame ( transa, 'N' ) ) then
!
!  Form B := alpha * b*A.
!
            if ( upper ) then
               do j = n, 1, -1
                  temp = alpha
                  if ( nounit ) then
                     temp = temp * a(j,j)
                  end if
                  do i = 1, m
                     b(i,j) = temp * b(i,j)
                  end do
                  do k = 1, j - 1
                     if ( a(k,j) /= zero ) then
                        temp = alpha * a(k,j)
                        do i = 1, m
                           b(i,j) = b(i,j) + temp * b(i,k)
                        end do
                     end if
                  end do
               end do
            else
               do j = 1, n
                  temp = alpha
                  if ( nounit ) then
                     temp = temp * a(j,j)
                  end if
                  do i = 1, m
                     b(i,j) = temp * b(i,j)
                  end do
                  do k = j + 1, n
                     if ( a(k,j) /= zero ) then
                        temp = alpha * a(k,j)
                        do i = 1, m
                           b(i,j) = b(i,j) + temp * b(i,k)
                        end do
                     end if
                  end do
               end do
            end if
         else
!
!  Form  B := alpha * b*A'   or   B := alpha*B*conjg( A' ).
!
            if ( upper ) then
               do k = 1, n
                  do j = 1, k - 1
                     if ( a(j,k) /= zero ) then
                        if ( noconj ) then
                           temp = alpha * a(j,k)
                        else
                           temp = alpha * conjg ( a(j,k) )
                        end if
                        do i = 1, m
                           b(i,j) = b(i,j) + temp * b(i,k)
                        end do
                     end if
                  end do
                  temp = alpha
                  if ( nounit ) then
                     if ( noconj ) then
                        temp = temp * a(k,k)
                     else
                        temp = temp*conjg ( a(k,k) )
                     end if
                  end if
                  if ( temp /= one ) then
                     do i = 1, m
                        b(i,k) = temp * b(i,k)
                     end do
                  end if
               end do
            else
               do k = n, 1, -1
                  do j = k + 1, n
                     if ( a(j,k) /= zero ) then
                        if ( noconj ) then
                           temp = alpha * a(j,k)
                        else
                           temp = alpha * conjg ( a(j,k) )
                        end if
                        do i = 1, m
                           b(i,j) = b(i,j) + temp * b(i,k)
                        end do
                     end if
                  end do
                  temp = alpha
                  if ( nounit ) then
                     if ( noconj ) then
                        temp = temp * a(k,k)
                     else
                        temp = temp * conjg ( a(k,k) )
                     end if
                  end if
                  if ( temp /= one ) then
                     do i = 1, m
                        b(i,k) = temp * b(i,k)
                     end do
                  end if
               end do
            end if
         end if
      end if

  return
end
subroutine ctrmv ( uplo, trans, diag, n, a, lda, x, incx )
!
!*******************************************************************************
!
!! CTRMV  performs one of the matrix-vector operations
!
!     x := A*x,   or   x := A'*x,   or   x := conjg( A' )*x,
!
!  where x is an n element vector and  A is an n by n unit, or non-unit,
!  upper or lower triangular matrix.
!
!  Parameters:
!
!  UPLO   - character.
!           On entry, UPLO specifies whether the matrix is an upper or
!           lower triangular matrix as follows:
!
!              UPLO = 'U' or 'u'   A is an upper triangular matrix.
!
!              UPLO = 'L' or 'l'   A is a lower triangular matrix.
!
!           Unchanged on exit.
!
!  TRANS  - character.
!           On entry, TRANS specifies the operation to be performed as
!           follows:
!
!              TRANS = 'N' or 'n'   x := A*x.
!
!              TRANS = 'T' or 't'   x := A'*x.
!
!              TRANS = 'C' or 'c'   x := conjg( A' )*x.
!
!           Unchanged on exit.
!
!  DIAG   - character.
!           On entry, DIAG specifies whether or not A is unit
!           triangular as follows:
!
!              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
!
!              DIAG = 'N' or 'n'   A is not assumed to be unit
!                                  triangular.
!
!           Unchanged on exit.
!
!  N      - integer.
!           On entry, N specifies the order of the matrix A.
!           N must be at least zero.
!           Unchanged on exit.
!
!  A      - complex          array of DIMENSION ( LDA, n ).
!           Before entry with  UPLO = 'U' or 'u', the leading n by n
!           upper triangular part of the array A must contain the upper
!           triangular matrix and the strictly lower triangular part of
!           A is not referenced.
!           Before entry with UPLO = 'L' or 'l', the leading n by n
!           lower triangular part of the array A must contain the lower
!           triangular matrix and the strictly upper triangular part of
!           A is not referenced.
!           Note that when  DIAG = 'U' or 'u', the diagonal elements of
!           A are not referenced either, but are assumed to be unity.
!           Unchanged on exit.
!
!  LDA    - integer.
!           On entry, LDA specifies the first dimension of A as declared
!           in the calling (sub) program. LDA must be at least
!           max( 1, n ).
!           Unchanged on exit.
!
!  X      - complex          array of dimension at least
!           ( 1 + ( n - 1 ) * abs( incx ) ).
!           Before entry, the incremented array X must contain the n
!           element vector x. On exit, X is overwritten with the
!           tranformed vector x.
!
!  incx   - integer.
!           On entry, incx specifies the increment for the elements of
!           X. incx must not be zero.
!           Unchanged on exit.
!
  implicit none
!
      integer            incx, lda, n
      character        diag, trans, uplo
      complex            a( lda, * ), x( * )
      complex            zero
      parameter        ( zero = ( 0.0e+0, 0.0e+0 ) )
      complex            temp
      integer            i, info, ix, j, jx, kx
      logical            noconj, nounit
      logical            lsame
!
!  Test the input.
!
      info = 0
      if   ( .not.lsame ( uplo , 'U' ) .and.  &
               .not.lsame ( uplo , 'L' )      ) then
         info = 1
      else if ( .not.lsame ( trans, 'N' ) .and.  &
               .not.lsame ( trans, 'T' ) .and.  &
               .not.lsame ( trans, 'C' )      ) then
         info = 2
      else if ( .not.lsame ( diag , 'U' ) .and.  &
               .not.lsame ( diag , 'N' )      ) then
         info = 3
      else if ( n<0 ) then
         info = 4
      else if ( lda<max ( 1, n ) ) then
         info = 6
      else if ( incx == 0 ) then
         info = 8
      end if
      if ( info/=0 ) then
         call xerbla ( 'ctrmv ', info )
         return
      end if
!
!  Quick return if possible.
!
      if ( n == 0 ) then
         return
      end if

      noconj = lsame ( trans, 'T' )
      nounit = 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 one pass through A.
!
      if ( lsame ( trans, 'N' ) ) then
!
!  Form  x := A*x.
!
         if ( lsame ( uplo, 'U' ) ) then
            if ( incx == 1 ) then
               do j = 1, n
                  if ( x(j) /= zero ) then
                     temp = x(j)
                     do i = 1, j - 1
                        x(i) = x(i) + temp * a(i,j)
                     end do
                     if ( nounit ) then
                        x(j) = x(j) * a(j,j)
                     end if
                  end if
               end do
            else
               jx = kx
               do j = 1, n
                  if ( x(jx) /= zero ) then
                     temp = x(jx)
                     ix = kx
                     do i = 1, j - 1
                        x(ix) = x(ix) + temp * a(i,j)
                        ix = ix + incx
                     end do
                     if ( nounit ) then
                        x(jx) = x(jx) * a(j,j)
                     end if
                  end if
                  jx = jx + incx
               end do
            end if
         else
            if ( incx == 1 ) then
               do j = n, 1, -1
                  if ( x(j) /= zero ) then
                     temp = x(j)
                     do i = n, j + 1, -1
                        x(i) = x(i) + temp * a(i,j)
                     end do
                     if ( nounit ) then
                        x(j) = x(j) * a(j,j)
                     end if
                  end if
               end do
            else
               kx = kx + ( n - 1 ) * incx
               jx = kx
               do j = n, 1, -1
                  if ( x(jx) /= zero ) then
                     temp = x(jx)
                     ix = kx
                     do i = n, j + 1, -1
                        x(ix) = x(ix) + temp * a(i,j)
                        ix = ix - incx
                     end do
                     if ( nounit ) then
                        x(jx) = x(jx) * a(j,j)
                     end if
                  end if
                  jx = jx - incx
               end do
            end if
         end if
      else
!
!  Form  x := A'*x  or  x := conjg( A' )*x.
!
         if ( lsame ( uplo, 'U' ) ) then
            if ( incx == 1 ) then
               do j = n, 1, -1
                  temp = x(j)
                  if ( noconj ) then
                     if ( nounit ) then
                        temp = temp * a(j,j)
                     end if
                     do i = j - 1, 1, -1
                        temp = temp + a(i,j) * x(i)
                     end do
                  else
                     if ( nounit ) then
                        temp = temp*conjg ( a(j,j) )
                     end if
                     do i = j - 1, 1, -1
                        temp = temp + conjg ( a(i,j) ) * x(i)
                     end do
                  end if
                  x(j) = temp
               end do
            else
               jx = kx + ( n - 1 ) * incx
               do j = n, 1, -1
                  temp = x(jx)
                  ix = jx
                  if ( noconj ) then
                     if ( nounit ) then
                        temp = temp * a(j,j)
                     end if
                     do i = j - 1, 1, -1
                        ix = ix - incx
                        temp = temp + a(i,j) * x(ix)
                     end do
                  else
                     if ( nounit ) then
                        temp = temp * conjg ( a(j,j) )
                     end if
                     do i = j - 1, 1, -1
                        ix = ix - incx
                        temp = temp + conjg ( a(i,j) ) * x(ix)
                     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)
                  if ( noconj ) then
                     if ( nounit ) then
                        temp = temp * a(j,j)
                     end if
                     do i = j + 1, n
                        temp = temp + a(i,j) * x(i)
                     end do
                  else
                     if ( nounit ) then
                        temp = temp * conjg ( a(j,j) )
                     end if
                     do i = j + 1, n
                        temp = temp + conjg ( a(i,j) ) * x(i)
                     end do
                  end if
                  x(j) = temp
               end do
            else
               jx = kx
               do j = 1, n
                  temp = x(jx)
                  ix = jx
                  if ( noconj ) then
                     if ( nounit ) then
                        temp = temp * a(j,j)
                     end if
                     do i = j + 1, n
                        ix = ix + incx
                        temp = temp + a(i,j) * x(ix)
                     end do
                  else
                     if ( nounit ) then
                        temp = temp * conjg ( a(j,j) )
                     end if
                     do i = j + 1, n
                        ix = ix + incx
                        temp = temp + conjg ( a(i,j) ) * x(ix)
                     end do
                  end if
                  x(jx) = temp
                  jx = jx + incx
               end do
            end if
         end if
      end if

  return
end
subroutine ctrsm ( side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb )
!
!*******************************************************************************
!
!! CTRSM  solves one of the matrix equations
!
!     op( A ) * x = alpha * b,   or   X*op( A ) = alpha*B,
!
!  where alpha is a scalar, X and B are m by n matrices, A is a unit, or
!  non-unit,  upper or lower triangular matrix and  op( A )  is one  of
!
!     op( A ) = A   or   op( A ) = A'   or   op( A ) = conjg( A' ).
!
!  The matrix X is overwritten on B.
!
!  Parameters:
!
!  SIDE   - character.
!           On entry, SIDE specifies whether op( A ) appears on the left
!           or right of X as follows:
!
!              SIDE = 'L' or 'l'   op( A ) * x = alpha * b.
!
!              SIDE = 'R' or 'r'   X*op( A ) = alpha * b.
!
!           Unchanged on exit.
!
!  UPLO   - character.
!           On entry, UPLO specifies whether the matrix A is an upper or
!           lower triangular matrix as follows:
!
!              UPLO = 'U' or 'u'   A is an upper triangular matrix.
!
!              UPLO = 'L' or 'l'   A is a lower triangular matrix.
!
!           Unchanged on exit.
!
!  TRANSA - character.
!           On entry, TRANSA specifies the form of op( A ) to be used in
!           the matrix multiplication as follows:
!
!              TRANSA = 'N' or 'n'   op( A ) = A.
!
!              TRANSA = 'T' or 't'   op( A ) = A'.
!
!              TRANSA = 'C' or 'c'   op( A ) = conjg( A' ).
!
!           Unchanged on exit.
!
!  DIAG   - character.
!           On entry, DIAG specifies whether or not A is unit triangular
!           as follows:
!
!              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
!
!              DIAG = 'N' or 'n'   A is not assumed to be unit
!                                  triangular.
!
!           Unchanged on exit.
!
!  M      - integer.
!           On entry, M specifies the number of rows of B. M must be at
!           least zero.
!           Unchanged on exit.
!
!  N      - integer.
!           On entry, N specifies the number of columns of B.  N must be
!           at least zero.
!           Unchanged on exit.
!
!  ALPHA  - complex         .
!           On entry,  ALPHA specifies the scalar  alpha. When  alpha is
!           zero then  A is not referenced and  B need not be set before
!           entry.
!           Unchanged on exit.
!
!  A      - complex          array of DIMENSION ( LDA, k ), where k is m
!           when  SIDE = 'L' or 'l'  and is  n  when  SIDE = 'R' or 'r'.
!           Before entry  with  UPLO = 'U' or 'u',  the  leading  k by k
!           upper triangular part of the array  A must contain the upper
!           triangular matrix and the strictly lower triangular part of
!           A is not referenced.
!           Before entry  with  UPLO = 'L' or 'l',  the  leading  k by k
!           lower triangular part of the array  A must contain the lower
!           triangular matrix and the strictly upper triangular part of
!           A is not referenced.
!           Note that when  DIAG = 'U' or 'u',  the diagonal elements of
!           A  are not referenced either,  but are assumed to be  unity.
!           Unchanged on exit.
!
!  LDA    - integer.
!           On entry, LDA specifies the first dimension of A as declared
!           in the calling (sub) program.  When  SIDE = 'L' or 'l'  then
!           LDA  must be at least  max( 1, m ),  when  SIDE = 'R' or 'r'
!           then LDA must be at least max( 1, n ).
!           Unchanged on exit.
!
!  B      - complex          array of DIMENSION ( LDB, n ).
!           Before entry,  the leading  m by n part of the array  B must
!           contain  the  right-hand  side  matrix B,  and  on exit  is
!           overwritten by the solution matrix X.
!
!  LDB    - integer.
!           On entry, LDB specifies the first dimension of B as declared
!           in  the  calling  (sub)  program.   LDB  must  be  at  least
!           max( 1, m ).
!           Unchanged on exit.
  implicit none
!
!     .. Scalar Arguments ..
      character        side, uplo, transa, diag
      integer            m, n, lda, ldb
      complex            alpha
!     .. Array Arguments ..
      complex            a( lda, * ), b( ldb, * )
      logical            lsame
      logical            lside, noconj, nounit, upper
      integer            i, info, j, k, nrowa
      complex            temp
      complex            one
      parameter        ( one  = ( 1.0e+0, 0.0e+0 ) )
      complex            zero
      parameter        ( zero = ( 0.0e+0, 0.0e+0 ) )
!
!  Test the input.
!
      lside  = lsame ( side  , 'L' )
      if ( lside ) then
         nrowa = m
      else
         nrowa = n
      end if
      noconj = lsame ( transa, 'T' )
      nounit = lsame ( diag  , 'N' )
      upper  = lsame ( uplo  , 'U' )
!
      info   = 0
      if (      ( .not.lside                ) .and.  &
               ( .not.lsame ( side  , 'R' ) )      ) then
         info = 1
      else if ( ( .not.upper                ) .and.  &
               ( .not.lsame ( uplo  , 'L' ) )      ) then
         info = 2
      else if ( ( .not.lsame ( transa, 'N' ) ) .and.  &
               ( .not.lsame ( transa, 'T' ) ) .and.  &
               ( .not.lsame ( transa, 'C' ) )      ) then
         info = 3
      else if ( ( .not.lsame ( diag  , 'U' ) ) .and.  &
               ( .not.lsame ( diag  , 'N' ) )      ) then
         info = 4
      else if ( m  <0               ) then
         info = 5
      else if ( n  <0               ) then
         info = 6
      else if ( lda<max ( 1, nrowa ) ) then
         info = 9
      else if ( ldb<max ( 1, m     ) ) then
         info = 11
      end if
      if ( info/=0 ) then
         call xerbla ( 'ctrsm ', info )
         return
      end if
!
!  Quick return if possible.
!
      if ( n == 0 ) then
         return
      end if
!
!  And when  alpha == zero.
!
      if ( alpha == zero ) then
         b(1:m,1:n) = zero
         return
      end if
!
!  Start the operations.
!
      if ( lside ) then
         if ( lsame ( transa, 'N' ) ) then
!
!  Form  B := alpha*inv( A ) * b.
!
            if ( upper ) then
               do j = 1, n
                  if ( alpha /= one ) then
                     do i = 1, m
                        b(i,j) = alpha * b(i,j)
                     end do
                  end if
                  do k = m, 1, -1
                     if ( b(k,j) /= zero ) then
                        if ( nounit ) then
                           b(k,j) = b(k,j) / a(k,k)
                        end if
                        do i = 1, k - 1
                           b(i,j) = b(i,j) - b(k,j) * a(i,k)
                        end do
                     end if
                  end do
               end do
            else
               do j = 1, n
                  if ( alpha /= one ) then
                     do i = 1, m
                        b(i,j) = alpha * b(i,j)
                     end do
                  end if
                  do k = 1, m
                     if ( b(k,j) /= zero ) then
                        if ( nounit ) then
                           b(k,j) = b(k,j) / a(k,k)
                        end if
                        do i = k + 1, m
                           b(i,j) = b(i,j) - b(k,j) * a(i,k)
                        end do
                     end if
                  end do
               end do
            end if
         else
!
!  Form  B := alpha*inv( A' ) * b
!  or    B := alpha*inv( conjg( A' ) ) * b.
!
            if ( upper ) then
               do j = 1, n
                  do i = 1, m
                     temp = alpha * b(i,j)
                     if ( noconj ) then
                        do k = 1, i - 1
                           temp = temp - a(k,i) * b(k,j)
                        end do
                        if ( nounit ) then
                           temp = temp / a(i,i)
                        end if
                     else
                        do k = 1, i - 1
                           temp = temp - conjg ( a(k,i) ) * b(k,j)
                        end do
                        if ( nounit ) then
                           temp = temp / conjg ( a(i,i) )
                        end if
                     end if
                     b(i,j) = temp
                  end do
               end do
            else
               do j = 1, n
                  do i = m, 1, -1
                     temp = alpha * b(i,j)
                     if ( noconj ) then
                        do k = i + 1, m
                           temp = temp - a(k,i) * b(k,j)
                        end do
                        if ( nounit ) then
                           temp = temp/a(i,i)
                        end if
                     else
                        do k = i + 1, m
                           temp = temp - conjg ( a(k,i) ) * b(k,j)
                        end do
                        if ( nounit ) then
                           temp = temp / conjg ( a(i,i) )
                        end if
                     end if
                     b(i,j) = temp
                  end do
               end do
            end if
         end if
      else
         if ( lsame ( transa, 'N' ) ) then
!
!  Form  B := alpha * b*inv( A ).
!
            if ( upper ) then
               do j = 1, n
                  if ( alpha /= one ) then
                     do i = 1, m
                        b(i,j) = alpha * b(i,j)
                     end do
                  end if
                  do k = 1, j - 1
                     if ( a(k,j) /= zero ) then
                        do i = 1, m
                           b(i,j) = b(i,j) - a(k,j) * b(i,k)
                        end do
                     end if
                  end do
                  if ( nounit ) then
                     temp = one / a(j,j)
                     do i = 1, m
                        b(i,j) = temp * b(i,j)
                     end do
                  end if
               end do
            else
               do j = n, 1, -1
                  if ( alpha/= one ) then
                     do i = 1, m
                        b(i,j) = alpha * b(i,j)
                     end do
                  end if
                  do k = j + 1, n
                     if ( a(k,j) /= zero ) then
                        do i = 1, m
                           b(i,j) = b(i,j) - a(k,j) * b(i,k)
                        end do
                     end if
                  end do
                  if ( nounit ) then
                     temp = one/a(j,j)
                     do i = 1, m
                       b(i,j) = temp * b(i,j)
                     end do
                  end if
               end do
            end if
         else
!
!  Form  B := alpha * b*inv( A' )
!  or    B := alpha * b*inv( conjg( A' ) ).
!
            if ( upper ) then
               do k = n, 1, -1
                  if ( nounit ) then
                     if ( noconj ) then
                        temp = one/a(k,k)
                     else
                        temp = one/conjg ( a(k,k) )
                     end if
                     do i = 1, m
                        b(i,k) = temp * b(i,k)
                     end do
                  end if
                  do j = 1, k - 1
                     if ( a(j,k) /= zero ) then
                        if ( noconj ) then
                           temp = a(j,k)
                        else
                           temp = conjg ( a(j,k) )
                        end if
                        do i = 1, m
                           b(i,j) = b(i,j) - temp * b(i,k)
                        end do
                     end if
                  end do
                  if ( alpha/= one ) then
                     do i = 1, m
                        b(i,k) = alpha * b(i,k)
                     end do
                  end if
               end do
            else
               do k = 1, n
                  if ( nounit ) then
                     if ( noconj ) then
                        temp = one/a(k,k)
                     else
                        temp = one/conjg ( a(k,k) )
                     end if
                     do i = 1, m
                        b(i,k) = temp * b(i,k)
                     end do
                  end if
                  do j = k + 1, n
                     if ( a(j,k) /= zero ) then
                        if ( noconj ) then
                           temp = a(j,k)
                        else
                           temp = conjg ( a(j,k) )
                        end if
                        do i = 1, m
                           b(i,j) = b(i,j) - temp * b(i,k)
                        end do
                     end if
                  end do
                  if ( alpha/= one ) then
                     do i = 1, m
                        b(i,k) = alpha * b(i,k)
                     end do
                  end if
               end do
            end if
         end if
      end if

  return
end
subroutine ctrsv ( uplo, trans, diag, n, a, lda, x, incx )
!
!*******************************************************************************
!
!! CTRSV  solves one of the systems of equations
!
!     A*x = b,   or   A'*x = b,   or   conjg( A' )*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 matrix.
!
!  No test for singularity or near-singularity is included in this
!  routine. Such tests must be performed before calling this routine.
!
!  Parameters:
!
!  UPLO   - character.
!           On entry, UPLO specifies whether the matrix is an upper or
!           lower triangular matrix as follows:
!
!              UPLO = 'U' or 'u'   A is an upper triangular matrix.
!
!              UPLO = 'L' or 'l'   A is a lower triangular matrix.
!
!           Unchanged on exit.
!
!  TRANS  - character.
!           On entry, TRANS specifies the equations to be solved as
!           follows:
!
!              TRANS = 'N' or 'n'   A*x = b.
!
!              TRANS = 'T' or 't'   A'*x = b.
!
!              TRANS = 'C' or 'c'   conjg( A' )*x = b.
!
!           Unchanged on exit.
!
!  DIAG   - character.
!           On entry, DIAG specifies whether or not A is unit
!           triangular as follows:
!
!              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
!
!              DIAG = 'N' or 'n'   A is not assumed to be unit
!                                  triangular.
!
!           Unchanged on exit.
!
!  N      - integer.
!           On entry, N specifies the order of the matrix A.
!           N must be at least zero.
!           Unchanged on exit.
!
!  A      - complex          array of DIMENSION ( LDA, n ).
!           Before entry with  UPLO = 'U' or 'u', the leading n by n
!           upper triangular part of the array A must contain the upper
!           triangular matrix and the strictly lower triangular part of
!           A is not referenced.
!           Before entry with UPLO = 'L' or 'l', the leading n by n
!           lower triangular part of the array A must contain the lower
!           triangular matrix and the strictly upper triangular part of
!           A is not referenced.
!           Note that when  DIAG = 'U' or 'u', the diagonal elements of
!           A are not referenced either, but are assumed to be unity.
!           Unchanged on exit.
!
!  LDA    - integer.
!           On entry, LDA specifies the first dimension of A as declared
!           in the calling (sub) program. LDA must be at least
!           max( 1, n ).
!           Unchanged on exit.
!
!  X      - complex          array of dimension at least
!           ( 1 + ( n - 1 ) * abs( incx ) ).
!           Before entry, the incremented array X must contain the n
!           element right-hand side vector b. On exit, X is overwritten
!           with the solution vector x.
!
!  incx   - integer.
!           On entry, incx specifies the increment for the elements of
!           X. incx must not be zero.
!           Unchanged on exit.
  implicit none
!
      integer            incx, lda, n
      character        diag, trans, uplo
      complex            a( lda, * ), x( * )
      complex            zero
      parameter        ( zero = ( 0.0e+0, 0.0e+0 ) )
      complex            temp
      integer            i, info, ix, j, jx, kx
      logical            noconj, nounit
      logical            lsame
!
!  Test the input.
!
      info = 0
      if   ( .not.lsame ( uplo , 'U' ) .and.  &
               .not.lsame ( uplo , 'L' )      ) then
         info = 1
      else if ( .not.lsame ( trans, 'N' ) .and.  &
               .not.lsame ( trans, 'T' ) .and.  &
               .not.lsame ( trans, 'C' )      ) then
         info = 2
      else if ( .not.lsame ( diag , 'U' ) .and.  &
               .not.lsame ( diag , 'N' )      ) then
         info = 3
      else if ( n<0 ) then
         info = 4
      else if ( lda<max ( 1, n ) ) then
         info = 6
      else if ( incx == 0 ) then
         info = 8
      end if
      if ( info/=0 ) then
         call xerbla ( 'ctrsv ', info )
         return
      end if
!
!  Quick return if possible.
!
      if ( n == 0 ) then
         return
      end if

      noconj = lsame ( trans, 'T' )
      nounit = 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 one pass through A.
!
      if ( lsame ( trans, 'N' ) ) then
!
!  Form  x := inv( A )*x.
!
         if ( lsame ( uplo, 'U' ) ) then
            if ( incx == 1 ) then
               do j = n, 1, -1
                  if ( x(j) /= zero ) then
                     if ( nounit ) then
                        x(j) = x(j)/a(j,j)
                     end if
                     temp = x(j)
                     do i = j - 1, 1, -1
                        x(i) = x(i) - temp * a(i,j)
                     end do
                  end if
               end do
            else
               jx = kx + ( n - 1 ) * incx
               do j = n, 1, -1
                  if ( x(jx) /= zero ) then
                     if ( nounit ) &
                        x(jx) = x(jx)/a(j,j)
                     temp = x(jx)
                     ix = jx
                     do i = j - 1, 1, -1
                        ix = ix - incx
                        x(ix) = x(ix) - temp * a(i,j)
                     end do
                  end if
                  jx = jx - incx
               end do
            end if
         else
            if ( incx == 1 ) then
               do j = 1, n
                  if ( x(j) /= zero ) then
                     if ( nounit ) then
                        x(j) = x(j) / a(j,j)
                     end if
                     temp = x(j)
                     do i = j + 1, n
                        x(i) = x(i) - temp * a(i,j)
                     end do
                  end if
               end do
            else
               jx = kx
               do j = 1, n
                  if ( x(jx) /= zero ) then
                     if ( nounit ) then
                        x(jx) = x(jx) / a(j,j)
                     end if
                     temp = x(jx)
                     ix = jx
                     do i = j + 1, n
                        ix = ix + incx
                        x(ix) = x(ix) - temp * a(i,j)
                     end do
                  end if
                  jx = jx + incx
               end do
            end if
         end if
      else
!
!  Form  x := inv( A' )*x  or  x := inv( conjg( A' ) )*x.
!
         if ( lsame ( uplo, 'U' ) ) then
            if ( incx == 1 ) then
               do j = 1, n
                  temp = x(j)
                  if ( noconj ) then
                     do i = 1, j - 1
                        temp = temp - a(i,j) * x(i)
                     end do
                     if ( nounit ) then
                        temp = temp/a(j,j)
                     end if
                  else
                     do i = 1, j - 1
                        temp = temp - conjg ( a(i,j) ) * x(i)
                     end do
                     if ( nounit ) then
                        temp = temp/conjg ( a(j,j) )
                     end if
                  end if
                  x(j) = temp
               end do
            else
               jx = kx
               do j = 1, n
                  ix = kx
                  temp = x(jx)
                  if ( noconj ) then
                     do i = 1, j - 1
                        temp = temp - a(i,j) * x(ix)
                        ix = ix + incx
                     end do
                     if ( nounit ) then
                        temp = temp/a(j,j)
                     end if
                  else
                     do i = 1, j - 1
                        temp = temp - conjg ( a(i,j) ) * x(ix)
                        ix = ix + incx
                     end do
                     if ( nounit ) &
                        temp = temp/conjg ( a(j,j) )
                  end if
                  x(jx) = temp
                  jx = jx + incx
               end do
            end if
         else
            if ( incx == 1 ) then
               do j = n, 1, -1
                  temp = x(j)
                  if ( noconj ) then
                     do i = n, j + 1, -1
                        temp = temp - a(i,j) * x(i)
                     end do
                     if ( nounit ) &
                        temp = temp/a(j,j)
                  else
                     do i = n, j + 1, -1
                        temp = temp - conjg ( a(i,j) ) * x(i)
                     end do
                     if ( nounit ) then
                        temp = temp/conjg ( a(j,j) )
                     end if
                  end if
                  x(j) = temp
               end do
            else
               kx = kx + ( n - 1 ) * incx
               jx = kx
               do j = n, 1, -1
                  ix = kx
                  temp = x(jx)
                  if ( noconj ) then
                     do i = n, j + 1, -1
                        temp = temp - a(i,j) * x(ix)
                        ix = ix - incx
                     end do
                     if ( nounit ) then
                        temp = temp/a(j,j)
                     end if
                  else
                     do i = n, j + 1, -1
                        temp = temp - conjg ( a(i,j) ) * x(ix)
                        ix = ix - incx
                     end do
                     if ( nounit ) then
                        temp = temp/conjg ( a(j,j) )
                     end if
                  end if
                  x(jx) = temp
                  jx = jx - incx
               end do
            end if
         end if
      end if

  return
end
function icamax ( n, x, incx )
!
!*******************************************************************************
!
!! ICAMAX finds the index of the element having maximum absolute value.
!
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, complex X(*), the vector.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Output, integer ICAMAX, the index of the element of maximum
!    absolute value.
!
  implicit none
!
  real cabs1
  integer i
  integer icamax
  integer incx
  integer ix
  integer n
  real smax
  complex x(*)
!
  icamax = 0
  if ( n < 1 .or. incx  <=  0 ) then
    return
  end if

  icamax = 1

  if ( n == 1 ) then
    return
  end if

  if ( incx /= 1 ) then

    ix = 1
    smax = cabs1 ( x(1) )
    ix = ix + incx
    do i = 2, n
      if ( cabs1 ( x(ix) ) > smax ) then
        icamax = i
        smax = cabs1 ( x(ix) )
      end if
      ix = ix + incx
    end do

  else

    smax = cabs1 ( x(1) )
    do i = 2, n
      if ( cabs1 ( x(i) ) > smax ) then
        icamax = i
        smax = cabs1 ( x(i) )
      end if
    end do

  end if

  return
end
subroutine i_swap ( i, j )
!
!*******************************************************************************
!
!! I_SWAP swaps two integer values.
!
!
!  Modified:
!
!    30 November 1998
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, integer I, J.  On output, the values of I and
!    J have been interchanged.
!
  integer i
  integer j
  integer k
!
  k = i
  i = j
  j = k

  return
end
function lsame ( ca, cb )
!
!*******************************************************************************
!
!! LSAME returns TRUE if CA is the same letter as CA, regardless of case.
!
!
!  -- LAPACK auxiliary routine (version 2.0) --
!     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
!     Courant Institute, Argonne National Lab, and Rice University
!     January 31, 1994
!
!
!  LSAME returns .TRUE. if CA is the same letter as CB regardless of
!  case.
!
!  Parameters:
!
!  CA      (input) character
!
!  CB      (input) character
!          CA and CB specify the single characters to be compared.
!
  implicit none
!
  character ca
  character cb
  integer inta
  integer intb
  logical lsame
  integer zcode
!
!  Test if the characters are equal.
!
  lsame = ca == cb

  if ( lsame ) then
    return
  end if
!
!  Now test for equivalence if both characters are alphabetic.
!
      zcode = ichar ( 'Z' )
!
!  Use 'Z' rather than 'A' so that ASCII can be detected on Prime
!  machines, on which ICHAR returns a value with bit 8 set.
!  ICHAR('A') on Prime machines returns 193 which is the same as
!  ICHAR('A') on an EBCDIC machine.
!
      inta = ICHAR( CA )
      INTB = ICHAR( CB )

      if ( ZCODE == 90 .or. ZCODE == 122 ) then
!
!  ASCII is assumed - ZCODE is the ASCII code of either lower or
!  upper case 'Z'.
!
         if ( inta >=97 .and. inta <= 122 ) inta = inta - 32
         if ( INTB>=97 .and. INTB <= 122 ) INTB = INTB - 32

      else if ( ZCODE == 233 .or. ZCODE == 169 ) then
!
!  EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or
!  upper case 'Z'.
!
         if ( inta >= 129 .and. inta <= 137 .or. &
             inta >= 145 .and. inta <= 153 .or. &
             inta >= 162 .and. inta <= 169 ) inta = inta + 64

         if ( INTB >= 129 .and. INTB <= 137 .or. &
             INTB >= 145 .and. INTB <= 153 .or. &
             INTB >= 162 .and. INTB <= 169 ) INTB = INTB + 64

      else if ( ZCODE == 218 .or. ZCODE == 250 ) then
!
!  ASCII is assumed, on Prime machines - ZCODE is the ASCII code
!  plus 128 of either lower or upper case 'Z'.
!
         if ( inta >=225 .and. inta <= 250 ) inta = inta - 32
         if ( INTB>=225 .and. INTB <= 250 ) INTB = INTB - 32
      end if

  LSAME = inta == INTB

  return
end
function scasum ( n, x, incx )
!
!*******************************************************************************
!
!! SCASUM takes the sum of the absolute values of a complex vector.
!
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, complex X(*), the vector.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Output, real SCASUM, the sum of the absolute values.
!
  implicit none
!
  integer incx
  integer n
  integer nincx
  real scasum
  complex x(*)
!
  scasum = 0.0e0

  if ( n <= 0 .or. incx <= 0 ) then
    return
  end if

  if ( incx == 1 ) then

    scasum = sum ( abs ( real ( x(1:n) ) ) + abs ( aimag ( x(1:n) ) ) )

  else

    nincx = n * incx

    scasum = sum ( abs ( real  ( x(1:nincx:incx) ) ) &
                 + abs ( aimag ( x(1:nincx:incx) ) ) )

  end if

  return
end
function scnrm2 ( n, x, incx )
!
!*******************************************************************************
!
!! SCNRM2 returns the euclidean norm of a complex vector.
!
!
!  Discussion:
!
!    SCNRM2 := sqrt ( sum ( conjg ( x(1:n) ) * x(1:n) ) )
!            = sqrt ( dot_product ( x(1:n), x(1:n) ) )
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, complex X(*), the vector.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Output, real SCNRM2, the norm of the vector.
!
  implicit none
!
  integer incx
  integer ix
  integer n
  real norm
  real, parameter :: one = 1.0e+00
  real scale
  real scnrm2
  real ssq
  real temp
  complex x(*)
  real, parameter :: zero = 0.0e+00
!
  if ( n < 1 .or. incx < 1 ) then

    norm  = zero

  else

    scale = zero
    ssq = one

    do ix = 1, 1 + ( n - 1 ) * incx, incx

      if ( real ( x(ix) ) /= zero ) then
        temp = abs ( real( x(ix) ) )
        if ( scale < temp ) then
          ssq = one + ssq * ( scale / temp )**2
          scale = temp
        else
          ssq = ssq + ( temp / scale )**2
        end if
      end if

      if ( aimag ( x(ix) ) /= zero ) then
        temp = abs ( aimag ( x(ix) ) )
        if ( scale < temp ) then
          ssq = one + ssq * ( scale / temp )**2
          scale = temp
        else
          ssq = ssq + ( temp / scale )**2
        end if

      end if

    end do

    norm  = scale * sqrt ( ssq )

  end if

  scnrm2 = norm

  return
end
subroutine xerbla ( srname, info )
!
!*******************************************************************************
!
!! XERBLA is an error handler for the LAPACK routines.
!
!
!  Discussion:
!
!    XERBLA is called by a LAPACK routine if an input parameter has an
!    invalid value.  A message is printed and execution stops.
!
!    Installers may consider modifying the STOP statement in order to
!    call system-specific exception-handling facilities.
!
!  Parameters:
!
!    Input, character ( len = 6 ) SRNAME,
!    the name of the routine which called XERBLA.
!
!    Input, integer INFO, the position of the invalid parameter in the
!    parameter list of the calling routine.
!
  implicit none
!
  integer info
  character (len = 6 ) srname
!
  write ( *, ' (a,i2,a)' ) 'on entry to ' // trim ( srname ) // &
    ' parameter number ', info, ' had an illegal value.'

  stop
end