[n,r,ldr,ipvt,diag,qtb,delta,par,x,sigma,wa1,wa2]=lmpar(n,r,ldr,ipvt,diag,qtb,delta,par,x,sigma,wa1,wa2);

function [n,r,ldr,ipvt,diag,qtb,delta,par,x,sigma,wa1,wa2]=lmpar(n,r,ldr,ipvt,diag,qtb,delta,par,x,sigma,wa1,wa2);
%***BEGIN PROLOGUE  LMPAR
%***SUBSIDIARY
%***PURPOSE  Subsidiary to SNLS1 and SNLS1E
%***LIBRARY   SLATEC
%***TYPE      SINGLE PRECISION (LMPAR-S, DMPAR-D)
%***AUTHOR  (UNKNOWN)
%***DESCRIPTION
%
%     Given an M by N matrix A, an N by N nonsingular DIAGONAL
%     matrix D, an M-vector B, and a positive number DELTA,
%     the problem is to determine a value for the parameter
%     PAR such that if X solves the system
%
%           A*X = B ,     SQRT(PAR)*D*X = 0 ,
%
%     in the least squares sense, and DXNORM is the Euclidean
%     norm of D*X, then either PAR is zero and
%
%           (DXNORM-DELTA) .LE. 0.1*DELTA ,
%
%     or PAR is positive and
%
%           ABS(DXNORM-DELTA) .LE. 0.1*DELTA .
%
%     This subroutine completes the solution of the problem
%     if it is provided with the necessary information from the
%     QR factorization, with column pivoting, of A. That is, if
%     A*P = Q*R, where P is a permutation matrix, Q has orthogonal
%     columns, and R is an upper triangular matrix with diagonal
%     elements of nonincreasing magnitude, then LMPAR expects
%     the full upper triangle of R, the permutation matrix P,
%     and the first N components of (Q TRANSPOSE)*B. On output
%     LMPAR also provides an upper triangular matrix S such that
%
%            T   T                   T
%           P *(A *A + PAR*D*D)*P = S *S .
%
%     S is employed within LMPAR and may be of separate interest.
%
%     Only a few iterations are generally needed for convergence
%     of the algorithm. If, however, the limit of 10 iterations
%     is reached, then the output PAR will contain the best
%     value obtained so far.
%
%     The subroutine statement is
%
%       subroutine LMPAR(N,R,LDR,IPVT,DIAG,QTB,DELTA,PAR,X,SIGMA,
%                        WA1,WA2)
%
%     where
%
%       N is a positive integer input variable set to the order of R.
%
%       R is an N by N array. On input the full upper triangle
%         must contain the full upper triangle of the matrix R.
%         On output the full upper triangle is unaltered, and the
%         strict lower triangle contains the strict upper triangle
%         (transposed) of the upper triangular matrix S.
%
%       LDR is a positive integer input variable not less than N
%         which specifies the leading dimension of the array R.
%
%       IPVT is an integer input array of length N which defines the
%         permutation matrix P such that A*P = Q*R. Column J of P
%         is column IPVT(J) of the identity matrix.
%
%       DIAG is an input array of length N which must contain the
%         diagonal elements of the matrix D.
%
%       QTB is an input array of length N which must contain the first
%         N elements of the vector (Q TRANSPOSE)*B.
%
%       DELTA is a positive input variable which specifies an upper
%         bound on the Euclidean norm of D*X.
%
%       PAR is a nonnegative variable. On input PAR contains an
%         initial estimate of the Levenberg-Marquardt parameter.
%         On output PAR contains the final estimate.
%
%       X is an output array of length N which contains the least
%         squares solution of the system A*X = B, SQRT(PAR)*D*X = 0,
%         for the output PAR.
%
%       SIGMA is an output array of length N which contains the
%         diagonal elements of the upper triangular matrix S.
%
%       WA1 and WA2 are work arrays of length N.
%
%***SEE ALSO  SNLS1, SNLS1E
%***ROUTINES CALLED  ENORM, QRSOLV, R1MACH
%***REVISION HISTORY  (YYMMDD)
%   800301  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890831  Modified array declarations.  (WRB)
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900326  Removed duplicate information from DESCRIPTION section.
%           (WRB)
%   900328  Added TYPE section.  (WRB)
%***end PROLOGUE  LMPAR
persistent dwarf dxnorm firstCall fp gnorm i iter j jm1 jp1 k l nsing p001 p1 parc parl paru summlv temp zero ; if isempty(firstCall),firstCall=1;end; 

ipvt_shape=size(ipvt);ipvt=reshape(ipvt,1,[]);
r_shape=size(r);r=reshape([r(:).',zeros(1,ceil(numel(r)./prod([ldr])).*prod([ldr])-numel(r))],ldr,[]);
diag_shape=size(diag);diag=reshape(diag,1,[]);
qtb_shape=size(qtb);qtb=reshape(qtb,1,[]);
x_shape=size(x);x=reshape(x,1,[]);
sigma_shape=size(sigma);sigma=reshape(sigma,1,[]);
wa1_shape=size(wa1);wa1=reshape(wa1,1,[]);
wa2_shape=size(wa2);wa2=reshape(wa2,1,[]);
if isempty(i), i=0; end;
if isempty(iter), iter=0; end;
if isempty(j), j=0; end;
if isempty(jm1), jm1=0; end;
if isempty(jp1), jp1=0; end;
if isempty(k), k=0; end;
if isempty(l), l=0; end;
if isempty(nsing), nsing=0; end;
if isempty(dxnorm), dxnorm=0; end;
if isempty(dwarf), dwarf=0; end;
if isempty(fp), fp=0; end;
if isempty(gnorm), gnorm=0; end;
if isempty(parc), parc=0; end;
if isempty(parl), parl=0; end;
if isempty(paru), paru=0; end;
if isempty(p1), p1=0; end;
if isempty(p001), p001=0; end;
if isempty(summlv), summlv=0; end;
if isempty(temp), temp=0; end;
if isempty(zero), zero=0; end;
if firstCall,   p1 =[1.0e-1];  end;
if firstCall,  p001 =[1.0e-3];  end;
if firstCall,  zero=[0.0e0];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  LMPAR
[dwarf ]=r1mach(1);
%
%     COMPUTE AND STORE IN X THE GAUSS-NEWTON DIRECTION. IF THE
%     JACOBIAN IS RANK-DEFICIENT, OBTAIN A LEAST SQUARES SOLUTION.
%
nsing = fix(n);
for j = 1 : n;
wa1(j) = qtb(j);
if( r(j,j)==zero && nsing==n )
nsing = fix(j - 1);
end;
if( nsing<n )
wa1(j) = zero;
end;
end; j = fix(n+1);
if( nsing>=1 )
for k = 1 : nsing;
j = fix(nsing - k + 1);
wa1(j) = wa1(j)./r(j,j);
temp = wa1(j);
jm1 = fix(j - 1);
if( jm1>=1 )
for i = 1 : jm1;
wa1(i) = wa1(i) - r(i,j).*temp;
end; i = fix(jm1+1);
end;
end; k = fix(nsing+1);
end;
for j = 1 : n;
l = fix(ipvt(j));
x(l) = wa1(j);
end; j = fix(n+1);
%
%     INITIALIZE THE ITERATION COUNTER.
%     EVALUATE THE FUNCTION AT THE ORIGIN, AND TEST
%     FOR ACCEPTANCE OF THE GAUSS-NEWTON DIRECTION.
%
iter = 0;
for j = 1 : n;
wa2(j) = diag(j).*x(j);
end; j = fix(n+1);
[dxnorm ,n,wa2]=enorm(n,wa2);
fp = dxnorm - delta;
if( fp<=p1.*delta )
%
%     TERMINATION.
%
if( iter==0 )
par = zero;
end;
else;
%
%     IF THE JACOBIAN IS NOT RANK DEFICIENT, THE NEWTON
%     STEP PROVIDES A LOWER BOUND, PARL, FOR THE ZERO OF
%     THE FUNCTION. OTHERWISE SET THIS BOUND TO ZERO.
%
parl = zero;
if( nsing>=n )
for j = 1 : n;
l = fix(ipvt(j));
wa1(j) = diag(l).*(wa2(l)./dxnorm);
end; j = fix(n+1);
for j = 1 : n;
summlv = zero;
jm1 = fix(j - 1);
if( jm1>=1 )
for i = 1 : jm1;
summlv = summlv + r(i,j).*wa1(i);
end; i = fix(jm1+1);
end;
wa1(j) =(wa1(j)-summlv)./r(j,j);
end; j = fix(n+1);
[temp ,n,wa1]=enorm(n,wa1);
parl =((fp./delta)./temp)./temp;
end;
%
%     CALCULATE AN UPPER BOUND, PARU, FOR THE ZERO OF THE FUNCTION.
%
for j = 1 : n;
summlv = zero;
for i = 1 : j;
summlv = summlv + r(i,j).*qtb(i);
end; i = fix(j+1);
l = fix(ipvt(j));
wa1(j) = summlv./diag(l);
end; j = fix(n+1);
[gnorm ,n,wa1]=enorm(n,wa1);
paru = gnorm./delta;
if( paru==zero )
paru = dwarf./min(delta,p1);
end;
%
%     IF THE INPUT PAR LIES OUTSIDE OF THE INTERVAL (PARL,PARU),
%     SET PAR TO THE CLOSER ENDPOINT.
%
par = max(par,parl);
par = min(par,paru);
if( par==zero )
par = gnorm./dxnorm;
end;
%
%     BEGINNING OF AN ITERATION.
%
while( true );
iter = fix(iter + 1);
%
%        EVALUATE THE FUNCTION AT THE CURRENT VALUE OF PAR.
%
if( par==zero )
par = max(dwarf,p001.*paru);
end;
temp = sqrt(par);
for j = 1 : n;
wa1(j) = temp.*diag(j);
end; j = fix(n+1);
[n,r,ldr,ipvt,wa1,qtb,x,sigma,wa2]=qrsolv(n,r,ldr,ipvt,wa1,qtb,x,sigma,wa2);
for j = 1 : n;
wa2(j) = diag(j).*x(j);
end; j = fix(n+1);
[dxnorm ,n,wa2]=enorm(n,wa2);
temp = fp;
fp = dxnorm - delta;
%
%        IF THE FUNCTION IS SMALL ENOUGH, ACCEPT THE CURRENT VALUE
%        OF PAR. ALSO TEST FOR THE EXCEPTIONAL CASES WHERE PARL
%        IS ZERO OR THE NUMBER OF ITERATIONS HAS REACHED 10.
%
if( abs(fp)<=p1.*delta || parl==zero && fp<=temp &&temp<zero || iter==10 )
break;
end;
%
%        COMPUTE THE NEWTON CORRECTION.
%
for j = 1 : n;
l = fix(ipvt(j));
wa1(j) = diag(l).*(wa2(l)./dxnorm);
end; j = fix(n+1);
for j = 1 : n;
wa1(j) = wa1(j)./sigma(j);
temp = wa1(j);
jp1 = fix(j + 1);
if( n>=jp1 )
for i = jp1 : n;
wa1(i) = wa1(i) - r(i,j).*temp;
end; i = fix(n+1);
end;
end; j = fix(n+1);
[temp ,n,wa1]=enorm(n,wa1);
parc =((fp./delta)./temp)./temp;
%
%        DEPENDING ON THE SIGN OF THE FUNCTION, UPDATE PARL OR PARU.
%
if( fp>zero )
parl = max(parl,par);
end;
if( fp<zero )
paru = min(paru,par);
end;
%
%        COMPUTE AN IMPROVED ESTIMATE FOR PAR.
%
par = max(parl,par+parc);
%
%        end OF AN ITERATION.
%
end;
if( iter==0 )
par = zero;
end;
end;
%
%     LAST CARD OF SUBROUTINE LMPAR.
%
ipvt_shape=zeros(ipvt_shape);ipvt_shape(:)=ipvt(1:numel(ipvt_shape));ipvt=ipvt_shape;
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_shape;
diag_shape=zeros(diag_shape);diag_shape(:)=diag(1:numel(diag_shape));diag=diag_shape;
qtb_shape=zeros(qtb_shape);qtb_shape(:)=qtb(1:numel(qtb_shape));qtb=qtb_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
sigma_shape=zeros(sigma_shape);sigma_shape(:)=sigma(1:numel(sigma_shape));sigma=sigma_shape;
wa1_shape=zeros(wa1_shape);wa1_shape(:)=wa1(1:numel(wa1_shape));wa1=wa1_shape;
wa2_shape=zeros(wa2_shape);wa2_shape(:)=wa2(1:numel(wa2_shape));wa2=wa2_shape;
end %subroutine lmpar
%DECK LPDP