function y = dspmv ( uplo, n, alpha, ap, x, incx, beta, y, incy )

%*****************************************************************************80
%
%% DSPMV computes y := alpha*A*x + beta*y for symmetric packed matrix A.
%
%  Discussion:
%
%    DSPMV 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.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 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 matrix A is being supplied in AP:
%    '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, real ALPHA, the first scalar multiplier.
%
%    Input, real AP((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.
%
%    Input, real X(1+(n-1)*abs(INCX)).
%    Before entry, the incremented array X must contain the n element 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+(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.
%
%    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 ( incx == 0 )
    info = 6;
  elseif ( incy == 0 )
    info = 9;
  end

  if ( info ~= 0 )
    fprintf ( 1, 'DSPMV - 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 AP
%  are accessed sequentially with one pass through AP.
%
%  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

  kk = 1;
  if ( uplo == 'u' || uplo == 'U' )
%
%  Form y when AP contains the upper triangle.
%
    if ( ( incx == 1 ) && ( incy == 1 ) )
      for j = 1 : n
        temp1 = alpha * x(j);
        temp2 = 0.0;
        k = kk;
        for i = 1 : j - 1
          y(i) = y(i) + temp1 * ap(k);
          temp2 = temp2 + ap(k) * x(i);
          k = k + 1;
        end
        y(j) = y(j) + temp1 * ap(kk+j-1) + alpha * temp2;
        kk = kk + j;
      end
    else
      jx = kx;
      jy = ky;
      for j = 1 : n
        temp1 = alpha * x(jx);
        temp2 = 0.0;
        ix = kx;
        iy = ky;
        for 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
        y(jy) = y(jy) + temp1 * ap(kk+j-1) + alpha * temp2;
        jx = jx + incx;
        jy = jy + incy;
        kk = kk + j;
      end
    end
  else
%
%  Form y when AP contains the lower triangle.
%
    if ( ( incx == 1 ) && ( incy == 1 ) )
      for j = 1 : n
        temp1 = alpha * x(j);
        temp2 = 0.0;
        y(j) = y(j) + temp1 * ap(kk);
        k = kk + 1;
        for i = j + 1 : n
          y(i) = y(i) + temp1 * ap(k);
          temp2 = temp2 + ap(k) * x(i);
          k = k + 1;
        end
        y(j) = y(j) + alpha * temp2;
        kk = kk + ( n - j + 1 );
      end
    else
      jx = kx;
      jy = ky;
      for j = 1 : n
        temp1 = alpha * x(jx);
        temp2 = 0.0;
        y(jy) = y(jy) + temp1 * ap(kk);
        ix = jx;
        iy = jy;
        for 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
        y(jy) = y(jy) + alpha * temp2;
        jx = jx + incx;
        jy = jy + incy;
        kk = kk + ( n - j + 1 );
      end
    end
  end

  return
end
