function y = dgbmv ( trans, m, n, kl, ku, alpha, a, lda, x, incx, beta, y, incy )

%*****************************************************************************80
%
%% DGBMV computes y: = alpha*A*x + beta*y for banded matrix A.
%
%  Discussion:
%
%    DGBMV performs one of the matrix-vector operations
%
%      y := alpha*A*x + beta*y,   or   y := alpha*transpose(A)*x + beta*y,
%
%    where alpha and beta are scalars, x and y are vectors and A is an
%    m by n band matrix, with kl sub-diagonals and ku super-diagonals.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 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 TRANS, the operation to be performed:
%    'n' or 'N'   y := alpha*A*x + beta*y.
%    't' or 'T'   y := alpha*A'*x + beta*y.
%    'c' or 'C'   y := alpha*A'*x + beta*y.
%
%    Input, integer M, the number of rows of the matrix A.
%    0 <= M.
%
%    Input, integer N, the number of columns of the matrix A.
%    0 <= N.
%
%    Input, integer KL, the number of subdiagonals.
%    0 <= KL.
%
%    Input, integer KU, the number of superdiagonals.
%    0 <= KU.
%
%    Input, real ALPHA, the first scalar multiplier.
%
%    Input, real A(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
%      end
%
%    Input, integer LDA, the first dimension of A as declared
%    in the calling routine.  KL+KU+1 <= LDA.
%
%    Input, real X(1+(n-1)*abs(INCX)) when TRANS = 'n' or 'N',
%    or real X((1+(m-1)*abs(INCX)) otherwise.
%    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.  When BETA is
%    supplied as zero then Y need not be set on input.
%
%    Input/output, real Y(1+(m-1)*abs(INCY)) if TRANS = 'n' or 'N',
%    real Y(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.
%
%    Input, integer INCY, the increment for the elements of
%    Y.  INCY must not be zero.
%
  x = x(:);
  y = y(:);
%
%  Test the input parameters.
%
  info = 0;
  if ( trans ~= 'n' && trans ~= 'N' && ...
       trans ~= 't' && trans ~= 'T' && ...
       trans ~= 'c' && trans ~= 'C' )
    info = 1;
  elseif ( m < 0 )
    info = 2;
  elseif ( n < 0 )
    info = 3;
  elseif ( kl < 0 )
    info = 4;
  elseif ( ku < 0 )
    info = 5;
  elseif ( lda < ( kl + ku + 1 ) )
    info = 8;
  elseif ( incx == 0 )
    info = 10;
  elseif ( incy == 0 )
    info = 13;
  end

  if ( info ~= 0 )
    printf ( 'DGBMV: Parameter %d had an illegal value.\n', info );
    return
  end
%
%  Quick return if possible.
%
  if ( m == 0 )
    return
  elseif ( n == 0 )
    return
  elseif ( alpha == 0.0 && beta == 1.0 )
    return
  end
%
%  Set  LENX  and  LENY, the lengths of the vectors x and y, and set
%  up the start points in  X  and  Y.
%
  if ( trans == 'N' || trans == 'n' )
    lenx = n;
    leny = m;
  else
    lenx = m;
    leny = n;
  end

  if ( 0 < incx )
    kx = 1;
  else
    kx = 1 - ( lenx - 1 ) * incx;
  end

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

  if ( alpha == 0.0 )
    return
  end

  kup1 = ku + 1;

  if ( trans == 'n' || trans == 'N' )
%
%  Form  y := alpha*A*x + y.
%
    jx = kx;
    if ( incy == 1 )
      for j = 1 : n
        if ( x(jx) ~= 0.0 )
          temp = alpha * x(jx);
          k = kup1 - j;
          for i = max ( 1, j - ku ) : min ( m, j + kl )
            y(i) = y(i) + temp * a(k+i,j);
          end
        end
        jx = jx + incx;
      end
    else
      for j = 1 : n
        if ( x(jx) ~= 0.0 )
          temp = alpha * x(jx);
          iy = ky;
          k = kup1 - j;
          for i = max ( 1, j - ku ) : min ( m, j + kl )
            y(iy) = y(iy) + temp * a(k+i,j);
            iy = iy + incy;
          end
        end
        jx = jx + incx;
        if ( ku < j )
         ky = ky + incy;
        end

      end
    end
  else
%
%  Form y := alpha*transpose(A)*x + y.
%
    jy = ky;
    if ( incx == 1 )
      for j = 1 : n
        temp = 0.0;
        k = kup1 - j
        for i = max ( 1, j - ku ) : min ( m, j + kl )
          temp = temp + a(k+i,j) * x(i);
        end
        y(jy) = y(jy) + alpha * temp;
        jy = jy + incy;
      end
    else
      for j = 1 : n
        temp = 0.0
        ix = kx
        k = kup1 - j
        for i = max ( 1, j - ku ) : min ( m, j + kl )
          temp = temp + a(k+i,j) * x(ix);
          ix = ix + incx;
        end
        y(jy) = y(jy) + alpha * temp;
        jy = jy + incy;

        if ( ku < j )
          kx = kx + incx;
        end

      end
    end
  end

  return
end