A New Nonlinear Equations Test Problem

[ivan-pi]

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]

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]

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.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 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 <#35 (comment)>, 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]

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

$ ./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.