LCOV - code coverage report
Current view: | top level - physics/clubb/src/CLUBB_core - matrix_operations.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 0 | 122 | 0.0 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 0 | 8 | 0.0 %

              Line data    Source code
           1             : !-----------------------------------------------------------------------
           2             : ! $Id$
           3             : !===============================================================================
           4             : module matrix_operations
           5             :
           6             :   implicit none
           7             :
           8             :
           9             :   public :: symm_covar_matrix_2_corr_matrix, Cholesky_factor, &
          10             :     row_mult_lower_tri_matrix, print_lower_triangular_matrix, &
          11             :     get_lower_triangular_matrix, set_lower_triangular_matrix, &
          12             :     mirror_lower_triangular_matrix
          13             :
          14             :   private :: Symm_matrix_eigenvalues
          15             :
          16             :   private ! Default scope
          17             :
          18             :   contains
          19             :
          20             : !-----------------------------------------------------------------------
          21           0 :   subroutine symm_covar_matrix_2_corr_matrix( ndim, covar, &
          22           0 :                                               corr )
          23             :
          24             : ! Description:
          25             : !   Convert a matrix of covariances in to a matrix of correlations.
          26             : !   This only does the computation the lower triangular portion of the
          27             : !   matrix.
          28             : ! References:
          29             : !   None
          30             : !-----------------------------------------------------------------------
          31             :
          32             :     use clubb_precision, only: &
          33             :         core_rknd ! double precision
          34             :
          35             :     implicit none
          36             :
          37             :     ! External
          38             :     intrinsic :: sqrt
          39             :
          40             :     ! Input Variables
          41             :     integer, intent(in) :: ndim
          42             :
          43             :     real( kind = core_rknd ), dimension(ndim,ndim), intent(in) :: &
          44             :       covar ! Covariance Matrix [units vary]
          45             :
          46             :     ! Output Variables
          47             :     real( kind = core_rknd ), dimension(ndim,ndim), intent(out) :: &
          48             :       corr ! Correlation Matrix [-]
          49             :
          50             :     ! Local Variables
          51             :     integer :: i, j
          52             :
          53             :     ! ---- Begin Code ----
          54             :
          55           0 :     corr = 0._core_rknd ! Initialize to 0
          56             :
          57           0 :     do i = 1, ndim
          58           0 :       do j = 1, i
          59           0 :         corr(i,j) = covar(i,j) / sqrt( covar(i,i) * covar(j,j) )
          60             :       end do
          61             :     end do
          62             :
          63           0 :     return
          64             :   end subroutine symm_covar_matrix_2_corr_matrix
          65             : !-----------------------------------------------------------------------
          66           0 :   subroutine row_mult_lower_tri_matrix( ndim, xvector, tmatrix_in, tmatrix_out )
          67             :
          68             : ! Description:
          69             : !   Do a row-wise multiply of the elements of a lower triangular matrix.
          70             : ! References:
          71             : !   None
          72             : !-----------------------------------------------------------------------
          73             :
          74             :     use clubb_precision, only: &
          75             :         core_rknd ! double precision
          76             :
          77             :     implicit none
          78             :
          79             :
          80             :     ! Input Variables
          81             :     integer, intent(in) :: ndim
          82             :
          83             :     real( kind = core_rknd ), dimension(ndim), intent(in) :: &
          84             :       xvector ! Factors to be multiplied across a row [units vary]
          85             :
          86             :     ! Input Variables
          87             :     real( kind = core_rknd ), dimension(ndim,ndim), intent(in) :: &
          88             :       tmatrix_in ! nxn matrix (usually a correlation matrix) [units vary]
          89             :
          90             :     ! Output Variables
          91             :     real( kind = core_rknd ), dimension(ndim,ndim), intent(inout) :: &
          92             :       tmatrix_out ! nxn matrix (usually a covariance matrix) [units vary]
          93             :
          94             :     ! Local Variables
          95             :     integer :: i, j
          96             :
          97             :     ! ---- Begin Code ----
          98             :
          99           0 :     do i = 1, ndim
         100           0 :       do j = 1, i
         101           0 :         tmatrix_out(i,j) = tmatrix_in(i,j) * xvector(i)
         102             :       end do
         103             :     end do
         104             :
         105           0 :     return
         106             :   end subroutine row_mult_lower_tri_matrix
         107             :
         108             : !-------------------------------------------------------------------------------
         109           0 :   subroutine Cholesky_factor( ndim, a_input, a_scaling, a_Cholesky, l_scaled )
         110             : !  Description:
         111             : !    Create a Cholesky factorization of a_input.
         112             : !    If the factorization fails we use a modified a_input matrix and attempt
         113             : !    to factorize again.
         114             : !
         115             : !  References:
         116             : !    <http://www.netlib.org/lapack/explore-html/a00868.html> dpotrf
         117             : !    <http://www.netlib.org/lapack/explore-html/a00860.html> dpoequ
         118             : !    <http://www.netlib.org/lapack/explore-html/a00753.html> dlaqsy
         119             : !-------------------------------------------------------------------------------
         120             :     use error_code, only: &
         121             :         clubb_at_least_debug_level ! Procedure
         122             :
         123             :     use constants_clubb, only: &
         124             :         fstderr ! Constant
         125             :
         126             :     use clubb_precision, only: &
         127             :         core_rknd
         128             :
         129             :     use lapack_interfaces, only: &
         130             :         lapack_potrf, &   ! Procedures
         131             :         lapack_poequ, &
         132             :         lapack_laqsy
         133             :
         134             :     implicit none
         135             :
         136             :     ! Constant Parameters
         137             :     integer, parameter :: itermax = 10 ! Max iterations of the modified method
         138             :
         139             :     real( kind = core_rknd), parameter :: d_coef = 0.1_core_rknd
         140             :        ! Coefficient applied if the decomposition doesn't work
         141             :
         142             :     ! Input Variables
         143             :     integer, intent(in) :: ndim
         144             :
         145             :     real( kind = core_rknd ), dimension(ndim,ndim), intent(in) :: a_input
         146             :
         147             :     ! Output Variables
         148             :     real( kind = core_rknd ), dimension(ndim), intent(out) :: a_scaling
         149             :
         150             :     real( kind = core_rknd ), dimension(ndim,ndim), intent(out) :: a_Cholesky
         151             :
         152             :     logical, intent(out) :: l_scaled
         153             :
         154             :     ! Local Variables
         155           0 :     real( kind = core_rknd ), dimension(ndim) :: a_eigenvalues
         156           0 :     real( kind = core_rknd ), dimension(ndim,ndim) ::  a_corr, a_scaled
         157             :
         158             :     real( kind = core_rknd ) :: tau, d_smallest
         159             :
         160             :     real( kind = core_rknd ) :: amax, scond
         161             :     integer :: info
         162             :     integer :: i, j, iter
         163             :
         164             :     character :: equed
         165             :
         166             :     ! ---- Begin code ----
         167             :
         168           0 :     a_scaled = a_input ! Copy input array into output array
         169             :
         170             : !   do i = 1, n
         171             : !     do j = 1, n
         172             : !       write(6,'(e10.3)',advance='no') a(i,j)
         173             : !     end do
         174             : !     write(6,*) ""
         175             : !   end do
         176             : !   pause
         177             :
         178           0 :     equed = 'N'
         179             :
         180             :     ! Compute scaling for a_input, using Lapack routine spoequ for single precision,
         181             :     ! or dpoequ for double precision
         182           0 :     call lapack_poequ( ndim, a_input, ndim, a_scaling, scond, amax, info )
         183             :
         184           0 :     if ( info == 0 ) then
         185             :       ! Apply scaling to a_input, using Lapack routine slaqsy for single precision,
         186             :       ! or dlaqsy for double precision
         187           0 :       call lapack_laqsy( 'Lower', ndim, a_scaled, ndim, a_scaling, scond, amax, equed )
         188             :     end if
         189             :
         190             :     ! Determine if scaling was necessary
         191           0 :     if ( equed == 'Y' ) then
         192           0 :       l_scaled = .true.
         193           0 :       a_Cholesky = a_scaled
         194             :     else
         195           0 :       l_scaled = .false.
         196           0 :       a_Cholesky = a_input
         197             :     end if
         198             :
         199           0 :     do iter = 1, itermax
         200             :
         201             :       ! Lapack Cholesky factorization, spotrf for single or dpotrf for double precision
         202           0 :       call lapack_potrf( 'Lower', ndim, a_Cholesky, ndim, info )
         203             :
         204           0 :       select case( info )
         205             :       case( :-1 )
         206             :         write(fstderr,*) "Cholesky_factor " // &
         207           0 :           " illegal value for argument ", -info
         208           0 :         error stop
         209             :       case( 0 )
         210             :         ! Success!
         211           0 :         if ( clubb_at_least_debug_level( 1 ) .and. iter > 1 ) then
         212           0 :           write(fstderr,*) "a_factored (worked)="
         213           0 :           do i = 1, ndim
         214           0 :             do j = 1, i
         215           0 :               write(fstderr,'(g10.3)',advance='no') a_Cholesky(i,j)
         216             :             end do
         217           0 :             write(fstderr,*) ""
         218             :           end do
         219             :         end if
         220           0 :         exit
         221             :       case( 1: )
         222           0 :         if ( clubb_at_least_debug_level( 1 ) ) then
         223             :           ! This shouldn't happen now that the s and t Mellor(chi/eta) elements have been
         224             :           ! modified to never be perfectly correlated, but it's here just in case.
         225             :           ! -dschanen 10 Sept 2010
         226           0 :           write(fstderr,*) "Cholesky_factor: leading minor of order ", &
         227           0 :             info, " is not positive definite."
         228           0 :           write(fstderr,*) "factorization failed."
         229           0 :           write(fstderr,*) "a_input="
         230           0 :           do i = 1, ndim
         231           0 :             do j = 1, i
         232           0 :               write(fstderr,'(g10.3)',advance='no') a_input(i,j)
         233             :             end do
         234           0 :             write(fstderr,*) ""
         235             :           end do
         236           0 :           write(fstderr,*) "a_Cholesky="
         237           0 :           do i = 1, ndim
         238           0 :             do j = 1, i
         239           0 :               write(fstderr,'(g10.3)',advance='no') a_Cholesky(i,j)
         240             :             end do
         241           0 :             write(fstderr,*) ""
         242             :           end do
         243             :         end if
         244             :
         245           0 :         if ( clubb_at_least_debug_level( 2 ) ) then
         246             :           call Symm_matrix_eigenvalues( ndim, a_input, & ! intent(in)
         247           0 :                                         a_eigenvalues ) ! intent(out)
         248           0 :           write(fstderr,*) "a_eigenvalues="
         249           0 :           do i = 1, ndim
         250           0 :             write(fstderr,'(g10.3)',advance='no') a_eigenvalues(i)
         251             :           end do
         252           0 :           write(fstderr,*) ""
         253             :
         254             :           call symm_covar_matrix_2_corr_matrix( ndim, a_input, & ! intent(in)
         255           0 :                                                 a_corr ) ! intent(out)
         256           0 :           write(fstderr,*) "a_correlations="
         257           0 :           do i = 1, ndim
         258           0 :             do j = 1, i
         259           0 :               write(fstderr,'(g10.3)',advance='no') a_corr(i,j)
         260             :             end do
         261           0 :             write(fstderr,*) ""
         262             :           end do
         263             :         end if
         264             :
         265           0 :         if ( iter == itermax ) then
         266           0 :           write(fstderr,*) "iteration =", iter, "itermax =", itermax
         267           0 :           write(fstderr,*) "Fatal error in Cholesky_factor"
         268           0 :         else if ( clubb_at_least_debug_level( 1 ) ) then
         269             :           ! Adding a STOP statement to prevent this problem from slipping under
         270             :           ! the rug.
         271           0 :           write(fstderr,*) "Fatal error in Cholesky_factor"
         272           0 :           write(fstderr,*) "Attempting to modify matrix to allow factorization."
         273             :         end if
         274             :
         275           0 :         if ( l_scaled ) then
         276           0 :           a_Cholesky = a_scaled
         277             :         else
         278           0 :           a_Cholesky = a_input
         279             :         end if
         280             :         ! The number used for tau here is case specific to the Sigma covariance
         281             :         ! matrix in the latin hypercube code and is not at all general.
         282             :         ! Tau should be a number that is small relative to the other diagonal
         283             :         ! elements of the matrix to have keep the error caused by modifying 'a' low.
         284             :         ! -dschanen 30 Aug 2010
         285           0 :         d_smallest = a_Cholesky(1,1)
         286           0 :         do i = 2, ndim
         287           0 :           if ( d_smallest > a_Cholesky(i,i) ) d_smallest = a_Cholesky(i,i)
         288             :         end do
         289             :         ! Use the smallest element * d_coef * iteration
         290           0 :         tau = d_smallest * d_coef * real( iter, kind=core_rknd )
         291             :
         292             : !       print *, "tau =", tau, "d_smallest = ", d_smallest
         293             :
         294           0 :         do i = 1, ndim
         295           0 :           do j = 1, ndim
         296           0 :             if ( i == j ) then
         297           0 :               a_Cholesky(i,j) = a_Cholesky(i,j) + tau ! Add tau to the diagonal
         298             :             else
         299           0 :               a_Cholesky(i,j) = a_Cholesky(i,j)
         300             :             end if
         301             :           end do
         302             :         end do
         303             :
         304           0 :         if ( clubb_at_least_debug_level( 2 ) ) then
         305             :           call Symm_matrix_eigenvalues( ndim, a_Cholesky, & ! intent(in)
         306           0 :                                         a_eigenvalues ) ! intent(out)
         307           0 :           write(fstderr,*) "a_modified eigenvalues="
         308           0 :           do i = 1, ndim
         309           0 :             write(fstderr,'(e10.3)',advance='no') a_eigenvalues(i)
         310             :           end do
         311           0 :           write(fstderr,*) ""
         312             :         end if
         313             :
         314             :       end select ! info
         315             :     end do ! 1..itermax
         316             :
         317           0 :     return
         318             :   end subroutine Cholesky_factor
         319             :
         320             : !----------------------------------------------------------------------
         321           0 :   subroutine Symm_matrix_eigenvalues( ndim, a_input, &
         322           0 :                                       a_eigenvalues )
         323             :
         324             : !   Description:
         325             : !     Computes the eigevalues of a_input
         326             : !
         327             : !   References:
         328             : !     None
         329             : !-----------------------------------------------------------------------
         330             :
         331             :     use constants_clubb, only: &
         332             :         fstderr ! Constant
         333             :
         334             :     use clubb_precision, only: &
         335             :         core_rknd ! double precision
         336             :
         337             :     use lapack_interfaces, only: &
         338             :         lapack_syev       ! Procedure
         339             :
         340             :     implicit none
         341             :
         342             :     ! Parameters
         343             :     integer, parameter :: &
         344             :       lwork = 180 ! This is the optimal value I obtained for an n of 5 -dschanen 31 Aug 2010
         345             :
         346             :     ! Input Variables
         347             :     integer, intent(in) :: ndim
         348             :
         349             :     real( kind = core_rknd ), dimension(ndim,ndim), intent(in) :: a_input
         350             :
         351             :     ! Output Variables
         352             :     real( kind = core_rknd ), dimension(ndim), intent(out) :: a_eigenvalues
         353             :
         354             :     ! Local Variables
         355           0 :     real( kind = core_rknd ), dimension(ndim,ndim) :: a_scratch
         356             :
         357             :     real( kind = core_rknd ), dimension(lwork) :: work
         358             :
         359             :     integer :: info
         360             : !   integer :: i, j
         361             :     ! ---- Begin code ----
         362             :
         363           0 :     a_scratch = a_input
         364             :
         365             : !   do i = 1, ndim
         366             : !     do j = 1, ndim
         367             : !       write(6,'(e10.3)',advance='no') a(i,j)
         368             : !     end do
         369             : !     write(6,*) ""
         370             : !   end do
         371             : !   pause
         372             :
         373             :     ! Lapack routine for computing eigenvalues and, optionally, eigenvectors. ssyev for
         374             :     ! single precision  or dsyev for double precision
         375             :     call lapack_syev( 'No eigenvectors', 'Lower', ndim, a_scratch, ndim, &
         376           0 :                       a_eigenvalues, work, lwork, info )
         377             :
         378           0 :     select case( info )
         379             :     case( :-1 )
         380             :       write(fstderr,*) "Symm_matrix_eigenvalues:" // &
         381           0 :         " illegal value for argument ", -info
         382             :     case( 0 )
         383             :       ! Success!
         384             :
         385             :     case( 1: )
         386           0 :       write(fstderr,*) "Symm_matrix_eigenvalues: Algorithm failed to converge."
         387             :     end select
         388             :
         389           0 :     return
         390             :   end subroutine Symm_matrix_eigenvalues
         391             : !-------------------------------------------------------------------------------
         392           0 :   subroutine set_lower_triangular_matrix( pdf_dim, index1, index2, xpyp, &
         393           0 :                                           matrix )
         394             : ! Description:
         395             : !   Set a value for the lower triangular portion of a matrix.
         396             : ! References:
         397             : !   None
         398             : !-------------------------------------------------------------------------------
         399             :
         400             :     use clubb_precision, only: &
         401             :       core_rknd ! user defined precision
         402             :
         403             :     implicit none
         404             :
         405             :     ! External
         406             :     intrinsic :: max, min
         407             :
         408             :     ! Input Variables
         409             :     integer, intent(in) :: &
         410             :       pdf_dim, & ! Number of variates
         411             :       index1, index2 ! Indices for 2 variates (the order doesn't matter)
         412             :
         413             :     real( kind = core_rknd ), intent(in) :: &
         414             :       xpyp ! Value for the matrix (usually a correlation or covariance) [units vary]
         415             :
         416             :     ! Input/Output Variables
         417             :     real( kind = core_rknd ), dimension(pdf_dim,pdf_dim), intent(inout) :: &
         418             :       matrix ! The lower triangular matrix
         419             :
         420             :     integer :: i,j
         421             :
         422             :     ! ---- Begin Code ----
         423             :
         424             :     ! Reverse these to set the values of upper triangular matrix
         425           0 :     i = max( index1, index2 )
         426           0 :     j = min( index1, index2 )
         427             :
         428           0 :     if( i > 0 .and. j > 0 ) then
         429           0 :       matrix(i,j) = xpyp
         430             :     end if
         431             :
         432           0 :     return
         433             :   end subroutine set_lower_triangular_matrix
         434             : !-------------------------------------------------------------------------------
         435             :
         436             : !-------------------------------------------------------------------------------
         437           0 :   subroutine get_lower_triangular_matrix( pdf_dim, index1, index2, matrix, &
         438             :                                           xpyp )
         439             : ! Description:
         440             : !   Returns a value from the lower triangular portion of a matrix.
         441             : ! References:
         442             : !   None
         443             : !-------------------------------------------------------------------------------
         444             :
         445             :     use clubb_precision, only: &
         446             :         core_rknd
         447             :
         448             :     implicit none
         449             :
         450             :     ! External
         451             :     intrinsic :: max, min
         452             :
         453             :     ! Input Variables
         454             :     integer, intent(in) :: &
         455             :       pdf_dim, & ! Number of variates
         456             :       index1, index2 ! Indices for 2 variates (the order doesn't matter)
         457             :
         458             :     ! Input/Output Variables
         459             :     real( kind = core_rknd ), dimension(pdf_dim,pdf_dim), intent(in) :: &
         460             :       matrix ! The covariance matrix
         461             :
         462             :     real( kind = core_rknd ), intent(out) :: &
         463             :       xpyp ! Value from the matrix (usually a correlation or covariance) [units vary]
         464             :
         465             :     integer :: i,j
         466             :
         467             :     ! ---- Begin Code ----
         468             :
         469             :     ! Reverse these to set the values of upper triangular matrix
         470           0 :     i = max( index1, index2 )
         471           0 :     j = min( index1, index2 )
         472             :
         473           0 :     xpyp = matrix(i,j)
         474             :
         475           0 :     return
         476             :   end subroutine get_lower_triangular_matrix
         477             :
         478             : !-----------------------------------------------------------------------
         479           0 :   subroutine print_lower_triangular_matrix( iunit, ndim, matrix )
         480             :
         481             : ! Description:
         482             : !   Print the values of lower triangular matrix to a file or console.
         483             :
         484             : ! References:
         485             : !   None
         486             : !-----------------------------------------------------------------------
         487             :
         488             :     use clubb_precision, only: &
         489             :         core_rknd ! Variable(s)
         490             :
         491             :     implicit none
         492             :
         493             :     ! Input Variables
         494             :     integer, intent(in) :: &
         495             :       iunit, & ! File I/O logical unit (usually 6 for stdout and 0 for stderr)
         496             :       ndim     ! Dimension of the matrix
         497             :
         498             :     real( kind = core_rknd ), dimension(ndim,ndim), intent(in) :: &
         499             :       matrix ! Lower triangular matrix [units vary]
         500             :
         501             :     ! Local Variables
         502             :     integer :: i, j
         503             :
         504             :     ! ---- Begin Code ----
         505             :
         506           0 :     do i = 1, ndim
         507           0 :       do j = 1, i
         508           0 :         write(iunit,fmt='(g15.6)',advance='no') matrix(i,j)
         509             :       end do
         510           0 :       write(iunit,fmt=*) "" ! newline
         511             :     end do
         512             :
         513           0 :     return
         514             :   end subroutine print_lower_triangular_matrix
         515             :
         516             :   !-----------------------------------------------------------------------
         517           0 :   subroutine mirror_lower_triangular_matrix( nvars, &
         518           0 :                                              matrix )
         519             :
         520             :   ! Description:
         521             :   !   Mirrors the elements of a lower triangular matrix to the upper
         522             :   !   triangle so that it is symmetric.
         523             :
         524             :   ! References:
         525             :   !   None
         526             :   !-----------------------------------------------------------------------
         527             :
         528             :     use clubb_precision, only: &
         529             :         core_rknd  ! Constant
         530             :
         531             :     implicit none
         532             :
         533             :     ! Input Variables
         534             :     integer, intent(in) :: &
         535             :       nvars ! Number of variables in each dimension of square matrix
         536             :
         537             :     ! Input/Output Variables
         538             :     real( kind = core_rknd ), dimension(nvars,nvars), intent(inout) :: &
         539             :       matrix  ! Lower triangluar square matrix
         540             :
         541             :     ! Local Variables
         542             :     integer :: row, col
         543             :
         544             :   !-----------------------------------------------------------------------
         545             :
         546             :     !----- Begin Code -----
         547             :
         548           0 :     if ( nvars > 1 ) then
         549             :
         550           0 :       do col=2, nvars
         551           0 :         do row=1, col-1
         552             :
         553           0 :           matrix(row,col) = matrix(col,row)
         554             :
         555             :         end do
         556             :       end do
         557             :
         558             :     end if ! nvars > 1
         559             :
         560           0 :     return
         561             :
         562             :   end subroutine mirror_lower_triangular_matrix
         563             :   !-----------------------------------------------------------------------
         564             :
         565             : end module matrix_operations

Generated by: LCOV version 1.14