function ap = dspr ( uplo, n, alpha, x, incx, ap )

%*****************************************************************************80
%
%% DSPR computes A := alpha*x*x' + A, where A is a symmetric packed matrix.
%
%  Discussion:
%
%    DSPR performs the symmetric rank 1 operation
%
%      A := alpha*x*x' + A,
%
%    where alpha is a real scalar, x is an n element vector and A is an
%    n by n symmetric matrix, supplied in packed form.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 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 scalar multiplier.
%
%    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/output, 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. 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.
%

%
%  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 = 5;
  end

  if ( info ~= 0 )
    fprintf ( 1, 'DSPR: Parameter %d had an illegal value.\n', info );
    return
  end
%
%  Quick return if possible.
%
  if ( ( n == 0 ) || ( alpha == 0.0 ) )
    return
  end
%
%  Set the start point in X if the increment is not unity.
%
  if ( incx <= 0 )
    kx = 1 - ( n - 1 ) * incx;
  elseif ( incx ~= 1 )
    kx = 1;
  end
%
%  Start the operations. In this version the elements of the array AP
%  are accessed sequentially with one pass through AP.
%
  kk = 1;
  if ( uplo == 'u' || uplo == 'U' )
%
%  Form A when upper triangle is stored in AP.
%
    if ( incx == 1 )
      for j = 1 : n
        if ( x(j) ~= 0.0 )
          temp = alpha * x(j);
          k = kk;
          for i = 1 : j
            ap(k) = ap(k) + x(i) * temp;
            k = k + 1;
          end
        end
        kk = kk + j;
      end
    else
      jx = kx;
      for j = 1 : n
        if ( x(jx) ~= 0.0 )
          temp = alpha * x(jx);
          ix = kx;
          for k = kk : kk + j - 1
            ap(k) = ap(k) + x(ix) * temp;
            ix = ix + incx;
          end
        end
        jx = jx + incx;
        kk = kk + j;
      end
    end
  else
%
%  Form A when lower triangle is stored in AP.
%
    if ( incx == 1 )
      for j = 1 : n
        if ( x(j) ~= 0.0 )
          temp = alpha * x(j);
          k = kk;
          for i = j : n
            ap(k) = ap(k) + x(i) * temp;
            k = k + 1;
          end
        end
        kk = kk + n - j + 1;
      end
    else
      jx = kx;
      for j = 1 : n
        if ( x(jx) ~= 0.0 )
          temp = alpha * x(jx);
          ix = jx;
          for k = kk : kk + n - j
            ap(k) = ap(k) + x(ix) * temp;
            ix = ix + incx;
          end
        end
        jx = jx + incx;
        kk = kk + n - j + 1;
      end
    end
  end

  return
end