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