*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