C ================================================================= C File: pedazos4.f C ================================================================= C ================================================================= C Module: Subroutines that define the problem C ================================================================= C Last update of any of the component of this module: C March 4, 2008. C Users are encouraged to download periodically updated versions of C this code at the COLLECTION home page: C C www.ime.usp.br/~egbirgin/collection/ C C and periodically updated versions of the TANGO Project solvers at C the TANGO home page: C C www.ime.usp.br/~egbirgin/tango/ C ================================================================= C Pedazos4 problem: Fit a continuous piecewise polynomial C ------------------------------------------------------- C We generate a table (t_i, y_i) with 30 equally spaced points C between 1 and 30. The first 10 points are on a line, the second C 10 on a parabola and the third 10 on a cubic. The overall C function is continuous at 10 and 20. We fit a line to the first C 10, a quadratic to the second 10 and a cubic to the third ten, C with the conditions that we preserve continuity at 10 and 20. So, C we have 9 unknowns and 2 equality constraints. The objective C function is quadratic and the constraints are linear. C ****************************************************************** C ****************************************************************** subroutine inip(n,x,l,u,m,lambda,equatn,linear,coded,checkder) implicit none C SCALAR ARGUMENTS integer m,n logical checkder C ARRAY ARGUMENTS logical coded(10),equatn(*),linear(*) double precision l(*),lambda(*),u(*),x(*) C COMMON ARRAYS double precision t(30),y(30) C LOCAL SCALARS integer i double precision drand,seed C COMMON BLOCKS common /probdata/ t,y C Number of variables n = 9 C Adjust a line, a parabola and a cubic do i = 1,10 t(i) = dfloat(i) y(i) = t(i) end do do i = 11,20 t(i) = dfloat(i) y(i) = 0.4d0 * ( t(i) - 15.0d0 ) ** 2 end do do i = 21,30 t(i) = dfloat(i) y(i) = - 10.0d0 * ( t(i) - 25.0d0 ) ** 3 / 125.0d0 end do C Perturbation seed = 17172937.0d0 do i = 1,30 y(i) = y(i) + + y(i) * 0.20d0 * ( drand(seed) * 2.0d0 - 1.0d0 ) end do C Initial point do i = 1,n x(i) = 0.0d0 end do C Lower and upper bounds do i = 1,n l(i) = - 1.0d+20 u(i) = 1.0d+20 end do C Number of constraints (equalities plus inequalities) m = 2 C Lagrange multipliers approximation. do i = 1,m lambda(i) = 0.0d0 end do C For each constraint i, set equatn(i) = .true. if it is an equality C constraint of the form c_i(x) = 0, and set equatn(i) = .false. if C it is an inequality constraint of the form c_i(x) <= 0. do i = 1,m equatn(i) = .true. end do C For each constraint i, set linear(i) = .true. if it is a linear C constraint, otherwise set linear(i) = .false. do i = 1,m linear(i) = .true. end do C Indicate which subroutines did you code. coded( 1) = .true. ! evalf coded( 2) = .true. ! evalg coded( 3) = .true. ! evalh coded( 4) = .true. ! evalc coded( 5) = .true. ! evaljac coded( 6) = .true. ! evalhc coded( 7) = .false. ! evalfc coded( 8) = .false. ! evalgjac coded( 9) = .false. ! evalhl coded(10) = .false. ! evalhlp C Set checkder = .true. if you code some derivatives and you would C like them to be tested by finite differences. It is highly C recommended. checkder = .false. end C ****************************************************************** C ****************************************************************** subroutine evalf(n,x,f,flag) implicit none C SCALAR ARGUMENTS integer flag,n double precision f C ARRAY ARGUMENTS double precision x(n) C COMMON ARRAYS double precision t(30),y(30) C LOCAL SCALARS integer i double precision r C COMMON BLOCKS common /probdata/ t,y flag = 0 f = 0.0d0 do i = 1,30 if ( t(i) .le. 10.0d0 ) then r = ( x(1) + x(2) * t(i) - y(i) ) ** 2 else if ( t(i) .le. 20.0d0 ) then r = ( x(3) + x(4) * t(i) + x(5) * t(i) ** 2 - y(i) ) ** 2 else r = ( x(6) + x(7) * t(i) + x(8) * t(i) ** 2 + + x(9) * t(i) ** 3 - y(i) ) ** 2 end if f = f + r end do end C ****************************************************************** C ****************************************************************** subroutine evalg(n,x,g,flag) implicit none C SCALAR ARGUMENTS integer flag,n C ARRAY ARGUMENTS double precision g(n),x(n) C COMMON ARRAYS double precision t(30),y(30) C LOCAL SCALARS integer i double precision r C COMMON BLOCKS common /probdata/ t,y flag = 0 do i = 1,n g(i) = 0.0d0 end do do i = 1,30 if ( t(i) .le. 10.0d0 ) then r = x(1) + x(2) * t(i) - y(i) g(1) = g(1) + 2.0d0 * r g(2) = g(2) + 2.0d0 * r * t(i) else if ( t(i) .le. 20.0d0 ) then r = x(3) + x(4) * t(i) + x(5) * t(i) ** 2 - y(i) g(3) = g(3) + 2.0d0 * r g(4) = g(4) + 2.0d0 * r * t(i) g(5) = g(5) + 2.0d0 * r * t(i) ** 2 else r = x(6) + x(7) * t(i) + x(8) * t(i) ** 2 + + x(9) * t(i) ** 3 - y(i) g(6) = g(6) + 2.0d0 * r g(7) = g(7) + 2.0d0 * r * t(i) g(8) = g(8) + 2.0d0 * r * t(i) ** 2 g(9) = g(9) + 2.0d0 * r * t(i) ** 3 end if end do end C ****************************************************************** C ****************************************************************** subroutine evalh(n,x,hlin,hcol,hval,hnnz,flag) implicit none C SCALAR ARGUMENTS integer flag,hnnz,n C ARRAY ARGUMENTS integer hcol(*),hlin(*) double precision hval(*),x(n) flag = 0 hnnz = 19 hlin(1) = 1 hcol(1) = 1 hval(1) = 20.0d0 hlin(2) = 2 hcol(2) = 1 hval(2) = 110.0d0 hlin(3) = 2 hcol(3) = 2 hval(3) = 770.0d0 hlin(4) = 3 hcol(4) = 3 hval(4) = 20.0d0 hlin(5) = 4 hcol(5) = 3 hval(5) = 310.0d0 hlin(6) = 4 hcol(6) = 4 hval(6) = 4970.0d0 hlin(7) = 5 hcol(7) = 3 hval(7) = 4970.0d0 hlin(8) = 5 hcol(8) = 4 hval(8) = 82150.0d0 hlin(9) = 5 hcol(9) = 5 hval(9) = 1394666.0d0 hlin(10) = 6 hcol(10) = 6 hval(10) = 20.0d0 hlin(11) = 7 hcol(11) = 6 hval(11) = 510.0d0 hlin(12) = 7 hcol(12) = 7 hval(12) = 13170.0d0 hlin(13) = 8 hcol(13) = 6 hval(13) = 13170.0d0 hlin(14) = 8 hcol(14) = 7 hval(14) = 344250.0d0 hlin(15) = 8 hcol(15) = 8 hval(15) = 9102666.0d0 hlin(16) = 9 hcol(16) = 6 hval(16) = 344250.0d0 hlin(17) = 9 hcol(17) = 7 hval(17) = 9102666.0d0 hlin(18) = 9 hcol(18) = 8 hval(18) = 243308250.0d0 hlin(19) = 9 hcol(19) = 9 hval(19) = 6.56895081d+09 end C ****************************************************************** C ****************************************************************** subroutine evalc(n,x,ind,c,flag) implicit none C SCALAR ARGUMENTS integer ind,flag,n double precision c C ARRAY ARGUMENTS double precision x(n) flag = 0 if ( ind .eq. 1 ) then c = ( x(1) + x(2) * 10.d0 ) - + ( x(3) + x(4) * 10.d0 + x(5) * 100.d0 ) else if( ind .eq. 2 ) then c = ( x(3) + x(4) * 20.d0 + x(5) * 400.d0 ) - + ( x(6) + x(7) * 20.d0 + x(8) * 400.d0 + x(9) * 8000.d0 ) end if end C ****************************************************************** C ****************************************************************** subroutine evaljac(n,x,ind,jcvar,jcval,jcnnz,flag) implicit none C SCALAR ARGUMENTS integer flag,ind,jcnnz,n C ARRAY ARGUMENTS integer jcvar(n) double precision x(n),jcval(n) flag = 0 if ( ind .eq. 1 ) then jcnnz = 5 jcvar(1) = 1 jcval(1) = 1.d0 jcvar(2) = 2 jcval(2) = 10.d0 jcvar(3) = 3 jcval(3) = - 1.d0 jcvar(4) = 4 jcval(4) = - 10.d0 jcvar(5) = 5 jcval(5) = - 100.d0 else if( ind .eq. 2 ) then jcnnz = 7 jcvar(1) = 3 jcval(1) = 1.d0 jcvar(2) = 4 jcval(2) = 20.d0 jcvar(3) = 5 jcval(3) = 400.d0 jcvar(4) = 6 jcval(4) = - 1.d0 jcvar(5) = 7 jcval(5) = - 20.d0 jcvar(6) = 8 jcval(6) = - 400.d0 jcvar(7) = 9 jcval(7) = - 8000.d0 end if end C ****************************************************************** C ****************************************************************** subroutine evalhc(n,x,ind,hclin,hccol,hcval,hcnnz,flag) implicit none C SCALAR ARGUMENTS integer flag,hcnnz,ind,n C ARRAY ARGUMENTS integer hccol(*),hclin(*) double precision hcval(*),x(n) flag = 0 hcnnz = 0 end C ****************************************************************** C ****************************************************************** subroutine evalhl(n,x,m,lambda,scalef,scalec,hllin,hlcol,hlval, +hlnnz,flag) implicit none C SCALAR ARGUMENTS integer flag,hlnnz,m,n double precision scalef C ARRAY ARGUMENTS integer hlcol(*),hllin(*) double precision hlval(*),lambda(m),scalec(m),x(n) flag = - 1 end C ****************************************************************** C ****************************************************************** subroutine evalhlp(n,x,m,lambda,scalef,scalec,p,hp,goth,flag) implicit none C SCALAR ARGUMENTS logical goth integer flag,m,n double precision scalef C ARRAY ARGUMENTS double precision hp(n),lambda(m),p(n),scalec(m),x(n) flag = - 1 end C ****************************************************************** C ****************************************************************** subroutine endp(n,x,l,u,m,lambda,equatn,linear) implicit none C SCALAR ARGUMENTS integer m,n C ARRAY ARGUMENTS logical equatn(m),linear(m) double precision l(n),lambda(m),u(n),x(n) C COMMON ARRAYS double precision t(30),y(30) C LOCAL SCALARS integer i double precision z,tt C COMMON BLOCKS common /probdata/ t,y write(*,*) 't, observed y, computed y' do i = 1,30 if ( t(i) .le. 10.0d0 ) then z = x(1) + x(2) * t(i) else if ( t(i) .le. 20.0d0 ) then z = x(3) + x(4) * t(i) + x(5) * t(i) ** 2 else z = x(6) + x(7) * t(i) + x(8) * t(i) ** 2 + + x(9) * t(i) ** 3 end if write(*,*) t(i),y(i),z end do tt = 1.0d0 10 if ( tt .le. 30.01d0 ) then if ( tt .le. 10.0d0 ) then z = x(1) + x(2) * tt else if ( tt .le. 20.0d0 ) then z = x(3) + x(4) * tt + x(5) * tt ** 2 else z = x(6) + x(7) * tt + x(8) * tt ** 2 + + x(9) * tt ** 3 end if write(*,*) tt,z tt = tt + 0.01d0 go to 10 end if end