001_Matrix_utilities.f90

Go to the documentation of this file.

!----------------------------------------------------------------------------------------
!
! This file is part of CosmicFish.
!
! Copyright (C) 2015-2016 by the CosmicFish authors
!
! The CosmicFish code is free software;
! You can use it, redistribute it, and/or modify it under the terms
! of the GNU General Public License as published by the Free Software Foundation;
! either version 3 of the License, or (at your option) any later version.
! The full text of the license can be found in the file LICENSE at
! the top level of the CosmicFish distribution.
!
!----------------------------------------------------------------------------------------


!----------------------------------------------------------------------------------------


module matrix_utilities

    use precision

    implicit none

    private

    public sym_matrix_inversion, matrix_determinant, matrix_eigenvalues_eigenvectors_sym

contains

    ! ---------------------------------------------------------------------------------------------
    subroutine sym_matrix_inversion(A, error_code)

        implicit none

        real(dl), dimension(:,:), intent(inout) :: A
        integer , intent(out) :: error_code

        real(dl), dimension(size(A,1)) :: work   ! work array for LAPACK
        integer , dimension(size(A,1)) :: ipiv   ! pivot indices
        integer :: n, info

        ! External procedures defined in LAPACK
        external dgetrf
        external dgetri

        error_code = 0
        n = size(a,1)

        ! dgetrf computes an LU factorization of a general M-by-N matrix A
        ! using partial pivoting with row interchanges.
        call dgetrf(n, n, a, n, ipiv, info)

        if (info /= 0) then
            error_code = 1
            write(*,*) 'Matrix is numerically singular!'
            return
        end if

        ! dgetri computes the inverse of a matrix using the LU factorization
        ! computed by dgetrf.
        call dgetri(n, a, n, ipiv, work, n, info)

        if (info /= 0) then
            error_code = 2
            write(*,*) 'Matrix inversion failed!'
            return
        end if

    end subroutine sym_matrix_inversion

    ! ---------------------------------------------------------------------------------------------
    subroutine matrix_determinant(A, det, error_code)

        implicit none

        real(dl), dimension(:,:), intent(in) :: A
        real(dl), intent(out) :: det
        integer , intent(out) :: error_code

        integer  :: N, i, info, AllocateStatus
        real(dl) :: sgn
        integer , dimension(:), allocatable :: ipiv   ! pivot indices
        ! External procedures defined in LAPACK
        external dgetrf

        n          = size(a,1)
        error_code = 0

        allocate( ipiv(n), stat = allocatestatus )
        if (allocatestatus /= 0) stop "matrix_determinant: allocation failed. Not enough memory."

        ipiv=0

        ! dgetrf computes an LU factorization of a general M-by-N matrix A
        ! using partial pivoting with row interchanges.
        call dgetrf(n, n, a, n, ipiv, info)

        if (info /= 0) then
            error_code = 1
            write(*,*) 'Matrix is numerically singular!'
            stop
        end if

        det = 1._dl

        do i=1, n
            det = det*a(i,i)
        end do

        sgn = 1._dl

        do i=1, n
            if(ipiv(i)/=i)then
                sgn=-sgn
            end if
        end do

        det=sgn*det

        return

    end subroutine matrix_determinant

    ! ---------------------------------------------------------------------------------------------
    subroutine matrix_eigenvalues_eigenvectors_sym(A, eigenvalues, eigenvectors, error_code)

        implicit none

        real(dl), dimension(:,:), intent(in)  :: A
        real(dl), dimension(:)  , intent(out) :: eigenvalues
        real(dl), dimension(:,:), intent(out) :: eigenvectors
        integer , intent(out) :: error_code

        eigenvectors = 0._dl
        eigenvalues  = 0._dl
        error_code   = 0

        stop 'matrix_eigenvalues_eigenvectors_sym not yet implemented.'

    end subroutine matrix_eigenvalues_eigenvectors_sym

    ! ---------------------------------------------------------------------------------------------

end module matrix_utilities