function y = dsbmv ( uplo, n, k, alpha, a, lda, x, incx, beta, y, incy )

%*****************************************************************************80
%
%% DSBMV computes y := alpha*A*x + beta*y for symmetric banded matrix A.
%
%  Discussion:
%
%    DSBMV 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.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 March 2017
%
%  Author:
%
%    Original FORTRAN77 version by Jack Dongarra,  Jeremy Du Croz,  
%    Sven Hammarling,  Richard Hanson.
%    This MATLAB version by John Burkardt.
%
%  Parameters:
%
%    Input, character UPLO, specifies whether the upper or lower
%    triangular part of the band matrix A is being supplied:
%    'u' or 'U': the upper triangular part of A is being supplied.
%    'l' or 'L': the lower triangular part of A is being supplied.
%
%    Input, integer N, the number of rows and columns of the 
%    matrix.  0 <= N.
%
%    Input, integer K, the number of super-diagonals of the
%    matrix A.  0 <= K.
%
%    Input, real ALPHA, the first scalar multiplier.
%
%    Input, real A(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
%      end
%    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
%      end
%
%    Input, integer LDA, the first dimension of A as declared
%    in the calling program. K + 1 <= LDA.
%
%    Input, real X(1+(n-1)*abs(INCX)).
%    Before entry, the incremented array X must contain the vector x.
%
%    Input, integer INCX, the increment for elements of X.
%    INCX must not be zero.
%
%    Input, real BETA, the second scalar multiplier.
%
%    Input/output, real Y(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.
%
%    Input, integer INCY, the increment for elements of Y.
%    INCY must not be zero.
%

%
%  Test the input parameters.
%
  info = 0;
  if ( uplo ~= 'u' && uplo ~= 'U' && uplo ~= 'l' && uplo ~= 'L' )
    info = 1;
  elseif ( n < 0 )
    info = 2;
  elseif ( k < 0 )
    info = 3;
  elseif ( lda < ( k + 1 ) )
    info = 6;
  elseif ( incx == 0 )
    info = 8;
  elseif ( incy == 0 )
    info = 11;
  end

  if ( info ~= 0 )
    printf ( 'DSBMV: illegal value for parameter %d\n', info );
    return
  end
%
%  Quick return if possible.
%
  if ( ( n == 0 ) || ( ( alpha == 0.0 ) && ( beta == 1.0 ) ) )
    return
  end
%
%  Set up the start points in  X  and  Y.
%
  if ( 0 < incx )
    kx = 1;
  else
    kx = 1 - ( n - 1 ) * incx;
  end

  if ( 0 < incy )
    ky = 1;
  else
    ky = 1 - ( n - 1 ) * incy;
  end
%
%  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 )
    if ( incy == 1 )
      if ( beta == 0.0 )
        for i = 1 : n
          y(i) = 0.0;
        end
      else
        for i = 1 : n
          y(i) = beta * y(i);
        end
      end
    else
      iy = ky;
      if ( beta == 0.0 )
        for i = 1 : n
          y(iy) = 0.0;
          iy = iy + incy;
        end
      else
        for i = 1 : n
          y(iy) = beta * y(iy);
          iy = iy + incy;
        end
      end
    end
  end

  if ( alpha == 0.0 )
    return
  end

  if ( uplo == 'u' || uplo == 'U' )
%
%  Form y when upper triangle of A is stored.
%
    kplus1 = k + 1;
    if ( ( incx == 1 ) && ( incy == 1 ) )
      for j = 1 : n
        temp1 = alpha * x(j);
        temp2 = 0.0;
        l = kplus1 - j;
        for 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
        y(j) = y(j) + temp1 * a(kplus1,j) + alpha * temp2;
      end
    else
      jx = kx;
      jy = ky;
      for j = 1 : n
        temp1 = alpha * x(jx);
        temp2 = 0.0;
        ix = kx;
        iy = ky;
        l = kplus1 - j;
        for 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
        y(jy) = y(jy) + temp1 * a(kplus1,j) + alpha * temp2;
        jx = jx + incx;
        jy = jy + incy;
        if ( k < j )
          kx = kx + incx;
          ky = ky + incy;
        end
      end
    end
  else
%
%  Form y when lower triangle of A is stored.
%
    if ( ( incx == 1 ) && ( incy == 1 ) )
      for j = 1 : n
        temp1 = alpha * x(j);
        temp2 = 0.0;
        y(j) = y(j) + temp1 * a(1,j);
        l = 1 - j;
        for 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
        y(j) = y(j) + alpha * temp2;
      end
    else
      jx = kx;
      jy = ky;
      for 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;
        for 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
        y(jy) = y(jy) + alpha * temp2;
        jx = jx + incx;
        jy = jy + incy;
      end
    end
  end

  return
end