C.hr snpcode3.f
C....*...1.........2.........3.........4.........5.........6.........7.*
C    
C     snpcode3.f
C
C     Copyright (C) 1990.
C
C     A. Ronald Gallant & George E. Tauchen 
C     P.O. Box 5513, Raleigh NC 27650-5513, USA   
C
C     Permission to use, copy, modify, and distribute this
C     software and its documentation for any purpose and
C     without fee is hereby granted, provided that the
C     above copyright notice appear in all copies and that
C     both that copyright notice and this permission notice
C     appear in supporting documentation.
C
C     This software is provided "as is" without any
C     expressed or implied warranty.
C
C....*...1.........2.........3.........4.........5.........6.........7.*
C.hr MULTIY
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C
C     MULTIY 1/26/88
C
C     PURPOSE: Fetch Ith multi-index for the Y polynomial
C     ARGUMENTS:
C          I - Integer pointer, between 1 and KY,
C              supplied by calling routine
C        IDX - Integer array of multi-indexes of length M,
C              returned by MULTIY.
C
C     WARNING: IF ZZMID IS CHANGED HERE IT MUST BE CHANGED IN
C               SUBROUTINES MULTIX, SETMYX, LNHYX AS WELL

C     WARNING: THE CALL TO SETMYX MUST PRECEDE USAGE OF THIS ROUTINE
C
      SUBROUTINE MULTIY(I,IDX)
      IMPLICIT REAL*8 (A-H,O-Z)
C
C----Declarations Shared between SETMYX and MULTIY,X -------------------
C
      INTEGER ZZMID
      PARAMETER ( ZZMID = 5000 )
      INTEGER ZM, ZML, ZKY, ZKX, ZMKY
      INTEGER MIDYX(ZZMID)
      COMMON /ZSVMYX/ MIDYX
      COMMON /ZSVPYX/ ZM, ZML, ZKY, ZKX, ZMKY
C
C     NOTE:
C        MIDYX(1..M*KY)            = Y MULTINDEXES STORED AS MxKY ARRAY
C        MIDYX(M*KY+1..M*KY+ML*KX) = X MULTINDEXES STORED AS MLxKX ARRAY
C        ZM = M, ZML = M*L, ZKY = KY, ZKX = KX, ZMKY = M*KY
C-----------------------------------------------------------------------
      DIMENSION IDX(1)
      SAVE
      NOUT=3
C
C ERROR TRAP
C
      IF ( (I .LT. 1) .OR. (I .GT. ZKY) ) THEN
         WRITE(NOUT,900) I
 900     FORMAT(//' *** MULTIY ERROR: INDEX OUT OF RANGE. ',
     1          //'     INDEX =  ', I15,
     2          //'     ABNORMAL RETURN'//)
         RETURN
      ENDIF
C
C FETCH Y MULTI-INDEX
C
      IB = (I-1)*ZM
      DO 10 K=1,ZM
         IDX(K) = MIDYX(IB+K)
 10   CONTINUE
      RETURN
      END
C.hr MULTIX
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C
C     MULTIX 1/26/88
C
C     PURPOSE: Fetch Ith multi-index for the X polynomial
C     ARGUMENTS:
C          I - Integer pointer, between 1 and KX,
C              supplied by calling routine
C        IDX - Integer array of multi-indexes of length M,
C              returned by MULTIY.
C
C     WARNING: IF ZZMID IS CHANGED HERE IT MUST BE CHANGED IN
C               SUBROUTINES MULTIY, SETMYX, LNHYX AS WELL
C
C     WARNING: THE CALL TO SETMYX MUST PRECEDE USAGE OF THIS ROUTINE
C
      SUBROUTINE MULTIX(I,IDX)
      IMPLICIT REAL*8 (A-H,O-Z)
C
C----Declarations Shared between SETMYX and MULTIY,X -------------------
C
      INTEGER ZZMID
      PARAMETER ( ZZMID = 5000 )
      INTEGER ZM, ZML, ZKY, ZKX, ZMKY
      INTEGER MIDYX(ZZMID)
      COMMON /ZSVMYX/ MIDYX
      COMMON /ZSVPYX/ ZM, ZML, ZKY, ZKX, ZMKY
C
C     NOTE:
C        MIDYX(1..M*KY)         =    Y MULTINDEXES STORED AS MxKY ARRAY
C        MIDYX(M*KY+1..M*KY+ML*KX) = X MULTINDEXES STORED AS MLxKX ARRAY
C        ZM = M, ZML = M*L, ZKY = KY, ZKX = KX, ZMKY = M*KY
C-----------------------------------------------------------------------
      DIMENSION IDX(1)
      SAVE
      NOUT=3
C
C ERROR TRAP
C
      IF ( (I .LT. 1) .OR. (I .GT. ZKX) ) THEN
         WRITE(NOUT,910) I
 910     FORMAT(//' *** MULTIX ERROR: INDEX OUT OF RANGE. ',
     1          //'     INDEX =  ', I15,
     2          //'     ABNORMAL RETURN'//)
         RETURN
      ENDIF
C
C FETCH X MULTI-INDEX
C
      IB = ZMKY + (I-1)*ZML
      DO 10 K=1,ZML
         IDX(K) = MIDYX(IB+K)
 10   CONTINUE
      RETURN
      END
C.hr SETMYX
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C
C     SETMYX  1/11/89
C
C     PURPOSE: Calculate KY, KX, and initialize the common blocks ZSVMYX
C              and ZSVPYX.  the block ZSVMYX holds the INTEGER array
C              MIDYX containing the multi-indexes for an SNP model.
C              Upon return, elements 1..M*KY of the array MIDYX are the
C              multi-indexes for the Y polynomial arranged as an MxKY
C              array; elements M*KY+1..M*KY+(M*L)*KX of MIDYX are
C              the multi-indexes for the X polynomial arranged as an
C              (M*L)xKX array.  the block ZSVPYX contains copies of
C              various constants as follows:
C                 ZM = M, ZML = M*L, ZKY = KY, ZKX = KX, ZMKY = M*KY
C              Thus, ZMKY+1 points to the start of the X multi-indexes.
C
C     USAGE:
C
C     CALL SETMYX(M,L,LPOLY,IDEGY,IDEGX,LINTY,LINTX,
C    &            KY,KX,IDUMMY,IWORK)
C
C     ARGUMENTS:
C           M - INTEGER, Dimension of the time series, supplied by the
C                 calling routine.
C           L - INTEGER, Lag length of X vector of the SNP model,
C                 supplied by calling routine.
C       LPOLY - INTEGER, Lag length of the X part of the polynomial,
C                 supplied by calling routine.  LPOLY cannot exceed L.
C       IDEGY - INTEGER, Degree of the Y polynomial, supplied by
C                 calling routine.
C       IDEGX - INTEGER, Degree of the X polynomial, supplied by
C                 calling routine.
C       LINTY - INTEGER, Limit on interactions included in the Y
C                 polynomial (See SUBROUTINE FLAGY), supplied by
C                 calling routine.
C       LINTX - INTEGER, Limit on interactions included in the X
C                 polynomial (See SUBROUTINE FLAGX), supplied by
C                 calling routine.
C          KY - INTEGER, Number of monomials in the Y polynomial
C                 characterized by M, IDEGY, and LINTY, returned.
C          KX - INTEGER, Number of monomials in the X polynomial
C                 characterized by M*L, IDEGX, and LINTX, returned.
C      IDUMMY - INTEGER, reserved for future use.
C       IWORK - An INTEGER work array of length 2*M*L.
C
C     WARNING:
C        IF ZZMID IS CHANGED HERE IT MUST BE CHANGED IN
C        SUBROUTINES MULTIY, MULTIX, LNHYX, EHYX AS WELL.
C
C     RESTRICTIONS:
C        0 <= LPOLY <= L
C        0 <= LINTY <= IDEGY
C        0 <= LINTX <= IDEGX
C
C
C     REMARKS: The routine works in tandem with TLOOP to generate all of
C        the multi-indexes.  For instance, generation of the Y multi-
C        indexes procedes as follows.  The indexes are generated by
C        degree, with the variable ITOT, ITOT=0,1,...,IDEGY, indicating
C        the degree of indexes currently be generated.  For each value
C        of ITOT, the routine TLOOP is repeatedly called.  The second
C        argument of TLOOP is an array of length M containing a multi-
C        index of degree ITOT; the first argument is an array of working
C        totals used by TLOOP.  (On the first call to TLOOP, its first
C        and second arguments must contains zeros in the first M-1
C        places  and ITOT in the Mth place.)  For each multi-index
C        returned by TLOOP, the routine FLAGY is called to see if that
C        multi-index meets the selection criterion. If it does, then
C        that multi-index is placed on a stack in the array MIDYX in the
C        common block ZZMID.  This continues until TLOOP sets IEND=1, in
C        which case all multi-indexes of degree ITOT have generated.
C           The X multi-indexes are generated in the same manner and
C        follow the Y multi-indexes on the stack.
C           If the routine determines that there is insufficient space
C        on the stack, then it stops placing them on the stack, but
C        continues generating and counting them to determine how much
C        larger the stack needs to be.
C
C     SUBROUTINES CALLED: TLOOP, FLAGY, FLAGX
C
      SUBROUTINE SETMYX(M,L,LPOLY,IDEGY,IDEGX,LINTY,LINTX,
     &                  KY,KX,IDUMMY,IWORK)
      IMPLICIT REAL*8 (A-H,O-Z)
C
C----Declarations Shared between SETMYX and MULTIY,X -------------------
C
      INTEGER ZZMID
      PARAMETER ( ZZMID = 5000 )
      INTEGER ZM, ZML, ZKY, ZKX, ZMKY
      INTEGER MIDYX(ZZMID)
      COMMON /ZSVMYX/ MIDYX
      COMMON /ZSVPYX/ ZM, ZML, ZKY, ZKX, ZMKY
C
C     NOTE:
C        MIDYX(1..M*KY)            = Y MULTINDEXES STORED AS MxKY ARRAY
C        MIDYX(M*KY+1..M*KY+ML*KX) = X MULTINDEXES STORED AS MLxKX ARRAY
C        ZM = M, ZML = M*L, ZKY = KY, ZKX = KX, ZMKY = M*KY
C-----------------------------------------------------------------------
      INTEGER IWORK(1)
      SAVE
      NOUT=3
C
C ERROR TRAP
C
      IF ( (M     .LT.   1) .OR. (L     .LT. 1    ) .OR.
     &     (LPOLY .LT.   0) .OR. (LPOLY .GT. L    ) .OR.
     &     (IDEGY .LT.   0) .OR. (IDEGX .LT. 0    ) .OR.
     &     (LINTY .LT.   0) .OR. (LINTY .GT. IDEGY) .OR.
     &     (LINTX .LT.   0) .OR. (LINTX .GT. IDEGX)   ) THEN
            WRITE(NOUT,900) M,L,LPOLY,IDEGY,IDEGX,LINTY,LINTX
 900     FORMAT(//' *** SETMYX ERROR: INVALID INPUT',/
     &            '      M = ', I5,/
     &            '      L = ', I5,/
     &            '  LPOLY = ', I5,/
     &            '  IDEGY = ', I5,/
     &            '  IDEGX = ', I5,//
     &            '  LINTY = ', I5,/
     &            '  LINTX = ', I5,//)
         STOP
      ENDIF
C
C START ROUTINE
C
      ML    = M*L
      ZM    = M
      ZML   = ML
      IOVER = 0
C
C  LOAD Y MULTI-INDEXES INTO ARRAY MIDYX
C
      IKS = 1
      IKV = IKS + M
      ICOUNT = 0
      DO 10 J=1,IDEGY+1
         ITOT = J - 1
         IEND = 0
         DO 15 K = 1, M
            IWORK(IKS-1+K) = 0
            IWORK(IKV-1+K) = 0
 15       CONTINUE
         IWORK(IKS-1+M) = ITOT
         IWORK(IKV-1+M) = ITOT
C
 20      IF ( IEND .EQ. 0 ) THEN
            LOC = ICOUNT*M
            CALL FLAGY(IWORK(IKV),M,IDEGY,LINTY,IFLAG,IDUMMY)
            IF (IFLAG .EQ. 0) THEN
               ICOUNT = ICOUNT + 1
               IF (LOC+M .LE. ZZMID) THEN
                  DO 30 K=1,M
                     MIDYX(LOC+K) = IWORK(IKV-1+K)
 30               CONTINUE
               ELSE
                  IOVER = 1
               ENDIF
            ENDIF
            CALL TLOOP(IWORK(IKS),IWORK(IKV),M,ITOT,IEND)
            GOTO 20
         ENDIF
 10   CONTINUE
      KY   = ICOUNT
      ZKY  = KY
      ZMKY = M*KY
C
C  LOAD X MULTI-INDEXES INTO ARRAY MIDYX
C
      IKS = 1
      IKV = IKS + ML
      ICOUNT = 0
      DO 60 J=1,IDEGX+1
         ITOT = J - 1
         IEND = 0
         DO 65 K = 1, ML
            IWORK(IKS-1+K) = 0
            IWORK(IKV-1+K) = 0
 65       CONTINUE
         IWORK(IKS-1+ML) = ITOT
         IWORK(IKV-1+ML) = ITOT
C
 70      IF ( IEND .EQ. 0 ) THEN
            LOC = ICOUNT*ML
            CALL FLAGX(IWORK(IKV),M,L,LPOLY,IDEGX,LINTX,IFLAG,IDUMMY)
            IF (IFLAG .EQ. 0) THEN
               ICOUNT = ICOUNT + 1
               IF (ZMKY+LOC+ML .LE. ZZMID) THEN
                  DO 80 K=1,ML
                     MIDYX(ZMKY+LOC+K) = IWORK(IKV-1+K)
 80               CONTINUE
               ELSE
                  IOVER = 1
               ENDIF
            ENDIF
            CALL TLOOP(IWORK(IKS),IWORK(IKV),ML,ITOT,IEND)
            GOTO 70
         ENDIF
 60   CONTINUE
      KX   = ICOUNT
      ZKX  = KX
C
C  CHECK IF SUFFICIENT SPACE HAD BEEN ALLOCATED TO STORE MULTI
C    INDEXES AND THEN EXIT
C
      IF (IOVER .EQ. 0) RETURN
C
C  IF WE GET HERE THERE WASN'T ENOUGH SPACE
C
      IREQ = M*KY + ML*KX
      WRITE(NOUT,910) ZZMID, IREQ, KY, KX
 910  FORMAT(' ***SETMYX ERROR: INSUFFICIENT SPACE FOR MULTINDEXES'/
     &       ' ALLOCATED (ZZMID) .........', I10/
     &       ' REQUIRED  (SET ZZMID TO)...', I10///
     &       ' KX ........................', I10/
     &       ' KY ........................', I10)
         STOP
      END
C.hr TLOOP
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
      SUBROUTINE TLOOP(KS,KV,M,ITOT,IEND)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION KS(M),KV(M)
      SAVE
C--------------------------------------C
C ERROR TRAPS                          C
C--------------------------------------C
      IER=0
      IEND=0
      NOUT=3
      IF (M .LT. 1) IER=1
      IF (M .GE. 1) THEN
         DO 5 L=1,M
            IF (KS(L) .LT. 0 .OR. KS(L) .GT. ITOT) IER=2
 5          IF (KV(L) .LT. 0 .OR. KV(L) .GT. ITOT) IER=3
      END IF
      IF ( IER .GT. 0 ) THEN
         WRITE(NOUT,900) IER, M
 900     FORMAT(//' *** TLOOP ERROR: BAD INPUT DATA ',
     1          //'     IER = ', I5,
     2          //'     NUMBER OF VARIABLES(M) = ', I15)
         IF ( M .GE. 1 ) THEN
            WRITE(NOUT,910)
 910        FORMAT(//'     ,   KS,   KV'//)
            DO 7 L=1,M
               WRITE(NOUT,920) L, KS(L), KV(L)
 920           FORMAT(3I5)
 7          CONTINUE
         ENDIF
      ENDIF
C--------------------------------------C
C START ROUTINE                        C
C--------------------------------------C
      IF (M .GE. 2) THEN
         KV(M-1) = KV(M-1) + 1
         IF (M .GE. 3) THEN
            DO 30 L=1,M-2
               IF (KV(M-L) .LT. ITOT-KS(M-L-1)+1) GO TO 40
               KV(M-L)=0
               KV(M-L-1) = KV(M-L-1) + 1
 30         CONTINUE
 40         CONTINUE
            KS(1)=KV(1)
            DO 50 L=2,M-1
 50         KS(L)=KS(L-1)+KV(L)
            KV(M)=ITOT-KS(M-1)
         ELSE
            KS(1)=KV(1)
            KV(2)=ITOT-KV(1)
         END IF
      END IF
      IF (M .EQ. 1 .OR. KV(1) .EQ. ITOT + 1) IEND=1
      RETURN
      END
C.HR  FLAGY
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C
C     1/16/89
C     FLAGY
C
C     PURPOSE:
C        Inspect a Y multi-index and set a flag indicating whether it
C        is to be screened out.
C
C     USAGE:
C     CALL FLAGY(ALPHA,M,IDEGY,LINTY,IFLAG,IDUMMY)
C
C     ARGUMENTS:
C
C     ALPHA  - INTEGER array of length M containing multi-index to be
C              inspected.
C     M      - INTEGER, length of Y, supplied by calling routine.
C     IDEGY  - INTEGER, degree of Y polynomial, supplied by
C                calling routine.
C     LINTY  - INTEGER, limit on iteractions, supplied by
C                calling routine.
C     IFLAG  - INTEGER, flag (0 means retain ALPHA, 1 means
C                             reject ALPHA), returned.
C     IDUMMY - INTEGER, reserved for future use.
C
C     REMARKS:
C        IFLAG will be set to 1, indicating rejection, if the degree
C        of the interaction exceeds IDEGY-LINTY; otherwise IFLAG is set
C        to zero.
C
      SUBROUTINE FLAGY(ALPHA,M,IDEGY,LINTY,IFLAG,IDUMMY)
      IMPLICIT REAL*8 (A-H,O-Z)
      INTEGER ALPHA(1)
      SAVE
      IFLAG = 0
      IDG = 0
      NGZ = 0
      DO 10 J=1,M
         IF (ALPHA(J) .GT. 0) THEN
            NGZ = NGZ + 1
            IDG = IDG + ALPHA(J)
         ENDIF
 10   CONTINUE
      IF (NGZ .GE. 2 .AND. IDG .GT. IDEGY-LINTY) IFLAG=1
      CONTINUE
      RETURN
      END
C.HR  FLAGX
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C
C
C     1/16/89
C     FLAGX
C
C     PURPOSE:
C        Inspect a X multi-index and set a flag indicating whether it
C        is to be screened out.
C
C     USAGE:
C     CALL FLAGX(BETA,M,L,LPOLY,IDEGX,LINTX,IFLAG,IDUMMY)
C
C     ARGUMENTS:
C
C     BETA  -  INTEGER array of length M*L containing multi-index to be
C                inspected.
C     M      - INTEGER, length of the time series, supplied by
C                calling routine.
C     L      - INTEGER, maximum lag length of the model (X is of
C                length M*L), supplied by calling routine.
C     LPOLY  - INTEGER, lag length of polynomial part of the model,
C                supplied by calling routine. LPOLY<=L.
C     IDEGX  - INTEGER, degree of X polynomial, supplied by
C                calling routine.
C     LINTX  - INTEGER, limit on iteractions, supplied by
C                calling routine.
C     IFLAG  - INTEGER, flag (0 means retain ALPHA, 1 means
C                             reject ALPHA), returned.
C     IDUMMY - INTEGER, reserved for future use.
C
C     REMARKS:
C        IFLAG will be set to 1, indicating rejection, if BETA
C        equals zero in the first M*(L-LPOLY) places and if the degree
C        of the interaction exceeds IDEGX-LINTX; otherwise IFLAG is set
C        to zero.
C
      SUBROUTINE FLAGX(BETA,M,L,LPOLY,IDEGX,LINTX,IFLAG,IDUMMY)
      IMPLICIT REAL*8 (A-H,O-Z)
      INTEGER BETA(1)
      SAVE
      IFLAG = 0
      DO 10 J=1,M*(L-LPOLY)
         IF (BETA(J) .GT. 0) THEN
            IFLAG = 1
            RETURN
         ENDIF
 10   CONTINUE
      IDG   = 0
      NGZ   = 0
      DO 20 J=M*(L-LPOLY)+1,M*L
         IF (BETA(J) .GT. 0) THEN
            NGZ = NGZ + 1
            IDG = IDG + BETA(J)
         ENDIF
 20   CONTINUE
      IF (NGZ .GE. 2 .AND. IDG .GT. IDEGX-LINTX) IFLAG=1
      RETURN
      END
C.hr SNPRST
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C     SUBROUTINE SNPRST 1/12/89
C
C     PURPOSE
C     MODIFY THETA FROM A SMALLER SNPRX MODEL TO A
C     LARGER SNPRX MODEL.
C
C     USAGE
C      CALL SNPRST(THETA,M,
C     &     LRX0,LPOLY0,IDEGY0,IDEGX0,LINTY0,LINTX0,
C     &     LRX1,LPOLY1,IDEGY1,IDEGX1,LINTY1,LINTX1,TWEAK,
C     &     IDUMMY,IWORK)
C
C     ARGUMENTS
C     THETA  - A REAL*8 ARRAY OF LENGTH LTHETA1.
C                THETA IS MODIFIED BY SNPRST:  ON ENTRY T0 SNPRST, THETA
C                IS ARRANGED FOR AN SNPRX(LPOLY0,LRX0,IDEGY0,IDEGX0)
C                MODEL; ON EXIT IT ARRANGED FOR AN
C                SNPRX(LPOLY1,LRX1,IDEGY1,IDEGX1) MODEL. THE POLYNOMIAL
C                PARAMETERS OF THE LARGER MODEL THAT ARE IMPLICITLY
C                ZERO IN THE LARGER MODEL ARE SET TO TWEAK, AS ARE
C                THE VAR AND RX PARAMETERS THAT ARE IMPLICITLY ZERO
C                IN THE LARGER MODEL.
C     M      - INTEGER, THE LENGTH OF Y, SUPPLIED BY CALLING ROUTINE.
C     LRX0   - INTEGER, LAG LENGTH OF RX PART OF THE SMALLER
C                MODEL, SUPPLIED BY CALLING ROUTINE.
C     LPOLY0 - INTEGER, LAG LENGTH OF THE X POLYNOMIAL IN THE
C                SMALLER MODEL, SUPPLIED BY CALLING ROUTINE.
C     IDEGY0 - INTEGER, DEGREE OF Y POLYNOMIAL OF THE SMALLER MODEL,
C                SUPPLIED BY CALLING ROUTINE.
C     IDEGX0 - INTEGER, DEGREE OF X POLYNOMIAL OF THE SMALLER MODEL,
C                SUPPLIED BY CALLING ROUTINE.
C     LINTY0 - INTEGER, HIGHEST DEGREE INTERACTION ALLOWED IN THE Y
C                 POLYNOMIAL OF THE SMALLER MODEL, SUPPLIED BY
C                 CALLING ROUTINE.
C     LINTX0 - INTEGER, HIGHEST DEGREE INTERACTION ALLOWED IN THE X
C                 POLYNOMIAL OF THE SMALLER MODEL, SUPPLIED BY
C                 CALLING ROUTINE.
C     LRX1   - INTEGER, LAG LENGTH OF RX PART OF THE LARGER MODEL
C                MODEL, SUPPLIED BY CALLING ROUTINE.
C     LPOLY1 - INTEGER, LAG LENGTH OF THE X POLYNOMIAL IN THE
C                LARGER MODEL, SUPPLIED BY CALLING ROUTINE.
C     IDEGY1 - INTEGER, DEGREE OF Y POLYNOMIAL OF THE LARGER MODEL,
C                SUPPLIED BY CALLING ROUTINE.
C     IDEGX1 - INTEGER, DEGREE OF X POLYNOMIAL OF THE LARGER MODEL,
C                SUPPLIED BY CALLING ROUTINE.
C     LINTY1 - INTEGER, HIGHEST DEGREE INTERACTION ALLOWED IN THE Y
C                 POLYNOMIAL OF THE LARGER MODEL, SUPPLIED BY
C                 CALLING ROUTINE.
C     LINTX1 - INTEGER, HIGHEST DEGREE INTERACTION ALLOWED IN THE X
C                 POLYNOMIAL OF THE LARGER MODEL, SUPPLIED BY
C                 CALLING ROUTINE.
C     TWEAK  - TWEAK VALUE FOR THOSE ELEMENTS OF THETA THAT ARE
C               IMPLICITLY SET TO ZERO IN THE SMALLER MODEL
C     IDUMMY - INTEGER, RESERVED FOR FUTURE USE.
C
C     IWORK  - INTEGER WORK ARRAY OF LENGTH AT LEAST
C                2*M*L1 + M*KY0 + M*L0*KX0 + KY0 + KX0 where
C                L0 = MAX0(LRX0,LPOLY0)
C                L1 = MAX0(LRX1,LPOLY1)
C
C     WHERE
C        LTHETA0 = KX0*KY0 + M + L0*M*M + (LRX0+1)*M*(M+1)/2
C        LTHETA1 = KX1*KY1 + M + L1*M*M + (LRX1+1)*M*(M+1)/2
C
C     SUBROUTINES CALLED
C        SETMYX, MULTIY, MULTIX
C
C     REMARKS
C        THE WORK ON THETA IS DONE IN PLACE.  IF THE MODELS ARE
C        OF THE SAME SIZE THEN THE ROUTINE HAS NO EFFECT ON THETA.
C
C     WARNING
C        THIS ROUTINE USES SETMYX TO GENERATE MULTI-INDEXES
C        THE CALLS TO SETMYX MIGHT NEED TO BE MODIFIED
C        IF THAT ROUTINE IS EVER CHANGED.
C
      SUBROUTINE SNPRST(THETA,M,
     &    LRX0,LPOLY0,IDEGY0,IDEGX0,LINTY0,LINTX0,
     &    LRX1,LPOLY1,IDEGY1,IDEGX1,LINTY1,LINTX1,TWEAK,
     &    IDUMMY,IWORK)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 THETA(1)
      INTEGER IWORK(1)
      COMMON /ZPRINT/ IPRINT
      SAVE
C
C  SET SOME CONSTANTS
C
      NOUT = 3
C
      L0   = MAX0(LRX0,LPOLY0)
      L1   = MAX0(LRX1,LPOLY1)
      ML0  = M*L0
      ML1  = M*L1
C
C   VERIFY PARAMETERS
C
      IF (IPRINT.EQ.0) THEN
        WRITE(NOUT,3000) M,L0,L1,
     &    LRX0,   LRX1,
     &    LPOLY0, LPOLY1,
     &    IDEGY0, IDEGY1,
     &    IDEGX0, IDEGX1,
     &    LINTY0, LINTY1,
     &    LINTX0, LINTX1, TWEAK
      ENDIF
 3000   FORMAT(//5X,'REINITIALIZING THETA FOR A RESTART'/
     &           5X,'----------------------------------'//
     & ' M.........', I3/
     & ' L0........', I3, ' L1........ ', I3//
     & ' LRX0......', I3, ' LRX1...... ', I3/
     & ' LPOLY0....', I3, ' LPOLY1.... ', I3/
     & ' IDEGY0....', I3, ' IDEGY1.... ', I3/
     & ' IDEGX0....', I3, ' IDEGX1.... ', I3/
     & ' LINTY0....', I3, ' LINTY1.... ', I3/
     & ' LINTX0....', I3, ' LINTX1.... ', I3/
     & ' TWEAK.....', F10.5/)
C
C  CALL SETMYX TO SETUP MULTI-INDEXES FOR THE SMALLER MODEL.
C
      LSETYX = 2*ML1
      ISETYX = 1
      CALL SETMYX(M,L0,LPOLY0,IDEGY0,IDEGX0,LINTY0,LINTX0,KY0,KX0,
     &            IDUMMY,IWORK(ISETYX))
C
C   DEFINE POINTERS TO MULTI-INDEXES OF THE SMALLER MODEL AND
C     SAVE THESE MULTI-INDEXES IN IWORK.
C
      IMY0   = ISETYX + LSETYX
      IMX0   = IMY0   + M*KY0
      ITOP   = IMX0   + ML0*KX0
      DO 110 IY=1,KY0
        CALL MULTIY(IY,IWORK(ITOP))
         DO 112 J=1,M
            IWORK(IMY0-1+J+(IY-1)*M) = IWORK(ITOP-1+J)
 112     CONTINUE
 110  CONTINUE
      DO 120 IX=1,KX0
         CALL MULTIX(IX,IWORK(ITOP))
         DO 122 J=1,ML0
            IWORK(IMX0-1+J+(IX-1)*ML0) = IWORK(ITOP-1+J)
 122     CONTINUE
 120  CONTINUE
C
C  CALL SETMYX TO RESET THE MULTI-INDEXES FOR THE LARGER MODEL.
C
      CALL SETMYX(M,L1,LPOLY1,IDEGY1,IDEGX1,LINTY1,LINTX1,KY1,KX1,
     &            IDUMMY,IWORK(ISETYX))
C
C  COMPUTE POINTERS TO ARRAYS THAT HOLD THE REMAPPINGS
C    OF THE Y AND X MULTI-INDEXS.
C
      IYD  = ITOP
      IXD  = IYD  +  KY0
      ITOP = IXD  +  KX0
C
C  DETERMINE THE MAPPINGS BETWEEN THE TWO SETS OF Y AND X MULTI-INDEXES.
C
C  FOR EACH Y MULTI-INDEX IN THE SMALLER MODEL, SCAN FOR A MATCH
C    IN THE LARGER MODEL.  FALL THROUGH LOOP 144 MEANS A
C    MATCH IS FOUND AND THE LOCATION IS SAVED.  FALL THROUGH
C    LOOP 142 MEANS NO MATCH WAS FOUND, WHICH IS FATAL.
C
      IALPHA = ISETYX
      DO 140 IY0=1,KY0
         DO 142 IY1=1,KY1
            CALL MULTIY(IY1,IWORK(IALPHA))
            DO 144 J=1,M
               IF ( IWORK(IMY0  -1+J+(IY0-1)*M) .NE.
     &              IWORK(IALPHA-1+J          )      ) GOTO 142
 144        CONTINUE
            IWORK(IYD-1+IY0) = IY1
            GOTO 140
 142     CONTINUE
         GOTO 980
 140  CONTINUE
C
C  FOR EACH X MULTI-INDEX IN THE SMALLER MODEL SCAN FOR A MATCH
C    IN THE LARGER MODEL.  FALL THROUGH LOOP 154 MEANS A
C    MATCH IS FOUND AND THE LOCATION IS SAVED.  FALL THROUGH
C    LOOP 152 MEANS NO MATCH WAS FOUND, WHICH IS A FATAL ERROR.
C
      IBETA = ISETYX
      DO 150 IX0=1,KX0
         DO 152 IX1=1,KX1
            CALL MULTIX(IX1,IWORK(IBETA))
            DO 153 J=1,(L1-LPOLY0)*M
               IF (IWORK(IBETA-1+J) .NE. 0) GOTO 152
 153        CONTINUE
            DO 154 J=1,LPOLY0*M
               IF ( IWORK(IBETA+ML1-J              ) .NE.
     &              IWORK(IMX0 +ML0-J + (IX0-1)*ML0)    ) GOTO 152
 154        CONTINUE
            IWORK(IXD-1+IX0) = IX1
            GOTO 150
 152     CONTINUE
         GOTO 990
 150  CONTINUE
C
C   SET SOME CONSTANTS AND THEN REARRANGE THETA, STARTING
C   AND THE END AND MOVING TOWARDS THE BEGINNING.
C
      NA0   = KX0*KY0
      NA1   = KX1*KY1
      NQ0   = M + M*ML0
      NQ1   = M + M*ML1
      LR    = (M*(M+1))/2
      NR0   = LR*(1+M*LRX0)
      NR1   = LR*(1+M*LRX1)
      NT0   = NA0 + NQ0 + NR0
      NT1   = NA1 + NQ1 + NR1
C
C   MOVE  RX(LRX1...1) AND PUT 0.0 IN THE EMPTY SPACE
C
      DO 21 I=1,NR0-LR
         THETA(NT1-I+1) = THETA(NT0-I+1)
 21   CONTINUE
      DO 22 I=1,NR1-NR0
         THETA(NA1+NQ1+LR+I) = 0.0
 22   CONTINUE
C
C   MOVE RX(0)
C
      DO 23 I=1,LR
         THETA(NA1+NQ1+LR-I+1) = THETA(NA0+NQ0+LR-I+1)
23    CONTINUE
C
C   MOVE  B AND PUT TWEAK IN EXTRA PLACE
C
      DO 24 I=1,M*ML0
         THETA(NA1+NQ1-I+1) = THETA(NA0+NQ0-I+1)
 24   CONTINUE
      DO 26 I=1,M*M*(L1-L0)
         THETA(NA1+M+I) = TWEAK
 26   CONTINUE
C
C   MOVE B0
C
      DO 28 I=1,M
         THETA(NA1+M-I+1) = THETA(NA0+M-I+1)
 28   CONTINUE
C
C  NOW PUT IN TWEAK AND MOVE THE POLYNOMIAL PARAMETERS
C
      DO 200 IXY=KX0*KY0+1,KX1*KY1
 200     THETA(IXY) = TWEAK
      DO 210 IY0=1,KY0
         IYY0 = KY0 - IY0  + 1
         IYY1 = IWORK(IYD-1+IYY0)
         DO 216 IX0=1,KX0
            IXX0 = KX0 - IX0  + 1
            IXX1 = IWORK(IXD-1+IXX0)
            IXY0 = IXX0 + (IYY0-1)*KX0
            IXY1 = IXX1 + (IYY1-1)*KX1
            THETA(IXY1) = THETA(IXY0)
            IF (IXY1 .GT. IXY0) THETA(IXY0) = TWEAK
 216     CONTINUE
 210  CONTINUE
C
C  DONE
C
      RETURN
C
C  HERE WE HANDLE THE SITUATIONS WHERE NO MATCHES WERE FOUND.
C  THE PROGRAM IS HALTED BECAUSE THE RESULTS WILL BE MEANINGLESS.
C
 980  CONTINUE
      WRITE(NOUT,3981) M,
     & LRX0,   LRX1,
     & LPOLY0, LPOLY1,
     & IDEGY0, IDEGY1,
     & IDEGX0, IDEGX1,
     & LINTY0, LINTY1,
     & LINTX0, LINTX1
      STOP
 990  CONTINUE
      WRITE(NOUT,3982) M,
     & LRX0,   LRX1,
     & LPOLY0, LPOLY1,
     & IDEGY0, IDEGY1,
     & IDEGX0, IDEGX1,
     & LINTY0, LINTY1,
     & LINTX0, LINTX1
       STOP
 3981 FORMAT(
     &/' ***ERROR: SNPRST: UNABLE TO FIND MATCHING Y MULTI-INDEX'/
     & ' M.........', I3/
     & ' LRX0......', I3, ' LRX1...... ', I3/
     & ' LPOLY0....', I3, ' LPOLY1.... ', I3/
     & ' IDEGY0....', I3, ' IDEGY1.... ', I3/
     & ' IDEGX0....', I3, ' IDEGX1.... ', I3/
     & ' LINTY0....', I3, ' LINTY1.... ', I3/
     & ' LINTX0....', I3, ' LINTX1.... ', I3/)
 3982 FORMAT(
     &/' ***ERROR: SNPRST: UNABLE TO FIND MATCHING X MULTI-INDEX'/
     & ' M.........', I3/
     & ' LRX0......', I3, ' LRX1...... ', I3/
     & ' LPOLY0....', I3, ' LPOLY1.... ', I3/
     & ' IDEGY0....', I3, ' IDEGY1.... ', I3/
     & ' IDEGX0....', I3, ' IDEGX1.... ', I3/
     & ' LINTY0....', I3, ' LINTY1.... ', I3/
     & ' LINTX0....', I3, ' LINTX1.... ', I3/)
      END
