1 subroutine SNQSOL(SNQFJ, N, X, FVEC, XTOL, IOPT, W, IDIMW) 2c Copyright (c) 1996 California Institute of Technology, Pasadena, CA. 3c ALL RIGHTS RESERVED. 4c Based on Government Sponsored Research NAS7-03001. 5c>> 2001-05-25 SNQSOL Krogh Minor change for making .f90 version. 6c>> 2000-12-01 SNQSOL Krogh Removed unused parameter P001. 7c>> 1996-05-16 SNQSOL Krogh Changes to use .C. and C%%. 8c>> 1996-03-30 SNQSOL Krogh Added external stmts. SIN => VSIN, etc. 9c>> 1994-11-02 SNQSOL Krogh Changes to use M77CON 10c>> 1992-04-27 SNQSOL CLL Deleted unreferenced stmt label. 11c>> 1992-04-07 CAO Extra comma in Print removed (error from VAX compile) 12c>> 1992-01-15 CLL 13c>> 1991-12-18 CLL & FTK Adding treatment of slow convergence to 0. 14c>> 1991-12-05 CLL & FTK Adding Option vector interface. 15c>> 1990-04-20 CLL@JPL Adapting code from Minpack for MATH77 16c 17c Solves a system of N nonlinear equations in N unknowns. 18c SNQSOL is the the user-interface subroutine. It calls SNQSL1 which 19c contains the top-level logic of the solution algorithm. 20c SNQSOL & SNQSL1 also need: 21c Other subroutines that are in this file: 22c SNQFDJ, SNQDOG, SNQQFM, SNQQRF, SNQUPD. 23c Other subprograms from the MATH77 library: SNRM2, SERV1, 24c [D/R]1MACH (Fortan 77 only), IERV1, & IERM1. 25C A user-provided subroutine: SNQFJ. 26c 27c Most of these subprograms are derived from MINPACK-1. 28c MINPACK-1, 1980, was developed by Jorge J. More', 29c Burton S. Garbow, and Kenneth E. Hillstrom, Argonne Nat'l Lab. 30c The MINPACK-1 code was obtained as FILE05 from MINPACK/EX from 31c Netlib, downloaded to JPL on Tue Feb 6 12:17:45 EST 1990. 32c 33c Old Name New Name 34c -------- -------- 35c HYBRJ1, HYBRD1 SNQSOL (Completely redesigned.) 36c HYBRJ, HYBRD SNQSL1 (Algorithm and code changes.) 37c DOGLEG SNQDOG 38c ENORM SNRM2 in BLAS and MATH77 39c FDJAC1 SNQFDJ 40c QFORM SNQQFM 41c QRFAC SNQQRF 42c R1MPYQ SNQAQ 43c R1UPDT SNQUPD 44c [D/S]PMPAR [D/R]1MACH in file amach.f (Fortran 77 only) 45c FCN SNQFJ 46c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 47c Arguments for SNQSOL 48c 49c call SNQSOL(SNQFJ, N, X, FVEC, XTOL, IOPT, W, IDIMW) 50c 51c SNQFJ Name of user-supplied subroutine. 52c 53c N [in] Problem size 54c X(N) [inout] Initial and final x-vector. 55c FVEC(N) [out] Final F values. 56c XTOL [in] Rel. Conv. tolerance on weighted X 57c IOPT() [inout] First 3 elements contain output values. 58c IOPT(1) = INFO. Output status. 59c IOPT(2) = NFEV. No. of F evals used. 60c IOPT(3) = NJEV. No. of evals of Jacobian. 61c 62c Ramaining elements in IOPT() select options. 63c 64c Option No. of Affected variables Affected variables 65c Number arguments in SNQSOL. in SNQSL1. 66c 1 0 HAVEJ HAVEJ 67c 2 1 DMODE, HAVED, W(4:3+N) HAVED, DIAG(1:N) 68c 3 1 NPRINT NPRINT 69c 4 1 MAXFEV MAXFEV 70c 5 2 ML, MU ML, MU 71c 6 0 W(1) EPSFCN 72c 7 0 W(2) FACTOR 73c 8 0 TRACE TRACE 74c 75c Functionality of options, listed by option numbers in square brackets. 76c [1] If set, user is not computing a Jacobian. 77c This subr sets HAVEJ = .false. 78c [2] Arg: DMODE = 1, 2, or 3. 79c 1. This subr sets DIAG() to all ones and HAVED = .true. 80c 2. User has set DIAG(). This subr sets HAVED = .true. 81c 3. This subr sets HAVED = .false. so SNQSL1 will set 82c DIAG() dynamically. 83c [3] Arg: NPRINT Print control. 84c [4] Arg: MAXFEV Limit on no. of F evals. 85c [5] Args: ML & MU Band structure. 86c [6] If set means EPSFCN has been set in W(1). 87c [7] If set means FACTOR has been set in W(2). 88c [8] If set, this subr sets TRACE = .true., else this subr 89c sets TRACE = .false. When TRACE is .true., SNQSL1 prints 90c detailed intermediate results. 91c 92c W() [inout] W(1) and W(2) may be used to pass EPSFCN and FACTOR 93c to the subroutine. W(3) contains TOLTST on return. 94c W( 4 : 3+(15*N + 3*N**2)/2 ) is used as work space. 95c 96c EPSFCN W(1) Error in F evals. Used in computing 97c approx derivs. 98c FACTOR W(2) Algorithm parameter. 99c TOLTST W(3) Output. Final value of quantity compared 100c with XTOL for convergence test. 101c DIAG(N) W(4:N+3) Scaling values. May be input or 102c computed. See option 2. 103c WA1(N) W() Work space of length N. 104c WA2(N) W() Work space of length N. 105c WA3(N) W() Work space of length N. 106c WA4(N) W() Work space of length N. 107c GNSTEP(N) W() Work space of length N. 108c QTF(N) W() Wrk space. At end has (Q**t)*F. 109c FJAC(N,N) W() Work space for Jacobian. At end has Q of 110c QR factorization. 111c R( (N + N**2)/2 ) W() Wrk spc. At end has Packed R of 112c QR factorization. 113c IDIMW [in] Dimension of W(). Require IDIMW .ge. 3+(15*N+3*N**2)/2 114c ------------------------------------------------------------------ 115c--S replaces "?": ?NQSOL,?NQSL1,?ERV1,?NQFJ,?NQDOG,?NRM2,?NQFDJ 116c--& ?NQQFM,?NQQRF,?NQAQ,?NQUPD 117c Also uses IERM1, IERV1 118c ------------------------------------------------------------------ 119 external R1MACH, SNQFJ 120 integer N, IOPT(*), IDIMW 121 real X(N), FVEC(N), XTOL, W(IDIMW) 122c 123 integer IWTOLT, IWDIAG, IWA1, IWA2, IWA3, IWA4, IWGNST 124 integer IWQTF, IWFJAC, IWR 125 parameter(IWTOLT = 3, IWDIAG = 4 ) 126 real R1MACH, EPSFCN, EPSMCH, FAC1, FACTOR 127 integer DMODE, J, JABS, K, NI, NPRINT, MAXF1, MAXFEV, ML, MU 128 logical JPOS, HAVEJ, HAVED, TRACE 129 parameter(FAC1 = 0.75e0, MAXF1 = 200) 130 save EPSMCH 131 data EPSMCH / 0.0e0 / 132c ------------------------------------------------------------------ 133c 134 if(EPSMCH .eq. 0.0e0) EPSMCH = R1MACH(4) 135 NI = N 136 IOPT(1) = 1 137 if (NI .le. 0) then 138 call IERM1('SNQSOL',IOPT(1),0,'Require N > 0','N',NI,'.') 139 go to 900 140 endif 141 if (IDIMW .lt. 3 + (NI*(15+3*NI))/2) then 142 call IERM1('SNQSOL',IOPT(1),0,'Require IDIMW .ge. NEED', 143 * 'IDIMW',IDIMW,',') 144 call IERV1('NEED =', 3 + (NI*(15+3*NI))/2,'.') 145 go to 900 146 endif 147c Set default values. 148 HAVEJ = .true. 149 DMODE = 1 150 NPRINT = 0 151 MAXFEV = MAXF1 * (NI + 1) 152 ML = NI - 1 153 MU = ML 154 EPSFCN = EPSMCH 155 FACTOR = FAC1 156 TRACE = .false. 157c 158c Loop on K beginning with K = 4 and 159c terminating when an option code, J, is zero. 160 K = 4 161 20 continue 162 J = IOPT(K) 163 JABS = abs(J) 164 JPOS = J .gt. 0 165 go to (40, 31, 32, 33, 34, 35, 36, 37, 38), JABS+1 166c 167c ANSI Standard Fortran 77 drops thru to here if JABS is 168c larger than 7. This is an error condition. 169c 170 call IERM1('SNQSOL',IOPT(1),0,'IOPT(K) must be in [-7..7]', 171 * 'K',K,',') 172 call IERV1('IOPT(K)',J,'.') 173 go to 900 174c 175 31 HAVEJ = .not. JPOS 176 K = K+1 177 go to 20 178c Option 2. Argument = 1, 2, or 3. Default = 1. 179c 1. This subr sets DIAG() to all ones. 180c 2. User has set DIAG(). 181c 3. Subr SNQSL1 sets DIAG() dynamically. 182 183 32 if( JPOS .and. IOPT(K+1) .eq. 2) then 184 DMODE = 2 185 elseif( JPOS .and. IOPT(K+1) .eq. 3) then 186 DMODE = 3 187 elseif(.not. JPOS .or. IOPT(K+1) .eq. 1) then 188 DMODE = 1 189 else 190c Error. 191 call IERM1('SNQSOL',IOPT(1),0,'Bad argument for Option 2.', 192 * 'Argument',IOPT(K+1),'.') 193 go to 900 194 endif 195 K = K+2 196 go to 20 197 33 if(JPOS) then 198 NPRINT = IOPT(K+1) 199 else 200 NPRINT = 0 201 endif 202 K = K+2 203 go to 20 204 34 if(JPOS) then 205 MAXFEV = IOPT(K+1) 206 else 207 MAXFEV = MAXF1 * (NI + 1) 208 endif 209 K = K+2 210 go to 20 211 35 if(JPOS) then 212 ML = IOPT(K+1) 213 MU = IOPT(K+2) 214 else 215 ML = NI+1 216 MU = ML 217 endif 218 K = K+3 219 go to 20 220 36 if(JPOS) then 221 EPSFCN = W(1) 222 else 223 EPSFCN = EPSMCH 224 endif 225 K = K+1 226 go to 20 227 37 If(JPOS) then 228 FACTOR = W(2) 229 else 230 FACTOR = FAC1 231 endif 232 K = K+1 233 go to 20 234 38 If(JPOS) then 235 TRACE = .true. 236 else 237 TRACE = .false. 238 endif 239 K = K+1 240 go to 20 241c End loop on K 242 40 continue 243c 244c Option 2. DMODE = 1, 2, or 3. 245c 1. This subr sets DIAG() to all ones. 246c 2. User has set DIAG(). 247c 3. Subr SNQSL1 sets DIAG() dynamically. 248 249 if(DMODE .eq. 1) then 250 HAVED = .true. 251 do 50 K = IWDIAG, IWDIAG+NI-1 252 W(K) = 1.0e0 253 50 continue 254 else 255 HAVED = DMODE .eq. 2 256 endif 257c 258 IWA1 = IWDIAG + NI 259 IWA2 = IWA1 + NI 260 IWA3 = IWA2 + NI 261 IWA4 = IWA3 + NI 262 IWGNST = IWA4 + NI 263 IWQTF = IWGNST + NI 264 IWFJAC = IWQTF + NI 265 IWR = IWFJAC + NI*NI 266c IWNEXT = IWR + (N * (N+1)) / 2 Next available loc in W(). 267c 268 call SNQSL1(SNQFJ, NI, X, FVEC, XTOL, 269 1 IOPT(1), IOPT(2), IOPT(3), 270 2 NPRINT, HAVEJ, MAXFEV, HAVED, ML, MU, 271 3 EPSFCN, FACTOR, TRACE, W(IWTOLT), W(IWDIAG), 272 4 W(IWA1), W(IWA2), W(IWA3), W(IWA4), W(IWGNST), W(IWQTF), 273 5 W(IWFJAC), W(IWR)) 274 return 275c Error return 276 900 continue 277 IOPT(2) = 0 278 IOPT(3) = 0 279 W(3) = 0.0e0 280 return 281 end 282c ================================================================== 283 subroutine SNQSL1(SNQFJ, N, X, FVEC, XTOL, 284 * INFO, NFEV, NJEV, 285 * NPRINT, HAVEJ, MAXFEV, HAVED, ML, MU, 286 * EPSFCN, FACTOR, TRACE, TOLTST, DIAG, 287 * WA1, WA2, WA3, WA4, GNSTEP, QTF, FJAC, R) 288c>> 1991-12-04 CLL 289c>> 1991-12-02 CLL 290c>> 1991-06-18 CLL@JPL Adapting code from Minpack for MATH77 291 292c 26 arguments. 293c Dimension of R() must be (N + N**2)/2. 294c Total space occupied by EPSFCN, FACTOR, and TOLTST through R is 295c 3 + (15*N + 3*N**2)/2 296 297 external SNQFJ 298 integer N, MAXFEV, NPRINT, INFO, NFEV, NJEV, ML, MU 299 logical HAVEJ, HAVED, TRACE 300 real XTOL, EPSFCN, FACTOR, TOLTST 301 real X(N), FVEC(N), FJAC(N,N), DIAG(N), R(*) 302 real QTF(N), WA1(N), WA2(N), WA3(N), WA4(N) 303 real GNSTEP(N) 304C ********** 305C 306C SUBROUTINE SNQSL1 307C 308C THE PURPOSE OF SNQSL1 IS TO FIND A ZERO OF A SYSTEM OF 309C N NONLINEAR FUNCTIONS IN N VARIABLES BY A MODIFICATION 310C OF THE POWELL HYBRID METHOD. THE USER MUST PROVIDE A 311C SUBROUTINE WHICH CALCULATES THE FUNCTIONS and THE JACOBIAN. 312C 313C ------------------------------------------------------------------ 314c Arguments 315c 316c SNQFJ is THE NAME OF THE USER-SUPPLIED SUBROUTINE WHICH 317c CALCULATES THE FUNCTIONS and THE JACOBIAN. SNQFJ MUST 318c BE DECLARED IN AN EXTERNAL STATEMENT IN THE USER 319c CALLING PROGRAM. SNQFJ will not be called with IFLAG = 2 320c if HAVEJ is .false. SNQFJ will not be called with IFLAG = 0 321c if NPRINT is <= 0. 322c SNQFJ is specified as follows: 323C 324c subroutine SNQFJ(N, X, FVEC, FJAC, IFLAG) 325c integer N, IFLAG 326c real X(N), FVEC(N), FJAC(N,N) 327c ---------- 328c if IFLAG = 0, Print X() and FVEC() and return. 329c IF IFLAG = 1 CALCULATE THE FUNCTIONS AT X AND 330c RETURN THIS VECTOR IN FVEC. DO NOT ALTER FJAC. 331c IF IFLAG = 2 CALCULATE THE JACOBIAN AT X AND 332c RETURN THIS MATRIX IN FJAC. DO NOT ALTER FVEC. 333c Set IFLAG to a negative value to force an immediate 334c termination of the solution procedure. Otherwise do not 335c alter IFLAG. 336c --------- 337c RETURN 338c END 339C 340c N is A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER 341c OF FUNCTIONS and VARIABLES. 342C 343c X is AN ARRAY OF LENGTH N. ON INPUT X MUST CONTAIN 344c AN INITIAL ESTIMATE OF THE SOLUTION VECTOR. ON OUTPUT X 345c CONTAINS THE FINAL ESTIMATE OF THE SOLUTION VECTOR. 346C 347c FVEC is AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS 348c THE FUNCTIONS EVALUATED AT THE OUTPUT X. 349C 350c XTOL is A NONNEGATIVE INPUT VARIABLE. TERMINATION 351c OCCURS WHEN THE RELATIVE ERROR BETWEEN TWO CONSECUTIVE 352c ITERATES is AT MOST XTOL. 353C 354c INFO [integer,out] If the user has terminated execution by setting 355c IFLAG negative in SNQFJ, INFO is set to IFLAG. 356c Otherwise, INFO is set as follows: 357C 358c INFO = 0 Successful termination. Radius of trust region has 359c been reduced to at most max(XTOL, machine precision). 360C 361c INFO = 1 IMPROPER INPUT PARAMETERS. 362C 363c INFO = 2 Number of calls to SNQFJ for function evaluations has 364c reached MAXFEV. 365C 366c INFO = 3 XTOL is TOO SMALL. NO FURTHER IMPROVEMENT IN 367c THE APPROXIMATE SOLUTION X is POSSIBLE. 368C 369c INFO = 4 Iteration is not making good progress, as 370c measured by the improvement through the last 371c five Jacobian evaluations. 372C 373c INFO = 5 Iteration is not making good progress, as 374c measured by the improvement through last 375c ten function evaluations. 376C 377c NFEV [out,integer] The number of calls to SNQFJ with IFLAG = 1. 378C 379c NJEV [out,integer] The number of evaluations of the Jacobian matrix. 380c If HAVEJ is .true. this will be the number of calls to SNQFJ with 381c IFLAG = 2. Otherwise it is the number of times the Jacobian has 382c been approximately computed by differencing. 383C 384c NPRINT [in, integer] Enables controlled printing of iterates if it 385c is positive. In this case, SNQFJ is called with IFLAG = 0 at the 386c beginning of the first iteration and every NPRINTth time a new X 387c vector is accepted as an improvement, and at termination. 388c On these calls the new best X and FVEC are made available for 389c printing. FVEC and FJAC should not be altered. 390c If NPRINT is not positive, no special calls to SNQFJ with 391c IFLAG = 0 will be made. 392C 393c HAVEJ [in, logical] True means the user subroutine SNQFJ contains 394c code for computing the Jacobian matrix, and false means it does 395c not. 396c 397c MAXFEV is A POSITIVE INTEGER INPUT VARIABLE. TERMINATION 398c OCCURS WHEN THE NUMBER OF CALLS TO SNQFJ WITH IFLAG = 1 399c HAS REACHED MAXFEV. 400C 401c HAVED = true means initial values of DIAG() are given by the 402c calling program. False means this subroutine must compute 403c initial values of DIAG(). It will set DIAG(j) = the euclidean 404c norm of column j, unless this is zero, in which case it will 405c set DIAG(j) = 0.0. 406C 407c ML and MU specify the band structure, if any, of the Jacobian 408c matrix. All nonzero elements of the Jacobian matrix lie 409c within the first ML subdiagonals, the main diagonal, and the 410c first MU superdiagonals. 411c ML and MU are only used when HAVEJ is .false. and are only useful 412c if ML+MU+1 < N. In this case they are used to 413c reduce the number of function evaluations in estimating 414c derivatives. If the Jacobian has no band structure set 415c ML = MU = N-1. 416C 417c EPSFCN is an input variable used in determining a suitable 418c step length for the forward-difference approximation. This 419c approximation assumes that the relative errors in the 420c functions are of the order of max(EPSFCN, Machine precision). 421C 422c FACTOR is a positive input variable used in determining the 423c initial step bound. This bound is set to the product of 424c FACTOR and the euclidean norm of DIAG*X if nonzero, or else 425c to FACTOR itself. In most cases FACTOR should lie in the 426c interval (0.1, 10.0). Default: FACTOR = 0.75. 427C 428c TRACE [in, logical] If true, this subr will print detailed 429c intermediate output. Otherwise it will not. 430c 431c TOLTST [out] Final value of quantity that is compared with 432c XTOL for convergence test. 433c 434c DIAG is an array of length N. If HAVED = false, 435c DIAG is internally set. If HAVED = true, DIAG() 436c MUST CONTAIN POSITIVE ENTRIES THAT SERVE AS 437c MULTIPLICATIVE SCALE FACTORS FOR THE VARIABLES. 438C 439c WA1, WA2, WA3, and WA4 are work arrays of length N. 440c 441c GNSTEP() [scratch] Work array of length N to save the 442c Gauss-Newton step vector computed in SNQDOG. 443c 444c QTF is AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS 445c THE VECTOR (Q TRANSPOSE)*FVEC. 446C 447c FJAC is AN OUTPUT N BY N ARRAY WHICH CONTAINS THE 448c ORTHOGONAL MATRIX Q PRODUCED BY THE QR FACTORIZATION 449c OF THE FINAL APPROXIMATE JACOBIAN. 450C 451c R is AN OUTPUT ARRAY OF LENGTH LR WHICH CONTAINS THE packed 452c UPPER TRIANGULAR MATRIX PRODUCED BY THE QR FACTORIZATION 453c OF THE FINAL APPROXIMATE JACOBIAN, STORED ROWWISE. 454C -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 455C SUBPROGRAMS CALLED 456C 457C USER-SUPPLIED ...... SNQFJ 458C 459C MINPACK-SUPPLIED ... SNQDOG,R1MACH,SNRM2,SNQFDJ, 460C SNQQFM,SNQQRF,SNQAQ,SNQUPD 461C 462C FORTRAN-SUPPLIED ... abs,max,min,mod 463C 464C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. 465C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE' 466c Argonne Reports: ANL-80-68 and ANL-80-74, 1980. 467C 468c 1991-12-09 CLL at JPL. Replacing integer argument MODE that had 469c values 2 or 1 with logical argument HAVED related to MODE by 470c HAVED = MODE .eq. 2. Thus the user must set HAVED = .true. when 471c supplying the DIAG() values, and .false. otherwise. 472C ********** 473c ------------------------------------------------------------------ 474c Description of some of the local variables. 475c 476c DELTA [flpt] Diameter of trust region. 477c HLIM0 [flpt] Upper limit on DELTA when working with a computed 478c Jacobian. 479c HLIM1 [flpt] Upper limit on DELTA when working with an updated 480c Jacobian. 481c JACT [integer] Can have values of COMPUT, UPDATE, or KEEP. 482c Initially set to COMPUT. At the beginning of the main loop 483c we either compute a new Jacobian, update the Jacobian, or keep 484c the old Jacobian, depending on the setting of JACT. 485c JACT0 [integer] Saves the value of JACT at the beginning of the 486c main loop. As JACT gets changed in the loop, JACT0 is still 487c available as a record of what it was at the beginning of the loop. 488c JEVAL [logical] Set true whenever the Jacobian is computed, and set 489c false when it is updated. 490c NBEST [integer] Counter, incremented each time an x-vector is 491c accepted as being a better approximation to the solution. Used 492c in connection with NPRINT to trigger calles to SNQFJ for printing. 493c NCFAIL [integer] Counts consecutive "failed" steps since the last 494c computation of the Jacobian. NCFAIL is set to 0 when the Jacobian 495c is computed or when a step "succeeds" in the sense that 496c RATIO .ge. 0.1. It is incremented when RATIO .lt. 0.1. 497c NLOOP [integer] Counter for main iteration loop. 498c NUPDAT [integer] Counts number of consecutive times the Jacobian 499c matrix is updated. 500c TRYZER [logical] Initially set to true. While true, the algorithm 501c will monitor X's to see if they seem to be all approaching zero. 502c If so will try setting them all to zero. If this gives an exactly 503c zero function vector then we are finished. If not, we set TRYZER 504c to false and restore X to its previous value (even if the function 505c value at X = 0 was an improvement) and omit any further testing 506c for X's approaching zero. (We tryed accepting the X reached by 507c this exceptional step if the function value was an improvement, 508c but in one test case this caused the algorithm to end at a local 509c nonzero minimum rather than finding a zero.) 510c ------------------------------------------------------------------ 511 external R1MACH, SNRM2 512 integer COMPUT, I, IFLAG, IWA(1), J, JACT, JACT0 513 integer KEEP, L, LDFJAC, LR 514 integer MSUM, NBEST, NCFAIL, NCSUC, NEXTPR 515 integer NLOOP, NSLOW1, NSLOW2, NUMNWT, NUPDAT, UPDATE 516 logical JEVAL, NEWX, NEWTOK, SING, TRYZER 517 real R1MACH,SNRM2 518 real ACTRED,DELTA,EPSMCH,FNORM,FNORM1, HLIM0, HLIM1 519 real ONE,PNORM, PRERED,P1,P5,P0001,RATIO 520 real SUM,TEMP,XNORM, ZERO 521 parameter(COMPUT = 1, UPDATE = 2, KEEP = 3) 522 parameter(ONE = 1.0e0, P1 = 0.1e0, P5 = 0.5e0) 523 parameter(P0001 = 0.0001e0, ZERO = 0.0e0) 524 save EPSMCH 525 data EPSMCH /0.0e0 / 526c ------------------------------------------------------------------ 527C Set EPSMCH to the machine precision. 528C 529 if(EPSMCH .eq. 0.0e0) EPSMCH = R1MACH(4) 530C 531C Initialize values of output arguments. 532C 533 INFO = 1 534 NFEV = 0 535 NJEV = 0 536 TOLTST = 0.0e0 537 TRYZER = .true. 538C 539C CHECK THE INPUT PARAMETERS FOR ERRORS. 540C We assume the condition N > 0 has already been checked in 541c the user-interface subroutine that called this one. 542 543 IF ( XTOL .lt. ZERO .or. MAXFEV .le. 0 544 * .or. FACTOR .le. ZERO ) then 545 call IERM1('SNQSL1',INFO,0, 546 * 'Require MAXFEV > 0, XTOL .gt. 0.0, FACTOR > 0.0', 547 * 'MAXFEV',MAXFEV,',') 548 call SERV1('XTOL',XTOL,',') 549 call SERV1('FACTOR',FACTOR,'.') 550 go to 300 551 endif 552 if( .not. HAVEJ .and. (ML .lt. 0 .or. MU .lt. 0)) then 553 call IERM1('SNQSL1',INFO,0, 554 * 'With HAVEJ false, require ML .ge. 0 and MU .ge. 0', 555 * 'ML',ML,',') 556 call IERV1('MU',MU,',') 557 go to 300 558 endif 559c HAVED = true means the user has set DIAG(). 560 IF ( HAVED ) then 561 DO 10 J = 1, N 562 IF (DIAG(J) .le. ZERO) then 563 call IERM1('SNQSL1',INFO,0, 564 * 'With HAVED = .true., require all DIAG(J) > 0.0', 565 * 'J',J,',') 566 call SERV1('DIAG(J)',DIAG(J),'.') 567 go to 300 568 endif 569 10 CONTINUE 570 endif 571c Initialize algorithm variables. 572 INFO = 0 573 JACT = COMPUT 574 LDFJAC = N 575 LR = (N*(N+1)) / 2 576 MSUM = min(ML + MU + 1, N) 577 NBEST = 1 578 NCSUC = 0 579 NEXTPR = 1 580 NLOOP = 0 581 NSLOW1 = 0 582 NSLOW2 = 0 583 NUMNWT = 0 584C 585C Evaluate the function at the starting point. 586C Calculate and test its norm. 587C 588 IFLAG = 1 589C%% (*snqfj)( n, x, fvec, fjac, &iflag ); 590 CALL SNQFJ(N, X, FVEC, FJAC, IFLAG) 591 NFEV = 1 592 IF (IFLAG .lt. 0) GO TO 300 593 FNORM = SNRM2(N,FVEC,1) 594 if(TRACE) then 595 print'(1x,i5,a/(6x,5g15.6))',NLOOP, 596 * ' Initial X:',(X(J),J=1,N) 597 print'(1x,5x,a,g15.6)', 598 * ' Initial FNORM:',FNORM 599 endif 600 if(FNORM .eq. 0.0e0) then 601 go to 300 602 endif 603C 604C Beginning of main loop. 605C 606 30 continue 607 NLOOP = NLOOP + 1 608 JACT0 = JACT 609C 610C Compute, Update, or Keep Jacobian, depending on JACT. 611C 612 if (JACT .eq. COMPUT) then 613 JEVAL = .TRUE. 614 NUPDAT = 0 615 NCFAIL = 0 616C 617C CALCULATE THE JACOBIAN MATRIX. 618C 619 if(TRACE) print'(1x,i5,a)',NLOOP, 620 * ' Computing new Jacobian matrix.' 621 NJEV = NJEV + 1 622 if(HAVEJ) then 623 IFLAG = 2 624C%% (*snqfj)( n, x, fvec, fjac, &iflag ); 625 CALL SNQFJ(N, X, FVEC, FJAC, IFLAG) 626 else 627 CALL SNQFDJ(SNQFJ,N,X,FVEC,FJAC,LDFJAC, 628 * IFLAG,ML,MU,EPSFCN,WA1, WA2) 629 NFEV = NFEV + MSUM 630 endif 631 IF (IFLAG .lt. 0) GO TO 300 632C 633C COMPUTE THE QR FACTORIZATION OF THE JACOBIAN. 634C 635 CALL SNQQRF(N,N,FJAC,LDFJAC, .false., IWA,1,WA1,WA2,WA3) 636C 637C On the first iteration and if HAVED is .false., scale according 638C to the norms of the columns of the initial Jacobian. 639C Also on the first iteration calculate the norm of the scaled X 640C and initialize the trust region diameter, DELTA. 641C 642 IF (NLOOP .eq. 1) then 643 IF ( .not. HAVED ) then 644 DO 40 J = 1, N 645 DIAG(J) = WA2(J) 646 IF (WA2(J) .eq. ZERO) DIAG(J) = ONE 647 40 CONTINUE 648 endif 649C 650 DO 60 J = 1, N 651 WA3(J) = DIAG(J)*X(J) 652 60 CONTINUE 653 XNORM = SNRM2(N,WA3,1) 654 DELTA = FACTOR*XNORM 655 IF (DELTA .eq. ZERO) DELTA = FACTOR 656 endif 657C 658C FORM (Q TRANSPOSE)*FVEC and STORE IN QTF. 659C 660 DO 80 I = 1, N 661 QTF(I) = FVEC(I) 662 80 CONTINUE 663 DO 120 J = 1, N 664 IF (FJAC(J,J) .eq. ZERO) GO TO 110 665 SUM = ZERO 666 DO 90 I = J, N 667 SUM = SUM + FJAC(I,J)*QTF(I) 668 90 CONTINUE 669 TEMP = -SUM/FJAC(J,J) 670 DO 100 I = J, N 671 QTF(I) = QTF(I) + FJAC(I,J)*TEMP 672 100 CONTINUE 673 110 CONTINUE 674 120 CONTINUE 675C 676C COPY THE TRIANGULAR FACTOR OF THE QR FACTORIZATION INTO R. 677c The diagonal elts come from WA1(). The strictly upper 678c triangular elts come from FJAC(,). The upper triangular matrix 679c will be stored, packed by rows, in R(). 680C 681 SING = .FALSE. 682 DO 150 J = 1, N 683 L = J 684 DO 130 I = 1, J-1 685 R(L) = FJAC(I,J) 686 L = L + N - I 687 130 CONTINUE 688 R(L) = WA1(J) 689 IF (WA1(J) .eq. ZERO) SING = .true. 690 150 CONTINUE 691C 692C ACCUMULATE THE ORTHOGONAL FACTOR IN FJAC. 693C 694 CALL SNQQFM(N,N,FJAC,LDFJAC,WA1) 695C 696C RESCALE IF NECESSARY. 697C 698 if ( .not. HAVED ) then 699 DO 160 J = 1, N 700 DIAG(J) = max(DIAG(J),WA2(J)) 701 160 CONTINUE 702 endif 703c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 704 elseif(JACT .eq. UPDATE) then 705 706C 707C CALCULATE THE RANK ONE MODIFICATION TO THE JACOBIAN 708C and UPDATE QTF IF NECESSARY. 709C 710 if(TRACE) print'(1x,i5,a)',NLOOP, 711 * ' Updating Jacobian matrix.' 712 NUPDAT = NUPDAT + 1 713 JEVAL = .FALSE. 714 DO 280 J = 1, N 715 SUM = ZERO 716 DO 270 I = 1, N 717 SUM = SUM + FJAC(I,J)*WA4(I) 718 270 CONTINUE 719 WA2(J) = (SUM - WA3(J))/PNORM 720 WA1(J) = DIAG(J)*((DIAG(J)*WA1(J))/PNORM) 721 IF (RATIO .ge. P0001) QTF(J) = SUM 722 280 CONTINUE 723C 724C COMPUTE THE QR FACTORIZATION OF THE UPDATED JACOBIAN. 725C 726 CALL SNQUPD(N,N,R,LR,WA1,WA2,WA3,SING) 727 CALL SNQAQ(N,N,FJAC,LDFJAC,WA2,WA3) 728 CALL SNQAQ(1,N,QTF,1,WA2,WA3) 729 else 730 if(TRACE) print'(1x,i5,a)',NLOOP, 731 * ' Keeping Jacobian matrix unchanged.' 732 endif 733C 734C Now have a new or updated or retained Jacobian matrix. 735C 736C IF REQUESTED, CALL SNQFJ TO ENABLE PRINTING OF ITERATES. 737C 738 if (NPRINT .gt. 0) then 739 if (NBEST .eq. NEXTPR) then 740 IFLAG = 0 741C%% (*snqfj)( n, x, fvec, fjac, &iflag ); 742 CALL SNQFJ(N, X, FVEC, FJAC, IFLAG) 743 IF (IFLAG .lt. 0) GO TO 300 744 NEXTPR = NEXTPR + NPRINT 745 endif 746 endif 747C 748C Determine the direction P, using a dogleg method, and 749c returning -P in WA1(). 750C 751 CALL SNQDOG(N,R,LR,DIAG,QTF,DELTA,WA1,NEWTOK,WA2,WA3, 752 * JACT0 .eq. KEEP, GNSTEP) 753c 754c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 755 if(TRYZER) then 756c NUMNWT counts number of consecutive 757c full Newton steps. 758 if(NEWTOK) then 759 NUMNWT = NUMNWT + 1 760 else 761 NUMNWT = 0 762 endif 763c 764c Test for convergence of some x components toward 0. 765c If this seems to be happening try setting such 766c components to 0. 767c 768 if(NUMNWT .ge. 5 .and. NCSUC .ge. 4) then 769 NUMNWT = 0 770 do 204 J = 1,N 771 WA2(J) = X(J) - WA1(J) 772 if(abs(WA2(J)) .le. 0.75e0 * abs(X(J)) ) then 773 WA2(J) = 0.0e0 774 else 775 go to 203 776 endif 777 204 continue 778 if(TRACE) print'(1x,i5,a)',NLOOP, 779 * ' Trial setting of X() to zero.' 780C 781C EVALUATE THE FUNCTION AT WA2() and CALCULATE ITS NORM. 782C 783 IFLAG = 1 784c%% (*snqfj)( n, wa2, wa4, fjac, &iflag ); 785 CALL SNQFJ(N, WA2, WA4, FJAC, IFLAG) 786 NFEV = NFEV + 1 787 IF (IFLAG .lt. 0) GO TO 300 788 FNORM1 = SNRM2(N,WA4,1) 789 if(TRACE) print'(1x,i5,a,g15.6)',NLOOP, 790 * ' FNORM1 = ', FNORM1 791 if(FNORM1 .eq. 0.0e0) then 792c 793C Accept new point as final solution. 794c Update X() and FVEC() and go to termination. 795C 796 INFO = 0 797 TOLTST = 0.0e0 798 do 201 J = 1, N 799 X(J) = WA2(J) 800 FVEC(J) = WA4(J) 801 201 continue 802 if(TRACE) print'(1x,i5,a,(6x,5g15.6))',NLOOP, 803 * ' Accepting X = all zeros.' 804 go to 300 805 else 806 TRYZER = .false. 807 endif 808 endif 809c The following "endif" matches "if(TRYZER)then" 810 endif 811c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 812C 813C STORE THE DIRECTION P and X + P. CALCULATE THE NORM OF P. 814C 815 203 continue 816 DO 200 J = 1, N 817 WA1(J) = -WA1(J) 818 WA2(J) = X(J) + WA1(J) 819 WA3(J) = DIAG(J)*WA1(J) 820 200 continue 821 822 PNORM = SNRM2(N,WA3,1) 823 if(TRACE) then 824 print'(1x,i5,a,/1x,5x,2g15.6)',NLOOP, 825 * ' DELTA PNORM',DELTA,PNORM 826 print'(6x,a/(6x,5g15.6))',' Trial X:',(WA2(J),J=1,N) 827 endif 828C 829C ON THE FIRST ITERATION, ADJUST THE INITIAL STEP BOUND. 830C 831 IF (NLOOP .eq. 1) then 832 DELTA = min(DELTA,PNORM) 833 HLIM0 = DELTA 834 HLIM1 = DELTA 835 endif 836C 837C EVALUATE THE FUNCTION AT X + P and CALCULATE ITS NORM. 838C 839 IFLAG = 1 840c%% (*snqfj)( n, wa2, wa4, fjac, &iflag ); 841 CALL SNQFJ(N, WA2, WA4, FJAC, IFLAG) 842 NFEV = NFEV + 1 843 IF (IFLAG .lt. 0) GO TO 300 844 FNORM1 = SNRM2(N,WA4,1) 845C 846C COMPUTE THE SCALED ACTUAL REDUCTION. 847C 848 ACTRED = -ONE 849 IF (FNORM1 .lt. FNORM) ACTRED = ONE - (FNORM1/FNORM)**2 850C 851C COMPUTE THE SCALED PREDICTED REDUCTION. 852C 853 L = 1 854 DO 220 I = 1, N 855 SUM = ZERO 856 DO 210 J = I, N 857 SUM = SUM + R(L)*WA1(J) 858 L = L + 1 859 210 CONTINUE 860 WA3(I) = QTF(I) + SUM 861 220 CONTINUE 862 TEMP = SNRM2(N,WA3,1) 863 PRERED = ZERO 864 IF (TEMP .lt. FNORM) PRERED = ONE - (TEMP/FNORM)**2 865C 866C COMPUTE THE RATIO OF THE ACTUAL TO THE PREDICTED 867C REDUCTION. 868C 869 RATIO = ZERO 870 IF (PRERED .gt. ZERO) RATIO = ACTRED/PRERED 871 if(TRACE) print'(1x,i5,a,/1x,5x,4g15.6)',NLOOP, 872 * ' FNORM1 ACTRED PRERED RATIO', 873 * FNORM1,ACTRED,PRERED,RATIO 874C 875c Analyze RATIO, NCSUC and JEVAL to decide on accepting or 876c rejecting the new X, and assigning new values to 877c NCSUC, JACT, and DELTA. 878c 879 if( RATIO .lt. 0.0000E0) then 880 NCSUC = 0 881 NCFAIL = NCFAIL + 1 882 NEWX = .false. 883 if(JEVAL) HLIM0 = min(HLIM0, 0.707107e0 * PNORM) 884 HLIM1 = min(HLIM1, 0.707107e0 * PNORM) 885 if( JEVAL .or. (NCFAIL .le. 1 .and. NUPDAT .le. 2)) then 886 JACT = KEEP 887 DELTA = 0.5e0 * PNORM 888 else 889 JACT = COMPUT 890 DELTA = HLIM0 891 endif 892 elseif( RATIO .lt. 0.1E0) then 893 NCSUC = 0 894 NCFAIL = NCFAIL + 1 895 NEWX = .true. 896 if(NCFAIL .le. 1 .and. NUPDAT .le. 2) then 897 JACT = UPDATE 898 DELTA = 0.5e0 * PNORM 899 else 900 JACT = COMPUT 901 DELTA = HLIM0 902 endif 903 else 904c Here we have RATIO .ge. 0.1 905 NCSUC = NCSUC + 1 906 NCFAIL = 0 907 NEWX = .true. 908 JACT = UPDATE 909 if(RATIO .lt. 0.5e0) then 910 if(NCSUC .ge. 5) 911 * HLIM1 = max(HLIM1, 1.414214e0 * PNORM) 912 if(NCSUC .ge. 2) 913 * DELTA = min(HLIM1, max(DELTA, 1.414214e0 * PNORM)) 914 elseif(RATIO .lt. 0.9e0) then 915 if(JACT0 .eq. COMPUT) 916 * HLIM0 = max(HLIM0, 1.414214e0 * PNORM) 917 if(NCSUC .ge. 4) 918 * HLIM1 = max(HLIM1, 1.414214e0 * PNORM) 919 if(NCSUC .ge. 2) 920 * DELTA = min(HLIM1, max(DELTA, 1.414214e0 * PNORM)) 921 elseif(RATIO .lt. 1.1e0) then 922 if(JACT0 .eq. COMPUT) 923 * HLIM0 = max(HLIM0, 2.0e0 * PNORM) 924 if(NCSUC .eq. 1) then 925 DELTA = 1.414214e0 * PNORM 926 else 927 DELTA = 2.0e0 * PNORM 928 endif 929 HLIM1 = max(HLIM1, DELTA) 930 endif 931 endif 932 HLIM0 = max(HLIM0, HLIM1) 933 if(TRACE) print'(1x,i5,a,a,/1x,5x,3i8,3g13.4)',NLOOP, 934 * ' NCSUC NCFAIL NUPDAT', 935 * ' DELTA HLIM0 HLIM1', 936 * NCSUC, NCFAIL, NUPDAT, DELTA,HLIM0, HLIM1 937C 938 if(NEWX) then 939c Accept new X, FVEC, and their norms. 940 DO 250 J = 1, N 941 X(J) = WA2(J) 942 WA2(J) = DIAG(J)*X(J) 943 FVEC(J) = WA4(J) 944 250 CONTINUE 945 XNORM = SNRM2(N,WA2,1) 946 FNORM = FNORM1 947 NBEST = NBEST + 1 948 if(TRACE) print'(1x,i5,a,g15.6)',NLOOP, 949 * ' Accepting new X with XNORM = ',XNORM 950 endif 951C 952C DETERMINE THE PROGRESS OF THE ITERATION. 953C 954 if( ACTRED .ge. 0.001e0) then 955 NSLOW1 = 0 956 else 957 NSLOW1 = NSLOW1 + 1 958 endif 959 if( ACTRED .ge. 0.1e0) then 960 NSLOW2 = 0 961 elseif( JACT0 .eq. COMPUT) then 962 NSLOW2 = NSLOW2 + 1 963 endif 964 if(TRACE) print'(1x,i5,a,/1x,5x,2(i11,4x))',NLOOP, 965 * ' NSLOW1 NSLOW2', 966 * NSLOW1, NSLOW2 967C 968C TEST FOR CONVERGENCE. 969C 970 IF (DELTA .le. XTOL*XNORM .or. FNORM .eq. ZERO) then 971 INFO = 0 972 if(TRACE) print'(1x,i5,a,/1x,5x,i14,g15.6)',NLOOP, 973 * ' INFO XNORM', INFO, XNORM 974 go to 295 975 endif 976C 977C TESTS FOR TERMINATION and STRINGENT TOLERANCES. 978C 979 IF (NFEV .ge. MAXFEV) INFO = 2 980 IF (P1*max(P1*DELTA,PNORM) .le. EPSMCH*XNORM) INFO = 3 981 IF (NSLOW2 .eq. 5) INFO = 4 982 IF (NSLOW1 .eq. 10) INFO = 5 983 IF (INFO .ne. 0) then 984 if(TRACE) print'(1x,i5,a,/1x,5x,i14,g15.6)',NLOOP, 985 * ' INFO XNORM', INFO, XNORM 986 call IERM1('SNQSL1',INFO, 0,'Unsuccessful termination.', 987 * 'INFO',INFO,'.') 988 go to 295 989 endif 990 go to 30 991C End of main loop. 992C 993c Come to following stmt when INFO has been set to 994c 2, 3, 4, or 5, or to 0 due to successful XTOL test. 995 295 continue 996 if(XNORM .ne. 0.0e0) then 997 TOLTST = DELTA / XNORM 998 else 999 TOLTST = DELTA 1000 endif 1001c 1002c Jump to following statement with IFLAG negative 1003c or INFO = 1 or INFO = 0 due to FNORM being zero. 1004c Here we have TOLTST = 0.0. 1005 300 continue 1006C 1007C TERMINATION, EITHER NORMAL OR USER IMPOSED. 1008C 1009 IF (IFLAG .lt. 0) INFO = IFLAG 1010 if(TRACE) print'(1x,i5,a,i3)',NLOOP, 1011 * ' Quitting with INFO = ',INFO 1012 IFLAG = 0 1013c%% if (nprint > 0) (*snqfj)( n, x, fvec, fjac, &iflag ); 1014 IF (NPRINT .gt. 0) CALL SNQFJ(N,X,FVEC,FJAC, IFLAG) 1015 if(INFO .lt. 0) then 1016 call IERM1('SNQSL1',INFO, 0, 1017 * 'Quitting because user code set IFLAG negative.', 1018 * 'IFLAG',INFO,'.') 1019 endif 1020 return 1021C 1022C Last line of subroutine SNQSL1. 1023C 1024 END 1025c ================================================================== 1026 subroutine SNQFDJ(SNQFJ,N,X,FVEC,FJAC,LDFJAC,IFLAG,ML,MU,EPSFCN, 1027 * WA1,WA2) 1028c>> 1991-12-04 CLL Changed arg list of user supplied subroutine. 1029c>> 1991-06-18 CLL@JPL Adapting code from Minpack for MATH77 1030 external SNQFJ 1031 integer N,LDFJAC,IFLAG,ML,MU 1032 real EPSFCN 1033 real X(N),FVEC(N),FJAC(LDFJAC,N),WA1(N),WA2(N) 1034C ********** 1035C 1036C SUBROUTINE SNQFDJ 1037C 1038C THIS SUBROUTINE COMPUTES A FORWARD-DIFFERENCE APPROXIMATION 1039C TO THE N BY N JACOBIAN MATRIX ASSOCIATED WITH A SPECIFIED 1040C PROBLEM OF N FUNCTIONS IN N VARIABLES. IF THE JACOBIAN HAS 1041C A BANDED FORM, THEN FUNCTION EVALUATIONS ARE SAVED BY ONLY 1042C APPROXIMATING THE NONZERO TERMS. 1043C 1044C THE SUBROUTINE STATEMENT IS 1045C 1046C subroutine SNQFDJ(SNQFJ,N,X,FVEC,FJAC,LDFJAC,IFLAG,ML,MU,EPSFCN, 1047C WA1,WA2) 1048C 1049C WHERE 1050C 1051C SNQFJ IS THE NAME OF THE USER-SUPPLIED SUBROUTINE WHICH 1052C CALCULATES THE FUNCTIONS. SNQFJ MUST BE DECLARED 1053C IN AN EXTERNAL STATEMENT IN THE USER CALLING 1054C PROGRAM, and SHOULD BE WRITTEN AS FOLLOWS. 1055C 1056C subroutine SNQFJ(N,X,FVEC,IFLAG) 1057C integer N,IFLAG 1058C real X(N),FVEC(N) 1059C ---------- 1060C CALCULATE THE FUNCTIONS AT X AND 1061C RETURN THIS VECTOR IN FVEC. 1062C ---------- 1063C RETURN 1064C END 1065C 1066C THE VALUE OF IFLAG SHOULD NOT BE CHANGED BY SNQFJ UNLESS 1067C THE USER WANTS TO TERMINATE EXECUTION OF SNQFDJ. 1068C IN THIS CASE SET IFLAG TO A NEGATIVE INTEGER. 1069C 1070C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER 1071C OF FUNCTIONS and VARIABLES. 1072C 1073C X IS AN INPUT ARRAY OF LENGTH N. 1074C 1075C FVEC IS AN INPUT ARRAY OF LENGTH N WHICH MUST CONTAIN THE 1076C FUNCTIONS EVALUATED AT X. 1077C 1078C FJAC IS AN OUTPUT N BY N ARRAY WHICH CONTAINS THE 1079C APPROXIMATION TO THE JACOBIAN MATRIX EVALUATED AT X. 1080C 1081C LDFJAC IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN N 1082C WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY FJAC. 1083C 1084C IFLAG IS AN INTEGER VARIABLE WHICH CAN BE USED TO TERMINATE 1085C THE EXECUTION OF SNQFDJ. SEE DESCRIPTION OF SNQFJ. 1086C 1087C ML IS A NONNEGATIVE INTEGER INPUT VARIABLE WHICH SPECIFIES 1088C THE NUMBER OF SUBDIAGONALS WITHIN THE BAND OF THE 1089C JACOBIAN MATRIX. IF THE JACOBIAN IS NOT BANDED, SET 1090C ML TO AT LEAST N - 1. 1091C 1092C MU IS A NONNEGATIVE INTEGER INPUT VARIABLE WHICH SPECIFIES 1093C THE NUMBER OF SUPERDIAGONALS WITHIN THE BAND OF THE 1094C JACOBIAN MATRIX. IF THE JACOBIAN IS NOT BANDED, SET 1095C MU TO AT LEAST N - 1. 1096C 1097C EPSFCN IS AN INPUT VARIABLE USED IN DETERMINING A SUITABLE 1098C STEP LENGTH FOR THE FORWARD-DIFFERENCE APPROXIMATION. THIS 1099C APPROXIMATION ASSUMES THAT THE RELATIVE ERRORS IN THE 1100C FUNCTIONS ARE OF THE ORDER OF EPSFCN. IF EPSFCN IS LESS 1101C THAN THE MACHINE PRECISION, IT IS ASSUMED THAT THE RELATIVE 1102C ERRORS IN THE FUNCTIONS ARE OF THE ORDER OF THE MACHINE 1103C PRECISION. 1104C 1105C WA1 and WA2 ARE WORK ARRAYS OF LENGTH N. IF ML + MU + 1 IS AT 1106C LEAST N, THEN THE JACOBIAN IS CONSIDERED DENSE, and WA2 IS 1107C NOT REFERENCED. 1108C 1109C SUBPROGRAMS CALLED 1110C 1111C MINPACK-SUPPLIED ... R1MACH 1112C 1113C FORTRAN-SUPPLIED ... abs,max,sqrt 1114C 1115C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. 1116C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE 1117C 1118C ********** 1119c ------------------------------------------------------------------ 1120 external R1MACH 1121 integer I,J,K,MSUM 1122 real EPS,EPSMCH,H,TEMP,ZERO 1123c++ CODE for ~.C. is active 1124 real DUMMY(1,1) 1125c++ CODE for .C. & (.N. == 'S') is inactive 1126c%% float *dummy; 1127c++ CODE for .C. & (.N. == 'D') is inactive 1128c%% double *dummy; 1129C++ End 1130 real R1MACH 1131 parameter(ZERO = 0.0e0) 1132C 1133C EPSMCH IS THE MACHINE PRECISION. 1134C 1135 EPSMCH = R1MACH(4) 1136C 1137 EPS = sqrt(max(EPSFCN,EPSMCH)) 1138 IFLAG = 1 1139 MSUM = ML + MU + 1 1140 IF (MSUM .lt. N) GO TO 40 1141C 1142C COMPUTATION OF DENSE APPROXIMATE JACOBIAN. 1143C 1144 DO 20 J = 1, N 1145 TEMP = X(J) 1146 H = EPS*abs(TEMP) 1147 IF (H .eq. ZERO) H = EPS 1148 X(J) = TEMP + H 1149c%% (*snqfj)( n, x, wa1, dummy, iflag ); 1150 CALL SNQFJ(N, X, WA1, DUMMY, IFLAG) 1151 IF (IFLAG .lt. 0) GO TO 30 1152 X(J) = TEMP 1153 DO 10 I = 1, N 1154 FJAC(I,J) = (WA1(I) - FVEC(I))/H 1155 10 CONTINUE 1156 20 CONTINUE 1157 30 CONTINUE 1158 GO TO 110 1159 40 CONTINUE 1160C 1161C COMPUTATION OF BANDED APPROXIMATE JACOBIAN. 1162C 1163 DO 90 K = 1, MSUM 1164 DO 60 J = K, N, MSUM 1165 WA2(J) = X(J) 1166 H = EPS*abs(WA2(J)) 1167 IF (H .eq. ZERO) H = EPS 1168 X(J) = WA2(J) + H 1169 60 CONTINUE 1170c%% (*snqfj)( n, x, wa1, dummy, iflag ); 1171 CALL SNQFJ(N, X, WA1, DUMMY, IFLAG) 1172 IF (IFLAG .lt. 0) GO TO 100 1173 DO 80 J = K, N, MSUM 1174 X(J) = WA2(J) 1175 H = EPS*abs(WA2(J)) 1176 IF (H .eq. ZERO) H = EPS 1177 DO 70 I = 1, N 1178 FJAC(I,J) = ZERO 1179 IF (I .ge. J - MU .and. I .le. J + ML) 1180 * FJAC(I,J) = (WA1(I) - FVEC(I))/H 1181 70 CONTINUE 1182 80 CONTINUE 1183 90 CONTINUE 1184 100 CONTINUE 1185 110 CONTINUE 1186 RETURN 1187C 1188C Last line of subroutine SNQFDJ. 1189C 1190 END 1191c ================================================================== 1192 subroutine SNQAQ(M,N,A,LDA,V,W) 1193 integer M,N,LDA 1194 real A(LDA,N),V(N),W(N) 1195C ********** 1196C 1197C SUBROUTINE SNQAQ 1198C 1199C GIVEN AN M BY N MATRIX A, THIS SUBROUTINE COMPUTES A*Q WHERE 1200C Q IS THE PRODUCT OF 2*(N - 1) TRANSFORMATIONS 1201C 1202C GV(N-1)*...*GV(1)*GW(1)*...*GW(N-1) 1203C 1204C and GV(I), GW(I) ARE GIVENS ROTATIONS IN THE (I,N) PLANE WHICH 1205C ELIMINATE ELEMENTS IN THE I-TH and N-TH PLANES, RESPECTIVELY. 1206C Q ITSELF IS NOT GIVEN, RATHER THE INFORMATION TO RECOVER THE 1207C GV, GW ROTATIONS IS SUPPLIED. 1208C 1209C THE SUBROUTINE STATEMENT IS 1210C 1211C subroutine SNQAQ(M,N,A,LDA,V,W) 1212C 1213C WHERE 1214C 1215C M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER 1216C OF ROWS OF A. 1217C 1218C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER 1219C OF COLUMNS OF A. 1220C 1221C A IS AN M BY N ARRAY. ON INPUT A MUST CONTAIN THE MATRIX 1222C TO BE POSTMULTIPLIED BY THE ORTHOGONAL MATRIX Q 1223C DESCRIBED ABOVE. ON OUTPUT A*Q HAS REPLACED A. 1224C 1225C LDA IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN M 1226C WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY A. 1227C 1228C V IS AN INPUT ARRAY OF LENGTH N. V(I) MUST CONTAIN THE 1229C INFORMATION NECESSARY TO RECOVER THE GIVENS ROTATION GV(I) 1230C DESCRIBED ABOVE. 1231C 1232C W IS AN INPUT ARRAY OF LENGTH N. W(I) MUST CONTAIN THE 1233C INFORMATION NECESSARY TO RECOVER THE GIVENS ROTATION GW(I) 1234C DESCRIBED ABOVE. 1235C 1236C SUBROUTINES CALLED 1237C 1238C FORTRAN-SUPPLIED ... abs,sqrt 1239C 1240C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. 1241C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE 1242C 1243C ********** 1244c ------------------------------------------------------------------ 1245 integer I,J,NMJ,NM1 1246 real VCOS,ONE,VSIN,TEMP 1247 parameter(ONE = 1.0e0) 1248C 1249C APPLY THE FIRST SET OF GIVENS ROTATIONS TO A. 1250C 1251 NM1 = N - 1 1252 IF (NM1 .lt. 1) GO TO 50 1253 DO 20 NMJ = 1, NM1 1254 J = N - NMJ 1255 IF (abs(V(J)) .gt. ONE) VCOS = ONE/V(J) 1256 IF (abs(V(J)) .gt. ONE) VSIN = sqrt(ONE-VCOS**2) 1257 IF (abs(V(J)) .le. ONE) VSIN = V(J) 1258 IF (abs(V(J)) .le. ONE) VCOS = sqrt(ONE-VSIN**2) 1259 DO 10 I = 1, M 1260 TEMP = VCOS*A(I,J) - VSIN*A(I,N) 1261 A(I,N) = VSIN*A(I,J) + VCOS*A(I,N) 1262 A(I,J) = TEMP 1263 10 CONTINUE 1264 20 CONTINUE 1265C 1266C APPLY THE SECOND SET OF GIVENS ROTATIONS TO A. 1267C 1268 DO 40 J = 1, NM1 1269 IF (abs(W(J)) .gt. ONE) VCOS = ONE/W(J) 1270 IF (abs(W(J)) .gt. ONE) VSIN = sqrt(ONE-VCOS**2) 1271 IF (abs(W(J)) .le. ONE) VSIN = W(J) 1272 IF (abs(W(J)) .le. ONE) VCOS = sqrt(ONE-VSIN**2) 1273 DO 30 I = 1, M 1274 TEMP = VCOS*A(I,J) + VSIN*A(I,N) 1275 A(I,N) = -VSIN*A(I,J) + VCOS*A(I,N) 1276 A(I,J) = TEMP 1277 30 CONTINUE 1278 40 CONTINUE 1279 50 CONTINUE 1280 RETURN 1281C 1282C Last line of subroutine SNQAQ. 1283C 1284 END 1285c ================================================================== 1286 subroutine SNQDOG(N,R,LR,DIAG,QTB,DELTA,X,NEWTOK,WA1,WA2, 1287 * SAMEJ, GNSTEP) 1288c>> 1992-01-03 CLL 1289 integer N,LR 1290 logical SAMEJ, NEWTOK 1291 real DELTA, GNSTEP(N) 1292 real R(LR),DIAG(N),QTB(N),X(N),WA1(N),WA2(N) 1293C ********** 1294C 1295C subroutine SNQDOG 1296C 1297c Given an M by N matrix A, an N by N nonsingular diagonal 1298c matrix D, an M-vector B, and a positive number DELTA, the 1299c problem is to determine the convex combination X of the 1300c gauss-newton and scaled gradient directions that minimizes 1301c (A*X - B) in the least squares sense, subject to the 1302c restriction that the euclidean norm of D*X be at most DELTA. 1303c 1304c This subroutine completes the solution of the problem 1305c if it is provided with the necessary information from the 1306c QR factorization of A. that is, if A = Q*R, where Q has 1307c orthogonal columns and R is an upper triangular matrix, 1308c then SNQDOG needs the full upper triangle of R and 1309c the first N components of (Q transpose)*B. 1310c 1311c The subroutine statement is 1312C 1313C subroutine SNQDOG(N,R,LR,DIAG,QTB,DELTA,X,WA1,WA2) 1314C 1315C where 1316C 1317c N is A POSITIVE INTEGER INPUT VARIABLE SET TO THE ORDER OF R. 1318C 1319c R() [in] An ARRAY OF LENGTH LR WHICH MUST CONTAIN THE UPPER 1320c TRIANGULAR MATRIX R STORED BY ROWS. 1321C 1322c LR is A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN 1323c (N*(N+1))/2. 1324C 1325c DIAG() [in] An ARRAY OF LENGTH N WHICH MUST CONTAIN THE 1326c DIAGONAL ELEMENTS OF THE MATRIX D. 1327C 1328c QTB() [in] An ARRAY OF LENGTH N WHICH MUST CONTAIN THE FIRST 1329c N ELEMENTS OF THE VECTOR (Q TRANSPOSE)*B. 1330C 1331c DELTA is a POSITIVE INPUT VARIABLE WHICH SPECIFIES AN UPPER 1332c BOUND ON THE EUCLIDEAN NORM OF D*X. 1333C 1334c X() [out] An ARRAY OF LENGTH N WHICH CONTAINS THE DESIRED 1335c CONVEX COMBINATION OF THE GAUSS-NEWTON DIRECTION and THE 1336c SCALED GRADIENT DIRECTION. 1337c 1338c NEWTOK [logical, out] True means the full Newton step was 1339c used. False means a modified, shorter, step was used. 1340c 1341c WA1() and WA2() are work arrays of length N. 1342c 1343c SAMEJ [logical, in] True means we have the same Jacobian matrix as 1344c on the previous call to this subr. The Gauss-Newton vector in 1345c GNSTEP() can be reused. 1346c 1347c GNSTEP() [inout] On return holds the Gauss-Newton vector. On entry 1348c with SAMEJ = .true., contains the GN vector from the previous call. 1349C ------------------------------------------------------------------ 1350C SUBPROGRAMS CALLED 1351C 1352c MINPACK-SUPPLIED ... R1MACH,SNRM2 1353C 1354C FORTRAN-SUPPLIED ... abs,max,min,sqrt 1355C 1356C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. 1357C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE 1358C 1359c ------------------------------------------------------------------ 1360 external R1MACH, SNRM2 1361 integer I,J,JJ,JP1,K,L, NP1 1362 real ALPHA,BNORM,EPSMCH,GNORM,ONE,QNORM,SGNORM,SUM 1363 real TEMP,ZERO 1364 real R1MACH,SNRM2 1365 parameter(ONE = 1.0e0, ZERO = 0.0e0) 1366 save EPSMCH 1367 data EPSMCH / 0.0e0 / 1368c ------------------------------------------------------------------ 1369C Set EPSMCH to the machine precision. 1370C 1371 if(EPSMCH .eq. 0.0e0) EPSMCH = R1MACH(4) 1372 if(.not. SAMEJ) then 1373C 1374C FIRST, CALCULATE THE GAUSS-NEWTON DIRECTION. 1375C 1376 NP1 = N+1 1377 JJ = (N*(N + 1))/2 + 1 1378 DO 50 K = 1, N 1379 J = NP1 - K 1380 JP1 = J + 1 1381 JJ = JJ - K 1382 L = JJ + 1 1383 SUM = ZERO 1384 DO 10 I = JP1, N 1385 SUM = SUM + R(L)*GNSTEP(I) 1386 L = L + 1 1387 10 CONTINUE 1388 TEMP = R(JJ) 1389 IF (TEMP .eq. ZERO) then 1390 L = J 1391 DO 30 I = 1, J-1 1392 TEMP = max(TEMP,abs(R(L))) 1393 L = L + N - I 1394 30 CONTINUE 1395 TEMP = EPSMCH*TEMP 1396 endif 1397 if (TEMP .eq. ZERO) then 1398 GNSTEP(J) = 0.0e0 1399 else 1400 GNSTEP(J) = (QTB(J) - SUM)/TEMP 1401 endif 1402 50 CONTINUE 1403 endif 1404C 1405C TEST WHETHER THE GAUSS-NEWTON DIRECTION is ACCEPTABLE. 1406C 1407 DO 60 J = 1, N 1408 WA1(J) = ZERO 1409 WA2(J) = DIAG(J)*GNSTEP(J) 1410 60 CONTINUE 1411 QNORM = SNRM2(N,WA2,1) 1412 NEWTOK = QNORM .le. DELTA 1413 if (NEWTOK) then 1414 do 65 J = 1,N 1415 X(J) = GNSTEP(J) 1416 65 continue 1417 go to 140 1418 endif 1419C 1420C THE GAUSS-NEWTON DIRECTION is NOT ACCEPTABLE. 1421C NEXT, CALCULATE THE SCALED GRADIENT DIRECTION. 1422C 1423 L = 1 1424 DO 80 J = 1, N 1425 TEMP = QTB(J) 1426 DO 70 I = J, N 1427 WA1(I) = WA1(I) + R(L)*TEMP 1428 L = L + 1 1429 70 CONTINUE 1430 WA1(J) = WA1(J)/DIAG(J) 1431 80 CONTINUE 1432C 1433C CALCULATE THE NORM OF THE SCALED GRADIENT and TEST FOR 1434C THE SPECIAL CASE IN WHICH THE SCALED GRADIENT is ZERO. 1435C 1436 GNORM = SNRM2(N,WA1,1) 1437 if (GNORM .eq. ZERO) then 1438 ALPHA = DELTA/QNORM 1439 do 85 J = 1, N 1440 X(J) = ALPHA*GNSTEP(J) 1441 85 continue 1442 go to 140 1443 endif 1444C 1445C CALCULATE THE POINT ALONG THE SCALED GRADIENT 1446C AT WHICH THE QUADRATIC is MINIMIZED. 1447C 1448 DO 90 J = 1, N 1449 WA1(J) = (WA1(J)/GNORM)/DIAG(J) 1450 90 CONTINUE 1451 L = 1 1452 DO 110 J = 1, N 1453 SUM = ZERO 1454 DO 100 I = J, N 1455 SUM = SUM + R(L)*WA1(I) 1456 L = L + 1 1457 100 CONTINUE 1458 WA2(J) = SUM 1459 110 CONTINUE 1460 TEMP = SNRM2(N,WA2,1) 1461 SGNORM = (GNORM/TEMP)/TEMP 1462C 1463C TEST WHETHER THE SCALED GRADIENT DIRECTION is ACCEPTABLE. 1464C 1465 ALPHA = ZERO 1466 if (SGNORM .lt. DELTA) then 1467C 1468C THE SCALED GRADIENT DIRECTION is NOT ACCEPTABLE. 1469C FINALLY, CALCULATE THE POINT ALONG THE dogleg 1470C AT WHICH THE QUADRATIC is MINIMIZED. 1471C 1472 BNORM = SNRM2(N,QTB,1) 1473 TEMP = (BNORM/GNORM)*(BNORM/QNORM)*(SGNORM/DELTA) 1474 TEMP = TEMP - (DELTA/QNORM)*(SGNORM/DELTA)**2 1475 * + sqrt((TEMP-(DELTA/QNORM))**2 1476 * +(ONE-(DELTA/QNORM)**2)*(ONE-(SGNORM/DELTA)**2)) 1477 ALPHA = ((DELTA/QNORM)*(ONE - (SGNORM/DELTA)**2))/TEMP 1478 endif 1479C 1480C FORM APPROPRIATE CONVEX COMBINATION OF THE GAUSS-NEWTON 1481C DIRECTION and THE SCALED GRADIENT DIRECTION. 1482C 1483 TEMP = (ONE - ALPHA)*min(SGNORM,DELTA) 1484 DO 130 J = 1, N 1485 X(J) = TEMP*WA1(J) + ALPHA*GNSTEP(J) 1486 130 CONTINUE 1487 140 CONTINUE 1488 RETURN 1489C 1490C Last line of subroutine SNQDOG. 1491C 1492 END 1493 subroutine SNQQFM(M,N,Q,LDQ,WA) 1494 integer M,N,LDQ 1495 real Q(LDQ,M),WA(M) 1496C ********** 1497C 1498C SUBROUTINE SNQQFM 1499C 1500C THIS SUBROUTINE PROCEEDS FROM THE COMPUTED QR FACTORIZATION OF 1501C AN M BY N MATRIX A TO ACCUMULATE THE M BY M ORTHOGONAL MATRIX 1502C Q FROM ITS FACTORED FORM. 1503C 1504C THE SUBROUTINE STATEMENT IS 1505C 1506C subroutine SNQQFM(M,N,Q,LDQ,WA) 1507C 1508C WHERE 1509C 1510C M is A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER 1511C OF ROWS OF A and THE ORDER OF Q. 1512C 1513C N is A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER 1514C OF COLUMNS OF A. 1515C 1516C Q is AN M BY M ARRAY. ON INPUT THE FULL LOWER TRAPEZOID IN 1517C THE FIRST MIN(M,N) COLUMNS OF Q CONTAINS THE FACTORED FORM. 1518C ON OUTPUT Q HAS BEEN ACCUMULATED INTO A SQUARE MATRIX. 1519C 1520C LDQ is A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN M 1521C WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY Q. 1522C 1523C WA is A WORK ARRAY OF LENGTH M. 1524C 1525C SUBPROGRAMS CALLED 1526C 1527C FORTRAN-SUPPLIED ... min 1528C 1529C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. 1530C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE 1531C 1532c ------------------------------------------------------------------ 1533 integer I,J,JM1,K,L,MINMN,NP1 1534 real ONE,SUM,TEMP,ZERO 1535 parameter(ONE = 1.0e0, ZERO = 0.0e0) 1536C 1537C ZERO OUT UPPER TRIANGLE OF Q IN THE FIRST MIN(M,N) COLUMNS. 1538C 1539 MINMN = min(M,N) 1540 IF (MINMN .lt. 2) GO TO 30 1541 DO 20 J = 2, MINMN 1542 JM1 = J - 1 1543 DO 10 I = 1, JM1 1544 Q(I,J) = ZERO 1545 10 CONTINUE 1546 20 CONTINUE 1547 30 CONTINUE 1548C 1549C INITIALIZE REMAINING COLUMNS TO THOSE OF THE IDENTITY MATRIX. 1550C 1551 NP1 = N + 1 1552 IF (M .lt. NP1) GO TO 60 1553 DO 50 J = NP1, M 1554 DO 40 I = 1, M 1555 Q(I,J) = ZERO 1556 40 CONTINUE 1557 Q(J,J) = ONE 1558 50 CONTINUE 1559 60 CONTINUE 1560C 1561C ACCUMULATE Q FROM ITS FACTORED FORM. 1562C 1563 DO 120 L = 1, MINMN 1564 K = MINMN - L + 1 1565 DO 70 I = K, M 1566 WA(I) = Q(I,K) 1567 Q(I,K) = ZERO 1568 70 CONTINUE 1569 Q(K,K) = ONE 1570 IF (WA(K) .eq. ZERO) GO TO 110 1571 DO 100 J = K, M 1572 SUM = ZERO 1573 DO 80 I = K, M 1574 SUM = SUM + Q(I,J)*WA(I) 1575 80 CONTINUE 1576 TEMP = SUM/WA(K) 1577 DO 90 I = K, M 1578 Q(I,J) = Q(I,J) - TEMP*WA(I) 1579 90 CONTINUE 1580 100 CONTINUE 1581 110 CONTINUE 1582 120 CONTINUE 1583 RETURN 1584C 1585C Last line of subroutine SNQQFM. 1586C 1587 END 1588 subroutine SNQQRF(M,N,A,LDA,PIVOT,IPVT,LIPVT,RDIAG,ACNORM,WA) 1589 integer M,N,LDA,LIPVT 1590 integer IPVT(LIPVT) 1591 logical PIVOT 1592 real A(LDA,N),RDIAG(N),ACNORM(N),WA(N) 1593C ********** 1594C 1595C SUBROUTINE SNQQRF 1596C 1597C THIS SUBROUTINE USES HOUSEHOLDER TRANSFORMATIONS WITH COLUMN 1598C PIVOTING (OPTIONAL) TO COMPUTE A QR FACTORIZATION OF THE 1599C M BY N MATRIX A. THAT IS, SNQQRF DETERMINES AN ORTHOGONAL 1600C MATRIX Q, A PERMUTATION MATRIX P, and AN UPPER TRAPEZOIDAL 1601C MATRIX R WITH DIAGONAL ELEMENTS OF NONINCREASING MAGNITUDE, 1602C SUCH THAT A*P = Q*R. THE HOUSEHOLDER TRANSFORMATION FOR 1603C COLUMN K, K = 1,2,...,MIN(M,N), is OF THE FORM 1604C 1605C T 1606C I - (1/U(K))*U*U 1607C 1608C WHERE U HAS ZEROS IN THE FIRST K-1 POSITIONS. THE FORM OF 1609C THIS TRANSFORMATION and THE METHOD OF PIVOTING FIRST 1610C APPEARED IN THE CORRESPONDING LINPACK SUBROUTINE. 1611C 1612C THE SUBROUTINE STATEMENT IS 1613C 1614C subroutine SNQQRF(M,N,A,LDA,PIVOT,IPVT,LIPVT,RDIAG,ACNORM,WA) 1615C 1616C WHERE 1617C 1618C M is A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER 1619C OF ROWS OF A. 1620C 1621C N is A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER 1622C OF COLUMNS OF A. 1623C 1624C A is AN M BY N ARRAY. ON INPUT A CONTAINS THE MATRIX FOR 1625C WHICH THE QR FACTORIZATION is TO BE COMPUTED. ON OUTPUT 1626C THE STRICT UPPER TRAPEZOIDAL PART OF A CONTAINS THE STRICT 1627C UPPER TRAPEZOIDAL PART OF R, and THE LOWER TRAPEZOIDAL 1628C PART OF A CONTAINS A FACTORED FORM OF Q (THE NON-TRIVIAL 1629C ELEMENTS OF THE U VECTORS DESCRIBED ABOVE). 1630C 1631C LDA is A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN M 1632C WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY A. 1633C 1634C PIVOT is A LOGICAL INPUT VARIABLE. IF PIVOT is SET TRUE, 1635C THEN COLUMN PIVOTING is ENFORCED. IF PIVOT is SET FALSE, 1636C THEN NO COLUMN PIVOTING is DONE. 1637C 1638C IPVT is AN INTEGER OUTPUT ARRAY OF LENGTH LIPVT. IPVT 1639C DEFINES THE PERMUTATION MATRIX P SUCH THAT A*P = Q*R. 1640C COLUMN J OF P is COLUMN IPVT(J) OF THE IDENTITY MATRIX. 1641C IF PIVOT is FALSE, IPVT is NOT REFERENCED. 1642C 1643C LIPVT is A POSITIVE INTEGER INPUT VARIABLE. IF PIVOT is FALSE, 1644C THEN LIPVT MAY BE AS SMALL AS 1. IF PIVOT is TRUE, THEN 1645C LIPVT MUST BE AT LEAST N. 1646C 1647C RDIAG is AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS THE 1648C DIAGONAL ELEMENTS OF R. 1649C 1650C ACNORM is AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS THE 1651C NORMS OF THE CORRESPONDING COLUMNS OF THE INPUT MATRIX A. 1652C IF THIS INFORMATION is NOT NEEDED, THEN ACNORM CAN COINCIDE 1653C WITH RDIAG. 1654C 1655C WA is A WORK ARRAY OF LENGTH N. IF PIVOT is FALSE, THEN WA 1656C CAN COINCIDE WITH RDIAG. 1657C 1658C SUBPROGRAMS CALLED 1659C 1660C MINPACK-SUPPLIED ... R1MACH,SNRM2 1661C 1662C FORTRAN-SUPPLIED ... max,sqrt,min 1663C 1664C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. 1665C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE 1666C 1667C ********** 1668c ------------------------------------------------------------------ 1669 external R1MACH, SNRM2 1670 integer I,J,JP1,K,KMAX,MINMN 1671 real AJNORM,EPSMCH,ONE,P05,SUM,TEMP,ZERO 1672 real R1MACH,SNRM2 1673 parameter(ONE = 1.0e0, P05 = 0.05e0, ZERO = 0.0e0) 1674C 1675C EPSMCH is THE MACHINE PRECISION. 1676C 1677 EPSMCH = R1MACH(4) 1678C 1679C COMPUTE THE INITIAL COLUMN NORMS and INITIALIZE SEVERAL ARRAYS. 1680C 1681 DO 10 J = 1, N 1682 ACNORM(J) = SNRM2(M,A(1,J),1) 1683 RDIAG(J) = ACNORM(J) 1684 WA(J) = RDIAG(J) 1685 IF (PIVOT) IPVT(J) = J 1686 10 CONTINUE 1687C 1688C REDUCE A TO R WITH HOUSEHOLDER TRANSFORMATIONS. 1689C 1690 MINMN = min(M,N) 1691 DO 110 J = 1, MINMN 1692 IF (.NOT.PIVOT) GO TO 40 1693C 1694C BRING THE COLUMN OF LARGEST NORM INTO THE PIVOT POSITION. 1695C 1696 KMAX = J 1697 DO 20 K = J, N 1698 IF (RDIAG(K) .gt. RDIAG(KMAX)) KMAX = K 1699 20 CONTINUE 1700 IF (KMAX .eq. J) GO TO 40 1701 DO 30 I = 1, M 1702 TEMP = A(I,J) 1703 A(I,J) = A(I,KMAX) 1704 A(I,KMAX) = TEMP 1705 30 CONTINUE 1706 RDIAG(KMAX) = RDIAG(J) 1707 WA(KMAX) = WA(J) 1708 K = IPVT(J) 1709 IPVT(J) = IPVT(KMAX) 1710 IPVT(KMAX) = K 1711 40 CONTINUE 1712C 1713C COMPUTE THE HOUSEHOLDER TRANSFORMATION TO REDUCE THE 1714C J-TH COLUMN OF A TO A MULTIPLE OF THE J-TH UNIT VECTOR. 1715C 1716 AJNORM = SNRM2(M-J+1,A(J,J),1) 1717 IF (AJNORM .eq. ZERO) GO TO 100 1718 IF (A(J,J) .lt. ZERO) AJNORM = -AJNORM 1719 DO 50 I = J, M 1720 A(I,J) = A(I,J)/AJNORM 1721 50 CONTINUE 1722 A(J,J) = A(J,J) + ONE 1723C 1724C APPLY THE TRANSFORMATION TO THE REMAINING COLUMNS 1725C and UPDATE THE NORMS. 1726C 1727 JP1 = J + 1 1728 IF (N .lt. JP1) GO TO 100 1729 DO 90 K = JP1, N 1730 SUM = ZERO 1731 DO 60 I = J, M 1732 SUM = SUM + A(I,J)*A(I,K) 1733 60 CONTINUE 1734 TEMP = SUM/A(J,J) 1735 DO 70 I = J, M 1736 A(I,K) = A(I,K) - TEMP*A(I,J) 1737 70 CONTINUE 1738 IF (.NOT.PIVOT .or. RDIAG(K) .eq. ZERO) GO TO 80 1739 TEMP = A(J,K)/RDIAG(K) 1740 RDIAG(K) = RDIAG(K)*sqrt(max(ZERO,ONE-TEMP**2)) 1741 IF (P05*(RDIAG(K)/WA(K))**2 .gt. EPSMCH) GO TO 80 1742 RDIAG(K) = SNRM2(M-J,A(JP1,K),1) 1743 WA(K) = RDIAG(K) 1744 80 CONTINUE 1745 90 CONTINUE 1746 100 CONTINUE 1747 RDIAG(J) = -AJNORM 1748 110 CONTINUE 1749 RETURN 1750C 1751C Last line of subroutine SNQQRF. 1752C 1753 END 1754 subroutine SNQUPD(M,N,S,LS,U,V,W,SING) 1755 integer M,N,LS 1756 logical SING 1757 real S(LS),U(M),V(N),W(M) 1758C ********** 1759C 1760C SUBROUTINE SNQUPD 1761C 1762C GIVEN AN M BY N LOWER TRAPEZOIDAL MATRIX S, AN M-VECTOR U, 1763C and AN N-VECTOR V, THE PROBLEM is TO DETERMINE AN 1764C ORTHOGONAL MATRIX Q SUCH THAT 1765C 1766C T 1767C (S + U*V )*Q 1768C 1769C is AGAIN LOWER TRAPEZOIDAL. 1770C 1771C THIS SUBROUTINE DETERMINES Q AS THE PRODUCT OF 2*(N - 1) 1772C TRANSFORMATIONS 1773C 1774C GV(N-1)*...*GV(1)*GW(1)*...*GW(N-1) 1775C 1776C WHERE GV(I), GW(I) ARE GIVENS ROTATIONS IN THE (I,N) PLANE 1777C WHICH ELIMINATE ELEMENTS IN THE I-TH and N-TH PLANES, 1778C RESPECTIVELY. Q ITSELF is NOT ACCUMULATED, RATHER THE 1779C INFORMATION TO RECOVER THE GV, GW ROTATIONS is RETURNED. 1780C 1781C THE SUBROUTINE STATEMENT IS 1782C 1783C subroutine SNQUPD(M,N,S,LS,U,V,W,SING) 1784C 1785C WHERE 1786C 1787C M is A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER 1788C OF ROWS OF S. 1789C 1790C N is A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER 1791C OF COLUMNS OF S. N MUST NOT EXCEED M. 1792C 1793C S is AN ARRAY OF LENGTH LS. ON INPUT S MUST CONTAIN THE LOWER 1794C TRAPEZOIDAL MATRIX S STORED BY COLUMNS. ON OUTPUT S CONTAINS 1795C THE LOWER TRAPEZOIDAL MATRIX PRODUCED AS DESCRIBED ABOVE. 1796C 1797C LS is A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN 1798C (N*(2*M-N+1))/2. 1799C 1800C U is AN INPUT ARRAY OF LENGTH M WHICH MUST CONTAIN THE 1801C VECTOR U. 1802C 1803C V is AN ARRAY OF LENGTH N. ON INPUT V MUST CONTAIN THE VECTOR 1804C V. ON OUTPUT V(I) CONTAINS THE INFORMATION NECESSARY TO 1805C RECOVER THE GIVENS ROTATION GV(I) DESCRIBED ABOVE. 1806C 1807C W is AN OUTPUT ARRAY OF LENGTH M. W(I) CONTAINS INFORMATION 1808C NECESSARY TO RECOVER THE GIVENS ROTATION GW(I) DESCRIBED 1809C ABOVE. 1810C 1811C SING is A LOGICAL OUTPUT VARIABLE. SING is SET TRUE IF ANY 1812C OF THE DIAGONAL ELEMENTS OF THE OUTPUT S ARE ZERO. OTHERWISE 1813C SING is SET FALSE. 1814C 1815C SUBPROGRAMS CALLED 1816C 1817C MINPACK-SUPPLIED ... R1MACH 1818C 1819C FORTRAN-SUPPLIED ... abs,sqrt 1820C 1821C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. 1822C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE, 1823C JOHN L. NAZARETH 1824C 1825C ********** 1826c ------------------------------------------------------------------ 1827 external R1MACH 1828 integer I,J,JJ,L,NMJ,NM1 1829 real VCOS,COTAN,GIANT,ONE,P5,P25,VSIN,VTAN,TAU,TEMP, 1830 * ZERO 1831 real R1MACH 1832 parameter(ONE = 1.0e0, P5 = 0.5e0, P25 = 0.25e0, ZERO = 0.0e0) 1833 save GIANT 1834 data GIANT / 0.0e0 / 1835C 1836C GIANT is THE LARGEST MAGNITUDE. 1837C 1838 if(GIANT .eq. 0.0e0) GIANT = R1MACH(2) 1839C 1840C INITIALIZE THE DIAGONAL ELEMENT POINTER. 1841C 1842 JJ = (N*(2*M - N + 1))/2 - (M - N) 1843C 1844C MOVE THE NONTRIVIAL PART OF THE LAST COLUMN OF S INTO W. 1845C 1846 L = JJ 1847 DO 10 I = N, M 1848 W(I) = S(L) 1849 L = L + 1 1850 10 CONTINUE 1851C 1852C ROTATE THE VECTOR V INTO A MULTIPLE OF THE N-TH UNIT VECTOR 1853C IN SUCH A WAY THAT A SPIKE is INTRODUCED INTO W. 1854C 1855 NM1 = N - 1 1856 IF (NM1 .lt. 1) GO TO 70 1857 DO 60 NMJ = 1, NM1 1858 J = N - NMJ 1859 JJ = JJ - (M - J + 1) 1860 W(J) = ZERO 1861 IF (V(J) .eq. ZERO) GO TO 50 1862C 1863C DETERMINE A GIVENS ROTATION WHICH ELIMINATES THE 1864C J-TH ELEMENT OF V. 1865C 1866 IF (abs(V(N)) .ge. abs(V(J))) GO TO 20 1867 COTAN = V(N)/V(J) 1868 VSIN = P5/sqrt(P25+P25*COTAN**2) 1869 VCOS = VSIN*COTAN 1870 TAU = ONE 1871 IF (abs(VCOS)*GIANT .gt. ONE) TAU = ONE/VCOS 1872 GO TO 30 1873 20 CONTINUE 1874 VTAN = V(J)/V(N) 1875 VCOS = P5/sqrt(P25+P25*VTAN**2) 1876 VSIN = VCOS*VTAN 1877 TAU = VSIN 1878 30 CONTINUE 1879C 1880C APPLY THE TRANSFORMATION TO V and STORE THE INFORMATION 1881C NECESSARY TO RECOVER THE GIVENS ROTATION. 1882C 1883 V(N) = VSIN*V(J) + VCOS*V(N) 1884 V(J) = TAU 1885C 1886C APPLY THE TRANSFORMATION TO S and EXTEND THE SPIKE IN W. 1887C 1888 L = JJ 1889 DO 40 I = J, M 1890 TEMP = VCOS*S(L) - VSIN*W(I) 1891 W(I) = VSIN*S(L) + VCOS*W(I) 1892 S(L) = TEMP 1893 L = L + 1 1894 40 CONTINUE 1895 50 CONTINUE 1896 60 CONTINUE 1897 70 CONTINUE 1898C 1899C ADD THE SPIKE FROM THE RANK 1 UPDATE TO W. 1900C 1901 DO 80 I = 1, M 1902 W(I) = W(I) + V(N)*U(I) 1903 80 CONTINUE 1904C 1905C ELIMINATE THE SPIKE. 1906C 1907 SING = .FALSE. 1908 IF (NM1 .lt. 1) GO TO 140 1909 DO 130 J = 1, NM1 1910 IF (W(J) .eq. ZERO) GO TO 120 1911C 1912C DETERMINE A GIVENS ROTATION WHICH ELIMINATES THE 1913C J-TH ELEMENT OF THE SPIKE. 1914C 1915 IF (abs(S(JJ)) .ge. abs(W(J))) GO TO 90 1916 COTAN = S(JJ)/W(J) 1917 VSIN = P5/sqrt(P25+P25*COTAN**2) 1918 VCOS = VSIN*COTAN 1919 TAU = ONE 1920 IF (abs(VCOS)*GIANT .gt. ONE) TAU = ONE/VCOS 1921 GO TO 100 1922 90 CONTINUE 1923 VTAN = W(J)/S(JJ) 1924 VCOS = P5/sqrt(P25+P25*VTAN**2) 1925 VSIN = VCOS*VTAN 1926 TAU = VSIN 1927 100 CONTINUE 1928C 1929C APPLY THE TRANSFORMATION TO S and REDUCE THE SPIKE IN W. 1930C 1931 L = JJ 1932 DO 110 I = J, M 1933 TEMP = VCOS*S(L) + VSIN*W(I) 1934 W(I) = -VSIN*S(L) + VCOS*W(I) 1935 S(L) = TEMP 1936 L = L + 1 1937 110 CONTINUE 1938C 1939C STORE THE INFORMATION NECESSARY TO RECOVER THE 1940C GIVENS ROTATION. 1941C 1942 W(J) = TAU 1943 120 CONTINUE 1944C 1945C TEST FOR ZERO DIAGONAL ELEMENTS IN THE OUTPUT S. 1946C 1947 IF (S(JJ) .eq. ZERO) SING = .TRUE. 1948 JJ = JJ + (M - J + 1) 1949 130 CONTINUE 1950 140 CONTINUE 1951C 1952C MOVE W BACK INTO THE LAST COLUMN OF THE OUTPUT S. 1953C 1954 L = JJ 1955 DO 150 I = N, M 1956 S(L) = W(I) 1957 L = L + 1 1958 150 CONTINUE 1959 IF (S(JJ) .eq. ZERO) SING = .TRUE. 1960 RETURN 1961C 1962C Last line of subroutine SNQUPD. 1963C 1964 END 1965