| [uplo,n,alpha,x,incx,a,lda]=cher(uplo,n,alpha,x,incx,a,lda); |
function [uplo,n,alpha,x,incx,a,lda]=cher(uplo,n,alpha,x,incx,a,lda);
%***BEGIN PROLOGUE CHER
%***PURPOSE Perform Hermitian rank 1 update of a complex Hermitian
% matrix.
%***LIBRARY SLATEC (BLAS)
%***CATEGORY D1B4
%***TYPE COMPLEX (SHER-S, DHER-D, CHER-C)
%***KEYWORDS LEVEL 2 BLAS, LINEAR ALGEBRA
%***AUTHOR Dongarra, J. J., (ANL)
% Du Croz, J., (NAG)
% Hammarling, S., (NAG)
% Hanson, R. J., (SNLA)
%***DESCRIPTION
%
% CHER 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.
%
% Parameters
% ==========
%
% UPLO - CHARACTER*1.
% On entry, UPLO specifies whether the upper or lower
% triangular part of the array A is to be referenced as
% follows:
%
% UPLO = 'U' or 'u' Only the upper triangular part of A
% is to be referenced.
%
% UPLO = 'L' or 'l' Only the lower triangular part of A
% is to be referenced.
%
% 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.
%
% A - COMPLEX array of DIMENSION ( LDA, n ).
% Before entry with UPLO = 'U' or 'u', the leading n by n
% upper triangular part of the array A must contain the upper
% triangular part of the hermitian matrix and the strictly
% lower triangular part of A is not referenced. On exit, the
% upper triangular part of the array A is overwritten by the
% upper triangular part of the updated matrix.
% Before entry with UPLO = 'L' or 'l', the leading n by n
% lower triangular part of the array A must contain the lower
% triangular part of the hermitian matrix and the strictly
% upper triangular part of A is not referenced. On exit, the
% lower triangular part of the array A 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.
%
% LDA - INTEGER.
% On entry, LDA specifies the first dimension of A as declared
% in the calling (sub) program. LDA must be at least
% max( 1, n ).
% Unchanged on exit.
%
%***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 CHER
% .. Scalar Arguments ..
% .. Array Arguments ..
persistent i info ix j jx kx temp zero ;
a_shape=size(a);a=reshape([a(:).',zeros(1,ceil(numel(a)./prod([lda])).*prod([lda])-numel(a))],lda,[]);
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(kx), kx=0; end;
% .. External Functions ..
% .. External Subroutines ..
% .. Intrinsic Functions ..
% intrinsic conjg , max , real ::;
%***FIRST EXECUTABLE STATEMENT CHER
%
% 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;
elseif( lda<max(1,n) ) ;
info = 7;
end;
if( info~=0 )
[dumvar1,info]=xerbla('CHER ',info);
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_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)) )
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_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 A are
% accessed sequentially with one pass through the triangular part
% of A.
%
if( lsame(uplo,'U') )
%
% Form A when A is stored in upper triangle.
%
if( incx==1 )
for j = 1 : n;
if( x(j)~=zero )
temp = alpha.*conj(x(j));
for i = 1 : j - 1;
a(i,j) = a(i,j) + x(i).*temp;
end; i = fix(j - 1+1);
a(j,j) = real(a(j,j)) + real(x(j).*temp);
else;
a(j,j) = real(a(j,j));
end;
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 i = 1 : j - 1;
a(i,j) = a(i,j) + x(ix).*temp;
ix = fix(ix + incx);
end; i = fix(j - 1+1);
a(j,j) = real(a(j,j)) + real(x(jx).*temp);
else;
a(j,j) = real(a(j,j));
end;
jx = fix(jx + incx);
end; j = fix(n+1);
end;
%
% Form A when A is stored in lower triangle.
%
elseif( incx==1 ) ;
for j = 1 : n;
if( x(j)~=zero )
temp = alpha.*conj(x(j));
a(j,j) = real(a(j,j)) + real(temp.*x(j));
for i = j + 1 : n;
a(i,j) = a(i,j) + x(i).*temp;
end; i = fix(n+1);
else;
a(j,j) = real(a(j,j));
end;
end; j = fix(n+1);
else;
jx = fix(kx);
for j = 1 : n;
if( x(jx)~=zero )
temp = alpha.*conj(x(jx));
a(j,j) = real(a(j,j)) + real(temp.*x(jx));
ix = fix(jx);
for i = j + 1 : n;
ix = fix(ix + incx);
a(i,j) = a(i,j) + x(ix).*temp;
end; i = fix(n+1);
else;
a(j,j) = real(a(j,j));
end;
jx = fix(jx + incx);
end; j = fix(n+1);
end;
%
%
% end of CHER .
%
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
end
%DECK CHERK
|
|