c  NEWTON.F  27 July 1992
c
      program main

c*********************************************************************72
c
cc MAIN is the main program for NEWTON.
c
c  Discussion:
c
c    NEWTON is A PROGRAM WHICH CARRIES OUT THE NEWTON ITERATION METHOD ON
c    A SYSTEM OF UP TO 10 EQUATIONS.
c
c
c  NEWTON CAN SOLVE TWO DISTINCT BUT RELATED GROUPS OF PROBLEMS.
c
c  THE FIRST TYPE OF PROBLEM INVOLVES N EQUATIONS IN N UNKNOWNS.
c  YOU SUPPLY A SINGLE APPROXIMATE GUESS, AND NEWTON ITERATION
c  IS USED TO FIND THE SOLUTION.
c
c  THE SECOND TYPE OF PROBLEM INVOLVES N-1 EQUATIONS IN N UNKNOWNS.
c  THIS IS KNOWN AS THE 'CONTINUATION PROCESS'.
c  IN THIS CASE, THE UNDERDETERMINED SYSTEM NOW DEFINES A
c  CURVE OF POINTS.  YOU SUPPLY AN APPROXIMATE POINT ON THE CURVE,
c  AND THEN USE NEWTON ITERATION TO GET ONTO THE CURVE.  NOW
c  YOU HAVE THE PROGRAM TAKE A STEP IN THE TANGENT DIRECTION
c  TO FIND ANOTHER APPROXIMATE POINT.  NEWTON ITERATION IS PERFORMED
c  ON THIS POINT TO GET BACK ON THE CURVE.  BY REPEATED APPLICATIONS
c  OF THIS PROCESS, YOU CAN FIND A SEQUENCE OF POINTS ALONG THE
c  CURVE.  NOTE THAT IN THIS CASE, THE PROGRAM ADDS AN EXTRA
c  EQUATION, SO THAT THE 'NUMBER OF EQUATIONS' MUST BE INCREASED
c  BY 1.
c
c  ALSO, BECAUSE OF THE CRUDE NATURE OF THIS ALGORITHM, THERE
c  MAY BE SOME PROBLEMS.  THERE IS NO STEPSIZE ALGORITHM IN THE CODE.
c  THAT MEANS YOU MUST CHOOSE A STEPSIZE H TO USE, WHICH MAY BE
c  TOO LARGE, CAUSING NEWTON TO FAIL.  ALSO, THE FORM OF THE EXTRA
c  EQUATION MAY NOT WORK WELL WITH SOME PROBLEMS.  AGAIN, YOU
c  WILL MOST LIKELY NOTICE THIS WHEN NEWTON DOES NOT CONVERGE WELL.
c
c
c  TO USE NEWTON, YOU MUST FIRST WRITE TWO SUBROUTINES WHICH
c  DEFINE YOUR PROBLEM.  THE FIRST OF THESE MUST BE NAMED
c  FUNK AND MUST HAVE THE FOLLOWING FORM
c
c     SUBROUTINE FUNK(FX,NEQN,X)
c     DIMENSION FX(NEQN)
c     DIMENSION X(NEQN)
c
c  FROM TIME TO TIME, THE PROGRAM WILL CALL THIS SUBROUTINE,
c  AND WILL GIVE YOU THE VALUE OF X.  IT EXPECTS YOU TO
c  ASSIGN VALUES TO FX AND RETURN.  IF YOU ARE USING THE CONTINUATION
c  OPTION, DO NOT BOTHER ASSIGNING A VALUE TO FX(NEQN).  THAT IS
c  DONE FOR YOU.
c
c  THE SECOND SUBROUTINE EVALUATES THE JACOBIAN MATRIX
c  OF PARTIAL DERIVATIVES OF THE FUNCTION DEFINED IN FUNK.
c  THE SUBROUTINE MUST HAVE THE FOLLOWING FORM
c
c     SUBROUTINE JACOB(FPRIME,NROW,NEQN,X)
c     DIMENSION FPRIME(NROW,NEQN)
c     DIMENSION X(NEQN)
c
c  FROM TIME TO TIME, THE PROGRAM WILL CALL JACOB, AND SUPPLY
c  THE VALUE OF X.  YOU MUST EVALUATE THE ENTRIES OF THE JACOBIAN
c  WHERE FPRIME(I,J)=D F(I)/D X(J).  NOTE THAT IF YOU DECIDE TO
c  USE THE FINITE DIFFERENCE JACOBIAN OPTION, YOU DO NOT HAVE
c  TO PERFORM ANY CALCULATIONS IN JACOB, BUT YOU STILL MUST
c  SUPPLY A SUBROUTINE OF THAT NAME.  IF YOU ARE USING THE CONTINUATION
c  OPTION, DO NOT BOTHER ASSIGNING VALUES TO THE LAST ROW OF
c  FPRIME.  THAT WILL BE DONE FOR YOU.
c
c
c  RUNNING NEWTON
c
c  ONCE YOU HAVE WRITTEN YOUR SUBROUTINES, AND SUPPOSING
c  THAT THEY ARE BOTH STORED IN ONE FILE NAMED 'MYPROG.FOR',
c  YOU CAN RUN NEWTON BY TYPING THE FOLLOWING COMMAND
c
c  EXECUTE MYPROG.FOR, NEWTON.FOR[??????,??????]
c
c  THE PROGRAM WILL FIRST ASK YOU WHETHER YOU WANT TO USE
c  THE STANDARD METHOD, OR THE CONTINUATION OPTION.
c  THEN IT WILL ASK YOU FOR THE NUMBER OF EQUATIONS
c  IN YOUR SYSTEM.  THEN IT WILL WANT THE CURRENT STARTING X.
c  AFTER THAT, THE PROGRAM WILL ASK YOU
c  AGAIN AND AGAIN WHAT THE NEXT COMMAND IS, AND WILL CARRY IT
c  OUT.  THE FOLLOWING COMMANDS ARE AVAILABLE TO YOU -
c
c  B = BEGIN NEWTON ITERATION
c  C = CONTINUE CURRENT NEWTON ITERATION
c  D = SET DAMPING OPTION
c  E = SET NUMBER OF EQUATIONS
c      (ABORTS CURRENT NEWTON ITERATION)
c  F = SET FINITE DIFFERENCE OPTION
c  G = CHECK JACOBIAN VERSUS FINITE DIFFERENCE JACOBIAN
c      (ABORTS CURRENT NEWTON ITERATION)
c  H = SET CONTINUATION STEPSIZES
c  I = SET STARTING VALUE OF PREDICTED NEXT CONTINUATION POINT.
c  J = PRINT JACOBIAN
c  K = DEFINE HOW OFTEN JACOBIAN IS EVALUATED
c  L = LIST COMMANDS
c  M = SET MAXIMUM NUMBER OF NEWTON ITERATIONS
c  N = SET THE NORM TO USE
c  O = CHANGE AMOUNT OF OUTPUT
c  P = PRINT X, F(X) (AND TAN(X) FOR CONTINUATION PROBLEMS)
c  S = STOP
c  T = SET ERROR TOLERANCES
c  V = PRINT VALUES OF CONTROL QUANTITIES
c  W = PRINT WORK STATISTICS
c  X = SET CURRENT X APPROXIMATION
c  Y = SET TARGET FOR CONTINUATION
c      (ABORTS CURRENT NEWTON ITERATION)
c
c  THE MEANING OF THESE COMMANDS IS
c
c
c  STARTUP COMMANDS
c
c
c  E = MEANS THAT YOU ARE DECLARING THE NUMBER OF EQUATIONS.
c  SINCE THIS QUESTION IS ASKED AT THE BEGINNING, YOU PROBABLY
c  SHOULD NOT NEED TO USE THIS COMMAND THEREAFTER.  A CLEVER
c  PROGRAMMER, HOWEVER, COULD RUN SEVERAL PROBLEMS IN ONE
c  RUN OF DIFFERENT SIZES USING THE SIZE AS A FLAG IN
c  THE FUNCTION AND JACOBIAN ROUTINES.
c
c  X = MEANS THAT YOU INPUT A POINT X FROM WHICH ALL
c  NEWTON ITERATIONS WILL BEGIN.  THE IMPROVED X PRODUCED
c  BY NEWTON ITERATION CAN BE USED TO PROCEED VIA THE 'C'
c  COMMAND, BUT ANY 'B' COMMAND CAUSES THE PROGRAM TO
c  BEGIN FROM THIS REFERENCE STARTING POINT.  YOU CAN CHANGE
c  X AT ANY TIME.
c
c  Y = MEANS THAT IN YOUR CONTINUATION PROCESS, YOU WOULD LIKE
c  THE PROGRAM TO SEARCH FOR A SPECIAL POINT.  THIS SPECIAL POINT
c  IS DEFINED BY THE VALUE OF A GIVEN INDEX.  IF YOU WERE TRACING
c  A CIRCLE OF RADIUS 9, YOU MIGHT WANT TO FIND THE POINT
c  (3,0).  TO DO THIS, YOU WOULD REQUEST THE POINT WITH INDEX 1
c  EQUAL TO 3, OR WITH INDEX 2 EQUAL TO 0.  NOTE HOWEVER
c  THAT THE PROGRAM ONLY CHECKS FOR THESE POINTS BY COMPARING
c  THE TWO MOST RECENT POINTS IT COMPUTED AND SEEING IF THE
c  GIVEN VALUE LIES BETWEEN THEM.  THUS, ON A CIRCLE, THE
c  TWO POINTS MIGHT BE (2.98, -0.1) AND (2.97, +0.2).  IN THAT
c  CASE, THE PROGRAM WILL NOTICE THAT THE POINT WITH X(2)=0
c  LIES BETWEEN THESE TWO POINTS, BUT IT WOULD NOT REALIZE
c  THAT A POINT WITH X(1)=3 LIES BETWEEN THEM.  SO BE CAREFUL
c  HOW YOU SPECIFY THIS TARGET POINT.
c
c
c  PRINT OUT COMMANDS
c
c
c  P = CAUSES THE PROGRAM TO PRINT THE VALUE OF THE
c  STARTING POINT AND THE FUNCTION F(X).
c
c  J = CAUSES THE PROGRAM TO PRINT THE VALUE OF THE JACOBIAN
c  MATRIX AT THE CURRENT POINT.  IF NO NEWTON STEPS HAVE
c  BEEN TAKEN, THIS IS THE STARTING POINT.  IF NEWTON STEPS
c  HAVE BEEN TAKEN, THE JACOBIAN IS DISPLAYED FOR THE
c  POINT OF LAST EVALUATION, WHICH IS THE POINT JUST
c  BEFORE THE CURRENT ITERATE, UNLESS THE MODIFIED NEWTON
c  METHOD HAS BEEN USED.
c
c
c  COMMANDS FOR THE CONTINUATION STEP
c
c
c  H = SET CONTINUATION STEPSIZE.  IF THE PROGRAM
c  HAS A POINT X, IT WILL COMPUTE THE TANGENT TO THE CURVE
c  AT X, TAN.  TO GET THE NEXT APPROXIMATE POINT ON THE
c  CURVE, IT WILL SET XN NEW = X + H*TAN.  THUS, FOR SMALL H,
c  THE NEW POINT SHOULD BE CLOSE TO THE CURVE, AND
c  FOR LARGE H, PROBABLY FAR OFF.  SINCE THIS POINT IS
c  THE STARTING POINT OF A NEWTON ITERATION, FAR AWAY POINTS
c  CAUSED BY LARGE H VALUES MAY CAUSE NEWTON TO FAIL.  YOU
c  MAY SET H NEGATIVE TO REVERSE DIRECTION.
c
c  I = TAKE CONTINUATION STEP.  THIS COMMAND CAUSES THE PROGRAM
c  TO SET THE STARTING POINT XN NEW = X + H*TANGENT.  IT ALSO
c  MEANS THAT THE PROGRAM WILL INFORM YOU OF THE EXTRA EQUATION,
c  WHICH WILL HAVE THE FORM FX(N)=X(IHOLD)-XHOLD.  THE PROGRAM
c  PICKS VALUES FOR IHOLD AND XHOLD EACH TIME.  IHOLD IS CALLED
c  THE CONTINUATION PARAMETER.  THE EFFECT OF THIS EQUATION
c  IS TO DEMAND THAT THE NEWTON ITERATES ALL FIX THE IHOLD-TH
c  INDEX OF X.  ONCE YOU HAVE ISSUED THE I COMMAND, IT IS UP
c  TO YOU TO USE THE B AND C COMMANDS TO GET A POINT
c  ON THE CURVE.  THUS A CONTINUATION PROCESS WILL INVOLVE
c  USING THE COMMANDS B, I, B, I, B  ETC.  OF COURSE YOU CAN
c  INSERT OTHER OUTPUT COMMANDS OR CHANGE VARIABLES AS YOU
c  GO ALONG.
c
c
c  COMMANDS THAT DEFINE THE NEWTON ITERATION
c
c
c  D = SET THE DAMPING OPTION.  NORMALLY, NEWTON TRIES
c  TO SET XN NEW=XN OLD+CORRECTION.  SOMETIMES, XN NEW IS NOT AN
c  IMPROVEMENT ON XN OLD.  IN THESE CASES, IF DAMPING IS NOT
c  USED, THEN THE PROGRAM IMMEDIATELY RETURNS TO YOU.  HOWEVER,
c  IT IS A FACT THAT IF 'DAMP' IS A SMALL ENOUGH NUMBER,
c  XN NEW=XN OLD+DAMP*CORRECTION WILL BE AN IMPROVEMENT ON XN OLD.
c  (THIS ASSUMES THAT WE ARE USING THE STANDARD NEWTON METHOD,
c  AND ARE NOT LAGGING THE JACOBIAN).
c  IF YOU SPECIFY DAMPING, THEN THE PROGRAM WILL ATTEMPT TO
c  SAVE THE DAY BY SEARCHING FOR A SMALL ENOUGH DAMP
c  TO REDUCE THE FUNCTION NORM OF XN NEW.  IN FACT, IT SETS
c  DAMP TO 1 UNTIL A FAILURE OCCURS, AND THEN REPEATEDLY
c  HALVES DAMP UNTIL THE POINT IS ACCEPTABLE.  WE ACTUALLY
c  ONLY ALLOW 5 HALVINGS BEFORE GIVING UP.
c
c  F = SET THE FINITE DIFFERENCE OPTION.  THE PROGRAM ASSUMES
c  THAT YOU HAVE ACTUALLY WRITTEN FORMULAS FOR THE JACOBIAN
c  IN YOUR JACOB SUBROUTINE.  IF YOU HAVEN'T, YOU CAN ASK
c  THE PROGRAM TO APPROXIMATE THE JACOBIAN USING FINITE
c  DIFFERENCES.  EVEN IF YOU HAVE, YOU CAN USE THIS OPTION
c  TO COMPARE THE RESULTS OF EXACT AND APPROXIMATE JACOBIANS.
c  THE FINITE DIFFERENCE OPTIONS INCLUDE BACKWARD, CENTERED,
c  AND FORWARD DIFFERENCES TO APPROXIMATE D F(I)/D X(J).
c
c  K = DEFINE HOW OFTEN THE JACOBIAN IS TO BE EVALUATED.
c  ALTHOUGH NEWTON'S METHOD GENERALLY EVALUATES THE JACOBIAN
c  AT EVERY STEP, IT IS CHEAPER (ESPECIALLY FOR LARGE PROBLEMS)
c  TO EVALUATE THE JACOBIAN ONLY EVERY SO OFTEN.  EVERY SO
c  OFTEN IS DEFINED TO BE 1 FOR THE ORIGINAL NEWTON'S METHOD,
c  BUT YOU CAN RESET THAT TO A HIGHER VALUE TO SEE WHAT
c  HAPPENS.
c
c  M = SET MAXIMUM NUMBER OF STEPS BEFORE RETURN.  TO KEEP
c  NEWTON FROM RUNNING WILD, THE PROGRAM ALLOWS IT TO TAKE
c  NO MORE THAN M STEPS BEFORE IT ASKS YOU FOR A NEW
c  COMMAND, WHICH MIGHT BE SIMPLY 'C' FOR CONTINUE IF YOU'RE
c  HAPPY WITH THE PROGRESS OF THE ITERATION.  RIGHT NOW,
c  THE NUMBER OF STEPS IS 10.  CHANGE THIS VALUE FOR
c  YOUR OWN TASTE, BUT DON'T MAKE IT TOO LARGE.
c
c  N = SET THE NORM.  THE ERROR TESTS FOR NEWTON'S METHOD
c  ARE BASED ON THE NORMS OF THE FUNCTION, THE POINT X AND
c  THE STEP SIZE (XN NEW-XN OLD).  THERE ARE THREE NORMS
c  AVAILABLE, THE MAXIMUM NORM, THE SUM OF THE ABSOLUTE VALUES,
c  AND THE SQUARE ROOT OF THE SUM OF SQUARES.
c
c  O = DEFINE HOW MUCH OUTPUT YOU WANT IN NEWTON ITERATION.
c  0 MEANS THAT THE NEWTON ITERATION WILL ONLY REPORT WHETHER
c  IT FAILED, SUCCEEDED, OR HASN'T REACHED A CONCLUSION YET.
c  1 MEANS THAT THE NORMS OF X, F(X) AND THE STEP
c  WILL BE PRINTED FOR EVERY STEP.
c  2 MEANS THAT THE ACTUAL VECTORS X, F(X) AND THE STEP WILL
c  BE PRINTED FOR EVERY STEP.
c
c  T = SET ERROR TOLERANCES.  YOU CAN TIGHTEN OR LOOSEN
c  THE ERROR TESTS THIS WAY.  TWO NUMBERS, ABSERR AND
c  RELERR, ARE USED.  THE ERROR TESTS ARE OF THE FOLLOWING FORM
c  WHERE XNORM1, FNORM1 AND SNORM1 ARE THE NORMS OF THE
c  CURRENT X, FX AND STEP, AND XNORM0, FNORM0 AND SNORM0 ARE
c  THE SAME FOR THE PREVIOUS STEP.
c
c  ACCEPTANCE TESTS-
c
c  STEP 1: ACCEPT IF FNORM1.LE.ABSERR
c  STEP N: ACCEPT IF FNORM1.LE.ABSERR
c                    AND SNORM1.LE.RELERR*(XNORM1+ABSERR)
c
c  REJECTION TESTS
c
c  STEP 1: REJECT IF FNORM1.GT.FNORM0+ABSERR
c  STEP N: REJECT IF FNORM1.GT.FNORM0+ABSERR
c          OR     IF SNORM1.GT.SNORM0+ABSERR
c
c
c  EXECUTION COMMANDS
c
c
c  B = BEGIN NEWTON ITERATION.  YOU MUST GIVE THIS COMMAND
c  AT LEAST ONCE BEFORE YOU CAN GIVE THE 'C' COMMAND.  THE
c  NEWTON ITERATION BEGINS FROM THE STARTING POINT YOU GAVE
c  AND TAKES UP TO THE MAXIMUM NUMBER OF STEPS.
c
c  C = CONTINUE NEWTON ITERATION.  NOT TO BE GIVEN UNLESS
c  THE B COMMAND HAS BEEN GIVEN AT LEAST ONCE.  THE NEWTON
c  ITERATION PICKS UP FROM THE LAST ITERATE OF THE ITERATION
c  THAT HAD PAUSED.
c
c
c  OTHER COMMANDS
c
c
c  G = CHECK JACOBIAN VERSUS FINITE DIFFERENCE JACOBIAN.
c  YOU SHOULD ONLY ISSUE THIS COMMAND WHEN YOU ARE NOT IN THE
c  MIDDLE OF A NEWTON ITERATION.  THE CODE WILL COMPARE YOUR
c  JACOBIAN WITH THE FINITE DIFFERENCE JACOBIAN.  YOU CAN
c  SEE THE ERRORS CAUSED BY DISCRETIZATION, AND YOU MAY ALSO
c  BE ABLE TO CATCH MISTAKES IN YOUR JACOBIAN ROUTINE.
c
c  L = LIST THE COMMANDS.  A SHORT LIST OF THE COMMANDS
c  IS GIVEN, IN CASE YOU FORGOT.
c
c  V = PRINT CONTROL QUANTITIES.  THE PROGRAM WILL PRINT
c  THE ERROR TOLERANCES, THE DAMPING OPTION, THE JACOBIAN
c  OPTION, THE NORM USED, THE MAXIMUM NUMBER OF NEWTON
c  STEPS ALLOWED, AND THE NUMBER OF NEWTON STEPS TAKEN BEFORE
c  THE JACOBIAN IS RE-EVALUATED.
c  IF THE CONTINUATION OPTION IS BEING USED, YOU WILL ALSO SEE
c  THE CURRENT STEPSIZE AND THE CURRENT PARAMETER.
c
c  W = PRINT WORK STATISTICS.  YOU CAN SEE HOW MANY TIMES
c  THE FUNCTION AND JACOBIAN WERE EVALUATED, HOW MANY TIMES
c  THE JACOBIAN WAS FACTORED AND SOLVED.
c
c  S  = STOP.
c
      INTEGER   MAXEQN
      PARAMETER (MAXEQN=10)

      CHARACTER COMMAN*1
      REAL      DELX(MAXEQN)
      REAL      FPFACT(MAXEQN,MAXEQN)
      REAL      FPRIME(MAXEQN,MAXEQN)
      REAL      FX(MAXEQN)
      REAL      FXBACK(MAXEQN)
      REAL      FXPLUS(MAXEQN)
      INTEGER   IPIVOT(MAXEQN)
      REAL      TAN(MAXEQN)
      REAL      TANOLD(MAXEQN)
      REAL      TANTMP(MAXEQN)
      REAL      XCNEW(MAXEQN)
      REAL      XNNEW(MAXEQN)
      REAL      XCOLD(MAXEQN)
      REAL      XDELT(MAXEQN)
      REAL      XNOLD(MAXEQN)
      REAL      XNPRED(MAXEQN)

      CALL SETCOM(ABSERR,DAMP,DELX,EPMACH,EPSQRT,FNORM0,FNORM1,
     1 FPFACT,FPRIME,FX,HSTEP,ICON,IDAMP,IDIF,IERROR,IHATE,
     2 IHOLD,IJAC,ILIKE,INORM,IPIVOT,IPRINT,ITARG,MAXEQN,MAXSTP,MODNEW,
     3 NEQN,NFACTS,NFPVAL,NFXVAL,NROW,NSOLVE,NSTP,RELERR,SNORM0,
     4 SNORM1,TAN,TARVAL,XNNEW,XHOLD,XCOLD,XNORM0,XNORM1,XNOLD,XNPRED)
c
c  Find out whether this is a continuation problem.
c
      CALL SETCON(ICON)
c
c  Get the number of equations
c
      CALL SETNEQ(ICON,MAXEQN,NEQN)
      CALL SETX(ICON,IHOLD,NCON,NEQN,NSTP,XCNEW,XCOLD,XHOLD,
     1 XNNEW,XNOLD,XNPRED)
      CALL SETPRM(EPMACH,EPSQRT,FPFACT,FPRIME,FXBACK,FXPLUS,
     & ICON,IDIF,IERROR,IHOLD,IJAC,IPIVOT,MAXEQN,MAXSTP,NEQN,
     & NFACTS,NFPVAL,NFXVAL,NROW,NSTP,TAN,XDELT,XNNEW,XHOLD)
c
c  Get the next command
c
   10 CONTINUE
      WRITE(*,*)' '
      write(*,1180)
      READ(*,1190)COMMAN
      CALL CHRCAP(COMMAN,1,1)
c
c  B = BEGIN NEWTON ITERATION OR CONTINUATION PROCESS
c
      IF(COMMAN.EQ.'B' ) then
        NSTP=0
        DO 30 I=1,NEQN
          XNNEW(I)=XNPRED(I)
   30     CONTINUE
        CALL NEWCON(ABSERR,DAMP,DELX,EPMACH,EPSQRT,FNORM0,FNORM1,
     1  FPFACT,FPRIME,FX,FXBACK,FXPLUS,HSTEP,ICON,IDAMP,IDIF,
     &  IERROR,IHATE,IHOLD,IJAC,ILIKE,INORM,IPIVOT,IPRINT,MAXEQN,
     &  MAXSTP,MODNEW,NEQN,NFACTS,NFPVAL,NFXVAL,NROW,NSOLVE,NSTP,
     &  RELERR,SNORM0,SNORM1,TAN,XNNEW,XHOLD,XCOLD,XDELT,XNORM0,
     &  XNORM1,XNOLD,XNPRED)
        IF(ILIKE.EQ.1 ) then
	  NSTP=0
          CALL SCOPY(NEQN,XNNEW,1,XCNEW,1)
          IF(ICON.EQ.1)write(*,1701)
	  end if
1701  FORMAT(' YOU MAY TAKE ANOTHER CONTINUATION STEP')
        IERROR=0
c
c  C = CONTINUE NEWTON ITERATION
c
      else if ( COMMAN.EQ.'C' ) then
        IF(NSTP.EQ.0 ) then
          write(*,1000)
          write(*,1010)
          GO TO 10
          end if
        CALL NEWCON(ABSERR,DAMP,DELX,EPMACH,EPSQRT,FNORM0,FNORM1,
     1  FPFACT,FPRIME,FX,HSTEP,ICON,IDAMP,IDIF,IERROR,IHATE,
     2  IHOLD,IJAC,ILIKE,INORM,IPIVOT,IPRINT,MAXEQN,MAXSTP,MODNEW,
     3  NEQN,NFACTS,NFPVAL,NFXVAL,NROW,NSOLVE,NSTP,RELERR,SNORM0,
     4  SNORM1,TAN,XNNEW,XHOLD,XCOLD,XNORM0,XNORM1,XNOLD,XNPRED)
        IF(ILIKE.EQ.1)NSTP=0
        IF(ILIKE.EQ.1)CALL SCOPY(NEQN,XNNEW,1,XCNEW,1)
        IF(ILIKE.EQ.1.AND.ICON.EQ.1)write(*,1701)
        IERROR=0
c
c  D = SET DAMPING OPTION
c
      else if ( COMMAN.EQ.'D' ) then
        WRITE(*,*)' '
        write(*,1110)
        write(*,1120)
        READ(*,*)IDAMP
c
c  E = SET NUMBER OF EQUATIONS
c
      else if ( COMMAN.EQ.'E' ) then
        IF(NSTP.GT.0)write(*,1020)
        NSTP=0
        CALL SETNEQ(ICON,MAXEQN,NEQN)
c
c  F = SET OR CHANGE FINITE DIFFERENCE OPTION
c
      else if ( COMMAN.EQ.'F' ) then
        CALL SETDIF(IDIF,IJAC)
c
c  G = CHECK FINITE DIFFERENCE JACOBIAN AND USER JACOBIAN
c
      else if ( COMMAN.EQ.'G' ) then
        IF(NSTP.NE.0)write(*,1020)
        NSTP=0
        CALL JACOB(FPFACT,NROW,NEQN,XNNEW)
        NFPVAL=NFPVAL+1
        CALL DIFJAC(EPMACH,EPSQRT,FPRIME,FXBACK,FXPLUS,
     &  ICON,IDIF,IHOLD,NEQN,NFXVAL,NROW,XDELT,XNNEW,XHOLD)
        CALL CHKJAC(EPSQRT,FPFACT,FPRIME,NROW,NEQN)
c
c  H = SET CONTINUATION STEPSIZES
c
      else if ( COMMAN.EQ.'H' ) then
        CALL SETSTP(HSTEP)
        IF(NSTP.EQ.0)GO TO 10
c       write(*,1030)
c       READ(5,1190)COMMAN
c       CALL CHRCAP(COMMAN,1,1)
c       IF(COMMAN.EQ.'Y')GO TO 100
c
c  I = TAKE CONTINUATION STEP
c
      else if ( COMMAN.EQ.'I' ) then
        CALL PREDIC(ABSERR,EPMACH,EPSQRT,FPFACT,FPRIME,FX,FXBACK,
     &  FXPLUS,HSTEP,ICON,IDIF,IERROR,IHOLD,IJAC,INORM,IPIVOT,ITARG,
     &  MAXEQN,MAXSTP,NEQN,NFACTS,NFPVAL,NFXVAL,NROW,NSOLVE,NSTP,TAN,
     &  TANOLD,TARVAL,XCNEW,XDELT,XHOLD,XCOLD,XNPRED)
c
c  J = PRINT JACOBIAN
c
      else if ( COMMAN.EQ.'J' ) then
        IF(NSTP.GT.0)write(*,1080)
        IF(NSTP.EQ.0)
     1  CALL GETJAC(EPMACH,EPSQRT,FPFACT,FPRIME,ICON,IDIF,IHOLD,
     2  IJAC,IPIVOT,NEQN,NFPVAL,NFXVAL,NROW,XNNEW,XHOLD)
        CALL PRIJAC(FPRIME,NROW,NEQN)
c
c  K = DEFINE FREQUENCY OF JACOBIAN EVALUATION
c
      else if ( COMMAN.EQ.'K' ) then
        write(*,1130)
        write(*,1140)
        READ(*,*)MODNEW
c
c  L = LIST COMMANDS
c
      else if ( COMMAN.EQ.'L' ) then
        WRITE(*,*)' '
        write(*,1210)
        IF(NSTP.GT.0)write(*,1220)
        write(*,1230)
        write(*,1240)
        write(*,1250)
        write(*,1260)
        IF(ICON.EQ.1)write(*,1270)
        IF(ICON.EQ.1)write(*,1280)
        write(*,1290)
        write(*,1300)
        write(*,1310)
        write(*,1320)
        write(*,1330)
        write(*,1340)
        write(*,1350)
        IF(ICON.EQ.1)write(*,1360)
        write(*,1370)
        write(*,1380)
        write(*,1390)
        write(*,1400)
        write(*,1410)
        IF(ICON.EQ.1)write(*,1420)
c
c  M = SET MAXIMUM NUMBER OF STEPS
c
      else if ( COMMAN.EQ.'M' ) then
        WRITE(*,*)' '
        write(*,1090)
        write(*,1100)
        READ(*,*)MAXSTP
        IF(MAXSTP.LT.1)MAXSTP=1
c
c  N = SET NORM
c
      else if ( COMMAN.EQ.'N' ) then
        CALL SETNRM(INORM)
c
c  O = SET AMOUNT OF OUTPUT
c
      else if ( COMMAN.EQ.'O' ) then
        WRITE(*,*)' '
        write(*,1150)
        write(*,1160)
        write(*,1170)
        READ(*,*)IPRINT
c
c  P = PRINT X, FX, (TANX)
c
      else if ( COMMAN.EQ.'P' ) then
        CALL PRIXFX(EPMACH,EPSQRT,FPFACT,FPRIME,FX,FXBACK,FXPLUS,
     &  ICON,IDIF,IERROR,IHOLD,IJAC,IPIVOT,MAXEQN,MAXSTP,NEQN,
     &  NFACTS,NFPVAL,NFXVAL,NROW,NSOLVE,NSTP,TAN,TANTMP,XDELT,
     &  XNNEW,XHOLD)
c
c  Q = CHANGE PARAMETERS
c
      else if ( COMMAN.EQ.'Q' ) then
        write(*,1040)IHOLD
        write(*,1050)NEQN
        READ(*,*)ITEMP
        IF(ITEMP.GE.1.AND.ITEMP.LE.NEQN)IHOLD=ITEMP
        XHOLD=XNNEW(IHOLD)
c
c  S = STOP
c
      else if ( COMMAN.EQ.'S' ) then
        WRITE(*,*)'NEWTON has been requested to stop.'
        STOP
c
c  T = ENTER ERROR TOLERANCES
c
      else if ( COMMAN.EQ.'T' ) then
        CALL SETERR(ABSERR,RELERR)
c
c  V = PRINT CURRENT PARAMETERS
c
      else if ( COMMAN.EQ.'V' ) then
        CALL PRIVAL(ABSERR,EPMACH,HSTEP,ICON,IDAMP,IDIF,IHOLD,
     1  IJAC,INORM,IPRINT,ITARG,MAXSTP,MODNEW,RELERR,TARVAL)
c
c  W = PRINT WORK STATISTICS
c
      else if ( COMMAN.EQ.'W' ) then
        CALL PRIWRK(NFACT,NFPVAL,NFXVAL,NSOLVE)
c
c  X = SET X
c
      else if ( COMMAN.EQ.'X' ) then
        IF(NSTP.GT.0)write(*,1020)
        NSTP=0
        CALL SETX(ICON,IHOLD,NCON,NEQN,NSTP,XCNEW,XCOLD,XHOLD,
     1  XNNEW,XNOLD,XNPRED)
        CALL SETPRM(EPMACH,EPSQRT,FPFACT,FPRIME,FXBACK,FXPLUS,
     &  ICON,IDIF,IERROR,IHOLD,IJAC,IPIVOT,MAXEQN,MAXSTP,NEQN,
     &  NFACTS,NFPVAL,NFXVAL,NROW,NSTP,TAN,XDELT,XNNEW,XHOLD)
c
c  Y = SET CONTINUATION TARGET
c
      else if ( COMMAN.EQ.'Y' ) then
        write(*,1060)NEQN
        READ(*,*)ITARG
        write(*,1070)ITARG
        READ(*,*)TARVAL
c
c  Unrecognized command
c
      else if ( COMMAN.NE.' ' ) then
        WRITE(*,*)'Your command was not recognized!'
	end if
c
      GO TO 10
c
 1000 FORMAT(' THE NEWTON ITERATION CANNOT BE CONTINUED NOW')
 1010 FORMAT(' YOU MAY TRY TO BEGIN INSTEAD')
 1020 FORMAT(' WARNING - THE CURRENT NEWTON ITERATION IS NOW ABORTED')
 1030 FORMAT(' WANT TO USE NEW STEPSIZE TO PICK NEW STARTING POINT?')
 1040 FORMAT(' CURRENT PARAMETER IS ',I6)
 1050 FORMAT(' ENTER NEW PARAMETER BETWEEN 1 AND ',I6)
 1060 FORMAT(' ENTER TARGET INDEX BETWEEN 1 AND ',I2)
 1070 FORMAT(' ENTER TARGET VALUE FOR X(',I2,')')
 1080 FORMAT(' THIS IS THE JACOBIAN USED FOR THE CURRENT ITERATE')
 1090 FORMAT(' WHAT IS THE MOST NUMBER OF STEPS OF NEWTON')
 1100 FORMAT(' ITERATION TO BE TAKEN BEFORE I MUST RETURN TO YOU?')
 1110 FORMAT(' SET THE DAMPING OPTION')
 1120 FORMAT(' 0=NO DAMPING, 1=DAMPING')
 1130 FORMAT(' HOW MANY NEWTON STEPS ARE TO BE TAKEN BEFORE')
 1140 FORMAT(' THE JACOBIAN IS UPDATED?')
 1150 FORMAT(' 0 = PRINT OUT ONLY RESULT OF ITERATION')
 1160 FORMAT(' 1 = PRINT NORMS OF X, FX AND STEP')
 1170 FORMAT(' 2 = PRINT X, FX AND STEP')
 1180 FORMAT(' COMMAND? L FOR LIST',$)
 1190 FORMAT(A1)
 1200 FORMAT(' ? I DID NOT UNDERSTAND THAT COMMAND')
 1210 FORMAT(' B = BEGIN NEWTON ITERATION')
 1220 FORMAT(' C = CONTINUE CURRENT NEWTON ITERATION')
 1230 FORMAT(' D = SET DAMPING OPTION')
 1240 FORMAT(' E = SET NUMBER OF EQUATIONS')
 1250 FORMAT(' F = SET FINITE DIFFERENCE OPTION')
 1260 FORMAT(' G = COMPARE JACOBIAN WITH FINITE DIFFERENCE JACOBIAN')
 1270 FORMAT(' H = DEFINE STEPSIZES FOR CONTINUATION')
 1280 FORMAT(' I = TAKE A CONTINUATION STEP')
 1290 FORMAT(' J = PRINT JACOBIAN')
 1300 FORMAT(' K = DEFINE HOW OFTEN JACOBIAN IS EVALUATED')
 1310 FORMAT(' L = LIST COMMANDS')
 1320 FORMAT(' M = SET MAXIMUM NUMBER OF STEPS')
 1330 FORMAT(' N = SET NORM TO USE')
 1340 FORMAT(' O = CHANGE AMOUNT OF OUTPUT')
 1350 FORMAT(' P = PRINT X, F(X)')
 1360 FORMAT(' Q = CHANGE PARAMETERS')
 1370 FORMAT(' S = STOP')
 1380 FORMAT(' T = SET ERROR TOLERANCES')
 1390 FORMAT(' V = PRINT VALUES OF CURRENT PARAMETERS')
 1400 FORMAT(' W = PRINT WORK STATISTICS')
 1410 FORMAT(' X = SET CURRENT X APPROXIMATION')
 1420 FORMAT(' Y = SET CONTINUATION TARGET VALUE')
      END
      SUBROUTINE ACCEPT(ABSERR,FNORM1,ILIKE,NSTP,RELERR,SNORM1,XNORM1)

c*********************************************************************72
c
cc ACCEPT ???
c
      ILIKE=0
      IF(NSTP.EQ.0 ) then
        IF(FNORM1.LE.ABSERR)ILIKE=1
      ELSE
        IF(FNORM1.LE.ABSERR.AND.SNORM1.LE.RELERR*(XNORM1+ABSERR))
     1  ILIKE=1
        end if
      RETURN
      END
      SUBROUTINE CHKJAC(EPSQRT,FPFACT,FPRIME,NROW,NEQN)

c*********************************************************************72
c
cc CHKJAC
c
      DIMENSION FPFACT(NROW,NEQN)
      DIMENSION FPRIME(NROW,NEQN)

      IFOUND=0
      DO 20 I=1,NEQN
        DO 10 J=1,NEQN
          DIF=FPFACT(I,J)-FPRIME(I,J)
          IF(ABS(DIF).LE.EPSQRT*(ABS(FPFACT(I,J))+EPSQRT+1.0))GO TO 10
          IFOUND=IFOUND+1
          IF(IFOUND.EQ.1 ) then
	    WRITE(*,*)' '
            write(*,1000)
            WRITE(*,*)' '
	    end if
          write(*,1010)I,J,FPFACT(I,J),FPRIME(I,J),DIF
   10     CONTINUE
   20   CONTINUE
      RETURN
 1000 FORMAT('   I  J  TRUE JACOBIAN   F-D JACOBIAN  DIFFERENCE')
 1010 FORMAT(1X,I3,I3,3G14.6)
      END
      SUBROUTINE CHRCAP(STRING,NWORDS,NCHPW)

c*********************************************************************72
c
cC CHRCAP ACCEPTS A VECTOR OF CHARACTERS, AND REPLACES ANY
c  LOWERCASE LETTERS BY UPPERCASE ONES.
c
      CHARACTER*(*) STRING(NWORDS)
      CHARACTER*1 INCHAR

      DO 20 I=1,NWORDS
        DO 10 J=1,NCHPW
          INCHAR=STRING(I)(J:J)
          IF('a'.GT.INCHAR)GO TO 10
          IF(INCHAR.GT.'z')GO TO 10
          ILITC=ICHAR(INCHAR)
          ILITA=ICHAR('a')
          IBIGA=ICHAR('A')
          IBIGC=IBIGA+(ILITC-ILITA)
          INCHAR=CHAR(IBIGC)
          STRING(I)(J:J)=INCHAR
   10     CONTINUE
   20   CONTINUE
      RETURN
      END
      SUBROUTINE DIFJAC(EPMACH,EPSQRT,FPRIME,FXBACK,FXPLUS,
     & ICON,IDIF,IHOLD,NEQN,NFXVAL,NROW,XDELT,XNNEW,XHOLD)

c*********************************************************************72
c
cC DIFJAC COMPUTES THE JACOBIAN OF A FUNCTION.
c
c  Discussion:
c
c    BACKWARD, CENTERED, OR FORWARD DIFFERENCES may be used.
c
c  Parameters:
c
c  FPRIME - AN ARRAY OF DIMENSION (NROW,NCOL), INTO WHICH THE
c           RESULTS OF THE COMPUTATION ARE TO BE STORED.
c
c  FXBACK - AN ARRAY OF DIMENSION NEQN, USED FOR THE COMPUTATION.
c           ON RETURN, IF IDIF.EQ.1 OR -1, FX CONTAINS F(X).
c
c  FXPLUS - AN ARRAY OF DIMENSION NEQN, USED FOR THE COMPUTATION.
c
c  IDIF   - BACKWARD, CENTERED, OR FORWARD DIFFERENCE SWITCH.
c           IDIF=-1, BACKWARD DIFFERENCES ARE USED.
c           IDIF= 0, CENTERED DIFFERENCES ARE USED.
c           IDIF=+1, FORWARD DIFFERENCES ARE USED.
c
c  NCOL   - THE NUMBER OF COLUMNS IN FPRIME.
c
c  NROW   - THE FORMAL FIRST DIMENSION OF FPRIME.
c
c  NEQN   - THE NUMBER OF VARIABLES.
c
c  X      - AN ARRAY OF DIMENSION NEQN, THE POINT AT WHICH
c           THE JACOBIAN IS TO BE EVALUATED.
c
c  XDELT  - AN ARRAY OF DIMENSION NEQN, USED AS
c           WORKSPACE BY THE PROGRAM.
c
      INTEGER NEQN
      INTEGER NROW

      DIMENSION FPRIME(NROW,NEQN)
      DIMENSION XNNEW(NEQN)
      DIMENSION FXBACK(NEQN)
      DIMENSION FXPLUS(NEQN)
      DIMENSION XDELT(NEQN)
c
c  EXCEPT FOR CENTERED DIFFERENCES, GET FX AT X
c
      IF(IDIF.EQ.1 ) then
        CALL GETFUN(FXBACK,ICON,IHOLD,NEQN,XNNEW,XHOLD)
        NFXVAL=NFXVAL+1
      end if

      IF(IDIF.EQ.-1 ) then
        CALL GETFUN(FXPLUS,ICON,IHOLD,NEQN,XNNEW,XHOLD)
        NFXVAL=NFXVAL+1
      end if
c
c  COMPUTATION OF APPROXIMATE JACOBIAN
c
      DO 60 J=1,NEQN
        IF(IDIF.EQ.-1)GO TO 20
c
c  FOR CENTERED OR FORWARD DIFFERENCES, INCREMENT UNRELATED X'S
c
        DO 10 I=1,NEQN
          XDELT(I)=XNNEW(I)
   10     CONTINUE
        XNORM=ABS(XNNEW(J))
        H=EPSQRT
        IF(XNORM.GE.1.0)H=EPSQRT*XNORM
        XDELT(J)=XDELT(J)+H
c
c  EVALUATE FUNCTION
c
        CALL GETFUN(FXPLUS,ICON,IHOLD,NEQN,XDELT,XHOLD)
        NFXVAL=NFXVAL+1
c
c  FOR BACKWARD OR CENTERED DIFFERENCES, DECREMENT UNRELATED X'S
c
        IF(IDIF.NE.0)GO TO 40
   20   CONTINUE
        DO I=1,NEQN
          XDELT(I)=XNNEW(I)
        end do
        XNORM=ABS(XNNEW(J))
        H=EPSQRT
        IF(XNORM.GE.1.0)H=EPSQRT*XNORM
        XDELT(J)=XDELT(J)-H
c
c  EVALUATE FUNCTION
c
        CALL GETFUN(FXBACK,ICON,IHOLD,NEQN,XDELT,XHOLD)
        NFXVAL=NFXVAL+1
   40   CONTINUE
c
c  COMPUTE JACOBIAN COLUMNS
c
        IF(IDIF.EQ.-1)H=XNNEW(J)-XDELT(J)
        IF(IDIF.EQ.0)H=2.0*(XNNEW(J)-XDELT(J))
        IF(IDIF.EQ.1)H=XDELT(J)-XNNEW(J)
        DO I=1,NEQN
          ENTRY=(FXPLUS(I)-FXBACK(I))/H
          FPRIME(I,J)=ENTRY
        end do
   60   CONTINUE

      RETURN
      END
      SUBROUTINE ERROR(IERROR,MAXEQN,MAXSTP,NEQN,NSTP)

c*********************************************************************72
c
cC ERROR prints out error messages.
c
      IF(IERROR.EQ.0)RETURN
      IF(IERROR.EQ.1)GO TO 10
      IF(IERROR.EQ.2)GO TO 20
      IF(IERROR.EQ.3)GO TO 30
      RETURN
   10 CONTINUE
      write(*,1000)NEQN
      write(*,1010)MAXEQN
      RETURN
   20 CONTINUE
      write(*,1020)NSTP
      write(*,1030)MAXSTP
      RETURN
   30 CONTINUE
      write(*,1040)
      RETURN
 1000 FORMAT(' ILLEGAL NUMBER OF EQUATIONS=',I6)
 1010 FORMAT(' IT SHOULD BE BETWEEN 1 AND ',I6)
 1020 FORMAT(' TOO MANY ITERATIONS WITHOUT CONVERGENCE=',I6)
 1030 FORMAT(' THE MAXIMUM ALLOWED IS ',I6)
 1040 FORMAT(' THE MATRIX IS SINGULAR')
      END
      SUBROUTINE FACTOR(FPFACT,FPRIME,ICON,IERROR,IHOLD,IPIVOT,MAXEQN,
     1 MAXSTP,NEQN,NFACTS,NROW,NSTP)

c*********************************************************************72
c
cC FACTOR factors the matrix.
c
      INTEGER NEQN
      INTEGER NROW

      DIMENSION FPFACT(NROW,NEQN)
      DIMENSION FPRIME(NROW,NEQN)
      DIMENSION IPIVOT(NEQN)
c
      IERROR=0
      NFACTS=NFACTS+1
c
c  LOAD JACOBIAN INTO FACTORED JACOBIAN STORAGE
c
      DO 20 I=1,MAXEQN
        DO 10 J=1,MAXEQN
          FPFACT(I,J)=FPRIME(I,J)
   10     CONTINUE
   20   CONTINUE
      CALL SGEFA(FPFACT,NROW,NEQN,IPIVOT,INFO,ICON,IHOLD)
      IF(INFO.NE.0)IERROR=3
      CALL ERROR(IERROR,MAXEQN,MAXSTP,NEQN,NSTP)
      RETURN
      END
      SUBROUTINE GETFUN(FX,ICON,IHOLD,NEQN,X,XHOLD)

c*********************************************************************72
c
cC GETFUN ???
c
      DIMENSION FX(NEQN)
      DIMENSION X(NEQN)

      CALL FUNK(FX,NEQN,X)
      IF(ICON.EQ.0)RETURN
      IF(IHOLD.NE.0)FX(NEQN)=X(IHOLD)-XHOLD
      RETURN
      END
      SUBROUTINE GETJAC(EPMACH,EPSQRT,FPFACT,FPRIME,FXBACK,FXPLUS,
     & ICON,IDIF,IHOLD,IJAC,IPIVOT,NEQN,NFPVAL,NFXVAL,NROW,XDELT,
     & XNNEW,XHOLD)

c*********************************************************************72
c
cc GETJAC
c
      INTEGER NEQN
      INTEGER NROW

      REAL    FPFACT(NROW,NEQN)
      REAL    FPRIME(NROW,NEQN)
      REAL    FXBACK(NEQN)
      REAL    FXPLUS(NEQN)
      INTEGER IPIVOT(NEQN)
      REAL    XDELT(NEQN)
      REAL    XNNEW(NEQN)

      IF(IJAC.EQ.0 ) then
        CALL JACOB(FPRIME,NROW,NEQN,XNNEW)
        IF(ICON.EQ.1 ) then
          DO I=1,NEQN
            FPRIME(NEQN,I)=0.0
          end do
          IF(IHOLD.NE.0)FPRIME(NEQN,IHOLD)=1.0
          end if
        NFPVAL=NFPVAL+1
      ELSE
        CALL DIFJAC(EPMACH,EPSQRT,FPRIME,FXBACK,FXPLUS,
     &  ICON,IDIF,IHOLD,NEQN,NFXVAL,NROW,XDELT,XNNEW,XHOLD)
        end if
      RETURN
      END
      SUBROUTINE GETNRM(DAMP,DELX,FNORM0,FNORM1,FX,INORM,NEQN,NSTP,
     1 SNORM0,SNORM1,XNNEW,XNORM0,XNORM1)

c*********************************************************************72
c
cc GETNRM
c
      DIMENSION DELX(NEQN)
      DIMENSION FX(NEQN)
      DIMENSION XNNEW(NEQN)

      XNORM0=XNORM1
      CALL NORM(INORM,NEQN,XNNEW,XNORM1)
      IF(DAMP.EQ.1.0)FNORM0=FNORM1
      CALL NORM(INORM,NEQN,FX,FNORM1)
      IF(NSTP.EQ.0)GO TO 10
      SNORM0=SNORM1
      CALL NORM(INORM,NEQN,DELX,SNORM1)
      RETURN
   10 CONTINUE
      SNORM0=0.0
      RETURN
      END
      SUBROUTINE GETTAN(EPMACH,EPSQRT,FPFACT,FPRIME,FXBACK,FXPLUS,
     & ICON,IDIF,IERROR,IHOLD,IJAC,IPIVOT,ITMAX,MAXEQN,MAXSTP,
     & NEQN,NFACTS,NFPVAL,NFXVAL,NROW,NSOLVE,NSTP,TAN,TANOLD,
     & XDELT,XNNEW,XHOLD)

c*********************************************************************72
c
cc GETTAN ???
c
      INTEGER NEQN
      INTEGER NROW

      DIMENSION FPFACT(NROW,NEQN)
      DIMENSION FPRIME(NROW,NEQN)
      REAL    FXBACK(NEQN)
      REAL    FXPLUS(NEQN)
      DIMENSION IPIVOT(NEQN)
      DIMENSION TAN(NEQN)
      DIMENSION TANOLD(NEQN)
      REAL    XDELT(NEQN)
      DIMENSION XNNEW(NEQN)
c
c  Set right hand side for tangent system
c
      DO 10 I=1,NEQN
        TAN(I)=0.0
   10   CONTINUE
      TAN(NEQN)=1.0
c
c  Get jacobian
c
      CALL GETJAC(EPMACH,EPSQRT,FPFACT,FPRIME,FXBACK,FXPLUS,
     & ICON,IDIF,IHOLD,IJAC,IPIVOT,NEQN,NFPVAL,NFXVAL,NROW,
     & XDELT,XNNEW,XHOLD)
c
c  Factor jacobian
c
      CALL FACTOR(FPFACT,FPRIME,ICON,IERROR,IHOLD,IPIVOT,MAXEQN,
     1 MAXSTP,NEQN,NFACTS,NROW,NSTP)
      IF(IERROR.NE.0 ) then
        WRITE(*,*)'GETTAN - FATAL ERROR!'
	WRITE(*,*)'Could not compute the tangent vector, because'
	WRITE(*,*)'the jacobian is singular!'
	RETURN
	end if
c
c  Solve the linear system for the tangent vector.
c
      CALL SGESL(FPFACT,NROW,NEQN,IPIVOT,TAN,0)
      NSOLVE=NSOLVE+1
c
c  Normalize the tangent vector
c
      INDEX=2
      CALL NORM(INDEX,NEQN,TAN,TANNRM)
      TANDOT=SDOT(NEQN,TAN,1,TANOLD,1)
      IF(TANDOT.LT.0.0)TANNRM=-TANNRM
      TANMAX=0.0
      ITMAX=0
      DO 20 I=1,NEQN
        TAN(I)=TAN(I)/TANNRM
        IF(ABS(TAN(I)).GT.TANMAX ) then
          TANMAX=ABS(TAN(I))
          ITMAX=I
          end if
   20   CONTINUE
      RETURN
      END
      SUBROUTINE NEWCON(ABSERR,DAMP,DELX,EPMACH,EPSQRT,FNORM0,FNORM1,
     1 FPFACT,FPRIME,FX,FXBACK,FXPLUS,HSTEP,ICON,IDAMP,IDIF,IERROR,
     & IHATE,IHOLD,IJAC,ILIKE,INORM,IPIVOT,IPRINT,MAXEQN,MAXSTP,
     & MODNEW,NEQN,NFACTS,NFPVAL,NFXVAL,NROW,NSOLVE,NSTP,RELERR,
     & SNORM0,SNORM1,TAN,XNNEW,XHOLD,XCOLD,XDELT,XNORM0,XNORM1,
     & XNOLD,XNPRED)

c*********************************************************************72
c
cc NEWCON ???
c
      INTEGER NEQN
      INTEGER NROW

      DIMENSION DELX(NEQN)
      DIMENSION FPFACT(NROW,NEQN)
      DIMENSION FPRIME(NROW,NEQN)
      DIMENSION FX(NEQN)
      REAL    FXBACK(NEQN)
      REAL    FXPLUS(NEQN)
      DIMENSION IPIVOT(NEQN)
      DIMENSION TAN(NEQN)
      DIMENSION XNNEW(NEQN)
      DIMENSION XCOLD(NEQN)
      REAL    XDELT(NEQN)
      DIMENSION XNOLD(NEQN)
      DIMENSION XNPRED(NEQN)

      ILIKE=0
      ISKIP=0
      IF(NSTP.EQ.0)IACCPT=0
      IHATE=0
      MSTP=0
      DAMP=1.0
c
c  Should we accept the initial X?
c
      IF(NSTP.EQ.0 ) then
        CALL GETFUN(FX,ICON,IHOLD,NEQN,XNNEW,XHOLD)
        NFXVAL=NFXVAL+1
        CALL GETNRM(DAMP,DELX,FNORM0,FNORM1,FX,INORM,NEQN,NSTP,
     1  SNORM0,SNORM1,XNNEW,XNORM0,XNORM1)
        CALL ACCEPT(ABSERR,FNORM1,ILIKE,NSTP,RELERR,SNORM1,XNORM1)
        IACCPT=IACCPT+1
        CALL OUTPUT(DAMP,DELX,FNORM1,FX,IDAMP,IHATE,ILIKE,IPRINT,
     1  ISKIP,NEQN,NSTP,SNORM1,XNNEW,XNORM1)
        IF(ILIKE.EQ.1)RETURN
	end if
c
c  NEWTON ITERATION LOOP
c
   10 CONTINUE
      NSTP=NSTP+1
      MSTP=MSTP+1
      IF(MSTP.LE.MAXSTP)GO TO 20
      write(*,1000)
      RETURN
   20 CONTINUE
c
c  IS IT TIME TO REEVALUATE THE JACOBIAN?
c
      ISKIP=0
      IF(MODNEW.EQ.0)GO TO 30
      IF(MODNEW.EQ.1)GO TO 30
      ISKIP=MOD(IACCPT,MODNEW)-1
      IF(ISKIP.EQ.0)GO TO 30
      ISKIP=1
      GO TO 40
   30 CONTINUE
c
c  Get jacobian
c
      CALL GETJAC(EPMACH,EPSQRT,FPFACT,FPRIME,FXBACK,FXPLUS,
     & ICON,IDIF,IHOLD,IJAC,IPIVOT,NEQN,NFPVAL,NFXVAL,NROW,XDELT,
     & XNNEW,XHOLD)
c
c  Factor jacobian
c
      CALL FACTOR(FPFACT,FPRIME,ICON,IERROR,IHOLD,IPIVOT,MAXEQN,
     1 MAXSTP,NEQN,NFACTS,NROW,NSTP)
      IF(IERROR.NE.0)RETURN
   40 CONTINUE
c
c  SOLVE NEWTON SYSTEM, UPDATE X, AND CHECK FOR ACCEPTANCE
c  OR REJECTION OF CURRENT ITERATE
c
      DO 50 I=1,NEQN
        DELX(I)=-FX(I)
   50   CONTINUE
      CALL SGESL(FPFACT,NROW,NEQN,IPIVOT,DELX,0)
      NSOLVE=NSOLVE+1
c
c  Damping loop begins here
c
      DAMP=1.0
   60 CONTINUE
      CALL UPDATE(DAMP,DELX,NEQN,XNNEW,XNOLD)
      CALL GETFUN(FX,ICON,IHOLD,NEQN,XNNEW,XHOLD)
      NFXVAL=NFXVAL+1
      CALL GETNRM(DAMP,DELX,FNORM0,FNORM1,FX,INORM,NEQN,NSTP,
     1 SNORM0,SNORM1,XNNEW,XNORM0,XNORM1)
      CALL ACCEPT(ABSERR,FNORM1,ILIKE,NSTP,RELERR,SNORM1,XNORM1)
      CALL REJECT(ABSERR,FNORM0,FNORM1,IHATE,NSTP,SNORM0,SNORM1)
      IF(IHATE.EQ.1.AND.IDAMP.EQ.1.AND.MSTP.LT.MAXSTP
     1 .AND.DAMP.GT.0.03125 ) then
        IHATE=0
        CALL OUTPUT(DAMP,DELX,FNORM1,FX,IDAMP,IHATE,ILIKE,IPRINT,
     1  ISKIP,NEQN,NSTP,SNORM1,XNNEW,XNORM1)
        DAMP=0.5*DAMP
        SNORM1=0.5*SNORM1
        DO 70 I=1,NEQN
          DELX(I)=0.5*DELX(I)
   70     CONTINUE
        ISKIP=1
        NSTP=NSTP+1
        MSTP=MSTP+1
        GO TO 60
        end if
c
      CALL OUTPUT(DAMP,DELX,FNORM1,FX,IDAMP,IHATE,ILIKE,IPRINT,
     1 ISKIP,NEQN,NSTP,SNORM1,XNNEW,XNORM1)
      IACCPT=IACCPT+1
      DAMP=1.0
      IF(ILIKE.EQ.1)RETURN
      IF(IHATE.EQ.1)RETURN
      GO TO 10
 1000 FORMAT(' ITERATION HALTED.  REACHED MAXIMUM NUMBER OF STEPS')
      END
      SUBROUTINE NORM(INORM,NVEC,VEC,VECNRM)

c*********************************************************************72
c
cc NORM ???
c
      integer nvec

      DIMENSION VEC(NVEC)

      VECNRM=0.0
      IF(NVEC.LE.0)RETURN
      IMAX=1
      VECNRM=ABS(VEC(1))
      IF(NVEC.EQ.1)RETURN
      IF(INORM.EQ.0)GO TO 10
      IF(INORM.EQ.1)GO TO 30
      IF(INORM.EQ.2)GO TO 50
c
c  MAXIMUM NORM
c
   10 CONTINUE
      DO 20 I=2,NVEC
        TEMP=ABS(VEC(I))
        IF(TEMP.LT.VECNRM)GO TO 20
        IMAX=I
        VECNRM=ABS(VEC(IMAX))
   20   CONTINUE
      RETURN
c
c  SUM OF ABSOLUTE VALUES
c
   30 CONTINUE
      VECNRM=0.0
      DO I=1,NVEC
        VECNRM=VECNRM+ABS(VEC(I))
      end do
      RETURN
c
c  SQUARE ROOT OF SUM OF SQUARES
c
   50 CONTINUE
      VECNRM=0.0
      DO 60 I=1,NVEC
        VECNRM=VECNRM+VEC(I)*VEC(I)
   60   CONTINUE
      VECNRM=SQRT(VECNRM)
      RETURN
      END
      SUBROUTINE OUTPUT(DAMP,DELX,FNORM1,FX,IDAMP,IHATE,ILIKE,IPRINT,
     1 ISKIP,NEQN,NSTP,SNORM1,XNNEW,XNORM1)

c*********************************************************************72
c
cc OUTPUT ???
c
      DIMENSION DELX(NEQN)
      DIMENSION FX(NEQN)
      DIMENSION XNNEW(NEQN)
      CHARACTER LABEL*3
c
      LABEL='OLD'
      IF(ISKIP.EQ.0)LABEL='NEW'
      IF(NSTP.EQ.0)WRITE(*,*)' '
      IF(IPRINT.LE.0)GO TO 30
      IF(NSTP.EQ.0.AND.IPRINT.EQ.1.AND.IDAMP.EQ.0)write(*,1000)
      IF(NSTP.EQ.0.AND.IPRINT.EQ.1.AND.IDAMP.NE.0)write(*,1010)
      IF(NSTP.EQ.0.AND.IPRINT.GE.2.AND.IDAMP.EQ.0)write(*,1020)
      IF(NSTP.EQ.0.AND.IPRINT.GE.2.AND.IDAMP.NE.0)write(*,1030)
      IF(NSTP.GT.0.AND.IPRINT.EQ.1.AND.IDAMP.EQ.0)
     1 write(*,1060)SNORM1,LABEL
      IF(NSTP.GT.0.AND.IPRINT.EQ.1.AND.IDAMP.NE.0)
     1 write(*,1060)SNORM1,LABEL,DAMP
      IF(NSTP.GT.0.AND.IPRINT.GE.2 ) then
        DO 10 I=1,NEQN
          IF(IDAMP.EQ.0.AND.I.EQ.1)write(*,1060)DELX(I),LABEL
          IF(IDAMP.EQ.1.AND.I.EQ.1)write(*,1060)DELX(I),LABEL,DAMP
          IF(IDAMP.EQ.0.AND.I.NE.1)write(*,1060)DELX(I)
          IF(IDAMP.EQ.1.AND.I.NE.1)write(*,1060)DELX(I)
   10     CONTINUE
        end if
      IF(IPRINT.EQ.1)write(*,1040)NSTP,XNORM1,FNORM1
      IF(IPRINT.GE.2 ) then
        DO 20 I=1,NEQN
          IF(I.EQ.1)write(*,1040)NSTP,XNNEW(I),FX(I)
          IF(I.NE.1)write(*,1050)     XNNEW(I),FX(I)
   20     CONTINUE
        end if
   30 CONTINUE
      IF(ILIKE.EQ.1.OR.IHATE.EQ.1)WRITE(*,*)' '
      IF(ILIKE.EQ.1)write(*,1070)
      IF(IHATE.EQ.1)write(*,1080)
      IF(IHATE.EQ.1)write(*,1090)
      RETURN
 1000 FORMAT('   STEP    X-NORM       FX-NORM        STEP-NORM'
     1 ,' JACOBIAN')
 1010 FORMAT('   STEP    X-NORM',7X,'FX-NORM',8X,'STEP-NORM',
     1 1X,'JACOBIAN',1X,'DAMPING')
 1020 FORMAT('   STEP',6X,'X',12X,'FX             STEP     JACOBIAN')
 1030 FORMAT('   STEP',6X,'X',12X,'FX             STEP     JACOBIAN'
     1 ,' DAMPING')
 1040 FORMAT(1X,I6,2G14.6)
 1050 FORMAT(1X,6X,2G14.6)
 1060 FORMAT(1X,6X,28X,G14.6,1X,A3,1X,G14.6)
 1070 FORMAT(' ITERATION HALTED.  POINT ACCEPTED')
 1080 FORMAT(' ITERATION HALTED.  A NORM INCREASED')
 1090 FORMAT(' YOU MUST DECIDE WHETHER TO GO ON OR NOT.')
      END
      SUBROUTINE PREDIC(ABSERR,EPMACH,EPSQRT,FPFACT,FPRIME,FX,FXBACK,
     & FXPLUS,HSTEP,ICON,IDIF,IERROR,IHOLD,IJAC,INORM,IPIVOT,ITARG,
     & MAXEQN,MAXSTP,NEQN,NFACTS,NFPVAL,NFXVAL,NROW,NSOLVE,NSTP,TAN,
     & TANOLD,TARVAL,XCNEW,XDELT,XHOLD,XCOLD,XNPRED)

c*********************************************************************72
c
cc PREDIC ???
c
      INTEGER NEQN
      INTEGER NROW
c
      DIMENSION FPFACT(NROW,NEQN)
      DIMENSION FPRIME(NROW,NEQN)
      DIMENSION FX(NEQN)
      REAL    FXBACK(NEQN)
      REAL    FXPLUS(NEQN)
      DIMENSION IPIVOT(NEQN)
      DIMENSION TAN(NEQN)
      DIMENSION TANOLD(NEQN)
      DIMENSION XCNEW(NEQN)
      DIMENSION XCOLD(NEQN)
      REAL    XDELT(NEQN)
      DIMENSION XNPRED(NEQN)
c
c  NO UPDATE IF STEP SIZE CHANGE
c
      IF(NSTP.GT.0)GO TO 20
c
c  REFUSE TO GO IF NORM TOO LARGE
c
      CALL GETFUN(FX,ICON,IHOLD,NEQN,XCNEW,XHOLD)
      CALL NORM(INORM,NEQN,FX,FXNORM)
      IF(FXNORM.GT.ABSERR ) then
        write(*,1010)
        write(*,1020)FXNORM
        write(*,1030)
        write(*,1040)
        RETURN
        end if
c
c  SAVE OLD POINT
c
      DO 10 I=1,NEQN
        XCOLD(I)=XCNEW(I)
        TANOLD(I)=TAN(I)
   10   CONTINUE
c
c  GET TANGENT VECTOR
c
   20 CONTINUE
      CALL GETTAN(EPMACH,EPSQRT,FPFACT,FPRIME,FXBACK,FXPLUS,
     & ICON,IDIF,IERROR,IHOLD,IJAC,IPIVOT,ITMAX,MAXEQN,MAXSTP,
     & NEQN,NFACTS,NFPVAL,NFXVAL,NROW,NSOLVE,NSTP,TAN,TANOLD,
     & XDELT,XCOLD,XHOLD)
      IF(IERROR.NE.0)RETURN
c
c  CHECK TO SEE IF WE CAN HIT A TARGET POINT BY REDUCING STEPSIZE
c
      HUSE=HSTEP
      IF(ITARG.EQ.0)GO TO 30
      XTOLD=XCOLD(ITARG)
      XTNEW=XCOLD(ITARG)+HUSE*TAN(ITARG)
      write(*,1701)XTOLD,XTNEW
1701  FORMAT(' XTOLD,XTNEW=',2G14.6)
      IF(TARVAL.LE.XTOLD.AND.TARVAL.LE.XTNEW)GO TO 30
      IF(TARVAL.GE.XTOLD.AND.TARVAL.GE.XTNEW)GO TO 30
      write(*,1000)
      ITMAX=ITARG
      HUSE=(TARVAL-XCOLD(ITARG))/TAN(ITARG)
   30 CONTINUE
      DO 40 I=1,NEQN
        XNPRED(I)=XCOLD(I)+HUSE*TAN(I)
   40   CONTINUE
c
c  UPDATE PARAMETER AND RIGHT HAND SIDE
c
      IHOLD=ITMAX
      XHOLD=XNPRED(IHOLD)
      WRITE(*,*)' '
      write(*,1050)HUSE
      write(*,1060)IHOLD
      write(*,1070)XHOLD
      write(*,1080)
      RETURN
 1000 FORMAT(' TARGET POINT SPOTTED')
 1010 FORMAT(' NO, I WILL NOT CONTINUE FROM THIS POINT')
 1020 FORMAT(' ITS NORM IS TOO LARGE, FXNORM=',G14.6)
 1030 FORMAT(' FIRST CALL NEWTON ON THIS POINT, THEN WE')
 1040 FORMAT(' CAN CONSIDER CONTINUATION')
 1050 FORMAT(' A STEP OF SIZE ',G14.6,' WILL BE TRIED')
 1060 FORMAT(' THE PARAMETER IS INDEX',I6)
 1070 FORMAT(' WHOSE VALUE WILL BE FIXED AT',G14.6)
 1080 FORMAT(' YOU MAY NOW BEGIN NEWTON ITERATION')
      END
      SUBROUTINE PRIJAC(FPRIME,NROW,NEQN)

c*********************************************************************72
c
cc PRIJAC prints the jacobian matrix.
c
      DIMENSION FPRIME(NROW,NEQN)
c
      WRITE(*,*)' '
      write(*,*)'Current jacobian matrix:'
      WRITE(*,*)' '
      DO 10 I=1,NEQN
        write(*,1010)(FPRIME(I,J),J=1,NEQN)
   10   CONTINUE
      RETURN
 1010 FORMAT(1X,5G14.6)
      END
      SUBROUTINE PRIVAL(ABSERR,EPMACH,HSTEP,ICON,IDAMP,IDIF,IHOLD,
     1 IJAC,INORM,IPRINT,ITARG,MAXSTP,MODNEW,RELERR,TARVAL)

c*********************************************************************72
c
cc PRIVAL ???
c
      WRITE(*,*)' '
      write(*,1040)ABSERR
      write(*,1050)RELERR
      write(*,1051)EPMACH
      IF(IJAC.EQ.0)write(*,1060)
      IF(IJAC.NE.0)write(*,1070)
      IF(IJAC.NE.0 ) then
        IF(IDIF.EQ.-1)write(*,1080)
        IF(IDIF.EQ.0)write(*,1090)
        IF(IDIF.EQ.1)write(*,1100)
        end if
      write(*,1110)IPRINT
      write(*,1120)MAXSTP
      write(*,1130)MODNEW
      IF(IDAMP.NE.0)write(*,1150)
      IF(IDAMP.EQ.0)write(*,1140)
      IF(INORM.EQ.0)write(*,1160)
      IF(INORM.EQ.1)write(*,1170)
      IF(INORM.EQ.2)write(*,1180)
      IF(ICON.EQ.0)RETURN
      WRITE(*,*)' '
      write(*,1000)
      write(*,1010)IHOLD
      write(*,1020)HSTEP
      IF(ITARG.NE.0)write(*,1031)ITARG,TARVAL
1031  FORMAT(' SEEKING POINT WITH INDEX',I3,' EQUAL TO',G14.6)
      RETURN
 1000 FORMAT(' CONTINUATION OPTION IS BEING USED')
 1010 FORMAT(' CURRENT PARAMETER INDEX IS ',I6)
 1020 FORMAT(' CURRENT STEPSIZE IS',G14.6)
 1040 FORMAT(' ABSOLUTE ERROR TOLERANCE=',G14.6)
 1050 FORMAT(' RELATIVE ERROR TOLERANCE=',G14.6)
1051  FORMAT(' RELATIVE MACHINE ACCURACY=',G14.6)
 1060 FORMAT(' USER SUPPLIES JACOBIAN')
 1070 FORMAT(' JACOBIAN COMPUTED VIA FINITE DIFFERENCES')
 1080 FORMAT(' BACKWARD DIFFERENCES USED')
 1090 FORMAT(' CENTERED DIFFERENCES USED')
 1100 FORMAT(' FORWARD DIFFERENCES USED')
 1110 FORMAT(' NEWTON OUTPUT (0-2) =',I6)
 1120 FORMAT(' MAXIMUM NUMBER OF STEPS BEFORE PAUSE=',I6)
 1130 FORMAT(' NUMBER OF STEPS BEFORE JACOBIAN REEVALUATED=',I6)
 1140 FORMAT(' DAMPING OPTION NOT BEING USED')
 1150 FORMAT(' DAMPING OPTION IS BEING USED')
 1160 FORMAT(' MAXIMUM VALUE NORM')
 1170 FORMAT(' SUM OF ABSOLUTE VALUES NORM')
 1180 FORMAT(' SQUARE ROOT OF SUM OF SQUARES NORM')
      END
      SUBROUTINE PRIXFX(EPMACH,EPSQRT,FPFACT,FPRIME,FX,FXBACK,
     & FXPLUS,ICON,IDIF,IERROR,IHOLD,IJAC,IPIVOT,MAXEQN,MAXSTP,
     & NEQN,NFACTS,NFPVAL,NFXVAL,NROW,NSOLVE,NSTP,TAN,TANTMP,
     & XDELT,XNNEW,XHOLD)

c*********************************************************************72
c
cc PRIXFX
c
      INTEGER NROW
      INTEGER NEQN
c
      DIMENSION FPFACT(NROW,NEQN)
      DIMENSION FPRIME(NROW,NEQN)
      DIMENSION FX(NEQN)
      REAL    FXBACK(NEQN)
      REAL    FXPLUS(NEQN)
      DIMENSION IPIVOT(NEQN)
      DIMENSION TAN(NEQN)
      DIMENSION TANTMP(NEQN)
      REAL    XDELT(NEQN)
      DIMENSION XNNEW(NEQN)
c
      WRITE(*,*)' '
      IF(ICON.EQ.0)write(*,1010)
      IF(ICON.EQ.1)write(*,1020)
      CALL GETFUN(FX,ICON,IHOLD,NEQN,XNNEW,XHOLD)
c
c  Get tangent, if continuation problem.
c
      IF(ICON.EQ.1)
     1 CALL GETTAN(EPMACH,EPSQRT,FPFACT,FPRIME,FXBACK,FXPLUS,
     & ICON,IDIF,IERROR,IHOLD,IJAC,IPIVOT,ITMAX,MAXEQN,MAXSTP,
     & NEQN,NFACTS,NFPVAL,NFXVAL,NROW,NSOLVE,NSTP,TANTMP,TAN,
     & XDELT,XNNEW,XHOLD)
      NFXVAL=NFXVAL+1
      DO 10 I=1,NEQN
        IF(ICON.EQ.0)write(*,1000)XNNEW(I),FX(I)
        IF(ICON.EQ.1)write(*,1000)XNNEW(I),FX(I),TAN(I)
   10   CONTINUE
      RETURN
 1000 FORMAT(1X,3G14.6)
 1010 FORMAT('        X             FX')
 1020 FORMAT('        X             FX            TAN')
      END
      SUBROUTINE PRIWRK(NFACTS,NFPVAL,NFXVAL,NSOLVE)

c*********************************************************************72
c
cc PRIWRK prints work statistics.
c
      WRITE(*,*)' '
      write(*,1000)
      write(*,1010)NFXVAL
      write(*,1020)NFPVAL
      write(*,1030)NFACTS
      write(*,1040)NSOLVE
      NFXVAL=0
      NFPVAL=0
      NFACTS=0
      NSOLVE=0
      RETURN
 1000 FORMAT(' WORK SINCE LAST PRINTOUT')
 1010 FORMAT(' FUNCTION EVALUATIONS=',I6)
 1020 FORMAT(' JACOBIAN EVALUATIONS=',I6)
 1030 FORMAT(' JACOBIAN FACTORIZATIONS=',I6)
 1040 FORMAT(' LINEAR SOLVES=',I6)
      END
      SUBROUTINE REJECT(ABSERR,FNORM0,FNORM1,IHATE,NSTP,SNORM0,SNORM1)

c*********************************************************************72
c
cc REJECT ???
c
      IF(NSTP.LE.0)RETURN
      IF(NSTP.EQ.1)GO TO 10
      IF(NSTP.GT.1)GO TO 20
   10 CONTINUE
      IF(FNORM1.GT.FNORM0+ABSERR)IHATE=1
      RETURN
   20 CONTINUE
      IF(FNORM1.GT.FNORM0+ABSERR)IHATE=1
      IF(SNORM1.GT.SNORM0+ABSERR)IHATE=1
      RETURN
      END
      SUBROUTINE SETCOM(ABSERR,DAMP,DELX,EPMACH,EPSQRT,FNORM0,FNORM1,
     1 FPFACT,FPRIME,FX,HSTEP,ICON,IDAMP,IDIF,IERROR,IHATE,
     2 IHOLD,IJAC,ILIKE,INORM,IPIVOT,IPRINT,ITARG,MAXEQN,MAXSTP,MODNEW,
     3 NEQN,NFACTS,NFPVAL,NFXVAL,NROW,NSOLVE,NSTP,RELERR,SNORM0,
     4 SNORM1,TAN,TARVAL,XNNEW,XHOLD,XCOLD,XNORM0,XNORM1,XNOLD,XNPRED)

c*********************************************************************72
c
cc SETCOM ???
c
      INTEGER MAXEQN
c
      DIMENSION DELX(MAXEQN)
      DIMENSION FPFACT(MAXEQN,MAXEQN)
      DIMENSION FPRIME(MAXEQN,MAXEQN)
      DIMENSION FX(MAXEQN)
      DIMENSION IPIVOT(MAXEQN)
      DIMENSION TAN(MAXEQN)
      DIMENSION XNNEW(MAXEQN)
      DIMENSION XCOLD(MAXEQN)
      DIMENSION XNOLD(MAXEQN)
      DIMENSION XNPRED(MAXEQN)
c
c  GET MACHINE NUMBERS
c
      EPMACH=1.0
   10 CONTINUE
      EPMACH=0.5*EPMACH
      IF(1.0+EPMACH.GT.1.0)GO TO 10
      EPMACH=EPMACH*2.0
c
c  INITIALIZE SCALARS
c
      ABSERR=0.00001
      DAMP=1.0
      EPSQRT=SQRT(EPMACH)
      FNORM0=0.0
      FNORM1=0.0
      HSTEP=0.25
      ICON=0
      IDAMP=0
      IDIF=0
      IERROR=0
      IHATE=0
      IHOLD=0
      IJAC=0
      ILIKE=0
      INORM=0
      IPRAM=0
      IPRINT=1
      ITARG=0
      MAXSTP=10
      MODNEW=1
      NEQN=0
      NFACTS=0
      NFPVAL=0
      NFXVAL=0
      NROW=MAXEQN
      NSOLVE=0
      NSTP=0
      RELERR=0.00001
      SNORM0=0.0
      SNORM1=0.0
      TARVAL=0.0
      XHOLD=0.0
      XNORM0=0.0
      XNORM1=0.0
c
c  INITIALIZE VECTORS
c
      DO 30 I=1,MAXEQN
        DELX(I)=0.0
        FX(I)=0.0
        IPIVOT(I)=0
        TAN(I)=0.0
        XNNEW(I)=0.0
        XNOLD(I)=0.0
        XNPRED(I)=0.0
        DO J=1,MAXEQN
          FPRIME(I,J)=0.0
        end do
   30   CONTINUE
      RETURN
      END
      SUBROUTINE SETCON(ICON)

c*********************************************************************72
c
cc SETCON ???
c
      CHARACTER ISAY*1
c
      write(*,*)'NEWTON can seek an N-dimensional point X'
      write(*,*)'satisfying N nonlinear equations.'
      write(*,*)'This is the standard option.'
      WRITE(*,*)' '
      write(*,1030)
      write(*,1040)
      write(*,1050)
      WRITE(*,*)' '
      write(*,1060)
      READ(*,1070)ISAY
      CALL CHRCAP(ISAY,1,1)
      IF(ISAY.EQ.'Y')ICON=1
      RETURN
 1030 FORMAT(' NEWTON CAN ALSO TRACE A CURVE OF POINTS')
 1040 FORMAT(' SATISFYING N EQUATIONS IN N+1 UNKNOWNS')
 1050 FORMAT(' THIS IS THE CONTINUATION OPTION')
 1060 FORMAT(' DO YOU WANT TO USE THE CONTINUATION OPTION?')
 1070 FORMAT(A1)
      END
      SUBROUTINE SETDIF(IDIF,IJAC)

c*********************************************************************72
c
cc SETDIF ???
c
      WRITE(*,*)' '
      write(*,1000)
      write(*,1010)
      write(*,1020)
      READ(*,*)IJAC
      IF(IJAC.EQ.0)RETURN
      IJAC=1
      write(*,1030)
      write(*,1040)
      write(*,1050)
      write(*,1060)
      READ(*,*)IDIF
      IF(IDIF.NE.-1.AND.IDIF.NE.1)IDIF=0
      RETURN
 1000 FORMAT(' TELL ME HOW TO EVALUATE THE JACOBIAN')
 1010 FORMAT(' 0 = USE THE SUBROUTINE YOU WROTE')
 1020 FORMAT(' 1 = USE FINITE DIFFERENCES')
 1030 FORMAT(' WHAT KIND OF FINITE DIFFERENCE?')
 1040 FORMAT(' -1 = BACKWARD DIFFERENCE')
 1050 FORMAT('  0 = CENTERED DIFFERENCE')
 1060 FORMAT('  1 = FORWARD DIFFERENCE')
      END
      SUBROUTINE SETERR(ABSERR,RELERR)

c*********************************************************************72
c
cc SETERR ???
c
      WRITE(*,*)' '
      write(*,1000)
      READ(*,*)ABSERR
      write(*,1010)
      READ(*,*)RELERR
      RETURN
 1000 FORMAT(' ENTER ABSOLUTE ERROR TOLERANCE')
 1010 FORMAT(' ENTER RELATIVE ERROR TOLERANCE')
      END
      SUBROUTINE SETNEQ(ICON,MAXEQN,NEQN)

c*********************************************************************72
c
cc SETNEQ ???
c
      WRITE(*,*)' '
      IF(ICON.EQ.1)write(*,1000)
      IF(ICON.EQ.1)write(*,1010)
      write(*,1020)MAXEQN
      READ(*,*)ITEMP
      IF(ITEMP.LE.0.OR.ITEMP.GT.MAXEQN)GO TO 10
      NEQN=ITEMP
      RETURN
   10 CONTINUE
      write(*,1030)
      RETURN
 1000 FORMAT(' WARNING - FOR CONTINUATION, WHICH YOU HAVE CHOSEN')
 1010 FORMAT(' INCREASE THE NUMBER OF EQUATIONS BY 1')
 1020 FORMAT(' ENTER NUMBER OF EQUATIONS UP TO',I6)
 1030 FORMAT(' UNACCEPTABLE NUMBER OF EQUATIONS')
      END
      SUBROUTINE SETNRM(INORM)

c*********************************************************************72
c
cc SETNRM ???
c
      WRITE(*,*)' '
      write(*,1000)
      write(*,1010)
      write(*,1020)
      write(*,1030)
      READ(*,*)INORM
      IF(INORM.LT.0.OR.INORM.GT.2)INORM=0
      RETURN
 1000 FORMAT(' PICK THE NORM TO USE FOR NEWTON')
 1010 FORMAT(' 0 = MAXIMUM ENTRY NORM')
 1020 FORMAT(' 1 = SUM OF ABSOLUTE VALUES')
 1030 FORMAT(' 2 = SQUARE ROOT OF SUM OF SQUARES')
      END
      SUBROUTINE SETPRM(EPMACH,EPSQRT,FPFACT,FPRIME,FXBACK,FXPLUS,
     & ICON,IDIF,IERROR,IHOLD,IJAC,IPIVOT,MAXEQN,MAXSTP,NEQN,NFACTS,
     & NFPVAL,NFXVAL,NROW,NSTP,TAN,XDELT,XNNEW,XHOLD)

c*********************************************************************72
c
cc SETPRM ???
c
      INTEGER NEQN
      INTEGER NROW

      DIMENSION FPFACT(NROW,NEQN)
      DIMENSION FPRIME(NROW,NEQN)
      REAL    FXBACK(NEQN)
      REAL    FXPLUS(NEQN)
      DIMENSION IPIVOT(NEQN)
      DIMENSION TAN(NEQN)
      REAL    XDELT(NEQN)
      DIMENSION XNNEW(NEQN)
c
      IF(ICON.NE.1)RETURN
c
c  Get the jacobian
c
      CALL GETJAC(EPMACH,EPSQRT,FPFACT,FPRIME,FXBACK,FXPLUS,
     & ICON,IDIF,IHOLD,IJAC,IPIVOT,NEQN,NFPVAL,NFXVAL,NROW,
     & XDELT,XNNEW,XHOLD)
c
c  Factor the jacobian
c
      CALL FACTOR(FPFACT,FPRIME,ICON,IERROR,IHOLD,IPIVOT,MAXEQN,
     1 MAXSTP,NEQN,NFACTS,NROW,NSTP)
      WRITE(*,*)' '
      write(*,1000)IHOLD
      XHOLD=XNNEW(IHOLD)
      DO 10 I=1,NEQN
        TAN(I)=0.0
   10   CONTINUE
      TAN(IHOLD)=1.0
      write(*,1010)XHOLD
      RETURN
 1000 FORMAT(' INITIAL PARAMETER IS INDEX=',I6)
 1010 FORMAT(' VALUE OF PARAMETER IS',G14.6)
      END
      SUBROUTINE SETSTP(HSTEP)

c*********************************************************************72
c
cc SETSTP ???
c
      WRITE(*,*)' '
      write(*,1000)
      READ(*,*,END=10,ERR=10)HSTEP
   10 CONTINUE
      RETURN
 1000 FORMAT(' ENTER STEPSIZE, OR S TO SKIP')
      END
      SUBROUTINE SETX(ICON,IHOLD,NCON,NEQN,NSTP,XCNEW,XCOLD,XHOLD,
     1 XNNEW,XNOLD,XNPRED)

c*********************************************************************72
c
cc SETX ???
c
      INTEGER NEQN
c
      DIMENSION XCNEW(NEQN)
      DIMENSION XCOLD(NEQN)
      DIMENSION XNNEW(NEQN)
      DIMENSION XNOLD(NEQN)
      DIMENSION XNPRED(NEQN)
c
c  NEW X POINT AUTOMATICALLY TERMINATES ANY CURRENT
c  CONTINUATION PROCESS
c
      NCON=0
      WRITE(*,*)' '
      write(*,1000)
      READ(*,*,END=20,ERR=20)(XNPRED(I),I=1,NEQN)
      DO 10 I=1,NEQN
        XCNEW(I)=XNPRED(I)
        XCOLD(I)=XNPRED(I)
        XNNEW(I)=XNPRED(I)
        XNOLD(I)=XNPRED(I)
   10   CONTINUE
   20 CONTINUE
      RETURN
 1000 FORMAT(' ENTER X VALUES')
      END
      SUBROUTINE UPDATE(DAMP,DELX,NEQN,XNNEW,XNOLD)

c*********************************************************************72
c
cc UPDATE ???
c
      INTEGER NEQN
c
      DIMENSION DELX(NEQN)
      DIMENSION XNNEW(NEQN)
      DIMENSION XNOLD(NEQN)
c
      IF(DAMP.EQ.1.0 ) then
        DO 10 I=1,NEQN
          XNOLD(I)=XNNEW(I)
   10     CONTINUE
        end if
      DO 20 I=1,NEQN
        XNNEW(I)=XNOLD(I)+DELX(I)
   20   CONTINUE
      RETURN
      END
      SUBROUTINE SGEFA(A,LDA,N,IPVT,INFO,ICON,IHOLD)

c*********************************************************************72
c
cC SGEFA FACTORS A MATRIX BY GAUSSIAN ELIMINATION.
c
c     SGEFA IS USUALLY CALLED BY SGECO, BUT IT CAN BE CALLED
c     DIRECTLY WITH A SAVING IN TIME IF  RCOND  IS NOT NEEDED.
c     (TIME FOR SGECO) = (1 + 9/N)*(TIME FOR SGEFA) .
c
c     ON ENTRY
c
c        A       DIMENSION(LDA, N)
c                THE MATRIX TO BE FACTORED.
c
c        LDA     INTEGER
c                THE LEADING DIMENSION OF THE ARRAY  A .
c
c        N       INTEGER
c                THE ORDER OF THE MATRIX  A .
c
c     ON RETURN
c
c        A       AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS
c                WHICH WERE USED TO OBTAIN IT.
c                THE FACTORIZATION CAN BE WRITTEN  A = L*U  WHERE
c                L  IS A PRODUCT OF PERMUTATION AND UNIT LOWER
c                TRIANGULAR MATRICES AND  U  IS UPPER TRIANGULAR.
c
c        IPVT    INTEGER(N)
c                AN INTEGER VECTOR OF PIVOT INDICES.
c
c        INFO    INTEGER
c                = 0  NORMAL VALUE.
c                = K  IF  U(K,K) .EQ. 0.0 .  THIS IS NOT AN ERROR
c                     CONDITION FOR THIS SUBROUTINE, BUT IT DOES
c                     INDICATE THAT SGESL OR SGEDI WILL DIVIDE BY ZERO
c                     IF CALLED.  USE  RCOND  IN SGECO FOR A RELIABLE
c                     INDICATION OF SINGULARITY.
c
c     LINPACK. THIS VERSION DATED 08/14/78 .
c     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
c
c     SUBROUTINES AND FUNCTIONS
c
c     BLAS SAXPY,SSCAL,ISAMAX
c
      INTEGER LDA,N,IPVT(1),INFO
      DIMENSION A(LDA,1)
      INTEGER ISAMAX,J,K,KP1,L,NM1
c
c
c     GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
c
      INFO = 0
      NM1=N-1
      IF (NM1 .LT. 1) GO TO 80
      DO 70 K = 1, NM1
         KP1 = K + 1
c
c  FIND L = PIVOT INDEX
c
         L=ISAMAX(N-K+1,A(K,K),1) + K - 1
         IPVT(K)=L
c
c  ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
c
         IF(A(L,K).EQ.0.0.AND.(ICON.EQ.0.OR.IHOLD.NE.0))GO TO 50
         IF(A(L,K).NE.0.0)GO TO 10
         L=N
         IPVT(K)=N
         A(NEQN,K)=1.0
         IHOLD=K
   10    CONTINUE
c
c  INTERCHANGE IF NECESSARY
c
            IF (L .EQ. K) GO TO 20
               T = A(L,K)
               A(L,K) = A(K,K)
               A(K,K) = T
   20       CONTINUE
c
c           COMPUTE MULTIPLIERS
c
            T = -1.0E0/A(K,K)
            CALL SSCAL(N-K,T,A(K+1,K),1)
c
c           ROW ELIMINATION WITH COLUMN INDEXING
c
            DO 40 J = KP1, N
               T = A(L,J)
               IF (L .EQ. K) GO TO 30
                  A(L,J) = A(K,J)
                  A(K,J) = T
   30          CONTINUE
               CALL SAXPY(N-K,T,A(K+1,K),1,A(K+1,J),1)
   40       CONTINUE
         GO TO 60
   50    CONTINUE
            INFO = K
   60    CONTINUE
   70 CONTINUE
   80 CONTINUE
      IPVT(N)=N
      IF (A(N,N).EQ.0.0) INFO = N
      IF(A(N,N).EQ.0.0.AND.ICON.EQ.1.AND.IHOLD.EQ.0 ) then
        A(N,N)=1.0
        IHOLD=N
        INFO=0
        end if
      RETURN
      END
      SUBROUTINE SGESL(A,LDA,N,IPVT,B,JOB)

c*********************************************************************72
c
cC SGESL SOLVES THE SYSTEM
c     A * X = B  OR  TRANS(A) * X = B
c     USING THE FACTORS COMPUTED BY SGECO OR SGEFA.
c
c     ON ENTRY
c
c        A       DIMENSION(LDA, N)
c                THE OUTPUT FROM SGECO OR SGEFA.
c
c        LDA     INTEGER
c                THE LEADING DIMENSION OF THE ARRAY  A .
c
c        N       INTEGER
c                THE ORDER OF THE MATRIX  A .
c
c        IPVT    INTEGER(N)
c                THE PIVOT VECTOR FROM SGECO OR SGEFA.
c
c        B       DIMENSION(N)
c                THE RIGHT HAND SIDE VECTOR.
c
c        JOB     INTEGER
c                = 0         TO SOLVE  A*X = B ,
c                = NONZERO   TO SOLVE  TRANS(A)*X = B  WHERE
c                            TRANS(A)  IS THE TRANSPOSE.
c
c     ON RETURN
c
c        B       THE SOLUTION VECTOR  X .
c
c     ERROR CONDITION
c
c        A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A
c        ZERO ON THE DIAGONAL.  TECHNICALLY THIS INDICATES SINGULARITY
c        BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER
c        SETTING OF LDA .  IT WILL NOT OCCUR IF THE SUBROUTINES ARE
c        CALLED CORRECTLY AND IF SGECO HAS SET RCOND .GT. 0.0
c        OR SGEFA HAS SET INFO .EQ. 0 .
c
c     TO COMPUTE  INVERSE(A) * C  WHERE  C  IS A MATRIX
c     WITH  P  COLUMNS
c           CALL SGECO(A,LDA,N,IPVT,RCOND,Z)
c           IF (RCOND IS TOO SMALL) GO TO ...
c           DO 10 J = 1, P
c              CALL SGESL(A,LDA,N,IPVT,C(1,J),0)
c        10 CONTINUE
c
c     LINPACK. THIS VERSION DATED 08/14/78 .
c     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
c
c     SUBROUTINES AND FUNCTIONS
c
c     BLAS SAXPY,SDOT
c
      INTEGER LDA,N,IPVT(1),JOB
      DIMENSION A(LDA,1),B(1)
      INTEGER K,KB,L,NM1
c
      NM1 = N - 1
      IF (JOB .NE. 0) GO TO 50
c
c        JOB = 0 , SOLVE  A * X = B
c        FIRST SOLVE  L*Y = B
c
         IF (NM1 .LT. 1) GO TO 30
         DO 20 K = 1, NM1
            L = IPVT(K)
            T = B(L)
            IF (L .EQ. K) GO TO 10
               B(L) = B(K)
               B(K) = T
   10       CONTINUE
            CALL SAXPY(N-K,T,A(K+1,K),1,B(K+1),1)
   20    CONTINUE
   30    CONTINUE
c
c        NOW SOLVE  U*X = Y
c
         DO 40 KB = 1, N
            K = N + 1 - KB
            B(K) = B(K)/A(K,K)
            T = -B(K)
            CALL SAXPY(K-1,T,A(1,K),1,B(1),1)
   40    CONTINUE
      GO TO 100
   50 CONTINUE
c
c        JOB = NONZERO, SOLVE  TRANS(A) * X = B
c        FIRST SOLVE  TRANS(U)*Y = B
c
         DO 60 K = 1, N
            T = SDOT(K-1,A(1,K),1,B(1),1)
            B(K) = (B(K) - T)/A(K,K)
   60    CONTINUE
c
c        NOW SOLVE TRANS(L)*X = Y
c
         IF (NM1 .LT. 1) GO TO 90
         DO 80 KB = 1, NM1
            K = N - KB
            B(K) = B(K) + SDOT(N-K,A(K+1,K),1,B(K+1),1)
            L = IPVT(K)
            IF (L .EQ. K) GO TO 70
               T = B(L)
               B(L) = B(K)
               B(K) = T
   70       CONTINUE
   80    CONTINUE
   90    CONTINUE
  100 CONTINUE
      RETURN
      END
      INTEGER FUNCTION ISAMAX(N,SX,INCX)

c*********************************************************************72
c
cC ISAMAX FINDS THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE.
c     JACK DONGARRA, LINPACK, 3/11/78.
c
      DIMENSION SX(1)
      INTEGER I,INCX,IX,N
c
      ISAMAX = 0
      IF( N .LT. 1 ) RETURN
      ISAMAX = 1
      IF(N.EQ.1)RETURN
      IF(INCX.EQ.1)GO TO 30
c
c        CODE FOR INCREMENT NOT EQUAL TO 1
c
      IX = 1
      SMAX = ABS(SX(1))
      IX = IX + INCX
      DO 20 I = 2,N
         IF(ABS(SX(IX)).LE.SMAX) GO TO 10
         ISAMAX = I
         SMAX = ABS(SX(IX))
   10    IX = IX + INCX
   20 CONTINUE
      RETURN
c
c        CODE FOR INCREMENT EQUAL TO 1
c
   30 SMAX = ABS(SX(1))
      DO 40 I = 2,N
         IF(ABS(SX(I)).LE.SMAX) GO TO 40
         ISAMAX = I
         SMAX = ABS(SX(I))
   40 CONTINUE
      RETURN
      END
      SUBROUTINE SAXPY(N,SA,SX,INCX,SY,INCY)

c*********************************************************************72
c
cC SAXPY CONSTANT TIMES A VECTOR PLUS A VECTOR.
c     USES UNROLLED LOOP FOR INCREMENTS EQUAL TO ONE.
c     JACK DONGARRA, LINPACK, 3/11/78.
c
      DIMENSION SX(1),SY(1)
      INTEGER I,INCX,INCY,IX,IY,M,MP1,N
c
      IF(N.LE.0)RETURN
      IF (SA .EQ. 0.0) RETURN
      IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
c
c        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
c          NOT EQUAL TO 1
c
      IX = 1
      IY = 1
      IF(INCX.LT.0)IX = (-N+1)*INCX + 1
      IF(INCY.LT.0)IY = (-N+1)*INCY + 1
      DO 10 I = 1,N
        SY(IY) = SY(IY) + SA*SX(IX)
        IX = IX + INCX
        IY = IY + INCY
   10 CONTINUE
      RETURN
c
c        CODE FOR BOTH INCREMENTS EQUAL TO 1
c
c
c        CLEAN-UP LOOP
c
   20 M = MOD(N,4)
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
        SY(I) = SY(I) + SA*SX(I)
   30 CONTINUE
      IF( N .LT. 4 ) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,4
        SY(I) = SY(I) + SA*SX(I)
        SY(I + 1) = SY(I + 1) + SA*SX(I + 1)
        SY(I + 2) = SY(I + 2) + SA*SX(I + 2)
        SY(I + 3) = SY(I + 3) + SA*SX(I + 3)
   50 CONTINUE
      RETURN
      END
      SUBROUTINE SCOPY(N,SX,INCX,SY,INCY)

c*********************************************************************72
c
cC SCOPY COPIES A VECTOR, X, TO A VECTOR, Y.
c     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO 1.
c     JACK DONGARRA, LINPACK, 3/11/78.
c
      DIMENSION SX(1),SY(1)
      INTEGER I,INCX,INCY,IX,IY,M,MP1,N
c
      IF(N.LE.0)RETURN
      IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
c
c        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
c          NOT EQUAL TO 1
c
      IX = 1
      IY = 1
      IF(INCX.LT.0)IX = (-N+1)*INCX + 1
      IF(INCY.LT.0)IY = (-N+1)*INCY + 1
      DO 10 I = 1,N
        SY(IY) = SX(IX)
        IX = IX + INCX
        IY = IY + INCY
   10 CONTINUE
      RETURN
c
c        CODE FOR BOTH INCREMENTS EQUAL TO 1
c
c
c        CLEAN-UP LOOP
c
   20 M = MOD(N,7)
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
        SY(I) = SX(I)
   30 CONTINUE
      IF( N .LT. 7 ) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,7
        SY(I) = SX(I)
        SY(I + 1) = SX(I + 1)
        SY(I + 2) = SX(I + 2)
        SY(I + 3) = SX(I + 3)
        SY(I + 4) = SX(I + 4)
        SY(I + 5) = SX(I + 5)
        SY(I + 6) = SX(I + 6)
   50 CONTINUE
      RETURN
      END
      FUNCTION SDOT(N,SX,INCX,SY,INCY)

c*********************************************************************72
c
cC SDOT FORMS THE DOT PRODUCT OF TWO VECTORS.
c     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
c     JACK DONGARRA, LINPACK, 3/11/78.
c
      DIMENSION SX(1),SY(1)
      INTEGER I,INCX,INCY,IX,IY,M,MP1,N
c
      STEMP = 0.0E0
      SDOT = 0.0E0
      IF(N.LE.0)RETURN
      IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
c
c        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
c          NOT EQUAL TO 1
c
      IX = 1
      IY = 1
      IF(INCX.LT.0)IX = (-N+1)*INCX + 1
      IF(INCY.LT.0)IY = (-N+1)*INCY + 1
      DO 10 I = 1,N
        STEMP = STEMP + SX(IX)*SY(IY)
        IX = IX + INCX
        IY = IY + INCY
   10 CONTINUE
      SDOT = STEMP
      RETURN
c
c        CODE FOR BOTH INCREMENTS EQUAL TO 1
c
c
c        CLEAN-UP LOOP
c
   20 M = MOD(N,5)
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
        STEMP = STEMP + SX(I)*SY(I)
   30 CONTINUE
      IF( N .LT. 5 ) GO TO 60
   40 MP1 = M + 1
      DO 50 I = MP1,N,5
        STEMP = STEMP + SX(I)*SY(I) + SX(I + 1)*SY(I + 1) +
     1   SX(I + 2)*SY(I + 2) + SX(I + 3)*SY(I + 3) + SX(I + 4)*SY(I + 4)
   50 CONTINUE
   60 SDOT = STEMP
      RETURN
      END
      SUBROUTINE SSCAL(N,SA,SX,INCX)

c*********************************************************************72
c
cC SSCAL SCALES A VECTOR BY A CONSTANT.
c     USES UNROLLED LOOPS FOR INCREMENT EQUAL TO 1.
c     JACK DONGARRA, LINPACK, 3/11/78.
c
      DIMENSION SX(1)
      INTEGER I,INCX,M,MP1,N,NINCX
c
      IF(N.LE.0)RETURN
      IF(INCX.EQ.1)GO TO 20
c
c        CODE FOR INCREMENT NOT EQUAL TO 1
c
      NINCX = N*INCX
      DO 10 I = 1,NINCX,INCX
        SX(I) = SA*SX(I)
   10 CONTINUE
      RETURN
c
c        CODE FOR INCREMENT EQUAL TO 1
c
c
c        CLEAN-UP LOOP
c
   20 M = MOD(N,5)
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
        SX(I) = SA*SX(I)
   30 CONTINUE
      IF( N .LT. 5 ) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,5
        SX(I) = SA*SX(I)
        SX(I + 1) = SA*SX(I + 1)
        SX(I + 2) = SA*SX(I + 2)
        SX(I + 3) = SA*SX(I + 3)
        SX(I + 4) = SA*SX(I + 4)
   50 CONTINUE
      RETURN
      END
      subroutine timestamp ( )

c*********************************************************************72
c
cc TIMESTAMP prints out the current YMDHMS date as a timestamp.
c
c  Discussion:
c
c    This FORTRAN77 version is made available for cases where the
c    FORTRAN90 version cannot be used.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    12 January 2007
c
c  Author:
c
c    John Burkardt
c
c  Parameters:
c
c    None
c
      implicit none

      character * ( 8 ) ampm
      integer d
      character * ( 8 ) date
      integer h
      integer m
      integer mm
      character * ( 9 ) month(12)
      integer n
      integer s
      character * ( 10 ) time
      integer y

      save month

      data month /
     &  'January  ', 'February ', 'March    ', 'April    ', 
     &  'May      ', 'June     ', 'July     ', 'August   ', 
     &  'September', 'October  ', 'November ', 'December ' /

      call date_and_time ( date, time )

      read ( date, '(i4,i2,i2)' ) y, m, d
      read ( time, '(i2,i2,i2,1x,i3)' ) h, n, s, mm

      if ( h .lt. 12 ) then
        ampm = 'AM'
      else if ( h .eq. 12 ) then
        if ( n .eq. 0 .and. s .eq. 0 ) then
          ampm = 'Noon'
        else
          ampm = 'PM'
        end if
      else
        h = h - 12
        if ( h .lt. 12 ) then
          ampm = 'PM'
        else if ( h .eq. 12 ) then
          if ( n .eq. 0 .and. s .eq. 0 ) then
            ampm = 'Midnight'
          else
            ampm = 'AM'
          end if
        end if
      end if

      write ( *, 
     &  '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) 
     &  d, month(m), y, h, ':', n, ':', s, '.', mm, ampm

      return
      end