 
C     ALGORITHM 634 COLLECTED ALGORITHMS FROM ACM.
C     ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.11, NO. 3,
C     SEPT., 1985, P. 201-217; 218-228.
C     PROGRAM GENERA
C
C ====================================================================
C DATA GENERATOR -- DATA GENERATOR -- DATA GENERATOR -- DATA GENERATOR
C ====================================================================
C
      INTEGER FITDEG,DIMEN,NFPOLS,NFPTS,NEPTS,NEPOLS,EVLDEG,TOPS
C
C     ***************
C     A  PRIMITIVE  DATA-GENERATING   PROGRAM   FOR   USE   WITH   THE
C     JEZIORANSKI-BARTELS    SUITE   OF   ROUTINES   FOR   MULTINOMIAL
C     LEAST-SQUARES FITTING.
C
C     THE OUTPUT FORMATS IN (GENDAT) AND IN  (GENEVL)  ARE  CONSISTENT
C     WITH  THE INPUT FORMATS TO BE FOUND IN THE SIMPLE DRIVER PROGRAM
C     WHICH IS INCLUDED WITH THE SUITE.
C
C     THE DATA  STATEMENTS BELOW  WILL CREATE  A PROBLEM REQUIRING THE
C     FIT OF A  3-VARIABLE MULTINOMIAL, CONSISTING OF  8 TERMS, TO 27
C     POINTS OF  GENERATED DATA  AND THEN  THE EVALUATION  OF THE FULL
C     RESULTING MULTINOMIAL AT 5 POINTS.
C
C     DATE LAST MODIFIED
C     ---- ---- --------
C     DECEMBER 14, 1984
C     ****************
C
      DATA FITDEG/-1/
      DATA EVLDEG/-1/
      DATA NFPOLS/8/
      DATA DIMEN/3/
      DATA NFPTS/27/
      DATA NEPTS/5/
      DATA NEPOLS/8/
      DATA TOPS/3/
C
      CALL GENDAT(DIMEN,NFPTS,FITDEG,NFPOLS,TOPS)
C
      CALL GENEVL(DIMEN,NEPTS,EVLDEG,NEPOLS,TOPS)
C
      STOP
      END
      SUBROUTINE GENDAT(DIMEN,NFPTS,FITDEG,NFPOLS,TOPS)
C
      INTEGER DIMEN,NFPTS,FITDEG,NFPOLS,I,J,K,COUNT,TOPS
      DOUBLE PRECISION X,Y,Z,W,F,XSTART,YSTART,ZSTART
      DOUBLE PRECISION XDEL,YDEL,ZDEL,WEIGHT
C
C     ***************
C     A SUBROUTINE TO GENERATE AND  PRINT  DATA  POINTS  FOR  THE
C     SPECIFIC FUNCTION
C
C                          SIN(X)
C               F(X,Y,Z) = ------ * COS(Y) + EXP(Z)
C                            X
C
C     ***************
C
      DATA XSTART /0.1D+00/
      DATA YSTART /0.1D+00/
      DATA ZSTART /0.1D+00/
      DATA XDEL   /0.2D+00/
      DATA YDEL   /0.2D+00/
      DATA ZDEL   /0.2D+00/
      DATA WEIGHT /1.0D+00/
C
      WRITE (6,6000) DIMEN,FITDEG,NFPOLS,NFPTS
      W = WEIGHT
      COUNT = 0
      X = XSTART
      DO 30 I = 1,TOPS
         Y = YSTART
         DO 20 J = 1,TOPS
            Z = ZSTART
            DO 10 K = 1,TOPS
               IF (COUNT .GE. NFPTS)  GO TO 40
               F = (DSIN(X)/X)*DCOS(Y) + DEXP(Z)
               WRITE (6,6010) W,X,Y,Z,F
               Z = Z + ZDEL
               COUNT = COUNT + 1
   10       CONTINUE
            Y = Y + YDEL
   20    CONTINUE
         X = X + XDEL
   30 CONTINUE
   40 CONTINUE
      RETURN
C
 6000 FORMAT(4I5)
 6010 FORMAT(5D14.6)
      END
      SUBROUTINE GENEVL(DIMEN,NEPTS,EVLDEG,NEPOLS,TOPS)
C
      INTEGER DIMEN,NEPTS,I,NEPOLS,EVLDEG,COUNT,TOPS
      DOUBLE PRECISION X,Y,Z,XSTART,YSTART,ZSTART
      DOUBLE PRECISION XDEL,YDEL,ZDEL
C
C     ***************
C     A SUBROUTINE TO GENERATE EVALUATION POINTS.
C     ***************
C
      DATA XSTART /0.15D+00/
      DATA YSTART /0.16D+00/
      DATA ZSTART /0.17D+00/
      DATA XDEL   /0.05D+00/
      DATA YDEL   /0.06D+00/
      DATA ZDEL   /0.07D+00/
C
      COUNT = 0
      WRITE (6,6000) EVLDEG,NEPOLS,NEPTS
      X = XSTART
      DO 30 I = 1,TOPS
         Y = YSTART
         DO 20 J = 1,TOPS
            Z = ZSTART
            DO 10 K = 1,TOPS
               IF (COUNT .GE. NEPTS)  GO TO 40
               WRITE (6,6010) X,Y,Z
               Z = Z + ZDEL
               COUNT = COUNT + 1
   10       CONTINUE
            Y = Y + YDEL
   20    CONTINUE
         X = X + XDEL
   30 CONTINUE
   40 CONTINUE
      RETURN
C
 6000 FORMAT(3I5)
 6010 FORMAT(3D14.6)
C
      END
C     PROGRAM DRIVER
C
C ====================================================================
C   SAMPLE DRIVER -- SAMPLE DRIVER -- SAMPLE DRIVER -- SAMPLE DRIVER
C ====================================================================
C
      INTEGER DIMEN,FITDEG,NFPOLS,NFPTS,EVLDEG,NEPOLS,NEPTS
      INTEGER ERROR,FIWKLN,FDWKLN,EDWKLN,IREQD,DREQD
      INTEGER FITIWK(89)
      DOUBLE PRECISION FITDWK(2201),FITVLS(125),FITCDS(375),WTS(125)
      DOUBLE PRECISION EVLDWK(125), EVLVLS(20), EVLCDS(60), RESIDS(125)
C
C     ***************
C     A  SIMPLE  DRIVER  PROGRAM  TO  READ  TEST  DATA  AND  CALL  THE
C     JEZIORANSKI-BARTELS  SUITE  OF MULTINOMIAL LEAST-SQUARES FITTING
C     ROUTINES.
C
C     THE DATA AND ARRAY DECLARATIONS  ARE  CONSISTENT  WITH  PROBLEMS
C     INVOLVING  UP  TO  3  VARIABLES,  125   FITTING  POINTS,  AND 20
C     EVALUATION POINTS USING 20 BASIS ORTHOGONAL MULTINOMIALS IN BOTH
C     THE  FIT AND THE EVALUATION.
C
C     DATE LAST MODIFIED
C     ---- ---- --------
C     DECEMBER 15, 1984
C     ****************
C
      DATA FDWKLN/2201/
      DATA EDWKLN/125/
      DATA FIWKLN/89/
C
C     ***************
C     READ FITTING DATA
C     ***************
C
      READ (5,5000) DIMEN,FITDEG,NFPOLS,NFPTS
      WRITE (6,6000) DIMEN,FITDEG,NFPOLS,NFPTS
      CALL INFIT(DIMEN,NFPTS,FITCDS,FITVLS,WTS)
C
C     ***************
C     COMPUTE THE LEAST-SQUARES FIT
C     ***************
C
      CALL CONSTR(DIMEN,FITDEG,NFPOLS,NFPTS,FITCDS,NFPTS,
     +            FITVLS,WTS,RESIDS,.TRUE.,ERROR,FITIWK,
     +            FITDWK,FIWKLN,FDWKLN,IREQD,DREQD)
C
C     ***************
C     PRINT RESIDUALS AT THE FITTING POINTS
C
C     THE USER COULD CHECK  IREQD,DREQD  AT THIS POINT TO SEE
C     THE NUMBER OF LOCATIONS WHICH WERE ACTUALLY REQUIRED IN
C     THE ARRAYS  FITIWK,FITDWK.
C     ***************
C
      WRITE (6,6010)
      CALL OUTRES(NFPTS,RESIDS)
C
C     ***************
C     READ POINTS OF EVALUATION
C     ***************
C
      READ (5,5000) EVLDEG,NEPOLS,NEPTS
      WRITE (6,6020) EVLDEG,NEPOLS,NEPTS
      CALL INEVL(DIMEN,NEPTS,EVLCDS)
C
C     ***************
C     EVALUATE THE FITTING MULTINOMIAL
C     ***************
C
      CALL EVAL(DIMEN,EVLDEG,NEPOLS,NEPTS,EVLCDS,EVLVLS,
     +          ERROR,FITIWK,FITDWK,FIWKLN,FDWKLN,EVLDWK,EDWKLN)
C
C     ***************
C     PRINT OUT THE ARRAY OF MULTINOMIAL VALUES
C     ***************
C
      CALL OUTEVL(DIMEN,NEPTS,NEPOLS,EVLCDS,EVLVLS)
C
      STOP
C
C     ***************
C     FORMATS
C     ***************
C
 5000 FORMAT(5I5)
 6000 FORMAT(//29H MULTINOMIAL FITTING PROBLEM.//
     +       10H DIMEN  = ,I5/
     +       10H FITDEG = ,I5/
     +       10H NFPOLS = ,I5/
     +       10H NFPTS  = ,I5)
 6010 FORMAT(//18H FITTING COMPLETE./13H RESIDUALS...
     +       //4X,1HI,5X,8HRESIDUAL)
 6020 FORMAT(//24H MULTINOMIAL EVALUATION.//
     +       10H EVLDEG = ,I5/
     +       10H NEPOLS = ,I5/
     +       10H NEPTS  = ,I5)
C
      END
      SUBROUTINE INEVL(DIMEN,NEPTS,EVLCDS)
C
      INTEGER DIMEN,NEPTS,I,J
      DOUBLE PRECISION EVLCDS(NEPTS,DIMEN)
C
C     ***************
C     SUBROUTINE TO READ THE ARRAY OF EVALUATION POINTS.
C     ***************
C
      DO 100 I=1,NEPTS
        READ (5,5000) (EVLCDS(I,J),J=1,DIMEN)
  100 CONTINUE
      RETURN
 5000 FORMAT(4D14.6)
      END
      SUBROUTINE INFIT(DIMEN,NFPTS,FITCDS,FITVLS,WTS)
C
      INTEGER NFPTS,DIMEN,I,J
      DOUBLE PRECISION FITCDS(NFPTS,DIMEN),FITVLS(NFPTS),WTS(NFPTS)
C
C     ***************
C     SUBROUTINE TO READ THE FITTING DATA.
C     ***************
C
      WRITE (6,6000)
      DO 100 I=1,NFPTS
        READ (5,5000) WTS(I),(FITCDS(I,J),J=1,DIMEN),FITVLS(I)
        WRITE (6,6010) I,WTS(I),(FITCDS(I,J),J=1,DIMEN),FITVLS(I)
  100 CONTINUE
      RETURN
 5000 FORMAT(5D14.6)
 6000 FORMAT(/8H DATA.../
     +       4X,1HI,5X,6HWEIGHT,19X,11HCOORDINATES,17X,11HDATA VALUES)
 6010 FORMAT(I5,5D14.6)
      END
      SUBROUTINE OUTEVL(DIMEN,NEPTS,NEPOLS,EVLCDS,EVLVLS)
C
      INTEGER DIMEN,NEPTS,NEPOLS,I,J
      DOUBLE PRECISION EVLCDS(NEPTS,DIMEN),EVLVLS(NEPTS)
C
C     ***************
C     SUBROUTINE TO PRINT OUT THE RESULTS OF THE EVALUATION.
C     ***************
C
      WRITE (6,6000)
      DO 100 I=1,NEPTS
        WRITE (6,6010) I,(EVLCDS(I,J),J=1,DIMEN),EVLVLS(I)
  100 CONTINUE
      RETURN
 6000 FORMAT(/8H DATA.../
     +       4X,1HI,22X,11HCOORDINATES,14X,5HVALUE)
 6010 FORMAT(I5,4D14.6)
      END
      SUBROUTINE OUTRES(NFPTS,RESIDS)
C
      INTEGER NFPTS,I
      DOUBLE PRECISION RESIDS(NFPTS)
C
C     ***************
C     SUBROUTINE TO PRINT OUT THE RESIDUALS FROM THE FIT.
C     ***************
C
      DO 100 I=1,NFPTS
        WRITE (6,6000) I,RESIDS(I)
  100 CONTINUE
      RETURN
C
 6000 FORMAT(I5,D14.6)
      END
C
C ====================================================================
C   CONSTRUCT FIT -- CONSTRUCT FIT -- CONSTRUCT FIT -- CONSTRUCT FIT
C ====================================================================
C
      SUBROUTINE CONSTR(DIMEN,FITDEG,NFPOLS,NFPTS,
     +                  FITCDS,NCROWS,FITVLS,WTS,
     +                  RESIDS,NEWFIT,ERROR,FITIWK,
     +                  FITDWK,FIWKLN,FDWKLN,IREQD,DREQD)
C
      INTEGER NFPOLS,ONPLYS,FITDEG,NFPTS,DIMEN,OLDEG,FIWKLN,FDWKLN
      INTEGER ERROR,IREQD,DREQD,OPSWID,OLALFL,INDSTT,P,DIMP1,NCROWS
      INTEGER NEWSTT,MAXSTT,ALFSTT,PSISTT,CSTT,SSQSTT,PSIWID,ALFL
      INTEGER FITIWK(FIWKLN)
      DOUBLE PRECISION FITDWK(FDWKLN),FITCDS(NCROWS,DIMEN)
      DOUBLE PRECISION FITVLS(NFPTS),RESIDS(NFPTS)
      DOUBLE PRECISION WTS(NFPTS)
      DOUBLE PRECISION SCALE
      LOGICAL NEWFIT
C
C     ***************
C     PURPOSE
C     -------
C
C     THIS SUBROUTINE CONSTRUCTS A LEAST-SQUARES  MULTINOMIAL  FIT  TO
C     GIVEN DATA USING A BASIS OF ORTHOGONAL MULTINOMIALS.
C
C     THE DATA FOR THE FIT IS GIVEN IN THE ARRAYS  FITCDS ,   FITVLS ,
C     AND   WTS .    FITCDS  IS A DOUBLE-PRECISION MATRIX, EACH ROW OF
C     WHICH CONTAINS  AN  OBSERVATION  POINT  (ORDERED  COLLECTION  OF
C     VARIABLE  VALUES).    FITVLS   IS  A  DOUBLE-PRECISION,  SINGLY-
C     INDEXED ARRAY,  EACH  ELEMENT  OF  WHICH  CONTAINS  AN  OBSERVED
C     FUNCTION VALUE CORRESPONDING TO AN  OBSERVATION POINT.   WTS  IS
C     A DOUBLE-PRECISION, SINGLY-INDEXED ARRAY, EACH ELEMENT OF  WHICH
C     IS A NONNEGATIVE WEIGHT FOR THE CORRESPONDING OBSERVATION.
C
C     THE FIT WHICH IS PRODUCED IS A MULTINOMIAL EXPRESSED IN THE FORM
C
C      C  PSI  (X ,...,X     ) +...+ C       PSI       (X ,...,X     )
C       1    1   1      DIMEN         NFPOLS    NFPOLS   1      DIMEN
C
C     WHERE THE VALUE OF  NFPOLS  WILL BE AS GIVEN (IF  FITDEG .LT. 0)
C     OR  AS  COMPUTED  BY  CONSTR  TO GIVE A FULL-DEGREE FIT (IN CASE
C      FITDEG  IS SPECIFIED .GE. 0).  THE ELEMENTS
C
C         PSI  (X ,...,X     )
C            K   1      DIMEN
C
C     FORM A BASIS FOR  THE  MULTINOMIALS  WHICH  IS  ORTHOGONAL  WITH
C     RESPECT TO THE WEIGHTS AND OBSERVATION POINTS.
C
C     THE EXTENT OF THE FIT CAN BE SPECIFIED IN ONE OF TWO WAYS.
C         IF THE PARAMETER FITDEG IS SET .GE. 0, THEN A COMPLETE BASIS
C         FOR THE MULTINOMIALS OF DEGREE =  FITDEG  WILL BE USED.  (AN
C         ERROR WILL BE  FLAGGED  IF  THIS  WILL  REQUIRE  MORE  BASIS
C         MULTINOMIALS  THAN  THE  NUMBER  OF  DATA  POINTS WHICH WERE
C         GIVEN.)
C         IF THE PARAMETER  FITDEG  IS .LT.  0, THEN  NFPOLS  WILL  BE
C         TAKEN AS THE COUNT OF THE NUMBER OF BASIS MULTINOMIALS TO BE
C         USED FOR A PARTIAL-DEGREE FIT.  (AN ERROR WILL BE FLAGGED IF
C          NFPOLS  .LT. 0.)
C
C     NOTE, THE CALL TO  CONSTR  WITH  NEWFIT  = .TRUE.  CAN  BE  MADE
C         WITH  THE  PARAMETERS  SET  FOR  THE  MAXIMUM  FIT  DESIRED.
C         SEVERAL SUBSEQUENT CALLS TO  CONSTR  WITH  NEWFIT  = .FALSE.
C         CAN  BE  MADE WITH SMALLER VALUES OF  FITDEG  OR  NFPOLS  AS
C         MAY BE DESIRED TO OBTAIN A PARTIAL FIT.
C
C     VARIABLES
C     ---------
C
C      DIMEN  -- (INTEGER) -- (PASSED)
C         THE NUMBER OF VARIABLES.
C      FITDEG  - (INTEGER) -- (PASSED/RETURNED)
C         IGNORED IF .LT. 0.
C         IF DEGREE .GE. 0 THEN  DEGREE  IS  CHECKED  AGAINST  NFPTS .
C         THE VALUE OF  DEGREE  WILL BE REDUCED IF THERE IS A BASIS OF
C         MULTINOMIALS, ALL OF  DEGREE  .LE. DEGREE ,  OF  CARDINALITY
C          NFPTS .  SEE  ERROR  BELOW.
C      NFPOLS  - (INTEGER) -- (PASSED/RETURNED)
C         IGNORED IF  DEGREE  .GE. 0.
C         IF  DEGREE .LT. 0 THEN THE VALUE OF NFPOLS WILL BE TAKEN  AS
C         THE SIZE OF THE BASIS OF MULTINOMIALS TO BE USED IN THE FIT.
C          NFPOLS  MUST SATISFY  NFPOLS  .LT. NFPTS  AND NFPOLS .GE. 1
C         SEE  ERROR  BELOW.
C      NFPTS  --- (INTEGER) -- (PASSED)
C         THE NUMBER OF DATA POINTS TO BE USED IN THE FIT.
C          NFPTS  MUST BE .GE. 1.  SEE  ERROR  BELOW.
C      FITCDS  -- (DOUBLE-PRECISION, 2-SUBSCRIPT ARRAY) -- (PASSED)
C          FITCDS (P,K) IS THE VALUE OF THE K-TH VARIABLE AT THE  P-TH
C         DATA POINT.
C      NCROWS  -- (INTEGER) -- (PASSED)
C         THE ROW DIMENSION  DECLARED  FOR   FITCDS   IN  THE  CALLING
C         PROGRAM.
C      FITVLS  -- (DOUBLE-PRECISION, 1-SUBSCRIPT ARRAY) -- (PASSED)
C          FITVLS (P) IS THE OBSERVED FUNCTION VALUE OF THE P-TH  DATA
C         POINT.
C      WTS  ----- (DOUBLE-PRECISION, 1-SUBSCRIPT ARRAY) -- (PASSED)
C          WTS (P) IS THE WEIGHT ATTACHED TO THE P-TH DATA POINT.
C      RESIDS  -- (DOUBLE-PRECISION, 1-SUBSCRIPT ARRAY) -- (RETURNED)
C          RESIDS (P) IS THE DIFFERENCE BETWEEN THE FITTED FUNCTION AT
C         POINT P AND  FITVLS (P).
C      NEWFIT  -- (LOGICAL) -- (PASSED)
C         A  LOGICAL FLAG.  IF  NEWFIT =.TRUE., THEN THIS IS THE FIRST
C         FIT TO BE CARRIED OUT WITH THE DATA TO BE FOUND IN THE OTHER
C         PARAMETERS TO   CONSTR ,  AND  SPACE  FOR A  FIT  IS  TO  BE
C         ALLOCATED. IF  NEWFIT =.FALSE., THEN A FIT OF ANOTHER DEGREE
C         CAN BE CONSTRUCTED IN THE SPACE ALLOCATED ON A PREVIOUS CALL
C         WITH THE SAME DATA, AND CERTAIN INITIALIZATION STEPS ARE BY-
C         PASSED.
C      ERROR  -- (INTEGER) -- (RETURNED)
C       0 IF  NFPOLS , DIMEN ,   DEGREE ,   NFPTS   AND   IWKLEN   ARE
C         VALID AND CONSISTENT WITH EACH OTHER.
C       1 IF  DEGREE  .GE. 0 BUT THERE IS AN INTERPOLATING MULTINOMIAL
C         OF SMALLER DEGREE OR IF DEGREE .LT. 0 AND NFPOLS .GT. NFPTS.
C       2 IF  DEGREE  .LT. 0 AND  NFPOLS  .LE. 0.
C       3 IF  NFPTS  .LT. 1 AND/OR  DIMEN  .LT. 1.
C       4 IF  IWKLEN  AND/OR  DWKLEN  IS TOO SMALL.  (SET  IWKLEN   TO
C         THE VALUE RETURNED IN  IREQD , AND SET  DWKLEN  TO THE VALUE
C         RETURNED IN  DREQD  TO RESOLVE THIS PROBLEM.)
C       5 NEWFIT  = .FALSE. BUT  ONPLYS  .GE.  NFPOLS.
C       6 ERROR RETURN FROM  INCDG .  THERE IS  NO  MORE ROOM  IN  THE
C         FITDWK  AND/OR  FITIWK  ARRAYS  TO INCLUDE MORE TERMS IN THE
C         FIT.
C      FITIWK  -- (INTEGER, 1-SUBSCRIPT ARRAY) -- (RETURNED)
C         AN INTEGER WORK ARRAY OF LENGTH  FIWKLN .  UPON RETURN  FROM
C         A CALL TO  CONSTR  WITH  NEWFIT  SET .TRUE., SOME  DIMENSION
C         AND ARRAY-LENGTH INFORMATION WILL BE INSERTED.  UPON  RETURN
C         FROM A CALL TO  CONSTR  WITH  NEWFIT  SET .FALSE.,  DETAILED
C         INDEXING INFORMATION (LOCATION OF COEFFICIENTS IN   FITDWK ,
C         ETC.) IS INSERTED.
C      FITDWK  -- (DOUBLE-PRECISION, 1-SUBSCRIPT ARRAY) -- (RETURNED)
C         A DOUBLE PRECISION ARRAY  OF LENGTH   FDWKLN .  UPON  RETURN
C         FROM  CONSTR  WITH  NEWFIT  SET .FALSE.,  THE  FULL  DETAILS
C         OF THE REQUESTED FIT (COEFFICIENTS, ETC.) WILL BE INSERTED.
C      FIWKLN  -- (INTEGER) -- (PASSED)
C         THE LENGTH OF THE ARRAY  FITIWK .
C      FDWKLN  -- (INTEGER) -- (PASSED)
C         THE LENGTH OF THE ARRAY  FITDWK .
C      IREQD  -- (INTEGER) -- (PASSED)
C         THE LENGTH WHICH THE ARRAY  FITIWK  REALLY NEEDS TO BE.
C      DREQD  -- (INTEGER) -- (PASSED)
C         THE LENGTH WHICH THE ARRAY  FITDWK  REALLY NEEDS TO BE.
C
C
C     NOTE, THE 10 AND 70  LOOPS  (I.E.  THE  LOOPS  FOR  SCALING  THE
C     RESIDUALS)  DEPEND  ON  THE  SCALING  SCHEME USED.  THE RESIDUAL
C     SCALING  MUST  BE  CONSISTENT  WITH  THAT  DEFINED BY   SCALPM ,
C      SCALDN , AND  SCALUP .
C
C      CONSTR   CALLS  ALLOT ,  RESTRT ,  INCDG , AND  GNRTP
C
C
C     DATE LAST MODIFIED
C     ---- ---- --------
C     OCTOBER 16, 1984
C     ****************
C
      DIMP1 = DIMEN + 1
C
C     ***************
C     TEST IF THIS IS A NEW FIT OR AN ADDITION TO THE OLD ONE.
C     ***************
C
      IF ( NEWFIT ) GO TO 20
C
C     ***************
C     THE SECTION BELOW IS FOR AN ADDITION TO A PREVIOUS FIT.
C     OBTAIN THE VALUES OF SAVED PARAMETERS.
C     ***************
C
         OLDEG = FITIWK(1)
         ONPLYS = FITIWK(2)
         OPSWID = FITIWK(3)
         OLALFL = FITIWK(4)
         SCALE = FITDWK(DIMP1)
         SCALE = SCALE * SCALE
C
         IF ( NFPOLS .GT. ONPLYS ) GO TO 10
            ERROR = 5
            RETURN
   10    CONTINUE
   20 CONTINUE
C
C     ***************
C     NEW FITTING PROBLEMS BEGIN HERE.  OLD FITTING PROBLEMS
C     PROCEED INTO THIS SECTION FROM ABOVE.
C     ***************
C
      CALL ALLOT(FITDEG,NFPOLS,NFPTS,DIMEN,FITIWK,FIWKLN,IREQD,DREQD,
     +           ERROR)
      IF ( ERROR .GE. 2 ) RETURN
C
      IF ( FDWKLN .GE. DREQD ) GO TO 30
         ERROR = 4
         RETURN
   30 CONTINUE
C
      PSIWID = FITIWK(3)
      ALFL = FITIWK(4)
      INDSTT = 1
      NEWSTT = 4 * NFPOLS + INDSTT
      MAXSTT = 1
      ALFSTT = MAXSTT + DIMP1
      CSTT = ALFSTT + ALFL
      SSQSTT = CSTT + NFPOLS
      PSISTT = SSQSTT + NFPOLS
C
      IF ( NEWFIT ) GO TO 50
C
C     ***************
C     FOR A  COMPLETELY NEW FIT,  SKIP TO 50.  IF THIS IS
C     AN ADDITION TO A PREVIOUS FIT, A PREVIOUSLY EXISTING
C     SCALING MUST BE RESTORED TO THE RESIDUALS.
C     ***************
C
         DO 40 P = 1,NFPTS
            RESIDS(P) = RESIDS(P) / SCALE
   40    CONTINUE
C
C     ***************
C     RE-ARRANGE AND SHUFFLE INFORMATION IN THE WORKING ARRAYS.
C     ***************
C
         CALL RESTRT(PSIWID,FITDWK,DREQD,ONPLYS,OPSWID,OLALFL,CSTT,
     +               SSQSTT,PSISTT,NFPTS,DIMEN)
C
C     ***************
C     PRODUCE THE AUGMENTED FIT.
C     ***************
C
         CALL INCDG(FITDEG,FITDWK(ALFSTT),FITDWK(PSISTT),FITIWK(INDSTT),
     +              FITIWK(NEWSTT),FITDWK(SSQSTT),FITCDS,
     +              NFPOLS,DIMEN,NFPTS,FITVLS,RESIDS,
     +              FITDWK(CSTT),PSIWID,WTS,ALFL,ONPLYS,OLDEG,ERROR)
         GO TO 60
C
   50    CONTINUE
C
C     ***************
C     THE CODE BELOW IS FOR A COMPLETELY NEW FIT.  ARRANGE PARAMETERS,
C     INDICES, AND STORAGE ALLOCATION IN THE WORK ARRAYS.  SCALE THE
C     DATA AND RESIDUALS, AND THEN PRODUCE THE FIT.
C     ***************
C
         CALL GNRTP(FITDEG,FITDWK(ALFSTT),
     +              FITDWK(PSISTT),FITIWK(INDSTT),
     +              FITIWK(NEWSTT),FITDWK(SSQSTT),
     +              FITCDS,NFPOLS,DIMEN,NFPTS,FITVLS,RESIDS,
     +              FITDWK(CSTT),PSIWID,WTS,ALFL,DIMP1,FITDWK(MAXSTT))
         SCALE = FITDWK(DIMEN + 1)
         SCALE = SCALE * SCALE
C
   60 CONTINUE
C
C     ***************
C     UNSCALE THE RESIDUALS FOR THE BENEFIT OF THE USER.
C     ***************
C
      DO 70 P = 1,NFPTS
        RESIDS(P) = RESIDS(P) * SCALE
   70 CONTINUE
      RETURN
      END
C
C ====================================================================
C     EVALUATE FIT -- EVALUATE FIT -- EVALUATE FIT -- EVALUATE FIT
C ====================================================================
C
      SUBROUTINE EVAL(DIMEN,EVLDEG,NEPOLS,NEPTS,EVLCDS,EVLVLS,
     +                ERROR,FITIWK,FITDWK,FIWKLN,FDWKLN,EVLDWK,EDWKLN)
C
      INTEGER FIWKLN,FDWKLN,NEPOLS,NEPTS,DIMEN,ERROR,MAXSTT,ALFSTT,CSTT
      INTEGER GBASIZ,ALFL,DIMP1,EVLDEG,TOP,BOT,CURDEG,EDWKLN
      INTEGER FITIWK(FIWKLN)
      DOUBLE PRECISION FITDWK(FDWKLN),EVLDWK(EDWKLN),EVLCDS(NEPTS,DIMEN)
      DOUBLE PRECISION EVLVLS(NEPTS)
C
C     ***************
C     PURPOSE
C     -------
C
C     THIS SUBROUTINE  EVALUATES  THE  LEAST-SQUARES  MULTINOMIAL  FIT
C     WHICH HAS BEEN PREVIOUSLY PRODUCED BY  CONSTR .  EITHER THE FULL
C     MULTINOMIAL AS PRODUCED MAY BE EVALUATED,  OR  ONLY  AN  INITIAL
C     SEGMENT THEREOF.  AS IN THE CASE WITH  CONSTR , IT IS POSSIBLE
C     (1) TO SPECIFY MULTINOMIALS OF A FULL GIVEN DEGREE, OR
C     (2) TO SPECIFY  THE  NUMBER  OF  ORTHOGONAL  BASIS  ELEMENTS  TO
C         ACHIEVE A PARTIAL-DEGREE FIT.
C
C     IN CASE (1), THE  DESIRED  DEGREE  IS  GIVEN  AS  THE  VALUE  OF
C      EVLDEG   (WHICH  MUST  BE  NONNEGATIVE AND NOT GREATER THAN THE
C     VALUE USED FOR  FITDEG  IN  CONSTR ), AND THE PARAMETER  NEPOLS
C     WILL  BE  SET  BY  EVAL  TO SPECIFY THE NUMBER OF BASIS ELEMENTS
C     REQUIRED. IF  EVLDEG  .LT.  FITDEG   IS  GIVEN,  THEN  ONLY  THE
C     INITIAL  PORTION OF THE FITTING MULTINOMIAL (OF DEGREE  EVLDEG )
C     WILL BE EVALUATED.
C
C     IN CASE (2),  EVLDEG  IS TO BE SET NEGATIVE, IN WHICH  CASE  THE
C     VALUE  OF   NEPOLS  (WHICH MUST BE POSITIVE AND NOT GREATER THAN
C     THE VALUE USED FOR   NFPOLS   IN   CONSTR )  WILL  BE  TAKEN  AS
C     DEFINING  THE  INITIAL  PORTION OF THE FITTING MULTINOMIAL TO BE
C     EVALUATED.
C
C     IF  NEPOLS  =  NFPOLS  (WITH  EVLDEG  .LT.  0),  OR   EVLDEG   =
C     FITDEG  (WITH   EVLDEG   .GT.  0),  THEN  THE  FULL  MULTINOMIAL
C     GENERATED BY  CONSTR  WILL BE EVALUATED.
C
C     THE  EVALUATION  WILL  TAKE  PLACE  FOR  EACH  OF   THE   POINTS
C     (COLLECTION  OF  VARIABLE  VALUES)  GIVEN AS A ROW OF THE MATRIX
C      EVLCDS .  THE  VALUES  PRODUCED  FROM  THE  FULL,  OR  PARTIAL,
C     MULTINOMIAL WILL BE PLACED IN THE ARRAY  EVLVLS .
C
C     VARIABLES
C     ---------
C
C      DIMEN  -- (INTEGER) -- (PASSED)
C         THE NUMBER OF VARIABLES.
C      EVLDEG  -- (INTEGER) -- (PASSED)
C         IF  EVLDEG  .LT. 0, THEN THIS PARAMETER WILL BE IGNORED.
C         IF  EVLDEG  .GE. 0, THEN THE VALUE OF  EVLDEG  MUST  SATISFY
C          EVLDEG  .LE. (THE  DEGREE  OF THE APPROXIMATING MULTINOMIAL
C         GENERATED IN  CONSTR ).  IN THIS CASE  EVLDEG  WILL  SPECIFY
C         THE DEGREE OF THE INITIAL PORTION OF THE FITTING MULTINOMIAL
C         TO BE EVALUATED.
C      NEPOLS  -- (INTEGER) -- (PASSED/RETURNED)
C         IF  EVLDEG  .GE. 0, THEN THIS PARAMETER WILL BE IGNORED.
C         IF EVLDEG .LT. 0, THEN THE PARTIAL MULTINOMIAL INVOLVING THE
C         FIRST   NEPOLS  ORTHOGONAL BASIS FUNCTIONS WILL BE EVALUATED
C         AT THE POINTS GIVEN BY  EVLCDS .  THE RESULTING VALUES  WILL
C         BE STORED IN  EVLVLS .
C         THE VALUE OF NEPOLS MUST BE .GE. 1 AND .LE. (THE SIZE OF THE
C         BASIS  GENERATED  IN   CONSTR ),  WHICH  WAS RETURNED AS THE
C         VALUE OF  NFPOLS .
C         NEPOLS  WILL BE CHANGED IF EVLDEG .GT. 0 TO GIVE THE SIZE OF
C         BASIS REQUIRED FOR THE MULTINOMIAL OF DEGREE  EVLDEG .
C      NEPTS  -- (INTEGER) -- (PASSED)
C         THE NUMBER OF EVALUATION POINTS.
C      EVLCDS  -- (INTEGER ) -- (PASSED)
C          EVLCDS (P,K) IS THE VALUE OF THE K-TH VARIABLE AT THE  P-TH
C         EVALUATION POINT.
C      EVLVLS  -- (INTEGER) -- (RETURNED)
C          EVLVLS (P) IS THE VALUE OF THE EVALUATED MULTINOMIAL AT THE
C         P-TH EVALUATION POINT.
C      ERROR  -- (INTEGER) -- (RETURNED)
C          0 .........  IF NO ERRORS
C         -1 .........  IF  NEPOLS  .GT.  NFPOLS  OR  NEPOLS  .LT. 1
C         -2 .........  IF  NEPTS  .LT. 1 OR  DIMEN  .LT. 1
C          NEPOLS  ...  IF  NEPOLS  .GT.  EDWKLN
C      FITIWK  -- (INTEGER, 1-SUBSCRIPT ARRAY) -- (PASSED)
C         THE INTEGER WORK ARRAY OF LENGTH  FIWKLN  THAT WAS  USED  IN
C          CONSTR .
C      FITDWK  -- (DOUBLE-PRECISION, 2-SUBSCRIPT ARRAY) -- (PASSED)
C         THE DOUBLE PRECISION WORK ARRAY OF LENGTH  FDWKLN  THAT  WAS
C         USED IN  CONSTR .
C      FIWKLN  -- (INTEGER) -- (PASSED)
C         THE LENGTH OF  FITIWK .
C      FDWKLN  -- (INTEGER) -- (PASSED)
C         THE LENGTH OF  FITDWK .
C      EVLDWK  -- (DOUBLE-PRECISION, 1-SUBSCRIPT ARRAY ) -- (RETURNED)
C         A WORK ARRAY OF LENGTH  NEPOLS  (OR LONGER).
C      EDWKLN  -- (INTEGER) -- (PASSED)
C         THE LENGTH OF  EVLDWK .
C
C     THE SUBROUTINE   EVALP   IS CALLED TO DO THE ACTUAL EVALUATION.
C
C     DATE LAST MODIFIED
C     ---- ---- --------
C     OCTOBER 16, 1984
C     ****************
C
C
C     ***************
C     SET UP INDEX POINTERS TO THE BEGINNING OF EACH ROW OF
C     THE TABLE -- THIS SETS THE BEGINNING POINT FOR EACH
C     FULL MULTINOMIAL DEGREE.
C     ***************
C
      ERROR = 0
      GBASIZ = FITIWK(2)
      IF ( EVLDEG .LT. 0 ) GO TO 20
         TOP = 1
         BOT = 1
         DO 10 CURDEG = 1,EVLDEG
            TOP = TOP * (DIMEN + CURDEG)
 10         BOT = BOT * CURDEG
         NEPOLS = TOP / BOT
         IF ( EVLDEG .EQ. 0 ) NEPOLS = 1
 20   IF ( NEPOLS .LE. GBASIZ .AND. NEPOLS.GE.1 ) GO TO 30
         ERROR = -1
         RETURN
 30   IF ( NEPTS .GE. 1 .AND. DIMEN .GE. 1 ) GO TO 40
         ERROR = -2
         RETURN
 40   IF ( NEPOLS .LE. EDWKLN ) GO TO 50
         ERROR = NEPOLS
         RETURN
 50   DIMP1 = DIMEN + 1
      ALFL = FITIWK(4)
      MAXSTT = 1
      ALFSTT = DIMP1 + MAXSTT
      CSTT = ALFSTT + ALFL
C
C     ***************
C     THE ACTUAL EVALUATION IS DONE INSIDE  EVALP .
C     ***************
C
      CALL EVALP(EVLCDS,FITDWK(CSTT),NEPTS,DIMEN,NEPOLS,FITDWK(ALFSTT),
     +           FITIWK,EVLDWK,EVLVLS,ALFL,FITDWK(MAXSTT),DIMP1)
      RETURN
      END
C
C ====================================================================
C     UTILITIES -- UTILITIES -- UTILITIES -- UTILITIES -- UTILITIES
C ====================================================================
C
      SUBROUTINE ALLOT(DEGREE,NPOLYS,NPTS,DIMEN,IWORK,IWKLEN,
     +                 IREQD,DREQD,ERROR)
C
      INTEGER IREQD,DREQD,ALFL,ERROR,NPOLYS,DEGREE,DIMEN,NPTS
      INTEGER NEWSTT,PSIWID,KMXBAS,STARTJ,KJP1D2,INDEX,IWKLEN
      INTEGER NPLYT4
      INTEGER IWORK(IWKLEN)
C
C     ***************
C     PURPOSE
C     -------
C
C      ALLOT  CHECKS FOR SUFFICIENCY THE DECLARED  DIMENSIONS  OF  THE
C     WORK  ARRAYS  USED BY THE SUBROUTINE  CONSTR .  VARIOUS SIZES OF
C     SUB-ARRAYS ARE COMPUTED AND REPORTED.
C
C     THIS ROUTINE IS CALLED BY  CONSTR .  IT IS NOT  CALLED  DIRECTLY
C     BY THE USER.
C
C     THIS ROUTINE CALLS   BASIZ   AND   TABLE   FOR  THE  SUBSTANTIVE
C     COMPUTATIONS.
C
C     VARIABLES
C     ---------
C
C      DEGREE  - (PASSED/RETURNED)
C         IGNORED IF .LT. 0.
C         IF  DEGREE  .GE. 0 THEN  DEGREE IS CHECKED  AGAINST   NPTS .
C         THE VALUE OF  DEGREE  WILL BE REDUCED IF THERE IS A BASIS OF
C         MULTINOMIALS, ALL OF  DEGREE .LE.  DEGREE ,  OF  CARDINALITY
C          NPTS
C      NPOLYS  - (PASSED/RETURNED)
C         IGNORED IF  DEGREE  .GE. 0.
C         IF DEGREE .LT. 0 THEN THE VALUE OF NPOLYS WILL BE TAKEN  AS
C         THE SIZE OF THE BASIS OF MULTINOMIALS TO BE USED IN THE FIT.
C         NPOLYS  MUST SATISFY NPOLYS .LT.  NPTS  AND  NPOLYS  .GE. 1
C      NPTS  --- (PASSED)
C         THE NUMBER OF DATA POINTS TO BE USED IN THE FIT.
C          NPTS  MUST BE .GE. 1.
C      DIMEN  -- (PASSED)
C         THE NUMBER OF VARIABLES.
C      IWORK  -- (RETURNED)
C         AN INTEGER WORK ARRAY OF LENGTH AT LEAST
C            IF  DEGREE  .GE. 0 THEN
C               4*BINOMIAL( DIMEN + DEGREE , DIMEN )
C                  +( DIMEN )*( DEGREE )
C            ELSE
C               4*BINOMIAL( DIMEN +D,D)+( DIMEN )*D
C         WHERE D IS THE MINIMUM CARDINALITY  OF  A  BASIS  OF  DEGREE
C          DEGREE  SUCH THAT
C            BINOMIAL( DIMEN +ABS( DEGREE ), DIMEN ) .GE.  NPOLYS
C      IWKLEN  - (PASSED)
C         THE LENGTH OF  IWORK
C      IREQD  -- (RETURNED)
C         THE SIZE OF THE INTEGER WORK ARRAY REQUIRED BY  CONSTR   FOR
C         THE FIT SPECIFIED BY THE 4 INPUT PARAMETERS.
C      DREQD  -- (RETURNED)
C         THE SIZE OF THE DOUBLE  PRECISION  WORK  ARRAY  REQUIRED  BY
C          CONSTR  FOR THE FIT SPECIFIED BY THE 4 INPUT PARAMETERS.
C      ERROR  -- (RETURNED)
C       0 IF  NPOLYS ,  DIMEN ,   DEGREE ,   NPTS   AND   IWKLEN   ARE
C         VALID AND CONSISTENT WITH EACH OTHER.
C       1 IF  DEGREE  .GE. 0 BUT THERE IS AN INTERPOLATING MULTINOMIAL
C         OF SMALLER DEGREE OR IF  DEGREE .LT. 0 AND  NPOLYS .GT. NPTS
C       2 IF  DEGREE  .LT. 0 AND  NPOLYS  .LE. 0
C       3 IF  NPTS  .LT. 1 AND/OR  DIMEN  .LT. 1
C       4 IF  IWKLEN  IS TOO SMALL (SET  IWKLEN  TO THE VALUE RETURNED
C         IN  IREQD  TO RESOLVE THIS PROBLEM)
C
C     NOTE THAT  DEGREE ,  NPOLYS ,  PSIWID  AND  ALFL   ARE  RETURNED
C     IN  IWORK (1-4), RESPECTIVELY.
C
C     DATE LAST MODIFIED
C     ---- ---- --------
C     DECEMBER 10, 1984
C     ****************
C
C     ***************
C     BASIZ  COMPUTES THE SIZE OF THE BASIS (AND AUXILIARY SIZES)
C     BASED PRIMARILY UPON THE DEGREE, NUMBER OF  FITTING POINTS,
C     AND THE DIMENSION.
C     ***************
C
      CALL BASIZ(DEGREE,NPTS,DIMEN,NPOLYS,ERROR)
      IF ( ERROR .GE. 2 ) RETURN
      IREQD = 4 * NPOLYS + DEGREE * DIMEN
      IF ( IWKLEN .GE. IREQD ) GO TO 5
         ERROR = 4
         RETURN
 5    NEWSTT = 4 * NPOLYS + 1
C
C     ***************
C     SET UP USEFUL INDEXING ARRAYS
C          IWORK(1) ,..., IWORK(NEWSTT-1)
C     AND
C          IWORK(NEWSTT ,..., IWORK(NEWSTT+DIMEN*DEGREE)
C     ***************
C
      CALL TABLE(DEGREE,DIMEN,NPOLYS,IWORK,IWORK(NEWSTT),ALFL)
      IWORK(1) = DEGREE
      IWORK(2) = NPOLYS
C
C     ***************
C     FORCE  ALFL  TO BE AT LEAST 1 SO THAT DIMENSION STATEMENTS
C     USING  ALFL  DO NOT BOMB.
C     ***************
C
      IF ( ALFL .GT. 1 ) ALFL = ALFL - 1
      IWORK(4) = ALFL
C
C     ***************
C     THE FOLLOWING IS A SECTION OF CODE FOR SETTING UP THE
C     STORAGE  MANAGEMENT OF THE  PSI  ARRAY.  THERE  IS  A
C     COMPLICATED DOVETAILING FORMULA USED TO PACK INFORMATION
C     INTO  PSI  WITHOUT LEAVING GAPS.
C
C     ARRAY          LENGTH
C     -----          ------
C      MAXABS         DIMEN  + 1
C      ALPHA          ALFL
C      C              NPOLYS
C      SUMSQS         NPOLYS
C
C     THE NUMBER OF COLUMNS IN  PSI ,  PSIWID , IS DETERMINED BY
C     PSIWID  =  NPOLYS  + 1 - (THE SMALLEST M SUCH THAT ALPHA(J,M)
C                                   IS NONZERO AND J .GE. NPOLYS)
C     THIS INSURES THAT IF THE USER EXTENDS THE BASIS, ALL THE  PSI
C     REQUIRED WILL CERTAINLY BE STORED
C
C     IF DEGREE( NPOLYS ) .LE. 2 THEN                      (CASE 1)
C       PSIWID  =  NPOLYS
C     ELSE
C       IF K =  DIMEN  THEN                                (CASE 2)
C         PSIWID  =  NPOLYS
C                    -  NEWKJ( 1 , DEGREE(NPOLYS)-1 )  + 1
C       ELSE
C            PSIWID  =  NPOLYS
C                      + 1
C                      - (
C                         THE SMALLER OF
C                           NEWKJ(K+1,DEGREE(NPOLYS)-2)    (CASE 3)
C                         AND
C                           INDEXS(3,NPOLYS)               (CASE 4)
C                        )
C
      IF ( DEGREE .GT. 2 ) GO TO 10
C
C     ***************
C     CASE 1
C     ***************
C
         PSIWID = NPOLYS
      GO TO 40
 10      NPLYT4 = 4 * NPOLYS
C
C     ***************
C     KMXBAS  IS K
C                 NPOLYS
C     ***************
C
         KMXBAS = IWORK(NPLYT4 - 2)
C
         IF ( KMXBAS .NE. DIMEN ) GO TO 20
C
C     ***************
C     CASE 2
C     ***************
C
            PSIWID = NPOLYS - IWORK(4 * NPOLYS - 1)
         GO TO 40
C
C     ***************
C     INDEX =  NEWKJ( K      + 1 , DEGREE(NPOLYS-2) )
C                      NPOLYS
C     ***************
C
 20         INDEX = NPLYT4 + (DEGREE - 3) * DIMEN + KMXBAS + 1
            KJP1D2 = IWORK(INDEX)
C
C     ***************
C     STARTJ  =  INDEXS(3,NPOLYS)
C     ***************
C
            STARTJ = IWORK(NPLYT4 - 1)
            IF ( STARTJ .GT. KJP1D2 ) GO TO 30
C
C     ***************
C     CASE 4
C     ***************
C
               PSIWID = NPOLYS - STARTJ + 1
            GO TO 40
C
C     ***************
C     CASE 3
C     ***************
C
 30            PSIWID = NPOLYS - KJP1D2 + 1
 40   IWORK(3) = PSIWID
      DREQD = 2 * NPOLYS + DIMEN + 1 + NPTS * PSIWID + ALFL
      RETURN
      END
      SUBROUTINE BASIZ(DEGREE,NPTS,DIMEN,NPOLYS,ERROR)
C
      INTEGER TOP,BOT,DEGREE,NPTS,DIMEN,NPOLYS,ERROR,I,ROWLEN
C
C     ***************
C     PURPOSE
C     -------
C
C     IF  DEGREE  .GE. 0 THEN
C       FIND THE SIZE OF A BASIS REQUIRED EITHER TO
C       1) APPROXIMATE THE DATA WITH A POLYNOMIAL OF DEGREE
C          GIVEN BY THE PARAMETER  DEGREE
C       OR TO
C       2) SPAN THE SPACE OF POLYNOMIALS OF DEGREE .LE. THE
C          SMALLEST DEGREE OF POLYNOMIAL WHICH INTERPOLATES THE
C          DATA.
C       IN CASE 1  ERROR  = 0.
C       IN CASE 2  ERROR  = 1.
C     ELSE
C       IF  NPOLYS  .GE. 1 THEN
C          IF NPOLYS .GT. NPTS THEN
C             SET  NPOLYS  =  NPTS , FIND THE SMALLEST DEGREE OF A
C             POLYNOMIAL  WHICH  INTERPOLATES  THE  DATA,  AND SET
C              ERROR  = 1.
C          ELSE
C             FIND THE LARGEST DEGREE  DEGREE  OF A POLYNOMIAL  IN
C             A  BASIS OF  NPOLYS  POLYNOMIALS GENERATED ACCORDING
C             TO OUR ORDERING AND SET  ERROR  = 0.
C       ELSE
C           ERROR  = 2
C
C     THIS SUBROUTINE IS CALLED BY  ALLOT .  IT IS NOT  CALLED  BY
C     THE USER DIRECTLY.
C
C     DATE LAST MODIFIED
C     ---- ---- --------
C     OCTOBER 16, 1984
C     ****************
C
      ERROR = 0
      IF ( NPTS .GE. 1 .AND. DIMEN .GE. 1 ) GO TO 10
         ERROR = 3
         RETURN
C
 10   CONTINUE
      IF ( DEGREE .LT. 0 ) GO TO 30
C
        ROWLEN = 1
        NPOLYS = 1
        TOP = DIMEN - 1
        BOT = 0
        IF ( DEGREE .LT. 1 )  GO TO 30
          DO 20 I=1,DEGREE
            TOP = TOP + 1
            BOT = BOT + 1
            ROWLEN = (ROWLEN*TOP)/BOT
            NPOLYS = NPOLYS + ROWLEN
   20     CONTINUE
C
   30 CONTINUE
      IF ( NPOLYS .GE. 1 ) GO TO 40
            ERROR = 2
            RETURN
   40 CONTINUE
      IF ( NPOLYS .LT. NPTS ) GO TO 50
            NPOLYS = NPTS
            ERROR = 1
   50 CONTINUE
      ROWLEN = 1
      I = 1
      DEGREE = 0
      TOP = DIMEN - 1
      BOT = 0
   60 CONTINUE
      IF ( I .GE. NPOLYS )  GO TO 70
          TOP = TOP + 1
          BOT = BOT + 1
          ROWLEN = (ROWLEN*TOP)/BOT
          I = I + ROWLEN
          DEGREE = DEGREE + 1
          IF ( I .LT. NPOLYS )  GO TO 60
   70 CONTINUE
      RETURN
      END
      SUBROUTINE EVALP(COORD,C,NEPTS,DIMEN,NPOLYS,ALPHA,INDEXS,
     +                 PSI,F,ALFL,MAXABS,DIMP1)
C
      INTEGER DIMEN,NEPTS,NPOLYS,ALFL,DIMP1
      INTEGER JM1,JPRIME,M,P,K,I,J,INDEX
      INTEGER INDEXS(4,NPOLYS)
      DOUBLE PRECISION ALPHA(ALFL),COORD(NEPTS,DIMEN),PSI(NPOLYS)
      DOUBLE PRECISION C(NPOLYS),F(NEPTS),MAXABS(DIMP1)
      DOUBLE PRECISION RUNTOT,RNTOT1
C
C     ***************
C     PURPOSE
C     -------
C
C     THIS SUBROUTINE PERFORMS THE MAIN WORK OF EVALUATING THE
C     FITTING MULTINOMIAL (OR THE INITIAL PORTION OF IT  WHICH
C     IS REQUESTED BY THE SETTING OF  NEPOLS ,  EVLDEG  IN THE
C     CALL TO SUBROUTINE  EVAL .
C
C     THIS SUBROUTINE IS CALLED BY  EVAL .  IT IS  NOT  CALLED
C     DIRECTLY BY THE USER.
C
C     SUBROUTINES  SCALDN  AND  SCALUP  ARE   CALLED  TO SCALE
C     AND TO UNSCALE VALUES.
C
C     THE BODY OF THIS SUBROUTINE FOLLOWS THE EXPLANATION
C     GIVEN IN
C          LEAST SQUARES FITTING USING
C          ORTHOGONAL MULTINOMIALS
C     BY
C          BARTELS AND JEZIORANSKI
C     IN
C          ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE
C
C     DATE LAST MODIFIED
C     ---- ---- --------
C     OCTOBER 16, 1984
C     ****************
C
C     ***************
C     SCALE DOWN THE INDEPENDENT CO-ORDINATES.
C     ***************
C
      DO 10 K = 1,DIMEN
 10      CALL SCALDN(COORD(1,K),NEPTS,MAXABS(K))
C
C     ***************
C     USE THE  BASIS FUNCTION COEFFICIENTS  C  AND RECURRENCE
C     COEFFICIENTS  ALPHA  TO EVALUATE THE FITTED MULTINOMIAL
C     AT THE EVALUATION CO-ORDINATES  COORD .
C     ***************
C
      PSI(1) = 1.0D+00
      DO 40 P = 1,NEPTS
         RNTOT1 = C(1)
         IF ( NPOLYS .EQ. 1 ) GO TO 40
          DO 30 J = 2,NPOLYS
              K = INDEXS(2,J)
              JPRIME = INDEXS(1,J)
              RUNTOT = COORD(P,K) * PSI(JPRIME)
              I = INDEXS(3,J)
              JM1 = J - 1
              DO 20 M = I,JM1
                 INDEX = INDEXS(4,J) + M - I
 20              RUNTOT = RUNTOT - PSI(M) * ALPHA(INDEX)
              PSI(J) = RUNTOT
 30           RNTOT1 = RNTOT1 + C(J) * PSI(J)
 40       F(P) = RNTOT1
C
C     ***************
C     SCALE UP THE DEPENDENT COORDINATES.
C     ***************
C
      CALL SCALUP(F,NEPTS,MAXABS(DIMP1))
      DO 50 K = 1,DIMEN
 50      CALL SCALUP(COORD(1,K),NEPTS,MAXABS(K))
C
      RETURN
      END
      SUBROUTINE GNRTP(DEGREE,ALPHA,PSI,INDEXS,
     +                 NEWKJ,SUMSQS,COORD,NPOLYS,
     +                 DIMEN,NPTS,F,Z,C,PSIWID,WEIGHT,
     +                 ALFL,DIMP1,MAXABS)
C
      INTEGER DEGREE,DIMEN,NPOLYS,NPTS,K,PSIWID,ALFL,P,STTDEG,ONPLYS
      INTEGER ERROR,DIMP1
      INTEGER INDEXS(4,NPOLYS),NEWKJ(DIMEN,DEGREE)
      DOUBLE PRECISION PSI(NPTS,PSIWID),ALPHA(ALFL),F(NPTS)
      DOUBLE PRECISION COORD(NPTS,DIMEN),MAXABS(DIMP1),WEIGHT(NPTS)
      DOUBLE PRECISION Z(NPTS),SUMSQS(NPOLYS),C(NPOLYS)
      DOUBLE PRECISION RUNTOT,RNTOT1
C
C     ***************
C     PURPOSE
C     -------
C
C     THE MULTINOMIAL FIT IS GENERATED INCREMENTALLY, A BASIS  ELEMENT
C     AT A TIME.  THIS SUBROUTINE STARTS THE PROCESS OFF BY SETTING UP
C     THE FIRST BASIS ELEMENT, SCALING THE  DATA,  FINDING  THE  FIRST
C     COEFFICIENT,  AND  INITIALIZING  THE  WORK ARRAY Z.  GNRTP  THEN
C     CALLS  INCDG  IF MORE THAN ONE BASIS ELEMENT IS REQUIRED.
C
C     THIS SUBROUTINE IS CALLED BY  CONSTR .  IT IS NOT CALLED BY  THE
C     USER.
C
C     THIS SUBROUTINE CALLS  SCALPM ,  SCALDN , AND  INCDG .
C
C     DATE LAST MODIFIED
C     ---- ---- --------
C     OCTOBER 16, 1984
C     ****************
C
C     ***************
C     SET UP THE SCALING.
C     ***************
C
      DO 10 K = 1,DIMEN
         CALL SCALPM(COORD(1,K),NPTS,MAXABS(K))
 10      CALL SCALDN(COORD(1,K),NPTS,MAXABS(K))
      CALL SCALPM(F,NPTS,MAXABS(DIMP1))
      CALL SCALDN(F,NPTS,MAXABS(DIMP1))
C
C     ***************
C      SUMSQS (1) = (1,1)
C     C  = (F,1) / (1,1)
C      1
C     ***************
C
      RUNTOT = 0.0D+00
      RNTOT1 = 0.0D+00
      DO 20 P = 1,NPTS
         PSI(P,1) = 1.0D+00
         RNTOT1 = RNTOT1 + WEIGHT(P)
 20      RUNTOT = RUNTOT + F(P) * WEIGHT(P)
      SUMSQS(1) = RNTOT1
      C(1) = RUNTOT / RNTOT1
C
C     ***************
C     Z = F - C
C              1
C     ***************
C
      DO 30 P = 1,NPTS
 30      Z(P) = F(P) - C(1)
C
      IF ( NPOLYS .EQ. 1 ) RETURN
      STTDEG = 1
      ONPLYS = 1
C
      CALL INCDG(DEGREE,ALPHA,PSI,INDEXS,NEWKJ,SUMSQS,
     +           COORD,NPOLYS,DIMEN,NPTS,F,Z,C,PSIWID,
     +           WEIGHT,ALFL,ONPLYS,STTDEG,ERROR)
      RETURN
      END
      SUBROUTINE INCDG(DEGREE,ALPHA,PSI,INDEXS,NEWKJ,
     +                 SUMSQS,COORD,NPOLYS,
     +                 DIMEN,NPTS,F,Z,C,PSIWID,WEIGHT,
     +                 ALFL,ONPLYS,STTDEG,ERROR)
C
      INTEGER JPRIME,P,J,CURDEG,KJ,KJP,L,JPM1,JM1
      INTEGER M,START,JINDEX,JPINDX,Q,J3,J1,J1MJ2,ERROR
      INTEGER J0MJ1,J1M1,STARTA,ONPLYS,ONPP1,STTDEG,INDEX1,INDEX2
      INTEGER DEGREE,NPOLYS,NPTS,DIMEN,PSIWID,ALFL
      DOUBLE PRECISION ALPHA(ALFL),COORD(NPTS,DIMEN),PSI(NPTS,PSIWID)
      DOUBLE PRECISION SUMSQS(NPOLYS),C(NPOLYS),F(NPTS),WEIGHT(NPTS)
      DOUBLE PRECISION Z(NPTS)
      INTEGER INDEXS(4,NPOLYS),NEWKJ(DIMEN,DEGREE)
      DOUBLE PRECISION RUNTOT,RNTOT1,RNTOT2
C
C     ***************
C     PURPOSE
C     -------
C
C     THE MULTINOMIAL FIT IS GENERATED INCREMENTALLY, A BASIS  ELEMENT
C     AT A TIME.  THIS SUBROUTINE CONTINUES THE PROCESS STARTED OFF BY
C      GNRTP .
C
C     THIS SUBROUTINE IS CALLED BY  GNRTP  AND NOT BY THE USER.
C
C
C     DATE LAST MODIFIED
C     ---- ---- --------
C     OCTOBER 16, 1984
C     ****************
C
      ERROR = 0
      IF ( ONPLYS .GE. 1 .AND. STTDEG .GE. 1 ) GO TO 10
         ERROR = 6
      RETURN
 10      IF ( INDEXS(2,ONPLYS) .EQ. DIMEN ) GO TO 20
            CURDEG = STTDEG
         GO TO 30
 20         CURDEG = STTDEG + 1
 30   ONPP1 = ONPLYS + 1
      DO 170 J = ONPP1,NPOLYS
         JPRIME = INDEXS(1,J)
         JINDEX = J - (J - 1) / PSIWID * PSIWID
         JPINDX = JPRIME - (JPRIME - 1) / PSIWID * PSIWID
         KJ = INDEXS(2,J)
         START = INDEXS(3,J)
         M = START
         STARTA = INDEXS(4,J) - START
         IF ( CURDEG .EQ. 1 ) GO TO 100
         KJP = INDEXS(2,JPRIME)
         J1 = NEWKJ(KJ,CURDEG - 1)
C
C     ***************
C     CALCULATE THOSE  ALPHA ( J , M ) THAT CAN BE CALCULATED FROM
C     PREVIOUSLY CALCULATED ALPHAS.
C     ***************
C
         IF ( KJ .LT. KJP ) GO TO 50
C
C     ***************
C     FIRST CALCULATE THOSE BETWEEN  JPP  AND THE END OF 2 ROWS BACK.
C     CALCULATE  ALPHA ( J , JPP )
C     ***************
C
         INDEX1 = INDEXS(4,J)
         ALPHA(INDEX1) = SUMSQS(JPRIME) / SUMSQS(START)
C
         M = START + 1
         J3 = NEWKJ(1,CURDEG - 1) - 1
         IF ( M .GT. J3 ) GO TO 50
C
C     ***************
C     CURDEG .GT. 2 IF CONTROL HAS PASSED THE BRANCHES IN THE 3-RD
C     PREVIOUS AND 8-TH PREVIOUS STATEMENTS.
C     ***************
C
            J1MJ2 = J1 - NEWKJ(KJ,CURDEG - 2)
C
            DO 40 L = M,J3
               Q = J1MJ2 + L
               INDEX1 = STARTA + L
               INDEX2 = INDEXS(4,Q) - INDEXS(3,Q) + JPRIME
 40            ALPHA(INDEX1) = ALPHA(INDEX2) * SUMSQS(JPRIME) /
     +         SUMSQS(L)
C
C     ***************
C     CALCULATE  ALPHA ( J , M ) FOR  M  BETWEEN THE 2
C     RANGES CALCULATED BEFORE USING
C
C        ALPHA ( J , L ) = (X  *  PSI  ,PSI )  / (PSI ,PSI )
C                             K      JP    L         L    L
C                              J
C     ***************
C
         M = J3 + 1
 50      IF ( JPRIME .EQ. J1 ) GO TO 100
            IF ( KJ .EQ. 1 ) GO TO 80
            J1M1 = J1 - 1
            DO 70 L = M,J1M1
               RUNTOT = 0.0D+00
               DO 60 P = 1,NPTS
                  INDEX1 = L - (L - 1) / PSIWID * PSIWID
 60               RUNTOT = RUNTOT + COORD(P,KJ) * PSI(P,JPINDX) *
     +            PSI(P,INDEX1) * WEIGHT(P)
               INDEX1 = STARTA + L
 70            ALPHA(INDEX1) = RUNTOT / SUMSQS(L)
C
C     ***************
C     CALCULATE  ALPHA ( J , M ) FOR  M  BETWEEN
C      NEWKJ ( KJ , CURDEG  - 1)  AND
C      JP  - 1.
C     ***************
C
 80         J0MJ1 = NEWKJ(KJ,CURDEG) - J1
            JPM1 = JPRIME - 1
            DO 90 L = J1,JPM1
               Q = J0MJ1 + L
               INDEX1 = STARTA + L
               INDEX2 = INDEXS(4,Q) - INDEXS(3,Q) + JPRIME
 90            ALPHA(INDEX1) = ALPHA(INDEX2) * SUMSQS(JPRIME) /
     +         SUMSQS(L)
            M = JPRIME
C
C     ***************
C     CALCULATE THE REMAINING  ALPHA ( J , M ) FROM
C
C      ALPHA ( J , L ) = (X   * PSI  ,PSI ) / (PSI ,PSI )
C                          K       JP    L        L    L
C                           J
C     ***************
C
 100     JM1 = J - 1
         DO 120 L = M,JM1
            RUNTOT = 0.0D+00
            DO 110 P = 1,NPTS
               INDEX1 = L - (L - 1) / PSIWID * PSIWID
 110           RUNTOT = RUNTOT + COORD(P,KJ) * PSI(P,JPINDX) *
     +         PSI(P,INDEX1) * WEIGHT(P)
            INDEX1 = STARTA + L
 120        ALPHA(INDEX1) = RUNTOT / SUMSQS(L)
C
C     ***************
C     NOW CALCULATE THE  PSI (P,J),  SUMSQS (J) AND  C (J) USING
C
C                          J-1
C     PSI  = X  * PSI   -  SUM   ALPHA(J,L) * PSI
C        J    K      JP    L=JPP                 L
C
C     SUMSQS  = (PSI ,PSI )
C           J       J    J
C
C     C  = (Z,PSI )
C      J         J
C     ***************
C
 130      RNTOT1 = 0.0D+00
          RNTOT2 = 0.0D+00
          DO 150 P = 1,NPTS
              RUNTOT = COORD(P,KJ) * PSI(P,JPINDX)
              JM1 = J - 1
              DO 140 L = START,JM1
                 INDEX1 = STARTA + L
                 INDEX2 = L - (L - 1) / PSIWID * PSIWID
 140             RUNTOT = RUNTOT - ALPHA(INDEX1) * PSI(P,INDEX2)
              PSI(P,JINDEX) = RUNTOT
              RNTOT1 = RNTOT1 + PSI(P,JINDEX) * PSI(P,JINDEX) *
     +                          WEIGHT(P)
 150          RNTOT2 = RNTOT2 + Z(P) * PSI(P,JINDEX) * WEIGHT(P)
          SUMSQS(J) = RNTOT1
          C(J) = RNTOT2 / RNTOT1
C
C     ***************
C     CALCULATE THE NEW  Z ( P ) AND THE NEW  SSRES  USING
C
C     Z = Z - C  * PSI
C              J      J
C     ***************
C
          DO 160 P = 1,NPTS
 160          Z(P) = Z(P) - C(J) * PSI(P,JINDEX)
 170      IF ( KJ .EQ. DIMEN ) CURDEG = CURDEG + 1
      RETURN
      END
      SUBROUTINE MOVE(OLDARR,NEWARR,START,OLDWID,NEWIDT,COLENG,ERROR)
C
      INTEGER START,OLDWID,NEWIDT,COLENG,BIG,BIG1,LILN,BIGN,I,J
      INTEGER ERROR,JO,JN,OLDWPS,BIG1PS,BIGP1,K
      DOUBLE PRECISION OLDARR(COLENG,OLDWID),NEWARR(COLENG,NEWIDT)
C
C     ***************
C     PURPOSE
C     -------
C
C     MOVE COLUMNS 1 THROUGH   OLDWID   OF   OLDARR   INTO  COLUMNS  1
C     THROUGH  NEWWID  OF  NEWARR  USING
C         ( START  + I) MOD  OLDWID
C     THROUGH
C         ( START  + I) MOD  NEWID
C     FOR
C         I = 0
C     THROUGH
C         I =  OLDWID  - 1
C     THE MOVEMENT STARTS FROM  THE  LARGEST  COLUM  OF   NEWARR   AND
C     PROCEEDS DOWNWARD TO THE SMALLEST (SO THAT  OLDARR  AND  NEWARR
C     CAN BE THE SAME ARRAY).
C
C
C     DATE LAST MODIFIED
C     ---- ---- --------
C     OCTOBER 16, 1984
C     ****************
C
      ERROR = 1
      IF ( OLDWID .LT. 1 .OR. NEWIDT .LT. 1 .OR. NEWIDT .LT. OLDWID)
     +         RETURN
      ERROR = 0
      BIG = START + OLDWID - 1
      BIGN = BIG - (BIG - 1) / NEWIDT * NEWIDT
      LILN = START - (START - 1) / NEWIDT * NEWIDT
      IF ( LILN .GT. BIGN ) GO TO 20
         OLDWPS = OLDWID + START
         DO 10 I = 1,OLDWID
            J = OLDWPS - I
            JO = J - (J - 1) / OLDWID * OLDWID
            JN = J - (J - 1) / NEWIDT * NEWIDT
            DO 10 K = 1,COLENG
 10            NEWARR(K,JN) = OLDARR(K,JO)
         RETURN
 20   BIG1 = NEWIDT - LILN + 1
      BIG1PS = BIG1 + START
      DO 30 I = 1,BIG1
         J = BIG1PS - I
         JO = J - (J - 1) / OLDWID * OLDWID
         JN = J - (J - 1) / NEWIDT * NEWIDT
         DO 30 K = 1,COLENG
 30         NEWARR(K,JN) = OLDARR(K,JO)
      BIGP1 = BIG + 1
      DO 40 I = 1,BIGN
         J = BIGP1 - I
         JO = J - (J - 1) / OLDWID * OLDWID
         JN = J - (J - 1) / NEWIDT * NEWIDT
         DO 40 K = 1,COLENG
 40         NEWARR(K,JN) = NEWARR(K,JO)
      RETURN
      END
      SUBROUTINE RESTRT(PSIWID,DWORK,DREQD,ONPLYS,OPSWID,OLALFL,CSTT,
     +                  SSQSTT,PSISTT,NPTS,DIMEN)
C
      INTEGER DREQD,ONPLYS,OPSWID,OLALFL,CSTT
      INTEGER OSSQST,OPSIST,PSIWID,NPTS,ERROR,DIMEN
      INTEGER I,J,SSQSTT,PSISTT,OCST,START,INDEX1,INDEX2
      DOUBLE PRECISION DWORK(DREQD)
C
C     ***************
C     PURPOSE
C     -------
C
C     THIS SUBROUTINE REARRANGES THE WORK  SPACE  (WITH THE HELP
C     OF SUBROUTINE  MOVE ) IN THE EVENT THAT A FIT OF INCREASED
C     DEGREE IS TO BE MADE.
C
C     CALLED INTERNALLY BY  CONSTR , NOT BY THE USER.
C
C     DATE LAST MODIFIED
C     ---- ---- --------
C     OCTOBER 16, 1984
C     ****************
C
      OCST = 2 + DIMEN + OLALFL
      OSSQST = OCST + ONPLYS
      OPSIST = OSSQST + ONPLYS
      START = ONPLYS - OPSWID
      CALL MOVE(DWORK(OPSIST),DWORK(PSISTT),START,OPSWID,PSIWID,NPTS,
     +          ERROR)
C
      DO 5 J = 1,ONPLYS
         I = ONPLYS - J
         INDEX1 = SSQSTT + I
         INDEX2 = OSSQST + I
 5       DWORK(INDEX1) = DWORK(INDEX2)
      DO 10 J = 1,ONPLYS
         I = ONPLYS - J
         INDEX1 = CSTT + I
         INDEX2 = OCST + I
 10      DWORK(INDEX1) = DWORK(INDEX2)
      RETURN
      END
      SUBROUTINE SCALPM(COORD,NPTS,MAXABS)
C
      INTEGER NPTS,P
      DOUBLE PRECISION MAXABS,A
      DOUBLE PRECISION COORD(NPTS)
C
C     ***************
C     PURPOSE
C     -------
C
C     FIND SCALING PARAMETER(S) FOR THE PROBLEM. IF THE SCALING SCHEME
C     IS CHANGED, ALL FOUR OF THE FOLLOWING WOULD HAVE TO BE CHANGED
C
C     1)  SCALPM  - FIND THE SCALING PARAMETERS
C     2)  SCALDN  - SCALE THE PROBLEM DATA
C     3)  SCALUP  - PERFORM THE INVERSE TRANSFORMATION TO  SCALDN
C     4) THE SCALING OF THE RESIDUALS IN  CONSTR
C
C     THIS SUBROUTINE IS CALLED BY  GNRTP .  IT IS NOT CALLED  BY  THE
C     USER.
C
C     THE SCALING WHICH  IT  DEFINES  MUST  BE  COORDINATED  WITH  THE
C     SCALING  OF RESIDUALS WHICH IS CARRIED OUT TOWARD THE END OF THE
C     SUBROUTINE  CONSTR .  THE SCALING DEFINED  BY  THIS  ROUTINE  IS
C     APPLIED  IN  THE  SECTION OF  CONSTR  JUST MENTIONED, AND BY THE
C     ROUTINES  SCALUP  AND  SCALDN .
C
C     DATE LAST MODIFIED
C     ---- ---- --------
C     OCTOBER 16, 1984
C     ****************
C
      MAXABS = 0.0D+00
      DO 5 P = 1,NPTS
         A = DABS(COORD(P))
 5       IF ( A .GT. MAXABS ) MAXABS = A
      RETURN
      END
      SUBROUTINE SCALDN(COORD,NPTS,MAXABS)
C
      INTEGER NPTS,P
      DOUBLE PRECISION MAXABS
      DOUBLE PRECISION COORD(NPTS)
C
C     ***************
C     PURPOSE
C     -------
C
C     CARRY OUT THE DATA-SCALING WHICH IS DEFINED  BY  THE  SUBROUTINE
C      SCALPM .
C
C     THIS SUBROUTINE IS CALLED BY  GNRTP  AND   EVALP .   IT  IS  NOT
C     CALLED BY THE USER.
C
C     THE SCALING WHICH THIS ROUTINE CARRIES OUT  MUST  BE  CONSISTENT
C     WITH  THE  SCALING  OF RESIDUALS WHICH IS CARRIED OUT TOWARD THE
C     END OF THE SUBROUTINE  CONSTR .
C
C     DATE LAST MODIFIED
C     ---- ---- --------
C     OCTOBER 16, 1984
C     ****************
C
      IF ( MAXABS .EQ. 0.0D+00 ) RETURN
         DO 5 P = 1,NPTS
 5          COORD(P) = COORD(P) / MAXABS
         RETURN
      END
      SUBROUTINE SCALUP(COORD,NPTS,MAXABS)
C
      INTEGER NPTS,P
      DOUBLE PRECISION MAXABS
      DOUBLE PRECISION COORD(NPTS)
C
C     ***************
C     PURPOSE
C     -------
C
C     REMOVE THE DATA-SCALING  WHICH  IS  DEFINED  BY  THE  SUBROUTINE
C      SCALPM  AND APPLIED BY THE SUBROUTINE  SCALDN .
C
C     THIS SUBROUTINE IS CALLED BY  EVALP .   IT IS  NOT CALLED BY THE
C     USER.
C
C     THE UNSCALING WHICH THIS ROUTINE CARRIES OUT MUST BE THE INVERSE
C     OF  THE SCALING OF RESIDUALS WHICH IS CARRIED OUT TOWARD THE END
C     OF THE SUBROUTINE  CONSTR .
C
C     DATE LAST MODIFIED
C     ---- ---- --------
C     OCTOBER 16, 1984
C     ****************
C
      IF ( MAXABS .EQ. 0.0D+00 ) RETURN
         DO 5 P = 1,NPTS
 5          COORD(P) = COORD(P) * MAXABS
         RETURN
      END
      SUBROUTINE TABLE(DEGREE,DIMEN,NPOLYS,INDEXS,NEWKJ,ALFLP1)
C
      INTEGER J,KJ,CURDEG,JPRIME,NWITHK,I,CURM1,RALEN,DIMM1,DIMM2
      INTEGER NPOLYS,DIMEN,DEGREE,ALFLP1,DIMP1
      INTEGER INDEXS(4,NPOLYS),NEWKJ(DIMEN,DEGREE)
C
C     ***************
C     PURPOSE
C     -------
C
C     TABULATE  JP  AND  KJ  FOR EACH  J
C
C     VARIABLES
C     ---------
C
C      ALFLP1  -- (INTEGER) -- (PASSED)
C         THE LENGTH REQUIRED FOR ARRAY  ALPHA , PLUS ONE
C      DEGREE  -- (INTEGER) -- (PASSED)
C         THE DEGREE OF THE POLYNOMIAL TO BE FITTED
C      DIMEN  -- (INTEGER) -- (PASSED)
C         NUMBER OF INDEPENDENT VARIABLES
C      INDEXS  -- (INTEGER, 2-SUBSCRIPT ARRAY) -- (RETURNED)
C          INDEXS (1, J )   IS    JP ,    INDEXS (2, J )   IS     KJ ,
C          INDEXS (3, J )  IS THE FIRST NONZERO RECURRENCE COEFFICIENT
C         IN  ALPHA  AND  INDEXS (4, J ) IS ITS LOCATION IN  ALPHA .
C      NEWKJ  -- (INTEGER, 2-SUBSCRIPT ARRAY) -- (RETURNED)
C          NEWKJ ( K , D ) IS THE FIRST MONOMIAL OF DEGREE  D   HAVING
C          KJ = K .
C      NPOLYS  -- (INTEGER) -- (PASSED)
C         NUMBER  OF  MONOMIALS  OF  DEGREE  .LE.  ORDER  IN   DIMEN
C         INDEPENDENT VARIABLES.
C
C     THIS SUBPROGRAM CAN BE CODED (EXCLUDING THE PART FOR CALCULATING
C      INDEXS (3, J )  AND   INDEXS (4, J )) MENTALLY MORE EFFICIENTLY
C     BUT COMPUTATIONALLY LESS EFFICIENTLY AS
C
C         J = 2
C         DO 5 KJ = 1,DIMEN
C           NEWKJ(KJ,1) = KJ + 1
C           INDEXS(1,J) = 1
C           INDEXS(2,J) = KJ
C           J = J + 1
C      5  CONTINUE
C         DO 10 CURDEG = 2,DEGREE
C           DO 10 KJ = 1,DIMEN
C             JPRIME = NEWKJ(KJ,CURDEG - 1)
C             NEWKJ(KJ,CURDEG) = J
C             NWITHK = COMB(DIMEN + CURDEG - KJ - 1,CURDEG - 1)
C             DO 10 I = 1,NWITHK
C               INDEXS(1,J) = JPRIME
C               INDEXS(2,J) = KJ
C               JPRIME = JPRIME + 1
C               J = J + 1
C     10  CONTINUE
C
C     WHERE COMB(N,KJ) IS N-FACTORIAL / ((N-KJ)-FACTORIAL * KJ-FACTORIAL
C     HERE WE MAKE USE OF THE RECURRENCE RELATIONS
C
C       COMB(DIMEN+CURDEG-2,CURDEG-1)
C
C                        (DIMEN+CURDEG-2)
C          = ----------------------------------------
C            (CURDEG-1)*COMB(DIMEN+CURDEG-3,CURDEG-2)
C
C     AND
C
C       COMB(DIMEN+CURDEG-KJ-1,CURDEG-1)
C
C                              (DIMEN-KJ+1)
C          = ----------------------------------------------
C            (DIMEN+CURDEG-KJ)*COMB(DIMEN+CURDEG-KJ,CURDEG-1)
C
C
C     DATE LAST MODIFIED
C     ---- ---- --------
C     OCTOBER 16, 1984
C     ****************
C
      ALFLP1 = 1
C
C     ***************
C     SET  INDEXS (4,1) TO 1 SO THAT  ALFL - INDEXS (4,1)+1 IS THE
C     NUMBER OF COLUMNS REQUIRED FOR  PSI  FOR  NPOLYS =1  ( ALFL
C     IS DEFINED IN THE MAINLINE TO BE  ALFLP1 -1 IF ALFLP1 .GT. 1
C     AND  ALFLP1  OTHERWISE.
C     ***************
C
      INDEXS(4,1) = 1
C
      IF ( NPOLYS .EQ. 1 ) RETURN
      J = 2
      DO 10 KJ = 1,DIMEN
         NEWKJ(KJ,1) = KJ + 1
         INDEXS(1,J) = 1
         INDEXS(2,J) = KJ
         INDEXS(3,J) = 1
         INDEXS(4,J) = ALFLP1
         ALFLP1 = ALFLP1 + J - 1
         IF ( J.EQ.NPOLYS ) RETURN
10       J = J + 1
      IF ( DEGREE .EQ. 1 ) RETURN
      RALEN = 1
      DIMM1 = DIMEN - 1
      DIMM2 = DIMEN - 2
      DIMP1 = DIMEN + 1
      DO 70 CURDEG = 2,DEGREE
         CURM1 = CURDEG - 1
         RALEN = RALEN * (DIMM2 + CURDEG) / CURM1
         NWITHK = RALEN
         KJ = 1
20       JPRIME = NEWKJ(KJ,CURM1)
         NEWKJ(KJ,CURDEG) = J
         IF ( KJ.EQ.DIMEN ) GO TO 60
            DO 50 I = 1,NWITHK
               INDEXS(1,J) = JPRIME
               INDEXS(2,J) = KJ
C
C     ***************
C     CALCULATE  INDEXS (3, J ),  INDEXS (4, J )
C     ***************
C
               IF ( KJ .LT. INDEXS(2,JPRIME) ) GO TO 30
                  INDEXS(3,J) = INDEXS(1,JPRIME)
                  GO TO 40
 30            INDEXS(3,J) = NEWKJ(1,CURDEG - 1)
 40            INDEXS(4,J) = ALFLP1
               ALFLP1 = ALFLP1 + J - INDEXS(3,J)
               IF ( J .EQ. NPOLYS ) RETURN
C
               JPRIME = JPRIME + 1
 50            J = J + 1
            KJ = KJ + 1
            NWITHK = NWITHK * (DIMP1 - KJ) / (DIMEN + CURDEG - KJ)
            GO TO 20
 60      INDEXS(1,J) = JPRIME
         INDEXS(2,J) = DIMEN
         INDEXS(3,J) = INDEXS(1,JPRIME)
         INDEXS(4,J) = ALFLP1
         ALFLP1 = ALFLP1 + J - INDEXS(3,J)
         IF ( J .EQ. NPOLYS ) RETURN
 70      J = J + 1
      RETURN
      END
 
