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 22:39:59 |  | 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   707206168 :   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   707206168 :     call dgtsv( n, nrhs, dl, d, du, b, ldb, info )
         194             :
         195   707206168 :   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   707206168 :   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 >33896*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