| [uplo,trans,n,k,alpha,a,lda,b,ldb,beta,c,ldc]=dsyr2k(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]=dsyr2k(uplo,trans,n,k,alpha,a,lda,b,ldb,beta,c,ldc);
%***BEGIN PROLOGUE DSYR2K
%***PURPOSE Perform one of the symmetric rank 2k operations.
%***LIBRARY SLATEC (BLAS)
%***CATEGORY D1B6
%***TYPE doubleprecision (SSYR2-S, DSYR2-D, CSYR2-C, DSYR2K-D)
%***KEYWORDS LEVEL 3 BLAS, LINEAR ALGEBRA
%***AUTHOR Dongarra, J., (ANL)
% Duff, I., (AERE)
% Du Croz, J., (NAG)
% Hammarling, S. (NAG)
%***DESCRIPTION
%
% DSYR2K performs one of the symmetric rank 2k operations
%
% C := alpha*A*B' + alpha*B*A' + beta*C,
%
% or
%
% C := alpha*A'*B + alpha*B'*A + beta*C,
%
% where alpha and beta are scalars, C is an n by n symmetric 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*B' + alpha*B*A' +
% beta*C.
%
% TRANS = 'T' or 't' C := alpha*A'*B + alpha*B'*A +
% beta*C.
%
% TRANS = 'C' or 'c' C := alpha*A'*B + alpha*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 = 'T' or 't' or '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 - doubleprecision.
% On entry, ALPHA specifies the scalar alpha.
% Unchanged on exit.
%
% A - doubleprecision 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 - doubleprecision 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 - doubleprecision.
% On entry, BETA specifies the scalar beta.
% Unchanged on exit.
%
% C - doubleprecision 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 DSYR2K
% .. 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 ..
% .. 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.0d+0; end;
if isempty(zero), zero=0.0d+0 ; end;
%***FIRST EXECUTABLE STATEMENT DSYR2K
%
% 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'))&& (~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('DSYR2K',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==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;
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*B' + alpha*B*A' + 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) ||(b(j,l)~=zero) )
temp1 = alpha.*b(j,l);
temp2 = alpha.*a(j,l);
for i = 1 : j;
c(i,j) = c(i,j) + a(i,l).*temp1 + b(i,l).*temp2;
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) ||(b(j,l)~=zero) )
temp1 = alpha.*b(j,l);
temp2 = alpha.*a(j,l);
for i = j : n;
c(i,j) = c(i,j) + a(i,l).*temp1 + b(i,l).*temp2;
end; i = fix(n+1);
end;
end; l = fix(k+1);
end; j = fix(n+1);
end;
%
% Form C := alpha*A'*B + alpha*B'*A + C.
%
elseif( uppermlv ) ;
for j = 1 : n;
for i = 1 : j;
temp1 = zero;
temp2 = zero;
for l = 1 : k;
temp1 = temp1 + a(l,i).*b(l,j);
temp2 = temp2 + b(l,i).*a(l,j);
end; l = fix(k+1);
if( beta==zero )
c(i,j) = alpha.*temp1 + alpha.*temp2;
else;
c(i,j) = beta.*c(i,j) + alpha.*temp1 + 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 + a(l,i).*b(l,j);
temp2 = temp2 + b(l,i).*a(l,j);
end; l = fix(k+1);
if( beta==zero )
c(i,j) = alpha.*temp1 + alpha.*temp2;
else;
c(i,j) = beta.*c(i,j) + alpha.*temp1 + alpha.*temp2;
end;
end; i = fix(n+1);
end; j = fix(n+1);
end;
%
%
% end of DSYR2K.
%
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 DSYR
|
|
