[side,uplo,m,n,alpha,a,lda,b,ldb,beta,c,ldc]=csymm(side,uplo,m,n,alpha,a,lda,b,ldb,beta,c,ldc);

function [side,uplo,m,n,alpha,a,lda,b,ldb,beta,c,ldc]=csymm(side,uplo,m,n,alpha,a,lda,b,ldb,beta,c,ldc);
%***BEGIN PROLOGUE  CSYMM
%***PURPOSE  Multiply a complex general matrix by a complex symmetric
%            matrix.
%***LIBRARY   SLATEC (BLAS)
%***CATEGORY  D1B6
%***TYPE      COMPLEX (SSYMM-S, DSYMM-D, CSYMM-C)
%***KEYWORDS  LEVEL 3 BLAS, LINEAR ALGEBRA
%***AUTHOR  Dongarra, J., (ANL)
%           Duff, I., (AERE)
%           Du Croz, J., (NAG)
%           Hammarling, S. (NAG)
%***DESCRIPTION
%
%  CSYMM  performs one of the matrix-matrix operations
%
%     C := alpha*A*B + beta*C,
%
%  or
%
%     C := alpha*B*A + beta*C,
%
%  where  alpha and beta are scalars, A is a symmetric matrix and  B and
%  C are m by n matrices.
%
%  Parameters
%  ==========
%
%  SIDE   - CHARACTER*1.
%           On entry,  SIDE  specifies whether  the  symmetric matrix  A
%           appears on the  left or right  in the  operation as follows:
%
%              SIDE = 'L' or 'l'   C := alpha*A*B + beta*C,
%
%              SIDE = 'R' or 'r'   C := alpha*B*A + beta*C,
%
%           Unchanged on exit.
%
%  UPLO   - CHARACTER*1.
%           On  entry,   UPLO  specifies  whether  the  uppermlv  or  lower
%           triangular  part  of  the  symmetric  matrix   A  is  to  be
%           referenced as follows:
%
%              UPLO = 'U' or 'u'   Only the uppermlv triangular part of the
%                                  symmetric matrix is to be referenced.
%
%              UPLO = 'L' or 'l'   Only the lower triangular part of the
%                                  symmetric matrix is to be referenced.
%
%           Unchanged on exit.
%
%  M      - INTEGER.
%           On entry,  M  specifies the number of rows of the matrix  C.
%           M  must be at least zero.
%           Unchanged on exit.
%
%  N      - INTEGER.
%           On entry, N specifies the number of columns of the matrix C.
%           N  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
%           m  when  SIDE = 'L' or 'l'  and is n  otherwise.
%           Before entry  with  SIDE = 'L' or 'l',  the  m by m  part of
%           the array  A  must contain the  symmetric matrix,  such that
%           when  UPLO = 'U' or 'u', the leading m by m uppermlv triangular
%           part of the array  A  must contain the uppermlv triangular part
%           of the  symmetric matrix and the  strictly  lower triangular
%           part of  A  is not referenced,  and when  UPLO = 'L' or 'l',
%           the leading  m by m  lower triangular part  of the  array  A
%           must  contain  the  lower triangular part  of the  symmetric
%           matrix and the  strictly uppermlv triangular part of  A  is not
%           referenced.
%           Before entry  with  SIDE = 'R' or 'r',  the  n by n  part of
%           the array  A  must contain the  symmetric matrix,  such that
%           when  UPLO = 'U' or 'u', the leading n by n uppermlv triangular
%           part of the array  A  must contain the uppermlv triangular part
%           of the  symmetric matrix and the  strictly  lower triangular
%           part of  A  is not referenced,  and when  UPLO = 'L' or 'l',
%           the leading  n by n  lower triangular part  of the  array  A
%           must  contain  the  lower triangular part  of the  symmetric
%           matrix and the  strictly uppermlv triangular part of  A  is not
%           referenced.
%           Unchanged on exit.
%
%  LDA    - INTEGER.
%           On entry, LDA specifies the first dimension of A as declared
%           in the  calling (sub) program. When  SIDE = 'L' or 'l'  then
%           LDA must be at least  max( 1, m ), otherwise  LDA must be at
%           least max( 1, n ).
%           Unchanged on exit.
%
%  B      - COMPLEX          array of DIMENSION ( LDB, n ).
%           Before entry, the leading  m by n part of the array  B  must
%           contain the matrix B.
%           Unchanged on exit.
%
%  LDB    - INTEGER.
%           On entry, LDB specifies the first dimension of B as declared
%           in  the  calling  (sub)  program.   LDB  must  be  at  least
%           max( 1, m ).
%           Unchanged on exit.
%
%  BETA   - COMPLEX         .
%           On entry,  BETA  specifies the scalar  beta.  When  BETA  is
%           supplied as zero then C need not be set on input.
%           Unchanged on exit.
%
%  C      - COMPLEX          array of DIMENSION ( LDC, n ).
%           Before entry, the leading  m by n  part of the array  C must
%           contain the matrix  C,  except when  beta  is zero, in which
%           case C need not be set on entry.
%           On exit, the array  C  is overwritten by the  m by n 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, m ).
%           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  CSYMM
%     .. Scalar Arguments ..
%     .. Array Arguments ..
persistent i info j k nrowa one temp1 temp2 uppermlv zero ; 

a_shape=size(a);a=reshape([a(:).',zeros(1,ceil(numel(a)./prod([lda])).*prod([lda])-numel(a))],lda,[]);
b_shape=size(b);b=reshape([b(:).',zeros(1,ceil(numel(b)./prod([ldb])).*prod([ldb])-numel(b))],ldb,[]);
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(k), k=0; end;
if isempty(nrowa), nrowa=0; end;
if isempty(temp1), temp1=0; end;
if isempty(temp2), temp2=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  CSYMM
%
%     Set NROWA as the number of rows of A.
%
if( lsame(side,'L') )
nrowa = fix(m);
else;
nrowa = fix(n);
end;
[uppermlv ,uplo]=lsame(uplo,'U');
%
%     Test the input parameters.
%
info = 0;
if( (~lsame(side,'L')) && (~lsame(side,'R')) )
info = 1;
elseif ( (~uppermlv) && (~lsame(uplo,'L')) ) ;
info = 2;
elseif( m<0 ) ;
info = 3;
elseif( n<0 ) ;
info = 4;
elseif( lda<max(1,nrowa) ) ;
info = 7;
elseif( ldb<max(1,m) ) ;
info = 9;
elseif( ldc<max(1,m) ) ;
info = 12;
end;
if( info~=0 )
[dumvar1,info]=xerbla('CSYMM ',info);
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
return;
end;
%
%     Quick return if possible.
%
if((m==0) ||(n==0) ||((alpha==zero) &&(beta==one)) )
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_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( beta==zero )
for j = 1 : n;
for i = 1 : m;
c(i,j) = zero;
end; i = fix(m+1);
end; j = fix(n+1);
else;
for j = 1 : n;
for i = 1 : m;
c(i,j) = beta.*c(i,j);
end; i = fix(m+1);
end; j = fix(n+1);
end;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
return;
end;
%
%     Start the operations.
%
if( ~(lsame(side,'L')) )
%
%        Form  C := alpha*B*A + beta*C.
%
for j = 1 : n;
temp1 = alpha.*a(j,j);
if( beta==zero )
for i = 1 : m;
c(i,j) = temp1.*b(i,j);
end; i = fix(m+1);
else;
for i = 1 : m;
c(i,j) = beta.*c(i,j) + temp1.*b(i,j);
end; i = fix(m+1);
end;
for k = 1 : j - 1;
if( uppermlv )
temp1 = alpha.*a(k,j);
else;
temp1 = alpha.*a(j,k);
end;
for i = 1 : m;
c(i,j) = c(i,j) + temp1.*b(i,k);
end; i = fix(m+1);
end; k = fix(j - 1+1);
for k = j + 1 : n;
if( uppermlv )
temp1 = alpha.*a(j,k);
else;
temp1 = alpha.*a(k,j);
end;
for i = 1 : m;
c(i,j) = c(i,j) + temp1.*b(i,k);
end; i = fix(m+1);
end; k = fix(n+1);
end; j = fix(n+1);
%
%        Form  C := alpha*A*B + beta*C.
%
elseif( uppermlv ) ;
for j = 1 : n;
for i = 1 : m;
temp1 = alpha.*b(i,j);
temp2 = zero;
for k = 1 : i - 1;
c(k,j) = c(k,j) + temp1.*a(k,i);
temp2 = temp2 + b(k,j).*a(k,i);
end; k = fix(i - 1+1);
if( beta==zero )
c(i,j) = temp1.*a(i,i) + alpha.*temp2;
else;
c(i,j) = beta.*c(i,j) + temp1.*a(i,i) + alpha.*temp2;
end;
end; i = fix(m+1);
end; j = fix(n+1);
else;
for j = 1 : n;
for i = m : -1: 1 ;
temp1 = alpha.*b(i,j);
temp2 = zero;
for k = i + 1 : m;
c(k,j) = c(k,j) + temp1.*a(k,i);
temp2 = temp2 + b(k,j).*a(k,i);
end; k = fix(m+1);
if( beta==zero )
c(i,j) = temp1.*a(i,i) + alpha.*temp2;
else;
c(i,j) = beta.*c(i,j) + temp1.*a(i,i) + alpha.*temp2;
end;
end; i = fix(1 -1);
end; j = fix(n+1);
end;
%
%
%     end of CSYMM .
%
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
end
%DECK CSYR2K