+++ /dev/null
-*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
git://source.jalview.org / jabaws.git / blobdiff