| [abd,lda,n,m,rcond,z,info]=spbco(abd,lda,n,m,rcond,z,info); |
function [abd,lda,n,m,rcond,z,info]=spbco(abd,lda,n,m,rcond,z,info);
%***BEGIN PROLOGUE SPBCO
%***PURPOSE Factor a real symmetric positive definite matrix stored in
% band form and estimate the condition number of the matrix.
%***LIBRARY SLATEC (LINPACK)
%***CATEGORY D2B2
%***TYPE SINGLE PRECISION (SPBCO-S, DPBCO-D, CPBCO-C)
%***KEYWORDS BANDED, CONDITION NUMBER, LINEAR ALGEBRA, LINPACK,
% MATRIX FACTORIZATION, POSITIVE DEFINITE
%***AUTHOR Moler, C. B., (U. of New Mexico)
%***DESCRIPTION
%
% SPBCO factors a real symmetric positive definite matrix
% stored in band form and estimates the condition of the matrix.
%
% If RCOND is not needed, SPBFA is slightly faster.
% To solve A*X = B , follow SPBCO by SPBSL.
% To compute INVERSE(A)*C , follow SPBCO by SPBSL.
% To compute DETERMINANT(A) , follow SPBCO by SPBDI.
%
% On Entry
%
% ABD REAL(LDA, N)
% the matrix to be factored. The columns of the upper
% triangle are stored in the columns of ABD and the
% diagonals of the upper triangle are stored in the
% rows of ABD . See the comments below for details.
%
% LDA INTEGER
% the leading dimension of the array ABD .
% LDA must be .GE. M + 1 .
%
% N INTEGER
% the order of the matrix A .
%
% M INTEGER
% the number of diagonals above the main diagonal.
% 0 .LE. M .LT. N .
%
% On Return
%
% ABD an upper triangular matrix R , stored in band
% form, so that A = TRANS(R)*R .
% If INFO .NE. 0 , the factorization is not complete.
%
% 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. If INFO .NE. 0 , RCOND is unchanged.
%
% Z REAL(N)
% a work vector whose contents are usually unimportant.
% If A is singular to working precision, then Z is
% an approximate null vector in the sense that
% NORM(A*Z) = RCOND*NORM(A)*NORM(Z) .
% If INFO .NE. 0 , Z is unchanged.
%
% INFO INTEGER
% = 0 for normal return.
% = K signals an error condition. The leading minor
% of order K is not positive definite.
%
% Band Storage
%
% If A is a symmetric positive definite band matrix,
% the following program segment will set up the input.
%
% M = (band width above diagonal)
% DO 20 J = 1, N
% I1 = MAX(1, J-M)
% DO 10 I = I1, J
% K = I-J+M+1
% ABD(K,J) = A(I,J)
% 10 CONTINUE
% 20 CONTINUE
%
% This uses M + 1 rows of A , except for the M by M
% upper left triangle, which is ignored.
%
% Example: If the original matrix is
%
% 11 12 13 0 0 0
% 12 22 23 24 0 0
% 13 23 33 34 35 0
% 0 24 34 44 45 46
% 0 0 35 45 55 56
% 0 0 0 46 56 66
%
% then N = 6 , M = 2 and ABD should contain
%
% * * 13 24 35 46
% * 12 23 34 45 56
% 11 22 33 44 55 66
%
%***REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
% Stewart, LINPACK Users' Guide, SIAM, 1979.
%***ROUTINES CALLED SASUM, SAXPY, SDOT, SPBFA, 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 SPBCO
persistent anorm ek i j j2 k kb kp1 l la lb lm mu s sm t wk wkm ynorm ;
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(i), i=0; end;
if isempty(j), j=0; end;
if isempty(j2), j2=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(lb), lb=0; end;
if isempty(lm), lm=0; end;
if isempty(mu), mu=0; end;
%
% FIND NORM OF A
%
%***FIRST EXECUTABLE STATEMENT SPBCO
for j = 1 : n;
l = fix(min(j,m+1));
mu = fix(max(m+2-j,1));
[z(j) ,l,abd(sub2ind(size(abd),mu,j):end)]=sasum(l,abd(sub2ind(size(abd),mu,j):end),1);
k = fix(j - l);
if( m>=mu )
for i = mu : m;
k = fix(k + 1);
z(k) = z(k) + abs(abd(i,j));
end; i = fix(m+1);
end;
end; j = fix(n+1);
anorm = 0.0e0;
for j = 1 : n;
anorm = max(anorm,z(j));
end; j = fix(n+1);
%
% FACTOR
%
[abd,lda,n,m,info]=spbfa(abd,lda,n,m,info);
if( info==0 )
%
% RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) .
% ESTIMATE = NORM(Z)/NORM(Y) WHERE A*Z = Y AND A*Y = E .
% THE COMPONENTS OF E ARE CHOSEN TO CAUSE MAXIMUM LOCAL
% GROWTH IN THE ELEMENTS OF W WHERE TRANS(R)*W = E .
% THE VECTORS ARE FREQUENTLY RESCALED TO AVOID OVERFLOW.
%
% SOLVE TRANS(R)*W = E
%
ek = 1.0e0;
for j = 1 : n;
z(j) = 0.0e0;
end; j = fix(n+1);
for k = 1 : n;
if( z(k)~=0.0e0 )
ek = (abs(ek).*sign(-z(k)));
end;
if( abs(ek-z(k))>abd(m+1,k) )
s = abd(m+1,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);
wk = wk./abd(m+1,k);
wkm = wkm./abd(m+1,k);
kp1 = fix(k + 1);
j2 = fix(min(k+m,n));
i = fix(m + 1);
if( kp1<=j2 )
for j = kp1 : j2;
i = fix(i - 1);
sm = sm + abs(z(j)+wkm.*abd(i,j));
z(j) = z(j) + wk.*abd(i,j);
s = s + abs(z(j));
end; j = fix(j2+1);
if( s<sm )
t = wkm - wk;
wk = wkm;
i = fix(m + 1);
for j = kp1 : j2;
i = fix(i - 1);
z(j) = z(j) + t.*abd(i,j);
end; j = fix(j2+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 R*Y = W
%
for kb = 1 : n;
k = fix(n + 1 - kb);
if( abs(z(k))>abd(m+1,k) )
s = abd(m+1,k)./abs(z(k));
[n,s,z]=sscal(n,s,z,1);
end;
z(k) = z(k)./abd(m+1,k);
lm = fix(min(k-1,m));
la = fix(m + 1 - lm);
lb = fix(k - lm);
t = -z(k);
[lm,t,abd(sub2ind(size(abd),la,k):end),dumvar4,z(sub2ind(size(z),max(lb,1)):end)]=saxpy(lm,t,abd(sub2ind(size(abd),la,k):end),1,z(sub2ind(size(z),max(lb,1)):end),1);
end; kb = fix(n+1);
s = 1.0e0./sasum(n,z,1);
[n,s,z]=sscal(n,s,z,1);
%
ynorm = 1.0e0;
%
% SOLVE TRANS(R)*V = Y
%
for k = 1 : n;
lm = fix(min(k-1,m));
la = fix(m + 1 - lm);
lb = fix(k - lm);
z(k) = z(k) - sdot(lm,abd(sub2ind(size(abd),la,k):end),1,z(sub2ind(size(z),max(lb,1)):end),1);
if( abs(z(k))>abd(m+1,k) )
s = abd(m+1,k)./abs(z(k));
[n,s,z]=sscal(n,s,z,1);
ynorm = s.*ynorm;
end;
z(k) = z(k)./abd(m+1,k);
end; k = fix(n+1);
s = 1.0e0./sasum(n,z,1);
[n,s,z]=sscal(n,s,z,1);
ynorm = s.*ynorm;
%
% SOLVE R*Z = W
%
for kb = 1 : n;
k = fix(n + 1 - kb);
if( abs(z(k))>abd(m+1,k) )
s = abd(m+1,k)./abs(z(k));
[n,s,z]=sscal(n,s,z,1);
ynorm = s.*ynorm;
end;
z(k) = z(k)./abd(m+1,k);
lm = fix(min(k-1,m));
la = fix(m + 1 - lm);
lb = fix(k - lm);
t = -z(k);
[lm,t,abd(sub2ind(size(abd),la,k):end),dumvar4,z(sub2ind(size(z),max(lb,1)):end)]=saxpy(lm,t,abd(sub2ind(size(abd),la,k):end),1,z(sub2ind(size(z),max(lb,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;
end;
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 SPBDI
|
|
