[uplo,n,alpha,x,incx,ap]=chpr(uplo,n,alpha,x,incx,ap);

function [uplo,n,alpha,x,incx,ap]=chpr(uplo,n,alpha,x,incx,ap);
%***BEGIN PROLOGUE  CHPR
%***PURPOSE  Perform the hermitian rank 1 operation.
%***LIBRARY   SLATEC (BLAS)
%***CATEGORY  D1B4
%***TYPE      COMPLEX (CHPR-C)
%***KEYWORDS  LEVEL 2 BLAS, LINEAR ALGEBRA
%***AUTHOR  Dongarra, J. J., (ANL)
%           Du Croz, J., (NAG)
%           Hammarling, S., (NAG)
%           Hanson, R. J., (SNLA)
%***DESCRIPTION
%
%  CHPR    performs the hermitian rank 1 operation
%
%     A := alpha*x*conjg( x') + A,
%
%  where alpha is a real scalar, x is an n element vector and A is an
%  n by n hermitian matrix, supplied in packed form.
%
%  Parameters
%  ==========
%
%  UPLO   - CHARACTER*1.
%           On entry, UPLO specifies whether the upper or lower
%           triangular part of the matrix A is supplied in the packed
%           array AP as follows:
%
%              UPLO = 'U' or 'u'   The upper triangular part of A is
%                                  supplied in AP.
%
%              UPLO = 'L' or 'l'   The lower triangular part of A is
%                                  supplied in AP.
%
%           Unchanged on exit.
%
%  N      - INTEGER.
%           On entry, N specifies the order of the matrix A.
%           N must be at least zero.
%           Unchanged on exit.
%
%  ALPHA  - REAL            .
%           On entry, ALPHA specifies the scalar alpha.
%           Unchanged on exit.
%
%  X      - COMPLEX          array of dimension at least
%           ( 1 + ( n - 1 )*abs( INCX ) ).
%           Before entry, the incremented array X must contain the n
%           element vector x.
%           Unchanged on exit.
%
%  INCX   - INTEGER.
%           On entry, INCX specifies the increment for the elements of
%           X. INCX must not be zero.
%           Unchanged on exit.
%
%  AP     - COMPLEX          array of DIMENSION at least
%           ( ( n*( n + 1 ) )/2 ).
%           Before entry with  UPLO = 'U' or 'u', the array AP must
%           contain the upper triangular part of the hermitian 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 hermitian 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.
%           Note that the imaginary parts of the diagonal elements need
%           not be set, they are assumed to be zero, and on exit they
%           are set to zero.
%
%***REFERENCES  Dongarra, J. J., Du Croz, J., Hammarling, S., and
%                 Hanson, R. J.  An extended set of Fortran basic linear
%                 algebra subprograms.  ACM TOMS, Vol. 14, No. 1,
%                 pp. 1-17, March 1988.
%***ROUTINES CALLED  LSAME, XERBLA
%***REVISION HISTORY  (YYMMDD)
%   861022  DATE WRITTEN
%   910605  Modified to meet SLATEC prologue standards.  Only comment
%           lines were modified.  (BKS)
%***end PROLOGUE  CHPR
%     .. Scalar Arguments ..
%     .. Array Arguments ..
persistent i info ix j jx k kk kx temp zero ; 

ap_shape=size(ap);ap=reshape(ap,1,[]);
x_shape=size(x);x=reshape(x,1,[]);
%     .. Parameters ..
if isempty(zero), zero=complex(0.0e+0,0.0e+0) ; end;
%     .. Local Scalars ..
if isempty(temp), temp=0; end;
if isempty(i), i=0; end;
if isempty(info), info=0; end;
if isempty(ix), ix=0; end;
if isempty(j), j=0; end;
if isempty(jx), jx=0; end;
if isempty(k), k=0; end;
if isempty(kk), kk=0; end;
if isempty(kx), kx=0; end;
%     .. External Functions ..
%     .. External Subroutines ..
%     .. Intrinsic Functions ..
% intrinsic conjg , real ::;
%***FIRST EXECUTABLE STATEMENT  CHPR
%
%     Test the input parameters.
%
info = 0;
if( ~lsame(uplo,'U') && ~lsame(uplo,'L') )
info = 1;
elseif( n<0 ) ;
info = 2;
elseif( incx==0 ) ;
info = 5;
end;
if( info~=0 )
[dumvar1,info]=xerbla('CHPR  ',info);
ap_shape=zeros(ap_shape);ap_shape(:)=ap(1:numel(ap_shape));ap=ap_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
return;
end;
%
%     Quick return if possible.
%
if((n==0) ||(alpha==real(zero)) )
ap_shape=zeros(ap_shape);ap_shape(:)=ap(1:numel(ap_shape));ap=ap_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
return;
end;
%
%     Set the start point in X if the increment is not unity.
%
if( incx<=0 )
kx = fix(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( lsame(uplo,'U') )
%
%        Form  A  when upper triangle is stored in AP.
%
if( incx==1 )
for j = 1 : n;
if( x(j)~=zero )
temp = alpha.*conj(x(j));
k = fix(kk);
for i = 1 : j - 1;
ap(k) = ap(k) + x(i).*temp;
k = fix(k + 1);
end; i = fix(j - 1+1);
ap(kk+j-1) = real(ap(kk+j-1)) + real(x(j).*temp);
else;
ap(kk+j-1) = real(ap(kk+j-1));
end;
kk = fix(kk + j);
end; j = fix(n+1);
else;
jx = fix(kx);
for j = 1 : n;
if( x(jx)~=zero )
temp = alpha.*conj(x(jx));
ix = fix(kx);
for k = kk : kk + j - 2;
ap(k) = ap(k) + x(ix).*temp;
ix = fix(ix + incx);
end; k = fix(kk + j - 2+1);
ap(kk+j-1) = real(ap(kk+j-1)) + real(x(jx).*temp);
else;
ap(kk+j-1) = real(ap(kk+j-1));
end;
jx = fix(jx + incx);
kk = fix(kk + j);
end; j = fix(n+1);
end;
%
%        Form  A  when lower triangle is stored in AP.
%
elseif( incx==1 ) ;
for j = 1 : n;
if( x(j)~=zero )
temp = alpha.*conj(x(j));
ap(kk) = real(ap(kk)) + real(temp.*x(j));
k = fix(kk + 1);
for i = j + 1 : n;
ap(k) = ap(k) + x(i).*temp;
k = fix(k + 1);
end; i = fix(n+1);
else;
ap(kk) = real(ap(kk));
end;
kk = fix(kk + n - j + 1);
end; j = fix(n+1);
else;
jx = fix(kx);
for j = 1 : n;
if( x(jx)~=zero )
temp = alpha.*conj(x(jx));
ap(kk) = real(ap(kk)) + real(temp.*x(jx));
ix = fix(jx);
for k = kk + 1 : kk + n - j;
ix = fix(ix + incx);
ap(k) = ap(k) + x(ix).*temp;
end; k = fix(kk + n - j+1);
else;
ap(kk) = real(ap(kk));
end;
jx = fix(jx + incx);
kk = fix(kk + n - j + 1);
end; j = fix(n+1);
end;
%
%
%     end of CHPR  .
%
ap_shape=zeros(ap_shape);ap_shape(:)=ap(1:numel(ap_shape));ap=ap_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
end
%DECK CHPSL