linalg: enable 80-bit extended precision for whole library, xdp (#839)

fortran-lang · Jul 4, 2024 · 4859081 · 4859081
2 parents 6d9d7fd + 12a6640
commit 4859081
Show file tree

Hide file tree

Showing 22 changed files with 25,332 additions and 23,752 deletions.
diff --git a/include/common.fypp b/include/common.fypp
@@ -60,6 +60,19 @@
 #:set CMPLX_INIT = CMPLX_INIT + ["w"]
 #:endif
 
+#! BLAS/LAPACK complex->real kind initial conversion
+#! Converts a BLAS/LAPACK complex kind initial to a real kind initial
+#!
+#! Args:
+#!     ci (character): Complex kind initial in ["c","z","y","w"]
+#!
+#! Returns:
+#!     Real kind initial in ["s","d","x","q"] or an empty string on invalid input
+#!
+#:def c2ri(cmplx)
+$:"s" if cmplx=="c" else "d" if cmplx=="z" else "x" if cmplx=="y" else "q" if cmplx=="w" else "ERROR"
+#:enddef
+
 #! Complex types to be considered during templating
 #:set CMPLX_TYPES = ["complex({})".format(k) for k in CMPLX_KINDS]
 

diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp
@@ -255,11 +255,9 @@ module stdlib_linalg
      !! or several (from a 2-d right-hand-side vector `b(:,:)`) systems.
      !! 
      !!@note The solution is based on LAPACK's generic LU decomposition based solvers `*GESV`.
-     !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).
      !!    
      #:for nd,ndsuf,nde in ALL_RHS
      #:for rk,rt,ri in RC_KINDS_TYPES
-     #:if rk!="xdp"
      module function stdlib_linalg_${ri}$_solve_${ndsuf}$(a,b,overwrite_a,err) result(x)
          !> Input matrix a[n,n]
          ${rt}$, intent(inout), target :: a(:,:)
@@ -280,7 +278,6 @@ module stdlib_linalg
          !> Result array/matrix x[n] or x[n,nrhs]
          ${rt}$, allocatable, target :: x${nd}$
      end function stdlib_linalg_${ri}$_pure_solve_${ndsuf}$
-     #:endif
      #:endfor
      #:endfor    
   end interface solve
@@ -306,11 +303,9 @@ module stdlib_linalg
      !! or several (from a 2-d right-hand-side vector `b(:,:)`) systems.
      !! 
      !!@note The solution is based on LAPACK's generic LU decomposition based solvers `*GESV`.
-     !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).
      !!        
      #:for nd,ndsuf,nde in ALL_RHS
      #:for rk,rt,ri in RC_KINDS_TYPES
-     #:if rk!="xdp"
      pure module subroutine stdlib_linalg_${ri}$_solve_lu_${ndsuf}$(a,b,x,pivot,overwrite_a,err)     
          !> Input matrix a[n,n]
          ${rt}$, intent(inout), target :: a(:,:)
@@ -325,7 +320,6 @@ module stdlib_linalg
          !> [optional] state return flag. On error if not requested, the code will stop
          type(linalg_state_type), optional, intent(out) :: err
      end subroutine stdlib_linalg_${ri}$_solve_lu_${ndsuf}$
-     #:endif
      #:endfor
      #:endfor    
   end interface solve_lu     
@@ -346,11 +340,9 @@ module stdlib_linalg
     !! Supported data types include `real` and `complex`.
     !! 
     !!@note The solution is based on LAPACK's singular value decomposition `*GELSD` methods.
-    !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).
     !! 
     #:for nd,ndsuf,nde in ALL_RHS
     #:for rk,rt,ri in RC_KINDS_TYPES
-    #:if rk!="xdp"
       module function stdlib_linalg_${ri}$_lstsq_${ndsuf}$(a,b,cond,overwrite_a,rank,err) result(x)
          !> Input matrix a[n,n]
          ${rt}$, intent(inout), target :: a(:,:)
@@ -367,7 +359,6 @@ module stdlib_linalg
          !> Result array/matrix x[n] or x[n,nrhs]
          ${rt}$, allocatable, target :: x${nd}$
       end function stdlib_linalg_${ri}$_lstsq_${ndsuf}$
-    #:endif
     #:endfor
     #:endfor
   end interface lstsq
@@ -389,11 +380,9 @@ module stdlib_linalg
     !! are provided, no internal memory allocations take place when using this interface.
     !! 
     !!@note The solution is based on LAPACK's singular value decomposition `*GELSD` methods.
-    !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).
     !! 
     #:for nd,ndsuf,nde in ALL_RHS
     #:for rk,rt,ri in RC_KINDS_TYPES
-    #:if rk!="xdp"
       module subroutine stdlib_linalg_${ri}$_solve_lstsq_${ndsuf}$(a,b,x,real_storage,int_storage,&
                         #{if rt.startswith('c')}#cmpl_storage,#{endif}#cond,singvals,overwrite_a,rank,err) 
          !> Input matrix a[n,n]
@@ -421,7 +410,6 @@ module stdlib_linalg
          !> [optional] state return flag. On error if not requested, the code will stop
          type(linalg_state_type), optional, intent(out) :: err
       end subroutine stdlib_linalg_${ri}$_solve_lstsq_${ndsuf}$
-    #:endif
     #:endfor
     #:endfor
   end interface solve_lstsq
@@ -442,7 +430,6 @@ module stdlib_linalg
     !! 
     #:for nd,ndsuf,nde in ALL_RHS
     #:for rk,rt,ri in RC_KINDS_TYPES
-    #:if rk!="xdp" 
       pure module subroutine stdlib_linalg_${ri}$_lstsq_space_${ndsuf}$(a,b,lrwork,liwork#{if rt.startswith('c')}#,lcwork#{endif}#)
          !> Input matrix a[m,n]
          ${rt}$, intent(in), target :: a(:,:)
@@ -451,7 +438,6 @@ module stdlib_linalg
          !> Size of the working space arrays                
          integer(ilp), intent(out) :: lrwork,liwork#{if rt.startswith('c')}#,lcwork#{endif}#         
       end subroutine stdlib_linalg_${ri}$_lstsq_space_${ndsuf}$
-    #:endif
     #:endfor
     #:endfor
   end interface lstsq_space
@@ -781,7 +767,6 @@ module stdlib_linalg
     !! It is possible to use partial storage [m,k] and [k,n], `k=min(m,n)`, choosing `full_matrices=.false.`.
     !! 
     !!@note The solution is based on LAPACK's singular value decomposition `*GESDD` methods.
-    !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).    
     !! 
     !!### Example
     !!
@@ -794,7 +779,6 @@ module stdlib_linalg
     !!```
     !!         
     #:for rk,rt,ri in RC_KINDS_TYPES
-    #:if rk!="xdp"
     module subroutine stdlib_linalg_svd_${ri}$(a,s,u,vt,overwrite_a,full_matrices,err)
      !!### Summary
      !! Compute singular value decomposition of a matrix \( A = U \cdot S \cdot \V^T \)
@@ -830,7 +814,6 @@ module stdlib_linalg
          !> [optional] state return flag. On error if not requested, the code will stop
          type(linalg_state_type), optional, intent(out) :: err
       end subroutine stdlib_linalg_svd_${ri}$
-    #:endif
     #:endfor     
   end interface svd
 
@@ -853,7 +836,6 @@ module stdlib_linalg
     !! singular values, with size [min(m,n)]. 
     !! 
     !!@note The solution is based on LAPACK's singular value decomposition `*GESDD` methods.
-    !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).    
     !! 
     !!### Example
     !!
@@ -866,7 +848,6 @@ module stdlib_linalg
     !!```
     !!   
     #:for rk,rt,ri in RC_KINDS_TYPES
-    #:if rk!="xdp"
     module function stdlib_linalg_svdvals_${ri}$(a,err) result(s)
      !!### Summary
      !! Compute singular values \(S \) from the singular-value decomposition of a matrix \( A = U \cdot S \cdot \V^T \).
@@ -890,7 +871,6 @@ module stdlib_linalg
          !> Array of singular values
          real(${rk}$), allocatable :: s(:)
       end function stdlib_linalg_svdvals_${ri}$
-    #:endif
     #:endfor             
   end interface svdvals