| [w,mdw,ma,mg,n,prgopt,x,rnorm,mode,ws,ip]=lsi(w,mdw,ma,mg,n,prgopt,x,rnorm,mode,ws,ip); |
function [w,mdw,ma,mg,n,prgopt,x,rnorm,mode,ws,ip]=lsi(w,mdw,ma,mg,n,prgopt,x,rnorm,mode,ws,ip);
%***BEGIN PROLOGUE LSI
%***SUBSIDIARY
%***PURPOSE Subsidiary to LSEI
%***LIBRARY SLATEC
%***TYPE SINGLE PRECISION (LSI-S, DLSI-D)
%***AUTHOR Hanson, R. J., (SNLA)
%***DESCRIPTION
%
% This is a companion subprogram to LSEI. The documentation for
% LSEI has complete usage instructions.
%
% Solve..
% AX = B, A MA by N (least squares equations)
% subject to..
%
% GX.GE.H, G MG by N (inequality constraints)
%
% Input..
%
% W(*,*) contains (A B) in rows 1,...,MA+MG, cols 1,...,N+1.
% (G H)
%
% MDW,MA,MG,N
% contain (resp) var. dimension of W(*,*),
% and matrix dimensions.
%
% PRGOPT(*),
% program option vector.
%
% OUTPUT..
%
% X(*),RNORM
%
% Solution vector(unless MODE=2), length of AX-B.
%
% MODE
% =0 Inequality constraints are compatible.
% =2 Inequality constraints contradictory.
%
% WS(*),
% Working storage of dimension K+N+(MG+2)*(N+7),
% where K=MAX(MA+MG,N).
% IP(MG+2*N+1)
% Integer working storage
%
%***ROUTINES CALLED H12, HFTI, LPDP, R1MACH, SASUM, SAXPY, SCOPY, SDOT,
% SSCAL, SSWAP
%***REVISION HISTORY (YYMMDD)
% 790701 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (WRB)
% 890618 Completely restructured and extensively revised (WRBRWC)
% 891214 Prologue converted to Version 4.0 format. (BAB)
% 900328 Added TYPE section. (WRB)
% 920422 Changed CALL to HFTI to include variable MA. (WRB)
%***end PROLOGUE LSI
persistent anorm cov fac first firstCall gam i j k key krank krm1 krp1 l last link m map1 mdlpdp minman n1 n2 n3 next np1 rb sclcov srelpr tau tol xnorm ; if isempty(firstCall),firstCall=1;end;
ip_shape=size(ip);ip=reshape(ip,1,[]);
prgopt_shape=size(prgopt);prgopt=reshape(prgopt,1,[]);
w_shape=size(w);w=reshape([w(:).',zeros(1,ceil(numel(w)./prod([mdw])).*prod([mdw])-numel(w))],mdw,[]);
ws_shape=size(ws);ws=reshape(ws,1,[]);
x_shape=size(x);x=reshape(x,1,[]);
%
%
if isempty(anorm), anorm=0; end;
if isempty(fac), fac=0; end;
if isempty(gam), gam=0; end;
if isempty(rb), rb=0; end;
if isempty(srelpr), srelpr=0; end;
if isempty(tau), tau=0; end;
if isempty(tol), tol=0; end;
if isempty(xnorm), xnorm=0; end;
if isempty(i), i=0; end;
if isempty(j), j=0; end;
if isempty(k), k=0; end;
if isempty(key), key=0; end;
if isempty(krank), krank=0; end;
if isempty(krm1), krm1=0; end;
if isempty(krp1), krp1=0; end;
if isempty(l), l=0; end;
if isempty(last), last=0; end;
if isempty(link), link=0; end;
if isempty(m), m=0; end;
if isempty(map1), map1=0; end;
if isempty(mdlpdp), mdlpdp=0; end;
if isempty(minman), minman=0; end;
if isempty(n1), n1=0; end;
if isempty(n2), n2=0; end;
if isempty(n3), n3=0; end;
if isempty(next), next=0; end;
if isempty(np1), np1=0; end;
if isempty(cov), cov=false; end;
if isempty(first), first=false; end;
if isempty(sclcov), sclcov=false; end;
%
if firstCall, first=[true]; end;
firstCall=0;
%
%***FIRST EXECUTABLE STATEMENT LSI
%
% Set the nominal tolerance used in the code.
%
if( first )
[ srelpr ]=r1mach(4);
end;
first = false;
tol = sqrt(srelpr);
%
mode = 0;
rnorm = 0.0e0;
m = fix(ma + mg);
np1 = fix(n + 1);
krank = 0;
if( n>0 && m>0 )
%
% To process option vector.
%
cov = false;
sclcov = true;
last = 1;
link = fix(prgopt(1));
%
while( link>1 );
key = fix(prgopt(last+1));
if( key==1 )
cov = prgopt(last+2)~=0.0e0;
end;
if( key==10 )
sclcov = prgopt(last+2)==0.0e0;
end;
if( key==5 )
tol = max(srelpr,prgopt(last+2));
end;
next = fix(prgopt(link));
last = fix(link);
link = fix(next);
end;
%
% Compute matrix norm of least squares equations.
%
anorm = 0.0e0;
for j = 1 : n;
anorm = max(anorm,sasum(ma,w(sub2ind(size(w),1,j):end),1));
end; j = fix(n+1);
%
% Set tolerance for HFTI( ) rank test.
%
tau = tol.*anorm;
%
% Compute Householder orthogonal decomposition of matrix.
%
[n,dumvar2,dumvar3,ws]=scopy(n,0.0e0,0,ws,1);
[ma,w(sub2ind(size(w),1,np1):end),dumvar3,ws]=scopy(ma,w(sub2ind(size(w),1,np1):end),1,ws,1);
k = fix(max(m,n));
minman = fix(min(ma,n));
n1 = fix(k + 1);
n2 = fix(n1 + n);
ma_orig=ma; [w,mdw,ma,n,ws,dumvar6,dumvar7,tau,krank,rnorm,dumvar11,dumvar12,ip]=hfti(w,mdw,ma,n,ws,ma,1,tau,krank,rnorm,ws(sub2ind(size(ws),max(n2,1)):end),ws(sub2ind(size(ws),max(n1,1)):end),ip); ma(dumvar6~=ma_orig)=dumvar6(dumvar6~=ma_orig); dumvar11i=find((ws(sub2ind(size(ws),max(n2,1)):end))~=(dumvar11));dumvar12i=find((ws(sub2ind(size(ws),max(n1,1)):end))~=(dumvar12)); ws(n2-1+dumvar11i)=dumvar11(dumvar11i); ws(n1-1+dumvar12i)=dumvar12(dumvar12i);
fac = 1.0e0;
gam = ma - krank;
if( krank<ma && sclcov )
fac = rnorm.^2./gam;
end;
%
% Reduce to LPDP and solve.
%
map1 = fix(ma + 1);
%
% Compute inequality rt-hand side for LPDP.
%
if( ma<m )
if( minman>0 )
for i = map1 : m;
w(i,np1) = w(i,np1) - sdot(n,w(sub2ind(size(w),i,1):end),mdw,ws,1);
end; i = fix(m+1);
%
% Apply permutations to col. of inequality constraint matrix.
%
for i = 1 : minman;
[mg,dumvar2,dumvar3,dumvar4]=sswap(mg,w(sub2ind(size(w),map1,i):end),1,w(sub2ind(size(w),map1,ip(i)):end),1); w(sub2ind(size(w),map1,i):end)=dumvar2; w(sub2ind(size(w),map1,ip(i)):end)=dumvar4;
end; i = fix(minman+1);
%
% Apply Householder transformations to constraint matrix.
%
if( krank>0 && krank<n )
for i = krank : -1: 1 ;
mdw_orig=mdw; [dumvar1,i,dumvar3,n,dumvar5,mdw,ws(n1+i-1),dumvar8,dumvar9,dumvar10,mg]=h12(2,i,krank+1,n,w(sub2ind(size(w),i,1):end),mdw,ws(n1+i-1),w(sub2ind(size(w),map1,1):end),mdw,1,mg); mdw(dumvar9~=mdw_orig)=dumvar9(dumvar9~=mdw_orig); w(sub2ind(size(w),i,1):end)=dumvar5; w(sub2ind(size(w),map1,1):end)=dumvar8;
end; i = fix(1 -1);
end;
%
% Compute permuted inequality constraint matrix times r-inv.
%
for i = map1 : m;
for j = 1 : krank;
w(i,j) =(w(i,j)-sdot(j-1,w(sub2ind(size(w),1,j):end),1,w(sub2ind(size(w),i,1):end),mdw))./w(j,j);
end; j = fix(krank+1);
end; i = fix(m+1);
end;
%
% Solve the reduced problem with LPDP algorithm,
% the least projected distance problem.
%
[w(sub2ind(size(w),map1,1):end),mdw,mg,krank,dumvar5,prgopt,x,xnorm,mdlpdp,ws(sub2ind(size(ws),max(n2,1)):end),ip(sub2ind(size(ip),max(n+1,1)):end)]=lpdp(w(sub2ind(size(w),map1,1):end),mdw,mg,krank,n-krank,prgopt,x,xnorm,mdlpdp,ws(sub2ind(size(ws),max(n2,1)):end),ip(sub2ind(size(ip),max(n+1,1)):end));
%
% Compute solution in original coordinates.
%
if( mdlpdp==1 )
for i = krank : -1: 1 ;
x(i) =(x(i)-sdot(krank-i,w(sub2ind(size(w),i,i+1):end),mdw,x(sub2ind(size(x),max(i+1,1)):end),1))./w(i,i);
end; i = fix(1 -1);
%
% Apply Householder transformation to solution vector.
%
if( krank<n )
for i = 1 : krank;
[dumvar1,i,dumvar3,n,w(sub2ind(size(w),i,1):end),mdw,ws(n1+i-1),x]=h12(2,i,krank+1,n,w(sub2ind(size(w),i,1):end),mdw,ws(n1+i-1),x,1,1,1);
end; i = fix(krank+1);
end;
%
% Repermute variables to their input order.
%
if( minman>0 )
for i = minman : -1: 1 ;
[dumvar1,dumvar2,dumvar3,dumvar4]=sswap(1,x(sub2ind(size(x),max(i,1)):end),1,x(sub2ind(size(x),max(ip(i),1)):end),1); dumvar2i=find((x(sub2ind(size(x),max(i,1)):end))~=(dumvar2));dumvar4i=find((x(sub2ind(size(x),max(ip(i),1)):end))~=(dumvar4)); x(i-1+dumvar2i)=dumvar2(dumvar2i); x(ip(i)-1+dumvar4i)=dumvar4(dumvar4i);
end; i = fix(1 -1);
%
% Variables are now in original coordinates.
% Add solution of unconstrained problem.
%
for i = 1 : n;
x(i) = x(i) + ws(i);
end; i = fix(n+1);
%
% Compute the residual vector norm.
%
rnorm = sqrt(rnorm.^2+xnorm.^2);
end;
else;
mode = 2;
end;
else;
[n,ws,dumvar3,x]=scopy(n,ws,1,x,1);
end;
%
% Compute covariance matrix based on the orthogonal decomposition
% from HFTI( ).
%
if( ~(~cov || krank<=0) )
krm1 = fix(krank - 1);
krp1 = fix(krank + 1);
%
% Copy diagonal terms to working array.
%
[krank,w,dumvar3,ws(sub2ind(size(ws),max(n2,1)):end)]=scopy(krank,w,mdw+1,ws(sub2ind(size(ws),max(n2,1)):end),1);
%
% Reciprocate diagonal terms.
%
for j = 1 : krank;
w(j,j) = 1.0e0./w(j,j);
end; j = fix(krank+1);
%
% Invert the upper triangular QR factor on itself.
%
if( krank>1 )
for i = 1 : krm1;
for j = i + 1 : krank;
w(i,j) = -sdot(j-i,w(sub2ind(size(w),i,i):end),mdw,w(sub2ind(size(w),i,j):end),1).*w(j,j);
end; j = fix(krank+1);
end; i = fix(krm1+1);
end;
%
% Compute the inverted factor times its transpose.
%
for i = 1 : krank;
for j = i : krank;
w(i,j) = sdot(krank+1-j,w(sub2ind(size(w),i,j):end),mdw,w(sub2ind(size(w),j,j):end),mdw);
end; j = fix(krank+1);
end; i = fix(krank+1);
%
% Zero out lower trapezoidal part.
% Copy upper triangular to lower triangular part.
%
if( krank<n )
for j = 1 : krank;
[j,dumvar2,dumvar3,dumvar4,mdw]=scopy(j,w(sub2ind(size(w),1,j):end),1,w(sub2ind(size(w),j,1):end),mdw); w(sub2ind(size(w),1,j):end)=dumvar2; w(sub2ind(size(w),j,1):end)=dumvar4;
end; j = fix(krank+1);
%
for i = krp1 : n;
[i,dumvar2,dumvar3,w(sub2ind(size(w),i,1):end),mdw]=scopy(i,0.0e0,0,w(sub2ind(size(w),i,1):end),mdw);
end; i = fix(n+1);
%
% Apply right side transformations to lower triangle.
%
n3 = fix(n2 + krp1);
for i = 1 : krank;
l = fix(n1 + i);
k = fix(n2 + i);
rb = ws(l-1).*ws(k-1);
%
% If RB.GE.0.0E0, transformation can be regarded as zero.
%
if( rb<0.0e0 )
rb = 1.0e0./rb;
%
% Store unscaled rank one Householder update in work array.
%
[n,dumvar2,dumvar3,ws(sub2ind(size(ws),max(n3,1)):end)]=scopy(n,0.0e0,0,ws(sub2ind(size(ws),max(n3,1)):end),1);
l = fix(n1 + i);
k = fix(n3 + i);
ws(k-1) = ws(l-1);
%
for j = krp1 : n;
ws(n3+j-1) = w(i,j);
end; j = fix(n+1);
%
for j = 1 : n;
ws(j) = rb.*(sdot(j-i,w(sub2ind(size(w),j,i):end),mdw,ws(sub2ind(size(ws),max(n3+i-1,1)):end),1)+sdot(n-j+1,w(sub2ind(size(w),j,j):end),1,ws(sub2ind(size(ws),max(n3+j-1,1)):end),1));
end; j = fix(n+1);
%
l = fix(n3 + i);
gam = 0.5e0.*rb.*sdot(n-i+1,ws(sub2ind(size(ws),max(l-1,1)):end),1,ws(sub2ind(size(ws),max(i,1)):end),1);
[dumvar1,gam,dumvar3,dumvar4,dumvar5]=saxpy(n-i+1,gam,ws(sub2ind(size(ws),max(l-1,1)):end),1,ws(sub2ind(size(ws),max(i,1)):end),1); dumvar3i=find((ws(sub2ind(size(ws),max(l-1,1)):end))~=(dumvar3));dumvar5i=find((ws(sub2ind(size(ws),max(i,1)):end))~=(dumvar5)); ws(l-1-1+dumvar3i)=dumvar3(dumvar3i); ws(i-1+dumvar5i)=dumvar5(dumvar5i);
for j = i : n;
for l = 1 : i - 1;
w(j,l) = w(j,l) + ws(n3+j-1).*ws(l);
end; l = fix(i - 1+1);
%
for l = i : j;
w(j,l) = w(j,l) + ws(j).*ws(n3+l-1) + ws(l).*ws(n3+j-1);
end; l = fix(j+1);
end; j = fix(n+1);
end;
end; i = fix(krank+1);
%
% Copy lower triangle to upper triangle to symmetrize the
% covariance matrix.
%
for i = 1 : n;
[i,dumvar2,mdw,dumvar4]=scopy(i,w(sub2ind(size(w),i,1):end),mdw,w(sub2ind(size(w),1,i):end),1); w(sub2ind(size(w),i,1):end)=dumvar2; w(sub2ind(size(w),1,i):end)=dumvar4;
end; i = fix(n+1);
end;
%
% Repermute rows and columns.
%
for i = minman : -1: 1 ;
k = fix(ip(i));
if( i~=k )
[dumvar1,dumvar2,dumvar3,dumvar4]=sswap(1,w(sub2ind(size(w),i,i):end),1,w(sub2ind(size(w),k,k):end),1); w(sub2ind(size(w),i,i):end)=dumvar2; w(sub2ind(size(w),k,k):end)=dumvar4;
[dumvar1,dumvar2,dumvar3,dumvar4]=sswap(i-1,w(sub2ind(size(w),1,i):end),1,w(sub2ind(size(w),1,k):end),1); w(sub2ind(size(w),1,i):end)=dumvar2; w(sub2ind(size(w),1,k):end)=dumvar4;
[dumvar1,dumvar2,mdw,dumvar4]=sswap(k-i-1,w(sub2ind(size(w),i,i+1):end),mdw,w(sub2ind(size(w),i+1,k):end),1); w(sub2ind(size(w),i,i+1):end)=dumvar2; w(sub2ind(size(w),i+1,k):end)=dumvar4;
mdw_orig=mdw; [dumvar1,dumvar2,mdw,dumvar4,dumvar5]=sswap(n-k,w(sub2ind(size(w),i,k+1):end),mdw,w(sub2ind(size(w),k,k+1):end),mdw); mdw(dumvar5~=mdw_orig)=dumvar5(dumvar5~=mdw_orig); w(sub2ind(size(w),i,k+1):end)=dumvar2; w(sub2ind(size(w),k,k+1):end)=dumvar4;
end;
end; i = fix(1 -1);
%
% Put in normalized residual sum of squares scale factor
% and symmetrize the resulting covariance matrix.
%
for j = 1 : n;
[j,fac,w(sub2ind(size(w),1,j):end)]=sscal(j,fac,w(sub2ind(size(w),1,j):end),1);
[j,dumvar2,dumvar3,dumvar4,mdw]=scopy(j,w(sub2ind(size(w),1,j):end),1,w(sub2ind(size(w),j,1):end),mdw); w(sub2ind(size(w),1,j):end)=dumvar2; w(sub2ind(size(w),j,1):end)=dumvar4;
end; j = fix(n+1);
end;
end;
%
ip(1) = fix(krank);
ip(2) = fix(n + max(m,n) +(mg+2).*(n+7));
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
prgopt_shape=zeros(prgopt_shape);prgopt_shape(:)=prgopt(1:numel(prgopt_shape));prgopt=prgopt_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
ws_shape=zeros(ws_shape);ws_shape(:)=ws(1:numel(ws_shape));ws=ws_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
end %subroutine lsi
%DECK LSOD
|
|
