NASA: Nonlocal Automated Sensitivity Analysis

Tracking Solutions of Parameterized Nonlinear
Systems Over Parameter Intervals

Last Updated: 22 April 2018

Site Maintained By:
Leigh Tesfatsion
Research Professor, and
Professor Emerita of Economics, Mathematics,
   and Electrical & Computer Engineering
Heady Hall 260
Iowa State University
Ames, Iowa 50011-1054
http://www2.econ.iastate.edu/tesfatsi/
tesfatsi AT iastate.edu

NASA Home Page

About NASA:

The Fortran program NASA (Nonlocal Automated Sensitivity Analysis), developed by Robert E. Kalaba and Leigh Tesfatsion, implements a complete system of ordinary differential equations (ODE) for the nonlocal sensitivity analysis of parameterized nonlinear systems using automated initialization, derivative evaluation, and solution generation procedures.

More precisely, suppose a nonlinear parameterized system H(x,b)=0 is given along with an initial value bStart for the scalar parameter b. The NASA program can then track the solution vector x(b) and sensitivity vector dx(b)/db, together with the adjoint A(b) and determinant D(b) of the Jacobian matrix J(b) = Hx(x(b),b), over any user-designated parameter path [bStart,bEnd] where the determinant D(b) remains nonzero.

This complete NASA ODE system was originally developed by Kalaba and Tesfatsion in "Nonlocal Automated Sensitivity Analysis" (pdf,1M), Computers and Mathematics With Applications 20 (1990), 53-65. This article constitutes the manual for the NASA program.

A summary report on the design and use of the NASA program is available online in both pdf (114K) and postscript (249K). A published version of this summary report can be found in C. A. Floudas and P. M. Pardalos (eds.), Encyclopedia of Optimization, Kluwer Academic Publishers, Vol. 4, 2001, 92-97.

For additional articles and introductory materials related to NASA, visit the NASA Home Page.

The NASA Program:

Software Release Disclaimer:

The software below is provided as-is, without warranty of any kind. It is released here as freeware by Robert E. Kalaba and Leigh Tesfatsion (the copyright holders) under the terms of the Artistic License Agreement.

// EXEC FORTVCLG                                                        00000010
//FORT.SYSIN DD *                                                       00000020
C                                                                       00000030
      PROGRAM  NASA                                                     00000040
C                                                                       00000050
C  NONLOCAL AUTOMATED SENSITIVITY ANALYSIS OF THE PARAMETERIZED SYSTEM  00000060
C  OF EQUATIONS PSI(X,ALPHA) = 0 [COMPUTERS AND MATHEMATICS WITH APPLI- 00000070
C  CATIONS, VOL. 20, NO. 2, PP. 53-65, 1990]                            00000080
C     BY ROBERT KALABA AND LEIGH TESFATSION                             00000090
C                                                                       00000100
C  LAST UPDATED: 5 DECEMBER 1989                                        00000110
C                                                                       00000120
      IMPLICIT COMPLEX (A-H,O-Z)                                        00000130
      REAL*4 VSS,TOL,RED                                                00000140
      DIMENSION  C(6),FVAL(6),Z(66),T(1110)                             00000150
      COMMON NVAR,NDIM,NSEC,N,NP1,NDIMT,MAXINT,INT,NDE                  00000160
      COMMON /FLAG/ C,FVAL,IFLGS,ALFA0                                  00000170
      COMMON /RSINIT/ VSS,NSPS,NP,MAXRS,TOL,NR,RED                      00000180
C                                                                       00000190
C  USER INPUTS FOR ALPHA INTEGRATION:                                   00000200
C                                                                       00000210
C  N IS THE NUMBER OF PSI COMPONENT FUNCTIONS, ALFA0 IS THE INITIAL     00000220
C  VALUE FOR THE PARAMETER ALPHA, H IS THE STEP SIZE FOR ALPHA, NSP     00000230
C  IS THE NUMBER OF ALPHA INTEGRATION STEPS TO BE TAKEN, NAP IS THE     00000240
C  NUMBER OF ALPHA INTEGRATION STEPS BETWEEN PRINTOUTS, MAXINT IS THE   00000250
C  TOTAL NUMBER OF INTEGRATION STEPS ALLOWED IN THE PROGRAM (INCLUDING  00000260
C  THE THETA-BETA INITIALIZATION PROCEDURE BELOW), C IS THE INITIAL     00000270
C  GUESS FOR X, AND NIT IS THE NUMBER OF THETA-BETA ITERATIONS TAKEN    00000280
C  TO IMPROVE ON C PRIOR TO THE ALPHA INTEGRATION.                      00000290
C                                                                       00000300
      N = 2                                                             00000310
      ALFA0 = CMPLX(0.5,0.0)                                            00000320
      H = CMPLX(0.001,0.0)                                              00000330
      NSP = 500                                                         00000340
      NAP = 50                                                          00000350
      MAXINT = 4000                                                     00000360
      C(1) = CMPLX(1.2,0.0)                                             00000370
      C(2) = CMPLX(1.1,0.0)                                             00000380
      NIT = 1                                                           00000390
C                                                                       00000400
C  N IS THE NUMBER OF COMPONENTS OF C AND FVAL WHICH ARE IN ACTUAL      00000410
C  USE, NDIM IS THE NUMBER OF COMPONENTS OF Z WHICH ARE IN ACTUAL       00000420
C  USE, AND NDIMT IS THE NUMBER OF COMPONENTS OF T WHICH ARE IN         00000430
C  ACTUAL USE. (T STORES THE COMPONENTS OF X, THE ADJOINT BY COLUMN,    00000440
C  AND THE DETERMINANT, PLUS THEIR DERIVATIVES AND WORKSPACE.)          00000450
      NDIM = 1 + (N+1) + (1 + (N+1))*(N+1)/2                            00000460
      NDIMT = 10*(1 + N  + N*N)                                         00000470
C                                                                       00000480
C  GETTING INITIAL CONDITIONS VIA A TWO-PHASE COMPLEX CONTINUATION:     00000490
C  FIRST SOLVE THE DERIVED SYSTEM OF EQUATIONS                          00000500
C     THETA*(PSI(X,ALFA0) - PSI(C,ALFA0)) + (1-THETA)*(X - C)  =  0     00000510
C  FOR THETA VARYING FROM 0 TO 1 TO GET ADJ AND DET OF THE JACOBIAN     00000520
C  OF PSI AT (C,ALFA0).  NEXT SOLVE THE DERIVED SYSTEM OF EQUATIONS     00000530
C     (PSI(X,ALFA0) - PSI(C,ALFA0)) + BETA*(PSI(C,ALFA0)) =  0          00000540
C  FOR BETA VARYING FROM 0 TO 1 TO GET APPROX SOLUTION X0 = X(ALFA0)    00000550
C  FOR PSI(X,ALFA0) = 0, PLUS THE ADJ AND DET OF THE JACOBIAN OF PSI    00000560
C  AT (X0,ALFA0).                                                       00000570
C                                                                       00000580
C  USER INPUTS FOR THE THETA AND BETA INTEGRATION PATHS THROUGH         00000590
C  THE COMPLEX SPIDER-WEB GRID:                                         00000600
C                                                                       00000610
C  VSS IS THE STEP SIZE ALONG SPOKES, NSPS IS THE TOTAL NUMBER OF       00000620
C  SPOKE-STEPS TO BE TAKEN (WITH VSS*NSPS = 1), NP IS THE NUMBER OF     00000630
C  SPOKE-STEPS BETWEEN DETAILED PRINTOUTS, AND MAXRS IS THE MAXIMUM     00000640
C  NUMBER OF ONE-DEGREE STEPS TO BE TAKEN ALONG ANY ONE RIM.            00000650
C                                                                       00000660
      VSS = 0.005                                                       00000670
      NSPS = 200                                                        00000680
      NP = 20                                                           00000690
      MAXRS = 360                                                       00000700
C                                                                       00000710
C  TOL IS THE MINIMUM TOLERABLE ABSOLUTE VALUE FOR THE JACOBIAN OF      00000720
C  THE DETERMINANT OF THE SYSTEM, NR IS THE NUMBER OF ONE-DEGREE RIM    00000730
C  STEPS ATTEMPTED BEFORE ANOTHER TOLERANCE TEST AFTER A TOLERANCE      00000740
C  TEST FAILURE ALONG A SPOKE, AND RED IS THE PERCENTAGE REDUCTION      00000750
C  FACTOR USED TO REDUCE TOL WHEN NR RIM-STEPS OR A SPOKE-STEP AWAY     00000760
C  FROM THE CURRENT POINT EACH RESULT IN A TOLERANCE TEST FAILURE, OR   00000770
C  THE NUMBER OF STEPS ALONG ANY ONE RIM REACHES MAXRS.                 00000780
C                                                                       00000790
      TOL = .007                                                        00000800
      NR  = 1                                                           00000810
      RED = .9                                                          00000820
C                                                                       00000830
C  THE REST OF THIS MAIN PROGRAM IS NOT PROBLEM-SPECIFIC.               00000840
C                                                                       00000850
C  OTHER VARIABLE DEFINITIONS AND INITIALIZATIONS                       00000860
      NP1 = N + 1                                                       00000870
      NDE = 1 + N + N*N                                                 00000880
      NVAR = NP1                                                        00000890
      NSEC = (NP1 +1)*NP1/2                                             00000900
      INT = 0                                                           00000910
C                                                                       00000920
C  THE VECTOR T WILL FIRST BE FILLED WITH THE INITIAL GUESS C FOR X,    00000930
C  THE ADJOINT I OF THE JACOBIAN OF (X - C) BY COLUMN, AND THE DET 1    00000940
C  OF THE JACOBIAN OF (X - C).                                          00000950
C                                                                       00000960
      LL = 1                                                            00000970
      DO 10 I = 1,N                                                     00000980
         T(LL) = C(LL)                                                  00000990
         LL = LL + 1                                                    00001000
   10 CONTINUE                                                          00001010
      DO 1000 IT = 1,NIT                                                00001020
      DO 20 J = 1,N                                                     00001030
         DO 30 I = 1,N                                                  00001040
            T(LL) = CMPLX(0.0,0.0)                                      00001050
            IF (J.EQ.I) T(LL) = CMPLX(1.0,0.0)                          00001060
            LL = LL + 1                                                 00001070
   30    CONTINUE                                                       00001080
   20 CONTINUE                                                          00001090
      T(LL) = CMPLX(1.0,0.0)                                            00001100
C                                                                       00001110
C  EVALUATION OF PSI(C,ALFA0)                                           00001120
      IFLGS = 2                                                         00001130
      ALFCUR = ALFA0                                                    00001140
      DO 40 K = 1,N                                                     00001150
           CALL FUN(K,C,ALFCUR,Z)                                       00001160
           FVAL(K) = Z(1)                                               00001170
   40 CONTINUE                                                          00001180
      WRITE (6,111)                                                     00001190
      WRITE (6,112) (C(K),K=1,N)                                        00001200
  111 FORMAT ('1',//1X,'HERE IS THE INITIAL GUESS C FOR X')             00001210
  112 FORMAT (1X,3X,4E12.3)                                             00001220
      WRITE (6,113)                                                     00001230
      WRITE (6,112) (FVAL(K),K=1,N)                                     00001240
  113 FORMAT (1X,'HERE ARE THE COMPONENT VALUES FOR PSI(C,ALFA0)')      00001250
C                                                                       00001260
C  THETA AND BETA INTEGRATIONS (IFLGS = 0,1)                            00001270
      IFLGS = 0                                                         00001280
      CALL INIT(T)                                                      00001290
      IFLGS = 1                                                         00001300
      CALL INIT(T)                                                      00001310
      IF (IT.LT.NIT) THEN                                               00001320
         LL = 1                                                         00001330
         DO 50 IJ = 1,N                                                 00001340
            C(LL) = T(LL)                                               00001350
            LL = LL + 1                                                 00001360
   50    CONTINUE                                                       00001370
      ENDIF                                                             00001380
 1000 CONTINUE                                                          00001390
C                                                                       00001400
C  ALPHA INTEGRATION (IFLGS = 2)                                        00001410
      IFLGS = 2                                                         00001420
      ALFCUR = ALFA0                                                    00001430
      BETA = CMPLX(1.0,0.0)                                             00001440
      THETA = CMPLX(1.0,0.0)                                            00001450
      CALL CINT1(T,NDE,ALFCUR,H)                                        00001460
      CALL OUTPUT(T,ALFCUR,BETA,THETA)                                  00001470
      JJ = NSP/NAP                                                      00001480
      DO 60 J = 1,JJ                                                    00001490
         DO 70 I = 1,NAP                                                00001500
            CALL CINT2(T,NDE,ALFCUR,H)                                  00001510
            INT = INT + 1                                               00001520
            IF (INT.GE.MAXINT) STOP 'ALPHA INTEGRATION'                 00001530
   70    CONTINUE                                                       00001540
         CALL OUTPUT(T,ALFCUR,BETA,THETA)                               00001550
   60 CONTINUE                                                          00001560
      WRITE (6,114) INT                                                 00001570
  114 FORMAT (1X,//1X,'THE TOTAL NUMBER OF INTEGRATION STEPS =',I7)     00001580
      STOP                                                              00001590
      END                                                               00001600
C                                                                       00001610
C                                                                       00001620
      SUBROUTINE PSIFUN(K,X,VCUR,VC,XX,PSIK)                            00001630
      IMPLICIT COMPLEX (A-H,O-Z)                                        00001640
      DIMENSION X1(66),X2(66),X3(66),X4(66),X5(66),X6(66)               00001650
      DIMENSION X(6),C(6),FVAL(6),VC(66),XX(6,66),PSIK(66),             00001660
     *ALFCUR(66)                                                        00001670
      DIMENSION U1(66),U2(66),U3(66),U4(66),U5(66)                      00001680
      COMMON NVAR,NDIM,NSEC,N,NP1,NDIMT,MAXINT,INT,NDE                  00001690
      COMMON /FLAG/ C,FVAL,IFLGS,ALFA0                                  00001700
C                                                                       00001710
C  THIS SUBROUTINE OBTAINS THE VALUE AND PARTIAL DERIVATIVES FOR THE    00001720
C  THE FOLLOWING PROBLEM-SPECIFIC COMPONENTS OF PSI:                    00001730
C                                                                       00001740
C         PSI1(X,ALPHA) = (1/2)(X1**(-1/2))*(X2**(1/3)) - ALPHA         00001750
C         PSI2(X,ALPHA) = (1/3)*(X1**(1/2))*(X2**(-2/3)) - (1/3)        00001760
C                                                                       00001770
C  NOT PROBLEM SPECIFIC:  THE INDEPENDENT VARIABLES X AND VCUR          00001780
C  ARE VECTORIZED AND ALPHA IS INITIALIZED                              00001790
C                                                                       00001800
      DO 1 I = 1,N                                                      00001810
         ZZ = X(I)                                                      00001820
         CALL VEC(I,ZZ,U1)                                              00001830
         DO 5 J = 1,NDIM                                                00001840
            XX(I,J) = U1(J)                                             00001850
    5    CONTINUE                                                       00001860
    1 CONTINUE                                                          00001870
      CALL VEC(NP1,VCUR,VC)                                             00001880
      IF (IFLGS.EQ.0.OR.IFLGS.EQ.1) THEN                                00001890
         CALL CON(ALFA0,U1)                                             00001900
         DO 10 J = 1,NDIM                                               00001910
            ALFCUR(J) = U1(J)                                           00001920
   10    CONTINUE                                                       00001930
      ELSEIF (IFLGS.EQ.2) THEN                                          00001940
         DO 15 J = 1,NDIM                                               00001950
            ALFCUR(J) = VC(J)                                           00001960
   15    CONTINUE                                                       00001970
      ENDIF                                                             00001980
C                                                                       00001990
C  PROBLEM-SPECIFIC STEPS START HERE                                    00002000
C                                                                       00002010
C  PUT THE N ROWS OF XX INTO VECTORS X1,...,XN                          00002020
      DO 20 J = 1,NDIM                                                  00002030
         X1(J) = XX(1,J)                                                00002040
         X2(J) = XX(2,J)                                                00002050
   20 CONTINUE                                                          00002060
C                                                                       00002070
C  SPECIFY ANY PROBLEM-SPECIFIC CONSTANTS USED TO DEFINE PSI            00002080
      P1 = CMPLX(0.5,0.0)                                               00002090
      P2 = CMPLX(-0.5,0.0)                                              00002100
      P3 = CMPLX(1.0,0.0)/CMPLX(3.0,0.0)                                00002110
      P4 = CMPLX(-2.0,0.0)/CMPLX(3.0,0.0)                               00002120
C                                                                       00002130
      IF (K.EQ.1) THEN                                                  00002140
C                                                                       00002150
C  GIVE THE PROBLEM-SPECIFIC LIST OF CALCULUS SUBROUTINES               00002160
C  FOR PSI1(X,ALFCUR) USING X = (X1,...,XN) AND PUT THE                 00002170
C  FINAL CALCULATION INTO PSIK                                          00002180
         CALL POW(X1,P2,U1)                                             00002190
         CALL CONVE(P1,U1,U2)                                           00002200
         CALL POW(X2,P3,U3)                                             00002210
         CALL MULT(U2,U3,U4)                                            00002220
         CALL SUB(U4,ALFCUR,PSIK)                                       00002230
         RETURN                                                         00002240
C                                                                       00002250
      ELSEIF (K.EQ.2) THEN                                              00002260
C                                                                       00002270
C  GIVE THE PROBLEM-SPECIFIC LIST OF CALCULUS SUBROUTINES               00002280
C  FOR PSI2(X,ALFCUR) USING X = (X1,...,XN) AND PUT THE                 00002290
C  FINAL CALCULATION INTO PSIK                                          00002300
         CALL POW(X1,P1,U1)                                             00002310
         CALL CONVE(P3,U1,U2)                                           00002320
         CALL POW(X2,P4,U3)                                             00002330
         CALL MULT(U2,U3,U4)                                            00002340
         CALL CON(P3,U5)                                                00002350
         CALL SUB(U4,U5,PSIK)                                           00002360
         RETURN                                                         00002370
C                                                                       00002380
      END IF                                                            00002390
      RETURN                                                            00002400
      END                                                               00002410
C                                                                       00002420
C                                                                       00002430
      SUBROUTINE OUTPUT(T,ALFCUR,BETA,THETA)                            00002440
      IMPLICIT COMPLEX (A-H,O-Z)                                        00002450
      DIMENSION T(1110),C(6),FVAL(6),X(6),Z(66)                         00002460
      DIMENSION PSI(6),PTHETA(6),PBETA(6)                               00002470
      COMMON NVAR,NDIM,NSEC,N,NP1,NDIMT,MAXINT,INT,NDE                  00002480
      COMMON /FLAG/ C,FVAL,IFLGS,ALFA0                                  00002490
      COMMON /RSINIT/ VSS,NSPS,NP,MAXRS,TOL,NR,RED                      00002500
C                                                                       00002510
C     THIS SUBROUTINE PRINTS OUTPUT FOR THE SYSTEM OF EQUATIONS         00002520
C     PSI(X,ALPHA) = 0, AND FOR THE THETA AND BETA CONTINUATIONS OF     00002530
C     PSI, WITH COMPONENT FUNCTIONS DEFINED AS IN SUBROUTINE PSIFUN     00002540
C     AND SUBROUTINE FUN.                                               00002550
C                                                                       00002560
      NDEP1 = NDE + 1                                                   00002570
      NDEPN = NDE + N                                                   00002580
      NDEMT = 2*NDE                                                     00002590
C                                                                       00002600
C  HERE ARE THE COMPONENTS OF PSI(X,ALFCUR)                             00002610
      DO 10 I = 1,N                                                     00002620
         X(I) = T(I)                                                    00002630
   10 CONTINUE                                                          00002640
      IFTEM = IFLGS                                                     00002650
      IFLGS = 2                                                         00002660
      DO 20 K = 1,N                                                     00002670
         CALL FUN(K,X,ALFCUR,Z)                                         00002680
         PSI(K) = Z(1)                                                  00002690
   20 CONTINUE                                                          00002700
      IFLGS = IFTEM                                                     00002710
C                                                                       00002720
C  HERE ARE THE COMPONENTS OF PTHETA(X,THETA) =                         00002730
C  THETA*(PSI(X,ALFA0) - PSI(C,ALFA0) + (1 - THETA)*(X - C)             00002740
      C1 = CMPLX(1.0,0.0)                                               00002750
      DO 30 K = 1,N                                                     00002760
         PTHETA(K) = THETA*(PSI(K)-FVAL(K))+(C1-THETA)*(T(K)-C(K))      00002770
   30 CONTINUE                                                          00002780
C                                                                       00002790
C  HERE ARE THE COMPONENTS OF PBETA(X,BETA) =                           00002800
C  (PSI(X,ALFA0) - PSI(C,ALFA0)) + BETA*PSI(C,ALFA0)                    00002810
      DO 40 K = 1,N                                                     00002820
         PBETA(K) = (PSI(K)-FVAL(K)) + BETA*FVAL(K)                     00002830
   40 CONTINUE                                                          00002840
C                                                                       00002850
      IF (IFLGS.EQ.0) THEN                                              00002860
         WRITE (6,100) THETA                                            00002870
  100    FORMAT (1X, //1X,'THETA =',2E12.3)                             00002880
         WRITE (6,110)                                                  00002890
  110    FORMAT (1X ,'HERE ARE THE COMPONENTS OF PTHETA')               00002900
         WRITE (6,120) (PTHETA(K),K=1,N)                                00002910
  120    FORMAT (1X,3X,4E18.6)                                          00002920
         WRITE (6,130)                                                  00002930
  130    FORMAT (1X,'HERE ARE THE COMPONENTS OF X AND THE SYSTEM DET')  00002940
         WRITE (6,120) (T(I),I=1,N) , T(NDE)                            00002950
         WRITE (6,140)                                                  00002960
  140    FORMAT (1X,'HERE ARE THE DERIVATIVES OF X AND THE DET ',       00002970
     *   'WRT THETA')                                                   00002980
         WRITE (6,120) (T(J),J=NDEP1,NDEPN), T(NDEMT)                   00002990
      ELSEIF (IFLGS.EQ.1) THEN                                          00003000
         WRITE (6,200) BETA                                             00003010
  200    FORMAT (1X,//1X,'BETA =',2E12.3)                               00003020
         WRITE (6,210)                                                  00003030
  210    FORMAT (1X,'HERE ARE THE COMPONENTS OF PBETA')                 00003040
         WRITE (6,120) (PBETA(K),K=1,N)                                 00003050
         WRITE (6,130)                                                  00003060
         WRITE (6,120) (T(I),I=1,N), T(NDE)                             00003070
         WRITE (6,240)                                                  00003080
  240    FORMAT (1X,'HERE ARE THE DERIVATIVES OF X AND THE DET ',       00003090
     *   'WRT BETA')                                                    00003100
         WRITE (6,120) (T(J),J=NDEP1,NDEPN),T(NDEMT)                    00003110
      ELSEIF (IFLGS.EQ.2) THEN                                          00003120
         WRITE (6,300) ALFCUR                                           00003130
  300    FORMAT (1X, //1X,'ALFA =',2E12.3)                              00003140
         WRITE (6,310)                                                  00003150
  310    FORMAT (1X, 'HERE ARE THE COMPONENTS OF PSI')                  00003160
         WRITE (6,120) (PSI(K),K=1,N)                                   00003170
         WRITE (6,130)                                                  00003180
         WRITE (6,120) (T(I),I=1,N), T(NDE)                             00003190
         WRITE (6,340)                                                  00003200
  340    FORMAT (1X,'HERE ARE THE DERIVATIVES OF X AND THE DET ',       00003210
     *   'WRT ALFA')                                                    00003220
         WRITE (6,120) (T(J),J=NDEP1,NDEPN), T(NDEMT)                   00003230
      ENDIF                                                             00003240
      RETURN                                                            00003250
      END                                                               00003260
C                                                                       00003270
      SUBROUTINE FUN(K,X,VCUR,Z)                                        00003280
      IMPLICIT COMPLEX (A-H,O-Z)                                        00003290
      DIMENSION X(6),Z(66),C(6),FVAL(6),VC(66),XX(6,66),PSIK(66),       00003300
     *XK(66)                                                            00003310
      DIMENSION U1(66),U2(66),U3(66),U4(66),U5(66),U6(66),U7(66),       00003320
     *U8(66)                                                            00003330
      COMMON NVAR,NDIM,NSEC,N,NP1,NDIMT,MAXINT,INT,NDE                  00003340
      COMMON /FLAG/ C,FVAL,IFLGS,ALFA0                                  00003350
C                                                                       00003360
C  THIS SUBROUTINE OBTAINS THE VALUE AND PARTIAL DERIVATIVES FOR        00003370
C  THE THETA AND BETA CONTINUATIONS PTHETA AND PBETA FOR ANY GIVEN      00003380
C  PROBLEM-SPECIFIC FUNCTION PSI GIVEN IN SUBROUTINE PSIFUN, WHERE      00003390
C  PTHETAK(X,THETA) = THETA*(PSIK(X,ALFA0)-PSIK(C,ALFA0))               00003400
C                     + (1 - THETA)*(XK - CK)                           00003410
C  PBETAK(X,BETA) = (PSIK(X,ALFA0)-PSIK(C,ALFA0) + BETA*PSIK(C,ALFA0)   00003420
C                                                                       00003430
      CALL PSIFUN(K,X,VCUR,VC,XX,PSIK)                                  00003440
      IF (IFLGS.EQ.2) THEN                                              00003450
         DO 1111 J = 1,NDIM                                             00003460
            Z(J) = PSIK(J)                                              00003470
 1111    CONTINUE                                                       00003480
         RETURN                                                         00003490
      ENDIF                                                             00003500
C                                                                       00003510
      FVK = FVAL(K)                                                     00003520
C                                                                       00003530
      IF (IFLGS.EQ.0) THEN                                              00003540
C                                                                       00003550
      B1 = CMPLX(1.0,0.0)                                               00003560
      CK = C(K)                                                         00003570
      DO 2222 I = 1,NDIM                                                00003580
        XK(I) = XX(K,I)                                                 00003590
 2222 CONTINUE                                                          00003600
C                                                                       00003610
C  HERE IS THETA*(PSIK(X,ALFA0)-PSIK(C,ALFA0)) + (1-THETA)*(XK-CK)      00003620
         CALL CON(FVK,U1)                                               00003630
         CALL SUB(PSIK,U1,U2)                                           00003640
         CALL MULT(VC,U2,U3)                                            00003650
         CALL CON(B1,U4)                                                00003660
         CALL SUB(U4,VC,U5)                                             00003670
         CALL CON(CK,U6)                                                00003680
         CALL SUB(XK,U6,U7)                                             00003690
         CALL MULT(U5,U7,U8)                                            00003700
         CALL ADD(U3,U8,Z)                                              00003710
         RETURN                                                         00003720
C                                                                       00003730
      ELSEIF (IFLGS.EQ.1) THEN                                          00003740
C                                                                       00003750
C  HERE IS (PSIK(X,ALFA0)-PSIK(C,ALFA0)) + BETA*PSIK(C,ALFA0)           00003760
         CALL CON(FVK,U1)                                               00003770
         CALL MULT(VC,U1,U2)                                            00003780
         CALL ADD(PSIK,U2,U3)                                           00003790
         CALL SUB(U3,U1,Z)                                              00003800
         RETURN                                                         00003810
C                                                                       00003820
      ENDIF                                                             00003830
      RETURN                                                            00003840
      END                                                               00003850
C                                                                       00003860
C                                                                       00003870
      SUBROUTINE INIT(T)                                                00003880
      IMPLICIT COMPLEX (A-H,O-Z)                                        00003890
      REAL*4 VSS,TOL,RED,RADIUS,ANGLE,DEGREE,PI,SADET,GSADET            00003900
      REAL*4 RNUM,MNUM,UANGLE,DANGLE,GANGLE,UADET,DADET                 00003910
      DIMENSION T(1110),C(6),FVAL(6),STNEW(1110),GSTNEW(1110)           00003920
      DIMENSION UTNEW(1110),DTNEW(1110)                                 00003930
      COMMON NVAR,NDIM,NSEC,N,NP1,NDIMT,MAXINT,INT,NDE                  00003940
      COMMON /FLAG/ C,FVAL,IFLGS,ALFA0                                  00003950
      COMMON /RSINIT/ VSS,NSPS,NP,MAXRS,TOL,NR,RED                      00003960
C                                                                       00003970
C  THIS SUBROUTINE CARRIES OUT THE INTEGRATIONS FOR THE THETA AND       00003980
C  BETA CONTINUATIONS TO OBTAIN THE INITIAL CONDITIONS NEEDED FOR       00003990
C  THE ALPHA INTEGRATION STARTING FROM ALPHA = ALFA0.                   00004000
C                                                                       00004010
      PI = 4.0*ATAN(1.0)                                                00004020
      DEGREE = PI/180.0                                                 00004030
C                                                                       00004040
C  VARIABLE INITIALIZATIONS                                             00004050
      NRIM = 0                                                          00004060
      IFDIR = 0                                                         00004070
      NRSTEP = 0                                                        00004080
      NSPOKE = 0                                                        00004090
      RADIUS = 1.0                                                      00004100
      ANGLE = 0.0                                                       00004110
      GSADET = 0.0                                                      00004120
      GANGLE = ANGLE                                                    00004130
C                                                                       00004140
      IF (IFLGS.EQ.0) THEN                                              00004150
         THETA = CMPLX(0.0,0.0)                                         00004160
         BETA = CMPLX(0.0,0.0)                                          00004170
         VCUR = THETA                                                   00004180
      ELSEIF (IFLGS.EQ.1) THEN                                          00004190
         DO 1 K = 1,N                                                   00004200
            T(K) = C(K)                                                 00004210
    1    CONTINUE                                                       00004220
         THETA = CMPLX(1.0,0.0)                                         00004230
         BETA = CMPLX(0.0,0.0)                                          00004240
         VCUR = BETA                                                    00004250
      ENDIF                                                             00004260
      HS = CMPLX(VSS,0.0)                                               00004270
      CALL CINT1(T,NDE,VCUR,HS)                                         00004280
      GSHS = HS                                                         00004290
      CALL OUTPUT(T,ALFA0,BETA,THETA)                                   00004300
      IF (IFLGS.EQ.0) THEN                                              00004310
         WRITE (6,211) NRIM,NSPOKE,NRSTEP,VCUR                          00004320
      ELSEIF (IFLGS.EQ.1) THEN                                          00004330
         WRITE (6,221) NRIM,NSPOKE,NRSTEP,VCUR                          00004340
      ENDIF                                                             00004350
      WRITE (6,231) TOL,T(NDE),RADIUS                                   00004360
C                                                                       00004370
C  FIRST TEST TO SEE WHETHER A STEP CAN BE MADE INWARD ALONG THE        00004380
C  CURRENT SPOKE AT THE CURRENT TOLERANCE LEVEL TOL, I.E., WHETHER      00004390
C  THE ABSOLUTE VALUE SADET OF THE DETERMINANT OF THE JACOBIAN          00004400
C  IS AT LEAST AS GREAT AS TOL WHEN THIS SPOKE STEP IS TAKEN.           00004410
C                                                                       00004420
 2000 IF (NRIM.LT.NSPS) THEN                                            00004430
         DO 10 IM = 1,NDIMT                                             00004440
            STNEW(IM) = T(IM)                                           00004450
   10    CONTINUE                                                       00004460
         SVNEW = VCUR                                                   00004470
         CALL CINT2(STNEW,NDE,SVNEW,HS)                                 00004480
         INT = INT + 1                                                  00004490
         IF (INT.GE.MAXINT) STOP 'INITIALIZATION 10'                    00004500
         DET = STNEW(NDE)                                               00004510
         SADET = CABS(DET)                                              00004520
C  IF SADET IS GREATER OR EQUAL TO TOL, MOVE ALONG THE CURRENT SPOKE    00004530
C  TO THE NEXT RIM.                                                     00004540
        IF (SADET.GE.TOL) THEN                                          00004550
          DO 20 ID = 1,NDIMT                                            00004560
            T(ID) = STNEW(ID)                                           00004570
   20     CONTINUE                                                      00004580
          VCUR = SVNEW                                                  00004590
          RADIUS = RADIUS - VSS                                         00004600
          NRIM = NRIM + 1                                               00004610
          IFDIR = 0                                                     00004620
          NRSTEP = 0                                                    00004630
          GSADET = 0.0                                                  00004640
          GSHS = HS                                                     00004650
          GANGLE = ANGLE                                                00004660
        ENDIF                                                           00004670
C  IF SADET IS LESS THAN TOL, AND THE NUMBER OF RIM-STEPS PREVIOUSLY    00004680
C  TAKEN ALONG THE CURRENT RIM IS LESS THAN MAXRS, TEST TO SEE IF NR    00004690
C  RIM-STEPS EITHER UPWARD OR DOWNWARD ALONG THE RIM ARE POSSIBLE AT    00004700
C  THE CURRENT TOLERANCE LEVEL.  IF NOT, MOVE TO THE TESTED SPOKE OR    00004710
C  OR RIM POINT WITH THE GREATEST ABSOLUTE VALUE ABS(DET) FOR THE       00004720
C  DET OF THE JACOBIAN AND REDUCE TOL TO RED*(GREATEST ABS(DET)).       00004730
        IF (SADET.LT.TOL.AND.NRSTEP.LT.MAXRS) THEN                      00004740
          IF (SADET.GT.GSADET) THEN                                     00004750
             GSADET = SADET                                             00004760
             GSHS = HS                                                  00004770
             GANGLE = ANGLE                                             00004780
             DO 30 LT = 1,NDIMT                                         00004790
                GSTNEW(LT) = STNEW(LT)                                  00004800
   30        CONTINUE                                                   00004810
             GSVNEW = SVNEW                                             00004820
          ENDIF                                                         00004830
          UADET = 0.0                                                   00004840
          IF (IFDIR.GE.0) THEN                                          00004850
            DO 40 LK = 1,NDIMT                                          00004860
                UTNEW(LK) = T(LK)                                       00004870
   40       CONTINUE                                                    00004880
            UVNEW = VCUR                                                00004890
            UANGLE = ANGLE                                              00004900
            DO 50 JJJ = 1,NR                                            00004910
             UANGLE = UANGLE + DEGREE                                   00004920
             RNUM = 1.0 - RADIUS*COS(UANGLE)                            00004930
             MNUM = RADIUS*SIN(UANGLE)                                  00004940
             UHS = CMPLX(RNUM,MNUM) - UVNEW                             00004950
             CALL CINT1(UTNEW,NDE,UVNEW,UHS)                            00004960
             CALL CINT2(UTNEW,NDE,UVNEW,UHS)                            00004970
             INT = INT + 1                                              00004980
             IF (INT.GE.MAXINT) STOP 'INITIALIZATION 50'                00004990
   50       CONTINUE                                                    00005000
            DET = UTNEW(NDE)                                            00005010
            UADET = CABS(DET)                                           00005020
          ENDIF                                                         00005030
          DADET = 0.0                                                   00005040
          IF (IFDIR.LE.0) THEN                                          00005050
            DO 60 IK = 1,NDIMT                                          00005060
             DTNEW(IK) = T(IK)                                          00005070
   60       CONTINUE                                                    00005080
            DVNEW = VCUR                                                00005090
            DANGLE = ANGLE                                              00005100
            DO 70 III = 1,NR                                            00005110
             DANGLE = DANGLE - DEGREE                                   00005120
             RNUM = 1.0 - RADIUS*COS(DANGLE)                            00005130
             MNUM = RADIUS*SIN(DANGLE)                                  00005140
             DHS = CMPLX(RNUM,MNUM) - DVNEW                             00005150
             CALL CINT1(DTNEW,NDE,DVNEW,DHS)                            00005160
             CALL CINT2(DTNEW,NDE,DVNEW,DHS)                            00005170
             INT = INT + 1                                              00005180
             IF (INT.GE.MAXINT) STOP 'INITIALIZATION 70'                00005190
   70       CONTINUE                                                    00005200
            DET = DTNEW(NDE)                                            00005210
            DADET = CABS(DET)                                           00005220
          ENDIF                                                         00005230
          IF (UADET.GT.SADET.AND.UADET.GE.DADET) THEN                   00005240
             IF (UADET.LT.TOL) TOL = RED*UADET                          00005250
             DO 80 LLK = 1,NDIMT                                        00005260
                T(LLK) = UTNEW(LLK)                                     00005270
   80        CONTINUE                                                   00005280
             VCUR = UVNEW                                               00005290
             ANGLE = UANGLE                                             00005300
             NRSTEP = NRSTEP + NR                                       00005310
             NSPOKE = NSPOKE + NR                                       00005320
             IFDIR = 1                                                  00005330
             RNUM = VSS*COS(ANGLE)                                      00005340
             MNUM = -VSS*SIN(ANGLE)                                     00005350
             HS = CMPLX(RNUM,MNUM)                                      00005360
             CALL CINT1(T,NDE,VCUR,HS)                                  00005370
          ENDIF                                                         00005380
          IF (DADET.GT.SADET.AND.DADET.GE.UADET) THEN                   00005390
             IF (DADET.LT.TOL) TOL = RED*DADET                          00005400
             DO 90 IIK = 1,NDIMT                                        00005410
                 T(IIK) = DTNEW(IIK)                                    00005420
   90        CONTINUE                                                   00005430
             VCUR = DVNEW                                               00005440
             ANGLE = DANGLE                                             00005450
             NRSTEP = NRSTEP + NR                                       00005460
             NSPOKE = NSPOKE - NR                                       00005470
             IFDIR = -1                                                 00005480
             RNUM = VSS*COS(ANGLE)                                      00005490
             MNUM = -VSS*SIN(ANGLE)                                     00005500
             HS = CMPLX(RNUM,MNUM)                                      00005510
             CALL CINT1(T,NDE,VCUR,HS)                                  00005520
          ENDIF                                                         00005530
          IF (SADET.GE.DADET.AND.SADET.GE.UADET) THEN                   00005540
             TOL = RED*SADET                                            00005550
             DO 100 ILK = 1,NDIMT                                       00005560
                T(ILK) = STNEW(ILK)                                     00005570
  100        CONTINUE                                                   00005580
             VCUR = SVNEW                                               00005590
             RADIUS = RADIUS - VSS                                      00005600
             NRIM = NRIM + 1                                            00005610
             IFDIR = 0                                                  00005620
             NRSTEP = 0                                                 00005630
             GSADET = 0.0                                               00005640
             GSHS = HS                                                  00005650
             GANGLE = ANGLE                                             00005660
          ENDIF                                                         00005670
        ENDIF                                                           00005680
C  IF MAXRS STEPS HAVE BEEN TAKEN AROUND ANY ONE RIM, DECREASE TOL      00005690
C  TO RED*(GREATEST ABSOLUTE VALUE OF DET ALONG NEXT INNER RIM)         00005700
C  AND MOVE TO THIS NEXT INNER RIM POINT.                               00005710
        IF (SADET.LT.TOL.AND.NRSTEP.GE.MAXRS) THEN                      00005720
          TOL = RED*GSADET                                              00005730
          DO 110 KKL = 1,NDIMT                                          00005740
             T(KKL) = GSTNEW(KKL)                                       00005750
  110     CONTINUE                                                      00005760
          VCUR = GSVNEW                                                 00005770
          RADIUS = RADIUS - VSS                                         00005780
          ANGLE = GANGLE                                                00005790
          NRIM = NRIM + 1                                               00005800
          IFDIR = 0                                                     00005810
          NRSTEP = 0                                                    00005820
          GSADET = 0.0                                                  00005830
          HS = GSHS                                                     00005840
          CALL CINT1(T,NDE,VCUR,HS)                                     00005850
        ENDIF                                                           00005860
        IF ((NRIM/NP)*NP.EQ.NRIM) THEN                                  00005870
          IF (IFLGS.EQ.0) THEN                                          00005880
             CALL OUTPUT(T,ALFA0,BETA,VCUR)                             00005890
             WRITE (6,211) NRIM,NSPOKE,NRSTEP,VCUR                      00005900
  211        FORMAT(1X, /1X,'NRIM =',I6,2X,'NSPOKE =',I6,2X,            00005910
     *       'NRSTEP =',I6,2X,'THETA =',2E12.3)                         00005920
          ELSEIF (IFLGS.EQ.1) THEN                                      00005930
             CALL OUTPUT(T,ALFA0,VCUR,THETA)                            00005940
             WRITE (6,221) NRIM,NSPOKE,NRSTEP,VCUR                      00005950
  221        FORMAT(1X, /1X,'NRIM =',I6,2X,'NSPOKE =',I6,2X,            00005960
     *       'NRSTEP =',I6,2X,'BETA =',2E12.3)                          00005970
          ENDIF                                                         00005980
          WRITE (6,231) TOL,T(NDE),RADIUS                               00005990
  231     FORMAT (1X, 'TOL =',E12.3,3X,'DET =',2E12.3,3X,'RADIUS =',    00006000
     *    F10.3)                                                        00006010
        ENDIF                                                           00006020
        GO TO 2000                                                      00006030
      ENDIF                                                             00006040
      WRITE (6,241) INT                                                 00006050
  241 FORMAT (1X,//1X, 'NUMBER OF INTEGRATION STEPS SO FAR =', I6)      00006060
      RETURN                                                            00006070
      END                                                               00006080
C                                                                       00006090
C                                                                       00006100
      SUBROUTINE DAUX(T,VCUR,H)                                         00006110
C                                                                       00006120
C     THIS DERIVATIVE AUXILIARY ROUTINE CALCULATES THE DERIVATIVES      00006130
C     DX/DVCUR, DADJ/DVCUR, AND DDET/DVCUR, AND STORES THEM IN T.       00006140
C                                                                       00006150
      IMPLICIT COMPLEX (A-H,O-Z)                                        00006160
      DIMENSION T(1110), A(10,10),PSUBA(10),V(10),DXSUBA(10),B(10,10),  00006170
     *D(10,10),E(10,10),F(10,10),ASUBA(10,10),X(6),G(10,10),            00006180
     *Z(66),DPSI(10,11),DDPSI(10,11,11),VF(10)                          00006190
      COMMON NVAR,NDIM,NSEC,N,NP1,NDIMT,MAXINT,INT,NDE                  00006200
C                                                                       00006210
C     DICTIONARY OF VALUES T(1) THROUGH T(NDE)                          00006220
C                                                                       00006230
      L = 1                                                             00006240
      DO 10 KJ=1,N                                                      00006250
      X(KJ)=T(L)                                                        00006260
      L=L+1                                                             00006270
   10 CONTINUE                                                          00006280
      DO 20 J=1,N                                                       00006290
      DO 30 I=1,N                                                       00006300
      A(I,J)=T(L)                                                       00006310
      L=L+1                                                             00006320
   30 CONTINUE                                                          00006330
   20 CONTINUE                                                          00006340
      DET=T(L)                                                          00006350
      L = L+1                                                           00006360
C                                                                       00006370
C  CALCULATE FIRST AND SECOND-ORDER PARTIAL DERIVATIVES OF PSI AND      00006380
C  STORE THEM IN DPSI AND DDPSI.                                        00006390
C                                                                       00006400
      DO 31 K=1,N                                                       00006410
      LL=1                                                              00006420
      CALL FUN(K,X,VCUR,Z)                                              00006430
      VF(K)=Z(LL)                                                       00006440
      LL=LL+1                                                           00006450
      DO 41 I=1,NP1                                                     00006460
      DPSI(K,I)=Z(LL)                                                   00006470
      LL=LL+1                                                           00006480
   41 CONTINUE                                                          00006490
      DO 51 I=1,NP1                                                     00006500
      DO 61 J=I,NP1                                                     00006510
      DDPSI(K,I,J)=Z(LL)                                                00006520
      DDPSI(K,J,I)=Z(LL)                                                00006530
      LL=LL+1                                                           00006540
   61 CONTINUE                                                          00006550
   51 CONTINUE                                                          00006560
   31 CONTINUE                                                          00006570
C                                                                       00006580
C  CALCULATE DX/DVCUR = -(ADJ/DET)*(DPSI/DVCUR)                         00006590
C                                                                       00006600
      DO 40 K=1,N                                                       00006610
      PSUBA(K)=DPSI(K,NP1)                                              00006620
   40 CONTINUE                                                          00006630
      CALL MATVEC(A,PSUBA,N,V)                                          00006640
      CC = CMPLX(-1.0,0.0)/DET                                          00006650
      CALL VECON(CC,V,N,DXSUBA)                                         00006660
C                                                                       00006670
C  CALCULATE B = D(JACOBIAN)/DVCUR                                      00006680
C                                                                       00006690
      DO 50 I=1,N                                                       00006700
      DO 60 J=1,N                                                       00006710
      SUM=DDPSI(I,J,NP1)                                                00006720
      DO 70 M=1,N                                                       00006730
      SUM=SUM+DDPSI(I,J,M)*DXSUBA(M)                                    00006740
   70 CONTINUE                                                          00006750
      B(I,J)=SUM                                                        00006760
   60 CONTINUE                                                          00006770
   50 CONTINUE                                                          00006780
C                                                                       00006790
C  CALCULATE DADJ/DVCUR = (A*TR(A*B) - A*B*A)/DET                       00006800
C                                                                       00006810
      CALL MULMAT(A,B,N,D)                                              00006820
      CALL TRACE(D,N,TR)                                                00006830
      CALL MULMAT(D,A,N,E)                                              00006840
      CALL CONMAT(TR,A,N,F)                                             00006850
      CALL SUBMAT(F,E,N,G)                                              00006860
      CC = -CC                                                          00006870
      CALL CONMAT(CC,G,N,ASUBA)                                         00006880
C                                                                       00006890
C  CALCULATE DDET/DVCUR = TR(A*B)                                       00006900
C                                                                       00006910
      DDETA=TR                                                          00006920
C                                                                       00006930
C  MAP THE DERIVATIVES DX/DVCUR, DADJ/DVCUR (BY COLUMN), AND            00006940
C  DDET/DVCUR INTO T IN ORDER, STARTING WITH T(NDE + 1).                00006950
C                                                                       00006960
      DO 80 J=1,N                                                       00006970
      T(L)=DXSUBA(J)                                                    00006980
      L=L+1                                                             00006990
   80 CONTINUE                                                          00007000
      DO 81 J=1,N                                                       00007010
      DO 82 I=1,N                                                       00007020
      T(L)=ASUBA(I,J)                                                   00007030
      L=L+1                                                             00007040
   82 CONTINUE                                                          00007050
   81 CONTINUE                                                          00007060
      T(L)=DDETA                                                        00007070
      RETURN                                                            00007080
      END                                                               00007090
C                                                                       00007100
C                                                                       00007110
      SUBROUTINE CINT1(T,N,X,H)                                         00007120
C     CALCULATES CONSTANTS NEEDED FOR INTEGRATION BY SUB CINT2          00007130
      IMPLICIT COMPLEX (A-H,O-Z)                                        00007140
      DIMENSION T(1110)                                                 00007150
      COMMON /NSEQ/ N2,N3,N4,N5,N6,N7,N8,N9,NN,KFLAG,IND                00007160
      COMMON /HSEQ/ H2,H4,H24,R                                         00007170
      NN=N                                                              00007180
      N2=N*2                                                            00007190
      N3=N*3                                                            00007200
      N4=N*4                                                            00007210
      N5=N*5                                                            00007220
      N6=N*6                                                            00007230
      N7=N*7                                                            00007240
      N8=N*8                                                            00007250
      N9=N*9                                                            00007260
      H2=H*CMPLX(0.5,0.0)                                               00007270
      H4=H2*CMPLX(0.5,0.0)                                              00007280
      H24=H/CMPLX(24.0,0.0)                                             00007290
      R=CMPLX(1.0,0.0)/CMPLX(6.0,0.0)                                   00007300
C     CALCULATE Y PRIME FOR INITIAL CONDITIONS                          00007310
      CALL DAUX(T,X,H)                                                  00007320
      DO 1 I=1,N                                                        00007330
      N9I=I+N9                                                          00007340
C     TEMPORARY STORAGE FOR Y                                           00007350
      N8I=I+N8                                                          00007360
      NNI=I+NN                                                          00007370
      T(N9I)=T(I)                                                       00007380
C     STORE Y PRIME AT N-3 FOR USE IN ADAMS-MOULTON INTEGRATION         00007390
    1 T(N8I)=T(NNI)                                                     00007400
      KFLAG=0                                                           00007410
      IND=0                                                             00007420
      RETURN                                                            00007430
      END                                                               00007440
C                                                                       00007450
C                                                                       00007460
      SUBROUTINE CINT2(T,N,X,H)                                         00007470
C     FOURTH-ORDER ADAMS-MOULTON INTEGRATION WITH A RUNGE-KUTTA START.  00007480
      IMPLICIT COMPLEX (A-H,O-Z)                                        00007490
      DIMENSION T(1110)                                                 00007500
      COMMON /NSEQ/ N2,N3,N4,N5,N6,N7,N8,N9,NN,KFLAG,IND                00007510
      COMMON /HSEQ/ H2,H4,H24,R                                         00007520
C                                                                       00007530
C     IND=FLAG FOR RUNGE-KUTTA INTEGRATION WHEN EQUAL OR LESS THAN 3.   00007540
C     KFLAG=FLAG FOR STORING PAST DERIVATIVES FOR USE IN ADAMS-MOULTON  00007550
C     INTEGRATION                                                       00007560
C       RUNGE-KUTTA INTEGRATION (6 STEPS AT H/2)                        00007570
C     NOTE-TWO STEPS OF RUNGE-KUTTA ARE DONE WITH EACH CALL TO THIS     00007580
C     SUBROUTINE, SO THAT PRINTOUT POINTS WILL BE AT STEPS OF H.        00007590
C                                                                       00007600
      IND=IND+1                                                         00007610
      KFLAG=KFLAG+1                                                     00007620
      IF (IND-3) 20,20,12                                               00007630
   20 DO 9 K=1,2                                                        00007640
C     CALCULATE K1                                                      00007650
      DO 1 I=1,N                                                        00007660
      N2I=I+N2                                                          00007670
      NNI=I+NN                                                          00007680
    1 T(N2I)=T(NNI)*H2                                                  00007690
C     CALCULATE K2                                                      00007700
C     STEP UP X                                                         00007710
      X=X+H4                                                            00007720
      DO 2 I=1,N                                                        00007730
      N9I=I+N9                                                          00007740
      N2I=I+N2                                                          00007750
    2 T(I)=T(N9I)+CMPLX(0.5,0.0)*T(N2I)                                 00007760
      CALL DAUX(T,X,H)                                                  00007770
C     STORE K2                                                          00007780
      DO 3 I=1,N                                                        00007790
      N3I=I+N3                                                          00007800
      NNI=I+NN                                                          00007810
    3 T(N3I)=T(NNI)*H2                                                  00007820
C     CALCULATE K3                                                      00007830
      DO 4 I=1,N                                                        00007840
      N9I=I+N9                                                          00007850
      N3I=I+N3                                                          00007860
    4 T(I)=T(N9I)+CMPLX(0.5,0.0)*T(N3I)                                 00007870
      CALL DAUX(T,X,H)                                                  00007880
C     STORE K3                                                          00007890
      DO 5 I=1,N                                                        00007900
      N4I=I+N4                                                          00007910
      NNI=I+NN                                                          00007920
    5 T(N4I)=T(NNI)*H2                                                  00007930
C     CALCULATE K4                                                      00007940
C     STEP UP X                                                         00007950
      X=X+H4                                                            00007960
      DO 6 I=1,N                                                        00007970
      N9I=I+N9                                                          00007980
      N4I=I+N4                                                          00007990
    6 T(I)=T(N9I)+T(N4I)                                                00008000
      CALL DAUX(T,X,H)                                                  00008010
C     STORE K4                                                          00008020
      DO 7 I=1,N                                                        00008030
      N5I=I+N5                                                          00008040
      NNI=I+NN                                                          00008050
    7 T(N5I)=T(NNI)*H2                                                  00008060
C     CALCULATE PREDICTED VALUE OF Y                                    00008070
      DO 8 I=1,N                                                        00008080
      N9I=I+N9                                                          00008090
      N2I=I+N2                                                          00008100
      N3I=I+N3                                                          00008110
      N4I=I+N4                                                          00008120
      N5I=I+N5                                                          00008130
      Y=T(N9I)+R*(T(N2I)+CMPLX(2.0,0.0)*T(N3I)                          00008140
     *+CMPLX(2.0,0.0)*T(N4I)+T(N5I))                                    00008150
C     STORE AS CURRENT VALUE OF Y                                       00008160
      T(I)=Y                                                            00008170
C     STORE Y IN TEMPORARY STORAGE                                      00008180
    8 T(N9I)=T(I)                                                       00008190
C     CALCULATE Y PRIME                                                 00008200
      CALL DAUX(T,X,H)                                                  00008210
C     STORE DERIVATIVES AS NEEDED FOR ADAMS-MOULTON                     00008220
      GO TO (10,11,9),KFLAG                                             00008230
C     STORE Y PRIME AT N-2                                              00008240
   10 DO 18 I=1,N                                                       00008250
      N7I=I+N7                                                          00008260
      NNI=I+NN                                                          00008270
   18 T(N7I)=T(NNI)                                                     00008280
      GO TO 9                                                           00008290
C     STORE Y PRIME AT N-1                                              00008300
   11 DO 19 I=1,N                                                       00008310
      N6I=I+N6                                                          00008320
      NNI=I+NN                                                          00008330
   19 T(N6I)=T(NNI)                                                     00008340
    9 CONTINUE                                                          00008350
      GO TO 17                                                          00008360
C     ADAMS-MOULTON INTEGRATION                                         00008370
C     STORE CURRENT Y AND Y PRIME IN TEMPORARY STORAGE                  00008380
   12 DO 13 I=1,N                                                       00008390
      N9I=I+N9                                                          00008400
      N2I=I+N2                                                          00008410
      NNI=I+NN                                                          00008420
      T(N9I)=T(I)                                                       00008430
   13 T(N2I)=T(NNI)                                                     00008440
C     CALCULATE PREDICTED VALUE OF Y                                    00008450
      DO 14 I=1,N                                                       00008460
      N2I=I+N2                                                          00008470
      N6I=I+N6                                                          00008480
      N7I=I+N7                                                          00008490
      N8I=I+N8                                                          00008500
      N9I=I+N9                                                          00008510
      YP=T(N9I)+H24*(CMPLX(55.0,0.0)*T(N2I)                             00008520
     *-CMPLX(59.0,0.0)*T(N6I)+CMPLX(37.0,0.0)*                          00008530
     *T(N7I)-CMPLX(9.0,0.0)*T(N8I))                                     00008540
C     STORE AS PREDICTED FUNCTIONAL VALUE OF Y AT X=X+H                 00008550
   14 T(I)=YP                                                           00008560
C     STEP UP X                                                         00008570
      X=X+H                                                             00008580
C     CALCULATE Y PRIME USING PREDICTED Y                               00008590
      CALL DAUX(T,X,H)                                                  00008600
C     CALCULATE CORRECTED Y                                             00008610
      DO 15 I=1,N                                                       00008620
      NNI=I+NN                                                          00008630
      N2I=I+N2                                                          00008640
      N6I=I+N6                                                          00008650
      N7I=I+N7                                                          00008660
      N9I=I+N9                                                          00008670
      YC=T(N9I)+H24*(CMPLX(9.0,0.0)*T(NNI)                              00008680
     *+CMPLX(19.0,0.0)*T(N2I)-CMPLX(5.0,0.0)*                           00008690
     *T(N6I)+T(N7I))                                                    00008700
C     STORE AS NEW CURRENT VALUE OF Y                                   00008710
   15 T(I)=YC                                                           00008720
C     CALCULATE Y PRIME TO BE USED IN NEW STEP                          00008730
      CALL DAUX(T,X,H)                                                  00008740
C     REARRANGE STORAGE OF PREVIOUS DERIVATIVES                         00008750
      DO 16 I=1,N                                                       00008760
      N8I=I+N8                                                          00008770
      N7I=I+N7                                                          00008780
      N6I=I+N6                                                          00008790
      N2I=I+N2                                                          00008800
C     Y PRIME (N-2) GOES TO (N-3)                                       00008810
      T(N8I)=T(N7I)                                                     00008820
C     Y PRIME (N-1) GOES TO (N-2)                                       00008830
      T(N7I)=T(N6I)                                                     00008840
C     Y PRIME (N) GOES TO (N-1)                                         00008850
   16 T(N6I)=T(N2I)                                                     00008860
   17 RETURN                                                            00008870
      END                                                               00008880
C                                                                       00008890
C  HERE ARE THE CALCULUS SUBROUTINES                                    00008900
C                                                                       00008910
      SUBROUTINE VEC(K,V,Z)                                             00008920
C     CALCULUS SUBROUTINE FOR VECTORIZING THE KTH INDEPENDENT VAR V.    00008930
      IMPLICIT COMPLEX (A-H,O-Z)                                        00008940
      DIMENSION Z(66)                                                   00008950
      COMMON NVAR,NDIM,NSEC,N,NP1,NDIMT,MAXINT,INT,NDE                  00008960
      DO 1 I=1,NDIM                                                     00008970
    1 Z(I)=CMPLX(0.0,0.0)                                               00008980
      Z(1)=V                                                            00008990
      Z(K+1)=CMPLX(1.0,0.0)                                             00009000
      RETURN                                                            00009010
      END                                                               00009020
C                                                                       00009030
      SUBROUTINE CON(C,Z)                                               00009040
C     CALCULUS SUBROUTINE FOR VECTORIZING A CONSTANT C.                 00009050
      IMPLICIT COMPLEX (A-H,O-Z)                                        00009060
      DIMENSION Z(66)                                                   00009070
      COMMON NVAR,NDIM,NSEC,N,NP1,NDIMT,MAXINT,INT,NDE                  00009080
      DO 1 I=1,NDIM                                                     00009090
    1 Z(I)=CMPLX(0.0,0.0)                                               00009100
      Z(1)=C                                                            00009110
      RETURN                                                            00009120
      END                                                               00009130
C                                                                       00009140
      SUBROUTINE RECP(X,Z)                                              00009150
C     CALCULUS SUBROUTINE FOR THE RECIPROCAL FUNCTION 1/X.              00009160
      IMPLICIT COMPLEX (A-H,O-Z)                                        00009170
      DIMENSION X(66),Z(66),F(5)                                        00009180
      F(1)=CMPLX(1.0,0.0)/X(1)                                          00009190
      F(2)=CMPLX(-1.0,0.0)/(X(1)*X(1))                                  00009200
      F(3)=CMPLX(2.0,0.0)/(X(1)*X(1)*X(1))                              00009210
      F(4)=CMPLX(-6.0,0.0)/(X(1)*X(1)*X(1)*X(1))                        00009220
      CALL DER(F,X,Z)                                                   00009230
      RETURN                                                            00009240
      END                                                               00009250
C                                                                       00009260
      SUBROUTINE COSS(X,Z)                                              00009270
C     CALCULUS SUBROUTINE FOR THE FUNCTION COS(X).                      00009280
      IMPLICIT COMPLEX (A-H,O-Z)                                        00009290
      DIMENSION X(66),Z(66),F(5)                                        00009300
      C=X(1)                                                            00009310
      F(1)=CCOS(C)                                                      00009320
      F(2)=-CSIN(C)                                                     00009330
      F(3)=-CCOS(C)                                                     00009340
      F(4)=CSIN(C)                                                      00009350
      CALL DER(F,X,Z)                                                   00009360
      RETURN                                                            00009370
      END                                                               00009380
C                                                                       00009390
      SUBROUTINE SINN(X,Y)                                              00009400
C     CALCULUS SUBROUTINE FOR THE FUNCTION SIN(X).                      00009410
      IMPLICIT COMPLEX (A-H,O-Z)                                        00009420
      DIMENSION X(66),Y(66),F(5)                                        00009430
      C=X(1)                                                            00009440
      F(1)=CSIN(C)                                                      00009450
      F(2)=CCOS(C)                                                      00009460
      F(3)=-CSIN(C)                                                     00009470
      F(4)=-CCOS(C)                                                     00009480
      CALL DER(F,X,Y)                                                   00009490
      RETURN                                                            00009500
      END                                                               00009510
C                                                                       00009520
      SUBROUTINE POW(X,C,Z)                                             00009530
C     CALCULUS SUBROUTINE FOR X RAISED TO A CONSTANT POWER C.           00009540
      IMPLICIT COMPLEX (A-H,O-Z)                                        00009550
      DIMENSION X(66),Z(66),F(5)                                        00009560
      C1=X(1)                                                           00009570
      C1=CLOG(C1)                                                       00009580
      C1=C1*C                                                           00009590
      C2=CEXP(C1)                                                       00009600
      F(1)=C2                                                           00009610
      F(2)=C*C2/X(1)                                                    00009620
      F(3)=(C-CMPLX(1.0,0.0))*F(2)/X(1)                                 00009630
      F(4)=(C-CMPLX(2.0,0.0))*F(3)/X(1)                                 00009640
      CALL DER(F,X,Z)                                                   00009650
      RETURN                                                            00009660
      END                                                               00009670
C                                                                       00009680
      SUBROUTINE EXPP(X,Z)                                              00009690
C     CALCULUS SUBROUTINE FOR THE EXPONENTIAL FUNCTION EXP(X).          00009700
      IMPLICIT COMPLEX (A-H,O-Z)                                        00009710
      DIMENSION X(66),Z(66),F(5)                                        00009720
      C=X(1)                                                            00009730
      F(1)=CEXP(C)                                                      00009740
      F(2)=F(1)                                                         00009750
      F(3)=F(1)                                                         00009760
      F(4)=F(1)                                                         00009770
      CALL DER(F,X,Z)                                                   00009780
      RETURN                                                            00009790
      END                                                               00009800
C                                                                       00009810
      SUBROUTINE XPY(X,Y,Z)                                             00009820
C     CALCULUS SUBROUTINE FOR X RAISED TO THE POWER Y.                  00009830
      IMPLICIT COMPLEX (A-H,O-Z)                                        00009840
      DIMENSION X(66),Y(66),Z(66),Z1(66),Z2(66)                         00009850
      CALL LOGG(X,Z1)                                                   00009860
      CALL MULT(Z1,Y,Z2)                                                00009870
      CALL EXPP(Z2,Z)                                                   00009880
      RETURN                                                            00009890
      END                                                               00009900
C                                                                       00009910
      SUBROUTINE CPX(C,X,Z)                                             00009920
C     CALCULUS SUBROUTINE FOR A CONSTANT C RAISED TO THE POWER X.       00009930
      IMPLICIT COMPLEX (A-H,O-Z)                                        00009940
      DIMENSION X(66),Z(66),F(5)                                        00009950
      C1=X(1)                                                           00009960
      C2=CLOG(C)                                                        00009970
      C3=C1*C2                                                          00009980
      F(1)=CEXP(C3)                                                     00009990
      F(2)=F(1)*C2                                                      00010000
      F(3)=F(2)*C2                                                      00010010
      F(4)=F(3)*C2                                                      00010020
      CALL DER(F,X,Z)                                                   00010030
      RETURN                                                            00010040
      END                                                               00010050
C                                                                       00010060
      SUBROUTINE LOGG(X,Z)                                              00010070
C     CALCULUS SUBROUTINE FOR THE FUNCTION LOG(X).                      00010080
      IMPLICIT COMPLEX (A-H,O-Z)                                        00010090
      DIMENSION X(66),Z(66),F(5)                                        00010100
      C=X(1)                                                            00010110
      F(1)=CLOG(C)                                                      00010120
      F(2)=CMPLX(1.0,0.0)/X(1)                                          00010130
      F(3)=-F(2)*F(2)                                                   00010140
      F(4)=CMPLX(2.0,0.0)/(C*C*C)                                       00010150
      CALL DER(F,X,Z)                                                   00010160
      RETURN                                                            00010170
      END                                                               00010180
C                                                                       00010190
      SUBROUTINE ADD(X,Y,Z)                                             00010200
C     CALCULUS SUBROUTINE FOR THE FUNCTION Z = X + Y.                   00010210
      IMPLICIT COMPLEX (A-H,O-Z)                                        00010220
      DIMENSION X(66),Y(66),Z(66)                                       00010230
      COMMON NVAR,NDIM,NSEC,N,NP1,NDIMT,MAXINT,INT,NDE                  00010240
      DO 1 I=1,NDIM                                                     00010250
      Z(I)=X(I)+Y(I)                                                    00010260
    1 CONTINUE                                                          00010270
      RETURN                                                            00010280
      END                                                               00010290
C                                                                       00010300
      SUBROUTINE SUB(X,Y,Z)                                             00010310
C     CALCULUS SUBROUTINE FOR THE FUNCTION Z = X - Y.                   00010320
      IMPLICIT COMPLEX (A-H,O-Z)                                        00010330
      DIMENSION X(66),Y(66),Z(66)                                       00010340
      COMMON NVAR,NDIM,NSEC,N,NP1,NDIMT,MAXINT,INT,NDE                  00010350
      DO 1 I=1,NDIM                                                     00010360
      Z(I)=X(I)-Y(I)                                                    00010370
    1 CONTINUE                                                          00010380
      RETURN                                                            00010390
      END                                                               00010400
C                                                                       00010410
      SUBROUTINE MULT(X,Y,Z)                                            00010420
C     CALCULUS SUBROUTINE FOR THE FUNCTION Z = X*Y.                     00010430
      IMPLICIT COMPLEX (A-H,O-Z)                                        00010440
      DIMENSION X(66),Y(66),Z(66)                                       00010450
      COMMON NVAR,NDIM,NSEC,N,NP1,NDIMT,MAXINT,INT,NDE                  00010460
      Z(1)=X(1)*Y(1)                                                    00010470
      DO 1 I=1,NVAR                                                     00010480
    1 Z(I+1)=X(1)*Y(I+1)+Y(1)*X(I+1)                                    00010490
      KK=NVAR+1                                                         00010500
      DO 2 I=1,NVAR                                                     00010510
      DO 3 J=I,NVAR                                                     00010520
      KK=KK+1                                                           00010530
    3 Z(KK)=X(KK)*Y(1)+X(I+1)*Y(J+1)+X(J+1)*Y(I+1)+X(1)*Y(KK)           00010540
    2 CONTINUE                                                          00010550
      RETURN                                                            00010560
      END                                                               00010570
C                                                                       00010580
      SUBROUTINE DIV(X,Y,Z)                                             00010590
C     CALCULUS SUBROUTINE FOR THE FUNCTION Z = X/Y.                     00010600
      IMPLICIT COMPLEX (A-H,O-Z)                                        00010610
      DIMENSION X(66),Y(66),Z(66),Z1(66)                                00010620
      CALL RECP(Y,Z1)                                                   00010630
      CALL MULT(X,Z1,Z)                                                 00010640
      RETURN                                                            00010650
      END                                                               00010660
C                                                                       00010670
      SUBROUTINE CONVE(C,X,Z)                                           00010680
C     CALCULUS SUBROUTINE FOR A CONSTANT C TIMES A VECTOR X.            00010690
      IMPLICIT COMPLEX (A-H,O-Z)                                        00010700
      DIMENSION X(66),Z(66)                                             00010710
      COMMON NVAR,NDIM,NSEC,N,NP1,NDIMT,MAXINT,INT,NDE                  00010720
      DO 1 I=1,NDIM                                                     00010730
      Z(I)=C*X(I)                                                       00010740
    1 CONTINUE                                                          00010750
      RETURN                                                            00010760
      END                                                               00010770
C                                                                       00010780
      SUBROUTINE DER(F,X,Z)                                             00010790
C     CALCULUS SUBROUTINE FOR THE FUNCTION OF A FUNCTION RULE:          00010800
C     GIVEN U = F(X), THIS SUBROUTINE OBTAINS DU = F'(X)*DX AND         00010810
C     DDU = F''(X)*DX*DX  +  F'(X)*DDX                                  00010820
      IMPLICIT COMPLEX (A-H,O-Z)                                        00010830
      DIMENSION F(5),X(66),Z(66)                                        00010840
      COMMON NVAR,NDIM,NSEC,N,NP1,NDIMT,MAXINT,INT,NDE                  00010850
      Z(1)=F(1)                                                         00010860
      DO 1 I=1,NVAR                                                     00010870
    1 Z(I+1)=F(2)*X(I+1)                                                00010880
      KK=NVAR+1                                                         00010890
      DO 2 I=1,NVAR                                                     00010900
      DO 3 J=I,NVAR                                                     00010910
      KK=KK+1                                                           00010920
    3 Z(KK)=F(3)*X(I+1)*X(J+1)+F(2)*X(KK)                               00010930
    2 CONTINUE                                                          00010940
      RETURN                                                            00010950
      END                                                               00010960
C                                                                       00010970
C  HERE ARE THE MATRIX SUBROUTINES                                      00010980
C                                                                       00010990
      SUBROUTINE VECON(C,V,N,W)                                         00011000
C     CALCULATES THE PRODUCT C*V OF A SCALAR C AND AN N-DIMENSIONAL     00011010
C     VECTOR V AND PUTS IT IN W.                                        00011020
      IMPLICIT COMPLEX (A-H,O-Z)                                        00011030
      DIMENSION V(10),W(10)                                             00011040
      DO 10 I=1,N                                                       00011050
      W(I)=C*V(I)                                                       00011060
   10 CONTINUE                                                          00011070
      RETURN                                                            00011080
      END                                                               00011090
C                                                                       00011100
      SUBROUTINE CONMAT(C,B,N,Z)                                        00011110
C     CALCULATES THE PRODUCT C*B OF A CONSTANT C TIMES AN N BY N MATRIX 00011120
C     B AND PUTS IT IS Z.                                               00011130
      IMPLICIT COMPLEX (A-H,O-Z)                                        00011140
      DIMENSION B(10,10),Z(10,10)                                       00011150
      DO 10 I=1,N                                                       00011160
      DO 20 J=1,N                                                       00011170
      Z(I,J)=B(I,J)*C                                                   00011180
   20 CONTINUE                                                          00011190
   10 CONTINUE                                                          00011200
      RETURN                                                            00011210
      END                                                               00011220
C                                                                       00011230
      SUBROUTINE MULMAT(A,B,N,D)                                        00011240
C     MULTIPLIES TWO N BY N MATRICES A AND B AND PUTS THE PRODUCT       00011250
C     IN D.                                                             00011260
      IMPLICIT COMPLEX (A-H,O-Z)                                        00011270
      DIMENSION A(10,10),B(10,10),D(10,10)                              00011280
      DO 10 I=1,N                                                       00011290
      DO 20 J=1,N                                                       00011300
      SUM=CMPLX(0.0,0.0)                                                00011310
      DO 30 K=1,N                                                       00011320
      SUM=SUM+A(I,K)*B(K,J)                                             00011330
   30 CONTINUE                                                          00011340
      D(I,J)=SUM                                                        00011350
   20 CONTINUE                                                          00011360
   10 CONTINUE                                                          00011370
      RETURN                                                            00011380
      END                                                               00011390
C                                                                       00011400
      SUBROUTINE SUBMAT(A,B,N,D)                                        00011410
C     SUBTRACTS TWO N BY N MATRICES A AND B AND PUTS THEIR DIFFERENCE   00011420
C     IN D.                                                             00011430
      IMPLICIT COMPLEX (A-H,O-Z)                                        00011440
      DIMENSION B(10,10),D(10,10),A(10,10)                              00011450
      DO 10 I=1,N                                                       00011460
      DO 20 J=1,N                                                       00011470
      D(I,J)=A(I,J)-B(I,J)                                              00011480
   20 CONTINUE                                                          00011490
   10 CONTINUE                                                          00011500
      RETURN                                                            00011510
      END                                                               00011520
C                                                                       00011530
      SUBROUTINE MATVEC(A,W,N,V)                                        00011540
C     CALCULATES THE PRODUCT A*W OF AN N BY N MATRIX A AND AN           00011550
C     N-DIMENSIONAL VECTOR W AND PUTS THE PRODUCT IN V.                 00011560
      IMPLICIT COMPLEX (A-H,O-Z)                                        00011570
      DIMENSION W(10),V(10),A(10,10)                                    00011580
      DO 10 I=1,N                                                       00011590
      SUM=CMPLX(0.0,0.0)                                                00011600
      DO 20 J=1,N                                                       00011610
      SUM=SUM+A(I,J)*W(J)                                               00011620
   20 CONTINUE                                                          00011630
      V(I)=SUM                                                          00011640
   10 CONTINUE                                                          00011650
      RETURN                                                            00011660
      END                                                               00011670
C                                                                       00011680
      SUBROUTINE TRACE(A,N,TR)                                          00011690
C     CALCULATES THE TRACE OF AN N BY N MATRIX A AND PUTS THE RESULT    00011700
C     IN TR                                                             00011710
      IMPLICIT COMPLEX (A-H,O-Z)                                        00011720
      DIMENSION A(10,10)                                                00011730
      SUM=CMPLX(0.0,0.0)                                                00011740
      DO 10 I=1,N                                                       00011750
      SUM=SUM+A(I,I)                                                    00011760
   10 CONTINUE                                                          00011770
      TR=SUM                                                            00011780
      RETURN                                                            00011790
      END                                                               00011800

Copyright © Leigh Tesfatsion. All Rights Reserved.