www.pudn.com > lm.zip > LGSTIC4.F90

MODULE data_storage 
! The data is stored in this module. 
! This is equivalent to a COMMON area in old Fortran. 
 
IMPLICIT NONE 
INTEGER, PARAMETER, PRIVATE  :: dp = SELECTED_REAL_KIND(12, 60) 
 
REAL (dp), SAVE  :: x(100), y(100) 
 
END MODULE data_storage 
 
 
 
PROGRAM logistic4 
! Fit the 4-parameter logistic curve: 
 
!     Y = A + C / [1 + exp{-B(X - D)}] 
 
! by unweighted least squares. 
 
! Latest revision - 25 October 2001 
 
USE Levenberg_Marquardt 
USE data_storage 
IMPLICIT NONE 
 
INTERFACE 
  SUBROUTINE fcn(m, n, x, fvec, fjac, iflag) 
    IMPLICIT NONE 
    INTEGER, PARAMETER         :: dp = SELECTED_REAL_KIND(12, 60) 
    INTEGER, INTENT(IN)        :: m, n 
    REAL (dp), INTENT(IN)      :: x(:) 
    REAL (dp), INTENT(IN OUT)  :: fvec(:) 
    REAL (dp), INTENT(OUT)     :: fjac(:,:) 
    INTEGER, INTENT(IN OUT)    :: iflag 
  END SUBROUTINE fcn 
END INTERFACE 
 
REAL (dp), ALLOCATABLE  :: fvec(:), fjac(:,:) 
INTEGER, PARAMETER      :: n = 4                 ! The number of parameters 
INTEGER                 :: info, iostatus, ipvt(n), m 
REAL (dp)               :: p(n), tol = 1.0E-04_dp 
 
! Read in the data, counting the number of cases. 
 
OPEN(UNIT=8, FILE='lgstic4.dat', STATUS='OLD') 
m = 1 
DO 
  READ(8, *, IOSTAT=iostatus) x(m), y(m) 
  IF (iostatus < 0) EXIT 
  IF (iostatus > 0) CYCLE    ! Skips cases with errors - e.g. a heading 
  m = m + 1 
END DO 
m = m - 1 
ALLOCATE( fvec(m), fjac(m,n) ) 
 
! Set starting values for parameters. 
! A = p(1), B = p(2), C = p(3), D = p(4) 
 
p(1) = MINVAL( y(1:m) ) 
p(3) = MAXVAL( y(1:m) ) - p(1) 
p(2) = 2.0_dp / ( MAXVAL( x(1:m) ) - MINVAL( x(1:m) ) ) 
p(4) = 0.5_dp * ( MAXVAL( x(1:m) ) + MINVAL( x(1:m) ) ) 
 
CALL lmder1(fcn, m, n, p, fvec, fjac, tol, info, ipvt) 
 
SELECT CASE (info) 
  CASE (:-1) 
    WRITE(*, *) 'Users FCN returned INFO = ', -info 
  CASE (0) 
    WRITE(*, *) 'Improper values for input parameters' 
  CASE (1:3) 
    WRITE(*, *) 'Convergence' 
    WRITE(*, '(a, 4g13.5)') ' Final values of A, B, C, D: ', p 
  CASE (4) 
    WRITE(*, *) 'Residuals orthogonal to the Jacobian' 
    WRITE(*, *) 'There may be an error in FCN' 
    WRITE(*, '(a, 4g13.5)') ' Final values of A, B, C, D: ', p 
  CASE (5) 
    WRITE(*, *) 'Too many calls to FCN' 
    WRITE(*, *) 'Either slow convergence, or an error in FCN' 
    WRITE(*, '(a, 4g13.5)') ' Final values of A, B, C, D: ', p 
  CASE (6:7) 
    WRITE(*, *) 'TOL was set too small' 
    WRITE(*, '(a, 4g13.5)') ' Final values of A, B, C, D: ', p 
  CASE DEFAULT 
    WRITE(*, *) 'INFO =', info, ' ???' 
END SELECT 
 
STOP 
END PROGRAM logistic4 
 
 
 
SUBROUTINE fcn(m, n, p, fvec, fjac, iflag) 
! Calculate either residuals or the Jacobian matrix. 
! A = p(1), B = p(2), C = p(3), D = p(4) 
! m = no. of cases, n = no. of parameters (4) 
 
USE data_storage 
IMPLICIT NONE 
INTEGER, PARAMETER         :: dp = SELECTED_REAL_KIND(12, 60) 
INTEGER, INTENT(IN)        :: m, n 
REAL (dp), INTENT(IN)      :: p(:) 
REAL (dp), INTENT(IN OUT)  :: fvec(:) 
REAL (dp), INTENT(OUT)     :: fjac(:,:) 
INTEGER, INTENT(IN OUT)    :: iflag 
 
! Local variables 
 
REAL (dp), PARAMETER  :: one = 1.0_dp 
INTEGER               :: i 
REAL (dp)             :: expntl, temp 
 
IF (iflag == 1) THEN 
  fvec = y(1:m) - p(1) - p(3)/(one + EXP(-p(2)*(x(1:m) - p(4)))) 
ELSE IF (iflag == 2) THEN 
  fjac(1:m,1) = -one 
  DO i = 1, m 
    expntl = EXP(-p(2)*(x(i) - p(4))) 
    temp = one / (one + expntl) 
    fjac(i,2) = - p(3) * temp**2 * (x(i) - p(4)) * expntl 
    fjac(i,3) = - temp 
    fjac(i,4) =   p(3) * temp**2 * expntl * p(2) 
  END DO 
END IF 
 
RETURN 
END SUBROUTINE fcn