*SODR SUBROUTINE SODR + (FCN, + N,M,NP,NQ, + BETA, + Y,LDY,X,LDX, + WE,LDWE,LD2WE,WD,LDWD,LD2WD, + JOB, + IPRINT,LUNERR,LUNRPT, + WORK,LWORK,IWORK,LIWORK, + INFO) C***BEGIN PROLOGUE SODR C***DATE WRITTEN 860529 (YYMMDD) C***REVISION DATE 920619 (YYMMDD) C***CATEGORY NO. G2E,I1B1 C***KEYWORDS ORTHOGONAL DISTANCE REGRESSION, C NONLINEAR LEAST SQUARES, C MEASUREMENT ERROR MODELS, C ERRORS IN VARIABLES C***AUTHOR BOGGS, PAUL T. C APPLIED AND COMPUTATIONAL MATHEMATICS DIVISION C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C BYRD, RICHARD H. C DEPARTMENT OF COMPUTER SCIENCE C UNIVERSITY OF COLORADO, BOULDER, CO 80309 C ROGERS, JANET E. C APPLIED AND COMPUTATIONAL MATHEMATICS DIVISION C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C BOULDER, CO 80303-3328 C SCHNABEL, ROBERT B. C DEPARTMENT OF COMPUTER SCIENCE C UNIVERSITY OF COLORADO, BOULDER, CO 80309 C AND C APPLIED AND COMPUTATIONAL MATHEMATICS DIVISION C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C BOULDER, CO 80303-3328 C***PURPOSE SINGLE PRECISION DRIVER ROUTINE FOR FINDING C THE WEIGHTED EXPLICIT OR IMPLICIT ORTHOGONAL DISTANCE C REGRESSION (ODR) OR ORDINARY LINEAR OR NONLINEAR LEAST C SQUARES (OLS) SOLUTION (SHORT CALL STATEMENT) C***DESCRIPTION C FOR DETAILS, SEE ODRPACK USER'S REFERENCE GUIDE. C***REFERENCES BOGGS, P. T., R. H. BYRD, J. R. DONALDSON, AND C R. B. SCHNABEL (1989), C "ALGORITHM 676 --- ODRPACK: SOFTWARE FOR WEIGHTED C ORTHOGONAL DISTANCE REGRESSION," C ACM TRANS. MATH. SOFTWARE., 15(4):348-364. C BOGGS, P. T., R. H. BYRD, J. E. ROGERS, AND C R. B. SCHNABEL (1992), C "USER'S REFERENCE GUIDE FOR ODRPACK VERSION 2.01, C SOFTWARE FOR WEIGHTED ORTHOGONAL DISTANCE REGRESSION," C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C INTERNAL REPORT NUMBER 92-4834. C BOGGS, P. T., R. H. BYRD, AND R. B. SCHNABEL (1987), C "A STABLE AND EFFICIENT ALGORITHM FOR NONLINEAR C ORTHOGONAL DISTANCE REGRESSION," C SIAM J. SCI. STAT. COMPUT., 8(6):1052-1078. C***ROUTINES CALLED SODCNT C***END PROLOGUE SODR C...SCALAR ARGUMENTS INTEGER + INFO,JOB,LDWD,LDWE,LDX,LDY,LD2WD,LD2WE,LIWORK,LWORK, + M,N,NDIGIT,NP,NQ C...ARRAY ARGUMENTS REAL + BETA(NP),WD(LDWD,LD2WD,M),WE(LDWE,LD2WE,NQ),WORK(LWORK), + X(LDX,M),Y(LDY,NQ) INTEGER + IWORK(LIWORK) C...SUBROUTINE ARGUMENTS EXTERNAL + FCN C...LOCAL SCALARS REAL + NEGONE,PARTOL,SSTOL,TAUFAC,ZERO INTEGER + IPRINT,LDIFX,LDSCLD,LDSTPD,LUNERR,LUNRPT,MAXIT LOGICAL + SHORT C...LOCAL ARRAYS REAL + SCLB(1),SCLD(1,1),STPB(1),STPD(1,1),WD1(1,1,1) INTEGER + IFIXB(1),IFIXX(1,1) C...EXTERNAL SUBROUTINES EXTERNAL + SODCNT C...DATA STATEMENTS DATA + NEGONE,ZERO + /-1.0E0,0.0E0/ C...ROUTINE NAMES USED AS SUBPROGRAM ARGUMENTS C FCN: THE USER-SUPPLIED SUBROUTINE FOR EVALUATING THE MODEL. C...VARIABLE DEFINITIONS (ALPHABETICALLY) C BETA: THE FUNCTION PARAMETERS. C IFIXB: THE VALUES DESIGNATING WHETHER THE ELEMENTS OF BETA ARE C FIXED AT THEIR INPUT VALUES OR NOT. C IFIXX: THE VALUES DESIGNATING WHETHER THE ELEMENTS OF X ARE C FIXED AT THEIR INPUT VALUES OR NOT. C INFO: THE VARIABLE DESIGNATING WHY THE COMPUTATIONS WERE STOPPED. C IPRINT: THE PRINT CONTROL VARIABLE. C IWORK: THE INTEGER WORK SPACE. C JOB: THE VARIABLE CONTROLLING PROBLEM INITIALIZATION AND C COMPUTATIONAL METHOD. C LDIFX: THE LEADING DIMENSION OF ARRAY IFIXX. C LDSCLD: THE LEADING DIMENSION OF ARRAY SCLD. C LDSTPD: THE LEADING DIMENSION OF ARRAY STPD. C LDWD: THE LEADING DIMENSION OF ARRAY WD. C LDWE: THE LEADING DIMENSION OF ARRAY WE. C LDX: THE LEADING DIMENSION OF ARRAY X. C LDY: THE LEADING DIMENSION OF ARRAY Y. C LD2WD: THE SECOND DIMENSION OF ARRAY WD. C LD2WE: THE SECOND DIMENSION OF ARRAY WE. C LIWORK: THE LENGTH OF VECTOR IWORK. C LUNERR: THE LOGICAL UNIT NUMBER FOR ERROR MESSAGES. C LUNRPT: THE LOGICAL UNIT NUMBER FOR COMPUTATION REPORTS. C LWORK: THE LENGTH OF VECTOR WORK. C M: THE NUMBER OF COLUMNS OF DATA IN THE EXPLANATORY VARIABLE. C MAXIT: THE MAXIMUM NUMBER OF ITERATIONS ALLOWED. C N: THE NUMBER OF OBSERVATIONS. C NEGONE: THE VALUE -1.0E0. C NDIGIT: THE NUMBER OF ACCURATE DIGITS IN THE FUNCTION RESULTS, AS C SUPPLIED BY THE USER. C NP: THE NUMBER OF FUNCTION PARAMETERS. C NQ: THE NUMBER OF RESPONSES PER OBSERVATION. C PARTOL: THE PARAMETER CONVERGENCE STOPPING TOLERANCE. C SCLB: THE SCALING VALUES FOR BETA. C SCLD: THE SCALING VALUES FOR DELTA. C STPB: THE RELATIVE STEP FOR COMPUTING FINITE DIFFERENCE C DERIVATIVES WITH RESPECT TO BETA. C STPD: THE RELATIVE STEP FOR COMPUTING FINITE DIFFERENCE C DERIVATIVES WITH RESPECT TO DELTA. C SHORT: THE VARIABLE DESIGNATING WHETHER THE USER HAS INVOKED C ODRPACK BY THE SHORT-CALL (SHORT=.TRUE.) OR THE LONG-CALL C (SHORT=.FALSE.). C SSTOL: THE SUM-OF-SQUARES CONVERGENCE STOPPING TOLERANCE. C TAUFAC: THE FACTOR USED TO COMPUTE THE INITIAL TRUST REGION C DIAMETER. C WD: THE DELTA WEIGHTS. C WD1: A DUMMY ARRAY USED WHEN WD(1,1,1)=0.0E0. C WE: THE EPSILON WEIGHTS. C WORK: THE REAL WORK SPACE. C X: THE EXPLANATORY VARIABLE. C Y: THE DEPENDENT VARIABLE. UNUSED WHEN THE MODEL IS IMPLICIT. C***FIRST EXECUTABLE STATEMENT SODR C INITIALIZE NECESSARY VARIABLES TO INDICATE USE OF DEFAULT VALUES IFIXB(1) = -1 IFIXX(1,1) = -1 LDIFX = 1 NDIGIT = -1 TAUFAC = NEGONE SSTOL = NEGONE PARTOL = NEGONE MAXIT = -1 STPB(1) = NEGONE STPD(1,1) = NEGONE LDSTPD = 1 SCLB(1) = NEGONE SCLD(1,1) = NEGONE LDSCLD = 1 SHORT = .TRUE. IF (WD(1,1,1).NE.ZERO) THEN CALL SODCNT + (SHORT, FCN, N,M,NP,NQ, BETA, Y,LDY,X,LDX, + WE,LDWE,LD2WE,WD,LDWD,LD2WD, IFIXB,IFIXX,LDIFX, + JOB,NDIGIT,TAUFAC, SSTOL,PARTOL,MAXIT, + IPRINT,LUNERR,LUNRPT, + STPB,STPD,LDSTPD, SCLB,SCLD,LDSCLD, + WORK,LWORK,IWORK,LIWORK, + INFO) ELSE WD1(1,1,1) = NEGONE CALL SODCNT + (SHORT, FCN, N,M,NP,NQ, BETA, Y,LDY,X,LDX, + WE,LDWE,LD2WE,WD1,1,1, IFIXB,IFIXX,LDIFX, + JOB,NDIGIT,TAUFAC, SSTOL,PARTOL,MAXIT, + IPRINT,LUNERR,LUNRPT, + STPB,STPD,LDSTPD, SCLB,SCLD,LDSCLD, + WORK,LWORK,IWORK,LIWORK, + INFO) END IF RETURN END *SODRC SUBROUTINE SODRC + (FCN, + N,M,NP,NQ, + BETA, + Y,LDY,X,LDX, + WE,LDWE,LD2WE,WD,LDWD,LD2WD, + IFIXB,IFIXX,LDIFX, + JOB,NDIGIT,TAUFAC, + SSTOL,PARTOL,MAXIT, + IPRINT,LUNERR,LUNRPT, + STPB,STPD,LDSTPD, + SCLB,SCLD,LDSCLD, + WORK,LWORK,IWORK,LIWORK, + INFO) C***BEGIN PROLOGUE SODRC C***DATE WRITTEN 860529 (YYMMDD) C***REVISION DATE 920619 (YYMMDD) C***CATEGORY NO. G2E,I1B1 C***KEYWORDS ORTHOGONAL DISTANCE REGRESSION, C NONLINEAR LEAST SQUARES, C MEASUREMENT ERROR MODELS, C ERRORS IN VARIABLES C***AUTHOR BOGGS, PAUL T. C APPLIED AND COMPUTATIONAL MATHEMATICS DIVISION C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C BYRD, RICHARD H. C DEPARTMENT OF COMPUTER SCIENCE C UNIVERSITY OF COLORADO, BOULDER, CO 80309 C ROGERS, JANET E. C APPLIED AND COMPUTATIONAL MATHEMATICS DIVISION C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C BOULDER, CO 80303-3328 C SCHNABEL, ROBERT B. C DEPARTMENT OF COMPUTER SCIENCE C UNIVERSITY OF COLORADO, BOULDER, CO 80309 C AND C APPLIED AND COMPUTATIONAL MATHEMATICS DIVISION C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C BOULDER, CO 80303-3328 C***PURPOSE SINGLE PRECISION DRIVER ROUTINE FOR FINDING C THE WEIGHTED EXPLICIT OR IMPLICIT ORTHOGONAL DISTANCE C REGRESSION (ODR) OR ORDINARY LINEAR OR NONLINEAR LEAST C SQUARES (OLS) SOLUTION (LONG CALL STATEMENT) C***DESCRIPTION C FOR DETAILS, SEE ODRPACK USER'S REFERENCE GUIDE. C***REFERENCES BOGGS, P. T., R. H. BYRD, J. R. DONALDSON, AND C R. B. SCHNABEL (1989), C "ALGORITHM 676 --- ODRPACK: SOFTWARE FOR WEIGHTED C ORTHOGONAL DISTANCE REGRESSION," C ACM TRANS. MATH. SOFTWARE., 15(4):348-364. C BOGGS, P. T., R. H. BYRD, J. E. ROGERS, AND C R. B. SCHNABEL (1992), C "USER'S REFERENCE GUIDE FOR ODRPACK VERSION 2.01, C SOFTWARE FOR WEIGHTED ORTHOGONAL DISTANCE REGRESSION," C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C INTERNAL REPORT NUMBER 92-4834. C BOGGS, P. T., R. H. BYRD, AND R. B. SCHNABEL (1987), C "A STABLE AND EFFICIENT ALGORITHM FOR NONLINEAR C ORTHOGONAL DISTANCE REGRESSION," C SIAM J. SCI. STAT. COMPUT., 8(6):1052-1078. C***ROUTINES CALLED SODCNT C***END PROLOGUE SODRC C...SCALAR ARGUMENTS REAL + PARTOL,SSTOL,TAUFAC INTEGER + INFO,IPRINT,JOB,LDIFX,LDSCLD,LDSTPD,LDWD,LDWE,LDX,LDY, + LD2WD,LD2WE,LIWORK,LUNERR,LUNRPT,LWORK,M,MAXIT,N,NDIGIT,NP,NQ C...ARRAY ARGUMENTS REAL + BETA(NP),SCLB(NP),SCLD(LDSCLD,M),STPB(NP),STPD(LDSTPD,M), + WD(LDWD,LD2WD,M),WE(LDWE,LD2WE,NQ),WORK(LWORK), + X(LDX,M),Y(LDY,NQ) INTEGER + IFIXB(NP),IFIXX(LDIFX,M),IWORK(LIWORK) C...SUBROUTINE ARGUMENTS EXTERNAL + FCN C...LOCAL SCALARS REAL + NEGONE,ZERO LOGICAL + SHORT C...LOCAL ARRAYS REAL + WD1(1,1,1) C...EXTERNAL SUBROUTINES EXTERNAL + SODCNT C...DATA STATEMENTS DATA + NEGONE,ZERO + /-1.0E0,0.0E0/ C...ROUTINE NAMES USED AS SUBPROGRAM ARGUMENTS C FCN: THE USER-SUPPLIED SUBROUTINE FOR EVALUATING THE MODEL. C...VARIABLE DEFINITIONS (ALPHABETICALLY) C BETA: THE FUNCTION PARAMETERS. C IFIXB: THE VALUES DESIGNATING WHETHER THE ELEMENTS OF BETA ARE C FIXED AT THEIR INPUT VALUES OR NOT. C IFIXX: THE VALUES DESIGNATING WHETHER THE ELEMENTS OF X ARE C FIXED AT THEIR INPUT VALUES OR NOT. C INFO: THE VARIABLE DESIGNATING WHY THE COMPUTATIONS WERE STOPPED. C IPRINT: THE PRINT CONTROL VARIABLE. C IWORK: THE INTEGER WORK SPACE. C JOB: THE VARIABLE CONTROLLING PROBLEM INITIALIZATION AND C COMPUTATIONAL METHOD. C LDIFX: THE LEADING DIMENSION OF ARRAY IFIXX. C LDSCLD: THE LEADING DIMENSION OF ARRAY SCLD. C LDSTPD: THE LEADING DIMENSION OF ARRAY STPD. C LDWD: THE LEADING DIMENSION OF ARRAY WD. C LDWE: THE LEADING DIMENSION OF ARRAY WE. C LDX: THE LEADING DIMENSION OF ARRAY X. C LDY: THE LEADING DIMENSION OF ARRAY Y. C LD2WD: THE SECOND DIMENSION OF ARRAY WD. C LD2WE: THE SECOND DIMENSION OF ARRAY WE. C LIWORK: THE LENGTH OF VECTOR IWORK. C LUNERR: THE LOGICAL UNIT NUMBER FOR ERROR MESSAGES. C LUNRPT: THE LOGICAL UNIT NUMBER FOR COMPUTATION REPORTS. C LWORK: THE LENGTH OF VECTOR WORK. C M: THE NUMBER OF COLUMNS OF DATA IN THE EXPLANATORY VARIABLE. C MAXIT: THE MAXIMUM NUMBER OF ITERATIONS ALLOWED. C N: THE NUMBER OF OBSERVATIONS. C NDIGIT: THE NUMBER OF ACCURATE DIGITS IN THE FUNCTION RESULTS, AS C SUPPLIED BY THE USER. C NP: THE NUMBER OF FUNCTION PARAMETERS. C NQ: THE NUMBER OF RESPONSES PER OBSERVATION. C PARTOL: THE PARAMETER CONVERGENCE STOPPING TOLERANCE. C SCLB: THE SCALING VALUES FOR BETA. C SCLD: THE SCALING VALUES FOR DELTA. C STPB: THE RELATIVE STEP FOR COMPUTING FINITE DIFFERENCE C DERIVATIVES WITH RESPECT TO BETA. C STPD: THE RELATIVE STEP FOR COMPUTING FINITE DIFFERENCE C DERIVATIVES WITH RESPECT TO DELTA. C SHORT: THE VARIABLE DESIGNATING WHETHER THE USER HAS INVOKED C ODRPACK BY THE SHORT-CALL (SHORT=.TRUE.) OR THE LONG-CALL C (SHORT=.FALSE.). C SSTOL: THE SUM-OF-SQUARES CONVERGENCE STOPPING TOLERANCE. C TAUFAC: THE FACTOR USED TO COMPUTE THE INITIAL TRUST REGION C DIAMETER. C WD: THE DELTA WEIGHTS. C WD1: A DUMMY ARRAY USED WHEN WD(1,1,1)=0.0E0. C WE: THE EPSILON WEIGHTS. C WORK: THE REAL WORK SPACE. C X: THE EXPLANATORY VARIABLE. C Y: THE DEPENDENT VARIABLE. UNUSED WHEN THE MODEL IS IMPLICIT. C***FIRST EXECUTABLE STATEMENT SODRC SHORT = .FALSE. IF (WD(1,1,1).NE.ZERO) THEN CALL SODCNT + (SHORT, FCN, N,M,NP,NQ, BETA, Y,LDY,X,LDX, + WE,LDWE,LD2WE,WD,LDWD,LD2WD, IFIXB,IFIXX,LDIFX, + JOB,NDIGIT,TAUFAC, SSTOL,PARTOL,MAXIT, + IPRINT,LUNERR,LUNRPT, + STPB,STPD,LDSTPD, SCLB,SCLD,LDSCLD, + WORK,LWORK,IWORK,LIWORK, + INFO) ELSE WD1(1,1,1) = NEGONE CALL SODCNT + (SHORT, FCN, N,M,NP,NQ, BETA, Y,LDY,X,LDX, + WE,LDWE,LD2WE,WD1,1,1, IFIXB,IFIXX,LDIFX, + JOB,NDIGIT,TAUFAC, SSTOL,PARTOL,MAXIT, + IPRINT,LUNERR,LUNRPT, + STPB,STPD,LDSTPD, SCLB,SCLD,LDSCLD, + WORK,LWORK,IWORK,LIWORK, + INFO) END IF RETURN END *SACCES SUBROUTINE SACCES + (N,M,NP,NQ,LDWE,LD2WE, + WORK,LWORK,IWORK,LIWORK, + ACCESS,ISODR, + JPVT,OMEGA,U,QRAUX,SD,VCV,WRK1,WRK2,WRK3,WRK4,WRK5,WRK6, + NNZW,NPP, + JOB,PARTOL,SSTOL,MAXIT,TAUFAC,ETA,NETA, + LUNRPT,IPR1,IPR2,IPR2F,IPR3, + WSS,RVAR,IDF, + TAU,ALPHA,NITER,NFEV,NJEV,INT2,OLMAVG, + RCOND,IRANK,ACTRS,PNORM,PRERS,RNORMS,ISTOP) C***BEGIN PROLOGUE SACCES C***REFER TO SODR,SODRC C***ROUTINES CALLED SIWINF,SWINF C***DATE WRITTEN 860529 (YYMMDD) C***REVISION DATE 920619 (YYMMDD) C***PURPOSE ACCESS OR STORE VALUES IN THE WORK ARRAYS C***END PROLOGUE SACESS C...SCALAR ARGUMENTS REAL + ACTRS,ALPHA,ETA,OLMAVG,PARTOL,PNORM,PRERS,RCOND, + RNORMS,RVAR,SSTOL,TAU,TAUFAC INTEGER + IDF,INT2,IPR1,IPR2,IPR2F,IPR3,IRANK,ISTOP,ISTOPI,JOB,JPVT, + LDWE,LD2WE,LIWORK,LUNRPT,LWORK,M,MAXIT,N,NETA,NFEV,NITER,NJEV, + NNZW,NP,NPP,NQ,OMEGA,QRAUX,SD,U,VCV, + WRK1,WRK2,WRK3,WRK4,WRK5,WRK6 LOGICAL + ACCESS,ISODR C...ARRAY ARGUMENTS REAL + WORK(LWORK),WSS(3) INTEGER + IWORK(LIWORK) C...LOCAL SCALARS INTEGER + ACTRSI,ALPHAI,BETACI,BETANI,BETASI,BETA0I, + DELTAI,DELTNI,DELTSI,DIFFI,EPSI, + EPSMAI,ETAI,FJACBI,FJACDI,FNI,FSI,IDFI,INT2I,IPRINI,IPRINT, + IRANKI,JOBI,JPVTI,LDTTI,LIWKMN,LUNERI,LUNRPI,LWKMN,MAXITI, + MSGB,MSGD,NETAI,NFEVI,NITERI,NJEVI,NNZWI,NPPI,NROWI, + NTOLI,OLMAVI,OMEGAI,PARTLI,PNORMI,PRERSI,QRAUXI,RCONDI, + RNORSI,RVARI,SDI,SI,SSFI,SSI,SSTOLI,TAUFCI,TAUI,TI,TTI,UI, + VCVI,WE1I,WRK1I,WRK2I,WRK3I,WRK4I,WRK5I,WRK6I,WRK7I, + WSSI,WSSDEI,WSSEPI,XPLUSI C...EXTERNAL SUBROUTINES EXTERNAL + SIWINF,SWINF C...VARIABLE DEFINITIONS (ALPHABETICALLY) C ACCESS: THE VARIABLE DESIGNATING WHETHER INFORMATION IS TO BE C ACCESSED FROM THE WORK ARRAYS (ACCESS=TRUE) OR STORED IN C THEM (ACCESS=FALSE). C ACTRS: THE SAVED ACTUAL RELATIVE REDUCTION IN THE SUM-OF-SQUARES. C ACTRSI: THE LOCATION IN ARRAY WORK OF VARIABLE ACTRS. C ALPHA: THE LEVENBERG-MARQUARDT PARAMETER. C ALPHAI: THE LOCATION IN ARRAY WORK OF VARIABLE ALPHA. C BETACI: THE STARTING LOCATION IN ARRAY WORK OF ARRAY BETAC. C BETANI: THE STARTING LOCATION IN ARRAY WORK OF ARRAY BETAN. C BETASI: THE STARTING LOCATION IN ARRAY WORK OF ARRAY BETAS. C BETA0I: THE STARTING LOCATION IN ARRAY WORK OF ARRAY BETA0. C DELTAI: THE STARTING LOCATION IN ARRAY WORK OF ARRAY DELTA. C DELTNI: THE STARTING LOCATION IN ARRAY WORK OF ARRAY DELTAN. C DELTSI: THE STARTING LOCATION IN ARRAY WORK OF ARRAY DELTAS. C DIFFI: THE STARTING LOCATION IN ARRAY WORK OF ARRAY DIFF. C EPSI: THE STARTING LOCATION IN ARRAY WORK OF ARRAY EPS. C EPSMAI: THE LOCATION IN ARRAY WORK OF VARIABLE EPSMAC. C ETA: THE RELATIVE NOISE IN THE FUNCTION RESULTS. C ETAI: THE LOCATION IN ARRAY WORK OF VARIABLE ETA. C FJACBI: THE STARTING LOCATION IN ARRAY WORK OF ARRAY FJACB. C FJACDI: THE STARTING LOCATION IN ARRAY WORK OF ARRAY FJACD. C FNI: THE STARTING LOCATION IN ARRAY WORK OF ARRAY FN. C FSI: THE STARTING LOCATION IN ARRAY WORK OF ARRAY FS. C IDF: THE DEGREES OF FREEDOM OF THE FIT, EQUAL TO THE NUMBER OF C OBSERVATIONS WITH NONZERO WEIGHTED DERIVATIVES MINUS THE C NUMBER OF PARAMETERS BEING ESTIMATED. C IDFI: THE STARTING LOCATION IN ARRAY IWORK OF VARIABLE IDF. C INT2: THE NUMBER OF INTERNAL DOUBLING STEPS. C INT2I: THE LOCATION IN ARRAY IWORK OF VARIABLE INT2. C IPR1: THE VALUE OF THE FOURTH DIGIT (FROM THE RIGHT) OF IPRINT, C WHICH CONTROLS THE INITIAL SUMMARY REPORT. C IPR2: THE VALUE OF THE THIRD DIGIT (FROM THE RIGHT) OF IPRINT, C WHICH CONTROLS THE ITERATION REPORTS. C IPR2F: THE VALUE OF THE SECOND DIGIT (FROM THE RIGHT) OF IPRINT, C WHICH CONTROLS THE FREQUENCY OF THE ITERATION REPORTS. C IPR3: THE VALUE OF THE FIRST DIGIT (FROM THE RIGHT) OF IPRINT, C WHICH CONTROLS THE FINAL SUMMARY REPORT. C IPRINI: THE LOCATION IN ARRAY IWORK OF VARIABLE IPRINT. C IPRINT: THE PRINT CONTROL VARIABLE. C IRANK: THE RANK DEFICIENCY OF THE JACOBIAN WRT BETA. C IRANKI: THE LOCATION IN ARRAY IWORK OF VARIABLE IRANK. C ISODR: THE VARIABLE DESIGNATING WHETHER THE SOLUTION IS TO BE C FOUND BY ODR (ISODR=TRUE) OR BY OLS (ISODR=FALSE). C ISTOP: THE VARIABLE DESIGNATING WHETHER THERE ARE PROBLEMS C COMPUTING THE FUNCTION AT THE CURRENT BETA AND DELTA. C ISTOPI: THE LOCATION IN ARRAY IWORK OF VARIABLE ISTOP. C IWORK: THE INTEGER WORK SPACE. C JOB: THE VARIABLE CONTROLING PROBLEM INITIALIZATION AND C COMPUTATIONAL METHOD. C JOBI: THE LOCATION IN ARRAY IWORK OF VARIABLE JOB. C JPVT: THE PIVOT VECTOR. C JPVTI: THE STARTING LOCATION IN ARRAY IWORK OF VARIABLE JPVT. C LDTTI: THE STARTING LOCATION IN ARRAY IWORK OF VARIABLE LDTT. C LDWE: THE LEADING DIMENSION OF ARRAY WE. C LD2WE: THE SECOND DIMENSION OF ARRAY WE. C LIWORK: THE LENGTH OF VECTOR IWORK. C LUNERI: THE LOCATION IN ARRAY IWORK OF VARIABLE LUNERR. C LUNERR: THE LOGICAL UNIT NUMBER USED FOR ERROR MESSAGES. C LUNRPI: THE LOCATION IN ARRAY IWORK OF VARIABLE LUNRPT. C LUNRPT: THE LOGICAL UNIT NUMBER USED FOR COMPUTATION REPORTS. C LWKMN: THE MINIMUM ACCEPTABLE LENGTH OF ARRAY WORK. C LWORK: THE LENGTH OF VECTOR WORK. C M: THE NUMBER OF COLUMNS OF DATA IN THE EXPLANATORY VARIABLE. C MAXIT: THE MAXIMUM NUMBER OF ITERATIONS ALLOWED. C MAXITI: THE LOCATION IN ARRAY IWORK OF VARIABLE MAXIT. C MSGB: THE STARTING LOCATION IN ARRAY IWORK OF ARRAY MSGB. C MSGD: THE STARTING LOCATION IN ARRAY IWORK OF ARRAY MSGD. C N: THE NUMBER OF OBSERVATIONS. C NETA: THE NUMBER OF ACCURATE DIGITS IN THE FUNCTION RESULTS. C NETAI: THE LOCATION IN ARRAY IWORK OF VARIABLE NETA. C NFEV: THE NUMBER OF FUNCTION EVALUATIONS. C NFEVI: THE LOCATION IN ARRAY IWORK OF VARIABLE NFEV. C NITER: THE NUMBER OF ITERATIONS TAKEN. C NITERI: THE LOCATION IN ARRAY IWORK OF VARIABLE NITER. C NJEV: THE NUMBER OF JACOBIAN EVALUATIONS. C NJEVI: THE LOCATION IN ARRAY IWORK OF VARIABLE NJEV. C NNZW: THE NUMBER OF NONZERO WEIGHTED OBSERVATIONS. C NNZWI: THE LOCATION IN ARRAY IWORK OF VARIABLE NNZW. C NP: THE NUMBER OF FUNCTION PARAMETERS. C NPP: THE NUMBER OF FUNCTION PARAMETERS ACTUALLY ESTIMATED. C NPPI: THE LOCATION IN ARRAY IWORK OF VARIABLE NPP. C NQ: THE NUMBER OF RESPONSES PER OBSERVATION. C NROWI: THE LOCATION IN ARRAY IWORK OF VARIABLE NROW. C NTOLI: THE LOCATION IN ARRAY IWORK OF VARIABLE NTOL. C OLMAVG: THE AVERAGE NUMBER OF LEVENBERG-MARQUARDT STEPS PER C ITERATION. C OLMAVI: THE LOCATION IN ARRAY WORK OF VARIABLE OLMAVG. C OMEGA: THE STARTING LOCATION IN ARRAY WORK OF ARRAY OMEGA. C OMEGAI: THE STARTING LOCATION IN ARRAY WORK OF ARRAY OMEGA. C PARTLI: THE LOCATION IN ARRAY WORK OF VARIABLE PARTOL. C PARTOL: THE PARAMETER CONVERGENCE STOPPING TOLERANCE. C PNORM: THE NORM OF THE SCALED ESTIMATED PARAMETERS. C PNORMI: THE LOCATION IN ARRAY WORK OF VARIABLE PNORM. C PRERS: THE SAVED PREDICTED RELATIVE REDUCTION IN THE C SUM-OF-SQUARES. C PRERSI: THE LOCATION IN ARRAY WORK OF VARIABLE PRERS. C QRAUX: THE STARTING LOCATION IN ARRAY WORK OF ARRAY QRAUX. C QRAUXI: THE STARTING LOCATION IN ARRAY WORK OF ARRAY QRAUX. C RCOND: THE APPROXIMATE RECIPROCAL CONDITION OF FJACB. C RCONDI: THE LOCATION IN ARRAY WORK OF VARIABLE RCOND. C RESTRT: THE VARIABLE DESIGNATING WHETHER THE CALL IS A RESTART C (RESTRT=TRUE) OR NOT (RESTRT=FALSE). C RNORMS: THE NORM OF THE SAVED WEIGHTED EPSILONS AND DELTAS. C RNORSI: THE LOCATION IN ARRAY WORK OF VARIABLE RNORMS. C RVAR: THE RESIDUAL VARIANCE, I.E. STANDARD DEVIATION SQUARED. C RVARI: THE LOCATION IN ARRAY WORK OF VARIABLE RVAR. C SCLB: THE SCALING VALUES USED FOR BETA. C SCLD: THE SCALING VALUES USED FOR DELTA. C SD: THE STARTING LOCATION IN ARRAY WORK OF ARRAY SD. C SDI: THE STARTING LOCATION IN ARRAY WORK OF ARRAY SD. C SHORT: THE VARIABLE DESIGNATING WHETHER THE USER HAS INVOKED C ODRPACK BY THE SHORT-CALL (SHORT=TRUE) OR THE LONG- C CALL (SHORT=FALSE). C SI: THE STARTING LOCATION IN ARRAY WORK OF ARRAY S. C SSFI: THE STARTING LOCATION IN ARRAY WORK OF ARRAY SSF. C SSI: THE STARTING LOCATION IN ARRAY WORK OF ARRAY SS. C SSTOL: THE SUM-OF-SQUARES CONVERGENCE STOPPING TOLERANCE. C SSTOLI: THE LOCATION IN ARRAY WORK OF VARIABLE SSTOL. C TAU: THE TRUST REGION DIAMETER. C TAUFAC: THE FACTOR USED TO COMPUTE THE INITIAL TRUST REGION C DIAMETER. C TAUFCI: THE LOCATION IN ARRAY WORK OF VARIABLE TAUFAC. C TAUI: THE LOCATION IN ARRAY WORK OF VARIABLE TAU. C TI: THE STARTING LOCATION IN ARRAY WORK OF ARRAY T. C TTI: THE STARTING LOCATION IN ARRAY WORK OF ARRAY TT. C U: THE STARTING LOCATION IN ARRAY WORK OF ARRAY U. C UI: THE STARTING LOCATION IN ARRAY WORK OF ARRAY U. C VCV: THE STARTING LOCATION IN ARRAY WORK OF ARRAY VCV. C VCVI: THE STARTING LOCATION IN ARRAY WORK OF ARRAY VCV. C WE1I: THE STARTING LOCATION IN ARRAY WORK OF ARRAY WE1. C WORK: THE REAL WORK SPACE. C WRK1: THE STARTING LOCATION IN ARRAY WORK OF ARRAY WRK1. C WRK1I: THE STARTING LOCATION IN ARRAY WORK OF ARRAY WRK1. C WRK2: THE STARTING LOCATION IN ARRAY WORK OF ARRAY WRK2. C WRK2I: THE STARTING LOCATION IN ARRAY WORK OF ARRAY WRK2. C WRK3: THE STARTING LOCATION IN ARRAY WORK OF ARRAY WRK3. C WRK3I: THE STARTING LOCATION IN ARRAY WORK OF ARRAY WRK3. C WRK4: THE STARTING LOCATION IN ARRAY WORK OF ARRAY WRK4. C WRK4I: THE STARTING LOCATION IN ARRAY WORK OF ARRAY WRK4. C WRK5: THE STARTING LOCATION IN ARRAY WORK OF ARRAY WRK5. C WRK5I: THE STARTING LOCATION IN ARRAY WORK OF ARRAY WRK5. C WRK6: THE STARTING LOCATION IN ARRAY WORK OF ARRAY WRK6. C WRK6I: THE STARTING LOCATION IN ARRAY WORK OF ARRAY WRK6. C WRK7I: THE STARTING LOCATION IN ARRAY WORK OF ARRAY WRK7. C WSS: THE SUM OF THE SQUARES OF THE WEIGHTED EPSILONS AND DELTAS, C THE SUM OF THE SQUARES OF THE WEIGHTED DELTAS, AND C THE SUM OF THE SQUARES OF THE WEIGHTED EPSILONS. C WSSI: THE STARTING LOCATION IN ARRAY WORK OF VARIABLE WSS(1). C WSSDEI: THE STARTING LOCATION IN ARRAY WORK OF VARIABLE WSS(2). C WSSEPI: THE STARTING LOCATION IN ARRAY WORK OF VARIABLE WSS(3). C XPLUSI: THE STARTING LOCATION IN ARRAY WORK OF ARRAY XPLUSD. C***FIRST EXECUTABLE STATEMENT SACCES C FIND STARTING LOCATIONS WITHIN INTEGER WORKSPACE CALL SIWINF(M,NP,NQ, + MSGB,MSGD,JPVTI,ISTOPI, + NNZWI,NPPI,IDFI, + JOBI,IPRINI,LUNERI,LUNRPI, + NROWI,NTOLI,NETAI, + MAXITI,NITERI,NFEVI,NJEVI,INT2I,IRANKI,LDTTI, + LIWKMN) C FIND STARTING LOCATIONS WITHIN REAL WORK SPACE CALL SWINF(N,M,NP,NQ,LDWE,LD2WE,ISODR, + DELTAI,EPSI,XPLUSI,FNI,SDI,VCVI, + RVARI,WSSI,WSSDEI,WSSEPI,RCONDI,ETAI, + OLMAVI,TAUI,ALPHAI,ACTRSI,PNORMI,RNORSI,PRERSI, + PARTLI,SSTOLI,TAUFCI,EPSMAI, + BETA0I,BETACI,BETASI,BETANI,SI,SSI,SSFI,QRAUXI,UI, + FSI,FJACBI,WE1I,DIFFI, + DELTSI,DELTNI,TI,TTI,OMEGAI,FJACDI, + WRK1I,WRK2I,WRK3I,WRK4I,WRK5I,WRK6I,WRK7I, + LWKMN) IF (ACCESS) THEN C SET STARTING LOCATIONS FOR WORK VECTORS JPVT = JPVTI OMEGA = OMEGAI QRAUX = QRAUXI SD = SDI VCV = VCVI U = UI WRK1 = WRK1I WRK2 = WRK2I WRK3 = WRK3I WRK4 = WRK4I WRK5 = WRK5I WRK6 = WRK6I C ACCESS VALUES FROM THE WORK VECTORS ACTRS = WORK(ACTRSI) ALPHA = WORK(ALPHAI) ETA = WORK(ETAI) OLMAVG = WORK(OLMAVI) PARTOL = WORK(PARTLI) PNORM = WORK(PNORMI) PRERS = WORK(PRERSI) RCOND = WORK(RCONDI) WSS(1) = WORK(WSSI) WSS(2) = WORK(WSSDEI) WSS(3) = WORK(WSSEPI) RVAR = WORK(RVARI) RNORMS = WORK(RNORSI) SSTOL = WORK(SSTOLI) TAU = WORK(TAUI) TAUFAC = WORK(TAUFCI) NETA = IWORK(NETAI) IRANK = IWORK(IRANKI) JOB = IWORK(JOBI) LUNRPT = IWORK(LUNRPI) MAXIT = IWORK(MAXITI) NFEV = IWORK(NFEVI) NITER = IWORK(NITERI) NJEV = IWORK(NJEVI) NNZW = IWORK(NNZWI) NPP = IWORK(NPPI) IDF = IWORK(IDFI) INT2 = IWORK(INT2I) C SET UP PRINT CONTROL VARIABLES IPRINT = IWORK(IPRINI) IPR1 = MOD(IPRINT,10000)/1000 IPR2 = MOD(IPRINT,1000)/100 IPR2F = MOD(IPRINT,100)/10 IPR3 = MOD(IPRINT,10) ELSE C STORE VALUES INTO THE WORK VECTORS WORK(ACTRSI) = ACTRS WORK(ALPHAI) = ALPHA WORK(OLMAVI) = OLMAVG WORK(PARTLI) = PARTOL WORK(PNORMI) = PNORM WORK(PRERSI) = PRERS WORK(RCONDI) = RCOND WORK(WSSI) = WSS(1) WORK(WSSDEI) = WSS(2) WORK(WSSEPI) = WSS(3) WORK(RVARI) = RVAR WORK(RNORSI) = RNORMS WORK(SSTOLI) = SSTOL WORK(TAUI) = TAU IWORK(IRANKI) = IRANK IWORK(ISTOPI) = ISTOP IWORK(NFEVI) = NFEV IWORK(NITERI) = NITER IWORK(NJEVI) = NJEV IWORK(IDFI) = IDF IWORK(INT2I) = INT2 END IF RETURN END *SESUBI SUBROUTINE SESUBI + (N,M,WD,LDWD,LD2WD,ALPHA,TT,LDTT,I,E) C***BEGIN PROLOGUE SESUBI C***REFER TO SODR,SODRC C***ROUTINES CALLED SZERO C***DATE WRITTEN 860529 (YYMMDD) C***REVISION DATE 920304 (YYMMDD) C***PURPOSE COMPUTE E = WD + ALPHA*TT**2 C***END PROLOGUE SESUBI C...SCALAR ARGUMENTS REAL + ALPHA INTEGER + LDTT,LDWD,LD2WD,M,N C...ARRAY ARGUMENTS REAL + E(M,M),TT(LDTT,M),WD(LDWD,LD2WD,M) C...LOCAL SCALARS REAL + ZERO INTEGER + I,J,J1,J2 C...EXTERNAL SUBROUTINES EXTERNAL + SZERO C...DATA STATEMENTS DATA + ZERO + /0.0E0/ C...VARIABLE DEFINITIONS (ALPHABETICALLY) C ALPHA: THE LEVENBERG-MARQUARDT PARAMETER. C E: THE VALUE OF THE ARRAY E = WD + ALPHA*TT**2 C I: AN INDEXING VARIABLE. C J: AN INDEXING VARIABLE. C J1: AN INDEXING VARIABLE. C J2: AN INDEXING VARIABLE. C LDWD: THE LEADING DIMENSION OF ARRAY WD. C LD2WD: THE SECOND DIMENSION OF ARRAY WD. C M: THE NUMBER OF COLUMNS OF DATA IN THE INDEPENDENT VARIABLE. C N: THE NUMBER OF OBSERVATIONS. C NP: THE NUMBER OF RESPONSES PER OBSERVATION. C TT: THE SCALING VALUES USED FOR DELTA. C WD: THE SQUARED DELTA WEIGHTS, D**2. C ZERO: THE VALUE 0.0E0. C***FIRST EXECUTABLE STATEMENT SESUBI C N.B. THE LOCATIONS OF WD AND TT ACCESSED DEPEND ON THE VALUE C OF THE FIRST ELEMENT OF EACH ARRAY AND THE LEADING DIMENSIONS C OF THE MULTIPLY SUBSCRIPTED ARRAYS. IF (N.EQ.0 .OR. M.EQ.0) RETURN IF (WD(1,1,1).GE.ZERO) THEN IF (LDWD.GE.N) THEN C THE ELEMENTS OF WD HAVE BEEN INDIVIDUALLY SPECIFIED IF (LD2WD.EQ.1) THEN C THE ARRAYS STORED IN WD ARE DIAGONAL CALL SZERO(M,M,E,M) DO 10 J=1,M E(J,J) = WD(I,1,J) 10 CONTINUE ELSE C THE ARRAYS STORED IN WD ARE FULL POSITIVE SEMIDEFINITE MATRICES DO 30 J1=1,M DO 20 J2=1,M E(J1,J2) = WD(I,J1,J2) 20 CONTINUE 30 CONTINUE END IF IF (TT(1,1).GT.ZERO) THEN IF (LDTT.GE.N) THEN DO 110 J=1,M E(J,J) = E(J,J) + ALPHA*TT(I,J)**2 110 CONTINUE ELSE DO 120 J=1,M E(J,J) = E(J,J) + ALPHA*TT(1,J)**2 120 CONTINUE END IF ELSE DO 130 J=1,M E(J,J) = E(J,J) + ALPHA*TT(1,1)**2 130 CONTINUE END IF ELSE C WD IS AN M BY M MATRIX IF (LD2WD.EQ.1) THEN C THE ARRAY STORED IN WD IS DIAGONAL CALL SZERO(M,M,E,M) DO 140 J=1,M E(J,J) = WD(1,1,J) 140 CONTINUE ELSE C THE ARRAY STORED IN WD IS A FULL POSITIVE SEMIDEFINITE MATRICES DO 160 J1=1,M DO 150 J2=1,M E(J1,J2) = WD(1,J1,J2) 150 CONTINUE 160 CONTINUE END IF IF (TT(1,1).GT.ZERO) THEN IF (LDTT.GE.N) THEN DO 210 J=1,M E(J,J) = E(J,J) + ALPHA*TT(I,J)**2 210 CONTINUE ELSE DO 220 J=1,M E(J,J) = E(J,J) + ALPHA*TT(1,J)**2 220 CONTINUE END IF ELSE DO 230 J=1,M E(J,J) = E(J,J) + ALPHA*TT(1,1)**2 230 CONTINUE END IF END IF ELSE C WD IS A DIAGONAL MATRIX WITH ELEMENTS ABS(WD(1,1,1)) CALL SZERO(M,M,E,M) IF (TT(1,1).GT.ZERO) THEN IF (LDTT.GE.N) THEN DO 310 J=1,M E(J,J) = ABS(WD(1,1,1)) + ALPHA*TT(I,J)**2 310 CONTINUE ELSE DO 320 J=1,M E(J,J) = ABS(WD(1,1,1)) + ALPHA*TT(1,J)**2 320 CONTINUE END IF ELSE DO 330 J=1,M E(J,J) = ABS(WD(1,1,1)) + ALPHA*TT(1,1)**2 330 CONTINUE END IF END IF RETURN END *SETAF SUBROUTINE SETAF + (FCN, + N,M,NP,NQ, + XPLUSD,BETA,EPSMAC,NROW, + PARTMP,PV0, + IFIXB,IFIXX,LDIFX, + ISTOP,NFEV,ETA,NETA, + WRK1,WRK2,WRK6,WRK7) C***BEGIN PROLOGUE SETAF C***REFER TO SODR,SODRC C***ROUTINES CALLED FCN C***DATE WRITTEN 860529 (YYMMDD) C***REVISION DATE 920619 (YYMMDD) C***PURPOSE COMPUTE NOISE AND NUMBER OF GOOD DIGITS IN FUNCTION RESULTS C (ADAPTED FROM STARPAC SUBROUTINE ETAFUN) C***END PROLOGUE SETAF C...SCALAR ARGUMENTS REAL + EPSMAC,ETA INTEGER + ISTOP,LDIFX,M,N,NETA,NFEV,NP,NQ,NROW C...ARRAY ARGUMENTS REAL + BETA(NP),PARTMP(NP),PV0(N,NQ), + WRK1(N,M,NQ),WRK2(N,NQ),WRK6(N,NP,NQ),WRK7(-2:2,NQ),XPLUSD(N,M) INTEGER + IFIXB(NP),IFIXX(LDIFX,M) C...SUBROUTINE ARGUMENTS EXTERNAL + FCN C...LOCAL SCALARS REAL + A,B,FAC,HUNDRD,ONE,P1,P2,P5,STP,TWO,ZERO INTEGER + J,K,L C...INTRINSIC FUNCTIONS INTRINSIC + ABS,INT,LOG10,MAX,SQRT C...DATA STATEMENTS DATA + ZERO,P1,P2,P5,ONE,TWO,HUNDRD + /0.0E0,0.1E0,0.2E0,0.5E0,1.0E0,2.0E0,1.0E2/ C...ROUTINE NAMES USED AS SUBPROGRAM ARGUMENTS C FCN: THE USER SUPPLIED SUBROUTINE FOR EVALUATING THE MODEL. C...VARIABLE DEFINITIONS (ALPHABETICALLY) C A: PARAMETERS OF THE LOCAL FIT. C B: PARAMETERS OF THE LOCAL FIT. C BETA: THE FUNCTION PARAMETERS. C EPSMAC: THE VALUE OF MACHINE PRECISION. C ETA: THE NOISE IN THE MODEL RESULTS. C FAC: A FACTOR USED IN THE COMPUTATIONS. C HUNDRD: THE VALUE 1.0E2. C IFIXB: THE VALUES DESIGNATING WHETHER THE ELEMENTS OF BETA ARE C FIXED AT THEIR INPUT VALUES OR NOT. C IFIXX: THE VALUES DESIGNATING WHETHER THE ELEMENTS OF X ARE C FIXED AT THEIR INPUT VALUES OR NOT. C ISTOP: THE VARIABLE DESIGNATING WHETHER THERE ARE PROBLEMS C COMPUTING THE FUNCTION AT THE CURRENT BETA AND DELTA. C J: AN INDEX VARIABLE. C K: AN INDEX VARIABLE. C L: AN INDEX VARIABLE. C LDIFX: THE LEADING DIMENSION OF ARRAY IFIXX. C M: THE NUMBER OF COLUMNS OF DATA IN THE EXPLANATORY VARIABLE. C N: THE NUMBER OF OBSERVATIONS. C NETA: THE NUMBER OF ACCURATE DIGITS IN THE MODEL RESULTS. C NFEV: THE NUMBER OF FUNCTION EVALUATIONS. C NP: THE NUMBER OF FUNCTION PARAMETERS. C NQ: THE NUMBER OF RESPONSES PER OBSERVATION. C NROW: THE ROW NUMBER AT WHICH THE DERIVATIVE IS TO BE CHECKED. C ONE: THE VALUE 1.0E0. C P1: THE VALUE 0.1E0. C P2: THE VALUE 0.2E0. C P5: THE VALUE 0.5E0. C PARTMP: THE MODEL PARAMETERS. C PV0: THE ORIGINAL PREDICTED VALUES. C STP: A SMALL VALUE USED TO PERTURB THE PARAMETERS. C WRK1: A WORK ARRAY OF (N BY M BY NQ) ELEMENTS. C WRK2: A WORK ARRAY OF (N BY NQ) ELEMENTS. C WRK6: A WORK ARRAY OF (N BY NP BY NQ) ELEMENTS. C WRK7: A WORK ARRAY OF (5 BY NQ) ELEMENTS. C XPLUSD: THE VALUES OF X + DELTA. C ZERO: THE VALUE 0.0E0. C***FIRST EXECUTABLE STATEMENT SETAF STP = HUNDRD*EPSMAC ETA = EPSMAC DO 40 J=-2,2 IF (J.EQ.0) THEN DO 10 L=1,NQ WRK7(J,L) = PV0(NROW,L) 10 CONTINUE ELSE DO 20 K=1,NP IF (IFIXB(1).LT.0) THEN PARTMP(K) = BETA(K) + J*STP*BETA(K) ELSE IF (IFIXB(K).NE.0) THEN PARTMP(K) = BETA(K) + J*STP*BETA(K) ELSE PARTMP(K) = BETA(K) END IF 20 CONTINUE ISTOP = 0 CALL FCN(N,M,NP,NQ, + N,M,NP, + PARTMP,XPLUSD, + IFIXB,IFIXX,LDIFX, + 003,WRK2,WRK6,WRK1,ISTOP) IF (ISTOP.NE.0) THEN RETURN ELSE NFEV = NFEV + 1 END IF DO 30 L=1,NQ WRK7(J,L) = WRK2(NROW,L) 30 CONTINUE END IF 40 CONTINUE DO 100 L=1,NQ A = ZERO B = ZERO DO 50 J=-2,2 A = A + WRK7(J,L) B = B + J*WRK7(J,L) 50 CONTINUE A = P2*A B = P1*B IF ((WRK7(0,L).NE.ZERO) .AND. + (ABS(WRK7(1,L)+WRK7(-1,L)).GT.HUNDRD*EPSMAC)) THEN FAC = ONE/ABS(WRK7(0,L)) ELSE FAC = ONE END IF DO 60 J=-2,2 WRK7(J,L) = ABS((WRK7(J,L)-(A+J*B))*FAC) ETA = MAX(WRK7(J,L),ETA) 60 CONTINUE 100 CONTINUE NETA = MAX(TWO,P5-LOG10(ETA)) RETURN END *SEVJAC SUBROUTINE SEVJAC + (FCN, + ANAJAC,CDJAC, + N,M,NP,NQ, + BETAC,BETA,STPB, + IFIXB,IFIXX,LDIFX, + X,LDX,DELTA,XPLUSD,STPD,LDSTPD, + SSF,TT,LDTT,NETA,FN, + STP,WRK1,WRK2,WRK3,WRK6, + FJACB,ISODR,FJACD,WE1,LDWE,LD2WE, + NJEV,NFEV,ISTOP,INFO) C***BEGIN PROLOGUE SEVJAC C***REFER TO SODR,SODRC C***ROUTINES CALLED FCN,SDOT,SIFIX,SJACCD,SJACFD,SWGHT,SUNPAC,SXPY C***DATE WRITTEN 860529 (YYMMDD) C***REVISION DATE 920304 (YYMMDD) C***PURPOSE COMPUTE THE WEIGHTED JACOBIANS WRT BETA AND DELTA C***END PROLOGUE SEVJAC C...SCALAR ARGUMENTS INTEGER + INFO,ISTOP,LDIFX,LDSTPD,LDTT,LDWE,LDX,LD2WE, + M,N,NETA,NFEV,NJEV,NP,NQ LOGICAL + ANAJAC,CDJAC,ISODR C...ARRAY ARGUMENTS REAL + BETA(NP),BETAC(NP),DELTA(N,M),FJACB(N,NP,NQ),FJACD(N,M,NQ), + FN(N,NQ),SSF(NP),STP(N),STPB(NP),STPD(LDSTPD,M),TT(LDTT,M), + WE1(LDWE,LD2WE,NQ),WRK1(N,M,NQ),WRK2(N,NQ),WRK3(NP), + WRK6(N,NP,NQ),X(LDX,M),XPLUSD(N,M) INTEGER + IFIXB(NP),IFIXX(LDIFX,M) C...SUBROUTINE ARGUMENTS EXTERNAL + FCN C...LOCAL SCALARS INTEGER + IDEVAL,J,K,K1,L REAL + ZERO LOGICAL + ERROR C...EXTERNAL SUBROUTINES EXTERNAL + SIFIX,SJACCD,SJACFD,SWGHT,SUNPAC,SXPY C...EXTERNAL FUNCTIONS REAL + SDOT EXTERNAL + SDOT C...DATA STATEMENTS DATA ZERO + /0.0E0/ C...ROUTINE NAMES USED AS SUBPROGRAM ARGUMENTS C FCN: THE USER-SUPPLIED SUBROUTINE FOR EVALUATING THE MODEL. C...VARIABLE DEFINITIONS (ALPHABETICALLY) C ANAJAC: THE VARIABLE DESIGNATING WHETHER THE JACOBIANS ARE C COMPUTED BY FINITE DIFFERENCES (ANAJAC=FALSE) OR NOT C (ANAJAC=TRUE). C BETA: THE FUNCTION PARAMETERS. C BETAC: THE CURRENT ESTIMATED VALUES OF THE UNFIXED BETA'S. C CDJAC: THE VARIABLE DESIGNATING WHETHER THE JACOBIANS ARE C COMPUTED BY CENTRAL DIFFERENCES (CDJAC=TRUE) OR BY FORWARD C DIFFERENCES (CDJAC=FALSE). C DELTA: THE ESTIMATED VALUES OF DELTA. C ERROR: THE VARIABLE DESIGNATING WHETHER ODRPACK DETECTED NONZERO C VALUES IN ARRAY DELTA IN THE OLS CASE, AND THUS WHETHER C THE USER MAY HAVE OVERWRITTEN IMPORTANT INFORMATION C BY COMPUTING FJACD IN THE OLS CASE. C FJACB: THE JACOBIAN WITH RESPECT TO BETA. C FJACD: THE JACOBIAN WITH RESPECT TO DELTA. C FN: THE PREDICTED VALUES OF THE FUNCTION AT THE CURRENT POINT. C IDEVAL: THE VARIABLE DESIGNATING WHAT COMPUTATIONS ARE TO BE C PERFORMED BY USER-SUPPLIED SUBROUTINE FCN. C IFIXB: THE VALUES DESIGNATING WHETHER THE ELEMENTS OF BETA ARE C FIXED AT THEIR INPUT VALUES OR NOT. C IFIXX: THE VALUES DESIGNATING WHETHER THE ELEMENTS OF DELTA ARE C FIXED AT THEIR INPUT VALUES OR NOT. C INFO: THE VARIABLE DESIGNATING WHY THE COMPUTATIONS WERE STOPPED. C ISTOP: THE VARIABLE DESIGNATING THAT THE USER WISHES THE C COMPUTATIONS STOPPED. C ISODR: THE VARIABLE DESIGNATING WHETHER THE SOLUTION IS BY ODR C (ISODR=TRUE) OR OLS (ISODR=FALSE). C J: AN INDEXING VARIABLE. C K: AN INDEXING VARIABLE. C K1: AN INDEXING VARIABLE. C L: AN INDEXING VARIABLE. C LDIFX: THE LEADING DIMENSION OF ARRAY IFIXX. C LDSTPD: THE LEADING DIMENSION OF ARRAY STPD. C LDTT: THE LEADING DIMENSION OF ARRAY TT. C LDWE: THE LEADING DIMENSION OF ARRAYS WE AND WE1. C LDX: THE LEADING DIMENSION OF ARRAY X. C LD2WE: THE SECOND DIMENSION OF ARRAYS WE AND WE1. C M: THE NUMBER OF COLUMNS OF DATA IN THE INDEPENDENT VARIABLE. C N: THE NUMBER OF OBSERVATIONS. C NETA: THE NUMBER OF ACCURATE DIGITS IN THE FUNCTION RESULTS. C NFEV: THE NUMBER OF FUNCTION EVALUATIONS. C NJEV: THE NUMBER OF JACOBIAN EVALUATIONS. C NP: THE NUMBER OF FUNCTION PARAMETERS. C NQ: THE NUMBER OF RESPONSES PER OBSERVATION. C SSF: THE SCALE USED FOR THE BETA'S. C STP: THE STEP USED FOR COMPUTING FINITE DIFFERENCE C DERIVATIVES WITH RESPECT TO DELTA. C STPB: THE RELATIVE STEP USED FOR COMPUTING FINITE DIFFERENCE C DERIVATIVES WITH RESPECT TO BETA. C STPD: THE RELATIVE STEP USED FOR COMPUTING FINITE DIFFERENCE C DERIVATIVES WITH RESPECT TO DELTA. C TT: THE SCALING VALUES USED FOR DELTA. C WE1: THE SQUARE ROOTS OF THE EPSILON WEIGHTS IN ARRAY WE. C WRK1: A WORK ARRAY OF (N BY M BY NQ) ELEMENTS. C WRK2: A WORK ARRAY OF (N BY NQ) ELEMENTS. C WRK3: A WORK ARRAY OF (NP) ELEMENTS. C WRK6: A WORK ARRAY OF (N BY NP BY NQ) ELEMENTS. C X: THE INDEPENDENT VARIABLE. C XPLUSD: THE VALUES OF X + DELTA. C ZERO: THE VALUE 0.0E0. C***FIRST EXECUTABLE STATEMENT SEVJAC C INSERT CURRENT UNFIXED BETA ESTIMATES INTO BETA CALL SUNPAC(NP,BETAC,BETA,IFIXB) C COMPUTE XPLUSD = X + DELTA CALL SXPY(N,M,X,LDX,DELTA,N,XPLUSD,N) C COMPUTE THE JACOBIAN WRT THE ESTIMATED BETAS (FJACB) AND C THE JACOBIAN WRT DELTA (FJACD) ISTOP = 0 IF (ISODR) THEN IDEVAL = 110 ELSE IDEVAL = 010 END IF IF (ANAJAC) THEN CALL FCN(N,M,NP,NQ, + N,M,NP, + BETA,XPLUSD, + IFIXB,IFIXX,LDIFX, + IDEVAL,WRK2,FJACB,FJACD, + ISTOP) IF (ISTOP.NE.0) THEN RETURN ELSE NJEV = NJEV+1 END IF C MAKE SURE FIXED ELEMENTS OF FJACD ARE ZERO IF (ISODR) THEN DO 10 L=1,NQ CALL SIFIX(N,M,IFIXX,LDIFX,FJACD(1,1,L),N,FJACD(1,1,L),N) 10 CONTINUE END IF ELSE IF (CDJAC) THEN CALL SJACCD(FCN, + N,M,NP,NQ, + BETA,X,LDX,DELTA,XPLUSD,IFIXB,IFIXX,LDIFX, + STPB,STPD,LDSTPD, + SSF,TT,LDTT,NETA,STP,WRK1,WRK2,WRK3,WRK6, + FJACB,ISODR,FJACD,NFEV,ISTOP) ELSE CALL SJACFD(FCN, + N,M,NP,NQ, + BETA,X,LDX,DELTA,XPLUSD,IFIXB,IFIXX,LDIFX, + STPB,STPD,LDSTPD, + SSF,TT,LDTT,NETA,FN,STP,WRK1,WRK2,WRK3,WRK6, + FJACB,ISODR,FJACD,NFEV,ISTOP) END IF IF (ISTOP.LT.0) THEN RETURN ELSE IF (.NOT.ISODR) THEN C TRY TO DETECT WHETHER THE USER HAS COMPUTED JFACD C WITHIN FCN IN THE OLS CASE ERROR = SDOT(N*M,DELTA,1,DELTA,1).NE.ZERO IF (ERROR) THEN INFO = 50300 RETURN END IF END IF C WEIGHT THE JACOBIAN WRT THE ESTIMATED BETAS IF (IFIXB(1).LT.0) THEN DO 20 K=1,NP CALL SWGHT(N,NQ,WE1,LDWE,LD2WE, + FJACB(1,K,1),N*NP,FJACB(1,K,1),N*NP) 20 CONTINUE ELSE K1 = 0 DO 30 K=1,NP IF (IFIXB(K).GE.1) THEN K1 = K1 + 1 CALL SWGHT(N,NQ,WE1,LDWE,LD2WE, + FJACB(1,K,1),N*NP,FJACB(1,K1,1),N*NP) END IF 30 CONTINUE END IF C WEIGHT THE JACOBIAN'S WRT DELTA AS APPROPRIATE IF (ISODR) THEN DO 40 J=1,M CALL SWGHT(N,NQ,WE1,LDWE,LD2WE, + FJACD(1,J,1),N*M,FJACD(1,J,1),N*M) 40 CONTINUE END IF RETURN END *SFCTR SUBROUTINE SFCTR(OKSEMI,A,LDA,N,INFO) C***BEGIN PROLOGUE SFCTR C***REFER TO SODR,SODRC C***ROUTINES CALLED SDOT C***DATE WRITTEN 910706 (YYMMDD) C***REVISION DATE 920619 (YYMMDD) C***PURPOSE FACTOR THE POSITIVE (SEMI)DEFINITE MATRIX A USING A C MODIFIED CHOLESKY FACTORIZATION C (ADAPTED FROM LINPACK SUBROUTINE SPOFA) C***REFERENCES DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W., C *LINPACK USERS GUIDE*, SIAM, 1979. C***END PROLOGUE SFCTR C...SCALAR ARGUMENTS INTEGER INFO,LDA,N LOGICAL OKSEMI C...ARRAY ARGUMENTS REAL A(LDA,N) C...LOCAL SCALARS REAL XI,S,T,TEN,ZERO INTEGER J,K C...EXTERNAL FUNCTIONS EXTERNAL SMPREC,SDOT REAL SMPREC,SDOT C...INTRINSIC FUNCTIONS INTRINSIC SQRT C...DATA STATEMENTS DATA + ZERO,TEN + /0.0E0,10.0E0/ C...VARIABLE DEFINITIONS (ALPHABETICALLY) C A: THE ARRAY TO BE FACTORED. UPON RETURN, A CONTAINS THE C UPPER TRIANGULAR MATRIX R SO THAT A = TRANS(R)*R C WHERE THE STRICT LOWER TRIANGLE IS SET TO ZERO C IF INFO .NE. 0 , THE FACTORIZATION IS NOT COMPLETE. C I: AN INDEXING VARIABLE. C INFO: AN IDICATOR VARIABLE, WHERE IF C INFO = 0 THEN FACTORIZATION WAS COMPLETED C INFO = K SIGNALS AN ERROR CONDITION. THE LEADING MINOR C OF ORDER K IS NOT POSITIVE (SEMI)DEFINITE. C J: AN INDEXING VARIABLE. C LDA: THE LEADING DIMENSION OF ARRAY A. C N: THE NUMBER OF ROWS AND COLUMNS OF DATA IN ARRAY A. C OKSEMI: THE INDICATING WHETHER THE FACTORED ARRAY CAN BE POSITIVE C SEMIDEFINITE (OKSEMI=TRUE) OR WHETHER IT MUST BE FOUND TO C BE POSITIVE DEFINITE (OKSEMI=FALSE). C TEN: THE VALUE 10.0E0. C XI: A VALUE USED TO TEST FOR NON POSITIVE SEMIDEFINITENESS. C ZERO: THE VALUE 0.0E0. C***FIRST EXECUTABLE STATEMENT SFCTR C SET RELATIVE TOLERANCE FOR DETECTING NON POSITIVE SEMIDEFINITENESS. XI = -TEN*SMPREC() C COMPUTE FACTORIZATION, STORING IN UPPER TRIANGULAR PORTION OF A DO 20 J=1,N INFO = J S = ZERO DO 10 K=1,J-1 IF (A(K,K).EQ.ZERO) THEN T = ZERO ELSE T = A(K,J) - SDOT(K-1,A(1,K),1,A(1,J),1) T = T/A(K,K) END IF A(K,J) = T S = S + T*T 10 CONTINUE S = A(J,J) - S C ......EXIT IF (A(J,J).LT.ZERO .OR. S.LT.XI*ABS(A(J,J))) THEN RETURN ELSE IF (.NOT.OKSEMI .AND. S.LE.ZERO) THEN RETURN ELSE IF (S.LE.ZERO) THEN A(J,J) = ZERO ELSE A(J,J) = SQRT(S) END IF 20 CONTINUE INFO = 0 C ZERO OUT LOWER PORTION OF A DO 40 J=2,N DO 30 K=1,J-1 A(J,K) = ZERO 30 CONTINUE 40 CONTINUE RETURN END *SFCTRW SUBROUTINE SFCTRW + (N,M,NQ,NPP, + ISODR, + WE,LDWE,LD2WE,WD,LDWD,LD2WD, + WRK0,WRK4, + WE1,NNZW,INFO) C***BEGIN PROLOGUE SFCTRW C***REFER TO SODR,SODRC C***ROUTINES CALLED SFCTR C***DATE WRITTEN 860529 (YYMMDD) C***REVISION DATE 920619 (YYMMDD) C***PURPOSE CHECK INPUT PARAMETERS, INDICATING ERRORS FOUND USING C NONZERO VALUES OF ARGUMENT INFO AS DESCRIBED IN THE C ODRPACK REFERENCE GUIDE C***END PROLOGUE SFCTRW C...SCALAR ARGUMENTS INTEGER + INFO,LDWD,LDWE,LD2WD,LD2WE, + M,N,NNZW,NPP,NQ LOGICAL + ISODR C...ARRAY ARGUMENTS REAL + WE(LDWE,LD2WE,NQ),WE1(LDWE,LD2WE,NQ),WD(LDWD,LD2WD,M), + WRK0(NQ,NQ),WRK4(M,M) C...LOCAL SCALARS REAL + ZERO INTEGER + I,INF,J,J1,J2,L,L1,L2 LOGICAL + NOTZRO C...EXTERNAL SUBROUTINES EXTERNAL + SFCTR C...DATA STATEMENTS DATA + ZERO + /0.0E0/ C...VARIABLE DEFINITIONS (ALPHABETICALLY) C I: AN INDEXING VARIABLE. C INFO: THE VARIABLE DESIGNATING WHY THE COMPUTATIONS WERE STOPPED. C ISODR: THE VARIABLE DESIGNATING WHETHER THE SOLUTION IS BY ODR C (ISODR=TRUE) OR BY OLS (ISODR=FALSE). C J: AN INDEXING VARIABLE. C J1: AN INDEXING VARIABLE. C J2: AN INDEXING VARIABLE. C L: AN INDEXING VARIABLE. C L1: AN INDEXING VARIABLE. C L2: AN INDEXING VARIABLE. C LAST: THE LAST ROW OF THE ARRAY TO BE ACCESSED. C LDWD: THE LEADING DIMENSION OF ARRAY WD. C LDWE: THE LEADING DIMENSION OF ARRAY WE. C LD2WD: THE SECOND DIMENSION OF ARRAY WD. C LD2WE: THE SECOND DIMENSION OF ARRAY WE. C M: THE NUMBER OF COLUMNS OF DATA IN THE EXPLANATORY VARIABLE. C N: THE NUMBER OF OBSERVATIONS. C NNZW: THE NUMBER OF NONZERO WEIGHTED OBSERVATIONS. C NOTZRO: THE VARIABLE DESIGNATING WHETHER A GIVEN COMPONENT OF THE C WEIGHT ARRAY WE CONTAINS A NONZERO ELEMENT (NOTZRO=FALSE) C OR NOT (NOTZRO=TRUE). C NPP: THE NUMBER OF FUNCTION PARAMETERS BEING ESTIMATED. C NQ: THE NUMBER OF RESPONSES PER OBSERVATIONS. C WE: THE (SQUARED) EPSILON WEIGHTS. C WE1: THE FACTORED EPSILON WEIGHTS, S.T. TRANS(WE1)*WE1 = WE. C WD: THE (SQUARED) DELTA WEIGHTS. C WRK0: A WORK ARRAY OF (NQ BY NQ) ELEMENTS. C WRK4: A WORK ARRAY OF (M BY M) ELEMENTS. C ZERO: THE VALUE 0.0E0. C***FIRST EXECUTABLE STATEMENT SFCTRW C CHECK EPSILON WEIGHTS, AND STORE FACTORIZATION IN WE1 IF (WE(1,1,1).LT.ZERO) THEN C WE CONTAINS A SCALAR WE1(1,1,1) = -SQRT(ABS(WE(1,1,1))) NNZW = N ELSE NNZW = 0 IF (LDWE.EQ.1) THEN IF (LD2WE.EQ.1) THEN C WE CONTAINS A DIAGONAL MATRIX DO 110 L=1,NQ IF (WE(1,1,L).GT.ZERO) THEN NNZW = N WE1(1,1,L) = SQRT(WE(1,1,L)) ELSE IF (WE(1,1,L).LT.ZERO) THEN INFO = 30010 GO TO 300 END IF 110 CONTINUE ELSE C WE CONTAINS A FULL NQ BY NQ SEMIDEFINITE MATRIX DO 130 L1=1,NQ DO 120 L2=L1,NQ WRK0(L1,L2) = WE(1,L1,L2) 120 CONTINUE 130 CONTINUE CALL SFCTR(.TRUE.,WRK0,NQ,NQ,INF) IF (INF.NE.0) THEN INFO = 30010 GO TO 300 ELSE DO 150 L1=1,NQ DO 140 L2=1,NQ WE1(1,L1,L2) = WRK0(L1,L2) 140 CONTINUE IF (WE1(1,L1,L1).NE.ZERO) THEN NNZW = N END IF 150 CONTINUE END IF END IF ELSE IF (LD2WE.EQ.1) THEN C WE CONTAINS AN ARRAY OF DIAGONAL MATRIX DO 220 I=1,N NOTZRO = .FALSE. DO 210 L=1,NQ IF (WE(I,1,L).GT.ZERO) THEN NOTZRO = .TRUE. WE1(I,1,L) = SQRT(WE(I,1,L)) ELSE IF (WE(I,1,L).LT.ZERO) THEN INFO = 30010 GO TO 300 END IF 210 CONTINUE IF (NOTZRO) THEN NNZW = NNZW + 1 END IF 220 CONTINUE ELSE C WE CONTAINS AN ARRAY OF FULL NQ BY NQ SEMIDEFINITE MATRICES DO 270 I=1,N DO 240 L1=1,NQ DO 230 L2=L1,NQ WRK0(L1,L2) = WE(I,L1,L2) 230 CONTINUE 240 CONTINUE CALL SFCTR(.TRUE.,WRK0,NQ,NQ,INF) IF (INF.NE.0) THEN INFO = 30010 GO TO 300 ELSE NOTZRO = .FALSE. DO 260 L1=1,NQ DO 250 L2=1,NQ WE1(I,L1,L2) = WRK0(L1,L2) 250 CONTINUE IF (WE1(I,L1,L1).NE.ZERO) THEN NOTZRO = .TRUE. END IF 260 CONTINUE END IF IF (NOTZRO) THEN NNZW = NNZW + 1 END IF 270 CONTINUE END IF END IF END IF C CHECK FOR A SUFFICIENT NUMBER OF NONZERO EPSILON WEIGHTS IF (NNZW.LT.NPP) THEN INFO = 30020 END IF C CHECK DELTA WEIGHTS 300 CONTINUE IF (.NOT.ISODR .OR. WD(1,1,1).LT.ZERO) THEN C PROBLEM IS NOT ODR, OR WD CONTAINS A SCALAR RETURN ELSE IF (LDWD.EQ.1) THEN IF (LD2WD.EQ.1) THEN C WD CONTAINS A DIAGONAL MATRIX DO 310 J=1,M IF (WD(1,1,J).LE.ZERO) THEN INFO = MAX(30001,INFO+1) RETURN END IF 310 CONTINUE ELSE C WD CONTAINS A FULL M BY M POSITIVE DEFINITE MATRIX DO 330 J1=1,M DO 320 J2=J1,M WRK4(J1,J2) = WD(1,J1,J2) 320 CONTINUE 330 CONTINUE CALL SFCTR(.FALSE.,WRK4,M,M,INF) IF (INF.NE.0) THEN INFO = MAX(30001,INFO+1) RETURN END IF END IF ELSE IF (LD2WD.EQ.1) THEN C WD CONTAINS AN ARRAY OF DIAGONAL MATRICES DO 420 I=1,N DO 410 J=1,M IF (WD(I,1,J).LE.ZERO) THEN INFO = MAX(30001,INFO+1) RETURN END IF 410 CONTINUE 420 CONTINUE ELSE C WD CONTAINS AN ARRAY OF FULL M BY M POSITIVE DEFINITE MATRICES DO 470 I=1,N DO 440 J1=1,M DO 430 J2=J1,M WRK4(J1,J2) = WD(I,J1,J2) 430 CONTINUE 440 CONTINUE CALL SFCTR(.FALSE.,WRK4,M,M,INF) IF (INF.NE.0) THEN INFO = MAX(30001,INFO+1) RETURN END IF 470 CONTINUE END IF END IF END IF RETURN END *SFLAGS SUBROUTINE SFLAGS + (JOB,RESTRT,INITD,DOVCV,REDOJ,ANAJAC,CDJAC,CHKJAC,ISODR,IMPLCT) C***BEGIN PROLOGUE SFLAGS C***REFER TO SODR,SODRC C***ROUTINES CALLED (NONE) C***DATE WRITTEN 860529 (YYMMDD) C***REVISION DATE 920304 (YYMMDD) C***PURPOSE SET FLAGS INDICATING CONDITIONS SPECIFIED BY JOB C***END PROLOGUE SFLAGS C...SCALAR ARGUMENTS INTEGER + JOB LOGICAL + ANAJAC,CDJAC,CHKJAC,DOVCV,IMPLCT,INITD,ISODR,REDOJ,RESTRT C...LOCAL SCALARS INTEGER + J C...INTRINSIC FUNCTIONS INTRINSIC + MOD C...VARIABLE DEFINITIONS (ALPHABETICALLY) C ANAJAC: THE VARIABLE DESIGNATING WHETHER THE JACOBIANS ARE COMPUTED C BY FINITE DIFFERENCES (ANAJAC=FALSE) OR NOT (ANAJAC=TRUE). C CDJAC: THE VARIABLE DESIGNATING WHETHER THE JACOBIANS ARE COMPUTED C BY CENTRAL DIFFERENCES (CDJAC=TRUE) OR BY FORWARD C DIFFERENCES (CDJAC=FALSE). C CHKJAC: THE VARIABLE DESIGNATING WHETHER THE USER-SUPPLIED C JACOBIANS ARE TO BE CHECKED (CHKJAC=TRUE) OR NOT C (CHKJAC=FALSE). C DOVCV: THE VARIABLE DESIGNATING WHETHER THE COVARIANCE MATRIX IS C TO BE COMPUTED (DOVCV=TRUE) OR NOT (DOVCV=FALSE). C IMPLCT: THE VARIABLE DESIGNATING WHETHER THE SOLUTION IS BY C IMPLICIT ODR (IMPLCT=TRUE) OR EXPLICIT ODR (IMPLCT=FALSE). C INITD: THE VARIABLE DESIGNATING WHETHER DELTA IS TO BE INITIALIZED C TO ZERO (INITD=TRUE) OR TO THE FIRST N BY M ELEMENTS OF C ARRAY WORK (INITD=FALSE). C ISODR: THE VARIABLE DESIGNATING WHETHER THE SOLUTION IS BY ODR C (ISODR=TRUE) OR BY OLS (ISODR=FALSE). C J: THE VALUE OF A SPECIFIC DIGIT OF JOB. C JOB: THE VARIABLE CONTROLING PROBLEM INITIALIZATION AND C COMPUTATIONAL METHOD. C REDOJ: THE VARIABLE DESIGNATING WHETHER THE JACOBIAN MATRIX IS TO C BE RECOMPUTED FOR THE COMPUTATION OF THE COVARIANCE MATRIX C (REDOJ=TRUE) OR NOT (REDOJ=FALSE). C RESTRT: THE VARIABLE DESIGNATING WHETHER THE CALL IS A RESTART C (RESTRT=TRUE) OR NOT (RESTRT=FALSE). C***FIRST EXECUTABLE STATEMENT SFLAGS IF (JOB.GE.0) THEN RESTRT= JOB.GE.10000 INITD = MOD(JOB,10000)/1000.EQ.0 J = MOD(JOB,1000)/100 IF (J.EQ.0) THEN DOVCV = .TRUE. REDOJ = .TRUE. ELSE IF (J.EQ.1) THEN DOVCV = .TRUE. REDOJ = .FALSE. ELSE DOVCV = .FALSE. REDOJ = .FALSE. END IF J = MOD(JOB,100)/10 IF (J.EQ.0) THEN ANAJAC = .FALSE. CDJAC = .FALSE. CHKJAC = .FALSE. ELSE IF (J.EQ.1) THEN ANAJAC = .FALSE. CDJAC = .TRUE. CHKJAC = .FALSE. ELSE IF (J.EQ.2) THEN ANAJAC = .TRUE. CDJAC = .FALSE. CHKJAC = .TRUE. ELSE ANAJAC = .TRUE. CDJAC = .FALSE. CHKJAC = .FALSE. END IF J = MOD(JOB,10) IF (J.EQ.0) THEN ISODR = .TRUE. IMPLCT = .FALSE. ELSE IF (J.EQ.1) THEN ISODR = .TRUE. IMPLCT = .TRUE. ELSE ISODR = .FALSE. IMPLCT = .FALSE. END IF ELSE RESTRT = .FALSE. INITD = .TRUE. DOVCV = .TRUE. REDOJ = .TRUE. ANAJAC = .FALSE. CDJAC = .FALSE. CHKJAC = .FALSE. ISODR = .TRUE. IMPLCT = .FALSE. END IF RETURN END *SHSTEP REAL FUNCTION SHSTEP + (ITYPE,NETA,I,J,STP,LDSTP) C***BEGIN PROLOGUE SHSTEP C***REFER TO SODR,SODRC C***ROUTINES CALLED (NONE) C***DATE WRITTEN 860529 (YYMMDD) C***REVISION DATE 920304 (YYMMDD) C***PURPOSE SET RELATIVE STEP SIZE FOR FINITE DIFFERENCE DERIVATIVES C***END PROLOGUE SHSTEP C...SCALAR ARGUMENTS INTEGER + I,ITYPE,J,LDSTP,NETA C...ARRAY ARGUMENTS REAL + STP(LDSTP,J) C...LOCAL SCALARS REAL + TEN,THREE,TWO,ZERO C...DATA STATEMENTS DATA + ZERO,TWO,THREE,TEN + /0.0E0,2.0E0,3.0E0,10.0E0/ C...VARIABLE DEFINITIONS (ALPHABETICALLY) C I: AN IDENTIFIER FOR SELECTING USER SUPPLIED STEP SIZES. C ITYPE: THE FINITE DIFFERENCE METHOD BEING USED, WHERE C ITYPE = 0 INDICATES FORWARD FINITE DIFFERENCES, AND C ITYPE = 1 INDICATES CENTRAL FINITE DIFFERENCES. C J: AN IDENTIFIER FOR SELECTING USER SUPPLIED STEP SIZES. C LDSTP: THE LEADING DIMENSION OF ARRAY STP. C NETA: THE NUMBER OF GOOD DIGITS IN THE FUNCTION RESULTS. C STP: THE STEP SIZE FOR THE FINITE DIFFERENCE DERIVATIVE. C TEN: THE VALUE 10.0E0. C THREE: THE VALUE 3.0E0. C TWO: THE VALUE 2.0E0. C ZERO: THE VALUE 0.0E0. C***FIRST EXECUTABLE STATEMENT SHSTEP C SET SHSTEP TO RELATIVE FINITE DIFFERENCE STEP SIZE IF (STP(1,1).LE.ZERO) THEN IF (ITYPE.EQ.0) THEN C USE DEFAULT FORWARD FINITE DIFFERENCE STEP SIZE SHSTEP = TEN**(-ABS(NETA)/TWO - TWO) ELSE C USE DEFAULT CENTRAL FINITE DIFFERENCE STEP SIZE SHSTEP = TEN**(-ABS(NETA)/THREE) END IF ELSE IF (LDSTP.EQ.1) THEN SHSTEP = STP(1,J) ELSE SHSTEP = STP(I,J) END IF RETURN END *SIFIX SUBROUTINE SIFIX + (N,M,IFIX,LDIFIX,T,LDT,TFIX,LDTFIX) C***BEGIN PROLOGUE SIFIX C***REFER TO SODR,SODRC C***ROUTINES CALLED (NONE) C***DATE WRITTEN 910612 (YYMMDD) C***REVISION DATE 920304 (YYMMDD) C***PURPOSE SET ELEMENTS OF T TO ZERO ACCORDING TO IFIX C***END PROLOGUE SIFIX C...SCALAR ARGUMENTS INTEGER + LDIFIX,LDT,LDTFIX,M,N C...ARRAY ARGUMENTS REAL + T(LDT,M),TFIX(LDTFIX,M) INTEGER + IFIX(LDIFIX,M) C...LOCAL SCALARS REAL + ZERO INTEGER + I,J C...INTRINSIC FUNCTIONS INTRINSIC + ABS C...DATA STATEMENTS DATA + ZERO + /0.0E0/ C...VARIABLE DEFINITIONS (ALPHABETICALLY) C I: AN INDEXING VARIABLE. C IFIX: THE ARRAY DESIGNATING WHETHER AN ELEMENT OF T IS TO BE C SET TO ZERO. C J: AN INDEXING VARIABLE. C LDT: THE LEADING DIMENSION OF ARRAY T. C LDIFIX: THE LEADING DIMENSION OF ARRAY IFIX. C LDTFIX: THE LEADING DIMENSION OF ARRAY TFIX. C M: THE NUMBER OF COLUMNS OF DATA IN THE ARRAY. C N: THE NUMBER OF ROWS OF DATA IN THE ARRAY. C T: THE ARRAY BEING SET TO ZERO ACCORDING TO THE ELEMENTS C OF IFIX. C TFIX: THE RESULTING ARRAY. C ZERO: THE VALUE 0.0E0. C***FIRST EXECUTABLE STATEMENT SIFIX IF (N.EQ.0 .OR. M.EQ.0) RETURN IF (IFIX(1,1).GE.ZERO) THEN IF (LDIFIX.GE.N) THEN DO 20 J=1,M DO 10 I=1,N IF (IFIX(I,J).EQ.0) THEN TFIX(I,J) = ZERO ELSE TFIX(I,J) = T(I,J) END IF 10 CONTINUE 20 CONTINUE ELSE DO 100 J=1,M IF (IFIX(1,J).EQ.0) THEN DO 30 I=1,N TFIX(I,J) = ZERO 30 CONTINUE ELSE DO 90 I=1,N TFIX(I,J) = T(I,J) 90 CONTINUE END IF 100 CONTINUE END IF END IF RETURN END *SINIWK SUBROUTINE SINIWK + (N,M,NP,WORK,LWORK,IWORK,LIWORK, + X,LDX,IFIXX,LDIFX,SCLD,LDSCLD, + BETA,SCLB, + SSTOL,PARTOL,MAXIT,TAUFAC, + JOB,IPRINT,LUNERR,LUNRPT, + EPSMAI,SSTOLI,PARTLI,MAXITI,TAUFCI, + JOBI,IPRINI,LUNERI,LUNRPI, + SSFI,TTI,LDTTI,DELTAI) C***BEGIN PROLOGUE SINIWK C***REFER TO SODR,SODRC C***ROUTINES CALLED SFLAGS,SMPREC,SSCLB,SSCLD,SZERO C***DATE WRITTEN 860529 (YYMMDD) C***REVISION DATE 920304 (YYMMDD) C***PURPOSE INITIALIZE WORK VECTORS AS NECESSARY C***END PROLOGUE SINIWK C...SCALAR ARGUMENTS REAL + PARTOL,SSTOL,TAUFAC INTEGER + DELTAI,EPSMAI,IPRINI,IPRINT,JOB,JOBI,LDIFX, + LDSCLD,LDTTI,LDX,LIWORK,LUNERI,LUNERR,LUNRPI,LUNRPT,LWORK,M, + MAXIT,MAXITI,N,NP,PARTLI,SSFI,SSTOLI,TAUFCI,TTI C...ARRAY ARGUMENTS REAL + BETA(NP),SCLB(NP),SCLD(LDSCLD,M),WORK(LWORK),X(LDX,M) INTEGER + IFIXX(LDIFX,M),IWORK(LIWORK) C...LOCAL SCALARS REAL + ONE,THREE,TWO,ZERO INTEGER + I,J LOGICAL + ANAJAC,CDJAC,CHKJAC,DOVCV,IMPLCT,INITD,ISODR,REDOJ,RESTRT C...EXTERNAL FUNCTIONS REAL + SMPREC EXTERNAL + SMPREC C...EXTERNAL SUBROUTINES EXTERNAL + SCOPY,SFLAGS,SSCLB,SSCLD,SZERO C...INTRINSIC FUNCTIONS INTRINSIC + MIN,SQRT C...DATA STATEMENTS DATA + ZERO,ONE,TWO,THREE + /0.0E0,1.0E0,2.0E0,3.0E0/ C...VARIABLE DEFINITIONS (ALPHABETICALLY) C ANAJAC: THE VARIABLE DESIGNATING WHETHER THE JACOBIANS ARE C COMPUTED BY FINITE DIFFERENCES (ANAJAC=FALSE) OR NOT C (ANAJAC=TRUE). C BETA: THE FUNCTION PARAMETERS. C CDJAC: THE VARIABLE DESIGNATING WHETHER THE JACOBIANS ARE C COMPUTED BY CENTRAL DIFFERENCES (CDJAC=TRUE) OR BY FORWARD C DIFFERENCES (CDJAC=FALSE). C CHKJAC: THE VARIABLE DESIGNATING WHETHER THE USER-SUPPLIED C JACOBIANS ARE TO BE CHECKED (CHKJAC=TRUE) OR NOT C (CHKJAC=FALSE). C DELTAI: THE STARTING LOCATION IN ARRAY WORK OF ARRAY DELTA. C DOVCV: THE VARIABLE DESIGNATING WHETHER THE COVARIANCE MATRIX IS C TO BE COMPUTED (DOVCV=TRUE) OR NOT (DOVCV=FALSE). C EPSMAI: THE LOCATION IN ARRAY WORK OF VARIABLE EPSMAC. C I: AN INDEXING VARIABLE. C IFIXX: THE VALUES DESIGNATING WHETHER THE ELEMENTS OF X ARE FIXED C AT THEIR INPUT VALUES OR NOT. C IMPLCT: THE VARIABLE DESIGNATING WHETHER THE SOLUTION IS BY C IMPLICIT ODR (IMPLCT=TRUE) OR EXPLICIT ODR (IMPLCT=FALSE). C INITD: THE VARIABLE DESIGNATING WHETHER DELTA IS TO BE INITIALIZED C TO ZERO (INITD=TRUE) OR TO THE VALUES IN THE FIRST N BY M C ELEMENTS OF ARRAY WORK (INITD=FALSE). C IPRINI: THE LOCATION IN ARRAY IWORK OF VARIABLE IPRINT. C IPRINT: THE PRINT CONTROL VARIABLE. C ISODR: THE VARIABLE DESIGNATING WHETHER THE SOLUTION IS BY ODR C (ISODR=TRUE) OR BY OLS (ISODR=FALSE). C IWORK: THE INTEGER WORK SPACE. C J: AN INDEXING VARIABLE. C JOB: THE VARIABLE CONTROLING PROBLEM INITIALIZATION AND C COMPUTATIONAL METHOD. C JOBI: THE LOCATION IN ARRAY IWORK OF VARIABLE JOB. C LDIFX: THE LEADING DIMENSION OF ARRAY IFIXX. C LDSCLD: THE LEADING DIMENSION OF ARRAY SCLD. C LDTTI: THE LEADING DIMENSION OF ARRAY TT. C LDX: THE LEADING DIMENSION OF ARRAY X. C LIWORK: THE LENGTH OF VECTOR IWORK. C LUNERI: THE LOCATION IN ARRAY IWORK OF VARIABLE LUNERR. C LUNERR: THE LOGICAL UNIT NUMBER USED FOR ERROR MESSAGES. C LUNRPI: THE LOCATION IN ARRAY IWORK OF VARIABLE LUNRPT. C LUNRPT: THE LOGICAL UNIT NUMBER USED FOR COMPUTATION REPORTS. C LWORK: THE LENGTH OF VECTOR WORK. C M: THE NUMBER OF COLUMNS OF DATA IN THE INDEPENDENT VARIABLE. C MAXIT: THE MAXIMUM NUMBER OF ITERATIONS ALLOWED. C MAXITI: THE LOCATION IN ARRAY IWORK OF VARIABLE MAXIT. C N: THE NUMBER OF OBSERVATIONS. C NP: THE NUMBER OF FUNCTION PARAMETERS. C ONE: THE VALUE 1.0E0. C PARTLI: THE LOCATION IN ARRAY WORK OF VARIABLE PARTOL. C PARTOL: THE PARAMETER CONVERGENCE STOPPING CRITERIA. C REDOJ: THE VARIABLE DESIGNATING WHETHER THE JACOBIAN MATRIX IS TO C BE RECOMPUTED FOR THE COMPUTATION OF THE COVARIANCE MATRIX C (REDOJ=TRUE) OR NOT (REDOJ=FALSE). C RESTRT: THE VARIABLE DESIGNATING WHETHER THE CALL IS A RESTART C (RESTRT=TRUE) OR NOT (RESTRT=FALSE). C SCLB: THE SCALING VALUES FOR BETA. C SCLD: THE SCALING VALUES FOR DELTA. C SSFI: THE STARTING LOCATION IN ARRAY WORK OF ARRAY SSF. C SSTOL: THE SUM-OF-SQUARES CONVERGENCE STOPPING CRITERIA. C SSTOLI: THE LOCATION IN ARRAY WORK OF VARIABLE SSTOL. C TAUFAC: THE FACTOR USED TO COMPUTE THE INITIAL TRUST REGION C DIAMETER. C TAUFCI: THE LOCATION IN ARRAY WORK OF VARIABLE TAUFAC. C THREE: THE VALUE 3.0E0. C TTI: THE STARTING LOCATION IN ARRAY WORK OF THE ARRAY TT. C TWO: THE VALUE 2.0E0. C WORK: THE REAL WORK SPACE. C X: THE INDEPENDENT VARIABLE. C ZERO: THE VALUE 0.0E0. C***FIRST EXECUTABLE STATEMENT SINIWK CALL SFLAGS(JOB,RESTRT,INITD,DOVCV,REDOJ, + ANAJAC,CDJAC,CHKJAC,ISODR,IMPLCT) C STORE VALUE OF MACHINE PRECISION IN WORK VECTOR WORK(EPSMAI) = SMPREC() C SET TOLERANCE FOR STOPPING CRITERIA BASED ON THE CHANGE IN THE C PARAMETERS (SEE ALSO SUBPROGRAM SODCNT) IF (PARTOL.LT.ZERO) THEN WORK(PARTLI) = WORK(EPSMAI)**(TWO/THREE) ELSE WORK(PARTLI) = MIN(PARTOL, ONE) END IF C SET TOLERANCE FOR STOPPING CRITERIA BASED ON THE CHANGE IN THE C SUM OF SQUARES OF THE WEIGHTED OBSERVATIONAL ERRORS IF (SSTOL.LT.ZERO) THEN WORK(SSTOLI) = SQRT(WORK(EPSMAI)) ELSE WORK(SSTOLI) = MIN(SSTOL, ONE) END IF C SET FACTOR FOR COMPUTING TRUST REGION DIAMETER AT FIRST ITERATION IF (TAUFAC.LE.ZERO) THEN WORK(TAUFCI) = ONE ELSE WORK(TAUFCI) = MIN(TAUFAC, ONE) END IF C SET MAXIMUM NUMBER OF ITERATIONS IF (MAXIT.LT.0) THEN IWORK(MAXITI) = 50 ELSE IWORK(MAXITI) = MAXIT END IF C STORE PROBLEM INITIALIZATION AND COMPUTATIONAL METHOD CONTROL C VARIABLE IF (JOB.LE.0) THEN IWORK(JOBI) = 0 ELSE IWORK(JOBI) = JOB END IF C SET PRINT CONTROL IF (IPRINT.LT.0) THEN IWORK(IPRINI) = 2001 ELSE IWORK(IPRINI) = IPRINT END IF C SET LOGICAL UNIT NUMBER FOR ERROR MESSAGES IF (LUNERR.LT.0) THEN IWORK(LUNERI) = 6 ELSE IWORK(LUNERI) = LUNERR END IF C SET LOGICAL UNIT NUMBER FOR COMPUTATION REPORTS IF (LUNRPT.LT.0) THEN IWORK(LUNRPI) = 6 ELSE IWORK(LUNRPI) = LUNRPT END IF C COMPUTE SCALING FOR BETA'S AND DELTA'S IF (SCLB(1).LE.ZERO) THEN CALL SSCLB(NP,BETA,WORK(SSFI)) ELSE CALL SCOPY(NP,SCLB,1,WORK(SSFI),1) END IF IF (ISODR) THEN IF (SCLD(1,1).LE.ZERO) THEN IWORK(LDTTI) = N CALL SSCLD(N,M,X,LDX,WORK(TTI),IWORK(LDTTI)) ELSE IF (LDSCLD.EQ.1) THEN IWORK(LDTTI) = 1 CALL SCOPY(M,SCLD(1,1),1,WORK(TTI),1) ELSE IWORK(LDTTI) = N DO 10 J=1,M CALL SCOPY(N,SCLD(1,J),1, + WORK(TTI+(J-1)*IWORK(LDTTI)),1) 10 CONTINUE END IF END IF END IF C INITIALIZE DELTA'S AS NECESSARY IF (ISODR) THEN IF (INITD) THEN CALL SZERO(N,M,WORK(DELTAI),N) ELSE IF (IFIXX(1,1).GE.0) THEN IF (LDIFX.EQ.1) THEN DO 20 J=1,M IF (IFIXX(1,J).EQ.0) THEN CALL SZERO(N,1,WORK(DELTAI+(J-1)*N),N) END IF 20 CONTINUE ELSE DO 40 J=1,M DO 30 I=1,N IF (IFIXX(I,J).EQ.0) THEN WORK(DELTAI-1+I+(J-1)*N) = ZERO END IF 30 CONTINUE 40 CONTINUE END IF END IF END IF ELSE CALL SZERO(N,M,WORK(DELTAI),N) END IF RETURN END *SIWINF SUBROUTINE SIWINF + (M,NP,NQ, + MSGBI,MSGDI,IFIX2I,ISTOPI, + NNZWI,NPPI,IDFI, + JOBI,IPRINI,LUNERI,LUNRPI, + NROWI,NTOLI,NETAI, + MAXITI,NITERI,NFEVI,NJEVI,INT2I,IRANKI,LDTTI, + LIWKMN) C***BEGIN PROLOGUE SIWINF C***REFER TO SODR,SODRC C***ROUTINES CALLED (NONE) C***DATE WRITTEN 860529 (YYMMDD) C***REVISION DATE 920304 (YYMMDD) C***PURPOSE SET STORAGE LOCATIONS WITHIN INTEGER WORK SPACE C***END PROLOGUE SIWINF C...SCALAR ARGUMENTS INTEGER + IDFI,INT2I,IPRINI,IRANKI,ISTOPI,JOBI,IFIX2I,LDTTI,LIWKMN, + LUNERI,LUNRPI,M,MAXITI,MSGBI,MSGDI,NETAI,NFEVI,NITERI,NJEVI, + NNZWI,NP,NPPI,NQ,NROWI,NTOLI C...VARIABLE DEFINITIONS (ALPHABETICALLY) C IDFI: THE LOCATION IN ARRAY IWORK OF VARIABLE IDF. C IFIX2I: THE STARTING LOCATION IN ARRAY IWORK OF ARRAY IFIX2. C INT2I: THE LOCATION IN ARRAY IWORK OF VARIABLE INT2. C IPRINI: THE LOCATION IN ARRAY IWORK OF VARIABLE IPRINT. C IRANKI: THE LOCATION IN ARRAY IWORK OF VARIABLE IRANK. C ISTOPI: THE LOCATION IN ARRAY IWORK OF VARIABLE ISTOP. C JOBI: THE LOCATION IN ARRAY IWORK OF VARIABLE JOB. C LDTTI: THE LOCATION IN ARRAY IWORK OF VARIABLE LDTT. C LIWKMN: THE MINIMUM ACCEPTABLE LENGTH OF ARRAY IWORK. C LUNERI: THE LOCATION IN ARRAY IWORK OF VARIABLE LUNERR. C LUNRPI: THE LOCATION IN ARRAY IWORK OF VARIABLE LUNRPT. C M: THE NUMBER OF COLUMNS OF DATA IN THE INDEPENDENT VARIABLE. C MAXITI: THE LOCATION IN ARRAY IWORK OF VARIABLE MAXIT. C MSGBI: THE STARTING LOCATION IN ARRAY IWORK OF ARRAY MSGB. C MSGDI: THE STARTING LOCATION IN ARRAY IWORK OF ARRAY MSGD. C NETAI: THE LOCATION IN ARRAY IWORK OF VARIABLE NETA. C NFEVI: THE LOCATION IN ARRAY IWORK OF VARIABLE NFEV. C NITERI: THE LOCATION IN ARRAY IWORK OF VARIABEL NITER. C NJEVI: THE LOCATION IN ARRAY IWORK OF VARIABLE NJEV. C NNZWI: THE LOCATION IN ARRAY IWORK OF VARIABLE NNZW. C NP: THE NUMBER OF FUNCTION PARAMETERS. C NPPI: THE LOCATION IN ARRAY IWORK OF VARIABLE NPP. C NQ: THE NUMBER OF RESPONSES PER OBSERVATION. C NROWI: THE LOCATION IN ARRAY IWORK OF VARIABLE NROW. C NTOLI: THE LOCATION IN ARRAY IWORK OF VARIABLE NTOL. C***FIRST EXECUTABLE STATEMENT SIWINF IF (NP.GE.1 .AND. M.GE.1) THEN MSGBI = 1 MSGDI = MSGBI + NQ*NP+1 IFIX2I = MSGDI + NQ*M+1 ISTOPI = IFIX2I + NP NNZWI = ISTOPI + 1 NPPI = NNZWI + 1 IDFI = NPPI + 1 JOBI = IDFI + 1 IPRINI = JOBI + 1 LUNERI = IPRINI + 1 LUNRPI = LUNERI + 1 NROWI = LUNRPI + 1 NTOLI = NROWI + 1 NETAI = NTOLI + 1 MAXITI = NETAI + 1 NITERI = MAXITI + 1 NFEVI = NITERI + 1 NJEVI = NFEVI + 1 INT2I = NJEVI + 1 IRANKI = INT2I + 1 LDTTI = IRANKI + 1 LIWKMN = LDTTI ELSE MSGBI = 1 MSGDI = 1 IFIX2I = 1 ISTOPI = 1 NNZWI = 1 NPPI = 1 IDFI = 1 JOBI = 1 IPRINI = 1 LUNERI = 1 LUNRPI = 1 NROWI = 1 NTOLI = 1 NETAI = 1 MAXITI = 1 NITERI = 1 NFEVI = 1 NJEVI = 1 INT2I = 1 IRANKI = 1 LDTTI = 1 LIWKMN = 1 END IF RETURN END *SJACCD SUBROUTINE SJACCD + (FCN, + N,M,NP,NQ, + BETA,X,LDX,DELTA,XPLUSD,IFIXB,IFIXX,LDIFX, + STPB,STPD,LDSTPD, + SSF,TT,LDTT,NETA,STP,WRK1,WRK2,WRK3,WRK6, + FJACB,ISODR,FJACD,NFEV,ISTOP) C***BEGIN PROLOGUE SJACCD C***REFER TO SODR,SODRC C***ROUTINES CALLED FCN,SHSTEP,SZERO C***DATE WRITTEN 860529 (YYMMDD) C***REVISION DATE 920619 (YYMMDD) C***PURPOSE COMPUTE CENTRAL DIFFERENCE APPROXIMATIONS TO THE C JACOBIAN WRT THE ESTIMATED BETAS AND WRT THE DELTAS C***END PROLOGUE SJACCD C...SCALAR ARGUMENTS INTEGER + ISTOP,LDIFX,LDSTPD,LDTT,LDX,M,N,NETA,NFEV,NP,NQ LOGICAL + ISODR C...ARRAY ARGUMENTS REAL + BETA(NP),DELTA(N,M),FJACB(N,NP,NQ),FJACD(N,M,NQ), + SSF(NP),STP(N),STPB(NP),STPD(LDSTPD,M),TT(LDTT,M), + WRK1(N,M,NQ),WRK2(N,NQ),WRK3(NP),WRK6(N,NP,NQ), + X(LDX,M),XPLUSD(N,M) INTEGER + IFIXB(NP),IFIXX(LDIFX,M) C...SUBROUTINE ARGUMENTS EXTERNAL + FCN C...LOCAL SCALARS REAL + BETAK,ONE,TYPJ,ZERO INTEGER + I,J,K,L LOGICAL + DOIT,SETZRO C...EXTERNAL SUBROUTINES EXTERNAL + SZERO C...EXTERNAL FUNCTIONS REAL + SHSTEP EXTERNAL + SHSTEP C...INTRINSIC FUNCTIONS INTRINSIC + ABS,MAX,SIGN,SQRT C...DATA STATEMENTS DATA + ZERO,ONE + /0.0E0,1.0E0/ C...ROUTINE NAMES USED AS SUBPROGRAM ARGUMENTS C FCN: THE USER SUPPLIED SUBROUTINE FOR EVALUATING THE MODEL. C...VARIABLE DEFINITIONS (ALPHABETICALLY) C BETA: THE FUNCTION PARAMETERS. C BETAK: THE K-TH FUNCTION PARAMETER. C DELTA: THE ESTIMATED ERRORS IN THE EXPLANATORY VARIABLES. C DOIT: THE VARIABLE DESIGNATING WHETHER THE DERIVATIVE WRT A GIVEN C BETA OR DELTA NEEDS TO BE COMPUTED (DOIT=TRUE) OR NOT C (DOIT=FALSE). C FJACB: THE JACOBIAN WITH RESPECT TO BETA. C FJACD: THE JACOBIAN WITH RESPECT TO DELTA. C I: AN INDEXING VARIABLE. C IFIXB: THE VALUES DESIGNATING WHETHER THE ELEMENTS OF BETA ARE C FIXED AT THEIR INPUT VALUES OR NOT. C IFIXX: THE VALUES DESIGNATING WHETHER THE ELEMENTS OF X ARE FIXED C AT THEIR INPUT VALUES OR NOT. C ISODR: THE VARIABLE DESIGNATING WHETHER THE SOLUTION IS BY ODR C (ISODR=TRUE) OR BY OLS (ISODR=FALSE). C ISTOP: THE VARIABLE DESIGNATING WHETHER THERE ARE PROBLEMS C COMPUTING THE FUNCTION AT THE CURRENT BETA AND DELTA. C J: AN INDEXING VARIABLE. C K: AN INDEXING VARIABLE. C L: AN INDEXING VARIABLE. C LDIFX: THE LEADING DIMENSION OF ARRAY IFIXX. C LDSTPD: THE LEADING DIMENSION OF ARRAY STPD. C LDTT: THE LEADING DIMENSION OF ARRAY TT. C LDX: THE LEADING DIMENSION OF ARRAY X. C M: THE NUMBER OF COLUMNS OF DATA IN THE EXPLANATORY VARIABLE. C N: THE NUMBER OF OBSERVATIONS. C NETA: THE NUMBER OF GOOD DIGITS IN THE FUNCTION RESULTS. C NFEV: THE NUMBER OF FUNCTION EVALUATIONS. C NP: THE NUMBER OF FUNCTION PARAMETERS. C ONE: THE VALUE 1.0E0. C SETZRO: THE VARIABLE DESIGNATING WHETHER THE DERIVATIVE WRT SOME C DELTA NEEDS TO BE SET TO ZERO (SETZRO=TRUE) OR NOT C (SETZRO=FALSE). C SSF: THE SCALING VALUES USED FOR BETA. C STP: THE STEP USED FOR COMPUTING FINITE DIFFERENCE C DERIVATIVES WITH RESPECT TO EACH DELTA. C STPB: THE RELATIVE STEP USED FOR COMPUTING FINITE DIFFERENCE C DERIVATIVES WITH RESPECT TO EACH BETA. C STPD: THE RELATIVE STEP USED FOR COMPUTING FINITE DIFFERENCE C DERIVATIVES WITH RESPECT TO EACH DELTA. C TT: THE SCALING VALUES USED FOR DELTA. C TYPJ: THE TYPICAL SIZE OF THE J-TH UNKNOWN BETA OR DELTA. C X: THE EXPLANATORY VARIABLE. C XPLUSD: THE VALUES OF X + DELTA. C WRK1: A WORK ARRAY OF (N BY M BY NQ) ELEMENTS. C WRK2: A WORK ARRAY OF (N BY NQ) ELEMENTS. C WRK3: A WORK ARRAY OF (NP) ELEMENTS. C WRK6: A WORK ARRAY OF (N BY NP BY NQ) ELEMENTS. C ZERO: THE VALUE 0.0E0. C***FIRST EXECUTABLE STATEMENT SJACCD C COMPUTE THE JACOBIAN WRT THE ESTIMATED BETAS DO 60 K=1,NP IF (IFIXB(1).GE.0) THEN IF (IFIXB(K).EQ.0) THEN DOIT = .FALSE. ELSE DOIT = .TRUE. END IF ELSE DOIT = .TRUE. END IF IF (.NOT.DOIT) THEN DO 10 L=1,NQ CALL SZERO(N,1,FJACB(1,K,L),N) 10 CONTINUE ELSE BETAK = BETA(K) IF (BETAK.EQ.ZERO) THEN IF (SSF(1).LT.ZERO) THEN TYPJ = ONE/ABS(SSF(1)) ELSE TYPJ = ONE/SSF(K) END IF ELSE TYPJ = ABS(BETAK) END IF WRK3(K) = BETAK + + SIGN(ONE,BETAK)*TYPJ*SHSTEP(1,NETA,1,K,STPB,1) WRK3(K) = WRK3(K) - BETAK BETA(K) = BETAK + WRK3(K) ISTOP = 0 CALL FCN(N,M,NP,NQ, + N,M,NP, + BETA,XPLUSD, + IFIXB,IFIXX,LDIFX, + 001,WRK2,WRK6,WRK1, + ISTOP) IF (ISTOP.NE.0) THEN RETURN ELSE NFEV = NFEV + 1 DO 30 L=1,NQ DO 20 I=1,N FJACB(I,K,L) = WRK2(I,L) 20 CONTINUE 30 CONTINUE END IF BETA(K) = BETAK - WRK3(K) ISTOP = 0 CALL FCN(N,M,NP,NQ, + N,M,NP, + BETA,XPLUSD, + IFIXB,IFIXX,LDIFX, + 001,WRK2,WRK6,WRK1, + ISTOP) IF (ISTOP.NE.0) THEN RETURN ELSE NFEV = NFEV + 1 END IF DO 50 L=1,NQ DO 40 I=1,N FJACB(I,K,L) = (FJACB(I,K,L)-WRK2(I,L))/(2*WRK3(K)) 40 CONTINUE 50 CONTINUE BETA(K) = BETAK END IF 60 CONTINUE C COMPUTE THE JACOBIAN WRT THE X'S IF (ISODR) THEN DO 220 J=1,M IF (IFIXX(1,1).LT.0) THEN DOIT = .TRUE. SETZRO = .FALSE. ELSE IF (LDIFX.EQ.1) THEN IF (IFIXX(1,J).EQ.0) THEN DOIT = .FALSE. ELSE DOIT = .TRUE. END IF SETZRO = .FALSE. ELSE DOIT = .FALSE. SETZRO = .FALSE. DO 100 I=1,N IF (IFIXX(I,J).NE.0) THEN DOIT = .TRUE. ELSE SETZRO = .TRUE. END IF 100 CONTINUE END IF IF (.NOT.DOIT) THEN DO 110 L=1,NQ CALL SZERO(N,1,FJACD(1,J,L),N) 110 CONTINUE ELSE DO 120 I=1,N IF (XPLUSD(I,J).EQ.ZERO) THEN IF (TT(1,1).LT.ZERO) THEN TYPJ = ONE/ABS(TT(1,1)) ELSE IF (LDTT.EQ.1) THEN TYPJ = ONE/TT(1,J) ELSE TYPJ = ONE/TT(I,J) END IF ELSE TYPJ = ABS(XPLUSD(I,J)) END IF STP(I) = XPLUSD(I,J) + + SIGN(ONE,XPLUSD(I,J)) + *TYPJ*SHSTEP(1,NETA,I,J,STPD,LDSTPD) STP(I) = STP(I) - XPLUSD(I,J) XPLUSD(I,J) = XPLUSD(I,J) + STP(I) 120 CONTINUE ISTOP = 0 CALL FCN(N,M,NP,NQ, + N,M,NP, + BETA,XPLUSD, + IFIXB,IFIXX,LDIFX, + 001,WRK2,WRK6,WRK1, + ISTOP) IF (ISTOP.NE.0) THEN RETURN ELSE NFEV = NFEV + 1 DO 140 L=1,NQ DO 130 I=1,N FJACD(I,J,L) = WRK2(I,L) 130 CONTINUE 140 CONTINUE END IF DO 150 I=1,N XPLUSD(I,J) = X(I,J) + DELTA(I,J) - STP(I) 150 CONTINUE ISTOP = 0 CALL FCN(N,M,NP,NQ, + N,M,NP, + BETA,XPLUSD, + IFIXB,IFIXX,LDIFX, + 001,WRK2,WRK6,WRK1, + ISTOP) IF (ISTOP.NE.0) THEN RETURN ELSE NFEV = NFEV + 1 END IF IF (SETZRO) THEN DO 180 I=1,N IF (IFIXX(I,J).EQ.0) THEN DO 160 L=1,NQ FJACD(I,J,L) = ZERO 160 CONTINUE ELSE DO 170 L=1,NQ FJACD(I,J,L) = (FJACD(I,J,L)-WRK2(I,L))/ + (2*STP(I)) 170 CONTINUE END IF 180 CONTINUE ELSE DO 200 L=1,NQ DO 190 I=1,N FJACD(I,J,L) = (FJACD(I,J,L)-WRK2(I,L))/ + (2*STP(I)) 190 CONTINUE 200 CONTINUE END IF DO 210 I=1,N XPLUSD(I,J) = X(I,J) + DELTA(I,J) 210 CONTINUE END IF 220 CONTINUE END IF RETURN END *SJACFD SUBROUTINE SJACFD + (FCN, + N,M,NP,NQ, + BETA,X,LDX,DELTA,XPLUSD,IFIXB,IFIXX,LDIFX, + STPB,STPD,LDSTPD, + SSF,TT,LDTT,NETA,FN,STP,WRK1,WRK2,WRK3,WRK6, + FJACB,ISODR,FJACD,NFEV,ISTOP) C***BEGIN PROLOGUE SJACFD C***REFER TO SODR,SODRC C***ROUTINES CALLED FCN,SHSTEP,SZERO C***DATE WRITTEN 860529 (YYMMDD) C***REVISION DATE 920619 (YYMMDD) C***PURPOSE COMPUTE FORWARD DIFFERENCE APPROXIMATIONS TO THE C JACOBIAN WRT THE ESTIMATED BETAS AND WRT THE DELTAS C***END PROLOGUE SJACFD C...SCALAR ARGUMENTS INTEGER + ISTOP,LDIFX,LDSTPD,LDTT,LDX,M,N,NETA,NFEV,NP,NQ LOGICAL + ISODR C...ARRAY ARGUMENTS REAL + BETA(NP),DELTA(N,M),FJACB(N,NP,NQ),FJACD(N,M,NQ),FN(N,NQ), + SSF(NP),STP(N),STPB(NP),STPD(LDSTPD,M),TT(LDTT,M), + WRK1(N,M,NQ),WRK2(N,NQ),WRK3(NP),WRK6(N,NP,NQ), + X(LDX,M),XPLUSD(N,M) INTEGER + IFIXB(NP),IFIXX(LDIFX,M) C...SUBROUTINE ARGUMENTS EXTERNAL + FCN C...LOCAL SCALARS REAL + BETAK,ONE,TYPJ,ZERO INTEGER + I,J,K,L LOGICAL + DOIT,SETZRO C...EXTERNAL SUBROUTINES EXTERNAL + SZERO C...EXTERNAL FUNCTIONS REAL + SHSTEP EXTERNAL + SHSTEP C...INTRINSIC FUNCTIONS INTRINSIC + ABS,MAX,SIGN,SQRT C...DATA STATEMENTS DATA + ZERO,ONE + /0.0E0,1.0E0/ C...ROUTINE NAMES USED AS SUBPROGRAM ARGUMENTS C FCN: THE USER SUPPLIED SUBROUTINE FOR EVALUATING THE MODEL. C...VARIABLE DEFINITIONS (ALPHABETICALLY) C BETA: THE FUNCTION PARAMETERS. C BETAK: THE K-TH FUNCTION PARAMETER. C DELTA: THE ESTIMATED ERRORS IN THE EXPLANATORY VARIABLES. C DOIT: THE VARIABLE DESIGNATING WHETHER THE DERIVATIVE WRT A C GIVEN BETA OR DELTA NEEDS TO BE COMPUTED (DOIT=TRUE) C OR NOT (DOIT=FALSE). C FJACB: THE JACOBIAN WITH RESPECT TO BETA. C FJACD: THE JACOBIAN WITH RESPECT TO DELTA. C FN: THE NEW PREDICTED VALUES FROM THE FUNCTION. C I: AN INDEXING VARIABLE. C IFIXB: THE VALUES DESIGNATING WHETHER THE ELEMENTS OF BETA ARE C FIXED AT THEIR INPUT VALUES OR NOT. C IFIXX: THE VALUES DESIGNATING WHETHER THE ELEMENTS OF X ARE C FIXED AT THEIR INPUT VALUES OR NOT. C ISODR: THE VARIABLE DESIGNATING WHETHER THE SOLUTION IS BY ODR C (ISODR=TRUE) OR BY OLS (ISODR=FALSE). C ISTOP: THE VARIABLE DESIGNATING WHETHER THERE ARE PROBLEMS C COMPUTING THE FUNCTION AT THE CURRENT BETA AND DELTA. C J: AN INDEXING VARIABLE. C K: AN INDEXING VARIABLE. C L: AN INDEXING VARIABLE. C LDIFX: THE LEADING DIMENSION OF ARRAY IFIXX. C LDSTPD: THE LEADING DIMENSION OF ARRAY STPD. C LDTT: THE LEADING DIMENSION OF ARRAY TT. C LDX: THE LEADING DIMENSION OF ARRAY X. C M: THE NUMBER OF COLUMNS OF DATA IN THE EXPLANATORY VARIABLE. C N: THE NUMBER OF OBSERVATIONS. C NETA: THE NUMBER OF GOOD DIGITS IN THE FUNCTION RESULTS. C NFEV: THE NUMBER OF FUNCTION EVALUATIONS. C NP: THE NUMBER OF FUNCTION PARAMETERS. C ONE: THE VALUE 1.0E0. C SETZRO: THE VARIABLE DESIGNATING WHETHER THE DERIVATIVE WRT SOME C DELTA NEEDS TO BE SET TO ZERO (SETZRO=TRUE) OR NOT C (SETZRO=FALSE). C SSF: THE SCALE USED FOR THE BETA'S. C STP: THE STEP USED FOR COMPUTING FINITE DIFFERENCE C DERIVATIVES WITH RESPECT TO DELTA. C STPB: THE RELATIVE STEP USED FOR COMPUTING FINITE DIFFERENCE C DERIVATIVES WITH RESPECT TO BETA. C STPD: THE RELATIVE STEP USED FOR COMPUTING FINITE DIFFERENCE C DERIVATIVES WITH RESPECT TO DELTA. C TT: THE SCALING VALUES USED FOR DELTA. C TYPJ: THE TYPICAL SIZE OF THE J-TH UNKNOWN BETA OR DELTA. C X: THE EXPLANATORY VARIABLE. C XPLUSD: THE VALUES OF X + DELTA. C WRK1: A WORK ARRAY OF (N BY M BY NQ) ELEMENTS. C WRK2: A WORK ARRAY OF (N BY NQ) ELEMENTS. C WRK3: A WORK ARRAY OF (NP) ELEMENTS. C WRK6: A WORK ARRAY OF (N BY NP BY NQ) ELEMENTS. C ZERO: THE VALUE 0.0E0. C***FIRST EXECUTABLE STATEMENT SJACFD C COMPUTE THE JACOBIAN WRT THE ESTIMATED BETAS DO 40 K=1,NP IF (IFIXB(1).GE.0) THEN IF (IFIXB(K).EQ.0) THEN DOIT = .FALSE. ELSE DOIT = .TRUE. END IF ELSE DOIT = .TRUE. END IF IF (.NOT.DOIT) THEN DO 10 L=1,NQ CALL SZERO(N,1,FJACB(1,K,L),N) 10 CONTINUE ELSE BETAK = BETA(K) IF (BETAK.EQ.ZERO) THEN IF (SSF(1).LT.ZERO) THEN TYPJ = ONE/ABS(SSF(1)) ELSE TYPJ = ONE/SSF(K) END IF ELSE TYPJ = ABS(BETAK) END IF WRK3(K) = BETAK + + SIGN(ONE,BETAK)*TYPJ*SHSTEP(0,NETA,1,K,STPB,1) WRK3(K) = WRK3(K) - BETAK BETA(K) = BETAK + WRK3(K) ISTOP = 0 CALL FCN(N,M,NP,NQ, + N,M,NP, + BETA,XPLUSD, + IFIXB,IFIXX,LDIFX, + 001,WRK2,WRK6,WRK1, + ISTOP) IF (ISTOP.NE.0) THEN RETURN ELSE NFEV = NFEV + 1 END IF DO 30 L=1,NQ DO 20 I=1,N FJACB(I,K,L) = (WRK2(I,L)-FN(I,L))/WRK3(K) 20 CONTINUE 30 CONTINUE BETA(K) = BETAK END IF 40 CONTINUE C COMPUTE THE JACOBIAN WRT THE X'S IF (ISODR) THEN DO 220 J=1,M IF (IFIXX(1,1).LT.0) THEN DOIT = .TRUE. SETZRO = .FALSE. ELSE IF (LDIFX.EQ.1) THEN IF (IFIXX(1,J).EQ.0) THEN DOIT = .FALSE. ELSE DOIT = .TRUE. END IF SETZRO = .FALSE. ELSE DOIT = .FALSE. SETZRO = .FALSE. DO 100 I=1,N IF (IFIXX(I,J).NE.0) THEN DOIT = .TRUE. ELSE SETZRO = .TRUE. END IF 100 CONTINUE END IF IF (.NOT.DOIT) THEN DO 110 L=1,NQ CALL SZERO(N,1,FJACD(1,J,L),N) 110 CONTINUE ELSE DO 120 I=1,N IF (XPLUSD(I,J).EQ.ZERO) THEN IF (TT(1,1).LT.ZERO) THEN TYPJ = ONE/ABS(TT(1,1)) ELSE IF (LDTT.EQ.1) THEN TYPJ = ONE/TT(1,J) ELSE TYPJ = ONE/TT(I,J) END IF ELSE TYPJ = ABS(XPLUSD(I,J)) END IF STP(I) = XPLUSD(I,J) + + SIGN(ONE,XPLUSD(I,J)) + *TYPJ*SHSTEP(0,NETA,I,J,STPD,LDSTPD) STP(I) = STP(I) - XPLUSD(I,J) XPLUSD(I,J) = XPLUSD(I,J) + STP(I) 120 CONTINUE ISTOP = 0 CALL FCN(N,M,NP,NQ, + N,M,NP, + BETA,XPLUSD, + IFIXB,IFIXX,LDIFX, + 001,WRK2,WRK6,WRK1, + ISTOP) IF (ISTOP.NE.0) THEN RETURN ELSE NFEV = NFEV + 1 DO 140 L=1,NQ DO 130 I=1,N FJACD(I,J,L) = WRK2(I,L) 130 CONTINUE 140 CONTINUE END IF IF (SETZRO) THEN DO 180 I=1,N IF (IFIXX(I,J).EQ.0) THEN DO 160 L=1,NQ FJACD(I,J,L) = ZERO 160 CONTINUE ELSE DO 170 L=1,NQ FJACD(I,J,L) = (FJACD(I,J,L)-FN(I,L))/STP(I) 170 CONTINUE END IF 180 CONTINUE ELSE DO 200 L=1,NQ DO 190 I=1,N FJACD(I,J,L) = (FJACD(I,J,L)-FN(I,L))/STP(I) 190 CONTINUE 200 CONTINUE END IF DO 210 I=1,N XPLUSD(I,J) = X(I,J) + DELTA(I,J) 210 CONTINUE END IF 220 CONTINUE END IF RETURN END *SJCK SUBROUTINE SJCK + (FCN, + N,M,NP,NQ, + BETA,XPLUSD, + IFIXB,IFIXX,LDIFX,STPB,STPD,LDSTPD, + SSF,TT,LDTT, + ETA,NETA,NTOL,NROW,ISODR,EPSMAC, + PV0,FJACB,FJACD, + MSGB,MSGD,DIFF,ISTOP,NFEV,NJEV, + WRK1,WRK2,WRK6) C***BEGIN PROLOGUE SJCK C***REFER TO SODR,SODRC C***ROUTINES CALLED FCN,SHSTEP,SJCKM C***DATE WRITTEN 860529 (YYMMDD) C***REVISION DATE 920619 (YYMMDD) C***PURPOSE DRIVER ROUTINE FOR THE DERIVATIVE CHECKING PROCESS C (ADAPTED FROM STARPAC SUBROUTINE DCKCNT) C***END PROLOGUE SJCK C...SCALAR ARGUMENTS REAL + EPSMAC,ETA INTEGER + ISTOP,LDIFX,LDSTPD,LDTT, + M,N,NETA,NFEV,NJEV,NP,NQ,NROW,NTOL LOGICAL + ISODR C...ARRAY ARGUMENTS REAL + BETA(NP),DIFF(NQ,NP+M),FJACB(N,NP,NQ),FJACD(N,M,NQ), + PV0(N,NQ),SSF(NP),STPB(NP),STPD(LDSTPD,M),TT(LDTT,M), + WRK1(N,M,NQ),WRK2(N,NQ),WRK6(N,NP,NQ),XPLUSD(N,M) INTEGER + IFIXB(NP),IFIXX(LDIFX,M),MSGB(1+NQ*NP),MSGD(1+NQ*M) C...SUBROUTINE ARGUMENTS EXTERNAL + FCN C...LOCAL SCALARS REAL + DIFFJ,H0,HC0,ONE,P5,PV,TOL,TYPJ,ZERO INTEGER + IDEVAL,J,LQ,MSGB1,MSGD1 LOGICAL + ISFIXD,ISWRTB C...EXTERNAL SUBROUTINES EXTERNAL + SJCKM C...EXTERNAL FUNCTIONS REAL + SHSTEP EXTERNAL + SHSTEP C...INTRINSIC FUNCTIONS INTRINSIC + ABS,INT,LOG10 C...DATA STATEMENTS DATA + ZERO,P5,ONE + /0.0E0,0.5E0,1.0E0/ C...ROUTINE NAMES USED AS SUBPROGRAM ARGUMENTS C FCN: THE USER SUPPLIED SUBROUTINE FOR EVALUATING THE MODEL. C...VARIABLE DEFINITIONS (ALPHABETICALLY) C BETA: THE FUNCTION PARAMETERS. C DIFF: THE RELATIVE DIFFERENCES BETWEEN THE USER SUPPLIED AND C FINITE DIFFERENCE DERIVATIVES FOR EACH DERIVATIVE CHECKED. C DIFFJ: THE RELATIVE DIFFERENCES BETWEEN THE USER SUPPLIED AND C FINITE DIFFERENCE DERIVATIVES FOR THE DERIVATIVE BEING C CHECKED. C EPSMAC: THE VALUE OF MACHINE PRECISION. C ETA: THE RELATIVE NOISE IN THE FUNCTION RESULTS. C FJACB: THE JACOBIAN WITH RESPECT TO BETA. C FJACD: THE JACOBIAN WITH RESPECT TO DELTA. C H0: THE INITIAL RELATIVE STEP SIZE FOR FORWARD DIFFERENCES. C HC0: THE INITIAL RELATIVE STEP SIZE FOR CENTRAL DIFFERENCES. C IDEVAL: THE VARIABLE DESIGNATING WHAT COMPUTATIONS ARE TO BE C PERFORMED BY USER SUPPLIED SUBROUTINE FCN. C IFIXB: THE VALUES DESIGNATING WHETHER THE ELEMENTS OF BETA ARE C FIXED AT THEIR INPUT VALUES OR NOT. C IFIXX: THE VALUES DESIGNATING WHETHER THE ELEMENTS OF X ARE C FIXED AT THEIR INPUT VALUES OR NOT. C ISFIXD: THE VARIABLE DESIGNATING WHETHER THE PARAMETER IS FIXED C (ISFIXD=TRUE) OR NOT (ISFIXD=FALSE). C ISTOP: THE VARIABLE DESIGNATING WHETHER THERE ARE PROBLEMS C COMPUTING THE FUNCTION AT THE CURRENT BETA AND DELTA. C ISODR: THE VARIABLE DESIGNATING WHETHER THE SOLUTION IS BY ODR C (ISODR=.TRUE.) OR BY OLS (ISODR=.FALSE.). C ISWRTB: THE VARIABLE DESIGNATING WHETHER THE DERIVATIVES WRT BETA C (ISWRTB=TRUE) OR DELTA (ISWRTB=FALSE) ARE BEING CHECKED. C J: AN INDEX VARIABLE. C LDIFX: THE LEADING DIMENSION OF ARRAY IFIXX. C LDSTPD: THE LEADING DIMENSION OF ARRAY STPD. C LDTT: THE LEADING DIMENSION OF ARRAY TT. C LQ: THE RESPONSE CURRENTLY BEING EXAMINED. C M: THE NUMBER OF COLUMNS OF DATA IN THE EXPLANATORY VARIABLE. C MSGB: THE ERROR CHECKING RESULTS FOR THE JACOBIAN WRT BETA. C MSGB1: THE ERROR CHECKING RESULTS FOR THE JACOBIAN WRT BETA. C MSGD: THE ERROR CHECKING RESULTS FOR THE JACOBIAN WRT DELTA. C MSGD1: THE ERROR CHECKING RESULTS FOR THE JACOBIAN WRT DELTA. C N: THE NUMBER OF OBSERVATIONS. C NETA: THE NUMBER OF RELIABLE DIGITS IN THE MODEL RESULTS, EITHER C SET BY THE USER OR COMPUTED BY SETAF. C NFEV: THE NUMBER OF FUNCTION EVALUATIONS. C NJEV: THE NUMBER OF JACOBIAN EVALUATIONS. C NP: THE NUMBER OF FUNCTION PARAMETERS. C NQ: THE NUMBER OF RESPONSES PER OBSERVATION. C NROW: THE ROW NUMBER OF THE EXPLANATORY VARIABLE ARRAY AT WHICH C THE DERIVATIVE IS CHECKED. C NTOL: THE NUMBER OF DIGITS OF AGREEMENT REQUIRED BETWEEN THE C NUMERICAL DERIVATIVES AND THE USER SUPPLIED DERIVATIVES. C ONE: THE VALUE 1.0E0. C P5: THE VALUE 0.5E0. C PV: THE SCALAR IN WHICH THE PREDICTED VALUE FROM THE MODEL FOR C ROW NROW IS STORED. C PV0: THE PREDICTED VALUES USING THE CURRENT PARAMETER ESTIMATES. C SSF: THE SCALING VALUES USED FOR BETA. C STPB: THE STEP SIZE FOR FINITE DIFFERENCE DERIVATIVES WRT BETA. C STPD: THE STEP SIZE FOR FINITE DIFFERENCE DERIVATIVES WRT DELTA. C TOL: THE AGREEMENT TOLERANCE. C TT: THE SCALING VALUES USED FOR DELTA. C TYPJ: THE TYPICAL SIZE OF THE J-TH UNKNOWN BETA OR DELTA. C WRK1: A WORK ARRAY OF (N BY M BY NQ) ELEMENTS. C WRK2: A WORK ARRAY OF (N BY NQ) ELEMENTS. C WRK6: A WORK ARRAY OF (N BY NP BY NQ) ELEMENTS. C XPLUSD: THE VALUES OF X + DELTA. C ZERO: THE VALUE 0.0E0. C***FIRST EXECUTABLE STATEMENT SJCK C SET TOLERANCE FOR CHECKING DERIVATIVES TOL = ETA**(0.25E0) NTOL = MAX(ONE,P5-LOG10(TOL)) C COMPUTE USER SUPPLIED DERIVATIVE VALUES ISTOP = 0 IF (ISODR) THEN IDEVAL = 110 ELSE IDEVAL = 010 END IF CALL FCN(N,M,NP,NQ, + N,M,NP, + BETA,XPLUSD, + IFIXB,IFIXX,LDIFX, + IDEVAL,WRK2,FJACB,FJACD, + ISTOP) IF (ISTOP.NE.0) THEN RETURN ELSE NJEV = NJEV + 1 END IF C CHECK DERIVATIVES WRT BETA FOR EACH RESPONSE OF OBSERVATION NROW MSGB1 = 0 MSGD1 = 0 DO 30 LQ=1,NQ C SET PREDICTED VALUE OF MODEL AT CURRENT PARAMETER ESTIMATES PV = PV0(NROW,LQ) ISWRTB = .TRUE. DO 10 J=1,NP IF (IFIXB(1).LT.0) THEN ISFIXD = .FALSE. ELSE IF (IFIXB(J).EQ.0) THEN ISFIXD = .TRUE. ELSE ISFIXD = .FALSE. END IF IF (ISFIXD) THEN MSGB(1+LQ+(J-1)*NQ) = -1 ELSE IF (BETA(J).EQ.ZERO) THEN IF (SSF(1).LT.ZERO) THEN TYPJ = ONE/ABS(SSF(1)) ELSE TYPJ = ONE/SSF(J) END IF ELSE TYPJ = ABS(BETA(J)) END IF H0 = SHSTEP(0,NETA,1,J,STPB,1) HC0 = H0 C CHECK DERIVATIVE WRT THE J-TH PARAMETER AT THE NROW-TH ROW CALL SJCKM(FCN, + N,M,NP,NQ, + BETA,XPLUSD, + IFIXB,IFIXX,LDIFX, + ETA,TOL,NROW,EPSMAC,J,LQ,TYPJ,H0,HC0, + ISWRTB,PV,FJACB(NROW,J,LQ), + DIFFJ,MSGB1,MSGB(2),ISTOP,NFEV, + WRK1,WRK2,WRK6) IF (ISTOP.NE.0) THEN MSGB(1) = -1 RETURN ELSE DIFF(LQ,J) = DIFFJ END IF END IF 10 CONTINUE C CHECK DERIVATIVES WRT X FOR EACH RESPONSE OF OBSERVATION NROW IF (ISODR) THEN ISWRTB = .FALSE. DO 20 J=1,M IF (IFIXX(1,1).LT.0) THEN ISFIXD = .FALSE. ELSE IF (LDIFX.EQ.1) THEN IF (IFIXX(1,J).EQ.0) THEN ISFIXD = .TRUE. ELSE ISFIXD = .FALSE. END IF ELSE ISFIXD = .FALSE. END IF IF (ISFIXD) THEN MSGD(1+LQ+(J-1)*NQ) = -1 ELSE IF (XPLUSD(NROW,J).EQ.ZERO) THEN IF (TT(1,1).LT.ZERO) THEN TYPJ = ONE/ABS(TT(1,1)) ELSE IF (LDTT.EQ.1) THEN TYPJ = ONE/TT(1,J) ELSE TYPJ = ONE/TT(NROW,J) END IF ELSE TYPJ = ABS(XPLUSD(NROW,J)) END IF H0 = SHSTEP(0,NETA,NROW,J,STPD,LDSTPD) HC0 = SHSTEP(1,NETA,NROW,J,STPD,LDSTPD) C CHECK DERIVATIVE WRT THE J-TH COLUMN OF DELTA AT ROW NROW CALL SJCKM(FCN, + N,M,NP,NQ, + BETA,XPLUSD, + IFIXB,IFIXX,LDIFX, + ETA,TOL,NROW,EPSMAC,J,LQ,TYPJ,H0,HC0, + ISWRTB,PV,FJACD(NROW,J,LQ), + DIFFJ,MSGD1,MSGD(2),ISTOP,NFEV, + WRK1,WRK2,WRK6) IF (ISTOP.NE.0) THEN MSGD(1) = -1 RETURN ELSE DIFF(LQ,NP+J) = DIFFJ END IF END IF 20 CONTINUE END IF 30 CONTINUE MSGB(1) = MSGB1 MSGD(1) = MSGD1 RETURN END *SJCKC SUBROUTINE SJCKC + (FCN, + N,M,NP,NQ, + BETA,XPLUSD,IFIXB,IFIXX,LDIFX, + ETA,TOL,NROW,EPSMAC,J,LQ,HC,ISWRTB, + FD,TYPJ,PVPSTP,STP0, + PV,D, + DIFFJ,MSG,ISTOP,NFEV, + WRK1,WRK2,WRK6) C***BEGIN PROLOGUE SJCKC C***REFER TO SODR,SODRC C***ROUTINES CALLED SJCKF,SPVB,SPVD C***DATE WRITTEN 860529 (YYMMDD) C***REVISION DATE 920619 (YYMMDD) C***PURPOSE CHECK WHETHER HIGH CURVATURE COULD BE THE CAUSE OF THE C DISAGREEMENT BETWEEN THE NUMERICAL AND ANALYTIC DERVIATIVES C (ADAPTED FROM STARPAC SUBROUTINE DCKCRV) C***END PROLOGUE SJCKC C...SCALAR ARGUMENTS REAL + D,DIFFJ,EPSMAC,ETA,FD,HC,PV,PVPSTP,STP0,TOL,TYPJ INTEGER + ISTOP,J,LDIFX,LQ,M,N,NFEV,NP,NQ,NROW LOGICAL + ISWRTB C...ARRAY ARGUMENTS REAL + BETA(NP),WRK1(N,M,NQ),WRK2(N,NQ),WRK6(N,NP,NQ),XPLUSD(N,M) INTEGER + IFIXB(NP),IFIXX(LDIFX,M),MSG(NQ,J) C...SUBROUTINE ARGUMENTS EXTERNAL + FCN C...LOCAL SCALARS REAL + CURVE,ONE,PVMCRV,PVPCRV,P01,STP,STPCRV,TEN,TWO C...EXTERNAL SUBROUTINES EXTERNAL + SJCKF,SPVB,SPVD C...INTRINSIC FUNCTIONS INTRINSIC + ABS,SIGN C...DATA STATEMENTS DATA + P01,ONE,TWO,TEN + /0.01E0,1.0E0,2.0E0,10.0E0/ C...ROUTINE NAMES USED AS SUBPROGRAM ARGUMENTS C FCN: THE USER SUPPLIED SUBROUTINE FOR EVALUATING THE MODEL. C...VARIABLE DEFINITIONS (ALPHABETICALLY) C BETA: THE FUNCTION PARAMETERS. C CURVE: A MEASURE OF THE CURVATURE IN THE MODEL. C D: THE DERIVATIVE WITH RESPECT TO THE JTH UNKNOWN PARAMETER. C DIFFJ: THE RELATIVE DIFFERENCES BETWEEN THE USER SUPPLIED AND C FINITE DIFFERENCE DERIVATIVES FOR THE DERIVATIVE BEING C CHECKED. C EPSMAC: THE VALUE OF MACHINE PRECISION. C ETA: THE RELATIVE NOISE IN THE MODEL C FD: THE FORWARD DIFFERENCE DERIVATIVE WRT THE JTH PARAMETER. C HC: THE RELATIVE STEP SIZE FOR CENTRAL FINITE DIFFERENCES. C IFIXB: THE VALUES DESIGNATING WHETHER THE ELEMENTS OF BETA ARE C FIXED AT THEIR INPUT VALUES OR NOT. C IFIXX: THE VALUES DESIGNATING WHETHER THE ELEMENTS OF X ARE C FIXED AT THEIR INPUT VALUES OR NOT. C ISTOP: THE VARIABLE DESIGNATING WHETHER THERE ARE PROBLEMS C COMPUTING THE FUNCTION AT THE CURRENT BETA AND DELTA. C ISWRTB: THE VARIABLE DESIGNATING WHETHER THE DERIVATIVES WRT BETA C (ISWRTB=TRUE) OR DELTA(ISWRTB=FALSE) ARE BEING CHECKED. C J: THE INDEX OF THE PARTIAL DERIVATIVE BEING EXAMINED. C LDIFX: THE LEADING DIMENSION OF ARRAY IFIXX. C LQ: THE RESPONSE CURRENTLY BEING EXAMINED. C M: THE NUMBER OF COLUMNS OF DATA IN THE EXPLANATORY VARIABLE. C MSG: THE ERROR CHECKING RESULTS. C N: THE NUMBER OF OBSERVATIONS. C NFEV: THE NUMBER OF FUNCTION EVALUATIONS. C NP: THE NUMBER OF FUNCTION PARAMETERS. C NQ: THE NUMBER OF RESPONSES PER OBSERVATION. C NROW: THE ROW NUMBER OF THE EXPLANATORY VARIABLE ARRAY AT WHICH C THE DERIVATIVE IS TO BE CHECKED. C ONE: THE VALUE 1.0E0. C PV: THE PREDICTED VALUE OF THE MODEL FOR ROW NROW . C PVMCRV: THE PREDICTED VALUE FOR ROW NROW OF THE MODEL C BASED ON THE CURRENT PARAMETER ESTIMATES FOR ALL BUT THE C JTH PARAMETER VALUE, WHICH IS BETA(J)-STPCRV. C PVPCRV: THE PREDICTED VALUE FOR ROW NROW OF THE MODEL C BASED ON THE CURRENT PARAMETER ESTIMATES FOR ALL BUT THE C JTH PARAMETER VALUE, WHICH IS BETA(J)+STPCRV. C PVPSTP: THE PREDICTED VALUE FOR ROW NROW OF THE MODEL C BASED ON THE CURRENT PARAMETER ESTIMATES FOR ALL BUT THE C JTH PARAMETER VALUE, WHICH IS BETA(J) + STP0. C P01: THE VALUE 0.01E0. C STP0: THE INITIAL STEP SIZE FOR THE FINITE DIFFERENCE DERIVATIVE. C STP: A STEP SIZE FOR THE FINITE DIFFERENCE DERIVATIVE. C STPCRV: THE STEP SIZE SELECTED TO CHECK FOR CURVATURE IN THE MODEL. C TEN: THE VALUE 10.0E0. C TOL: THE AGREEMENT TOLERANCE. C TWO: THE VALUE 2.0E0. C TYPJ: THE TYPICAL SIZE OF THE J-TH UNKNOWN BETA OR DELTA. C WRK1: A WORK ARRAY OF (N BY M BY NQ) ELEMENTS. C WRK2: A WORK ARRAY OF (N BY NQ) ELEMENTS. C WRK6: A WORK ARRAY OF (N BY NP BY NQ) ELEMENTS. C XPLUSD: THE VALUES OF X + DELTA. C***FIRST EXECUTABLE STATEMENT SJCKC IF (ISWRTB) THEN C PERFORM CENTRAL DIFFERENCE COMPUTATIONS FOR DERIVATIVES WRT BETA STPCRV = (HC*TYPJ*SIGN(ONE,BETA(J))+BETA(J)) - BETA(J) CALL SPVB(FCN, + N,M,NP,NQ, + BETA,XPLUSD,IFIXB,IFIXX,LDIFX, + NROW,J,LQ,STPCRV, + ISTOP,NFEV,PVPCRV, + WRK1,WRK2,WRK6) IF (ISTOP.NE.0) THEN RETURN END IF CALL SPVB(FCN, + N,M,NP,NQ, + BETA,XPLUSD,IFIXB,IFIXX,LDIFX, + NROW,J,LQ,-STPCRV, + ISTOP,NFEV,PVMCRV, + WRK1,WRK2,WRK6) IF (ISTOP.NE.0) THEN RETURN END IF ELSE C PERFORM CENTRAL DIFFERENCE COMPUTATIONS FOR DERIVATIVES WRT DELTA STPCRV = (HC*TYPJ*SIGN(ONE,XPLUSD(NROW,J))+XPLUSD(NROW,J)) - + XPLUSD(NROW,J) CALL SPVD(FCN, + N,M,NP,NQ, + BETA,XPLUSD,IFIXB,IFIXX,LDIFX, + NROW,J,LQ,STPCRV, + ISTOP,NFEV,PVPCRV, + WRK1,WRK2,WRK6) IF (ISTOP.NE.0) THEN RETURN END IF CALL SPVD(FCN, + N,M,NP,NQ, + BETA,XPLUSD,IFIXB,IFIXX,LDIFX, + NROW,J,LQ,-STPCRV, + ISTOP,NFEV,PVMCRV, + WRK1,WRK2,WRK6) IF (ISTOP.NE.0) THEN RETURN END IF END IF C ESTIMATE CURVATURE BY SECOND DERIVATIVE OF MODEL CURVE = ABS((PVPCRV-PV)+(PVMCRV-PV)) / (STPCRV*STPCRV) CURVE = CURVE + + ETA*(ABS(PVPCRV)+ABS(PVMCRV)+TWO*ABS(PV)) / (STPCRV**2) C CHECK IF FINITE PRECISION ARITHMETIC COULD BE THE CULPRIT. CALL SJCKF(FCN, + N,M,NP,NQ, + BETA,XPLUSD,IFIXB,IFIXX,LDIFX, + ETA,TOL,NROW,J,LQ,ISWRTB, + FD,TYPJ,PVPSTP,STP0,CURVE,PV,D, + DIFFJ,MSG,ISTOP,NFEV, + WRK1,WRK2,WRK6) IF (ISTOP.NE.0) THEN RETURN END IF IF (MSG(LQ,J).EQ.0) THEN RETURN END IF C CHECK IF HIGH CURVATURE COULD BE THE PROBLEM. STP = TWO*MAX(TOL*ABS(D)/CURVE,EPSMAC) IF (STP.LT.ABS(TEN*STP0)) THEN STP = MIN(STP,P01*ABS(STP0)) END IF IF (ISWRTB) THEN C PERFORM COMPUTATIONS FOR DERIVATIVES WRT BETA STP = (STP*SIGN(ONE,BETA(J)) + BETA(J)) - BETA(J) CALL SPVB(FCN, + N,M,NP,NQ, + BETA,XPLUSD,IFIXB,IFIXX,LDIFX, + NROW,J,LQ,STP, + ISTOP,NFEV,PVPSTP, + WRK1,WRK2,WRK6) IF (ISTOP.NE.0) THEN RETURN END IF ELSE C PERFORM COMPUTATIONS FOR DERIVATIVES WRT DELTA STP = (STP*SIGN(ONE,XPLUSD(NROW,J)) + XPLUSD(NROW,J)) - + XPLUSD(NROW,J) CALL SPVD(FCN, + N,M,NP,NQ, + BETA,XPLUSD,IFIXB,IFIXX,LDIFX, + NROW,J,LQ,STP, + ISTOP,NFEV,PVPSTP, + WRK1,WRK2,WRK6) IF (ISTOP.NE.0) THEN RETURN END IF END IF C COMPUTE THE NEW NUMERICAL DERIVATIVE FD = (PVPSTP-PV)/STP DIFFJ = MIN(DIFFJ,ABS(FD-D)/ABS(D)) C CHECK WHETHER THE NEW NUMERICAL DERIVATIVE IS OK IF (ABS(FD-D).LE.TOL*ABS(D)) THEN MSG(LQ,J) = 0 C CHECK IF FINITE PRECISION MAY BE THE CULPRIT (FUDGE FACTOR = 2) ELSE IF (ABS(STP*(FD-D)).LT.TWO*ETA*(ABS(PV)+ABS(PVPSTP)) + + CURVE*(EPSMAC*TYPJ)**2) THEN MSG(LQ,J) = 5 END IF RETURN END *SJCKF SUBROUTINE SJCKF + (FCN, + N,M,NP,NQ, + BETA,XPLUSD,IFIXB,IFIXX,LDIFX, + ETA,TOL,NROW,J,LQ,ISWRTB, + FD,TYPJ,PVPSTP,STP0,CURVE,PV,D, + DIFFJ,MSG,ISTOP,NFEV, + WRK1,WRK2,WRK6) C***BEGIN PROLOGUE SJCKF C***REFER TO SODR,SODRC C***ROUTINES CALLED SPVB,SPVD C***DATE WRITTEN 860529 (YYMMDD) C***REVISION DATE 920619 (YYMMDD) C***PURPOSE CHECK WHETHER FINITE PRECISION ARITHMETIC COULD BE THE C CAUSE OF THE DISAGREEMENT BETWEEN THE DERIVATIVES C (ADAPTED FROM STARPAC SUBROUTINE DCKFPA) C***END PROLOGUE SJCKF C...SCALAR ARGUMENTS REAL + CURVE,D,DIFFJ,ETA,FD,PV,PVPSTP,STP0,TOL,TYPJ INTEGER + ISTOP,J,LDIFX,LQ,M,N,NFEV,NP,NQ,NROW LOGICAL + ISWRTB C...ARRAY ARGUMENTS REAL + BETA(NP),WRK1(N,M,NQ),WRK2(N,NQ),WRK6(N,NP,NQ),XPLUSD(N,M) INTEGER + IFIXB(NP),IFIXX(LDIFX,M),MSG(NQ,J) C...SUBROUTINE ARGUMENTS EXTERNAL + FCN C...LOCAL SCALARS REAL + HUNDRD,ONE,P1,STP,TWO LOGICAL + LARGE C...EXTERNAL SUBROUTINES EXTERNAL + SPVB,SPVD C...INTRINSIC FUNCTIONS INTRINSIC + ABS,SIGN C...DATA STATEMENTS DATA + P1,ONE,TWO,HUNDRD + /0.1E0,1.0E0,2.0E0,100.0E0/ C...ROUTINE NAMES USED AS SUBPROGRAM ARGUMENTS C FCN: THE USER SUPPLIED SUBROUTINE FOR EVALUATING THE MODEL. C...VARIABLE DEFINITIONS (ALPHABETICALLY) C BETA: THE FUNCTION PARAMETERS. C CURVE: A MEASURE OF THE CURVATURE IN THE MODEL. C D: THE DERIVATIVE WITH RESPECT TO THE JTH UNKNOWN PARAMETER. C DIFFJ: THE RELATIVE DIFFERENCES BETWEEN THE USER SUPPLIED AND C FINITE DIFFERENCE DERIVATIVES FOR THE DERIVATIVE BEING C CHECKED. C ETA: THE RELATIVE NOISE IN THE MODEL C FD: THE FORWARD DIFFERENCE DERIVATIVE WRT THE JTH PARAMETER. C HUNDRD: THE VALUE 100.0E0. C IFIXB: THE VALUES DESIGNATING WHETHER THE ELEMENTS OF BETA ARE C FIXED AT THEIR INPUT VALUES OR NOT. C IFIXX: THE VALUES DESIGNATING WHETHER THE ELEMENTS OF X ARE C FIXED AT THEIR INPUT VALUES OR NOT. C ISTOP: THE VARIABLE DESIGNATING WHETHER THERE ARE PROBLEMS C COMPUTING THE FUNCTION AT THE CURRENT BETA AND DELTA. C ISWRTB: THE VARIABLE DESIGNATING WHETHER THE DERIVATIVES WRT BETA C (ISWRTB=TRUE) OR DELTA(ISWRTB=FALSE) ARE BEING CHECKED. C J: THE INDEX OF THE PARTIAL DERIVATIVE BEING EXAMINED. C LARGE: THE VALUE DESIGNATING WHETHER THE RECOMMENDED INCREASE IN C THE STEP SIZE WOULD BE GREATER THAN TYPJ. C LDIFX: THE LEADING DIMENSION OF ARRAY IFIXX. C LQ: THE RESPONSE CURRENTLY BEING EXAMINED. C M: THE NUMBER OF COLUMNS OF DATA IN THE EXPLANATORY VARIABLE. C MSG: THE ERROR CHECKING RESULTS. C N: THE NUMBER OF OBSERVATIONS. C NFEV: THE NUMBER OF FUNCTION EVALUATIONS. C NP: THE NUMBER OF FUNCTION PARAMETERS. C NQ: THE NUMBER OF RESPONSES PER OBSERVATION. C NROW: THE ROW NUMBER OF THE EXPLANATORY VARIABLE ARRAY AT WHICH C THE DERIVATIVE IS TO BE CHECKED. C ONE: THE VALUE 1.0E0. C PV: THE PREDICTED VALUE FOR ROW NROW . C PVPSTP: THE PREDICTED VALUE FOR ROW NROW OF THE MODEL C BASED ON THE CURRENT PARAMETER ESTIMATES FOR ALL BUT THE C JTH PARAMETER VALUE, WHICH IS BETA(J) + STP0. C P1: THE VALUE 0.1E0. C STP0: THE STEP SIZE FOR THE FINITE DIFFERENCE DERIVATIVE. C TOL: THE AGREEMENT TOLERANCE. C TWO: THE VALUE 2.0E0. C TYPJ: THE TYPICAL SIZE OF THE J-TH UNKNOWN BETA OR DELTA. C WRK1: A WORK ARRAY OF (N BY M BY NQ) ELEMENTS. C WRK2: A WORK ARRAY OF (N BY NQ) ELEMENTS. C WRK6: A WORK ARRAY OF (N BY NP BY NQ) ELEMENTS. C XPLUSD: THE VALUES OF X + DELTA. C***FIRST EXECUTABLE STATEMENT SJCKF C FINITE PRECISION ARITHMETIC COULD BE THE PROBLEM. C TRY A LARGER STEP SIZE BASED ON ESTIMATE OF CONDITION ERROR STP = ETA*(ABS(PV)+ABS(PVPSTP))/(TOL*ABS(D)) IF (STP.GT.ABS(P1*STP0)) THEN STP = MAX(STP,HUNDRD*ABS(STP0)) END IF IF (STP.GT.TYPJ) THEN STP = TYPJ LARGE = .TRUE. ELSE LARGE = .FALSE. END IF IF (ISWRTB) THEN C PERFORM COMPUTATIONS FOR DERIVATIVES WRT BETA STP = (STP*SIGN(ONE,BETA(J))+BETA(J)) - BETA(J) CALL SPVB(FCN, + N,M,NP,NQ, + BETA,XPLUSD,IFIXB,IFIXX,LDIFX, + NROW,J,LQ,STP, + ISTOP,NFEV,PVPSTP, + WRK1,WRK2,WRK6) ELSE C PERFORM COMPUTATIONS FOR DERIVATIVES WRT DELTA STP = (STP*SIGN(ONE,XPLUSD(NROW,J)) + XPLUSD(NROW,J)) - + XPLUSD(NROW,J) CALL SPVD(FCN, + N,M,NP,NQ, + BETA,XPLUSD,IFIXB,IFIXX,LDIFX, + NROW,J,LQ,STP, + ISTOP,NFEV,PVPSTP, + WRK1,WRK2,WRK6) END IF IF (ISTOP.NE.0) THEN RETURN END IF FD = (PVPSTP-PV)/STP DIFFJ = MIN(DIFFJ,ABS(FD-D)/ABS(D)) C CHECK FOR AGREEMENT IF ((ABS(FD-D)).LE.TOL*ABS(D)) THEN C FORWARD DIFFERENCE QUOTIENT AND ANALYTIC DERIVATIVES AGREE. MSG(LQ,J) = 0 ELSE IF ((ABS(FD-D).LE.ABS(TWO*CURVE*STP)) .OR. LARGE) THEN C CURVATURE MAY BE THE CULPRIT (FUDGE FACTOR = 2) IF (LARGE) THEN MSG(LQ,J) = 4 ELSE MSG(LQ,J) = 5 END IF END IF RETURN END *SJCKM SUBROUTINE SJCKM + (FCN, + N,M,NP,NQ, + BETA,XPLUSD,IFIXB,IFIXX,LDIFX, + ETA,TOL,NROW,EPSMAC,J,LQ,TYPJ,H0,HC0, + ISWRTB,PV,D, + DIFFJ,MSG1,MSG,ISTOP,NFEV, + WRK1,WRK2,WRK6) C***BEGIN PROLOGUE SJCKM C***REFER TO SODR,SODRC C***ROUTINES CALLED SJCKC,SJCKZ,SPVB,SPVD C***DATE WRITTEN 860529 (YYMMDD) C***REVISION DATE 920619 (YYMMDD) C***PURPOSE CHECK USER SUPPLIED ANALYTIC DERIVATIVES AGAINST NUMERICAL C DERIVATIVES C (ADAPTED FROM STARPAC SUBROUTINE DCKMN) C***END PROLOGUE SJCKM C...SCALAR ARGUMENTS REAL + D,DIFFJ,EPSMAC,ETA,H0,HC0,PV,TOL,TYPJ INTEGER + ISTOP,J,LDIFX,LQ,M,MSG1,N,NFEV,NP,NQ,NROW LOGICAL + ISWRTB C...ARRAY ARGUMENTS REAL + BETA(NP),WRK1(N,M,NQ),WRK2(N,NQ),WRK6(N,NP,NQ),XPLUSD(N,M) INTEGER + IFIXB(NP),IFIXX(LDIFX,M),MSG(NQ,J) C...SUBROUTINE ARGUMENTS EXTERNAL + FCN C...LOCAL SCALARS REAL + BIG,FD,H,HC,H1,HC1,HUNDRD,ONE,PVPSTP,P01,P1,STP0, + TEN,THREE,TOL2,TWO,ZERO INTEGER + I C...EXTERNAL SUBROUTINES EXTERNAL + SJCKC,SJCKZ,SPVB,SPVD C...INTRINSIC FUNCTIONS INTRINSIC + ABS,MAX,SIGN,SQRT C...DATA STATEMENTS DATA + ZERO,P01,P1,ONE,TWO,THREE,TEN,HUNDRD + /0.0E0,0.01E0,0.1E0,1.0E0,2.0E0,3.0E0,1.0E1,1.0E2/ DATA + BIG,TOL2 + /1.0E19,5.0E-2/ C...ROUTINE NAMES USED AS SUBPROGRAM ARGUMENTS C FCN: THE USER SUPPLIED SUBROUTINE FOR EVALUATING THE MODEL. C...VARIABLE DEFINITIONS (ALPHABETICALLY) C BETA: THE FUNCTION PARAMETERS. C BIG: A BIG VALUE, USED TO INITIALIZE DIFFJ. C D: THE DERIVATIVE WITH RESPECT TO THE JTH UNKNOWN PARAMETER. C DIFFJ: THE RELATIVE DIFFERENCES BETWEEN THE USER SUPPLIED AND C FINITE DIFFERENCE DERIVATIVES FOR THE DERIVATIVE BEING C CHECKED. C EPSMAC: THE VALUE OF MACHINE PRECISION. C ETA: THE RELATIVE NOISE IN THE FUNCTION RESULTS. C FD: THE FORWARD DIFFERENCE DERIVATIVE WRT THE JTH PARAMETER. C H: THE RELATIVE STEP SIZE FOR FORWARD DIFFERENCES. C H0: THE INITIAL RELATIVE STEP SIZE FOR FORWARD DIFFERENCES. C H1: THE DEFAULT RELATIVE STEP SIZE FOR FORWARD DIFFERENCES. C HC: THE RELATIVE STEP SIZE FOR CENTRAL DIFFERENCES. C HC0: THE INITIAL RELATIVE STEP SIZE FOR CENTRAL DIFFERENCES. C HC1: THE DEFAULT RELATIVE STEP SIZE FOR CENTRAL DIFFERENCES. C HUNDRD: THE VALUE 100.0E0. C IFIXB: THE VALUES DESIGNATING WHETHER THE ELEMENTS OF BETA ARE C FIXED AT THEIR INPUT VALUES OR NOT. C IFIXX: THE VALUES DESIGNATING WHETHER THE ELEMENTS OF X ARE C FIXED AT THEIR INPUT VALUES OR NOT. C ISTOP: THE VARIABLE DESIGNATING WHETHER THERE ARE PROBLEMS C COMPUTING THE FUNCTION AT THE CURRENT BETA AND DELTA. C ISWRTB: THE VARIABLE DESIGNATING WHETHER THE DERIVATIVES WRT BETA C (ISWRTB=TRUE) OR DELTAS (ISWRTB=FALSE) ARE BEING CHECKED. C J: THE INDEX OF THE PARTIAL DERIVATIVE BEING EXAMINED. C LDIFX: THE LEADING DIMENSION OF ARRAY IFIXX. C LQ: THE RESPONSE CURRENTLY BEING EXAMINED. C M: THE NUMBER OF COLUMNS OF DATA IN THE EXPLANATORY VARIABLE. C MSG: THE ERROR CHECKING RESULTS. C MSG1: THE ERROR CHECKING RESULTS SUMMARY. C N: THE NUMBER OF OBSERVATIONS. C NFEV: THE NUMBER OF FUNCTION EVALUATIONS. C NP: THE NUMBER OF FUNCTION PARAMETERS. C NQ: THE NUMBER OF RESPONSES PER OBSERVATION. C NROW: THE ROW NUMBER OF THE EXPLANATORY VARIABLE ARRAY AT WHICH C THE DERIVATIVE IS TO BE CHECKED. C ONE: THE VALUE 1.0E0. C PV: THE PREDICTED VALUE FROM THE MODEL FOR ROW NROW . C PVPSTP: THE PREDICTED VALUE FOR ROW NROW OF THE MODEL C USING THE CURRENT PARAMETER ESTIMATES FOR ALL BUT THE JTH C PARAMETER VALUE, WHICH IS BETA(J) + STP0. C P01: THE VALUE 0.01E0. C P1: THE VALUE 0.1E0. C STP0: THE INITIAL STEP SIZE FOR THE FINITE DIFFERENCE DERIVATIVE. C TEN: THE VALUE 10.0E0. C THREE: THE VALUE 3.0E0. C TWO: THE VALUE 2.0E0. C TOL: THE AGREEMENT TOLERANCE. C TOL2: A MINIMUM AGREEMENT TOLERANCE. C TYPJ: THE TYPICAL SIZE OF THE J-TH UNKNOWN BETA OR DELTA. C WRK1: A WORK ARRAY OF (N BY M BY NQ) ELEMENTS. C WRK2: A WORK ARRAY OF (N BY NQ) ELEMENTS. C WRK6: A WORK ARRAY OF (N BY NP BY NQ) ELEMENTS. C XPLUSD: THE VALUES OF X + DELTA. C ZERO: THE VALUE 0.0E0. C***FIRST EXECUTABLE STATEMENT SJCKM C CALCULATE THE JTH PARTIAL DERIVATIVE USING FORWARD DIFFERENCE C QUOTIENTS AND DECIDE IF IT AGREES WITH USER SUPPLIED VALUES H1 = SQRT(ETA) HC1 = ETA**(ONE/THREE) MSG(LQ,J) = 7 DIFFJ = BIG DO 10 I=1,3 IF (I.EQ.1) THEN C TRY INITIAL RELATIVE STEP SIZE H = H0 HC = HC0 ELSE IF (I.EQ.2) THEN C TRY LARGER RELATIVE STEP SIZE H = MAX(TEN*H1, MIN(HUNDRD*H0, ONE)) HC = MAX(TEN*HC1,MIN(HUNDRD*HC0,ONE)) ELSE IF (I.EQ.3) THEN C TRY SMALLER RELATIVE STEP SIZE H = MIN(P1*H1, MAX(P01*H,TWO*EPSMAC)) HC = MIN(P1*HC1,MAX(P01*HC,TWO*EPSMAC)) END IF IF (ISWRTB) THEN C PERFORM COMPUTATIONS FOR DERIVATIVES WRT BETA STP0 = (H*TYPJ*SIGN(ONE,BETA(J))+BETA(J)) - BETA(J) CALL SPVB(FCN, + N,M,NP,NQ, + BETA,XPLUSD,IFIXB,IFIXX,LDIFX, + NROW,J,LQ,STP0, + ISTOP,NFEV,PVPSTP, + WRK1,WRK2,WRK6) ELSE C PERFORM COMPUTATIONS FOR DERIVATIVES WRT DELTA STP0 = (H*TYPJ*SIGN(ONE,XPLUSD(NROW,J))+XPLUSD(NROW,J)) + - XPLUSD(NROW,J) CALL SPVD(FCN, + N,M,NP,NQ, + BETA,XPLUSD,IFIXB,IFIXX,LDIFX, + NROW,J,LQ,STP0, + ISTOP,NFEV,PVPSTP, + WRK1,WRK2,WRK6) END IF IF (ISTOP.NE.0) THEN RETURN END IF FD = (PVPSTP-PV)/STP0 C CHECK FOR AGREEMENT IF (ABS(FD-D).LE.TOL*ABS(D)) THEN C NUMERICAL AND ANALYTIC DERIVATIVES AGREE C SET RELATIVE DIFFERENCE FOR DERIVATIVE CHECKING REPORT IF ((D.EQ.ZERO) .OR. (FD.EQ.ZERO)) THEN DIFFJ = ABS(FD-D) ELSE DIFFJ = ABS(FD-D)/ABS(D) END IF C SET MSG FLAG. IF (D.EQ.ZERO) THEN C JTH ANALYTIC AND NUMERICAL DERIVATIVES ARE BOTH ZERO. MSG(LQ,J) = 1 ELSE C JTH ANALYTIC AND NUMERICAL DERIVATIVES ARE BOTH NONZERO. MSG(LQ,J) = 0 END IF ELSE C NUMERICAL AND ANALYTIC DERIVATIVES DISAGREE. CHECK WHY IF ((D.EQ.ZERO) .OR. (FD.EQ.ZERO)) THEN CALL SJCKZ(FCN, + N,M,NP,NQ, + BETA,XPLUSD,IFIXB,IFIXX,LDIFX, + NROW,EPSMAC,J,LQ,ISWRTB, + TOL,D,FD,TYPJ,PVPSTP,STP0,PV, + DIFFJ,MSG,ISTOP,NFEV, + WRK1,WRK2,WRK6) ELSE CALL SJCKC(FCN, + N,M,NP,NQ, + BETA,XPLUSD,IFIXB,IFIXX,LDIFX, + ETA,TOL,NROW,EPSMAC,J,LQ,HC,ISWRTB, + FD,TYPJ,PVPSTP,STP0,PV,D, + DIFFJ,MSG,ISTOP,NFEV, + WRK1,WRK2,WRK6) END IF IF (MSG(LQ,J).LE.2) THEN GO TO 20 END IF END IF 10 CONTINUE C SET SUMMARY FLAG TO INDICATE QUESTIONABLE RESULTS 20 CONTINUE IF ((MSG(LQ,J).GE.7) .AND. (DIFFJ.LE.TOL2)) MSG(LQ,J) = 6 IF ((MSG(LQ,J).GE.1) .AND. (MSG(LQ,J).LE.6)) THEN MSG1 = MAX(MSG1,1) ELSE IF (MSG(LQ,J).GE.7) THEN MSG1 = 2 END IF RETURN END *SJCKZ SUBROUTINE SJCKZ + (FCN, + N,M,NP,NQ, + BETA,XPLUSD,IFIXB,IFIXX,LDIFX, + NROW,EPSMAC,J,LQ,ISWRTB, + TOL,D,FD,TYPJ,PVPSTP,STP0,PV, + DIFFJ,MSG,ISTOP,NFEV, + WRK1,WRK2,WRK6) C***BEGIN PROLOGUE SJCKZ C***REFER TO SODR,SODRC C***ROUTINES CALLED SPVB,SPVD C***DATE WRITTEN 860529 (YYMMDD) C***REVISION DATE 920619 (YYMMDD) C***PURPOSE RECHECK THE DERIVATIVES IN THE CASE WHERE THE FINITE C DIFFERENCE DERIVATIVE DISAGREES WITH THE ANALYTIC C DERIVATIVE AND THE ANALYTIC DERIVATIVE IS ZERO C (ADAPTED FROM STARPAC SUBROUTINE DCKZRO) C***END PROLOGUE SJCKZ C...SCALAR ARGUMENTS REAL + D,DIFFJ,EPSMAC,FD,PV,PVPSTP,STP0,TOL,TYPJ INTEGER + ISTOP,J,LDIFX,LQ,M,N,NFEV,NP,NQ,NROW LOGICAL + ISWRTB C...ARRAY ARGUMENTS REAL + BETA(NP),WRK1(N,M,NQ),WRK2(N,NQ),WRK6(N,NP,NQ),XPLUSD(N,M) INTEGER + IFIXB(NP),IFIXX(LDIFX,M),MSG(NQ,J) C...SUBROUTINE ARGUMENTS EXTERNAL + FCN C...LOCAL SCALARS REAL + CD,ONE,PVMSTP,THREE,TWO,ZERO C...EXTERNAL SUBROUTINES EXTERNAL + SPVB,SPVD C...INTRINSIC FUNCTIONS INTRINSIC + ABS,MIN C...DATA STATEMENTS DATA + ZERO,ONE,TWO,THREE + /0.0E0,1.0E0,2.0E0,3.0E0/ C...ROUTINE NAMES USED AS SUBPROGRAM ARGUMENTS C FCN: THE USER SUPPLIED SUBROUTINE FOR EVALUATING THE MODEL. C...VARIABLE DEFINITIONS (ALPHABETICALLY) C BETA: THE FUNCTION PARAMETERS. C CD: THE CENTRAL DIFFERENCE DERIVATIVE WRT THE JTH PARAMETER. C D: THE DERIVATIVE WITH RESPECT TO THE JTH UNKNOWN PARAMETER. C DIFFJ: THE RELATIVE DIFFERENCES BETWEEN THE USER SUPPLIED AND C FINITE DIFFERENCE DERIVATIVES FOR THE DERIVATIVE BEING C CHECKED. C EPSMAC: THE VALUE OF MACHINE PRECISION. C FD: THE FORWARD DIFFERENCE DERIVATIVE WRT THE JTH PARAMETER. C IFIXB: THE VALUES DESIGNATING WHETHER THE ELEMENTS OF BETA ARE C FIXED AT THEIR INPUT VALUES OR NOT. C IFIXX: THE VALUES DESIGNATING WHETHER THE ELEMENTS OF X ARE C FIXED AT THEIR INPUT VALUES OR NOT. C ISTOP: THE VARIABLE DESIGNATING WHETHER THERE ARE PROBLEMS C COMPUTING THE FUNCTION AT THE CURRENT BETA AND DELTA. C ISWRTB: THE VARIABLE DESIGNATING WHETHER THE DERIVATIVES WRT BETA C (ISWRTB=TRUE) OR X (ISWRTB=FALSE) ARE BEING CHECKED. C J: THE INDEX OF THE PARTIAL DERIVATIVE BEING EXAMINED. C LDIFX: THE LEADING DIMENSION OF ARRAY IFIXX. C LQ: THE RESPONSE CURRENTLY BEING EXAMINED. C M: THE NUMBER OF COLUMNS OF DATA IN THE EXPLANATORY VARIABLE. C MSG: THE ERROR CHECKING RESULTS. C N: THE NUMBER OF OBSERVATIONS. C NFEV: THE NUMBER OF FUNCTION EVALUATIONS. C NP: THE NUMBER OF FUNCTION PARAMETERS. C NQ: THE NUMBER OF RESPONSES PER OBSERVATION. C NROW: THE ROW NUMBER OF THE EXPLANATORY VARIABLE ARRAY AT WHICH C THE DERIVATIVE IS TO BE CHECKED. C ONE: THE VALUE 1.0E0. C PV: THE PREDICTED VALUE FROM THE MODEL FOR ROW NROW . C PVMSTP: THE PREDICTED VALUE FOR ROW NROW OF THE MODEL C USING THE CURRENT PARAMETER ESTIMATES FOR ALL BUT THE C JTH PARAMETER VALUE, WHICH IS BETA(J) - STP0. C PVPSTP: THE PREDICTED VALUE FOR ROW NROW OF THE MODEL C USING THE CURRENT PARAMETER ESTIMATES FOR ALL BUT THE C JTH PARAMETER VALUE, WHICH IS BETA(J) + STP0. C STP0: THE INITIAL STEP SIZE FOR THE FINITE DIFFERENCE DERIVATIVE. C THREE: THE VALUE 3.0E0. C TWO: THE VALUE 2.0E0. C TOL: THE AGREEMENT TOLERANCE. C TYPJ: THE TYPICAL SIZE OF THE J-TH UNKNOWN BETA OR DELTA. C WRK1: A WORK ARRAY OF (N BY M BY NQ) ELEMENTS. C WRK2: A WORK ARRAY OF (N BY NQ) ELEMENTS. C WRK6: A WORK ARRAY OF (N BY NP BY NQ) ELEMENTS. C XPLUSD: THE VALUES OF X + DELTA. C ZERO: THE VALUE 0.0E0. C***FIRST EXECUTABLE STATEMENT SJCKZ C RECALCULATE NUMERICAL DERIVATIVE USING CENTRAL DIFFERENCE AND STEP C SIZE OF 2*STP0 IF (ISWRTB) THEN C PERFORM COMPUTATIONS FOR DERIVATIVES WRT BETA CALL SPVB(FCN, + N,M,NP,NQ, + BETA,XPLUSD,IFIXB,IFIXX,LDIFX, + NROW,J,LQ,-STP0, + ISTOP,NFEV,PVMSTP, + WRK1,WRK2,WRK6) ELSE C PERFORM COMPUTATIONS FOR DERIVATIVES WRT DELTA CALL SPVD(FCN, + N,M,NP,NQ, + BETA,XPLUSD,IFIXB,IFIXX,LDIFX, + NROW,J,LQ,-STP0, + ISTOP,NFEV,PVMSTP, + WRK1,WRK2,WRK6) END IF IF (ISTOP.NE.0) THEN RETURN END IF CD = (PVPSTP-PVMSTP)/(TWO*STP0) DIFFJ = MIN(ABS(CD-D),ABS(FD-D)) C CHECK FOR AGREEMENT IF (DIFFJ.LE.TOL*ABS(D)) THEN C FINITE DIFFERENCE AND ANALYTIC DERIVATIVES NOW AGREE. IF (D.EQ.ZERO) THEN MSG(LQ,J) = 1 ELSE MSG(LQ,J) = 0 END IF ELSE IF (DIFFJ*TYPJ.LE.ABS(PV*EPSMAC**(ONE/THREE))) THEN C DERIVATIVES ARE BOTH CLOSE TO ZERO MSG(LQ,J) = 2 ELSE C DERIVATIVES ARE NOT BOTH CLOSE TO ZERO MSG(LQ,J) = 3 END IF RETURN END *SODCHK SUBROUTINE SODCHK + (N,M,NP,NQ, + ISODR,ANAJAC,IMPLCT, + IFIXB, + LDX,LDIFX,LDSCLD,LDSTPD,LDWE,LD2WE,LDWD,LD2WD, + LDY, + LWORK,LWKMN,LIWORK,LIWKMN, + SCLB,SCLD,STPB,STPD, + INFO) C***BEGIN PROLOGUE SODCHK C***REFER TO SODR,SODRC C***ROUTINES CALLED (NONE) C***DATE WRITTEN 860529 (YYMMDD) C***REVISION DATE 920619 (YYMMDD) C***PURPOSE CHECK INPUT PARAMETERS, INDICATING ERRORS FOUND USING C NONZERO VALUES OF ARGUMENT INFO C***END PROLOGUE SODCHK C...SCALAR ARGUMENTS INTEGER + INFO,LDIFX,LDSCLD,LDSTPD,LDWD,LDWE,LDX,LDY,LD2WD,LD2WE, + LIWKMN,LIWORK,LWKMN,LWORK,M,N,NP,NQ LOGICAL + ANAJAC,IMPLCT,ISODR C...ARRAY ARGUMENTS REAL + SCLB(NP),SCLD(LDSCLD,M),STPB(NP),STPD(LDSTPD,M) INTEGER + IFIXB(NP) C...LOCAL SCALARS INTEGER + I,J,K,LAST,NPP C...VARIABLE DEFINITIONS (ALPHABETICALLY)