program main !*****************************************************************************80 ! !! DQED_PRB4 tests DQED. ! ! Discussion: ! ! This program illustrates the use of the Hanson-Krogh nonlinear least ! squares solver DQED. ! ! The variables X(1:6) represent the dietary proportion of 6 prey ! animals in the diet of an alligator. ! ! The percentages of isotopic Carbon 13 and Nitrogen 15 are known ! for each of these prey animals, as well as for the alligator. ! ! The problem is to find values for the variables x(1:6) ! so that the carbon and nitrogen levels in the alligator can be ! explained purely in terms of diet. In other words, we know ! C1 through C6 and CA, N1 through N6 and NA, and we seek X1 through ! X6 so that, as closely as possible, it is true that: ! ! (X1) ! (X2) ! ( C1 C2 C3 C4 C5 C6 ) * (X3) = (CA NA) ! ( N1 N2 N3 N4 N5 N6 ) (X4) ! (X5) ! (X6) ! ! We have constraints on our variables, namely, that each X is nonnegative, ! and that the X's sum to 1. ! ! The impetus for this problem, and the field data, came from ! ! Matt Aresco, ! Biology Department, ! Florida State University, ! aresco@bio.fsu.edu ! ! Modified: ! ! 21 October 2002 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Laurel Saito, Brett Johnson, John Bartholow, Blair Hanna, ! Assessing Ecosystem Effects of Reservoir Operations Using Food ! Web-Energy Transfer and Water Quality Models, ! Ecosystems, ! Volume 4, Number 2, pages 105-125, 2001. ! implicit none integer ( kind = 4 ), parameter :: mcon = 1 integer ( kind = 4 ), parameter :: mequa = 2 integer ( kind = 4 ), parameter :: npmax = 5 integer ( kind = 4 ), parameter :: nvars = 6 integer ( kind = 4 ), parameter :: nall = mcon + 2 * nvars + npmax + 1 integer ( kind = 4 ), parameter :: nplus = 3 * nall + 2 integer ( kind = 4 ), parameter :: liwork = 3 * mcon + 9 * nvars + 4 * npmax + nall + 11 integer ( kind = 4 ), parameter :: ldfj = mcon + mequa integer ( kind = 4 ), parameter :: lwork = nall * nall + 4 * nall + nvars * npmax & + 33 * nvars + mequa * npmax + nvars * npmax + 13 * npmax & + 9 * mcon + 26 + nplus real ( kind = 8 ) bl(nvars+mcon) real ( kind = 8 ) bu(nvars+mcon) real ( kind = 8 ), dimension ( nvars ) :: c = (/ & -27.60D+00, -25.40D+00, -25.63D+00, -27.17D+00, -24.81D+00, -26.22D+00 /) real ( kind = 8 ) :: c_total = -25.34 real ( kind = 8 ) fj(ldfj,nvars+1) real ( kind = 8 ) fnorm integer ( kind = 4 ) igo integer ( kind = 4 ) ind(nvars+mcon) integer ( kind = 4 ) iopt(24) integer ( kind = 4 ) iwork(liwork) real ( kind = 8 ), dimension ( nvars ) :: n = (/ & 2.79D+00, 6.93D+00, 8.21D+00, 8.34D+00, 1.97D+00, 4.36D+00 /) real ( kind = 8 ) :: n_total = 5.10D+00 integer ( kind = 4 ) niters real ( kind = 8 ) ropt(1) external dqedev real ( kind = 8 ) work(lwork) real ( kind = 8 ) x(nvars) call timestamp ( ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DQED_PRB4' write ( *, '(a)' ) ' FORTRAN90 version' write ( *, '(a)' ) ' A test for DQED' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' This example uses reverse communication.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The stable isotope problem for the alligator.' write ( *, '(a)' ) ' We assume an alligator has 6 kinds of prey, and ' write ( *, '(a)' ) ' that the concentrations of isotopic C13 and N15 ' write ( *, '(a)' ) ' are known in the prey and in the alligator.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' We want to determine dietary proportions X1 through X6' write ( *, '(a)' ) ' so that the isotope values in the alligator can ' write ( *, '(a)' ) ' be explained' write ( *, '(a)' ) ' in terms of how much of what prey it eats.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Our model:' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' (C_All) = (C1 C2 C3 C4 C5 C6 ) * (X1 X2 X3 X4 X5 X6)''' write ( *, '(a)' ) ' (N_All) = (N1 N2 N3 N4 N5 N6 )' ! ! Define the bounding constraints on the variables. ! ! 1: 0.0 <= X(1) ! ind(1) = 1 bl(1) = 0.0D+00 bu(1) = 0.0D+00 ! ! 2: 0.0 <= X(2) ! ind(2) = 1 bl(2) = 0.0D+00 bu(2) = 0.0D+00 ! ! 3: 0.0 <= X(3) ! ind(3) = 1 bl(3) = 0.0D+00 bu(3) = 0.0D+00 ! ! 4: 0.0 <= X(4) ! ind(4) = 1 bl(4) = 0.0D+00 bu(4) = 0.0D+00 ! ! 5: 0.0 <= X(5) ! ind(5) = 1 bl(5) = 0.0D+00 bu(5) = 0.0D+00 ! ! 6: 0.0 <= X(6) ! ind(6) = 1 bl(6) = 0.0D+00 bu(6) = 0.0D+00 ! ! 7: the linear constraint 1 = X(1) + X(2) + X(3) + X(4) + X(5) + X(6) ! ind(7) = 3 bl(7) = 1.0D+00 bu(7) = 1.0D+00 ! ! Define the initial values of the variables. ! We don't know anything at all, so all variables are set equal. ! x(1:nvars) = 1.0D+00 / real ( nvars, kind = 8 ) ! ! Tell how much storage we gave the solver. ! iwork(1) = lwork iwork(2) = liwork ! ! Initialize the call counter. ! niters = 0 ! ! Use reverse commumication to evaluate the derivatives. ! iopt(1) = 16 iopt(2) = 1 ! ! No more options. ! iopt(3) = 99 do call dqed ( dqedev, mequa, nvars, mcon, ind, bl, bu, x, fj, ldfj, fnorm, & igo, iopt, ropt, iwork, work ) if ( 1 < igo ) then exit end if ! ! Count function evaluations. ! niters = niters + 1 ! ! Set the value of the constraint. ! fj(1,nvars+1) = sum ( x(1:nvars) ) ! ! Set the value of the residual functions. ! fj(mcon+1,nvars+1) = dot_product ( x(1:nvars), c(1:nvars) ) - c_total fj(mcon+2,nvars+1) = dot_product ( x(1:nvars), n(1:nvars) ) - n_total ! ! If IGO is nonzero, compute the derivatives. ! if ( igo == 0 ) then cycle end if ! ! Partial derivatives of the constraint equation. ! fj(1,1:nvars) = 1.0D+00 ! ! Partial derivatives of the least squares equations. ! fj(mcon+1,1:nvars) = c(1:nvars) fj(mcon+2,1:nvars) = n(1:nvars) end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Our computed answer X:' write ( *, '(a)' ) ' ' write ( *, '(g14.6)' ) x(1:nvars) write ( *, '(a)' ) ' ' write ( *, '(a,g14.6)' ) ' The residual of the fitted data, FNORM = ', fnorm write ( *, '(a)' ) ' ' write ( *, '(a,g14.6)' ) ' Sum of variables = ', sum ( x(1:nvars) ) write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) ' Number of model evaluations NITERS = ', niters write ( *, '(a,i6)' ) ' DQED output flag IGO = ', igo ! ! Terminate. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DQED_PRB4' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) ' ' call timestamp ( ) stop end subroutine timestamp ( ) !*****************************************************************************80 ! !! timestamp() prints the current YMDHMS date as a time stamp. ! ! Example: ! ! 31 May 2001 9:45:54.872 AM ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 15 August 2021 ! ! Author: ! ! John Burkardt ! implicit none character ( len = 8 ) ampm integer d integer h integer m integer mm character ( len = 9 ), parameter, dimension(12) :: month = (/ & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' /) integer n integer s integer values(8) integer y call date_and_time ( values = values ) y = values(1) m = values(2) d = values(3) h = values(5) n = values(6) s = values(7) mm = values(8) if ( h < 12 ) then ampm = 'AM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Noon' else ampm = 'PM' end if else h = h - 12 if ( h < 12 ) then ampm = 'PM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Midnight' else ampm = 'AM' end if end if end if write ( *, '(i2.2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm ) return end