function isamax ( n, x, incx )
!
!*******************************************************************************
!
!! ISAMAX finds the index of the vector element of maximum absolute value.
!
!
!  Modified:
!
!    08 April 1999
!
!  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.
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, real X(*), the vector to be examined.
!
!    Input, integer INCX, the increment between successive entries of SX.
!
!    Output, integer ISAMAX, the index of the element of SX of maximum
!    absolute value.
!
  implicit none
!
  integer i
  integer incx
  integer isamax
  integer ix
  integer n
  real samax
  real x(*)
!
  if ( n <= 0 ) then

    isamax = 0

  else if ( n == 1 ) then

    isamax = 1

  else if ( incx == 1 ) then

    isamax = 1
    samax = abs ( x(1) )

    do i = 2, n

      if ( abs ( x(i) ) > samax ) then
        isamax = i
        samax = abs ( x(i) )
      end if

    end do

  else

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

    isamax = 1
    samax = abs ( x(ix) )

    ix = ix + incx

    do i = 2, n
      if ( abs ( x(ix) ) > samax ) then
        isamax = i
        samax = abs ( x(ix) )
      end if
      ix = ix + incx
    end do

  end if

  return
end
function lsame ( ca, cb )
!
!*******************************************************************************
!
!! LSAME tests if CA is the same letter as CB regardless of case.
!
!
!  Discussion:
!
!    This version of the routine is only correct for ASCII code.
!    Installers must modify the routine for other character-codes.
!
!    For EBCDIC systems the constant IOFF must be changed to -64.
!
!  Author:
!
!    Richard Hanson, Sandia National Labs.
!    Jeremy du croz, NAG Central Office.
!
!  Parameters:
!
!    Input, character CA, CB, characters to be compared.  CB is assumed
!    to be an upper case letter.
!
!    Output, logical LSAME, is TRUE if CA is either the same as CB, or
!    is the equivalent lower case letter.
!
  implicit none
!
  character ca
  character cb
  integer, parameter :: ioff = 32
  logical lsame
!
!  Test for equality.
!
  lsame = ca == cb
!
!  Now test for equivalence
!
  if ( .not. lsame ) then
    lsame = ichar(ca) - ioff == ichar(cb)
  end if

  return
end
function sasum ( n, x, incx )
!
!*******************************************************************************
!
!! SASUM sums the absolute values of the entries of a vector.
!
!
!  Modified:
!
!    15 February 2001
!
!  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.
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, real X(*), the vector to be examined.
!
!    Input, integer INCX, the increment between successive entries of X.
!    INCX must not be negative.
!
!    Output, real SASUM, the sum of the absolute values of X.
!
  implicit none
!
  integer incx
  integer n
  real sasum
  real x(*)
!
  sasum = sum ( abs ( x(1:1+(n-1)*incx:incx) ) )

  return
end
subroutine saxpy ( n, sa, x, incx, y, incy )
!
!*******************************************************************************
!
!! SAXPY adds a constant times one vector to another.
!
!
!  Modified:
!
!    08 April 1999
!
!  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.
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, real SA, the multiplier.
!
!    Input, real X(*), the vector to be scaled and added to Y.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Input/output, real Y(*), the vector to which a multiple of X is to
!    be added.
!
!    Input, integer INCY, the increment between successive entries of Y.
!
  implicit none
!
  integer i
  integer incx
  integer incy
  integer ix
  integer iy
  integer n
  real sa
  real x(*)
  real y(*)
!
  if ( n <= 0 ) then

  else if ( sa == 0.0E+00 ) then

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

    y(1:n) = y(1:n) + sa * x(1:n)

  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
      y(iy) = y(iy) + sa * x(ix)
      ix = ix + incx
      iy = iy + incy
    end do

  end if

  return
end
subroutine scopy ( n, x, incx, y, incy )
!
!*******************************************************************************
!
!! SCOPY copies one real vector into another.
!
!
!  Modified:
!
!    08 April 1999
!
!  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.
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, real X(*), the vector to be copied into Y.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Output, real Y(*), the copy of X.
!
!    Input, integer INCY, the increment between successive elements of Y.
!
  implicit none
!
  integer i
  integer incx
  integer incy
  integer ix
  integer iy
  integer n
  real x(*)
  real y(*)
!
  if ( n <= 0 ) then

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

    y(1:n) = x(1:n)

  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
      y(iy) = x(ix)
      ix = ix + incx
      iy = iy + incy
    end do

  end if

  return
end
function sdot ( n, x, incx, y, incy )
!
!*******************************************************************************
!
!! SDOT forms the dot product of two vectors.
!
!
!  Modified:
!
!    02 June 2000
!
!  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.
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vectors.
!
!    Input, real X(*), one of the vectors to be multiplied.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Input, real Y(*), one of the vectors to be multiplied.
!
!    Input, integer INCY, the increment between successive elements of Y.
!
!    Output, real SDOT, the dot product of X and Y.
!
  implicit none
!
  integer i
  integer incx
  integer incy
  integer ix
  integer iy
  integer n
  real sdot
  real stemp
  real x(*)
  real y(*)
!
  if ( n <= 0 ) then

    sdot = 0.0E+00

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

    sdot = dot_product ( x(1:n), y(1:n) )

  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

    stemp = 0.0E+00
    do i = 1, n
      stemp = stemp + x(ix) * y(iy)
      ix = ix + incx
      iy = iy + incy
    end do

    sdot = stemp

  end if

  return
end
subroutine sgbmv ( trans, m, n, kl, ku, alpha, a, lda, x, incx, beta, y, incy )
!
!*******************************************************************************
!
!! SGBMV computes Y := ALPHA * A * X + BETA * Y or Y := ALPHA * A' * X + BETA * Y.
!
!
!  Discussion:
!
!    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.
!
!  Author:
!
!    Jack Dongarra, Argonne National Lab.
!    Jeremy du Croz, NAG Central Office.
!    Sven Hammarling, NAG Central Office.
!    Richard Hanson, Sandia National Labs.
!
!  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 * 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 0.
!           unchanged on exit.
!
!  n - integer.
!           on entry, n specifies the number of columns of the matrix a.
!           n must be at least 0.
!           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 - real            .
!           on entry, alpha specifies the scalar alpha.
!           unchanged on exit.
!
!  a - real             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 - real             array of dimension at least
!           ( 1 + ( n - 1 ) * abs( incx ) ) when trans = 'N' or 'N'
!           and at least
!           ( 1 + ( m - 1 ) * abs( incx ) ) otherwise.
!           before entry, the incremented array x must contain the
!           vector x.
!           unchanged on exit.
!
!  incx - integer.
!           on entry, incx specifies the increment for the elements of
!           x. incx must not be 0.
!           unchanged on exit.
!
!  beta - real            .
!           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 - real             array of dimension at least
!           ( 1 + ( m - 1 ) * abs( incy ) ) when trans = 'N' or 'N'
!           and at least
!           ( 1 + ( n - 1 ) * abs( incy ) ) otherwise.
!           before entry, 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 0.
!           unchanged on exit.
!
  implicit none
!
  integer lda
!
  real a(lda,*)
  real alpha
  real beta
  integer i
  integer incx
  integer incy
  integer info
  integer ix
  integer iy
  integer j
  integer jx
  integer jy
  integer k
  integer kl
  integer ku
  integer kup1
  integer kx
  integer ky
  integer lenx
  integer leny
  logical, external :: lsame
  integer m
  integer n
  character trans
  real temp
  real x(*)
  external xerbla
  real y(*)
!
!  Test the input parameters.
!
  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 ( 'sgbmv ', info )
     return
  end if
!
!  Quick return if possible.
!
  if ( m == 0 ) then
    return
  else if ( n == 0 ) then
    return
  else if ( alpha == 0.0 .and. beta == 1.0 ) then
    return
  end if
!
!  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 /= 1.0 ) then
     if ( incy == 1 ) then
        if ( beta == 0.0 ) then
           y(1:leny) = 0.0
        else
           y(1:leny) = beta * y(1:leny)
        end if
     else
        iy = ky
        if ( beta == 0.0 ) then
           do i = 1, leny
              y(iy) = 0.0
              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 == 0.0 ) 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) /= 0.0 ) 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) /= 0.0 ) 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.
!
     jy = ky
     if ( incx == 1 ) then
        do j = 1, n
           temp = 0.0
           k = kup1 - j
           do i = max ( 1, j - ku ), min ( m, j + kl )
              temp = temp + a( k + i,j) * x(i)
           end do
           y(jy) = y(jy) + alpha * temp
           jy = jy + incy
        end do
     else
        do j = 1, n
           temp = 0.0
           ix = kx
           k = kup1 - j
           do i = max ( 1, j - ku ), min ( m, j + kl )
              temp = temp + a( k + i,j) * x(ix)
              ix = ix + incx
           end do
           y(jy) = y(jy) + alpha * temp
           jy = jy + incy
           if ( j>ku ) then
             kx = kx + incx
           end if
        end do
     end if
  end if

  return
end
subroutine sgemm ( transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, &
  ldc )
!
!*******************************************************************************
!
!! SGEMM performs one of the matrix-matrix operations
!
!     c := alpha*op( a )*op( b ) + beta*c,
!
!  where  op( x ) is one of
!
!     op( x ) = x   or   op( x ) = x',
!
!  alpha and beta are scalars, and a, b and c are matrices, with op( a )
!  an m by k matrix,  op( b )  a  k by n matrix and  c an m by n matrix.
!
!  Author:
!
!    Jack Dongarra, Argonne National Laboratory.
!    Iain Duff, AERE Harwell.
!    Jeremy du Croz, Numerical Algorithms Group LTD.
!    Sven Hammarling, Numerical Algorithms Group LTD.
!
!  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 ) = 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 ) = 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  0.
!           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 0.
!           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  0.
!           unchanged on exit.
!
!  alpha  - real            .
!           on entry, alpha specifies the scalar alpha.
!           unchanged on exit.
!
!  a      - real             array of dimension ( lda, ka ), where ka is
!           k  when  transa = 'n' or 'n',  and is  m  otherwise.
!           before entry with  transa = 'n' or 'n',  the leading  m by k
!           part of the array  a  must contain the matrix  a,  otherwise
!           the leading  k by m  part of the array  a  must contain  the
!           matrix a.
!           unchanged on exit.
!
!  lda    - 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      - real             array of dimension ( ldb, kb ), where kb is
!           n  when  transb = 'n' or 'n',  and is  k  otherwise.
!           before entry with  transb = 'n' or 'n',  the leading  k by n
!           part of the array  b  must contain the matrix  b,  otherwise
!           the leading  n by k  part of the array  b  must contain  the
!           matrix b.
!           unchanged on exit.
!
!  ldb    - 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   - real            .
!           on entry,  beta  specifies the scalar  beta.  when  beta  is
!           supplied as 0.0 then c need not be set on input.
!           unchanged on exit.
!
!  c      - real             array of dimension ( ldc, n ).
!           before entry, the leading  m by n  part of the array  c must
!           contain the matrix  c,  except when  beta  is 0.0, 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, transb
  integer            m, n, k, lda, ldb, ldc
  real               alpha, beta
!     .. array arguments ..
  real               a( lda, * ), b( ldb, * ), c( ldc, * )
!     ..
!     .. external functions ..
  logical            lsame
  external           lsame
!     .. external subroutines ..
  external           xerbla
!     .. intrinsic functions ..
  intrinsic          max
!     .. local scalars ..
  logical            nota, notb
  integer            i, info, j, l, ncola, nrowa, nrowb
  real               temp
!
!     set  nota  and  notb  as  true if  a  and  b  respectively are not
!     transposed 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' )
  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 parameters.
!
  info = 0
  if ( ( .not.nota                 ).and. &
       ( .not.lsame( transa, 'C' ) ).and. &
       ( .not.lsame( transa, 'T' ) )      ) then
     info = 1
  else if ( ( .not.notb                 ).and. &
           ( .not.lsame( transb, 'C' ) ).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( 'sgemm ', info )
     return
  end if
!
!     quick return if possible.
!
  if ( ( m==0 ).or.( n==0 ).or. &
     ( ( ( alpha== 0.0  ).or.( k==0 ) ).and.( beta== 1.0  ) ) ) then
    return
  end if
!
!     and if  alpha== 0.0 .
!
  if ( alpha== 0.0  ) then
     if ( beta== 0.0  ) then
        do j = 1, n
           do i = 1, m
              c( i, j ) = 0.0
           end do
        end do
     else
        do j = 1, n
           do i = 1, m
              c( i, j ) = beta*c( i, j )
           end do
        end do
     end if
     return
  end if
!
!  start the operations.
!
  if ( notb ) then
     if ( nota ) then
!
!  form  c := alpha*a*b + beta*c.
!
        do j = 1, n
           if ( beta== 0.0  ) then
              do i = 1, m
                 c( i, j ) = 0.0
              end do
           else if ( beta/= 1.0  ) then
              do i = 1, m
                 c( i, j ) = beta*c( i, j )
              end do
           end if
           do l = 1, k
              if ( b( l, j )/= 0.0  ) 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
!
!  form  c := alpha*a'*b + beta*c
!
        do j = 1, n
           do i = 1, m
              temp = 0.0
              do l = 1, k
                 temp = temp + a( l, i )*b( l, j )
              end do
              if ( beta== 0.0  ) 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
!
!  form  c := alpha*a*b' + beta*c
!
        do j = 1, n
           if ( beta== 0.0  ) then
              do i = 1, m
                 c( i, j ) = 0.0
              end do
           else if ( beta/= 1.0  ) then
              do i = 1, m
                 c( i, j ) = beta*c( i, j )
              end do
           end if
           do l = 1, k
              if ( b( j, l )/= 0.0  ) 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
     else
!
!  form  c := alpha*a'*b' + beta*c
!
        do j = 1, n
           do i = 1, m
              temp = 0.0
              do l = 1, k
                 temp = temp + a( l, i )*b( j, l )
              end do
              if ( beta== 0.0  ) 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 sgemv ( trans, m, n, alpha, a, lda, x, incx, beta, y, incy )
!
!*******************************************************************************
!
!! SGEMV performs y := alpha * A * x + beta*y or y := alpha * A'*x + beta*y,
!
!  where alpha and beta are scalars, x and y are vectors and a is an
!  m by n matrix.
!
!  Author:
!
!    Jack Dongarra, Argonne National Lab.
!    Jeremy du Croz, NAG Central Office.
!    Sven Hammarling, NAG Central Office.
!    Richard Hanson, Sandia National Labs.
!
!  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 * 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 0.
!           unchanged on exit.
!
!  n - integer.
!           on entry, n specifies the number of columns of the matrix a.
!           n must be at least 0.
!           unchanged on exit.
!
!  alpha - real            .
!           on entry, alpha specifies the scalar alpha.
!           unchanged on exit.
!
!  a - real             array of dimension ( lda, n ).
!           before entry, the leading m by n part of the array a must
!           contain the matrix of coefficients.
!           unchanged on exit.
!
!  lda - 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 - real             array of dimension at least
!           ( 1 + ( n - 1 ) * abs( incx ) ) when trans = 'N' or 'N'
!           and at least
!           ( 1 + ( m - 1 ) * abs( incx ) ) otherwise.
!           before entry, the incremented array x must contain the
!           vector x.
!           unchanged on exit.
!
!  incx - integer.
!           on entry, incx specifies the increment for the elements of
!           x. incx must not be 0.
!           unchanged on exit.
!
!  beta - real            .
!           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 - real             array of dimension at least
!           ( 1 + ( m - 1 ) * abs( incy ) ) when trans = 'N' or 'N'
!           and at least
!           ( 1 + ( n - 1 ) * abs( incy ) ) otherwise.
!           before entry with beta non-zero, the incremented array y
!           must contain the vector y. on exit, y is overwritten by the
!           updated vector y.
!
!  incy - integer.
!           on entry, incy specifies the increment for the elements of
!           y. incy must not be 0.
!           unchanged on exit.
!
  implicit none
!
  real               alpha, beta
  integer            incx, incy, lda, m, n
  character        trans
  real               a( lda, * ), x( * ), y( * )
  real               temp
  integer            i, info, ix, iy, j, jx, jy, kx, ky, lenx, leny
!
  logical            lsame
  external           lsame
!
  external           xerbla
!
!     test the input parameters.
!
  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 ( 'sgemv ', info )
     return
  end if
!
!  Quick return if possible.
!
  if ( ( m == 0 ).or.( n== 0 ).or. &
    ( ( alpha == 0.0 ).and.( beta== 1.0 ) ) ) then
    return
  end if
!
!  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 /= 1.0 ) then
     if ( incy == 1 ) then
        if ( beta == 0.0 ) then
           do i = 1, leny
              y(i) = 0.0
           end do
        else
           do i = 1, leny
              y(i) = beta * y(i)
           end do
        end if
     else
        iy = ky
        if ( beta == 0.0 ) then
           do i = 1, leny
              y(iy) = 0.0
              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 == 0.0 ) 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) /= 0.0 ) 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) /= 0.0 ) 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.
!
     jy = ky
     if ( incx == 1 ) then
        do j = 1, n
           temp = 0.0
           do i = 1, m
              temp = temp + a(i,j) * x(i)
           end do
           y(jy) = y(jy) + alpha * temp
           jy = jy + incy
        end do
     else
        do j = 1, n
           temp = 0.0
           ix = kx
           do i = 1, m
              temp = temp + a(i,j) * x(ix)
              ix = ix + incx
           end do
           y(jy) = y(jy) + alpha * temp
           jy = jy + incy
        end do
     end if
  end if

  return
end
subroutine sger ( m, n, alpha, x, incx, y, incy, a, lda )
!
!*******************************************************************************
!
!! SGER performs the rank 1 operation A := A + alpha * x*y'.
!
!
!  Discussion:
!
!    ALPHA is a scalar, 
!    X is an M element vector, 
!    Y is an N element vector,
!    A is an M by N matrix.
!
!  Author:
!
!    Jack Dongarra, Argonne National Lab.
!    Jeremy du Croz, NAG Central Office.
!    Sven Hammarling, NAG Central Office.
!    Richard Hanson, Sandia National Labs.
!
!  Parameters:
!
!    Input, integer M, the number of rows in A.  M must be at least 0.
!
!  n - integer.
!           on entry, n specifies the number of columns of the matrix a.
!           n must be at least 0.
!           unchanged on exit.
!
!  alpha - real            .
!           on entry, alpha specifies the scalar alpha.
!           unchanged on exit.
!
!  x - real             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 0.
!           unchanged on exit.
!
!  y - real             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 0.
!           unchanged on exit.
!
!  a - real             array of dimension ( lda, n ).
!           before entry, the leading m by n part of the array a must
!           contain the matrix of coefficients. 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
!
  real               alpha
  integer            incx, incy, lda, m, n
!     .. array arguments ..
  real               a( lda, * ), x( * ), y( * )
  real               temp
  integer            i, info, ix, j, jy, kx
!     .. external subroutines ..
  external           xerbla
!
!  Test the input parameters.
!
  info = 0

  if     ( m < 0 ) then
     info = 1
  else if ( n < 0 ) then
     info = 2
  else if ( incx == 0 ) then
     info = 5
  else if ( incy == 0 ) then
     info = 7
  else if ( lda < max ( 1, m ) ) then
     info = 9
  end if

  if ( info /= 0 ) then
     call xerbla ( 'sger  ', info )
     return
  end if
!
!  Quick returns.
!
  if ( m == 0 .or. n== 0 .or. alpha == 0.0 ) 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) /= 0.0 ) 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) /= 0.0 ) 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
function snrm2 ( n, x, incx )
!
!*******************************************************************************
!
!! SNRM2 computes the Euclidean norm of a vector.
!
!
!  Discussion:
!
!    The original SNRM2 algorithm is accurate but written in a bizarre,
!    unreadable and obsolete format.  This version goes for clarity.
!
!  Modified:
!
!    01 June 2000
!
!  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.
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, real X(*), the vector whose norm is to be computed.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Output, real SNRM2, the Euclidean norm of X.
!
  implicit none
!
  integer i
  integer incx
  integer ix
  integer n
  real samax
  real snrm2
  real stemp
  real x(*)
  real xmax
!
  if ( n <= 0 ) then

    snrm2 = 0.0E+00

  else

    xmax = samax ( n, x, incx )

    if ( xmax == 0.0E+00 ) then

      snrm2 = 0.0E+00

    else

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

      stemp = 0.0E+00
      do i = 1, n
        stemp = stemp + ( x(ix) / xmax )**2
        ix = ix + incx
      end do

      snrm2 = xmax * sqrt ( stemp )

    end if

  end if

  return
end
subroutine srot ( n, x, incx, y, incy, c, s )
!
!*******************************************************************************
!
!! SROT applies a plane rotation.
!
!
!  Modified:
!
!    08 April 1999
!
!  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.
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vectors.
!
!    Input/output, real X(*), one of the vectors to be rotated.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Input/output, real Y(*), one of the vectors to be rotated.
!
!    Input, integer INCY, the increment between successive elements of Y.
!
!    Input, real C, S, parameters (presumably the cosine and sine of
!    some angle) that define a plane rotation.
!
  implicit none
!
  real c
  integer i
  integer incx
  integer incy
  integer ix
  integer iy
  integer n
  real s
  real stemp
  real x(*)
  real y(*)
!
  if ( n <= 0 ) then

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

    do i = 1, n
      stemp = c * x(i) + s * y(i)
      y(i) = c * y(i) - s * x(i)
      x(i) = stemp
    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
      stemp = c * x(ix) + s * y(iy)
      y(iy) = c * y(iy) - s * x(ix)
      x(ix) = stemp
      ix = ix + incx
      iy = iy + incy
    end do

  end if

  return
end
subroutine srotg ( sa, sb, c, s )
!
!*******************************************************************************
!
!! SROTG constructs a Givens plane rotation.
!
!
!  Modified:
!
!    08 April 1999
!
!  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.
!
!  Parameters:
!
!    Input/output, real SA, SB, ...
!
!    Output, real C, S, ...
!
  implicit none
!
  real c
  real r
  real roe
  real s
  real sa
  real sb
  real scale
  real z
!
  if ( abs ( sa ) > abs ( sb ) ) then
    roe = sa
  else
    roe = sb
  end if

  scale = abs ( sa ) + abs ( sb )

  if ( scale == 0.0E+00 ) then
    c = 1.0E+00
    s = 0.0E+00
    r = 0.0E+00
  else
    r = scale * sqrt ( ( sa / scale )**2 + ( sb / scale )**2 )
    r = sign ( 1.0E+00, roe ) * r
    c = sa / r
    s = sb / r
  end if

  if ( abs ( c ) > 0.0E+00 .and. abs ( c ) <= s ) then
    z = 1.0E+00 / c
  else
    z = s
  end if

  sa = r
  sb = z

  return
end
subroutine ssbmv ( uplo, n, k, alpha, a, lda, x, incx, beta, y, incy )
!
!*******************************************************************************
!
!! SSBMV 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 symmetric band matrix, with k super-diagonals.
!
!  Author:
!
!    Jack Dongarra, Argonne National Lab.
!    Jeremy du Croz, NAG Central Office.
!    Sven Hammarling, NAG Central Office.
!    Richard Hanson, Sandia National Labs.
!
!  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 0.
!           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 - real            .
!           on entry, alpha specifies the scalar alpha.
!           unchanged on exit.
!
!  a - real             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 symmetric 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 symmetric 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 symmetric 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 symmetric 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
!
!           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 - real             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 0.
!           unchanged on exit.
!
!  beta - real            .
!           on entry, beta specifies the scalar beta.
!           unchanged on exit.
!
!  y - real             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 0.
!           unchanged on exit.
!
  implicit none
!
  real               alpha, beta
  integer            incx, incy, k, lda, n
  character        uplo
!     .. array arguments ..
  real               a( lda, * ), x( * ), y( * )
  real               temp1, temp2
  integer            i, info, ix, iy, j, jx, jy, kplus1, kx, ky, l
!     .. external functions ..
  logical            lsame
  external           lsame
!     .. external subroutines ..
  external           xerbla
!
!  test the input parameters.
!
  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 ( 'ssbmv ', info )
     return
  end if
!
!  Quick returns.
!
  if ( ( n == 0 ).or.( ( alpha== 0.0 ).and.( beta== 1.0 ) ) ) 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 /= 1.0 ) then
     if ( incy == 1 ) then
        if ( beta == 0.0 ) then
           do i = 1, n
              y(i) = 0.0
           end do
        else
           do i = 1, n
              y(i) = beta*y(i)
           end do
        end if
     else
        iy = ky
        if ( beta == 0.0 ) then
           do i = 1, n
              y(iy) = 0.0
              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 == 0.0 ) 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 = 0.0
           l = kplus1 - j
           do i = max ( 1, j - k ), j - 1
              y(i) = y(i) + temp1 * a( l + i,j)
              temp2 = temp2 + a( l + i,j) * x(i)
           end do
           y(j) = y(j) + temp1 * a( kplus1,j) + alpha * temp2
        end do
     else
        jx = kx
        jy = ky
        do j = 1, n
           temp1 = alpha * x(jx)
           temp2 = 0.0
           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 + a( l + i,j) * x(ix)
              ix = ix + incx
              iy = iy + incy
           end do
           y(jy) = y(jy) + temp1 * 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 = 0.0
           y(j) = y(j) + temp1 * 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 + 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 = 0.0
           y(jy) = y(jy) + temp1 * 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 + 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 sscal ( n, sa, x, incx )
!
!*******************************************************************************
!
!! SSCAL scales a vector by a constant.
!
!
!  Modified:
!
!    08 April 1999
!
!  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.
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, real SA, the multiplier.
!
!    Input/output, real X(*), the vector to be scaled.
!
!    Input, integer INCX, the increment between successive entries of X.
!
  implicit none
!
  integer i
  integer incx
  integer ix
  integer m
  integer n
  real sa
  real x(*)
!
  if ( n <= 0 ) then

  else if ( incx == 1 ) then

    m = mod ( n, 5 )

    x(1:m) = sa * x(1:m)

    do i = m+1, n, 5
      x(i)   = sa * x(i)
      x(i+1) = sa * x(i+1)
      x(i+2) = sa * x(i+2)
      x(i+3) = sa * x(i+3)
      x(i+4) = sa * x(i+4)
    end do

  else

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

    do i = 1, n
      x(ix) = sa * x(ix)
      ix = ix + incx
    end do

  end if

  return
end
subroutine sspmv ( uplo, n, alpha, ap, x, incx, beta, y, incy )
!
!*******************************************************************************
!
!! SSPMV 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 symmetric matrix, supplied in packed form.
!
!  Author:
!
!    Jack Dongarra, Argonne National Lab.
!    Jeremy du Croz, NAG Central Office.
!    Sven Hammarling, NAG Central Office.
!    Richard Hanson, Sandia National Labs.
!
!  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 0.
!           unchanged on exit.
!
!  alpha - real            .
!           on entry, alpha specifies the scalar alpha.
!           unchanged on exit.
!
!  ap - real             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 symmetric 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 symmetric 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.
!           unchanged on exit.
!
!  x - real             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 0.
!           unchanged on exit.
!
!  beta - real            .
!           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 - real             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 0.
!           unchanged on exit.
!
  implicit none
!
  real               alpha, beta
  integer            incx, incy, n
  character        uplo
!     .. array arguments ..
  real               ap( * ), x( * ), y( * )
!     ..
!
!     .. local scalars ..
  real               temp1, temp2
  integer            i, info, ix, iy, j, jx, jy, k, kk, kx, ky
!     .. external functions ..
  logical            lsame
  external           lsame
!     .. external subroutines ..
  external           xerbla
!
!     test the input parameters.
!
  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 ( 'sspmv ', info )
     return
  end if
!
!  Quick returns.
!
  if ( n == 0 .or. ( ( alpha== 0.0 ).and.( beta== 1.0 ) ) ) 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 ap
!  are accessed sequentially with one pass through ap.
!
!  First form  y := beta*y.
!
  if ( beta /= 1.0 ) then
     if ( incy == 1 ) then
        if ( beta == 0.0 ) then
           do i = 1, n
              y(i) = 0.0
           end do
        else
           do i = 1, n
              y(i) = beta * y(i)
           end do
        end if
     else
        iy = ky
        if ( beta == 0.0 ) then
           do i = 1, n
              y(iy) = 0.0
              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 == 0.0 ) 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 = 0.0
           k = kk
           do i = 1, j - 1
              y(i) = y(i) + temp1 * ap(k)
              temp2 = temp2 + ap(k) * x(i)
              k = k + 1
           end do
           y(j) = y(j) + temp1 * 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 = 0.0
           ix = kx
           iy = ky
           do k = kk, kk + j - 2
              y(iy) = y(iy) + temp1 * ap(k)
              temp2 = temp2 + ap(k) * x(ix)
              ix = ix + incx
              iy = iy + incy
           end do
           y(jy) = y(jy) + temp1 * 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 = 0.0
           y(j) = y(j) + temp1 * ap(kk)
           k = kk + 1
           do i = j + 1, n
              y(i) = y(i) + temp1 * ap(k)
              temp2 = temp2 + 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 = 0.0
           y(jy) = y(jy) + temp1 * 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 + 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 sspr ( uplo, n, alpha, x, incx, ap )
!
!*******************************************************************************
!
!! SSPR performs the symmetric rank 1 operation A := A + alpha * x*x'.
!
!
!  Discussion:
!
!    ALPHA is a real scalar, X is an N element vector and A is an
!    N by N symmetric matrix, supplied in packed form.
!
!  Author:
!
!    Jack Dongarra, Argonne National Lab.
!    Jeremy du Croz, NAG Central Office.
!    Sven Hammarling, NAG Central Office.
!    Richard Hanson, Sandia National Labs.
!
!  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 0.
!           unchanged on exit.
!
!  alpha - real            .
!           on entry, alpha specifies the scalar alpha.
!           unchanged on exit.
!
!  x - real             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 0.
!           unchanged on exit.
!
!  ap - real             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 symmetric 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 symmetric 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.
!
!
  implicit none
!
  real               alpha
  integer            incx, n
  character        uplo
!     .. array arguments ..
  real               ap( * ), x( * )
!     ..
!
!     .. local scalars ..
  real               temp
  integer            i, info, ix, j, jx, k, kk, kx
!     .. external functions ..
  logical            lsame
  external           lsame
!     .. external subroutines ..
  external           xerbla
!
!  test the input parameters.
!
  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 ( 'sspr  ', info )
     return
  end if
!
!  Quick return if possible.
!
  if ( ( n == 0 ).or.( alpha== 0.0 ) ) 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 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) /= 0.0 ) then
              temp = alpha * x(j)
              k = kk
              do i = 1, j
                 ap(k) = ap( k ) + x(i) * temp
                 k = k + 1
              end do
           end if
           kk = kk + j
        end do
     else
        jx = kx
        do j = 1, n
           if ( x(jx) /= 0.0 ) then
              temp = alpha * x(jx)
              ix = kx
              do k = kk, kk + j - 1
                 ap(k) = ap( k ) + x(ix) * temp
                 ix = ix + incx
              end do
           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) /= 0.0 ) then
              temp = alpha * x(j)
              k = kk
              do i = j, n
                 ap(k) = ap( k ) + x(i) * temp
                 k = k + 1
              end do
           end if
           kk = kk + n - j + 1
        end do
     else
        jx = kx
        do j = 1, n
           if ( x(jx) /= 0.0 ) then
              temp = alpha * x(jx)
              ix = jx
              do k = kk, kk + n - j
                 ap(k) = ap( k ) + x(ix) * temp
                 ix = ix + incx
              end do
           end if
           jx = jx + incx
           kk = kk + n - j + 1
        end do
     end if
  end if

  return
end
subroutine sspr2 ( uplo, n, alpha, x, incx, y, incy, ap )
!
!*******************************************************************************
!
!! SSPR2 performs the symmetric rank 2 operation A := A + alpha * x*y' + alpha * y*x',
!
!  where alpha is a scalar, x and y are n element vectors and a is an
!  n by n symmetric matrix, supplied in packed form.
!
!
!  Author:
!
!    Jack Dongarra, Argonne National Lab.
!    Jeremy du Croz, NAG Central Office.
!    Sven Hammarling, NAG Central Office.
!    Richard Hanson, Sandia National Labs.
!
!  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 0.
!           unchanged on exit.
!
!  alpha - real            .
!           on entry, alpha specifies the scalar alpha.
!           unchanged on exit.
!
!  x - real             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 0.
!           unchanged on exit.
!
!  y - real             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 0.
!           unchanged on exit.
!
!  ap - real             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 symmetric 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 symmetric 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.
!
  implicit none
!
  real               alpha
  integer            incx, incy, n
  character        uplo
!     .. array arguments ..
  real               ap( * ), x( * ), y( * )
  real               temp1, temp2
  integer            i, info, ix, iy, j, jx, jy, k, kk, kx, ky
!     .. external functions ..
  logical            lsame
  external           lsame
!     .. external subroutines ..
  external           xerbla
!
!  Test the input parameters.
!
  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 ( 'sspr2 ', info )
     return
  end if
!
!  Quick return if possible.
!
  if ( ( n == 0 ).or.( alpha== 0.0 ) ) 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 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) /= 0.0 ).or.( y(j)/=0.0 ) ) then
              temp1 = alpha * y(j)
              temp2 = alpha * x(j)
              k = kk
              do i = 1, j
                 ap(k) = ap( k ) + x(i) * temp1 + y(i) * temp2
                 k = k + 1
              end do
           end if
           kk = kk + j
        end do
     else
        do j = 1, n
           if ( ( x(jx) /= 0.0 ).or.( y(jy)/=0.0 ) ) then
              temp1 = alpha * y(jy)
              temp2 = alpha * x(jx)
              ix = kx
              iy = ky
              do k = kk, kk + j - 1
                 ap(k) = ap( k ) + x(ix) * temp1 + y(iy) * temp2
                 ix = ix + incx
                 iy = iy + incy
              end do
           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) /= 0.0 ).or.( y(j)/= 0.0 ) ) then
              temp1 = alpha * y(j)
              temp2 = alpha * x(j)
              k = kk
              do i = j, n
                 ap(k) = ap( k ) + x(i) * temp1 + y(i) * temp2
                 k = k + 1
              end do
           end if
           kk = kk + n - j + 1
        end do
     else
        do j = 1, n
           if ( ( x(jx) /= 0.0 ).or.( y(jy)/= 0.0  ) ) then
              temp1 = alpha * y(jy)
              temp2 = alpha * x(jx)
              ix = jx
              iy = jy
              do k = kk, kk + n - j
                 ap(k) = ap( k ) + x(ix) * temp1 + y(iy) * temp2
                 ix = ix + incx
                 iy = iy + incy
              end do
           end if
           jx = jx + incx
           jy = jy + incy
           kk = kk + n - j + 1
        end do
     end if
  end if

  return
end
subroutine sswap ( n, x, incx, y, incy )
!
!*******************************************************************************
!
!! SSWAP interchanges two vectors.
!
!
!  Modified:
!
!    08 April 1999
!
!  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.
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vectors.
!
!    Input/output, real X(*), one of the vectors to swap.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Input/output, real Y(*), one of the vectors to swap.
!
!    Input, integer INCY, the increment between successive elements of Y.
!
  implicit none
!
  integer i
  integer incx
  integer incy
  integer ix
  integer iy
  integer m
  integer n
  real temp
  real x(*)
  real y(*)
!
  if ( n <= 0 ) then

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

    m = mod ( n, 3 )

    do i = 1, m
      temp = x(i)
      x(i) = y(i)
      y(i) = temp
    end do

    do i = m+1, n, 3

      temp = x(i)
      x(i) = y(i)
      y(i) = temp

      temp = x(i+1)
      x(i+1) = y(i+1)
      y(i+1) = temp

      temp = x(i+2)
      x(i+2) = y(i+2)
      y(i+2) = temp

    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
      temp = x(ix)
      x(ix) = y(iy)
      y(iy) = temp
      ix = ix + incx
      iy = iy + incy
    end do

  end if

  return
end
subroutine ssymm ( side, uplo, m, n, alpha, a, lda, b, ldb, beta, c, ldc )
!
!*******************************************************************************
!
!! SSYMM 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.
!
!
!  Author:
!
!    Jack Dongarra, Argonne National Laboratory.
!    Iain Duff, AERE Harwell.
!    Jeremy du Croz, Numerical Algorithms Group LTD.
!    Sven Hammarling, Numerical Algorithms Group LTD.
!
!  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 0.
!           unchanged on exit.
!
!  n      - integer.
!           on entry, n specifies the number of columns of the matrix c.
!           n  must be at least 0.
!           unchanged on exit.
!
!  alpha  - real            .
!           on entry, alpha specifies the scalar alpha.
!           unchanged on exit.
!
!  a      - real             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      - real             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   - real            .
!           on entry,  beta  specifies the scalar  beta.  when  beta  is
!           supplied as 0.0 then c need not be set on input.
!           unchanged on exit.
!
!  c      - real             array of dimension ( ldc, n ).
!           before entry, the leading  m by n  part of the array  c must
!           contain the matrix  c,  except when  beta  is 0.0, 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
  real               alpha, beta
!     .. array arguments ..
  real               a( lda, * ), b( ldb, * ), c( ldc, * )
!     ..
!
!     .. external functions ..
  logical            lsame
  external           lsame
!     .. external subroutines ..
  external           xerbla
!     .. local scalars ..
  logical            upper
  integer            i, info, j, k, nrowa
  real               temp1, temp2
!
!  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 parameters.
!
  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( 'ssymm ', info )
     return
  end if
!
!  Quick return if possible.
!
  if ( ( m==0 ).or.( n==0 ).or. &
      ( ( alpha== 0.0  ).and.( beta== 1.0  ) ) ) then
    return
  end if
!
!  and when  alpha== 0.0 .
!
  if ( alpha== 0.0  ) then
     if ( beta== 0.0  ) then
        do j = 1, n
           do i = 1, m
              c( i, j ) = 0.0
           end do
        end do
     else
        do j = 1, n
           do i = 1, m
              c( i, j ) = beta*c( i, j )
           end do
        end do
     end if
     return
  end if
!
!  start the operations.
!
  if ( 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 = 0.0
              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== 0.0  ) 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 = 0.0
              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== 0.0  ) 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== 0.0  ) 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 ssymv ( uplo, n, alpha, a, lda, x, incx, beta, y, incy )
!
!*******************************************************************************
!
!! SSYMV 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 symmetric matrix.
!
!
!  Author:
!
!    Jack Dongarra, Argonne National Lab.
!    Jeremy du Croz, NAG Central Office.
!    Sven Hammarling, NAG Central Office.
!    Richard Hanson, Sandia National Labs.
!
!  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 0.
!           unchanged on exit.
!
!  alpha - real            .
!           on entry, alpha specifies the scalar alpha.
!           unchanged on exit.
!
!  a - real             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 symmetric 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 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. lda must be at least
!           max ( 1, n ).
!           unchanged on exit.
!
!  x - real             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 0.
!           unchanged on exit.
!
!  beta - real            .
!           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 - real             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 0.
!           unchanged on exit.
!
  implicit none
!
  real               alpha, beta
  integer            incx, incy, lda, n
  character        uplo
!     .. array arguments ..
  real               a( lda, * ), x( * ), y( * )
!     ..
!
!     .. parameters ..
!     .. local scalars ..
  real               temp1, temp2
  integer            i, info, ix, iy, j, jx, jy, kx, ky
!     .. external functions ..
  logical            lsame
  external           lsame
!     .. external subroutines ..
  external           xerbla
!
!     test the input parameters.
!
  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 ( 'ssymv ', info )
     return
  end if
!
!  Quick return if possible.
!
  if ( n == 0 ) then
    return
  else if ( alpha == 0.0 .and. beta== 1.0 ) 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 /= 1.0 ) then
     if ( incy == 1 ) then
        if ( beta == 0.0 ) then
           do i = 1, n
              y(i) = 0.0
           end do
        else
           do i = 1, n
              y(i) = beta*y(i)
           end do
        end if
     else
        iy = ky
        if ( beta == 0.0 ) then
           do i = 1, n
              y(iy) = 0.0
              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 == 0.0 ) then
    return
  end if

  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 = 0.0
           do i = 1, j - 1
              y(i) = y(i) + temp1 * a(i,j)
              temp2 = temp2 + a(i,j) * x(i)
           end do
           y(j) = y(j) + temp1 * a(j,j) + alpha * temp2
        end do
     else
        jx = kx
        jy = ky
        do j = 1, n
           temp1 = alpha * x(jx)
           temp2 = 0.0
           ix = kx
           iy = ky
           do i = 1, j - 1
              y(iy) = y(iy) + temp1 * a(i,j)
              temp2 = temp2 + a(i,j) * x(ix)
              ix = ix + incx
              iy = iy + incy
           end do
           y(jy) = y(jy) + temp1 * 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 = 0.0
           y(j) = y(j) + temp1 * a(j,j)
           do i = j + 1, n
              y(i) = y(i) + temp1 * a(i,j)
              temp2 = temp2 + 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 = 0.0
           y(jy) = y(jy) + temp1 * 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 + 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 ssyr ( uplo, n, alpha, x, incx, a, lda )
!
!*******************************************************************************
!
!! SSYR performs the symmetric rank 1 operation A := A + alpha * x*x',
!
!  where alpha is a real scalar, x is an n element vector and a is an
!  n by n symmetric matrix.
!
!
!  Author:
!
!    Jack Dongarra, Argonne National Lab.
!    Jeremy du Croz, NAG Central Office.
!    Sven Hammarling, NAG Central Office.
!    Richard Hanson, Sandia National Labs.
!
!  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 0.
!           unchanged on exit.
!
!  alpha - real            .
!           on entry, alpha specifies the scalar alpha.
!           unchanged on exit.
!
!  x - real             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 0.
!           unchanged on exit.
!
!  a - real             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 symmetric 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 symmetric 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.
!
!  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
!     .. array arguments ..
  real               a( lda, * ), x( * )
!     ..
!
!     .. parameters ..
!     .. local scalars ..
  real               temp
  integer            i, info, ix, j, jx, kx
!     .. external functions ..
  logical            lsame
  external           lsame
!     .. external subroutines ..
  external           xerbla
!
!  test the input parameters.
!
  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 ( 'ssyr  ', info )
     return
  end if
!
!  Quick return if possible.
!
  if ( ( n == 0 ).or.( alpha== 0.0 ) ) 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) /= 0.0 ) then
              temp = alpha * x(j)
              do i = 1, j
                 a(i,j) = a(i,j) + x(i) * temp
              end do
           end if
        end do
     else
        jx = kx
        do j = 1, n
           if ( x(jx) /= 0.0 ) then
              temp = alpha * x(jx)
              ix = kx
              do i = 1, j
                 a(i,j) = a(i,j) + x(ix) * temp
                 ix = ix + incx
              end do
           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) /= 0.0 ) then
              temp = alpha * x(j)
              do i = j, n
                 a(i,j) = a(i,j) + x(i) * temp
              end do
           end if
        end do
     else
        jx = kx
        do j = 1, n
           if ( x(jx) /= 0.0 ) then
              temp = alpha * x(jx)
              ix = jx
              do i = j, n
                 a(i,j) = a(i,j) + x(ix) * temp
                 ix = ix + incx
              end do
           end if
           jx = jx + incx
        end do
     end if
  end if

  return
end
subroutine ssyr2 ( uplo, n, alpha, x, incx, y, incy, a, lda )
!
!*******************************************************************************
!
!! SSYR2 performs A := A + alpha * x*y' + alpha * y*x',
!
!  where alpha is a scalar, x and y are n element vectors and a is an n
!  by n symmetric matrix.
!
!
!  Author:
!
!    Jack Dongarra, Argonne National Lab.
!    Jeremy du Croz, NAG Central Office.
!    Sven Hammarling, NAG Central Office.
!    Richard Hanson, Sandia National Labs.
!
!  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 0.
!           unchanged on exit.
!
!  alpha - real            .
!           on entry, alpha specifies the scalar alpha.
!           unchanged on exit.
!
!  x - real             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 0.
!           unchanged on exit.
!
!  y - real             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 0.
!           unchanged on exit.
!
!  a - real             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 symmetric 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 symmetric 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.
!
!  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, incy, lda, n
  character        uplo
!     .. array arguments ..
  real               a( lda, * ), x( * ), y( * )
!     .. parameters ..
!     .. local scalars ..
  real               temp1, temp2
  integer            i, info, ix, iy, j, jx, jy, kx, ky
!     .. external functions ..
  logical            lsame
  external           lsame
!     .. external subroutines ..
  external           xerbla
!
!  test the input parameters.
!
  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 ( 'ssyr2 ', info )
     return
  end if
!
!  Quick return if possible.
!
  if ( ( n == 0 ).or.( alpha== 0.0 ) ) 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) /= 0.0 ).or.( y(j)/= 0.0  ) ) then
              temp1 = alpha * y(j)
              temp2 = alpha * x(j)
              do i = 1, j
                 a(i,j) = a(i,j) + x(i) * temp1 + y(i) * temp2
              end do
           end if
        end do
     else
        do j = 1, n
           if ( ( x(jx) /= 0.0 ).or.( y(jy)/= 0.0  ) ) then
              temp1 = alpha * y(jy)
              temp2 = alpha * x(jx)
              ix = kx
              iy = ky
              do i = 1, j
                 a(i,j) = a(i,j) + x(ix) * temp1 + y(iy) * temp2
                 ix = ix + incx
                 iy = iy + incy
              end do
           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) /= 0.0 ).or.( y(j)/= 0.0  ) ) then
              temp1 = alpha * y(j)
              temp2 = alpha * x(j)
              do i = j, n
                 a(i,j) = a(i,j) + x(i) * temp1 + y(i) * temp2
              end do
           end if
        end do
     else
        do j = 1, n
           if ( ( x(jx) /= 0.0 ).or.( y(jy)/= 0.0  ) ) then
              temp1 = alpha * y(jy)
              temp2 = alpha * x(jx)
              ix = jx
              iy = jy
              do i = j, n
                 a(i,j) = a(i,j) + x(ix) * temp1 + y(iy) * temp2
                 ix = ix + incx
                 iy = iy + incy
              end do
           end if
           jx = jx + incx
           jy = jy + incy
        end do
     end if
  end if

  return
end
subroutine ssyr2k ( uplo, trans, n, k, alpha, a, lda, b, ldb, beta, c, ldc )
!
!*******************************************************************************
!
!! SSYR2K 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.
!
!
!  Author:
!
!    Jack Dongarra, Argonne National Laboratory.
!    Iain Duff, AERE Harwell.
!    Jeremy du Croz, Numerical Algorithms Group LTD.
!    Sven Hammarling, Numerical Algorithms Group LTD.
!
!  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.
!
!              trans = 'c' or 'C'   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 0.
!           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' or 'c' or 'C',  k  specifies  the  number
!           of rows of the matrices  a and b.  k must be at least  0.
!           unchanged on exit.
!
!  alpha  - real            .
!           on entry, alpha specifies the scalar alpha.
!           unchanged on exit.
!
!  a      - real             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      - real             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      - real             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
  real               alpha, beta
!     .. array arguments ..
  real               a( lda, * ), b( ldb, * ), c( ldc, * )
!     ..
!
!     .. external functions ..
  logical            lsame
  external           lsame
!     .. external subroutines ..
  external           xerbla
!     .. local scalars ..
  logical            upper
  integer            i, info, j, l, nrowa
  real               temp1, temp2
!
!  test the input parameters.
!
  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' ) ).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( 'ssyr2k', info )
     return
  end if
!
!  quick return if possible.
!
  if ( ( n==0 ).or. &
    ( ( ( alpha== 0.0  ).or.( k==0 ) ).and.( beta== 1.0  ) ) ) then
    return
  end if
!
!  and when  alpha== 0.0 .
!
  if ( alpha== 0.0  ) then
     if ( upper ) then
        if ( beta== 0.0  ) then
           do j = 1, n
              do i = 1, j
                 c( i, j ) = 0.0
              end do
           end do
        else
           do j = 1, n
              do i = 1, j
                 c( i, j ) = beta*c( i, j )
              end do
           end do
        end if
     else
        if ( beta== 0.0  ) then
           do j = 1, n
              do i = j, n
                 c( i, j ) = 0.0
              end do
           end do
        else
           do j = 1, n
              do i = j, n
                 c( i, j ) = beta*c( i, j )
              end do
           end do
        end if
     end if
     return
  end if
!
!  start the operations.
!
  if ( lsame( trans, 'N' ) ) then
!
!  form  c := alpha*a*b' + alpha*b*a' + c.
!
     if ( upper ) then
        do j = 1, n
           if ( beta== 0.0  ) then
              do i = 1, j
                 c( i, j ) = 0.0
              end do
           else if ( beta/= 1.0  ) then
              do i = 1, j
                 c( i, j ) = beta*c( i, j )
              end do
           end if
           do l = 1, k
              if ( a(j,l) /= 0.0 .or. b(j,l) /= 0.0 ) 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== 0.0  ) then
              do i = j, n
                 c( i, j ) = 0.0
              end do
           else if ( beta/= 1.0  ) then
              do i = j, n
                 c( i, j ) = beta*c( i, j )
              end do
           end if
           do l = 1, k
              if ( a( j, l )/= 0.0 .or. b( j, l )/= 0.0 ) 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 = 0.0
              temp2 = 0.0
              do l = 1, k
                 temp1 = temp1 + a( l, i )*b( l, j )
                 temp2 = temp2 + b( l, i )*a( l, j )
              end do
              if ( beta== 0.0  ) 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 = 0.0
              temp2 = 0.0
              do l = 1, k
                 temp1 = temp1 + a( l, i )*b( l, j )
                 temp2 = temp2 + b( l, i )*a( l, j )
              end do
              if ( beta== 0.0  ) 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 ssyrk ( uplo, trans, n, k, alpha, a, lda, beta, c, ldc )
!
!*******************************************************************************
!
!! SSYRK 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.
!
!
!  Author:
!
!    Jack Dongarra, Argonne National Laboratory.
!    Iain Duff, AERE Harwell.
!    Jeremy du Croz, Numerical Algorithms Group LTD.
!    Sven Hammarling, Numerical Algorithms Group LTD.
!
!  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.
!
!              trans = 'c' or 'C'   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 0.
!           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' or 'c' or 'C',  k  specifies  the  number
!           of rows of the matrix  a.  k must be at least 0.
!           unchanged on exit.
!
!  alpha  - real            .
!           on entry, alpha specifies the scalar alpha.
!           unchanged on exit.
!
!  a      - real             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      - real             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, ldc
  real               alpha, beta
!     .. array arguments ..
  real               a( lda, * ), c( ldc, * )
!     ..
!     .. external functions ..
  logical            lsame
  external           lsame
!     .. external subroutines ..
  external           xerbla
!     .. local scalars ..
  logical            upper
  integer            i, info, j, l, nrowa
  real               temp
!
!  test the input parameters.
!
  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' ) ).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( 'ssyrk ', info )
     return
  end if
!
!  quick return if possible.
!
  if ( ( n==0 ).or. &
    ( ( ( alpha== 0.0  ).or.( k==0 ) ).and.( beta== 1.0  ) ) ) then
    return
  end if
!
!  and when  alpha== 0.0 .
!
  if ( alpha== 0.0  ) then
     if ( upper ) then
        if ( beta== 0.0  ) then
           do j = 1, n
              do i = 1, j
                 c( i, j ) = 0.0
              end do
           end do
        else
           do j = 1, n
              do i = 1, j
                 c( i, j ) = beta*c( i, j )
              end do
           end do
        end if
     else
        if ( beta== 0.0  ) then
           do j = 1, n
              do i = j, n
                 c( i, j ) = 0.0
              end do
           end do
        else
           do j = 1, n
              do i = j, n
                 c( i, j ) = beta*c( i, j )
              end do
           end do
        end if
     end if
     return
  end if
!
!  Start the operations.
!
  if ( lsame( trans, 'N' ) ) then
!
!  form  c := alpha*a*a' + beta*c.
!
     if ( upper ) then
        do j = 1, n
           if ( beta== 0.0  ) then
              do i = 1, j
                 c( i, j ) = 0.0
              end do
           else if ( beta/= 1.0  ) then
              do i = 1, j
                 c( i, j ) = beta*c( i, j )
              end do
           end if
           do l = 1, k
              if ( a( j, l )/= 0.0  ) 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== 0.0  ) then
              do i = j, n
                 c( i, j ) = 0.0
              end do
           else if ( beta/= 1.0  ) then
              do i = j, n
                 c( i, j ) = beta*c( i, j )
              end do
           end if
           do l = 1, k
              if ( a( j, l )/= 0.0  ) 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 = 0.0
              do l = 1, k
                 temp = temp + a( l, i )*a( l, j )
              end do
              if ( beta== 0.0  ) 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 = 0.0
              do l = 1, k
                 temp = temp + a( l, i )*a( l, j )
              end do
              if ( beta== 0.0  ) 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 stbmv ( uplo, trans, diag, n, k, a, lda, x, incx )
!
!*******************************************************************************
!
!! STBMV performs x := A * x or x := 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.
!
!
!  Author:
!
!    Jack Dongarra, Argonne National Lab.
!    Jeremy du Croz, NAG Central Office.
!    Sven Hammarling, NAG Central Office.
!    Richard Hanson, Sandia National Labs.
!
!  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 := 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 0.
!           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 - real             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 - real             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 0.
!           unchanged on exit.
!
  implicit none
!
  integer            incx, k, lda, n
  character        diag, trans, uplo
!     .. array arguments ..
  real               a( lda, * ), x( * )
!     ..
!     .. parameters ..
!     .. local scalars ..
  real               temp
  integer            i, info, ix, j, jx, kplus1, kx, l
  logical            nounit
!     .. external functions ..
  logical            lsame
  external           lsame
!     .. external subroutines ..
  external           xerbla
!
!  test the input parameters.
!
  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 ( 'stbmv ', info )
     return
  end if
!
!  Quick return if possible.
!
  if ( n == 0 ) then
    return
  end if

  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) /= 0.0 ) 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) /= 0.0 ) 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 ) then
                kx = kx + incx
              end if

           end do
        end if
     else
        if ( incx == 1 ) then
           do j = n, 1, -1
              if ( x(j) /= 0.0 ) 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) /= 0.0 ) 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 ) then
                   x(jx) = x(jx) * a(1,j)
                 end if

              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.
!
     if ( lsame ( uplo, 'U' ) ) then
        kplus1 = k + 1
        if ( incx == 1 ) then
           do j = n, 1, -1
              temp = x(j)
              l = kplus1 - j

              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
              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 ( 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
              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 ( 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
              x(j) = temp
           end do
        else
           jx = kx
           do j = 1, n
              temp = x(jx)
              kx = kx + incx
              ix = kx
              l = 1 - j

              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(ix)
                 ix = ix + incx
              end do
              x(jx) = temp
              jx = jx + incx
           end do
        end if
     end if
  end if

  return
end
subroutine stbsv ( uplo, trans, diag, n, k, a, lda, x, incx )
!
!*******************************************************************************
!
!! STBSV solves A * x = b or 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.
!
!
!  Author:
!
!    Jack Dongarra, Argonne National Lab.
!    Jeremy du Croz, NAG Central Office.
!    Sven Hammarling, NAG Central Office.
!    Richard Hanson, Sandia National Labs.
!
!  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'   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 0.
!           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 - real             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 - real             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 0.
!           unchanged on exit.
!
  implicit none
!
  integer            incx, k, lda, n
  character        diag, trans, uplo
!     .. array arguments ..
  real               a( lda, * ), x( * )
!     ..
!
!     .. parameters ..
!     .. local scalars ..
  real               temp
  integer            i, info, ix, j, jx, kplus1, kx, l
  logical            nounit
!     .. external functions ..
  logical            lsame
  external           lsame
!     .. external subroutines ..
  external           xerbla
!
!  test the input parameters.
!
  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 ( 'stbsv ', info )
     return
  end if
!
!  Quick return if possible.
!
  if ( n == 0 ) then
    return
  end if

  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) /= 0.0 ) 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) /= 0.0 ) then
                 ix = kx
                 l = kplus1 - j
                 if ( nounit ) then
                   x(jx) = x(jx) / a(kplus1,j)
                 end if
                 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) /= 0.0 ) 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) /= 0.0 ) 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.
!
     if ( lsame ( uplo, 'U' ) ) then
        kplus1 = k + 1
        if ( incx == 1 ) then
           do j = 1, n
              temp = x(j)
              l = kplus1 - j
              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
              x(j) = temp
           end do
        else
           jx = kx
           do j = 1, n
              temp = x(jx)
              ix = kx
              l = kplus1 - j
              do i = max ( 1, j - k ), j - 1
                 temp = temp - a( l + i,j) * x(ix)
                 ix = ix + incx
              end do
              if ( nounit ) then
                temp = temp/a( kplus1,j)
              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
              do i = min ( n, j + k ), j + 1, -1
                 temp = temp - a( l + i,j) * x(i)
              end do
              if ( nounit ) then
                temp = temp/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
              do i = min ( n, j + k ), j + 1, -1
                 temp = temp - a( l + i,j) * x(ix)
                 ix = ix - incx
              end do
              if ( nounit ) then
                temp = temp/a( 1,j)
              end if
              x(jx) = temp
              jx = jx - incx
              if ( ( n - j ) >= k ) then
                kx = kx - incx
              end if
           end do
        end if
     end if
  end if

  return
end
subroutine stpmv ( uplo, trans, diag, n, ap, x, incx )
!
!*******************************************************************************
!
!! STPMV performs x := A * x or x := 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.
!
!
!  Author:
!
!    Jack Dongarra, Argonne National Lab.
!    Jeremy du Croz, NAG Central Office.
!    Sven Hammarling, NAG Central Office.
!    Richard Hanson, Sandia National Labs.
!
!  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 := 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 0.
!           unchanged on exit.
!
!  ap - real             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 - real             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 0.
!           unchanged on exit.
!
  implicit none
!
  integer            incx, n
  character        diag, trans, uplo
!     .. array arguments ..
  real               ap( * ), x( * )
!     ..
!
!     .. parameters ..
!     .. local scalars ..
  real               temp
  integer            i, info, ix, j, jx, k, kk, kx
  logical            nounit
!     .. external functions ..
  logical            lsame
  external           lsame
!     .. external subroutines ..
  external           xerbla
!
!  test the input parameters.
!
  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 ( 'stpmv ', info )
     return
  end if
!
!  Quick return if possible.
!
  if ( n == 0 ) then
    return
  end if

  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) /= 0.0 ) 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) /= 0.0 ) 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) /= 0.0 ) 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 ) then
                   x(j) = x(j) * ap( kk - n + j )
                 end if
              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) /= 0.0 ) 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 ) then
                   x(jx) = x(jx) * ap( kk - n + j )
                 end if
              end if
              jx = jx - incx
              kk = kk - ( n - j + 1 )
           end do
        end if
     end if
  else
!
!  Form  x := a' * x.
!
     if ( lsame ( uplo, 'U' ) ) then
        kk = ( n*( n + 1 ) )/2
        if ( incx == 1 ) then
           do j = n, 1, -1
              temp = x(j)
              if ( nounit ) then
                temp = temp * ap(kk)
              end if
              k = kk - 1
              do i = j - 1, 1, -1
                 temp = temp + ap(k) * x(i)
                 k = k - 1
              end do
              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 ( 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
              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)
              if ( nounit ) then
                temp = temp * ap(kk)
              end if
              k = kk + 1
              do i = j + 1, n
                 temp = temp + ap(k) * x(i)
                 k = k + 1
              end do
              x(j) = temp
              kk = kk + ( n - j + 1 )
           end do
        else
           jx = kx
           do j = 1, n
              temp = x(jx)
              ix = jx
              if ( nounit ) then
                temp = temp * ap(kk)
              end if
              do k = kk + 1, kk + n - j
                 ix = ix + incx
                 temp = temp + ap(k) * x(ix)
              end do
              x(jx) = temp
              jx = jx + incx
              kk = kk + ( n - j + 1 )
           end do
        end if
     end if
  end if

  return
end
subroutine stpsv ( uplo, trans, diag, n, ap, x, incx )
!
!*******************************************************************************
!
!! STPSV solves A * x = b or 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.
!
!  Author:
!
!    Jack Dongarra, Argonne National Lab.
!    Jeremy du Croz, NAG Central Office.
!    Sven Hammarling, NAG Central Office.
!    Richard Hanson, Sandia National Labs.
!
!  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'   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 0.
!           unchanged on exit.
!
!  ap - real             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 - real             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 0.
!           unchanged on exit.
!
  implicit none
!
  integer            incx, n
  character        diag, trans, uplo
!     .. array arguments ..
  real               ap( * ), x( * )
!     ..
!
!     .. parameters ..
!     .. local scalars ..
  real               temp
  integer            i, info, ix, j, jx, k, kk, kx
  logical            nounit
!     .. external functions ..
  logical            lsame
  external           lsame
!     .. external subroutines ..
  external           xerbla
!
!  test the input parameters.
!
  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 ( 'stpsv ', info )
     return
  end if
!
!  Quick return if possible.
!
  if ( n == 0 ) then
    return
  end if

  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) /= 0.0 ) then
                 if ( nounit ) then
                   x(j) = x(j)/ap(kk)
                 end if
                 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) /= 0.0 ) then
                 if ( nounit ) then
                   x(jx) = x(jx) / ap(kk)
                 end if
                 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) /= 0.0 ) then
                 if ( nounit ) then
                   x(j) = x(j) / ap(kk)
                 end if
                 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) /= 0.0 ) then
                 if ( nounit ) then
                   x(jx) = x(jx) / ap(kk)
                 end if
                 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.
!
     if ( lsame ( uplo, 'U' ) ) then
        kk = 1
        if ( incx == 1 ) then
           do j = 1, n
              temp = x(j)
              k = kk
              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
              x(j) = temp
              kk = kk + j
           end do
        else
           jx = kx
           do j = 1, n
              temp = x(jx)
              ix = kx
              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
              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
              do i = n, j + 1, -1
                 temp = temp - ap(k) * x(i)
                 k = k - 1
              end do
              if ( nounit ) then
                temp = temp / 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
              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
              x(jx) = temp
              jx = jx - incx
              kk = kk - ( n - j + 1 )
           end do
        end if
     end if
  end if

  return
end
subroutine strmm ( side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb )
!
!*******************************************************************************
!
!! STRMM 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'.
!
!
!  Author:
!
!    Jack Dongarra, Argonne National Laboratory.
!    Iain Duff, AERE Harwell.
!    Jeremy du Croz, Numerical Algorithms Group LTD.
!    Sven Hammarling, Numerical Algorithms Group LTD.
!
!  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 ) = 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 0.
!           unchanged on exit.
!
!  n      - integer.
!           on entry, n specifies the number of columns of b.  n must be
!           at least 0.
!           unchanged on exit.
!
!  alpha  - real            .
!           on entry,  alpha specifies the scalar  alpha. when  alpha is
!           0.0 then  a is not referenced and  b need not be set before
!           entry.
!           unchanged on exit.
!
!  a      - real             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      - real             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
!
  character        side, uplo, transa, diag
  integer            m, n, lda, ldb
  real               alpha
!     .. array arguments ..
  real               a( lda, * ), b( ldb, * )
!     ..
!
!     .. external functions ..
  logical            lsame
  external           lsame
!     .. external subroutines ..
  external           xerbla
!     .. local scalars ..
  logical            lside, nounit, upper
  integer            i, info, j, k, nrowa
  real               temp
!
!  test the input parameters.
!
  lside  = lsame( side  , 'L' )
  if ( lside ) then
     nrowa = m
  else
     nrowa = n
  end if
  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( 'strmm ', info )
     return
  end if
!
!  quick return if possible.
!
  if ( n==0 ) then
    return
  end if
!
!  and when  alpha== 0.0 .
!
  if ( alpha== 0.0  ) then
     do j = 1, n
        do i = 1, m
           b( i, j ) = 0.0
        end do
     end do
     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 )/= 0.0  ) 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 )/= 0.0  ) 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*b*a'.
!
        if ( upper ) then
           do j = 1, n
              do i = m, 1, -1
                 temp = b( i, j )
                 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
                 b( i, j ) = alpha*temp
              end do
           end do
        else
           do j = 1, n
              do i = 1, m
                 temp = b( i, j )
                 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
                 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 )/= 0.0  ) 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 )/= 0.0  ) 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'.
!
        if ( upper ) then
           do k = 1, n
              do j = 1, k - 1
                 if ( a( j, k )/= 0.0  ) then
                    temp = alpha*a( j, k )
                    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
                 temp = temp*a( k, k )
              end if
              if ( temp/= 1.0  ) 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 )/= 0.0  ) then
                    temp = alpha*a( j, k )
                    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
                 temp = temp*a( k, k )
              end if
              if ( temp/= 1.0  ) 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 strmv ( uplo, trans, diag, n, a, lda, x, incx )
!
!*******************************************************************************
!
!! STRMV performs x := A * x  or x := 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.
!
!
!  Author:
!
!    Jack Dongarra, Argonne National Lab.
!    Jeremy du Croz, NAG Central Office.
!    Sven Hammarling, NAG Central Office.
!    Richard Hanson, Sandia National Labs.
!
!  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 := 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 0.
!           unchanged on exit.
!
!  a - real             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 - real             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 0.
!           unchanged on exit.
!
  implicit none
!
  integer            incx, lda, n
  character        diag, trans, uplo
!     .. array arguments ..
  real               a( lda, * ), x( * )
!     ..
!
!     .. local scalars ..
  real               temp
  integer            i, info, ix, j, jx, kx
  logical            nounit
!     .. external functions ..
  logical            lsame
  external           lsame
!     .. external subroutines ..
  external           xerbla
!
!  test the input parameters.
!
  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 ( 'strmv ', info )
     return
  end if
!
!  Quick return if possible.
!
  if ( n == 0 ) then
    return
  end if

  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) /= 0.0 ) 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) /= 0.0 ) 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) /= 0.0 ) 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) /= 0.0 ) 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.
!
     if ( lsame ( uplo, 'U' ) ) then
        if ( incx == 1 ) then
           do j = n, 1, -1
              temp = x(j)
              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
              x(j) = temp
           end do
        else
           jx = kx + ( n - 1 ) * incx
           do j = n, 1, -1
              temp = x(jx)
              ix = jx
              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
              x(jx) = temp
              jx = jx - incx
           end do
        end if
     else
        if ( incx == 1 ) then
           do j = 1, n
              temp = x(j)
              if ( nounit ) then
                temp = temp * a(j,j)
              end if
              do i = j + 1, n
                 temp = temp + a(i,j) * x(i)
              end do
              x(j) = temp
           end do
        else
           jx = kx
           do j = 1, n
              temp = x(jx)
              ix = jx
              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
              x(jx) = temp
              jx = jx + incx
           end do
        end if
     end if
  end if

  return
end
subroutine strsm ( side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb )
!
!*******************************************************************************
!
!! STRSM 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'.
!
!  the matrix x is overwritten on b.
!
!
!  Author:
!
!    Jack Dongarra, Argonne National Laboratory.
!    Iain Duff, AERE Harwell.
!    Jeremy du Croz, Numerical Algorithms Group LTD.
!    Sven Hammarling, Numerical Algorithms Group LTD.
!
!  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 ) = 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 0.
!           unchanged on exit.
!
!  n      - integer.
!           on entry, n specifies the number of columns of b.  n must be
!           at least 0.
!           unchanged on exit.
!
!  alpha  - real            .
!           on entry,  alpha specifies the scalar  alpha. when  alpha is
!           0.0 then  a is not referenced and  b need not be set before
!           entry.
!           unchanged on exit.
!
!  a      - real             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      - real             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
!
  character        side, uplo, transa, diag
  integer            m, n, lda, ldb
  real               alpha
!     .. array arguments ..
  real               a( lda, * ), b( ldb, * )
!     ..
!
!     .. external functions ..
  logical            lsame
  external           lsame
!     .. external subroutines ..
  external           xerbla
!     .. local scalars ..
  logical            lside, nounit, upper
  integer            i, info, j, k, nrowa
  real               temp
!
!  test the input parameters.
!
  lside  = lsame( side  , 'L' )
  if ( lside ) then
     nrowa = m
  else
     nrowa = n
  end if
  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( 'strsm ', info )
     return
  end if
!
!  quick return if possible.
!
  if ( n==0 ) then
    return
  end if
!
!  and when  alpha== 0.0 .
!
  if ( alpha== 0.0  ) then
     do j = 1, n
        do i = 1, m
           b( i, j ) = 0.0
        end do
     end do
     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/= 1.0  ) 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 )/= 0.0  ) 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/= 1.0  ) then
                 do i = 1, m
                    b( i, j ) = alpha*b( i, j )
                 end do
              end if
              do k = 1, m
                 if ( b( k, j )/= 0.0  ) 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.
!
        if ( upper ) then
           do j = 1, n
              do i = 1, m
                 temp = alpha*b( i, j )
                 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
                 b( i, j ) = temp
              end do
           end do
        else
           do j = 1, n
              do i = m, 1, -1
                 temp = alpha*b( i, j )
                 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
                 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/= 1.0  ) 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 )/= 0.0  ) 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 = 1.0 / 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/= 1.0  ) 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 )/= 0.0  ) 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 = 1.0/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' ).
!
        if ( upper ) then
           do k = n, 1, -1
              if ( nounit ) then
                 temp = 1.0/a( k, k )
                 do i = 1, m
                    b( i, k ) = temp*b( i, k )
                 end do
              end if
              do j = 1, k - 1
                 if ( a( j, k )/= 0.0  ) then
                    temp = a( j, k )
                    do i = 1, m
                       b( i, j ) = b( i, j ) - temp*b( i, k )
                    end do
                 end if
              end do
              if ( alpha/= 1.0  ) 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
                 temp = 1.0/a( k, k )
                 do i = 1, m
                    b( i, k ) = temp*b( i, k )
                 end do
              end if

              do j = k + 1, n
                 if ( a( j, k )/= 0.0  ) then
                    temp = a( j, k )
                    do i = 1, m
                       b( i, j ) = b( i, j ) - temp*b( i, k )
                    end do
                 end if
              end do

              if ( alpha/= 1.0  ) 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 strsv ( uplo, trans, diag, n, a, lda, x, incx )
!
!*******************************************************************************
!
!! STRSV solves A * x = b or 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.
!
!  Author:
!
!    Jack Dongarra, Argonne National Lab.
!    Jeremy du Croz, NAG Central Office.
!    Sven Hammarling, NAG Central Office.
!    Richard Hanson, Sandia National Labs.
!
!  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'   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 0.
!           unchanged on exit.
!
!  a - real             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 - real             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 0.
!           unchanged on exit.
!
  implicit none
!
  integer            incx, lda, n
  character        diag, trans, uplo
!
  real               a( lda, * ), x( * )
  real               temp
  integer            i, info, ix, j, jx, kx
  logical            nounit
!     .. external functions ..
  logical            lsame
  external           lsame
!     .. external subroutines ..
  external           xerbla
!
!  test the input parameters.
!
  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 ( 'strsv ', info )
     return
  end if
!
!  Quick return if possible.
!
  if ( n == 0 ) then
    return
  end if

  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) /= 0.0 ) 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) /= 0.0 ) then
                 if ( nounit ) then
                   x(jx) = x(jx)/a(j,j)
                 end if
                 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) /= 0.0 ) 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) /= 0.0 ) 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.
!
     if ( lsame ( uplo, 'U' ) ) then
        if ( incx == 1 ) then
           do j = 1, n
              temp = x(j)
              do i = 1, j - 1
                 temp = temp - a(i,j) * x(i)
              end do
              if ( nounit ) then
                temp = temp/a(j,j)
              end if
              x(j) = temp
           end do
        else
           jx = kx
           do j = 1, n
              temp = x(jx)
              ix = kx
              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
              x(jx) = temp
              jx = jx + incx
           end do
        end if
     else
        if ( incx == 1 ) then
           do j = n, 1, -1
              temp = x(j)
              do i = n, j + 1, -1
                 temp = temp - a(i,j) * x(i)
              end do
              if ( nounit ) then
                temp = temp/a(j,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
              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
              x(jx) = temp
              jx = jx - incx
           end do
        end if
     end if
  end if

  return
end
subroutine xerbla ( srname, info )
!
!*******************************************************************************
!
!! XERBLA is an error handler for the level 2 BLAS routines.
!
!
!  Discussion:
!
!    This routine is called by the level 2 blas routines if an input 
!    parameter is invalid.
!
!    Installers should consider modifying the stop statement in order to
!    call system-specific exception-handling facilities.
!
!  Parameters:
!
!    srname - character*6.
!           on entry, srname specifies the name of the routine which
!           called xerbla.
!
!    info - integer.
!           on entry, info specifies the position of the invalid
!           parameter in the parameter-list of the calling routine.
!
  implicit none
!
  integer info
  character ( len = 6 ) srname
!
  write ( *, 99999 ) srname, info

  stop

99999 format ( ' ** on entry to ', a6, ' parameter number ', i2, &
              ' had an illegal value' )

end