SUBROUTINE n2fb(n, p, x, b, calcr, iv, liv, lv, v) ! *** MINIMIZE A NONLINEAR SUM OF SQUARES USING RESIDUAL VALUES ONLY.. ! *** THIS AMOUNTS TO N2G WITHOUT THE SUBROUTINE PARAMETER CALCJ. USE nlsol IMPLICIT NONE ! *** PARAMETERS *** INTEGER, INTENT(IN) :: n INTEGER, INTENT(IN) :: p REAL (dp), INTENT(IN OUT) :: x(:) REAL (dp), INTENT(IN) :: b(:,:) ! b(2,p) INTEGER, INTENT(IN OUT) :: iv(:) INTEGER, INTENT(IN) :: liv INTEGER, INTENT(IN) :: 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 rn2gb(b, d, dr, iv, liv, lv, n, nd, n1, n2, p, r, rd, v, x) USE nlsol IMPLICIT NONE REAL (dp), INTENT(IN) :: b(:,:) ! b(2,p) 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 rn2gb END INTERFACE !----------------------------- DISCUSSION ---------------------------- ! THIS AMOUNTS TO SUBROUTINE NL2SNO (REF. 1) MODIFIED TO HANDLE ! SIMPLE BOUNDS ON THE VARIABLES... ! B(1,I) .LE. X(I) .LE. B(2,I), I = 1(1)P. ! THE PARAMETERS FOR N2FB ARE THE SAME AS THOSE FOR N2GB ! (WHICH SEE), EXCEPT THAT CALCJ IS OMITTED. INSTEAD OF CALLING ! CALCJ TO OBTAIN THE JACOBIAN MATRIX OF R AT X, N2FB COMPUTES ! AN APPROXIMATION TO IT BY FINITE (FORWARD) DIFFERENCES -- SEE ! V(DLTFDJ) BELOW. N2FB DOES NOT COMPUTE A COVARIANCE MATRIX. ! THE NUMBER OF EXTRA CALLS ON CALCR USED IN COMPUTING THE JACO- ! BIAN APPROXIMATION ARE NOT INCLUDED IN THE FUNCTION EVALUATION ! COUNT IV(NFCALL), BUT ARE RECORDED IN IV(NGCALL) INSTEAD. ! V(DLTFDJ)... V(43) HELPS CHOOSE THE STEP SIZE USED WHEN COMPUTING THE ! FINITE-DIFFERENCE JACOBIAN MATRIX. FOR DIFFERENCES IN- ! VOLVING X(I), THE STEP SIZE FIRST TRIED IS ! V(DLTFDJ) * MAX(ABS(X(I)), 1/D(I)), ! WHERE D IS THE CURRENT SCALE VECTOR (SEE REF. 1). (IF ! THIS STEP IS TOO BIG, I.E., IF CALCR SETS NF TO 0, THEN ! SMALLER STEPS ARE TRIED UNTIL THE STEP SIZE IS SHRUNK BE- ! LOW 1000 * MACHEP, WHERE MACHEP IS THE UNIT ROUNDOFF. ! DEFAULT = MACHEP**0.5. ! *** REFERENCE *** ! 1. DENNIS, J.E., GAY, D.M., AND WELSCH, R.E. (1981), AN ADAPTIVE ! NONLINEAR LEAST-SQUARES ALGORITHM, ACM TRANS. MATH. ! SOFTWARE, VOL. 7, NO. 3. ! *** GENERAL *** ! CODED BY DAVID M. GAY. !+++++++++++++++++++++++++++ DECLARATIONS +++++++++++++++++++++++++++ ! *** EXTERNAL SUBROUTINES *** ! EXTERNAL ivset, rn2gb, v7scp ! IVSET.... PROVIDES DEFAULT IV AND V INPUT COMPONENTS. ! RN2GB... CARRIES OUT OPTIMIZATION ITERATIONS. ! N2RDP... PRINTS REGRESSION DIAGNOSTICS. ! V7SCP... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR. ! *** LOCAL VARIABLES *** INTEGER :: d1, dk, dr1, i, iv1, j1k, k, n1, n2, nf, ng, rd1, r1, rn REAL (dp) :: h1, h0, t, xk, xk1 REAL (dp), ALLOCATABLE :: dr(:,:) ! *** IV AND V COMPONENTS *** INTEGER, PARAMETER :: d=27, j=70, r=61, regd0=82 REAL (dp), PARAMETER :: hlim=0.1_dp, negpt5=-0.5_dp !--------------------------------- BODY ------------------------------ IF (iv(1) == 0) CALL ivset(1, iv, liv, lv, v) iv(covreq) = 0 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 rn2gb(b, 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) dr1 = iv(j) r1 = iv(r) rn = r1 + n - 1 rd1 = iv(regd0) 20 CALL rn2gb(b, 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 999 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 FINITE-DIFFERENCE APPROXIMATION TO DR = GRAD. OF R *** ! *** INITIALIZE D IF NECESSARY *** 50 IF (iv(mode) < 0 .AND. v(dinit) == zero) v(d1:d1+p-1) = one j1k = dr1 dk = d1 ng = iv(ngcall) - 1 IF (iv(1) == -1) iv(ngcov) = iv(ngcov) - 1 DO k = 1, p IF (b(1,k) >= b(2,k)) GO TO 110 xk = x(k) h1 = v(dltfdj) * MAX( ABS(xk), one/v(dk)) h0 = h1 dk = dk + 1 t = negpt5 xk1 = xk + h1 IF (xk - h1 >= b(1,k)) GO TO 60 t = -t IF (xk1 > b(2,k)) GO TO 80 60 IF (xk1 <= b(2,k)) GO TO 70 t = -t h1 = -h1 xk1 = xk + h1 IF (xk1 < b(1,k)) GO TO 80 70 x(k) = xk1 nf = iv(nfgcal) CALL calcr (n, p, x, nf, v(j1k:)) ng = ng + 1 IF (nf > 0) GO TO 90 h1 = t * h1 xk1 = xk + h1 IF ( ABS(h1/h0) >= hlim) GO TO 70 80 iv(toobig) = 1 iv(ngcall) = ng GO TO 20 90 x(k) = xk iv(ngcall) = ng DO i = r1, rn dr(i-r1+1,k) = (v(j1k) - v(i)) / h1 j1k = j1k + 1 END DO CYCLE ! *** SUPPLY A ZERO DERIVATIVE FOR CONSTANT COMPONENTS... 110 v(j1k:j1k+n-1) = zero j1k = j1k + n END DO GO TO 20 999 DEALLOCATE( dr ) RETURN ! *** LAST CARD OF N2FB FOLLOWS *** END SUBROUTINE n2fb