program main !*****************************************************************************80 ! !! MAIN is the main program for DQED_PRB1. ! ! Discussion: ! ! DQED_PRB1 tests the DQED library. ! ! This routine illustrates the use of DQED, the Hanson-Krogh nonlinear least ! squares solver, by solving a certain heart dipole moment equation. ! ! Modified: ! ! 11 September 2002 ! ! Reference: ! ! John Dennis, David Gay, Phuong Vu. ! A new nonlinear equations test problem, ! Rice University Department of Mathematics ! Report 83-16, 6/83, 6/85 ! implicit none integer ( kind = 4 ), parameter :: ldfj = 8 integer ( kind = 4 ), parameter :: liwa = 150 integer ( kind = 4 ), parameter :: lwa = 2000 integer ( kind = 4 ), parameter :: nvars = 8 real ( kind = 8 ) bl(nvars+2) real ( kind = 8 ) bu(nvars+2) real ( kind = 8 ) df(nvars) external dqedhd real ( kind = 8 ) fj(ldfj,nvars+1) integer ( kind = 4 ) i integer ( kind = 4 ) ind(nvars+2) integer ( kind = 4 ) iopt(28) integer ( kind = 4 ) ios integer ( kind = 4 ) iwa(liwa) integer ( kind = 4 ) mode real ( kind = 8 ) ropt(1) real ( kind = 8 ) sigma(nvars) character ( len = 20 ) title real ( kind = 8 ) wa(lwa) ! ! To have the jacobian checked, set WANT_PR to true. ! logical, parameter :: want_pr = .false. real ( kind = 8 ) x(nvars) real ( kind = 8 ) xsave(nvars) real ( kind = 8 ) y(nvars) common /sigma/ sigma common /mode/ mode call timestamp ( ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DQED_PRB1' write ( *, '(a)' ) ' FORTRAN90 version' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' A set of tests for DQED, which can solve' write ( *, '(a)' ) ' bounded and constrained linear least squares problems' write ( *, '(a)' ) ' and systems of nonlinear equations.' do ! ! Read the problem title. ! read ( *, '(a)', iostat = ios ) title if ( ios /= 0 ) then exit end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) ! ! Read in the sums of values and initial parameter values. ! read ( *, *, iostat = ios ) sigma(1:nvars), xsave(1:nvars) if ( ios /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Unexpected end of file or other error.' exit end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'The input SIGMA vector:' write ( *, '(a)' ) ' ' do i = 1, nvars write ( *, '(i6,g14.6)' ) i, sigma(i) end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'The initial estimate for X:' write ( *, '(a)' ) ' ' do i = 1, nvars write ( *, '(i6,g14.6)' ) i, xsave(i) end do ! ! Test the partial derivative computation. ! if ( want_pr ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Test the partial derivative computation:' write ( *, '(a)' ) ' ' call dpchek ( df, dqedhd, fj, iopt, ldfj, nvars, ropt, xsave, y ) end if ! ! Test 1, with all equations normal. ! mode = 0 x(1:nvars) = xsave(1:nvars) call test01 ( bl, bu, fj, ind, iopt, iwa, ldfj, liwa, lwa, mode, nvars, & ropt, sigma, wa, x ) ! ! Test 1, with first two linear equations used as constraints. ! mode = 1 x(1:nvars) = xsave(1:nvars) call test01 ( bl, bu, fj, ind, iopt, iwa, ldfj, liwa, lwa, mode, nvars, & ropt, sigma, wa, x ) ! ! Test 2, with all equations normal. ! mode = 0 x(1:nvars) = xsave(1:nvars) call test02 ( bl, bu, fj, ind, iopt, iwa, ldfj, liwa, lwa, mode, nvars, & ropt, sigma, wa, x ) ! ! Test 2, with first two linear equations used as constraints. ! mode = 1 x(1:nvars) = xsave(1:nvars) call test02 ( bl, bu, fj, ind, iopt, iwa, ldfj, liwa, lwa, mode, nvars, & ropt, sigma, wa, x ) end do ! ! Terminate. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DQED_PRB1' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) ' ' call timestamp ( ) stop end subroutine test01 ( bl, bu, fj, ind, iopt, iwa, ldfj, liwa, lwa, mode, & nvars, ropt, sigma, wa, x ) !*****************************************************************************80 ! !! TEST01 uses an analytic jacobian. ! ! Discussion: ! ! The linear equations will not be constraints if MODE = 0. ! Convergence is normally faster if the linear equations are used ! as constraints, MODE = 1. ! implicit none integer ( kind = 4 ) ldfj integer ( kind = 4 ) liwa integer ( kind = 4 ) lwa integer ( kind = 4 ) nvars real ( kind = 8 ) bl(nvars+2) real ( kind = 8 ) bu(nvars+2) external dqedhd real ( kind = 8 ) fj(ldfj,nvars+1) real ( kind = 8 ) fnorm integer ( kind = 4 ) igo integer ( kind = 4 ) ind(nvars+2) integer ( kind = 4 ) iopt(28) integer ( kind = 4 ) iwa(liwa) integer ( kind = 4 ) mcon integer ( kind = 4 ) mequa integer ( kind = 4 ) mode real ( kind = 8 ) ropt(1) real ( kind = 8 ) sigma(nvars) real ( kind = 8 ) wa(lwa) real ( kind = 8 ) x(nvars) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST01' write ( *, '(a)' ) ' Use an analytic jacobian.' write ( *, '(a,i6)' ) ' MODE = ', mode ! ! Tell how much storage the solver has. ! iwa(1) = lwa iwa(2) = liwa ! ! Set up print option. (not used here). ! iopt(1) = -1 iopt(2) = 1 ! ! Set up for linear model without any quadratic terms. (not used). ! iopt(3) = -14 iopt(4) = 1 ! ! Do not allow convergence to be claimed on small steps. ! iopt(5) = 17 iopt(6) = 1 iopt(7) = 15 ! ! Allow up to NVARS quadratic model terms. ! iopt(8) = nvars ! ! Change condition number for quadratic model degree. ! iopt(9) = 10 iopt(10) = 1 ropt(1) = 1.0D+04 ! ! No more options. ! iopt(11) = 99 ! ! MODE = 0, no constraints. ! if ( mode == 0 ) then mcon = 0 ! ! MODE = 1, there are constraints. ! ! (The first two equations are linear, and can be used as constraints). ! else mcon = 2 ind(nvars+1) = 3 ind(nvars+2) = 3 bl(nvars+1) = sigma(1) bu(nvars+1) = sigma(1) bl(nvars+2) = sigma(2) bu(nvars+2) = sigma(2) end if mequa = nvars - mcon ! ! All variables are otherwise free. ! ind(1:nvars) = 4 call dqed ( dqedhd, mequa, nvars, mcon, ind, bl, bu, x, fj, ldfj, fnorm, & igo, iopt, ropt, iwa, wa ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Computed minimizing X:' write ( *, '(a)' ) ' ' write ( *, '(g14.6)' ) x(1:nvars) write ( *, '(a)' ) ' ' write ( *, '(a,g14.6)' ) ' Residual after the fit = ',fnorm write ( *, '(a,i6)' ) ' DQED output flag IGO = ',igo return end subroutine test02 ( bl, bu, fj, ind, iopt, iwa, ldfj, liwa, lwa, mode, & nvars, ropt, sigma, wa, x ) !*****************************************************************************80 ! !! TEST02 uses a finite-difference approximated jacobian. ! ! Discussion: ! ! The linear equations will not be constraints if MODE = 0. ! Convergence is normally faster if the linear equations are used ! as constraints, MODE = 1. ! implicit none integer ( kind = 4 ) ldfj integer ( kind = 4 ) liwa integer ( kind = 4 ) lwa integer ( kind = 4 ) nvars real ( kind = 8 ) bl(nvars+2) real ( kind = 8 ) bu(nvars+2) real ( kind = 8 ) fj(ldfj,nvars+1) external fjaprx real ( kind = 8 ) fnorm integer ( kind = 4 ) igo integer ( kind = 4 ) ind(nvars+2) integer ( kind = 4 ) iopt(28) integer ( kind = 4 ) iwa(liwa) integer ( kind = 4 ) mcon integer ( kind = 4 ) mequa integer ( kind = 4 ) mode real ( kind = 8 ) ropt(1) real ( kind = 8 ) sigma(nvars) real ( kind = 8 ) wa(lwa) real ( kind = 8 ) x(nvars) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST02' write ( *, '(a)' ) ' Use an approximate jacobian.' write ( *, '(a,i6)' ) ' MODE = ', mode ! ! Tell how much storage the solver has. ! iwa(1) = lwa iwa(2) = liwa ! ! Set up print option. (not used here). ! iopt(1) = -1 iopt(2) = 1 ! ! Set up for linear model without any quadratic terms. (not used). ! iopt(3) = -14 iopt(4) = 1 ! ! Do not allow convergence to be claimed on small steps. ! iopt(5) = 17 iopt(6) = 1 iopt(7) = 15 ! ! Allow up to NVARS quadratic model terms. ! iopt(8) = nvars ! ! Change condition number for quadratic model degree. ! iopt(9) = 10 iopt(10) = 1 ropt(1) = 10000.0D+00 ! ! No more options. ! iopt(11) = 99 ! ! MODE = 0, no constraints. ! if ( mode == 0 ) then mcon = 0 ! ! MODE = 1, there are constraints. ! (The first two equations are linear, and can be used as ! constraints instead). ! else mcon = 2 ind(nvars+1) = 3 ind(nvars+2) = 3 bl(nvars+1) = sigma(1) bu(nvars+1) = sigma(1) bl(nvars+2) = sigma(2) bu(nvars+2) = sigma(2) end if mequa = nvars - mcon ! ! All variables are otherwise free. ! ind(1:nvars) = 4 call dqed ( fjaprx, mequa, nvars, mcon, ind, bl, bu, x, fj, ldfj, fnorm, & igo, iopt, ropt, iwa, wa ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Computed minimizing X:' write ( *, '(a)' ) ' ' write ( *, '(g14.6)' ) x(1:nvars) write ( *, '(a)' ) ' ' write ( *, '(a,g14.6)' ) ' Residual after the fit = ', fnorm write ( *, '(a,i6)' ) ' DQED output flag IGO = ', igo return end subroutine dqedhd ( x, fj, ldfj, igo, iopt, ropt ) !*****************************************************************************80 ! !! DQEDHD evaluates functions and derivatives for DQED. ! ! Discussion: ! ! The user problem has MCON constraint functions, ! MEQUA least squares equations, and involves NVARS ! unknown variables. ! ! When this subprogram is entered, the general (near) ! linear constraint partial derivatives, the derivatives ! for the least squares equations, and the associated ! function values are placed into the array FJ. ! all partials and functions are evaluated at the point ! in X. then the subprogram returns to the calling ! program unit. ! ! Typically one could do the following steps: ! ! if ( igo /= 0 ) then ! place the partials of the i-th constraint function with respect to ! variable j in the array fj(i,j), i = 1,...,mcon, j=1,...,nvars. ! ! place the values of the i-th constraint equation into fj(i,nvars+1). ! ! if ( igo /= 0 ) then ! place the partials of the i-th least squares equation with respect ! to variable j in the array fj(i,j), i = 1,...,mequa, j = 1,...,nvars. ! ! place the value of the i-th least squares equation into fj(i,nvars+1). ! implicit none integer ( kind = 4 ) ldfj integer ( kind = 4 ), parameter :: nvars = 8 real ( kind = 8 ) fj(ldfj,nvars+1) integer ( kind = 4 ) igo integer ( kind = 4 ) iopt(*) integer ( kind = 4 ) mcon integer ( kind = 4 ) mequa integer ( kind = 4 ) mode real ( kind = 8 ) ropt(*) real ( kind = 8 ) x(nvars) common /mode/ mode if ( mode == 0 ) then mcon = 0 else mcon = 2 end if mequa = nvars - mcon if ( igo /= 0 ) then call jack ( fj, ldfj, nvars, x ) end if call func ( fj(1,9), iopt, mcon, mequa, nvars, ropt, x ) return end subroutine fjaprx ( x, fj, ldfj, igo, iopt, ropt ) !*****************************************************************************80 ! !! FJAPRX uses DIFCEN to approximate the jacobian matrix. ! implicit none integer ( kind = 4 ) ldfj integer ( kind = 4 ), parameter :: nvars = 8 real ( kind = 8 ) fj(ldfj,nvars+1) external func real ( kind = 8 ) fx(8) integer ( kind = 4 ) igo integer ( kind = 4 ) iopt(*) integer ( kind = 4 ) mcon integer ( kind = 4 ) mequa integer ( kind = 4 ) mode real ( kind = 8 ) ropt(*) real ( kind = 8 ) x(nvars) common /mode/ mode if ( mode == 0 ) then mcon = 0 else mcon = 2 end if mequa = nvars - mcon call difcen ( fj, func, fx, iopt, ldfj, mcon, mequa, nvars, ropt, x ) call func ( fj(1,nvars+1), iopt, mcon, mequa, nvars, ropt, x ) return end subroutine func ( fx, iopt, mcon, mequa, nvars, ropt, x ) !*****************************************************************************80 ! !! FUNC evaluates the functions. ! ! Discussion: ! ! The system of equations has the form: ! ! x(1)+x(2) = sigma(1) ! ! x(3)+x(4) = sigma(2) ! ! x(1)*x(5)+x(2)*x(6)-x(3)*x(7)-x(4)*x(8) = sigma(3) ! ! x(1)*x(7)+x(2)*x(8)+x(3)*x(5)+x(4)*x(6) = sigma(4) ! ! x(1)*(x(5)**2-x(7)**2)+x(2)*(x(6)**2-x(8)**2)+x(3)*(-2.0*x(5)*x(7)) ! +x(4)*(-2.0*x(6)*x(8)) = sigma(5) ! ! x(1)*(2.0*x(5)*x(7))+x(2)*(2.0*x(6)*x(8))+x(3)*(x(5)**2-x(7)**2) ! +x(4)*(x(6)**2-x(8)**2) = sigma(6) ! ! x(1)*(x(5)*(x(5)**2-3.0*x(7)**2))+x(2)*(x(6)*(x(6)**2-3.0*x(8)**2)) ! +x(3)*(x(7)*(x(7)**2-3.0*x(5)**2))+x(4)*(x(8)*(x(8)**2-3.0*x(6)**2)) ! = sigma(7) ! ! x(1)*(-x(7)*(x(7)**2-3.0*x(5)**2))+x(2)*(-x(8)*(x(8)**2-3.0*x(6)**2)) ! +x(3)*(x(5)*(x(5)**2-3.0*x(7)**2))+x(4)*(x(6)*(x(6)**2-3.0*x(8)**2)) ! = sigma(8) ! implicit none integer ( kind = 4 ) mcon integer ( kind = 4 ) mequa integer ( kind = 4 ) nvars real ( kind = 8 ) fx(mequa+mcon) integer ( kind = 4 ) iopt(*) integer ( kind = 4 ) mode real ( kind = 8 ) ropt(*) real ( kind = 8 ) sigma(8) real ( kind = 8 ) x(nvars) common /sigma/ sigma common /mode/ mode ! if ( mode == 0 ) then fx(1) = x(1) + x(2) - sigma(1) else fx(1) = x(1) + x(2) end if if ( mode == 0 ) then fx(2) = x(3) + x(4) - sigma(2) else fx(2) = x(3) + x(4) end if fx(3) = x(1)*x(5) + x(2)*x(6) - x(3)*x(7) - x(4)*x(8) - sigma(3) fx(4) = x(1)*x(7) + x(2)*x(8) + x(3)*x(5) + x(4)*x(6) - sigma(4) fx(5) = x(1)*(x(5)**2-x(7)**2) + x(2)*(x(6)**2-x(8)**2) & +x(3)*(-2.0D+00*x(5)*x(7)) + x(4)*(-2.0D+00*x(6)*x(8)) - sigma(5) fx(6) = x(1)*(2.0D+00*x(5)*x(7)) + x(2)*(2.0D+00*x(6)*x(8)) & + x(3)*(x(5)**2-x(7)**2) & + x(4)*(x(6)**2-x(8)**2) - sigma(6) fx(7) = x(1)*(x(5)*(x(5)**2-3.0*x(7)**2)) & + x(2)*(x(6)*(x(6)**2-3.0D+00*x(8)**2)) & + x(3)*(x(7)*(x(7)**2-3.0D+00*x(5)**2)) & + x(4)*(x(8)*(x(8)**2-3.0D+00*x(6)**2)) & - sigma(7) fx(8) = x(1)*(-x(7)*(x(7)**2-3.0D+00*x(5)**2)) & + x(2)*(-x(8)*(x(8)**2-3.0D+00*x(6)**2)) & + x(3)*(x(5)*(x(5)**2-3.0D+00*x(7)**2)) & + x(4)*(x(6)*(x(6)**2-3.0D+00*x(8)**2)) - sigma(8) return end subroutine jack ( fj, ldfj, nvars, x ) !*****************************************************************************80 ! !! JACK evaluates the partial derivatives of the functions. ! ! Parameters: ! ! Output, real ( kind = 8 ) FJ(LDFJ,NVARS), the MCON+MEQUA by NVARS ! matrix of partial derivatives. ! ! Input, integer ( kind = 4 ) LDFJ, the leading dimension of FJ, which must be ! at least MCON+MEQUA. ! ! Input, integer ( kind = 4 ) NVARS, the number of variables. ! ! Input, real X(NVARS), the point at which the partial derivatives ! are to be evaluated. ! implicit none integer ( kind = 4 ) ldfj integer ( kind = 4 ) nvars real ( kind = 8 ) fj(ldfj,nvars) real ( kind = 8 ) x(nvars) ! ! Equation #1: ! ! x(1)+x(2) = sigma(1) ! fj(1,1) = 1.0D+00 fj(1,2) = 1.0D+00 fj(1,3) = 0.0D+00 fj(1,4) = 0.0D+00 fj(1,5) = 0.0D+00 fj(1,6) = 0.0D+00 fj(1,7) = 0.0D+00 fj(1,8) = 0.0D+00 ! ! Equation #2: ! ! x(3)+x(4) = sigma(2) ! fj(2,1) = 0.0D+00 fj(2,2) = 0.0D+00 fj(2,3) = 1.0D+00 fj(2,4) = 1.0D+00 fj(2,5) = 0.0D+00 fj(2,6) = 0.0D+00 fj(2,7) = 0.0D+00 fj(2,8) = 0.0D+00 ! ! Equation #3: ! x(1)*x(5)+x(2)*x(6)-x(3)*x(7)-x(4)*x(8) = sigma(3) ! fj(3,1) = x(5) fj(3,2) = x(6) fj(3,3) = -x(7) fj(3,4) = -x(8) fj(3,5) = x(1) fj(3,6) = x(2) fj(3,7) = -x(3) fj(3,8) = -x(4) ! ! Equation #4: ! x(1)*x(7)+x(2)*x(8)+x(3)*x(5)+x(4)*x(6) = sigma(4) ! fj(4,1) = x(7) fj(4,2) = x(8) fj(4,3) = x(5) fj(4,4) = x(6) fj(4,5) = x(3) fj(4,6) = x(4) fj(4,7) = x(1) fj(4,8) = x(2) ! ! Equation #5: ! x(1)*(x(5)**2-x(7)**2)+x(2)*(x(6)**2-x(8)**2)+x(3)*(-2.0*x(5)*x(7)) ! +x(4)*(-2.0*x(6)*x(8)) = sigma(5) ! fj(5,1) = x(5)**2-x(7)**2 fj(5,2) = x(6)**2-x(8)**2 fj(5,3) = -2.0D+00 *x(5)*x(7) fj(5,4) = -2.0D+00 *x(6)*x(8) fj(5,5) = 2.0D+00 *(x(1)*x(5)-x(3)*x(7)) fj(5,6) = 2.0D+00 *(x(2)*x(6)-x(4)*x(8)) fj(5,7) = -2.0D+00 *(x(1)*x(7)+x(3)*x(5)) fj(5,8) = -2.0D+00 *(x(2)*x(8)+x(4)*x(6)) ! ! Equation #6: ! x(1)*(2.0*x(5)*x(7))+x(2)*(2.0*x(6)*x(8))+x(3)*(x(5)**2-x(7)**2) ! +x(4)*(x(6)**2-x(8)**2) = sigma(6) ! fj(6,1) = 2.0D+00 *x(5)*x(7) fj(6,2) = 2.0D+00 *x(6)*x(8) fj(6,3) = x(5)**2-x(7)**2 fj(6,4) = x(6)**2-x(8)**2 fj(6,5) = 2.0D+00 *(x(1)*x(7)+x(3)*x(5)) fj(6,6) = 2.0D+00 *(x(2)*x(8)+x(4)*x(6)) fj(6,7) = 2.0D+00 *(x(1)*x(5)-x(3)*x(7)) fj(6,8) = 2.0D+00 *(x(2)*x(6)-x(4)*x(8)) ! ! Equation #7: ! x(1)*(x(5)*(x(5)**2-3.0*x(7)**2))+x(2)*(x(6)*(x(6)**2-3.0*x(8)**2)) ! +x(3)*(x(7)*(x(7)**2-3.0*x(5)**2))+x(4)*(x(8)*(x(8)**2-3.0*x(6)**2)) ! = sigma(7) ! fj(7,1) = x(5)*(x(5)**2-3.0D+00*x(7)**2) fj(7,2) = x(6)*(x(6)**2-3.0D+00*x(8)**2) fj(7,3) = x(7)*(x(7)**2-3.0D+00*x(5)**2) fj(7,4) = x(8)*(x(8)**2-3.0D+00*x(6)**2) fj(7,5) = 3.0D+00*(x(1)*(x(5)**2-x(7)**2)+x(3)*(-2.0D+00*x(5)*x(7))) fj(7,6) = 3.0D+00*(x(2)*(x(6)**2-x(8)**2)+x(4)*(-2.0D+00*x(6)*x(8))) fj(7,7) = -3.0D+00*(x(1)*(2.0*x(5)*x(7))+x(3)*(x(5)**2-x(7)**2)) fj(7,8) = -3.0D+00*(x(2)*(2.0*x(6)*x(8))+x(4)*(x(6)**2-x(8)**2)) ! ! Equation #8: ! x(1)*(-x(7)*(x(7)**2-3.0*x(5)**2))+x(2)*(-x(8)*(x(8)**2-3.0*x(6)**2)) ! +x(3)*(x(5)*(x(5)**2-3.0*x(7)**2))+x(4)*(x(6)*(x(6)**2-3.0*x(8)**2)) ! = sigma(8) ! fj(8,1) = -x(7)*(x(7)**2-3.0D+00*x(5)**2) fj(8,2) = -x(8)*(x(8)**2-3.0D+00*x(6)**2) fj(8,3) = x(5)*(x(5)**2-3.0D+00*x(7)**2) fj(8,4) = x(6)*(x(6)**2-3.0D+00*x(8)**2) fj(8,5) = 3.0D+00*(x(1)*(2.0D+00*x(5)*x(7))+x(3)*(x(5)**2-x(7)**2)) fj(8,6) = 3.0D+00*(x(2)*(2.0D+00*x(6)*x(8))+x(4)*(x(6)**2-x(8)**2)) fj(8,7) = 3.0D+00*(x(1)*(x(5)**2-x(7)**2)+x(3)*(-2.0D+00*x(5)*x(7))) fj(8,8) = 3.0D+00*(x(2)*(x(6)**2-x(8)**2)+x(4)*(-2.0D+00*x(6)*x(8))) return end