*DECK LMPAR SUBROUTINE LMPAR (N, R, LDR, IPVT, DIAG, QTB, DELTA, PAR, X, + SIGMA, WA1, WA2) C***BEGIN PROLOGUE LMPAR C***SUBSIDIARY C***PURPOSE Subsidiary to SNLS1 and SNLS1E C***LIBRARY SLATEC C***TYPE SINGLE PRECISION (LMPAR-S, DMPAR-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, C the problem is to determine a value for the parameter C PAR such that if X solves the system C C A*X = B , SQRT(PAR)*D*X = 0 , C C in the least squares sense, and DXNORM is the Euclidean C norm of D*X, then either PAR is zero and C C (DXNORM-DELTA) .LE. 0.1*DELTA , C C or PAR is positive and C C ABS(DXNORM-DELTA) .LE. 0.1*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, with column pivoting, of A. That is, if C A*P = Q*R, where P is a permutation matrix, Q has orthogonal C columns, and R is an upper triangular matrix with diagonal C elements of nonincreasing magnitude, then LMPAR expects C the full upper triangle of R, the permutation matrix P, C and the first N components of (Q TRANSPOSE)*B. On output C LMPAR also provides an upper triangular matrix S such that C C T T T C P *(A *A + PAR*D*D)*P = S *S . C C S is employed within LMPAR and may be of separate interest. C C Only a few iterations are generally needed for convergence C of the algorithm. If, however, the limit of 10 iterations C is reached, then the output PAR will contain the best C value obtained so far. C C The subroutine statement is C C SUBROUTINE LMPAR(N,R,LDR,IPVT,DIAG,QTB,DELTA,PAR,X,SIGMA, C 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 N by N array. On input the full upper triangle C must contain the full upper triangle of the matrix R. C On output the full upper triangle is unaltered, and the C strict lower triangle contains the strict upper triangle C (transposed) of the upper triangular matrix S. C C LDR is a positive integer input variable not less than N C which specifies the leading dimension of the array R. C C IPVT is an integer input array of length N which defines the C permutation matrix P such that A*P = Q*R. Column J of P C is column IPVT(J) of the identity matrix. 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 PAR is a nonnegative variable. On input PAR contains an C initial estimate of the Levenberg-Marquardt parameter. C On output PAR contains the final estimate. C C X is an output array of length N which contains the least C squares solution of the system A*X = B, SQRT(PAR)*D*X = 0, C for the output PAR. C C SIGMA is an output array of length N which contains the C diagonal elements of the upper triangular matrix S. C C WA1 and WA2 are work arrays of length N. C C***SEE ALSO SNLS1, SNLS1E C***ROUTINES CALLED ENORM, QRSOLV, R1MACH 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 LMPAR INTEGER N,LDR INTEGER IPVT(*) REAL DELTA,PAR REAL R(LDR,*),DIAG(*),QTB(*),X(*),SIGMA(*),WA1(*),WA2(*) INTEGER I,ITER,J,JM1,JP1,K,L,NSING REAL DXNORM,DWARF,FP,GNORM,PARC,PARL,PARU,P1,P001,SUM,TEMP,ZERO REAL R1MACH,ENORM SAVE P1, P001, ZERO DATA P1,P001,ZERO /1.0E-1,1.0E-3,0.0E0/ C***FIRST EXECUTABLE STATEMENT LMPAR DWARF = R1MACH(1) C C COMPUTE AND STORE IN X THE GAUSS-NEWTON DIRECTION. IF THE C JACOBIAN IS RANK-DEFICIENT, OBTAIN A LEAST SQUARES SOLUTION. C NSING = N DO 10 J = 1, N WA1(J) = QTB(J) IF (R(J,J) .EQ. ZERO .AND. NSING .EQ. N) NSING = J - 1 IF (NSING .LT. N) WA1(J) = ZERO 10 CONTINUE IF (NSING .LT. 1) GO TO 50 DO 40 K = 1, NSING J = NSING - K + 1 WA1(J) = WA1(J)/R(J,J) TEMP = WA1(J) JM1 = J - 1 IF (JM1 .LT. 1) GO TO 30 DO 20 I = 1, JM1 WA1(I) = WA1(I) - R(I,J)*TEMP 20 CONTINUE 30 CONTINUE 40 CONTINUE 50 CONTINUE DO 60 J = 1, N L = IPVT(J) X(L) = WA1(J) 60 CONTINUE C C INITIALIZE THE ITERATION COUNTER. C EVALUATE THE FUNCTION AT THE ORIGIN, AND TEST C FOR ACCEPTANCE OF THE GAUSS-NEWTON DIRECTION. C ITER = 0 DO 70 J = 1, N WA2(J) = DIAG(J)*X(J) 70 CONTINUE DXNORM = ENORM(N,WA2) FP = DXNORM - DELTA IF (FP .LE. P1*DELTA) GO TO 220 C C IF THE JACOBIAN IS NOT RANK DEFICIENT, THE NEWTON C STEP PROVIDES A LOWER BOUND, PARL, FOR THE ZERO OF C THE FUNCTION. OTHERWISE SET THIS BOUND TO ZERO. C PARL = ZERO IF (NSING .LT. N) GO TO 120 DO 80 J = 1, N L = IPVT(J) WA1(J) = DIAG(L)*(WA2(L)/DXNORM) 80 CONTINUE DO 110 J = 1, N SUM = ZERO JM1 = J - 1 IF (JM1 .LT. 1) GO TO 100 DO 90 I = 1, JM1 SUM = SUM + R(I,J)*WA1(I) 90 CONTINUE 100 CONTINUE WA1(J) = (WA1(J) - SUM)/R(J,J) 110 CONTINUE TEMP = ENORM(N,WA1) PARL = ((FP/DELTA)/TEMP)/TEMP 120 CONTINUE C C CALCULATE AN UPPER BOUND, PARU, FOR THE ZERO OF THE FUNCTION. C DO 140 J = 1, N SUM = ZERO DO 130 I = 1, J SUM = SUM + R(I,J)*QTB(I) 130 CONTINUE L = IPVT(J) WA1(J) = SUM/DIAG(L) 140 CONTINUE GNORM = ENORM(N,WA1) PARU = GNORM/DELTA IF (PARU .EQ. ZERO) PARU = DWARF/MIN(DELTA,P1) C C IF THE INPUT PAR LIES OUTSIDE OF THE INTERVAL (PARL,PARU), C SET PAR TO THE CLOSER ENDPOINT. C PAR = MAX(PAR,PARL) PAR = MIN(PAR,PARU) IF (PAR .EQ. ZERO) PAR = GNORM/DXNORM C C BEGINNING OF AN ITERATION. C 150 CONTINUE ITER = ITER + 1 C C EVALUATE THE FUNCTION AT THE CURRENT VALUE OF PAR. C IF (PAR .EQ. ZERO) PAR = MAX(DWARF,P001*PARU) TEMP = SQRT(PAR) DO 160 J = 1, N WA1(J) = TEMP*DIAG(J) 160 CONTINUE CALL QRSOLV(N,R,LDR,IPVT,WA1,QTB,X,SIGMA,WA2) DO 170 J = 1, N WA2(J) = DIAG(J)*X(J) 170 CONTINUE DXNORM = ENORM(N,WA2) TEMP = FP FP = DXNORM - DELTA C C IF THE FUNCTION IS SMALL ENOUGH, ACCEPT THE CURRENT VALUE C OF PAR. ALSO TEST FOR THE EXCEPTIONAL CASES WHERE PARL C IS ZERO OR THE NUMBER OF ITERATIONS HAS REACHED 10. C IF (ABS(FP) .LE. P1*DELTA 1 .OR. PARL .EQ. ZERO .AND. FP .LE. TEMP 2 .AND. TEMP .LT. ZERO .OR. ITER .EQ. 10) GO TO 220 C C COMPUTE THE NEWTON CORRECTION. C DO 180 J = 1, N L = IPVT(J) WA1(J) = DIAG(L)*(WA2(L)/DXNORM) 180 CONTINUE DO 210 J = 1, N WA1(J) = WA1(J)/SIGMA(J) TEMP = WA1(J) JP1 = J + 1 IF (N .LT. JP1) GO TO 200 DO 190 I = JP1, N WA1(I) = WA1(I) - R(I,J)*TEMP 190 CONTINUE 200 CONTINUE 210 CONTINUE TEMP = ENORM(N,WA1) PARC = ((FP/DELTA)/TEMP)/TEMP C C DEPENDING ON THE SIGN OF THE FUNCTION, UPDATE PARL OR PARU. C IF (FP .GT. ZERO) PARL = MAX(PARL,PAR) IF (FP .LT. ZERO) PARU = MIN(PARU,PAR) C C COMPUTE AN IMPROVED ESTIMATE FOR PAR. C PAR = MAX(PARL,PAR+PARC) C C END OF AN ITERATION. C GO TO 150 220 CONTINUE C C TERMINATION. C IF (ITER .EQ. 0) PAR = ZERO RETURN C C LAST CARD OF SUBROUTINE LMPAR. C END