lincoa_test Subroutine

public subroutine lincoa_test()

Test problem for lincoa.

Calculate the tetrahedron of least volume that encloses the points   (XP(J),YP(J),ZP(J)), J=1,2,...,NP. Our method requires the origin   to be strictly inside the convex hull of these points. There are   twelve variables that define the four faces of each tetrahedron   that is considered. Each face has the form ALPHA*X+BETA*Y+GAMMA*Z=1,   the variables X(3K-2), X(3K-1) and X(3K) being the values of ALPHA,   BETA and GAMMA for the K-th face, K=1,2,3,4. Let the set T contain   all points in three dimensions that can be reached from the origin   without crossing a face. Because the volume of T may be infinite,   the objective function is the smaller of FMAX and the volume of T,   where FMAX is set to an upper bound on the final volume initially.   There are 4*NP linear constraints on the variables, namely that each   of the given points (XP(J),YP(J),ZP(J)) shall be in T. Let XS = min   XP(J), YS = min YP(J), ZS = min ZP(J) and SS = max XP(J)+YP(J)+ZP(J),   where J runs from 1 to NP. The initial values of the variables are   X(1)=1/XS, X(5)=1/YS, X(9)=1/ZS, X(2)=X(3)=X(4)=X(6)=X(7)=X(8)=0   and X(10)=X(11)=X(12)=1/SS, which satisfy the linear constraints,   and which provide the bound FMAX=(SS-XS-YS-ZS)**3/6. Other details   of the test calculation are given below, including the choice of   the data points (XP(J),YP(J),ZP(J)), J=1,2,...,NP. The smaller final   value of the objective function in the case NPT=35 shows that the   problem has local minima.

Arguments

None

Calls

proc~~lincoa_test~~CallsGraph proc~lincoa_test lincoa_test proc~lincoa lincoa proc~lincoa_test->proc~lincoa proc~lincob lincob proc~lincoa->proc~lincob proc~prelim prelim proc~lincob->proc~prelim proc~qmstep qmstep proc~lincob->proc~qmstep proc~trstep trstep proc~lincob->proc~trstep proc~update~2 update proc~lincob->proc~update~2 proc~prelim->proc~update~2 proc~getact getact proc~trstep->proc~getact

Source Code

    subroutine lincoa_test ()

        implicit none

        real(wp) :: xp (50), yp (50), zp (50), a (12, 200), b (200), x (12)
        integer :: ia,n,np,j,iw,iprint,jcase,k,i,maxfun,npt,m
        real(wp) :: sumx,sumy,sumz,theta,fmax,rhobeg,rhoend,ss,xs,ys,zs
    !
    !     Set some constants.
    !
        real(wp),parameter :: one  = 1.0_wp
        real(wp),parameter :: two  = 2.0_wp
        real(wp),parameter :: zero = 0.0_wp
        real(wp),parameter :: pi   = 4.0_wp * atan(one)

        ia = 12
        n = 12
    !
    !     Set the data points.
    !
        np = 50
        sumx = zero
        sumy = zero
        sumz = zero
        do j = 1, np
            theta = real (j-1, wp) * pi / real (np-1, wp)
            xp (j) = cos (theta) * cos (two*theta)
            sumx = sumx + xp (j)
            yp (j) = sin (theta) * cos (two*theta)
            sumy = sumy + yp (j)
            zp (j) = sin (two*theta)
            sumz = sumz + zp (j)
        end do
        sumx = sumx / real (np, wp)
        sumy = sumy / real (np, wp)
        sumz = sumz / real (np, wp)
        do j = 1, np
            xp (j) = xp (j) - sumx
            yp (j) = yp (j) - sumy
            zp (j) = zp (j) - sumz
        end do
    !
    !     Set the linear constraints.
    !
        m = 4 * np
        do k = 1, m
            b (k) = one
            do i = 1, n
                a (i, k) = zero
            end do
        end do
        do j = 1, np
            do i = 1, 4
                k = 4 * j + i - 4
                iw = 3 * i
                a (iw-2, k) = xp (j)
                a (iw-1, k) = yp (j)
                a (iw, k) = zp (j)
            end do
        end do
    !
    !     Set the initial vector of variables. The JCASE=1,6 loop gives six
    !       different choices of NPT when LINCOA is called.
    !
        xs = zero
        ys = zero
        zs = zero
        ss = zero
        do j = 1, np
            xs = min (xs, xp(j))
            ys = min (ys, yp(j))
            zs = min (zs, zp(j))
            ss = max (ss, xp(j)+yp(j)+zp(j))
        end do
        fmax = (ss-xs-ys-zs) ** 3 / 6.0_wp
        do jcase = 1, 6
            do i = 2, 8
                x (i) = zero
            end do
            x (1) = one / xs
            x (5) = one / ys
            x (9) = one / zs
            x (10) = one / ss
            x (11) = one / ss
            x (12) = one / ss
    !
    !     Call of LINCOA, which provides the printing given at the end of this
    !       note.
    !
            npt = 5 * jcase + 10
            rhobeg = 1.0_wp
            rhoend = 1.0e-6_wp
            iprint = 1
            maxfun = 10000
            print 70, npt, rhoend
70          format (/ / 4 x, 'Output from LINCOA with  NPT =', i4, '  and  RHOEND =', 1 &
           & pd12.4)
            call lincoa(n,npt,m,a,ia,b,x,rhobeg,rhoend,iprint,maxfun,calfun)
        end do

    contains

        subroutine calfun (n, x, f)

            implicit none

            integer :: n
            real(wp) :: x (*)
            real(wp) :: f

            real(wp) :: v12,v13,v14,v23,v24,v34,del1,del2,del3,del4,temp

            f = fmax
            v12 = x (1) * x (5) - x (4) * x (2)
            v13 = x (1) * x (8) - x (7) * x (2)
            v14 = x (1) * x (11) - x (10) * x (2)
            v23 = x (4) * x (8) - x (7) * x (5)
            v24 = x (4) * x (11) - x (10) * x (5)
            v34 = x (7) * x (11) - x (10) * x (8)
            del1 = v23 * x (12) - v24 * x (9) + v34 * x (6)
            if (del1 <= zero) return
            del2 = - v34 * x (3) - v13 * x (12) + v14 * x (9)
            if (del2 <= zero) return
            del3 = - v14 * x (6) + v24 * x (3) + v12 * x (12)
            if (del3 <= zero) return
            del4 = - v12 * x (9) + v13 * x (6) - v23 * x (3)
            if (del4 <= zero) return
            temp = (del1+del2+del3+del4) ** 3 / (del1*del2*del3*del4)
            f = min (temp/6.0_wp, fmax)

        end subroutine calfun

    end subroutine lincoa_test