(new) Added pocon routine for testing new approach

14NGiestas · Nov 28, 2024 · 88d86fb · 88d86fb
1 parent bc60c03
commit 88d86fb
Show file tree

Hide file tree

Showing 4 changed files with 108 additions and 28 deletions.
diff --git a/README.md b/README.md
@@ -100,6 +100,7 @@ Please note that this project is experimental, is missing a test suite and error
 ### BLAS
 #### Level 1
 Most of BLAS level 1 routines can be replaced by intrinsincs and other features in modern fortran.
+<details>
 
 | done? | name   | description                                             | modern alternative |
 | ----- | ------ | ------------------------------------------------------- | ------------------ |
@@ -117,9 +118,22 @@ Most of BLAS level 1 routines can be replaced by intrinsincs and other features
 | :+1:  | rotmg  | Generate modified Givens plane rotation of points       | |
 |       | scal   | Vector-scalar product                                   | `a*x + b` |
 | :+1:  | swap   | Vector-vector swap                                      | |
+</details>
+
+#### Level 1 - Utils / Extensions
+<details>
+
+| done? | name  | description                                              |  modern alternatives | obs |
+| ----- | ----- | -------------------------------------------------------- | ------------------- | --- |
+| :+1:  | iamax | Index of the maximum absolute value element of a vector  | [maxval](https://gcc.gnu.org/onlinedocs/gfortran/MAXVAL.html), [maxloc](https://gcc.gnu.org/onlinedocs/gfortran/MAXLOC.html) | |
+| :+1:  | iamin | Index of the minimum absolute value element of a vector  | [minval](https://gcc.gnu.org/onlinedocs/gfortran/MINVAL.html), [minloc](https://gcc.gnu.org/onlinedocs/gfortran/MINLOC.html) | |
+|  :(   | lamch | Determines precision machine parameters.                 | [huge](https://gcc.gnu.org/onlinedocs/gfortran/intrinsic-procedures/huge.html), [tiny](https://gcc.gnu.org/onlinedocs/gfortran/intrinsic-procedures/tiny.html), [epsilon](https://gcc.gnu.org/onlinedocs/gfortran/intrinsic-procedures/epsilon.html) | Obs: had to add a parameter so fortran can distinguish between the single and double precision with the same interface. For values of cmach see: [lamch](https://www.netlib.org/lapack//explore-html/d4/d86/group__lamch.html)|
+</details>
 
 #### Level 2
 
+<details>
+
 | done? | name | description                                                              |
 | ----- | ---- | ------------------------------------------------------------------------ |
 | :+1:  | gbmv | Matrix-vector product using a general band matrix                        |
@@ -147,9 +161,12 @@ Most of BLAS level 1 routines can be replaced by intrinsincs and other features
 | :+1:  | tpsv | Solution of a linear system of equations with a triangular packed matrix |
 | :+1:  | trmv | Matrix-vector product using a triangular matrix                          |
 | :+1:  | trsv | Solution of a linear system of equations with a triangular matrix        |
+</details>
 
 #### Level 3
 
+<details>
+
 | done? | name  | description                                                                                            |
 | ----- | ----- | ------------------------------------------------------------------------------------------------------ |
 | :+1:  | gemm  | Computes a matrix-matrix product with general matrices.                                                |
@@ -162,34 +179,52 @@ Most of BLAS level 1 routines can be replaced by intrinsincs and other features
 | :+1:  | trmm  | Computes a matrix-matrix product where one input matrix is triangular and one input matrix is general. |
 | :+1:  | trsm  | Solves a triangular matrix equation (forward or backward solve).                                       |
 
-#### Utils / Extensions
-#### Level 1
-Here are some extensions that may be useful.
-Again, BLAS level 1 routines can be replaced by intrinsincs and other features in modern fortran.
-
-| done? | name  | description                                              |  modern alternatives | obs |
-| ----- | ----- | -------------------------------------------------------- | ------------------- | --- |
-| :+1:  | iamax | Index of the maximum absolute value element of a vector  | [maxval](https://gcc.gnu.org/onlinedocs/gfortran/MAXVAL.html), [maxloc](https://gcc.gnu.org/onlinedocs/gfortran/MAXLOC.html) | |
-| :+1:  | iamin | Index of the minimum absolute value element of a vector  | [minval](https://gcc.gnu.org/onlinedocs/gfortran/MINVAL.html), [minloc](https://gcc.gnu.org/onlinedocs/gfortran/MINLOC.html) | |
-|  :(   | lamch | Determines precision machine parameters.                 | [huge](https://gcc.gnu.org/onlinedocs/gfortran/intrinsic-procedures/huge.html), [tiny](https://gcc.gnu.org/onlinedocs/gfortran/intrinsic-procedures/tiny.html), [epsilon](https://gcc.gnu.org/onlinedocs/gfortran/intrinsic-procedures/epsilon.html) | Obs: had to add a parameter so fortran can distinguish between the single and double precision with the same interface. For values of cmach see: [lamch](https://www.netlib.org/lapack//explore-html/d4/d86/group__lamch.html)|
-
-### LAPACK
-#### Linear Equation Routines
-
-| done? | name  | description                                                                                                                                      |
-| ----- | ----- | ------------------------------------------------------------------------------------------------------------------------------------------------ |
-| :+1:  | geqrf | Computes the QR factorization of a general m-by-n matrix.                                                                                        |
-| :+1:  | gerqf | Computes the RQ factorization of a general m-by-n matrix.                                                                                        |
-| :+1:  | getrf | Computes the LU factorization of a general m-by-n matrix.                                                                                        |
-| :+1:  | getri | Computes the inverse of an LU-factored general matrix.                                                                                           |
-| :+1:  | getrs | Solves a system of linear equations with an LU-factored square coefficient matrix, with multiple right-hand sides.                               |
-| :+1:  | hetrf | Computes the Bunch-Kaufman factorization of a complex Hermitian matrix.                                                                          |
-|       | orgqr | Generates the real orthogonal matrix Q of the QR factorization formed by geqrf.                                                                  |
-|       | ormqr | Multiplies a real matrix by the orthogonal matrix Q of the QR factorization formed by geqrf.                                                     |
-|       | ormrq | Multiplies a real matrix by the orthogonal matrix Q of the RQ factorization formed by gerqf.                                                     |
-| :+1:  | potrf | Computes the Cholesky factorization of a symmetric (Hermitian) positive-definite matrix.                                                         |
-| :+1:  | potri | Computes the inverse of a Cholesky-factored symmetric (Hermitian) positive-definite matrix.                                                      |
+</details>
+
+### LAPACK :warning:
+
+- Lapack is really huge, so I'm going to focus on getting the improved f77 interfaces ready first.
+  Anything I end up using I'm going to implement.
+
+#### Linear solve, $AX = B$
+<details>
+<!-- ##### LU: General matrix, driver -->
+
+<!-- ##### LU: computational routines (factor, cond, etc.) -->
+
+<!-- ##### Cholesky: Hermitian/symmetric positive definite matrix, driver -->
+
+##### Cholesky: computational routines (factor, cond, etc.)
+| done?| name  | description               |
+| ---- | ----- | ------------------------- |
+| f77  | pocon | condition number estimate |
+
+<!-- ##### LDL: Hermitian/symmetric indefinite matrix, driver -->
+
+<!-- ##### LDL: computational routines (factor, cond, etc.) -->
+
+<!-- ##### Triangular computational routines (solve, cond, etc.) -->
+
+<!-- ##### Auxiliary routines -->
+</details>
+
+##### Orthogonal/unitary factors (QR, CS, etc.)
+<details>
+
+| done? | name  | description  |
+| ----- | ----- | ------------ |
+| :+1:  | geqrf | Computes the QR factorization of a general m-by-n matrix. |
+| :+1:  | gerqf | Computes the RQ factorization of a general m-by-n matrix. |
+| :+1:  | getrf | Computes the LU factorization of a general m-by-n matrix. |
+| :+1:  | getri | Computes the inverse of an LU-factored general matrix.    |
+| :+1:  | getrs | Solves a system of linear equations with an LU-factored square coefficient matrix, with multiple right-hand sides. |
+| :+1:  | hetrf | Computes the Bunch-Kaufman factorization of a complex Hermitian matrix. |
+| :+1:  | potrf | Computes the Cholesky factorization of a symmetric (Hermitian) positive-definite matrix.     |
+| :+1:  | potri | Computes the inverse of a Cholesky-factored symmetric (Hermitian) positive-definite matrix.  |
 | :+1:  | potrs | Solves a system of linear equations with a Cholesky-factored symmetric (Hermitian) positive-definite coefficient matrix, with multiple right-hand sides.  |
+|      | orgqr | Generates the real orthogonal matrix Q of the QR factorization formed by geqrf. |
+|      | ormqr | Multiplies a real matrix by the orthogonal matrix Q of the QR factorization formed by geqrf. |
+|      | ormrq | Multiplies a real matrix by the orthogonal matrix Q of the RQ factorization formed by gerqf. |
 |      | sytrf | Computes the Bunch-Kaufman factorization of a symmetric matrix.                        |
 |      | trtrs | Solves a system of linear equations with a triangular coefficient matrix, with multiple right-hand sides. |
 |      | ungqr | Generates the complex unitary matrix Q of the QR factorization formed by geqrf.  |
@@ -199,12 +234,12 @@ Again, BLAS level 1 routines can be replaced by intrinsincs and other features i
 #### Singular Value and Eigenvalue Problem Routines
 | done?| name  | description             |
 | ---- | ----- | ----------------------- |
-|      | gebrd | Reduces a general matrix to bidiagonal form.     |
 | :+1: | gesvd | Computes the singular value decomposition of a general rectangular matrix.  |
 | :+1: | heevd | Computes all eigenvalues and, optionally, all eigenvectors of a complex Hermitian matrix using divide and conquer algorithm. |
 | :+1: | hegvd | Computes all eigenvalues and, optionally, all eigenvectors of a complex generalized Hermitian definite eigenproblem using divide and conquer algorithm. |
 | f77  | heevr | Computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices. |
 | f77  | heevx | Computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices. |
+|      | gebrd | Reduces a general matrix to bidiagonal form.     |
 |      | hetrd | Reduces a complex Hermitian matrix to tridiagonal form. |
 |      | orgbr | Generates the real orthogonal matrix Q or PT determined by gebrd. |
 |      | orgtr | Generates the real orthogonal matrix Q determined by sytrd. |

diff --git a/common.fpp b/common.fpp
@@ -19,6 +19,13 @@
     'z':   'REAL64',  &
 }
 
+#:set TYPE_AND_KIND_TO_PREFIX = { &
+    'real(REAL32)':    's',       &
+    'real(REAL64)':    'd',       &
+    'complex(REAL32)': 'c',       &
+    'complex(REAL64)': 'z',       &
+}
+
 #! Defines a optional variable, creating local corresponding variable by default
 #:def optional(dtype, intent, *args)
 #:for variable in args

diff --git a/src/f77/lapack.fpp b/src/f77/lapack.fpp
@@ -13,6 +13,7 @@
 #:include "src/f77/lapack/heevd.fypp"
 #:include "src/f77/lapack/potrf_potri.fypp"
 #:include "src/f77/lapack/potrs.fypp"
+#:include "src/f77/lapack/pocon.fypp"
 #:endmute
 module f77_lapack
 use iso_fortran_env
@@ -32,6 +33,7 @@ $:f77_interface('?gesvd',  DEFAULT_TYPES, gesvd)
 $:f77_interface('?potrf',  DEFAULT_TYPES, potrf_potri)
 $:f77_interface('?potri',  DEFAULT_TYPES, potrf_potri)
 $:f77_interface('?potrs',  DEFAULT_TYPES, potrs)
+$:f77_interface('?pocon',  DEFAULT_TYPES, pocon)
 
 ! Other Auxiliary Routines
 $:f77_interface('?lartg',  DEFAULT_TYPES, aux_lartg)

diff --git a/src/f77/lapack/pocon.fypp b/src/f77/lapack/pocon.fypp
@@ -0,0 +1,36 @@
+#:mute
+
+!subroutine ?pocon (
+!   character                      uplo,
+!   integer                           n,
+!   type(wp), dimension( lda, * )     a,
+!   integer                         lda,
+!   real(wp)                      anorm,
+!   real(wp)                      rcond,
+!   type(wp),    dimension( * )    work,
+!   real(wp),    dimension( * )   rwork,
+!   integer                        info 
+!)
+#:def pocon(NAME, TYPE, KIND)
+#:set TYPE_AND_KIND = TYPE.replace('wp',KIND)
+#:set PREFIX        = TYPE_AND_KIND_TO_PREFIX.get(TYPE_AND_KIND).upper()
+!> ${NAME}$ estimates the reciprocal of the condition number (in the
+!> 1-norm) of a ${TYPE_AND_KIND}$ Hermitian positive definite matrix using the
+!> Cholesky factorization A = U**H*U or A = L*L**H computed by ${PREFIX}$POTRF.
+!> An estimate is obtained for norm(inv(A)), and the reciprocal of the
+!> condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
+pure subroutine ${NAME}$(uplo, n, a, lda, anorm, rcond, work, rwork, info)
+    import :: ${KIND}$
+@:parameter(integer, wp=${KIND}$)
+@:args(character,     in,    uplo)
+@:args(integer,       in,    n, lda)
+@:args(${TYPE}$,      inout, a(lda,*))
+@:args(${REAL_TYPE}$, in,    anorm)
+@:args(${REAL_TYPE}$, out,   rcond)
+@:args(${TYPE}$,      inout, work(*))
+@:args(${REAL_TYPE}$, inout, rwork(*))
+@:args(integer,       out,   info)
+end subroutine
+#:enddef
+
+#:endmute