moments.f90

Go to the documentation of this file.

!********************************************************************************
!
!********************************************************************************

MODULE moments_module

  IMPLICIT NONE

  INTEGER :: n_nodes

  INTEGER, PARAMETER :: max_nmom = 20

  INTEGER :: n_mom

CONTAINS

  !******************************************************************************
  !
  !******************************************************************************

  SUBROUTINE wheeler_algorithm(m,xi,w)

    IMPLICIT NONE

    REAL*8, DIMENSION(n_mom), INTENT(IN) :: m
    REAL*8, DIMENSION(n_nodes), INTENT(OUT) :: xi
    REAL*8, DIMENSION(n_nodes), INTENT(OUT) :: w


    REAL*8, DIMENSION(n_nodes,n_nodes) :: jacobi
    REAL*8, DIMENSION(n_nodes) :: d
    REAL*8, DIMENSION(n_nodes-1) :: e

    REAL*8, DIMENSION(n_nodes) :: a , b
    REAL*8, DIMENSION(n_nodes+1, 2*n_nodes) :: sigma

    REAL*8, DIMENSION(n_nodes,n_nodes) :: evec

    INTEGER :: i , l , k

    REAL*8, DIMENSION(2*n_nodes-2) :: work
    INTEGER :: info
    CHARACTER*1 :: jobz
    INTEGER :: ldz


    DO l=1,2*n_nodes-2

       sigma(1,l+1) = 0.0

    END DO

    DO l=0,2*n_nodes-1

       sigma(2,l+1) = m(l+1)

    END DO

    !
    ! compute coefficients for Jacobi matrix 
    !


    a(1) = m(2) / m(1)
    b(1) = 0.d0


    DO k=1,n_nodes-1

       DO l=k,2*n_nodes-k-1

          sigma(k+2,l+1) = sigma(k+1,l+2) - a(k) * sigma(k+1,l+1) - b(k) *      &
               sigma(k,l+1)

          a(k+1) = -sigma(k+1,k+1) / sigma(k+1,k) + sigma(k+2,k+2) /            &
               sigma(k+2,k+1)

          b(k+1) = sigma(k+2,k+1) / sigma(k+1,k)

       END DO

    END DO

    !
    ! compute Jacobi matrix
    !

    DO i=1,n_nodes

       jacobi(i,i) = a(i)

       d(i) = jacobi(i,i)

    END DO

    DO i=1,n_nodes-1

       jacobi(i,i+1) = -(abs(b(i+1)))**0.5
       jacobi(i+1,i) = -(abs(b(i+1)))**0.5

       e(i) = jacobi(i,i+1)

    END DO

    !
    ! compute eigenvalues and eigenvectors
    !

    jobz = 'V'    ! compute the eigenvectors

    ldz = n_nodes

    CALL dstev(jobz, n_nodes, d, e, evec, ldz, work, info)

    !
    ! return weights
    !

    DO i=1,n_nodes

       w(i) = evec(1,i)**2 * m(1)

    END DO

    !
    ! return abscissas
    !

    DO i=1,n_nodes

       xi(i) = d(i)

    END DO

  END SUBROUTINE wheeler_algorithm

  !******************************************************************************
  !
  !******************************************************************************

  SUBROUTINE moments_correction(m,iter)

    IMPLICIT NONE

    REAL*8, DIMENSION(:), INTENT(INOUT) :: m
    INTEGER, INTENT(OUT) :: iter

    INTEGER :: iter_max

    REAL*8, DIMENSION(n_mom,n_mom) :: bk , bkk

    REAL*8, DIMENSION(n_mom,n_mom) :: d

    INTEGER :: i , j , k

    LOGICAL :: check_corr

    INTEGER :: k_star

    REAL*8, DIMENSION(n_mom) :: cos_quad_alfa

    REAL*8 :: lnck

    iter_max = 100

    !
    ! define the unitary correction matrix
    !
    DO i = 1,n_mom

       DO j = 1,n_mom

          bkk(i,j) = 0.0

       END DO

    END DO

    DO k = 1,n_mom

       DO i = 1,n_mom

          IF ( i.EQ.k ) THEN 

             bkk(i,1) = 1.d0

          ELSE

             bkk(i,1) = 0.0

          END IF

       END DO

       DO j=2,4

          DO i=1,n_mom-j+1

             bkk(i,j) = bkk(i+1,j-1) - bkk(i,j-1)

          END DO

       END DO

       DO i=1,n_mom

          bk(k,i) = bkk(i,4)

       END DO

    END DO

    check_corr = .true.
    iter = 0

    ! start the correction loop

    DO WHILE ( ( check_corr ) .AND. ( iter .LT. iter_max ) )

       check_corr = .false.

       DO i = 1,n_mom

          d(i,1) = log(m(i))

       END DO

       DO j = 2,n_mom+1

          DO i = 1,n_mom-j+1

             d(i,j) = d(i+1,j-1) - d(i,j-1)

             IF ( ( j .EQ. 3 ) .AND. ( d(i,j) .LT. 0.d0 ) ) THEN 

                check_corr = .true.

             END IF

          END DO
       END DO

       IF ( check_corr ) THEN

          iter = iter + 1
          k_star = 1

          DO k = 1,n_mom

             cos_quad_alfa(k) = ( ( dot_product( bk(k,:) , d(:,4) ) ) /         &
                  ( sqrt( dot_product( bk(k,:) , bk(k,:) ) ) *                  &
                  sqrt( dot_product( d(:,4) , d(:,4) ) ) ) ) **2.d0

             IF ( cos_quad_alfa(k) .GE. cos_quad_alfa(k_star) ) THEN

                k_star = k

             END IF

          END DO

          lnck = -( dot_product( bk(k_star,:) , d(:,4) ) ) /                    &
               dot_product( bk(k_star,:) , bk(k_star,:) )

          m(k_star) = dexp(lnck) * m(k_star)

       END IF

    END DO

  END SUBROUTINE moments_correction

  !******************************************************************************
  !
  !******************************************************************************

  REAL*8 FUNCTION beta_function(z,w)

    IMPLICIT NONE

    REAL*8, INTENT(IN):: w,z

    REAL*8:: gamln_z , gamln_w , gamln_zw

    gamln_z = gammln(z)
    gamln_w = gammln(w)
    gamln_zw = gammln(z+w)

    beta_function = dexp( gamln_z + gamln_w - gamln_zw )

    RETURN

  END FUNCTION beta_function

  !******************************************************************************
  !
  !******************************************************************************

  FUNCTION gammln(xx)

    IMPLICIT NONE

    REAL*8 :: gammln,xx

    INTEGER :: j

    REAL*8 :: ser , tmp , x , y , cof(14)

    cof(1) = 57.1562356658629235
    cof(2) = -59.5979603554754912
    cof(3) = 14.1360979747417471
    cof(4) = -0.491913816097620199
    cof(5) = 0.339946499848118887d-4
    cof(6) = 0.465236289270485756d-4
    cof(7) = -0.983744753048795646d-4
    cof(8) = 0.158088703224912494d-3
    cof(9) = -0.210264441724104883d-3
    cof(10) = 0.217439618115212643d-3
    cof(11) = -0.164318106536763890d-3
    cof(12) = 0.844182239838527433d-4
    cof(13) = -0.261908384015814087d-4
    cof(14) = .368991826595316234d-5

    IF (xx .LE. 0) THEN
       
       WRITE(*,*) "bad arg in gammln"
       RETURN

    END IF

    x = xx
    y = x

    tmp = x + 5.2421875d0
    tmp = ( x + 0.5d0 ) * log(tmp) - tmp
    
    ser = 0.999999999999997092

    DO j=1,14

       y = y+1.d0
       ser = ser + cof(j)/y

    END DO

    gammln = tmp + log(2.5066282746310005*ser/x)

    RETURN

  END FUNCTION gammln

  !******************************************************************************
  !
  !******************************************************************************


  SUBROUTINE coefficient(nmax,c)
    IMPLICIT NONE

    ! This example evaluates the binomial coefficients c(n,m) and 
    ! stores the values in the matrix (array) c(n,m)                            
    ! declarations
    INTEGER :: n,m
    INTEGER, INTENT(IN) :: nmax
    ! specify the dimensions of the two dimensional array 
    ! (matrix) c(n,m).
    INTEGER, DIMENSION(0:nmax,0:nmax), INTENT(OUT) :: c

    ! initialize array. Note how all elements are set to zero 
    ! by this one statement. Fortran90 has the facility to treat 
    ! arrays as single objects. 
    c = 0

    ! use recurrence relation c(n,m)=c(n-1,m-1) + c(n-1,m) 
    ! where m<n to generate elements

    c(0,0) = 1
    c(1,0) = 1
    c(1,1) = 1

    DO  n = 2,nmax

       c(n,0) = 1
       c(n,n) = 1

       DO  m = 1,n-1

          c(n,m) = c(n-1,m-1) + c(n-1,m)

       END DO

    END DO

    RETURN

  END SUBROUTINE coefficient

END MODULE moments_module