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.
!
!  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.
!
!
!  level 3 blas routine.
!
!  -- written on 8-february-1989.
!     jack dongarra, argonne national laboratory.
!     iain duff, aere harwell.
!     jeremy du croz, numerical algorithms group ltd.
!     sven hammarling, numerical algorithms group ltd.
!
  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 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.
!
!  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.
!
!
!  level 3 blas routine.
!
!  -- written on 8-february-1989.
!     jack dongarra, argonne national laboratory.
!     iain duff, aere harwell.
!     jeremy du croz, numerical algorithms group ltd.
!     sven hammarling, numerical algorithms group ltd.
!
  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 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.
!
!  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.
!
!
!  level 3 blas routine.
!
!  -- written on 8-february-1989.
!     jack dongarra, argonne national laboratory.
!     iain duff, aere harwell.
!     jeremy du croz, numerical algorithms group ltd.
!     sven hammarling, numerical algorithms group ltd.
!
!     .. scalar arguments ..
  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 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.
!
!  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.
!
!
!  level 3 blas routine.
!
!
!  -- written on 8-february-1989.
!     jack dongarra, argonne national laboratory.
!     iain duff, aere harwell.
!     jeremy du croz, numerical algorithms group ltd.
!     sven hammarling, numerical algorithms group ltd.
!
!     .. scalar arguments ..
  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 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'.
!
!  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.
!
!
!  level 3 blas routine.
!
!  -- written on 8-february-1989.
!     jack dongarra, argonne national laboratory.
!     iain duff, aere harwell.
!     jeremy du croz, numerical algorithms group ltd.
!     sven hammarling, numerical algorithms group ltd.
!
!     .. scalar arguments ..
  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 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.
!
!  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.
!
!
!  level 3 blas routine.
!
!
!  -- written on 8-february-1989.
!     jack dongarra, argonne national laboratory.
!     iain duff, aere harwell.
!     jeremy du croz, numerical algorithms group ltd.
!     sven hammarling, numerical algorithms group ltd.
!
!     .. scalar arguments ..
  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