[a,mda,m,n1,n2,prgopt,x,wnorm,mode,ws,is]=lpdp(a,mda,m,n1,n2,prgopt,x,wnorm,mode,ws,is);

function [a,mda,m,n1,n2,prgopt,x,wnorm,mode,ws,is]=lpdp(a,mda,m,n1,n2,prgopt,x,wnorm,mode,ws,is);
persistent fac firstCall i iw ix j l modew n np1 one rnorm sc ynorm zero ; if isempty(firstCall),firstCall=1;end; 

if isempty(i), i=0; end;
if isempty(iw), iw=0; end;
if isempty(ix), ix=0; end;
if isempty(j), j=0; end;
if isempty(l), l=0; end;
if isempty(modew), modew=0; end;
if isempty(n), n=0; end;
if isempty(np1), np1=0; end;
%***BEGIN PROLOGUE  LPDP
%***SUBSIDIARY
%***PURPOSE  Subsidiary to LSEI
%***LIBRARY   SLATEC
%***TYPE      SINGLE PRECISION (LPDP-S, DLPDP-D)
%***AUTHOR  Hanson, R. J., (SNLA)
%           Haskell, K. H., (SNLA)
%***DESCRIPTION
%
%     DIMENSION A(MDA,N+1),PRGOPT(*),X(N),WS((M+2)*(N+7)),IS(M+N+1),
%     where N=N1+N2.  This is a slight overestimate for WS(*).
%
%     Determine an N1-vector W, and
%               an N2-vector Z
%     which minimizes the Euclidean length of W
%     subject to G*W+H*Z .GE. Y.
%     This is the least projected distance problem, LPDP.
%     The matrices G and H are of respective
%     dimensions M by N1 and M by N2.
%
%     Called by subprogram LSI( ).
%
%     The matrix
%                (G H Y)
%
%     occupies rows 1,...,M and cols 1,...,N1+N2+1 of A(*,*).
%
%     The solution (W) is returned in X(*).
%                  (Z)
%
%     The value of MODE indicates the status of
%     the computation after returning to the user.
%
%          MODE=1  The solution was successfully obtained.
%
%          MODE=2  The inequalities are inconsistent.
%
%***SEE ALSO  LSEI
%***ROUTINES CALLED  SCOPY, SDOT, SNRM2, SSCAL, WNNLS
%***REVISION HISTORY  (YYMMDD)
%   790701  DATE WRITTEN
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900328  Added TYPE section.  (WRB)
%   910408  Updated the AUTHOR section.  (WRB)
%***end PROLOGUE  LPDP
%
%     SUBROUTINES CALLED
%
%     WNNLS         SOLVES A NONNEGATIVELY CONSTRAINED LINEAR LEAST
%                   SQUARES PROBLEM WITH LINEAR EQUALITY CONSTRAINTS.
%                   PART OF THIS PACKAGE.
%
%++
%     SDOT,         SUBROUTINES FROM THE BLAS PACKAGE.
%     SSCAL,SNRM2,  SEE TRANS. MATH. SOFT., VOL. 5, NO. 3, P. 308.
%     SCOPY
%
a_shape=size(a);a=reshape([a(:).',zeros(1,ceil(numel(a)./prod([mda])).*prod([mda])-numel(a))],mda,[]);
prgopt_shape=size(prgopt);prgopt=reshape(prgopt,1,[]);
ws_shape=size(ws);ws=reshape(ws,1,[]);
x_shape=size(x);x=reshape(x,1,[]);
is_shape=size(is);is=reshape(is,1,[]);
if isempty(fac), fac=0; end;
if isempty(one), one=0; end;
if isempty(rnorm), rnorm=0; end;
if isempty(sc), sc=0; end;
if isempty(ynorm), ynorm=0; end;
if isempty(zero), zero=0; end;
if firstCall,   zero =[0.0e0];  end;
if firstCall,  one=[1.0e0];  end;
if firstCall,  fac=[0.1e0];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  LPDP
n = fix(n1 + n2);
mode = 1;
if( m>0 )
np1 = fix(n + 1);
%
%     SCALE NONZERO ROWS OF INEQUALITY MATRIX TO HAVE LENGTH ONE.
for i = 1 : m;
[sc ,n,a(sub2ind(size(a),i,1):end),mda]=snrm2(n,a(sub2ind(size(a),i,1):end),mda);
if( sc~=zero )
sc = one./sc;
[np1,sc,a(sub2ind(size(a),i,1):end),mda]=sscal(np1,sc,a(sub2ind(size(a),i,1):end),mda);
end;
end; i = fix(m+1);
%
%     SCALE RT.-SIDE VECTOR TO HAVE LENGTH ONE (OR ZERO).
[ynorm ,m,a(sub2ind(size(a),1,np1):end)]=snrm2(m,a(sub2ind(size(a),1,np1):end),1);
if( ynorm~=zero )
sc = one./ynorm;
[m,sc,a(sub2ind(size(a),1,np1):end)]=sscal(m,sc,a(sub2ind(size(a),1,np1):end),1);
end;
%
%     SCALE COLS OF MATRIX H.
j = fix(n1 + 1);
while( j<=n );
[sc ,m,a(sub2ind(size(a),1,j):end)]=snrm2(m,a(sub2ind(size(a),1,j):end),1);
if( sc~=zero )
sc = one./sc;
end;
[m,sc,a(sub2ind(size(a),1,j):end)]=sscal(m,sc,a(sub2ind(size(a),1,j):end),1);
x(j) = sc;
j = fix(j + 1);
end;
if( n1>0 )
%
%     COPY TRANSPOSE OF (H G Y) TO WORK ARRAY WS(*).
iw = 0;
for i = 1 : m;
%
%     MOVE COL OF TRANSPOSE OF H INTO WORK ARRAY.
[n2,a(sub2ind(size(a),i,n1+1):end),mda,ws(sub2ind(size(ws),max(iw+1,1)):end)]=scopy(n2,a(sub2ind(size(a),i,n1+1):end),mda,ws(sub2ind(size(ws),max(iw+1,1)):end),1);
iw = fix(iw + n2);
%
%     MOVE COL OF TRANSPOSE OF G INTO WORK ARRAY.
[n1,a(sub2ind(size(a),i,1):end),mda,ws(sub2ind(size(ws),max(iw+1,1)):end)]=scopy(n1,a(sub2ind(size(a),i,1):end),mda,ws(sub2ind(size(ws),max(iw+1,1)):end),1);
iw = fix(iw + n1);
%
%     MOVE COMPONENT OF VECTOR Y INTO WORK ARRAY.
ws(iw+1) = a(i,np1);
iw = fix(iw + 1);
end; i = fix(m+1);
ws(iw+1) = zero;
[n,dumvar2,dumvar3,dumvar2]=scopy(n,ws(sub2ind(size(ws),max(iw+1,1)):end),0,ws(sub2ind(size(ws),max(iw+1,1)):end),1);   dumvar2i=find((ws(sub2ind(size(ws),max(iw+1,1)):end))~=(dumvar2));   ws(iw+1-1+dumvar2i)=dumvar2(dumvar2i); 
iw = fix(iw + n);
ws(iw+1) = one;
iw = fix(iw + 1);
%
%     SOLVE EU=F SUBJECT TO (TRANSPOSE OF H)U=0, U.GE.0.  THE
%     MATRIX E = TRANSPOSE OF (G Y), AND THE (N+1)-VECTOR
%     F = TRANSPOSE OF (0,...,0,1).
ix = fix(iw + 1);
iw = fix(iw + m);
%
%     DO NOT CHECK LENGTHS OF WORK ARRAYS IN THIS USAGE OF WNNLS( ).
is(1) = 0;
is(2) = 0;
[ws,np1,n2,dumvar4,m,dumvar6,prgopt,dumvar8,rnorm,modew,is,dumvar12]=wnnls(ws,np1,n2,np1-n2,m,0,prgopt,ws(sub2ind(size(ws),max(ix,1)):end),rnorm,modew,is,ws(sub2ind(size(ws),max(iw+1,1)):end));   dumvar8i=find((ws(sub2ind(size(ws),max(ix,1)):end))~=(dumvar8));dumvar12i=find((ws(sub2ind(size(ws),max(iw+1,1)):end))~=(dumvar12));   ws(ix-1+dumvar8i)=dumvar8(dumvar8i); ws(iw+1-1+dumvar12i)=dumvar12(dumvar12i); 
%
%     COMPUTE THE COMPONENTS OF THE SOLN DENOTED ABOVE BY W.
sc = one - sdot(m,a(sub2ind(size(a),1,np1):end),1,ws(sub2ind(size(ws),max(ix,1)):end),1);
if( one+fac.*abs(sc)==one || rnorm<=zero )
mode = 2;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
prgopt_shape=zeros(prgopt_shape);prgopt_shape(:)=prgopt(1:numel(prgopt_shape));prgopt=prgopt_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;
is_shape=zeros(is_shape);is_shape(:)=is(1:numel(is_shape));is=is_shape;
return;
else;
sc = one./sc;
for j = 1 : n1;
x(j) = sc.*sdot(m,a(sub2ind(size(a),1,j):end),1,ws(sub2ind(size(ws),max(ix,1)):end),1);
end; j = fix(n1+1);
%
%     COMPUTE THE VECTOR Q=Y-GW.  OVERWRITE Y WITH THIS VECTOR.
for i = 1 : m;
a(i,np1) = a(i,np1) - sdot(n1,a(sub2ind(size(a),i,1):end),mda,x,1);
end; i = fix(m+1);
end;
end;
if( n2>0 )
%
%     COPY TRANSPOSE OF (H Q) TO WORK ARRAY WS(*).
iw = 0;
for i = 1 : m;
[n2,a(sub2ind(size(a),i,n1+1):end),mda,ws(sub2ind(size(ws),max(iw+1,1)):end)]=scopy(n2,a(sub2ind(size(a),i,n1+1):end),mda,ws(sub2ind(size(ws),max(iw+1,1)):end),1);
iw = fix(iw + n2);
ws(iw+1) = a(i,np1);
iw = fix(iw + 1);
end; i = fix(m+1);
ws(iw+1) = zero;
[n2,dumvar2,dumvar3,dumvar2]=scopy(n2,ws(sub2ind(size(ws),max(iw+1,1)):end),0,ws(sub2ind(size(ws),max(iw+1,1)):end),1);   dumvar2i=find((ws(sub2ind(size(ws),max(iw+1,1)):end))~=(dumvar2));   ws(iw+1-1+dumvar2i)=dumvar2(dumvar2i); 
iw = fix(iw + n2);
ws(iw+1) = one;
iw = fix(iw + 1);
ix = fix(iw + 1);
iw = fix(iw + m);
%
%     SOLVE RV=S SUBJECT TO V.GE.0.  THE MATRIX R =(TRANSPOSE
%     OF (H Q)), WHERE Q=Y-GW.  THE (N2+1)-VECTOR S =(TRANSPOSE
%     OF (0,...,0,1)).
%
%     DO NOT CHECK LENGTHS OF WORK ARRAYS IN THIS USAGE OF WNNLS( ).
is(1) = 0;
is(2) = 0;
[ws,dumvar2,dumvar3,dumvar4,m,dumvar6,prgopt,dumvar8,rnorm,modew,is,dumvar12]=wnnls(ws,n2+1,0,n2+1,m,0,prgopt,ws(sub2ind(size(ws),max(ix,1)):end),rnorm,modew,is,ws(sub2ind(size(ws),max(iw+1,1)):end));   dumvar8i=find((ws(sub2ind(size(ws),max(ix,1)):end))~=(dumvar8));dumvar12i=find((ws(sub2ind(size(ws),max(iw+1,1)):end))~=(dumvar12));   ws(ix-1+dumvar8i)=dumvar8(dumvar8i); ws(iw+1-1+dumvar12i)=dumvar12(dumvar12i); 
%
%     COMPUTE THE COMPONENTS OF THE SOLN DENOTED ABOVE BY Z.
sc = one - sdot(m,a(sub2ind(size(a),1,np1):end),1,ws(sub2ind(size(ws),max(ix,1)):end),1);
if( one+fac.*abs(sc)==one || rnorm<=zero )
mode = 2;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
prgopt_shape=zeros(prgopt_shape);prgopt_shape(:)=prgopt(1:numel(prgopt_shape));prgopt=prgopt_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;
is_shape=zeros(is_shape);is_shape(:)=is(1:numel(is_shape));is=is_shape;
return;
else;
sc = one./sc;
for j = 1 : n2;
l = fix(n1 + j);
x(l) = sc.*sdot(m,a(sub2ind(size(a),1,l):end),1,ws(sub2ind(size(ws),max(ix,1)):end),1).*x(l);
end; j = fix(n2+1);
end;
end;
%
%     ACCOUNT FOR SCALING OF RT.-SIDE VECTOR IN SOLUTION.
[n,ynorm,x]=sscal(n,ynorm,x,1);
[wnorm ,n1,x]=snrm2(n1,x,1);
else;
if( n>0 )
x(1) = zero;
x_orig=x;    [n,x,dumvar3,dumvar4]=scopy(n,x,0,x,1);    x(dumvar4~=x_orig)=dumvar4(dumvar4~=x_orig);
end;
wnorm = zero;
end;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
prgopt_shape=zeros(prgopt_shape);prgopt_shape(:)=prgopt(1:numel(prgopt_shape));prgopt=prgopt_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;
is_shape=zeros(is_shape);is_shape(:)=is(1:numel(is_shape));is=is_shape;
end %subroutine lpdp
%DECK LSAME