qmomk_tools_mod.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: QMOMK_TOOLS                                            C
!  Purpose: Helper functions for QMOMK implementation                  C
!                                                                      C
!  Author: Alberto Passalacqua                        Date:            C
!  Reviewer:                                          Date:            C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

#include "version.inc"

MODULE qmomk_tools

  IMPLICIT NONE

  PRIVATE

  PUBLIC :: det3
  PUBLIC :: inv3
  PUBLIC :: diag3
  PUBLIC :: cholesky3
  PUBLIC :: transpose3
  PUBLIC :: multiplymatrix3

CONTAINS

  SUBROUTINE det3 (A, det)

    IMPLICIT NONE

    DOUBLE PRECISION, INTENT(IN), DIMENSION(3,3) :: A
    DOUBLE PRECISION, INTENT(OUT) :: det

    det = a(1,1)*a(2,2)*a(3,3) + a(1,2)*a(2,3)*a(3,1) +  a(1,3)*a(2,1)*a(2,3) &
         - a(3,1)*a(2,2)*a(1,3) - a(3,2)*a(2,3)*a(1,1) - a(3,3)*a(2,1)*a(1,2)

  END SUBROUTINE det3

  SUBROUTINE inv3 (A, B)

    IMPLICIT none

    DOUBLE PRECISION, INTENT(IN), DIMENSION(3,3) :: A
    DOUBLE PRECISION, iNTENT(OUT), DIMENSION(3,3) :: B

    DOUBLE PRECISION :: detA;

    CALL det3(a, deta)

    IF (deta == 0.d0) THEN
       print *,'QMOMK: Null determinant in matrix'
       error_stop
    END IF

    b(1,1) = (a(3,3)*a(2,2) - a(3,2)*a(2,3))/deta
    b(1,2) = -(a(3,3)*a(1,2) - a(3,2)*a(1,3))/deta
    b(1,3) = (a(2,3)*a(1,2) - a(2,2)*a(1,3))/deta

    b(2,1) = -(a(3,3)*a(2,1) - a(3,1)*a(2,3))/deta
    b(2,2) = (a(3,3)*a(1,1) - a(3,1)*a(1,3))/deta
    b(2,3) = -(a(2,3)*a(1,1) - a(2,1)*a(1,3))/deta

    b(3,1) = (a(3,2)*a(2,1) - a(3,1)*a(2,2))/deta
    b(3,2) = -(a(3,2)*a(1,1) - a(3,1)*a(1,2))/deta
    b(3,3) = (a(2,2)*a(1,1) - a(2,1)*a(1,2))/deta
  END SUBROUTINE inv3

  SUBROUTINE diag3 (a, b, c, mat)

    IMPLICIT NONE

    INTEGER :: i, j
    DOUBLE PRECISION, INTENT(IN) :: a, b, c
    DOUBLE PRECISION, INTENT(OUT), DIMENSION(3,3) :: mat

    DO i = 1, 3
       DO j = 1, 3
          mat(i,j) = 0.d0
       END DO
    END DO

    mat(1,1) = a
    mat(2,2) = b
    mat(3,3) = c
  END SUBROUTINE diag3

  SUBROUTINE cholesky3 (A, L)

    IMPLICIT NONE

    DOUBLE PRECISION, INTENT(IN), DIMENSION(3,3) :: A
    DOUBLE PRECISION, INTENT(OUT), DIMENSION(3,3) :: L

    INTEGER :: i, j

    DO i = 1, 3
       DO j = 1, 3
          l(i,j) = 0.d0
       END DO
    END DO

    l(1,1) = sqrt(a(1,1))
    IF (l(1,1) == 0.d0) THEN
       print *,'Impossible to find Cholesky decomposition.'
       error_stop
    ELSE
       l(2,1) = a(2,1)/l(1,1)
       l(2,2) = sqrt(a(2,2) - (l(2,1))**2)
       l(3,1) = a(3,1)/l(1,1)

       IF (l(2,2) ==  0.d0) THEN
          print *,'Impossible to find Cholesky decomposition.'
          error_stop
       ELSE
          l(3,2) = (a(3,2) - l(3,1)*l(2,1))/l(2,2)
          l(3,3) = sqrt(a(3,3) - (l(3,1))**2 - (l(3,2))**2)
       END IF
    END IF
  END SUBROUTINE cholesky3

  SUBROUTINE transpose3 (A, T)

    IMPLICIT NONE

    DOUBLE PRECISION, INTENT(IN), DIMENSION(3,3) :: A
    DOUBLE PRECISION, INTENT(OUT), DIMENSION(3,3) :: T

    INTEGER :: i, j

    DO i = 1, 3
       DO j = 1, 3
          t(i,j) = a(j,i)
       END DO
    END DO

  END SUBROUTINE transpose3

  SUBROUTINE multiplymatrix3 (A, B, P)

    IMPLICIT NONE

    DOUBLE PRECISION, INTENT(IN), DIMENSION(3,3) :: A
    DOUBLE PRECISION, INTENT(IN), DIMENSION(3,3) :: B
    DOUBLE PRECISION, INTENT(OUT), DIMENSION(3,3) :: P

    INTEGER :: i, j, k
    DOUBLE PRECISION ::sum

    sum = 0.d0

    DO i = 1, 3
       DO j = 1, 3
          DO k = 1, 3
             sum = sum + a(i,k)*b(k,j)
          END DO
          p(i,j) = sum
          sum = 0.d0
       END DO
    END DO
  END SUBROUTINE multiplymatrix3

END MODULE qmomk_tools