MODULE common_refnum
IMPLICIT NONE

! COMMON /refnum/ nprob,nfev,njev
INTEGER, SAVE :: nprob, nfev, njev

END MODULE common_refnum



PROGRAM test_lmder
 
! Code converted using TO_F90 by Alan Miller
! Date: 1999-12-09  Time: 23:25:02
 
!  **********

!  THIS PROGRAM TESTS CODES FOR THE LEAST-SQUARES SOLUTION OF
!  M NONLINEAR EQUATIONS IN N VARIABLES. IT CONSISTS OF A DRIVER
!  AND AN INTERFACE SUBROUTINE FCN. THE DRIVER READS IN DATA,
!  CALLS THE NONLINEAR LEAST-SQUARES SOLVER, AND FINALLY PRINTS
!  OUT INFORMATION ON THE PERFORMANCE OF THE SOLVER. THIS IS
!  ONLY A SAMPLE DRIVER, MANY OTHER DRIVERS ARE POSSIBLE. THE
!  INTERFACE SUBROUTINE FCN IS NECESSARY TO TAKE INTO ACCOUNT THE
!  FORMS OF CALLING SEQUENCES USED BY THE FUNCTION AND JACOBIAN
!  SUBROUTINES IN THE VARIOUS NONLINEAR LEAST-SQUARES SOLVERS.

!  SUBPROGRAMS CALLED

!    USER-SUPPLIED ...... FCN

!    MINPACK-SUPPLIED ... DPMPAR,ENORM,INITPT,LMDER1,SSQFCN

!    FORTRAN-SUPPLIED ... SQRT

!  ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980.
!  BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE

!  **********

USE Levenberg_Marquardt
USE common_refnum
IMPLICIT NONE
INTEGER   :: i, ic, info, k, m, n, ntries
INTEGER   :: iwa(40), ma(60), na(60), nf(60), nj(60), np(60), nx(60)
REAL (dp) :: factor, fnorm1, fnorm2, tol
REAL (dp) :: fjac(65,40), fnm(60), fvec(65), x(40)

!     LOGICAL INPUT UNIT IS ASSUMED TO BE NUMBER 5.
!     LOGICAL OUTPUT UNIT IS ASSUMED TO BE NUMBER 6.

INTEGER, PARAMETER :: nread = 5, nwrite = 6

! EXTERNAL fcn

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

  SUBROUTINE ssqfcn (m, n, x, fvec, nprob)
    USE Levenberg_Marquardt
    IMPLICIT NONE
    INTEGER, INTENT(IN)     :: m, n
    REAL (dp), INTENT(IN)   :: x(:)
    REAL (dp), INTENT(OUT)  :: fvec(:)
    INTEGER, INTENT(IN)     :: nprob
  END SUBROUTINE ssqfcn

  SUBROUTINE ssqjac (m, n, x, fjac, nprob)
    USE Levenberg_Marquardt
    IMPLICIT NONE
    INTEGER, INTENT(IN)     :: m, n, nprob
    REAL (dp), INTENT(IN)   :: x(:)
    REAL (dp), INTENT(OUT)  :: fjac(:,:)
  END SUBROUTINE ssqjac
END INTERFACE

REAL (dp), PARAMETER :: one = 1.0_dp, ten = 10.0_dp

tol = SQRT( EPSILON(one) )
ic = 0
10 READ (nread,50) nprob, n, m, ntries
IF (nprob <= 0) GO TO 30
factor = one
DO  k=1,ntries
  ic = ic + 1
  CALL initpt (n, x, nprob, factor)
  CALL ssqfcn (m, n, x, fvec, nprob)
  fnorm1 = enorm(m, fvec)
  WRITE (nwrite,60) nprob, n, m
  nfev = 0
  njev = 0
  CALL lmder1 (fcn, m, n, x, fvec, fjac, tol, info, iwa)
  CALL ssqfcn (m, n, x, fvec, nprob)
  fnorm2 = enorm(m, fvec)
  np(ic) = nprob
  na(ic) = n
  ma(ic) = m
  nf(ic) = nfev
  nj(ic) = njev
  nx(ic) = info
  fnm(ic) = fnorm2
  WRITE (nwrite,70) fnorm1, fnorm2, nfev, njev, info, x(1:n)
  factor = ten*factor
END DO
GO TO 10

30 WRITE (nwrite,80) ic
WRITE (nwrite,90)
DO  i=1,ic
  WRITE (nwrite,100) np(i), na(i), ma(i), nf(i), nj(i), nx(i), fnm(i)
END DO
STOP

50 FORMAT (4I5)
60 FORMAT (////'      PROBLEM', i5, '      DIMENSIONS', 2I5//)
70 FORMAT ('      INITIAL L2 NORM OF THE RESIDUALS', g15.7//  &
           '      FINAL L2 NORM OF THE RESIDUALS  ', g15.7//  &
           '      NUMBER OF FUNCTION EVALUATIONS  ', i10//  &
           '      NUMBER OF JACOBIAN EVALUATIONS  ', i10//  &
           '      EXIT PARAMETER', t39, i10//  &
           '      FINAL APPROXIMATE SOLUTION'// (t6, 5g15.7))
80 FORMAT (' SUMMARY OF ', i3, ' CALLS TO LMDER1'/)
90 FORMAT (' NPROB   N    M   NFEV  NJEV  INFO  FINAL L2 NORM'/)
100 FORMAT (3I5, 3I6, ' ', g15.7)


CONTAINS


SUBROUTINE initpt (n, x, nprob, factor)
INTEGER, INTENT(IN)     :: n, nprob
REAL (dp), INTENT(IN)   :: factor
REAL (dp), INTENT(OUT)  :: x(:)
!     **********

!     SUBROUTINE INITPT

!     THIS SUBROUTINE SPECIFIES THE STANDARD STARTING POINTS FOR THE
!     FUNCTIONS DEFINED BY SUBROUTINE SSQFCN. THE SUBROUTINE RETURNS
!     IN X A MULTIPLE (FACTOR) OF THE STANDARD STARTING POINT. FOR
!     THE 11TH FUNCTION THE STANDARD STARTING POINT IS ZERO, SO IN
!     THIS CASE, IF FACTOR IS NOT UNITY, THEN THE SUBROUTINE RETURNS
!     THE VECTOR  X(J) = FACTOR, J=1,...,N.

!     THE SUBROUTINE STATEMENT IS

!       SUBROUTINE INITPT(N,X,NPROB,FACTOR)

!     WHERE

!       N IS A POSITIVE INTEGER INPUT VARIABLE.

!       X IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS THE STANDARD
!         STARTING POINT FOR PROBLEM NPROB MULTIPLIED BY FACTOR.

!       NPROB IS A POSITIVE INTEGER INPUT VARIABLE WHICH DEFINES THE
!         NUMBER OF THE PROBLEM. NPROB MUST NOT EXCEED 18.

!       FACTOR IS AN INPUT VARIABLE WHICH SPECIFIES THE MULTIPLE OF
!         THE STANDARD STARTING POINT.  IF FACTOR IS UNITY, NO
!         MULTIPLICATION IS PERFORMED.

!     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980.
!     BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE

!     **********
INTEGER   :: j
REAL (dp) :: h
REAL (dp), PARAMETER :: zero = 0.0_dp, half = 0.5_dp, one = 1.0_dp,  &
                        two = 2.0_dp, three = 3.0_dp, five = 5.0_dp, &
                        seven = 7.0_dp, ten = 10.0_dp, twenty = 20.0_dp, &
                        twntf = 25.0_dp
REAL (dp), PARAMETER :: c1 = 1.2_dp, c2 = 0.25_dp, c3 = 0.39_dp,    &
                        c4 = 0.415_dp, c5 = 0.02_dp, c6 = 4000._dp, &
                        c7 = 250._dp, c8 = 0.3_dp, c9 = 0.4_dp,     &
                        c10 = 1.5_dp, c11 = 0.01_dp, c12 = 1.3_dp,  &
                        c13 = 0.65_dp, c14 = 0.7_dp, c15 = 0.6_dp,  &
                        c16 = 4.5_dp, c17 = 5.5_dp

!     SELECTION OF INITIAL POINT.

SELECT CASE ( nprob )
  CASE (    1:3)   !     LINEAR FUNCTION - FULL RANK OR RANK 1.

    x(1:n) = one

  CASE (    4)     !     ROSENBROCK FUNCTION.

    x(1) = -c1
    x(2) = one

  CASE (    5)     !     HELICAL VALLEY FUNCTION.

    x(1) = -one
    x(2) = zero
    x(3) = zero

  CASE (    6)     !     POWELL SINGULAR FUNCTION.

    x(1) = three
    x(2) = -one
    x(3) = zero
    x(4) = one

  CASE (    7)     !     FREUDENSTEIN AND ROTH FUNCTION.

x(1) = half
x(2) = -two

  CASE (    8)     !     BARD FUNCTION.

x(1) = one
x(2) = one
x(3) = one

  CASE (    9)     !     KOWALIK AND OSBORNE FUNCTION.

x(1) = c2
x(2) = c3
x(3) = c4
x(4) = c3

  CASE (   10)     !     MEYER FUNCTION.

    x(1) = c5
    x(2) = c6
    x(3) = c7

  CASE (   11)     !     WATSON FUNCTION.

    x(1:n) = zero

  CASE (   12)     !     BOX 3-DIMENSIONAL FUNCTION.

    x(1) = zero
    x(2) = ten
    x(3) = twenty

  CASE (   13)     !     JENNRICH AND SAMPSON FUNCTION.

    x(1) = c8
    x(2) = c9

  CASE (   14)     !     BROWN AND DENNIS FUNCTION.

    x(1) = twntf
    x(2) = five
    x(3) = -five
    x(4) = -one

  CASE (   15)     !     CHEBYQUAD FUNCTION.

    h = one / DBLE(n+1)
    DO  j=1,n
      x(j) = DBLE(j)*h
    END DO

  CASE (   16)     !     BROWN ALMOST-LINEAR FUNCTION.

    x(1:n) = half

  CASE (   17)     !     OSBORNE 1 FUNCTION.

    x(1) = half
    x(2) = c10
    x(3) = -one
    x(4) = c11
    x(5) = c5

  CASE (   18)     !     OSBORNE 2 FUNCTION.

    x(1) = c12
    x(2) = c13
    x(3) = c13
    x(4) = c14
    x(5) = c15
    x(6) = three
    x(7) = five
    x(8) = seven
    x(9) = two
    x(10) = c16
    x(11) = c17

END SELECT

!     COMPUTE MULTIPLE OF INITIAL POINT.

IF (factor == one) GO TO 250
IF (nprob == 11) GO TO 230
x(1:n) = factor*x(1:n)
GO TO 250

230 x(1:n) = factor

250 RETURN

!     LAST CARD OF SUBROUTINE INITPT.

END SUBROUTINE initpt

!     LAST CARD OF DRIVER.

END PROGRAM test_lmder



SUBROUTINE fcn (m, n, x, fvec, fjac, iflag)
USE common_refnum
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

INTERFACE
  SUBROUTINE ssqfcn (m, n, x, fvec, nprob)
    USE Levenberg_Marquardt
    IMPLICIT NONE
    INTEGER, INTENT(IN)     :: m, n
    REAL (dp), INTENT(IN)   :: x(:)
    REAL (dp), INTENT(OUT)  :: fvec(:)
    INTEGER, INTENT(IN)     :: nprob
  END SUBROUTINE ssqfcn

  SUBROUTINE ssqjac (m, n, x, fjac, nprob)
    USE Levenberg_Marquardt
    IMPLICIT NONE
    INTEGER, INTENT(IN)     :: m, n, nprob
    REAL (dp), INTENT(IN)   :: x(:)
    REAL (dp), INTENT(OUT)  :: fjac(:,:)
  END SUBROUTINE ssqjac
END INTERFACE

!  **********

!  THE CALLING SEQUENCE OF FCN SHOULD BE IDENTICAL TO THE
!  CALLING SEQUENCE OF THE FUNCTION SUBROUTINE IN THE NONLINEAR
!  LEAST-SQUARES SOLVER. FCN SHOULD ONLY CALL THE TESTING
!  FUNCTION AND JACOBIAN SUBROUTINES SSQFCN AND SSQJAC WITH
!  THE APPROPRIATE VALUE OF PROBLEM NUMBER (NPROB).

!  SUBPROGRAMS CALLED

!    MINPACK-SUPPLIED ... SSQFCN,SSQJAC

!  ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980.
!  BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE

!  **********

IF (iflag == 1) CALL ssqfcn (m, n, x, fvec, nprob)
IF (iflag == 2) CALL ssqjac (m, n, x, fjac, nprob)
IF (iflag == 1) nfev = nfev + 1
IF (iflag == 2) njev = njev + 1
RETURN

!     LAST CARD OF INTERFACE SUBROUTINE FCN.

END SUBROUTINE fcn



SUBROUTINE ssqjac (m, n, x, fjac, nprob)
USE Levenberg_Marquardt
IMPLICIT NONE
INTEGER, INTENT(IN)     :: m, n, nprob
REAL (dp), INTENT(IN)   :: x(:)
REAL (dp), INTENT(OUT)  :: fjac(:,:)

!  **********

!  SUBROUTINE SSQJAC

!  THIS SUBROUTINE DEFINES THE JACOBIAN MATRICES OF EIGHTEEN
!  NONLINEAR LEAST SQUARES PROBLEMS. THE PROBLEM DIMENSIONS ARE
!  AS DESCRIBED IN THE PROLOGUE COMMENTS OF SSQFCN.

!  THE SUBROUTINE STATEMENT IS

!    SUBROUTINE SSQJAC(M,N,X,FJAC,LDFJAC,NPROB)

!  WHERE

!    M AND N ARE POSITIVE INTEGER INPUT VARIABLES. N MUST NOT EXCEED M.

!    X IS AN INPUT ARRAY OF LENGTH N.

!    FJAC IS AN M BY N OUTPUT ARRAY WHICH CONTAINS THE JACOBIAN
!      MATRIX OF THE NPROB FUNCTION EVALUATED AT X.

!    LDFJAC IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN M
!      WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY FJAC.

!    NPROB IS A POSITIVE INTEGER VARIABLE WHICH DEFINES THE
!      NUMBER OF THE PROBLEM. NPROB MUST NOT EXCEED 18.

!  SUBPROGRAMS CALLED

!    FORTRAN-SUPPLIED ... DATAN,DCOS,DEXP,DSIN,SQRT

!  ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980.
!  BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE

!  **********

INTEGER   :: i, j, k, mm1, nm1
REAL (dp) :: div, dx, prod, s2, temp, ti, tmp1, tmp2, tmp3, tmp4, tpi
REAL (dp), PARAMETER :: v(11) = (/ 4.0D0, 2.0D0, 1.0D0, 5.0D-1, 2.5D-1, &
                                   1.67D-1, 1.25D-1, 1.0D-1, 8.33D-2,  &
                                   7.14D-2, 6.25D-2 /)

REAL (dp), PARAMETER :: zero = 0.0_dp, one = 1.0_dp, two = 2.0_dp,  &
                        three = 3.0_dp, four = 4.0_dp, five = 5.0_dp, &
                        eight = 8.0_dp, ten = 10.0_dp, c14 = 14.0_dp, &
                        c20 = 20.0_dp, c29 = 29.0_dp, c45 = 45.0_dp,  &
                        c100 = 100.0_dp

!     JACOBIAN ROUTINE SELECTOR.

SELECT CASE ( nprob )
  CASE (    1)     !     LINEAR FUNCTION - FULL RANK.

    temp = two / DBLE(m)
    fjac(1:m,1:n) = -temp
    DO  j=1,n
      fjac(j,j) = fjac(j,j) + one
    END DO

  CASE (    2)     !     LINEAR FUNCTION - RANK 1.

    DO  j=1,n
      DO  i=1,m
        fjac(i,j) = i * j
      END DO
    END DO

  CASE (    3)     !     LINEAR FUNCTION - RANK 1 WITH ZERO COLUMNS AND ROWS.

    fjac(1:m,1:n) = zero
    nm1 = n-1
    mm1 = m-1
    DO  j=2,nm1
      DO  i=2,mm1
        fjac(i,j) = (i-1) * j
      END DO
    END DO

  CASE (    4)     !     ROSENBROCK FUNCTION.

    fjac(1,1) = -c20*x(1)
    fjac(1,2) = ten
    fjac(2,1) = -one
    fjac(2,2) = zero

  CASE (    5)     !     HELICAL VALLEY FUNCTION.

    tpi = eight*ATAN(one)
    temp = x(1)**2 + x(2)**2
    tmp1 = tpi*temp
    tmp2 = SQRT(temp)
    fjac(1,1) = c100*x(2)/tmp1
    fjac(1,2) = -c100*x(1)/tmp1
    fjac(1,3) = ten
    fjac(2,1) = ten*x(1)/tmp2
    fjac(2,2) = ten*x(2)/tmp2
    fjac(2,3) = zero
    fjac(3,1) = zero
    fjac(3,2) = zero
    fjac(3,3) = one