1C PCON61.FOR 27 February 1991 2C 3 SUBROUTINE PITCON(DF,FPAR,FX,IERROR,IPAR,IWORK,LIW,NVAR,RWORK, 4 *LRW,XR,SLNAME) 5C 6C*********************************************************************** 7C 8C This is the interface program between the user and the PITCON package. 9C 10C 11C A. Introduction 12C 13C 14C PCON61 is version 6.1 of PITCON, the University of Pittsburgh continuation 15C program. 16C 17C PITCON was written by 18C 19C Professor Werner C Rheinboldt and John Burkardt, 20C Department of Mathematics and Statistics 21C University of Pittsburgh, 22C Pittsburgh, Pennsylvania, 15260, USA. 23C 24C E-Mail: wcrhein@vms.cis.pitt.edu 25C burkardt@psc.edu 26C 27C The original work on this package was partially supported by the National 28C Science Foundation under grants MCS-78-05299 and MCS-83-09926. 29C 30C PITCON computes a sequence of solution points along a one dimensional 31C manifold of a system of nonlinear equations F(X)=0 involving NVAR-1 32C equations and an NVAR dimensional unknown vector X. 33C 34C The operation of PITCON is somewhat analogous to that of an initial value 35C ODE solver. In particular, the user must begin the computation by 36C specifying an approximate initial solution, and subsequent points returned 37C by PITCON lie on the curve which passes through this initial point and is 38C implicitly defined by F(X)=0. The extra degree of freedom in the system is 39C analogous to the role of the independent variable in a differential 40C equations. 41C 42C However, PITCON does not try to solve the algebraic problem by turning it 43C into a differential equation system. Unlike differential equations, the 44C solution curve may bend and switch back in any direction, and there may be 45C many solutions for a fixed value of one of the variables. Accordingly, 46C PITCON is not required to parametrize the implicitly defined curve with a 47C fixed parameter. Instead, at each step, PITCON selects a suitable variable 48C as the current parameter and then determines the other variables as 49C functions of it. This allows PITCON to go around relatively sharp bends. 50C Moreover, if the equations were actually differentiated - that is, replaced 51C by some linearization - this would introduce an inevitable "drift" away from 52C the true solution curve. Errors at previous steps would be compounded in a 53C way that would make later solution points much less reliable than earlier 54C ones. Instead, PITCON solves the algebraic equations explicitly and each 55C solution has to pass an acceptance test in an iterative solution process 56C with tolerances provided by the user. 57C 58C PITCON is only designed for systems with one degree of freedom. However, 59C it may be used on systems with more degrees of freedom if the user reduces 60C the available degrees of freedom by the introduction of suitable constraints 61C that are added to the set of nonlinear equations. In this sense, PITCON may 62C be used to investigate the equilibrium behavior of physical systems with 63C several degrees of freedom. 64C 65C Program options include the ability to search for solutions for which a 66C given component has a specified value. Another option is a search for a 67C limit or turning point with respect to a given component; that is, of a 68C point where this particular solution component has a local extremum. 69C 70C Another feature of the program is the use of two work arrays, IWORK and 71C RWORK. All information required for continuing any interrupted computation 72C is saved in these two arrays. 73C 74C 75C B. PITCON Calling Sequence 76C 77C 78C SUBROUTINE PITCON(DF,FPAR,FX,IERROR,IPAR,IWORK,LIW,NVAR,RWORK,LRW,XR,SLVNAM) 79C 80C On the first call, PITCON expects a point XR and a routine FX defining a 81C nonlinear function F. As noted earlier, XR and FX specify a curve which 82C PITCON is to trace. On the first call, PITCON simply verifies that F(XR)=0, 83C and if this is not the case, the program attempts to correct XR to a new 84C value satisfying the equation. 85C 86C On subsequent calls, PITCON assumes that the input vector XR contains the 87C point which had been computed on the previous call. It also assumes that 88C the work arrays IWORK and RWORK contain the results of the prior 89C calculations. PITCON estimates an appropriate stepsize, computes the 90C tangent direction to the curve at the given input point, and calculates a 91C predicted new point on the curve. A form of Newton's method is used to 92C correct this point so that it lies on the curve. If the iteration is 93C successful, the code returns with a new point XR. Otherwise, the stepsize 94C may be reduced, and the calculation retried. 95C 96C Aside from its ability to produce successive points on the solution curve, 97C PITCON may be asked to search for "target points" or "limit points". 98C Target points are solution vectors for which a certain component has a 99C specified value. One might ask for all solutions for which XR(17)=4.0, for 100C instance. Limit points occur when the curve turns back in a given 101C direction, and have the property that the corresponding component of the 102C tangent vector vanishes there. 103C 104C If the user has asked for the computation of target or limit points, then 105C PITCON will usually proceed as described earlier, producing a new 106C continuation point on each return. But if a target or limit point is 107C found, such a point is returned as the value of XR, temporarily interrupting 108C the usual form of the computation. 109C 110C 111C C. Overview of PITCON parameters: 112C 113C 114C DF Input, EXTERNAL DF, routine for evaluating the Jacobian of F. 115C FPAR Input/output, REAL FPAR(*), user defined parameter array. 116C FX Input, EXTERNAL FX, routine for evaluating the function F. 117C IERROR Output, INTEGER IERROR, error return flag. 118C IPAR Input/output, INTEGER IPAR(*), user defined parameter array. 119C IWORK Input/output, INTEGER IWORK(LIW). Communication and work array. 120C LIW Input, INTEGER LIW, the dimension of IWORK. 121C NVAR Input, INTEGER NVAR, number of variables, the dimension of X. 122C RWORK Input/output, REAL RWORK(LRW), Communication and work space array. 123C LRW Input, INTEGER LRW, the dimension of RWORK. 124C XR Input/output, REAL XR(NVAR), the current solution point. 125C SLVNAM Input, EXTERNAL SLVNAM, the solver to use on linear systems. 126C 127C 128C D. Details of PITCON parameters: 129C 130C 131C DF Input, EXTERNAL DF, the name of the Jacobian evaluation routine. 132C This name must be declared EXTERNAL in the calling program. 133C 134C DF is not needed if the finite difference option is used 135C (IWORK(9)=1 or 2). In that case, only a dummy name is needed for DF. 136C 137C Otherwise, the user must write a routine which evaluates the 138C Jacobian matrix of the function FX at a given point X and stores it 139C in the FJAC array in accordance with the format used by the solver 140C specified in SLVNAM. 141C 142C In the simplest case, when the full matrix solverDENSLV solver 143C provided with the package is used, DF must store D F(I)/D X(J) into 144C FJAC(I,J). 145C 146C The array to contain the Jacobian will be zeroed out before DF is 147C called, so that only nonzero elements need to be stored. DF must 148C have the form: 149C 150C SUBROUTINE DF(NVAR,FPAR,IPAR,X,FJAC,IERROR) 151C 152C NVAR Input, INTEGER NVAR, number of variables. 153C 154C FPAR Input, REAL FPAR(*), vector for passing real parameters. 155C This vector is not used by the program, and is only provided 156C for the user's convenience. 157C 158C IPAR Input, INTEGER IPAR(*), vector for passing integer 159C parameters. This vector is not used by the program, and is 160C only provided for the user's convenience. 161C 162C X Input, REAL X(NVAR), the point at which the Jacobian is 163C desired. 164C 165C FJAC Output, REAL FJAC(*), the array containing the Jacobian. 166C 167C If DENSLV is the solver: FJAC must be dimensioned 168C FJAC(NVAR,NVAR) as shown above, and DF sets 169C FJAC(I,J)=D F(I)/DX(J). 170C 171C If BANSLV is the solver: the main portion of the Jacobian, 172C rows and columns 1 through NVAR-1, is assumed to be a banded 173C matrix in the standard LINPACK form with lower bandwidth ML 174C and upper bandwidth MU. However, the final column of the 175C Jacobian is allowed to be full. 176C 177C BANSLV will pass to DF the beginning of the storage for 178C FJAC, but it is probably best not to doubly dimension FJAC 179C inside of DF, since it is a "hybrid" object. The first 180C portion of it is a (2*ML+MU+1, NEQN) array, followed by a 181C single column of length NEQN (the last column of the 182C Jacobian). Thus the simplest approach is to declare FJAC to 183C be a vector, and then then to store values as follows: 184C 185C If J is less than NVAR, then 186C if I-J .LE. ML and J-I .LE. MU, 187C set K=(2*ML+MU)*J + I - ML 188C set FJAC(K)=D F(I)/DX(J). 189C else 190C do nothing, index is outside the band 191C endif 192C Else if J equals NVAR, then 193C set K=(2*ML+MU+1)*(NVAR-1)+I, 194C set FJAC(K)=D F(I)/DX(J). 195C endif. 196C 197C IERROR Output, INTEGER IERROR, error return flag. DF should set 198C this to 0 for normal return, nonzero if trouble. 199C 200C FPAR Input/output, REAL FPAR(*), a user defined parameter array. 201C 202C FPAR is not used in any way by PITCON. It is provided for the user's 203C convenience. It is passed to DF and FX, and hence may be used to 204C transmit information between the user calling program and these user 205C subprograms. The dimension of FPAR and its contents are up to the 206C user. Internally, the program declares DIMENSION FPAR(*) but never 207C references its contents. 208C 209C FX Input, EXTERNAL FX, the name of the routine which evaluates the 210C function. 211C 212C FX computes the value of the nonlinear function. This name must be 213C declared EXTERNAL in the calling program. FX should evaluate the 214C NVAR-1 function components at the input point X, and store the result 215C in the vector FVEC. An augmenting equation will be stored in entry 216C NVAR of FVEC by the PITCON program. 217C 218C FX should have the following form: 219C 220C SUBROUTINE FX(NVAR,FPAR,IPAR,X,FVEC,IERROR) 221C 222C NVAR Input, INTEGER NVAR, number of variables. 223C 224C FPAR Input/output, REAL FPAR(*), array of user parameters. 225C 226C IPAR Input/output, INTEGER IPAR(*), array of user parameters. 227C 228C X Input, REAL X(NVAR), the point at which function evaluation 229C is required. 230C 231C FVEC Output, REAL FVEC(NVAR), the value of the function at point 232C X. Only the first NVAR-1 entries of FVEC are to be set by 233C the routine. PITCON sets the final value itself. 234C 235C IERROR Output, INTEGER IERROR, error flag. FX should set this to 0 236C if there are no problems, or to 2 if there is a problem. 237C 238C IERROR Output, INTEGER IERROR, error return flag. 239C 240C On return from PITCON, a nonzero value of IERROR is a warning of some 241C problem which may be minor, serious, or fatal. 242C 243C 0, No errors occurred. 244C 245C 1, Insufficient storage provided in RWORK and IWORK, or NVAR is less 246C than 2. This is a fatal error, which occurs on the first call to 247C PITCON. 248C 249C 2, A user defined error condition occurred in the FX or DF 250C subroutines. PITCON treats this as a fatal error. 251C 252C 3, A numerically singular matrix was encountered. Continuation 253C cannot proceed without some redefinition of the problem. This is 254C a fatal error. 255C 256C 4, Unsuccessful corrector iteration. Loosening the tolerances 257C RWORK(1) and RWORK(2), or decreasing the maximum stepsize RWORK(4) 258C might help. This is a fatal error. 259C 260C 5, Too many corrector steps. The corrector iteration was proceeding 261C properly, but too slowly. Increase number of Newton steps 262C IWORK(17), increase the error tolerances RWORK(1) or RWORK(2), or 263C decrease RWORK(4). This is a fatal error. 264C 265C 6, Null tangent vector. A serious error which indicates a data 266C problem or singularity in the nonlinear system. This is a fatal 267C error. 268C 269C 7, Root finder failed while searching for a limit point. 270C This is a warning. It means that the target or limit point 271C computation has failed, but the main computation (computing the 272C continuation curve itself) may continue. 273C 274C 8, Limit point iteration took too many steps. This is a warning 275C error. It means that the limit point computation has failed, but 276C the main computation (computing the continuation curve itself) may 277C continue. 278C 279C 9, Undiagnosed error condition. This is a fatal error. 280C 281C IPAR Input/output, INTEGER IPAR(*), user defined parameter array. 282C 283C IPAR is not used in any way by PITCON. It is provided for the user's 284C convenience in transmitting parameters between the calling program 285C and the user routines FX and DF. IPAR is declared in the PITCON 286C program and passed through it to FX and DF, but otherwise ignored. 287C Note, however, that if BANSLV is used for the solver routine, then 288C IPAR(1) must contain the lower bandwidth, and IPAR(2) the upper 289C bandwidth of the Jacobian matrix. 290C 291C IWORK Input/output, INTEGER IWORK(LIW). Communication and workspace array. 292C 293C The specific allocation of IWORK is described in the section devoted 294C to the work arrays. Some elements of IWORK MUST be set by the user, 295C others may be set to change the way PITCON works. 296C 297C LIW Input, INTEGER LIW, the dimension of IWORK. 298C 299C The minimum acceptable value of LIW depends on the solver chosen, 300C but for either DENSLV or BANSLV, setting LIW=29+NVAR is sufficient. 301C 302C NVAR Input, INTEGER NVAR, the number of variables, the dimension of X. 303C 304C This is, of course, one greater than the number of equations or 305C functions. NVAR must be at least 2. 306C 307C RWORK Input/output, REAL RWORK(LRW), work array. 308C 309C The specific allocation of RWORK is described in the section 310C devoted to the work arrays. 311C 312C LRW Input, INTEGER LRW, the dimension of RWORK. 313C 314C The minimum acceptable value depends heavily on the solver options. 315C There is storage required for scalars, vectors, and the Jacobian 316C array. The minimum acceptable value of LRW is the sum of three 317C corresponding numbers. 318C 319C For DENSLV with user-supplied Jacobian, 320C 321C LRW=29 + 4*NVAR + NVAR*NVAR. 322C 323C For DENSLV with internally approximated Jacobian, 324C 325C LRW=29 + 6*NVAR + NVAR*NVAR. 326C 327C For BANSLV, with a Jacobian matrix with upper bandwidth MU and lower 328C bandwidth ML, and NBAND=2*ML+MU+1, with user supplied Jacobian, 329C 330C LRW=29 + 6*NVAR + (NVAR-1)*NBAND. 331C 332C For BANSLV with internally approximated Jacobian, 333C 334C LRW=29 + 9*NVAR + (NVAR-1)*NBAND. 335C 336C XR Input/output, REAL XR(NVAR), the current solution point. 337C 338C On the first call, the user should set XR to a starting point which 339C at least approximately satisfies F(XR)=0. The user need never 340C update XR again. 341C 342C Thereafter, on each return from the program with IERROR=0, XR will 343C hold the most recently computed point, whether a continuation, target 344C or limit point. 345C 346C SLVNAM Input, EXTERNAL SLVNAM, the name of the solver to use on linear 347C systems. 348C 349C The linear systems have the form A*x=b, where A is the augmented 350C Jacobian matrix. A will be square, and presumably nonsingular. 351C The routine must return the value of the solution x. 352C 353C Two possible choices for SLVNAM are "DENSLV" and "BANSLV", which are 354C the names of routines provided with the package. DENSLV is 355C appropriate for a full storage jacobian, and BANSLV for a jacobian 356C which is banded except for the last column. 357C 358C The advanced user may study the source code for these two routines 359C and write an equivalent solver more suited to a given problem. 360C 361C 362C E. The Integer Work Array IWORK 363C 364C 365C Input to the program includes the setting of some of the entries in IWORK. 366C Some of this input is optional. The user input section of IWORK involves 367C entries 1 through 8, and, possibly also 17 and 29. 368C 369C IWORK(1) must be set by the user. All other entries have default values. 370C 371C 372C IWORK(1) On first call only, the user must set IWORK(1)=0. 373C Thereafter, the program sets IWORK(1) before return to 374C explain what kind of point is being returned. This return 375C code is: 376C 377C 1 return of corrected starting point. 378C 2 return of continuation point. 379C 3 return of target point. 380C 4 return of limit point. 381C 382C NOTE: At any time, PITCON may be called with a negative 383C value of IWORK(1). This requests a check of the 384C jacobian routine against a finite difference approximation. 385C The program will call the jacobian routine, and 386C then estimate the jacobian. If IWORK(1)=-1, then it will 387C print out the value of the maximum difference, and the row 388C and column of the jacobian in which it appears. Otherwise, 389C the program will print out the entire matrix 390C FP(I,J)-DEL(J)F(I), where DEL(J)F(I) represents the finite 391C difference approximation. 392C 393C Before a call with negative IWORK(1), the current value of 394C IWORK(1) should be saved, and then restored to the previous 395C value after the call, in order to resume calculation. 396C 397C IWORK(1) does not have a default value. The user MUST set 398C it. 399C 400C IWORK(2) The component of the current continuation point XR which is 401C to be used as the continuation parameter. On first call, 402C the program is willing to use the index NVAR as a default, 403C but the user should set this value if better information is 404C available. 405C 406C After the first call, the program sets this value for each 407C step automatically unless the user prevents this by setting 408C the parameterization option IWORK(3) to a non-zero valus. 409C Note that a poor choice of IWORK(2) may cause the algorithm 410C to fail. IWORK(2) defaults to NVAR on the first step. 411C 412C IWORK(3) Parameterization option. The program would prefer to be 413C free to choose a new local parameter from step to step. 414C The value of IWORK(3) allows or prohibits this action. 415C IWORK(3)=0 allows the program to vary the parameter, 416C IWORK(3)=1 forces the program to use whatever the contents 417C of IWORK(2) are, which will not be changed from the user's 418C input or the default. The default is IWORK(3)=0. 419C 420C IWORK(4) Newton iteration Jacobian update option. 421C 0, the Jacobian is reevaluated at every step of the 422C Newton iteration. This is costly, but may result in 423C fewer Newton steps and fewer Newton iteration rejections. 424C 1, the Jacobian is evaluated only on the first and 425C IWORK(17)-th steps of the Newton process. 426C 2, the Jacobian is evaluated only when absolutely 427C necessary, namely, at the first step, and when the 428C process fails. This option is most suitable for problems 429C with mild nonlinearities. 430C 431C The default is IWORK(4)=0. 432C 433C IWORK(5) Target point index. If IWORK(5) is not zero, it is presumed 434C to be the component index between 1 and NVAR for which 435C target points are sought. In this case, the value of 436C RWORK(7) is assumed to be the target value. The program 437C will monitor every new continuation point, and if it finds 438C that a target point may lie between the new point and the 439C previous point, it will compute this target point and 440C return. This target point is defined by the property that 441C its component with the index prescribed in IWORK(5) will 442C have the value given in RWORK(7). For a given problem there 443C may be zero, one, or many target points requested. 444C The default of IWORK(5) is 0. 445C 446C IWORK(6) Limit point index. If IWORK(6) is nonzero, then the program 447C will search for limit points with respect to the component 448C with index IWORK(6); that is, of points for which the 449C IWORK(6)-th variable has a local extremum, or equivalently 450C where the IWORK(6)-th component of the tangent vector is 451C zero. The default of IWORK(6) is zero. 452C 453C IWORK(7) Control of the amount of intermediate output produced by the 454C program. IWORK(7) may have a value between 0 and 3. 455C For IWORK(7) = 0 there is no intermediate output while for 456C IWORK(7) = 3 the most intermediate output is produced. 457C The default is 1. 458C 459C IWORK(8) FORTRAN unit to which output is to be written. The 460C default value is site dependent but should be the standard 461C output device. 462C The default is 6 on the Cray, Vax and PC, or 9 on the 463C Macintosh. 464C 465C IWORK(9) Control of the Jacobian option specifying whether the user 466C has supplied a Jacobian routine, or wants the program 467C to approximate the Jacobian. 468C 0, the user has supplied the Jacobian. 469C 1, program is to use forward difference approximation. 470C 2, program is to use central difference approximation. 471C IWORK(9) defaults to 0. 472C 473C IWORK(10) State indicator of the progress of the program. 474C The values are: 475C 0, start up with unchecked starting point. 476C 1, first step. Corrected starting point available. 477C 2, two successive continuation points available, as well 478C as the tangent vector at the oldest of them. 479C 3, two successive continuation points available, as well 480C as the tangent vector at the newest of them. 481C 482C IWORK(11) Index of the last computed target point. This is used to 483C avoid repeated computation of a target point. If a target 484C point has been found, then the target index IWORK(5) is 485C copied into IWORK(11). 486C 487C IWORK(12) Second best choice for the local parameterization index. 488C This index may be tried if the first choice causes poor 489C performance in the Newton corrector. 490C 491C IWORK(13) Beginning location in IWORK of unused integer work space 492C available for use by the solver. 493C 494C IWORK(14) LIW, the user declared dimension of the array IWORK. 495C 496C IWORK(15) Beginning location in RWORK of unused real work space 497C available for use by the solver. 498C 499C IWORK(16) LRW, the user declared dimension of RWORK. 500C 501C IWORK(17) Maximum number of corrector steps allowed during one run 502C of the Newton process in which the Jacobian is updated at 503C every step. If the Jacobian is only evaluated at 504C the beginning of the Newton iteration then 2*IWORK(17) steps 505C are allowed. 506C IWORK(17) must be greater than 0. It defaults to 10. 507C 508C IWORK(18) Number of stepsize reductions that were needed for 509C producing the last continuation point. 510C 511C IWORK(19) Total number of calls to the user Jacobian routine DF. 512C 513C IWORK(20) Total number of calls to the matrix factorization routine. 514C If DENSLV is the chose solver then factorization is done by 515C the LINPACK routine SGEFA. If BANSLV is the solver, the 516C LINPACK routine SGBFA will be used. 517C 518C IWORK(21) Total number of calls to the back-substitution routine. 519C If DENSLV is the chosen solver, then back substitution is 520C done by the LINPACK routine SGESL. If BANSLV is used, then 521C the LINPACK routine SGBSL will be used. 522C 523C IWORK(22) Total number of calls to the user function routine FX. 524C 525C IWORK(23) Total number of steps taken in limit point iterations. 526C Each step involves determining an approximate limit point 527C and applying a Newton iteration to correct it. 528C 529C IWORK(24) Total number of Newton corrector steps used during the 530C computation of target points. 531C 532C IWORK(25) Total number of Newton steps taken during the correction 533C of a starting point or the continuation points. 534C 535C IWORK(26) Total number of predictor stepsize-reductions needed 536C since the start of the continuation procesds. 537C 538C IWORK(27) Total number of calls to the program. This also 539C corresponds to the number of points computed. 540C 541C IWORK(28) Total number of Newton steps taken during current iteration. 542C 543C IWORK(30) and on are reserved for use by the linear equation solver, 544C and typically are used for pivoting. 545C 546C 547C F. The Real Work Array RWORK 548C 549C 550C Input to the program includes the setting of some of the entries in RWORK. 551C Some of this input is optional. The user input section of RWORK involves 552C entries 1 through 7 and possibly 20. All entries of RWORK have default 553C values. 554C 555C 556C RWORK(1) Absolute error tolerance. This value is used mainly during 557C the Newton iteration. RWORK(1) defaults to SQRT(EPMACH) 558C where EPMACH is the machine relative precision stored in 559C RWORK(8). 560C 561C RWORK(2) Relative error tolerance. This value is used mainly during 562C the Newton iteration. RWORK(2) defaults to SQRT(EPMACH) 563C where EPMACH is the machine relative precision stored in 564C RWORK(8). 565C 566C RWORK(3) Minimum allowable predictor stepsize. If failures of 567C the Newton correction force the stepsize down to this level, 568C then the program will give up. The default value is 569C SQRT(EPMACH). 570C 571C RWORK(4) Maximum allowable predictor step. Too generous a value 572C may cause erratic behavior of the program. The default 573C value is SQRT(NVAR). 574C 575C RWORK(5) Predictor stepsize. On first call, it should be set by 576C the user. Thereafter it is set by the program. 577C RWORK(5) should be positive. In order to travel in the 578C negative direction, see RWORK(6). 579C The default initial value equals 0.5*(RWORK(3)+RWORK(4)). 580C 581C RWORK(6) The local continuation direction, which is either +1.0 582C or -1.0 . This asserts that the program is moving in the 583C direction of increasing or decreasing values of the local 584C continuation variable, whose index is in IWORK(2). On first 585C call, the user must choose IWORK(2). Therefore, by setting 586C RWORK(6), the user may also specify whether the program is 587C to move initially to increase or decrease the variable whose 588C index is IWORK(2). 589C RWORK(6) defaults to +1. 590C 591C RWORK(7) A target value. It is only used if a target index 592C has been specified through IWORK(5). In that case, solution 593C points with the IWORK(5) component equal to RWORK(7) are 594C to be computed. The code will return each time it finds such 595C a point. RWORK(7) does not have a default value. The 596C program does not set it, and it is not referenced unless 597C IWORK(5) has been set. 598C 599C RWORK(8) EPMACH, the value of the machine precision. The computer 600C can distinguish 1.0+EPMACH from 1.0, but it cannot 601C distinguish 1.0+(EPMACH/2) from 1.0. This number is used 602C when estimating a reasonable accuracy request on a given 603C computer. PITCON computes a value for EPMACH internally. 604C 605C RWORK(9) Not currently used. 606C 607C RWORK(10) A minimum angle used in the steplength computation, 608C equal to 2.0*ARCCOS(1-EPMACH). 609C 610C RWORK(11) Estimate of the angle between the tangent vectors at the 611C last two continuation points. 612C 613C RWORK(12) The pseudo-arclength coordinate of the previous continuation 614C pointl; that is, the sum of the Euclidean distances between 615C all computed continuation points beginning with the start 616C point. Thus each new point should have a larger coordinate, 617C except for target and limit points which lie between the two 618C most recent continuation points. 619C 620C RWORK(13) Estimate of the pseudo-arclength coordinate of the current 621C continuation point. 622C 623C RWORK(14) Estimate of the pseudoarclength coordinate of the current 624C limit or target point, if any. 625C 626C RWORK(15) Size of the correction of the most recent continuation 627C point; that is, the maximum norm of the distance between the 628C predicted point and the accepted corrected point. 629C 630C RWORK(16) Estimate of the curvature between the last two 631C continuation points. 632C 633C RWORK(17) Sign of the determinant of the augmented matrix at the 634C last continuation point whose tangent vector has been 635C calculated. 636C 637C RWORK(18) Not currently used. 638C 639C RWORK(19) Not currently used. 640C 641C RWORK(20) Maximum growth factor for the predictor stepsize based 642C on the previous secant stepsize. The stepsize algorithm 643C will produce a suggested step that is no less that the 644C previous secant step divided by this factor, and no greater 645C than the previous secant step multiplied by that factor. 646C RWORK(20) defaults to 3. 647C 648C RWORK(21) The (Euclidean) secant distance between the last two 649C computed continuation points. 650C 651C RWORK(22) The previous value of RWORK(21). 652C 653C RWORK(23) A number judging the quality of the Newton corrector 654C convergence at the last continuation point. 655C 656C RWORK(24) Value of the component of the current tangent vector 657C corresponding to the current continuation index. 658C 659C RWORK(25) Value of the component of the previous tangent vector 660C corresponding to the current continuation index. 661C 662C RWORK(26) Value of the component of the current tangent vector 663C corresponding to the limit index in IWORK(6). 664C 665C RWORK(27) Value of the component of the previous tangent vector 666C corresponding to the limit index in IWORK(6). 667C 668C RWORK(28) Value of RWORK(7) when the last target point was 669C computed. 670C 671C RWORK(29) Sign of the determinant of the augmented matrix at the 672C previous continuation point whose tangent vector has been 673C calculated. 674C 675C RWORK(30) through RWORK(30+4*NVAR-1) are used by the program to hold 676C an old and new continuation point, a tangent vector and a 677C work vector. Subsequent entries of RWORK are used by the 678C linear solver. 679C 680C 681C G. Programming Notes 682C 683C 684C The minimal input and user routines required to apply the program are as 685C follows: 686C Write a function routine FX of the form described above. 687C Use DENSLV as the linear equation solver by setting SLVNAM to DENSLV. 688C Skip writing a Jacobian routine by using the finite difference option. 689C Pass the name of FX as the Jacobian name as well. 690C Declare the name of the function FX as EXTERNAL. 691C Set NVAR in accordance with your problem. 692C 693C Then: 694C 695C Dimension the vector IWORK to the size LIW=29+NVAR. 696C Dimension the vector RWORK to the size LRW=29+NVAR*(NVAR+6). 697C Dimension IPAR(1) and FPAR(1) as required in the function routine. 698C Dimension XR(NVAR) and set it to an approximate solution of F(XR)=0. 699C 700C Set the work arrays as follows: 701C 702C Initialize IWORK to 0 and RWORK to 0.0. 703C 704C Set IWORK(1)=0 (Problem startup) 705C Set IWORK(7)=3 (Maximum internally generated output) 706C Set IWORK(9)=1 (Forward difference Jacobian) 707C 708C Now call the program repeatedly, and never change any of its arguments. 709C Check IERROR to decide whether the code is working satisfactorily. 710C Print out the vector XR to see the current solution point. 711C 712C The most obvious input to try to set appropriately after some practice 713C would be the error tolerances RWORK(1) and RWORK(2), the minimum, maximum 714C and initial stepsizes RWORK(3), RWORK(4) and RWORK(5), and the initial 715C continuation index IWORK(2). 716C 717C For speed and efficiency, a Jacobian routine should be written. It can be 718C checked by comparing its results with the finite difference Jacobian. 719C 720C For a particular problem, the target and limit point input can be very 721C useful. For instance, in the case of a discretized boundary value problem 722C with a real parameter it may be desirable to compare the computed solutions 723C for different discretization-dimensions and equal values of the parameter. 724C For this the target option can be used with the prescribed values of the 725C parameter. Limit points usually are of importance in connection with 726C stability considerations. 727C 728C 729C H. References 730C 731C 732C 1. 733C Werner Rheinboldt, 734C Solution Field of Nonlinear Equations and Continuation Methods, 735C SIAM Journal of Numerical Analysis, 736C Volume 17, 1980, pages 221-237. 737C 738C 2. 739C Cor den Heijer and Werner Rheinboldt, 740C On Steplength Algorithms for a Class of Continuation Methods, 741C SIAM Journal of Numerical Analysis, 742C Volume 18, 1981, pages 925-947. 743C 744C 3. 745C Werner Rheinboldt, 746C Numerical Analysis of Parametrized Nonlinear Equations 747C John Wiley and Sons, New York, 1986 748C 749C 4. 750C Werner Rheinboldt and John Burkardt, 751C A Locally Parameterized Continuation Process, 752C ACM Transactions on Mathematical Software, 753C Volume 9, Number 2, June 1983, pages 215-235. 754C 755C 5. 756C Werner Rheinboldt and John Burkardt, 757C Algorithm 596, A Program for a Locally Parameterized Continuation Process, 758C ACM Transactions on Mathematical Software, 759C Volume 9, Number 2, June 1983, Pages 236-241. 760C 761C 6. 762C J J Dongarra, J R Bunch, C B Moler and G W Stewart, 763C LINPACK User's Guide, 764C Society for Industrial and Applied Mathematics, 765C Philadelphia, 1979. 766C 767C 7. 768C Richard Brent, 769C Algorithms for Minimization without Derivatives, 770C Prentice Hall, 1973. 771C 772C 8. 773C Tony Chan, 774C Deflated Decomposition of Solutions of Nearly Singular Systems, 775C Technical Report 225, 776C Computer Science Department, 777C Yale University, 778C New Haven, Connecticut, 06520, 779C 1982. 780C 781C 782C I. Converting between REAL and DOUBLE PRECISION versions 783C 784C 785C PITCON was written to use single precision floating point arithmetic, and 786C for most problems, this has sufficed, even on machines for which the 787C default REAL variable is stored in 32 bits. However, if PITCON cannot 788C solve a problem, the failure is sometimes caused by poor scaling of the 789C functions, a near-singular jacobian, or similar conditions that are 790C sometimes helped by using higher precision. 791C 792C The current version of the code is fairly easy to convert. The type of 793C each variable is declared in every routine. Generic FORTRAN functions 794C are used. All constants are declared in PARAMETER statements at the 795C beginning of each routine. 796C 797C To create a DOUBLE PRECISION version from a REAL version, a user should 798C make the following global substitutions: 799C 800C Replace by 801C 802C REAL DOUBLE PRECISION or REAL*8 803C 804C ISAMAX IDAMAX 805C SAXPY DAXPY 806C SCOPY DCOPY 807C SDOT DDOT 808C SGBDI DGBDI 809C SGBFA DGBFA 810C SGBSL DGBSL 811C SGEDI DGEDI 812C SGEFA DGEFA 813C SGESL DGESL 814C SNRM2 DNRM2 815C SSCAL DSCAL 816C SSWAP DSWAP 817C 818C and to convert from DOUBLE PRECISION to REAL, the changes should be 819C reversed. 820C 821C The same changes can be made to the sample programs. In most cases, this 822C should be sufficient to create proper double precision programs. 823C 824C 825C J. A sample calling program 826C 827C 828C The following is a sample program that calls PITCON to handle the 829C Freudenstein-Roth function. 830C 831C 832C C PCPRB1.FOR 22 February 1991 833C C 834C PROGRAM PCPRB1 835C C 836C C The Freudenstein-Roth function. 837C C 838C C For a more complicated version of this example, see PCPRB8. 839C C 840C C Reference 841C C 842C C F Freudenstein, B Roth, 843C C Numerical Solutions of Nonlinear Equations, 844C C Journal of the Association for Computing Machinery, 845C C Volume 10, 1963, Pages 550-556. 846C C 847C C The function F(X) is of the form 848C C 849C C FX(1) = X1 - X2**3 + 5*X2**2 - 2*X2 - 13 + 34*(X3-1) 850C C FX(2) = X1 + X2**3 + X2**2 - 14*X2 - 29 + 10*(X3-1) 851C C 852C C Starting from the point (15,-2,0), the program is required to produce 853C C solution points along the curve until it reaches a solution point 854C C (*,*,1). It also may be requested to look for limit points in the 855C C first or third components. 856C C 857C C The correct value of the solution at X3=1 is (5,4,1). 858C C 859C C Limit points in the first variable occur at: 860C C 861C C (14.28309, -1.741377, 0.2585779) 862C C (61.66936, 1.983801, -0.6638797) 863C C 864C C Limit points for the third variable occur at: 865C C 866C C (20.48586, -0.8968053, 0.5875873) 867C C (61.02031, 2.230139, -0.6863528) 868C C 869C C 870C INTEGER LIW 871C INTEGER LRW 872C INTEGER NVAR 873C C 874C PARAMETER (NVAR=3) 875C PARAMETER (LIW=NVAR+29) 876C PARAMETER (LRW=29+(6+NVAR)*NVAR) 877C C 878C EXTERNAL FXROTH 879C EXTERNAL DFROTH 880C EXTERNAL DENSLV 881C EXTERNAL PITCON 882C C 883C REAL FPAR(1) 884C INTEGER I 885C INTEGER IERROR 886C INTEGER IPAR(1) 887C INTEGER IWORK(LIW) 888C INTEGER J 889C CHARACTER NAME*12 890C REAL RWORK(LRW) 891C REAL XR(NVAR) 892C C 893C C Set work arrays to zero: 894C C 895C DO 10 I=1,LIW 896C IWORK(I)=0 897C 10 CONTINUE 898C DO 20 I=1,LRW 899C RWORK(I)=0.0 900C 20 CONTINUE 901C C 902C C Set some entries of work arrays. 903C C 904C C IWORK(1)=0 ; This is a startup 905C C IWORK(2)=2 ; Use X(2) for initial parameter 906C C IWORK(3)=0 ; Program may choose parameter index 907C C IWORK(4)=0 ; Update jacobian every newton step 908C C IWORK(5)=3 ; Seek target values for X(3) 909C C IWORK(6)=1 ; Seek limit points in X(1) 910C C IWORK(7)=1 ; Control amount of output. 911C C IWORK(8)=6 ; Output unit 912C C IWORK(9)=2 ; Jacobian choice. 913C C 914C IWORK(1)=0 915C IWORK(2)=2 916C IWORK(3)=0 917C IWORK(4)=0 918C IWORK(5)=3 919C IWORK(6)=1 920C IWORK(7)=1 921C IWORK(8)=6 922C IWORK(9)=0 923C C 924C C RWORK(1)=0.00001; Absolute error tolerance 925C C RWORK(2)=0.00001; Relative error tolerance 926C C RWORK(3)=0.01 ; Minimum stepsize 927C C RWORK(4)=20.0 ; Maximum stepsize 928C C RWORK(5)=0.3 ; Starting stepsize 929C C RWORK(6)=1.0 ; Starting direction 930C C RWORK(7)=1.0 ; Target value (Seek solution with X(3)=1) 931C C 932C RWORK(1)=0.00001 933C RWORK(2)=0.00001 934C RWORK(3)=0.01 935C RWORK(4)=20.0 936C RWORK(5)=0.3 937C RWORK(6)=1.0 938C RWORK(7)=1.0 939C C 940C C Set starting point. 941C C 942C XR(1)=15.0 943C XR(2)=-2.0 944C XR(3)=0.0 945C C 946C WRITE(6,*)' ' 947C WRITE(6,*)'PCPRB1:' 948C WRITE(6,*)' ' 949C WRITE(6,*)'PITCON sample program.' 950C WRITE(6,*)'Freudenstein-Roth function' 951C WRITE(6,*)'2 equations, 3 variables.' 952C WRITE(6,*)' ' 953C WRITE(6,*)'Step Type of Point '// 954C *'X(1) X(2) X(3)' 955C WRITE(6,*)' ' 956C I=0 957C NAME='Start point ' 958C WRITE(6,'(1X,I3,2X,A12,2X,3G14.6)')I,NAME,(XR(J),J=1,NVAR) 959C C 960C DO 30 I=1,30 961C C 962C CALL PITCON(DFROTH,FPAR,FXROTH,IERROR,IPAR,IWORK,LIW, 963C * NVAR,RWORK,LRW,XR,DENSLV) 964C C 965C IF(IWORK(1).EQ.1)THEN 966C NAME='Corrected ' 967C ELSEIF(IWORK(1).EQ.2)THEN 968C NAME='Continuation ' 969C ELSEIF(IWORK(1).EQ.3)THEN 970C NAME='Target point ' 971C ELSEIF(IWORK(1).EQ.4)THEN 972C NAME='Limit point ' 973C ENDIF 974C C 975C WRITE(6,'(1X,I3,2X,A12,2X,3G14.6)')I,NAME,(XR(J),J=1,NVAR) 976C C 977C IF(IWORK(1).EQ.3)THEN 978C WRITE(6,*)' ' 979C WRITE(6,*)'We have reached the point we wanted.' 980C WRITE(6,*)'The code may stop now.' 981C STOP 982C ENDIF 983C C 984C IF(IERROR.NE.0)THEN 985C WRITE(6,*)' ' 986C WRITE(6,*)'An error occurred.' 987C WRITE(6,*)'We will terminate the computation now.' 988C STOP 989C ENDIF 990C C 991C 30 CONTINUE 992C C 993C WRITE(6,*)' ' 994C WRITE(6,*)'PITCON did not reach the point of interest.' 995C STOP 996C END 997C SUBROUTINE FXROTH(NVAR,FPAR,IPAR,X,F,IERROR) 998C C 999C C Evaluate the function at X. 1000C C 1001C C ( X1 - ((X2-5.0)*X2+2.0)*X2 - 13.0 + 34.0*(X3-1.0) ) 1002C C ( X1 + ((X2+1.0)*X2-14.0)*X2 - 29.0 + 10.0*(X3-1.0) ) 1003C C 1004C REAL F(*) 1005C REAL FPAR(*) 1006C INTEGER IERROR 1007C INTEGER IPAR(*) 1008C REAL X(NVAR) 1009C C 1010C F(1)=X(1) 1011C * -((X(2)-5.0)*X(2)+2.0)*X(2) 1012C * -13.0 1013C * +34.0*(X(3)-1.0) 1014C C 1015C F(2)=X(1) 1016C * +((X(2)+1.0)*X(2)-14.0)*X(2) 1017C * -29.0 1018C * +10.0*(X(3)-1.0) 1019C C 1020C RETURN 1021C END 1022C SUBROUTINE DFROTH(NVAR,FPAR,IPAR,X,FJAC,IERROR) 1023C C 1024C C Evaluate the Jacobian: 1025C C 1026C C ( 1.0 (-3.0*X(2)+10.0)*X(2)-2.0 34.0 ) 1027C C ( 1.0 (3.0*X(2)+2.0)*X(2)-14.0 10.0 ) 1028C C 1029C REAL FJAC(NVAR,NVAR) 1030C REAL FPAR(*) 1031C INTEGER IERROR 1032C INTEGER IPAR(*) 1033C REAL X(NVAR) 1034C C 1035C FJAC(1,1)=1.0 1036C FJAC(1,2)=(-3.0*X(2)+10.0)*X(2)-2.0 1037C FJAC(1,3)=34.0 1038C C 1039C FJAC(2,1)=1.0 1040C FJAC(2,2)=(3.0*X(2)+2.0)*X(2)-14.0 1041C FJAC(2,3)=10.0 1042C C 1043C RETURN 1044C END 1045C 1046C 1047C Here is the output from the run of this sample problem 1048C 1049C 1050C PCPRB1: 1051C 1052C PITCON sample program. 1053C Freudenstein-Roth function 1054C 2 equations, 3 variables. 1055C 1056C Step Type of Point X(1) X(2) X(3) 1057C 1058C 0 Start point 15.0000 -2.00000 0.000000E+00 1059C 1060C PITCON 6.1 1061C University of Pittsburgh Continuation Code 1062C 1063C Last modification on 25 February 1991 1064C 1065C 1 Corrected 15.0000 -2.00000 0.000000E+00 1066C 2 Continuation 14.7105 -1.94205 0.653814E-01 1067C 3 Continuation 14.2846 -1.72915 0.268742 1068C 4 Limit point 14.2831 -1.74138 0.258578 1069C 5 Continuation 16.9061 -1.20941 0.546845 1070C 6 Continuation 24.9179 -0.599064 0.555136 1071C 7 Continuation 44.8783 0.487670 0.595261E-01 1072C 8 Continuation 60.0889 1.57585 -0.542365 1073C 9 Continuation -11.1752 4.23516 1.55667 1074C 10 Target point 5.00000 4.00000 1.00000 1075C 1076C We have reached the point we wanted. 1077C The code may stop now. 1078C FORTRAN STOP 1079C 1080 DOUBLE PRECISION EIGHT 1081 DOUBLE PRECISION HALF 1082 DOUBLE PRECISION HUNDRD 1083 DOUBLE PRECISION ONE 1084 DOUBLE PRECISION THREE 1085 DOUBLE PRECISION TWO 1086 DOUBLE PRECISION ZERO 1087C 1088 PARAMETER (EIGHT=8.0) 1089 PARAMETER (HALF=0.5) 1090 PARAMETER (HUNDRD=100.0) 1091 PARAMETER (ONE=1.0) 1092 PARAMETER (THREE=3.0) 1093 PARAMETER (TWO=2.0) 1094 PARAMETER (ZERO=0.0) 1095C 1096 EXTERNAL COQUAL 1097 EXTERNAL CORECT 1098 EXTERNAL DF 1099 EXTERNAL FX 1100 EXTERNAL IDAMAX 1101 EXTERNAL REPS 1102 EXTERNAL ROOT 1103 EXTERNAL DDOT 1104 EXTERNAL SETSTP 1105 EXTERNAL SLNAME 1106 EXTERNAL DNRM2 1107 EXTERNAL DSCAL 1108 EXTERNAL START 1109 EXTERNAL TANGNT 1110C 1111 INTRINSIC ABS 1112 INTRINSIC ATAN2 1113 INTRINSIC MAX 1114 INTRINSIC SIGN 1115 INTRINSIC SQRT 1116C 1117 INTEGER LIW 1118 INTEGER LRW 1119 INTEGER NVAR 1120C 1121 DOUBLE PRECISION A 1122 DOUBLE PRECISION ATCIPC 1123 DOUBLE PRECISION ATCJPC 1124 DOUBLE PRECISION B 1125 DOUBLE PRECISION DETS 1126 DOUBLE PRECISION DIRLPC 1127 DOUBLE PRECISION FA 1128 DOUBLE PRECISION FB 1129 DOUBLE PRECISION FPAR(*) 1130 DOUBLE PRECISION HTAN 1131 INTEGER I 1132 INTEGER ICALL 1133 INTEGER ICRIT 1134 INTEGER IERROR 1135 INTEGER IFLAG 1136 INTEGER IHOLD 1137 INTEGER IMITL 1138 INTEGER IPAR(*) 1139 INTEGER IPC 1140 INTEGER IPL 1141 INTEGER IDAMAX 1142 INTEGER IT 1143 INTEGER ITMP 1144 INTEGER IWORK(LIW) 1145 INTEGER IWRITE 1146 INTEGER JOB 1147 INTEGER JPC 1148 INTEGER LIM 1149 INTEGER LOUNIT 1150 INTEGER LPC 1151 INTEGER LTC 1152 INTEGER LTC1 1153 INTEGER LWK 1154 INTEGER LWK1 1155 INTEGER LXC 1156 INTEGER LXC1 1157 INTEGER LXF 1158 INTEGER LXF1 1159 INTEGER MLSTEP 1160 INTEGER MODSAV 1161 DOUBLE PRECISION REPS 1162 DOUBLE PRECISION RWORK(LRW) 1163 DOUBLE PRECISION SCIPL 1164 DOUBLE PRECISION DDOT 1165 DOUBLE PRECISION SKALE 1166 DOUBLE PRECISION SN 1167 DOUBLE PRECISION SNL 1168 DOUBLE PRECISION DNRM2 1169 DOUBLE PRECISION STEPX 1170 DOUBLE PRECISION TCIPC 1171 DOUBLE PRECISION TCIPL 1172 DOUBLE PRECISION TCOS 1173 DOUBLE PRECISION TEMP 1174 DOUBLE PRECISION TLIPC 1175 DOUBLE PRECISION TMAX 1176 DOUBLE PRECISION TMAX2 1177 DOUBLE PRECISION TMP 1178 DOUBLE PRECISION TMP1 1179 DOUBLE PRECISION TSIN 1180 DOUBLE PRECISION TSN 1181 DOUBLE PRECISION XABS 1182 DOUBLE PRECISION XDIF 1183 DOUBLE PRECISION XLOW 1184 DOUBLE PRECISION XR(NVAR) 1185 DOUBLE PRECISION XUP 1186C 1187 SAVE ICALL 1188C 1189 DATA ICALL /0/ 1190C 1191 ICALL=ICALL+1 1192C 1193C********************************************************************** 1194C 1195C 1. Preparations. 1196C 1197C Check entries of IWORK and RWORK, compute some constants. 1198C For negative values of IWORK(1), compare user and finite difference 1199C jacobian and return. 1200C For IWORK(1)=0, set up starting data, check starting point and return. 1201C For positive IWORK(1), prepare to compute next point. 1202C 1203C********************************************************************** 1204C 1205 IERROR=0 1206 IF(IWORK(8).LT.1.OR.IWORK(8).GT.99)IWORK(8)=6 1207 LOUNIT=IWORK(8) 1208 IF(IWORK(7).LT.0.OR.IWORK(7).GT.3)IWORK(7)=1 1209 IWRITE=IWORK(7) 1210C 1211C Initialization for first call 1212C 1213 IF(IWORK(1).EQ.0.OR.ICALL.EQ.1)THEN 1214 IF(IWRITE.GE.1)THEN 1215 WRITE(LOUNIT,*)' ' 1216 WRITE(LOUNIT,*)'PITCON 6.1' 1217 WRITE(LOUNIT,*)'University of Pittsburgh Continuation Code' 1218 WRITE(LOUNIT,*)' ' 1219 WRITE(LOUNIT,*)'Last modification on 27 February 1991' 1220 WRITE(LOUNIT,*)' ' 1221 ENDIF 1222 DO 10 I=11,17 1223 RWORK(I)=ZERO 122410 CONTINUE 1225 DO 20 I=21,29 1226 RWORK(I)=ZERO 122720 CONTINUE 1228 IWORK(10)=0 1229 IWORK(11)=0 1230 IWORK(12)=0 1231 DO 30 I=18,28 1232 IWORK(I)=0 123330 CONTINUE 1234C 1235C This is the only place REPS is called. It might be replaced by 1236C a call to the SLATEC machine parameter function R1MACH: 1237C 1238C RWORK(8)=R1MACH(4) 1239C 1240C or by simply setting RWORK(8) to the "machine epsilon", usually 1241C taken as the smallest power of 2, EPMACH, such that (1.0+EPMACH) is 1242C greater than 1.0, but (1.0+(EPMACH/2.0)) equals 1.0. 1243C 1244 RWORK(8)=REPS() 1245 IF(IWRITE.GE.2)WRITE(LOUNIT,919)RWORK(8) 1246 TCOS=SQRT(ONE-RWORK(8)) 1247 TSIN=SQRT(RWORK(8)) 1248 RWORK(10)=TWO*ATAN2(TSIN,TCOS) 1249 IF(NVAR.LE.1)THEN 1250 IERROR=1 1251 IF(IWRITE.GE.1)WRITE(LOUNIT,*) 1252 * 'PITCON - Number of variables must be at least 2.' 1253 GO TO 900 1254 ENDIF 1255 ENDIF 1256C 1257 LXC=29 1258 LXC1=LXC+1 1259 LXF=LXC+NVAR 1260 LXF1=LXF+1 1261 LTC=LXF+NVAR 1262 LTC1=LTC+1 1263 LWK=LTC+NVAR 1264 LWK1=LWK+1 1265 IWORK(13)=30 1266 IWORK(14)=LIW 1267 IWORK(15)=LWK+NVAR+1 1268 IWORK(16)=LRW 1269 IF(LIW.LT.IWORK(13))THEN 1270 WRITE(LOUNIT,*)'PITCON - Input value of LIW=',LIW 1271 WRITE(LOUNIT,*)' Minimal acceptable value=',IWORK(13) 1272 IERROR=1 1273 GO TO 900 1274 ELSEIF(LRW.LT.IWORK(15))THEN 1275 WRITE(LOUNIT,*)'PITCON - Input value of LRW=',LRW 1276 WRITE(LOUNIT,*)' Minimal acceptable value=',IWORK(15) 1277 IERROR=1 1278 GO TO 900 1279 ENDIF 1280C 1281C Check entries of IWORK 1282C 1283 IF(IWORK(2).LT.1.OR.IWORK(2).GT.NVAR)IWORK(2)=NVAR 1284 IPC=IWORK(2) 1285 IF(IWORK(3).NE.1)IWORK(3)=0 1286 IF(IWORK(4).LT.0.OR.IWORK(4).GT.2)IWORK(4)=0 1287 IF(IWORK(5).LT.1.OR.IWORK(5).GT.NVAR)IWORK(5)=0 1288 IT=IWORK(5) 1289 IF(IWORK(6).LT.1.OR.IWORK(6).GT.NVAR)IWORK(6)=0 1290 LIM=IWORK(6) 1291 IF(IWORK(9).LT.0.OR.IWORK(9).GT.2)IWORK(9)=0 1292 IF(IWORK(17).LT.1)IWORK(17)=10 1293C 1294C Check entries of RWORK 1295C 1296 IF(RWORK(1).LE.ZERO)RWORK(1)=SQRT(RWORK(8)) 1297 IF(RWORK(2).LE.ZERO)RWORK(2)=SQRT(RWORK(8)) 1298 IF(RWORK(3).LT.SQRT(RWORK(8)))RWORK(3)=SQRT(RWORK(8)) 1299 IF(RWORK(4).LT.RWORK(3))THEN 1300 TEMP=NVAR 1301 RWORK(4)=MAX(RWORK(3),SQRT(TEMP)) 1302 ENDIF 1303 IF(RWORK(5).LT.ZERO)THEN 1304 RWORK(5)=-RWORK(5) 1305 RWORK(6)=-RWORK(6) 1306 ENDIF 1307 IF(RWORK(5).LT.RWORK(3).OR.RWORK(5).GT.RWORK(4))THEN 1308 RWORK(5)=HALF*(RWORK(3)+RWORK(4)) 1309 ENDIF 1310 HTAN=RWORK(5) 1311 IF(RWORK(6).NE.(-ONE))RWORK(6)=ONE 1312 IF(RWORK(20).LT.ONE.OR.RWORK(20).GT.HUNDRD)RWORK(20)=THREE 1313C 1314C Check user Jacobian routine versus finite difference calculation. 1315C 1316 IF(IWORK(1).LT.0)THEN 1317 JOB=3 1318 CALL SLNAME(DETS,FX,DF,FPAR,IERROR,IPC,IPAR,IWORK,LIW, 1319 1 JOB,NVAR,RWORK,LRW,XR,RWORK(LWK1)) 1320 GO TO 900 1321 ENDIF 1322C 1323C Starting point check 1324C 1325 IF(IWORK(1).EQ.0)THEN 1326 CALL START(DF,FPAR,FX,IERROR,IPAR,IPC,IWORK,IWRITE,LIW, 1327 * LOUNIT,LRW,NVAR,RWORK,RWORK(LTC1),RWORK(LWK1),RWORK(LXC1), 1328 * RWORK(LXF1),XR,SLNAME) 1329 GO TO 900 1330 ENDIF 1331C 1332C*********************************************************************** 1333C 1334C 2. Target point 1335C 1336C If (IT.NE.0) target points are sought. Check to see if target component 1337C IT has value XIT lying between XC(IT) and XF(IT). If so, get linearly 1338C interpolated starting point, and use Newton's method to get target point. 1339C 1340C*********************************************************************** 1341C 1342 IF(IT.EQ.0.OR.IWORK(10).LE.1)GO TO 300 1343 IF( (IWORK(1).EQ.3) .AND. (RWORK(7).EQ.RWORK(28)) .AND. 1344 * (IT.EQ.IWORK(11)) ) GO TO 300 1345 MODSAV=IWORK(4) 1346210 CONTINUE 1347 XLOW=RWORK(LXC+IT) 1348 XUP=RWORK(LXF+IT) 1349 IF(RWORK(7).LT.XLOW.AND.RWORK(7).LT.XUP)GO TO 300 1350 IF(RWORK(7).GT.XLOW.AND.RWORK(7).GT.XUP)GO TO 300 1351 IF(IWRITE.GE.2)WRITE(LOUNIT,*) 1352 *'PITCON - Attempt correction of approximate target point.' 1353C 1354C Approximate the target point using the bracketing solutions. 1355C 1356 IF(XLOW.NE.XUP)THEN 1357 SKALE=(RWORK(7)-XLOW)/(XUP-XLOW) 1358 ELSE 1359 SKALE=ONE 1360 ENDIF 1361 CALL DCOPY(NVAR,RWORK(LXF1),1,XR,1) 1362 CALL DSCAL(NVAR,SKALE,XR,1) 1363 SKALE=ONE-SKALE 1364 CALL DAXPY(NVAR,SKALE,RWORK(LXC1),1,XR,1) 1365 XR(IT)=RWORK(7) 1366C 1367C Call CORECT to compute the exact target point, holding index IT fixed. 1368C 1369 ICRIT=0 1370 CALL CORECT(DF,FPAR,FX,IERROR,IT,IPAR,IWORK, 1371 1 NVAR,RWORK,TMP,RWORK(LWK1),XR,LRW,LIW,ICRIT,SLNAME) 1372 IWORK(24)=IWORK(24)+IWORK(28) 1373 IF(IERROR.NE.0.AND.IWORK(4).GT.0)THEN 1374 IERROR=0 1375 IWORK(4)=IWORK(4)-1 1376 IF(IWRITE.GE.1)WRITE(LOUNIT,1080)IWORK(4) 1377 GO TO 210 1378 ENDIF 1379 IWORK(4)=MODSAV 1380 IF(IERROR.NE.0)THEN 1381 WRITE(LOUNIT,*)'PITCON - Target point calculation failed.' 1382 GO TO 900 1383 ENDIF 1384C 1385C Record the values of IT and XIT, and compute the arclength to the target 1386C point. 1387C 1388 IWORK(1)=3 1389 IWORK(11)=IT 1390 RWORK(28)=RWORK(7) 1391 SKALE=-ONE 1392 CALL DCOPY(NVAR,XR,1,RWORK(LWK1),1) 1393 CALL DAXPY(NVAR,SKALE,RWORK(LXC1),1,RWORK(LWK1),1) 1394 RWORK(14)=RWORK(12)+DNRM2(NVAR,RWORK(LWK1),1) 1395 IWORK(27)=IWORK(27)+1 1396 GO TO 900 1397C 1398C*********************************************************************** 1399C 1400C 3. Tangent and local continuation parameter calculation. 1401C 1402C Unless the tangent and limit point calculations were already performed, 1403C (because the loop was interrupted for a limit point calculation), 1404C set up and solve the equation for the tangent vector. 1405C 1406C Force the tangent vector to be of unit length, and try to preserve 1407C the "sense" or "direction" of the curve by forcing the IPL-th component 1408C of the tangent vector to agree in sign with the IPL-th component of the 1409C previous secant vector. (On the first step, however, we have to use 1410C the user's input direction to choose the sign). 1411C 1412C Set the local continuation parameter IPC as the index of the 1413C component of greatest magnitude in the tangent vector, unless a 1414C limit point appears to be coming in that direction, and another 1415C choice is available. 1416C 1417C*********************************************************************** 1418C 1419 300 CONTINUE 1420 IF(IWORK(10).EQ.3) GO TO 600 1421C 1422C STORE OLD TANGENT IN WORK ARRAY RWORK(LWK1)-RWORK(LWK+NVAR), COMPUTE 1423C NEW TANGENT FOR CURRENT POINT XF IN RWORK(LXF1)-RWORK(LXF+NVAR). 1424C THE TANGENT IS CALLED TC AND STORED IN RWORK(LTC1)-RWORK(LTC+NVAR). 1425C NOTE THAT FOR A NORMAL RETURN THE CURRENT POINT XF WILL BE STORED 1426C AS XC IN RWORK(LXC1)-RWORK(LXC+NVAR) AND HENCE THAT TC WILL THEN 1427C CORRESPOND TO THE POINT XC 1428C 1429 IPL=IPC 1430 CALL DCOPY(NVAR,RWORK(LTC1),1,RWORK(LWK1),1) 1431 RWORK(29)=RWORK(17) 1432 CALL TANGNT(RWORK(17),FX,DF,FPAR,IERROR,IPC,IPAR,IWORK,NVAR, 1433 1 RWORK,RWORK(LTC1),RWORK(LXF1),LIW,LRW,SLNAME) 1434 IF(IERROR.NE.0)THEN 1435 IF(IWRITE.GE.1)WRITE(LOUNIT,*) 1436 * 'PITCON - Tangent calculation failed.' 1437 GO TO 900 1438 ENDIF 1439C 1440C FIND LARGEST AND SECOND LARGEST COMPONENTS OF TANGENT 1441C 1442 IPC=IDAMAX(NVAR,RWORK(LTC1),1) 1443 TMAX=RWORK(LTC+IPC) 1444 RWORK(LTC+IPC)=ZERO 1445 JPC=IDAMAX(NVAR,RWORK(LTC1),1) 1446 RWORK(LTC+IPC)=TMAX 1447 TMAX=ABS(TMAX) 1448 TMAX2=ABS(RWORK(LTC+JPC)) 1449 IF(JPC.EQ.0)JPC=IPC 1450 IF(IWORK(3).EQ.1)JPC=IPC 1451C 1452C ADJUST SIGN OF THE TANGENT VECTOR. 1453C COMPARE THE SIGN OF THE COMPONENT TC(IPL) WITH THE 1454C SIGN OF XF(IPL)-XC(IPL). IF A TARGET OR LIMIT POINT 1455C HAS BEEN COMPUTED BETWEEN XC AND XF THEN USE XF(IPL)-XR(IPL) 1456C INSTEAD. AT THE FIRST STEP COMPARE WITH THE USER INPUT DIRECTION. 1457C IF THESE SIGNS DIFFER, CHANGE THE SIGN OF THE TANGENT VECTOR 1458C AND THE SIGN OF THE DETERMINANT. 1459C 1460 SCIPL=RWORK(6) 1461 TCIPL=RWORK(LTC+IPL) 1462 IF(IWORK(10).GT.1)THEN 1463 SCIPL=RWORK(LXF+IPL) 1464 TMP=RWORK(LXC+IPL) 1465 IF(IWORK(1).EQ.3.OR.IWORK(1).EQ.4)TMP=XR(IPL) 1466 SCIPL=SCIPL-TMP 1467 ENDIF 1468 IF(SIGN(ONE,TCIPL).NE.SIGN(ONE,SCIPL))THEN 1469 SKALE=-ONE 1470 CALL DSCAL(NVAR,SKALE,RWORK(LTC1),1) 1471 RWORK(17)=-RWORK(17) 1472 ENDIF 1473C 1474C Unless we are computing a starting point, record the new state. 1475C 1476 IF(IWORK(10).GT.1)IWORK(10)=3 1477C 1478C THE LOCATION OF THE LARGEST COMPONENT OF THE TANGENT VECTOR 1479C WILL BE USED FOR THE LOCAL PARAMETERIZATION OF THE CURVE UNLESS 1480C A LIMIT POINT IN IPC APPEARS TO BE COMING. 1481C TO CHECK THIS, WE COMPARE TCIPC=TC(IPC) AND THE 1482C SECOND LARGEST COMPONENT TCJPC=TC(JPC). IF TCJPC IS NO LESS 1483C THAN .1 OF TCIPC, AND TC(JPC) IS LARGER THAN TL(JPC), 1484C WHEREAS TC(IPC) IS LESS THAN TL(IPC), WE WILL RESET THE 1485C LOCAL PARAMETER IPC TO JPC. 1486C BUT IF NOT ALLOWED TO SWITCH PARAMETERS, IGNORE THE CHECK 1487C 1488 IF(IWORK(3).NE.1)THEN 1489 ATCIPC=ABS(RWORK(LTC+IPC)) 1490 ATCJPC=ABS(RWORK(LTC+JPC)) 1491 IF(JPC.NE.IPC.AND.IWORK(10).GT.1)THEN 1492 TLIPC=RWORK(LWK+IPC) 1493 TCIPC=RWORK(LTC+IPC) 1494 TMP=ABS(RWORK(LWK+JPC)) 1495 IF((SIGN(ONE,TCIPC).EQ.SIGN(ONE,TLIPC)) 1496 1 .AND.(ATCIPC.LT.ABS(TLIPC)) 1497 2 .AND.(ATCJPC.GE.MAX(0.1*ATCIPC,TMP)))THEN 1498 IF(IWRITE.GE.3)WRITE(LOUNIT,350)IPC 1499 ITMP=IPC 1500 IPC=JPC 1501 JPC=ITMP 1502 ENDIF 1503 ENDIF 1504 ENDIF 1505 IF(IWRITE.GE.3.AND.IWORK(3).EQ.0)THEN 1506 WRITE(LOUNIT,360)IPC 1507 WRITE(LOUNIT,370)JPC 1508 ENDIF 1509C 1510C RECORD VALUES OF THE COMPONENT TCIPC=TC(IPC) OF 1511C NEW TANGENT VECTOR, AS WELL AS OF THE CORRESPONDING COMPONENT 1512C TLIPC=WK(IPC) OF THE OLD TANGENT VECTOR. 1513C SET SIGN OF THE DETERMINANT. 1514C FOR USE IN THE LIMIT POINT COMPUTATION RECORD THE 1515C COMPONENT TC(LIM) OF THE NEW TANGENT VECTOR 1516C 1517 IWORK(2)=IPC 1518 IWORK(12)=JPC 1519 RWORK(24)=RWORK(LTC+IPC) 1520 RWORK(25)=RWORK(LWK+IPC) 1521 RWORK(6)=RWORK(17) 1522 IF(LIM.GT.0)THEN 1523 RWORK(27)=RWORK(26) 1524 RWORK(26)=RWORK(LTC+LIM) 1525 IF(IWRITE.GE.2)WRITE(LOUNIT,1020)RWORK(26) 1526 ENDIF 1527C 1528C COMPUTE ALPHLC, THE ANGLE BETWEEN OLD TANGENT AND TANGENT AT XC 1529C 1530 IF(IWORK(10).LE.1) GO TO 600 1531 TCOS=DDOT(NVAR,RWORK(LWK1),1,RWORK(LTC1),1) 1532 IF(TCOS.GT.ONE)TCOS=ONE 1533 IF(TCOS.LT.(-ONE))TCOS=-ONE 1534 TSIN=SQRT(ONE-TCOS*TCOS) 1535 RWORK(11)=ATAN2(TSIN,TCOS) 1536 IF((IWORK(10).LE.1).OR.(RWORK(17).EQ.RWORK(29)))GO TO 400 1537 IF(IWRITE.GE.2)WRITE(LOUNIT,*) 1538 *'PITCON - Possible bifurcation point detected.' 1539C 1540C*********************************************************************** 1541C 1542C 4. Limit point check. 1543C 1544C Skip this section if LIM=0. 1545C 1546C Otherwise, user has requested a search for limit points in a given 1547C index by setting LIM to a nonzero value. 1548C 1549C Compare LIM-th components of previous and current tangent vectors. 1550C If a sign difference occurs, we assume a limit point has been passed. 1551C Attempt to compute a point XR between the previous and current points, 1552C for which the LIM-th component of the tangent is zero. 1553C 1554C This search will be guided by a rootfinder. The scalar function 1555C to be zeroed out is the LIM-th tangent vector component. 1556C 1557C*********************************************************************** 1558C 1559 400 CONTINUE 1560 IF((LIM.LE.0).OR. 1561 * (IWORK(10).LE.1).OR. 1562 * (SIGN(ONE,RWORK(26)).EQ.SIGN(ONE,RWORK(27)))) GO TO 600 1563 IF(IWRITE.GE.2)WRITE(LOUNIT,*) 1564 *'PITCON - Attempt correction of approximate limit point.' 1565 MODSAV=IWORK(4) 1566410 CONTINUE 1567C 1568C SEE IF EITHER ENDPOINT CAN BE ACCEPTED AS LIMIT POINT. 1569C IF XC IS LIMIT POINT, WORK ALREADY CONTAINS TANGENT AT XC. 1570C 1571 IF(ABS(RWORK(27)).LE.HALF*RWORK(1)) THEN 1572 CALL DCOPY(NVAR,RWORK(LXC1),1,XR,1) 1573 GO TO 490 1574 ENDIF 1575 IF(ABS(RWORK(26)).LE.HALF*RWORK(1)) THEN 1576 CALL DCOPY(NVAR,RWORK(LXF1),1,XR,1) 1577 CALL DCOPY(NVAR,RWORK(LTC1),1,RWORK(LWK1),1) 1578 GO TO 490 1579 ENDIF 1580C 1581C If interval is extremely small, simply assign one 1582C endpoint of the interval as the answer. 1583C 1584 XDIF=ABS(RWORK(LXC+LIM)-RWORK(LXF+LIM)) 1585 XABS=MAX(ABS(RWORK(LXC+LIM)),ABS(RWORK(LXF+LIM))) 1586 IF(XDIF.LE.EIGHT*RWORK(8)*(ONE+XABS)) THEN 1587 IF(ABS(RWORK(27)).GT.ABS(RWORK(26))) THEN 1588 CALL DCOPY(NVAR,RWORK(LXF1),1,XR,1) 1589 CALL DCOPY(NVAR,RWORK(LTC1),1,RWORK(LWK1),1) 1590 ELSE 1591 CALL DCOPY(NVAR,RWORK(LXC1),1,XR,1) 1592 ENDIF 1593 GO TO 490 1594 ENDIF 1595C 1596C BEGIN ROOT-FINDING ITERATION WITH INTERVAL (0,1) AND 1597C FUNCTION VALUES TLLIM, TCLIM. 1598C 1599 A=ZERO 1600 FA=RWORK(27) 1601 B=ONE 1602 FB=RWORK(26) 1603C 1604C FIND LPC, INDEX OF MAXIMUM ENTRY OF SECANT (EXCEPT THAT 1605C LPC MUST NOT EQUAL LIM) AND SAVE SIGN IN DIRLPC 1606C SO THAT NEW TANGENTS MAY BE PROPERLY SIGNED. 1607C 1608 TEMP=RWORK(LXC+LIM) 1609 RWORK(LXC+LIM)=RWORK(LXF+LIM) 1610 CALL DCOPY(NVAR,RWORK(LXF1),1,XR,1) 1611 SKALE=-ONE 1612 CALL DAXPY(NVAR,SKALE,RWORK(LXC1),1,XR,1) 1613 LPC=IDAMAX(NVAR,XR,1) 1614 RWORK(LXC+LIM)=TEMP 1615 DIRLPC=SIGN(ONE,RWORK(LXF+LPC)-RWORK(LXC+LPC)) 1616C 1617C SET FIRST APPROXIMATION TO ROOT TO WHICHEVER ENDPOINT 1618C HAS SMALLEST LIM-TH COMPONENT OF TANGENT 1619C 1620 IF(ABS(RWORK(26)).GE.ABS(RWORK(27))) THEN 1621 SN=ZERO 1622 TSN=RWORK(27) 1623 CALL DCOPY(NVAR,RWORK(LXC1),1,XR,1) 1624 ELSE 1625 SN=ONE 1626 TSN=RWORK(26) 1627 CALL DCOPY(NVAR,RWORK(LXF1),1,XR,1) 1628 ENDIF 1629C 1630C CALL ROOTFINDER FOR APPROXIMATE ROOT SN. USE LINEAR COMBINATION 1631C OF PREVIOUS ROOT SNL, AND ONE OF 0.0 AND 1.0 TO GET A STARTING 1632C POINT FOR CORRECTION. RETURN TO CURVE ON LINE X(LPC)=CONSTANT, 1633C COMPUTE TANGENT THERE, AND SET FUNCTION VALUE TO TANGENT(LIM) 1634C 1635 MLSTEP=25 1636 IMITL=0 1637 IF(IWRITE.GE.2)THEN 1638 TMP=ZERO 1639 WRITE(LOUNIT,1100)TMP,RWORK(27) 1640 TMP=ONE 1641 WRITE(LOUNIT,1100)TMP,RWORK(26) 1642 ENDIF 1643441 CONTINUE 1644 SNL=SN 1645 CALL ROOT(A,FA,B,FB,SN,TSN,IMITL,IFLAG,IERROR,RWORK(8)) 1646 IWORK(23)=IWORK(23)+1 1647 IF(IERROR.NE.0)THEN 1648 IF(IWRITE.GE.1)WRITE(LOUNIT,*) 1649 * 'PITCON - Bad interval for rootfinder' 1650 GO TO 900 1651 ENDIF 1652 IF(IFLAG.EQ.(-1).OR.IFLAG.EQ.0)GO TO 490 1653C 1654C FIND WHETHER SN LIES IN (0.0,SNL) OR (SNL,1.0). 1655C USE APPROPRIATE LINEAR COMBINATION TO GET STARTING POINT. 1656C 1657 IF(SN.LE.SNL) THEN 1658C 1659C SET X(SN)=(SNL-SN)/(SNL-0.0)*X(0.0)+(SN-0.0)/(SNL-0.0)*X(SNL) 1660C 1661 IF(SNL.LE.ZERO)THEN 1662 SKALE=ZERO 1663 ELSE 1664 SKALE=SN/SNL 1665 ENDIF 1666 IF(SKALE.LE.ZERO)SKALE=ZERO 1667 IF(SKALE.GE.ONE)SKALE=ONE 1668 CALL DSCAL(NVAR,SKALE,XR,1) 1669 SKALE=ONE-SKALE 1670 CALL DAXPY(NVAR,SKALE,RWORK(LXC1),1,XR,1) 1671 ELSE 1672C 1673C SET X(SN)=(SN-SNL)/(1.0-SNL)*X(1.0)+(1.0-SN)/(1.0-SNL)*X(SNL) 1674C 1675 IF(SNL.GE.ONE)THEN 1676 SKALE=ZERO 1677 ELSE 1678 SKALE=(ONE-SN)/(ONE-SNL) 1679 ENDIF 1680 IF(SKALE.LE.ZERO)SKALE=ZERO 1681 IF(SKALE.GE.ONE)SKALE=ONE 1682 CALL DSCAL(NVAR,SKALE,XR,1) 1683 SKALE=ONE-SKALE 1684 CALL DAXPY(NVAR,SKALE,RWORK(LXF1),1,XR,1) 1685 ENDIF 1686C 1687582 CONTINUE 1688 ICRIT=0 1689 CALL CORECT(DF,FPAR,FX,IERROR,LPC,IPAR,IWORK, 1690 1 NVAR,RWORK,TMP,RWORK(LWK1),XR,LRW,LIW,ICRIT,SLNAME) 1691 IF(IERROR.NE.0.AND.IWORK(4).GT.0)THEN 1692 IWORK(4)=IWORK(4)-1 1693 CALL DCOPY(NVAR,RWORK(LXC1),1,XR,1) 1694 SKALE=ONE-SN 1695 CALL DSCAL(NVAR,SKALE,XR,1) 1696 SKALE=SN 1697 CALL DAXPY(NVAR,SKALE,RWORK(LXF1),1,XR,1) 1698 IF(IWRITE.GE.1)WRITE(LOUNIT,1090)IWORK(4) 1699 GO TO 582 1700 ENDIF 1701 IWORK(4)=MODSAV 1702 IF(IERROR.NE.0)THEN 1703 WRITE(LOUNIT,*) 1704 * 'PITCON - Corrector failed on limit point.' 1705 GO TO 900 1706 ENDIF 1707 CALL TANGNT(TEMP,FX,DF,FPAR,IERROR,LPC,IPAR,IWORK,NVAR, 1708 1 RWORK,RWORK(LWK1),XR,LIW,LRW,SLNAME) 1709 IF(IERROR.NE.0)THEN 1710 IF(IWRITE.GE.1)WRITE(LOUNIT,*) 1711 * 'PITCON - Tangent error in limit point calculation.' 1712 GO TO 900 1713 ENDIF 1714C 1715C Adjust the sign of the tangent vector so the LPC-th component 1716C has the same sign as the LPC-th component of the secant. 1717C 1718 IF(DIRLPC.NE.SIGN(ONE,RWORK(LWK+LPC))) THEN 1719 SKALE=-ONE 1720 CALL DSCAL(NVAR,SKALE,RWORK(LWK1),1) 1721 ENDIF 1722C 1723C SEE IF WE CAN ACCEPT THE NEW POINT BECAUSE TANGENT(LIM) IS SMALL 1724C OR MUST STOP BECAUSE 25 ITERATIONS TAKEN, 1725C OR IF WE MUST GO ON. 1726C 1727 TSN=RWORK(LWK+LIM) 1728 IF(IWRITE.GE.2)WRITE(LOUNIT,1100)SN,TSN 1729 IF(ABS(TSN).LE.RWORK(1))GO TO 490 1730 IF(IMITL.LT.MLSTEP)GO TO 441 1731C 1732C Limit point iteration has failed. We set an error flag, 1733C but we return the partial results of the abortive computation 1734C anyway. 1735C 1736 IERROR=8 1737 IF(IWRITE.GE.1)WRITE(LOUNIT,*) 1738 *'PITCON - Too many steps in limit point calculation.' 1739C 1740C Limit point iteration is over. 1741C Compute and store information. 1742C 1743 490 CONTINUE 1744 CALL DCOPY(NVAR,XR,1,RWORK(LWK1),1) 1745 SKALE=-ONE 1746 CALL DAXPY(NVAR,SKALE,RWORK(LXC1),1,RWORK(LWK1),1) 1747 RWORK(14)=RWORK(12)+DNRM2(NVAR,RWORK(LWK1),1) 1748 IWORK(27)=IWORK(27)+1 1749 IWORK(1)=4 1750 GO TO 900 1751C 1752C*********************************************************************** 1753C 1754C 5. Compute next predictor step length, HTAN. 1755C 1756C*********************************************************************** 1757C 1758 600 CONTINUE 1759 IF(IWORK(10).GT.1)CALL SETSTP(HTAN,IWORK,LIW,LRW,RWORK) 1760C 1761C*********************************************************************** 1762C 1763C 6. Continuation step 1764C 1765C Our current data is the current point XC, its tangent vector TC, and 1766C a steplength HTAN. We predict the location of the next point on the 1767C curve using the Euler approximation XR=XC+HTAN*TC. 1768C 1769C Newton iteration is applied to this point, to force it to lie on the 1770C curve. In order to make the system square, an augmenting equation 1771C is added to the system, specifying that XR(IPC)=XC(IPC)+HTAN*TC(IPC). 1772C (The right hand side is a constant.) 1773C 1774C If the Newton correction process fails, the stepsize is reduced and 1775C prediction and correction retried. Failure will most likely be 1776C signaled by repeated step reductions, until the minimum allowable 1777C stepsize is reached. If this occurs, PITCON has failed, and cannot 1778C proceed along the curve any more. 1779C 1780C*********************************************************************** 1781C 1782 700 CONTINUE 1783 IWORK(18)=0 1784 IPC=IWORK(2) 1785 JPC=IWORK(12) 1786 IHOLD=IPC 1787 MODSAV=IWORK(4) 1788710 CONTINUE 1789 IF(IWRITE.GE.2)WRITE(LOUNIT,715)HTAN 1790 CALL DCOPY(NVAR,RWORK(LXF1),1,XR,1) 1791 CALL DAXPY(NVAR,HTAN,RWORK(LTC1),1,XR,1) 1792 740 CONTINUE 1793 IF(IWRITE.GE.2.AND.IWORK(3).EQ.0)WRITE(LOUNIT,725)IHOLD 1794 ICRIT=0 1795 CALL CORECT(DF,FPAR,FX,IERROR,IHOLD,IPAR,IWORK, 1796 1 NVAR,RWORK,STEPX,RWORK(LWK1),XR,LRW,LIW,ICRIT,SLNAME) 1797 IWORK(25)=IWORK(25)+IWORK(28) 1798 IF(IERROR.EQ.0)GO TO 800 1799C 1800C Only VERY fatal errors should abort the correction process. 1801C Abort this process only after kicking and screaming. 1802C 1803 IF(IERROR.EQ.2)THEN 1804 WRITE(LOUNIT,*)'PITCON - Fatal error during Newton correction' 1805 WRITE(LOUNIT,*)' of a continutation point!' 1806 GO TO 900 1807 ENDIF 1808C 1809C FOR CASES WHERE CORRECTOR FAILED TO CONVERGE, 1810C AND THE VALUE OF IHOLD WAS NOT JPC, 1811C AND THE USER HAS NOT FIXED THE CONTINUATION PARAMETER 1812C KEEP THE SAME STEPSIZE BUT USE JPC INSTEAD OF IPC. 1813C 1814 IF(IWORK(3).NE.1.AND.IHOLD.NE.JPC.AND.JPC.NE.0)THEN 1815 IHOLD=JPC 1816 IF(IERROR.EQ.5)GO TO 740 1817 GO TO 710 1818 ENDIF 1819C 1820C IF JPC FAILS AS IHOLD VALUE, RETURN TO IPC. 1821C 1822 IHOLD=IPC 1823C 1824C IF WE ARE USING SOME FORM OF MODIFIED NEWTON, BUT 1825C WE HAVE REACHED MINIMUM STEPSIZE OR 1826C HAVE HAD TWO FAILURES IN A ROW ON THIS STEP, 1827C USE A BETTER NEWTON METHOD 1828C 1829 IF(IWORK(4).GT.0.AND. 1830 * (HTAN.LT.RWORK(20)*RWORK(3).OR.IWORK(18).GE.2) )THEN 1831 IWORK(4)=IWORK(4)-1 1832 WRITE(LOUNIT,1130)IWORK(4) 1833 IF(IERROR.NE.5)GO TO 710 1834 GO TO 740 1835 ENDIF 1836C 1837C NO CONVERGENCE, TRY A SMALLER STEPSIZE 1838C 1839 IF(HTAN.LT.RWORK(20)*RWORK(3))THEN 1840 IERROR=4 1841 IF(IWRITE.GE.1)WRITE(LOUNIT,*) 1842 * 'PITCON - Stepsize fell below minimum.' 1843 GO TO 900 1844 ENDIF 1845 HTAN=HTAN/RWORK(20) 1846 IWORK(18)=IWORK(18)+1 1847 IF(IERROR.NE.5)GO TO 710 1848 SKALE=ONE/RWORK(20) 1849 CALL DSCAL(NVAR,SKALE,XR,1) 1850 SKALE=ONE-SKALE 1851 CALL DAXPY(NVAR,SKALE,RWORK(LXF1),1,XR,1) 1852 GO TO 740 1853C 1854C*********************************************************************** 1855C 1856C 7. Successful step. Update information. 1857C 1858C*********************************************************************** 1859C 1860 800 CONTINUE 1861 IWORK(1)=2 1862 IF(IWORK(3).NE.1)IWORK(2)=IHOLD 1863 IWORK(4)=MODSAV 1864 IWORK(10)=2 1865 IWORK(26)=IWORK(26)+IWORK(18) 1866 IWORK(27)=IWORK(27)+1 1867 RWORK(5)=HTAN 1868 RWORK(22)=RWORK(21) 1869 IF(IWRITE.GE.1.AND.IWORK(18).GT.0)WRITE(LOUNIT,755)IWORK(18) 1870 CALL DCOPY(NVAR,XR,1,RWORK(LWK1),1) 1871 SKALE=-ONE 1872 CALL DAXPY(NVAR,SKALE,RWORK(LXF1),1,RWORK(LWK1),1) 1873 RWORK(21)=DNRM2(NVAR,RWORK(LWK1),1) 1874 RWORK(12)=RWORK(13) 1875 RWORK(13)=RWORK(12)+RWORK(21) 1876 RWORK(14)=RWORK(13) 1877C 1878C UPDATE VALUES OF COMPUTED POINTS. STORE OLD POINT IN 1879C RWORK(LXC1)-RWORK(LXC+NVAR) AND NEW POINT IN RWORK(LXF1)- 1880C RWORK(LXF+NVAR). COMPUTE AND STORE CORDIS, THE MAXIMUM NORM 1881C OF THE CORRECTOR DISTANCE. 1882C 1883 RWORK(15)=ZERO 1884 IF(IWORK(28).NE.0)THEN 1885 DO 810 I=1,NVAR 1886 TMP1=XR(I)-(RWORK(LXF+I)+HTAN*RWORK(LTC+I)) 1887 IF(ABS(TMP1).GT.RWORK(15))RWORK(15)=ABS(TMP1) 1888 810 CONTINUE 1889 ENDIF 1890 CALL DCOPY(NVAR,RWORK(LXF1),1,RWORK(LXC1),1) 1891 CALL DCOPY(NVAR,XR,1,RWORK(LXF1),1) 1892C 1893C COMPUTE CONVERGENCE QUALITY 1894C 1895 CALL COQUAL(STEPX,IWORK,LIW,RWORK,LRW) 1896C 1897C********************************************************************** 1898C 1899C 8. All returns to calling program should exit here. 1900C 1901C********************************************************************** 1902C 1903 900 CONTINUE 1904 IF(IERROR.EQ.0)THEN 1905 ELSEIF(IERROR.EQ.1)THEN 1906 WRITE(LOUNIT,*)'PITCON - FATAL ERROR.' 1907 WRITE(LOUNIT,*)' Storage or data error.' 1908 ELSEIF(IERROR.EQ.2)THEN 1909 WRITE(LOUNIT,*)'PITCON - FATAL ERROR.' 1910 WRITE(LOUNIT,*) 1911 * ' User-defined function or jacobian error.' 1912 ELSEIF(IERROR.EQ.3)THEN 1913 WRITE(LOUNIT,*)'PITCON - FATAL ERROR.' 1914 WRITE(LOUNIT,*)' Solver failed.' 1915 ELSEIF(IERROR.EQ.4)THEN 1916 WRITE(LOUNIT,*)'PITCON - FATAL ERROR.' 1917 WRITE(LOUNIT,*)' Iteration failed.' 1918 ELSEIF(IERROR.EQ.5)THEN 1919 WRITE(LOUNIT,*)'PITCON - FATAL ERROR.' 1920 WRITE(LOUNIT,*)' Too many corrector steps.' 1921 ELSEIF(IERROR.EQ.6)THEN 1922 WRITE(LOUNIT,*)'PITCON - FATAL ERROR.' 1923 WRITE(LOUNIT,*)' Zero tangent vector.' 1924 ELSEIF(IERROR.EQ.7)THEN 1925 WRITE(LOUNIT,*)'PITCON - Nonfatal error.' 1926 WRITE(LOUNIT,*)' Root finder failed on limit point.' 1927 ELSEIF(IERROR.EQ.8)THEN 1928 WRITE(LOUNIT,*)'PITCON - Nonfatal error.' 1929 WRITE(LOUNIT,*)' Too many limit point steps.' 1930 ELSE 1931 IERROR=9 1932 WRITE(LOUNIT,*)'PITCON - FATAL ERROR.' 1933 WRITE(LOUNIT,*)' Unidentified error.' 1934 ENDIF 1935 RETURN 1936350 FORMAT(' PITCON - Expect limit point in index ',I5) 1937360 FORMAT(' PITCON - Continuation index: First choice',I5) 1938370 FORMAT(' PITCON - Second choice',I5) 1939715 FORMAT(' PITCON - Predictor using stepsize ',G14.6) 1940725 FORMAT(' PITCON - Corrector fixing index',I5) 1941755 FORMAT(' PITCON - Predictor stepsize was reduced ',I6,' times.') 1942919 FORMAT(' PITCON - Machine epsilon=',G14.6) 19431020 FORMAT(' PITCON - Tangent vector has limit component =',G14.6) 19441080 FORMAT(' PITCON - Retry target computation with IWORK(4)=',I6) 19451090 FORMAT(' PITCON - Retry limit computation with IWORK(4)=',I6) 19461100 FORMAT(' PITCON - For S=',G14.6,' TAN(X(S))(LIM)=',G14.6) 19471130 FORMAT(' PITCON - Retrying step with IWORK(4)=',I6) 1948 END 1949 SUBROUTINE COQUAL(STEPX,IWORK,LIW,RWORK,LRW) 1950C 1951C*********************************************************************** 1952C 1953C COQUAL computes the factor QUAL which is based on the 'quality' 1954C (number of steps, last step/total step) of the corrector iteration. 1955C This factor is used as part of the stepsize algorithm. 1956C 1957C See the paper "On Steplength Algorithms for a Class of Continuation 1958C Methods", listed in the documentation. 1959C 1960C STEPX = Input, size of the last correction. 1961C CORDIS = Input, distance between first and last point. 1962C MAXCOR = Input, maximum number of corrections allowed. 1963C MODNEW = Input, defines type of corrector process used. 1964C NCOR = Input, actual number of corrections used. 1965C QUAL = Output, the quality factor, between 1/8 and 8, 1966C with 1/8 poor, 1 average, 8 superior. Used to 1967C modify the next stepsize. 1968C 1969 DOUBLE PRECISION EIGHT 1970 DOUBLE PRECISION ONE 1971 DOUBLE PRECISION THGIE 1972C 1973 PARAMETER (EIGHT=8.0) 1974 PARAMETER (ONE=1.0) 1975 PARAMETER (THGIE=0.125) 1976C 1977 INTEGER LIW 1978 INTEGER LRW 1979C 1980 DOUBLE PRECISION BASE 1981 DOUBLE PRECISION BOT 1982 DOUBLE PRECISION CORDIS 1983 DOUBLE PRECISION EPSATE 1984 DOUBLE PRECISION ESAB 1985 DOUBLE PRECISION EXPO 1986 INTEGER IBOT 1987 INTEGER ITOP 1988 INTEGER IWORK(LIW) 1989 INTEGER LOUNIT 1990 INTEGER MAXCOR 1991 INTEGER MODNEW 1992 INTEGER NAVE 1993 INTEGER NMAX 1994 DOUBLE PRECISION QUAL 1995 DOUBLE PRECISION RWORK(LRW) 1996 DOUBLE PRECISION STEPX 1997 DOUBLE PRECISION TERM 1998 DOUBLE PRECISION TEST 1999 DOUBLE PRECISION TOP 2000C 2001 CORDIS=RWORK(15) 2002 EPSATE=EIGHT*RWORK(8) 2003 MODNEW=IWORK(4) 2004 MAXCOR=IWORK(17) 2005C 2006C SEE IF MINIMAL NUMBER OF STEPS TAKEN 2007C 2008 IF((IWORK(28).LE.1).OR.(CORDIS.LE.EPSATE)) THEN 2009 QUAL=EIGHT 2010 GO TO 88 2011 ENDIF 2012C 2013C SEE IF AVERAGE NUMBER OF STEPS TAKEN 2014C 2015 IF(MODNEW.EQ.0)THEN 2016 NAVE=(MAXCOR-1)/2 2017 ELSE 2018 NAVE=MAXCOR 2019 ENDIF 2020 IF(IWORK(28).EQ.NAVE)THEN 2021 QUAL=ONE 2022 GO TO 88 2023 ENDIF 2024C 2025C SEE IF MAXIMUM NUMBER OF STEPS TAKEN 2026C 2027 IF(MODNEW.EQ.0)THEN 2028 NMAX=MAXCOR 2029 ELSE 2030 NMAX=2*MAXCOR 2031 ENDIF 2032 IF(IWORK(28).GE.NMAX)THEN 2033 QUAL=THGIE 2034 GO TO 88 2035 ENDIF 2036C 2037C COMPUTE QUAL 2038C 2039C FOR MODNEW=0, LET W=(STEPX/CORDIS), 2040C IEXP=1/(2**(NCOR-1)-1) 2041C U=W**IEXP 2042C JEXP=2**(NCOR-NAVE) 2043C THEN QUAL=(U+1+(1/U)) / (U**JEXP+1+(1/U**JEXP)) 2044C 2045C OTHERWISE, EXP=(NCOR-NAVE)/(NCOR-1) 2046C W=(STEPX/CORDIS) 2047C QUAL=W**EXP 2048C 2049 IF(MODNEW.EQ.0)THEN 2050 ITOP=2**(IWORK(28)-NAVE) 2051 IBOT=2**(IWORK(28)-1)-1 2052 EXPO=ONE/FLOAT(IBOT) 2053 BASE=(STEPX/CORDIS)**EXPO 2054 TOP=BASE+ONE+(ONE/BASE) 2055 TERM=BASE**ITOP 2056 BOT=TERM+ONE+(ONE/TERM) 2057 QUAL=TOP/BOT 2058 IF(QUAL.GT.EIGHT)QUAL=EIGHT 2059 IF(QUAL.LT.THGIE)QUAL=THGIE 2060 GO TO 88 2061 ENDIF 2062C 2063C TRY TO ANTICIPATE UNDERFLOW AND OVERFLOW 2064C 2065 ITOP=IWORK(28)-NAVE 2066 IBOT=IWORK(28)-1 2067 EXPO=FLOAT(IBOT)/FLOAT(ITOP) 2068 TEST=EIGHT**EXPO 2069 BASE=STEPX/CORDIS 2070 ESAB=CORDIS/STEPX 2071 IF((IWORK(28).LT.NAVE.AND.TEST.GT.BASE) 2072 1 .OR.(IWORK(28).GT.NAVE.AND.TEST.LT.BASE))THEN 2073 QUAL=EIGHT 2074 GO TO 88 2075 ENDIF 2076 IF((IWORK(28).LT.NAVE.AND.TEST.GT.ESAB) 2077 1 .OR.(IWORK(28).GT.NAVE.AND.TEST.LT.ESAB))THEN 2078 QUAL=THGIE 2079 GO TO 88 2080 ENDIF 2081C 2082C COMPUTE QUAL 2083C 2084 EXPO=FLOAT(ITOP)/FLOAT(IBOT) 2085 QUAL=BASE**EXPO 2086 IF(QUAL.GT.EIGHT)QUAL=EIGHT 2087 IF(QUAL.LT.THGIE)QUAL=THGIE 208888 CONTINUE 2089 RWORK(23)=QUAL 2090 IF(IWORK(7).GE.3)THEN 2091 LOUNIT=IWORK(8) 2092 WRITE(LOUNIT,1000)QUAL 2093 ENDIF 2094 RETURN 20951000 FORMAT(' COQUAL - Corrector convergence quality factor=',G14.6) 2096 END 2097 SUBROUTINE CORECT(DF,FPAR,FX,IERROR,IHOLD,IPAR,IWORK, 2098 1 NVAR,RWORK,STEPX,WK,XR,LRW,LIW,ICRIT,SLNAME) 2099C 2100C*********************************************************************** 2101C 2102C CORECT PERFORMS THE CORRECTION OF AN APPROXIMATE 2103C SOLUTION OF THE EQUATION F(XR)=0.0. FOR THE CORRECTION EITHER 2104C NEWTON'S METHOD OR THE CHORD NEWTON METHOD IS USED. IN THE LATTER 2105C CASE THE JACOBIAN IS EVALUATED ONLY AT THE STARTING POINT. 2106C IF B IS THE VALUE OF XR(IHOLD) FOR THE INPUT STARTING POINT, 2107C THEN THE AUGMENTING EQUATION IS XR(IHOLD)=B, THAT IS, THE IHOLD-TH 2108C COMPONENT OF XR IS TO BE HELD FIXED. THE AUGMENTED SYSTEM TO BE 2109C SOLVED IS THEN DFA(XR,IHOLD)*(-DELX)=FA(XR) 2110C 2111C ACCEPTANCE AND REJECTION CRITERIA- 2112C 2113C LET NCOR DENOTE THE CURRENT ITERATION INDEX; THAT IS, NCOR=0 FOR 2114C THE PREDICTED POINT WHICH SERVES AS STARTING POINT. 2115C After each iteration, the current point will be accepted 2116C if any of the following conditions hold: 2117C 2118C STRONG ACCEPTANCE CRITERION 2119C 2120C 1. FNRM.LE.ABSERR AND STEPX.LE.(ABSERR+RELERR*XNRM)) 2121C 2122C WEAK ACCEPTANCE CRITERIA 2123C 2124C 2. (NCOR.EQ.0) AND 2125C FNRM.LE.(.5*ABSERR) 2126C 3. FNRM.LE.EPSATE OR STEPX.LE.EPSATE 2127C 4. (NCOR.GE.2) AND 2128C (FNRM+FNRML).LE.ABSERR AND 2129C STEPX.LE.8*(ABSERR+RELERR*XNRM) 2130C 5. (NCOR.GE.2) AND 2131C FNRM.LE.8.0*ABSERR AND 2132C (STEPX+STEPXL).LE.(ABSERR+RELERR*XNRM) 2133C 2134C AFTER EACH STEP OF THE ITERATION, THE PROCESS IS TO BE ABORTED IF 2135C ANY OF THE FOLLOWING CONDITIONS HOLDS 2136C 2137C 1. FNRM.GT.(FMP*FNRML+ABSERR) 2138C 2. (NCOR.GE.2) AND (ICRIT.EQ.0) AND 2139C (STEPX.GT.(FMP*STEPXL+ABSERR)) 2140C 3. NCOR.GE.MAXCOR (NUMBER OF ITERATIONS EXCEEDED). 2141C 2142C HERE FMP=2.0 FOR NCOR.EQ.1, FMP=1.05 FOR NCOR.GT.1 2143C 2144C ERROR CONDITIONS 2145C IERROR=1 DATA OR STORAGE ERROR 2146C IERROR=2 ERROR IN FUNCTION OR DERIVATIVE CALL 2147C IERROR=3 SOLVER FAILED 2148C IERROR=4 UNSUCCESSFUL ITERATION 2149C IERROR=5 TOO MANY CORRECTOR STEPS 2150C 2151 DOUBLE PRECISION EIGHT 2152 DOUBLE PRECISION HALF 2153 DOUBLE PRECISION ONE 2154 DOUBLE PRECISION ONEFIV 2155 DOUBLE PRECISION TWO 2156 DOUBLE PRECISION ZERO 2157C 2158 PARAMETER (EIGHT=8.0) 2159 PARAMETER (HALF=0.5) 2160 PARAMETER (ONE=1.0) 2161 PARAMETER (ONEFIV=1.05) 2162 PARAMETER (TWO=2.0) 2163 PARAMETER (ZERO=0.0) 2164C 2165 EXTERNAL DF 2166 EXTERNAL FX 2167 EXTERNAL IDAMAX 2168 EXTERNAL DAXPY 2169 EXTERNAL SLNAME 2170 EXTERNAL DNRM2 2171C 2172 INTRINSIC ABS 2173C 2174 INTEGER LIW 2175 INTEGER LRW 2176 INTEGER NVAR 2177C 2178 DOUBLE PRECISION ABSERR 2179 DOUBLE PRECISION DETS 2180 DOUBLE PRECISION EPSATE 2181 DOUBLE PRECISION FMP 2182 DOUBLE PRECISION FNRM 2183 DOUBLE PRECISION FNRML 2184 DOUBLE PRECISION FPAR(*) 2185 INTEGER I 2186 INTEGER ICRIT 2187 INTEGER IERROR 2188 INTEGER IFMAX 2189 INTEGER IHOLD 2190 INTEGER IPAR(*) 2191 INTEGER IDAMAX 2192 INTEGER IWORK(LIW) 2193 INTEGER IWRITE 2194 INTEGER IXMAX 2195 INTEGER JOB 2196 INTEGER KSMAX 2197 INTEGER LOUNIT 2198 INTEGER MAXCOR 2199 INTEGER MAXNEW 2200 INTEGER MODNEW 2201 DOUBLE PRECISION RELERR 2202 DOUBLE PRECISION RWORK(LRW) 2203 DOUBLE PRECISION SKALE 2204 DOUBLE PRECISION DNRM2 2205 DOUBLE PRECISION STEPX 2206 DOUBLE PRECISION STEPXL 2207 DOUBLE PRECISION TLSTEP 2208 DOUBLE PRECISION WK(NVAR) 2209 DOUBLE PRECISION XNRM 2210 DOUBLE PRECISION XR(NVAR) 2211 DOUBLE PRECISION XVALUE 2212C 2213C Initialize. 2214C 2215 ABSERR=RWORK(1) 2216 RELERR=RWORK(2) 2217 EPSATE=EIGHT*RWORK(8) 2218 MODNEW=IWORK(4) 2219 IWRITE=IWORK(7) 2220 LOUNIT=IWORK(8) 2221 MAXCOR=IWORK(17) 2222 IERROR=0 2223 IWORK(28)=0 2224 MAXNEW=MAXCOR 2225 IF(MODNEW.NE.0)MAXNEW=2*MAXCOR 2226 FMP=TWO 2227 STEPX=ZERO 2228 XVALUE=XR(IHOLD) 2229C 2230C STORE INITIAL FUNCTION VALUE IN THE WORK VECTOR WK 2231C 2232 CALL FX(NVAR,FPAR,IPAR,XR,WK,IERROR) 2233 IWORK(22)=IWORK(22)+1 2234 IF(IERROR.NE.0)THEN 2235 WRITE(LOUNIT,*) 2236 * 'CORECT - Error flag from user function routine.' 2237 RETURN 2238 ENDIF 2239 WK(NVAR)=XR(IHOLD)-XVALUE 2240 IFMAX=IDAMAX(NVAR,WK,1) 2241 FNRM=DNRM2(NVAR-1,WK,1) 2242 IXMAX=IDAMAX(NVAR,XR,1) 2243 XNRM=DNRM2(NVAR,XR,1) 2244 IF(IWRITE.GE.2)WRITE(LOUNIT,70)IWORK(28),FNRM,IFMAX 2245 IF(FNRM.LE.HALF*ABSERR)RETURN 2246C 2247C NEWTON ITERATION LOOP BEGINS 2248C 2249 DO 100 I=1,MAXNEW 2250 IWORK(28)=I 2251C 2252C SOLVE SYSTEM FPRIME*(-DELX)=FX FOR NEWTON DIRECTION DELX 2253C 2254 JOB=0 2255 IF(I.NE.1.AND.I.NE.MAXCOR.AND.MODNEW.EQ.1)JOB=1 2256 IF(MODNEW.EQ.2)JOB=1 2257 CALL SLNAME(DETS,FX,DF,FPAR,IERROR,IHOLD,IPAR,IWORK,LIW, 2258 1 JOB,NVAR,RWORK,LRW,XR,WK) 2259 IF(IERROR.NE.0)THEN 2260 WRITE(LOUNIT,12)IERROR 2261 RETURN 2262 ENDIF 2263C 2264C UPDATE X, COMPUTE SIZE OF STEP AND OF X 2265C 2266 SKALE=-ONE 2267 CALL DAXPY(NVAR,SKALE,WK,1,XR,1) 2268 STEPXL=STEPX 2269 KSMAX=IDAMAX(NVAR,WK,1) 2270 STEPX=ABS(WK(KSMAX)) 2271 IXMAX=IDAMAX(NVAR,XR,1) 2272 XNRM=DNRM2(NVAR,XR,1) 2273C 2274C Compute F(X) and take its norm. 2275C 2276 CALL FX(NVAR,FPAR,IPAR,XR,WK,IERROR) 2277 IWORK(22)=IWORK(22)+1 2278 IF(IERROR.NE.0) THEN 2279 WRITE(LOUNIT,*) 2280 * 'CORECT - Error flag from user function routine.' 2281 RETURN 2282 ENDIF 2283 WK(NVAR)=XR(IHOLD)-XVALUE 2284 FNRML=FNRM 2285 IFMAX=IDAMAX(NVAR,WK,1) 2286 FNRM=DNRM2(NVAR-1,WK,1) 2287 IF(IWRITE.GE.2) THEN 2288 WRITE(LOUNIT,80)IWORK(28),STEPX,KSMAX 2289 WRITE(LOUNIT,70)IWORK(28),FNRM,IFMAX 2290 ENDIF 2291C 2292C Check for strong acceptance of function and stepsize. 2293C 2294 TLSTEP=ABSERR+RELERR*XNRM 2295 IF(FNRM.LE.ABSERR.AND.STEPX.LE.TLSTEP)RETURN 2296C 2297C Check for weak acceptance of function and stepsize. 2298C 2299 IF(FNRM.LE.EPSATE.OR.STEPX.LE.EPSATE)RETURN 2300 IF(IWORK(28).GT.1)THEN 2301 IF((FNRM+FNRML).LE.ABSERR.AND.STEPX.LE.EIGHT*TLSTEP)RETURN 2302 IF(FNRM.LE.EIGHT*ABSERR.AND.(STEPX+STEPXL).LE.TLSTEP)RETURN 2303 ENDIF 2304C 2305C Decide if iteration should be aborted 2306C 2307 IF(IWORK(28).GT.1)THEN 2308 IF(ICRIT.LT.1.AND.STEPX.GT.(FMP*STEPXL+ABSERR)) THEN 2309 IERROR=4 2310 IF(IWRITE.GE.2)WRITE(LOUNIT,*) 2311 * 'CORECT - Size of correction not decreasing.' 2312 RETURN 2313 ENDIF 2314 ENDIF 2315 IF(ICRIT.LT.2.AND.FNRM.GT.(FMP*FNRML+ABSERR)) THEN 2316 IERROR=4 2317 IF(IWRITE.GE.2)WRITE(LOUNIT,*) 2318 * 'CORECT - Residual is not decreasing.' 2319 RETURN 2320 ENDIF 2321 FMP=ONEFIV 2322 100 CONTINUE 2323C 2324C Reached maximum number of steps without acceptance or rejection. 2325C 2326 IERROR=5 2327 IF(IWRITE.GE.2)WRITE(LOUNIT,*)'CORECT - Convergence too slow.' 2328 RETURN 232970 FORMAT(' CORECT - Residual ',I6,'=',G14.6,' I=',I6) 233080 FORMAT(' CORECT - Step ',I6,15X,G14.6,' I=',I6) 233112 FORMAT(' CORECT - Error flag=',I6,' from solver.') 2332 END 2333 SUBROUTINE ROOT(A,FA,B,FB,U,FU,KOUNT,IFLAG,IERROR,EPMACH) 2334C 2335C*********************************************************************** 2336C 2337C ROOT SEEKS A ROOT OF THE EQUATION F(X)=0.0, 2338C GIVEN A STARTING INTERVAL (A,B) ON WHICH F CHANGES SIGN. 2339C ON FIRST CALL TO ROOT, THE INTERVAL AND FUNCTION VALUES FA AND FB 2340C ARE INPUT AND AN APPROXIMATION U FOR THE ROOT IS RETURNED. 2341C BEFORE EACH SUBSEQUENT CALL, THE USER EVALUATES FU=F(U), AND THE 2342C PROGRAM TRIES TO RETURN A BETTER APPROXIMATION U. 2343C 2344C SEE THE BOOK BY BRENT LISTED in the documentation. 2345C 2346C VARIABLES 2347C 2348C A = IS ONE ENDPOINT OF AN INTERVAL IN WHICH F CHANGES SIGN. 2349C FA = THE VALUE OF F(A). THE USER MUST EVALUATE F(A) BEFORE 2350C FIRST CALL ONLY. THEREAFTER THE PROGRAM SETS FA. 2351C B = IS THE OTHER ENDPOINT OF THE INTERVAL IN WHICH 2352C F CHANGES SIGN. NOTE THAT THE PROGRAM WILL RETURN 2353C IMMEDIATELY WITH AN ERROR FLAG IF FB*FA.GT.0.0. 2354C FB = THE VALUE OF F(B). THE USER MUST EVALUATE F(B) BEFORE 2355C FIRST CALL ONLY. THERAFTER THE PROGRAM SETS FB. 2356C U = ON FIRST CALL, U SHOULD NOT BE SET BY THE USER. 2357C ON SUBSEQUENT CALLS, U SHOULD NOT BE CHANGED 2358C FROM ITS OUTPUT VALUE, THE CURRENT APPROXIMANT 2359C TO THE ROOT. 2360C FU = ON FIRST CALL, FU SHOULD NOT BE SET BY THE USER. 2361C THEREAFTER, THE USER SHOULD EVALUATE THE FUNCTION 2362C AT THE OUTPUT VALUE U, AND RETURN FU=F(U). 2363C KOUNT = A COUNTER FOR THE NUMBER OF CALLS TO ROOT. KOUNT 2364C SHOULD BE SET TO ZERO ON THE FIRST CALL FOR A GIVEN 2365C ROOT PROBLEM. 2366C IFLAG = PROGRAM RETURN FLAG 2367C IFLAG=-1 THE CURRENT BRACKETING INTERVAL WITH 2368C ENDPOINTS STORED IN A AND B IS LESS THAN 2369C 4*EPMACH*ABS(U)+EPMACH 2370C HENCE U SHOULD BE ACCEPTED AS THE ROOT. 2371C THE FUNCTION VALUE F(U) IS STORED IN FU. 2372C IFLAG= 0 THE INPUT VALUE FU IS EXACTLY ZERO, AND 2373C U SHOULD BE ACCEPTED AS THE ROOT. 2374C IFLAG>0 THE CURRENT APPROXIMATION TO THE ROOT IS 2375C CONTAINED IN U. IF A BETTER APPROXIMATION IS 2376C DESIRED, SET FU=F(U) AND CALL ROOT AGAIN. 2377C THE VALUE OF IFLAG INDICATES 2378C THE METHOD THAT WAS USED TO PRODUCE U. 2379C 2380C IFLAG= 1 BISECTION WAS USED. 2381C IFLAG= 2 LINEAR INTERPOLATION (SECANT METHOD). 2382C IFLAG= 3 INVERSE QUADRATIC INTERPOLATION. 2383C 2384C IERROR= GLOBAL ERROR FLAG. 2385C IERROR=0 NO ERROR FOUND 2386C IERROR=7 ON FIRST CALL, FA*FB.GT.0.0. AND HENCE 2387C THE GIVEN INTERVAL IS UNACCEPTABLE. 2388C EPMACH= THE RELATIVE MACHINE PRECISION 2389C 2390 DOUBLE PRECISION EIGHT 2391 DOUBLE PRECISION HALF 2392 DOUBLE PRECISION ONE 2393 DOUBLE PRECISION ONEP5 2394 DOUBLE PRECISION TWO 2395 DOUBLE PRECISION ZERO 2396C 2397 PARAMETER (EIGHT=8.0) 2398 PARAMETER (HALF=0.5) 2399 PARAMETER (ONE=1.0) 2400 PARAMETER (ONEP5=1.5) 2401 PARAMETER (TWO=2.0) 2402 PARAMETER (ZERO=0.0) 2403C 2404 INTRINSIC ABS 2405 INTRINSIC SIGN 2406C 2407 DOUBLE PRECISION A 2408 DOUBLE PRECISION B 2409 DOUBLE PRECISION EPMACH 2410 DOUBLE PRECISION FA 2411 DOUBLE PRECISION FB 2412 DOUBLE PRECISION FU 2413 DOUBLE PRECISION HALFUB 2414 INTEGER IERROR 2415 INTEGER IFLAG 2416 INTEGER KOUNT 2417 DOUBLE PRECISION P 2418 DOUBLE PRECISION Q 2419 DOUBLE PRECISION R 2420 DOUBLE PRECISION S 2421 DOUBLE PRECISION SDEL1 2422 DOUBLE PRECISION SDEL2 2423 DOUBLE PRECISION SDEL3 2424 DOUBLE PRECISION SDEL4 2425 DOUBLE PRECISION STEP 2426 DOUBLE PRECISION TOLER 2427 DOUBLE PRECISION U 2428C 2429C SEGMENT 1 FIRST CALL HANDLED SPECIALLY. DO BOOKKEEPING. 2430C 2431 IF(KOUNT.LE.0)THEN 2432 IF((FA.GT.ZERO.AND.FB.GT.ZERO) 2433 1 .OR.(FA.LT.ZERO.AND.FB.LT.ZERO))THEN 2434 IERROR=7 2435 KOUNT=0 2436 RETURN 2437 ENDIF 2438 KOUNT=1 2439 SDEL1=TWO*ABS(B-A) 2440 SDEL2=TWO*SDEL1 2441 SDEL3=TWO*SDEL2 2442 U=B 2443 B=A 2444 FU=FB 2445 FB=FA 2446 ELSE 2447C 2448C INCREMENT COUNTER AND CHECK WHETHER F(U)=0. 2449C 2450 KOUNT=KOUNT+1 2451 IF(FU.EQ.ZERO)THEN 2452 IFLAG=0 2453 RETURN 2454 ENDIF 2455C 2456C IF FU AND FB HAVE THE SAME SIGN OVERWRITE B WITH A 2457C 2458 IF(SIGN(ONE,FU).EQ.SIGN(ONE,FB))THEN 2459 B=A 2460 FB=FA 2461 ENDIF 2462 ENDIF 2463C 2464C SEGMENT 2 REARRANGE POINTS AND FUNCTION VALUES IF 2465C NECESSARY SO THAT ABS(FU).LT.ABS(FB) 2466C 2467 IF(ABS(FB).LT.ABS(FU)) THEN 2468 A=U 2469 U=B 2470 B=A 2471 FA=FU 2472 FU=FB 2473 FB=FA 2474 ENDIF 2475C 2476C SEGMENT 3 CHECK FOR ACCEPTANCE BECAUSE OF SMALL INTERVAL 2477C CURRENT CHANGE IN SIGN INTERVAL IS (B,U) OR (U,B). 2478C 2479 TOLER=TWO*EPMACH*ABS(U)+EPMACH 2480 HALFUB=HALF*(B-U) 2481 SDEL4=SDEL3 2482 SDEL3=SDEL2 2483 SDEL2=SDEL1 2484 SDEL1=ABS(B-U) 2485 IF(ABS(HALFUB).LE.TOLER) THEN 2486 IFLAG=-1 2487 A=U 2488 FA=FU 2489 RETURN 2490 ENDIF 2491C 2492C SEGMENT 4 COMPUTE NEW APPROXIMANT TO ROOT OF THE FORM 2493C U(NEW)=U(OLD)+STEP. 2494C METHODS AVAILABLE ARE LINEAR INTERPOLATION 2495C INVERSE QUADRATIC INTERPOLATION 2496C AND BISECTION. 2497C 2498 IF(ABS(FU).GE.ABS(FA)) THEN 2499 IFLAG=1 2500 STEP=HALFUB 2501 GO TO 80 2502 ENDIF 2503C 2504C ATTEMPT LINEAR INTERPOLATION IF ONLY TWO POINTS AVAILABLE 2505C COMPUTE P AND Q FOR APPROXIMATION U(NEW)=U(OLD)+P/Q 2506C 2507 IF(A.EQ.B) THEN 2508 IFLAG=2 2509 S=FU/FA 2510 P=TWO*HALFUB*S 2511 Q=ONE-S 2512C 2513C ATTEMPT INVERSE QUADRATIC INTERPOLATION IF THREE POINTS AVAILABLE 2514C COMPUTE P AND Q FOR APPROXIMATION U(NEW)=U(OLD)+P/Q 2515C 2516 ELSE 2517 IFLAG=3 2518 S=FU/FA 2519 Q=FA/FB 2520 R=FU/FB 2521 P=S*(TWO*HALFUB*Q*(Q-R)-(U-A)*(R-ONE)) 2522 Q=(Q-ONE)*(R-ONE)*(S-ONE) 2523 ENDIF 2524C 2525C CORRECT THE SIGNS OF P AND Q 2526C 2527 IF(P.GT.ZERO)Q=-Q 2528 P=ABS(P) 2529C 2530C IF P/Q IS TOO LARGE, GO BACK TO BISECTION 2531C 2532 IF((EIGHT*SDEL1.GT.SDEL4) 2533 1 .OR.(P.GE.ONEP5*ABS(HALFUB*Q)-ABS(TOLER*Q))) THEN 2534 IFLAG=1 2535 STEP=HALFUB 2536 GO TO 80 2537 ENDIF 2538 STEP=P/Q 2539C 2540C SEGMENT 5 VALUE OF STEP HAS BEEN COMPUTED. 2541C UPDATE INFORMATION A =U, FA =FU, U =U+STEP. 2542C CHANGE IN SIGN INTERVAL IS NOW (A,B) OR (B,A). 2543C 254480 CONTINUE 2545 A=U 2546 FA=FU 2547 IF(ABS(STEP).LE.TOLER) STEP=SIGN(TOLER,HALFUB) 2548 U=U+STEP 2549 RETURN 2550 END 2551 SUBROUTINE SETSTP(HTAN,IWORK,LIW,LRW,RWORK) 2552C 2553C*********************************************************************** 2554C 2555C Predictor stepsize computation. 2556C 2557C THE FORMULAS UNDERLYING THE ALGORITHM ARE 2558C 2559C ALPHLC=MAXIMUM OF ARCCOS(TL,TC) AND ALFMIN=2*ARCCOS(1-EPMACH) 2560C HSEC =NORM(XC-XF) 2561C HSECL =NORM(XL-XC) 2562C ABSNLC=ABS(SIN(.5*ALPHLC)) 2563C CURV =LAST VALUE OF THE CURVATURE 2564C CURVN =2*ABSNLC/HSEC 2565C CORDIS=CORRECTOR DISTANCE TO CURRENT CONTINUATION POINT. 2566C FORCE CORDIS TO LIE BETWEEN .01*HSEC AND HSEC. 2567C UNLESS (CORDIS.EQ.0.0), BECAUSE THE PREDICTED POINT WAS 2568C ACCEPTED. IN SUCH A CASE, SET HTAN=HFACT*HSEC 2569C INSTEAD OF USING FIRST ESTIMATE FOR HTAN. 2570C 2571C THEN 2572C 2573C CURVX=CURVN+HSEC*(CURVN-CURV)/(HSEC+HSECL) 2574C A SIMPLER FORMULA IS USED IF WE DO NOT HAVE DATA AT TWO PREVIOUS 2575C POINTS. 2576C 2577C IF (IWORK(10).GE.2) CURVX MUST BE AT LEAST AS LARGE AS 2578C THE MAXIMUM OF 0.001 AND .01/HSEC 2579C 2580C FIRST ESTIMATE (UNLESS (CORDIS.EQ.0.0) ) 2581C 2582C HTAN =SQRT(2*QUAL*CORDIS/CURVX) 2583C 2584C ADJUSTED VALUE 2585C 2586C HTAN=HTAN*(1.0+HTAN*(TC(IPC)-TL(IPC))/(2*HSEC*TC(IPC))) 2587C 2588C READJUSTMENT AND TRUNCATIONS 2589C 2590C IF STEPSIZE REDUCTION OCCURRED DURING LAST CORRECTOR PROCESS, 2591C HTAN IS FORCED TO BE LESS THAN (HFACT-1)*HSEC/2. 2592C 2593C HTAN MUST LIE BETWEEN (HSEC/HFACT) AND (HSEC*HFACT). 2594C HTAN IS ALWAYS FORCED TO LIE BETWEEN USER SPECIFIED 2595C MINIMUM AND MAXIMUM STEPS. 2596C 2597 DOUBLE PRECISION HALF 2598 DOUBLE PRECISION HUNDTH 2599 DOUBLE PRECISION ONE 2600 DOUBLE PRECISION THOUTH 2601 DOUBLE PRECISION TWO 2602 DOUBLE PRECISION ZERO 2603C 2604 PARAMETER (HALF=0.5) 2605 PARAMETER (HUNDTH=0.01) 2606 PARAMETER (ONE=1.0) 2607 PARAMETER (THOUTH=0.001) 2608 PARAMETER (TWO=2.0) 2609 PARAMETER (ZERO=0.0) 2610C 2611 INTRINSIC ABS 2612 INTRINSIC MAX 2613 INTRINSIC MIN 2614 INTRINSIC SIN 2615 INTRINSIC SQRT 2616C 2617 INTEGER LIW 2618 INTEGER LRW 2619C 2620 DOUBLE PRECISION CURV 2621 DOUBLE PRECISION CURVX 2622 DOUBLE PRECISION HTAN 2623 INTEGER IWORK(LIW) 2624 DOUBLE PRECISION RWORK(LRW) 2625 DOUBLE PRECISION TEMP 2626 DOUBLE PRECISION TMP 2627C 2628 RWORK(11)=MAX(RWORK(10),RWORK(11)) 2629C 2630C COMPUTE NEW CURVATURE DATA 2631C 2632 CURV=RWORK(16) 2633 RWORK(16)=TWO*ABS(SIN(HALF*RWORK(11)))/RWORK(21) 2634 IF(CURV.EQ.ZERO)CURV=RWORK(16) 2635 CURVX=RWORK(16) 2636 IF(RWORK(22).NE.ZERO)CURVX= 2637 1 RWORK(16)+RWORK(21)*(RWORK(16)-CURV)/(RWORK(21)+RWORK(22)) 2638 TMP=MAX(THOUTH,(HUNDTH)/RWORK(21)) 2639 CURVX=MAX(CURVX,TMP) 2640C 2641C IF CONVERGENCE DISTANCE IS ZERO, SET ESTIMATE TO HFACT*HSEC. 2642C ELSE TRUNCATE QUAL*CORDIS TO LIE BETWEEN .01*HSEC AND HSEC. 2643C 2644 HTAN=RWORK(20)*RWORK(21) 2645 IF(RWORK(15).NE.ZERO)THEN 2646 TEMP=MAX(RWORK(23)*RWORK(15),HUNDTH*RWORK(21)) 2647 TEMP=MIN(TEMP,RWORK(21)) 2648 HTAN=SQRT(TWO*TEMP/CURVX) 2649 ENDIF 2650C 2651C Adjust the step to account for estimated curvature in the direction 2652C of the parameter. 2653C 2654 IF(IWORK(18).GT.0) 2655 * HTAN=MIN(HTAN,HALF*(RWORK(20)-ONE)*RWORK(21)) 2656 IF(IWORK(3).NE.1)THEN 2657 TEMP=ONE+(ONE-RWORK(25)/RWORK(24))*(HALF*HTAN)/RWORK(21) 2658 HTAN=HTAN*TEMP 2659 ENDIF 2660C 2661C Enforce restrictions on growth and size of the step. 2662C 2663 HTAN=MAX(HTAN,RWORK(21)/RWORK(20)) 2664 HTAN=MIN(HTAN,RWORK(21)*RWORK(20)) 2665 HTAN=MAX(HTAN,RWORK(3)) 2666 HTAN=MIN(HTAN,RWORK(4)) 2667 RETURN 2668 END 2669 SUBROUTINE START(DF,FPAR,FX,IERROR,IPAR,IPC,IWORK,IWRITE,LIW, 2670 *LOUNIT,LRW,NVAR,RWORK,TC,WORK1,XC,XF,XR,SLNAME) 2671C 2672C*********************************************************************** 2673C 2674C On first call, make sure that the input point satisfies the 2675C nonlinear equations. If not, call CORECT and attempt to force this. 2676C 2677 DOUBLE PRECISION ONE 2678 DOUBLE PRECISION ZERO 2679C 2680 PARAMETER (ONE=1.0) 2681 PARAMETER (ZERO=0.0) 2682C 2683 EXTERNAL COQUAL 2684 EXTERNAL CORECT 2685 EXTERNAL DF 2686 EXTERNAL FX 2687 EXTERNAL IDAMAX 2688 EXTERNAL DAXPY 2689 EXTERNAL DCOPY 2690 EXTERNAL SLNAME 2691C 2692 INTRINSIC ABS 2693C 2694 INTEGER LIW 2695 INTEGER LRW 2696 INTEGER NVAR 2697C 2698 DOUBLE PRECISION DETS 2699 DOUBLE PRECISION FPAR(*) 2700 INTEGER I 2701 INTEGER ICRIT 2702 INTEGER IERROR 2703 INTEGER IMAX 2704 INTEGER IPAR(*) 2705 INTEGER IPC 2706 INTEGER IDAMAX 2707 INTEGER IWORK(LIW) 2708 INTEGER IWRITE 2709 INTEGER JOB 2710 INTEGER LOUNIT 2711 INTEGER MODSAV 2712 DOUBLE PRECISION RWORK(LRW) 2713 DOUBLE PRECISION SKALE 2714 DOUBLE PRECISION STEPX 2715 DOUBLE PRECISION TC(NVAR) 2716 DOUBLE PRECISION WORK1(NVAR) 2717 DOUBLE PRECISION XC(NVAR) 2718 DOUBLE PRECISION XF(NVAR) 2719 DOUBLE PRECISION XR(NVAR) 2720C 2721C If user is requesting that Jacobian be used as long as possible, 2722C then go ahead, generate and factor the first one now. 2723C 2724 IF(IWORK(4).EQ.2)THEN 2725 JOB=2 2726 CALL SLNAME(DETS,FX,DF,FPAR,IERROR,IPC,IPAR,IWORK,LIW, 2727 1 JOB,NVAR,RWORK,LRW,XR,WORK1) 2728 RWORK(17)=DETS 2729 IF(IERROR.NE.0)THEN 2730 WRITE(LOUNIT,*) 2731 * 'START - Could not factor initial jacobian.' 2732 RETURN 2733 ENDIF 2734 ENDIF 2735 IF(IWRITE.GE.2)WRITE(LOUNIT,1040)IPC 27361040 FORMAT(' START - Correct initial point, fixing index ',I5) 2737C 2738C Set up a pseudo-tangent vector and other data for first time, 2739C check that starting point satisfies the equations. 2740C 2741 DO 50 I=1,NVAR 2742 TC(I)=ZERO 274350 CONTINUE 2744 CALL DCOPY(NVAR,XR,1,XC,1) 2745 TC(IPC)=ONE 2746 MODSAV=IWORK(4) 2747 ICRIT=1 2748C 274960 CONTINUE 2750 CALL DCOPY(NVAR,XC,1,XR,1) 2751 CALL CORECT(DF,FPAR,FX,IERROR,IPC,IPAR,IWORK,NVAR,RWORK,STEPX, 2752 *WORK1,XR,LRW,LIW,ICRIT,SLNAME) 2753 IWORK(25)=IWORK(25)+IWORK(28) 2754C 2755C If possible, retry with ICRIT=2. 2756C 2757 IF(IERROR.NE.0.AND.ICRIT.EQ.1)THEN 2758 IF(IWRITE.GE.1)WRITE(LOUNIT,*) 2759 * 'START - Retry starting point correction' 2760 ICRIT=2 2761 GO TO 60 2762 ELSEIF(IERROR.NE.0)THEN 2763 ICRIT=1 2764 ENDIF 2765C 2766 IF(IERROR.NE.0.AND.IWORK(4).GT.0)THEN 2767 IWORK(4)=IWORK(4)-1 2768 IERROR=0 27691110 FORMAT(' START - Retrying starting point with IWORK(4)=',I6) 2770 IF(IWRITE.GE.1)WRITE(LOUNIT,1110)IWORK(4) 2771 GO TO 60 2772 ENDIF 2773 IWORK(4)=MODSAV 2774 IF(IERROR.NE.0)THEN 2775 WRITE(LOUNIT,*) 2776 * 'START - Starting point correction failed.' 2777 RETURN 2778 ENDIF 2779 SKALE=-ONE 2780 CALL DAXPY(NVAR,SKALE,XR,1,XC,1) 2781 IMAX=IDAMAX(NVAR,XC,1) 2782 RWORK(15)=ABS(XC(IMAX)) 2783 CALL DCOPY(NVAR,XR,1,XC,1) 2784 CALL DCOPY(NVAR,XR,1,XF,1) 2785 CALL COQUAL(STEPX,IWORK,LIW,RWORK,LRW) 2786 RWORK(14)=RWORK(13) 2787 IWORK(27)=IWORK(27)+1 2788 IWORK(10)=1 2789 IWORK(1)=1 2790 RETURN 2791 END 2792 SUBROUTINE TANGNT(DETSN,FX,DF,FPAR,IERROR,IP,IPAR,IWORK, 2793 1 NVAR,RWORK,TAN,XR,LIW,LRW,SLNAME) 2794C 2795C*********************************************************************** 2796C 2797C TANGNT COMPUTES A UNIT TANGENT VECTOR TO THE SOLUTION 2798C CURVE OF THE UNDERDETERMINED NONLINEAR SYSTEM FX = 0. THE 2799C TANGENT VECTOR TAN IS THE SOLUTION OF THE LINEAR SYSTEM 2800C 2801C DFA(X,IPL)*TAN = E(NVAR) 2802C 2803C HERE, E(I) ALWAYS DENOTES THE I-TH BASIS NATURAL BASIS VECTOR 2804C IN NVAR-SPACE; THAT IS, THE VECTOR WITH A 1 IN THE I-TH 2805C COMPONENT AND ZEROS ELSEWHERE. THEN DFA(X,IPL) IS THE 2806C AUGMENTED JACOBIAN WHOSE FIRST NVAR-1 ROWS ARE DFX/DX (X), THE 2807C DERIVATIVE OF FX EVELUATED AT X, AND WHOSE LAST ROW IS 2808C THE TRANSPOSE OF E(IPL). 2809C 2810C THE TANGENT VECTOR IS THEN NORMALIZED, BUT THE ADJUSTMENT 2811C OF ITS SIGN IS PERFORMED OUTSIDE THE ROUTINE. 2812C 2813C ERROR RETURNS IERROR=1 DATA OR STORAGE ERROR 2814C IERROR=2 ERROR IN DERIVATIVE CALL 2815C IERROR=3 SOLVER FAILED 2816C IERROR=6 TANGENT VECTOR ZERO 2817C 2818 DOUBLE PRECISION ONE 2819 DOUBLE PRECISION ZERO 2820C 2821 PARAMETER (ONE=1.0) 2822 PARAMETER (ZERO=0.0) 2823C 2824 EXTERNAL FX 2825 EXTERNAL DF 2826 EXTERNAL SLNAME 2827 EXTERNAL DNRM2 2828 EXTERNAL DSCAL 2829C 2830 INTEGER LIW 2831 INTEGER LRW 2832 INTEGER NVAR 2833C 2834 DOUBLE PRECISION DETSN 2835 DOUBLE PRECISION FPAR(*) 2836 INTEGER I 2837 INTEGER IERROR 2838 INTEGER IP 2839 INTEGER IPAR(*) 2840 INTEGER IWORK(LIW) 2841 INTEGER JOB 2842 DOUBLE PRECISION RWORK(LRW) 2843 DOUBLE PRECISION SKALE 2844 DOUBLE PRECISION DNRM2 2845 DOUBLE PRECISION TAN(NVAR) 2846 DOUBLE PRECISION TNORM 2847 DOUBLE PRECISION XR(NVAR) 2848C 2849C SET RIGHT HAND SIDE OF TANGENT SYSTEM 2850C 2851 DO 10 I=1,NVAR 2852 TAN(I)=ZERO 2853 10 CONTINUE 2854 TAN(NVAR)=ONE 2855C 2856C CALL SOLVER 2857C 2858 JOB=0 2859 IF(IWORK(4).EQ.2)JOB=1 2860 CALL SLNAME(DETSN,FX,DF,FPAR,IERROR,IP,IPAR,IWORK,LIW, 2861 1JOB,NVAR,RWORK,LRW,XR,TAN) 2862 IF(IERROR.NE.0)RETURN 2863C 2864C Normalize the tangent vector. 2865C 2866 TNORM=DNRM2(NVAR,TAN,1) 2867 IF(TNORM.EQ.ZERO)THEN 2868 IERROR=6 2869 ELSE 2870 SKALE=ONE/TNORM 2871 CALL DSCAL(NVAR,SKALE,TAN,1) 2872 IERROR=0 2873 ENDIF 2874 RETURN 2875 END 2876 SUBROUTINE DENSLV(DETS,FX,DF,FPAR,IERROR,IPC,IPAR, 2877 1 IWORK,LIW,JOB,NVAR,RWORK,LRW,X,Y) 2878C 2879C*********************************************************************** 2880C 2881C DENSLV SOLVES THE NVAR*NVAR LINEAR SYSTEM 2882C 2883C ( DF(X) ) 2884C (1) ( ) * Y = Y 2885C ( E(IPC) ) 2886C 2887C WHERE DF(X) IS THE JACOBIAN OF THE GIVEN FUNCTION FX AT X, AND 2888C E(IPC) IS THE IP-TH NATURAL BASIS VECTOR; THAT IS, THE LAST ROW 2889C CONSISTS OF ALL ZEROS EXCEPT FOR A 1 IN COLUMN IPC. 2890C 2891C THE MATRIX IS STORED IN DENSE FORM IN THE ONE-DIMENSIONAL ARRAY 2892C RWORK(LRBEG)-RWORK(LRLST) WHERE LRLST=LRBEG+NVAR*NVAR - 1 AND 2893C LRBEG IS STORED IN IWORK(15). A DATA ERROR IS CAUSED IF LRLST.GT.LRW 2894C 2895C THE ROUTINE STORES THE PIVOT ARRAY IN IWORK(IWORK(13))) THROUGH 2896C IWORK(LILST) WHERE LILST=IWORK(13)+NVAR - 1. 2897C THE ROUTINE RETURNS WITH A DATA ERROR IF LILST.GT.LIW. 2898C 2899C SINCE THE MATRIX OF THE SYSTEM (1) IS STORED IN DENSE FORM THE USER 2900C PROGRAM DF MAY HAVE THE FORM 2901C 2902C SUBROUTINE DF(NVAR,FPAR,IPAR,X,A,IERR) 2903C DIMENSION A(NVAR,NVAR), X(NVAR), FPAR(*),IPAR(*) 2904C ..... 2905C A(I,J)=DF(I)/DX(J) , I=1,...,NVAR-1,J=1,...,NVAR 2906C ..... 2907C 2908C WHERE DF(I)/DX(J) DENOTES THE DERIVATIVE OF THE I-TH COMPONENT 2909C OF THE FUNCTION FX WITH RESPECT TO THE J-TH VARIABLE. THE USER 2910C NEED NOT SUPPLY THE LAST ROW OF THE MATRIX OF THE SYSTEM (1). 2911C 2912C THE LINPACK SUBROUTINES DGEFA, DGEDI, AND DGESL ARE USED. 2913C 2914C LIST OF VARIABLES 2915C 2916C DETS THE SIGN OF THE DETERMINANT OF THE MATRIX IN (1). 2917C FX AN EXTERNAL ROUTINE, THE NAME OF THE USER SUPPLIED 2918C SUBROUTINE WHICH COMPUTES THE NVAR-1 DIMENSIONAL 2919C FUNCTION. 2920C 2921C DF AN EXTERNAL ROUTINE, THE NAME OF THE USER SUPPLIED 2922C SUBROUTINE WHICH COMPUTES THE JACOBIAN OF THE FUNCTION. 2923C FPAR A REAL ARRAY FOR THE TRANSMISSION OF PARAMETERS TO 2924C DENSLV AND DF. 2925C IERROR ERROR FLAG 2926C = 0 NORMAL RETURN 2927C = 1 DATA OR STORAGE ERROR 2928C = 2 ERROR IN DERIVATIVE ROUTINE DF 2929C = 3 NUMERICALLY SINGULAR MATRIX DETECTED. 2930C IPC THE LOCATION OF THE 1 IN THE LAST ROW OF THE MATRIX 2931C IPAR AN INTEGER ARRAY FOR TRANSMISSION OF PARAMETERS TO DENSLV 2932C AND DF. 2933C IWORK INTEGER WORK ARRAY USED TO STORE THE PIVOT ARRAY 2934C (SEE ABOVE) 2935C LIW DIMENSION OF IWORK IN THE CALLING PROGRAM 2936C JOB WORK REQUEST SWITCH 2937C = 0 DECOMPOSE THE MATRIX, SOLVE THE SYSTEM, AND 2938C COMPUTE THE SIGN OF THE DETERMINANT. 2939C = 1 SOLVE THE SYSTEM, ASSUMING THAT THE MATRIX HAS 2940C BEEN PREVIOUSLY DECOMPOSED. 2941C 2 = Factor system, compute determinant. 2942C 3 = Check jacobian matrix. Call user jacobian routine, 2943C multiply by -1.0, add finite difference jacobian, 2944C print largest entry. 2945C 2946C NVAR THE DIMENSION OF THE SYSTEM. 2947C RWORK REAL WORK ARRAY USED TO STORE THE MATRIX (SEE ABOVE). 2948C LRW THE DIMENSION OF RWORK IN THE CALLING PROGRAM. 2949C X A REAL VECTOR OF LENGTH NVAR, THE POINT AT WHICH THE 2950C JACOBIAN IS TO BE EVALUATED. 2951C Y A REAL VECTOR OF LENGTH NVAR. ON INPUT Y CONTAINS THE 2952C RIGHT HAND SIDE OF THE LINEAR SYSTEM (1), ON RETURN 2953C WITH IERROR = 0 THE SOLUTION OF THE SYSTEM IS GIVEN IN Y. 2954C 2955 DOUBLE PRECISION ONE 2956 DOUBLE PRECISION ZERO 2957C 2958 PARAMETER (ONE=1.0) 2959 PARAMETER (ZERO=0.0) 2960C 2961 EXTERNAL DENJAC 2962 EXTERNAL DF 2963 EXTERNAL FX 2964 EXTERNAL IDAMAX 2965 EXTERNAL DGEDI 2966 EXTERNAL DGEFA 2967 EXTERNAL DGESL 2968 EXTERNAL DSCAL 2969C 2970 INTRINSIC MOD 2971 INTRINSIC SQRT 2972C 2973 INTEGER LIW 2974 INTEGER LRW 2975 INTEGER NVAR 2976C 2977 DOUBLE PRECISION DET(2) 2978 DOUBLE PRECISION DETS 2979 DOUBLE PRECISION EPS 2980 DOUBLE PRECISION FPAR(*) 2981 INTEGER I 2982 INTEGER IERROR 2983 INTEGER IPAR(*) 2984 INTEGER IPC 2985 INTEGER IDAMAX 2986 INTEGER IWORK(LIW) 2987 INTEGER J 2988 INTEGER JAC 2989 INTEGER JACK 2990 INTEGER JOB 2991 INTEGER JOB2 2992 INTEGER K 2993 INTEGER LDF 2994 INTEGER LFXM 2995 INTEGER LFXP 2996 INTEGER LILST 2997 INTEGER LOUNIT 2998 INTEGER LPIV 2999 INTEGER LRLST 3000 INTEGER NDIM 3001 DOUBLE PRECISION RWORK(LRW) 3002 DOUBLE PRECISION SKALE 3003 DOUBLE PRECISION X(NVAR) 3004 DOUBLE PRECISION Y(NVAR) 3005C 3006C CHECK INPUT VALUES OF NVAR, AND IPC. 3007C DEPENDING ON VALUE OF JOB, EITHER SET UP 3008C AUGMENTED JACOBIAN, DECOMPOSE INTO L-U FACTORS, AND GET DETERMINANT, 3009C OR USE CURRENT FACTORED JACOBIAN WITH NEW RIGHT HAND SIDE. 3010C 3011 LOUNIT=IWORK(8) 3012 IERROR=0 3013 LPIV=IWORK(13) 3014 LDF=IWORK(15) 3015 LILST=LPIV+NVAR-1 3016 JAC=IWORK(9) 3017 IF(JAC.NE.0.AND.JOB.EQ.3)THEN 3018 IERROR=4 3019 WRITE(LOUNIT,*)'DENSLV - Error! Jacobian check requested' 3020 WRITE(LOUNIT,*)' but no user Jacobian routine.' 3021 RETURN 3022 ENDIF 3023 IF(JAC.EQ.0.AND.JOB.NE.3)THEN 3024 LRLST=LDF-1+NVAR*NVAR 3025 ELSE 3026 LFXP=LDF+NVAR*NVAR 3027 LFXM=LFXP+NVAR 3028 LRLST=LDF-1+NVAR*NVAR+2*NVAR 3029 ENDIF 3030 IF((LILST.GT.LIW).OR.(LRLST.GT.LRW))THEN 3031 IERROR=1 3032 WRITE(LOUNIT,1040)LILST,LIW 3033 WRITE(LOUNIT,1050)LRLST,LRW 3034 RETURN 3035 ENDIF 3036 IF(JOB.EQ.1)GO TO 10 3037C 3038C IF JOB IS 0 OR 2, THEN WE MUST EVALUATE THE JACOBIAN, USING THE USER 3039C ROUTINE OR FINITE DIFFERENCES, AND THEN FACTOR IT AND GET THE 3040C SIGN OF THE DETERMINANT. 3041C 3042 NDIM=NVAR*NVAR 3043 DO 5 I=1,NDIM 3044 RWORK(LDF+I-1)=ZERO 30455 CONTINUE 3046 IF(JAC.EQ.0)THEN 3047 CALL DF(NVAR,FPAR,IPAR,X,RWORK(LDF),IERROR) 3048 IWORK(19)=IWORK(19)+1 3049 RWORK(LDF+IPC*NVAR-1)=ONE 3050 ENDIF 3051 IF(JOB.EQ.3)THEN 3052 SKALE=-ONE 3053 CALL DSCAL(NDIM,SKALE,RWORK(LDF),1) 3054 ENDIF 3055 IF(JAC.EQ.1.OR.JAC.EQ.2.OR.JOB.EQ.3)THEN 3056 JACK=JAC 3057 IF(JOB.EQ.3)JACK=2 3058 EPS=SQRT(SQRT(RWORK(8))) 3059 CALL DENJAC(EPS,FPAR,RWORK(LDF),FX,IERROR,IPAR,IPC,IWORK, 3060 * JACK,LIW,LOUNIT,NVAR,X,RWORK(LFXP),RWORK(LFXM)) 3061 ENDIF 3062 IF(IERROR.NE.0)RETURN 3063C 3064C If JOB=3, print out largest entry of matrix. 3065C 3066 IF(JOB.EQ.3)THEN 3067 K=IDAMAX(NDIM,RWORK(LDF),1) 3068 I=MOD(K-1,NVAR)+1 3069 J=(K-I)/NVAR+1 3070 WRITE(LOUNIT,1070)RWORK(LDF+K-1),I,J 3071 IF(IWORK(1).EQ.(-2))THEN 3072 WRITE(LOUNIT,*)' ' 3073 WRITE(LOUNIT,*)'DENSLV - Entire difference matrix:' 3074 WRITE(LOUNIT,*)' ' 3075 DO 99 I=1,NVAR 3076 DO 98 J=1,NVAR 3077 K=LDF+(J-1)*NVAR+I-1 3078 WRITE(LOUNIT,1080)RWORK(K),I,J 307998 CONTINUE 3080 WRITE(LOUNIT,*)' ' 308199 CONTINUE 3082 ENDIF 3083 RETURN 3084 ENDIF 3085C 3086C DECOMPOSE MATRIX 3087C 3088 CALL DGEFA(RWORK(LDF),NVAR,NVAR,IWORK(LPIV),IERROR) 3089 IWORK(20)=IWORK(20)+1 3090 IF(IERROR.NE.0)THEN 3091 WRITE(LOUNIT,1492)IERROR 3092 IERROR=3 3093 RETURN 3094 ENDIF 3095C 3096C COMPUTE DETERMINANT 3097C 3098 JOB2=10 3099 CALL DGEDI(RWORK(LDF),NVAR,NVAR,IWORK(LPIV),DET,Y,JOB2) 3100C 3101C RECORD THE SIGN OF THE DETERMINANT 3102C 3103 DETS=ZERO 3104 IF(DET(1).GT.ZERO)THEN 3105 DETS=ONE 3106 ELSEIF(DET(1).LT.ZERO)THEN 3107 DETS=-ONE 3108 ENDIF 3109 IF(JOB.EQ.2)RETURN 3110C 3111C USING L-U FACTORED MATRIX, SOLVE SYSTEM USING FORWARD-BACKWARD 3112C ELIMINATION, AND OVERWRITE RIGHT HAND SIDE WITH SOLUTION 3113C 311410 CONTINUE 3115 CALL DGESL(RWORK(LDF),NVAR,NVAR,IWORK(LPIV),Y,0) 3116 IWORK(21)=IWORK(21)+1 3117 RETURN 31181040 FORMAT(' DENSLV - Need LIW=',I6,', have LIW=',I6) 31191050 FORMAT(' DENSLV - Need LRW=',I6,', have LRW=',I6) 31201070 FORMAT(' DENSLV - Maximum difference between user and estimated'/ 3121 * ' jacobian is ',G14.6,' row ',I6,' column ',I6) 31221080 FORMAT(1X,G14.6,' =FP(I,J)-DELF(I,J), I, J=',2I6) 31231492 FORMAT(' DENSLV - Zero pivot, DGEFA returns INFO=',I6) 3124 END 3125 SUBROUTINE DENJAC(EPS,FPAR,FPRIME,FX,IERROR,IPAR,IPC,IWORK, 3126 * JAC,LIW,LOUNIT,NVAR,X,WORK1,WORK2) 3127C 3128C************************************************************************ 3129C 3130C DENJAC ESTIMATES THE JACOBIAN FPRIME OF THE FUNCTION FX ASSUMING 3131C THAT FULL (DENSE) STORAGE IS USED. IT IS CALLED AUTOMATICALLY 3132C BY DENSLV WHEN JAC=1 OR 2. 3133C 3134C EPS - A TOLERANCE TO USE FOR SHIFTING THE X VALUES. 3135C EPS=SQRT(MACHINE EPSILON) FOR THE VAX, OR 3136C EPS=SQRT(SQRT(MACHINE EPSILON)) FOR THE PC 3137C HAVE BEEN FOUND SUITABLE. 3138C 3139C FPAR - USER REAL PARAMETER VECTOR, PASSED TO FX BUT NOT OTHERWISE 3140C REFERENCED. 3141C 3142C FPRIME - REAL(NVAR,NVAR), THE ARRAY INTO WHICH THE JACOBIAN WILL 3143C BE STORED, INCLUDING THE LAST AUXILLIARY ROW. 3144C 3145C FX - EXTERNAL NAME OF USER ROUTINE DEFINING FUNCTION, OF THE 3146C FORM SUBROUTINE FX(NVAR,FPAR,IPAR,X,F,IERROR). 3147C 3148C IERROR - ERROR RETURN FLAG. IERROR NONZERO MEANS THERE WAS A PROBLEM 3149C IN THE USER ROUTINE FX OR IN DENJAC. 3150C 3151C IPAR - USER INTEGER PARAMETER VECTOR, PASSED TO FX BUT NOT OTHERWISE 3152C REFERENCED. 3153C 3154C IPC - INDEX OF CURRENT CONTINUATION PARAMETER, DETERMINES FORM OF 3155C LAST ROW OF JACOBIAN. 3156C 3157C IWORK - INTEGER WORK VECTOR OF DIMENSION LIW. ONLY REQUIRED IN ORDER 3158C TO UPDATE IWORK(22), WHICH COUNTS FUNCTION EVALUATIONS. 3159C 3160C JAC - USER DEFINED JACOBIAN OPTION. 0 MEANS USER SUPPLIES JACOBIAN. 3161C 1 MEANS USE FORWARD, 2 MEANS USE CENTRAL DIFFERENCES TO 3162C ESTIMATE JACOBIAN FROM FX. 3163C 3164C LIW - DIMENSION OF IWORK. 3165C 3166C LOUNIT - OUTPUT UNIT. 3167C 3168C NVAR - DIMENSION OF X, FPRIME, WORK1, WORK2. THE NUMBER OF VARIABLES. 3169C 3170C X - POINT AT WHICH JACOBIAN IS TO BE ESTIMATED. 3171C 3172C WORK1 - WORK VECTOR OF DIMENSION NVAR. 3173C 3174C WORK2 - WORK VECTOR OF DIMENSION NVAR. 3175C 3176 DOUBLE PRECISION ONE 3177 DOUBLE PRECISION ZERO 3178C 3179 PARAMETER (ONE=1.0) 3180 PARAMETER (ZERO=0.0) 3181C 3182 EXTERNAL FX 3183 EXTERNAL DAXPY 3184 EXTERNAL DSCAL 3185C 3186 INTRINSIC ABS 3187C 3188 INTEGER LIW 3189 INTEGER NVAR 3190C 3191 DOUBLE PRECISION DELM 3192 DOUBLE PRECISION DELP 3193 DOUBLE PRECISION EPS 3194 DOUBLE PRECISION FPAR(*) 3195 DOUBLE PRECISION FPRIME(NVAR,NVAR) 3196 INTEGER IPAR(*) 3197 INTEGER IERROR 3198 INTEGER IPC 3199 INTEGER IWORK(LIW) 3200 INTEGER J 3201 INTEGER JAC 3202 INTEGER LOUNIT 3203 DOUBLE PRECISION SKALE 3204 DOUBLE PRECISION X(NVAR) 3205 DOUBLE PRECISION XSAVE 3206 DOUBLE PRECISION WORK1(NVAR) 3207 DOUBLE PRECISION WORK2(NVAR) 3208C 3209 IF(JAC.EQ.1)THEN 3210 CALL FX(NVAR,FPAR,IPAR,X,WORK2,IERROR) 3211 DELM=ZERO 3212 IWORK(22)=IWORK(22)+1 3213 IF(IERROR.NE.0)RETURN 3214 ENDIF 3215C 3216 DO 20 J=1,NVAR 3217 XSAVE=X(J) 3218 DELP=EPS*(ONE+ABS(X(J))) 3219 X(J)=X(J)+DELP 3220 CALL FX(NVAR,FPAR,IPAR,X,WORK1,IERROR) 3221 IWORK(22)=IWORK(22)+1 3222 IF(IERROR.NE.0)RETURN 3223 IF(JAC.EQ.2)THEN 3224 DELM=-DELP 3225 X(J)=XSAVE+DELM 3226 CALL FX(NVAR,FPAR,IPAR,X,WORK2,IERROR) 3227 IWORK(22)=IWORK(22)+1 3228 IF(IERROR.NE.0)RETURN 3229 ENDIF 3230 X(J)=XSAVE 3231C 3232C Compute DFDX(*,J)= (F(X+)-F(X-))/DELX. 3233C Note that we contrive to ADD this quantity to DFDX, 3234C This makes it possible to use this routine also to check 3235C the accuracy of the user's jacobian routine. 3236C 3237 SKALE=-ONE 3238 CALL DAXPY(NVAR-1,SKALE,WORK2,1,WORK1,1) 3239 SKALE=ONE/(DELP-DELM) 3240 CALL DSCAL(NVAR-1,SKALE,WORK1,1) 3241 SKALE=ONE 3242 CALL DAXPY(NVAR-1,SKALE,WORK1,1,FPRIME(1,J),1) 324320 CONTINUE 3244C 3245 FPRIME(NVAR,IPC)=FPRIME(NVAR,IPC)+ONE 3246 RETURN 3247 END 3248 SUBROUTINE BANSLV(DETS,FX,DF,FPAR,IERROR,IPC,IPAR,IWORK,LIW, 3249 1 JOB,NVAR,RWORK,LRW,X,Y) 3250C 3251C*********************************************************************** 3252C 3253C BANSLV SOLVES A LINEAR SYSTEM OF THE FORM 3254C 3255C ( DFDY DFDZ ) 3256C ( ) * X=B 3257C ( E(IPC) ) 3258C 3259C WHERE DFDY IS FORMALLY AN (NVAR-1) BY (NVAR-1) MATRIX WHICH IS 3260C TO BE STORED IN LINPACK GENERAL BAND STORAGE, DFDZ IS A COLUMN 3261C VECTOR, AND E(IPC) IS AN NVAR-DIMENSIONAL ROW VECTOR WHICH IS 0 3262C EXCEPT FOR A 1 IN THE IPC-TH POSITION. IT IS ASSUMED THAT THE 3263C FULL MATRIX ABOVE REPRESENTS THE JACOBIAN OF A SET OF NONLINEAR 3264C EQUATIONS, WITH THE LAST EQUATION BEING X(IPC)=CONSTANT. 3265C 3266C THE PURPOSE OF THIS ROUTINE IS TO SOLVE THE FULL SYSTEM, WHILE 3267C TAKING ADVANTAGE OF SPECIAL PROPERTIES OF THE (NVAR-1) SUBSYSTEM 3268C (IN THIS CASE BANDEDNESS), WHICH WOULD BE LOST IF THE FULL 3269C SYSTEM WAS SOLVED DIRECTLY. 3270C 3271C SEE THE PAPER BY CHAN LISTED ABOVE FOR DETAILS. 3272C 3273C SUBROUTINE ARGUMENTS 3274C 3275C DETS - THE SIGN OF THE DETERMINANT OF THE FULL MATRIX. 3276C 3277C FX - AN EXTERNAL ROUTINE, THE NAME OF THE USER SUPPLIED 3278C SUBROUTINE WHICH COMPUTES THE NVAR-1 DIMENSIONAL 3279C FUNCTION. 3280C 3281C DF - AN EXTERNAL ROUTINE, THE NAME OF THE USER SUPPLIED 3282C SUBROUTINE WHICH COMPUTES THE JACOBIAN OF THE NVAR-1 3283C NONLINEAR FUNCTIONS. THE LAST, AUGMENTING FUNCTION, 3284C IS HANDLED BY THIS ROUTINE. BECAUSE OF THE SPECIAL 3285C STORAGE ARRANGEMENTS, GREAT CARE MUST BE TAKEN TO 3286C STORE INFORMATION PROPERLY. 3287C 3288C LET MU AND ML BE THE UPRER AND LOWER BANDWIDTHS. 3289C THESE VALUES MUST BE STORED IN IPAR(*) AND IPAR(2), 3290C RESPECTIVELY. SET NBAND=2*ML+MU+1. 3291C THE INDEX K OF RWORK IN WHICH WE ARE TO STORE THE 3292C (I,J) ENTRY OF THE JACOBIAN IS DETERMINED AS FOLLOWS: 3293C 3294C IF(J.EQ.NVAR) K=(NVAR-1)*NBAND+I 3295C 3296C OTHERWISE 3297C 3298C IF((J-I.LE.ML).AND.(I-J.LE.MU)) K=I+J*(NBAND-1)-ML 3299C 3300C FOR ALL OTHER VALUES OF I AND J THE ENTRY IS ZERO AND 3301C NOT STORED. 3302C 3303C THE FORM OF THE CALLING SEQUENCE OF DF IS 3304C SUBROUTINE DF(NVAR,FPAR,IPAR,X,RWORK,IERR) 3305C DIMENSION FPAR(*),IPAR(*),X(NVAR),RWORK(1) 3306C 3307C SET IERR NONZERO IF EXECUTION IS TO BE HALTED FOR ANY REASON 3308C 3309C FPAR - A REAL ARRAY FOR TRANSMISSION OF PARAMETERS FROM THE 3310C USER CALLING PROGRAM THROUGH THE CONTINUATION CODE TO 3311C USER SUBROUTINES. THE CONTINUATION CODE ASSUMES A 3312C DIMENSION OF 1 FOR FPAR, AND NEVER ACCESSES ANY ENTRIES. 3313C 3314C IERROR - ERROR RETURN FLAG. 3315C IERROR= 0, NO ERRORS DETECTED. 3316C IERROR= 1, DATA OR STORAGE ERROR. 3317C ILLEGAL VALUES OF NVAR, IPC, MU, ML, OR 3318C INSUFFICIENT STORAGE IN RWORK OR IWORK. 3319C IERROR= 2, ERROR RETURN FROM DF. 3320C THE USER ROUTINE DF HAS SET AN ERROR FLAG. 3321C IERROR= 3, NUMERICALLY SINGULAR MATRIX DETECTED. 3322C EITHER THE SUBMATRIX AS MODIFIED BY THIS 3323C ROUTINE IS SINGULAR, AND THE DECOMPOSITION 3324C OR THE CHOICE OF IPC IS INAPPROPRIATE, 3325C OR THE FULL SYSTEM MAY ACTUALLY BE SINGULAR. 3326C 3327C IPC - AN INDEX DETERMINED BY THE CONTINUATION CODE, WHICH 3328C DEFINES THE AUGMENTING EQUATION, AND HENCE THE FORM 3329C OF THE AUGMENTING ROW OF THE JACOBIAN. THE FULL 3330C ALGORITHM EMPLOYED HERE IS REALLY ONLY REQUIRED 3331C WHEN IPC IS NOT EQUAL TO NVAR. OTHERWISE, THE SOLUTION 3332C IS QUITE EASY. 3333C 3334C IPAR - AN INTEGER ARRAY SIMILAR TO RPAR, EXCEPT THAT LOCATIONS 3335C 1 AND 2 OF IPAR ARE REFERENCED BY BANSLV. IPAR(*)=ML, 3336C THE LOWER BANDWIDTH OF THE JACOBIAN, AND IPAR(2)=MU, 3337C THE UPPER BANDWIDTH. 3338C 3339C IWORK - AN INTEGER WORK ARRAY, USED BY THE CONTINUATION CODE 3340C TO STORE STATISTICS, POINTERS, AND THE PIVOT VECTOR. 3341C NOTE THAT IWORK(13) STORES THE FIRST ENTRY IN IWORK 3342C DEVOTED TO THE PIVOT VECTOR, IWORK(15) STORES THE 3343C FIRST ELEMENT OF RWORK DEVOTED TO THE JACOBIAN SUBSYSTEM. 3344C IWORK(20) COUNTS MATRIX FACTORIZATIONS, AND IWORK(21) 3345C COUNTS BACK-SOLVES. 3346C 3347C LIW - DIMENSION OF IWORK IN THE CALLING PROGRAM. 3348C 3349C JOB - WORK CONTROL SWITCH. 3350C JOB=0, DECOMPOSE MATRIX, GET DETERMINANT, SOLVE SYSTEM. 3351C JOB=1, SOLVE A NEW SYSTEM WITH SAME OLD MATRIX. 3352C JOB=2, DECOMPOSE, COMPUTE DETERMINANT. 3353C 3, Check jacobian matrix. Call user jacobian routine, 3354C multiply by -1.0, add finite difference jacobian, 3355C print largest entry. 3356C 3357C NVAR - THE NUMBER OF VARIABLES X, THE FORMAL DIMENSION OF 3358C THE FULL LINEAR SYSTEM. 3359C 3360C RWORK - A REAL ARRAY USED FOR STORING THE MATRIX AND THE TWO 3361C AUXILIARY VECTORS NEEDED. 3362C 3363C LRW - THE DIMENSION OF RWORK IN THE CALLING PROGRAM. 3364C 3365C X - A REAL VECTOR OF LENGTH NVAR, THE POINT AT WHICH 3366C THE FULL LINEAR SYSTEM IS TO BE EVALUATED. 3367C 3368C Y - A REAL VECTOR OF LENGTH NVAR, WHICH ON INPUT IS 3369C THE RIGHT HAND SIDE OF THE LINEAR SYSTEM TO BE SOLVED, 3370C AND ON SUCCESSFUL RETURN WITH IERROR=0, IS THE 3371C SOLUTION OF THAT LINEAR SYSTEM. 3372C 3373C 3374 DOUBLE PRECISION ONE 3375 DOUBLE PRECISION ZERO 3376C 3377 PARAMETER (ONE=1.0) 3378 PARAMETER (ZERO=0.0) 3379C 3380 EXTERNAL BANJAC 3381 EXTERNAL DF 3382 EXTERNAL FX 3383 EXTERNAL IDAMAX 3384 EXTERNAL DAXPY 3385 EXTERNAL DDOT 3386 EXTERNAL DGBDI 3387 EXTERNAL DGBFA 3388 EXTERNAL DGBSL 3389 EXTERNAL DSCAL 3390 EXTERNAL DSWAP 3391C 3392 INTRINSIC MAX 3393 INTRINSIC MIN 3394 INTRINSIC SQRT 3395C 3396 INTEGER LIW 3397 INTEGER LRW 3398 INTEGER NVAR 3399C 3400 DOUBLE PRECISION AK 3401 DOUBLE PRECISION DET(2) 3402 DOUBLE PRECISION DETS 3403 DOUBLE PRECISION EPS 3404 DOUBLE PRECISION FPAR(*) 3405 INTEGER I 3406 INTEGER IERROR 3407 INTEGER IPAR(*) 3408 INTEGER IPC 3409 INTEGER IRL 3410 INTEGER IRU 3411 INTEGER IDAMAX 3412 INTEGER ITEMP 3413 INTEGER IWORK(LIW) 3414 INTEGER J 3415 INTEGER JAC 3416 INTEGER JACK 3417 INTEGER JOB 3418 INTEGER JTEMP 3419 INTEGER K 3420 INTEGER LDFL 3421 INTEGER LDFX 3422 INTEGER LDX 3423 INTEGER LFXM 3424 INTEGER LFXP 3425 INTEGER LILST 3426 INTEGER LOUNIT 3427 INTEGER LPIV 3428 INTEGER LRLST 3429 INTEGER LROWIP 3430 INTEGER MBAND 3431 INTEGER ML 3432 INTEGER MU 3433 INTEGER NBAND 3434 INTEGER NDIM 3435 INTEGER NEQN 3436 INTEGER NSWAP 3437 DOUBLE PRECISION RWORK(LRW) 3438 DOUBLE PRECISION DDOT 3439 DOUBLE PRECISION SKALE 3440 DOUBLE PRECISION TEMP 3441 DOUBLE PRECISION X(NVAR) 3442 DOUBLE PRECISION Y(NVAR) 3443C 3444C CHECK INPUT DATA 3445C 3446 IERROR=0 3447 LOUNIT=IWORK(8) 3448 ML=IPAR(1) 3449 MU=IPAR(2) 3450 NEQN=NVAR-1 3451 IF(ML.LT.0.OR.ML.GT.NEQN)THEN 3452 IERROR=1 3453 WRITE(LOUNIT,1020)ML 3454 RETURN 3455 ENDIF 3456 IF(MU.LT.0.OR.MU.GT.NEQN)THEN 3457 IERROR=1 3458 WRITE(LOUNIT,1030)MU 3459 RETURN 3460 ENDIF 3461 MBAND=ML+MU+1 3462 NBAND=MBAND+ML 3463 LPIV=IWORK(13) 3464 LILST=LPIV+NEQN-1 3465 LDFX=IWORK(15) 3466 LDFL=LDFX+NBAND*NEQN 3467 LROWIP=LDFL+NVAR 3468 JAC=IWORK(9) 3469 IF(JAC.EQ.0.AND.JOB.NE.3)THEN 3470 LRLST=LROWIP+NVAR-1 3471 ELSE 3472 LFXP=LROWIP+NVAR 3473 LFXM=LFXP+NVAR 3474 LDX=LFXM+NVAR 3475 LRLST=LDX+NVAR-1 3476 ENDIF 3477 NDIM=NEQN*NBAND+NVAR+NVAR-1 3478 IF(LILST.GT.LIW.OR.LRLST.GT.LRW)THEN 3479 WRITE(LOUNIT,1040)LILST,LIW 3480 WRITE(LOUNIT,1050)LRLST,LRW 3481 IERROR=1 3482 RETURN 3483 ENDIF 3484 IF(JOB.EQ.1)THEN 3485 IF(IPC.EQ.NVAR)THEN 3486 GO TO 70 3487 ELSE 3488 GO TO 40 3489 ENDIF 3490 ENDIF 3491C 3492C JOB=0 OR 2 MEANS WE MUST EVALUATE THE JACOBIAN, EITHER WITH THE USER'S 3493C ROUTINE OR BY FINITE DIFFERENCES, FACTOR IT AND GET THE SIGN 3494C OF THE DETERMINANT 3495C 3496 DO 10 I=1,NDIM 3497 RWORK(LDFX+I-1)=ZERO 349810 CONTINUE 3499 IF(JAC.EQ.0)THEN 3500 CALL DF(NVAR,FPAR,IPAR,X,RWORK(LDFX),IERROR) 3501 IWORK(19)=IWORK(19)+1 3502 IF(IERROR.NE.0)RETURN 3503 RWORK(LROWIP-1+IPC)=ONE 3504 ENDIF 3505 IF(JOB.EQ.3)THEN 3506 SKALE=-ONE 3507 CALL DSCAL(NDIM,SKALE,RWORK(LDFX),1) 3508 ENDIF 3509 IF(JAC.EQ.1.OR.JAC.EQ.2.OR.JOB.EQ.3)THEN 3510 JACK=JAC 3511 IF(JOB.EQ.3)JACK=2 3512 EPS=SQRT(SQRT(RWORK(8))) 3513 CALL BANJAC(EPS,RWORK(LDFL),FPAR,RWORK(LDFX),RWORK(LROWIP), 3514 * FX,IERROR,IPAR,IPC,IWORK,JACK,LIW,LOUNIT,NBAND,NEQN,NVAR,X, 3515 * RWORK(LDX),RWORK(LFXP),RWORK(LFXM)) 3516 ENDIF 3517C 3518C Figure out formulas for I and J 3519C 3520 IF(JOB.EQ.3)THEN 3521 K=IDAMAX(NDIM,RWORK(LDFX),1) 3522 AK=RWORK(LDFX+K-1) 3523 IF(K.LE.NEQN*NBAND)THEN 3524 J=((K-1)/NBAND)+1 3525 I=K-(J-1)*NBAND+J-ML-MU-1 3526 ELSEIF(K.LE.NEQN*NBAND+NEQN)THEN 3527 I=K-NEQN*NBAND 3528 J=NVAR 3529 ELSE 3530 I=NVAR 3531 J=K-NEQN*NBAND-NEQN 3532 ENDIF 3533 WRITE(LOUNIT,1070)AK,I,J 3534 IF(IWORK(1).EQ.-2)THEN 3535 WRITE(LOUNIT,*)' ' 3536 WRITE(LOUNIT,*)'BANSLV - Entire difference matrix:' 3537 WRITE(LOUNIT,*)' ' 3538 DO 99 I=1,NVAR 3539 DO 98 J=1,NVAR 3540 IF(J.EQ.NVAR)THEN 3541 K=LDFX-1+(NVAR-1)*NBAND+I 3542 WRITE(LOUNIT,1080)RWORK(K),I,J 3543 ELSEIF((J-I.LE.ML).AND.(I-J.LE.MU))THEN 3544 K=LDFX-1+I+J*(NBAND-1)-ML 3545 WRITE(LOUNIT,1080)RWORK(K),I,J 3546 ENDIF 354798 CONTINUE 3548 WRITE(LOUNIT,*)' ' 354999 CONTINUE 3550 ENDIF 3551 RETURN 3552 ENDIF 3553 IF(IPC.EQ.NVAR)GO TO 60 3554C 3555C SWITCH THE NVAR-TH AND IPC-TH ROWS. 3556C 3557 IRL=MAX(1,IPC-ML) 3558 IRU=MIN(NEQN,IPC+MU) 3559 NSWAP=IRU+1-IRL 3560 ITEMP=LDFX-ML-1+IPC+IRL*(NBAND-1) 3561 JTEMP=NBAND-1 3562 CALL DSWAP(NSWAP,RWORK(LROWIP-1+IRL),1,RWORK(ITEMP),JTEMP) 3563C 3564C DECOMPOSE THE SUBMATRIX AND OBTAIN DETERMINANT SIGN 3565C 3566 CALL DGBFA(RWORK(LDFX),NBAND,NEQN,ML,MU,IWORK(LPIV),IERROR) 3567 IWORK(20)=IWORK(20)+1 3568 IF(IERROR.NE.0)THEN 3569 WRITE(LOUNIT,1000)IERROR,IPC 3570 IERROR=3 3571 RETURN 3572 ENDIF 3573 CALL DGBDI(RWORK(LDFX),NBAND,NEQN,ML,MU,IWORK(LPIV),DET) 3574 DETS=ZERO 3575 IF(DET(1).GT.ZERO)THEN 3576 DETS=-ONE 3577 ELSEIF(DET(1).LT.ZERO)THEN 3578 DETS=ONE 3579 ENDIF 3580C 3581C SET THE RIGHT HAND SIDE OF THE AUXILLIARY SYSTEM TO 3582C THE LAST COLUMN OF THE JACOBIAN, MINUS E(IPC). SHUFFLE THE 3583C IPC-TH AND NVAR-TH ENTRIES OF THIS RIGHT HAND SIDE 3584C TO REFLECT THE PIVOTING OF THE EQUATIONS, THEN SOLVE. 3585C 3586 TEMP=RWORK(LDFL+IPC-1)-ONE 3587 RWORK(LDFL+IPC-1)=RWORK(LDFL+NVAR-1) 3588 RWORK(LDFL+NVAR-1)=TEMP 3589 CALL DGBSL(RWORK(LDFX),NBAND,NEQN,ML,MU,IWORK(LPIV),RWORK(LDFL) 3590 * ,0) 3591 IWORK(21)=IWORK(21)+1 3592C 3593C SOLVE FOR LAST ENTRY OF AUXILLIARY SOLUTION 3594C 3595 RWORK(LDFL+NVAR-1)=RWORK(LDFL+NVAR-1) 3596 * -DDOT(NEQN,RWORK(LROWIP),1,RWORK(LDFL),1) 3597C 3598C ADJUST DETERMINANT SIGN 3599C 3600 IF(ONE+RWORK(LDFL+NVAR-1).EQ.ZERO)THEN 3601 IERROR=3 3602 WRITE(LOUNIT,*)'BANSLV - Algorithm fails, DENOM=0.0' 3603 RETURN 3604 ENDIF 3605 IF(ONE+RWORK(LDFL+NVAR-1).LT.ZERO)DETS=-DETS 3606 IF(JOB.EQ.2)RETURN 3607C 3608C SOLVE SYSTEM 3609C 361040 CONTINUE 3611 IF(IPC.EQ.NVAR)GO TO 70 3612C 3613C MODIFY RIGHT HAND SIDE OF MAIN SYSTEM 3614C 3615 TEMP=Y(IPC) 3616 Y(IPC)=Y(NVAR) 3617 Y(NVAR)=TEMP 3618C 3619C SOLVE SUBSYSTEM 3620C 3621 CALL DGBSL(RWORK(LDFX),NBAND,NEQN,ML,MU,IWORK(LPIV),Y,0) 3622 IWORK(21)=IWORK(21)+1 3623C 3624C SOLVE FOR LAST ENTRY OF MAIN SOLUTION 3625C 3626 Y(NVAR)=Y(NVAR)-DDOT(NEQN,RWORK(LROWIP),1,Y,1) 3627C 3628C CORRECT MAIN SOLUTION WITH MULTIPLE OF SUBSOLUTION 3629C 3630 SKALE=-Y(NVAR)/(ONE+RWORK(LDFL+NVAR-1)) 3631 CALL DAXPY(NVAR,SKALE,RWORK(LDFL),1,Y,1) 3632 RETURN 3633C 3634C FACTOR MATRIX FOR THE SPECIAL CASE IPC=NVAR 3635C 3636 60 CONTINUE 3637 CALL DGBFA(RWORK(LDFX),NBAND,NEQN,ML,MU,IWORK(LPIV),IERROR) 3638 IWORK(20)=IWORK(20)+1 3639 IF(IERROR.NE.0)THEN 3640 WRITE(LOUNIT,1000)IERROR,IPC 3641 IERROR=3 3642 RETURN 3643 ENDIF 3644 CALL DGBDI(RWORK(LDFX),NBAND,NEQN,ML,MU,IWORK(LPIV),DET) 3645 DETS=ZERO 3646 IF(DET(1).LT.ZERO)THEN 3647 DETS=-ONE 3648 ELSEIF(DET(1).GT.ZERO)THEN 3649 DETS=ONE 3650 ENDIF 3651 IF(JOB.EQ.2)RETURN 3652C 3653C SOLVE SYSTEM, FOR SPECIAL CASE IPC=NVAR 3654C 3655 70 CONTINUE 3656 SKALE=-Y(NVAR) 3657 CALL DAXPY(NEQN,SKALE,RWORK(LDFL),1,Y,1) 3658 CALL DGBSL(RWORK(LDFX),NBAND,NEQN,ML,MU,IWORK(LPIV),Y,0) 3659 IWORK(21)=IWORK(21)+1 3660 RETURN 36611000 FORMAT(' BANSLV - Zero pivot, DGBFA returns INFO=',I6,' IPC=',I6) 36621020 FORMAT(' BANSLV - Illegal ML=',I6) 36631030 FORMAT(' BANSLV - Illegal MU=',I6) 36641040 FORMAT(' BANSLV - Need LIW=',I6,', have LIW=',I6) 36651050 FORMAT(' BANSLV - Need LRW=',I6,', have LRW=',I6) 36661070 FORMAT(' BANSLV - Maximum difference between user and estimated'/ 3667 * ' jacobian is ',G14.6,' row ',I6,' column ',I6) 36681080 FORMAT(1X,G14.6,' =FP(I,J)-DELF(I,J), I, J=',2I6) 3669 END 3670 SUBROUTINE BANJAC(EPS,FCOL,FPAR,FPRIME,FROW,FX,IERROR,IPAR,IPC, 3671 * IWORK,JAC,LIW,LOUNIT,NBAND,NEQN,NVAR,X,XTEMP,WORK1,WORK2) 3672C 3673C********************************************************************** 3674C 3675C BANJAC estimates the jacobian matrix FPRIME of the function FX, 3676C using forward or central finite differences. BANJAC is called by 3677C BANSLV when the user has specified the jacobian option as 1 or 2. 3678C 3679C EPS Input, a tolerance used for shifting the X values. A value 3680C of the square root of the machine precision is usually 3681C appropriate. 3682C 3683C FCOL Output, REAL FCOL(NEQN), the last column of the approximate 3684C jacobian, which is allowed to be "full". This comprises 3685C matrix entries FPRIME(1,NVAR) through FPRIME(NEQN,NVAR). 3686C 3687C FPAR Input/output, REAL FPAR(*), user parameter vector, not 3688C touched by this routine, but passed on to user routines. 3689C 3690C FPRIME Output, REAL FPRIME(NBAND,NEQN), is the array into which the 3691C the banded portion of the computed jacobian will be stored. 3692C The LINPACK general band format is used, assigning entry (I,J) 3693C to FPRIME(I-J+ML+MU+1,J), where ML and MU are the lower and 3694C upper half bandwidths respectively. 3695C 3696C FROW Output, REAL FROW(NVAR), storage for the last (augmenting) row 3697C of the jacobian, which will be all zero except for a 1 in 3698C location IPC. 3699C 3700C FX Input, EXTERNAL FX, the name of the routine which evaluates the 3701C function. 3702C 3703C FX computes the value of the nonlinear function. This name must be 3704C declared EXTERNAL in the calling program. FX should evaluate the 3705C NVAR-1 function components at the input point X, and store the result 3706C in the vector FVEC. An augmenting equation will be stored in entry 3707C NVAR of FVEC by the PITCON program. 3708C 3709C FX should have the following form: 3710C 3711C SUBROUTINE FX(NVAR,FPAR,IPAR,X,FVEC,IERROR) 3712C 3713C NVAR Input, INTEGER NVAR, number of variables. 3714C 3715C FPAR Input/output, REAL FPAR(*), array of user parameters. 3716C 3717C IPAR Input/output, INTEGER IPAR(*), array of user parameters. 3718C 3719C X Input, REAL X(NVAR), the point at which function evaluation 3720C is required. 3721C 3722C FVEC Output, REAL FVEC(NVAR), the value of the function at point 3723C X. Only the first NVAR-1 entries of FVEC are to be set by 3724C the routine. PITCON sets the final value itself. 3725C 3726C IERROR Output, INTEGER IERROR, error flag. FX should set this to 0 3727C if there are no problems, or to 2 if there is a problem. 3728C 3729C IERROR Output, INTEGER IERROR, error flag. A nonzero value means that 3730C there was an error in the user routine FX, or in BANJAC itself. 3731C In either case, the jacobian has not been computed. 3732C 3733C IPAR Input, INTEGER IPAR(*), a user parameter vector passed to FX. 3734C However, because this is a problem with a banded jacobian, entries 3735C IPAR(1) and IPAR(2) are read by this routine. IPAR(1) contains 3736C ML, the lower half bandwidth of the jacobian, and IPAR(2) contains 3737C MU, the upper half bandwidth of the jacobian. 3738C 3739C IPC Input, INTEGER IPC, the index of the current continuation parameter, 3740C which is needed to determine the form of FROW. 3741C 3742C IWORK Input, INTEGER IWORK(LIW), work and statistics vector. Only 3743C required here so that we can count the number of function 3744C evaluations. 3745C 3746C JAC Input, INTEGER JAC, the user requested jacobian option. For 3747C our purposes, the only two values of interest are: 3748C 3749C 1 = estimate jacobian with forward differences, 3750C 2 = estimate jacobian with central differences (twice the work) 3751C 3752C LIW Input, INTEGER LIW, the dimension of IWORK. 3753C 3754C LOUNIT Input, INTEGER LOUNIT, the FORTRAN output unit for messages. 3755C 3756C NBAND Input, INTEGER NBAND, the first dimension of the jacobian matrix 3757C FPRIME, NBAND=ML+MU+1. 3758C 3759C NEQN Input, INTEGER NEQN, the number of equations, equal to NVAR-1. 3760C 3761C NVAR Input, INTEGER NVAR, the number of variables. 3762C 3763C X Input, REAL X(NVAR), the point at which the jacobian is desired. 3764C 3765C XTEMP, 3766C WORK1, 3767C WORK2 Work arrays, REAL XTEMP(NVAR), WORK1(NVAR), WORK2(NVAR). 3768C 3769 DOUBLE PRECISION ONE 3770 DOUBLE PRECISION TWO 3771C 3772 PARAMETER (ONE=1.0) 3773 PARAMETER (TWO=2.0) 3774C 3775 EXTERNAL FX 3776 EXTERNAL DAXPY 3777 EXTERNAL DCOPY 3778 EXTERNAL DSCAL 3779C 3780 INTRINSIC ABS 3781 INTRINSIC MAX 3782 INTRINSIC MIN 3783C 3784 INTEGER LIW 3785 INTEGER NBAND 3786 INTEGER NEQN 3787 INTEGER NVAR 3788C 3789 DOUBLE PRECISION EPS 3790 DOUBLE PRECISION FCOL(NEQN) 3791 DOUBLE PRECISION FPAR(*) 3792 DOUBLE PRECISION FPRIME(NBAND,NEQN) 3793 DOUBLE PRECISION FROW(NVAR) 3794 INTEGER IBAND 3795 INTEGER IERROR 3796 INTEGER IHI 3797 INTEGER ILO 3798 INTEGER IPAR(2) 3799 INTEGER IPC 3800 INTEGER IROW 3801 INTEGER IWORK(LIW) 3802 INTEGER J 3803 INTEGER JAC 3804 INTEGER KCALL 3805 INTEGER LOUNIT 3806 INTEGER MBAND 3807 INTEGER ML 3808 INTEGER MU 3809 DOUBLE PRECISION SKALE 3810 DOUBLE PRECISION X(NVAR) 3811 DOUBLE PRECISION XJAC 3812 DOUBLE PRECISION XTEMP(NVAR) 3813 DOUBLE PRECISION WORK1(NVAR) 3814 DOUBLE PRECISION WORK2(NVAR) 3815C 3816 ML=IPAR(1) 3817 MU=IPAR(2) 3818 MBAND=ML+MU+1 3819 IF(JAC.EQ.1)THEN 3820 CALL FX(NVAR,FPAR,IPAR,X,WORK2,IERROR) 3821 IWORK(22)=IWORK(22)+1 3822 IF(IERROR.NE.0)RETURN 3823 ENDIF 3824 XJAC=ONE 3825 IF(JAC.EQ.2)XJAC=TWO 3826 DO 40 KCALL=1,MBAND 3827 CALL DCOPY(NVAR,X,1,XTEMP,1) 3828 DO 10 J=KCALL,NEQN,MBAND 3829 XTEMP(J)=X(J)+EPS*(ONE+ABS(X(J))) 383010 CONTINUE 3831 CALL FX(NVAR,FPAR,IPAR,XTEMP,WORK1,IERROR) 3832 IWORK(22)=IWORK(22)+1 3833 IF(IERROR.NE.0)RETURN 3834 IF(JAC.EQ.2)THEN 3835 CALL DCOPY(NVAR,X,1,XTEMP,1) 3836 DO 20 J=KCALL,NEQN,MBAND 3837 XTEMP(J)=X(J)-EPS*(ONE+ABS(X(J))) 383820 CONTINUE 3839 CALL FX(NVAR,FPAR,IPAR,XTEMP,WORK2,IERROR) 3840 IWORK(22)=IWORK(22)+1 3841 IF(IERROR.NE.0)RETURN 3842 ENDIF 3843 DO 30 J=KCALL,NEQN,MBAND 3844 ILO=MAX(1,J-MU) 3845 IHI=MIN(NEQN,J+ML) 3846 IROW=ILO-J+ML+MU+1 3847 IBAND=IHI-ILO+1 3848 SKALE=-ONE 3849 CALL DAXPY(IBAND,SKALE,WORK2(ILO),1,WORK1(ILO),1) 3850 SKALE=ONE/(XJAC*EPS*(ONE+ABS(X(J)))) 3851 CALL DSCAL(IBAND,SKALE,WORK1(ILO),1) 3852 SKALE=ONE 3853 CALL DAXPY(IBAND,SKALE,WORK1(ILO),1,FPRIME(IROW,J),1) 385430 CONTINUE 385540 CONTINUE 3856C 3857C Compute last column of jacobian, rows 1 to NEQN 3858C 3859 CALL DCOPY(NVAR,X,1,XTEMP,1) 3860 XTEMP(NVAR)=X(NVAR)+EPS*(ONE+ABS(X(NVAR))) 3861 CALL FX(NVAR,FPAR,IPAR,XTEMP,WORK1,IERROR) 3862 IWORK(22)=IWORK(22)+1 3863 IF(IERROR.NE.0)RETURN 3864 IF(JAC.EQ.2)THEN 3865 XTEMP(NVAR)=X(NVAR)-EPS*(ONE+ABS(X(NVAR))) 3866 CALL FX(NVAR,FPAR,IPAR,XTEMP,WORK2,IERROR) 3867 IWORK(22)=IWORK(22)+1 3868 IF(IERROR.NE.0)RETURN 3869 ENDIF 3870 SKALE=-ONE 3871 CALL DAXPY(NEQN,SKALE,WORK2,1,WORK1,1) 3872 SKALE=ONE/(XJAC*EPS*(ONE+ABS(X(NVAR)))) 3873 CALL DSCAL(NEQN,SKALE,WORK1,1) 3874 SKALE=ONE 3875 CALL DAXPY(NEQN,SKALE,WORK1,1,FCOL,1) 3876C 3877C Do last row, J=1,NVAR 3878C 3879 FROW(IPC)=FROW(IPC)+ONE 3880 RETURN 3881 END 3882 FUNCTION REPS() 3883C 3884C******************************************************************** 3885C 3886C Compute the machine relative precision. 3887C 3888C If you would prefer a simpler computation, you can replace this routine 3889C by a line that reads "REPS=1.0E-7" for example, or "REPS=R1MACH(4)" 3890C if you have the PORT/SLATEC machine constant routines. 3891C 3892 DOUBLE PRECISION HALF 3893 DOUBLE PRECISION ONE 3894 DOUBLE PRECISION TWO 3895 DOUBLE PRECISION ZERO 3896C 3897 PARAMETER (HALF=0.5) 3898 PARAMETER (ONE=1.0) 3899 PARAMETER (TWO=2.0) 3900 PARAMETER (ZERO=0.0) 3901C 3902 DOUBLE PRECISION EPS 3903 INTEGER I 3904 DOUBLE PRECISION RADD 3905 DOUBLE PRECISION RADDOLD 3906 DOUBLE PRECISION REPS 3907 DOUBLE PRECISION RMANT 3908 DOUBLE PRECISION RMANTOLD 3909 DOUBLE PRECISION TEMP1 3910 DOUBLE PRECISION TEMP2 3911C 3912 RMANT=ONE 3913 RADD=ONE 3914 DO 10 I=1,100 3915 RADDOLD=RADD 3916 RADD=HALF*RADD 3917 RMANTOLD=RMANT 3918 RMANT=RMANT+RADD 3919 IF((RMANT.EQ.RMANTOLD).OR. 3920 * (RMANT-RADD.NE.RMANTOLD))GO TO 20 3921 IF(RMANT.EQ.TWO)THEN 3922 RADDOLD=RADD 3923 GO TO 20 3924 ENDIF 392510 CONTINUE 3926C 3927C Loop doesn't terminate after 100 steps! 3928C We'll just use that value. 3929C 393020 CONTINUE 3931C 3932 RMANT=RMANTOLD 3933 EPS=RADDOLD 3934C 3935C Additional requirement on the value of EPS is that 3936C (1.0+EPS)-EPS.EQ.(1.0-EPS)+EPS.EQ.1.0) 3937C 393830 CONTINUE 3939 TEMP1=ONE+EPS 3940 TEMP2=ONE-EPS 3941 IF(TEMP1-EPS.NE.ONE.OR.TEMP2+EPS.NE.ONE)THEN 3942 EPS=TWO*EPS 3943 GO TO 30 3944 ENDIF 3945C 3946 REPS=EPS 3947 RETURN 3948 END 3949