search

fortran-lang / minpack

Modernized Minpack: for solving nonlinear equations and nonlinear least squares problems
https://fortran-lang.github.io/minpack/
Other
90 stars 18 forks source link

A New Nonlinear Equations Test Problem #35

Open ivan-pi opened 2 years ago

ivan-pi commented 2 years ago

Just found this report which tested MINPACK and NL2SOL on a problem originating from a heart model:

A New Nonlinear Equations Test Problem, J.E. Dennis, Jr., David M. Gay, and Phuong Ahn Vu, Technical Report 83-1.P, December 1985.

A Google search gives a PDF copy at the following address: https://scholarship.rice.edu/bitstream/handle/1911/101557/TR83-16.pdf

Aside from the report, Fortran source code for the problem can also be found in opt/dqed.f on Netlib.

ivan-pi commented 2 years ago

Here's a (hopefully) faithful reproduction of the callback from the report:

      SUBROUTINE FCN(N, X, FVEC, IFLAG)
C
C    This is an example of a subroutine to evaluate the "heart model"
C    function to use with HYBRD of MINPACK
C
C   IPICK = 1 ---> full problem.
C   IPICK = 2 ---> reduced problem.
C
      INTEGER N, IFLAG
      DOUBLE PRECISION X(N), FVEC(N)
C
      INTEGER IPICK
      DOUBLE PRECISION SIGMAX, SIGMAY, SIGMAA, SIGMAB, SIGMAC, SIGMAD,
     1                 SIGMAE, SIGMAF
      COMMON  /SIGMA/  SIGMAX, SIGMAY, SIGMAA, SIGMAB, SIGMAC, SIGMAD,
     1                 SIGMAE, SIGMAF, IPICK
C
      DOUBLE PRECISION AT, ATV, AV, B, BU, BUW, BW, CT, CTV, CV, D, DU,
     1            DUW, DW, T2M3V2, T2MV2, U2M3W2, U2MW2, V2M3T2, W2M3U2
C
      IF (IPICK .NE. 1) GO TO 10
C
      T2MV2  = X(5)**2 - X(7)**2
      U2MW2  = X(6)**2 - X(8)**2
      T2M3V2 = X(5)**2 - 3.0D0 * X(7)**2
      V2M3T2 = X(7)**2 - 3.0D0 * X(5)**2
      U2M3W2 = X(6)**2 - 3.0D0 * X(8)**2
      W2M3U2 = X(8)**2 - 3.0D0 * X(6)**2
      CTV = X(3) * X(5) * X(7)
      DUW = X(4) * X(6) * X(8)
      ATV = X(1) * X(5) * X(7)
      BUW = X(2) * X(6) * X(8)
      AT  = X(1) * X(5)
      BU  = X(2) * X(6)
      CV  = X(3) * X(7)
      DW  = X(4) * X(8)
      CT  = X(3) * X(5)
      AV  = X(1) * X(7)
      DU  = X(4) * X(6)
      BW  = X(2) * X(8)
C
      FVEC(1) = X(1) + X(2) - SIGMAX
      FVEC(2) = X(3) + X(4) - SIGMAY
      FVEC(3) = AT + BU - CV - DW - SIGMAA
      FVEC(4) = AV + BW + CT + DU - SIGMAB
      FVEC(5) = X(1)*T2MV2 - 2.0D0*CTV + X(2)*U2MW2 - 2.0D0*DUW - SIGMAC
      FVEC(6) = X(3)*T2MV2 + 2.0D0*ATV + X(4)*U2MW2 + 2.0D0*BUW - SIGMAD
      FVEC(7) = AT*T2M3V2 + CV*V2M3T2 + BU*U2M3W2 + DW*W2M3U2 - SIGMAE
      FVEC(8) = CT*T2M3V2 - AV*V2M3T2 + DU*U2M3W2 - BW*W2M3U2 - SIGMAF
      GO TO 999
C
 10   T2MV2  = X(3)**2 - X(5)**2
      U2MW2  = X(4)**2 - X(6)**2
      T2M3V2 = X(3)**2 - 3.0D0 * X(5)**2
      V2M3T2 = X(5)**2 - 3.0D0 * X(3)**2
      U2M3W2 = X(4)**2 - 3.0D0 * X(6)**2
      W2M3U2 = X(6)**2 - 3.0D0 * X(4)**2
      B = SIGMAX - X(1)
      D = SIGMAY - X(2)
      CTV = X(2) * X(3) * X(5)
      DUW = D * X(4) * X(6)
      ATV = X(1) * X(3) * X(5)
      BUW = B * X(4) * X(6)
      AT  = X(1) * X(3)
      BU  = B * X(4)
      CV  = X(2) * X(5)
      DW  = D * X(6)
      CT  = X(2) * X(3)
      AV  = X(1) * X(5)
      DU  = D * X(4)
      BW  = B * X(6)
C
      FVEC(1) = AT + BU - CV - DW - SIGMAA
      FVEC(2) = AV + BW + CT + DU - SIGMAB
      FVEC(3) = X(1)*T2MV2 - 2.0D0*CTV + B*U2MW2 - 2.0D0*DUW - SIGMAC
      FVEC(4) = X(2)*T2MV2 + 2.0D0*ATV + D*U2MW2 + 2.0D0*BUW - SIGMAD
      FVEC(5) = AT*T2M3V2 + CV*V2M3T2 + BU*U2M3W2 + DW*W2M3U2 - SIGMAE
      FVEC(6) = CT*T2M3V2 - AV*V2M3T2 + DU*U2M3W2 - BW*W2M3U2 - SIGMAF
 999  RETURN
C
      END

I've read the whole subroutine out loud side-by-side with the original. I've also checked it with implicit none to make sure all variables are defined. The code is perfectly legal <=F2008. The COMMON block is obsolescent since F2018.

arjenmarkus commented 2 years ago

The COMMON blocks allow you to set the parameters of the minimisation problem "from the outside". Since Fortran 90 we can use module variables instead, thereby deleting the COMMON blocks. Then an auxiliary routine that simply sets the problem-specific parameters and you have a nice test case, I'd say. (Although I do not quite understand the model ;))

Op di 15 feb. 2022 om 10:58 schreef Ivan Pribec @.***>:

Here's a (hopefully) faithful reproduction of the callback from the report:

  SUBROUTINE FCN(N, X, FVEC, IFLAG)CC    This is an example of a subroutine to evaluate the "heart model"C    function to use with HYBRD of MINPACKCC   IPICK = 1 ---> full problem.C   IPICK = 2 ---> reduced problem.C
  INTEGER N, IFLAG
  DOUBLE PRECISION X(N), FVEC(N)C
  INTEGER IPICK
  DOUBLE PRECISION SIGMAX, SIGMAY, SIGMAA, SIGMAB, SIGMAC, SIGMAD,
 1                 SIGMAE, SIGMAF
  COMMON  /SIGMA/  SIGMAX, SIGMAY, SIGMAA, SIGMAB, SIGMAC, SIGMAD,
 1                 SIGMAE, SIGMAF, IPICKC
  DOUBLE PRECISION AT, ATV, AV, B, BU, BUW, BW, CT, CTV, CV, D, DU,
 1            DUW, DW, T2M3V2, T2MV2, U2M3W2, U2MW2, V2M3T2, W2M3U2C
  IF (IPICK .NE. 1) GO TO 10C
  T2MV2  = X(5)**2 - X(7)**2
  U2MW2  = X(6)**2 - X(8)**2
  T2M3V2 = X(5)**2 - 3.0D0 * X(7)**2
  V2M3T2 = X(7)**2 - 3.0D0 * X(5)**2
  U2M3W2 = X(6)**2 - 3.0D0 * X(8)**2
  W2M3U2 = X(8)**2 - 3.0D0 * X(6)**2
  CTV = X(3) * X(5) * X(7)
  DUW = X(4) * X(6) * X(8)
  ATV = X(1) * X(5) * X(7)
  BUW = X(2) * X(6) * X(8)
  AT  = X(1) * X(5)
  BU  = X(2) * X(6)
  CV  = X(3) * X(7)
  DW  = X(4) * X(8)
  CT  = X(3) * X(5)
  AV  = X(1) * X(7)
  DU  = X(4) * X(6)
  BW  = X(2) * X(8)C
  FVEC(1) = X(1) + X(2) - SIGMAX
  FVEC(2) = X(3) + X(4) - SIGMAY
  FVEC(3) = AT + BU - CV - DW - SIGMAA
  FVEC(4) = AV + BW + CT + DU - SIGMAB
  FVEC(5) = X(1)*T2MV2 - 2.0D0*CTV + X(2)*U2MW2 - 2.0D0*DUW - SIGMAC
  FVEC(6) = X(3)*T2MV2 + 2.0D0*ATV + X(4)*U2MW2 + 2.0D0*BUW - SIGMAD
  FVEC(7) = AT*T2M3V2 + CV*V2M3T2 + BU*U2M3W2 + DW*W2M3U2 - SIGMAE
  FVEC(8) = CT*T2M3V2 - AV*V2M3T2 + DU*U2M3W2 - BW*W2M3U2 - SIGMAF
  GO TO 999C

10 T2MV2 = X(3)2 - X(5)2 U2MW2 = X(4)2 - X(6)2 T2M3V2 = X(3)2 - 3.0D0 X(5)2 V2M3T2 = X(5)2 - 3.0D0 X(3)2 U2M3W2 = X(4)2 - 3.0D0 X(6)2 W2M3U2 = X(6)2 - 3.0D0 X(4)2 B = SIGMAX - X(1) D = SIGMAY - X(2) CTV = X(2) X(3) X(5) DUW = D X(4) X(6) ATV = X(1) X(3) X(5) BUW = B X(4) X(6) AT = X(1) X(3) BU = B X(4) CV = X(2) X(5) DW = D X(6) CT = X(2) X(3) AV = X(1) X(5) DU = D X(4) BW = B X(6)C FVEC(1) = AT + BU - CV - DW - SIGMAA FVEC(2) = AV + BW + CT + DU - SIGMAB FVEC(3) = X(1)T2MV2 - 2.0D0CTV + BU2MW2 - 2.0D0DUW - SIGMAC FVEC(4) = X(2)T2MV2 + 2.0D0ATV + DU2MW2 + 2.0D0BUW - SIGMAD FVEC(5) = ATT2M3V2 + CVV2M3T2 + BUU2M3W2 + DWW2M3U2 - SIGMAE FVEC(6) = CTT2M3V2 - AVV2M3T2 + DUU2M3W2 - BWW2M3U2 - SIGMAF 999 RETURNC END

I've read the whole subroutine out loud side-by-side with the original. I've also checked it with implicit none to make sure all variables are defined. The code is perfectly legal <=F2008. The COMMON block is obsolescent since F2018.

— Reply to this email directly, view it on GitHub https://github.com/fortran-lang/minpack/issues/35#issuecomment-1040075224, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAN6YR3UTJQKME7RAQ6FVFLU3IPTLANCNFSM5N7PKH2A . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

You are receiving this because you are subscribed to this thread.Message ID: @.***>

ivan-pi commented 2 years ago

37 years later, the results match exactly with the values reported in the original document :rocket:.

$ ./heart_test
 Experiment file: exp791129.txt
sigmax =      4.850000000E-01
sigmay =     -1.900000000E-03
sigmaa =     -5.810000000E-02
sigmab =      1.500000000E-02
sigmac =      1.050000000E-01
sigmad =      4.060000000E-02
sigmae =      1.670000000E-01
sigmaf =     -3.990000000E-01
        x                 xt        
    -6.321349025E-03    -6.321349025E-03
     4.913213490E-01     4.913213490E-01
    -1.998156408E-03    -1.998156408E-03
     9.815640840E-05     9.815640840E-05
     1.226569755E-01     1.226569755E-01
    -1.003153205E-01    -1.003153205E-01
    -4.023517593E+00    -4.023517593E+00
    -2.071785527E-02    -2.071785527E-02
norm  ( 0.5 F(x)^T . F(X) ) =      2.883244015E-17

I'll see to submit a test driver soon.