API tools faq
paste
Advertisement
Guest User

Algorithm 634 - CONSTR and EVAL - Routines for fitting mult

a guest
Jun 13th, 2018
241
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Fortran 55.33 KB | None | 0 0
raw download clone embed print report
  1. C     ALGORITHM 634 COLLECTED ALGORITHMS FROM ACM.
  2. C     ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.11, NO. 3,
  3. C     SEPT., 1985, P. 201-217; 218-228.
  4. C     PROGRAM GENERA
  5. C
  6. C ====================================================================
  7. C DATA GENERATOR -- DATA GENERATOR -- DATA GENERATOR -- DATA GENERATOR
  8. C ====================================================================
  9. C
  10.       INTEGER FITDEG,DIMEN,NFPOLS,NFPTS,NEPTS,NEPOLS,EVLDEG,TOPS
  11. C
  12. C     ***************
  13. C     A  PRIMITIVE  DATA-GENERATING   PROGRAM   FOR   USE   WITH   THE
  14. C     JEZIORANSKI-BARTELS    SUITE   OF   ROUTINES   FOR   MULTINOMIAL
  15. C     LEAST-SQUARES FITTING.
  16. C
  17. C     THE OUTPUT FORMATS IN (GENDAT) AND IN  (GENEVL)  ARE  CONSISTENT
  18. C     WITH  THE INPUT FORMATS TO BE FOUND IN THE SIMPLE DRIVER PROGRAM
  19. C     WHICH IS INCLUDED WITH THE SUITE.
  20. C
  21. C     THE DATA  STATEMENTS BELOW  WILL CREATE  A PROBLEM REQUIRING THE
  22. C     FIT OF A  3-VARIABLE MULTINOMIAL, CONSISTING OF  8 TERMS, TO 27
  23. C     POINTS OF  GENERATED DATA  AND THEN  THE EVALUATION  OF THE FULL
  24. C     RESULTING MULTINOMIAL AT 5 POINTS.
  25. C
  26. C     DATE LAST MODIFIED
  27. C     ---- ---- --------
  28. C     DECEMBER 14, 1984
  29. C     ****************
  30. C
  31.       DATA FITDEG/-1/
  32.       DATA EVLDEG/-1/
  33.       DATA NFPOLS/8/
  34.       DATA DIMEN/3/
  35.       DATA NFPTS/27/
  36.       DATA NEPTS/5/
  37.       DATA NEPOLS/8/
  38.       DATA TOPS/3/
  39. C
  40.       CALL GENDAT(DIMEN,NFPTS,FITDEG,NFPOLS,TOPS)
  41. C
  42.       CALL GENEVL(DIMEN,NEPTS,EVLDEG,NEPOLS,TOPS)
  43. C
  44.       STOP
  45.       END
  46.       SUBROUTINE GENDAT(DIMEN,NFPTS,FITDEG,NFPOLS,TOPS)
  47. C
  48.       INTEGER DIMEN,NFPTS,FITDEG,NFPOLS,I,J,K,COUNT,TOPS
  49.       DOUBLE PRECISION X,Y,Z,W,F,XSTART,YSTART,ZSTART
  50.       DOUBLE PRECISION XDEL,YDEL,ZDEL,WEIGHT
  51. C
  52. C     ***************
  53. C     A SUBROUTINE TO GENERATE AND  PRINT  DATA  POINTS  FOR  THE
  54. C     SPECIFIC FUNCTION
  55. C
  56. C                          SIN(X)
  57. C               F(X,Y,Z) = ------ * COS(Y) + EXP(Z)
  58. C                            X
  59. C
  60. C     ***************
  61. C
  62.       DATA XSTART /0.1D+00/
  63.       DATA YSTART /0.1D+00/
  64.       DATA ZSTART /0.1D+00/
  65.       DATA XDEL   /0.2D+00/
  66.       DATA YDEL   /0.2D+00/
  67.       DATA ZDEL   /0.2D+00/
  68.       DATA WEIGHT /1.0D+00/
  69. C
  70.       WRITE (6,6000) DIMEN,FITDEG,NFPOLS,NFPTS
  71.       W = WEIGHT
  72.       COUNT = 0
  73.       X = XSTART
  74.       DO 30 I = 1,TOPS
  75.          Y = YSTART
  76.          DO 20 J = 1,TOPS
  77.             Z = ZSTART
  78.             DO 10 K = 1,TOPS
  79.                IF (COUNT .GE. NFPTS)  GO TO 40
  80.                F = (DSIN(X)/X)*DCOS(Y) + DEXP(Z)
  81.                WRITE (6,6010) W,X,Y,Z,F
  82.                Z = Z + ZDEL
  83.                COUNT = COUNT + 1
  84.    10       CONTINUE
  85.             Y = Y + YDEL
  86.    20    CONTINUE
  87.          X = X + XDEL
  88.    30 CONTINUE
  89.    40 CONTINUE
  90.       RETURN
  91. C
  92.  6000 FORMAT(4I5)
  93.  6010 FORMAT(5D14.6)
  94.       END
  95.       SUBROUTINE GENEVL(DIMEN,NEPTS,EVLDEG,NEPOLS,TOPS)
  96. C
  97.       INTEGER DIMEN,NEPTS,I,NEPOLS,EVLDEG,COUNT,TOPS
  98.       DOUBLE PRECISION X,Y,Z,XSTART,YSTART,ZSTART
  99.       DOUBLE PRECISION XDEL,YDEL,ZDEL
  100. C
  101. C     ***************
  102. C     A SUBROUTINE TO GENERATE EVALUATION POINTS.
  103. C     ***************
  104. C
  105.       DATA XSTART /0.15D+00/
  106.       DATA YSTART /0.16D+00/
  107.       DATA ZSTART /0.17D+00/
  108.       DATA XDEL   /0.05D+00/
  109.       DATA YDEL   /0.06D+00/
  110.       DATA ZDEL   /0.07D+00/
  111. C
  112.       COUNT = 0
  113.       WRITE (6,6000) EVLDEG,NEPOLS,NEPTS
  114.       X = XSTART
  115.       DO 30 I = 1,TOPS
  116.          Y = YSTART
  117.          DO 20 J = 1,TOPS
  118.             Z = ZSTART
  119.             DO 10 K = 1,TOPS
  120.                IF (COUNT .GE. NEPTS)  GO TO 40
  121.                WRITE (6,6010) X,Y,Z
  122.                Z = Z + ZDEL
  123.                COUNT = COUNT + 1
  124.    10       CONTINUE
  125.             Y = Y + YDEL
  126.    20    CONTINUE
  127.          X = X + XDEL
  128.    30 CONTINUE
  129.    40 CONTINUE
  130.       RETURN
  131. C
  132.  6000 FORMAT(3I5)
  133.  6010 FORMAT(3D14.6)
  134. C
  135.       END
  136. C     PROGRAM DRIVER
  137. C
  138. C ====================================================================
  139. C   SAMPLE DRIVER -- SAMPLE DRIVER -- SAMPLE DRIVER -- SAMPLE DRIVER
  140. C ====================================================================
  141. C
  142.       INTEGER DIMEN,FITDEG,NFPOLS,NFPTS,EVLDEG,NEPOLS,NEPTS
  143.       INTEGER ERROR,FIWKLN,FDWKLN,EDWKLN,IREQD,DREQD
  144.       INTEGER FITIWK(89)
  145.       DOUBLE PRECISION FITDWK(2201),FITVLS(125),FITCDS(375),WTS(125)
  146.       DOUBLE PRECISION EVLDWK(125), EVLVLS(20), EVLCDS(60), RESIDS(125)
  147. C
  148. C     ***************
  149. C     A  SIMPLE  DRIVER  PROGRAM  TO  READ  TEST  DATA  AND  CALL  THE
  150. C     JEZIORANSKI-BARTELS  SUITE  OF MULTINOMIAL LEAST-SQUARES FITTING
  151. C     ROUTINES.
  152. C
  153. C     THE DATA AND ARRAY DECLARATIONS  ARE  CONSISTENT  WITH  PROBLEMS
  154. C     INVOLVING  UP  TO  3  VARIABLES,  125   FITTING  POINTS,  AND 20
  155. C     EVALUATION POINTS USING 20 BASIS ORTHOGONAL MULTINOMIALS IN BOTH
  156. C     THE  FIT AND THE EVALUATION.
  157. C
  158. C     DATE LAST MODIFIED
  159. C     ---- ---- --------
  160. C     DECEMBER 15, 1984
  161. C     ****************
  162. C
  163.       DATA FDWKLN/2201/
  164.       DATA EDWKLN/125/
  165.       DATA FIWKLN/89/
  166. C
  167. C     ***************
  168. C     READ FITTING DATA
  169. C     ***************
  170. C
  171.       READ (5,5000) DIMEN,FITDEG,NFPOLS,NFPTS
  172.       WRITE (6,6000) DIMEN,FITDEG,NFPOLS,NFPTS
  173.       CALL INFIT(DIMEN,NFPTS,FITCDS,FITVLS,WTS)
  174. C
  175. C     ***************
  176. C     COMPUTE THE LEAST-SQUARES FIT
  177. C     ***************
  178. C
  179.       CALL CONSTR(DIMEN,FITDEG,NFPOLS,NFPTS,FITCDS,NFPTS,
  180.      +            FITVLS,WTS,RESIDS,.TRUE.,ERROR,FITIWK,
  181.      +            FITDWK,FIWKLN,FDWKLN,IREQD,DREQD)
  182. C
  183. C     ***************
  184. C     PRINT RESIDUALS AT THE FITTING POINTS
  185. C
  186. C     THE USER COULD CHECK  IREQD,DREQD  AT THIS POINT TO SEE
  187. C     THE NUMBER OF LOCATIONS WHICH WERE ACTUALLY REQUIRED IN
  188. C     THE ARRAYS  FITIWK,FITDWK.
  189. C     ***************
  190. C
  191.       WRITE (6,6010)
  192.       CALL OUTRES(NFPTS,RESIDS)
  193. C
  194. C     ***************
  195. C     READ POINTS OF EVALUATION
  196. C     ***************
  197. C
  198.       READ (5,5000) EVLDEG,NEPOLS,NEPTS
  199.       WRITE (6,6020) EVLDEG,NEPOLS,NEPTS
  200.       CALL INEVL(DIMEN,NEPTS,EVLCDS)
  201. C
  202. C     ***************
  203. C     EVALUATE THE FITTING MULTINOMIAL
  204. C     ***************
  205. C
  206.       CALL EVAL(DIMEN,EVLDEG,NEPOLS,NEPTS,EVLCDS,EVLVLS,
  207.      +          ERROR,FITIWK,FITDWK,FIWKLN,FDWKLN,EVLDWK,EDWKLN)
  208. C
  209. C     ***************
  210. C     PRINT OUT THE ARRAY OF MULTINOMIAL VALUES
  211. C     ***************
  212. C
  213.       CALL OUTEVL(DIMEN,NEPTS,NEPOLS,EVLCDS,EVLVLS)
  214. C
  215.       STOP
  216. C
  217. C     ***************
  218. C     FORMATS
  219. C     ***************
  220. C
  221.  5000 FORMAT(5I5)
  222.  6000 FORMAT(//29H MULTINOMIAL FITTING PROBLEM.//
  223.      +       10H DIMEN  = ,I5/
  224.      +       10H FITDEG = ,I5/
  225.      +       10H NFPOLS = ,I5/
  226.      +       10H NFPTS  = ,I5)
  227.  6010 FORMAT(//18H FITTING COMPLETE./13H RESIDUALS...
  228.      +       //4X,1HI,5X,8HRESIDUAL)
  229.  6020 FORMAT(//24H MULTINOMIAL EVALUATION.//
  230.      +       10H EVLDEG = ,I5/
  231.      +       10H NEPOLS = ,I5/
  232.      +       10H NEPTS  = ,I5)
  233. C
  234.       END
  235.       SUBROUTINE INEVL(DIMEN,NEPTS,EVLCDS)
  236. C
  237.       INTEGER DIMEN,NEPTS,I,J
  238.       DOUBLE PRECISION EVLCDS(NEPTS,DIMEN)
  239. C
  240. C     ***************
  241. C     SUBROUTINE TO READ THE ARRAY OF EVALUATION POINTS.
  242. C     ***************
  243. C
  244.       DO 100 I=1,NEPTS
  245.         READ (5,5000) (EVLCDS(I,J),J=1,DIMEN)
  246.   100 CONTINUE
  247.       RETURN
  248.  5000 FORMAT(4D14.6)
  249.       END
  250.       SUBROUTINE INFIT(DIMEN,NFPTS,FITCDS,FITVLS,WTS)
  251. C
  252.       INTEGER NFPTS,DIMEN,I,J
  253.       DOUBLE PRECISION FITCDS(NFPTS,DIMEN),FITVLS(NFPTS),WTS(NFPTS)
  254. C
  255. C     ***************
  256. C     SUBROUTINE TO READ THE FITTING DATA.
  257. C     ***************
  258. C
  259.       WRITE (6,6000)
  260.       DO 100 I=1,NFPTS
  261.         READ (5,5000) WTS(I),(FITCDS(I,J),J=1,DIMEN),FITVLS(I)
  262.         WRITE (6,6010) I,WTS(I),(FITCDS(I,J),J=1,DIMEN),FITVLS(I)
  263.   100 CONTINUE
  264.       RETURN
  265.  5000 FORMAT(5D14.6)
  266.  6000 FORMAT(/8H DATA.../
  267.      +       4X,1HI,5X,6HWEIGHT,19X,11HCOORDINATES,17X,11HDATA VALUES)
  268.  6010 FORMAT(I5,5D14.6)
  269.       END
  270.       SUBROUTINE OUTEVL(DIMEN,NEPTS,NEPOLS,EVLCDS,EVLVLS)
  271. C
  272.       INTEGER DIMEN,NEPTS,NEPOLS,I,J
  273.       DOUBLE PRECISION EVLCDS(NEPTS,DIMEN),EVLVLS(NEPTS)
  274. C
  275. C     ***************
  276. C     SUBROUTINE TO PRINT OUT THE RESULTS OF THE EVALUATION.
  277. C     ***************
  278. C
  279.       WRITE (6,6000)
  280.       DO 100 I=1,NEPTS
  281.         WRITE (6,6010) I,(EVLCDS(I,J),J=1,DIMEN),EVLVLS(I)
  282.   100 CONTINUE
  283.       RETURN
  284.  6000 FORMAT(/8H DATA.../
  285.      +       4X,1HI,22X,11HCOORDINATES,14X,5HVALUE)
  286.  6010 FORMAT(I5,4D14.6)
  287.       END
  288.       SUBROUTINE OUTRES(NFPTS,RESIDS)
  289. C
  290.       INTEGER NFPTS,I
  291.       DOUBLE PRECISION RESIDS(NFPTS)
  292. C
  293. C     ***************
  294. C     SUBROUTINE TO PRINT OUT THE RESIDUALS FROM THE FIT.
  295. C     ***************
  296. C
  297.       DO 100 I=1,NFPTS
  298.         WRITE (6,6000) I,RESIDS(I)
  299.   100 CONTINUE
  300.       RETURN
  301. C
  302.  6000 FORMAT(I5,D14.6)
  303.       END
  304. C
  305. C ====================================================================
  306. C   CONSTRUCT FIT -- CONSTRUCT FIT -- CONSTRUCT FIT -- CONSTRUCT FIT
  307. C ====================================================================
  308. C
  309.       SUBROUTINE CONSTR(DIMEN,FITDEG,NFPOLS,NFPTS,
  310.      +                  FITCDS,NCROWS,FITVLS,WTS,
  311.      +                  RESIDS,NEWFIT,ERROR,FITIWK,
  312.      +                  FITDWK,FIWKLN,FDWKLN,IREQD,DREQD)
  313. C
  314.       INTEGER NFPOLS,ONPLYS,FITDEG,NFPTS,DIMEN,OLDEG,FIWKLN,FDWKLN
  315.       INTEGER ERROR,IREQD,DREQD,OPSWID,OLALFL,INDSTT,P,DIMP1,NCROWS
  316.       INTEGER NEWSTT,MAXSTT,ALFSTT,PSISTT,CSTT,SSQSTT,PSIWID,ALFL
  317.       INTEGER FITIWK(FIWKLN)
  318.       DOUBLE PRECISION FITDWK(FDWKLN),FITCDS(NCROWS,DIMEN)
  319.       DOUBLE PRECISION FITVLS(NFPTS),RESIDS(NFPTS)
  320.       DOUBLE PRECISION WTS(NFPTS)
  321.       DOUBLE PRECISION SCALE
  322.       LOGICAL NEWFIT
  323. C
  324. C     ***************
  325. C     PURPOSE
  326. C     -------
  327. C
  328. C     THIS SUBROUTINE CONSTRUCTS A LEAST-SQUARES  MULTINOMIAL  FIT  TO
  329. C     GIVEN DATA USING A BASIS OF ORTHOGONAL MULTINOMIALS.
  330. C
  331. C     THE DATA FOR THE FIT IS GIVEN IN THE ARRAYS  FITCDS ,   FITVLS ,
  332. C     AND   WTS .    FITCDS  IS A DOUBLE-PRECISION MATRIX, EACH ROW OF
  333. C     WHICH CONTAINS  AN  OBSERVATION  POINT  (ORDERED  COLLECTION  OF
  334. C     VARIABLE  VALUES).    FITVLS   IS  A  DOUBLE-PRECISION,  SINGLY-
  335. C     INDEXED ARRAY,  EACH  ELEMENT  OF  WHICH  CONTAINS  AN  OBSERVED
  336. C     FUNCTION VALUE CORRESPONDING TO AN  OBSERVATION POINT.   WTS  IS
  337. C     A DOUBLE-PRECISION, SINGLY-INDEXED ARRAY, EACH ELEMENT OF  WHICH
  338. C     IS A NONNEGATIVE WEIGHT FOR THE CORRESPONDING OBSERVATION.
  339. C
  340. C     THE FIT WHICH IS PRODUCED IS A MULTINOMIAL EXPRESSED IN THE FORM
  341. C
  342. C      C  PSI  (X ,...,X     ) +...+ C       PSI       (X ,...,X     )
  343. C       1    1   1      DIMEN         NFPOLS    NFPOLS   1      DIMEN
  344. C
  345. C     WHERE THE VALUE OF  NFPOLS  WILL BE AS GIVEN (IF  FITDEG .LT. 0)
  346. C     OR  AS  COMPUTED  BY  CONSTR  TO GIVE A FULL-DEGREE FIT (IN CASE
  347. C      FITDEG  IS SPECIFIED .GE. 0).  THE ELEMENTS
  348. C
  349. C         PSI  (X ,...,X     )
  350. C            K   1      DIMEN
  351. C
  352. C     FORM A BASIS FOR  THE  MULTINOMIALS  WHICH  IS  ORTHOGONAL  WITH
  353. C     RESPECT TO THE WEIGHTS AND OBSERVATION POINTS.
  354. C
  355. C     THE EXTENT OF THE FIT CAN BE SPECIFIED IN ONE OF TWO WAYS.
  356. C         IF THE PARAMETER FITDEG IS SET .GE. 0, THEN A COMPLETE BASIS
  357. C         FOR THE MULTINOMIALS OF DEGREE =  FITDEG  WILL BE USED.  (AN
  358. C         ERROR WILL BE  FLAGGED  IF  THIS  WILL  REQUIRE  MORE  BASIS
  359. C         MULTINOMIALS  THAN  THE  NUMBER  OF  DATA  POINTS WHICH WERE
  360. C         GIVEN.)
  361. C         IF THE PARAMETER  FITDEG  IS .LT.  0, THEN  NFPOLS  WILL  BE
  362. C         TAKEN AS THE COUNT OF THE NUMBER OF BASIS MULTINOMIALS TO BE
  363. C         USED FOR A PARTIAL-DEGREE FIT.  (AN ERROR WILL BE FLAGGED IF
  364. C          NFPOLS  .LT. 0.)
  365. C
  366. C     NOTE, THE CALL TO  CONSTR  WITH  NEWFIT  = .TRUE.  CAN  BE  MADE
  367. C         WITH  THE  PARAMETERS  SET  FOR  THE  MAXIMUM  FIT  DESIRED.
  368. C         SEVERAL SUBSEQUENT CALLS TO  CONSTR  WITH  NEWFIT  = .FALSE.
  369. C         CAN  BE  MADE WITH SMALLER VALUES OF  FITDEG  OR  NFPOLS  AS
  370. C         MAY BE DESIRED TO OBTAIN A PARTIAL FIT.
  371. C
  372. C     VARIABLES
  373. C     ---------
  374. C
  375. C      DIMEN  -- (INTEGER) -- (PASSED)
  376. C         THE NUMBER OF VARIABLES.
  377. C      FITDEG  - (INTEGER) -- (PASSED/RETURNED)
  378. C         IGNORED IF .LT. 0.
  379. C         IF DEGREE .GE. 0 THEN  DEGREE  IS  CHECKED  AGAINST  NFPTS .
  380. C         THE VALUE OF  DEGREE  WILL BE REDUCED IF THERE IS A BASIS OF
  381. C         MULTINOMIALS, ALL OF  DEGREE  .LE. DEGREE ,  OF  CARDINALITY
  382. C          NFPTS .  SEE  ERROR  BELOW.
  383. C      NFPOLS  - (INTEGER) -- (PASSED/RETURNED)
  384. C         IGNORED IF  DEGREE  .GE. 0.
  385. C         IF  DEGREE .LT. 0 THEN THE VALUE OF NFPOLS WILL BE TAKEN  AS
  386. C         THE SIZE OF THE BASIS OF MULTINOMIALS TO BE USED IN THE FIT.
  387. C          NFPOLS  MUST SATISFY  NFPOLS  .LT. NFPTS  AND NFPOLS .GE. 1
  388. C         SEE  ERROR  BELOW.
  389. C      NFPTS  --- (INTEGER) -- (PASSED)
  390. C         THE NUMBER OF DATA POINTS TO BE USED IN THE FIT.
  391. C          NFPTS  MUST BE .GE. 1.  SEE  ERROR  BELOW.
  392. C      FITCDS  -- (DOUBLE-PRECISION, 2-SUBSCRIPT ARRAY) -- (PASSED)
  393. C          FITCDS (P,K) IS THE VALUE OF THE K-TH VARIABLE AT THE  P-TH
  394. C         DATA POINT.
  395. C      NCROWS  -- (INTEGER) -- (PASSED)
  396. C         THE ROW DIMENSION  DECLARED  FOR   FITCDS   IN  THE  CALLING
  397. C         PROGRAM.
  398. C      FITVLS  -- (DOUBLE-PRECISION, 1-SUBSCRIPT ARRAY) -- (PASSED)
  399. C          FITVLS (P) IS THE OBSERVED FUNCTION VALUE OF THE P-TH  DATA
  400. C         POINT.
  401. C      WTS  ----- (DOUBLE-PRECISION, 1-SUBSCRIPT ARRAY) -- (PASSED)
  402. C          WTS (P) IS THE WEIGHT ATTACHED TO THE P-TH DATA POINT.
  403. C      RESIDS  -- (DOUBLE-PRECISION, 1-SUBSCRIPT ARRAY) -- (RETURNED)
  404. C          RESIDS (P) IS THE DIFFERENCE BETWEEN THE FITTED FUNCTION AT
  405. C         POINT P AND  FITVLS (P).
  406. C      NEWFIT  -- (LOGICAL) -- (PASSED)
  407. C         A  LOGICAL FLAG.  IF  NEWFIT =.TRUE., THEN THIS IS THE FIRST
  408. C         FIT TO BE CARRIED OUT WITH THE DATA TO BE FOUND IN THE OTHER
  409. C         PARAMETERS TO   CONSTR ,  AND  SPACE  FOR A  FIT  IS  TO  BE
  410. C         ALLOCATED. IF  NEWFIT =.FALSE., THEN A FIT OF ANOTHER DEGREE
  411. C         CAN BE CONSTRUCTED IN THE SPACE ALLOCATED ON A PREVIOUS CALL
  412. C         WITH THE SAME DATA, AND CERTAIN INITIALIZATION STEPS ARE BY-
  413. C         PASSED.
  414. C      ERROR  -- (INTEGER) -- (RETURNED)
  415. C       0 IF  NFPOLS , DIMEN ,   DEGREE ,   NFPTS   AND   IWKLEN   ARE
  416. C         VALID AND CONSISTENT WITH EACH OTHER.
  417. C       1 IF  DEGREE  .GE. 0 BUT THERE IS AN INTERPOLATING MULTINOMIAL
  418. C         OF SMALLER DEGREE OR IF DEGREE .LT. 0 AND NFPOLS .GT. NFPTS.
  419. C       2 IF  DEGREE  .LT. 0 AND  NFPOLS  .LE. 0.
  420. C       3 IF  NFPTS  .LT. 1 AND/OR  DIMEN  .LT. 1.
  421. C       4 IF  IWKLEN  AND/OR  DWKLEN  IS TOO SMALL.  (SET  IWKLEN   TO
  422. C         THE VALUE RETURNED IN  IREQD , AND SET  DWKLEN  TO THE VALUE
  423. C         RETURNED IN  DREQD  TO RESOLVE THIS PROBLEM.)
  424. C       5 NEWFIT  = .FALSE. BUT  ONPLYS  .GE.  NFPOLS.
  425. C       6 ERROR RETURN FROM  INCDG .  THERE IS  NO  MORE ROOM  IN  THE
  426. C         FITDWK  AND/OR  FITIWK  ARRAYS  TO INCLUDE MORE TERMS IN THE
  427. C         FIT.
  428. C      FITIWK  -- (INTEGER, 1-SUBSCRIPT ARRAY) -- (RETURNED)
  429. C         AN INTEGER WORK ARRAY OF LENGTH  FIWKLN .  UPON RETURN  FROM
  430. C         A CALL TO  CONSTR  WITH  NEWFIT  SET .TRUE., SOME  DIMENSION
  431. C         AND ARRAY-LENGTH INFORMATION WILL BE INSERTED.  UPON  RETURN
  432. C         FROM A CALL TO  CONSTR  WITH  NEWFIT  SET .FALSE.,  DETAILED
  433. C         INDEXING INFORMATION (LOCATION OF COEFFICIENTS IN   FITDWK ,
  434. C         ETC.) IS INSERTED.
  435. C      FITDWK  -- (DOUBLE-PRECISION, 1-SUBSCRIPT ARRAY) -- (RETURNED)
  436. C         A DOUBLE PRECISION ARRAY  OF LENGTH   FDWKLN .  UPON  RETURN
  437. C         FROM  CONSTR  WITH  NEWFIT  SET .FALSE.,  THE  FULL  DETAILS
  438. C         OF THE REQUESTED FIT (COEFFICIENTS, ETC.) WILL BE INSERTED.
  439. C      FIWKLN  -- (INTEGER) -- (PASSED)
  440. C         THE LENGTH OF THE ARRAY  FITIWK .
  441. C      FDWKLN  -- (INTEGER) -- (PASSED)
  442. C         THE LENGTH OF THE ARRAY  FITDWK .
  443. C      IREQD  -- (INTEGER) -- (PASSED)
  444. C         THE LENGTH WHICH THE ARRAY  FITIWK  REALLY NEEDS TO BE.
  445. C      DREQD  -- (INTEGER) -- (PASSED)
  446. C         THE LENGTH WHICH THE ARRAY  FITDWK  REALLY NEEDS TO BE.
  447. C
  448. C
  449. C     NOTE, THE 10 AND 70  LOOPS  (I.E.  THE  LOOPS  FOR  SCALING  THE
  450. C     RESIDUALS)  DEPEND  ON  THE  SCALING  SCHEME USED.  THE RESIDUAL
  451. C     SCALING  MUST  BE  CONSISTENT  WITH  THAT  DEFINED BY   SCALPM ,
  452. C      SCALDN , AND  SCALUP .
  453. C
  454. C      CONSTR   CALLS  ALLOT ,  RESTRT ,  INCDG , AND  GNRTP
  455. C
  456. C
  457. C     DATE LAST MODIFIED
  458. C     ---- ---- --------
  459. C     OCTOBER 16, 1984
  460. C     ****************
  461. C
  462.       DIMP1 = DIMEN + 1
  463. C
  464. C     ***************
  465. C     TEST IF THIS IS A NEW FIT OR AN ADDITION TO THE OLD ONE.
  466. C     ***************
  467. C
  468.       IF ( NEWFIT ) GO TO 20
  469. C
  470. C     ***************
  471. C     THE SECTION BELOW IS FOR AN ADDITION TO A PREVIOUS FIT.
  472. C     OBTAIN THE VALUES OF SAVED PARAMETERS.
  473. C     ***************
  474. C
  475.          OLDEG = FITIWK(1)
  476.          ONPLYS = FITIWK(2)
  477.          OPSWID = FITIWK(3)
  478.          OLALFL = FITIWK(4)
  479.          SCALE = FITDWK(DIMP1)
  480.          SCALE = SCALE * SCALE
  481. C
  482.          IF ( NFPOLS .GT. ONPLYS ) GO TO 10
  483.             ERROR = 5
  484.             RETURN
  485.    10    CONTINUE
  486.    20 CONTINUE
  487. C
  488. C     ***************
  489. C     NEW FITTING PROBLEMS BEGIN HERE.  OLD FITTING PROBLEMS
  490. C     PROCEED INTO THIS SECTION FROM ABOVE.
  491. C     ***************
  492. C
  493.       CALL ALLOT(FITDEG,NFPOLS,NFPTS,DIMEN,FITIWK,FIWKLN,IREQD,DREQD,
  494.      +           ERROR)
  495.       IF ( ERROR .GE. 2 ) RETURN
  496. C
  497.       IF ( FDWKLN .GE. DREQD ) GO TO 30
  498.          ERROR = 4
  499.          RETURN
  500.    30 CONTINUE
  501. C
  502.       PSIWID = FITIWK(3)
  503.       ALFL = FITIWK(4)
  504.       INDSTT = 1
  505.       NEWSTT = 4 * NFPOLS + INDSTT
  506.       MAXSTT = 1
  507.       ALFSTT = MAXSTT + DIMP1
  508.       CSTT = ALFSTT + ALFL
  509.       SSQSTT = CSTT + NFPOLS
  510.       PSISTT = SSQSTT + NFPOLS
  511. C
  512.       IF ( NEWFIT ) GO TO 50
  513. C
  514. C     ***************
  515. C     FOR A  COMPLETELY NEW FIT,  SKIP TO 50.  IF THIS IS
  516. C     AN ADDITION TO A PREVIOUS FIT, A PREVIOUSLY EXISTING
  517. C     SCALING MUST BE RESTORED TO THE RESIDUALS.
  518. C     ***************
  519. C
  520.          DO 40 P = 1,NFPTS
  521.             RESIDS(P) = RESIDS(P) / SCALE
  522.    40    CONTINUE
  523. C
  524. C     ***************
  525. C     RE-ARRANGE AND SHUFFLE INFORMATION IN THE WORKING ARRAYS.
  526. C     ***************
  527. C
  528.          CALL RESTRT(PSIWID,FITDWK,DREQD,ONPLYS,OPSWID,OLALFL,CSTT,
  529.      +               SSQSTT,PSISTT,NFPTS,DIMEN)
  530. C
  531. C     ***************
  532. C     PRODUCE THE AUGMENTED FIT.
  533. C     ***************
  534. C
  535.          CALL INCDG(FITDEG,FITDWK(ALFSTT),FITDWK(PSISTT),FITIWK(INDSTT),
  536.      +              FITIWK(NEWSTT),FITDWK(SSQSTT),FITCDS,
  537.      +              NFPOLS,DIMEN,NFPTS,FITVLS,RESIDS,
  538.      +              FITDWK(CSTT),PSIWID,WTS,ALFL,ONPLYS,OLDEG,ERROR)
  539.          GO TO 60
  540. C
  541.    50    CONTINUE
  542. C
  543. C     ***************
  544. C     THE CODE BELOW IS FOR A COMPLETELY NEW FIT.  ARRANGE PARAMETERS,
  545. C     INDICES, AND STORAGE ALLOCATION IN THE WORK ARRAYS.  SCALE THE
  546. C     DATA AND RESIDUALS, AND THEN PRODUCE THE FIT.
  547. C     ***************
  548. C
  549.          CALL GNRTP(FITDEG,FITDWK(ALFSTT),
  550.      +              FITDWK(PSISTT),FITIWK(INDSTT),
  551.      +              FITIWK(NEWSTT),FITDWK(SSQSTT),
  552.      +              FITCDS,NFPOLS,DIMEN,NFPTS,FITVLS,RESIDS,
  553.      +              FITDWK(CSTT),PSIWID,WTS,ALFL,DIMP1,FITDWK(MAXSTT))
  554.          SCALE = FITDWK(DIMEN + 1)
  555.          SCALE = SCALE * SCALE
  556. C
  557.    60 CONTINUE
  558. C
  559. C     ***************
  560. C     UNSCALE THE RESIDUALS FOR THE BENEFIT OF THE USER.
  561. C     ***************
  562. C
  563.       DO 70 P = 1,NFPTS
  564.         RESIDS(P) = RESIDS(P) * SCALE
  565.    70 CONTINUE
  566.       RETURN
  567.       END
  568. C
  569. C ====================================================================
  570. C     EVALUATE FIT -- EVALUATE FIT -- EVALUATE FIT -- EVALUATE FIT
  571. C ====================================================================
  572. C
  573.       SUBROUTINE EVAL(DIMEN,EVLDEG,NEPOLS,NEPTS,EVLCDS,EVLVLS,
  574.      +                ERROR,FITIWK,FITDWK,FIWKLN,FDWKLN,EVLDWK,EDWKLN)
  575. C
  576.       INTEGER FIWKLN,FDWKLN,NEPOLS,NEPTS,DIMEN,ERROR,MAXSTT,ALFSTT,CSTT
  577.       INTEGER GBASIZ,ALFL,DIMP1,EVLDEG,TOP,BOT,CURDEG,EDWKLN
  578.       INTEGER FITIWK(FIWKLN)
  579.       DOUBLE PRECISION FITDWK(FDWKLN),EVLDWK(EDWKLN),EVLCDS(NEPTS,DIMEN)
  580.       DOUBLE PRECISION EVLVLS(NEPTS)
  581. C
  582. C     ***************
  583. C     PURPOSE
  584. C     -------
  585. C
  586. C     THIS SUBROUTINE  EVALUATES  THE  LEAST-SQUARES  MULTINOMIAL  FIT
  587. C     WHICH HAS BEEN PREVIOUSLY PRODUCED BY  CONSTR .  EITHER THE FULL
  588. C     MULTINOMIAL AS PRODUCED MAY BE EVALUATED,  OR  ONLY  AN  INITIAL
  589. C     SEGMENT THEREOF.  AS IN THE CASE WITH  CONSTR , IT IS POSSIBLE
  590. C     (1) TO SPECIFY MULTINOMIALS OF A FULL GIVEN DEGREE, OR
  591. C     (2) TO SPECIFY  THE  NUMBER  OF  ORTHOGONAL  BASIS  ELEMENTS  TO
  592. C         ACHIEVE A PARTIAL-DEGREE FIT.
  593. C
  594. C     IN CASE (1), THE  DESIRED  DEGREE  IS  GIVEN  AS  THE  VALUE  OF
  595. C      EVLDEG   (WHICH  MUST  BE  NONNEGATIVE AND NOT GREATER THAN THE
  596. C     VALUE USED FOR  FITDEG  IN  CONSTR ), AND THE PARAMETER  NEPOLS
  597. C     WILL  BE  SET  BY  EVAL  TO SPECIFY THE NUMBER OF BASIS ELEMENTS
  598. C     REQUIRED. IF  EVLDEG  .LT.  FITDEG   IS  GIVEN,  THEN  ONLY  THE
  599. C     INITIAL  PORTION OF THE FITTING MULTINOMIAL (OF DEGREE  EVLDEG )
  600. C     WILL BE EVALUATED.
  601. C
  602. C     IN CASE (2),  EVLDEG  IS TO BE SET NEGATIVE, IN WHICH  CASE  THE
  603. C     VALUE  OF   NEPOLS  (WHICH MUST BE POSITIVE AND NOT GREATER THAN
  604. C     THE VALUE USED FOR   NFPOLS   IN   CONSTR )  WILL  BE  TAKEN  AS
  605. C     DEFINING  THE  INITIAL  PORTION OF THE FITTING MULTINOMIAL TO BE
  606. C     EVALUATED.
  607. C
  608. C     IF  NEPOLS  =  NFPOLS  (WITH  EVLDEG  .LT.  0),  OR   EVLDEG   =
  609. C     FITDEG  (WITH   EVLDEG   .GT.  0),  THEN  THE  FULL  MULTINOMIAL
  610. C     GENERATED BY  CONSTR  WILL BE EVALUATED.
  611. C
  612. C     THE  EVALUATION  WILL  TAKE  PLACE  FOR  EACH  OF   THE   POINTS
  613. C     (COLLECTION  OF  VARIABLE  VALUES)  GIVEN AS A ROW OF THE MATRIX
  614. C      EVLCDS .  THE  VALUES  PRODUCED  FROM  THE  FULL,  OR  PARTIAL,
  615. C     MULTINOMIAL WILL BE PLACED IN THE ARRAY  EVLVLS .
  616. C
  617. C     VARIABLES
  618. C     ---------
  619. C
  620. C      DIMEN  -- (INTEGER) -- (PASSED)
  621. C         THE NUMBER OF VARIABLES.
  622. C      EVLDEG  -- (INTEGER) -- (PASSED)
  623. C         IF  EVLDEG  .LT. 0, THEN THIS PARAMETER WILL BE IGNORED.
  624. C         IF  EVLDEG  .GE. 0, THEN THE VALUE OF  EVLDEG  MUST  SATISFY
  625. C          EVLDEG  .LE. (THE  DEGREE  OF THE APPROXIMATING MULTINOMIAL
  626. C         GENERATED IN  CONSTR ).  IN THIS CASE  EVLDEG  WILL  SPECIFY
  627. C         THE DEGREE OF THE INITIAL PORTION OF THE FITTING MULTINOMIAL
  628. C         TO BE EVALUATED.
  629. C      NEPOLS  -- (INTEGER) -- (PASSED/RETURNED)
  630. C         IF  EVLDEG  .GE. 0, THEN THIS PARAMETER WILL BE IGNORED.
  631. C         IF EVLDEG .LT. 0, THEN THE PARTIAL MULTINOMIAL INVOLVING THE
  632. C         FIRST   NEPOLS  ORTHOGONAL BASIS FUNCTIONS WILL BE EVALUATED
  633. C         AT THE POINTS GIVEN BY  EVLCDS .  THE RESULTING VALUES  WILL
  634. C         BE STORED IN  EVLVLS .
  635. C         THE VALUE OF NEPOLS MUST BE .GE. 1 AND .LE. (THE SIZE OF THE
  636. C         BASIS  GENERATED  IN   CONSTR ),  WHICH  WAS RETURNED AS THE
  637. C         VALUE OF  NFPOLS .
  638. C         NEPOLS  WILL BE CHANGED IF EVLDEG .GT. 0 TO GIVE THE SIZE OF
  639. C         BASIS REQUIRED FOR THE MULTINOMIAL OF DEGREE  EVLDEG .
  640. C      NEPTS  -- (INTEGER) -- (PASSED)
  641. C         THE NUMBER OF EVALUATION POINTS.
  642. C      EVLCDS  -- (INTEGER ) -- (PASSED)
  643. C          EVLCDS (P,K) IS THE VALUE OF THE K-TH VARIABLE AT THE  P-TH
  644. C         EVALUATION POINT.
  645. C      EVLVLS  -- (INTEGER) -- (RETURNED)
  646. C          EVLVLS (P) IS THE VALUE OF THE EVALUATED MULTINOMIAL AT THE
  647. C         P-TH EVALUATION POINT.
  648. C      ERROR  -- (INTEGER) -- (RETURNED)
  649. C          0 .........  IF NO ERRORS
  650. C         -1 .........  IF  NEPOLS  .GT.  NFPOLS  OR  NEPOLS  .LT. 1
  651. C         -2 .........  IF  NEPTS  .LT. 1 OR  DIMEN  .LT. 1
  652. C          NEPOLS  ...  IF  NEPOLS  .GT.  EDWKLN
  653. C      FITIWK  -- (INTEGER, 1-SUBSCRIPT ARRAY) -- (PASSED)
  654. C         THE INTEGER WORK ARRAY OF LENGTH  FIWKLN  THAT WAS  USED  IN
  655. C          CONSTR .
  656. C      FITDWK  -- (DOUBLE-PRECISION, 2-SUBSCRIPT ARRAY) -- (PASSED)
  657. C         THE DOUBLE PRECISION WORK ARRAY OF LENGTH  FDWKLN  THAT  WAS
  658. C         USED IN  CONSTR .
  659. C      FIWKLN  -- (INTEGER) -- (PASSED)
  660. C         THE LENGTH OF  FITIWK .
  661. C      FDWKLN  -- (INTEGER) -- (PASSED)
  662. C         THE LENGTH OF  FITDWK .
  663. C      EVLDWK  -- (DOUBLE-PRECISION, 1-SUBSCRIPT ARRAY ) -- (RETURNED)
  664. C         A WORK ARRAY OF LENGTH  NEPOLS  (OR LONGER).
  665. C      EDWKLN  -- (INTEGER) -- (PASSED)
  666. C         THE LENGTH OF  EVLDWK .
  667. C
  668. C     THE SUBROUTINE   EVALP   IS CALLED TO DO THE ACTUAL EVALUATION.
  669. C
  670. C     DATE LAST MODIFIED
  671. C     ---- ---- --------
  672. C     OCTOBER 16, 1984
  673. C     ****************
  674. C
  675. C
  676. C     ***************
  677. C     SET UP INDEX POINTERS TO THE BEGINNING OF EACH ROW OF
  678. C     THE TABLE -- THIS SETS THE BEGINNING POINT FOR EACH
  679. C     FULL MULTINOMIAL DEGREE.
  680. C     ***************
  681. C
  682.       ERROR = 0
  683.       GBASIZ = FITIWK(2)
  684.       IF ( EVLDEG .LT. 0 ) GO TO 20
  685.          TOP = 1
  686.          BOT = 1
  687.          DO 10 CURDEG = 1,EVLDEG
  688.             TOP = TOP * (DIMEN + CURDEG)
  689.  10         BOT = BOT * CURDEG
  690.          NEPOLS = TOP / BOT
  691.          IF ( EVLDEG .EQ. 0 ) NEPOLS = 1
  692.  20   IF ( NEPOLS .LE. GBASIZ .AND. NEPOLS.GE.1 ) GO TO 30
  693.          ERROR = -1
  694.          RETURN
  695.  30   IF ( NEPTS .GE. 1 .AND. DIMEN .GE. 1 ) GO TO 40
  696.          ERROR = -2
  697.          RETURN
  698.  40   IF ( NEPOLS .LE. EDWKLN ) GO TO 50
  699.          ERROR = NEPOLS
  700.          RETURN
  701.  50   DIMP1 = DIMEN + 1
  702.       ALFL = FITIWK(4)
  703.       MAXSTT = 1
  704.       ALFSTT = DIMP1 + MAXSTT
  705.       CSTT = ALFSTT + ALFL
  706. C
  707. C     ***************
  708. C     THE ACTUAL EVALUATION IS DONE INSIDE  EVALP .
  709. C     ***************
  710. C
  711.       CALL EVALP(EVLCDS,FITDWK(CSTT),NEPTS,DIMEN,NEPOLS,FITDWK(ALFSTT),
  712.      +           FITIWK,EVLDWK,EVLVLS,ALFL,FITDWK(MAXSTT),DIMP1)
  713.       RETURN
  714.       END
  715. C
  716. C ====================================================================
  717. C     UTILITIES -- UTILITIES -- UTILITIES -- UTILITIES -- UTILITIES
  718. C ====================================================================
  719. C
  720.       SUBROUTINE ALLOT(DEGREE,NPOLYS,NPTS,DIMEN,IWORK,IWKLEN,
  721.      +                 IREQD,DREQD,ERROR)
  722. C
  723.       INTEGER IREQD,DREQD,ALFL,ERROR,NPOLYS,DEGREE,DIMEN,NPTS
  724.       INTEGER NEWSTT,PSIWID,KMXBAS,STARTJ,KJP1D2,INDEX,IWKLEN
  725.       INTEGER NPLYT4
  726.       INTEGER IWORK(IWKLEN)
  727. C
  728. C     ***************
  729. C     PURPOSE
  730. C     -------
  731. C
  732. C      ALLOT  CHECKS FOR SUFFICIENCY THE DECLARED  DIMENSIONS  OF  THE
  733. C     WORK  ARRAYS  USED BY THE SUBROUTINE  CONSTR .  VARIOUS SIZES OF
  734. C     SUB-ARRAYS ARE COMPUTED AND REPORTED.
  735. C
  736. C     THIS ROUTINE IS CALLED BY  CONSTR .  IT IS NOT  CALLED  DIRECTLY
  737. C     BY THE USER.
  738. C
  739. C     THIS ROUTINE CALLS   BASIZ   AND   TABLE   FOR  THE  SUBSTANTIVE
  740. C     COMPUTATIONS.
  741. C
  742. C     VARIABLES
  743. C     ---------
  744. C
  745. C      DEGREE  - (PASSED/RETURNED)
  746. C         IGNORED IF .LT. 0.
  747. C         IF  DEGREE  .GE. 0 THEN  DEGREE IS CHECKED  AGAINST   NPTS .
  748. C         THE VALUE OF  DEGREE  WILL BE REDUCED IF THERE IS A BASIS OF
  749. C         MULTINOMIALS, ALL OF  DEGREE .LE.  DEGREE ,  OF  CARDINALITY
  750. C          NPTS
  751. C      NPOLYS  - (PASSED/RETURNED)
  752. C         IGNORED IF  DEGREE  .GE. 0.
  753. C         IF DEGREE .LT. 0 THEN THE VALUE OF NPOLYS WILL BE TAKEN  AS
  754. C         THE SIZE OF THE BASIS OF MULTINOMIALS TO BE USED IN THE FIT.
  755. C         NPOLYS  MUST SATISFY NPOLYS .LT.  NPTS  AND  NPOLYS  .GE. 1
  756. C      NPTS  --- (PASSED)
  757. C         THE NUMBER OF DATA POINTS TO BE USED IN THE FIT.
  758. C          NPTS  MUST BE .GE. 1.
  759. C      DIMEN  -- (PASSED)
  760. C         THE NUMBER OF VARIABLES.
  761. C      IWORK  -- (RETURNED)
  762. C         AN INTEGER WORK ARRAY OF LENGTH AT LEAST
  763. C            IF  DEGREE  .GE. 0 THEN
  764. C               4*BINOMIAL( DIMEN + DEGREE , DIMEN )
  765. C                  +( DIMEN )*( DEGREE )
  766. C            ELSE
  767. C               4*BINOMIAL( DIMEN +D,D)+( DIMEN )*D
  768. C         WHERE D IS THE MINIMUM CARDINALITY  OF  A  BASIS  OF  DEGREE
  769. C          DEGREE  SUCH THAT
  770. C            BINOMIAL( DIMEN +ABS( DEGREE ), DIMEN ) .GE.  NPOLYS
  771. C      IWKLEN  - (PASSED)
  772. C         THE LENGTH OF  IWORK
  773. C      IREQD  -- (RETURNED)
  774. C         THE SIZE OF THE INTEGER WORK ARRAY REQUIRED BY  CONSTR   FOR
  775. C         THE FIT SPECIFIED BY THE 4 INPUT PARAMETERS.
  776. C      DREQD  -- (RETURNED)
  777. C         THE SIZE OF THE DOUBLE  PRECISION  WORK  ARRAY  REQUIRED  BY
  778. C          CONSTR  FOR THE FIT SPECIFIED BY THE 4 INPUT PARAMETERS.
  779. C      ERROR  -- (RETURNED)
  780. C       0 IF  NPOLYS ,  DIMEN ,   DEGREE ,   NPTS   AND   IWKLEN   ARE
  781. C         VALID AND CONSISTENT WITH EACH OTHER.
  782. C       1 IF  DEGREE  .GE. 0 BUT THERE IS AN INTERPOLATING MULTINOMIAL
  783. C         OF SMALLER DEGREE OR IF  DEGREE .LT. 0 AND  NPOLYS .GT. NPTS
  784. C       2 IF  DEGREE  .LT. 0 AND  NPOLYS  .LE. 0
  785. C       3 IF  NPTS  .LT. 1 AND/OR  DIMEN  .LT. 1
  786. C       4 IF  IWKLEN  IS TOO SMALL (SET  IWKLEN  TO THE VALUE RETURNED
  787. C         IN  IREQD  TO RESOLVE THIS PROBLEM)
  788. C
  789. C     NOTE THAT  DEGREE ,  NPOLYS ,  PSIWID  AND  ALFL   ARE  RETURNED
  790. C     IN  IWORK (1-4), RESPECTIVELY.
  791. C
  792. C     DATE LAST MODIFIED
  793. C     ---- ---- --------
  794. C     DECEMBER 10, 1984
  795. C     ****************
  796. C
  797. C     ***************
  798. C     BASIZ  COMPUTES THE SIZE OF THE BASIS (AND AUXILIARY SIZES)
  799. C     BASED PRIMARILY UPON THE DEGREE, NUMBER OF  FITTING POINTS,
  800. C     AND THE DIMENSION.
  801. C     ***************
  802. C
  803.       CALL BASIZ(DEGREE,NPTS,DIMEN,NPOLYS,ERROR)
  804.       IF ( ERROR .GE. 2 ) RETURN
  805.       IREQD = 4 * NPOLYS + DEGREE * DIMEN
  806.       IF ( IWKLEN .GE. IREQD ) GO TO 5
  807.          ERROR = 4
  808.          RETURN
  809.  5    NEWSTT = 4 * NPOLYS + 1
  810. C
  811. C     ***************
  812. C     SET UP USEFUL INDEXING ARRAYS
  813. C          IWORK(1) ,..., IWORK(NEWSTT-1)
  814. C     AND
  815. C          IWORK(NEWSTT ,..., IWORK(NEWSTT+DIMEN*DEGREE)
  816. C     ***************
  817. C
  818.       CALL TABLE(DEGREE,DIMEN,NPOLYS,IWORK,IWORK(NEWSTT),ALFL)
  819.       IWORK(1) = DEGREE
  820.       IWORK(2) = NPOLYS
  821. C
  822. C     ***************
  823. C     FORCE  ALFL  TO BE AT LEAST 1 SO THAT DIMENSION STATEMENTS
  824. C     USING  ALFL  DO NOT BOMB.
  825. C     ***************
  826. C
  827.       IF ( ALFL .GT. 1 ) ALFL = ALFL - 1
  828.       IWORK(4) = ALFL
  829. C
  830. C     ***************
  831. C     THE FOLLOWING IS A SECTION OF CODE FOR SETTING UP THE
  832. C     STORAGE  MANAGEMENT OF THE  PSI  ARRAY.  THERE  IS  A
  833. C     COMPLICATED DOVETAILING FORMULA USED TO PACK INFORMATION
  834. C     INTO  PSI  WITHOUT LEAVING GAPS.
  835. C
  836. C     ARRAY          LENGTH
  837. C     -----          ------
  838. C      MAXABS         DIMEN  + 1
  839. C      ALPHA          ALFL
  840. C      C              NPOLYS
  841. C      SUMSQS         NPOLYS
  842. C
  843. C     THE NUMBER OF COLUMNS IN  PSI ,  PSIWID , IS DETERMINED BY
  844. C     PSIWID  =  NPOLYS  + 1 - (THE SMALLEST M SUCH THAT ALPHA(J,M)
  845. C                                   IS NONZERO AND J .GE. NPOLYS)
  846. C     THIS INSURES THAT IF THE USER EXTENDS THE BASIS, ALL THE  PSI
  847. C     REQUIRED WILL CERTAINLY BE STORED
  848. C
  849. C     IF DEGREE( NPOLYS ) .LE. 2 THEN                      (CASE 1)
  850. C       PSIWID  =  NPOLYS
  851. C     ELSE
  852. C       IF K =  DIMEN  THEN                                (CASE 2)
  853. C         PSIWID  =  NPOLYS
  854. C                    -  NEWKJ( 1 , DEGREE(NPOLYS)-1 )  + 1
  855. C       ELSE
  856. C            PSIWID  =  NPOLYS
  857. C                      + 1
  858. C                      - (
  859. C                         THE SMALLER OF
  860. C                           NEWKJ(K+1,DEGREE(NPOLYS)-2)    (CASE 3)
  861. C                         AND
  862. C                           INDEXS(3,NPOLYS)               (CASE 4)
  863. C                        )
  864. C
  865.       IF ( DEGREE .GT. 2 ) GO TO 10
  866. C
  867. C     ***************
  868. C     CASE 1
  869. C     ***************
  870. C
  871.          PSIWID = NPOLYS
  872.       GO TO 40
  873.  10      NPLYT4 = 4 * NPOLYS
  874. C
  875. C     ***************
  876. C     KMXBAS  IS K
  877. C                 NPOLYS
  878. C     ***************
  879. C
  880.          KMXBAS = IWORK(NPLYT4 - 2)
  881. C
  882.          IF ( KMXBAS .NE. DIMEN ) GO TO 20
  883. C
  884. C     ***************
  885. C     CASE 2
  886. C     ***************
  887. C
  888.             PSIWID = NPOLYS - IWORK(4 * NPOLYS - 1)
  889.          GO TO 40
  890. C
  891. C     ***************
  892. C     INDEX =  NEWKJ( K      + 1 , DEGREE(NPOLYS-2) )
  893. C                      NPOLYS
  894. C     ***************
  895. C
  896.  20         INDEX = NPLYT4 + (DEGREE - 3) * DIMEN + KMXBAS + 1
  897.             KJP1D2 = IWORK(INDEX)
  898. C
  899. C     ***************
  900. C     STARTJ  =  INDEXS(3,NPOLYS)
  901. C     ***************
  902. C
  903.             STARTJ = IWORK(NPLYT4 - 1)
  904.             IF ( STARTJ .GT. KJP1D2 ) GO TO 30
  905. C
  906. C     ***************
  907. C     CASE 4
  908. C     ***************
  909. C
  910.                PSIWID = NPOLYS - STARTJ + 1
  911.             GO TO 40
  912. C
  913. C     ***************
  914. C     CASE 3
  915. C     ***************
  916. C
  917.  30            PSIWID = NPOLYS - KJP1D2 + 1
  918.  40   IWORK(3) = PSIWID
  919.       DREQD = 2 * NPOLYS + DIMEN + 1 + NPTS * PSIWID + ALFL
  920.       RETURN
  921.       END
  922.       SUBROUTINE BASIZ(DEGREE,NPTS,DIMEN,NPOLYS,ERROR)
  923. C
  924.       INTEGER TOP,BOT,DEGREE,NPTS,DIMEN,NPOLYS,ERROR,I,ROWLEN
  925. C
  926. C     ***************
  927. C     PURPOSE
  928. C     -------
  929. C
  930. C     IF  DEGREE  .GE. 0 THEN
  931. C       FIND THE SIZE OF A BASIS REQUIRED EITHER TO
  932. C       1) APPROXIMATE THE DATA WITH A POLYNOMIAL OF DEGREE
  933. C          GIVEN BY THE PARAMETER  DEGREE
  934. C       OR TO
  935. C       2) SPAN THE SPACE OF POLYNOMIALS OF DEGREE .LE. THE
  936. C          SMALLEST DEGREE OF POLYNOMIAL WHICH INTERPOLATES THE
  937. C          DATA.
  938. C       IN CASE 1  ERROR  = 0.
  939. C       IN CASE 2  ERROR  = 1.
  940. C     ELSE
  941. C       IF  NPOLYS  .GE. 1 THEN
  942. C          IF NPOLYS .GT. NPTS THEN
  943. C             SET  NPOLYS  =  NPTS , FIND THE SMALLEST DEGREE OF A
  944. C             POLYNOMIAL  WHICH  INTERPOLATES  THE  DATA,  AND SET
  945. C              ERROR  = 1.
  946. C          ELSE
  947. C             FIND THE LARGEST DEGREE  DEGREE  OF A POLYNOMIAL  IN
  948. C             A  BASIS OF  NPOLYS  POLYNOMIALS GENERATED ACCORDING
  949. C             TO OUR ORDERING AND SET  ERROR  = 0.
  950. C       ELSE
  951. C           ERROR  = 2
  952. C
  953. C     THIS SUBROUTINE IS CALLED BY  ALLOT .  IT IS NOT  CALLED  BY
  954. C     THE USER DIRECTLY.
  955. C
  956. C     DATE LAST MODIFIED
  957. C     ---- ---- --------
  958. C     OCTOBER 16, 1984
  959. C     ****************
  960. C
  961.       ERROR = 0
  962.       IF ( NPTS .GE. 1 .AND. DIMEN .GE. 1 ) GO TO 10
  963.          ERROR = 3
  964.          RETURN
  965. C
  966.  10   CONTINUE
  967.       IF ( DEGREE .LT. 0 ) GO TO 30
  968. C
  969.         ROWLEN = 1
  970.         NPOLYS = 1
  971.         TOP = DIMEN - 1
  972.         BOT = 0
  973.         IF ( DEGREE .LT. 1 )  GO TO 30
  974.           DO 20 I=1,DEGREE
  975.             TOP = TOP + 1
  976.             BOT = BOT + 1
  977.             ROWLEN = (ROWLEN*TOP)/BOT
  978.             NPOLYS = NPOLYS + ROWLEN
  979.    20     CONTINUE
  980. C
  981.    30 CONTINUE
  982.       IF ( NPOLYS .GE. 1 ) GO TO 40
  983.             ERROR = 2
  984.             RETURN
  985.    40 CONTINUE
  986.       IF ( NPOLYS .LT. NPTS ) GO TO 50
  987.             NPOLYS = NPTS
  988.             ERROR = 1
  989.    50 CONTINUE
  990.       ROWLEN = 1
  991.       I = 1
  992.       DEGREE = 0
  993.       TOP = DIMEN - 1
  994.       BOT = 0
  995.    60 CONTINUE
  996.       IF ( I .GE. NPOLYS )  GO TO 70
  997.           TOP = TOP + 1
  998.           BOT = BOT + 1
  999.           ROWLEN = (ROWLEN*TOP)/BOT
  1000.           I = I + ROWLEN
  1001.           DEGREE = DEGREE + 1
  1002.           IF ( I .LT. NPOLYS )  GO TO 60
  1003.    70 CONTINUE
  1004.       RETURN
  1005.       END
  1006.       SUBROUTINE EVALP(COORD,C,NEPTS,DIMEN,NPOLYS,ALPHA,INDEXS,
  1007.      +                 PSI,F,ALFL,MAXABS,DIMP1)
  1008. C
  1009.       INTEGER DIMEN,NEPTS,NPOLYS,ALFL,DIMP1
  1010.       INTEGER JM1,JPRIME,M,P,K,I,J,INDEX
  1011.       INTEGER INDEXS(4,NPOLYS)
  1012.       DOUBLE PRECISION ALPHA(ALFL),COORD(NEPTS,DIMEN),PSI(NPOLYS)
  1013.       DOUBLE PRECISION C(NPOLYS),F(NEPTS),MAXABS(DIMP1)
  1014.       DOUBLE PRECISION RUNTOT,RNTOT1
  1015. C
  1016. C     ***************
  1017. C     PURPOSE
  1018. C     -------
  1019. C
  1020. C     THIS SUBROUTINE PERFORMS THE MAIN WORK OF EVALUATING THE
  1021. C     FITTING MULTINOMIAL (OR THE INITIAL PORTION OF IT  WHICH
  1022. C     IS REQUESTED BY THE SETTING OF  NEPOLS ,  EVLDEG  IN THE
  1023. C     CALL TO SUBROUTINE  EVAL .
  1024. C
  1025. C     THIS SUBROUTINE IS CALLED BY  EVAL .  IT IS  NOT  CALLED
  1026. C     DIRECTLY BY THE USER.
  1027. C
  1028. C     SUBROUTINES  SCALDN  AND  SCALUP  ARE   CALLED  TO SCALE
  1029. C     AND TO UNSCALE VALUES.
  1030. C
  1031. C     THE BODY OF THIS SUBROUTINE FOLLOWS THE EXPLANATION
  1032. C     GIVEN IN
  1033. C          LEAST SQUARES FITTING USING
  1034. C          ORTHOGONAL MULTINOMIALS
  1035. C     BY
  1036. C          BARTELS AND JEZIORANSKI
  1037. C     IN
  1038. C          ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE
  1039. C
  1040. C     DATE LAST MODIFIED
  1041. C     ---- ---- --------
  1042. C     OCTOBER 16, 1984
  1043. C     ****************
  1044. C
  1045. C     ***************
  1046. C     SCALE DOWN THE INDEPENDENT CO-ORDINATES.
  1047. C     ***************
  1048. C
  1049.       DO 10 K = 1,DIMEN
  1050.  10      CALL SCALDN(COORD(1,K),NEPTS,MAXABS(K))
  1051. C
  1052. C     ***************
  1053. C     USE THE  BASIS FUNCTION COEFFICIENTS  C  AND RECURRENCE
  1054. C     COEFFICIENTS  ALPHA  TO EVALUATE THE FITTED MULTINOMIAL
  1055. C     AT THE EVALUATION CO-ORDINATES  COORD .
  1056. C     ***************
  1057. C
  1058.       PSI(1) = 1.0D+00
  1059.       DO 40 P = 1,NEPTS
  1060.          RNTOT1 = C(1)
  1061.          IF ( NPOLYS .EQ. 1 ) GO TO 40
  1062.           DO 30 J = 2,NPOLYS
  1063.               K = INDEXS(2,J)
  1064.               JPRIME = INDEXS(1,J)
  1065.               RUNTOT = COORD(P,K) * PSI(JPRIME)
  1066.               I = INDEXS(3,J)
  1067.               JM1 = J - 1
  1068.               DO 20 M = I,JM1
  1069.                  INDEX = INDEXS(4,J) + M - I
  1070.  20              RUNTOT = RUNTOT - PSI(M) * ALPHA(INDEX)
  1071.               PSI(J) = RUNTOT
  1072.  30           RNTOT1 = RNTOT1 + C(J) * PSI(J)
  1073.  40       F(P) = RNTOT1
  1074. C
  1075. C     ***************
  1076. C     SCALE UP THE DEPENDENT COORDINATES.
  1077. C     ***************
  1078. C
  1079.       CALL SCALUP(F,NEPTS,MAXABS(DIMP1))
  1080.       DO 50 K = 1,DIMEN
  1081.  50      CALL SCALUP(COORD(1,K),NEPTS,MAXABS(K))
  1082. C
  1083.       RETURN
  1084.       END
  1085.       SUBROUTINE GNRTP(DEGREE,ALPHA,PSI,INDEXS,
  1086.      +                 NEWKJ,SUMSQS,COORD,NPOLYS,
  1087.      +                 DIMEN,NPTS,F,Z,C,PSIWID,WEIGHT,
  1088.      +                 ALFL,DIMP1,MAXABS)
  1089. C
  1090.       INTEGER DEGREE,DIMEN,NPOLYS,NPTS,K,PSIWID,ALFL,P,STTDEG,ONPLYS
  1091.       INTEGER ERROR,DIMP1
  1092.       INTEGER INDEXS(4,NPOLYS),NEWKJ(DIMEN,DEGREE)
  1093.       DOUBLE PRECISION PSI(NPTS,PSIWID),ALPHA(ALFL),F(NPTS)
  1094.       DOUBLE PRECISION COORD(NPTS,DIMEN),MAXABS(DIMP1),WEIGHT(NPTS)
  1095.       DOUBLE PRECISION Z(NPTS),SUMSQS(NPOLYS),C(NPOLYS)
  1096.       DOUBLE PRECISION RUNTOT,RNTOT1
  1097. C
  1098. C     ***************
  1099. C     PURPOSE
  1100. C     -------
  1101. C
  1102. C     THE MULTINOMIAL FIT IS GENERATED INCREMENTALLY, A BASIS  ELEMENT
  1103. C     AT A TIME.  THIS SUBROUTINE STARTS THE PROCESS OFF BY SETTING UP
  1104. C     THE FIRST BASIS ELEMENT, SCALING THE  DATA,  FINDING  THE  FIRST
  1105. C     COEFFICIENT,  AND  INITIALIZING  THE  WORK ARRAY Z.  GNRTP  THEN
  1106. C     CALLS  INCDG  IF MORE THAN ONE BASIS ELEMENT IS REQUIRED.
  1107. C
  1108. C     THIS SUBROUTINE IS CALLED BY  CONSTR .  IT IS NOT CALLED BY  THE
  1109. C     USER.
  1110. C
  1111. C     THIS SUBROUTINE CALLS  SCALPM ,  SCALDN , AND  INCDG .
  1112. C
  1113. C     DATE LAST MODIFIED
  1114. C     ---- ---- --------
  1115. C     OCTOBER 16, 1984
  1116. C     ****************
  1117. C
  1118. C     ***************
  1119. C     SET UP THE SCALING.
  1120. C     ***************
  1121. C
  1122.       DO 10 K = 1,DIMEN
  1123.          CALL SCALPM(COORD(1,K),NPTS,MAXABS(K))
  1124.  10      CALL SCALDN(COORD(1,K),NPTS,MAXABS(K))
  1125.       CALL SCALPM(F,NPTS,MAXABS(DIMP1))
  1126.       CALL SCALDN(F,NPTS,MAXABS(DIMP1))
  1127. C
  1128. C     ***************
  1129. C      SUMSQS (1) = (1,1)
  1130. C     C  = (F,1) / (1,1)
  1131. C      1
  1132. C     ***************
  1133. C
  1134.       RUNTOT = 0.0D+00
  1135.       RNTOT1 = 0.0D+00
  1136.       DO 20 P = 1,NPTS
  1137.          PSI(P,1) = 1.0D+00
  1138.          RNTOT1 = RNTOT1 + WEIGHT(P)
  1139.  20      RUNTOT = RUNTOT + F(P) * WEIGHT(P)
  1140.       SUMSQS(1) = RNTOT1
  1141.       C(1) = RUNTOT / RNTOT1
  1142. C
  1143. C     ***************
  1144. C     Z = F - C
  1145. C              1
  1146. C     ***************
  1147. C
  1148.       DO 30 P = 1,NPTS
  1149.  30      Z(P) = F(P) - C(1)
  1150. C
  1151.       IF ( NPOLYS .EQ. 1 ) RETURN
  1152.       STTDEG = 1
  1153.       ONPLYS = 1
  1154. C
  1155.       CALL INCDG(DEGREE,ALPHA,PSI,INDEXS,NEWKJ,SUMSQS,
  1156.      +           COORD,NPOLYS,DIMEN,NPTS,F,Z,C,PSIWID,
  1157.      +           WEIGHT,ALFL,ONPLYS,STTDEG,ERROR)
  1158.       RETURN
  1159.       END
  1160.       SUBROUTINE INCDG(DEGREE,ALPHA,PSI,INDEXS,NEWKJ,
  1161.      +                 SUMSQS,COORD,NPOLYS,
  1162.      +                 DIMEN,NPTS,F,Z,C,PSIWID,WEIGHT,
  1163.      +                 ALFL,ONPLYS,STTDEG,ERROR)
  1164. C
  1165.       INTEGER JPRIME,P,J,CURDEG,KJ,KJP,L,JPM1,JM1
  1166.       INTEGER M,START,JINDEX,JPINDX,Q,J3,J1,J1MJ2,ERROR
  1167.       INTEGER J0MJ1,J1M1,STARTA,ONPLYS,ONPP1,STTDEG,INDEX1,INDEX2
  1168.       INTEGER DEGREE,NPOLYS,NPTS,DIMEN,PSIWID,ALFL
  1169.       DOUBLE PRECISION ALPHA(ALFL),COORD(NPTS,DIMEN),PSI(NPTS,PSIWID)
  1170.       DOUBLE PRECISION SUMSQS(NPOLYS),C(NPOLYS),F(NPTS),WEIGHT(NPTS)
  1171.       DOUBLE PRECISION Z(NPTS)
  1172.       INTEGER INDEXS(4,NPOLYS),NEWKJ(DIMEN,DEGREE)
  1173.       DOUBLE PRECISION RUNTOT,RNTOT1,RNTOT2
  1174. C
  1175. C     ***************
  1176. C     PURPOSE
  1177. C     -------
  1178. C
  1179. C     THE MULTINOMIAL FIT IS GENERATED INCREMENTALLY, A BASIS  ELEMENT
  1180. C     AT A TIME.  THIS SUBROUTINE CONTINUES THE PROCESS STARTED OFF BY
  1181. C      GNRTP .
  1182. C
  1183. C     THIS SUBROUTINE IS CALLED BY  GNRTP  AND NOT BY THE USER.
  1184. C
  1185. C
  1186. C     DATE LAST MODIFIED
  1187. C     ---- ---- --------
  1188. C     OCTOBER 16, 1984
  1189. C     ****************
  1190. C
  1191.       ERROR = 0
  1192.       IF ( ONPLYS .GE. 1 .AND. STTDEG .GE. 1 ) GO TO 10
  1193.          ERROR = 6
  1194.       RETURN
  1195.  10      IF ( INDEXS(2,ONPLYS) .EQ. DIMEN ) GO TO 20
  1196.             CURDEG = STTDEG
  1197.          GO TO 30
  1198.  20         CURDEG = STTDEG + 1
  1199.  30   ONPP1 = ONPLYS + 1
  1200.       DO 170 J = ONPP1,NPOLYS
  1201.          JPRIME = INDEXS(1,J)
  1202.          JINDEX = J - (J - 1) / PSIWID * PSIWID
  1203.          JPINDX = JPRIME - (JPRIME - 1) / PSIWID * PSIWID
  1204.          KJ = INDEXS(2,J)
  1205.          START = INDEXS(3,J)
  1206.          M = START
  1207.          STARTA = INDEXS(4,J) - START
  1208.          IF ( CURDEG .EQ. 1 ) GO TO 100
  1209.          KJP = INDEXS(2,JPRIME)
  1210.          J1 = NEWKJ(KJ,CURDEG - 1)
  1211. C
  1212. C     ***************
  1213. C     CALCULATE THOSE  ALPHA ( J , M ) THAT CAN BE CALCULATED FROM
  1214. C     PREVIOUSLY CALCULATED ALPHAS.
  1215. C     ***************
  1216. C
  1217.          IF ( KJ .LT. KJP ) GO TO 50
  1218. C
  1219. C     ***************
  1220. C     FIRST CALCULATE THOSE BETWEEN  JPP  AND THE END OF 2 ROWS BACK.
  1221. C     CALCULATE  ALPHA ( J , JPP )
  1222. C     ***************
  1223. C
  1224.          INDEX1 = INDEXS(4,J)
  1225.          ALPHA(INDEX1) = SUMSQS(JPRIME) / SUMSQS(START)
  1226. C
  1227.          M = START + 1
  1228.          J3 = NEWKJ(1,CURDEG - 1) - 1
  1229.          IF ( M .GT. J3 ) GO TO 50
  1230. C
  1231. C     ***************
  1232. C     CURDEG .GT. 2 IF CONTROL HAS PASSED THE BRANCHES IN THE 3-RD
  1233. C     PREVIOUS AND 8-TH PREVIOUS STATEMENTS.
  1234. C     ***************
  1235. C
  1236.             J1MJ2 = J1 - NEWKJ(KJ,CURDEG - 2)
  1237. C
  1238.             DO 40 L = M,J3
  1239.                Q = J1MJ2 + L
  1240.                INDEX1 = STARTA + L
  1241.                INDEX2 = INDEXS(4,Q) - INDEXS(3,Q) + JPRIME
  1242.  40            ALPHA(INDEX1) = ALPHA(INDEX2) * SUMSQS(JPRIME) /
  1243.      +         SUMSQS(L)
  1244. C
  1245. C     ***************
  1246. C     CALCULATE  ALPHA ( J , M ) FOR  M  BETWEEN THE 2
  1247. C     RANGES CALCULATED BEFORE USING
  1248. C
  1249. C        ALPHA ( J , L ) = (X  *  PSI  ,PSI )  / (PSI ,PSI )
  1250. C                             K      JP    L         L    L
  1251. C                              J
  1252. C     ***************
  1253. C
  1254.          M = J3 + 1
  1255.  50      IF ( JPRIME .EQ. J1 ) GO TO 100
  1256.             IF ( KJ .EQ. 1 ) GO TO 80
  1257.             J1M1 = J1 - 1
  1258.             DO 70 L = M,J1M1
  1259.                RUNTOT = 0.0D+00
  1260.                DO 60 P = 1,NPTS
  1261.                   INDEX1 = L - (L - 1) / PSIWID * PSIWID
  1262.  60               RUNTOT = RUNTOT + COORD(P,KJ) * PSI(P,JPINDX) *
  1263.      +            PSI(P,INDEX1) * WEIGHT(P)
  1264.                INDEX1 = STARTA + L
  1265.  70            ALPHA(INDEX1) = RUNTOT / SUMSQS(L)
  1266. C
  1267. C     ***************
  1268. C     CALCULATE  ALPHA ( J , M ) FOR  M  BETWEEN
  1269. C      NEWKJ ( KJ , CURDEG  - 1)  AND
  1270. C      JP  - 1.
  1271. C     ***************
  1272. C
  1273.  80         J0MJ1 = NEWKJ(KJ,CURDEG) - J1
  1274.             JPM1 = JPRIME - 1
  1275.             DO 90 L = J1,JPM1
  1276.                Q = J0MJ1 + L
  1277.                INDEX1 = STARTA + L
  1278.                INDEX2 = INDEXS(4,Q) - INDEXS(3,Q) + JPRIME
  1279.  90            ALPHA(INDEX1) = ALPHA(INDEX2) * SUMSQS(JPRIME) /
  1280.      +         SUMSQS(L)
  1281.             M = JPRIME
  1282. C
  1283. C     ***************
  1284. C     CALCULATE THE REMAINING  ALPHA ( J , M ) FROM
  1285. C
  1286. C      ALPHA ( J , L ) = (X   * PSI  ,PSI ) / (PSI ,PSI )
  1287. C                          K       JP    L        L    L
  1288. C                           J
  1289. C     ***************
  1290. C
  1291.  100     JM1 = J - 1
  1292.          DO 120 L = M,JM1
  1293.             RUNTOT = 0.0D+00
  1294.             DO 110 P = 1,NPTS
  1295.                INDEX1 = L - (L - 1) / PSIWID * PSIWID
  1296.  110           RUNTOT = RUNTOT + COORD(P,KJ) * PSI(P,JPINDX) *
  1297.      +         PSI(P,INDEX1) * WEIGHT(P)
  1298.             INDEX1 = STARTA + L
  1299.  120        ALPHA(INDEX1) = RUNTOT / SUMSQS(L)
  1300. C
  1301. C     ***************
  1302. C     NOW CALCULATE THE  PSI (P,J),  SUMSQS (J) AND  C (J) USING
  1303. C
  1304. C                          J-1
  1305. C     PSI  = X  * PSI   -  SUM   ALPHA(J,L) * PSI
  1306. C        J    K      JP    L=JPP                 L
  1307. C
  1308. C     SUMSQS  = (PSI ,PSI )
  1309. C           J       J    J
  1310. C
  1311. C     C  = (Z,PSI )
  1312. C      J         J
  1313. C     ***************
  1314. C
  1315.  130      RNTOT1 = 0.0D+00
  1316.           RNTOT2 = 0.0D+00
  1317.           DO 150 P = 1,NPTS
  1318.               RUNTOT = COORD(P,KJ) * PSI(P,JPINDX)
  1319.               JM1 = J - 1
  1320.               DO 140 L = START,JM1
  1321.                  INDEX1 = STARTA + L
  1322.                  INDEX2 = L - (L - 1) / PSIWID * PSIWID
  1323.  140             RUNTOT = RUNTOT - ALPHA(INDEX1) * PSI(P,INDEX2)
  1324.               PSI(P,JINDEX) = RUNTOT
  1325.               RNTOT1 = RNTOT1 + PSI(P,JINDEX) * PSI(P,JINDEX) *
  1326.      +                          WEIGHT(P)
  1327.  150          RNTOT2 = RNTOT2 + Z(P) * PSI(P,JINDEX) * WEIGHT(P)
  1328.           SUMSQS(J) = RNTOT1
  1329.           C(J) = RNTOT2 / RNTOT1
  1330. C
  1331. C     ***************
  1332. C     CALCULATE THE NEW  Z ( P ) AND THE NEW  SSRES  USING
  1333. C
  1334. C     Z = Z - C  * PSI
  1335. C              J      J
  1336. C     ***************
  1337. C
  1338.           DO 160 P = 1,NPTS
  1339.  160          Z(P) = Z(P) - C(J) * PSI(P,JINDEX)
  1340.  170      IF ( KJ .EQ. DIMEN ) CURDEG = CURDEG + 1
  1341.       RETURN
  1342.       END
  1343.       SUBROUTINE MOVE(OLDARR,NEWARR,START,OLDWID,NEWIDT,COLENG,ERROR)
  1344. C
  1345.       INTEGER START,OLDWID,NEWIDT,COLENG,BIG,BIG1,LILN,BIGN,I,J
  1346.       INTEGER ERROR,JO,JN,OLDWPS,BIG1PS,BIGP1,K
  1347.       DOUBLE PRECISION OLDARR(COLENG,OLDWID),NEWARR(COLENG,NEWIDT)
  1348. C
  1349. C     ***************
  1350. C     PURPOSE
  1351. C     -------
  1352. C
  1353. C     MOVE COLUMNS 1 THROUGH   OLDWID   OF   OLDARR   INTO  COLUMNS  1
  1354. C     THROUGH  NEWWID  OF  NEWARR  USING
  1355. C         ( START  + I) MOD  OLDWID
  1356. C     THROUGH
  1357. C         ( START  + I) MOD  NEWID
  1358. C     FOR
  1359. C         I = 0
  1360. C     THROUGH
  1361. C         I =  OLDWID  - 1
  1362. C     THE MOVEMENT STARTS FROM  THE  LARGEST  COLUM  OF   NEWARR   AND
  1363. C     PROCEEDS DOWNWARD TO THE SMALLEST (SO THAT  OLDARR  AND  NEWARR
  1364. C     CAN BE THE SAME ARRAY).
  1365. C
  1366. C
  1367. C     DATE LAST MODIFIED
  1368. C     ---- ---- --------
  1369. C     OCTOBER 16, 1984
  1370. C     ****************
  1371. C
  1372.       ERROR = 1
  1373.       IF ( OLDWID .LT. 1 .OR. NEWIDT .LT. 1 .OR. NEWIDT .LT. OLDWID)
  1374.      +         RETURN
  1375.       ERROR = 0
  1376.       BIG = START + OLDWID - 1
  1377.       BIGN = BIG - (BIG - 1) / NEWIDT * NEWIDT
  1378.       LILN = START - (START - 1) / NEWIDT * NEWIDT
  1379.       IF ( LILN .GT. BIGN ) GO TO 20
  1380.          OLDWPS = OLDWID + START
  1381.          DO 10 I = 1,OLDWID
  1382.             J = OLDWPS - I
  1383.             JO = J - (J - 1) / OLDWID * OLDWID
  1384.             JN = J - (J - 1) / NEWIDT * NEWIDT
  1385.             DO 10 K = 1,COLENG
  1386.  10            NEWARR(K,JN) = OLDARR(K,JO)
  1387.          RETURN
  1388.  20   BIG1 = NEWIDT - LILN + 1
  1389.       BIG1PS = BIG1 + START
  1390.       DO 30 I = 1,BIG1
  1391.          J = BIG1PS - I
  1392.          JO = J - (J - 1) / OLDWID * OLDWID
  1393.          JN = J - (J - 1) / NEWIDT * NEWIDT
  1394.          DO 30 K = 1,COLENG
  1395.  30         NEWARR(K,JN) = OLDARR(K,JO)
  1396.       BIGP1 = BIG + 1
  1397.       DO 40 I = 1,BIGN
  1398.          J = BIGP1 - I
  1399.          JO = J - (J - 1) / OLDWID * OLDWID
  1400.          JN = J - (J - 1) / NEWIDT * NEWIDT
  1401.          DO 40 K = 1,COLENG
  1402.  40         NEWARR(K,JN) = NEWARR(K,JO)
  1403.       RETURN
  1404.       END
  1405.       SUBROUTINE RESTRT(PSIWID,DWORK,DREQD,ONPLYS,OPSWID,OLALFL,CSTT,
  1406.      +                  SSQSTT,PSISTT,NPTS,DIMEN)
  1407. C
  1408.       INTEGER DREQD,ONPLYS,OPSWID,OLALFL,CSTT
  1409.       INTEGER OSSQST,OPSIST,PSIWID,NPTS,ERROR,DIMEN
  1410.       INTEGER I,J,SSQSTT,PSISTT,OCST,START,INDEX1,INDEX2
  1411.       DOUBLE PRECISION DWORK(DREQD)
  1412. C
  1413. C     ***************
  1414. C     PURPOSE
  1415. C     -------
  1416. C
  1417. C     THIS SUBROUTINE REARRANGES THE WORK  SPACE  (WITH THE HELP
  1418. C     OF SUBROUTINE  MOVE ) IN THE EVENT THAT A FIT OF INCREASED
  1419. C     DEGREE IS TO BE MADE.
  1420. C
  1421. C     CALLED INTERNALLY BY  CONSTR , NOT BY THE USER.
  1422. C
  1423. C     DATE LAST MODIFIED
  1424. C     ---- ---- --------
  1425. C     OCTOBER 16, 1984
  1426. C     ****************
  1427. C
  1428.       OCST = 2 + DIMEN + OLALFL
  1429.       OSSQST = OCST + ONPLYS
  1430.       OPSIST = OSSQST + ONPLYS
  1431.       START = ONPLYS - OPSWID
  1432.       CALL MOVE(DWORK(OPSIST),DWORK(PSISTT),START,OPSWID,PSIWID,NPTS,
  1433.      +          ERROR)
  1434. C
  1435.       DO 5 J = 1,ONPLYS
  1436.          I = ONPLYS - J
  1437.          INDEX1 = SSQSTT + I
  1438.          INDEX2 = OSSQST + I
  1439.  5       DWORK(INDEX1) = DWORK(INDEX2)
  1440.       DO 10 J = 1,ONPLYS
  1441.          I = ONPLYS - J
  1442.          INDEX1 = CSTT + I
  1443.          INDEX2 = OCST + I
  1444.  10      DWORK(INDEX1) = DWORK(INDEX2)
  1445.       RETURN
  1446.       END
  1447.       SUBROUTINE SCALPM(COORD,NPTS,MAXABS)
  1448. C
  1449.       INTEGER NPTS,P
  1450.       DOUBLE PRECISION MAXABS,A
  1451.       DOUBLE PRECISION COORD(NPTS)
  1452. C
  1453. C     ***************
  1454. C     PURPOSE
  1455. C     -------
  1456. C
  1457. C     FIND SCALING PARAMETER(S) FOR THE PROBLEM. IF THE SCALING SCHEME
  1458. C     IS CHANGED, ALL FOUR OF THE FOLLOWING WOULD HAVE TO BE CHANGED
  1459. C
  1460. C     1)  SCALPM  - FIND THE SCALING PARAMETERS
  1461. C     2)  SCALDN  - SCALE THE PROBLEM DATA
  1462. C     3)  SCALUP  - PERFORM THE INVERSE TRANSFORMATION TO  SCALDN
  1463. C     4) THE SCALING OF THE RESIDUALS IN  CONSTR
  1464. C
  1465. C     THIS SUBROUTINE IS CALLED BY  GNRTP .  IT IS NOT CALLED  BY  THE
  1466. C     USER.
  1467. C
  1468. C     THE SCALING WHICH  IT  DEFINES  MUST  BE  COORDINATED  WITH  THE
  1469. C     SCALING  OF RESIDUALS WHICH IS CARRIED OUT TOWARD THE END OF THE
  1470. C     SUBROUTINE  CONSTR .  THE SCALING DEFINED  BY  THIS  ROUTINE  IS
  1471. C     APPLIED  IN  THE  SECTION OF  CONSTR  JUST MENTIONED, AND BY THE
  1472. C     ROUTINES  SCALUP  AND  SCALDN .
  1473. C
  1474. C     DATE LAST MODIFIED
  1475. C     ---- ---- --------
  1476. C     OCTOBER 16, 1984
  1477. C     ****************
  1478. C
  1479.       MAXABS = 0.0D+00
  1480.       DO 5 P = 1,NPTS
  1481.          A = DABS(COORD(P))
  1482.  5       IF ( A .GT. MAXABS ) MAXABS = A
  1483.       RETURN
  1484.       END
  1485.       SUBROUTINE SCALDN(COORD,NPTS,MAXABS)
  1486. C
  1487.       INTEGER NPTS,P
  1488.       DOUBLE PRECISION MAXABS
  1489.       DOUBLE PRECISION COORD(NPTS)
  1490. C
  1491. C     ***************
  1492. C     PURPOSE
  1493. C     -------
  1494. C
  1495. C     CARRY OUT THE DATA-SCALING WHICH IS DEFINED  BY  THE  SUBROUTINE
  1496. C      SCALPM .
  1497. C
  1498. C     THIS SUBROUTINE IS CALLED BY  GNRTP  AND   EVALP .   IT  IS  NOT
  1499. C     CALLED BY THE USER.
  1500. C
  1501. C     THE SCALING WHICH THIS ROUTINE CARRIES OUT  MUST  BE  CONSISTENT
  1502. C     WITH  THE  SCALING  OF RESIDUALS WHICH IS CARRIED OUT TOWARD THE
  1503. C     END OF THE SUBROUTINE  CONSTR .
  1504. C
  1505. C     DATE LAST MODIFIED
  1506. C     ---- ---- --------
  1507. C     OCTOBER 16, 1984
  1508. C     ****************
  1509. C
  1510.       IF ( MAXABS .EQ. 0.0D+00 ) RETURN
  1511.          DO 5 P = 1,NPTS
  1512.  5          COORD(P) = COORD(P) / MAXABS
  1513.          RETURN
  1514.       END
  1515.       SUBROUTINE SCALUP(COORD,NPTS,MAXABS)
  1516. C
  1517.       INTEGER NPTS,P
  1518.       DOUBLE PRECISION MAXABS
  1519.       DOUBLE PRECISION COORD(NPTS)
  1520. C
  1521. C     ***************
  1522. C     PURPOSE
  1523. C     -------
  1524. C
  1525. C     REMOVE THE DATA-SCALING  WHICH  IS  DEFINED  BY  THE  SUBROUTINE
  1526. C      SCALPM  AND APPLIED BY THE SUBROUTINE  SCALDN .
  1527. C
  1528. C     THIS SUBROUTINE IS CALLED BY  EVALP .   IT IS  NOT CALLED BY THE
  1529. C     USER.
  1530. C
  1531. C     THE UNSCALING WHICH THIS ROUTINE CARRIES OUT MUST BE THE INVERSE
  1532. C     OF  THE SCALING OF RESIDUALS WHICH IS CARRIED OUT TOWARD THE END
  1533. C     OF THE SUBROUTINE  CONSTR .
  1534. C
  1535. C     DATE LAST MODIFIED
  1536. C     ---- ---- --------
  1537. C     OCTOBER 16, 1984
  1538. C     ****************
  1539. C
  1540.       IF ( MAXABS .EQ. 0.0D+00 ) RETURN
  1541.          DO 5 P = 1,NPTS
  1542.  5          COORD(P) = COORD(P) * MAXABS
  1543.          RETURN
  1544.       END
  1545.       SUBROUTINE TABLE(DEGREE,DIMEN,NPOLYS,INDEXS,NEWKJ,ALFLP1)
  1546. C
  1547.       INTEGER J,KJ,CURDEG,JPRIME,NWITHK,I,CURM1,RALEN,DIMM1,DIMM2
  1548.       INTEGER NPOLYS,DIMEN,DEGREE,ALFLP1,DIMP1
  1549.       INTEGER INDEXS(4,NPOLYS),NEWKJ(DIMEN,DEGREE)
  1550. C
  1551. C     ***************
  1552. C     PURPOSE
  1553. C     -------
  1554. C
  1555. C     TABULATE  JP  AND  KJ  FOR EACH  J
  1556. C
  1557. C     VARIABLES
  1558. C     ---------
  1559. C
  1560. C      ALFLP1  -- (INTEGER) -- (PASSED)
  1561. C         THE LENGTH REQUIRED FOR ARRAY  ALPHA , PLUS ONE
  1562. C      DEGREE  -- (INTEGER) -- (PASSED)
  1563. C         THE DEGREE OF THE POLYNOMIAL TO BE FITTED
  1564. C      DIMEN  -- (INTEGER) -- (PASSED)
  1565. C         NUMBER OF INDEPENDENT VARIABLES
  1566. C      INDEXS  -- (INTEGER, 2-SUBSCRIPT ARRAY) -- (RETURNED)
  1567. C          INDEXS (1, J )   IS    JP ,    INDEXS (2, J )   IS     KJ ,
  1568. C          INDEXS (3, J )  IS THE FIRST NONZERO RECURRENCE COEFFICIENT
  1569. C         IN  ALPHA  AND  INDEXS (4, J ) IS ITS LOCATION IN  ALPHA .
  1570. C      NEWKJ  -- (INTEGER, 2-SUBSCRIPT ARRAY) -- (RETURNED)
  1571. C          NEWKJ ( K , D ) IS THE FIRST MONOMIAL OF DEGREE  D   HAVING
  1572. C          KJ = K .
  1573. C      NPOLYS  -- (INTEGER) -- (PASSED)
  1574. C         NUMBER  OF  MONOMIALS  OF  DEGREE  .LE.  ORDER  IN   DIMEN
  1575. C         INDEPENDENT VARIABLES.
  1576. C
  1577. C     THIS SUBPROGRAM CAN BE CODED (EXCLUDING THE PART FOR CALCULATING
  1578. C      INDEXS (3, J )  AND   INDEXS (4, J )) MENTALLY MORE EFFICIENTLY
  1579. C     BUT COMPUTATIONALLY LESS EFFICIENTLY AS
  1580. C
  1581. C         J = 2
  1582. C         DO 5 KJ = 1,DIMEN
  1583. C           NEWKJ(KJ,1) = KJ + 1
  1584. C           INDEXS(1,J) = 1
  1585. C           INDEXS(2,J) = KJ
  1586. C           J = J + 1
  1587. C      5  CONTINUE
  1588. C         DO 10 CURDEG = 2,DEGREE
  1589. C           DO 10 KJ = 1,DIMEN
  1590. C             JPRIME = NEWKJ(KJ,CURDEG - 1)
  1591. C             NEWKJ(KJ,CURDEG) = J
  1592. C             NWITHK = COMB(DIMEN + CURDEG - KJ - 1,CURDEG - 1)
  1593. C             DO 10 I = 1,NWITHK
  1594. C               INDEXS(1,J) = JPRIME
  1595. C               INDEXS(2,J) = KJ
  1596. C               JPRIME = JPRIME + 1
  1597. C               J = J + 1
  1598. C     10  CONTINUE
  1599. C
  1600. C     WHERE COMB(N,KJ) IS N-FACTORIAL / ((N-KJ)-FACTORIAL * KJ-FACTORIAL
  1601. C     HERE WE MAKE USE OF THE RECURRENCE RELATIONS
  1602. C
  1603. C       COMB(DIMEN+CURDEG-2,CURDEG-1)
  1604. C
  1605. C                        (DIMEN+CURDEG-2)
  1606. C          = ----------------------------------------
  1607. C            (CURDEG-1)*COMB(DIMEN+CURDEG-3,CURDEG-2)
  1608. C
  1609. C     AND
  1610. C
  1611. C       COMB(DIMEN+CURDEG-KJ-1,CURDEG-1)
  1612. C
  1613. C                              (DIMEN-KJ+1)
  1614. C          = ----------------------------------------------
  1615. C            (DIMEN+CURDEG-KJ)*COMB(DIMEN+CURDEG-KJ,CURDEG-1)
  1616. C
  1617. C
  1618. C     DATE LAST MODIFIED
  1619. C     ---- ---- --------
  1620. C     OCTOBER 16, 1984
  1621. C     ****************
  1622. C
  1623.       ALFLP1 = 1
  1624. C
  1625. C     ***************
  1626. C     SET  INDEXS (4,1) TO 1 SO THAT  ALFL - INDEXS (4,1)+1 IS THE
  1627. C     NUMBER OF COLUMNS REQUIRED FOR  PSI  FOR  NPOLYS =1  ( ALFL
  1628. C     IS DEFINED IN THE MAINLINE TO BE  ALFLP1 -1 IF ALFLP1 .GT. 1
  1629. C     AND  ALFLP1  OTHERWISE.
  1630. C     ***************
  1631. C
  1632.       INDEXS(4,1) = 1
  1633. C
  1634.       IF ( NPOLYS .EQ. 1 ) RETURN
  1635.       J = 2
  1636.       DO 10 KJ = 1,DIMEN
  1637.          NEWKJ(KJ,1) = KJ + 1
  1638.          INDEXS(1,J) = 1
  1639.          INDEXS(2,J) = KJ
  1640.          INDEXS(3,J) = 1
  1641.          INDEXS(4,J) = ALFLP1
  1642.          ALFLP1 = ALFLP1 + J - 1
  1643.          IF ( J.EQ.NPOLYS ) RETURN
  1644. 10       J = J + 1
  1645.       IF ( DEGREE .EQ. 1 ) RETURN
  1646.       RALEN = 1
  1647.       DIMM1 = DIMEN - 1
  1648.       DIMM2 = DIMEN - 2
  1649.       DIMP1 = DIMEN + 1
  1650.       DO 70 CURDEG = 2,DEGREE
  1651.          CURM1 = CURDEG - 1
  1652.          RALEN = RALEN * (DIMM2 + CURDEG) / CURM1
  1653.          NWITHK = RALEN
  1654.          KJ = 1
  1655. 20       JPRIME = NEWKJ(KJ,CURM1)
  1656.          NEWKJ(KJ,CURDEG) = J
  1657.          IF ( KJ.EQ.DIMEN ) GO TO 60
  1658.             DO 50 I = 1,NWITHK
  1659.                INDEXS(1,J) = JPRIME
  1660.                INDEXS(2,J) = KJ
  1661. C
  1662. C     ***************
  1663. C     CALCULATE  INDEXS (3, J ),  INDEXS (4, J )
  1664. C     ***************
  1665. C
  1666.                IF ( KJ .LT. INDEXS(2,JPRIME) ) GO TO 30
  1667.                   INDEXS(3,J) = INDEXS(1,JPRIME)
  1668.                   GO TO 40
  1669.  30            INDEXS(3,J) = NEWKJ(1,CURDEG - 1)
  1670.  40            INDEXS(4,J) = ALFLP1
  1671.                ALFLP1 = ALFLP1 + J - INDEXS(3,J)
  1672.                IF ( J .EQ. NPOLYS ) RETURN
  1673. C
  1674.                JPRIME = JPRIME + 1
  1675.  50            J = J + 1
  1676.             KJ = KJ + 1
  1677.             NWITHK = NWITHK * (DIMP1 - KJ) / (DIMEN + CURDEG - KJ)
  1678.             GO TO 20
  1679.  60      INDEXS(1,J) = JPRIME
  1680.          INDEXS(2,J) = DIMEN
  1681.          INDEXS(3,J) = INDEXS(1,JPRIME)
  1682.          INDEXS(4,J) = ALFLP1
  1683.          ALFLP1 = ALFLP1 + J - INDEXS(3,J)
  1684.          IF ( J .EQ. NPOLYS ) RETURN
  1685.  70      J = J + 1
  1686.       RETURN
  1687.       END
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!