website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_f/slatec/lmpar.f

   1 *DECK LMPAR
   2       SUBROUTINE LMPAR (N, R, LDR, IPVT, DIAG, QTB, DELTA, PAR, X,
   3      +   SIGMA, WA1, WA2)
   4 C***BEGIN PROLOGUE  LMPAR
   5 C***SUBSIDIARY
   6 C***PURPOSE  Subsidiary to SNLS1 and SNLS1E
   7 C***LIBRARY   SLATEC
   8 C***TYPE      SINGLE PRECISION (LMPAR-S, DMPAR-D)
   9 C***AUTHOR  (UNKNOWN)
  10 C***DESCRIPTION
  11 C
  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
  16 C
  17 C           A*X = B ,     SQRT(PAR)*D*X = 0 ,
  18 C
  19 C     in the least squares sense, and DXNORM is the Euclidean
  20 C     norm of D*X, then either PAR is zero and
  21 C
  22 C           (DXNORM-DELTA) .LE. 0.1*DELTA ,
  23 C
  24 C     or PAR is positive and
  25 C
  26 C           ABS(DXNORM-DELTA) .LE. 0.1*DELTA .
  27 C
  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
  37 C
  38 C            T   T                   T
  39 C           P *(A *A + PAR*D*D)*P = S *S .
  40 C
  41 C     S is employed within LMPAR and may be of separate interest.
  42 C
  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.
  47 C
  48 C     The subroutine statement is
  49 C
  50 C       SUBROUTINE LMPAR(N,R,LDR,IPVT,DIAG,QTB,DELTA,PAR,X,SIGMA,
  51 C                        WA1,WA2)
  52 C
  53 C     where
  54 C
  55 C       N is a positive integer input variable set to the order of R.
  56 C
  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.
  62 C
  63 C       LDR is a positive integer input variable not less than N
  64 C         which specifies the leading dimension of the array R.
  65 C
  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.
  69 C
  70 C       DIAG is an input array of length N which must contain the
  71 C         diagonal elements of the matrix D.
  72 C
  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.
  75 C
  76 C       DELTA is a positive input variable which specifies an upper
  77 C         bound on the Euclidean norm of D*X.
  78 C
  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.
  82 C
  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,
  85 C         for the output PAR.
  86 C
  87 C       SIGMA is an output array of length N which contains the
  88 C         diagonal elements of the upper triangular matrix S.
  89 C
  90 C       WA1 and WA2 are work arrays of length N.
  91 C
  92 C***SEE ALSO  SNLS1, SNLS1E
  93 C***ROUTINES CALLED  ENORM, QRSOLV, R1MACH
  94 C***REVISION HISTORY  (YYMMDD)
  95 C   800301  DATE WRITTEN
  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.
 100 C           (WRB)
 101 C   900328  Added TYPE section.  (WRB)
 102 C***END PROLOGUE  LMPAR
 103       INTEGER N,LDR
 104       INTEGER IPVT(*)
 105       REAL DELTA,PAR
 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
 109       REAL R1MACH,ENORM
 110       SAVE P1, P001, ZERO
 111       DATA P1,P001,ZERO /1.0E-1,1.0E-3,0.0E0/
 112 C***FIRST EXECUTABLE STATEMENT  LMPAR
 113       DWARF = R1MACH(1)
 114 C
 115 C     COMPUTE AND STORE IN X THE GAUSS-NEWTON DIRECTION. IF THE
 116 C     JACOBIAN IS RANK-DEFICIENT, OBTAIN A LEAST SQUARES SOLUTION.
 117 C
 118       NSING = N
 119       DO 10 J = 1, N
 120          WA1(J) = QTB(J)
 121          IF (R(J,J) .EQ. ZERO .AND. NSING .EQ. N) NSING = J - 1
 122          IF (NSING .LT. N) WA1(J) = ZERO
 123    10    CONTINUE
 124       IF (NSING .LT. 1) GO TO 50
 125       DO 40 K = 1, NSING
 126          J = NSING - K + 1
 127          WA1(J) = WA1(J)/R(J,J)
 128          TEMP = WA1(J)
 129          JM1 = J - 1
 130          IF (JM1 .LT. 1) GO TO 30
 131          DO 20 I = 1, JM1
 132             WA1(I) = WA1(I) - R(I,J)*TEMP
 133    20       CONTINUE
 134    30    CONTINUE
 135    40    CONTINUE
 136    50 CONTINUE
 137       DO 60 J = 1, N
 138          L = IPVT(J)
 139          X(L) = WA1(J)
 140    60    CONTINUE
 141 C
 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.
 145 C
 146       ITER = 0
 147       DO 70 J = 1, N
 148          WA2(J) = DIAG(J)*X(J)
 149    70    CONTINUE
 150       DXNORM = ENORM(N,WA2)
 151       FP = DXNORM - DELTA
 152       IF (FP .LE. P1*DELTA) GO TO 220
 153 C
 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.
 157 C
 158       PARL = ZERO
 159       IF (NSING .LT. N) GO TO 120
 160       DO 80 J = 1, N
 161          L = IPVT(J)
 162          WA1(J) = DIAG(L)*(WA2(L)/DXNORM)
 163    80    CONTINUE
 164       DO 110 J = 1, N
 165          SUM = ZERO
 166          JM1 = J - 1
 167          IF (JM1 .LT. 1) GO TO 100
 168          DO 90 I = 1, JM1
 169             SUM = SUM + R(I,J)*WA1(I)
 170    90       CONTINUE
 171   100    CONTINUE
 172          WA1(J) = (WA1(J) - SUM)/R(J,J)
 173   110    CONTINUE
 174       TEMP = ENORM(N,WA1)
 175       PARL = ((FP/DELTA)/TEMP)/TEMP
 176   120 CONTINUE
 177 C
 178 C     CALCULATE AN UPPER BOUND, PARU, FOR THE ZERO OF THE FUNCTION.
 179 C
 180       DO 140 J = 1, N
 181          SUM = ZERO
 182          DO 130 I = 1, J
 183             SUM = SUM + R(I,J)*QTB(I)
 184   130       CONTINUE
 185          L = IPVT(J)
 186          WA1(J) = SUM/DIAG(L)
 187   140    CONTINUE
 188       GNORM = ENORM(N,WA1)
 189       PARU = GNORM/DELTA
 190       IF (PARU .EQ. ZERO) PARU = DWARF/MIN(DELTA,P1)
 191 C
 192 C     IF THE INPUT PAR LIES OUTSIDE OF THE INTERVAL (PARL,PARU),
 193 C     SET PAR TO THE CLOSER ENDPOINT.
 194 C
 195       PAR = MAX(PAR,PARL)
 196       PAR = MIN(PAR,PARU)
 197       IF (PAR .EQ. ZERO) PAR = GNORM/DXNORM
 198 C
 199 C     BEGINNING OF AN ITERATION.
 200 C
 201   150 CONTINUE
 202          ITER = ITER + 1
 203 C
 204 C        EVALUATE THE FUNCTION AT THE CURRENT VALUE OF PAR.
 205 C
 206          IF (PAR .EQ. ZERO) PAR = MAX(DWARF,P001*PARU)
 207          TEMP = SQRT(PAR)
 208          DO 160 J = 1, N
 209             WA1(J) = TEMP*DIAG(J)
 210   160       CONTINUE
 211          CALL QRSOLV(N,R,LDR,IPVT,WA1,QTB,X,SIGMA,WA2)
 212          DO 170 J = 1, N
 213             WA2(J) = DIAG(J)*X(J)
 214   170       CONTINUE
 215          DXNORM = ENORM(N,WA2)
 216          TEMP = FP
 217          FP = DXNORM - DELTA
 218 C
 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.
 222 C
 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
 226 C
 227 C        COMPUTE THE NEWTON CORRECTION.
 228 C
 229          DO 180 J = 1, N
 230             L = IPVT(J)
 231             WA1(J) = DIAG(L)*(WA2(L)/DXNORM)
 232   180       CONTINUE
 233          DO 210 J = 1, N
 234             WA1(J) = WA1(J)/SIGMA(J)
 235             TEMP = WA1(J)
 236             JP1 = J + 1
 237             IF (N .LT. JP1) GO TO 200
 238             DO 190 I = JP1, N
 239                WA1(I) = WA1(I) - R(I,J)*TEMP
 240   190          CONTINUE
 241   200       CONTINUE
 242   210       CONTINUE
 243          TEMP = ENORM(N,WA1)
 244          PARC = ((FP/DELTA)/TEMP)/TEMP
 245 C
 246 C        DEPENDING ON THE SIGN OF THE FUNCTION, UPDATE PARL OR PARU.
 247 C
 248          IF (FP .GT. ZERO) PARL = MAX(PARL,PAR)
 249          IF (FP .LT. ZERO) PARU = MIN(PARU,PAR)
 250 C
 251 C        COMPUTE AN IMPROVED ESTIMATE FOR PAR.
 252 C
 253          PAR = MAX(PARL,PAR+PARC)
 254 C
 255 C        END OF AN ITERATION.
 256 C
 257          GO TO 150
 258   220 CONTINUE
 259 C
 260 C     TERMINATION.
 261 C
 262       IF (ITER .EQ. 0) PAR = ZERO
 263       RETURN
 264 C
 265 C     LAST CARD OF SUBROUTINE LMPAR.
 266 C
 267       END