cnmn07 Subroutine

private subroutine cnmn07(me, ii, xbar, eps, x1, y1, x2, y2, x3, y3)

Routine to find first xbar>=epscorresponding to a real zero of a one-dimensional function by polynomial interpolation.

BY G. N. VANDERPLAATS, APRIL, 1972.

If required zero on y does not exits, or the function is ill-conditioned, xbar = eps-1.0 will be returned as an error indicator. if desired interpolation is ill-conditioned, a lower order interpolation, consistant with input data, will be attempted and ii will be changed accordingly.

Type Bound

conmin_class

Arguments

Type IntentOptional Attributes Name
class(conmin_class), intent(inout) :: me
integer, intent(inout) :: ii

CALCULATION CONTROL:

  1. 2-POINT LINEAR INTERPOLATION, GIVEN X1, Y1, X2 AND Y2.
  2. 3-POINT QUADRATIC INTERPOLATION, GIVEN X1, Y1, X2, Y2, X3 AND Y3.
real(kind=wp), intent(out) :: xbar
real(kind=wp), intent(in) :: eps

may be negative
real(kind=wp), intent(in) :: x1
real(kind=wp), intent(in) :: y1
real(kind=wp), intent(in) :: x2
real(kind=wp), intent(in) :: y2
real(kind=wp), intent(in) :: x3
real(kind=wp), intent(in) :: y3

Called by

proc~~cnmn07~~CalledByGraph proc~cnmn07 conmin_class%cnmn07 proc~cnmn06 conmin_class%cnmn06 proc~cnmn06->proc~cnmn07 proc~conmin conmin_class%conmin proc~conmin->proc~cnmn06

Source Code

    subroutine cnmn07(me, ii, xbar, eps, x1, y1, x2, y2, x3, y3)

        !!  Routine to find first `xbar>=eps `corresponding to a real zero
        !!  of a one-dimensional function by polynomial interpolation.
        !!
        !!  BY G. N. VANDERPLAATS, APRIL, 1972.
        !!
        !!  If required zero on `y` does not exits, or the function is
        !!  ill-conditioned, `xbar = eps-1.0` will be returned as an error indicator.
        !!  if desired interpolation is ill-conditioned, a lower order
        !!  interpolation, consistant with input data, will be attempted and
        !!  ii will be changed accordingly.

        class(conmin_class), intent(inout) :: me
        integer, intent(inout)     :: ii  !! CALCULATION CONTROL:
                                          !!
                                          !! 1.  2-POINT LINEAR INTERPOLATION, GIVEN X1, Y1, X2 AND Y2.
                                          !! 2.  3-POINT QUADRATIC INTERPOLATION, GIVEN X1, Y1, X2, Y2, X3 AND Y3.
        real(wp), intent(out)     :: xbar
        real(wp), intent(in)      :: eps !! may be negative
        real(wp), intent(in)      :: x1, y1, x2, y2, x3, y3

        real(wp)  :: aa, bac, bb, cc, dy, qq, x21, x31, x32, xb2, xbar1, yy
        integer   :: jj

        xbar1 = eps - 1.0_wp
        xbar = xbar1
        jj = 0
        x21 = x2 - x1
        if (abs(x21) < me%small) return
        if (ii == 2) then
            ! ------------------------------------------------------------------
            !                 ii=2: 3-point quadratic interpolation
            ! ------------------------------------------------------------------
            jj = 1
            x31 = x3 - x1
            x32 = x3 - x2
            qq = x21*x31*x32
            if (abs(qq) < me%small) return
            aa = (y1*x32 - y2*x31 + y3*x21)/qq
            if (abs(aa) >= me%small) then
                bb = (y2 - y1)/x21 - aa*(x1 + x2)
                cc = y1 - x1*(aa*x1 + bb)
                bac = bb*bb - 4.0_wp*aa*cc
                if (bac >= 0.0_wp) then
                    bac = sqrt(bac)
                    aa = 0.5_wp/aa
                    xbar = aa*(bac - bb)
                    xb2 = -aa*(bac + bb)
                    if (xbar < eps) xbar = xb2
                    if (xb2 < xbar .and. xb2 > eps) xbar = xb2
                    if (xbar < eps) xbar = xbar1
                    return
                end if
            end if
        end if

        ! ------------------------------------------------------------------
        !                  ii=1: 2-point linear interpolation
        ! ------------------------------------------------------------------
        ii = 1
        yy = y1*y2
        if (jj /= 0 .and. yy >= 0.0_wp) then
            ! interpolate between x2 and x3.
            dy = y3 - y2
            if (abs(dy) >= me%small) then
                xbar = x2 + y2*(x2 - x3)/dy
                if (xbar < eps) xbar = xbar1
                return
            end if
        end if
        dy = y2 - y1
        ! interpolate between x1 and x2.
        if (abs(dy) < me%small) return
        xbar = x1 + y1*(x1 - x2)/dy
        if (xbar < eps) xbar = xbar1

    end subroutine cnmn07