| [abd,lda,n,ml,mu,ipvt,info]=sgbfa(abd,lda,n,ml,mu,ipvt,info); |
function [abd,lda,n,ml,mu,ipvt,info]=sgbfa(abd,lda,n,ml,mu,ipvt,info);
%***BEGIN PROLOGUE SGBFA
%***PURPOSE Factor a band matrix using Gaussian elimination.
%***LIBRARY SLATEC (LINPACK)
%***CATEGORY D2A2
%***TYPE SINGLE PRECISION (SGBFA-S, DGBFA-D, CGBFA-C)
%***KEYWORDS BANDED, LINEAR ALGEBRA, LINPACK, MATRIX FACTORIZATION
%***AUTHOR Moler, C. B., (U. of New Mexico)
%***DESCRIPTION
%
% SGBFA factors a real band matrix by elimination.
%
% SGBFA is usually called by SBGCO, but it can be called
% directly with a saving in time if RCOND is not needed.
%
% On Entry
%
% ABD REAL(LDA, N)
% contains the matrix in band storage. The columns
% of the matrix are stored in the columns of ABD and
% the diagonals of the matrix are stored in rows
% ML+1 through 2*ML+MU+1 of ABD .
% See the comments below for details.
%
% LDA INTEGER
% the leading dimension of the array ABD .
% LDA must be .GE. 2*ML + MU + 1 .
%
% N INTEGER
% the order of the original matrix.
%
% ML INTEGER
% number of diagonals below the main diagonal.
% 0 .LE. ML .LT. N .
%
% MU INTEGER
% number of diagonals above the main diagonal.
% 0 .LE. MU .LT. N .
% More efficient if ML .LE. MU .
% On Return
%
% ABD an upper triangular matrix in band storage and
% the multipliers which were used to obtain it.
% The factorization can be written A = L*U , where
% L is a product of permutation and unit lower
% triangular matrices and U is upper triangular.
%
% IPVT INTEGER(N)
% an integer vector of pivot indices.
%
% INFO INTEGER
% = 0 normal value.
% = K if U(K,K) .EQ. 0.0 . This is not an error
% condition for this subroutine, but it does
% indicate that SGBSL will divide by zero if
% called. Use RCOND in SBGCO for a reliable
% indication of singularity.
%
% Band Storage
%
% If A is a band matrix, the following program segment
% will set up the input.
%
% ML = (band width below the diagonal)
% MU = (band width above the diagonal)
% M = ML + MU + 1
% DO 20 J = 1, N
% I1 = MAX(1, J-MU)
% I2 = MIN(N, J+ML)
% DO 10 I = I1, I2
% K = I - J + M
% ABD(K,J) = A(I,J)
% 10 CONTINUE
% 20 CONTINUE
%
% This uses rows ML+1 through 2*ML+MU+1 of ABD .
% In addition, the first ML rows in ABD are used for
% elements generated during the triangularization.
% The total number of rows needed in ABD is 2*ML+MU+1 .
% The ML+MU by ML+MU upper left triangle and the
% ML by ML lower right triangle are not referenced.
%
%***REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
% Stewart, LINPACK Users' Guide, SIAM, 1979.
%***ROUTINES CALLED ISAMAX, SAXPY, SSCAL
%***REVISION HISTORY (YYMMDD)
% 780814 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (WRB)
% 890831 Modified array declarations. (WRB)
% 890831 REVISION DATE from Version 3.2
% 891214 Prologue converted to Version 4.0 format. (BAB)
% 900326 Removed duplicate information from DESCRIPTION section.
% (WRB)
% 920501 Reformatted the REFERENCES section. (WRB)
%***end PROLOGUE SGBFA
persistent i i0 j j0 j1 ju jz k kp1 l lm m mm nm1 t ;
ipvt_shape=size(ipvt);ipvt=reshape(ipvt,1,[]);
abd_shape=size(abd);abd=reshape([abd(:).',zeros(1,ceil(numel(abd)./prod([lda])).*prod([lda])-numel(abd))],lda,[]);
%
if isempty(t), t=0; end;
if isempty(i), i=0; end;
if isempty(i0), i0=0; end;
if isempty(j), j=0; end;
if isempty(ju), ju=0; end;
if isempty(jz), jz=0; end;
if isempty(j0), j0=0; end;
if isempty(j1), j1=0; end;
if isempty(k), k=0; end;
if isempty(kp1), kp1=0; end;
if isempty(l), l=0; end;
if isempty(lm), lm=0; end;
if isempty(m), m=0; end;
if isempty(mm), mm=0; end;
if isempty(nm1), nm1=0; end;
%
%***FIRST EXECUTABLE STATEMENT SGBFA
m = fix(ml + mu + 1);
info = 0;
%
% ZERO INITIAL FILL-IN COLUMNS
%
j0 = fix(mu + 2);
j1 = fix(min(n,m) - 1);
if( j1>=j0 )
for jz = j0 : j1;
i0 = fix(m + 1 - jz);
for i = i0 : ml;
abd(i,jz) = 0.0e0;
end; i = fix(ml+1);
end; jz = fix(j1+1);
end;
jz = fix(j1);
ju = 0;
%
% GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
%
nm1 = fix(n - 1);
if( nm1>=1 )
for k = 1 : nm1;
kp1 = fix(k + 1);
%
% ZERO NEXT FILL-IN COLUMN
%
jz = fix(jz + 1);
if( jz<=n )
if( ml>=1 )
for i = 1 : ml;
abd(i,jz) = 0.0e0;
end; i = fix(ml+1);
end;
end;
%
% FIND L = PIVOT INDEX
%
lm = fix(min(ml,n-k));
l = fix(isamax(lm+1,abd(sub2ind(size(abd),m,k):end),1) + m - 1);
ipvt(k) = fix(l + k - m);
%
% ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
%
if( abd(l,k)==0.0e0 )
info = fix(k);
else;
%
% INTERCHANGE IF NECESSARY
%
if( l~=m )
t = abd(l,k);
abd(l,k) = abd(m,k);
abd(m,k) = t;
end;
%
% COMPUTE MULTIPLIERS
%
t = -1.0e0./abd(m,k);
[lm,t,abd(sub2ind(size(abd),m+1,k):end)]=sscal(lm,t,abd(sub2ind(size(abd),m+1,k):end),1);
%
% ROW ELIMINATION WITH COLUMN INDEXING
%
ju = fix(min(max(ju,mu+ipvt(k)),n));
mm = fix(m);
if( ju>=kp1 )
for j = kp1 : ju;
l = fix(l - 1);
mm = fix(mm - 1);
t = abd(l,j);
if( l~=mm )
abd(l,j) = abd(mm,j);
abd(mm,j) = t;
end;
[lm,t,dumvar3,dumvar4,dumvar5]=saxpy(lm,t,abd(sub2ind(size(abd),m+1,k):end),1,abd(sub2ind(size(abd),mm+1,j):end),1); abd(sub2ind(size(abd),m+1,k):end)=dumvar3; abd(sub2ind(size(abd),mm+1,j):end)=dumvar5;
end; j = fix(ju+1);
end;
end;
end; k = fix(nm1+1);
end;
ipvt(n) = fix(n);
if( abd(m,n)==0.0e0 )
info = fix(n);
end;
ipvt_shape=zeros(ipvt_shape);ipvt_shape(:)=ipvt(1:numel(ipvt_shape));ipvt=ipvt_shape;
abd_shape=zeros(abd_shape);abd_shape(:)=abd(1:numel(abd_shape));abd=abd_shape;
end
%DECK SGBMV
|
|
