*DECK DDOGLG
SUBROUTINE DDOGLG (N, R, LR, DIAG, QTB, DELTA, X, WA1, WA2)
C***BEGIN PROLOGUE DDOGLG
C***SUBSIDIARY
C***PURPOSE Subsidiary to DNSQ and DNSQE
C***LIBRARY SLATEC
C***TYPE DOUBLE PRECISION (DOGLEG-S, DDOGLG-D)
C***AUTHOR (UNKNOWN)
C***DESCRIPTION
C
C Given an M by N matrix A, an N by N nonsingular diagonal
C matrix D, an M-vector B, and a positive number DELTA, the
C problem is to determine the convex combination X of the
C Gauss-Newton and scaled gradient directions that minimizes
C (A*X - B) in the least squares sense, subject to the
C restriction that the Euclidean norm of D*X be at most DELTA.
C
C This subroutine completes the solution of the problem
C if it is provided with the necessary information from the
C QR factorization of A. That is, if A = Q*R, where Q has
C orthogonal columns and R is an upper triangular matrix,
C then DDOGLG expects the full upper triangle of R and
C the first N components of (Q transpose)*B.
C
C The subroutine statement is
C
C SUBROUTINE DDOGLG(N,R,LR,DIAG,QTB,DELTA,X,WA1,WA2)
C
C where
C
C N is a positive integer input variable set to the order of R.
C
C R is an input array of length LR which must contain the upper
C triangular matrix R stored by rows.
C
C LR is a positive integer input variable not less than
C (N*(N+1))/2.
C
C DIAG is an input array of length N which must contain the
C diagonal elements of the matrix D.
C
C QTB is an input array of length N which must contain the first
C N elements of the vector (Q transpose)*B.
C
C DELTA is a positive input variable which specifies an upper
C bound on the Euclidean norm of D*X.
C
C X is an output array of length N which contains the desired
C convex combination of the Gauss-Newton direction and the
C scaled gradient direction.
C
C WA1 and WA2 are work arrays of length N.
C
C***SEE ALSO DNSQ, DNSQE
C***ROUTINES CALLED D1MACH, DENORM
C***REVISION HISTORY (YYMMDD)
C 800301 DATE WRITTEN
C 890531 Changed all specific intrinsics to generic. (WRB)
C 890831 Modified array declarations. (WRB)
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 900326 Removed duplicate information from DESCRIPTION section.
C (WRB)
C 900328 Added TYPE section. (WRB)
C***END PROLOGUE DDOGLG
DOUBLE PRECISION D1MACH,DENORM
INTEGER I, J, JJ, JP1, K, L, LR, N
DOUBLE PRECISION ALPHA, BNORM, DELTA, DIAG(*), EPSMCH, GNORM,
1 ONE, QNORM, QTB(*), R(*), SGNORM, SUM, TEMP, WA1(*),
2 WA2(*), X(*), ZERO
SAVE ONE, ZERO
DATA ONE,ZERO /1.0D0,0.0D0/
C
C EPSMCH IS THE MACHINE PRECISION.
C
C***FIRST EXECUTABLE STATEMENT DDOGLG
EPSMCH = D1MACH(4)
C
C FIRST, CALCULATE THE GAUSS-NEWTON DIRECTION.
C
JJ = (N*(N + 1))/2 + 1
DO 50 K = 1, N
J = N - K + 1
JP1 = J + 1
JJ = JJ - K
L = JJ + 1
SUM = ZERO
IF (N .LT. JP1) GO TO 20
DO 10 I = JP1, N
SUM = SUM + R(L)*X(I)
L = L + 1
10 CONTINUE
20 CONTINUE
TEMP = R(JJ)
IF (TEMP .NE. ZERO) GO TO 40
L = J
DO 30 I = 1, J
TEMP = MAX(TEMP,ABS(R(L)))
L = L + N - I
30 CONTINUE
TEMP = EPSMCH*TEMP
IF (TEMP .EQ. ZERO) TEMP = EPSMCH
40 CONTINUE
X(J) = (QTB(J) - SUM)/TEMP
50 CONTINUE
C
C TEST WHETHER THE GAUSS-NEWTON DIRECTION IS ACCEPTABLE.
C
DO 60 J = 1, N
WA1(J) = ZERO
WA2(J) = DIAG(J)*X(J)
60 CONTINUE
QNORM = DENORM(N,WA2)
IF (QNORM .LE. DELTA) GO TO 140
C
C THE GAUSS-NEWTON DIRECTION IS NOT ACCEPTABLE.
C NEXT, CALCULATE THE SCALED GRADIENT DIRECTION.
C
L = 1
DO 80 J = 1, N
TEMP = QTB(J)
DO 70 I = J, N
WA1(I) = WA1(I) + R(L)*TEMP
L = L + 1
70 CONTINUE
WA1(J) = WA1(J)/DIAG(J)
80 CONTINUE
C
C CALCULATE THE NORM OF THE SCALED GRADIENT AND TEST FOR
C THE SPECIAL CASE IN WHICH THE SCALED GRADIENT IS ZERO.
C
GNORM = DENORM(N,WA1)
SGNORM = ZERO
ALPHA = DELTA/QNORM
IF (GNORM .EQ. ZERO) GO TO 120
C
C CALCULATE THE POINT ALONG THE SCALED GRADIENT
C AT WHICH THE QUADRATIC IS MINIMIZED.
C
DO 90 J = 1, N
WA1(J) = (WA1(J)/GNORM)/DIAG(J)
90 CONTINUE
L = 1
DO 110 J = 1, N
SUM = ZERO
DO 100 I = J, N
SUM = SUM + R(L)*WA1(I)
L = L + 1
100 CONTINUE
WA2(J) = SUM
110 CONTINUE
TEMP = DENORM(N,WA2)
SGNORM = (GNORM/TEMP)/TEMP
C
C TEST WHETHER THE SCALED GRADIENT DIRECTION IS ACCEPTABLE.
C
ALPHA = ZERO
IF (SGNORM .GE. DELTA) GO TO 120
C
C THE SCALED GRADIENT DIRECTION IS NOT ACCEPTABLE.
C FINALLY, CALCULATE THE POINT ALONG THE DOGLEG
C AT WHICH THE QUADRATIC IS MINIMIZED.
C
BNORM = DENORM(N,QTB)
TEMP = (BNORM/GNORM)*(BNORM/QNORM)*(SGNORM/DELTA)
TEMP = TEMP - (DELTA/QNORM)*(SGNORM/DELTA)**2
1 + SQRT((TEMP-(DELTA/QNORM))**2
2 +(ONE-(DELTA/QNORM)**2)*(ONE-(SGNORM/DELTA)**2))
ALPHA = ((DELTA/QNORM)*(ONE - (SGNORM/DELTA)**2))/TEMP
120 CONTINUE
C
C FORM APPROPRIATE CONVEX COMBINATION OF THE GAUSS-NEWTON
C DIRECTION AND THE SCALED GRADIENT DIRECTION.
C
TEMP = (ONE - ALPHA)*MIN(SGNORM,DELTA)
DO 130 J = 1, N
X(J) = TEMP*WA1(J) + ALPHA*X(J)
130 CONTINUE
140 CONTINUE
RETURN
C
C LAST CARD OF SUBROUTINE DDOGLG.
C
END