| [uplo,trans,n,k,alpha,a,lda,b,ldb,beta,c,ldc]=cher2k(uplo,trans,n,k,alpha,a,lda,b,ldb,beta,c,ldc); |
function [uplo,trans,n,k,alpha,a,lda,b,ldb,beta,c,ldc]=cher2k(uplo,trans,n,k,alpha,a,lda,b,ldb,beta,c,ldc);
%***BEGIN PROLOGUE CHER2K
%***PURPOSE Perform Hermitian rank 2k update of a complex.
%***LIBRARY SLATEC (BLAS)
%***CATEGORY D1B6
%***TYPE COMPLEX (SHER2-S, DHER2-D, CHER2-C, CHER2K-C)
%***KEYWORDS LEVEL 3 BLAS, LINEAR ALGEBRA
%***AUTHOR Dongarra, J., (ANL)
% Duff, I., (AERE)
% Du Croz, J., (NAG)
% Hammarling, S. (NAG)
%***DESCRIPTION
%
% CHER2K performs one of the hermitian rank 2k operations
%
% C := alpha*A*conjg( B' ) + conjg( alpha )*B*conjg( A' ) + beta*C,
%
% or
%
% C := alpha*conjg( A' )*B + conjg( alpha )*conjg( B' )*A + beta*C,
%
% where alpha and beta are scalars with beta real, C is an n by n
% hermitian matrix and A and B are n by k matrices in the first case
% and k by n matrices 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*conjg( B' ) +
% conjg( alpha )*B*conjg( A' ) +
% beta*C.
%
% TRANS = 'C' or 'c' C := alpha*conjg( A' )*B +
% conjg( alpha )*conjg( B' )*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 matrices A and B, and on entry with
% TRANS = 'C' or 'c', K specifies the number of rows of the
% matrices A and B. 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.
%
% B - COMPLEX array of DIMENSION ( LDB, kb ), where kb 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 B must contain the matrix B, otherwise
% the leading k 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. When TRANS = 'N' or 'n'
% then LDB must be at least max( 1, n ), otherwise LDB must
% be at least max( 1, k ).
% Unchanged on exit.
%
% BETA - REAL .
% 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 hermitian 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 hermitian 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.
% 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.
%
% 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 CHER2K
% .. Scalar Arguments ..
% .. Array Arguments ..
persistent i info j l 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 ..
% intrinsic conjg , max , real ::;
% .. 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(temp1), temp1=0; end;
if isempty(temp2), temp2=0; end;
% .. Parameters ..
if isempty(one), one=1.0e+0 ; end;
if isempty(zero), zero=complex(0.0e+0,0.0e+0) ; end;
%***FIRST EXECUTABLE STATEMENT CHER2K
%
% 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,'C')) );
info = 2;
elseif( n<0 ) ;
info = 3;
elseif( k<0 ) ;
info = 4;
elseif( lda<max(1,nrowa) ) ;
info = 7;
elseif( ldb<max(1,nrowa) ) ;
info = 9;
elseif( ldc<max(1,n) ) ;
info = 12;
end;
if( info~=0 )
[dumvar1,info]=xerbla('CHER2K',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((n==0) ||(((alpha==zero) ||(k==0)) &&(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( uppermlv )
if( beta==real(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 - 1;
c(i,j) = beta.*c(i,j);
end; i = fix(j - 1+1);
c(j,j) = beta.*real(c(j,j));
end; j = fix(n+1);
end;
elseif( beta==real(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;
c(j,j) = beta.*real(c(j,j));
for i = j + 1 : 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;
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(trans,'N') )
%
% Form C := alpha*A*conjg( B' ) + conjg( alpha )*B*conjg( A' ) +
% C.
%
if( uppermlv )
for j = 1 : n;
if( beta==real(zero) )
for i = 1 : j;
c(i,j) = zero;
end; i = fix(j+1);
elseif( beta~=one ) ;
for i = 1 : j - 1;
c(i,j) = beta.*c(i,j);
end; i = fix(j - 1+1);
c(j,j) = beta.*real(c(j,j));
end;
for l = 1 : k;
if((a(j,l)~=zero) ||(b(j,l)~=zero) )
temp1 = alpha.*conj(b(j,l));
temp2 = conj(alpha.*a(j,l));
for i = 1 : j - 1;
c(i,j) = c(i,j) + a(i,l).*temp1 + b(i,l).*temp2;
end; i = fix(j - 1+1);
c(j,j) = real(c(j,j)) + real(a(j,l).*temp1+b(j,l).*temp2);
end;
end; l = fix(k+1);
end; j = fix(n+1);
else;
for j = 1 : n;
if( beta==real(zero) )
for i = j : n;
c(i,j) = zero;
end; i = fix(n+1);
elseif( beta~=one ) ;
for i = j + 1 : n;
c(i,j) = beta.*c(i,j);
end; i = fix(n+1);
c(j,j) = beta.*real(c(j,j));
end;
for l = 1 : k;
if((a(j,l)~=zero) ||(b(j,l)~=zero) )
temp1 = alpha.*conj(b(j,l));
temp2 = conj(alpha.*a(j,l));
for i = j + 1 : n;
c(i,j) = c(i,j) + a(i,l).*temp1 + b(i,l).*temp2;
end; i = fix(n+1);
c(j,j) = real(c(j,j)) + real(a(j,l).*temp1+b(j,l).*temp2);
end;
end; l = fix(k+1);
end; j = fix(n+1);
end;
%
% Form C := alpha*conjg( A' )*B + conjg( alpha )*conjg( B' )*A +
% C.
%
elseif( uppermlv ) ;
for j = 1 : n;
for i = 1 : j;
temp1 = zero;
temp2 = zero;
for l = 1 : k;
temp1 = temp1 + conj(a(l,i)).*b(l,j);
temp2 = temp2 + conj(b(l,i)).*a(l,j);
end; l = fix(k+1);
if( i==j )
if( beta==real(zero) )
c(j,j) = real(alpha.*temp1+conj(alpha).*temp2);
else;
c(j,j) = beta.*real(c(j,j))+ real(alpha.*temp1+conj(alpha).*temp2);
end;
elseif( beta==real(zero) ) ;
c(i,j) = alpha.*temp1 + conj(alpha).*temp2;
else;
c(i,j) = beta.*c(i,j) + alpha.*temp1 + conj(alpha).*temp2;
end;
end; i = fix(j+1);
end; j = fix(n+1);
else;
for j = 1 : n;
for i = j : n;
temp1 = zero;
temp2 = zero;
for l = 1 : k;
temp1 = temp1 + conj(a(l,i)).*b(l,j);
temp2 = temp2 + conj(b(l,i)).*a(l,j);
end; l = fix(k+1);
if( i==j )
if( beta==real(zero) )
c(j,j) = real(alpha.*temp1+conj(alpha).*temp2);
else;
c(j,j) = beta.*real(c(j,j))+ real(alpha.*temp1+conj(alpha).*temp2);
end;
elseif( beta==real(zero) ) ;
c(i,j) = alpha.*temp1 + conj(alpha).*temp2;
else;
c(i,j) = beta.*c(i,j) + alpha.*temp1 + conj(alpha).*temp2;
end;
end; i = fix(n+1);
end; j = fix(n+1);
end;
%
%
% end of CHER2K.
%
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 CHER
|
|
