| [uplo,trans,diag,n,k,a,lda,x,incx]=stbmv(uplo,trans,diag,n,k,a,lda,x,incx); |
function [uplo,trans,diag,n,k,a,lda,x,incx]=stbmv(uplo,trans,diag,n,k,a,lda,x,incx);
%***BEGIN PROLOGUE STBMV
%***PURPOSE Multiply a real vector by a real triangular band matrix.
%***LIBRARY SLATEC (BLAS)
%***CATEGORY D1B4
%***TYPE SINGLE PRECISION (STBMV-S, DTBMV-D, CTBMV-C)
%***KEYWORDS LEVEL 2 BLAS, LINEAR ALGEBRA
%***AUTHOR Dongarra, J. J., (ANL)
% Du Croz, J., (NAG)
% Hammarling, S., (NAG)
% Hanson, R. J., (SNLA)
%***DESCRIPTION
%
% STBMV performs one of the matrix-vector operations
%
% x := A*x, or x := A'*x,
%
% where x is an n element vector and A is an n by n unit, or non-unit,
% upper or lower triangular band matrix, with ( k + 1) diagonals.
%
% Parameters
% ==========
%
% UPLO - CHARACTER*1.
% On entry, UPLO specifies whether the matrix is an upper or
% lower triangular matrix as follows:
%
% UPLO = 'U' or 'u' A is an upper triangular matrix.
%
% UPLO = 'L' or 'l' A is a lower triangular matrix.
%
% Unchanged on exit.
%
% TRANS - CHARACTER*1.
% On entry, TRANS specifies the operation to be performed as
% follows:
%
% TRANS = 'N' or 'n' x := A*x.
%
% TRANS = 'T' or 't' x := A'*x.
%
% TRANS = 'C' or 'c' x := A'*x.
%
% Unchanged on exit.
%
% DIAG - CHARACTER*1.
% On entry, DIAG specifies whether or not A is unit
% triangular as follows:
%
% DIAG = 'U' or 'u' A is assumed to be unit triangular.
%
% DIAG = 'N' or 'n' A is not assumed to be unit
% triangular.
%
% Unchanged on exit.
%
% N - INTEGER.
% On entry, N specifies the order of the matrix A.
% N must be at least zero.
% Unchanged on exit.
%
% K - INTEGER.
% On entry with UPLO = 'U' or 'u', K specifies the number of
% super-diagonals of the matrix A.
% On entry with UPLO = 'L' or 'l', K specifies the number of
% sub-diagonals of the matrix A.
% K must satisfy 0 <= K.
% Unchanged on exit.
%
% A - REAL array of DIMENSION ( LDA, n ).
% Before entry with UPLO = 'U' or 'u', the leading ( k + 1 )
% by n part of the array A must contain the upper triangular
% band part of the matrix of coefficients, supplied column by
% column, with the leading diagonal of the matrix in row
% ( k + 1 ) of the array, the first super-diagonal starting at
% position 2 in row k, and so on. The top left k by k triangle
% of the array A is not referenced.
% The following program segment will transfer an upper
% triangular band matrix from conventional full matrix storage
% to band storage:
%
% DO 20, J = 1, N
% M = K + 1 - J
% DO 10, I = MAX( 1, J - K ), J
% A( M + I, J ) = matrix( I, J )
% 10 CONTINUE
% 20 CONTINUE
%
% Before entry with UPLO = 'L' or 'l', the leading ( k + 1 )
% by n part of the array A must contain the lower triangular
% band part of the matrix of coefficients, supplied column by
% column, with the leading diagonal of the matrix in row 1 of
% the array, the first sub-diagonal starting at position 1 in
% row 2, and so on. The bottom right k by k triangle of the
% array A is not referenced.
% The following program segment will transfer a lower
% triangular band matrix from conventional full matrix storage
% to band storage:
%
% DO 20, J = 1, N
% M = 1 - J
% DO 10, I = J, MIN( N, J + K )
% A( M + I, J ) = matrix( I, J )
% 10 CONTINUE
% 20 CONTINUE
%
% Note that when DIAG = 'U' or 'u' the elements of the array A
% corresponding to the diagonal elements of the matrix are not
% referenced, but are assumed to be unity.
% Unchanged on exit.
%
% LDA - INTEGER.
% On entry, LDA specifies the first dimension of A as declared
% in the calling (sub) program. LDA must be at least
% ( k + 1 ).
% Unchanged on exit.
%
% X - REAL array of dimension at least
% ( 1 + ( n - 1 )*abs( INCX ) ).
% Before entry, the incremented array X must contain the n
% element vector x. On exit, X is overwritten with the
% transformed vector x.
%
% INCX - INTEGER.
% On entry, INCX specifies the increment for the elements of
% X. INCX must not be zero.
% 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 STBMV
% .. Scalar Arguments ..
% .. Array Arguments ..
persistent i info ix j jx kplus1 kx l nounit 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=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(kplus1), kplus1=0; end;
if isempty(kx), kx=0; end;
if isempty(l), l=0; end;
if isempty(nounit), nounit=false; end;
% .. External Functions ..
% .. External Subroutines ..
% .. Intrinsic Functions ..
%***FIRST EXECUTABLE STATEMENT STBMV
%
% Test the input parameters.
%
info = 0;
if( ~lsame(uplo,'U') && ~lsame(uplo,'L') )
info = 1;
elseif ( ~lsame(trans,'N') && ~lsame(trans,'T') &&~lsame(trans,'C') ) ;
info = 2;
elseif ( ~lsame(diag,'U') && ~lsame(diag,'N') ) ;
info = 3;
elseif( n<0 ) ;
info = 4;
elseif( k<0 ) ;
info = 5;
elseif( lda<(k+1) ) ;
info = 7;
elseif( incx==0 ) ;
info = 9;
end;
if( info~=0 )
[dumvar1,info]=xerbla('STBMV ',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 )
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;
%
[nounit ,diag]=lsame(diag,'N');
%
% Set up the start point in X if the increment is not unity. This
% will be ( N - 1 )*INCX too small for descending loops.
%
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 A.
%
if( lsame(trans,'N') )
%
% Form x := A*x.
%
if( lsame(uplo,'U') )
kplus1 = fix(k + 1);
if( incx==1 )
for j = 1 : n;
if( x(j)~=zero )
temp = x(j);
l = fix(kplus1 - j);
for i = max(1,j-k) : j - 1;
x(i) = x(i) + temp.*a(l+i,j);
end; i = fix(j - 1+1);
if( nounit )
x(j) = x(j).*a(kplus1,j);
end;
end;
end; j = fix(n+1);
else;
jx = fix(kx);
for j = 1 : n;
if( x(jx)~=zero )
temp = x(jx);
ix = fix(kx);
l = fix(kplus1 - j);
for i = max(1,j-k) : j - 1;
x(ix) = x(ix) + temp.*a(l+i,j);
ix = fix(ix + incx);
end; i = fix(j - 1+1);
if( nounit )
x(jx) = x(jx).*a(kplus1,j);
end;
end;
jx = fix(jx + incx);
if( j>k )
kx = fix(kx + incx);
end;
end; j = fix(n+1);
end;
elseif( incx==1 ) ;
for j = n : -1: 1 ;
if( x(j)~=zero )
temp = x(j);
l = fix(1 - j);
for i = min(n,j+k) : -1: j + 1 ;
x(i) = x(i) + temp.*a(l+i,j);
end; i = fix(j + 1 -1);
if( nounit )
x(j) = x(j).*a(1,j);
end;
end;
end; j = fix(1 -1);
else;
kx = fix(kx +(n-1).*incx);
jx = fix(kx);
for j = n : -1: 1 ;
if( x(jx)~=zero )
temp = x(jx);
ix = fix(kx);
l = fix(1 - j);
for i = min(n,j+k) : -1: j + 1 ;
x(ix) = x(ix) + temp.*a(l+i,j);
ix = fix(ix - incx);
end; i = fix(j + 1 -1);
if( nounit )
x(jx) = x(jx).*a(1,j);
end;
end;
jx = fix(jx - incx);
if((n-j)>=k )
kx = fix(kx - incx);
end;
end; j = fix(1 -1);
end;
%
% Form x := A'*x.
%
elseif ( lsame(uplo,'U') ) ;
kplus1 = fix(k + 1);
if( incx==1 )
for j = n : -1: 1 ;
temp = x(j);
l = fix(kplus1 - j);
if( nounit )
temp = temp.*a(kplus1,j);
end;
for i = j - 1 : -1: max(1,j-k) ;
temp = temp + a(l+i,j).*x(i);
end; i = fix(max(1,j-k) -1);
x(j) = temp;
end; j = fix(1 -1);
else;
kx = fix(kx +(n-1).*incx);
jx = fix(kx);
for j = n : -1: 1 ;
temp = x(jx);
kx = fix(kx - incx);
ix = fix(kx);
l = fix(kplus1 - j);
if( nounit )
temp = temp.*a(kplus1,j);
end;
for i = j - 1 : -1: max(1,j-k) ;
temp = temp + a(l+i,j).*x(ix);
ix = fix(ix - incx);
end; i = fix(max(1,j-k) -1);
x(jx) = temp;
jx = fix(jx - incx);
end; j = fix(1 -1);
end;
elseif( incx==1 ) ;
for j = 1 : n;
temp = x(j);
l = fix(1 - j);
if( nounit )
temp = temp.*a(1,j);
end;
for i = j + 1 : min(n,j+k);
temp = temp + a(l+i,j).*x(i);
end; i = fix(min(n,j+k)+1);
x(j) = temp;
end; j = fix(n+1);
else;
jx = fix(kx);
for j = 1 : n;
temp = x(jx);
kx = fix(kx + incx);
ix = fix(kx);
l = fix(1 - j);
if( nounit )
temp = temp.*a(1,j);
end;
for i = j + 1 : min(n,j+k);
temp = temp + a(l+i,j).*x(ix);
ix = fix(ix + incx);
end; i = fix(min(n,j+k)+1);
x(jx) = temp;
jx = fix(jx + incx);
end; j = fix(n+1);
end;
%
%
% end of STBMV .
%
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 STBSV
|
|
