| [a,mda,m,n,b,mdb,nb,tau,krank,rnorm,h,g,ip]=dhfti(a,mda,m,n,b,mdb,nb,tau,krank,rnorm,h,g,ip); |
function [a,mda,m,n,b,mdb,nb,tau,krank,rnorm,h,g,ip]=dhfti(a,mda,m,n,b,mdb,nb,tau,krank,rnorm,h,g,ip);
%***BEGIN PROLOGUE DHFTI
%***PURPOSE Solve a least squares problem for banded matrices using
% sequential accumulation of rows of the data matrix.
% Exactly one right-hand side vector is permitted.
%***LIBRARY SLATEC
%***CATEGORY D9
%***TYPE doubleprecision (HFTI-S, DHFTI-D)
%***KEYWORDS CURVE FITTING, LEAST SQUARES
%***AUTHOR Lawson, C. L., (JPL)
% Hanson, R. J., (SNLA)
%***DESCRIPTION
%
% DIMENSION A(MDA,N),(B(MDB,NB) or B(M)),RNORM(NB),H(N),G(N),IP(N)
%
% This subroutine solves a linear least squares problem or a set of
% linear least squares problems having the same matrix but different
% right-side vectors. The problem data consists of an M by N matrix
% A, an M by NB matrix B, and an absolute tolerance parameter TAU
% whose usage is described below. The NB column vectors of B
% represent right-side vectors for NB distinct linear least squares
% problems.
%
% This set of problems can also be written as the matrix least
% squares problem
%
% AX = B,
%
% where X is the N by NB solution matrix.
%
% Note that if B is the M by M identity matrix, then X will be the
% pseudo-inverse of A.
%
% This subroutine first transforms the augmented matrix (A B) to a
% matrix (R C) using premultiplying Householder transformations with
% column interchanges. All subdiagonal elements in the matrix R are
% zero and its diagonal elements satisfy
%
% ABS(R(I,I)).GE.ABS(R(I+1,I+1)),
%
% I = 1,...,L-1, where
%
% L = MIN(M,N).
%
% The subroutine will compute an integer, KRANK, equal to the number
% of diagonal terms of R that exceed TAU in magnitude. Then a
% solution of minimum Euclidean length is computed using the first
% KRANK rows of (R C).
%
% To be specific we suggest that the user consider an easily
% computable matrix norm, such as, the maximum of all column sums of
% magnitudes.
%
% Now if the relative uncertainty of B is EPS, (norm of uncertainty/
% norm of B), it is suggested that TAU be set approximately equal to
% EPS*(norm of A).
%
% The user must dimension all arrays appearing in the call list..
% A(MDA,N),(B(MDB,NB) or B(M)),RNORM(NB),H(N),G(N),IP(N). This
% permits the solution of a range of problems in the same array
% space.
%
% The entire set of parameters for DHFTI are
%
% INPUT.. All TYPE REAL variables are doubleprecision
%
% A(*,*),MDA,M,N The array A(*,*) initially contains the M by N
% matrix A of the least squares problem AX = B.
% The first dimensioning parameter of the array
% A(*,*) is MDA, which must satisfy MDA.GE.M
% Either M.GE.N or M.LT.N is permitted. There
% is no restriction on the rank of A. The
% condition MDA.LT.M is considered an error.
%
% B(*),MDB,NB If NB = 0 the subroutine will perform the
% orthogonal decomposition but will make no
% references to the array B(*). If NB.GT.0
% the array B(*) must initially contain the M by
% NB matrix B of the least squares problem AX =
% B. If NB.GE.2 the array B(*) must be doubly
% subscripted with first dimensioning parameter
% MDB.GE.MAX(M,N). If NB = 1 the array B(*) may
% be either doubly or singly subscripted. In
% the latter case the value of MDB is arbitrary
% but it should be set to some valid integer
% value such as MDB = M.
%
% The condition of NB.GT.1.AND.MDB.LT. MAX(M,N)
% is considered an error.
%
% TAU Absolute tolerance parameter provided by user
% for pseudorank determination.
%
% H(*),G(*),IP(*) Arrays of working space used by DHFTI.
%
% OUTPUT.. All TYPE REAL variables are doubleprecision
%
% A(*,*) The contents of the array A(*,*) will be
% modified by the subroutine. These contents
% are not generally required by the user.
%
% B(*) On return the array B(*) will contain the N by
% NB solution matrix X.
%
% KRANK Set by the subroutine to indicate the
% pseudorank of A.
%
% RNORM(*) On return, RNORM(J) will contain the Euclidean
% norm of the residual vector for the problem
% defined by the J-th column vector of the array
% B(*,*) for J = 1,...,NB.
%
% H(*),G(*) On return these arrays respectively contain
% elements of the pre- and post-multiplying
% Householder transformations used to compute
% the minimum Euclidean length solution.
%
% IP(*) Array in which the subroutine records indices
% describing the permutation of column vectors.
% The contents of arrays H(*),G(*) and IP(*)
% are not generally required by the user.
%
%***REFERENCES C. L. Lawson and R. J. Hanson, Solving Least Squares
% Problems, Prentice-Hall, Inc., 1974, Chapter 14.
%***ROUTINES CALLED D1MACH, DH12, XERMSG
%***REVISION HISTORY (YYMMDD)
% 790101 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (WRB)
% 891006 Cosmetic changes to prologue. (WRB)
% 891006 REVISION DATE from Version 3.2
% 891214 Prologue converted to Version 4.0 format. (BAB)
% 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
% 901005 Replace usage of DDIFF with usage of D1MACH. (RWC)
% 920501 Reformatted the REFERENCES section. (WRB)
%***end PROLOGUE DHFTI
persistent dzero factor firstCall gt hmax i ii iopt ip1 j jb jj k kp1 l ldiag lmax nerr releps sm sm1 szero tmp ; if isempty(firstCall),firstCall=1;end;
if isempty(i), i=0; end;
if isempty(ii), ii=0; end;
if isempty(iopt), iopt=0; end;
ip_shape=size(ip);ip=reshape(ip,1,[]);
if isempty(ip1), ip1=0; end;
if isempty(j), j=0; end;
if isempty(jb), jb=0; end;
if isempty(jj), jj=0; end;
if isempty(k), k=0; end;
if isempty(kp1), kp1=0; end;
if isempty(l), l=0; end;
if isempty(ldiag), ldiag=0; end;
if isempty(lmax), lmax=0; end;
if isempty(nerr), nerr=0; end;
if isempty(gt), gt=0; end;
if isempty(dzero), dzero=0; end;
if isempty(factor), factor=0; end;
if isempty(hmax), hmax=0; end;
if isempty(releps), releps=0; end;
if isempty(sm), sm=0; end;
if isempty(sm1), sm1=0; end;
if isempty(szero), szero=0; end;
if isempty(tmp), tmp=0; end;
a_shape=size(a);a=reshape([a(:).',zeros(1,ceil(numel(a)./prod([mda])).*prod([mda])-numel(a))],mda,[]);
b_shape=size(b);b=reshape([b(:).',zeros(1,ceil(numel(b)./prod([mdb])).*prod([mdb])-numel(b))],mdb,[]);
h_shape=size(h);h=reshape(h,1,[]);
g_shape=size(g);g=reshape(g,1,[]);
rnorm_shape=size(rnorm);rnorm=reshape(rnorm,1,[]);
if firstCall, releps=[0.0d0]; end;
firstCall=0;
% BEGIN BLOCK PERMITTING ...EXITS TO 360
%***FIRST EXECUTABLE STATEMENT DHFTI
gt=0;
if( releps==0.0d0 )
[ releps ]=d1mach(4);
end;
szero = 0.0d0;
dzero = 0.0d0;
factor = 0.001d0;
%
k = 0;
ldiag = fix(min(m,n));
if( ldiag>0 )
% BEGIN BLOCK PERMITTING ...EXITS TO 130
% BEGIN BLOCK PERMITTING ...EXITS TO 120
if( mda<m )
nerr = 1;
iopt = 2;
[dumvar1,dumvar2,dumvar3,nerr,iopt]=xermsg('SLATEC','DHFTI','MDA.LT.M, PROBABLE ERROR.',nerr,iopt);
% ...............EXIT
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
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;
h_shape=zeros(h_shape);h_shape(:)=h(1:numel(h_shape));h=h_shape;
g_shape=zeros(g_shape);g_shape(:)=g(1:numel(g_shape));g=g_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
return;
%
elseif( nb<=1 || max(m,n)<=mdb ) ;
%
for j = 1 : ldiag;
gt=0;
% BEGIN BLOCK PERMITTING ...EXITS TO 70
if( j~=1 )
%
% UPDATE SQUARED COLUMN LENGTHS AND FIND LMAX
% ..
lmax = fix(j);
for l = j : n;
h(l) = h(l) - a(j-1,l).^2;
if( h(l)>h(lmax) )
lmax = fix(l);
end;
end; l = fix(n+1);
% ......EXIT
if( factor.*h(lmax)>hmax.*releps )
gt=1;
end;
end;
%
% COMPUTE SQUARED COLUMN LENGTHS AND FIND LMAX
% ..
if(gt==0)
lmax = fix(j);
for l = j : n;
h(l) = 0.0d0;
for i = j : m;
h(l) = h(l) + a(i,l).^2;
end; i = fix(m+1);
if( h(l)>h(lmax) )
lmax = fix(l);
end;
end; l = fix(n+1);
hmax = h(lmax);
end;
gt=0;
% ..
% LMAX HAS BEEN DETERMINED
%
% DO COLUMN INTERCHANGES IF NEEDED.
% ..
ip(j) = fix(lmax);
if( ip(j)~=j )
for i = 1 : m;
tmp = a(i,j);
a(i,j) = a(i,lmax);
a(i,lmax) = tmp;
end; i = fix(m+1);
h(lmax) = h(j);
end;
%
% COMPUTE THE J-TH TRANSFORMATION AND APPLY IT TO A
% AND B.
% ..
[dumvar1,j,dumvar3,m,dumvar5,dumvar6,h(j),dumvar8,dumvar9,mda]=dh12(1,j,j+1,m,a(sub2ind(size(a),1,j):end),1,h(j),a(sub2ind(size(a),1,j+1):end),1,mda,n-j); a(sub2ind(size(a),1,j):end)=dumvar5; a(sub2ind(size(a),1,j+1):end)=dumvar8;
[dumvar1,j,dumvar3,m,a(sub2ind(size(a),1,j):end),dumvar6,h(j),b,dumvar9,mdb,nb]=dh12(2,j,j+1,m,a(sub2ind(size(a),1,j):end),1,h(j),b,1,mdb,nb);
end; j = fix(ldiag+1);
%
% DETERMINE THE PSEUDORANK, K, USING THE TOLERANCE,
% TAU.
% ..
gt=0;
for j = 1 : ldiag;
% ......EXIT
if( abs(a(j,j))<=tau )
k = fix(j - 1);
gt=1;
break;
end;
end;
% ......EXIT
if(gt==0)
k = fix(ldiag);
end;
kp1 = fix(k + 1);
%
% COMPUTE THE NORMS OF THE RESIDUAL VECTORS.
%
if( nb>=1 )
for jb = 1 : nb;
tmp = szero;
if( m>=kp1 )
for i = kp1 : m;
tmp = tmp + b(i,jb).^2;
end; i = fix(m+1);
end;
rnorm(jb) = sqrt(tmp);
end; jb = fix(nb+1);
end;
% SPECIAL FOR PSEUDORANK = 0
if( k>0 )
%
% IF THE PSEUDORANK IS LESS THAN N COMPUTE HOUSEHOLDER
% DECOMPOSITION OF FIRST K ROWS.
% ..
if( k~=n )
for ii = 1 : k;
i = fix(kp1 - ii);
mda_orig=mda; [dumvar1,i,kp1,n,dumvar5,mda,g(i),a,dumvar9]=dh12(1,i,kp1,n,a(sub2ind(size(a),i,1):end),mda,g(i),a,mda,1,i-1); mda(dumvar9~=mda_orig)=dumvar9(dumvar9~=mda_orig); a(sub2ind(size(a),i,1):end)=dumvar5;
end; ii = fix(k+1);
end;
%
%
if( nb>=1 )
for jb = 1 : nb;
%
% SOLVE THE K BY K TRIANGULAR SYSTEM.
% ..
for l = 1 : k;
sm = dzero;
i = fix(kp1 - l);
ip1 = fix(i + 1);
if( k>=ip1 )
for j = ip1 : k;
sm = sm + a(i,j).*b(j,jb);
end; j = fix(k+1);
end;
sm1 = sm;
b(i,jb) =(b(i,jb)-sm1)./a(i,i);
end; l = fix(k+1);
%
% COMPLETE COMPUTATION OF SOLUTION VECTOR.
% ..
if( k~=n )
for j = kp1 : n;
b(j,jb) = szero;
end; j = fix(n+1);
for i = 1 : k;
[dumvar1,i,kp1,n,a(sub2ind(size(a),i,1):end),mda,g(i),b(sub2ind(size(b),1,jb):end),dumvar9,mdb]=dh12(2,i,kp1,n,a(sub2ind(size(a),i,1):end),mda,g(i),b(sub2ind(size(b),1,jb):end),1,mdb,1);
end; i = fix(k+1);
end;
%
% RE-ORDER THE SOLUTION VECTOR TO COMPENSATE FOR THE
% COLUMN INTERCHANGES.
% ..
for jj = 1 : ldiag;
j = fix(ldiag + 1 - jj);
if( ip(j)~=j )
l = fix(ip(j));
tmp = b(l,jb);
b(l,jb) = b(j,jb);
b(j,jb) = tmp;
end;
end; jj = fix(ldiag+1);
end; jb = fix(nb+1);
end;
elseif( nb>=1 ) ;
for jb = 1 : nb;
for i = 1 : n;
b(i,jb) = szero;
end; i = fix(n+1);
end; jb = fix(nb+1);
end;
else;
nerr = 2;
iopt = 2;
[dumvar1,dumvar2,dumvar3,nerr,iopt]=xermsg('SLATEC','DHFTI','MDB.LT.MAX(M,N).AND.NB.GT.1. PROBABLE ERROR.',nerr,iopt);
% ...............EXIT
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
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;
h_shape=zeros(h_shape);h_shape(:)=h(1:numel(h_shape));h=h_shape;
g_shape=zeros(g_shape);g_shape(:)=g(1:numel(g_shape));g=g_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
return;
end;
end;
% ..
% THE SOLUTION VECTORS, X, ARE NOW
% IN THE FIRST N ROWS OF THE ARRAY B().
%
krank = fix(k);
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
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;
h_shape=zeros(h_shape);h_shape(:)=h(1:numel(h_shape));h=h_shape;
g_shape=zeros(g_shape);g_shape(:)=g(1:numel(g_shape));g=g_shape;
rnorm_shape=zeros(rnorm_shape);rnorm_shape(:)=rnorm(1:numel(rnorm_shape));rnorm=rnorm_shape;
end %subroutine dhfti
%DECK DHKSEQ
|
|
