LCOV - code coverage report
Current view: top level - physics/clubb/src/CLUBB_core - lapack_interfaces.F90 (source / functions) Hit Total Coverage
Test: coverage.info Lines: 8 62 12.9 %
Date: 2024-12-17 17:57:11 Functions: 3 20 15.0 %

          Line data    Source code
       1             : !-----------------------------------------------------------------------
       2             : ! $Id$
       3             : !===============================================================================
       4             : module lapack_interfaces
       5             :     
       6             :   ! This module acts as an interface to Lapack routines. It's main purpose
       7             :   ! is to make interfaces available to clubb that can handle both 
       8             :   ! single and doulbe precision. This may be compiled along with Lapack
       9             :   ! source code, or along with a linked Lapack library such as MKL. 
      10             :     
      11             :   implicit none
      12             :   
      13             :   public :: lapack_gbsv, lapack_gbsvx, &
      14             :             lapack_gtsv, lapack_gtsvx, &
      15             :             lapack_isnan, lapack_potrf, &
      16             :             lapack_poequ, lapack_laqsy, &
      17             :             lapack_syev, lapack_trmv
      18             :             
      19             :   private :: &
      20             :     dgbsv_wrap, sgbsv_wrap, &
      21             :     dgbsvx_wrap, sgbsvx_wrap, &
      22             :     dgtsv_wrap, sgtsv_wrap, &
      23             :     dgtsvx_wrap, sgtsvx_wrap, &
      24             :     disnan_wrap, sisnan_wrap, &
      25             :     dpotrf_wrap, spotrf_wrap, &
      26             :     dpoequ_wrap, spoequ_wrap, &
      27             :     dlaqsy_wrap, slaqsy_wrap, &
      28             :     dsyev_wrap, ssyev_wrap, &
      29             :     dtrmv_wrap, strmv_wrap
      30             :     
      31             :   ! Interface for Lapack general band solver, single or double precision
      32             :   interface lapack_gbsv
      33             :       module procedure dgbsv_wrap
      34             :       module procedure sgbsv_wrap
      35             :   end interface
      36             : 
      37             :   ! Interface for Lapack general band solver, expert version, single or double precision
      38             :   interface lapack_gbsvx
      39             :       module procedure dgbsvx_wrap
      40             :       module procedure sgbsvx_wrap
      41             :   end interface
      42             :   
      43             :   ! Interface for Lapack tridiagonal matrix solver, single or double precision
      44             :   interface lapack_gtsv
      45             :       module procedure dgtsv_wrap
      46             :       module procedure sgtsv_wrap
      47             :   end interface
      48             :   
      49             :   ! Interface for Lapack tridiagonal matrix solver, expert version, single or double precision
      50             :   interface lapack_gtsvx
      51             :       module procedure dgtsvx_wrap
      52             :       module procedure sgtsvx_wrap
      53             :   end interface
      54             :   
      55             :   ! Interface for Lapack nan check, single or double precision
      56             :   interface lapack_isnan
      57             :       module procedure disnan_wrap
      58             :       module procedure sisnan_wrap
      59             :   end interface
      60             :   
      61             :   ! Interface for Lapack's Cholesky factorization of a real symmetric positive definite
      62             :   ! matrix, single or double precision
      63             :   interface lapack_potrf
      64             :       module procedure dpotrf_wrap
      65             :       module procedure spotrf_wrap
      66             :   end interface
      67             :   
      68             :   ! Interface for Lapack routine to compute row and column scalings intended to 
      69             :   ! equilibriate a symmetric positive definite matrix, single or doulbe precision
      70             :   interface lapack_poequ
      71             :       module procedure dpoequ_wrap
      72             :       module procedure spoequ_wrap
      73             :   end interface
      74             :   
      75             :   ! Interface for Lapack routine to equilibriate a symmetric matrix, single or double precision
      76             :   interface lapack_laqsy
      77             :       module procedure dlaqsy_wrap
      78             :       module procedure slaqsy_wrap
      79             :   end interface
      80             :   
      81             :   ! Interface for Lapack routine to compute all eigenvalues and, optionally, eigenvectors
      82             :   ! of a real symmetric matrix, single or double precision
      83             :   interface lapack_syev
      84             :       module procedure dsyev_wrap
      85             :       module procedure ssyev_wrap
      86             :   end interface
      87             :   
      88             :   ! Interface for Lapack routines to performe one of the following matrix-vector operations
      89             :   !     x := A*x,   or   x := A**T*x,
      90             :   ! where A is an upper or lower triangular matrix, single or double precision
      91             :   interface lapack_trmv
      92             :       module procedure dtrmv_wrap
      93             :       module procedure strmv_wrap
      94             :   end interface
      95             :       
      96             :       
      97             :   private ! Set Default Scope
      98             :   
      99             :   contains
     100             :       
     101             :   ! ==================== General Band Solver Wrappers ====================
     102             :   
     103             :   ! Double precision wrapper
     104   280519200 :   subroutine dgbsv_wrap( n, kl, ku, nrhs, ab, ldab, ipiv, b, ldb, info )
     105             :       
     106             :       implicit none
     107             :       
     108             :       external :: dgbsv
     109             :       
     110             :       integer            info, kl, ku, ldab, ldb, n, nrhs
     111             :       integer            ipiv( * )
     112             :       double precision   ab( ldab, * ), b( ldb, * )
     113             :       
     114   280519200 :       call dgbsv( n, kl, ku, nrhs, ab, ldab, ipiv, b, ldb, info )
     115             :       
     116   280519200 :   end subroutine dgbsv_wrap
     117             :   
     118             :   ! Single precision wrapper
     119           0 :   subroutine sgbsv_wrap( n, kl, ku, nrhs, ab, ldab, ipiv, b, ldb, info )
     120             :       
     121             :       implicit none
     122             :       
     123             :       external :: sgbsv
     124             :       
     125             :       integer   info, kl, ku, ldab, ldb, n, nrhs
     126             :       integer   ipiv( * )
     127             :       real      ab( ldab, * ), b( ldb, * )
     128             :       
     129           0 :       call sgbsv( n, kl, ku, nrhs, ab, ldab, ipiv, b, ldb, info )
     130             :       
     131           0 :   end subroutine sgbsv_wrap
     132             :   
     133             :   
     134             :   ! ==================== Band Solver Expert Wrappers ====================
     135             : 
     136             :   ! Double precision wrapper
     137           0 :   subroutine dgbsvx_wrap( fact, trans, n, kl, ku, nrhs, ab, ldab, afb, &
     138           0 :                           ldafb, ipiv, equed, r, c, b, ldb, x, ldx, &
     139             :                           rcond, ferr, berr, work, iwork, info )
     140             :     implicit none
     141             :     
     142             :     external :: dgbsvx
     143             :  
     144             :     character          equed, fact, trans
     145             :     integer            info, kl, ku, ldab, ldafb, ldb, ldx, n, nrhs
     146             :     double precision   rcond
     147             :     integer            ipiv( * ), iwork( * )
     148             :     double precision   ab( ldab, * ), afb( ldafb, * ), b( ldb, * ), &
     149             :                        berr( * ), c( * ), ferr( * ), r( * ), &
     150             :                        work( * ), x( ldx, * )
     151             :                        
     152             :     call dgbsvx( fact, trans, n, kl, ku, nrhs, ab, ldab, afb, &
     153             :                  ldafb, ipiv, equed, r, c, b, ldb, x, ldx, &
     154           0 :                  rcond, ferr, berr, work, iwork, info )
     155             :                  
     156           0 :   end subroutine dgbsvx_wrap
     157             :   
     158             :   ! Single precision wrapper
     159           0 :   subroutine sgbsvx_wrap( fact, trans, n, kl, ku, nrhs, ab, ldab, afb, &
     160           0 :                           ldafb, ipiv, equed, r, c, b, ldb, x, ldx, &
     161             :                           rcond, ferr, berr, work, iwork, info )                  
     162             :     implicit none
     163             :     
     164             :     external :: sgbsvx
     165             :  
     166             :     character   equed, fact, trans
     167             :     integer     info, kl, ku, ldab, ldafb, ldb, ldx, n, nrhs
     168             :     real        rcond
     169             :     integer     ipiv( * ), iwork( * )
     170             :     real        ab( ldab, * ), afb( ldafb, * ), b( ldb, * ), &
     171             :                 berr( * ), c( * ), ferr( * ), r( * ), &
     172             :                 work( * ), x( ldx, * )
     173             :                        
     174             :     call sgbsvx( fact, trans, n, kl, ku, nrhs, ab, ldab, afb, &
     175             :                  ldafb, ipiv, equed, r, c, b, ldb, x, ldx, &
     176           0 :                  rcond, ferr, berr, work, iwork, info )
     177             :                  
     178           0 :   end subroutine sgbsvx_wrap
     179             :   
     180             :   
     181             :   ! ==================== Tridiagonal Solver Wrappers ====================
     182             : 
     183             :   ! Double precision wrapper
     184   707664624 :   subroutine dgtsv_wrap( n, nrhs, dl, d, du, b, ldb, info )
     185             :       
     186             :     implicit none
     187             :     
     188             :     external :: dgtsv
     189             :       
     190             :     integer            info, ldb, n, nrhs
     191             :     double precision   b( ldb, * ), d( * ), dl( * ), du( * )
     192             :     
     193   707664624 :     call dgtsv( n, nrhs, dl, d, du, b, ldb, info )
     194             :     
     195   707664624 :   end subroutine dgtsv_wrap
     196             :     
     197             :   ! Single precision wrapper
     198           0 :   subroutine sgtsv_wrap( n, nrhs, dl, d, du, b, ldb, info )
     199             :       
     200             :     implicit none
     201             :     
     202             :     external :: sgtsv
     203             :       
     204             :     integer     info, ldb, n, nrhs
     205             :     real        b( ldb, * ), d( * ), dl( * ), du( * )
     206             :     
     207           0 :     call sgtsv( n, nrhs, dl, d, du, b, ldb, info )
     208             :     
     209           0 :   end subroutine sgtsv_wrap
     210             :   
     211             :   
     212             :   ! ==================== Tridiagonal Solver Expert Wrappers ====================
     213             : 
     214             :   ! Double precision wrapper
     215           0 :   subroutine dgtsvx_wrap( fact, trans, n, nrhs, dl, d, du, dlf, df, duf, &
     216           0 :                           du2, ipiv, b, ldb, x, ldx, rcond, ferr, berr, &
     217             :                           work, iwork, info )
     218             :     implicit none
     219             :     
     220             :     external :: dgtsvx
     221             :   
     222             :     character          fact, trans
     223             :     integer            info, ldb, ldx, n, nrhs
     224             :     double precision   rcond
     225             :     integer            ipiv( * ), iwork( * )
     226             :     double precision   b( ldb, * ), berr( * ), d( * ), df( * ), &
     227             :                        dl( * ), dlf( * ), du( * ), du2( * ), duf( * ), &
     228             :                        ferr( * ), work( * ), x( ldx, * )
     229             :                        
     230             :     call dgtsvx( fact, trans, n, nrhs, dl, d, du, dlf, df, duf, &
     231             :                  du2, ipiv, b, ldb, x, ldx, rcond, ferr, berr, &
     232           0 :                  work, iwork, info )
     233             :                  
     234           0 :   end subroutine dgtsvx_wrap
     235             :   
     236             :   ! Single precision wrapper
     237           0 :   subroutine sgtsvx_wrap( fact, trans, n, nrhs, dl, d, du, dlf, df, duf, &
     238           0 :                           du2, ipiv, b, ldb, x, ldx, rcond, ferr, berr, &
     239             :                           work, iwork, info )
     240             :     implicit none
     241             :     
     242             :     external :: sgtsvx
     243             :  
     244             :     character   fact, trans
     245             :     integer     info, ldb, ldx, n, nrhs
     246             :     real        rcond
     247             :     integer     ipiv( * ), iwork( * )
     248             :     real        b( ldb, * ), berr( * ), d( * ), df( * ), &
     249             :                 dl( * ), dlf( * ), du( * ), du2( * ), duf( * ), &
     250             :                 ferr( * ), work( * ), x( ldx, * )
     251             :                        
     252             :     call sgtsvx( fact, trans, n, nrhs, dl, d, du, dlf, df, duf, &
     253             :                  du2, ipiv, b, ldb, x, ldx, rcond, ferr, berr, &
     254           0 :                  work, iwork, info )
     255             :                  
     256           0 :   end subroutine sgtsvx_wrap
     257             :   
     258             :   
     259             :   ! ==================== NaN Check Wrappers ====================
     260             : 
     261             :   ! Double precision wrapper
     262             :   !-----------------------------------------------------------------------
     263   707664624 :   logical function disnan_wrap( ndim, nrhs, variable )
     264             : 
     265             :   ! Description:
     266             :   !   Check for NaN values in a variable using the LAPACK subroutines
     267             : 
     268             :   ! References:
     269             :   !   <http://www.netlib.org/lapack/single/sisnan.f>
     270             :   !   <http://www.netlib.org/lapack/double/disnan.f>
     271             :   !-----------------------------------------------------------------------
     272             :   
     273             :     implicit none
     274             :       
     275             : #ifdef NO_LAPACK_ISNAN /* Used for older LAPACK libraries that don't have sisnan/disnan */
     276             : 
     277             :     integer, intent(in) :: &
     278             :       ndim, & ! Size of variable
     279             :       nrhs    ! Number of right hand sides
     280             : 
     281             :     double precision, dimension(ndim,nrhs), intent(in) :: &
     282             :       variable ! Variable to check
     283             : 
     284 >21837*10^7 :       disnan_wrap = any( variable(:,1:nrhs) /= variable(:,1:nrhs) )
     285             : #else
     286             :   
     287             :     logical, external :: &
     288             :       disnan    ! Procedure
     289             : 
     290             :     integer, intent(in) :: &
     291             :       ndim, & ! Size of variable
     292             :       nrhs    ! Number of right hand sides
     293             : 
     294             :     double precision, dimension(ndim,nrhs), intent(in) :: &
     295             :       variable ! Variable to check
     296             : 
     297             :     integer :: k, j
     298             : 
     299             :     ! ---- Begin Code ----
     300             : 
     301             :     disnan_wrap = .false.
     302             : 
     303             :     do k = 1, ndim
     304             :       do j = 1, nrhs
     305             :             
     306             :         ! Lapack NaN check function, sisnan for single precision or disnan for double precision
     307             :         disnan_wrap = disnan( variable(k,j) )
     308             :           
     309             :         if ( disnan_wrap ) exit
     310             :       end do
     311             :       if ( disnan_wrap ) exit
     312             :     end do
     313             :       
     314             : #endif /* NO_LAPACK_ISNAN */
     315             : 
     316             :     return
     317             :   end function disnan_wrap
     318             :   
     319             :   ! Single precision wrapper
     320             :   !-----------------------------------------------------------------------
     321           0 :   logical function sisnan_wrap( ndim, nrhs, variable )
     322             : 
     323             :   ! Description:
     324             :   !   Check for NaN values in a variable using the LAPACK subroutines
     325             : 
     326             :   ! References:
     327             :   !   <http://www.netlib.org/lapack/single/sisnan.f>
     328             :   !   <http://www.netlib.org/lapack/double/disnan.f>
     329             :   !-----------------------------------------------------------------------
     330             : 
     331             :     implicit none
     332             :       
     333             : #ifdef NO_LAPACK_ISNAN /* Used for older LAPACK libraries that don't have sisnan/disnan */
     334             : 
     335             :     integer, intent(in) :: &
     336             :       ndim, & ! Size of variable
     337             :       nrhs    ! Number of right hand sides
     338             : 
     339             :     real, dimension(ndim,nrhs), intent(in) :: &
     340             :       variable ! Variable to check
     341             : 
     342           0 :       sisnan_wrap = any( variable(:,1:nrhs) /= variable(:,1:nrhs) )
     343             : #else
     344             :   
     345             :     logical, external :: &
     346             :       sisnan    ! Procedure
     347             : 
     348             :     integer, intent(in) :: &
     349             :       ndim, & ! Size of variable
     350             :       nrhs    ! Number of right hand sides
     351             : 
     352             :     real, dimension(ndim,nrhs), intent(in) :: &
     353             :       variable ! Variable to check
     354             : 
     355             :     integer :: k, j
     356             : 
     357             :     ! ---- Begin Code ----
     358             : 
     359             :     sisnan_wrap = .false.
     360             : 
     361             :     do k = 1, ndim
     362             :       do j = 1, nrhs
     363             :             
     364             :         ! Lapack NaN check function, sisnan for single precision or disnan for double precision
     365             :         sisnan_wrap = sisnan( variable(k,j) )
     366             :           
     367             :         if ( sisnan_wrap ) exit
     368             :       end do
     369             :       if ( sisnan_wrap ) exit
     370             :     end do
     371             :       
     372             : #endif /* NO_LAPACK_ISNAN */
     373             : 
     374             :     return
     375             :   end function sisnan_wrap
     376             :  
     377             :  
     378             :   ! ==================== Cholesky Factorization Wrappers ====================
     379             : 
     380             :   ! Double precision wrapper
     381           0 :   subroutine dpotrf_wrap( uplo, n, a, lda, info )
     382             :       
     383             :     implicit none
     384             :     
     385             :     external :: dpotrf
     386             :       
     387             :     character          uplo
     388             :     integer            info, lda, n
     389             :     double precision   a( lda, * )
     390             : 
     391           0 :     call dpotrf( uplo, n, a, lda, info )
     392             :     
     393           0 :   end subroutine dpotrf_wrap
     394             :   
     395             :   ! Single precision wrapper
     396           0 :   subroutine spotrf_wrap( uplo, n, a, lda, info )
     397             :       
     398             :     implicit none
     399             :     
     400             :     external :: spotrf
     401             :       
     402             :     character   uplo
     403             :     integer     info, lda, n
     404             :     real        a( lda, * )
     405             : 
     406           0 :     call spotrf( uplo, n, a, lda, info )
     407             :     
     408           0 :   end subroutine spotrf_wrap
     409             :   
     410             :   
     411             :   ! ==================== Equilibrium Scaling Calculation Wrappers ====================
     412             : 
     413             :   ! Double precision wrapper
     414           0 :   subroutine dpoequ_wrap( n, a, lda, s, scond, amax, info )
     415             :       
     416             :     implicit none
     417             :     
     418             :     external :: dpoequ
     419             :       
     420             :     integer   info, lda, n
     421             :     double precision      amax, scond
     422             :     double precision      a( lda, * ), s( * )
     423             : 
     424           0 :     call dpoequ( n, a, lda, s, scond, amax, info )
     425             :   
     426           0 :   end subroutine dpoequ_wrap
     427             :   
     428             :   ! Single precision wrapper
     429           0 :   subroutine spoequ_wrap( n, a, lda, s, scond, amax, info )
     430             :       
     431             :     implicit none
     432             :     
     433             :     external :: spoequ
     434             :       
     435             :     integer   info, lda, n
     436             :     real      amax, scond
     437             :     real      a( lda, * ), s( * )
     438             : 
     439           0 :     call spoequ( n, a, lda, s, scond, amax, info )
     440             :   
     441           0 :   end subroutine spoequ_wrap
     442             :   
     443             :   
     444             :   ! ==================== Matrix Equilibriator Wrappers ====================
     445             : 
     446             :   ! Double precision wrapper
     447           0 :   subroutine dlaqsy_wrap( uplo, n, a, lda, s, scond, amax, equed )
     448             :       
     449             :     implicit none
     450             :     
     451             :     external :: dlaqsy
     452             :       
     453             :     character          equed, uplo
     454             :     integer            lda, n
     455             :     double precision   amax, scond
     456             :     double precision   a( lda, * ), s( * )
     457             :   
     458           0 :     call dlaqsy( uplo, n, a, lda, s, scond, amax, equed )
     459             :   
     460           0 :   end subroutine dlaqsy_wrap
     461             :     
     462             :   ! Single precision wrapper
     463           0 :   subroutine slaqsy_wrap( uplo, n, a, lda, s, scond, amax, equed )
     464             :       
     465             :     implicit none
     466             :     
     467             :     external :: slaqsy
     468             :       
     469             :     character   equed, uplo
     470             :     integer     lda, n
     471             :     real        amax, scond
     472             :     real        a( lda, * ), s( * )
     473             :   
     474           0 :     call slaqsy( uplo, n, a, lda, s, scond, amax, equed )
     475             :     
     476           0 :   end subroutine slaqsy_wrap
     477             :   
     478             :   
     479             :   ! ==================== Eigenvalue/vector Calculation Wrappers ====================
     480             : 
     481             :   ! Double precision wrapper
     482           0 :   subroutine dsyev_wrap( jobz, uplo, n, a, lda, w, work, lwork, info )
     483             :       
     484             :     implicit none
     485             :     
     486             :     external :: dsyev
     487             :       
     488             :     character          jobz, uplo
     489             :     integer            info, lda, lwork, n
     490             :     double precision   a( lda, * ), w( * ), work( * )
     491             :     
     492           0 :     call dsyev( jobz, uplo, n, a, lda, w, work, lwork, info )
     493             :       
     494           0 :   end subroutine dsyev_wrap
     495             :   
     496             :   ! Single precision wrapper
     497           0 :   subroutine ssyev_wrap( jobz, uplo, n, a, lda, w, work, lwork, info )
     498             :      
     499             :     implicit none
     500             :     
     501             :     external :: ssyev
     502             :       
     503             :     character   jobz, uplo
     504             :     integer     info, lda, lwork, n
     505             :     real        a( lda, * ), w( * ), work( * )
     506             :     
     507           0 :     call ssyev( jobz, uplo, n, a, lda, w, work, lwork, info )
     508             :       
     509           0 :   end subroutine ssyev_wrap
     510             :   
     511             :   
     512             :   ! ==================== Matrix Operations Wrappers ====================
     513             :   
     514             :   ! Double precision wrapper
     515           0 :   subroutine dtrmv_wrap( uplo, trans, diag, n, a, lda, x, incx)
     516             :      
     517             :     implicit none
     518             :     
     519             :     external :: dtrmv
     520             :       
     521             :     integer             incx,lda,n
     522             :     character           diag,trans,uplo             
     523             :     double precision    a(lda,*),x(*)
     524             : 
     525           0 :     call dtrmv( uplo, trans, diag, n, a, lda, x, incx)
     526             :     
     527           0 :   end subroutine dtrmv_wrap
     528             :   
     529             :   ! Single precision wrapper
     530           0 :   subroutine strmv_wrap( uplo, trans, diag, n, a, lda, x, incx)
     531             :      
     532             :     implicit none
     533             :     
     534             :     external :: strmv
     535             :       
     536             :     integer     incx,lda,n
     537             :     character   diag,trans,uplo             
     538             :     real        a(lda,*),x(*)
     539             : 
     540           0 :     call strmv( uplo, trans, diag, n, a, lda, x, incx)
     541             :     
     542           0 :   end subroutine strmv_wrap
     543             :           
     544             : end module lapack_interfaces

Generated by: LCOV version 1.14