[uplo,trans,n,k,alpha,a,lda,beta,c,ldc]=csyrk(uplo,trans,n,k,alpha,a,lda,beta,c,ldc);

function [uplo,trans,n,k,alpha,a,lda,beta,c,ldc]=csyrk(uplo,trans,n,k,alpha,a,lda,beta,c,ldc);
%***BEGIN PROLOGUE  CSYRK
%***PURPOSE  Perform symmetric rank k update of a complex symmetric
%            matrix.
%***LIBRARY   SLATEC (BLAS)
%***CATEGORY  D1B6
%***TYPE      COMPLEX (SSYRK-S, DSYRK-D, CSYRK-C)
%***KEYWORDS  LEVEL 3 BLAS, LINEAR ALGEBRA
%***AUTHOR  Dongarra, J., (ANL)
%           Duff, I., (AERE)
%           Du Croz, J., (NAG)
%           Hammarling, S. (NAG)
%***DESCRIPTION
%
%  CSYRK  performs one of the symmetric rank k operations
%
%     C := alpha*A*A' + beta*C,
%
%  or
%
%     C := alpha*A'*A + beta*C,
%
%  where  alpha and beta  are scalars,  C is an  n by n symmetric matrix
%  and  A  is an  n by k  matrix in the first case and a  k by n  matrix
%  in the second case.
%
%  Parameters
%  ==========
%
%  UPLO   - CHARACTER*1.
%           On  entry,   UPLO  specifies  whether  the  uppermlv  or  lower
%           triangular  part  of the  array  C  is to be  referenced  as
%           follows:
%
%              UPLO = 'U' or 'u'   Only the  uppermlv triangular part of  C
%                                  is to be referenced.
%
%              UPLO = 'L' or 'l'   Only the  lower triangular part of  C
%                                  is to be referenced.
%
%           Unchanged on exit.
%
%  TRANS  - CHARACTER*1.
%           On entry,  TRANS  specifies the operation to be performed as
%           follows:
%
%              TRANS = 'N' or 'n'   C := alpha*A*A' + beta*C.
%
%              TRANS = 'T' or 't'   C := alpha*A'*A + beta*C.
%
%           Unchanged on exit.
%
%  N      - INTEGER.
%           On entry,  N specifies the order of the matrix C.  N must be
%           at least zero.
%           Unchanged on exit.
%
%  K      - INTEGER.
%           On entry with  TRANS = 'N' or 'n',  K  specifies  the number
%           of  columns   of  the   matrix   A,   and  on   entry   with
%           TRANS = 'T' or 't',  K  specifies  the number of rows of the
%           matrix A.  K must be at least zero.
%           Unchanged on exit.
%
%  ALPHA  - COMPLEX         .
%           On entry, ALPHA specifies the scalar alpha.
%           Unchanged on exit.
%
%  A      - COMPLEX          array of DIMENSION ( LDA, ka ), where ka is
%           k  when  TRANS = 'N' or 'n',  and is  n  otherwise.
%           Before entry with  TRANS = 'N' or 'n',  the  leading  n by k
%           part of the array  A  must contain the matrix  A,  otherwise
%           the leading  k by n  part of the array  A  must contain  the
%           matrix A.
%           Unchanged on exit.
%
%  LDA    - INTEGER.
%           On entry, LDA specifies the first dimension of A as declared
%           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
%           then  LDA must be at least  max( 1, n ), otherwise  LDA must
%           be at least  max( 1, k ).
%           Unchanged on exit.
%
%  BETA   - COMPLEX         .
%           On entry, BETA specifies the scalar beta.
%           Unchanged on exit.
%
%  C      - COMPLEX          array of DIMENSION ( LDC, n ).
%           Before entry  with  UPLO = 'U' or 'u',  the leading  n by n
%           uppermlv triangular part of the array C must contain the uppermlv
%           triangular part  of the  symmetric matrix  and the strictly
%           lower triangular part of C is not referenced.  On exit, the
%           uppermlv triangular part of the array  C is overwritten by the
%           uppermlv triangular part of the updated matrix.
%           Before entry  with  UPLO = 'L' or 'l',  the leading  n by n
%           lower triangular part of the array C must contain the lower
%           triangular part  of the  symmetric matrix  and the strictly
%           uppermlv triangular part of C is not referenced.  On exit, the
%           lower triangular part of the array  C is overwritten by the
%           lower triangular part of the updated matrix.
%
%  LDC    - INTEGER.
%           On entry, LDC specifies the first dimension of C as declared
%           in  the  calling  (sub)  program.   LDC  must  be  at  least
%           max( 1, n ).
%           Unchanged on exit.
%
%***REFERENCES  Dongarra, J., Du Croz, J., Duff, I., and Hammarling, S.
%                 A set of level 3 basic linear algebra subprograms.
%                 ACM TOMS, Vol. 16, No. 1, pp. 1-17, March 1990.
%***ROUTINES CALLED  LSAME, XERBLA
%***REVISION HISTORY  (YYMMDD)
%   890208  DATE WRITTEN
%   910605  Modified to meet SLATEC prologue standards.  Only comment
%           lines were modified.  (BKS)
%***end PROLOGUE  CSYRK
%     .. Scalar Arguments ..
%     .. Array Arguments ..
persistent i info j l nrowa one temp uppermlv zero ; 

a_shape=size(a);a=reshape([a(:).',zeros(1,ceil(numel(a)./prod([lda])).*prod([lda])-numel(a))],lda,[]);
c_shape=size(c);c=reshape([c(:).',zeros(1,ceil(numel(c)./prod([ldc])).*prod([ldc])-numel(c))],ldc,[]);
%     .. External Functions ..
%     .. External Subroutines ..
%     .. Intrinsic Functions ..
%     .. Local Scalars ..
if isempty(uppermlv), uppermlv=false; end;
if isempty(i), i=0; end;
if isempty(info), info=0; end;
if isempty(j), j=0; end;
if isempty(l), l=0; end;
if isempty(nrowa), nrowa=0; end;
if isempty(temp), temp=0; end;
%     .. Parameters ..
if isempty(one), one=complex(1.0e+0,0.0e+0) ; end;
if isempty(zero), zero=complex(0.0e+0,0.0e+0) ; end;
%***FIRST EXECUTABLE STATEMENT  CSYRK
%
%     Test the input parameters.
%
if( lsame(trans,'N') )
nrowa = fix(n);
else;
nrowa = fix(k);
end;
[uppermlv ,uplo]=lsame(uplo,'U');
%
info = 0;
if( (~uppermlv) && (~lsame(uplo,'L')) )
info = 1;
elseif ( (~lsame(trans,'N')) && (~lsame(trans,'T')) );
info = 2;
elseif( n<0 ) ;
info = 3;
elseif( k<0 ) ;
info = 4;
elseif( lda<max(1,nrowa) ) ;
info = 7;
elseif( ldc<max(1,n) ) ;
info = 10;
end;
if( info~=0 )
[dumvar1,info]=xerbla('CSYRK ',info);
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
return;
end;
%
%     Quick return if possible.
%
if((n==0) ||(((alpha==zero) ||(k==0)) &&(beta==one)) )
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
return;
end;
%
%     And when  alpha==zero.
%
if( alpha==zero )
if( uppermlv )
if( beta==zero )
for j = 1 : n;
for i = 1 : j;
c(i,j) = zero;
end; i = fix(j+1);
end; j = fix(n+1);
else;
for j = 1 : n;
for i = 1 : j;
c(i,j) = beta.*c(i,j);
end; i = fix(j+1);
end; j = fix(n+1);
end;
elseif( beta==zero ) ;
for j = 1 : n;
for i = j : n;
c(i,j) = zero;
end; i = fix(n+1);
end; j = fix(n+1);
else;
for j = 1 : n;
for i = j : n;
c(i,j) = beta.*c(i,j);
end; i = fix(n+1);
end; j = fix(n+1);
end;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
return;
end;
%
%     Start the operations.
%
if( lsame(trans,'N') )
%
%        Form  C := alpha*A*A' + beta*C.
%
if( uppermlv )
for j = 1 : n;
if( beta==zero )
for i = 1 : j;
c(i,j) = zero;
end; i = fix(j+1);
elseif( beta~=one ) ;
for i = 1 : j;
c(i,j) = beta.*c(i,j);
end; i = fix(j+1);
end;
for l = 1 : k;
if( a(j,l)~=zero )
temp = alpha.*a(j,l);
for i = 1 : j;
c(i,j) = c(i,j) + temp.*a(i,l);
end; i = fix(j+1);
end;
end; l = fix(k+1);
end; j = fix(n+1);
else;
for j = 1 : n;
if( beta==zero )
for i = j : n;
c(i,j) = zero;
end; i = fix(n+1);
elseif( beta~=one ) ;
for i = j : n;
c(i,j) = beta.*c(i,j);
end; i = fix(n+1);
end;
for l = 1 : k;
if( a(j,l)~=zero )
temp = alpha.*a(j,l);
for i = j : n;
c(i,j) = c(i,j) + temp.*a(i,l);
end; i = fix(n+1);
end;
end; l = fix(k+1);
end; j = fix(n+1);
end;
%
%        Form  C := alpha*A'*A + beta*C.
%
elseif( uppermlv ) ;
for j = 1 : n;
for i = 1 : j;
temp = zero;
for l = 1 : k;
temp = temp + a(l,i).*a(l,j);
end; l = fix(k+1);
if( beta==zero )
c(i,j) = alpha.*temp;
else;
c(i,j) = alpha.*temp + beta.*c(i,j);
end;
end; i = fix(j+1);
end; j = fix(n+1);
else;
for j = 1 : n;
for i = j : n;
temp = zero;
for l = 1 : k;
temp = temp + a(l,i).*a(l,j);
end; l = fix(k+1);
if( beta==zero )
c(i,j) = alpha.*temp;
else;
c(i,j) = alpha.*temp + beta.*c(i,j);
end;
end; i = fix(n+1);
end; j = fix(n+1);
end;
%
%
%     end of CSYRK .
%
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
end
%DECK CTAN