1 SUBROUTINE NLEQ2(N,FCN,JAC,X,XSCAL,RTOL,IOPT,IERR, 2 $LIWK,IWK,LRWK,RWK) 3C* Begin Prologue NLEQ2 4 INTEGER N 5 EXTERNAL FCN,JAC 6 DOUBLE PRECISION X(N),XSCAL(N) 7 DOUBLE PRECISION RTOL 8 INTEGER IOPT(50) 9 INTEGER IERR 10 INTEGER LIWK 11 INTEGER IWK(LIWK) 12 INTEGER LRWK 13 DOUBLE PRECISION RWK(LRWK) 14C ------------------------------------------------------------ 15C 16C* Title 17C 18C Numerical solution of nonlinear (NL) equations (EQ) 19C especially designed for numerically sensitive problems. 20C 21C* Written by U. Nowak, L. Weimann 22C* Purpose Solution of systems of highly nonlinear equations 23C* Method Damped affine invariant Newton method with rank- 24C strategy (see references below) 25C* Category F2a. - Systems of nonlinear equations 26C* Keywords Nonlinear equations, Newton methods 27C* Version 2.3 28C* Revision September 1991 29C* Latest Change July 2000 30C* Library CodeLib 31C* Code Fortran 77, Double Precision 32C* Environment Standard Fortran 77 environment on PC's, 33C workstations and hosts. 34C* Copyright (c) Konrad-Zuse-Zentrum fuer 35C Informationstechnik Berlin (ZIB) 36C Takustrasse 7, D-14195 Berlin-Dahlem 37C phone : + 49/30/84185-0 38C fax : + 49/30/84185-125 39C* Contact Lutz Weimann 40C ZIB, Division Scientific Computing, 41C Department Scientific Software 42C phone : + 49/30/84185-185 43C fax : + 49/30/84185-107 44C e-mail: weimann@zib.de 45C 46C* References: 47C 48C /1/ P. Deuflhard: 49C Newton Techniques for Highly Nonlinear Problems - 50C Theory and Algorithms. 51C Academic press Inc. (To be published) 52C 53C /2/ U. Nowak, L. Weimann: 54C A Family of Newton Codes for Systems of Highly Nonlinear 55C Equations - Algorithm, Implementation, Application. 56C ZIB, Technical Report TR 90-10 (December 1990) 57C 58C --------------------------------------------------------------- 59C 60C* Licence 61C You may use or modify this code for your own non commercial 62C purposes for an unlimited time. 63C In any case you should not deliver this code without a special 64C permission of ZIB. 65C In case you intend to use the code commercially, we oblige you 66C to sign an according licence agreement with ZIB. 67C 68C* Warranty 69C This code has been tested up to a certain level. Defects and 70C weaknesses, which may be included in the code, do not establish 71C any warranties by ZIB. ZIB does not take over any liabilities 72C which may follow from acquisition or application of this code. 73C 74C* Software status 75C This code is under care of ZIB and belongs to ZIB software class 1. 76C 77C ------------------------------------------------------------ 78C 79C* Summary: 80C ======== 81C Damped Newton-algorithm with rank strategy for systems of 82C highly nonlinear equations - damping strategy due to Ref.(1). 83C 84C (The iteration is done by subroutine N2INT currently. NLEQ2 85C itself does some house keeping and builds up workspace.) 86C 87C Jacobian approximation by numerical differences or user 88C supplied subroutine JAC. 89C 90C The numerical solution of the arising linear equations is 91C done by means of the subroutines DECCON and SOLCON (QR de- 92C composition with subcondition estimation, rank decision and 93C computation of the rank-deficient pseudoinverse) . 94C For special purposes these routines may be substituted. 95C 96C This is a driver routine for the core solver N2INT. 97C 98C ------------------------------------------------------------ 99C 100C* Parameters list description (* marks inout parameters) 101C ====================================================== 102C 103C* External subroutines (to be supplied by the user) 104C ================================================= 105C 106C (Caution: Arguments declared as (input) must not 107C be altered by the user subroutines ! ) 108C 109C FCN(N,X,F,IFAIL) Ext Function subroutine 110C N Int Number of vector components (input) 111C X(N) Dble Vector of unknowns (input) 112C F(N) Dble Vector of function values (output) 113C IFAIL Int FCN evaluation-failure indicator. (output) 114C On input: Has always value 0 (zero). 115C On output: Indicates failure of FCN eval- 116C uation, if having a value <= 2. 117C If <0: NLEQ2 will be terminated with 118C error code = 82, and IFAIL stored 119C to IWK(23). 120C If =1: A new trial Newton iterate will 121C computed, with the damping factor 122C reduced to it's half. 123C If =2: A new trial Newton iterate will 124C computed, with the damping factor 125C reduced by a reduct. factor, which 126C must be output through F(1) by FCN, 127C and it's value must be >0 and < 1. 128C Note, that if IFAIL = 1 or 2, additional 129C conditions concerning the damping factor, 130C e.g. the minimum damping factor or the 131C bounded damping strategy may also influ- 132C ence the value of the reduced damping 133C factor. 134C 135C 136C JAC(N,LDJAC,X,DFDX,IFAIL) 137C Ext Jacobian matrix subroutine 138C N Int Number of vector components (input) 139C LDJAC Int Leading dimension of Jacobian array 140C (input) 141C X(N) Dble Vector of unknowns (input) 142C DFDX(LDJAC,N) Dble DFDX(i,k): partial derivative of 143C I-th component of FCN with respect 144C to X(k) (output) 145C IFAIL Int JAC evaluation-failure indicator. 146C (output) 147C Has always value 0 (zero) on input. 148C Indicates failure of JAC evaluation 149C and causes termination of NLEQ2, 150C if set to a negative value on output 151C 152C 153C* Input parameters of NLEQ2 154C ========================= 155C 156C N Int Number of unknowns 157C * X(N) Dble Initial estimate of the solution 158C * XSCAL(N) Dble User scaling (lower threshold) of the 159C iteration vector X(N) 160C * RTOL Dble Required relative precision of 161C solution components - 162C RTOL.GE.EPMACH*TEN*N 163C * IOPT(50) Int Array of run-time options. Set to zero 164C to get default values (details see below) 165C 166C* Output parameters of NLEQ2 167C ========================== 168C 169C * X(N) Dble Solution values ( or final values, 170C respectively ) 171C * XSCAL(N) Dble After return with IERR.GE.0, it contains 172C the latest internal scaling vector used 173C After return with IERR.EQ.-1 in onestep- 174C mode it contains a possibly adapted 175C (as described below) user scaling vector: 176C If (XSCAL(I).LT. SMALL) XSCAL(I) = SMALL , 177C If (XSCAL(I).GT. GREAT) XSCAL(I) = GREAT . 178C For SMALL and GREAT, see section machine 179C constants below and regard note 1. 180C * RTOL Dble Finally achieved (relative) accuracy 181C The estimated absolute error of component i 182C of x_out is approximately given by 183C abs_err(i) = RTOL * XSCAL_out(i) , 184C where (approximately) 185C XSCAL_out(i) = 186C max(abs(X_out(i)),XSCAL_in(i)). 187C IERR Int Return value parameter 188C =-1 sucessfull completion of one iteration 189C step, subsequent iterations are needed 190C to get a solution. (stepwise mode only) 191C = 0 successfull completion of iteration 192C > 0 see list of error messages below 193C 194C Note 1. 195C The machine dependent values SMALL, GREAT and EPMACH are 196C gained from calls of the machine constants function D1MACH. 197C As delivered, this function is adapted to use constants 198C suitable for all machines with IEEE arithmetic. If you use 199C another type of machine, you have to change the DATA state- 200C ments for IEEE arithmetic in D1MACH into comments and to 201C uncomment the set of DATA statements suitable for your machine. 202C 203C* Workspace parameters of NLEQ2 204C ============================= 205C 206C LIWK Int Declared dimension of integer 207C workspace. 208C Required minimum (for standard linear system 209C solver) : N+52 210C * IWK(LIWK) Int Integer Workspace 211C LRWK Int Declared dimension of real workspace. 212C Required minimum (for standard linear system 213C solver and Jacobian computed by numerical 214C approximation - if the Jacobian is computed 215C by a user subroutine JAC, decrease the 216C expression noted below by N): 217C (N+NBROY+15)*N+61 218C NBROY = Maximum number of Broyden steps 219C (Default: if Broyden steps are enabled, e.g. 220C IOPT(32)=1 - 221C NBROY=MAX(N,10) 222C else (if IOPT(32)=0) - 223C NBROY=0 ; 224C see equally named IWK-field below) 225C * RWK(LRWK) Dble Real Workspace 226C 227C Note 2a. A test on sufficient workspace is made. If this 228C test fails, IERR is set to 10 and an error-message 229C is issued from which the minimum of required 230C workspace size can be obtained. 231C 232C Note 2b. The first 50 elements of IWK and RWK are partially 233C used as input for internal algorithm parameters (for 234C details, see below). In order to set the default values 235C of these parameters, the fields must be set to zero. 236C Therefore, it's recommanded always to initialize the 237C first 50 elements of both workspaces to zero. 238C 239C* Options IOPT: 240C ============= 241C 242C Pos. Name Default Meaning 243C 244C 1 QSUCC 0 =0 (.FALSE.) initial call: 245C NLEQ2 is not yet initialized, i.e. this is 246C the first call for this nonlinear system. 247C At successfull return with MODE=1, 248C QSUCC is set to 1. 249C =1 (.TRUE.) successive call: 250C NLEQ2 is initialized already and is now 251C called to perform one or more following 252C Newton-iteration steps. 253C ATTENTION: 254C Don't destroy the contents of 255C IOPT(i) for 1 <= i <= 50 , 256C IWK(j) for 1 <= j < NIWKFR and 257C RWK(k) for 1 <= k < NRWKFR. 258C (Nevertheless, some of the options, e.g. 259C FCMIN, SIGMA, MPR..., can be modified 260C before successive calls.) 261C 2 MODE 0 =0 Standard mode initial call: 262C Return when the required accuracy for the 263C iteration vector is reached. User defined 264C parameters are evaluated and checked. 265C Standard mode successive call: 266C If NLEQ2 was called previously with MODE=1, 267C it performs all remaining iteration steps. 268C =1 Stepwise mode: 269C Return after one Newton iteration step. 270C 3 JACGEN 0 Method of Jacobian generation 271C =0 Standard method is JACGEN=2 272C =1 User supplied subroutine JAC will be 273C called to generate Jacobian matrix 274C =2 Jacobian approximation by numerical 275C differentation (see subroutine N2JAC) 276C =3 Jacobian approximation by numerical 277C differentation with feedback control 278C (see subroutine N2JCF) 279C 4..8 Reserved 280C 9 ISCAL 0 Determines how to scale the iterate-vector: 281C =0 The user supplied scaling vector XSCAL is 282C used as a (componentwise) lower threshold 283C of the current scaling vector 284C =1 The vector XSCAL is always used as the 285C current scaling vector 286C 10 Reserved 287C 11 MPRERR 0 Print error messages 288C =0 No output 289C =1 Error messages 290C =2 Warnings additionally 291C =3 Informal messages additionally 292C 12 LUERR 6 Logical unit number for error messages 293C 13 MPRMON 0 Print iteration Monitor 294C =0 No output 295C =1 Standard output 296C =2 Summary iteration monitor additionally 297C =3 Detailed iteration monitor additionally 298C =4,5,6 Outputs with increasing level addi- 299C tional increasing information for code 300C testing purposes. Level 6 produces 301C in general extremely large output! 302C 14 LUMON 6 Logical unit number for iteration monitor 303C 15 MPRSOL 0 Print solutions 304C =0 No output 305C =1 Initial values and solution values 306C =2 Intermediate iterates additionally 307C 16 LUSOL 6 Logical unit number for solutions 308C 17..18 Reserved 309C 19 MPRTIM 0 Output level for the time monitor 310C = 0 : no time measurement and no output 311C = 1 : time measurement will be done and 312C summary output will be written - 313C regard note 4a. 314C 20 LUTIM 6 Logical output unit for time monitor 315C 21..30 Reserved 316C 31 NONLIN 3 Problem type specification 317C =1 Linear problem 318C Warning: If specified, no check will be 319C done, if the problem is really linear, and 320C NLEQ2 terminates unconditionally after one 321C Newton-iteration step. 322C =2 Mildly nonlinear problem 323C =3 Highly nonlinear problem 324C =4 Extremely nonlinear problem 325C 32 QRANK1 0 =0 (.FALSE.) Rank-1 updates by Broyden- 326C approximation are inhibited. 327C =1 (.TRUE.) Rank-1 updates by Broyden- 328C approximation are allowed. 329C 33..34 Reserved 330C 35 QNSCAL 0 Inhibit automatic row scaling: 331C =0 (.FALSE.) Automatic row scaling of 332C the linear system is activ: 333C Rows i=1,...,N will be divided by 334C max j=1,...,N (abs(a(i,j))) 335C =1 (.TRUE.) No row scaling of the linear 336C system. Recommended only for well row- 337C scaled nonlinear systems. 338C 36..37 Reserved 339C 38 IBDAMP Bounded damping strategy switch: 340C =0 The default switch takes place, dependent 341C on the setting of NONLIN (=IOPT(31)): 342C NONLIN = 0,1,2,3 -> IBDAMP = off , 343C NONLIN = 4 -> IBDAMP = on 344C =1 means always IBDAMP = on 345C =2 means always IBDAMP = off 346C 39 IORMON Convergence order monitor 347C =0 Standard option is IORMON=2 348C =1 Convergence order is not checked, 349C the iteration will be always proceeded 350C until the solution has the required 351C precision RTOL (or some error condition 352C occured) 353C =2 Use additional 'weak stop' criterion: 354C Convergence order is monitored 355C and termination due to slowdown of the 356C convergence may occur. 357C =3 Use additional 'hard stop' criterion: 358C Convergence order is monitored 359C and termination due to superlinear 360C convergence slowdown may occur. 361C In case of termination due to convergence 362C slowdown, the warning code IERR=4 will be 363C set. 364C In cases, where the Newton iteration con- 365C verges but superlinear convergence order has 366C never been detected, the warning code IERR=5 367C is returned. 368C 40..45 Reserved 369C 46..50 User options (see note 4b) 370C 371C Note 3: 372C If NLEQ2 terminates with IERR=2 (maximum iterations) 373C or IERR=3 (small damping factor), you may try to continue 374C the iteration by increasing NITMAX or decreasing FCMIN 375C (see RWK) and setting QSUCC to 1. 376C 377C Note 4a: 378C The integrated time monitor calls the machine dependent 379C subroutine SECOND to get the current time stamp in form 380C of a real number (Single precision). As delivered, this 381C subroutine always return 0.0 as time stamp value. Refer 382C to the compiler- or library manual of the FORTRAN compiler 383C which you currently use to find out how to get the current 384C time stamp on your machine. 385C 386C Note 4b: 387C The user options may be interpreted by the user replacable 388C routines N2SOUT, N2FACT, N2SOLV - the distributed version 389C of N2SOUT currently uses IOPT(46) as follows: 390C 0 = standard plotdata output (may be postprocessed by a user- 391C written graphical program) 392C 1 = plotdata output is suitable as input to the graphical 393C package GRAZIL (based on GKS), which has been developed 394C at ZIB. 395C 396C 397C* Optional INTEGER input/output in IWK: 398C ======================================= 399C 400C Pos. Name Meaning 401C 402C 1 NITER IN/OUT Number of Newton-iterations 403C 2 reserved 404C 3 NCORR IN/OUT Number of corrector steps 405C 4 NFCN IN/OUT Number of FCN-evaluations 406C 5 NJAC IN/OUT Number of Jacobian generations or 407C JAC-calls 408C 6 reserved 409C 7 reserved 410C 8 NFCNJ IN/OUT Number of FCN-evaluations for Jacobian 411C approximation 412C 9 NREJR1 IN/OUT Number of rejected Newton iteration steps 413C done with a rank-1 approximated Jacobian 414C 10..11 Reserved 415C 12 IDCODE IN/OUT Output: The 8 decimal digits program identi- 416C fication number ppppvvvv, consisting of the 417C program code pppp and the version code vvvv. 418C Input: If containing a negative number, 419C it will only be overwritten by the identi- 420C fication number, immediately followed by 421C a return to the calling program. 422C 13..15 Reserved 423C 16 NIWKFR OUT First element of IWK which is free to be used 424C as workspace between Newton iteration steps. 425C For standard linear solver: N+53 426C 17 NRWKFR OUT First element of RWK which is free to be used 427C as workspace between Newton iteration steps. 428C For standard linear solver and numerically 429C approximated Jacobian computed by the 430C expression: (N+9+NBROY)*N+62 431C If the Jacobian is computed by a user routine 432C JAC, subtract N in this expression. 433C 18 LIWKA OUT Length of IWK currently required 434C 19 LRWKA OUT Length of RWK currently required 435C 20..22 Reserved 436C 23 IFAIL OUT Set in case of failure of N2FACT (IERR=80), 437C N2SOLV (IERR=81), FCN (IERR=82) or JAC(IERR=83) 438C to the nonzero IFAIL value returned by the 439C routine indicating the failure . 440C 24..30 Reserved 441C 31 NITMAX IN Maximum number of permitted iteration 442C steps (default: 50) 443C 32 IRANK IN Initial rank of the Jacobian 444C (at the iteration starting point) 445C =0 : full rank N 446C =k with min_rank <= k < N : 447C deficient rank assumed, 448C where min_rank = max (1,N-max(N/10,10)) 449C 33 NEW IN/OUT Count of consecutive rank-1 updates 450C 34..35 Reserved 451C 36 NBROY IN Maximum number of possible consecutive 452C iterative Broyden steps. The total real 453C workspace needed (RWK) depends on this value 454C (see LRWK above). 455C Default is N (see parameter N), 456C if MSTOR=0 (=IOPT(4)), 457C and ML+MU+1 (=IOPT(6)+IOPT(7)+1), if MSTOR=1 458C (but minimum is always 10) - 459C provided that Broyden is allowed. 460C If Broyden is inhibited, NBROY is always set to 461C zero. 462C 37..50 Reserved 463C 464C* Optional REAL input/output in RWK: 465C ==================================== 466C 467C Pos. Name Meaning 468C 469C 1..16 Reserved 470C 17 CONV OUT The achieved relative accuracy after the 471C current step 472C 18 SUMX OUT Natural level (not Normx of printouts) 473C of the current iterate, i.e. Sum(DX(i)**2), 474C where DX = scaled Newton correction. 475C 19 DLEVF OUT Standard level (not Normf of printouts) 476C of the current iterate, i.e. Norm2(F(X)), 477C where F = nonlinear problem function. 478C 20 FCBND IN Bounded damping strategy restriction factor 479C (Default is 10) 480C 21 FCSTRT IN Damping factor for first Newton iteration - 481C overrides option NONLIN, if set (see note 5) 482C 22 FCMIN IN Minimal allowed damping factor (see note 5) 483C 23 SIGMA IN Broyden-approximation decision parameter 484C Required choice: SIGMA.GE.1. Increasing this 485C parameter make it less probable that the algo- 486C rith performs rank-1 updates. 487C Rank-1 updates are inhibited, if 488C SIGMA.GT.1/FCMIN is set. (see note 5) 489C 24 SIGMA2 IN Decision parameter about increasing damping 490C factor to corrector if predictor is small. 491C Required choice: SIGMA2.GE.1. Increasing this 492C parameter make it less probable that the algo- 493C rith performs rank-1 updates. 494C 25 COND IN Maximum permitted subcondition for rank- 495C decision by linear solver. 496C (Default: 1/epmach, epmach: relative 497C machine precision) 498C 26 AJDEL IN Jacobian approximation without feedback: 499C Relative pertubation for components 500C (Default: sqrt(epmach*10), epmach: relative 501C machine precision) 502C 27 AJMIN IN Jacobian approximation without feedback: 503C Threshold value (Default: 0.0d0) 504C The absolute pertubation for component k is 505C computed by 506C DELX := AJDEL*max(abs(Xk),AJMIN) 507C 28 ETADIF IN Jacobian approximation with feedback: 508C Target value for relative pertubation ETA of X 509C (Default: 1.0d-6) 510C 29 ETAINI IN Jacobian approximation with feedback: 511C Initial value for denominator differences 512C (Default: 1.0d-6) 513C 30..50 Reserved 514C 515C Note 5: 516C The default values of the internal parameters may be obtained 517C from the monitor output with at least IOPT field MPRMON set to 2 518C and by initializing the corresponding RWK-fields to zero. 519C 520C* Error and warning messages: 521C =========================== 522C 523C 1 Termination at stationary point (rank deficient Jacobian 524C and termination criterion fulfilled) 525C 2 Termination after NITMAX iterations ( as indicated by 526C input parameter NITMAX=IWK(31) ) 527C 3 Termination, since damping factor became to small and 528C rank of Jacobian matrix is already zero 529C 4 Warning: Superlinear or quadratic convergence slowed down 530C near the solution. 531C Iteration has been stopped therefore with an approximation 532C of the solution not such accurate as requested by RTOL, 533C because possibly the RTOL requirement may be too stringent 534C (i.e. the nonlinear problem is ill-conditioned) 535C 5 Warning: Iteration stopped with termination criterion 536C (using RTOL as requested precision) satisfied, but no 537C superlinear or quadratic convergence has been indicated yet. 538C Therefore, possibly the error estimate for the solution may 539C not match good enough the really achieved accuracy. 540C 10 Integer or real workspace too small 541C 20 Bad input to dimensional parameter N 542C 21 Nonpositive value for RTOL supplied 543C 22 Negative scaling value via vector XSCAL supplied 544C 30 One or more fields specified in IOPT are invalid 545C (for more information, see error-printout) 546C 80 Error signalled by linear solver routine N2FACT, 547C for more detailed information see IFAIL-value 548C stored to IWK(23) 549C (not used by standard routine N2FACT) 550C 81 Error signalled by linear solver routine N2SOLV, 551C for more detailed information see IFAIL-value 552C stored to IWK(23) 553C (not used by standard routine N2SOLV) 554C 82 Error signalled by user routine FCN (Nonzero value 555C returned via IFAIL-flag; stored to IWK(23) ) 556C 83 Error signalled by user routine JAC (Nonzero value 557C returned via IFAIL-flag; stored to IWK(23) ) 558C 559C Note 6 : in case of failure: 560C - use non-standard options 561C - use another initial guess 562C - or reformulate model 563C - or apply continuation techniques 564C 565C* Machine dependent constants used: 566C ================================= 567C 568C DOUBLE PRECISION EPMACH in NLEQ2, N2PCHK, N2INT 569C DOUBLE PRECISION GREAT in N2PCHK 570C DOUBLE PRECISION SMALL in N2PCHK, N2INT, N2SCAL 571C 572C* Subroutines called: N2PCHK, N2INT, D1MACH 573C 574C ------------------------------------------------------------ 575C* End Prologue 576C 577C* Summary of changes: 578C =================== 579C 580C 2.2.1 91, June 3 Time monitor included 581C 2.2.2 91, June 3 Bounded damping strategy implemented 582C 2.2.3 91, July 26 AJDEL, AJMIN as RWK-options for JACGEN.EQ.2, 583C ETADIF, ETAINI as RWK-opts. for JACGEN.EQ.3 584C FCN-count changed for anal. Jacobian 585C 2.2.4 91, August 16 Convergence order monitor included 586C 2.2.5 91, August 19 Standard Broyden updates replaced by 587C iterative Broyden 588C 91, Sept. Rank strategy modified 589C DECCON with new fail exit, for the case that 590C the square root of a negative number will 591C appear during pseudo inverse computation. 592C (Occured, although theoretical impossible!) 593C 2.2.6 91, Sept. 17 Damping factor reduction by FCN-fail imple- 594C mented 595C 2.3 91, Dec. 20 New Release for CodeLib 596C 00, July 12 RTOL output-value bug fixed 597C 598C ------------------------------------------------------------ 599C 600C PARAMETER (IRWKI=xx, LRWKI=yy) 601C IRWKI: Start position of internally used RWK part 602C LRWKI: Length of internally used RWK part 603C (current values see parameter statement below) 604C 605C INTEGER L4,L41,L5,L51,L6,L61,L62,L63,L7,L71,L9,L11,L12,L121, 606C L13,L14,L20 607C Starting positions in RWK of formal array parameters of internal 608C routine N1INT (dynamically determined in driver routine NLEQ1, 609C dependent on N and options setting) 610C 611C Further RWK positions (only internally used) 612C 613C Position Name Meaning 614C 615C IRWKI FCKEEP Damping factor of previous successfull iter. 616C IRWKI+1 FCA Previous damping factor 617C IRWKI+2 FCPRI A priori estimate of damping factor 618C IRWKI+3 DMYCOR Number My of latest corrector damping factor 619C (kept for use in rank-1 decision criterium) 620C IRWKI+(4..LRWKI-1) Free 621C 622C Internal arrays stored in RWK (see routine N2INT for descriptions) 623C 624C Position Array Type Remarks 625C 626C L4 QA(N,N) Perm Arrays QA and DXSAVE need never to 627C L4 DXSAVE(N,NBROY) be kept the same time and therefore 628C Perm may be stored to the same workspace 629C part 630C L41 A(N,N) Perm 631C L5 DX(N) Perm 632C L51 DXQ(N) Perm 633C L6 XA(N) Perm 634C L61 F(N) Perm 635C L62 FW(N) Perm 636C L63 XWA(N) Perm 637C L7 FA(N) Perm 638C L71 ETA(N) Perm Only used for JACGEN=IOPT(3)=3 639C L8 Perm Start position of array workspace 640C needed for linear solver 641C L9 XW(N) Temp 642C L11 DXQA(N) Temp 643C L111 QU(N) Temp 644C L12 T1(N) Temp 645C L121 T2(N) Temp 646C L13 T3(N) Temp Not yet used or even reserved 647C (for future band mode implementat.) 648C 649C 650 EXTERNAL D1MACH,N2INT 651 INTRINSIC DBLE 652 INTEGER IRWKI, LRWKI 653 PARAMETER (IRWKI=51, LRWKI=10) 654 DOUBLE PRECISION ONE 655 PARAMETER (ONE=1.0D0) 656 DOUBLE PRECISION TEN 657 PARAMETER (TEN=1.0D1) 658 DOUBLE PRECISION ZERO 659 PARAMETER (ZERO=0.0D0) 660 INTEGER IRANK,NITMAX,LUERR,LUMON,LUSOL,MPRERR,MPRMON, 661 $MPRSOL,M1,M2,NRWKFR,NRFRIN,NRW,NIWKFR,NIFRIN,NIW,NONLIN,JACGEN 662 INTEGER L4,L41,L5,L51,L6,L61,L62,L63,L7,L71,L8,L9,L11,L12,L121, 663 $ L13,L20 664 DOUBLE PRECISION COND,D1MACH,FC,FCMIN,PERCI,PERCR 665 DOUBLE PRECISION EPMACH 666 LOGICAL QINIMO,QRANK1,QFCSTR,QSUCC,QBDAMP 667 CHARACTER CHGDAT*20, PRODCT*8 668C Which version ? 669 LOGICAL QVCHK 670 INTEGER IVER 671 PARAMETER( IVER=21122301 ) 672C 673C Version: 2.3 Latest change: 674C ----------------------------------------- 675C 676 DATA CHGDAT /'July 12, 2000 '/ 677 DATA PRODCT /'NLEQ2 '/ 678C* Begin 679 EPMACH = D1MACH(3) 680 IERR = 0 681 QVCHK = IWK(12).LT.0 682 IWK(12) = IVER 683 IF (QVCHK) RETURN 684C Print error messages? 685 MPRERR = IOPT(11) 686 LUERR = IOPT(12) 687 IF (LUERR .EQ. 0) THEN 688 LUERR = 6 689 IOPT(12)=LUERR 690 ENDIF 691C Print iteration monitor? 692 MPRMON = IOPT(13) 693 LUMON = IOPT(14) 694 IF (LUMON .LE. 0 .OR. LUMON .GT. 99) THEN 695 LUMON = 6 696 IOPT(14)=LUMON 697 ENDIF 698C Print intermediate solutions? 699 MPRSOL = IOPT(15) 700 LUSOL = IOPT(16) 701 IF (LUSOL .EQ. 0) THEN 702 LUSOL = 6 703 IOPT(16)=LUSOL 704 ENDIF 705C Print time summary statistics? 706 MPRTIM = IOPT(19) 707 LUTIM = IOPT(20) 708 IF (LUTIM .EQ. 0) THEN 709 LUTIM = 6 710 IOPT(20)=LUTIM 711 ENDIF 712 QSUCC = IOPT(1).EQ.1 713 QINIMO = MPRMON.GE.1.AND..NOT.QSUCC 714C Print NLEQ2 heading lines 715 IF(QINIMO)THEN 71610000 FORMAT(' N L E Q 2 ***** V e r s i o n ', 717 $ '2 . 3 ***',//,1X,'Newton-Method ', 718 $ 'for the solution of nonlinear systems',//) 719 WRITE(LUMON,10000) 720 ENDIF 721C Check input parameters and options 722 CALL N2PCHK(N,X,XSCAL,RTOL,IOPT,IERR,LIWK,IWK,LRWK,RWK) 723C Exit, if any parameter error was detected till here 724 IF (IERR.NE.0) RETURN 725 M1=N 726 M2=N 727 JACGEN=IOPT(3) 728 IF (JACGEN.EQ.0) JACGEN=2 729 IOPT(3)=JACGEN 730 QRANK1=IOPT(32).EQ.1 731 IF (QRANK1) THEN 732 NBROY=IWK(36) 733 IF (NBROY.EQ.0) NBROY=MAX(M2,10) 734 IWK(36)=NBROY 735 ELSE 736 NBROY=0 737 ENDIF 738C WorkSpace: RWK 739 L4=IRWKI+LRWKI 740 L41=L4+NBROY*N 741 L5=L41+M1*N 742 L51=L5+N 743 L6=L51+N 744 L61=L6+N 745 L62=L61+N 746 L63=L62+N 747 L7=L63+N 748 L71=L7+N 749 IF (JACGEN.NE.3) THEN 750 L8=L71 751 ELSE 752 L8=L71+N 753 ENDIF 754 NRWKFR = L8 755 L9=LRWK+1-N 756 L11=L9-N 757 L111=L11-N 758 L12=L111-N 759 L121=L12-N 760C L13 : Work array T3, for future band mode implementation 761 L13=L121 762 NRW=NRWKFR+LRWK-L13+1 763C End WorkSpace at NRW 764C WorkSpace: IWK 765 L20=51 766 NIWKFR = L20 767 NIW = NIWKFR-1 768 NIFRIN = NIWKFR 769 NRFRIN = NRWKFR 770 LIWL=N+2 771 LRWL=2*N+1 772 IF (QRANK1) THEN 773 NRWKFR=NRWKFR+LRWL 774 NIWKFR=NIWKFR+LIWL 775 ENDIF 776 NRW = NRW+LRWL 777 NIW = NIW+LIWL 778C End WorkSpace at NIW 779 IWK(16) = NIWKFR 780 IWK(17) = NRWKFR 781C 782 IF(NRW.GT.LRWK.OR.NIW.GT.LIWK)THEN 783 IERR=10 784 ELSE 785 IF(QINIMO)THEN 786 PERCR = DBLE(NRW)/DBLE(LRWK)*100.0D0 787 PERCI = DBLE(NIW)/DBLE(LIWK)*100.0D0 788C Print statistics concerning workspace usage 78910050 FORMAT(' Real Workspace declared as ',I9, 790 $ ' is used up to ',I9,' (',F5.1,' percent)',//, 791 $ ' Integer Workspace declared as ',I9, 792 $ ' is used up to ',I9,' (',F5.1,' percent)',//) 793 WRITE(LUMON,10050)LRWK,NRW,PERCR,LIWK,NIW,PERCI 794 ENDIF 795 IF(QINIMO)THEN 79610051 FORMAT(/,' N =',I4,//,' Prescribed relative ', 797 $ 'precision',D10.2,/) 798 WRITE(LUMON,10051)N,RTOL 79910052 FORMAT(' The Jacobian is supplied by ',A) 800 IF (JACGEN.EQ.1) THEN 801 WRITE(LUMON,10052) 'a user subroutine' 802 ELSE IF (JACGEN.EQ.2) THEN 803 WRITE(LUMON,10052) 804 $ 'numerical differentiation (without feedback strategy)' 805 ELSE IF (JACGEN.EQ.3) THEN 806 WRITE(LUMON,10052) 807 $ 'numerical differentiation (feedback strategy included)' 808 ENDIF 80910057 FORMAT(' Automatic row scaling of the Jacobian is ',A,/) 810 IF (IOPT(35).EQ.1) THEN 811 WRITE(LUMON,10057) 'inhibited' 812 ELSE 813 WRITE(LUMON,10057) 'allowed' 814 ENDIF 815 ENDIF 816 NONLIN=IOPT(31) 817 IF (IOPT(38).EQ.0) QBDAMP = NONLIN.EQ.4 818 IF (IOPT(38).EQ.1) QBDAMP = .TRUE. 819 IF (IOPT(38).EQ.2) QBDAMP = .FALSE. 820 IF (QBDAMP) THEN 821 IF (RWK(20).LT.ONE) RWK(20) = TEN 822 ENDIF 82310064 FORMAT(' Rank-1 updates are ',A) 824 IF (QINIMO) THEN 825 IF (QRANK1) THEN 826 WRITE(LUMON,10064) 'allowed' 827 ELSE 828 WRITE(LUMON,10064) 'inhibited' 829 ENDIF 83010065 FORMAT(' Problem is specified as being ',A) 831 IF (NONLIN.EQ.1) THEN 832 WRITE(LUMON,10065) 'linear' 833 ELSE IF (NONLIN.EQ.2) THEN 834 WRITE(LUMON,10065) 'mildly nonlinear' 835 ELSE IF (NONLIN.EQ.3) THEN 836 WRITE(LUMON,10065) 'highly nonlinear' 837 ELSE IF (NONLIN.EQ.4) THEN 838 WRITE(LUMON,10065) 'extremely nonlinear' 839 ENDIF 84010066 FORMAT(' Bounded damping strategy is ',A,:,/, 841 $ ' Bounding factor is ',D10.3) 842 IF (QBDAMP) THEN 843 WRITE(LUMON,10066) 'active', RWK(20) 844 ELSE 845 WRITE(LUMON,10066) 'off' 846 ENDIF 847 ENDIF 848C Maximum permitted number of iteration steps 849 NITMAX=IWK(31) 850 IF (NITMAX.LE.0) NITMAX=50 851 IWK(31)=NITMAX 85210068 FORMAT(' Maximum permitted number of iteration steps : ', 853 $ I6) 854 IF (QINIMO) WRITE(LUMON,10068) NITMAX 855C Initial damping factor for highly nonlinear problems 856 QFCSTR=RWK(21).GT.ZERO 857 IF (.NOT.QFCSTR) THEN 858 RWK(21)=1.0D-2 859 IF (NONLIN.EQ.4) RWK(21)=1.0D-4 860 ENDIF 861C Minimal permitted damping factor 862 IF (RWK(22).LE.ZERO) THEN 863 RWK(22)=1.0D-4 864 IF (NONLIN.EQ.4) RWK(22)=1.0D-8 865 ENDIF 866 FCMIN=RWK(22) 867C Rank1 decision parameter SIGMA 868 IF (RWK(23).LT.ONE) RWK(23)=3.0D0 869 IF (.NOT.QRANK1) RWK(23)=10.0D0/FCMIN 870C Decision parameter about increasing too small predictor 871C to greater corrector value 872 IF (RWK(24).LT.ONE) RWK(24)=10.0D0/FCMIN 873C Starting value of damping factor (FCMIN.LE.FC.LE.1.0) 874 IF(NONLIN.LE.2.AND..NOT.QFCSTR)THEN 875C for linear or mildly nonlinear problems 876 FC = ONE 877 ELSE 878C for highly or extremely nonlinear problems 879 FC = RWK(21) 880 ENDIF 881 RWK(21)=FC 882C Initial rank 883 IRANK = IWK(32) 884 IF (IRANK.LE.0.OR.IRANK.GT.N) IWK(32) = N 885C Maximum permitted sub condition number of matrix A 886 COND = RWK(25) 887 IF (COND.LT.ONE) COND = ONE/EPMACH 888 RWK(25) = COND 889 IF (MPRMON.GE.2.AND..NOT.QSUCC) THEN 89010069 FORMAT(//,' Internal parameters:',//, 891 $ ' Starting value for damping factor FCSTART = ',D9.2,/, 892 $ ' Minimum allowed damping factor FCMIN = ',D9.2,/, 893 $ ' Rank-1 updates decision parameter SIGMA = ',D9.2,/, 894 $ ' Initial Jacobian pseudo-rank IRANK =',I6,/, 895 $ ' Maximum permitted subcondition COND = ',D9.2) 896 WRITE(LUMON,10069) RWK(21),FCMIN,RWK(23),IWK(32),COND 897 ENDIF 898C Store lengths of currently required workspaces 899 IWK(18) = NIFRIN-1 900 IWK(19) = NRFRIN-1 901C 902C Initialize and start time measurements monitor 903C 904 IF ( IOPT(1).EQ.0 .AND. MPRTIM.NE.0 ) THEN 905 CALL MONINI (' NLEQ2',LUTIM) 906 CALL MONDEF (0,'NLEQ2') 907 CALL MONDEF (1,'FCN') 908 CALL MONDEF (2,'Jacobi') 909 CALL MONDEF (3,'Lin-Fact') 910 CALL MONDEF (4,'Lin-Sol') 911 CALL MONDEF (5,'Output') 912 CALL MONSRT () 913 ENDIF 914C 915C 916 IERR=-1 917C If IERR is unmodified on exit, successive steps are required 918C to complete the Newton iteration 919 IF (NBROY.EQ.0) NBROY=1 920 CALL N2INT(N,FCN,JAC,X,XSCAL,RTOL,NITMAX,NONLIN,IWK(32),IOPT, 921 $ IERR,LRWK,RWK,NRFRIN,LRWL,LIWK,IWK,NIFRIN,LIWL,M1,M2,NBROY, 922 $ RWK(L4),RWK(L41),RWK(L4),RWK(L5),RWK(L51),RWK(L6),RWK(L63), 923 $ RWK(L61), 924 $ RWK(L7),RWK(L71),RWK(L9),RWK(L62),RWK(L11),RWK(L111),RWK(L12), 925 $ RWK(L121),RWK(L13),RWK(21),RWK(22),RWK(23),RWK(24),RWK(IRWKI+1), 926 $ RWK(IRWKI),RWK(IRWKI+2),COND,RWK(IRWKI+3),RWK(17),RWK(18), 927 $ RWK(19),MPRERR,MPRMON,MPRSOL,LUERR,LUMON,LUSOL,IWK(1),IWK(3), 928 $ IWK(4),IWK(5),IWK(8),IWK(9),IWK(33),QBDAMP) 929C 930 IF (MPRTIM.NE.0.AND.IERR.NE.-1.AND.IERR.NE.10) CALL MONEND 931C 932C Free workspaces, so far not used between steps 933 IWK(16) = NIWKFR 934 IWK(17) = NRWKFR 935 ENDIF 936C Print statistics 937 IF (MPRMON.GE.1.AND.IERR.NE.-1.AND.IERR.NE.10) THEN 93810080 FORMAT(/, ' ****** Statistics * ', A8, ' *******', /, 939 $ ' *** Newton iterations : ', I7,' ***', /, 940 $ ' *** Corrector steps : ', I7,' ***', /, 941 $ ' *** Rejected rk-1 st. : ', I7,' ***', /, 942 $ ' *** Jacobian eval. : ', I7,' ***', /, 943 $ ' *** Function eval. : ', I7,' ***', /, 944 $ ' *** ... for Jacobian : ', I7,' ***', /, 945 $ ' *************************************', /) 946 WRITE (LUMON,10080) PRODCT,IWK(1),IWK(3),IWK(9),IWK(5), 947 $ IWK(4),IWK(8) 948 ENDIF 949C Print workspace requirements, if insufficient 950 IF (IERR.EQ.10) THEN 95110090 FORMAT(///,20('*'),'Workspace Error',20('*')) 952 IF (MPRERR.GE.1) WRITE(LUERR,10090) 953 IF(NRW.GT.LRWK)THEN 95410091 FORMAT(/,' Real Workspace dimensioned as',1X,I9, 955 $ 1X,'must be enlarged at least up to ', 956 $ I9,//) 957 IF (MPRERR.GE.1) WRITE(LUERR,10091)LRWK,NRFRIN-1 958 ENDIF 959 IF(NIW.GT.LIWK)THEN 96010092 FORMAT(/,' Integer Workspace dimensioned as ', 961 $ I9,' must be enlarged at least up ', 962 $ 'to ',I9,//) 963 IF (MPRERR.GE.1) WRITE(LUERR,10092)LIWK,NIFRIN-1 964 ENDIF 965 ENDIF 966C End of subroutine NLEQ2 967 RETURN 968 END 969C 970 SUBROUTINE N2PCHK(N,X,XSCAL,RTOL,IOPT,IERR,LIWK,IWK,LRWK,RWK) 971C* Begin Prologue N2PCHK 972 INTEGER N 973 DOUBLE PRECISION X(N),XSCAL(N) 974 DOUBLE PRECISION RTOL 975 INTEGER IOPT(50) 976 INTEGER IERR 977 INTEGER LIWK 978 INTEGER IWK(LIWK) 979 INTEGER LRWK 980 DOUBLE PRECISION RWK(LRWK) 981C ------------------------------------------------------------ 982C 983C* Summary : 984C 985C N 2 P C H K : Checking of input parameters and options 986C for NLEQ2. 987C 988C* Parameters: 989C =========== 990C 991C See parameter description in driver routine. 992C 993C* Subroutines called: D1MACH 994C 995C* Machine dependent constants used: 996C ================================= 997C 998C EPMACH = relative machine precision 999C GREAT = squareroot of maxreal divided by 10 1000C SMALL = squareroot of "smallest positive machine number 1001C divided by relative machine precision" 1002 DOUBLE PRECISION EPMACH,GREAT,SMALL 1003C 1004C ------------------------------------------------------------ 1005C* End Prologue 1006C 1007 EXTERNAL D1MACH 1008 INTRINSIC DBLE 1009 DOUBLE PRECISION ONE 1010 PARAMETER (ONE=1.0D0) 1011 DOUBLE PRECISION TEN 1012 PARAMETER (TEN=1.0D1) 1013 DOUBLE PRECISION ZERO 1014 PARAMETER (ZERO=0.0D0) 1015C 1016 PARAMETER (NUMOPT=50) 1017 INTEGER IOPTL(NUMOPT),IOPTU(NUMOPT) 1018 DOUBLE PRECISION D1MACH,TOLMIN,TOLMAX,DEFSCL 1019C 1020 DATA IOPTL /0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,1,0,0,0,1, 1021 $ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 1022 $ 0,0,0,0,0,0,0,0,0,0, 1023 $ -9999999,-9999999,-9999999,-9999999,-9999999/ 1024 DATA IOPTU /1,1,3,0,0,0,0,0,1,0,3,99,6,99,3,99,0,0,1,99, 1025 $ 0,0,0,0,0,0,0,0,0,0,4,1,0,0,1, 1026 $ 0,0,2,3,0,0,0,0,0,0, 1027 $ 9999999,9999999,9999999,9999999,9999999/ 1028C 1029 EPMACH = D1MACH(3) 1030 GREAT = DSQRT(D1MACH(2)/TEN) 1031 SMALL = D1MACH(6) 1032 IERR = 0 1033C Print error messages? 1034 MPRERR = IOPT(11) 1035 LUERR = IOPT(12) 1036 IF (LUERR .LE. 0 .OR. LUERR .GT. 99) THEN 1037 LUERR = 6 1038 IOPT(12)=LUERR 1039 ENDIF 1040C 1041C Checking dimensional parameter N 1042 IF ( N.LE.0 ) THEN 1043 IF (MPRERR.GE.1) WRITE(LUERR,10011) N 104410011 FORMAT(/,' Error: Bad input to dimensional parameter N supplied' 1045 $ ,/,8X,'choose N positive, your input is: N = ',I5) 1046 IERR = 20 1047 ENDIF 1048C 1049C Problem type specification by user 1050 NONLIN=IOPT(31) 1051 IF (NONLIN.EQ.0) NONLIN=3 1052 IOPT(31)=NONLIN 1053C 1054C Checking and conditional adaption of the user-prescribed RTOL 1055 IF (RTOL.LE.ZERO) THEN 1056 IF (MPRERR.GE.1) 1057 $ WRITE(LUERR,'(/,A)') ' Error: Nonpositive RTOL supplied' 1058 IERR = 21 1059 ELSE 1060 TOLMIN = EPMACH*TEN*DBLE(N) 1061 IF(RTOL.LT.TOLMIN) THEN 1062 RTOL = TOLMIN 1063 IF (MPRERR.GE.2) 1064 $ WRITE(LUERR,10012) 'increased ','smallest',RTOL 1065 ENDIF 1066 TOLMAX = 1.0D-1 1067 IF(RTOL.GT.TOLMAX) THEN 1068 RTOL = TOLMAX 1069 IF (MPRERR.GE.2) 1070 $ WRITE(LUERR,10012) 'decreased ','largest',RTOL 1071 ENDIF 107210012 FORMAT(/,' Warning: User prescribed RTOL ',A,'to ', 1073 $ 'reasonable ',A,' value RTOL = ',D11.2) 1074 ENDIF 1075C 1076C Test user prescribed accuracy and scaling on proper values 1077 IF (N.LE.0) RETURN 1078 IF (NONLIN.GE.3) THEN 1079 DEFSCL = RTOL 1080 ELSE 1081 DEFSCL = ONE 1082 ENDIF 1083 DO 10 I=1,N 1084 IF (XSCAL(I).LT.ZERO) THEN 1085 IF (MPRERR.GE.1) THEN 1086 WRITE(LUERR,10013) I 108710013 FORMAT(/,' Error: Negative value in XSCAL(',I5,') supplied') 1088 ENDIF 1089 IERR = 22 1090 ENDIF 1091 IF (XSCAL(I).EQ.ZERO) XSCAL(I) = DEFSCL 1092 IF ( XSCAL(I).GT.ZERO .AND. XSCAL(I).LT.SMALL ) THEN 1093 IF (MPRERR.GE.2) THEN 1094 WRITE(LUERR,10014) I,XSCAL(I),SMALL 109510014 FORMAT(/,' Warning: XSCAL(',I5,') = ',D9.2,' too small, ', 1096 $ 'increased to',D9.2) 1097 ENDIF 1098 XSCAL(I) = SMALL 1099 ENDIF 1100 IF (XSCAL(I).GT.GREAT) THEN 1101 IF (MPRERR.GE.2) THEN 1102 WRITE(LUERR,10015) I,XSCAL(I),GREAT 110310015 FORMAT(/,' Warning: XSCAL(',I5,') = ',D9.2,' too big, ', 1104 $ 'decreased to',D9.2) 1105 ENDIF 1106 XSCAL(I) = GREAT 1107 ENDIF 110810 CONTINUE 1109C Checks options 1110 DO 20 I=1,30 1111 IF (IOPT(I).LT.IOPTL(I) .OR. IOPT(I).GT.IOPTU(I)) THEN 1112 IERR=30 1113 IF (MPRERR.GE.1) THEN 1114 WRITE(LUERR,20001) I,IOPT(I),IOPTL(I),IOPTU(I) 111520001 FORMAT(' Invalid option specified: IOPT(',I2,')=',I12,';', 1116 $ /,3X,'range of permitted values is ',I8,' to ',I8) 1117 ENDIF 1118 ENDIF 111920 CONTINUE 1120C End of subroutine N2PCHK 1121 RETURN 1122 END 1123C 1124 SUBROUTINE N2INT(N,FCN,JAC,X,XSCAL,RTOL,NITMAX,NONLIN,IRANK,IOPT, 1125 $IERR,LRWK,RWK,NRWKFR,LRWL,LIWK,IWK,NIWKFR,LIWL,M1,M2,NBROY, 1126 $QA,A,DXSAVE,DX,DXQ,XA,XWA,F,FA,ETA,XW,FW,DXQA,QU,T1,T2,T3,FC, 1127 $FCMIN,SIGMA,SIGMA2,FCA,FCKEEP,FCPRI,COND,DMYCOR,CONV,SUMX,DLEVF, 1128 $MPRERR,MPRMON,MPRSOL,LUERR,LUMON,LUSOL,NITER,NCORR,NFCN,NJAC, 1129 $NFCNJ,NREJR1,NEW,QBDAMP) 1130C* Begin Prologue N2INT 1131 INTEGER N 1132 EXTERNAL FCN,JAC 1133 DOUBLE PRECISION X(N),XSCAL(N) 1134 DOUBLE PRECISION RTOL 1135 INTEGER NITMAX,NONLIN,IRANK 1136 INTEGER IOPT(50) 1137 INTEGER IERR 1138 INTEGER LRWK 1139 DOUBLE PRECISION RWK(LRWK) 1140 INTEGER NRWKFR,LRWL,LIWK 1141 INTEGER IWK(LIWK) 1142 INTEGER NIWKFR,LIWL,M1,M2,NBROY 1143 DOUBLE PRECISION QA(M2,N),A(M1,N),DXSAVE(N,NBROY) 1144 DOUBLE PRECISION DX(N),DXQ(N),XA(N),XWA(N),F(N),FA(N),ETA(N) 1145 DOUBLE PRECISION XW(N),FW(N),DXQA(N),QU(N),T1(N),T2(N),T3(N) 1146 DOUBLE PRECISION FC,FCMIN,SIGMA,SIGMA2,FCA,FCKEEP,COND,CONV,SUMX, 1147 $ DLEVF,FCPRI,DMYCOR 1148 INTEGER MPRERR,MPRMON,MPRSOL,LUERR,LUMON,LUSOL,NITER, 1149 $NCORR,NFCN,NJAC,NFCNJ,NREJR1,NEW 1150 LOGICAL QBDAMP 1151C ------------------------------------------------------------ 1152C 1153C* Summary : 1154C 1155C N 2 I N T : Core routine for NLEQ2 . 1156C Damped Newton-algorithm with rank-strategy for systems of 1157C highly nonlinear equations especially designed for 1158C numerically sensitive problems. 1159C 1160C* Parameters: 1161C =========== 1162C 1163C N,FCN,JAC,X,XSCAL,RTOL 1164C See parameter description in driver routine 1165C 1166C NITMAX Int Maximum number of allowed iterations 1167C NONLIN Int Problem type specification 1168C (see IOPT-field NONLIN) 1169C IRANK Int Initially proposed (in) and final (out) rank 1170C of Jacobian 1171C IOPT Int See parameter description in driver routine 1172C IERR Int See parameter description in driver routine 1173C LRWK Int Length of real workspace 1174C RWK(LRWK) Dble Real workspace array 1175C NRWKFR Int First free position of RWK on exit 1176C LIWK Int Length of integer workspace 1177C IWK(LIWK) Int Integer workspace array 1178C NIWKFR Int First free position of IWK on exit 1179C M1 Int Leading dimension of Jacobian array A 1180C for full case Jacobian: N 1181C (other matrix types are not yet implemented) 1182C M2 Int Leading dimension of Jacobian array QA 1183C for full case Jacobian: N 1184C NBROY Int Maximum number of possible consecutive 1185C iterative Broyden steps. (See IWK(36)) 1186C QA(M2,N) Dble Holds the originally computed Jacobian 1187C or the pseudo inverse in case of rank- 1188C deficiency 1189C A(M1,N) Dble Holds the Jacobian matrix (decomposed form 1190C after call of linear decomposition 1191C routine) 1192C DXSAVE(X,NBROY) 1193C Dble Used to save the quasi Newton corrections of 1194C all previously done consecutive Broyden 1195C steps. 1196C DX(N) Dble Current Newton correction 1197C DXQ(N) Dble Simplified Newton correction J(k-1)*X(k) 1198C XA(N) Dble Previous Newton iterate 1199C XWA(N) Dble Scaling factors used for latest decomposed 1200C Jacobian for column scaling - may differ 1201C from XW, if Broyden updates are performed 1202C F(N) Dble Function (FCN) value of current iterate 1203C FA(N) Dble Function (FCN) value of previous iterate 1204C ETA(N) Dble Jacobian approximation: updated scaled 1205C denominators 1206C XW(N) Dble Scaling factors for iteration vector 1207C FW(N) Dble Scaling factors for rows of the system 1208C DXQA(N) Dble Previous Newton correction 1209C QU(N) Dble Savespace for right hand side belonging 1210C to upper triangular linear system 1211C T1(N) Dble Workspace for linear solvers and internal 1212C subroutines 1213C T2(N) Dble Workspace array for internal subroutines 1214C T3(N) Dble Workspace array for internal subroutines 1215C FC Dble Current Newton iteration damping factor. 1216C FCMIN Dble Minimum permitted damping factor. If 1217C FC becomes smaller than this value, one 1218C of the following may occur: 1219C a. Recomputation of the Jacobian 1220C matrix by means of difference 1221C approximation (instead of Rank1 1222C update), if Rank1 - update 1223C previously was used 1224C b. Rank reduction of Jacobi 1225C matrix , if difference 1226C approximation was used previously 1227C and Rank(A).NE.0 1228C c. Fail exit otherwise 1229C SIGMA Dble Decision parameter for rank1-updates. 1230C SIGMA2 Dble Decision parameter for damping factor 1231C increasing to corrector value 1232C FCA Dble Previous Newton iteration damping factor. 1233C FCKEEP Dble Keeps the damping factor as it is at start 1234C of iteration step. 1235C COND Dble Maximum permitted subcondition for rank- 1236C decision by linear solver. 1237C CONV Dble Scaled maximum norm of the Newton- 1238C correction. Passed to RWK-field on output. 1239C SUMX Dble Square of the natural level (see equal- 1240C named IOPT-output field) 1241C DLEVF Dble Square of the standard level (see equal- 1242C named IOPT-output field) 1243C MPRERR,MPRMON,MPRSOL,LUERR,LUMON,LUSOL : 1244C See description of equal named IOPT-fields 1245C in the driver subroutine 1246C NITER,NCORR,NFCN,NJAC,NFCNJ,NREJR1,NEW : 1247C See description of equal named IWK-fields 1248C in the driver subroutine 1249C QBDAMP Logic Flag, that indicates, whether bounded damping 1250C strategy is active: 1251C .true. = bounded damping strategy is active 1252C .false. = normal damping strategy is active 1253C 1254C 1255C* Internal double variables 1256C ========================= 1257C 1258C AJDEL See RWK(26) (num. diff. without feedback) 1259C AJMIN See RWK(27) (num. diff. without feedback) 1260C COND1 Gets the subcondition of the linear system 1261C as estimated by the linear solver (N2FACT) 1262C CONVA Holds the previous value of CONV . 1263C DEL Gets the projection defect in case of rank- 1264C deficiency. 1265C DMUE Temporary value used during computation of damping 1266C factors predictor. 1267C EPDIFF sqrt(10*epmach) (num. diff. with feedback) 1268C ETADIF See description of RWK(28) (num. diff. with feedback) 1269C ETAINI Initial value for all ETA-components (num. diff. fb.) 1270C ETAMAX Maximum allowed pertubation (num. diff. with feedback) 1271C ETAMIN Minimum allowed pertubation (num. diff. with feedback) 1272C FCDNM Used to compute the denominator of the damping 1273C factor FC during computation of it's predictor, 1274C corrector and aposteriori estimate (in the case of 1275C performing a Rank1 update) . 1276C FCK2 Aposteriori estimate of FC. 1277C FCMIN2 FCMIN**2 . Used for FC-predictor computation. 1278C FCMINH DSQRT(FCMIN). 1279C Used in rank decision logical expression. 1280C FCNUMP Gets the numerator of the predictor formula for FC. 1281C FCNMP2 Temporary used for predictor numerator computation. 1282C FCNUMK Gets the numerator of the corrector computation 1283C of FC . 1284C SENS1 Gets the sensitivity of the Jacobian as 1285C estimated by the linear solver N2FACT. 1286C SUMXA Natural level of the previous iterate. 1287C TH Temporary variable used during corrector- and 1288C aposteriori computations of FC. 1289C 1290C* Internal integer variables 1291C ========================== 1292C 1293C IFAIL Gets the return value from subroutines called from 1294C N2INT (N2FACT, N2SOLV, FCN, JAC) 1295C ISCAL Holds the scaling option from the IOPT-field ISCAL 1296C MODE Matrix storage mode (see IOPT-field MODE) 1297C NRED Count of successive corrector steps 1298C NILUSE Gets the amount of IWK used by the linear solver 1299C NRLUSE Gets the amount of RWK used by the linear solver 1300C NIWLA Index of first element of IWK provided to the 1301C linear solver 1302C NRWLA Index of first element of RWK provided to the 1303C linear solver 1304C LIWL Holds the maximum amount of integer workspace 1305C available to the linear solver 1306C LRWL Holds the maximum amount of real workspace 1307C available to the linear solver 1308C 1309C* Internal logical variables 1310C ========================== 1311C 1312C QGENJ Jacobian updating technique flag: 1313C =.TRUE. : Call of analytical subroutine JAC or 1314C numerical differentiation 1315C =.FALSE. : rank1- (Broyden-) update 1316C QINISC Iterate initial-scaling flag: 1317C =.TRUE. : at first call of N2SCAL 1318C =.FALSE. : at successive calls of N2SCAL 1319C QSUCC See description of IOPT-field QSUCC. 1320C QJCRFR Jacobian refresh flag: 1321C set to .TRUE. if damping factor gets too small 1322C and Jacobian was computed by rank1-update. 1323C Indicates, that the Jacobian needs to be recomputed 1324C by subroutine JAC or numerical differentation. 1325C QLINIT Initialization state of linear solver workspace: 1326C =.FALSE. : Not yet initialized 1327C =.TRUE. : Initialized - N2FACT has been called at 1328C least one time. 1329C QREPET Operation mode flag for linear solver: 1330C =.FALSE. : Normal operation (full rank matrix) 1331C =.TRUE. : Special operation in rank deficient case: 1332C Compute rank-deficient pseudo-inverse, 1333C partial recomputation when solving the 1334C linear system again prescribing a lower 1335C rank as before. 1336C QNEXT Set to .TRUE. to indicate success of the current 1337C Newton-step, i.e. : sucessfull monotonicity-test. 1338C 1339C QREDU Set to .TRUE. to indicate that rank-reduction (or 1340C refreshment of the Jacobian) is needed - if the 1341C computed damping factor gets too small. 1342C QSCALE Holds the value of .NOT.QNSCAL. See description 1343C of IOPT-field QNSCAL. 1344C 1345C* Subroutines called: 1346C =================== 1347C 1348C N2FACT, N2SOLV, N2JAC, N2JCF, N2LVLS, N2PRJN, 1349C N2SCRF, N2SOUT, N2PRV1, N2PRV2, N2SCAL, 1350C MONON, MONOFF 1351C 1352C* Functions called: 1353C ================= 1354C 1355C D1MACH, WNORM 1356C 1357C* Machine constants used 1358C ====================== 1359C 1360 DOUBLE PRECISION EPMACH,SMALL 1361C 1362C ------------------------------------------------------------ 1363C* End Prologue 1364 EXTERNAL N2FACT, N2SOLV, N2JAC, N2JCF, N2LVLS, N2PRJN, 1365 $ N2SCRF, N2SOUT, N2PRV1, N2PRV2, N2SCAL, 1366 $ MONON, MONOFF, D1MACH, WNORM 1367 INTRINSIC DSQRT,DMIN1 1368 DOUBLE PRECISION ONE 1369 PARAMETER (ONE=1.0D0) 1370 DOUBLE PRECISION ZERO 1371 PARAMETER (ZERO=0.0D0) 1372 DOUBLE PRECISION HALF 1373 PARAMETER (HALF=0.5D0) 1374 DOUBLE PRECISION TEN 1375 PARAMETER (TEN=10.0D0) 1376 INTEGER IFAIL,ILOOP,ISCAL,K,MODE,NRED,NILUSE,NRLUSE,NIWLA, 1377 $NRWLA,L1,JACGEN,MINRNK 1378 DOUBLE PRECISION AJDEL,AJMIN,ALFA1,ALFA2,ALFA,BETA, 1379 $COND1,CONVA,DLEVXA,DEL,DMYPRI,D1MACH,DXANRM,DXNRM,WNORM,EPDIFF, 1380 $ETAMIN,ETAMAX,ETAINI,ETADIF,FCDNM,FCBND,FCBH,FCK2,FCMIN2,FCMINH, 1381 $FCNUMP,FCCOR,FCNMP2,FCNUMK,FCREDU,DLEVFN,SENS1,SUMXA,SUM1,SUM2, 1382 $S1,TH,RSMALL,APREC 1383 LOGICAL QGENJ,QINISC,QSUCC,QJCRFR,QLINIT,QNEXT,QRANK1,QREPET, 1384 $ QREDU,QREP,QSCALE,QMIXIO 1385CWEI 1386 INTRINSIC DLOG 1387 DOUBLE PRECISION CLIN0,CLIN1,CALPHA,CALPHK,ALPHAE,ALPHAK,ALPHAA, 1388 $ SUMXA0,SUMXA1,SUMXA2,SUMXTE,FCMON,DLOG 1389 INTEGER ICONV, IORMON 1390 LOGICAL QMSTOP 1391 SAVE CLIN0,CLIN1,CALPHA,ALPHAE,ALPHAK,ALPHAA,SUMXA0,SUMXA1,SUMXA2, 1392 $ ICONV,QMSTOP 1393C 1394 EPMACH = D1MACH(3) 1395 SMALL = D1MACH(6) 1396C* Begin 1397C ---------------------------------------------------------- 1398C 1 Initialization 1399C ---------------------------------------------------------- 1400C 1.1 Control-flags and -integers 1401 QSUCC = IOPT(1).EQ.1 1402 QSCALE = .NOT. IOPT(35).EQ.1 1403 QRANK1 = IOPT(32).EQ.1 1404 IORMON = IOPT(39) 1405 IF (IORMON.EQ.0) IORMON=2 1406 ISCAL = IOPT(9) 1407 MODE = IOPT(2) 1408 JACGEN = IOPT(3) 1409 QMIXIO = LUMON.EQ.LUSOL .AND. MPRMON.NE.0 .AND. MPRSOL.NE.0 1410 MPRTIM = IOPT(19) 1411C ---------------------------------------------------------- 1412C 1.2 Derivated dimensional parameters 1413 MINRNK = MAX0(1,N-MAX0(INT(FLOAT(N)/10.0),10)) 1414C ---------------------------------------------------------- 1415C 1.3 Derivated internal parameters 1416 FCMIN2 = FCMIN*FCMIN 1417 FCMINH = DSQRT(FCMIN) 1418 TOLMIN = DSQRT(TEN*EPMACH) 1419 RSMALL = DSQRT(TEN*RTOL) 1420C ---------------------------------------------------------- 1421C 1.4 Adaption of input parameters, if necessary 1422 IF(FC.LT.FCMIN) FC = FCMIN 1423 IF(FC.GT.ONE) FC = ONE 1424C ---------------------------------------------------------- 1425C 1.5 Initial preparations 1426 QJCRFR = .FALSE. 1427 QLINIT = .FALSE. 1428 QREPET = .FALSE. 1429 IFAIL = 0 1430 FCBND = ZERO 1431 IF (QBDAMP) FCBND = RWK(20) 1432C ---------------------------------------------------------- 1433C 1.5.1 Numerical differentation related initializations 1434 IF (JACGEN.EQ.2) THEN 1435 AJDEL = RWK(26) 1436 IF (AJDEL.LE.SMALL) AJDEL = DSQRT(EPMACH*TEN) 1437 AJMIN = RWK(27) 1438 ELSE IF (JACGEN.EQ.3) THEN 1439 ETADIF = RWK(28) 1440 IF (ETADIF .LE. SMALL) ETADIF = 1.0D-6 1441 ETAINI = RWK(29) 1442 IF (ETAINI .LE. SMALL) ETAINI = 1.0D-6 1443 EPDIFF = DSQRT(EPMACH*TEN) 1444 ETAMAX = DSQRT(EPDIFF) 1445 ETAMIN = EPDIFF*ETAMAX 1446 ENDIF 1447C ---------------------------------------------------------- 1448C 1.5.2 Miscellaneous preparations of first iteration step 1449 IF (.NOT.QSUCC) THEN 1450 NITER = 0 1451 NCORR = 0 1452 NREJR1 = 0 1453 NFCN = 0 1454 NJAC = 0 1455 NFCNJ = 0 1456 QGENJ = .TRUE. 1457 QINISC = .TRUE. 1458 FCKEEP = FC 1459 FCA = FC 1460 FCPRI = FC 1461 FCK2 = FC 1462 CONV = ZERO 1463 IF (JACGEN.EQ.3) THEN 1464 DO 1520 L1=1,N 1465 ETA(L1)=ETAINI 14661520 CONTINUE 1467 ENDIF 1468 DO 1521 L1=1,N 1469 XA(L1)=X(L1) 14701521 CONTINUE 1471CWEI 1472 ICONV = 0 1473 ALPHAE = ZERO 1474 SUMXA1 = ZERO 1475 SUMXA0 = ZERO 1476 CLIN0 = ZERO 1477 QMSTOP = .FALSE. 1478C ------------------------------------------------------ 1479C 1.6 Print monitor header 1480 IF(MPRMON.GE.2 .AND. .NOT.QMIXIO)THEN 148116003 FORMAT(///,2X,66('*')) 1482 WRITE(LUMON,16003) 148316004 FORMAT(/,8X,'It',7X,'Normf ',10X,'Normx ',8X, 1484 $ 'Damp.Fct.',3X,'New',6X,'Rank',8X,'Cond') 1485 WRITE(LUMON,16004) 1486 ENDIF 1487C -------------------------------------------------------- 1488C 1.7 Startup step 1489C -------------------------------------------------------- 1490C 1.7.1 Computation of the residual vector 1491 IF (MPRTIM.NE.0) CALL MONON(1) 1492 CALL FCN(N,X,F,IFAIL) 1493 IF (MPRTIM.NE.0) CALL MONOFF(1) 1494 NFCN = NFCN+1 1495C Exit, if ... 1496 IF (IFAIL.NE.0) THEN 1497 IERR = 82 1498 GOTO 4299 1499 ENDIF 1500 ELSE 1501 QINISC = .FALSE. 1502 ENDIF 1503C 1504C Main iteration loop 1505C =================== 1506C 1507C Repeat 15082 CONTINUE 1509C -------------------------------------------------------- 1510C 2 Startup of iteration step 1511 IF (.NOT.QJCRFR) THEN 1512C ------------------------------------------------------ 1513C 2.1 Scaling of variables X(N) 1514 CALL N2SCAL(N,X,XA,XSCAL,XW,ISCAL,QINISC,IOPT,LRWK,RWK) 1515 QINISC = .FALSE. 1516 IF(NITER.NE.0)THEN 1517C Preliminary pseudo-rank 1518 IRANKA = IRANK 1519C ---------------------------------------------------- 1520C 2.2 Aposteriori estimate of damping factor 1521 DO 2200 L1=1,N 1522 DXQA(L1)=DXQ(L1) 15232200 CONTINUE 1524 FCNUMP = ZERO 1525 DO 2201 L1=1,N 1526 FCNUMP=FCNUMP+(DX(L1)/XW(L1))**2 15272201 CONTINUE 1528 TH = FC-ONE 1529 FCDNM = ZERO 1530 DO 2202 L1=1,N 1531 FCDNM=FCDNM+((DXQA(L1)+TH*DX(L1))/XW(L1))**2 15322202 CONTINUE 1533C ---------------------------------------------------- 1534C 2.2.2 Decision criterion for Jacobian updating 1535C technique: 1536C QGENJ.EQ..TRUE. numerical differentation, 1537C QGENJ.EQ..FALSE. rank1 updating 1538 QGENJ = .TRUE. 1539 IF (FC.EQ.FCPRI) THEN 1540 QGENJ = FC.LT.ONE.OR.FCA.LT.ONE.OR.DMYCOR.LE.FC*SIGMA 1541 $ .OR. .NOT.QRANK1 .OR. NEW+2.GT.NBROY 1542 FCA = FC 1543 ELSE 1544 DMYCOR = FCA*FCA*HALF*DSQRT(FCNUMP/FCDNM) 1545 IF (NONLIN.LE.3) THEN 1546 FCCOR = DMIN1(ONE,DMYCOR) 1547 ELSE 1548 FCCOR = DMIN1(ONE,HALF*DMYCOR) 1549 ENDIF 1550 FCA = DMAX1(DMIN1(FC,FCCOR),FCMIN) 1551C$Test-begin 1552 IF (MPRMON.GE.5) THEN 1553 WRITE(LUMON,22201) FCCOR, FC, DMYCOR, FCNUMP, 1554 $ FCDNM 155522201 FORMAT (/, ' +++ aposteriori estimate +++', /, 1556 $ ' FCCOR = ', D18.10, ' FC = ', D18.10, /, 1557 $ ' DMYCOR = ', D18.10, ' FCNUMP = ', D18.10, /, 1558 $ ' FCDNM = ', D18.10, /, 1559 $ ' ++++++++++++++++++++++++++++', /) 1560 ENDIF 1561C$Test-end 1562 ENDIF 1563 FCK2 = FCA 1564C ------------------------------------------------------ 1565C 2.2.1 Computation of the numerator of damping 1566C factor predictor 1567 FCNMP2 = ZERO 1568 DO 221 L1=1,N 1569 FCNMP2=FCNMP2+(DXQA(L1)/XW(L1))**2 1570221 CONTINUE 1571 FCNUMP = FCNUMP*FCNMP2 1572 ENDIF 1573 ENDIF 1574 QJCRFR =.FALSE. 1575C -------------------------------------------------------- 1576C 2.3 Jacobian matrix (stored to array A(M1,N)) 1577C -------------------------------------------------------- 1578C 2.3.1 Jacobian generation by routine JAC or 1579C difference approximation (If QGENJ.EQ..TRUE.) 1580C - or - 1581C Rank-1 update of Jacobian (If QGENJ.EQ..FALSE.) 1582 IF(QGENJ)THEN 1583 NEW = 0 1584 IF (JACGEN.EQ.1) THEN 1585 IF (MPRTIM.NE.0) CALL MONON(2) 1586 CALL JAC(N,M1,X,A,IFAIL) 1587 IF (MPRTIM.NE.0) CALL MONOFF(2) 1588 ELSE 1589 IF (MPRTIM.NE.0) CALL MONON(2) 1590 IF (JACGEN.EQ.3) 1591 $ CALL N2JCF(FCN,N,M1,X,F,A,XW,ETA,ETAMIN,ETAMAX, 1592 $ ETADIF,CONV,NFCNJ,T1,IFAIL) 1593 IF (JACGEN.EQ.2) 1594 $ CALL N2JAC(FCN, N, M1, X, F, A, XW, AJDEL, AJMIN, 1595 $ NFCNJ, T1, IFAIL) 1596 IF (MPRTIM.NE.0) CALL MONOFF(2) 1597 ENDIF 1598 NJAC = NJAC + 1 1599C Exit, If ... 1600 IF (JACGEN.EQ.1 .AND. IFAIL.LT.0) THEN 1601 IERR = 83 1602 GOTO 4299 1603 ENDIF 1604 IF (JACGEN.NE.1 .AND. IFAIL.NE.0) THEN 1605 IERR = 82 1606 GOTO 4299 1607 ENDIF 1608 ELSE 1609 NEW = NEW+1 1610 ENDIF 1611 IF ( NEW.EQ.0 ) THEN 1612C ------------------------------------------------------ 1613C 2.3.2 Save scaling values 1614 DO 232 L1=1,N 1615 XWA(L1) = XW(L1) 1616232 CONTINUE 1617C -------------------------------------------------------- 1618C 2.4 Prepare solution of the linear system 1619C -------------------------------------------------------- 1620C 2.4.1 internal column scaling of matrix A 1621 DO 2410 K=1,N 1622 S1 =-XW(K) 1623 DO 2412 L1=1,N 1624 A(L1,K)=A(L1,K)*S1 16252412 CONTINUE 16262410 CONTINUE 1627C ------------------------------------------------------ 1628C 2.4.2 Row scaling of matrix A 1629 IF (QSCALE) THEN 1630 CALL N2SCRF(N,N,A,FW) 1631 ELSE 1632 DO 242 K=1,N 1633 FW(K)=ONE 1634242 CONTINUE 1635 ENDIF 1636 ENDIF 1637C -------------------------------------------------------- 1638C 2.4.3 Save and scale values of F(N) 1639 DO 243 L1=1,N 1640 FA(L1)=F(L1) 1641 T1(L1)=F(L1)*FW(L1) 1642243 CONTINUE 1643 IRANKA = IRANK 1644 IF (NITER.NE.0) IRANK = N 1645 QNEXT = .FALSE. 1646C -------------------------------------------------------- 1647C 3 Central part of iteration step 1648C 1649C Pseudo-rank reduction loop 1650C ========================== 1651C DO (Until) 16523 CONTINUE 1653C -------------------------------------------------------- 1654C 3.1 Solution of the linear system 1655C -------------------------------------------------------- 1656C 3.1.1 Decomposition of (N,N)-matrix A 1657 IF (.NOT.QLINIT) THEN 1658 NIWLA = NIWKFR 1659 NRWLA = NRWKFR 1660 ENDIF 1661 IF (NEW.EQ.0) THEN 1662 COND1 = COND 1663 IF (QREPET) THEN 1664 IWK(NIWLA) = 1 1665 ELSE 1666 IWK(NIWLA) = 0 1667 ENDIF 1668 IF (MPRTIM.NE.0) CALL MONON(3) 1669 CALL N2FACT(N,M1,N,ML,MU,A,QA,COND1,IRANK,IOPT,IFAIL,LIWL, 1670 $ IWK(NIWLA),NILUSE,LRWL,RWK(NRWLA),NRLUSE) 1671 IF (MPRTIM.NE.0) CALL MONOFF(3) 1672C Exit Repeat If ... 1673 IF(IFAIL.NE.0) THEN 1674 IERR = 80 1675 GOTO 4299 1676 ENDIF 1677 IF (.NOT.QLINIT) THEN 1678 NIWKFR = NIWKFR+NILUSE 1679 NRWKFR = NRWKFR+NRLUSE 1680C Store lengths of currently required workspaces 1681 IWK(18) = NIWKFR-1 1682 IWK(19) = NRWKFR-1 1683 ENDIF 1684 SENS1 = RWK(NRWLA) 1685 ENDIF 1686 QLINIT = .TRUE. 1687C -------------------------------------------------------- 1688C 3.1.2 Solution of linear (N,N)-system 1689 IF(NEW.EQ.0) THEN 1690 IF (MPRTIM.NE.0) CALL MONON(4) 1691 CALL N2SOLV(N,M1,N,ML,MU,A,QA,T1,T2,IRANK,IOPT,IFAIL,LIWL, 1692 $ IWK(NIWLA),IDUMMY,LRWL,RWK(NRWLA),IDUMMY) 1693 IF (MPRTIM.NE.0) CALL MONOFF(4) 1694C Exit Repeat If ... 1695 IF(IFAIL.NE.0) THEN 1696 IERR = 81 1697 GOTO 4299 1698 ENDIF 1699 IF(.NOT.QREPET.AND.IRANK.NE.0)THEN 1700 DO 312 L1=1,N 1701 QU(L1)=T1(L1) 1702312 CONTINUE 1703 ENDIF 1704 ELSE 1705 ALFA1=ZERO 1706 ALFA2=ZERO 1707 DO 3121 I=1,N 1708 ALFA1=ALFA1+DX(I)*DXQ(I)/XW(I)**2 1709 ALFA2=ALFA2+DX(I)**2/XW(I)**2 17103121 CONTINUE 1711 ALFA=ALFA1/ALFA2 1712 BETA=ONE-ALFA 1713 DO 3122 I=1,N 1714 T2(I)=(DXQ(I)+(FCA-ONE)*ALFA*DX(I))/BETA 17153122 CONTINUE 1716 IF(NEW.EQ.1) THEN 1717 DO 3123 I=1,N 1718 DXSAVE(I,1)=DX(I) 17193123 CONTINUE 1720 ENDIF 1721 DO 3124 I=1,N 1722 DXSAVE(I,NEW+1)=T2(I) 1723 DX(I)=T2(I) 1724 T2(I)=T2(I)/XW(I) 17253124 CONTINUE 1726 ENDIF 1727C -------------------------------------------------------- 1728C 3.2 Evaluation of scaled natural level function SUMX 1729C scaled maximum error norm CONV 1730C evaluation of (scaled) standard level function 1731C DLEVF ( DLEVF only, if MPRMON.GE.2 ) 1732C and computation of ordinary Newton corrections DX( 1733C N) 1734 CALL N2LVLS(N,T2,XW,F,DX,CONV,SUMX,DLEVF,MPRMON,NEW.EQ.0) 1735 DO 32 L1=1,N 1736 XA(L1)=X(L1) 173732 CONTINUE 1738 SUMXA = SUMX 1739 DLEVXA = DSQRT(SUMXA/DBLE(FLOAT(N))) 1740 CONVA = CONV 1741 DXANRM = WNORM(N,DX,XW) 1742C -------------------------------------------------------- 1743C 3.3 A - priori estimate of damping factor FC 1744 QREDU = .FALSE. 1745 IF(NITER.NE.0.AND.NONLIN.NE.1)THEN 1746CWei; IF(NEW.EQ.0.AND.(IRANK.EQ.N.OR.IRANKA.EQ.N).OR. 1747CWei; * QREPET)THEN 1748 IF(NEW.EQ.0.OR.QREPET)THEN 1749C ------------------------------------------------------ 1750C 3.3.1 Computation of the denominator of a-priori 1751C estimate 1752 FCDNM = ZERO 1753 DO 331 L1=1,N 1754 FCDNM=FCDNM+((DX(L1)-DXQ(L1))/XW(L1))**2 1755331 CONTINUE 1756 IF(IRANK.NE.N)THEN 1757C -------------------------------------------- 1758C 3.3.2 Rank-deficient case (if previous rank 1759C was full) computation of the projected 1760C denominator of a-priori estimate 1761 DO 332 L1=1,N 1762 T1(L1)=DXQ(L1)/XW(L1) 1763332 CONTINUE 1764C Norm of projection of reduced component T1(N) 1765 CALL N2PRJN(N,IRANK,DEL,T1,RWK(NRWLA+1),T2,QA, 1766 $ IWK(NIWLA+2)) 1767 FCDNM = FCDNM-DEL 1768 ENDIF 1769 FCDNM = FCDNM*SUMX 1770C ------------------------------------------------------ 1771C 3.3.3 New damping factor 1772 IF(FCDNM.GT.FCNUMP*FCMIN2 .OR. 1773 $ (NONLIN.EQ.4 .AND. FCA**2*FCNUMP .LT. 4.0D0*FCDNM)) 1774 $ THEN 1775 DMYPRI = FCA*DSQRT(FCNUMP/FCDNM) 1776 FCPRI = DMIN1(DMYPRI,ONE) 1777 IF (NONLIN.EQ.4) FCPRI = DMIN1(HALF*DMYPRI,ONE) 1778 ELSE 1779 FCPRI = ONE 1780C$Test-begin 1781 DMYPRI = -1.0D0 1782C$Test-end 1783 ENDIF 1784C$Test-begin 1785 IF (MPRMON.GE.5) THEN 1786 WRITE(LUMON,33201) FCPRI, FC, FCA, DMYPRI, FCNUMP, 1787 $ FCDNM 178833201 FORMAT (/, ' +++ apriori estimate +++', /, 1789 $ ' FCPRI = ', D18.10, ' FC = ', D18.10, /, 1790 $ ' FCA = ', D18.10, ' DMYPRI = ', D18.10, /, 1791 $ ' FCNUMP = ', D18.10, ' FCDNM = ', D18.10, /, 1792 $ ' ++++++++++++++++++++++++', /) 1793 ENDIF 1794C$Test-end 1795 FC = FCPRI 1796 IF ( IRANK.EQ.N .OR. IRANK.LE.MINRNK ) 1797 $ FC=DMAX1(FC,FCMIN) 1798 IF (QBDAMP) THEN 1799 FCBH = FCA*FCBND 1800 IF (FC.GT.FCBH) THEN 1801 FC = FCBH 1802 IF (MPRMON.GE.4) 1803 $ WRITE(LUMON,*)' *** incr. rest. act. (a prio) ***' 1804 ENDIF 1805 FCBH = FCA/FCBND 1806 IF (FC.LT.FCBH) THEN 1807 FC = FCBH 1808 IF (MPRMON.GE.4) 1809 $ WRITE(LUMON,*)' *** decr. rest. act. (a prio) ***' 1810 ENDIF 1811 ENDIF 1812 ENDIF 1813 QREDU = FC.LT.FCMIN 1814 ENDIF 1815 QREPET = .FALSE. 1816CWEI 1817 IF (IORMON.GE.2) THEN 1818 SUMXA2=SUMXA1 1819 SUMXA1=SUMXA0 1820 SUMXA0=DLEVXA 1821 IF (SUMXA0.EQ.ZERO) SUMXA0=SMALL 1822C Check convergence rates (linear and superlinear) 1823C ICONV : Convergence indicator 1824C =0: No convergence indicated yet 1825C =1: Damping factor is 1.0d0 1826C =2: Superlinear convergence detected (alpha >=1.2) 1827C =3: Quadratic convergence detected (alpha > 1.8) 1828 FCMON = DMIN1(FC,FCMON) 1829 IF (FCMON.LT.ONE) THEN 1830 ICONV = 0 1831 ALPHAE = ZERO 1832 ENDIF 1833 IF (FCMON.EQ.ONE .AND. ICONV.EQ.0) ICONV=1 1834 IF (NITER.GE.1) THEN 1835 CLIN1 = CLIN0 1836 CLIN0 = SUMXA0/SUMXA1 1837 ENDIF 1838 IF (ICONV.GE.1.AND.NITER.GE.2) THEN 1839 ALPHAK = ALPHAE 1840 ALPHAE = ZERO 1841 IF (CLIN1.LE.0.95D0) ALPHAE = DLOG(CLIN0)/DLOG(CLIN1) 1842 IF (ALPHAK.NE.ZERO) ALPHAK =0.5D0*(ALPHAE+ALPHAK) 1843 ALPHAA = DMIN1(ALPHAK,ALPHAE) 1844 CALPHK = CALPHA 1845 CALPHA = ZERO 1846 IF (ALPHAE.NE.ZERO) CALPHA = SUMXA1/SUMXA2**ALPHAE 1847 SUMXTE = DSQRT(CALPHA*CALPHK)*SUMXA1**ALPHAK-SUMXA0 1848 IF (ALPHAA.GE.1.2D0 .AND. ICONV.EQ.1) ICONV = 2 1849 IF (ALPHAA.GT.1.8D0) ICONV = 3 1850 IF (MPRMON.GE.4) WRITE(LUMON,32001) ICONV, ALPHAE, 1851 $ CALPHA, CLIN0, ALPHAK, SUMXTE 185232001 FORMAT(' ** ICONV: ',I1,' ALPHA: ',D9.2, 1853 $ ' CONST-ALPHA: ',D9.2,' CONST-LIN: ',D9.2,' **', 1854 $ /,' **',11X,'ALPHA-POST: ',D9.2,' CHECK: ',D9.2, 1855 $ 25X,'**') 1856 IF ( ICONV.GE.2 .AND. ALPHAA.LT.0.9D0 ) THEN 1857 IF (IORMON.EQ.3) THEN 1858 IERR = 4 1859 GOTO 4299 1860 ELSE 1861 QMSTOP = .TRUE. 1862 ENDIF 1863 ENDIF 1864 ENDIF 1865 ENDIF 1866 FCMON = FC 1867C 1868 IF (MPRMON.GE.2) THEN 1869 IF (MPRTIM.NE.0) CALL MONON(5) 1870 CALL N2PRV1(DLEVF,DLEVXA,FCKEEP,NITER,NEW,IRANK,MPRMON, 1871 $ LUMON,QMIXIO,COND1) 1872 IF (MPRTIM.NE.0) CALL MONOFF(5) 1873 ENDIF 1874 IF(.NOT.QREDU)THEN 1875C -------------------------------------------------------- 1876C 3.4 Save natural level for later computations of 1877C corrector and print iterate 1878 FCNUMK = SUMX 1879 NRED = 0 1880 QREP = .FALSE. 1881C QREP = ITER .GT. ITMAX or QREP = ITER .GT. 0 1882C Damping-factor reduction loop 1883C ================================ 1884C DO (Until) 188534 CONTINUE 1886C ------------------------------------------------------ 1887C 3.5 Preliminary new iterate 1888 DO 35 L1=1,N 1889 X(L1)=XA(L1)+DX(L1)*FC 189035 CONTINUE 1891C ----------------------------------------------------- 1892C 3.5.2 Exit, if problem is specified as being linear 1893C Exit Repeat If ... 1894 IF( NONLIN.EQ.1 )THEN 1895 IERR = 0 1896 GOTO 4299 1897 ENDIF 1898C ------------------------------------------------------ 1899C 3.6.1 Computation of the residual vector 1900 IF (MPRTIM.NE.0) CALL MONON(1) 1901 CALL FCN(N,X,F,IFAIL) 1902 IF (MPRTIM.NE.0) CALL MONOFF(1) 1903 NFCN = NFCN+1 1904C Exit, if ... 1905 IF (IFAIL.LT.0) THEN 1906 IERR = 82 1907 GOTO 4299 1908 ENDIF 1909 IF(IFAIL.EQ.1 .OR. IFAIL.EQ.2) THEN 1910 IF (IFAIL.EQ.1) THEN 1911 FCREDU = HALF 1912 ELSE 1913 FCREDU = F(1) 1914C Exit, If ... 1915 IF (FCREDU.LE.0 .OR. FCREDU.GE.1) THEN 1916 IERR = 83 1917 GOTO 4299 1918 ENDIF 1919 ENDIF 1920 IF (MPRMON.GE.2) THEN 192136101 FORMAT(8X,I2,' FCN could not be evaluated ', 1922 $ 8X,F7.5,4X,I2,6X,I4) 1923 WRITE(LUMON,36101)NITER,FC,NEW,IRANK 1924 ENDIF 1925 FCH = FC 1926 FC = FCREDU*FC 1927 IF (FCH.GT.FCMIN) FC = DMAX1(FC,FCMIN) 1928 IF (QBDAMP) THEN 1929 FCBH = FCH/FCBND 1930 IF (FC.LT.FCBH) THEN 1931 FC = FCBH 1932 IF (MPRMON.GE.4) WRITE(LUMON,*) 1933 $ ' *** decr. rest. act. (FCN redu.) ***' 1934 ENDIF 1935 ENDIF 1936 IF (FC.LT.FCMIN) THEN 1937 IERR = 3 1938 GOTO 4299 1939 ENDIF 1940C Break DO (Until) ... 1941 GOTO 3109 1942 ENDIF 1943 DO 361 L1=1,N 1944 T1(L1)=F(L1)*FW(L1) 1945361 CONTINUE 1946C ------------------------------------------------------ 1947C 3.6.2 Solution of linear (N,N)-system 1948 IF (QREPET) THEN 1949 IWK(NIWLA) = 1 1950 ELSE 1951 IWK(NIWLA) = 0 1952 ENDIF 1953 IF (MPRTIM.NE.0) CALL MONON(4) 1954 CALL N2SOLV(N,M1,N,ML,MU,A,QA,T1,T2,IRANK,IOPT,IFAIL, 1955 $ LIWL,IWK(NIWLA),IDUMMY,LRWL,RWK(NRWLA), 1956 $ IDUMMY) 1957 IF (MPRTIM.NE.0) CALL MONOFF(4) 1958C Exit Repeat If ... 1959 IF(IFAIL.NE.0) THEN 1960 IERR = 81 1961 GOTO 4299 1962 ENDIF 1963 IF(NEW.GT.0) THEN 1964 DO 3630 I=1,N 1965 DXQ(I) = T2(I)*XWA(I) 19663630 CONTINUE 1967 DO 363 ILOOP=1,NEW 1968 SUM1=ZERO 1969 SUM2=ZERO 1970 DO 3631 I=1,N 1971 SUM1=SUM1+(DXQ(I)*DXSAVE(I,ILOOP))/ XW(I)**2 1972 SUM2=SUM2+(DXSAVE(I,ILOOP)/XW(I))**2 19733631 CONTINUE 1974 BETA=SUM1/SUM2 1975 DO 3632 I=1,N 1976 DXQ(I)=DXQ(I)+BETA*DXSAVE(I,ILOOP+1) 1977 T2(I) = DXQ(I)/XW(I) 19783632 CONTINUE 1979363 CONTINUE 1980 ENDIF 1981C ------------------------------------------------------ 1982C 3.6.3 Evaluation of scaled natural level function 1983C SUMX 1984C scaled maximum error norm CONV and evaluation 1985C of (scaled) standard level function DLEVF 1986 CALL N2LVLS(N,T2,XW,F,DXQ,CONV,SUMX,DLEVFN,MPRMON, 1987 $ NEW.EQ.0) 1988 DXNRM = WNORM(N,DXQ,XW) 1989C ------------------------------------------------------ 1990C 3.6.4 Convergence test 1991C Exit Repeat If ... 1992 IF ( DXNRM.LE.RTOL .AND. DXANRM.LE.RSMALL .AND. 1993 $ FC.EQ.ONE ) THEN 1994 IERR = 0 1995 GOTO 4299 1996 ENDIF 1997C 1998 FCA = FC 1999C ---------------------------------------------------- 2000C 3.6.5 Evaluation of reduced damping factor 2001 TH = FCA-ONE 2002 FCDNM = ZERO 2003 DO 39 L1=1,N 2004 FCDNM=FCDNM+((DXQ(L1)+TH*DX(L1))/XW(L1))**2 200539 CONTINUE 2006 IF (FCDNM.NE.ZERO) THEN 2007 DMYCOR = FCA*FCA*HALF*DSQRT(FCNUMK/FCDNM) 2008 ELSE 2009 DMYCOR = 1.0D+35 2010 ENDIF 2011 IF (NONLIN.LE.3) THEN 2012 FCCOR = DMIN1(ONE,DMYCOR) 2013 ELSE 2014 FCCOR = DMIN1(ONE,HALF*DMYCOR) 2015 ENDIF 2016C$Test-begin 2017 IF (MPRMON.GE.5) THEN 2018 WRITE(LUMON,39001) FCCOR, FC, DMYCOR, FCNUMK, 2019 $ FCDNM, FCA 202039001 FORMAT (/, ' +++ corrector computation +++', /, 2021 $ ' FCCOR = ', D18.10, ' FC = ', D18.10, /, 2022 $ ' DMYCOR = ', D18.10, ' FCNUMK = ', D18.10, /, 2023 $ ' FCDNM = ', D18.10, ' FCA = ', D18.10, /, 2024 $ ' +++++++++++++++++++++++++++++', /) 2025 ENDIF 2026C$Test-end 2027C ------------------------------------------------------ 2028C 3.7 Natural monotonicity test 2029 IF(SUMX.GT.SUMXA)THEN 2030C ---------------------------------------------------- 2031C 3.8 Output of iterate 2032 IF(MPRMON.GE.3) THEN 2033 IF (MPRTIM.NE.0) CALL MONON(5) 2034 CALL N2PRV2(DLEVFN,DSQRT(SUMX/DBLE(FLOAT(N))),FC, 2035 $ NITER,MPRMON,LUMON,QMIXIO,'*') 2036 IF (MPRTIM.NE.0) CALL MONOFF(5) 2037 ENDIF 2038 IF (QMSTOP) THEN 2039 IERR = 4 2040 GOTO 4299 2041 ENDIF 2042 FC = DMIN1(FCCOR,HALF*FC) 2043 IF ((IRANK.EQ.N .OR. IRANK.EQ.MINRNK) .AND. 2044 $ FCA.GT.FCMIN) 2045 $ FC=DMAX1(FC,FCMIN) 2046 IF (QBDAMP) THEN 2047 FCBH = FCA/FCBND 2048 IF (FC.LT.FCBH) THEN 2049 FC = FCBH 2050 IF (MPRMON.GE.4) WRITE(LUMON,*) 2051 $ ' *** decr. rest. act. (a post) ***' 2052 ENDIF 2053 ENDIF 2054CWEI 2055 FCMON = FC 2056C 2057C$Test-begin 2058 IF (MPRMON.GE.5) THEN 2059 WRITE(LUMON,39002) FC 206039002 FORMAT (/, ' +++ corrector setting 1 +++', /, 2061 $ ' FC = ', D18.10, /, 2062 $ ' +++++++++++++++++++++++++++', /) 2063 ENDIF 2064C$Test-end 2065 QREP = .TRUE. 2066 NCORR = NCORR+1 2067 NRED = NRED+1 2068C ---------------------------------------------------- 2069C 3.10 If damping factor is too small: 2070C Refresh Jacobian,if current Jacobian was computed 2071C by a Rank1-update, else reduce pseudo-rank 2072 QREDU = FC.LT.FCMIN.OR.NEW.GT.0.AND.NRED.GT.1 2073 ELSE 2074 IF (.NOT.QREP .AND. FCCOR.GT.SIGMA2*FC) THEN 2075 IF(MPRMON.GE.3) THEN 2076 IF (MPRTIM.NE.0) CALL MONON(5) 2077 CALL N2PRV2(DLEVFN,DSQRT(SUMX/DBLE(FLOAT(N))),FC, 2078 $ NITER,MPRMON,LUMON,QMIXIO,'+') 2079 IF (MPRTIM.NE.0) CALL MONOFF(5) 2080 ENDIF 2081 FC = FCCOR 2082C$Test-begin 2083 IF (MPRMON.GE.5) THEN 2084 WRITE(LUMON,39003) FC 208539003 FORMAT (/, ' +++ corrector setting 2 +++', /, 2086 $ ' FC = ', D18.10, /, 2087 $ ' +++++++++++++++++++++++++++', /) 2088 ENDIF 2089C$Test-end 2090 QREP = .TRUE. 2091 ELSE 2092 QNEXT = .TRUE. 2093 ENDIF 2094 ENDIF 20953109 CONTINUE 2096 IF(.NOT.(QNEXT.OR.QREDU)) GOTO 34 2097C UNTIL ( expression - negated above) 2098C End of damping-factor reduction loop 2099C ======================================= 2100 ENDIF 2101 IF(QREDU)THEN 2102C ------------------------------------------------------ 2103C 3.11 Restore former values for repeting step 2104C step 2105 NREJR1 = NREJR1+1 2106 DO 3111 L1=1,N 2107 X(L1)=XA(L1) 21083111 CONTINUE 2109 DO 3112 L1=1,N 2110 F(L1)=FA(L1) 21113112 CONTINUE 2112 DO 3113 L1=1,N 2113 DXQ(L1)=DXQA(L1) 21143113 CONTINUE 2115 IF(MPRMON.GE.2)THEN 211631130 FORMAT(8X,I2,' Not accepted damping ', 2117 $ 'factor',8X,F7.5,4X,I2,6X,I4) 2118 WRITE(LUMON,31130)NITER,FC,NEW,IRANK 2119 ENDIF 2120 FC = FCKEEP 2121 FCA = FCK2 2122 IF(NITER.EQ.0)THEN 2123 FC = FCMIN 2124 ENDIF 2125 IF(NEW.GT.0)THEN 2126 QGENJ = .TRUE. 2127 QJCRFR = .TRUE. 2128 QREDU = .FALSE. 2129 ELSE 2130C ------------------------------------------------ 2131C 3.12 Pseudo-rank reduction 2132 QREPET = .TRUE. 2133 DO 42 L1=1,N 2134 T1(L1)=QU(L1) 213542 CONTINUE 2136 IRANK = IRANK-1 2137 IF(IRANK.LT.MINRNK)THEN 2138 IERR = 3 2139 GOTO 4299 2140 ENDIF 2141 ENDIF 2142 ENDIF 2143 IF(.NOT.(.NOT.QREDU)) GOTO 3 2144C UNTIL ( expression - negated above) 2145C 2146C End of pseudo-rank reduction loop 2147C ================================= 2148 IF (QNEXT) THEN 2149C ------------------------------------------------------ 2150C 4 Preparations to start the following iteration step 2151C ------------------------------------------------------ 2152C 4.1 Print values 2153 IF(MPRMON.GE.3) THEN 2154 IF (MPRTIM.NE.0) CALL MONON(5) 2155 CALL N2PRV2(DLEVFN,DSQRT(SUMX/DBLE(FLOAT(N))),FC,NITER+1, 2156 $ MPRMON,LUMON,QMIXIO,'*') 2157 IF (MPRTIM.NE.0) CALL MONOFF(5) 2158 ENDIF 2159C Print the natural level of the current iterate and return 2160C it in one-step mode 2161 SUMX = SUMXA 2162 IF(MPRSOL.GE.2.AND.NITER.NE.0) THEN 2163 IF (MPRTIM.NE.0) CALL MONON(5) 2164 CALL N2SOUT(N,XA,2,IOPT,RWK,LRWK,IWK,LIWK,MPRSOL,LUSOL) 2165 IF (MPRTIM.NE.0) CALL MONOFF(5) 2166 ELSE IF(MPRSOL.GE.1.AND.NITER.EQ.0)THEN 2167 IF (MPRTIM.NE.0) CALL MONON(5) 2168 CALL N2SOUT(N,XA,1,IOPT,RWK,LRWK,IWK,LIWK,MPRSOL,LUSOL) 2169 IF (MPRTIM.NE.0) CALL MONOFF(5) 2170 ENDIF 2171 NITER = NITER+1 2172 DLEVF = DLEVFN 2173C Exit Repeat If ... 2174 IF(NITER.GE.NITMAX)THEN 2175 IERR = 2 2176 GOTO 4299 2177 ENDIF 2178 FCKEEP = FC 2179C ------------------------------------------------------ 2180C 4.2 Return, if in one-step mode 2181C Exit Subroutine If ... 2182 IF (MODE.EQ.1) THEN 2183 IWK(18)=NIWLA-1 2184 IWK(19)=NRWLA-1 2185 IOPT(1)=1 2186 RETURN 2187 ENDIF 2188 ENDIF 2189 GOTO 2 2190C End Repeat 21914299 CONTINUE 2192C 2193C End of main iteration loop 2194C ========================== 2195C ---------------------------------------------------------- 2196C 9 Exits 2197C ---------------------------------------------------------- 2198C 9.1 Solution exit 2199 APREC = -1.0D0 2200C 2201 IF(IERR.EQ.0 .OR. IERR.EQ.4)THEN 2202 IF (NONLIN.NE.1) THEN 2203 IF ( IERR.EQ.0 ) THEN 2204 APREC = DSQRT(SUMX/DBLE(FLOAT(N))) 2205 DO 91 L1=1,N 2206 X(L1)=X(L1)+DXQ(L1) 220791 CONTINUE 2208 ELSE 2209 APREC = DSQRT(SUMXA/DBLE(FLOAT(N))) 2210 IF (ALPHAA.GT.ZERO .AND. IORMON.EQ.3) THEN 2211 DO 92 L1=1,N 2212 X(L1)=X(L1)+DX(L1) 221392 CONTINUE 2214 ENDIF 2215 ENDIF 2216 IF(IRANK.LT.N)THEN 2217 IERR = 1 2218 ENDIF 2219C Print final monitor output 2220 IF(MPRMON.GE.2) THEN 2221 IF (IERR.EQ.0) THEN 2222 IF (MPRTIM.NE.0) CALL MONON(5) 2223 CALL N2PRV2(DLEVFN,DSQRT(SUMX/DBLE(FLOAT(N))),FC, 2224 $ NITER+1,MPRMON,LUMON,QMIXIO,'*') 2225 IF (MPRTIM.NE.0) CALL MONOFF(5) 2226 ELSE IF (IORMON.EQ.3) THEN 2227 IF (MPRTIM.NE.0) CALL MONON(5) 2228 CALL N2PRV1(DLEVFN,DSQRT(SUMXA/DBLE(FLOAT(N))),FC, 2229 $ NITER,NEW,IRANK,MPRMON,LUMON,QMIXIO,COND1) 2230 IF (MPRTIM.NE.0) CALL MONOFF(5) 2231 ENDIF 2232 ENDIF 2233 IF ( IORMON.GE.2 ) THEN 2234 IF ( ICONV.LE.1 .AND. ALPHAE .NE. ZERO 2235 $ .AND. ALPHAK .NE. ZERO ) IERR = 5 2236 ENDIF 2237C 2238 IF(MPRMON.GE.1.AND. IERR.NE.1) THEN 223991001 FORMAT(///' Solution of nonlinear system ', 2240 $ 'of equations obtained within ',I3, 2241 $ ' iteration steps',//,' Achieved relative accuracy',D10.3) 2242 IF (IERR.EQ.4) THEN 2243 WRITE(LUMON,91001) NITER,APREC 2244 ELSE 2245 WRITE(LUMON,91001) NITER+1,APREC 2246 ENDIF 2247 ENDIF 2248 ELSE 2249 IF(MPRMON.GE.1) THEN 225091002 FORMAT(///' Solution of linear system ', 2251 $ 'of equations obtained by NLEQ2',//,' No estimate ', 2252 $ 'available for the achieved relative accuracy') 2253 WRITE(LUMON,91002) 2254 ENDIF 2255 ENDIF 2256 ENDIF 2257C ---------------------------------------------------------- 2258C 9.2 Fail exit messages 2259C ---------------------------------------------------------- 2260C 9.2.1 Termination at stationary point 2261 IF(IERR.EQ.1.AND.MPRERR.GE.1)THEN 226292101 FORMAT(/,' Iteration terminates at stationary point',/) 2263 WRITE(LUERR,92101) 2264 ENDIF 2265C ---------------------------------------------------------- 2266C 9.2.2 Termination after more than NITMAX iterations 2267 IF(IERR.EQ.2.AND.MPRERR.GE.1)THEN 226892201 FORMAT(/,' Iteration terminates after NITMAX ', 2269 $ '=',I3,' Iteration steps') 2270 WRITE(LUERR,92201)NITMAX 2271 ENDIF 2272C ---------------------------------------------------------- 2273C 9.2.3 Newton method fails to converge 2274 IF(IERR.EQ.3.AND.MPRERR.GE.1)THEN 227592301 FORMAT(/,' Newton method fails to converge') 2276 WRITE(LUERR,92301) 2277 ENDIF 2278CWEI 2279C ---------------------------------------------------------- 2280C 9.2.4.1 Superlinear convergence slowed down 2281 IF(IERR.EQ.4.AND.MPRERR.GE.1)THEN 228292401 FORMAT(/,' Warning: Monotonicity test failed after ',A, 2283 $ ' convergence was already checked;',/, 2284 $ ' RTOL requirement may be too stringent',/) 228592402 FORMAT(/,' Warning: ',A,' convergence slowed down;',/, 2286 $ ' RTOL requirement may be too stringent',/) 2287 IF (QMSTOP) THEN 2288 IF (ICONV.EQ.2) WRITE(LUERR,92401) 'superlinear' 2289 IF (ICONV.EQ.3) WRITE(LUERR,92401) 'quadratic' 2290 ELSE 2291 IF (ICONV.EQ.2) WRITE(LUERR,92402) 'superlinear' 2292 IF (ICONV.EQ.3) WRITE(LUERR,92402) 'quadratic' 2293 ENDIF 2294 ENDIF 2295C ---------------------------------------------------------- 2296C 9.2.4.2 Convergence criterion satisfied before superlinear 2297C convergence has been established 2298 IF(IERR.EQ.5.AND.MPRERR.GE.1)THEN 229992410 FORMAT(/,' Warning: No quadratic or superlinear convergence ', 2300 $ 'established yet',/, 2301 $ 10X,'your solution may perhaps may be less accurate ', 2302 $ /,10X,'as indicated by the standard error estimate') 2303 WRITE(LUERR,92410) 2304 ENDIF 2305C ---------------------------------------------------------- 2306C 9.2.5 Error exit due to linear solver routine N2FACT 2307 IF(IERR.EQ.80.AND.MPRERR.GE.1)THEN 230892501 FORMAT(/,' Error ',I5,' signalled by linear solver N2FACT') 2309 WRITE(LUERR,92501) IFAIL 2310 ENDIF 2311C ---------------------------------------------------------- 2312C 9.2.6 Error exit due to linear solver routine N2SOLV 2313 IF(IERR.EQ.81.AND.MPRERR.GE.1)THEN 231492601 FORMAT(/,' Error ',I5,' signalled by linear solver N2SOLV') 2315 WRITE(LUERR,92601) IFAIL 2316 ENDIF 2317C ---------------------------------------------------------- 2318C 9.2.7 Error exit due to fail of user function FCN 2319 IF(IERR.EQ.82.AND.MPRERR.GE.1)THEN 232092701 FORMAT(/,' Error ',I5,' signalled by user function FCN') 2321 WRITE(LUERR,92701) IFAIL 2322 ENDIF 2323C ---------------------------------------------------------- 2324C 9.2.8 Error exit due to fail of user function JAC 2325 IF(IERR.EQ.83.AND.MPRERR.GE.1)THEN 232692801 FORMAT(/,' Error ',I5,' signalled by user function JAC') 2327 WRITE(LUERR,92801) IFAIL 2328 ENDIF 2329 IF(IERR.GE.80.AND.IERR.LE.83) IWK(23) = IFAIL 2330 IF ((IERR.EQ.82.OR.IERR.EQ.83).AND.NITER.LE.1.AND.MPRERR.GE.1) 2331 $ THEN 2332 WRITE (LUERR,92810) 233392810 FORMAT(' Try to find a better initial guess for the solution') 2334 ENDIF 2335C ---------------------------------------------------------- 2336C 9.3 Common exit 2337 IF (MPRERR.GE.3.AND.IERR.NE.0.AND.IERR.NE.4.AND.NONLIN.NE.1) 2338 $ THEN 233993100 FORMAT(/,' Achieved relative accuracy',D10.3,2X) 2340 WRITE(LUERR,93100)CONVA 2341 APREC = CONVA 2342 ENDIF 2343 IF(MPRMON.GE.1)THEN 234493001 FORMAT(/,3X,'Subcondition ( 1,',I4,')',1X,D10.3,2X, 2345 $ /,3X,'Sensitivity ( 1,',I4,')',1X,D10.3,2X,/) 2346 WRITE(LUMON,93001)IRANK,COND1,IRANK,SENS1 2347 ENDIF 2348 RTOL = APREC 2349 SUMX = SUMXA 2350 IF(MPRSOL.GE.2.AND.NITER.NE.0) THEN 2351 IF (MPRTIM.NE.0) CALL MONON(5) 2352 CALL N2SOUT(N,XA,2,IOPT,RWK,LRWK,IWK,LIWK,MPRSOL,LUSOL) 2353 IF (MPRTIM.NE.0) CALL MONOFF(5) 2354 ELSE IF(MPRSOL.GE.1.AND.NITER.EQ.0)THEN 2355 IF (MPRTIM.NE.0) CALL MONON(5) 2356 CALL N2SOUT(N,XA,1,IOPT,RWK,LRWK,IWK,LIWK,MPRSOL,LUSOL) 2357 IF (MPRTIM.NE.0) CALL MONOFF(5) 2358 ENDIF 2359 IF (IERR.NE.4) NITER = NITER+1 2360 DLEVF = DLEVFN 2361 IF(MPRSOL.GE.1)THEN 2362C Print Solution or final iteration vector 2363 IF(IERR.EQ.0)THEN 2364 MODEFI = 3 2365 ELSE 2366 MODEFI = 4 2367 ENDIF 2368 IF (MPRTIM.NE.0) CALL MONON(5) 2369 CALL N2SOUT(N,X,MODEFI,IOPT,RWK,LRWK,IWK,LIWK,MPRSOL,LUSOL) 2370 IF (MPRTIM.NE.0) CALL MONOFF(5) 2371 ENDIF 2372C Return the latest internal scaling to XSCAL 2373 DO 93 I=1,N 2374 XSCAL(I)=XW(I) 237593 CONTINUE 2376C End of exits 2377C End of subroutine N2INT 2378 RETURN 2379 END 2380C 2381 SUBROUTINE N2SCAL(N,X,XA,XSCAL,XW,ISCAL,QINISC,IOPT,LRWK,RWK) 2382C* Begin Prologue SCALE 2383 INTEGER N 2384 DOUBLE PRECISION X(N),XSCAL(N),XA(N),XW(N) 2385 INTEGER ISCAL 2386 LOGICAL QINISC 2387 INTEGER IOPT(50),LRWK 2388 DOUBLE PRECISION RWK(LRWK) 2389C ------------------------------------------------------------ 2390C 2391C* Summary : 2392C 2393C S C A L E : To be used in connection with NLEQ2 . 2394C Computation of the internal scaling vector XW used for the 2395C Jacobian matrix, the iterate vector and it's related 2396C vectors - especially for the solution of the linear system 2397C and the computations of norms to avoid numerical overflow. 2398C 2399C* Input parameters 2400C ================ 2401C 2402C N Int Number of unknowns 2403C X(N) Dble Current iterate 2404C XA(N) Dble Previous iterate 2405C XSCAL(N) Dble User scaling passed from parameter XSCAL 2406C of interface routine NLEQ2 2407C ISCAL Int Option ISCAL passed from IOPT-field 2408C (for details see description of IOPT-fields) 2409C QINISC Logical = .TRUE. : Initial scaling 2410C = .FALSE. : Subsequent scaling 2411C IOPT(50) Int Options array passed from NLEQ2 parameter list 2412C LRWK Int Length of real workspace 2413C RWK(LRWK) Dble Real workspace (see description above) 2414C 2415C* Output parameters 2416C ================= 2417C 2418C XW(N) Dble Scaling vector computed by this routine 2419C All components must be positive. The follow- 2420C ing relationship between the original vector 2421C X and the scaled vector XSCAL holds: 2422C XSCAL(I) = X(I)/XW(I) for I=1,...N 2423C 2424C* Subroutines called: D1MACH 2425C 2426C* Machine constants used 2427C ====================== 2428C 2429 DOUBLE PRECISION SMALL 2430C 2431C ------------------------------------------------------------ 2432C* End Prologue 2433 EXTERNAL D1MACH 2434 INTRINSIC DABS,DMAX1 2435 DOUBLE PRECISION D1MACH,HALF 2436 PARAMETER (HALF=0.5D0) 2437 INTEGER MPRMON,LUMON 2438 SMALL = D1MACH(6) 2439C* Begin 2440 DO 1 L1=1,N 2441 IF (ISCAL.EQ.1) THEN 2442 XW(L1) = XSCAL(L1) 2443 ELSE 2444 XW(L1)=DMAX1(XSCAL(L1),(DABS(X(L1))+DABS(XA(L1)))*HALF,SMALL) 2445 ENDIF 24461 CONTINUE 2447C$Test-Begin 2448 MPRMON = IOPT(13) 2449 IF (MPRMON.GE.6) THEN 2450 LUMON = IOPT(14) 2451 WRITE(LUMON,*) ' ' 2452 WRITE(LUMON,*) ' ++++++++++++++++++++++++++++++++++++++++++' 2453 WRITE(LUMON,*) ' X-components Scaling-components ' 2454 WRITE(LUMON,10) (X(L1), XW(L1), L1=1,N) 245510 FORMAT(' ',D18.10,' ',D18.10) 2456 WRITE(LUMON,*) ' ++++++++++++++++++++++++++++++++++++++++++' 2457 WRITE(LUMON,*) ' ' 2458 ENDIF 2459C$Test-End 2460C End of subroutine N2SCAL 2461 RETURN 2462 END 2463C 2464 SUBROUTINE N2SCRF(M,N,A,FW) 2465C* Begin Prologue SCROWF 2466 INTEGER M,N 2467 DOUBLE PRECISION A(M,N),FW(M) 2468C ------------------------------------------------------------ 2469C 2470C* Summary : 2471C 2472C S C R O W F : Row Scaling of a (M,N)-matrix in full storage 2473C mode 2474C 2475C* Input parameters (* marks inout parameters) 2476C =========================================== 2477C 2478C M Int Number of rows of the matrix 2479C N Int Number of columns of the matrix 2480C * A(M,N) Dble Matrix to be scaled 2481C 2482C* Output parameters 2483C ================= 2484C 2485C FW(M) Dble Row scaling factors - FW(i) contains 2486C the factor by which the i-th row of A 2487C has been multiplied 2488C 2489C ------------------------------------------------------------ 2490C* End Prologue 2491 INTRINSIC DABS 2492 DOUBLE PRECISION ONE 2493 PARAMETER (ONE=1.0D0) 2494 DOUBLE PRECISION ZERO 2495 PARAMETER (ZERO=0.0D0) 2496 INTEGER J,K 2497 DOUBLE PRECISION S1,S2 2498C* Begin 2499 DO 1 K=1,M 2500 S1=ZERO 2501 DO 2 J=1,N 2502 S2=DABS(A(K,J)) 2503 IF (S2.GT.S1) S1=S2 25042 CONTINUE 2505 IF (S1.GT.ZERO) THEN 2506 S1=ONE/S1 2507 FW(K)=S1 2508 DO 3 J=1,N 2509 A(K,J)=A(K,J)*S1 25103 CONTINUE 2511 ELSE 2512 FW(K)=ONE 2513 ENDIF 25141 CONTINUE 2515C End of subroutine N1SCRF 2516 RETURN 2517 END 2518C 2519 SUBROUTINE N2FACT(N,LDA,LDAINV,ML,MU,A,AINV,COND,IRANK,IOPT, 2520 $IFAIL,LIWK,IWK,LAIWK,LRWK,RWK,LARWK) 2521C* Begin Prologue FACT 2522 INTEGER N,LDA,LDAINV,ML,MU 2523 DOUBLE PRECISION A(LDA,N),AINV(LDAINV,N) 2524 DOUBLE PRECISION COND 2525 INTEGER IRANK 2526 INTEGER IOPT(50) 2527 INTEGER IFAIL 2528 INTEGER LIWK 2529 INTEGER IWK(LIWK) 2530 INTEGER LAIWK,LRWK 2531 DOUBLE PRECISION RWK(LRWK) 2532 INTEGER LARWK 2533C ------------------------------------------------------------ 2534C 2535C* Summary : 2536C 2537C F A C T : Call linear algebra subprogram for factorization of 2538C a (N,N)-matrix with rank decision and casual compu- 2539C tation of the rank deficient pseudo-inverse matrix 2540C 2541C* Input parameters (* marks inout parameters) 2542C =========================================== 2543C 2544C N Int Order of the linear system 2545C LDA Int Leading dimension of the matrix array A 2546C LDAINV Int Leading dimension of the matrix array AINV 2547C ML Int Lower bandwidth of the matrix (only for 2548C banded systems) 2549C MU Int Upper bandwidth of the matrix (only for 2550C banded systems) 2551C * A(LDA,N) Dble Matrix storage. 2552C * COND Dble On Input, COND holds the maximum permitted 2553C subcondition for the prescribed rank 2554C On Output, it holds the estimated subcon- 2555C dition of A 2556C IOPT(50) Int Option vector passed from NLEQ2 2557C 2558C* Output parameters 2559C ================= 2560C 2561C AINV(LDAINV,N) Dble If matrix A is rank deficient, this array 2562C holds the pseudo-inverse of A 2563C IFAIL Int Error indicator returned by this routine: 2564C = 0 matrix decomposition successfull 2565C =10 supplied (integer) workspace too small 2566C 2567C* Workspace parameters 2568C ==================== 2569C 2570C LIWK Int Length of integer workspace passed to this 2571C routine (In) 2572C IWK(LIWK) Int Integer Workspace supplied for this routine 2573C LAIWK Int Length of integer Workspace used by this 2574C routine (out) 2575C LRWK Int Length of real workspace passed to this 2576C routine (In) 2577C RWK(LRWK) Dble Real Workspace supplied for this routine 2578C LARWK Int Length of real Workspace used by this 2579C routine (out) 2580C 2581C* Subroutines called: DECCON 2582C 2583C ------------------------------------------------------------ 2584C* End Prologue 2585 EXTERNAL DECCON 2586 INTRINSIC DABS 2587 DOUBLE PRECISION ONE 2588 PARAMETER (ONE=1.0D0) 2589 DOUBLE PRECISION ZERO 2590 PARAMETER (ZERO=0.0D0) 2591 INTEGER IREPET,MCON 2592C* Begin 2593 MPRERR=IOPT(11) 2594 LUERR=IOPT(12) 2595 LAIWK = N+2 2596 LARWK = 2*N+1 2597 IF (LIWK.GE.LAIWK.AND.LRWK.GE.LARWK) THEN 2598 MCON = 0 2599 IREPET = -IWK(1) 2600 IF (IREPET.EQ.0) IWK(2) = MCON 2601 CALL DECCON(A,LDA,N,MCON,N,N,IWK(2),IRANK,COND,RWK(2),IWK(3), 2602 $ IREPET,AINV,RWK(N+2),IFAIL) 2603 IF (IFAIL.EQ.-2 .AND. MPRERR.GT.0) WRITE(LUERR,10001) 260410001 FORMAT(1X, 2605 $ 'DECCON failed to compute rank-deficient QR-decomposition', 2606 $ /) 2607 IF(IRANK.NE.0)THEN 2608 COND = DABS(RWK(2)/RWK(IRANK+1)) 2609 RWK(1) = DABS(RWK(2)) 2610 ELSE 2611 COND = ONE 2612 RWK(1) = ZERO 2613 ENDIF 2614 ELSE 2615 IFAIL = 10 261610002 FORMAT(/,' Insuffient workspace for linear solver,', 2617 $ ' at least needed more needed : ',/, 2618 $ ' ',A,' workspace : ',I4) 2619 IF (LIWK.LT.LAIWK.AND.MPRERR.GT.0) 2620 $ WRITE(LUERR,10002) 'Integer',LAIWK-LIWK 2621 IF (LRWK.LT.LARWK.AND.MPRERR.GT.0) 2622 $ WRITE(LUERR,10002) 'Double',LARWK-LRWK 2623 ENDIF 2624 RETURN 2625 END 2626C 2627 SUBROUTINE N2SOLV(N,LDA,LDAINV,ML,MU,A,AINV,B,Z,IRANK,IOPT, 2628 $IFAIL,LIWK,IWK,LAIWK,LRWK,RWK,LARWK) 2629C* Begin Prologue SOLVE 2630 INTEGER N,LDA,LDAINV,ML,MU 2631 DOUBLE PRECISION A(LDA,N),AINV(LDAINV,N) 2632 DOUBLE PRECISION B(N),Z(N) 2633 INTEGER IRANK 2634 INTEGER IOPT(50) 2635 INTEGER IFAIL 2636 INTEGER LIWK 2637 INTEGER IWK(LIWK) 2638 INTEGER LRWK,LAIWK 2639 DOUBLE PRECISION RWK(LRWK) 2640 INTEGER LARWK 2641C ------------------------------------------------------------ 2642C 2643C* Summary : 2644C 2645C S O L V E : Call linear algebra subprogram for solution of 2646C the linear system A*Z = B 2647C 2648C* Parameters 2649C ========== 2650C 2651C N,LDA,LDAINV,ML,MU,A,AINV,IRANK,IOPT,IFAIL,LIWK,IWK,LAIWK,LRWK, 2652C RWK,LARWK : 2653C See description for subroutine N2FACT. 2654C B Dble In: Right hand side of the linear system 2655C Out: Rhs. transformed to the upper trian- 2656C gular part of the linear system 2657C Z Dble Out: Solution of the linear system 2658C 2659C Subroutines called: SOLCON 2660C 2661C ------------------------------------------------------------ 2662C* End Prologue 2663 EXTERNAL SOLCON 2664 INTEGER IREPET 2665C* Begin 2666 MCON = 0 2667 IREPET = -IWK(1) 2668 CALL SOLCON(A,LDA,N,MCON,N,N,Z,B,IWK(2),IRANK,RWK(2),IWK(3), 2669 $ IREPET,AINV,RWK(N+2)) 2670 IFAIL = 0 2671 RETURN 2672 END 2673C 2674 SUBROUTINE N2LVLS(N,DX1,XW,F,DXQ,CONV,SUMX,DLEVF,MPRMON,QDSCAL) 2675C* Begin Prologue LEVELS 2676 INTEGER N,MPRMON 2677 DOUBLE PRECISION DX1(N),XW(N),F(N),DXQ(N) 2678 DOUBLE PRECISION CONV,SUMX,DLEVF 2679 LOGICAL QDSCAL 2680C ------------------------------------------------------------ 2681C 2682C* Summary : 2683C 2684C L E V E L S : To be used in connection with NLEQ2 . 2685C provides descaled solution, error norm and level functions 2686C 2687C* Input parameters (* marks inout parameters) 2688C =========================================== 2689C 2690C N Int Number of parameters to be estimated 2691C DX1(N) Dble array containing the scaled Newton 2692C correction 2693C XW(N) Dble Array containing the scaling values 2694C F(N) Dble Array containing the residuum 2695C 2696C* Output parameters 2697C ================= 2698C 2699C DXQ(N) Dble Array containing the descaled Newton 2700C correction 2701C CONV Dble Scaled maximum norm of the Newton 2702C correction 2703C SUMX Dble Scaled natural level function value 2704C DLEVF Dble Standard level function value (only 2705C if needed for print) 2706C MPRMON Int Print information parameter (see 2707C driver routine NLEQ2 ) 2708C QDSCAL Logic .TRUE., if descaling of DX1 required, 2709C else .FALSE. 2710C 2711C ------------------------------------------------------------ 2712C* End Prologue 2713 INTRINSIC DABS 2714 DOUBLE PRECISION ZERO 2715 PARAMETER (ZERO=0.0D0) 2716 INTEGER L1 2717 DOUBLE PRECISION S1 2718C* Begin 2719 IF (QDSCAL) THEN 2720C ------------------------------------------------------------ 2721C 1.2 Descaling of solution DX1 ( stored to DXQ ) 2722 DO 12 L1=1,N 2723 DXQ(L1)=DX1(L1)*XW(L1) 272412 CONTINUE 2725 ENDIF 2726C ------------------------------------------------------------ 2727C 2 Evaluation of scaled natural level function SUMX and 2728C scaled maximum error norm CONV 2729 CONV = ZERO 2730 DO 20 L1=1,N 2731 S1 = DABS(DX1(L1)) 2732 IF(S1.GT.CONV) CONV=S1 273320 CONTINUE 2734 SUMX = ZERO 2735 DO 21 L1=1,N 2736 SUMX = SUMX+DX1(L1)**2 273721 CONTINUE 2738C ------------------------------------------------------------ 2739C 3 Evaluation of (scaled) standard level function DLEVF (only 2740C if needed for print) 2741 IF(MPRMON.GE.2)THEN 2742 DLEVF = ZERO 2743 DO 3 L1=1,N 2744 DLEVF = DLEVF+F(L1)**2 27453 CONTINUE 2746 DLEVF = DSQRT(DLEVF/DBLE(FLOAT(N))) 2747 ENDIF 2748C End of subroutine N2LVLS 2749 RETURN 2750 END 2751C 2752 SUBROUTINE N2JAC (FCN, N, LDA, X, FX, A, YSCAL, AJDEL, AJMIN, 2753 $ NFCN, FU, IFAIL) 2754C* Begin Prologue N2JAC 2755 EXTERNAL FCN 2756 INTEGER N, LDA 2757 DOUBLE PRECISION X(N), FX(N), A(LDA,N), YSCAL(N), AJDEL, AJMIN 2758 INTEGER NFCN 2759 DOUBLE PRECISION FU(N) 2760 INTEGER IFAIL 2761C 2762C --------------------------------------------------------------------- 2763C 2764C* Title 2765C 2766C Evaluation of a dense Jacobian matrix using finite difference 2767C approximation adapted for use in nonlinear systems solver NLEQ2 2768C 2769C* Environment Fortran 77 2770C Double Precision 2771C Sun 3/60, Sun OS 2772C* Latest Revision January 1991 2773C 2774C 2775C* Parameter list description 2776C -------------------------- 2777C 2778C* External subroutines (to be supplied by the user) 2779C ------------------------------------------------- 2780C 2781C FCN Ext FCN (N, X, FX, IFAIL) 2782C Subroutine in order to provide right-hand 2783C side of first-order differential equations 2784C N Int Number of rows and columns of the Jacobian 2785C X(N) Dble The current scaled iterates 2786C FX(N) Dble Array containing FCN(X) 2787C IFAIL Int Return code 2788C Whenever a negative value is returned by FCN 2789C routine N2JAC is terminated immediately. 2790C 2791C 2792C* Input parameters (* marks inout parameters) 2793C ---------------- 2794C 2795C N Int Number of rows and columns of the Jacobian 2796C LDA Int Leading Dimension of array A 2797C X(N) Dble Array containing the current scaled 2798C iterate 2799C FX(N) Dble Array containing FCN(X) 2800C YSCAL(N) Dble Array containing the scaling factors 2801C AJDEL Dble Perturbation of component k: abs(Y(k))*AJDEL 2802C AJMIN Dble Minimum perturbation is AJMIN*AJDEL 2803C NFCN Int * FCN - evaluation count 2804C 2805C* Output parameters (* marks inout parameters) 2806C ----------------- 2807C 2808C A(LDA,N) Dble Array to contain the approximated 2809C Jacobian matrix ( dF(i)/dx(j)in A(i,j)) 2810C NFCN Int * FCN - evaluation count adjusted 2811C IFAIL Int Return code non-zero if Jacobian could not 2812C be computed. 2813C 2814C* Workspace parameters 2815C -------------------- 2816C 2817C FU(N) Dble Array to contain FCN(x+dx) for evaluation of 2818C the numerator differences 2819C 2820C* Called 2821C ------ 2822C 2823 INTRINSIC DABS, DMAX1, DSIGN 2824C --------------------------------------------------------------------- 2825C 2826C* End Prologue 2827C 2828C* Local variables 2829C --------------- 2830C 2831 INTEGER I, K 2832 DOUBLE PRECISION U, W 2833C 2834C* Begin 2835C 2836 IFAIL = 0 2837 DO 1 K = 1,N 2838 W = X(K) 2839 U = DSIGN(DMAX1(DABS(X(K)),AJMIN,YSCAL(K))*AJDEL, X(K)) 2840 X(K) = W + U 2841C 2842 CALL FCN (N, X, FU, IFAIL) 2843 NFCN = NFCN + 1 2844 IF (IFAIL .NE. 0) GOTO 99 2845C 2846 X(K) = W 2847 DO 11 I = 1,N 2848 A(I,K) = (FU(I) - FX(I)) / U 2849 11 CONTINUE 2850 1 CONTINUE 2851C 285299 CONTINUE 2853 RETURN 2854C 2855C 2856C* End of N2JAC 2857C 2858 END 2859 SUBROUTINE N2JCF (FCN, N, LDA, X, FX, A, YSCAL, ETA, ETAMIN, 2860 $ ETAMAX, ETADIF, CONV, NFCN, FU, IFAIL) 2861C* Begin Prologue N2JCF 2862 EXTERNAL FCN 2863 INTEGER N, LDA 2864 DOUBLE PRECISION X(N), FX(N), A(LDA,N), YSCAL(N), ETA(N), 2865 $ ETAMIN, ETAMAX, ETADIF, CONV 2866 INTEGER NFCN 2867 DOUBLE PRECISION FU(N) 2868 INTEGER IFAIL 2869C 2870C --------------------------------------------------------------------- 2871C 2872C* Title 2873C 2874C Approximation of dense Jacobian matrix for nonlinear systems solver 2875C NLEQ2 with feed-back control of discretization and rounding errors 2876C 2877C* Environment Fortran 77 2878C Double Precision 2879C Sun 3/60, Sun OS 2880C* Latest Revision May 1990 2881C 2882C 2883C* Parameter list description 2884C -------------------------- 2885C 2886C* External subroutines (to be supplied by the user) 2887C ------------------------------------------------- 2888C 2889C FCN Ext FCN (N, X, FX, IFAIL) 2890C Subroutine in order to provide right-hand 2891C side of first-order differential equations 2892C N Int Number of rows and columns of the Jacobian 2893C X(N) Dble The current scaled iterates 2894C FX(N) Dble Array containing FCN(X) 2895C IFAIL Int Return code 2896C Whenever a negative value is returned by FCN 2897C routine N2JCF is terminated immediately. 2898C 2899C 2900C* Input parameters (* marks inout parameters) 2901C ---------------- 2902C 2903C N Int Number of rows and columns of the Jacobian 2904C LDA Int Leading dimension of A (LDA .GE. N) 2905C X(N) Dble Array containing the current scaled 2906C iterate 2907C FX(N) Dble Array containing FCN(X) 2908C YSCAL(N) Dble Array containing the scaling factors 2909C ETA(N) Dble * Array containing the scaled denominator 2910C differences 2911C ETAMIN Dble Minimum allowed scaled denominator 2912C ETAMAX Dble Maximum allowed scaled denominator 2913C ETADIF Dble DSQRT (1.1*EPMACH) 2914C EPMACH = machine precision 2915C CONV Dble Maximum norm of last (unrelaxed) Newton correction 2916C NFCN Int * FCN - evaluation count 2917C 2918C* Output parameters (* marks inout parameters) 2919C ----------------- 2920C 2921C A(LDA,N) Dble Array to contain the approximated 2922C Jacobian matrix ( dF(i)/dx(j)in A(i,j)) 2923C ETA(N) Dble * Scaled denominator differences adjusted 2924C NFCN Int * FCN - evaluation count adjusted 2925C IFAIL Int Return code non-zero if Jacobian could not 2926C be computed. 2927C 2928C* Workspace parameters 2929C -------------------- 2930C 2931C FU(N) Dble Array to contain FCN(x+dx) for evaluation of 2932C the numerator differences 2933C 2934C* Called 2935C ------ 2936C 2937 INTRINSIC DABS, DMAX1, DMIN1, DSIGN, DSQRT 2938C 2939C* Constants 2940C --------- 2941C 2942 DOUBLE PRECISION SMALL2, ZERO 2943 PARAMETER (SMALL2 = 0.1D0, 2944 $ ZERO = 0.0D0) 2945C 2946C --------------------------------------------------------------------- 2947C 2948C* End Prologue 2949C 2950C* Local variables 2951C --------------- 2952C 2953 INTEGER I, K, IS 2954 DOUBLE PRECISION FHI, HG, U, SUMD, W 2955 LOGICAL QFINE 2956C 2957C* Begin 2958C 2959 DO 1 K = 1,N 2960 IS = 0 2961C DO (Until) 2962 11 CONTINUE 2963 W = X(K) 2964 U = DSIGN (ETA(K)*YSCAL(K), X(K)) 2965 X(K) = W + U 2966 CALL FCN (N, X, FU, IFAIL) 2967 NFCN = NFCN + 1 2968C Exit, If ... 2969 IF (IFAIL .NE. 0) GOTO 99 2970 X(K) = W 2971 SUMD = ZERO 2972 DO 111 I = 1,N 2973 HG = DMAX1 (DABS (FX(I)), DABS (FU(I))) 2974 FHI = FU(I) - FX(I) 2975 IF (HG .NE. ZERO) SUMD = SUMD + (FHI/HG)**2 2976 A(I,K) = FHI / U 2977 111 CONTINUE 2978 SUMD = DSQRT (SUMD / DBLE(N)) 2979 QFINE = .TRUE. 2980 IF (SUMD .NE. ZERO .AND. IS .EQ. 0)THEN 2981 ETA(K) = DMIN1 (ETAMAX, 2982 $ DMAX1 (ETAMIN, DSQRT (ETADIF / SUMD)*ETA(K))) 2983 IS = 1 2984 QFINE = CONV .LT. SMALL2 .OR. SUMD .GE. ETAMIN 2985 ENDIF 2986 IF (.NOT.(QFINE)) GOTO 11 2987C UNTIL ( expression - negated above) 2988 1 CONTINUE 2989C 2990C Exit from DO-loop 2991 99 CONTINUE 2992C 2993 RETURN 2994C 2995C* End of subroutine N2JCF 2996C 2997 END 2998 SUBROUTINE N2PRJN(N,IRANK,DEL,U,D,V,QE,PIVOT) 2999C* Begin Prologue PRJCTN 3000 INTEGER IRANK,N 3001 INTEGER PIVOT(N) 3002 DOUBLE PRECISION DEL 3003 DOUBLE PRECISION QE(N,N) 3004 DOUBLE PRECISION U(N),D(N),V(N) 3005C ------------------------------------------------------------ 3006C 3007C* Summary : 3008C 3009C P R J C T N : 3010C To be used in connection with either DECOMP/SOLVE or 3011C DECCON/SOLCON . 3012C Provides the projection to the appropriate subspace in case 3013C of rank - reduction 3014C 3015C* Input parameters (* marks inout parameters) 3016C =========================================== 3017C 3018C N Int Number of parameters to be estimated 3019C IRANK Pseudo rank of decomposed Jacobian 3020C matrix 3021C U(N) Dble Scaled Newton correction 3022C D(N) Dble Diagonal elements of upper 3023C triangular matrix 3024C QE(N,N) Dble Part of pseudoinverse Jacobian 3025C matrix ( see QA of DECCON ) 3026C PIVOT(N) Dble Pivot vector resulting from matrix 3027C decomposition (DECCON) 3028C V(N) Dble Real work array 3029C 3030C* Output parameters 3031C ================= 3032C 3033C DEL Dble Defekt 3034C 3035C ------------------------------------------------------------ 3036C* End Prologue 3037 DOUBLE PRECISION ZERO 3038 PARAMETER (ZERO=0.0D0) 3039 INTEGER L1,I,IRK1 3040 DOUBLE PRECISION S,SH 3041C* Begin 3042 DO 1 I=1,N 3043 V(I)=U(PIVOT(I)) 30441 CONTINUE 3045 IRK1 = IRANK+1 3046 DEL = ZERO 3047 DO 2 I=IRK1,N 3048 SH = 0.0 3049 DO 21 L1=1,I-1 3050 SH = SH+QE(L1,I)*V(L1) 305121 CONTINUE 3052 S =(V(I)-SH)/D(I) 3053 DEL = S*S+DEL 3054 V(I)=S 30552 CONTINUE 3056C End of subroutine N2PRJN 3057 RETURN 3058 END 3059C 3060 SUBROUTINE N2PRV1(DLEVF,DLEVX,FC,NITER,NEW,IRANK,MPRMON,LUMON, 3061 $ QMIXIO,COND1) 3062C* Begin Prologue N2PRV1 3063 DOUBLE PRECISION DLEVF,DLEVX,FC,COND1 3064 INTEGER NITER,NEW,IRANK,MPRMON,LUMON 3065 LOGICAL QMIXIO 3066C ------------------------------------------------------------ 3067C 3068C* Summary : 3069C 3070C N 2 P R V 1 : Printing of intermediate values (Type 1 routine) 3071C 3072C Parameters 3073C ========== 3074C 3075C DLEVF, DLEVX See descr. of internal double variables of N2INT 3076C FC,NITER,NEW,IRANK,MPRMON,LUMON,COND1 3077C See parameter descr. of subroutine N2INT 3078C QMIXIO Logical = .TRUE. , if LUMON.EQ.LUSOL 3079C = .FALSE. , if LUMON.NE.LUSOL 3080C 3081C ------------------------------------------------------------ 3082C* End Prologue 3083C Print Standard - and natural level 3084 IF(QMIXIO)THEN 30851 FORMAT(2X,66('*')) 3086 WRITE(LUMON,1) 30872 FORMAT(8X,'It',7X,'Normf ',10X,'Normx ',20X,'New',6X,'Rank', 3088 $ 8X,'Cond') 3089 IF (MPRMON.GE.3) WRITE(LUMON,2) 30903 FORMAT(8X,'It',7X,'Normf ',10X,'Normx ',8X,'Damp.Fct.',3X,'New', 3091 $ 6X,'Rank',8X,'Cond') 3092 IF (MPRMON.EQ.2) WRITE(LUMON,3) 3093 ENDIF 30944 FORMAT(6X,I4,5X,D10.3,2X,4X,D10.3,17X,I2,6X,I4,2X,D10.3) 3095 IF (MPRMON.GE.3.OR.NITER.EQ.0) 3096 $ WRITE(LUMON,4) NITER,DLEVF,DLEVX,NEW,IRANK,COND1 30975 FORMAT(6X,I4,5X,D10.3,6X,D10.3,6X,F7.5,4X,I2,6X,I4,2X,D10.3) 3098 IF (MPRMON.EQ.2.AND.NITER.NE.0) 3099 $ WRITE(LUMON,5) NITER,DLEVF,DLEVX,FC,NEW,IRANK,COND1 3100 IF(QMIXIO)THEN 31016 FORMAT(2X,66('*')) 3102 WRITE(LUMON,6) 3103 ENDIF 3104C End of subroutine N2PRV1 3105 RETURN 3106 END 3107C 3108 SUBROUTINE N2PRV2(DLEVF,DLEVX,FC,NITER,MPRMON,LUMON,QMIXIO, 3109 $ CMARK) 3110C* Begin Prologue N2PRV2 3111 DOUBLE PRECISION DLEVF,DLEVX,FC 3112 INTEGER NITER,MPRMON,LUMON 3113 LOGICAL QMIXIO 3114 CHARACTER*1 CMARK 3115C ------------------------------------------------------------ 3116C 3117C* Summary : 3118C 3119C N 2 P R V 2 : Printing of intermediate values (Type 2 routine) 3120C 3121C* Parameters 3122C ========== 3123C 3124C DLEVF, DLEVX See descr. of internal double variables of N2INT 3125C FC,NITER,MPRMON,LUMON 3126C See parameter descr. of subroutine N2INT 3127C QMIXIO Logical = .TRUE. , if LUMON.EQ.LUSOL 3128C = .FALSE. , if LUMON.NE.LUSOL 3129C CMARK Char*1 Marker character to be printed before DLEVX 3130C 3131C ------------------------------------------------------------ 3132C* End Prologue 3133C Print Standard - and natural level, and damping 3134C factor 3135 IF(QMIXIO)THEN 31361 FORMAT(2X,66('*')) 3137 WRITE(LUMON,1) 31382 FORMAT(8X,'It',7X,'Normf ',10X,'Normx ',8X,'Damp.Fct.') 3139 WRITE(LUMON,2) 3140 ENDIF 31413 FORMAT(6X,I4,5X,D10.3,4X,A1,1X,D10.3,6X,F7.5) 3142 WRITE(LUMON,3)NITER,DLEVF,CMARK,DLEVX,FC 3143 IF(QMIXIO)THEN 31444 FORMAT(2X,66('*')) 3145 WRITE(LUMON,4) 3146 ENDIF 3147C End of subroutine N2PRV2 3148 RETURN 3149 END 3150C 3151 SUBROUTINE N2SOUT(N,X,MODE,IOPT,RWK,NRW,IWK,NIW,MPRINT,LUOUT) 3152C* Begin Prologue SOLOUT 3153 INTEGER N 3154 DOUBLE PRECISION X(N) 3155 INTEGER NRW 3156 INTEGER MODE 3157 INTEGER IOPT(50) 3158 DOUBLE PRECISION RWK(NRW) 3159 INTEGER NIW 3160 INTEGER IWK(NIW) 3161 INTEGER MPRINT,LUOUT 3162C ------------------------------------------------------------ 3163C 3164C* Summary : 3165C 3166C S O L O U T : Printing of iterate (user customizable routine) 3167C 3168C* Input parameters 3169C ================ 3170C 3171C N Int Number of equations/unknowns 3172C X(N) Dble iterate vector 3173C MODE =1 This routine is called before the first 3174C Newton iteration step 3175C =2 This routine is called with an intermedi- 3176C ate iterate X(N) 3177C =3 This is the last call with the solution 3178C vector X(N) 3179C =4 This is the last call with the final, but 3180C not solution vector X(N) 3181C IOPT(50) Int The option array as passed to the driver 3182C routine (elements 46 to 50 may be used 3183C for user options) 3184C MPRINT Int Solution print level 3185C (see description of IOPT-field MPRINT) 3186C LUOUT Int the solution print unit 3187C (see description of see IOPT-field LUSOL) 3188C 3189C 3190C* Workspace parameters 3191C ==================== 3192C 3193C NRW, RWK, NIW, IWK see description in driver routine 3194C 3195C* Use of IOPT by this routine 3196C =========================== 3197C 3198C Field 46: =0 Standard output 3199C =1 GRAZIL suitable output 3200C 3201C ------------------------------------------------------------ 3202C* End Prologue 3203 LOGICAL QGRAZ,QNORM 3204C* Begin 3205 QNORM = IOPT(46).EQ.0 3206 QGRAZ = IOPT(46).EQ.1 3207 IF(QNORM) THEN 32081 FORMAT(' ',A,' data:',/) 3209 IF (MODE.EQ.1) THEN 3210101 FORMAT(' Start data:',/,' N =',I5,//, 3211 $ ' Format: iteration-number, (x(i),i=1,...N), ', 3212 $ 'Normf , Normx ',/) 3213 WRITE(LUOUT,101) N 3214 WRITE(LUOUT,1) 'Initial' 3215 ELSE IF (MODE.EQ.3) THEN 3216 WRITE(LUOUT,1) 'Solution' 3217 ELSE IF (MODE.EQ.4) THEN 3218 WRITE(LUOUT,1) 'Final' 3219 ENDIF 32202 FORMAT(' ',I5) 3221C WRITE NITER 3222 WRITE(LUOUT,2) IWK(1) 32233 FORMAT((12X,3(D18.10,1X))) 3224 WRITE(LUOUT,3)(X(L1),L1=1,N) 3225C WRITE DLEVF, DLEVX 3226 WRITE(LUOUT,3) RWK(19),DSQRT(RWK(18)/DBLE(FLOAT(N))) 3227 IF(MODE.EQ.1.AND.MPRINT.GE.2) THEN 3228 WRITE(LUOUT,1) 'Intermediate' 3229 ELSE IF(MODE.GE.3) THEN 3230 WRITE(LUOUT,1) 'End' 3231 ENDIF 3232 ENDIF 3233 IF(QGRAZ) THEN 3234 IF(MODE.EQ.1) THEN 323510 FORMAT('&name com',I3.3,:,255(7(', com',I3.3,:),/)) 3236 WRITE(LUOUT,10)(I,I=1,N+2) 323715 FORMAT('&def com',I3.3,:,255(7(', com',I3.3,:),/)) 3238 WRITE(LUOUT,15)(I,I=1,N+2) 323916 FORMAT(6X,': X=1, Y=',I3) 3240 WRITE(LUOUT,16) N+2 3241 ENDIF 324220 FORMAT('&data ',I5) 3243C WRITE NITER 3244 WRITE(LUOUT,20) IWK(1) 324521 FORMAT((6X,4(D18.10))) 3246 WRITE(LUOUT,21)(X(L1),L1=1,N) 3247C WRITE DLEVF, DLEVX 3248 WRITE(LUOUT,21) RWK(19),DSQRT(RWK(18)/DBLE(FLOAT(N))) 3249 IF(MODE.GE.3) THEN 325030 FORMAT('&wktype 3111',/,'&atext x ''iter''') 3251 WRITE(LUOUT,30) 325235 FORMAT('&vars = com',I3.3,/,'&atext y ''x',I3,'''', 3253 $ /,'&run') 3254 WRITE(LUOUT,35) (I,I,I=1,N) 325536 FORMAT('&vars = com',I3.3,/,'&atext y ''',A,'''', 3256 $ /,'&run') 3257 WRITE(LUOUT,36) N+1,'Normf ',N+2,'Normx ' 3258C39 FORMAT('&stop') 3259C WRITE(LUOUT,39) 3260 ENDIF 3261 ENDIF 3262C End of subroutine N2SOUT 3263 RETURN 3264 END 3265C 3266 DOUBLE PRECISION FUNCTION WNORM(N,Z,XW) 3267 INTEGER N 3268 DOUBLE PRECISION Z(N), XW(N) 3269C ------------------------------------------------------------ 3270C 3271C* Summary : 3272C 3273C E N O R M : Return the norm to be used in exit (termination) 3274C criteria 3275C 3276C* Input parameters 3277C ================ 3278C 3279C N Int Number of equations/unknowns 3280C Z(N) Dble The vector, of which the norm is to be computed 3281C XW(N) Dble The scaling values of Z(N) 3282C 3283C* Output 3284C ====== 3285C 3286C WNORM(N,Z,XW) Dble The mean square root norm of Z(N) subject 3287C to the scaling values in XW(N): 3288C = Sqrt( Sum(1,...N)((Z(I)/XW(I))**2) / N ) 3289C 3290C ------------------------------------------------------------ 3291C* End Prologue 3292 INTEGER I 3293 DOUBLE PRECISION S 3294C* Begin 3295 S = 0.0D0 3296 DO 10 I=1,N 3297 S = S + ( Z(I)/XW(I) ) ** 2 329810 CONTINUE 3299 WNORM = DSQRT( S / DBLE(FLOAT(N)) ) 3300C End of function WNORM 3301 RETURN 3302 END 3303C 3304C 3305C* Group Linear Solver subroutines (Code DECCON/SOLCON) 3306C 3307 SUBROUTINE DECCON(A,NROW,NCOL,MCON,M,N,IRANKC,IRANK,COND,D,PIVOT, 3308 *KRED,AH,V,IERR) 3309C* Begin Prologue DECCON 3310 INTEGER IRANKC,IRANK,MCON 3311 INTEGER M,N,NROW,NCOL,KRED 3312 INTEGER PIVOT(NCOL) 3313 DOUBLE PRECISION COND 3314 DOUBLE PRECISION A(NROW,NCOL),AH(NCOL,NCOL) 3315 DOUBLE PRECISION D(NCOL),V(NCOL) 3316 INTEGER IERR 3317C ------------------------------------------------------------ 3318C 3319C* Title 3320C 3321C* Deccon - Constrained Least Squares QR-Decomposition 3322C 3323C* Written by P. Deuflhard, U. Nowak, L. Weimann 3324C* Purpose Solution of least squares problems, optionally 3325C with equality constraints. 3326C* Method Constrained Least Squares QR-Decomposition 3327C (see references below) 3328C* Category D9b1. - Singular, overdetermined or 3329C underdetermined systems of linear 3330C equations, generalized inverses. 3331C Constrained Least Squares solution 3332C* Keywords Linear Least Square Problems, constrained, 3333C QR-decomposition, pseudo inverse. 3334C* Version 1.1 3335C* Revision March 1989 3336C* Latest Change January 1991 3337C* Library CodeLib 3338C* Code Fortran 77, Double Precision 3339C* Environment Standard Fortran 77 environment on PC's, 3340C workstations and hosts. 3341C* Copyright (c) Konrad Zuse Zentrum fuer 3342C Informationstechnik Berlin 3343C Heilbronner Str. 10, D-1000 Berlin 31 3344C phone 0049+30+89604-0, 3345C telefax 0049+30+89604-125 3346C* Contact Lutz Weimann 3347C ZIB, Numerical Software Development 3348C phone: 0049+30+89604-185 ; 3349C e-mail: 3350C RFC822 notation: weimann@sc.zib-berlin.de 3351C X.400: C=de;A=dbp;P=zib-berlin;OU=sc;S=Weimann 3352C 3353C* References: 3354C =========== 3355C 3356C /1/ P.Deuflhard, V.Apostolescu: 3357C An underrelaxed Gauss-Newton method for equality 3358C constrained nonlinear least squares problems. 3359C Lecture Notes Control Inform. Sci. vol. 7, p. 3360C 22-32 (1978) 3361C /2/ P.Deuflhard, W.Sautter: 3362C On rank-deficient pseudoinverses. 3363C J. Lin. Alg. Appl. vol. 29, p. 91-111 (1980) 3364C 3365C* Related Programs: SOLCON 3366C 3367C --------------------------------------------------------------- 3368C 3369C* Licence 3370C You may use or modify this code for your own non commercial 3371C purposes for an unlimited time. 3372C In any case you should not deliver this code without a special 3373C permission of ZIB. 3374C In case you intend to use the code commercially, we oblige you 3375C to sign an according licence agreement with ZIB. 3376C 3377C* Warranty 3378C This code has been tested up to a certain level. Defects and 3379C weaknesses, which may be included in the code, do not establish 3380C any warranties by ZIB. ZIB does not take over any liabilities 3381C which may follow from acquisition or application of this code. 3382C 3383C* Software status 3384C This code is under care of ZIB and belongs to ZIB software class 1. 3385C 3386C ------------------------------------------------------------ 3387C 3388C* Summary: 3389C ======== 3390C Constrained QR-decomposition of (M,N)-system with 3391C computation of pseudoinverse in case of rank-defeciency . 3392C First MCON rows belong to equality constraints. 3393C 3394C ------------------------------------------------------------ 3395C 3396C* Parameters list description (* marks inout parameters) 3397C ====================================================== 3398C 3399C* Input parameters 3400C ================ 3401C 3402C A(NROW,NCOL) Dble Array holding the (M,N)-Matrix to be 3403C decomposed 3404C NROW Int Declared number of rows of array A 3405C NCOL Int Declared number of columns of array A and 3406C rows and columns of array AH 3407C MCON Int Number of equality constraints (MCON.LE.N) 3408C Internally reduced if equality constraints 3409C are linearly dependent 3410C M Int Current number of rows of matrix A 3411C N Int Current number of columns of matrix A 3412C * IRANKC Int Prescribed maximum pseudo-rank of 3413C constrained part of matrix A (IRANKC.LE.MCON) 3414C * IRANK Int Prescribed maximum pseudo-rank of matrix A 3415C (IRANK.LE.N) 3416C * COND Dble Permitted upper bound of DABS(D(1)/D(IRANK)) 3417C and of DABS(D(IRANK+1)/D(IRANK)) 3418C (subcondition numbers of A) 3419C KRED Int Type of operation 3420C >=0 Householder triangularization 3421C (build up pseudo-inverse,if IRANK.LT.N) 3422C < 0 Reduction of pseudo-rank of matrix A, 3423C skipping Householder triangularization, of 3424C build-up new pseudo-inverse 3425C 3426C* Output parameters 3427C ================= 3428C 3429C A(NROW,NCOL) Dble Array holding the (M,N)-output consisting 3430C of the transformed matrix in the upper 3431C right triangle and the performed House- 3432C holder transf. in the lower left triangle. 3433C * IRANKC Int New pseudo-rank of constrained part of 3434C matrix A, determined so that 3435C DABS(D(1)/D(IRANKC))<COND 3436C * IRANK Int New pseudo-rank of matrix A, determined 3437C so that DABS(D(IRANKC+1)/D(IRANK)) < COND 3438C D(IRANK) Dble Diagonal elements of upper triangular matr. 3439C PIVOT(N) Int Index vector storing permutation of columns 3440C due to pivoting 3441C * COND Dble The maximum of the two sub-condition 3442C numbers belonging to the constrained 3443C part and the remaining part of A. 3444C (in case of rank reduction: 3445C sub-condition number which led to 3446C rank reduction) 3447C AH(NCOL,NCOL) Dble In case of rank-defect used to compute the 3448C pseudo-inverse (currently used will be an 3449C (N,N)-part of this array) 3450C IERR Int Error indicator: 3451C = 0 : DECCON computations are successfull. 3452C =-2 : Numerically negative diagonal element 3453C encountered during computation of 3454C pseudo inverse - due to extremely bad 3455C conditioned Matrix A. DECCON is 3456C unable to continue rank-reduction. 3457C 3458C* Workspace parameters 3459C ==================== 3460C 3461C V(N) Dble Workspace array 3462C 3463C* Subroutines called: D1MACH 3464C 3465C* Machine constants used 3466C ====================== 3467C 3468C EPMACH = relative machine precision 3469 DOUBLE PRECISION EPMACH 3470C 3471C ------------------------------------------------------------ 3472C* End Prologue 3473 EXTERNAL D1MACH 3474 INTRINSIC DABS,DSQRT 3475 DOUBLE PRECISION ZERO 3476 PARAMETER (ZERO=0.0D0) 3477 DOUBLE PRECISION ONE 3478 PARAMETER (ONE=1.0D0) 3479 DOUBLE PRECISION REDUCE 3480 PARAMETER (REDUCE=0.05D0) 3481 INTEGER L1 3482 DOUBLE PRECISION S1 3483 INTEGER I,II,IRANKH,IRK1,I1,J,JD,JJ,K,K1,LEVEL,MH 3484 DOUBLE PRECISION DD,D1MACH,H,HMAX,S,SH,SMALLS,T 3485C* Begin 3486C -------------------------------------------------------------- 3487C 1 Initialization 3488 EPMACH = D1MACH(3) 3489 SMALLS = DSQRT(EPMACH*10.0D0) 3490 IF(IRANK.GT.N) IRANK = N 3491 IF(IRANK.GT.M) IRANK = M 3492C -------------------------------------------------------------- 3493C 1.1 Special case M=1 and N=1 3494 IF(M.EQ.1.AND.N.EQ.1)THEN 3495 PIVOT(1)=1 3496 D(1)=A(1,1) 3497 COND = ONE 3498 RETURN 3499 ENDIF 3500 IF(KRED.GE.0)THEN 3501C ------------------------------------------------------------ 3502C 1.1 Initialize pivot-array 3503 DO 11 J=1,N 3504 PIVOT(J)=J 350511 CONTINUE 3506C ------------------------------------------------------------ 3507C 2. Constrained Householder triangularization 3508 JD = 1 3509 IRANC1 = IRANKC + 1 3510 MH = MCON 3511 IRANKH = IRANKC 3512 IF(MH.EQ.0) THEN 3513 IRANKH = IRANK 3514 MH = M 3515 ENDIF 3516 IRK1 = IRANK 3517 DO 2 K=1,IRK1 35182000 LEVEL = 1 3519 IF(K.NE.N)THEN 3520 K1 = K+1 3521C DO (Until) 352220 CONTINUE 3523 IF(JD.NE.0)THEN 3524 DO 201 J=K,N 3525 S = ZERO 3526 DO 2011 L1=K,MH 3527 S = S+A(L1,J)**2 35282011 CONTINUE 3529 D(J)=S 3530201 CONTINUE 3531 ENDIF 3532C ------------------------------------------------------ 3533C 2.1 Column pivoting 3534 S1 = D(K) 3535 JJ = K 3536 DO 21 L1=K,N 3537 IF(D(L1).GT.S1) THEN 3538 S1=D(L1) 3539 JJ = L1 3540 ENDIF 354121 CONTINUE 3542 H = D(JJ) 3543 IF(JD.EQ.1) HMAX = H/(COND* REDUCE) 3544 JD = 0 3545 IF(H.LT.HMAX) JD = 1 3546 IF(.NOT.(H.GE.HMAX)) GOTO 20 3547C UNTIL ( expression - negated above) 3548 IF(JJ.NE.K)THEN 3549C ------------------------------------------------------ 3550C 2.2 Column interchange 3551 I = PIVOT(K) 3552 PIVOT(K)=PIVOT(JJ) 3553 PIVOT(JJ)=I 3554 D(JJ)=D(K) 3555 DO 221 L1=1,M 3556 S1=A(L1,JJ) 3557 A(L1,JJ)=A(L1,K) 3558 A(L1,K)=S1 3559221 CONTINUE 3560 ENDIF 3561 ENDIF 3562 H = ZERO 3563 DO 222 L1=K,MH 3564 H = H+A(L1,K)**2 3565222 CONTINUE 3566 T = DSQRT(H) 3567C ---------------------------------------------------------- 3568C 2.3.0 A-priori test on pseudo-rank 3569 IF ( K.EQ.1 .OR. K.EQ.IRANC1 ) DD = T/COND 3570 IF(T.LE.DD .OR. K.GT.IRANKH)THEN 3571C ------------------------------------------------------ 3572C 2.3.1 Rank reduction 3573 IRANKH = K-1 3574 IF (MH.NE.MCON) THEN 3575 IRANK = IRANKH 3576 IF (IRANKC.EQ.IRANK) THEN 3577 LEVEL = 4 3578 ELSE 3579 LEVEL = 3 3580 ENDIF 3581 ELSE 3582 IRANKC = IRANKH 3583 IF (IRANKC.NE.MCON) THEN 3584 MH = M 3585 IRANKH = IRANK 3586 JD = 1 3587 GOTO 2000 3588 ENDIF 3589 ENDIF 3590 ENDIF 3591 IF (LEVEL.EQ.1) THEN 3592C ------------------------------------------------------ 3593C 2.4 Householder step 3594 S = A(K,K) 3595 T = -DSIGN(T,S) 3596 D(K)=T 3597 A(K,K)=S-T 3598 IF(K.NE.N)THEN 3599 T = ONE/(H-S*T) 3600 DO 24 J=K1,N 3601 S = ZERO 3602 DO 241 L1=K,MH 3603 S = S+A(L1,K)*A(L1,J) 3604241 CONTINUE 3605 S = S*T 3606 S1 =-S 3607 DO 242 L1=K,M 3608 A(L1,J) = A(L1,J)+A(L1,K)*S1 3609242 CONTINUE 3610 D(J)=D(J)-A(K,J)**2 361124 CONTINUE 3612 IF(K.EQ.IRANKC)THEN 3613 MH = M 3614 JD = 1 3615 IRANKH=IRANK 3616 IF (K.EQ.IRK1) LEVEL=3 3617 ENDIF 3618 ELSE 3619 LEVEL = 4 3620 ENDIF 3621 ENDIF 3622C Exit Do 2 If ... 3623 IF(LEVEL.GT.1) GOTO 2999 36242 CONTINUE 3625C ENDDO 36262999 CONTINUE 3627 ELSE 3628 K = -1 3629 LEVEL = 3 3630 ENDIF 3631C -------------------------------------------------------------- 3632C 3 Rank-deficient pseudo-inverse 3633 IF(LEVEL.EQ.3)THEN 3634 IRK1 = IRANK+1 3635 DO 3 J=IRK1,N 3636 DO 31 II=1,IRANK 3637 I = IRK1-II 3638 S = A(I,J) 3639 IF(II.NE.1)THEN 3640 SH = ZERO 3641 DO 311 L1=I1,IRANK 3642 SH=SH+A(I,L1)*V(L1) 3643311 CONTINUE 3644 S = S-SH 3645 ENDIF 3646 I1 = I 3647 V(I)=S/D(I) 3648 AH(I,J)=V(I) 364931 CONTINUE 3650 DO 32 I=IRK1,J 3651 S = ZERO 3652 DO 321 L1=1,I-1 3653 S = S+AH(L1,I)*V(L1) 3654321 CONTINUE 3655 IF(I.NE.J)THEN 3656 V(I)=-S/D(I) 3657 AH(I,J)=-V(I) 3658 ENDIF 365932 CONTINUE 3660 IF (S.GT.-ONE) THEN 3661 D(J)=DSQRT(S+ONE) 3662 ELSE 3663 IERR=-2 3664 GOTO 999 3665 ENDIF 36663 CONTINUE 3667 ENDIF 3668C -------------------------------------------------------------- 3669C 9 Exit 3670 IF (IRANKC.NE.0) THEN 3671 SH = D(IRANKC) 3672 IF(SH.NE.ZERO) SH = DABS(D(1)/SH) 3673 ELSE 3674 SH=ZERO 3675 ENDIF 3676 IF (K.EQ.IRANK) T = D(IRANK) 3677 IF (IRANKC+1.LE.IRANK .AND. T.NE.ZERO) THEN 3678 S = DABS(D(IRANKC+1)/T) 3679 ELSE 3680 S = ZERO 3681 ENDIF 3682 IF (S.GT.SH) THEN 3683 COND=S 3684 ELSE 3685 COND=SH 3686 ENDIF 3687 IERR=0 3688999 RETURN 3689 END 3690 SUBROUTINE SOLCON(A,NROW,NCOL,MCON,M,N,X,B,IRANKC,IRANK,D,PIVOT, 3691 *KRED,AH,V) 3692C* Begin Prologue SOLCON 3693 DOUBLE PRECISION A(NROW,NCOL),AH(NCOL,NCOL) 3694 DOUBLE PRECISION X(NCOL),B(NROW),D(NCOL),V(NCOL) 3695 INTEGER NROW,NCOL,MCON,M,N,IRANKC,IRANK,KRED 3696 INTEGER PIVOT(NCOL) 3697C ------------------------------------------------------------ 3698C 3699C* Summary 3700C ======= 3701C 3702C Best constrained linear least squares solution of (M,N)- 3703C system . First MCON rows comprise MCON equality constraints. 3704C To be used in connection with subroutine DECCON 3705C References: See DECCON 3706C Related Programs: DECCON 3707C 3708C* Parameters: 3709C =========== 3710C 3711C* Input parameters (* marks inout parameters) 3712C =========================================== 3713C 3714C A(M,N), NROW, NCOL, M, N, MCON, IRANKC, IRANK, 3715C D(N), PIVOT(N), AH(N,N), KRED 3716C See input- respective output-parameters 3717C description of subroutine DECCON 3718C * B(M) Dble Right-hand side of linear system, if 3719C KRED.GE.0 3720C Right-hand side of upper linear system, 3721C if KRED.LT.0 3722C 3723C* Output parameters 3724C ================= 3725C 3726C X(N) Dble Best LSQ-solution of linear system 3727C B(M) Dble Right-hand of upper trigular system 3728C (transformed right-hand side of linear 3729C system) 3730C 3731C* Workspace parameters 3732C ==================== 3733C 3734C V(N) Dble Workspace array 3735C 3736C ------------------------------------------------------------ 3737C* End Prologue 3738 DOUBLE PRECISION ZERO 3739 PARAMETER (ZERO=0.0D0) 3740 INTEGER L1,L2 3741 INTEGER I,II,I1,IRANC1,IRK1,J,JJ,J1,MH 3742 DOUBLE PRECISION S,SH 3743C* Begin 3744C ------------------------------------------------------------ 3745C 1 Solution for pseudo-rank zero 3746 IF(IRANK.EQ.0)THEN 3747 DO 11 L1=1,N 3748 X(L1)=ZERO 374911 CONTINUE 3750 RETURN 3751 ENDIF 3752 IF ((IRANK.LE.IRANKC).AND.(IRANK.NE.N)) THEN 3753 IRANC1 = IRANKC + 1 3754 DO 12 L1=IRANC1,N 3755 V(L1) = ZERO 375612 CONTINUE 3757 ENDIF 3758 IF(KRED.GE.0.AND.(M.NE.1.OR.N.NE.1))THEN 3759C ---------------------------------------------------------- 3760C 2 Constrained householder transformations of right-hand side 3761 MH = MCON 3762 IF(IRANKC.EQ.0) MH = M 3763 DO 21 J=1,IRANK 3764 S = ZERO 3765 DO 211 L1=J,MH 3766 S = S+A(L1,J)*B(L1) 3767211 CONTINUE 3768 S = S/(D(J)*A(J,J)) 3769 DO 212 L1=J,M 3770 B(L1)=B(L1)+A(L1,J)*S 3771212 CONTINUE 3772 IF(J.EQ.IRANKC) MH = M 377321 CONTINUE 3774 ENDIF 3775C ------------------------------------------------------------ 3776C 3 Solution of upper triangular system 3777 IRK1 = IRANK+1 3778 DO 31 II=1,IRANK 3779 I = IRK1-II 3780 I1 = I+1 3781 S = B(I) 3782 IF(II.NE.1)THEN 3783 SH = ZERO 3784 DO 311 L1=I1,IRANK 3785 SH=SH+A(I,L1)*V(L1) 3786311 CONTINUE 3787 S = S-SH 3788 ENDIF 3789 V(I)=S/D(I) 379031 CONTINUE 3791 IF((IRANK.NE.N).AND.(IRANK.NE.IRANKC))THEN 3792C ---------------------------------------------------------- 3793C 3.2 Computation of the best constrained least squares- 3794C solution 3795 DO 321 J=IRK1,N 3796 S = ZERO 3797 DO 3211 L1=1,J-1 3798 S = S+AH(L1,J)*V(L1) 37993211 CONTINUE 3800 V(J)=-S/D(J) 3801321 CONTINUE 3802 DO 322 JJ=1,N 3803 J = N-JJ+1 3804 S = ZERO 3805 IF(JJ.NE.1)THEN 3806 DO 3221 L1=J1,N 3807 S=S+AH(J,L1)*V(L1) 38083221 CONTINUE 3809 ENDIF 3810 IF(JJ.NE.1.AND.J.LE.IRANK)THEN 3811 V(J)=V(J)-S 3812 ELSE 3813 J1 = J 3814 V(J)=-(S+V(J))/D(J) 3815 ENDIF 3816322 CONTINUE 3817 ENDIF 3818C ------------------------------------------------------------ 3819C 4 Back-permutation of solution components 3820 DO 4 L1=1,N 3821 L2 = PIVOT(L1) 3822 X(L2) = V(L1) 38234 CONTINUE 3824 RETURN 3825 END 3826C* End package 3827C 3828C 3829C* Group Time monitor package 3830C 3831C* Begin Prologue 3832C ------------------------------------------------------------ 3833C 3834C* Title 3835C 3836C Monitor - A package for making multiple time measurements and 3837C summary statistics 3838C 3839C* Written by U. Nowak, L. Weimann 3840C* Version 1.0 3841C* Revision January 1991 3842C* Latest Change January 1991 3843C* Library CodeLib 3844C* Code Fortran 77, Double Precision 3845C* Environment Standard Fortran 77 environment on PC's, 3846C* Copyright (c) Konrad Zuse Zentrum fuer 3847C Informationstechnik Berlin 3848C Heilbronner Str. 10, D-1000 Berlin 31 3849C phone 0049+30+89604-0, 3850C telefax 0049+30+89604-125 3851C* Contact Lutz Weimann 3852C ZIB, Numerical Software Development 3853C phone: 0049+30+89604-185 ; 3854C e-mail: 3855C RFC822 notation: weimann@sc.zib-berlin.de 3856C X.400: C=de;A=dbp;P=zib-berlin;OU=sc;S=Weimann 3857C 3858C --------------------------------------------------------------- 3859C 3860C* Licence 3861C You may use or modify this code for your own non commercial 3862C purposes for an unlimited time. 3863C In any case you should not deliver this code without a special 3864C permission of ZIB. 3865C In case you intend to use the code commercially, we oblige you 3866C to sign an according licence agreement with ZIB. 3867C 3868C* Warranty 3869C This code has been tested up to a certain level. Defects and 3870C weaknesses, which may be included in the code, do not establish 3871C any warranties by ZIB. ZIB does not take over any liabilities 3872C which may follow from aquisition or application of this code. 3873C 3874C* Software status 3875C This code is under care of ZIB and belongs to ZIB software class 1. 3876C 3877C --------------------------------------------------------------- 3878C 3879C* Summary: 3880C 3881C Monitor is a package for generating time and summary statistics 3882C about the execution of multiple program parts of any program. 3883C Nested measurements of program parts are possible. 3884C ------------------------------------------------------------ 3885C 3886C* Usage: 3887C 3888C The usage of Monitor is naturally divided into three phases: 3889C 1. the initialization and setup phase before the start of 3890C the program or subroutines package to be measured; 3891C 2. the run phase of the program to be measured; 3892C 3. the final evaluation call. 3893C 3894C The phase 1 must start with exactly one call of the subroutine 3895C MONINI, which passes a title string and a logical unit for 3896C later statistics output and possible error messages to the 3897C package. This call follows a number of calls of the subroutine 3898C MONDEF, where each call associates an identification string 3899C to a positive integer number, called the measurement index 3900C - up to maxtab, where maxtab is a package constant. Multiple 3901C measurement indices may be used for measurements of multiple 3902C program parts. The index 0 must also be associated with some 3903C identification string, and corresponds to all parts of the 3904C measured program from the measurement start call till the final 3905C evaluation call, which are not associated with specific positive 3906C measurement indices. After all necessary MONDEF calls are done, 3907C the measurements are started at begin of the program to be 3908C measured by a parameterless call of MONSRT. 3909C In phase 2, each program part to be measured must be immediately 3910C preceeded by a call of the subroutine MONON with the associated 3911C measurement index, and must be immediately followed by a call of 3912C the subroutine MONOFF with the same measurement index. Measure- 3913C ments of nested program parts are possible, and nesting is allowed 3914C up to the number mnest, where mnest is a package constant. 3915C Calling MONOFF without a preceeding MONON call with the same 3916C measurement index, or calling one of these subroutines with a 3917C measurement index not previously defined by a MONDEF call causes 3918C an error stop of the program. 3919C Finally at the end of the program to be measured, the parameter- 3920C less call of the subroutine MONEND closes all measurements and 3921C prints the summary statistics. 3922C As delivered, maxtab has a value 20 and mnest a value 10, but 3923C both constants may be increased, if needed, to any possible 3924C integer value, by simply changing it's values in the first 3925C parameter statement of the subroutine MONTOR below. 3926C 3927C* Subroutines and their parameters: 3928C ================================= 3929C 3930C MONINI(CIDENT,LUMON) : Initialize Monitor 3931C CIDENT char*20 Identification string for the total measurement 3932C ( printed in summary ) 3933C LUMON int The logical unit for printing out the summary 3934C 3935C MONDEF(MESIND,CIDMES) : Define one measurement index 3936C MESIND int >=1 : measurement index for a specific part 3937C = 0 : measurement index for all remaining parts 3938C (i.e. not belonging to parts with 3939C index >=1) 3940C CIDMES char*15 Identification string for the part associated 3941C with MESIND ( printed in summary ) 3942C 3943C MONSRT : Start measurements 3944C (no parameters) 3945C 3946C MONON(MESIND) : Start measurement of a specific part 3947C MESIND int >=1 : measurement index for a specific part 3948C 3949C MONOFF(MESIND) : Stop measurement of a specific part 3950C MESIND int >=1 : measurement index for a specific part 3951C 3952C MONEND : Finish measurements and print summary 3953C (no parameters) 3954C 3955C 3956C* Example: 3957C Calling sequence: 3958C 3959C CALL MONINI (' Example',6) 3960C CALL MONDEF (0,'Solver') 3961C CALL MONDEF (1,'User function') 3962C CALL MONDEF (2,'User matrix') 3963C CALL MONSRT () 3964C ... 3965C program to be measured (part without specific measurement index) 3966C ... 3967C 1 CONTINUE 3968C ... 3969C CALL MONON (2) 3970C ... user matrix code ... 3971C CALL MONOFF(2) 3972C ... 3973C program to be measured (part without specific measurement index) 3974C ... 3975C CALL MONON (1) 3976C ... user function code ... 3977C CALL MONOFF(1) 3978C ... 3979C program to be measured (part without specific measurement index) 3980C ... 3981C IF (no termination) GOTO 1 3982C ... 3983C CALL MONEND () 3984C ------------------------------------------------------------ 3985C 3986 SUBROUTINE MONTOR 3987 PARAMETER(MAXTAB=20,MNEST=10) 3988 CHARACTER*15 NAME(MAXTAB),NAME0 3989 CHARACTER*20 TEXT 3990 CHARACTER*(*) TEXTH 3991 CHARACTER*(*) NAMEH 3992 REAL SEC(MAXTAB),ASEC(MAXTAB),PC1(MAXTAB),PC2(MAXTAB) 3993 INTEGER COUNT(MAXTAB),INDACT(MNEST) 3994 LOGICAL QON(MAXTAB) 3995 INTEGER IOUNIT 3996C 3997 SAVE SEC,COUNT,ASEC,PC1,PC2,INDXO,TIME1,TIME0,MAXIND,NAME 3998 SAVE SEC0,NAME0,TEXT,MONI,QON,IONCNT,INDACT 3999C 4000C 4001 DATA MONI/6/ , INFO/1/ , IGRAPH/1/ 4002C 4003 RETURN 4004C 4005C initialize monitor 4006C 4007 ENTRY MONINI (TEXTH,IOUNIT) 4008C 4009 MONI=IOUNIT 4010 MAXIND=0 4011 TEXT=TEXTH 4012 DO 100 I=1,MAXTAB 4013 SEC(I)=0. 4014 ASEC(I)=0. 4015 COUNT(I)=0 4016 QON(I)=.FALSE. 4017100 CONTINUE 4018 DO 105 I=1,MNEST 4019 INDACT(I)=0 4020105 CONTINUE 4021C 4022 SEC0=0. 4023 IONCNT=0 4024 RETURN 4025C 4026C define one monitor entry 4027C 4028 ENTRY MONDEF(INDX,NAMEH) 4029 IF(INDX.LT.0 .OR. INDX.GT.MAXTAB) GOTO 1190 4030 IF (INDX.GT.MAXIND) MAXIND=INDX 4031 IF (INDX.GT.0) THEN 4032 IF (COUNT(INDX).GT.0) GOTO 1290 4033 ENDIF 4034 IF (INDX.EQ.0) THEN 4035 NAME0 = NAMEH 4036 ELSE 4037 NAME(INDX) = NAMEH 4038 ENDIF 4039 RETURN 4040C 4041C start monitor measurements 4042C 4043 ENTRY MONSRT() 4044 CALL SECOND (TIME1) 4045C 4046C if(igraph.gt.0) call gmini(maxind,name0,name) 4047C 4048 RETURN 4049C 4050C start one measurement 4051C 4052 ENTRY MONON (INDX) 4053 IF(INDX.GT.MAXIND.OR.INDX.LE.0) GOTO 1010 4054 IF (QON(INDX)) GOTO 1030 4055 CALL SECOND(ASEC(INDX)) 4056 QON(INDX)=.TRUE. 4057 IF (IONCNT.EQ.0) THEN 4058 SEC0=SEC0+ASEC(INDX)-TIME1 4059 ELSE 4060 INDXO=INDACT(IONCNT) 4061 SEC(INDXO)=SEC(INDXO)+ASEC(INDX)-ASEC(INDXO) 4062 ENDIF 4063 IONCNT=IONCNT+1 4064 INDACT(IONCNT)=INDX 4065 IF(INFO.GT.1) WRITE(MONI,*) ' enter',NAME(INDX),ASEC(INDX) 4066C 4067C if(igraph.gt.0) call gmon(indx,sec0) 4068C 4069 RETURN 4070C 4071C stop one measurement 4072C 4073 ENTRY MONOFF (INDX) 4074 IF(INDX.GT.MAXIND.OR.INDX.LE.0) GOTO 1010 4075 IF (.NOT. QON(INDX)) GOTO 1040 4076 CALL SECOND(TIME2) 4077 QON(INDX)=.FALSE. 4078 SEC(INDX)=SEC(INDX)+TIME2-ASEC(INDX) 4079 COUNT(INDX)=COUNT(INDX)+1 4080 IONCNT=IONCNT-1 4081 IF (IONCNT.EQ.0) THEN 4082 TIME1=TIME2 4083 ELSE 4084 ASEC(INDACT(IONCNT))=TIME2 4085 ENDIF 4086 IF(INFO.GT.1) WRITE(MONI,*) ' exit ',NAME(INDX),TIME2 4087C 4088C if(igraph.gt.0) call gmoff(indx,sec(indx)) 4089C 4090 RETURN 4091C 4092C terminate monitor and print statistics 4093C 4094 ENTRY MONEND 4095 CALL SECOND (TIME0) 4096 SEC0=SEC0+TIME0-TIME1 4097C 4098 SUM=1.E-10 4099 DO 200 I=1,MAXIND 4100 SUM=SUM+SEC(I) 4101 IF(COUNT(I).LE.0) GOTO 200 4102 ASEC(I)=SEC(I)/FLOAT(COUNT(I)) 4103200 CONTINUE 4104 SUM0=SUM+SEC0 4105C 4106 DO 250 I=1,MAXIND 4107 PC1(I)=100.*SEC(I)/SUM0 4108 PC2(I)=100.*SEC(I)/SUM 4109250 CONTINUE 4110 PC10=100.*SEC0/SUM0 4111 PC20=100.*SEC0/SUM 4112C 4113 WRITE(MONI,9500) 4114 WRITE(MONI,9510) 4115 WRITE(MONI,9505) 41169500 FORMAT(///) 41179510 FORMAT(1X,75('#')) 41189505 FORMAT(' #',73X,'#') 4119 WRITE(MONI,9505) 4120 WRITE(MONI,9512) TEXT 41219512 FORMAT(' # Results from time monitor program for: ',A29,2X,'#') 4122 WRITE(MONI,9505) 4123 WRITE(MONI,9514) SUM0,SUM 41249514 FORMAT(' # Total time:',F11.3,5X,'Sum of parts:',F11.3,19X,'#') 4125 WRITE(MONI,9505) 4126 WRITE(MONI,9520) 41279520 FORMAT(' # ',2X,'name',12X,'calls',7X,'time',4X,'av-time', 4128 1 4X,'% total',6X,'% sum #') 4129C 4130 I0=1 4131 WRITE(MONI,9550) NAME0,I0,SEC0,SEC0,PC10,PC20 41329550 FORMAT(' # ',A15,I8,F11.3,F11.4,F11.2,F11.2,' #') 4133C 4134 DO 300 I=1,MAXIND 4135 WRITE(MONI,9550) NAME(I),COUNT(I),SEC(I),ASEC(I),PC1(I),PC2(I) 4136300 CONTINUE 4137C 4138C 4139 WRITE(MONI,9505) 4140 WRITE(MONI,9510) 4141 WRITE(MONI,9500) 4142C 4143C 4144C IF(IGRAPH.GT.0) CALL GMEND 4145C 4146 RETURN 4147C 4148C error exits 4149C 41501010 CONTINUE 4151 WRITE(MONI,9010) INDX 41529010 FORMAT(/,' error in subroutine monon or monoff',/, 4153 $ ' indx out of range indx=',I4) 4154 GOTO 1111 4155C 41561020 CONTINUE 4157 WRITE(MONI,9020) INDX 41589020 FORMAT(/,' error in subroutine monoff',/,' indx out of range',/, 4159 1 ' indx=',I4) 4160 GOTO 1111 4161C 41621030 CONTINUE 4163 WRITE(MONI,9030) INDX 41649030 FORMAT(/,' error in subroutine monon',/, 4165 $ ' measurement is already running for ', 4166 1 ' indx=',I4) 4167 GOTO 1111 4168C 41691040 CONTINUE 4170 WRITE(MONI,9040) INDX 41719040 FORMAT(/,' error in subroutine monoff',/, 4172 $ ' measurement has never been activated for ', 4173 1 ' indx=',I4) 4174 GOTO 1111 4175C 41761190 CONTINUE 4177 WRITE(MONI,9190) MAXTAB,INDX 41789190 FORMAT(/,' error in subroutine mondef',/,' indx gt ',I4,/, 4179 1 ' indx=',I4) 4180 GOTO 1111 4181C 41821290 CONTINUE 4183 WRITE(MONI,9290) INDX 41849290 FORMAT(/,' error in subroutine mondef',/,' indx = ',I4, 4185 1 ' already in use' ) 4186 GOTO 1111 4187C 41881111 STOP 4189C 4190C end subroutine monitor 4191C 4192 END 4193C 4194C* Group Machine dependent subroutines and functions 4195C 4196 SUBROUTINE SECOND(TIME) 4197 REAL TIME 4198 TIME=0.0E0 4199 RETURN 4200 END 4201 DOUBLE PRECISION FUNCTION D1MACH(I) 4202C 4203C DOUBLE-PRECISION MACHINE CONSTANTS 4204C 4205C D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE. 4206C 4207C D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE. 4208C 4209C D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING. 4210C 4211C D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING. 4212C 4213C D1MACH( 5) = LOG10(B) 4214C 4215C TO ALTER THIS FUNCTION FOR A PARTICULAR ENVIRONMENT, 4216C THE DESIRED SET OF DATA STATEMENTS SHOULD BE ACTIVATED BY 4217C REMOVING THE C FROM COLUMN 1. 4218C ON RARE MACHINES A STATIC STATEMENT MAY NEED TO BE ADDED. 4219C (BUT PROBABLY MORE SYSTEMS PROHIBIT IT THAN REQUIRE IT.) 4220C 4221C FOR IEEE-ARITHMETIC MACHINES (BINARY STANDARD), ONE OF THE FIRST 4222C TWO SETS OF CONSTANTS BELOW SHOULD BE APPROPRIATE. 4223C 4224C WHERE POSSIBLE, DECIMAL, OCTAL OR HEXADECIMAL CONSTANTS ARE USED 4225C TO SPECIFY THE CONSTANTS EXACTLY. SOMETIMES THIS REQUIRES USING 4226C EQUIVALENT INTEGER ARRAYS. IF YOUR COMPILER USES HALF-WORD 4227C INTEGERS BY DEFAULT (SOMETIMES CALLED INTEGER*2), YOU MAY NEED TO 4228C CHANGE INTEGER TO INTEGER*4 OR OTHERWISE INSTRUCT YOUR COMPILER 4229C TO USE FULL-WORD INTEGERS IN THE NEXT 5 DECLARATIONS. 4230C 4231 INTEGER SMALL(2) 4232 INTEGER LARGE(2) 4233 INTEGER RIGHT(2) 4234 INTEGER DIVER(2) 4235 INTEGER LOG10(2) 4236C 4237 DOUBLE PRECISION DMACH(5) 4238C 4239 EQUIVALENCE (DMACH(1),SMALL(1)) 4240 EQUIVALENCE (DMACH(2),LARGE(1)) 4241 EQUIVALENCE (DMACH(3),RIGHT(1)) 4242 EQUIVALENCE (DMACH(4),DIVER(1)) 4243 EQUIVALENCE (DMACH(5),LOG10(1)) 4244C 4245C MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES, SUCH AS THE AT&T 4246C 3B SERIES AND MOTOROLA 68000 BASED MACHINES (E.G. SUN 3 AND AT&T 4247C PC 7300), IN WHICH THE MOST SIGNIFICANT BYTE IS STORED FIRST. 4248C 4249C DATA SMALL(1),SMALL(2) / 1048576, 0 / 4250C DATA LARGE(1),LARGE(2) / 2146435071, -1 / 4251C DATA RIGHT(1),RIGHT(2) / 1017118720, 0 / 4252C DATA DIVER(1),DIVER(2) / 1018167296, 0 / 4253C DATA LOG10(1),LOG10(2) / 1070810131, 1352628735 / 4254C 4255C MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES AND 8087-BASED 4256C MICROS, SUCH AS THE IBM PC AND AT&T 6300, IN WHICH THE LEAST 4257C SIGNIFICANT BYTE IS STORED FIRST. 4258C 4259 DATA SMALL(1),SMALL(2) / 0, 1048576 / 4260 DATA LARGE(1),LARGE(2) / -1, 2146435071 / 4261 DATA RIGHT(1),RIGHT(2) / 0, 1017118720 / 4262 DATA DIVER(1),DIVER(2) / 0, 1018167296 / 4263 DATA LOG10(1),LOG10(2) / 1352628735, 1070810131 / 4264C 4265C MACHINE CONSTANTS FOR AMDAHL MACHINES. 4266C 4267C DATA SMALL(1),SMALL(2) / 1048576, 0 / 4268C DATA LARGE(1),LARGE(2) / 2147483647, -1 / 4269C DATA RIGHT(1),RIGHT(2) / 856686592, 0 / 4270C DATA DIVER(1),DIVER(2) / 873463808, 0 / 4271C DATA LOG10(1),LOG10(2) / 1091781651, 1352628735 / 4272C 4273C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM. 4274C 4275C DATA SMALL(1) / ZC00800000 / 4276C DATA SMALL(2) / Z000000000 / 4277C 4278C DATA LARGE(1) / ZDFFFFFFFF / 4279C DATA LARGE(2) / ZFFFFFFFFF / 4280C 4281C DATA RIGHT(1) / ZCC5800000 / 4282C DATA RIGHT(2) / Z000000000 / 4283C 4284C DATA DIVER(1) / ZCC6800000 / 4285C DATA DIVER(2) / Z000000000 / 4286C 4287C DATA LOG10(1) / ZD00E730E7 / 4288C DATA LOG10(2) / ZC77800DC0 / 4289C 4290C MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM. 4291C 4292C DATA SMALL(1) / O1771000000000000 / 4293C DATA SMALL(2) / O0000000000000000 / 4294C 4295C DATA LARGE(1) / O0777777777777777 / 4296C DATA LARGE(2) / O0007777777777777 / 4297C 4298C DATA RIGHT(1) / O1461000000000000 / 4299C DATA RIGHT(2) / O0000000000000000 / 4300C 4301C DATA DIVER(1) / O1451000000000000 / 4302C DATA DIVER(2) / O0000000000000000 / 4303C 4304C DATA LOG10(1) / O1157163034761674 / 4305C DATA LOG10(2) / O0006677466732724 / 4306C 4307C MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS. 4308C 4309C DATA SMALL(1) / O1771000000000000 / 4310C DATA SMALL(2) / O7770000000000000 / 4311C 4312C DATA LARGE(1) / O0777777777777777 / 4313C DATA LARGE(2) / O7777777777777777 / 4314C 4315C DATA RIGHT(1) / O1461000000000000 / 4316C DATA RIGHT(2) / O0000000000000000 / 4317C 4318C DATA DIVER(1) / O1451000000000000 / 4319C DATA DIVER(2) / O0000000000000000 / 4320C 4321C DATA LOG10(1) / O1157163034761674 / 4322C DATA LOG10(2) / O0006677466732724 / 4323C 4324C MACHINE CONSTANTS FOR FTN4 ON THE CDC 6000/7000 SERIES. 4325C 4326C DATA SMALL(1) / 00564000000000000000B / 4327C DATA SMALL(2) / 00000000000000000000B / 4328C 4329C DATA LARGE(1) / 37757777777777777777B / 4330C DATA LARGE(2) / 37157777777777777774B / 4331C 4332C DATA RIGHT(1) / 15624000000000000000B / 4333C DATA RIGHT(2) / 00000000000000000000B / 4334C 4335C DATA DIVER(1) / 15634000000000000000B / 4336C DATA DIVER(2) / 00000000000000000000B / 4337C 4338C DATA LOG10(1) / 17164642023241175717B / 4339C DATA LOG10(2) / 16367571421742254654B / 4340C 4341C MACHINE CONSTANTS FOR FTN5 ON THE CDC 6000/7000 SERIES. 4342C 4343C DATA SMALL(1) / O"00564000000000000000" / 4344C DATA SMALL(2) / O"00000000000000000000" / 4345C 4346C DATA LARGE(1) / O"37757777777777777777" / 4347C DATA LARGE(2) / O"37157777777777777774" / 4348C 4349C DATA RIGHT(1) / O"15624000000000000000" / 4350C DATA RIGHT(2) / O"00000000000000000000" / 4351C 4352C DATA DIVER(1) / O"15634000000000000000" / 4353C DATA DIVER(2) / O"00000000000000000000" / 4354C 4355C DATA LOG10(1) / O"17164642023241175717" / 4356C DATA LOG10(2) / O"16367571421742254654" / 4357C 4358C MACHINE CONSTANTS FOR CONVEX C-1 4359C 4360C DATA SMALL(1),SMALL(2) / '00100000'X, '00000000'X / 4361C DATA LARGE(1),LARGE(2) / '7FFFFFFF'X, 'FFFFFFFF'X / 4362C DATA RIGHT(1),RIGHT(2) / '3CC00000'X, '00000000'X / 4363C DATA DIVER(1),DIVER(2) / '3CD00000'X, '00000000'X / 4364C DATA LOG10(1),LOG10(2) / '3FF34413'X, '509F79FF'X / 4365C 4366C MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3. 4367C 4368C DATA SMALL(1) / 201354000000000000000B / 4369C DATA SMALL(2) / 000000000000000000000B / 4370C 4371C DATA LARGE(1) / 577767777777777777777B / 4372C DATA LARGE(2) / 000007777777777777776B / 4373C 4374C DATA RIGHT(1) / 376434000000000000000B / 4375C DATA RIGHT(2) / 000000000000000000000B / 4376C 4377C DATA DIVER(1) / 376444000000000000000B / 4378C DATA DIVER(2) / 000000000000000000000B / 4379C 4380C DATA LOG10(1) / 377774642023241175717B / 4381C DATA LOG10(2) / 000007571421742254654B / 4382C 4383C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200 4384C 4385C NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING LINE - 4386C STATIC DMACH(5) 4387C 4388C DATA SMALL/20K,3*0/,LARGE/77777K,3*177777K/ 4389C DATA RIGHT/31420K,3*0/,DIVER/32020K,3*0/ 4390C DATA LOG10/40423K,42023K,50237K,74776K/ 4391C 4392C MACHINE CONSTANTS FOR THE HARRIS SLASH 6 AND SLASH 7 4393C 4394C DATA SMALL(1),SMALL(2) / '20000000, '00000201 / 4395C DATA LARGE(1),LARGE(2) / '37777777, '37777577 / 4396C DATA RIGHT(1),RIGHT(2) / '20000000, '00000333 / 4397C DATA DIVER(1),DIVER(2) / '20000000, '00000334 / 4398C DATA LOG10(1),LOG10(2) / '23210115, '10237777 / 4399C 4400C MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES. 4401C 4402C DATA SMALL(1),SMALL(2) / O402400000000, O000000000000 / 4403C DATA LARGE(1),LARGE(2) / O376777777777, O777777777777 / 4404C DATA RIGHT(1),RIGHT(2) / O604400000000, O000000000000 / 4405C DATA DIVER(1),DIVER(2) / O606400000000, O000000000000 / 4406C DATA LOG10(1),LOG10(2) / O776464202324, O117571775714 / 4407C 4408C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, 4409C THE XEROX SIGMA 5/7/9 AND THE SEL SYSTEMS 85/86. 4410C 4411C DATA SMALL(1),SMALL(2) / Z00100000, Z00000000 / 4412C DATA LARGE(1),LARGE(2) / Z7FFFFFFF, ZFFFFFFFF / 4413C DATA RIGHT(1),RIGHT(2) / Z33100000, Z00000000 / 4414C DATA DIVER(1),DIVER(2) / Z34100000, Z00000000 / 4415C DATA LOG10(1),LOG10(2) / Z41134413, Z509F79FF / 4416C 4417C MACHINE CONSTANTS FOR THE INTERDATA 8/32 4418C WITH THE UNIX SYSTEM FORTRAN 77 COMPILER. 4419C 4420C FOR THE INTERDATA FORTRAN VII COMPILER REPLACE 4421C THE Z'S SPECIFYING HEX CONSTANTS WITH Y'S. 4422C 4423C DATA SMALL(1),SMALL(2) / Z'00100000', Z'00000000' / 4424C DATA LARGE(1),LARGE(2) / Z'7EFFFFFF', Z'FFFFFFFF' / 4425C DATA RIGHT(1),RIGHT(2) / Z'33100000', Z'00000000' / 4426C DATA DIVER(1),DIVER(2) / Z'34100000', Z'00000000' / 4427C DATA LOG10(1),LOG10(2) / Z'41134413', Z'509F79FF' / 4428C 4429C MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR). 4430C 4431C DATA SMALL(1),SMALL(2) / "033400000000, "000000000000 / 4432C DATA LARGE(1),LARGE(2) / "377777777777, "344777777777 / 4433C DATA RIGHT(1),RIGHT(2) / "113400000000, "000000000000 / 4434C DATA DIVER(1),DIVER(2) / "114400000000, "000000000000 / 4435C DATA LOG10(1),LOG10(2) / "177464202324, "144117571776 / 4436C 4437C MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR). 4438C 4439C DATA SMALL(1),SMALL(2) / "000400000000, "000000000000 / 4440C DATA LARGE(1),LARGE(2) / "377777777777, "377777777777 / 4441C DATA RIGHT(1),RIGHT(2) / "103400000000, "000000000000 / 4442C DATA DIVER(1),DIVER(2) / "104400000000, "000000000000 / 4443C DATA LOG10(1),LOG10(2) / "177464202324, "047674776746 / 4444C 4445C MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING 4446C 32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). 4447C 4448C DATA SMALL(1),SMALL(2) / 8388608, 0 / 4449C DATA LARGE(1),LARGE(2) / 2147483647, -1 / 4450C DATA RIGHT(1),RIGHT(2) / 612368384, 0 / 4451C DATA DIVER(1),DIVER(2) / 620756992, 0 / 4452C DATA LOG10(1),LOG10(2) / 1067065498, -2063872008 / 4453C 4454C DATA SMALL(1),SMALL(2) / O00040000000, O00000000000 / 4455C DATA LARGE(1),LARGE(2) / O17777777777, O37777777777 / 4456C DATA RIGHT(1),RIGHT(2) / O04440000000, O00000000000 / 4457C DATA DIVER(1),DIVER(2) / O04500000000, O00000000000 / 4458C DATA LOG10(1),LOG10(2) / O07746420232, O20476747770 / 4459C 4460C MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING 4461C 16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). 4462C 4463C DATA SMALL(1),SMALL(2) / 128, 0 / 4464C DATA SMALL(3),SMALL(4) / 0, 0 / 4465C 4466C DATA LARGE(1),LARGE(2) / 32767, -1 / 4467C DATA LARGE(3),LARGE(4) / -1, -1 / 4468C 4469C DATA RIGHT(1),RIGHT(2) / 9344, 0 / 4470C DATA RIGHT(3),RIGHT(4) / 0, 0 / 4471C 4472C DATA DIVER(1),DIVER(2) / 9472, 0 / 4473C DATA DIVER(3),DIVER(4) / 0, 0 / 4474C 4475C DATA LOG10(1),LOG10(2) / 16282, 8346 / 4476C DATA LOG10(3),LOG10(4) / -31493, -12296 / 4477C 4478C DATA SMALL(1),SMALL(2) / O000200, O000000 / 4479C DATA SMALL(3),SMALL(4) / O000000, O000000 / 4480C 4481C DATA LARGE(1),LARGE(2) / O077777, O177777 / 4482C DATA LARGE(3),LARGE(4) / O177777, O177777 / 4483C 4484C DATA RIGHT(1),RIGHT(2) / O022200, O000000 / 4485C DATA RIGHT(3),RIGHT(4) / O000000, O000000 / 4486C 4487C DATA DIVER(1),DIVER(2) / O022400, O000000 / 4488C DATA DIVER(3),DIVER(4) / O000000, O000000 / 4489C 4490C DATA LOG10(1),LOG10(2) / O037632, O020232 / 4491C DATA LOG10(3),LOG10(4) / O102373, O147770 / 4492C 4493C MACHINE CONSTANTS FOR THE PRIME 50 SERIES SYSTEMS 4494C WTIH 32-BIT INTEGERS AND 64V MODE INSTRUCTIONS, 4495C SUPPLIED BY IGOR BRAY. 4496C 4497C DATA SMALL(1),SMALL(2) / :10000000000, :00000100001 / 4498C DATA LARGE(1),LARGE(2) / :17777777777, :37777677775 / 4499C DATA RIGHT(1),RIGHT(2) / :10000000000, :00000000122 / 4500C DATA DIVER(1),DIVER(2) / :10000000000, :00000000123 / 4501C DATA LOG10(1),LOG10(2) / :11504046501, :07674600177 / 4502C 4503C MACHINE CONSTANTS FOR THE SEQUENT BALANCE 8000 4504C 4505C DATA SMALL(1),SMALL(2) / $00000000, $00100000 / 4506C DATA LARGE(1),LARGE(2) / $FFFFFFFF, $7FEFFFFF / 4507C DATA RIGHT(1),RIGHT(2) / $00000000, $3CA00000 / 4508C DATA DIVER(1),DIVER(2) / $00000000, $3CB00000 / 4509C DATA LOG10(1),LOG10(2) / $509F79FF, $3FD34413 / 4510C 4511C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. 4512C 4513C DATA SMALL(1),SMALL(2) / O000040000000, O000000000000 / 4514C DATA LARGE(1),LARGE(2) / O377777777777, O777777777777 / 4515C DATA RIGHT(1),RIGHT(2) / O170540000000, O000000000000 / 4516C DATA DIVER(1),DIVER(2) / O170640000000, O000000000000 / 4517C DATA LOG10(1),LOG10(2) / O177746420232, O411757177572 / 4518C 4519C MACHINE CONSTANTS FOR THE VAX UNIX F77 COMPILER 4520C 4521C DATA SMALL(1),SMALL(2) / 128, 0 / 4522C DATA LARGE(1),LARGE(2) / -32769, -1 / 4523C DATA RIGHT(1),RIGHT(2) / 9344, 0 / 4524C DATA DIVER(1),DIVER(2) / 9472, 0 / 4525C DATA LOG10(1),LOG10(2) / 546979738, -805796613 / 4526C 4527C MACHINE CONSTANTS FOR THE VAX-11 WITH 4528C FORTRAN IV-PLUS COMPILER 4529C 4530C DATA SMALL(1),SMALL(2) / Z00000080, Z00000000 / 4531C DATA LARGE(1),LARGE(2) / ZFFFF7FFF, ZFFFFFFFF / 4532C DATA RIGHT(1),RIGHT(2) / Z00002480, Z00000000 / 4533C DATA DIVER(1),DIVER(2) / Z00002500, Z00000000 / 4534C DATA LOG10(1),LOG10(2) / Z209A3F9A, ZCFF884FB / 4535C 4536C MACHINE CONSTANTS FOR VAX/VMS VERSION 2.2 4537C 4538C DATA SMALL(1),SMALL(2) / '80'X, '0'X / 4539C DATA LARGE(1),LARGE(2) / 'FFFF7FFF'X, 'FFFFFFFF'X / 4540C DATA RIGHT(1),RIGHT(2) / '2480'X, '0'X / 4541C DATA DIVER(1),DIVER(2) / '2500'X, '0'X / 4542C DATA LOG10(1),LOG10(2) / '209A3F9A'X, 'CFF884FB'X / 4543C 4544C/6S 4545 IF (I .LT. 1 .OR. I .GT. 6) GOTO 999 4546 IF (I .LE. 5 ) THEN 4547 D1MACH = DMACH(I) 4548 ELSE IF (I .EQ. 6) THEN 4549C D1MACH = DSQRT(DMACH(1)/DMACH(3)) 4550 D1MACH = 4.94D-32 4551 ENDIF 4552 RETURN 4553 999 WRITE(6,1999) I 4554 1999 FORMAT(' D1MACH - I OUT OF BOUNDS',I10) 4555 STOP 4556 END 4557