Code covered by the BSD License

Highlights fromslatec

from slatec by Ben Barrowes
The slatec library converted into matlab functions.

[abd,lda,n,ml,mu,ipvt,rcond,z]=sgbco(abd,lda,n,ml,mu,ipvt,rcond,z);
```function [abd,lda,n,ml,mu,ipvt,rcond,z]=sgbco(abd,lda,n,ml,mu,ipvt,rcond,z);
%***BEGIN PROLOGUE  SGBCO
%***PURPOSE  Factor a band matrix by Gaussian elimination and
%            estimate the condition number of the matrix.
%***LIBRARY   SLATEC (LINPACK)
%***CATEGORY  D2A2
%***TYPE      SINGLE PRECISION (SGBCO-S, DGBCO-D, CGBCO-C)
%***KEYWORDS  BANDED, CONDITION NUMBER, LINEAR ALGEBRA, LINPACK,
%             MATRIX FACTORIZATION
%***AUTHOR  Moler, C. B., (U. of New Mexico)
%***DESCRIPTION
%
%     SBGCO factors a real band matrix by Gaussian
%     elimination and estimates the condition of the matrix.
%
%     If  RCOND  is not needed, SGBFA is slightly faster.
%     To solve  A*X = B , follow SBGCO by SGBSL.
%     To compute  INVERSE(A)*C , follow SBGCO by SGBSL.
%     To compute  DETERMINANT(A) , follow SBGCO by SGBDI.
%
%     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.
%
%        RCOND   REAL
%                an estimate of the reciprocal condition of  A .
%                For the system  A*X = B , relative perturbations
%                in  A  and  B  of size  EPSILON  may cause
%                relative perturbations in  X  of size  EPSILON/RCOND .
%                If  RCOND  is so small that the logical expression
%                           1.0 + RCOND .EQ. 1.0
%                is truemlv, then  A  may be singular to working
%                precision.  In particular,  RCOND  is zero  if
%                exact singularity is detected or the estimate
%                underflows.
%
%        Z       REAL(N)
%                a work vector whose contents are usually unimportant.
%                If  A  is close to a singular matrix, then  Z  is
%                an approximate null vector in the sense that
%                NORM(A*Z) = RCOND*NORM(A)*NORM(Z) .
%
%     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.
%
%     Example:  If the original matrix is
%
%           11 12 13  0  0  0
%           21 22 23 24  0  0
%            0 32 33 34 35  0
%            0  0 43 44 45 46
%            0  0  0 54 55 56
%            0  0  0  0 65 66
%
%      then  N = 6, ML = 1, MU = 2, LDA .GE. 5  and ABD should contain
%
%            *  *  *  +  +  +  , * = not used
%            *  * 13 24 35 46  , + = used for pivoting
%            * 12 23 34 45 56
%           11 22 33 44 55 66
%           21 32 43 54 65  *
%
%***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
%                 Stewart, LINPACK Users' Guide, SIAM, 1979.
%***ROUTINES CALLED  SASUM, SAXPY, SDOT, SGBFA, 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  SGBCO
persistent anorm ek info is j ju k kb kp1 l la lm lz m mm s sm t wk wkm ynorm ;

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,[]);
z_shape=size(z);z=reshape(z,1,[]);
%
if isempty(ek), ek=0; end;
if isempty(t), t=0; end;
if isempty(wk), wk=0; end;
if isempty(wkm), wkm=0; end;
if isempty(anorm), anorm=0; end;
if isempty(s), s=0; end;
if isempty(sm), sm=0; end;
if isempty(ynorm), ynorm=0; end;
if isempty(is), is=0; end;
if isempty(info), info=0; end;
if isempty(j), j=0; end;
if isempty(ju), ju=0; end;
if isempty(k), k=0; end;
if isempty(kb), kb=0; end;
if isempty(kp1), kp1=0; end;
if isempty(l), l=0; end;
if isempty(la), la=0; end;
if isempty(lm), lm=0; end;
if isempty(lz), lz=0; end;
if isempty(m), m=0; end;
if isempty(mm), mm=0; end;
%
%     COMPUTE 1-NORM OF A
%
%***FIRST EXECUTABLE STATEMENT  SGBCO
anorm = 0.0e0;
l = fix(ml + 1);
is = fix(l + mu);
for j = 1 : n;
anorm = max(anorm,sasum(l,abd(sub2ind(size(abd),is,j):end),1));
if( is>ml+1 )
is = fix(is - 1);
end;
if( j<=mu )
l = fix(l + 1);
end;
if( j>=n-ml )
l = fix(l - 1);
end;
end; j = fix(n+1);
%
%     FACTOR
%
[abd,lda,n,ml,mu,ipvt,info]=sgbfa(abd,lda,n,ml,mu,ipvt,info);
%
%     RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) .
%     ESTIMATE = NORM(Z)/NORM(Y) WHERE  A*Z = Y  AND  TRANS(A)*Y = E .
%     TRANS(A)  IS THE TRANSPOSE OF A .  THE COMPONENTS OF  E  ARE
%     CHOSEN TO CAUSE MAXIMUM LOCAL GROWTH IN THE ELEMENTS OF W  WHERE
%     TRANS(U)*W = E .  THE VECTORS ARE FREQUENTLY RESCALED TO AVOID
%     OVERFLOW.
%
%     SOLVE TRANS(U)*W = E
%
ek = 1.0e0;
for j = 1 : n;
z(j) = 0.0e0;
end; j = fix(n+1);
m = fix(ml + mu + 1);
ju = 0;
for k = 1 : n;
if( z(k)~=0.0e0 )
ek = (abs(ek).*sign(-z(k)));
end;
if( abs(ek-z(k))>abs(abd(m,k)) )
s = abs(abd(m,k))./abs(ek-z(k));
[n,s,z]=sscal(n,s,z,1);
ek = s.*ek;
end;
wk = ek - z(k);
wkm = -ek - z(k);
s = abs(wk);
sm = abs(wkm);
if( abd(m,k)==0.0e0 )
wk = 1.0e0;
wkm = 1.0e0;
else;
wk = wk./abd(m,k);
wkm = wkm./abd(m,k);
end;
kp1 = fix(k + 1);
ju = fix(min(max(ju,mu+ipvt(k)),n));
mm = fix(m);
if( kp1<=ju )
for j = kp1 : ju;
mm = fix(mm - 1);
sm = sm + abs(z(j)+wkm.*abd(mm,j));
z(j) = z(j) + wk.*abd(mm,j);
s = s + abs(z(j));
end; j = fix(ju+1);
if( s<sm )
t = wkm - wk;
wk = wkm;
mm = fix(m);
for j = kp1 : ju;
mm = fix(mm - 1);
z(j) = z(j) + t.*abd(mm,j);
end; j = fix(ju+1);
end;
end;
z(k) = wk;
end; k = fix(n+1);
s = 1.0e0./sasum(n,z,1);
[n,s,z]=sscal(n,s,z,1);
%
%     SOLVE TRANS(L)*Y = W
%
for kb = 1 : n;
k = fix(n + 1 - kb);
lm = fix(min(ml,n-k));
if( k<n )
z(k) = z(k) + sdot(lm,abd(sub2ind(size(abd),m+1,k):end),1,z(sub2ind(size(z),max(k+1,1)):end),1);
end;
if( abs(z(k))>1.0e0 )
s = 1.0e0./abs(z(k));
[n,s,z]=sscal(n,s,z,1);
end;
l = fix(ipvt(k));
t = z(l);
z(l) = z(k);
z(k) = t;
end; kb = fix(n+1);
s = 1.0e0./sasum(n,z,1);
[n,s,z]=sscal(n,s,z,1);
%
ynorm = 1.0e0;
%
%     SOLVE L*V = Y
%
for k = 1 : n;
l = fix(ipvt(k));
t = z(l);
z(l) = z(k);
z(k) = t;
lm = fix(min(ml,n-k));
if( k<n )
[lm,t,abd(sub2ind(size(abd),m+1,k):end),dumvar4,z(sub2ind(size(z),max(k+1,1)):end)]=saxpy(lm,t,abd(sub2ind(size(abd),m+1,k):end),1,z(sub2ind(size(z),max(k+1,1)):end),1);
end;
if( abs(z(k))>1.0e0 )
s = 1.0e0./abs(z(k));
[n,s,z]=sscal(n,s,z,1);
ynorm = s.*ynorm;
end;
end; k = fix(n+1);
s = 1.0e0./sasum(n,z,1);
[n,s,z]=sscal(n,s,z,1);
ynorm = s.*ynorm;
%
%     SOLVE  U*Z = W
%
for kb = 1 : n;
k = fix(n + 1 - kb);
if( abs(z(k))>abs(abd(m,k)) )
s = abs(abd(m,k))./abs(z(k));
[n,s,z]=sscal(n,s,z,1);
ynorm = s.*ynorm;
end;
if( abd(m,k)~=0.0e0 )
z(k) = z(k)./abd(m,k);
end;
if( abd(m,k)==0.0e0 )
z(k) = 1.0e0;
end;
lm = fix(min(k,m) - 1);
la = fix(m - lm);
lz = fix(k - lm);
t = -z(k);
[lm,t,abd(sub2ind(size(abd),la,k):end),dumvar4,z(sub2ind(size(z),max(lz,1)):end)]=saxpy(lm,t,abd(sub2ind(size(abd),la,k):end),1,z(sub2ind(size(z),max(lz,1)):end),1);
end; kb = fix(n+1);
%     MAKE ZNORM = 1.0
s = 1.0e0./sasum(n,z,1);
[n,s,z]=sscal(n,s,z,1);
ynorm = s.*ynorm;
%
if( anorm~=0.0e0 )
rcond = ynorm./anorm;
end;
if( anorm==0.0e0 )
rcond = 0.0e0;
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;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
end
%DECK SGBDI
```

Contact us