1 PSCDOC:PCON60.DOC 17 February 1988 PCON60, The University of Pittsburgh Continuation Program 2 Introduction PCON60 is a version of PITCON, the University of Pittsburgh continuation program, rewritten for use on the Cray and other vector machines. Versions for IBM PC and the Apple Macintosh are also available. PITCON was written by Professor Werner C Rheinboldt and John Burkardt, of the Department of Mathematics and Statistics at the University of Pittsburgh, Pittsburgh, Pennsylvania, 15260, USA. Work on this package was partially supported by the National Science Foundation under grants MCS-78-05299 and MCS-83-09926. PITCON computes successive solution points along a one dimensional manifold of a system of nonlinear equations F(X)=0 where there are NVAR variables X and NVAR-1 functions F. In many ways, PITCON solves a nonlinear algebraic system with one degree of freedom much like an initial value ODE solver handles a differential equation. For instance, the user must begin the computation by specifying an approximate initial solution, and subsequent solutions returned by PITCON lie on the curve through the initial point and implicitly defined by F(X)=0. The extra degree of freedom in the system is analogous to the role of time in differential equations, except that at each step, PITCON is free to choose any of of the variables to temporarily play the parameterizing role that time plays for ODE's. This allows PITCON to easily handle singularities such as limit points that would defeat a more limited program. While PITCON formally may be applied to any nonlinear system with one degree of freedom, it is typically used to investigate the equilibrium behavior of physical systems with several degrees of freedom. In any single run of PITCON, all but one of the degrees of freedom are held fixed. Examples of this sort of problem include the behavior of a loaded beam, where the magnitude of the load is varied, the attitude of an aircraft as the pilot alters the controls, or the composition of a chemical system as the temperature is varied. Program options include the ability to report all solutions discovered for which a given component of the solution has a specified value. Unlike differential equations, there may be several solutions for a given value of the 'parameter'. Another option will search for limit points in a given component. Limit points are solution points at which the given component has reached a (local) extreme value. Viewed as part of the curve F(X)=0, the tangent vector is perpendicular to the given coordinate direction. Another feature of the program is the use of two work arrays, IWORK and RWORK. All the information required to carry out further computations is saved in these two vectors. Should the user save these arrays, and restore them later, the interrupted computation can proceed. 2 Calling Sequence SUBROUTINE PITCON(DF,FPAR,FX,IERROR,IPAR,IWORK,LIW, NVAR,RWORK,LRW,XR) 3 DF The name of the routine which computes the Jacobian of the nonlinear function. This name must be declared EXTERNAL in the calling program. DF must evaluate the Jacobian matrix at a given point X and store it in an array in accordance with the format used by the solver. For the simplest case (using DENSLV), DF stores D F(I)/D X(J) into FJAC(I,J). DF is not needed if the finite difference option is used (IWORK(9)=1 or 2). In that case, only a dummy name is required for DF. The array to contain the Jacobian is zeroed out before DF is called, so only nonzero elements need to be stored. DF must have the form: SUBROUTINE DF(NVAR,FPAR,IPAR,X,FJAC,IERROR) DIMENSION FPAR(1),IPAR(1),X(NVAR) DIMENSION FJAC(NVAR,NVAR) NVAR Number of variables. FPAR Vector for passing real user parameters. IPAR Vector for passing integer user parameters. X Point at which Jacobian is desired. FJAC Array in which Jacobian is to be put. If DENSLV is the solver: FJAC must be dimensioned (NVAR,NVAR) as shown above, and DF sets FJAC(I,J)=D F(I)/DX(J). If BANSLV is the solver: Then the main portion of the Jacobian, rows and columns 1 through NVAR-1, is assumed to be a banded matrix in the standard LINPACK form with lower bandwidth ML and upper bandwidth MU. However, the final column of the Jacobian is assumed to be full. BANSLV will pass to DF the beginning of the storage for FJAC, but it is probably best not to doubly dimension FJAC since the first portion of it is a (2*ML+MU+1, NEQN) array, followed by a single column of length NEQN. Thus the simplest approach is to declare FJAC of dimension 1, and then to store values as follows: If J is less than NVAR, then if I-J .LE. ML and J-I .LE. MU, set K=(2*ML+MU)*J + I - ML and set FJAC(K)=D F(I)/DX(J). If J equals NVAR, then set K=(2*ML+MU+1)*(NVAR-1)+I, and set FJAC(K)=D F(I)/DX(J). IERROR Error return flag. DF should set this to 0 for normal return, nonzero if trouble. 3 FPAR A real array provided for the user's convenience. It is passed to DF and FX and may be used by the user to transmit information between his calling program and those subprograms. The dimension of FPAR and its contents are up to the user. Internally, the program declares DIMENSION FPAR(1) but never references it. 3 FX The name of the routine which computes the value of the nonlinear function. This name must be declared EXTERNAL in the calling program. FX should evaluate the NVAR-1 function components at the input point X, and store the result in FVEC. (An augmenting equation will be stored in row NVAR by the program). SUBROUTINE FX(NVAR,FPAR,IPAR,X,FVEC,IERROR) DIMENSION FPAR(1),IPAR(1),X(NVAR),FVEC(NVAR) NVAR Number of variables. FPAR Array of real user parameters. IPAR Array of integer user parameters. X Point at which function evaluation is required. FVEC Value of function at point X. This is to be computed by FX. IERROR FX should set this to 0 if there are no problems, or to 2 if there is a problem. 3 IERROR Program error return flag. On return from the program, a nonzero value of IERROR is a warning of some problem which may be minor, serious, or fatal. 0, No errors occurred. 1, Insufficient storage provided in RWORK and IWORK, or NVAR is less than 2. 2, User defined error condition reported through FX or DF subroutines. 3, A numerically singular matrix was encountered. Continuation cannot proceed without some redefinition of the problem. 4, Unsuccessful corrector iteration. Loosening the tolerances RWORK(1) and RWORK(2), or decreasing the maximum stepsize RWORK(4) might help. 5, Too many corrector steps. The corrector iteration was proceeding properly, but too slowly. Increase number of Newton steps IWORK(17), increase the error tolerances RWORK(1) or RWORK(2), or decrease RWORK(4). 6, Null tangent vector. A serious error which indicates a data problem or singularity in the nonlinear system. 7, Root finder failed while searching for a limit point. 8, Limit point iteration took too many steps. 9, Undiagnosed error condition. 3 IPAR Integer array for the user's convenience in transmitting parameters between the calling program and the user routines FX and DF. IPAR is declared of DIMENSION (1) in the program and otherwise ignored. Note, however, that if BANSLV is used for the solver routine, then IPAR(1) must contain the lower bandwidth, and IPAR(2) the upper bandwidth of the Jacobian matrix (ignoring the final column, which may be full). 3 IWORK Integer work array used inside the package. IWORK is presumed to be dimensioned LIW in the calling program. The specific allocation of IWORK is described in the section devoted to the work arrays. 3 LIW The dimension of IWORK set in the calling program. The minimum acceptable value of LIW depends on the solver chosen, but for either DENSLV or BANSLV this means that LIW=29+NVAR is acceptable. 3 NVAR The number of variables X. This is, of course, one greater than the number of equations. NVAR must be at least 2. 3 RWORK Real work array used in the program. RWORK is dimensioned LRW. The specific allocation of RWORK is described in the section devoted to the work arrays. 3 LRW The dimension of RWORK as set in the calling program. The minimum acceptable value depends heavily on the solver options chosen. For DENSLV with user-supplied Jacobian, LRW=29+NVAR*(NVAR+4). For DENSLV with internally approximated Jacobian, LRW=29+NVAR*(NVAR+6). For BANSLV, with a Jacobian matrix with upper bandwidth MU and lower bandwidth ML, and NBAND=2*ML+MU+1, with user supplied Jacobian, LRW=29+NVAR*(NBAND+6)-NBAND For BANSLV with internally approximated Jacobian, LRW=29+NVAR*(NBAND+9)-NBAND. 3 XR A vector of length NVAR. On the first call, XR should contain a starting point which at least approximately satisfies F(XR)=0. Thereafter, on return from the program with IERROR=0, XR will hold the most recently computed point, whether continuation, target or limit point. 2 Work Arrays Input to the program includes the setting of some of the entries in IWORK and RWORK. Some of this input is optional. The user input section of IWORK is entries 1 through 8, 17 and 29. For RWORK it is entries 1 through 7 and entry 20. IWORK(1) must be set by the user. All other entries have default values. 3 IWORK IWORK(1) On first call only, the user must set IWORK(1)=0. Thereafter, the program sets IWORK(1) before return to explain what kind of point is being returned. This return code is: 1 return of corrected starting point. 2 return of continuation point. 3 return of target point. 4 return of limit point. NOTE: At any time, you may call PCON60 with IWORK(1) negative. This indicates that you wish to check your jacobian routine against a finite difference approximation. The program will call your jacobian routine, and then estimate the jacobian, and print out the value of the maximum difference, and the row and column of the jacobian in which it appears. IWORK(2) The component of the current continuation point (never target or limit point) which is to be used as the continuation parameter. On first call, the program is willing to use index NVAR as a default, but the user should set this value if more is known. After the first call, the program sets this value for each step automatically. Note that a poor choice may cause the algorithm to fail. IWORK(3) Parameterization option. The program would prefer to be free to choose a new continuation index from step to step. The value of IWORK(3) allows or prohibits this action. IWORK(3)=0 allows the program to vary the parameter, IWORK(3)=1 forces the program to use whatever the contents of IWORK(2) are, which will not be changed from the user's input or the default. The default is IWORK(3)=0. IWORK(4) Newton Jacobian update option. If IWORK(4)=0, the Jacobian is reevaluated at every step of the Newton iteration. This is costly, but accurate. If IWORK(4)=1, the Jacobian is evaluated only on the first and IWORK(17)-th steps of the Newton iteration. If IWORK(4)=2, the Jacobian is evaluated only when absolutely necessary: the first step, and when a Newton iteration fails. This option is most suitable for problems with mild nonlinearities, such as discretizations of two point boundary value problems. The default is IWORK(4)=0. IWORK(5) Target point index. If IWORK(5) is not zero, it is presumed to be the component index between 1 and NVAR for which target points are sought. In this case, the value of RWORK(7) is assumed to be the target value. The program will monitor every new continuation point, and if it finds that a target point lies between the new and previous points, it will compute the target point and return. This target point will have the property that the IWORK(5) component of the solution will have the value RWORK(7). For a given problem there may be zero, one or many target points returned. The default of IWORK(5) is 0. IWORK(6) is the limit point index. If it is nonzero, then the program will search for limit points in the specified component of the solution vector, that is, points whose IWORK(6)-th component is locally extremal, or equivalently where the IWORK(6)-th component of the tangent vector is zero. The default for IWORK(6) is zero. IWORK(7) controls the amount of output produced by the program. A value of zero means no output is produced. 3 produces the most output. The default is 1. IWORK(8) is the FORTRAN unit to which output is to be written. The default value is site dependent but should be the standard output device. The default is 6 on the Cray, Vax and PC, or 9 on the Macintosh. IWORK(9) contains the Jacobian option, which tells the program whether the user has supplied a Jacobian routine, or wants the program to approximate the Jacobian. 0, the user has supplied the Jacobian. 1, program is to use forward difference approximation. 2, program is to use central difference approximation. IWORK(9) defaults to 0. IWORK(10) is a state indicator by which the program keeps track of its progress. Values include: 0, start up. Unchecked starting point available. 1, first step. Corrected starting point available. 2, two successive continuation points available, and the tangent at the oldest of these. 3, two successive continuation points available, and the tangent at the newest of these. IWORK(11) is used to help the program avoid computing a target point twice. If a target point has been found, then the target index IWORK(5) is copied into IWORK(11). IWORK(12) contains the second best choice for the local parameterization. This index might be used if the first choice causes poor performance in the Newton correction iteration. IWORK(13) contains the beginning location of unused integer work space, which may be used by the solver. IWORK(14) equals LIW, the user declared dimension of the array IWORK. IWORK(15) is the beginning location in RWORK of the unused real storage which may be used by the linear equation solver. IWORK(16) is LRW, the user declared dimension of RWORK. IWORK(17) contains the maximum number of corrector steps that may be taken in a single Newton iteration in which the Jacobian is updated at every step. If the Jacobian is only evaluated at the beginning of the Newton iteration then 2*IWORK(17) steps may be taken. IWORK(17) must be greater than 0. It defaults to 10. IWORK(18) contains the number of stepsize reductions required before convergence which produced the last continuation point. IWORK(19) counts the total number of calls to the user Jacobian routine DF. IWORK(20) counts the number of calls to factor the matrix. If DENSLV is the solver chosen by IWORK(29) then factorization is done by the Linpack routine SGEFA. Otherwise, SGBFA will be used. IWORK(21) counts the number of calls to the back-substitution routine. If DENSLV is the solver chosen by IWORK(29) then back substitution is done by the Linpack routine SGELS. Otherwise SGBSL will be used. IWORK(22) counts the number of calls to the user function routine FX. IWORK(23) counts the number of steps taken in the limit point iteration. Each step involves setting an approximate limit point and using Newton iteration to correct it. IWORK(24) counts the number of Newton corrector steps used during the computation of target points. IWORK(25) counts the number of Newton steps taken during the correction of a starting point or the continuation points. IWORK(26) counts the number of times the predictor stepsize had to be reduced before the Newton correction would converge. IWORK(27) counts the number of calls to the program. This also corresponds to the number of points computed. IWORK(28) is used to count the number of Newton steps taken during the current iteration. IWORK(29) is user input, the solver chosen. A value of 0 chooses the full storage solver DENSLV, while a value of 1 chooses the banded solver BANSLV. The default is 0. IWORK(30) and on are reserved for use by the linear equation solver, and typically are used for pivoting. 3 RWORK RWORK(1) contains the absolute error tolerance. This value is used mainly during the Newton iteration, where the size of the residual is compared against an error tolerance which is the sum of a relative term plus RWORK(1). RWORK(1) defaults to SQRT(EPMACH) where EPMACH is the machine relative precision. RWORK(2) contains the relative error tolerance. This value is used mainly during the Newton iteration, where the size of the residual is compared against an error tolerance which is the sum of the relative error tolerance times the norm of the current solution plus an absolute error tolerance. RWORK(2) defaults to SQRT(EPMACH) where EPMACH is the machine relative precision. RWORK(3) is the minimum allowable predictor stepsize. If failures of the Newton correction force the stepsize down to this level, then the program will give up. The default value is SQRT(EPMACH). RWORK(4) contains the maximum allowable predictor step. Too generous a value may cause erratic behavior of the program. The default value is SQRT(NVAR). RWORK(5) contains the predictor stepsize. On first call, it should be set by the user. Thereafter it is set by the program. RWORK(5) should be positive. In order to travel in the negative direction, see RWORK(6). The default initial value is 0.5*(RWORK(3)+RWORK(4)). RWORK(6) contains the local continuation direction, which is either +1 or -1. This simply asserts that the program is moving in the direction of increasing or decreasing values of the local continuation variable, whose index is in IWORK(2). On first call, the user must choose IWORK(2). Therefore, by setting RWORK(6), the user may also specify whether the program is to move initially to increase or decrease the variable whose index is IWORK(2). RWORK(6) defaults to +1. RWORK(7) contains the target value. It is only used if a target index has been specified through IWORK(5). In that case, solution points whose IWORK(5) component is equal to RWORK(7) are sought, and the code will return every time it finds such a point. RWORK(7) does not have a default value. The program does not set it, and it is not referenced unless IWORK(5) has been set. RWORK(8) Not currently used. RWORK(9) Not currently used. RWORK(10) contains a minimum angle used in the steplength computation, equal to 2.0*ARCCOS(1-EPMACH). RWORK(11) is an estimate of the angle between the tangent vectors at the last two continuation points. RWORK(12) is an estimate of the pseudoarclength coordinate of the previous continuation point. The pseudoarclength is a rough measure of the distance of the point from the starting point, measured along the curve we are following. Thus each new point should have a larger coordinate, except for target and limit points which lie between the two most recent continuation points. RWORK(13) is an estimate of the pseudoarclength coordinate of the current continuation point. RWORK(14) is an estimate of the pseudoarclength coordinate of the current limit or target point, if any. RWORK(15) records the size of the correction of the most recent continuation point. That is, it is the maximum norm of the difference between the initial predicted value and the final corrected value. RWORK(16) is an estimate of the curvature between the last two continuation points. RWORK(17) is the sign of the determinant of the augmented matrix at the last continuation point whose tangent vector has been calculated. RWORK(18) is not used. RWORK(19) is not used. RWORK(20) is the maximum growth factor for the predictor stepsize based on the previous secant stepsize. The stepsize algorithm will produce a suggested step that is no less that the previous secant step divided by the factor, and no greater than the previous secant step multiplied by that factor. RWORK(20) defaults to 3 RWORK(21) is the (Euclidean) secant distance between the last two computed continuation points. RWORK(22) is the previous value of RWORK(21). RWORK(23) is a judgment of the quality of the Newton corrector convergence at the last continuation point. RWORK(24) is the value of the component of the current tangent vector in the current continuation index. RWORK(25) is the value of the component of the previous tangent vector in the current continuation index. RWORK(26) is the value of the component of the current tangent vector in the IWORK(6) index (limit index). RWORK(27) is the value of the component of the previous tangent vector in the IWORK(6) index (limit index). RWORK(28) is the value of RWORK(7) when the last target point was computed. RWORK(29) is the sign of the determinant of the augmented matrix at the previous continuation point whose tangent vector has been calculated. RWORK(30) through RWORK(30+4*NVAR-1) are used by the program to hold an old and new continuation point, a tangent vector and a work vector. Subsequent entries of RWORK are used by the linear solver. 2 Programming Notes The minimal input and extra routines required to use the program are as follows: Write a function routine FX of the form described above. Use DENSLV (see IWORK(29) ) for the solver. Skip writing a Jacobian routine by using the finite difference option, and instead pass the name of FX as the Jacobian name as well. Declare the name of the function FX as EXTERNAL. For input, set NVAR in accordance with your problem. Then: Dimension the vector IWORK of size LIW=29+NVAR. Dimension RWORK of size LRW=29+NVAR*(NVAR+6). Dimension IPAR(1) and FPAR(1). Dimension XR(NVAR) and set it to an approximate solution of F(XR)=0. Set the following entries of IWORK and RWORK: IWORK(1)=0 (Problem startup) IWORK(7)=3 (Maximum internally generated output) IWORK(9)=1 (Forward difference Jacobian) Now call the program repeatedly, and never change any of its arguments. Check the value of IERROR to decide whether the code is working satisfactorily. Print out the vector XR to see the current solution point. The most obvious input to try to set appropriately after some practice would be the error tolerances RWORK(1) and RWORK(2), the minimum, maximum and initial stepsizes RWORK(3), RWORK(4) and RWORK(5), and the initial continuation index IWORK(2). For speed and efficiency, you're going to have to write your own Jacobian routine. Fortunately, you can check it by comparing results with the finite difference Jacobian. For a particular problem, the target and limit point input can be very useful. Say you have a discretized problem and you want to compare solutions for 6, 11, 21, 41 and 81 points. Clearly, to see if the solutions converge for more points, you want to compare the solutions at the same value of the parameter or load. You can only do so with the target option. Limit points can also be of importance - in structural problems they can represent the loss of stability. 2 Sample Problem What follows is a sample main program, function and Jacobian which have been used with PCON60. PARAMETER (NVAR=3) PARAMETER (LIW=NVAR+29) PARAMETER (LRW=29+(6+NVAR)*NVAR) EXTERNAL FXROTH,FPROTH DIMENSION FPAR(1),IPAR(1),IWORK(LIW),RWORK(LRW),XR(NVAR) LOUNIT=20 OPEN(UNIT=LOUNIT,FILE='PCPRB1.LPT',STATUS='NEW') DO 10 I=1,LIW IWORK(I)=0 10 CONTINUE DO 20 I=1,LRW RWORK(I)=0.0 20 CONTINUE C C IWORK(1)=0 ; This is a startup C IWORK(2)=2 ; Use index 2 for first parameter C IWORK(3)=0 ; Program may choose index C IWORK(4)=0 ; Update Jacobian every newton step C IWORK(5)=3 ; Seek target values for index 3 C IWORK(6)=1 ; Seek limit points in index 1 C IWORK(7)=0 ; No internal output C IWORK(8)=20; Output unit to use C IWORK(9)=0 ; Use user's Jacobian routine C IWORK(29)=0; Full storage solver C IWORK(1)=0 IWORK(2)=2 IWORK(3)=0 IWORK(4)=0 IWORK(5)=3 IWORK(6)=1 IWORK(7)=0 IWORK(8)=LOUNIT IWORK(9)=0 IWORK(29)=0 C C RWORK(1)=0.0001 ; Absolute error tolerance C RWORK(2)=0.0001 ; Relative error tolerance C RWORK(3)=0.01 ; Minimum stepsize C RWORK(4)=20.0 ; Maximum stepsize C RWORK(5)=0.3 ; Starting stepsize C RWORK(6)=1.0 ; Starting direction C RWORK(7)=1.0 ; Target value (Seek solution with X(3)=1) C RWORK(1)=0.0001 RWORK(2)=0.0001 RWORK(3)=0.01 RWORK(4)=20.0 RWORK(5)=0.3 RWORK(6)=1.0 RWORK(7)=1.0 C C SET STARTING POINT C XR(1)=15.0 XR(2)=-2.0 XR(3)=0.0 WRITE(LOUNIT,1000) WRITE(LOUNIT,1010) I=0 WRITE(LOUNIT,1030)I WRITE(LOUNIT,1040)(XR(J),J=1,NVAR) C C TAKE ANOTHER STEP C DO 30 I=1,30 CALL PITCON(FPROTH,FPAR,FXROTH,IERROR,IPAR,IWORK,LIW, * NVAR,RWORK,LRW,XR) IF(IWORK(1).EQ.1)WRITE(LOUNIT,1080)I IF(IWORK(1).EQ.2)WRITE(LOUNIT,1050)I IF(IWORK(1).EQ.3)WRITE(LOUNIT,1060)I IF(IWORK(1).EQ.4)WRITE(LOUNIT,1070)I WRITE(LOUNIT,1040)(XR(J),J=1,NVAR) IF(IWORK(1).EQ.3)GO TO 40 IF(IERROR.NE.0)THEN WRITE(LOUNIT,1090)IERROR GO TO 40 ENDIF 30 CONTINUE 40 CONTINUE CLOSE(UNIT=LOUNIT) STOP 1000 FORMAT('1FREUDENSTEIN-ROTH FUNCTION') 1010 FORMAT(' TEST EXAMPLE FOR CONTINUATION PROGRAM') 1030 FORMAT(' STARTING POINT, STEP',I6) 1040 FORMAT(1X,3G14.6) 1050 FORMAT(' CONTINUATION POINT, STEP',I6) 1060 FORMAT(' TARGET POINT, STEP',I6) 1070 FORMAT(' LIMIT POINT, STEP',I6) 1080 FORMAT(' CORRECTED STARTING POINT, STEP',I6) 1090 FORMAT(' ITERATION HALTED, ERROR FLAG=',I6) END SUBROUTINE FXROTH(NVAR,FPAR,IPAR,X,FX,IERROR) DIMENSION FPAR(1),FX(NVAR),IPAR(1),X(NVAR) FX(1)=X(1)-((X(2)-5.0)*X(2)+2.0)*X(2)-13.0+34.0*(X(3)-1.0) FX(2)=X(1)+((X(2)+1.0)*X(2)-14.0)*X(2)-29.0+10.0*(X(3)-1.0) RETURN END SUBROUTINE FPROTH(NVAR,FPAR,IPAR,X,FJAC,IERROR) DIMENSION FPAR(1),FJAC(NVAR,NVAR),IPAR(1),X(NVAR) FJAC(1,1)=1.0 FJAC(1,2)=(-3.0*X(2)+10.0)*X(2)-2.0 FJAC(1,3)=34.0 FJAC(2,1)=1.0 FJAC(2,2)=(3.0*X(2)+2.0)*X(2)-14.0 FJAC(2,3)=10.0 RETURN END The output from this program would be: FREUDENSTEIN-ROTH FUNCTION TEST EXAMPLE FOR CONTINUATION PROGRAM STARTING POINT, STEP 0 15.0000 -2.00000 0.000000E+00 CORRECTED STARTING POINT, STEP 1 15.0000 -2.00000 0.000000E+00 CONTINUATION POINT, STEP 2 14.7105 -1.94205 0.653814E-01 CONTINUATION POINT, STEP 3 14.2846 -1.72915 0.268742 LIMIT POINT, STEP 4 14.2831 -1.74137 0.258583 CONTINUATION POINT, STEP 5 16.9061 -1.20941 0.546845 CONTINUATION POINT, STEP 6 24.9179 -0.599064 0.555136 CONTINUATION POINT, STEP 7 44.8783 0.487669 0.595272E-01 CONTINUATION POINT, STEP 8 60.0889 1.57585 -0.542365 CONTINUATION POINT, STEP 9 -11.1747 4.23515 1.55665 TARGET POINT, STEP 10 5.00000 4.00000 1.00000 2 VMS Usage The PCON60 document is available in the file PSCDOC:PCON60.DOC Online help is available by typing "HELP PCON60" To compile your main calling program, link it to PCON60 and run it, use commands like the following: $ FORTRAN MYPROG $ LINK MYPROG, PSC$LIB:PCON60/LIB $ RUN MYPROG 2 COS Usage The PCON60 document is available in the directory PSCDOC. To attach it and make a copy which may be moved to your own account or printed, use commands like: ATTACH,DN=PCON60,ID=PSCDOC. COPYD,I=PCON60,O=MYCOPY. DISPOSE,DN=MYCOPY,DC=ST,TEXT='PCON60.DOC'. Similarly, the source code for PCON60 is available in the directory PSCSRC, and could be copied by commands like: ATTACH,DN=PCON60,ID=PSCSRC. COPYD,I=PCON60,O=MYCOPY. DISPOSE,DN=MYCOPY,DC=ST,TEXT='PCON60.FOR'. To compile your program, link it to PCON60 and run it under COS, use commands like FETCH,DN=MYPROG,TEXT='MYPROG.FOR'. CFT,L=0,I=MYPROG,ON=AS. ATTACH,DN=PCON60,ID=PUBLIC. LDR,MAP=OFF,SET=ZERO,LIB=PCON60. 2 UNICOS Usage The PCON60 document is available in the directory /usr/local/doc/pcon60.doc To compile your program, link it with PCON60 and run it under UNICOS, use commands like cft myprog segldr myprog.o /usr/local/lib/libpcon60.o mv a.out myprog myprog 2 Personal Computer Versions PCON60 is available for the IBM PC and compatible machines as well as the Apple Macintosh. Information concerning these implementations may be obtained by contacting John Burkardt at: Pittsburgh Supercomputing Center Mellon Institute 4400 Fifth Avenue Pittsburgh, Pennsylvania 15213 2 References 1. Werner Rheinboldt, Solution Field of Nonlinear Equations and Continuation Methods, SIAM Journal of Numerical Analysis, 17, 1980, pages 221-237. 2. Cor den Heijer and Werner Rheinboldt, On Steplength Algorithms for a Class of Continuation Methods, SIAM Journal of Numerical Analysis, 18, 1981, pages 925-947. 3. Werner Rheinboldt, Numerical Analysis of Continuation Methods for Nonlinear Structural Problems, Computers and Structures, 13, 1981, pages 103-114. 4. Werner Rheinboldt and John Burkardt, A Locally Parameterized Continuation Process, ACM Transactions on Mathematical Software, Volume 9, Number 2, June 1983, pages 215-235. 5. Werner Rheinboldt and John Burkardt, Algorithm 596, A Program for a Locally Parameterized Continuation Process, ACM Transactions on Mathematical Software, Volume 9, Number 2, June 1983, Pages 236-241. 6. J J Dongarra, J R Bunch, C B Moler and G W Stewart, LINPACK User's Guide, Society for Industrial and Applied Mathematics, Philadelphia, 1979. 7. Richard Brent, Algorithms for Minimization without Derivatives, Prentice Hall, 1973. 8. Tony Chan, Deflated Decomposition of Solutions of Nearly Singular Systems, Technical Report 225, Computer Science Department, Yale University, New Haven, Connecticut, 06520, 1982.