--- /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