common.f90

一个基于打靶法的最优控制求解软件 求解过程中采用参数延续算法

!!****************************************
!!* Simplicial package v1.4              *
!!* Auxilliary subroutines and functions *
!!* Author: Pierre Martinon              *
!!* ENSEEIHT-IRIT, UMR CNRS 5505         *
!!* CMAP-INRIA FUTURS, UMR CNRS 7641     *
!!* 12/2006                              *
!!****************************************

module CommonFuns

contains

  !!********************
  !!* Matrix inversion *
  !!********************
  Subroutine Invert(dim,mat,invmat,ifail,verbose)
    implicit none
    integer, intent(in) :: dim, verbose
    integer, intent(out) :: ifail
    real(kind=8), dimension(dim,dim), intent(in) :: mat
    real(kind=8), dimension(dim,dim), intent(out) :: invmat

    !Local declarations
    integer :: lwa, info
    integer, dimension(dim) :: ipiv
    real(kind=8), dimension(dim*dim) :: wa

    ifail = 0
    lwa = dim * dim
    ipiv = 0
    info = 0
    wa = 0d0

    invmat = mat
    call dgetrf(dim, dim, invmat, dim, ipiv, info)
    if (info /= 0) then 
       if (verbose > 1) write (0,*) 'ERROR: DGETRF >>> info =', info
       !write (0,*) mat
       ifail = 2
    end if

    call dgetri(dim, invmat, dim, ipiv, wa, lwa, info)
    if (info /= 0) then
       if (verbose > 1) write (0,*) 'ERROR: DGETRI >>> info =', info
       ifail = 3
    end if


  end subroutine Invert


  Subroutine Ecart(n,tab,moy,dist,sigma)
    implicit none
    integer, intent(in) :: n
    integer, dimension(n), intent(in) :: tab
    real(kind=8), intent(out) :: moy, sigma
    real(kind=8), dimension(n), intent(out) :: dist

    !Local declarations
    integer :: i

    moy = 0d0
    do i=1,n
       moy = moy + tab(i)
    end do
    moy = moy / n

    dist = tab - moy

    sigma = sqrt(sum(dist*dist)) / moy


  end Subroutine Ecart


  !!**********************
  !!* min tab (minloc)   *
  !!**********************
  function mintab(dim,tab) result(loc)
    implicit none
    integer, intent(in) :: dim
    real(kind=8), dimension(dim), intent(in) :: tab
    integer :: loc

    !Local declarations
    integer :: i
    real(kind=8) :: minval

    minval = tab(1)
    loc = 1
    do i=2,dim
       if (tab(i) < minval) then
          loc = i
          minval = tab(i)
       end if
    end do

  end function mintab


!!$  function minval(dim,tab) result(val)
!!$    implicit none
!!$    integer, intent(in) :: dim
!!$    real(kind=8), dimension(dim), intent(in) :: tab
!!$    real(kind=8) :: val
!!$
!!$    !Local declarations
!!$    integer :: i
!!$
!!$    val = tab(1)
!!$    do i=2,dim
!!$       if (tab(i) < val) then
!!$          val = tab(i)
!!$       end if
!!$    end do
!!$
!!$  end function minval

  !!**********************
  !!* factorial function *
  !!**********************
  function factorial(n) result(fact)
    implicit none
    integer, intent(in) :: n
    integer :: fact

    !Local declarations
    integer :: i

    fact = 1
    do i = 2, n
       fact = fact * i
    end do

  end function factorial


  !!************************************
  !!* Data scaling internal subroutine *
  !!************************************
  Subroutine SubScale(dim,value,scaling,scalemode)
    implicit none
    integer, intent(in) :: dim, scalemode
    real(kind=8), dimension(dim), intent(inout) :: value
    real(kind=8), dimension(dim), intent(inout) :: scaling

    !Local variables
    integer :: i
    real(kind=8) :: lowerscale, upperscale, scaleps

    scaleps = 1d-8


    !value: unscaled value to be scaled to [lower,upper] in abs val
    !scaling: corresponding scaling
    lowerscale = 1d-1
    upperscale = 1d0
    
    !write(0,*) 'Scaling Mode',scalemode
    !write(0,*) 'incoming',value, scaling

    do i=1,dim
       if (abs(value(i)) <= scaleps) then
          scaling(i) = 1d0
       else
          do while (abs(value(i)) >= upperscale)
             value(i) = value(i) / 1d1
             scaling(i) = scaling(i) / 1d1
             !write(0,*) 'Composante',i,'valeur',value(i),'scal',scaling(i)
          end do
          do while (abs(value(i)) < lowerscale)
             value(i) = value(i) * 1d1
             scaling(i) = scaling(i) * 1d1
             !write(0,*) 'Composante',i,'valeur',value(i),'scal',scaling(i)
          end do
       end if
    end do

    !write(0,*) 'outgoing',value, scaling

  end subroutine Subscale

end module CommonFuns