*DODR SUBROUTINE DODR + (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 DODR 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 DOUBLE 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 DODCNT C***END PROLOGUE DODR C...SCALAR ARGUMENTS INTEGER + INFO,JOB,LDWD,LDWE,LDX,LDY,LD2WD,LD2WE,LIWORK,LWORK, + M,N,NDIGIT,NP,NQ C...ARRAY ARGUMENTS DOUBLE PRECISION + 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 DOUBLE PRECISION + NEGONE,PARTOL,SSTOL,TAUFAC,ZERO INTEGER + IPRINT,LDIFX,LDSCLD,LDSTPD,LUNERR,LUNRPT,MAXIT LOGICAL + SHORT C...LOCAL ARRAYS DOUBLE PRECISION + SCLB(1),SCLD(1,1),STPB(1),STPD(1,1),WD1(1,1,1) INTEGER + IFIXB(1),IFIXX(1,1) C...EXTERNAL SUBROUTINES EXTERNAL + DODCNT C...DATA STATEMENTS DATA + NEGONE,ZERO + /-1.0D0,0.0D0/ 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.0D0. 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.0D0. C WE: THE EPSILON WEIGHTS. C WORK: THE DOUBLE PRECISION WORK SPACE. C X: THE EXPLANATORY VARIABLE. C Y: THE DEPENDENT VARIABLE. UNUSED WHEN THE MODEL IS IMPLICIT. C***FIRST EXECUTABLE STATEMENT DODR 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 DODCNT + (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 DODCNT + (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 *DODRC SUBROUTINE DODRC + (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 DODRC 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 DOUBLE 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 DODCNT C***END PROLOGUE DODRC C...SCALAR ARGUMENTS DOUBLE PRECISION + 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 DOUBLE PRECISION + 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 DOUBLE PRECISION + NEGONE,ZERO LOGICAL + SHORT C...LOCAL ARRAYS DOUBLE PRECISION + WD1(1,1,1) C...EXTERNAL SUBROUTINES EXTERNAL + DODCNT C...DATA STATEMENTS DATA + NEGONE,ZERO + /-1.0D0,0.0D0/ 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.0D0. C WE: THE EPSILON WEIGHTS. C WORK: THE DOUBLE PRECISION WORK SPACE. C X: THE EXPLANATORY VARIABLE. C Y: THE DEPENDENT VARIABLE. UNUSED WHEN THE MODEL IS IMPLICIT. C***FIRST EXECUTABLE STATEMENT DODRC SHORT = .FALSE. IF (WD(1,1,1).NE.ZERO) THEN CALL DODCNT + (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 DODCNT + (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 *DACCES SUBROUTINE DACCES + (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 DACCES C***REFER TO DODR,DODRC C***ROUTINES CALLED DIWINF,DWINF C***DATE WRITTEN 860529 (YYMMDD) C***REVISION DATE 920619 (YYMMDD) C***PURPOSE ACCESS OR STORE VALUES IN THE WORK ARRAYS C***END PROLOGUE DACESS C...SCALAR ARGUMENTS DOUBLE PRECISION + 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 DOUBLE PRECISION + 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 + DIWINF,DWINF 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 DOUBLE PRECISION 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 DACCES C FIND STARTING LOCATIONS WITHIN INTEGER WORKSPACE CALL DIWINF(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 DOUBLE PRECISION WORK SPACE CALL DWINF(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 *DETAF SUBROUTINE DETAF + (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 DETAF C***REFER TO DODR,DODRC 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 DETAF C...SCALAR ARGUMENTS DOUBLE PRECISION + EPSMAC,ETA INTEGER + ISTOP,LDIFX,M,N,NETA,NFEV,NP,NQ,NROW C...ARRAY ARGUMENTS DOUBLE PRECISION + 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 DOUBLE PRECISION + 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.0D0,0.1D0,0.2D0,0.5D0,1.0D0,2.0D0,1.0D2/ 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.0D2. 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.0D0. C P1: THE VALUE 0.1D0. C P2: THE VALUE 0.2D0. C P5: THE VALUE 0.5D0. 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.0D0. C***FIRST EXECUTABLE STATEMENT DETAF 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 *DFCTR SUBROUTINE DFCTR(OKSEMI,A,LDA,N,INFO) C***BEGIN PROLOGUE DFCTR C***REFER TO DODR,DODRC C***ROUTINES CALLED DDOT 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 DPOFA) C***REFERENCES DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W., C *LINPACK USERS GUIDE*, SIAM, 1979. C***END PROLOGUE DFCTR C...SCALAR ARGUMENTS INTEGER INFO,LDA,N LOGICAL OKSEMI C...ARRAY ARGUMENTS DOUBLE PRECISION A(LDA,N) C...LOCAL SCALARS DOUBLE PRECISION XI,S,T,TEN,ZERO INTEGER J,K C...EXTERNAL FUNCTIONS EXTERNAL DMPREC,DDOT DOUBLE PRECISION DMPREC,DDOT C...INTRINSIC FUNCTIONS INTRINSIC SQRT C...DATA STATEMENTS DATA + ZERO,TEN + /0.0D0,10.0D0/ 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.0D0. C XI: A VALUE USED TO TEST FOR NON POSITIVE SEMIDEFINITENESS. C ZERO: THE VALUE 0.0D0. C***FIRST EXECUTABLE STATEMENT DFCTR C SET RELATIVE TOLERANCE FOR DETECTING NON POSITIVE SEMIDEFINITENESS. XI = -TEN*DMPREC() 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) - DDOT(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 *DFCTRW SUBROUTINE DFCTRW + (N,M,NQ,NPP, + ISODR, + WE,LDWE,LD2WE,WD,LDWD,LD2WD, + WRK0,WRK4, + WE1,NNZW,INFO) C***BEGIN PROLOGUE DFCTRW C***REFER TO DODR,DODRC C***ROUTINES CALLED DFCTR 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 DFCTRW C...SCALAR ARGUMENTS INTEGER + INFO,LDWD,LDWE,LD2WD,LD2WE, + M,N,NNZW,NPP,NQ LOGICAL + ISODR C...ARRAY ARGUMENTS DOUBLE PRECISION + WE(LDWE,LD2WE,NQ),WE1(LDWE,LD2WE,NQ),WD(LDWD,LD2WD,M), + WRK0(NQ,NQ),WRK4(M,M) C...LOCAL SCALARS DOUBLE PRECISION + ZERO INTEGER + I,INF,J,J1,J2,L,L1,L2 LOGICAL + NOTZRO C...EXTERNAL SUBROUTINES EXTERNAL + DFCTR C...DATA STATEMENTS DATA + ZERO + /0.0D0/ 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.0D0. C***FIRST EXECUTABLE STATEMENT DFCTRW 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 DFCTR(.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 DFCTR(.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 DFCTR(.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 DFCTR(.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 *DJACCD SUBROUTINE DJACCD + (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 DJACCD C***REFER TO DODR,DODRC C***ROUTINES CALLED FCN,DHSTEP,DZERO 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 DJACCD C...SCALAR ARGUMENTS INTEGER + ISTOP,LDIFX,LDSTPD,LDTT,LDX,M,N,NETA,NFEV,NP,NQ LOGICAL + ISODR C...ARRAY ARGUMENTS DOUBLE PRECISION + 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 DOUBLE PRECISION + BETAK,ONE,TYPJ,ZERO INTEGER + I,J,K,L LOGICAL + DOIT,SETZRO C...EXTERNAL SUBROUTINES EXTERNAL + DZERO C...EXTERNAL FUNCTIONS DOUBLE PRECISION + DHSTEP EXTERNAL + DHSTEP C...INTRINSIC FUNCTIONS INTRINSIC + ABS,MAX,SIGN,SQRT C...DATA STATEMENTS DATA + ZERO,ONE + /0.0D0,1.0D0/ 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.0D0. 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.0D0. C***FIRST EXECUTABLE STATEMENT DJACCD 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 DZERO(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*DHSTEP(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 DZERO(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*DHSTEP(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 *DJACFD SUBROUTINE DJACFD + (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 DJACFD C***REFER TO DODR,DODRC C***ROUTINES CALLED FCN,DHSTEP,DZERO 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 DJACFD C...SCALAR ARGUMENTS INTEGER + ISTOP,LDIFX,LDSTPD,LDTT,LDX,M,N,NETA,NFEV,NP,NQ LOGICAL + ISODR C...ARRAY ARGUMENTS DOUBLE PRECISION + 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 DOUBLE PRECISION + BETAK,ONE,TYPJ,ZERO INTEGER + I,J,K,L LOGICAL + DOIT,SETZRO C...EXTERNAL SUBROUTINES EXTERNAL + DZERO C...EXTERNAL FUNCTIONS DOUBLE PRECISION + DHSTEP EXTERNAL + DHSTEP C...INTRINSIC FUNCTIONS INTRINSIC + ABS,MAX,SIGN,SQRT C...DATA STATEMENTS DATA + ZERO,ONE + /0.0D0,1.0D0/ 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.0D0. 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.0D0. C***FIRST EXECUTABLE STATEMENT DJACFD 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 DZERO(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*DHSTEP(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 DZERO(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*DHSTEP(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 *DJCK SUBROUTINE DJCK + (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 DJCK C***REFER TO DODR,DODRC C***ROUTINES CALLED FCN,DHSTEP,DJCKM 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 DJCK C...SCALAR ARGUMENTS DOUBLE PRECISION + EPSMAC,ETA INTEGER + ISTOP,LDIFX,LDSTPD,LDTT, + M,N,NETA,NFEV,NJEV,NP,NQ,NROW,NTOL LOGICAL + ISODR C...ARRAY ARGUMENTS DOUBLE PRECISION + 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 DOUBLE PRECISION + DIFFJ,H0,HC0,ONE,P5,PV,TOL,TYPJ,ZERO INTEGER + IDEVAL,J,LQ,MSGB1,MSGD1 LOGICAL + ISFIXD,ISWRTB C...EXTERNAL SUBROUTINES EXTERNAL + DJCKM C...EXTERNAL FUNCTIONS DOUBLE PRECISION + DHSTEP EXTERNAL + DHSTEP C...INTRINSIC FUNCTIONS INTRINSIC + ABS,INT,LOG10 C...DATA STATEMENTS DATA + ZERO,P5,ONE + /0.0D0,0.5D0,1.0D0/ 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 DETAF. 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.0D0. C P5: THE VALUE 0.5D0. 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.0D0. C***FIRST EXECUTABLE STATEMENT DJCK C SET TOLERANCE FOR CHECKING DERIVATIVES TOL = ETA**(0.25D0) 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 = DHSTEP(0,NETA,1,J,STPB,1) HC0 = H0 C CHECK DERIVATIVE WRT THE J-TH PARAMETER AT THE NROW-TH ROW CALL DJCKM(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 = DHSTEP(0,NETA,NROW,J,STPD,LDSTPD) HC0 = DHSTEP(1,NETA,NROW,J,STPD,LDSTPD) C CHECK DERIVATIVE WRT THE J-TH COLUMN OF DELTA AT ROW NROW CALL DJCKM(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 *DJCKC SUBROUTINE DJCKC + (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 DJCKC C***REFER TO DODR,DODRC C***ROUTINES CALLED DJCKF,DPVB,DPVD 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 DJCKC C...SCALAR ARGUMENTS DOUBLE PRECISION + 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 DOUBLE PRECISION + 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 DOUBLE PRECISION + CURVE,ONE,PVMCRV,PVPCRV,P01,STP,STPCRV,TEN,TWO C...EXTERNAL SUBROUTINES EXTERNAL + DJCKF,DPVB,DPVD C...INTRINSIC FUNCTIONS INTRINSIC + ABS,SIGN C...DATA STATEMENTS DATA + P01,ONE,TWO,TEN + /0.01D0,1.0D0,2.0D0,10.0D0/ 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.0D0. 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.01D0. 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.0D0. C TOL: THE AGREEMENT TOLERANCE. C TWO: THE VALUE 2.0D0. 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 DJCKC IF (ISWRTB) THEN C PERFORM CENTRAL DIFFERENCE COMPUTATIONS FOR DERIVATIVES WRT BETA STPCRV = (HC*TYPJ*SIGN(ONE,BETA(J))+BETA(J)) - BETA(J) CALL DPVB(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 DPVB(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 DPVD(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 DPVD(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 DJCKF(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 DPVB(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 DPVD(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 *DJCKF SUBROUTINE DJCKF + (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 DJCKF C***REFER TO DODR,DODRC C***ROUTINES CALLED DPVB,DPVD 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 DJCKF C...SCALAR ARGUMENTS DOUBLE PRECISION + 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 DOUBLE PRECISION + 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 DOUBLE PRECISION + HUNDRD,ONE,P1,STP,TWO LOGICAL + LARGE C...EXTERNAL SUBROUTINES EXTERNAL + DPVB,DPVD C...INTRINSIC FUNCTIONS INTRINSIC + ABS,SIGN C...DATA STATEMENTS DATA + P1,ONE,TWO,HUNDRD + /0.1D0,1.0D0,2.0D0,100.0D0/ 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.0D0. 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.0D0. 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.1D0. C STP0: THE STEP SIZE FOR THE FINITE DIFFERENCE DERIVATIVE. C TOL: THE AGREEMENT TOLERANCE. C TWO: THE VALUE 2.0D0. 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 DJCKF 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 DPVB(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 DPVD(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 *DJCKM SUBROUTINE DJCKM + (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 DJCKM C***REFER TO DODR,DODRC C***ROUTINES CALLED DJCKC,DJCKZ,DPVB,DPVD 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 DJCKM C...SCALAR ARGUMENTS DOUBLE PRECISION + 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 DOUBLE PRECISION + 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 DOUBLE PRECISION + BIG,FD,H,HC,H1,HC1,HUNDRD,ONE,PVPSTP,P01,P1,STP0, + TEN,THREE,TOL2,TWO,ZERO INTEGER + I C...EXTERNAL SUBROUTINES EXTERNAL + DJCKC,DJCKZ,DPVB,DPVD C...INTRINSIC FUNCTIONS INTRINSIC + ABS,MAX,SIGN,SQRT C...DATA STATEMENTS DATA + ZERO,P01,P1,ONE,TWO,THREE,TEN,HUNDRD + /0.0D0,0.01D0,0.1D0,1.0D0,2.0D0,3.0D0,1.0D1,1.0D2/ DATA + BIG,TOL2 + /1.0D19,5.0D-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.0D0. 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.0D0. 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.01D0. C P1: THE VALUE 0.1D0. C STP0: THE INITIAL STEP SIZE FOR THE FINITE DIFFERENCE DERIVATIVE. C TEN: THE VALUE 10.0D0. C THREE: THE VALUE 3.0D0. C TWO: THE VALUE 2.0D0. 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.0D0. C***FIRST EXECUTABLE STATEMENT DJCKM 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 DPVB(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 DPVD(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 DJCKZ(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 DJCKC(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 *DJCKZ SUBROUTINE DJCKZ + (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 DJCKZ C***REFER TO DODR,DODRC C***ROUTINES CALLED DPVB,DPVD 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 DJCKZ C...SCALAR ARGUMENTS DOUBLE PRECISION + 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 DOUBLE PRECISION + 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 DOUBLE PRECISION + CD,ONE,PVMSTP,THREE,TWO,ZERO C...EXTERNAL SUBROUTINES EXTERNAL + DPVB,DPVD C...INTRINSIC FUNCTIONS INTRINSIC + ABS,MIN C...DATA STATEMENTS DATA + ZERO,ONE,TWO,THREE + /0.0D0,1.0D0,2.0D0,3.0D0/ 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.0D0. 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.0D0. C TWO: THE VALUE 2.0D0. 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.0D0. C***FIRST EXECUTABLE STATEMENT DJCKZ 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 DPVB(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 DPVD(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 *DODCHK SUBROUTINE DODCHK + (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 DODCHK C***REFER TO DODR,DODRC 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 DODCHK 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 DOUBLE PRECISION + 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) C ANAJAC: THE VARIABLE DESIGNATING WHETHER THE JACOBIANS ARE C COMPUTED BY FINITE DIFFERENCES (ANAJAC=FALSE) OR NOT C (ANAJAC=TRUE). 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 IMPLCT: THE VARIABLE DESIGNATING WHETHER THE SOLUTION IS BY C IMPLICIT ODR (IMPLCT=TRUE) OR EXPLICIT ODR (IMPLCT=FALSE). 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 K: AN INDEXING VARIABLE. C LAST: THE LAST ROW OF THE ARRAY TO BE ACCESSED. 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 X. C LD2WD: THE SECOND DIMENSION OF ARRAY WD. C LD2WE: THE SECOND DIMENSION OF ARRAY WE. C LIWKMN: THE MINIMUM ACCEPTABLE LENGTH OF ARRAY IWORK. C LIWORK: THE LENGTH OF VECTOR IWORK. 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 N: THE NUMBER OF OBSERVATIONS. C NP: THE NUMBER OF FUNCTION PARAMETERS. C NPP: THE NUMBER OF FUNCTION PARAMETERS BEING ESTIMATED. C NQ: THE NUMBER OF RESPONSES PER OBSERVATIONS. C SCLB: THE SCALING VALUES FOR BETA. C SCLD: THE SCALING VALUE FOR DELTA. C STPB: THE STEP FOR THE FINITE DIFFERENCE DERIVITIVE WRT BETA. C STPD: THE STEP FOR THE FINITE DIFFERENCE DERIVITIVE WRT DELTA. C***FIRST EXECUTABLE STATEMENT DODCHK C FIND ACTUAL NUMBER OF PARAMETERS BEING ESTIMATED IF (NP.LE.0 .OR. IFIXB(1).LT.0) THEN NPP = NP ELSE NPP = 0 DO 10 K=1,NP IF (IFIXB(K).NE.0) THEN NPP = NPP + 1 END IF 10 CONTINUE END IF C CHECK PROBLEM SPECIFICATION PARAMETERS IF (N.LE.0 .OR. + M.LE.0 .OR. + (NPP.LE.0 .OR. NPP.GT.N) .OR. + (NQ.LE.0)) THEN INFO = 10000 IF (N.LE.0) THEN INFO = INFO + 1000 END IF IF (M.LE.0) THEN INFO = INFO + 100 END IF IF (NPP.LE.0 .OR. NPP.GT.N) THEN INFO = INFO + 10 END IF IF (NQ.LE.0) THEN INFO = INFO + 1 END IF RETURN END IF C CHECK DIMENSION SPECIFICATION PARAMETERS IF ((.NOT.IMPLCT .AND. LDY.LT.N) .OR. + (LDX.LT.N) .OR. + (LDWE.NE.1 .AND. LDWE.LT.N) .OR. + (LD2WE.NE.1 .AND. LD2WE.LT.NQ) .OR. + (ISODR .AND. (LDWD.NE.1 .AND. LDWD.LT.N)) .OR. + (ISODR .AND. (LD2WD.NE.1 .AND. LD2WD.LT.M)) .OR. + (ISODR .AND. (LDIFX.NE.1 .AND. LDIFX.LT.N)) .OR. + (ISODR .AND. (LDSTPD.NE.1 .AND. LDSTPD.LT.N)) .OR. + (ISODR .AND. (LDSCLD.NE.1 .AND. LDSCLD.LT.N)) .OR. + (LWORK.LT.LWKMN) .OR. + (LIWORK.LT.LIWKMN)) THEN INFO = 20000 IF (.NOT.IMPLCT .AND. LDY.LT.N) THEN INFO = INFO + 1000 END IF IF (LDX.LT.N) THEN INFO = INFO + 2000 END IF IF ((LDWE.NE.1 .AND. LDWE.LT.N) .OR. + (LD2WE.NE.1 .AND. LD2WE.LT.NQ)) THEN INFO = INFO + 100 END IF IF (ISODR .AND. ((LDWD.NE.1 .AND. LDWD.LT.N) .OR. + (LD2WD.NE.1 .AND. LD2WD.LT.M))) THEN INFO = INFO + 200 END IF IF (ISODR .AND. (LDIFX.NE.1 .AND. LDIFX.LT.N)) THEN INFO = INFO + 10 END IF IF (ISODR .AND. (LDSTPD.NE.1 .AND. LDSTPD.LT.N)) THEN INFO = INFO + 20 END IF IF (ISODR .AND. (LDSCLD.NE.1 .AND. LDSCLD.LT.N)) THEN INFO = INFO + 40 END IF IF (LWORK.LT.LWKMN) THEN INFO = INFO + 1 END IF IF (LIWORK.LT.LIWKMN) THEN INFO = INFO + 2 END IF RETURN END IF C CHECK DELTA SCALING IF (ISODR .AND. SCLD(1,1).GT.0) THEN IF (LDSCLD.GE.N) THEN LAST = N ELSE LAST = 1 END IF DO 120 J=1,M DO 110 I=1,LAST IF (SCLD(I,J).LE.0) THEN INFO = 30200 GO TO 130 END IF 110 CONTINUE 120 CONTINUE END IF 130 CONTINUE C CHECK BETA SCALING IF (SCLB(1).GT.0) THEN DO 210 K=1,NP IF (SCLB(K).LE.0) THEN IF (INFO.EQ.0) THEN INFO = 30100 ELSE INFO = INFO + 100 END IF GO TO 220 END IF 210 CONTINUE END IF 220 CONTINUE C CHECK DELTA FINITE DIFFERENCE STEP SIZES IF (ANAJAC .AND. ISODR .AND. STPD(1,1).GT.0) THEN IF (LDSTPD.GE.N) THEN LAST = N ELSE LAST = 1 END IF DO 320 J=1,M DO 310 I=1,LAST IF (STPD(I,J).LE.0) THEN IF (INFO.EQ.0) THEN INFO = 32000 ELSE INFO = INFO + 2000 END IF GO TO 330 END IF 310 CONTINUE 320 CONTINUE END IF 330 CONTINUE C CHECK BETA FINITE DIFFERENCE STEP SIZES IF (ANAJAC .AND. STPB(1).GT.0) THEN DO 410 K=1,NP IF (STPB(K).LE.0) THEN IF (INFO.EQ.0) THEN INFO = 31000 ELSE INFO = INFO + 1000 END IF GO TO 420 END IF 410 CONTINUE END IF 420 CONTINUE RETURN END *DODDRV SUBROUTINE DODDRV + (SHORT,HEAD,FSTITR,PRTPEN, + 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, + MAXIT1,TSTIMP, INFO) C***BEGIN PROLOGUE DODDRV C***REFER TO DODR,DODRC C***ROUTINES CALLED FCN,DCOPY,DDOT,DETAF,DFCTRW,DFLAGS, C DINIWK,DIWINF,DJCK,DNRM2,DODCHK,DODMN, C DODPER,DPACK,DSETN,DUNPAC,DWGHT,DWINF,DXMY,DXPY C***DATE WRITTEN 860529 (YYMMDD) C***REVISION DATE 920619 (YYMMDD) C***PURPOSE PERFORM ERROR CHECKING AND INITIALIZATION, AND BEGIN C PROCEDURE FOR PERFORMING ORTHOGONAL DISTANCE REGRESSION C (ODR) OR ORDINARY LINEAR OR NONLINEAR LEAST SQUARES (OLS) C***END PROLOGUE DODDRV C...SCALAR ARGUMENTS DOUBLE PRECISION + PARTOL,SSTOL,TAUFAC,TSTIMP INTEGER + INFO,IPRINT,JOB,LDIFX,LDSCLD,LDSTPD,LDWD,LDWE,LDX,LDY, + LD2WD,LD2WE,LIWORK,LUNERR,LUNRPT,LWORK,M,MAXIT,MAXIT1, + N,NDIGIT,NP,NQ LOGICAL + FSTITR,HEAD,PRTPEN,SHORT C...ARRAY ARGUMENTS DOUBLE PRECISION + BETA(NP),SCLB(NP),SCLD(LDSCLD,M),STPB(NP),STPD(LDSTPD,M), + WE(LDWE,LD2WE,NQ),WD(LDWD,LD2WD,M),WORK(LWORK), + X(LDX,M),Y(LDY,NQ) INTEGER + IFIXB(NP),IFIXX(LDIFX,M),IWORK(LIWORK) C...SUBROUTINE ARGUMENTS EXTERNAL + FCN C...LOCAL SCALARS DOUBLE PRECISION + EPSMAC,ETA,P5,ONE,TEN,ZERO INTEGER + ACTRSI,ALPHAI,BETACI,BETANI,BETASI,BETA0I,DELTAI,DELTNI,DELTSI, + DIFFI,EPSMAI,ETAI,FI,FJACBI,FJACDI,FNI,FSI,I,IDFI,INT2I,IPRINI, + IRANKI,ISTOP,ISTOPI,JOBI,JPVTI,K,LDTT,LDTTI,LIWKMN, + LUNERI,LUNRPI,LWKMN,LWRK,MAXITI,MSGB,MSGD,NETA,NETAI, + NFEV,NFEVI,NITERI,NJEV,NJEVI,NNZW,NNZWI,NPP,NPPI,NROW,NROWI, + NTOL,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,WRK, + WSSI,WSSDEI,WSSEPI,XPLUSI LOGICAL + ANAJAC,CDJAC,CHKJAC,DOVCV,IMPLCT,INITD,ISODR,REDOJ,RESTRT C...EXTERNAL FUNCTIONS DOUBLE PRECISION + DDOT,DNRM2 EXTERNAL + DDOT,DNRM2 C...EXTERNAL SUBROUTINES EXTERNAL + DCOPY,DETAF,DFCTRW,DFLAGS,DINIWK,DIWINF,DJCK,DODCHK, + DODMN,DODPER,DPACK,DSETN,DUNPAC,DWGHT,DWINF,DXMY,DXPY C...DATA STATEMENTS DATA + ZERO,P5,ONE,TEN + /0.0D0,0.5D0,1.0D0,10.0D0/ C...ROUTINE NAMES USED AS SUBPROGRAM ARGUMENTS C FCN: THE USER SUPPLIED SUBROUTINE FOR EVALUATING THE MODEL. C...VARIABLE DEFINITIONS (ALPHABETICALLY) C ACTRSI: THE LOCATION IN ARRAY WORK OF VARIABLE ACTRS. C ALPHAI: THE LOCATION IN ARRAY WORK OF VARIABLE ALPHA. 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 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 CDJAC: THE VARIABLE DESIGNATING WHETHER THE JACOBIANS ARE C COMPUTED BY CENTRAL DIFFERENCES (CDJAC=TRUE) OR 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 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 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 ETA: THE RELATIVE NOISE IN THE FUNCTION RESULTS. C ETAI: THE LOCATION IN ARRAY WORK OF VARIABLE ETA. C FI: THE STARTING LOCATION IN ARRAY WORK OF ARRAY F. 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 FSTITR: THE VARIABLE DESIGNATING WHETHER THIS IS THE FIRST C ITERATION (FSTITR=TRUE) OR NOT (FSTITR=FALSE). C HEAD: THE VARIABLE DESIGNATING WHETHER THE HEADING IS TO BE C PRINTED (HEAD=TRUE) OR NOT (HEAD=FALSE). C I: AN INDEX VARIABLE. C IDFI: THE LOCATION IN ARRAY IWORK OF VARIABLE IDF. 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 IMPLCT: THE VARIABLE DESIGNATING WHETHER THE SOLUTION IS BY C IMPLICIT ODR (IMPLCT=TRUE) OR EXPLICIT ODR (IMPLCT=FALSE). C INFO: THE VARIABLE DESIGNATING WHY THE COMPUTATIONS WERE STOPPED. 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 INT2I: THE IN ARRAY IWORK OF VARIABLE INT2. C IPRINI: THE LOCATION IN ARRAY IWORK OF VARIABLE IPRINT. C IPRINT: THE PRINT CONTROL VARIABLE. C IRANKI: THE LOCATION IN ARRAY IWORK OF VARIABLE IRANK. 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 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 JPVTI: THE STARTING LOCATION IN ARRAY IWORK OF ARRAY JPVT. C K: AN INDEX VARIABLE. 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 LDTT: THE LEADING DIMENSION OF ARRAY TT. C LDTTI: THE LOCATION IN ARRAY IWORK OF VARIABLE LDTT. 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 LIWKMN: THE MINIMUM ACCEPTABLE LENGTH OF ARRAY IWORK. 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 LWRK: THE LENGTH OF VECTOR WRK. C M: THE NUMBER OF COLUMNS OF DATA IN THE EXPLANATORY VARIABLE. C MAXIT: THE MAXIMUM NUMBER OF ITERATIONS ALLOWED. C MAXIT1: FOR IMPLICIT MODELS, THE ITERATIONS ALLOWED FOR THE NEXT C PENALTY PARAMETER VALUE. 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 NDIGIT: THE NUMBER OF ACCURATE DIGITS IN THE FUNCTION RESULTS, AS C SUPPLIED BY THE USER. 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 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 OBSERVATIONAL ERROR WEIGHTS. C NNZWI: THE LOCATION IN ARRAY IWORK OF VARIABLE NNZW. C NP: THE NUMBER OF FUNCTION PARAMETERS. C NPP: THE NUMBER OF FUNCTION PARAMETERS BEING ESTIMATED. C NPPI: THE LOCATION IN ARRAY IWORK OF VARIABLE NPP. C NQ: THE NUMBER OF RESPONSES PER OBSERVATION. C NROW: THE ROW NUMBER AT WHICH THE DERIVATIVE IS TO BE CHECKED. C NROWI: THE LOCATION IN ARRAY IWORK OF VARIABLE NROW. C NTOL: THE NUMBER OF DIGITS OF AGREEMENT REQUIRED BETWEEN THE C NUMERICAL DERIVATIVES AND THE USER SUPPLIED DERIVATIVES, C SET BY DJCK. C NTOLI: THE LOCATION IN ARRAY IWORK OF VARIABLE NTOL. C OLMAVI: THE LOCATION IN ARRAY WORK OF VARIABLE OLMAVG. C OMEGAI: THE STARTING LOCATION IN ARRAY WORK OF ARRAY OMEGA. C ONE: THE VALUE 1.0D0. 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 PRERSI: THE LOCATION IN ARRAY WORK OF VARIABLE PRERS. C PRTPEN: THE VARIABLE DESIGNATING WHETHER THE PENALTY PARAMETER IS C TO BE PRINTED IN THE ITERATION REPORT (PRTPEN=TRUE) OR NOT C (PRTPEN=FALSE). C P5: THE VALUE 0.5D0. C QRAUXI: THE STARTING LOCATION IN ARRAY WORK OF ARRAY QRAUX. C RCONDI: THE LOCATION IN ARRAY WORK OF VARIABLE RCOND. 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 RNORSI: THE LOCATION IN ARRAY WORK OF VARIABLE RNORMS. C RVARI: THE LOCATION IN ARRAY WORK OF VARIABLE RVAR. C SCLB: THE SCALING VALUES FOR BETA. C SCLD: THE SCALING VALUES FOR DELTA. 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-CALL C (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 STPB: THE STEP SIZE FOR FINITE DIFFERENCE DERIVATIVES WRT BETA. C STPD: THE STEP SIZE FOR FINITE DIFFERENCE DERIVATIVES WRT DELTA. 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 TEN: THE VALUE 10.0D0. C TI: THE STARTING LOCATION IN ARRAY WORK OF ARRAY T. C TSTIMP: THE RELATIVE CHANGE IN THE PARAMETERS BETWEEN THE INITIAL C VALUES AND THE SOLUTION. C TTI: THE STARTING LOCATION IN ARRAY WORK OF ARRAY TT. C UI: THE STARTING LOCATION IN ARRAY WORK OF ARRAY U. C VCVI: THE STARTING LOCATION IN ARRAY WORK OF ARRAY VCV. C WD: THE DELTA WEIGHTS. C WE: THE EPSILON WEIGHTS. C WE1I: THE STARTING LOCATION IN ARRAY WORK OF ARRAY WE1. C WORK: THE DOUBLE PRECISION WORK SPACE. C WRK: THE STARTING LOCATION IN ARRAY WORK OF ARRAY WRK, C EQUIVALENCED TO WRK1 AND WRK2. C WRK1I: THE STARTING LOCATION IN ARRAY WORK OF ARRAY WRK1. C WRK2I: THE STARTING LOCATION IN ARRAY WORK OF ARRAY WRK2. C WRK3I: THE STARTING LOCATION IN ARRAY WORK OF ARRAY WRK3. C WRK4I: THE STARTING LOCATION IN ARRAY WORK OF ARRAY WRK4. C WRK5I: THE STARTING LOCATION IN ARRAY WORK OF ARRAY WRK5. C WRK6I: THE STARTING LOCATION IN ARRAY WORK OF ARRAY WRK6. C WRK7I: THE STARTING LOCATION IN ARRAY WORK OF ARRAY WRK7. C WSSI: THE LOCATION IN ARRAY WORK OF VARIABLE WSS. C WSSDEI: THE LOCATION IN ARRAY WORK OF VARIABLE WSSDEL. C WSSEPI: THE LOCATION IN ARRAY WORK OF VARIABLE WSSEPS. C X: THE EXPLANATORY VARIABLE. C XPLUSI: THE STARTING LOCATION IN ARRAY WORK OF ARRAY XPLUSD. C Y: THE DEPENDENT VARIABLE. UNUSED WHEN THE MODEL IS IMPLICIT. C ZERO: THE VALUE 0.0D0. C***FIRST EXECUTABLE STATEMENT DODDRV C INITIALIZE NECESSARY VARIABLES CALL DFLAGS(JOB,RESTRT,INITD,DOVCV,REDOJ, + ANAJAC,CDJAC,CHKJAC,ISODR,IMPLCT) C SET STARTING LOCATIONS WITHIN INTEGER WORKSPACE C (INVALID VALUES OF M, NP AND/OR NQ ARE HANDLED REASONABLY BY DIWINF) CALL DIWINF(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 SET STARTING LOCATIONS WITHIN DOUBLE PRECISION WORK SPACE C (INVALID VALUES OF N, M, NP, NQ, LDWE AND/OR LD2WE C ARE HANDLED REASONABLY BY DWINF) CALL DWINF(N,M,NP,NQ,LDWE,LD2WE,ISODR, + DELTAI,FI,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 (ISODR) THEN WRK = WRK1I LWRK = N*M*NQ + N*NQ ELSE WRK = WRK2I LWRK = N*NQ END IF C UPDATE THE PENALTY PARAMETERS C (WE(1,1,1) IS NOT A USER SUPPLIED ARRAY IN THIS CASE) IF (RESTRT .AND. IMPLCT) THEN WE(1,1,1) = MAX(WORK(WE1I)**2,ABS(WE(1,1,1))) WORK(WE1I) = -SQRT(ABS(WE(1,1,1))) END IF IF (RESTRT) THEN C RESET MAXIMUM NUMBER OF ITERATIONS IF (MAXIT.GE.0) THEN IWORK(MAXITI) = IWORK(NITERI) + MAXIT ELSE IWORK(MAXITI) = IWORK(NITERI) + 10 END IF IF (IWORK(NITERI).LT.IWORK(MAXITI)) THEN INFO = 0 END IF IF (JOB.GE.0) IWORK(JOBI) = JOB IF (IPRINT.GE.0) IWORK(IPRINI) = IPRINT IF (PARTOL.GE.ZERO .AND. PARTOL.LT.ONE) WORK(PARTLI) = PARTOL IF (SSTOL.GE.ZERO .AND. SSTOL.LT.ONE) WORK(SSTOLI) = SSTOL WORK(OLMAVI) = WORK(OLMAVI)*IWORK(NITERI) IF (IMPLCT) THEN CALL DCOPY(N*NQ,WORK(FNI),1,WORK(FI),1) ELSE CALL DXMY(N,NQ,WORK(FNI),N,Y,LDY,WORK(FI),N) END IF CALL DWGHT(N,NQ,WORK(WE1I),LDWE,LD2WE,WORK(FI),N,WORK(FI),N) WORK(WSSEPI) = DDOT(N*NQ,WORK(FI),1,WORK(FI),1) WORK(WSSI) = WORK(WSSEPI) + WORK(WSSDEI) ELSE C PERFORM ERROR CHECKING INFO = 0 CALL DODCHK(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) IF (INFO.GT.0) THEN GO TO 50 END IF C INITIALIZE WORK VECTORS AS NECESSARY DO 10 I=N*M+N*NQ+1,LWORK WORK(I) = ZERO 10 CONTINUE DO 20 I=1,LIWORK IWORK(I) = 0 20 CONTINUE CALL DINIWK(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) IWORK(MSGB) = -1 IWORK(MSGD) = -1 WORK(TAUI) = -WORK(TAUFCI) C SET UP FOR PARAMETER ESTIMATION - C PULL BETA'S TO BE ESTIMATED AND CORRESPONDING SCALE VALUES C AND STORE IN WORK(BETACI) AND WORK(SSI), RESPECTIVELY CALL DPACK(NP,IWORK(NPPI),WORK(BETACI),BETA,IFIXB) CALL DPACK(NP,IWORK(NPPI),WORK(SSI),WORK(SSFI),IFIXB) NPP = IWORK(NPPI) C CHECK THAT WD IS POSITIVE DEFINITE AND WE IS POSITIVE SEMIDEFINITE, C SAVING FACTORIZATION OF WE, AND COUNTING NUMBER OF NONZERO WEIGHTS CALL DFCTRW(N,M,NQ,NPP, + ISODR, + WE,LDWE,LD2WE,WD,LDWD,LD2WD, + WORK(WRK2I),WORK(WRK4I), + WORK(WE1I),NNZW,INFO) IWORK(NNZWI) = NNZW IF (INFO.NE.0) THEN GO TO 50 END IF C EVALUATE THE PREDICTED VALUES AND C WEIGHTED EPSILONS AT THE STARTING POINT CALL DUNPAC(NP,WORK(BETACI),BETA,IFIXB) CALL DXPY(N,M,X,LDX,WORK(DELTAI),N,WORK(XPLUSI),N) ISTOP = 0 CALL FCN(N,M,NP,NQ, + N,M,NP, + BETA,WORK(XPLUSI), + IFIXB,IFIXX,LDIFX, + 002,WORK(FNI),WORK(WRK6I),WORK(WRK1I), + ISTOP) IWORK(ISTOPI) = ISTOP IF (ISTOP.EQ.0) THEN IWORK(NFEVI) = IWORK(NFEVI) + 1 IF (IMPLCT) THEN CALL DCOPY(N*NQ,WORK(FNI),1,WORK(FI),1) ELSE CALL DXMY(N,NQ,WORK(FNI),N,Y,LDY,WORK(FI),N) END IF CALL DWGHT(N,NQ,WORK(WE1I),LDWE,LD2WE,WORK(FI),N,WORK(FI),N) ELSE INFO = 52000 GO TO 50 END IF C COMPUTE NORM OF THE INITIAL ESTIMATES CALL DWGHT(NPP,1,WORK(SSI),NPP,1,WORK(BETACI),NPP, + WORK(WRK),NPP) IF (ISODR) THEN CALL DWGHT(N,M,WORK(TTI),IWORK(LDTTI),1,WORK(DELTAI),N, + WORK(WRK+NPP),N) WORK(PNORMI) = DNRM2(NPP+N*M,WORK(WRK),1) ELSE WORK(PNORMI) = DNRM2(NPP,WORK(WRK),1) END IF C COMPUTE SUM OF SQUARES OF THE WEIGHTED EPSILONS AND WEIGHTED DELTAS WORK(WSSEPI) = DDOT(N*NQ,WORK(FI),1,WORK(FI),1) IF (ISODR) THEN CALL DWGHT(N,M,WD,LDWD,LD2WD,WORK(DELTAI),N,WORK(WRK),N) WORK(WSSDEI) = DDOT(N*M,WORK(DELTAI),1,WORK(WRK),1) ELSE WORK(WSSDEI) = ZERO END IF WORK(WSSI) = WORK(WSSEPI) + WORK(WSSDEI) C SELECT FIRST ROW OF X + DELTA THAT CONTAINS NO ZEROS NROW = -1 CALL DSETN(N,M,WORK(XPLUSI),N,NROW) IWORK(NROWI) = NROW C SET NUMBER OF GOOD DIGITS IN FUNCTION RESULTS EPSMAC = WORK(EPSMAI) IF (NDIGIT.LT.2) THEN IWORK(NETAI) = -1 NFEV = IWORK(NFEVI) CALL DETAF(FCN, + N,M,NP,NQ, + WORK(XPLUSI),BETA,EPSMAC,NROW, + WORK(BETANI),WORK(FNI), + IFIXB,IFIXX,LDIFX, + ISTOP,NFEV,ETA,NETA, + WORK(WRK1I),WORK(WRK2I),WORK(WRK6I),WORK(WRK7I)) IWORK(ISTOPI) = ISTOP IWORK(NFEVI) = NFEV IF (ISTOP.NE.0) THEN INFO = 53000 IWORK(NETAI) = 0 WORK(ETAI) = ZERO GO TO 50 ELSE IWORK(NETAI) = -NETA WORK(ETAI) = ETA END IF ELSE IWORK(NETAI) = MIN(NDIGIT,INT(P5-LOG10(EPSMAC))) WORK(ETAI) = MAX(EPSMAC,TEN**(-NDIGIT)) END IF C CHECK DERIVATIVES IF NECESSARY IF (CHKJAC .AND. ANAJAC) THEN NTOL = -1 NFEV = IWORK(NFEVI) NJEV = IWORK(NJEVI) NETA = IWORK(NETAI) LDTT = IWORK(LDTTI) ETA = WORK(ETAI) EPSMAC = WORK(EPSMAI) CALL DJCK(FCN, + N,M,NP,NQ, + BETA,WORK(XPLUSI), + IFIXB,IFIXX,LDIFX,STPB,STPD,LDSTPD, + WORK(SSFI),WORK(TTI),LDTT, + ETA,NETA,NTOL,NROW,ISODR,EPSMAC, + WORK(FNI),WORK(FJACBI),WORK(FJACDI), + IWORK(MSGB),IWORK(MSGD),WORK(DIFFI), + ISTOP,NFEV,NJEV, + WORK(WRK1I),WORK(WRK2I),WORK(WRK6I)) IWORK(ISTOPI) = ISTOP IWORK(NFEVI) = NFEV IWORK(NJEVI) = NJEV IWORK(NTOLI) = NTOL IF (ISTOP.NE.0) THEN INFO = 54000 ELSE IF (IWORK(MSGB).NE.0 .OR. IWORK(MSGD).NE.0) THEN INFO = 40000 END IF ELSE C INDICATE USER SUPPLIED DERIVATIVES WERE NOT CHECKED IWORK(MSGB) = -1 IWORK(MSGD) = -1 END IF C PRINT APPROPRIATE ERROR MESSAGES 50 IF ((INFO.NE.0) .OR. (IWORK(MSGB).NE.-1)) THEN IF (LUNERR.NE.0 .AND. IPRINT.NE.0) THEN CALL DODPER + (INFO,LUNERR,SHORT, + N,M,NP,NQ, + LDSCLD,LDSTPD,LDWE,LD2WE,LDWD,LD2WD, + LWKMN,LIWKMN, + WORK(FJACBI),WORK(FJACDI), + WORK(DIFFI),IWORK(MSGB),ISODR,IWORK(MSGD), + WORK(XPLUSI),IWORK(NROWI),IWORK(NETAI),IWORK(NTOLI)) END IF C SET INFO TO REFLECT ERRORS IN THE USER SUPPLIED JACOBIANS IF (INFO.EQ.40000) THEN IF (IWORK(MSGB).EQ.2 .OR. IWORK(MSGD).EQ.2) THEN IF (IWORK(MSGB).EQ.2) THEN INFO = INFO + 1000 END IF IF (IWORK(MSGD).EQ.2) THEN INFO = INFO + 100 END IF ELSE INFO = 0 END IF END IF IF (INFO.NE.0) THEN RETURN END IF END IF END IF C SAVE THE INITIAL VALUES OF BETA CALL DCOPY(NP,BETA,1,WORK(BETA0I),1) C FIND LEAST SQUARES SOLUTION CALL DCOPY(N*NQ,WORK(FNI),1,WORK(FSI),1) LDTT = IWORK(LDTTI) CALL DODMN(HEAD,FSTITR,PRTPEN, + FCN, N,M,NP,NQ, JOB, BETA,Y,LDY,X,LDX, + WE,WORK(WE1I),LDWE,LD2WE,WD,LDWD,LD2WD, + IFIXB,IFIXX,LDIFX, + WORK(BETACI),WORK(BETANI),WORK(BETASI),WORK(SI), + WORK(DELTAI),WORK(DELTNI),WORK(DELTSI), + WORK(TI),WORK(FI),WORK(FNI),WORK(FSI), + WORK(FJACBI),IWORK(MSGB),WORK(FJACDI),IWORK(MSGD), + WORK(SSFI),WORK(SSI),WORK(TTI),LDTT, + STPB,STPD,LDSTPD, + WORK(XPLUSI),WORK(WRK),LWRK, + WORK,LWORK,IWORK,LIWORK,INFO) MAXIT1 = IWORK(MAXITI) - IWORK(NITERI) TSTIMP = ZERO DO 100 K=1,NP IF (BETA(K).EQ.ZERO) THEN TSTIMP = MAX(TSTIMP, + ABS(BETA(K)-WORK(BETA0I-1+K))/WORK(SSFI-1+K)) ELSE TSTIMP = MAX(TSTIMP, + ABS(BETA(K)-WORK(BETA0I-1+K))/ABS(BETA(K))) END IF 100 CONTINUE RETURN END *DODLM SUBROUTINE DODLM + (N,M,NP,NQ,NPP, + F,FJACB,FJACD, + WD,LDWD,LD2WD,SS,TT,LDTT,DELTA, + ALPHA2,TAU,EPSFCN,ISODR, + TFJACB,OMEGA,U,QRAUX,JPVT, + S,T,NLMS,RCOND,IRANK, + WRK1,WRK2,WRK3,WRK4,WRK5,WRK,LWRK,ISTOPC) C***BEGIN PROLOGUE DODLM C***REFER TO DODR,DODRC C***ROUTINES CALLED DDOT,DNRM2,DODSTP,DSCALE,DWGHT C***DATE WRITTEN 860529 (YYMMDD) C***REVISION DATE 920619 (YYMMDD) C***PURPOSE COMPUTE LEVENBERG-MARQUARDT PARAMETER AND STEPS S AND T C USING ANALOG OF THE TRUST-REGION LEVENBERG-MARQUARDT C ALGORITHM C***END PROLOGUE DODLM C...SCALAR ARGUMENTS DOUBLE PRECISION + ALPHA2,EPSFCN,RCOND,TAU INTEGER + IRANK,ISTOPC,LDTT,LDWD,LD2WD,LWRK,M,N,NLMS,NP,NPP,NQ LOGICAL + ISODR C...ARRAY ARGUMENTS DOUBLE PRECISION + DELTA(N,M),F(N,NQ),FJACB(N,NP,NQ),FJACD(N,M,NQ), + OMEGA(NQ,NQ),QRAUX(NP),S(NP),SS(NP), + T(N,M),TFJACB(N,NQ,NP),TT(LDTT,M),U(NP),WD(LDWD,LD2WD,M), + WRK(LWRK),WRK1(N,NQ,M),WRK2(N,NQ),WRK3(NP),WRK4(M,M),WRK5(M) INTEGER + JPVT(NP) C...LOCAL SCALARS DOUBLE PRECISION + ALPHA1,ALPHAN,BOT,P001,P1,PHI1,PHI2,SA,TOP,ZERO INTEGER + I,IWRK,J,K,L LOGICAL + FORVCV C...EXTERNAL FUNCTIONS DOUBLE PRECISION + DDOT,DNRM2 EXTERNAL + DDOT,DNRM2 C...EXTERNAL SUBROUTINES EXTERNAL + DODSTP,DSCALE,DWGHT C...INTRINSIC FUNCTIONS INTRINSIC + ABS,MAX,MIN,SQRT C...DATA STATEMENTS DATA + ZERO,P001,P1 + /0.0D0,0.001D0,0.1D0/ C...VARIABLE DEFINITIONS (ALPHABETICALLY) C ALPHAN: THE NEW LEVENBERG-MARQUARDT PARAMETER. C ALPHA1: THE PREVIOUS LEVENBERG-MARQUARDT PARAMETER.