SUBROUTINE n2g(n, p, x, calcr, calcj, iv, liv, lv, v) ! *** VERSION OF NL2SOL THAT CALLS RN2G *** ! i.e. it requires programmed derivatives from routine calcj USE nlsol IMPLICIT NONE ! *** PARAMETERS *** INTEGER, INTENT(IN) :: n, p REAL (dp), INTENT(IN OUT) :: x(:) INTEGER, INTENT(IN OUT) :: iv(:) INTEGER, INTENT(IN) :: liv, lv REAL (dp), INTENT(IN OUT) :: v(:) INTERFACE SUBROUTINE calcr(n, p, x, nf, r) USE nlsol IMPLICIT NONE INTEGER, INTENT(IN) :: n, p INTEGER, INTENT(IN OUT) :: nf REAL (dp), INTENT(IN) :: x(:) REAL (dp), INTENT(OUT) :: r(:) END SUBROUTINE calcr SUBROUTINE calcj(n, p, x, nf, j) USE nlsol IMPLICIT NONE INTEGER, INTENT(IN) :: n, p INTEGER, INTENT(IN OUT) :: nf REAL (dp), INTENT(IN) :: x(:) REAL (dp), INTENT(OUT) :: j(:,:) ! j(n,p) END SUBROUTINE calcj SUBROUTINE rn2g(d, dr, iv, liv, lv, n, nd, n1, n2, p, r, rd, v, x) USE nlsol IMPLICIT NONE REAL (dp), INTENT(IN OUT) :: d(:) REAL (dp), INTENT(IN OUT) :: dr(:,:) ! dr(nd,p) INTEGER, INTENT(IN OUT) :: iv(:) INTEGER, INTENT(IN) :: liv INTEGER, INTENT(IN) :: lv INTEGER, INTENT(IN) :: n INTEGER, INTENT(IN) :: nd INTEGER, INTENT(IN OUT) :: n1 INTEGER, INTENT(IN OUT) :: n2 INTEGER, INTENT(IN) :: p REAL (dp), INTENT(IN OUT) :: r(:) REAL (dp), INTENT(IN OUT) :: rd(:) REAL (dp), INTENT(OUT) :: v(:) REAL (dp), INTENT(IN OUT) :: x(:) END SUBROUTINE rn2g END INTERFACE ! *** PARAMETER USAGE *** ! N....... TOTAL NUMBER OF RESIDUALS. ! P....... NUMBER OF PARAMETERS (COMPONENTS OF X) BEING ESTIMATED. ! X....... PARAMETER VECTOR BEING ESTIMATED (INPUT = INITIAL GUESS, ! OUTPUT = BEST VALUE FOUND). ! CALCR... SUBROUTINE FOR COMPUTING RESIDUAL VECTOR. ! CALCJ... SUBROUTINE FOR COMPUTING JACOBIAN MATRIX = MATRIX OF FIRST ! PARTIALS OF THE RESIDUAL VECTOR. ! IV...... INTEGER VALUES ARRAY. ! LIV..... LENGTH OF IV (SEE DISCUSSION BELOW). ! LV...... LENGTH OF V (SEE DISCUSSION BELOW). ! V....... FLOATING-POINT VALUES ARRAY. ! *** DISCUSSION *** ! NOTE... NL2SOL (MENTIONED BELOW) IS A CODE FOR SOLVING NONLINEAR ! LEAST-SQUARES PROBLEMS. IT IS DESCRIBED IN ACM TRANS. MATH. SOFTWARE, ! VOL. 7 (1981), PP. 369-383 (AN ADAPTIVE NONLINEAR LEAST-SQUARES ALGORITHM, ! BY J.E. DENNIS, D.M. GAY, AND R.E. WELSCH). ! LIV GIVES THE LENGTH OF IV. IT MUST BE AT LEAST 82+P. IF NOT, THEN ! N2G RETURNS WITH IV(1) = 15. WHEN N2G RETURNS, THE MINIMUM ACCEPTABLE ! VALUE OF LIV IS STORED IN IV(LASTIV) = IV(44), (PROVIDED THAT LIV >= 44). ! LV GIVES THE LENGTH OF V. THE MINIMUM VALUE FOR LV IS ! LV0 = 105 + P*(N + 2*P + 17) + 2*N. IF LV IS SMALLER THAN THIS, THEN ! N2G RETURNS WITH IV(1) = 16. WHEN N2G RETURNS, THE MINIMUM ACCEPTABLE ! VALUE OF LV IS STORED IN IV(LASTV) = IV(45) (PROVIDED LIV >= 45). ! RETURN CODES AND CONVERGENCE TOLERANCES ARE THE SAME AS FOR ! NL2SOL, WITH SOME SMALL EXTENSIONS... IV(1) = 15 MEANS LIV WAS ! TOO SMALL. IV(1) = 16 MEANS LV WAS TOO SMALL. ! THERE ARE TWO NEW V INPUT COMPONENTS... V(LMAXS) = V(36) AND ! V(SCTOL) = V(37) SERVE WHERE V(LMAX0) AND V(RFCTOL) FORMERLY DID ! IN THE SINGULAR CONVERGENCE TEST -- SEE THE NL2SOL DOCUMENTATION. ! *** DEFAULT VALUES *** ! DEFAULT VALUES ARE PROVIDED BY SUBROUTINE IVSET, RATHER THAN ! DFAULT. THE CALLING SEQUENCE IS... ! CALL IVSET(1, IV, LIV, LV, V) ! THE FIRST PARAMETER IS AN INTEGER 1. IF LIV AND LV ARE LARGE ENOUGH FOR ! IVSET, THEN IVSET SETS IV(1) TO 12. OTHERWISE IT SETS IV(1) TO 15 OR 16. ! CALLING N2G WITH IV(1) = 0 CAUSES ALL DEFAULT VALUES TO BE USED FOR THE ! INPUT COMPONENTS OF IV AND V. ! IF YOU FIRST CALL IVSET, THEN SET IV(1) TO 13 AND CALL N2G, THEN ! STORAGE ALLOCATION ONLY WILL BE PERFORMED. IN PARTICULAR, IV(D) = IV(27), ! IV(J) = IV(70), AND IV(R) = IV(61) WILL BE SET TO THE FIRST SUBSCRIPT IN ! V OF THE SCALE VECTOR, THE JACOBIAN MATRIX, AND THE RESIDUAL VECTOR ! RESPECTIVELY, PROVIDED LIV AND LV ARE LARGE ENOUGH. IF SO, THEN N2G ! RETURNS WITH IV(1) = 14. WHEN CALLED WITH IV(1) = 14, N2G ASSUMES THAT ! STORAGE HAS BEEN ALLOCATED, AND IT BEGINS THE MINIMIZATION ALGORITHM. ! *** SCALE VECTOR *** ! ONE DIFFERENCE WITH NL2SOL IS THAT THE SCALE VECTOR D IS STORED IN V, ! STARTING AT A DIFFERENT SUBSCRIPT. THE STARTING SUBSCRIPT VALUE IS STILL ! STORED IN IV(D) = IV(27). THE DISCUSSION OF DEFAULT VALUES ABOVE TELLS ! HOW TO HAVE IV(D) SET BEFORE THE ALGORITHM IS STARTED. ! *** REGRESSION DIAGNOSTICS *** ! IF IV(RDREQ) SO DICTATES, THEN ESTIMATES ARE COMPUTED OF THE INFLUENCE ! EACH RESIDUAL COMPONENT HAS ON THE FINAL PARAMETER ESTIMATE X. ! THE GENERAL IDEA IS THAT ONE MAY WISH TO EXAMINE RESIDUAL COMPONENTS (AND ! THE DATA BEHIND THEM) FOR WHICH THE INFLUENCE ESTIMATE IS SIGNIFICANTLY ! LARGER THAN MOST OF THE OTHER INFLUENCE ESTIMATES. THESE ESTIMATES, ! HEREAFTER CALLED REGRESSION DIAGNOSTICS, ARE ONLY COMPUTED IF ! IV(RDREQ) = 2 OR 3. IN THIS CASE, FOR I = 1(1)N, ! SQRT( G(I)**T * H(I)**-1 * G(I) ) ! IS COMPUTED AND STORED IN V, STARTING AT V(IV(REGD)), WHERE RDREQ = 57 ! AND REGD = 67. HERE G(I) STANDS FOR THE GRADIENT RESULTING WHEN THE I-TH ! OBSERVATION IS DELETED AND H(I) STANDS FOR AN APPROXIMATION TO THE ! CORRESPONDING HESSIAN AT X, THE SOLUTION CORRESPONDING TO ALL OBSERVATIONS. ! (THIS APPROXIMATION IS OBTAINED BY SUBTRACTING THE FIRST-ORDER ! CONTRIBUTION OF THE I-TH OBSERVATION TO THE HESSIAN FROM A FINITE- ! DIFFERENCE HESSIAN APPROXIMATION. IF H IS INDEFINITE, THEN IV(REGD) IS ! SET TO -1. IF H(I) IS INDEFINITE, THEN -1 IS RETURNED AS THE DIAGNOSTIC ! FOR OBSERVATION I. IF NO DIAGNOSTICS ARE COMPUTED, PERHAPS BECAUSE ! OF A FAILURE TO CONVERGE, THEN IV(REGD) = 0 IS RETURNED.) ! PRINTING OF THE REGRESSION DIAGNOSTICS IS CONTROLLED BY ! IV(COVPRT) = IV(14)... IF IV(COVPRT) = 3, THEN BOTH THE COVARIANCE MATRIX ! AND THE REGRESSION DIAGNOSTICS ARE PRINTED. IV(COVPRT) = 2 CAUSES ONLY ! THE REGRESSION DIAGNOSTICS TO BE PRINTED, IV(COVPRT) = 1 CAUSES ONLY THE ! COVARIANCE MATRIX TO BE PRINTED, AND IV(COVPRT) = 0 CAUSES NEITHER TO BE ! PRINTED. ! RDREQ = 57 AND REGD = 67. ! *** GENERAL *** ! CODED BY DAVID M. GAY. !+++++++++++++++++++++++++++ DECLARATIONS +++++++++++++++++++++++++++ ! *** EXTERNAL SUBROUTINES *** ! EXTERNAL ivset, rn2g, n2rdp ! IVSET.... PROVIDES DEFAULT IV AND V INPUT COMPONENTS. ! RN2G... CARRIES OUT OPTIMIZATION ITERATIONS. ! N2RDP... PRINTS REGRESSION DIAGNOSTICS. ! *** NO INTRINSIC FUNCTIONS *** ! *** LOCAL VARIABLES *** INTEGER :: d1, iv1, n1, n2, nf, r1, rd1 REAL (dp), ALLOCATABLE :: dr(:,:) ! *** IV COMPONENTS *** INTEGER, PARAMETER :: d=27, j=70, r=61, regd0=82 !--------------------------------- BODY ------------------------------ IF (iv(1) == 0) CALL ivset(1, iv, liv, lv, v) iv1 = iv(1) IF (iv1 == 14) GO TO 10 IF (iv1 > 2 .AND. iv1 < 12) GO TO 10 IF (iv1 == 12) iv(1) = 13 IF (iv(1) == 13) iv(vneed) = iv(vneed) + p + n*(p+2) ALLOCATE( dr(n,p) ) CALL rn2g(x, dr, iv, liv, lv, n, n, n1, n2, p, v, v, v, x) IF (iv(1) /= 14) GO TO 999 ! *** STORAGE ALLOCATION *** iv(d) = iv(nextv) iv(r) = iv(d) + p iv(regd0) = iv(r) + n iv(j) = iv(regd0) + n iv(nextv) = iv(j) + n*p IF (iv1 == 13) GO TO 999 10 d1 = iv(d) r1 = iv(r) rd1 = iv(regd0) 20 CALL rn2g(v(d1:), dr, iv, liv, lv, n, n, n1, n2, p, v(r1:), v(rd1:), v, x) IF (iv(1) == 2) THEN GO TO 50 ELSE IF (iv(1) > 2) THEN GO TO 60 END IF ! *** NEW FUNCTION VALUE (R VALUE) NEEDED *** nf = iv(nfcall) CALL calcr(n, p, x, nf, v(r1:)) IF (nf > 0) GO TO 40 iv(toobig) = 1 GO TO 20 40 IF (iv(1) > 0) GO TO 20 ! *** COMPUTE DR = GRADIENT OF R COMPONENTS *** 50 CALL calcj(n, p, x, iv(nfgcal), dr) IF (iv(nfgcal) == 0) iv(toobig) = 1 GO TO 20 ! *** INDICATE WHETHER THE REGRESSION DIAGNOSTIC ARRAY WAS COMPUTED ! *** AND PRINT IT IF SO REQUESTED... 60 IF (iv(regd) > 0) iv(regd) = rd1 CALL n2rdp(iv, v(rd1:), n, v) DEALLOCATE( dr ) 999 RETURN ! *** LAST LINE OF N2G FOLLOWS *** END SUBROUTINE n2g