2 SUBROUTINE LMPAR (N, R, LDR, IPVT, DIAG, QTB, DELTA, PAR, X,
4 C***BEGIN PROLOGUE LMPAR
6 C***PURPOSE Subsidiary to SNLS1 and SNLS1E
8 C***TYPE SINGLE PRECISION (LMPAR-S, DMPAR-D)
12 C Given an M by N matrix A, an N by N nonsingular DIAGONAL
13 C matrix D, an M-vector B, and a positive number DELTA,
14 C the problem is to determine a value for the parameter
15 C PAR such that if X solves the system
17 C A*X = B , SQRT(PAR)*D*X = 0 ,
19 C in the least squares sense, and DXNORM is the Euclidean
20 C norm of D*X, then either PAR is zero and
22 C (DXNORM-DELTA) .LE. 0.1*DELTA ,
24 C or PAR is positive and
26 C ABS(DXNORM-DELTA) .LE. 0.1*DELTA .
28 C This subroutine completes the solution of the problem
29 C if it is provided with the necessary information from the
30 C QR factorization, with column pivoting, of A. That is, if
31 C A*P = Q*R, where P is a permutation matrix, Q has orthogonal
32 C columns, and R is an upper triangular matrix with diagonal
33 C elements of nonincreasing magnitude, then LMPAR expects
34 C the full upper triangle of R, the permutation matrix P,
35 C and the first N components of (Q TRANSPOSE)*B. On output
36 C LMPAR also provides an upper triangular matrix S such that
39 C P *(A *A + PAR*D*D)*P = S *S .
41 C S is employed within LMPAR and may be of separate interest.
43 C Only a few iterations are generally needed for convergence
44 C of the algorithm. If, however, the limit of 10 iterations
45 C is reached, then the output PAR will contain the best
46 C value obtained so far.
48 C The subroutine statement is
50 C SUBROUTINE LMPAR(N,R,LDR,IPVT,DIAG,QTB,DELTA,PAR,X,SIGMA,
55 C N is a positive integer input variable set to the order of R.
57 C R is an N by N array. On input the full upper triangle
58 C must contain the full upper triangle of the matrix R.
59 C On output the full upper triangle is unaltered, and the
60 C strict lower triangle contains the strict upper triangle
61 C (transposed) of the upper triangular matrix S.
63 C LDR is a positive integer input variable not less than N
64 C which specifies the leading dimension of the array R.
66 C IPVT is an integer input array of length N which defines the
67 C permutation matrix P such that A*P = Q*R. Column J of P
68 C is column IPVT(J) of the identity matrix.
70 C DIAG is an input array of length N which must contain the
71 C diagonal elements of the matrix D.
73 C QTB is an input array of length N which must contain the first
74 C N elements of the vector (Q TRANSPOSE)*B.
76 C DELTA is a positive input variable which specifies an upper
77 C bound on the Euclidean norm of D*X.
79 C PAR is a nonnegative variable. On input PAR contains an
80 C initial estimate of the Levenberg-Marquardt parameter.
81 C On output PAR contains the final estimate.
83 C X is an output array of length N which contains the least
84 C squares solution of the system A*X = B, SQRT(PAR)*D*X = 0,
87 C SIGMA is an output array of length N which contains the
88 C diagonal elements of the upper triangular matrix S.
90 C WA1 and WA2 are work arrays of length N.
92 C***SEE ALSO SNLS1, SNLS1E
93 C***ROUTINES CALLED ENORM, QRSOLV, R1MACH
94 C***REVISION HISTORY (YYMMDD)
96 C 890531 Changed all specific intrinsics to generic. (WRB)
97 C 890831 Modified array declarations. (WRB)
98 C 891214 Prologue converted to Version 4.0 format. (BAB)
99 C 900326 Removed duplicate information from DESCRIPTION section.
101 C 900328 Added TYPE section. (WRB)
102 C***END PROLOGUE LMPAR
106 REAL R(LDR,*),DIAG(*),QTB(*),X(*),SIGMA(*),WA1(*),WA2(*)
107 INTEGER I,ITER,J,JM1,JP1,K,L,NSING
108 REAL DXNORM,DWARF,FP,GNORM,PARC,PARL,PARU,P1,P001,SUM,TEMP,ZERO
111 DATA P1,P001,ZERO /1.0E-1,1.0E-3,0.0E0/
112 C***FIRST EXECUTABLE STATEMENT LMPAR
115 C COMPUTE AND STORE IN X THE GAUSS-NEWTON DIRECTION. IF THE
116 C JACOBIAN IS RANK-DEFICIENT, OBTAIN A LEAST SQUARES SOLUTION.
121 IF (R(J,J) .EQ. ZERO .AND. NSING .EQ. N) NSING = J - 1
122 IF (NSING .LT. N) WA1(J) = ZERO
124 IF (NSING .LT. 1) GO TO 50
127 WA1(J) = WA1(J)/R(J,J)
130 IF (JM1 .LT. 1) GO TO 30
132 WA1(I) = WA1(I) - R(I,J)*TEMP
142 C INITIALIZE THE ITERATION COUNTER.
143 C EVALUATE THE FUNCTION AT THE ORIGIN, AND TEST
144 C FOR ACCEPTANCE OF THE GAUSS-NEWTON DIRECTION.
148 WA2(J) = DIAG(J)*X(J)
150 DXNORM = ENORM(N,WA2)
152 IF (FP .LE. P1*DELTA) GO TO 220
154 C IF THE JACOBIAN IS NOT RANK DEFICIENT, THE NEWTON
155 C STEP PROVIDES A LOWER BOUND, PARL, FOR THE ZERO OF
156 C THE FUNCTION. OTHERWISE SET THIS BOUND TO ZERO.
159 IF (NSING .LT. N) GO TO 120
162 WA1(J) = DIAG(L)*(WA2(L)/DXNORM)
167 IF (JM1 .LT. 1) GO TO 100
169 SUM = SUM + R(I,J)*WA1(I)
172 WA1(J) = (WA1(J) - SUM)/R(J,J)
175 PARL = ((FP/DELTA)/TEMP)/TEMP
178 C CALCULATE AN UPPER BOUND, PARU, FOR THE ZERO OF THE FUNCTION.
183 SUM = SUM + R(I,J)*QTB(I)
190 IF (PARU .EQ. ZERO) PARU = DWARF/MIN(DELTA,P1)
192 C IF THE INPUT PAR LIES OUTSIDE OF THE INTERVAL (PARL,PARU),
193 C SET PAR TO THE CLOSER ENDPOINT.
197 IF (PAR .EQ. ZERO) PAR = GNORM/DXNORM
199 C BEGINNING OF AN ITERATION.
204 C EVALUATE THE FUNCTION AT THE CURRENT VALUE OF PAR.
206 IF (PAR .EQ. ZERO) PAR = MAX(DWARF,P001*PARU)
209 WA1(J) = TEMP*DIAG(J)
211 CALL QRSOLV(N,R,LDR,IPVT,WA1,QTB,X,SIGMA,WA2)
213 WA2(J) = DIAG(J)*X(J)
215 DXNORM = ENORM(N,WA2)
219 C IF THE FUNCTION IS SMALL ENOUGH, ACCEPT THE CURRENT VALUE
220 C OF PAR. ALSO TEST FOR THE EXCEPTIONAL CASES WHERE PARL
221 C IS ZERO OR THE NUMBER OF ITERATIONS HAS REACHED 10.
223 IF (ABS(FP) .LE. P1*DELTA
224 1 .OR. PARL .EQ. ZERO .AND. FP .LE. TEMP
225 2 .AND. TEMP .LT. ZERO .OR. ITER .EQ. 10) GO TO 220
227 C COMPUTE THE NEWTON CORRECTION.
231 WA1(J) = DIAG(L)*(WA2(L)/DXNORM)
234 WA1(J) = WA1(J)/SIGMA(J)
237 IF (N .LT. JP1) GO TO 200
239 WA1(I) = WA1(I) - R(I,J)*TEMP
244 PARC = ((FP/DELTA)/TEMP)/TEMP
246 C DEPENDING ON THE SIGN OF THE FUNCTION, UPDATE PARL OR PARU.
248 IF (FP .GT. ZERO) PARL = MAX(PARL,PAR)
249 IF (FP .LT. ZERO) PARU = MIN(PARU,PAR)
251 C COMPUTE AN IMPROVED ESTIMATE FOR PAR.
253 PAR = MAX(PARL,PAR+PARC)
255 C END OF AN ITERATION.
262 IF (ITER .EQ. 0) PAR = ZERO
265 C LAST CARD OF SUBROUTINE LMPAR.