[n,r,lr,diag,qtb,delta,x,wa1,wa2]=ddoglg(n,r,lr,diag,qtb,delta,x,wa1,wa2);

function [n,r,lr,diag,qtb,delta,x,wa1,wa2]=ddoglg(n,r,lr,diag,qtb,delta,x,wa1,wa2);
%***BEGIN PROLOGUE  DDOGLG
%***SUBSIDIARY
%***PURPOSE  Subsidiary to DNSQ and DNSQE
%***LIBRARY   SLATEC
%***TYPE      doubleprecision (DOGLEG-S, DDOGLG-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 the convex combination X of the
%     Gauss-Newton and scaled gradient directions that minimizes
%     (A*X - B) in the least squares sense, subject to the
%     restriction that the Euclidean norm of D*X be at most DELTA.
%
%     This subroutine completes the solution of the problem
%     if it is provided with the necessary information from the
%     QR factorization of A. That is, if A = Q*R, where Q has
%     orthogonal columns and R is an upper triangular matrix,
%     then DDOGLG expects the full upper triangle of R and
%     the first N components of (Q transpose)*B.
%
%     The subroutine statement is
%
%       subroutine DDOGLG(N,R,LR,DIAG,QTB,DELTA,X,WA1,WA2)
%
%     where
%
%       N is a positive integer input variable set to the order of R.
%
%       R is an input array of length LR which must contain the upper
%         triangular matrix R stored by rows.
%
%       LR is a positive integer input variable not less than
%         (N*(N+1))/2.
%
%       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.
%
%       X is an output array of length N which contains the desired
%         convex combination of the Gauss-Newton direction and the
%         scaled gradient direction.
%
%       WA1 and WA2 are work arrays of length N.
%
%***SEE ALSO  DNSQ, DNSQE
%***ROUTINES CALLED  D1MACH, DENORM
%***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  DDOGLG
persistent alpha bnorm epsmch firstCall gnorm i j jj jp1 k l one qnorm sgnorm summlv temp zero ; if isempty(firstCall),firstCall=1;end; 

if isempty(i), i=0; end;
if isempty(j), j=0; end;
if isempty(jj), jj=0; end;
if isempty(jp1), jp1=0; end;
if isempty(k), k=0; end;
if isempty(l), l=0; end;
if isempty(alpha), alpha=0; end;
if isempty(bnorm), bnorm=0; end;
diag_shape=size(diag);diag=reshape(diag,1,[]);
if isempty(epsmch), epsmch=0; end;
if isempty(gnorm), gnorm=0; end;
if isempty(one), one=0; end;
if isempty(qnorm), qnorm=0; end;
qtb_shape=size(qtb);qtb=reshape(qtb,1,[]);
r_shape=size(r);r=reshape(r,1,[]);
if isempty(sgnorm), sgnorm=0; end;
if isempty(summlv), summlv=0; end;
if isempty(temp), temp=0; end;
wa1_shape=size(wa1);wa1=reshape(wa1,1,[]);
wa2_shape=size(wa2);wa2=reshape(wa2,1,[]);
x_shape=size(x);x=reshape(x,1,[]);
if isempty(zero), zero=0; end;
if firstCall,   one =[1.0d0];  end;
if firstCall,  zero=[0.0d0];  end;
firstCall=0;
%
%     EPSMCH IS THE MACHINE PRECISION.
%
%***FIRST EXECUTABLE STATEMENT  DDOGLG
[epsmch ]=d1mach(4);
%
%     FIRST, CALCULATE THE GAUSS-NEWTON DIRECTION.
%
jj =fix(fix((n.*(n+1))./2) + 1);
for k = 1 : n;
j = fix(n - k + 1);
jp1 = fix(j + 1);
jj = fix(jj - k);
l = fix(jj + 1);
summlv = zero;
if( n>=jp1 )
for i = jp1 : n;
summlv = summlv + r(l).*x(i);
l = fix(l + 1);
end; i = fix(n+1);
end;
temp = r(jj);
if( temp==zero )
l = fix(j);
for i = 1 : j;
temp = max(temp,abs(r(l)));
l = fix(l + n - i);
end; i = fix(j+1);
temp = epsmch.*temp;
if( temp==zero )
temp = epsmch;
end;
end;
x(j) =(qtb(j)-summlv)./temp;
end; k = fix(n+1);
%
%     TEST WHETHER THE GAUSS-NEWTON DIRECTION IS ACCEPTABLE.
%
for j = 1 : n;
wa1(j) = zero;
wa2(j) = diag(j).*x(j);
end; j = fix(n+1);
[qnorm ,n,wa2]=denorm(n,wa2);
if( qnorm>delta )
%
%     THE GAUSS-NEWTON DIRECTION IS NOT ACCEPTABLE.
%     NEXT, CALCULATE THE SCALED GRADIENT DIRECTION.
%
l = 1;
for j = 1 : n;
temp = qtb(j);
for i = j : n;
wa1(i) = wa1(i) + r(l).*temp;
l = fix(l + 1);
end; i = fix(n+1);
wa1(j) = wa1(j)./diag(j);
end; j = fix(n+1);
%
%     CALCULATE THE NORM OF THE SCALED GRADIENT AND TEST FOR
%     THE SPECIAL CASE IN WHICH THE SCALED GRADIENT IS ZERO.
%
[gnorm ,n,wa1]=denorm(n,wa1);
sgnorm = zero;
alpha = delta./qnorm;
if( gnorm~=zero )
%
%     CALCULATE THE POINT ALONG THE SCALED GRADIENT
%     AT WHICH THE QUADRATIC IS MINIMIZED.
%
for j = 1 : n;
wa1(j) =(wa1(j)./gnorm)./diag(j);
end; j = fix(n+1);
l = 1;
for j = 1 : n;
summlv = zero;
for i = j : n;
summlv = summlv + r(l).*wa1(i);
l = fix(l + 1);
end; i = fix(n+1);
wa2(j) = summlv;
end; j = fix(n+1);
[temp ,n,wa2]=denorm(n,wa2);
sgnorm =(gnorm./temp)./temp;
%
%     TEST WHETHER THE SCALED GRADIENT DIRECTION IS ACCEPTABLE.
%
alpha = zero;
if( sgnorm<delta )
%
%     THE SCALED GRADIENT DIRECTION IS NOT ACCEPTABLE.
%     FINALLY, CALCULATE THE POINT ALONG THE DOGLEG
%     AT WHICH THE QUADRATIC IS MINIMIZED.
%
[bnorm ,n,qtb]=denorm(n,qtb);
temp =(bnorm./gnorm).*(bnorm./qnorm).*(sgnorm./delta);
temp = temp -(delta./qnorm).*(sgnorm./delta).^2 + sqrt((temp-(delta./qnorm)).^2+(one-(delta./qnorm).^2).*(one-(sgnorm./delta).^2));
alpha =((delta./qnorm).*(one-(sgnorm./delta).^2))./temp;
end;
end;
%
%     FORM APPROPRIATE CONVEX COMBINATION OF THE GAUSS-NEWTON
%     DIRECTION AND THE SCALED GRADIENT DIRECTION.
%
temp =(one-alpha).*min(sgnorm,delta);
for j = 1 : n;
x(j) = temp.*wa1(j) + alpha.*x(j);
end; j = fix(n+1);
end;
%
%     LAST CARD OF SUBROUTINE DDOGLG.
%
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;
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_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;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
end
%DECK DDOT