From 61bf91c330951d0ee87658d93cc3fa5e448f8874 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 27 Nov 2024 17:29:56 +0100 Subject: [PATCH] write specs --- doc/specs/stdlib_linalg.md | 136 ++++++++++++++++++++++++++++++++++++ src/stdlib_linalg.fypp | 6 +- src/stdlib_linalg_pinv.fypp | 10 +-- 3 files changed, 144 insertions(+), 8 deletions(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 2f03f0fad..d507c006e 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -1459,6 +1459,142 @@ If `err` is not present, exceptions trigger an `error stop`. {!example/linalg/example_inverse_function.f90!} ``` +## `pinv` - Moore-Penrose pseudo-inverse of a matrix. + +### Status + +Experimental + +### Description + +This function computes the Moore-Penrose pseudo-inverse of a `real` or `complex` matrix. +The pseudo-inverse, \( A^{+} \), generalizes the matrix inverse and satisfies the conditions: +- \( A \cdot A^{+} \cdot A = A \) +- \( A^{+} \cdot A \cdot A^{+} = A^{+} \) +- \( (A \cdot A^{+})^T = A \cdot A^{+} \) +- \( (A^{+} \cdot A)^T = A^{+} \cdot A \) + +The computation is based on singular value decomposition (SVD). Singular values below a relative +tolerance threshold \( \text{rtol} \cdot \sigma_{\max} \), where \( \sigma_{\max} \) is the largest +singular value, are treated as zero. + +### Syntax + +`b ` [[stdlib_linalg(module):pinv(interface)]] `(a, [, rtol, err])` + +### Arguments + +`a`: Shall be a rank-2, `real` or `complex` array of shape `[m, n]` containing the coefficient matrix. +It is an `intent(in)` argument. + +`rtol` (optional): Shall be a scalar `real` value specifying the relative tolerance for singular value cutoff. +If `rtol` is not provided, the default relative tolerance is \( \text{rtol} = \text{max}(m, n) \cdot \epsilon \), +where \( \epsilon \) is the machine precision for the element type of `a`. It is an `intent(in)` argument. + +`err` (optional): Shall be a `type(linalg_state_type)` value. It is an `intent(out)` argument. + +### Return value + +Returns an array value of the same type, kind, and rank as `a` with shape `[n, m]`, that contains the pseudo-inverse matrix \( A^{+} \). + +Raises `LINALG_ERROR` if the underlying SVD did not converge. +Raises `LINALG_VALUE_ERROR` if `pinva` and `a` have degenerate or incompatible sizes. +If `err` is not present, exceptions trigger an `error stop`. + +### Example + +```fortran +{!example/linalg/example_pseudoinverse.f90!} +``` + +## `pseudoinvert` - Moore-Penrose pseudo-inverse of a matrix + +### Status + +Experimental + +### Description + +This subroutine computes the Moore-Penrose pseudo-inverse of a `real` or `complex` matrix. +The pseudo-inverse \( A^{+} \) is a generalization of the matrix inverse and satisfies the following properties: +- \( A \cdot A^{+} \cdot A = A \) +- \( A^{+} \cdot A \cdot A^{+} = A^{+} \) +- \( (A \cdot A^{+})^T = A \cdot A^{+} \) +- \( (A^{+} \cdot A)^T = A^{+} \cdot A \) + +The computation is based on singular value decomposition (SVD). Singular values below a relative +tolerance threshold \( \text{rtol} \cdot \sigma_{\max} \), where \( \sigma_{\max} \) is the largest +singular value, are treated as zero. + +On return, matrix `pinva` `[n, m]` will store the pseudo-inverse of `a` `[m, n]`. + +### Syntax + +`call ` [[stdlib_linalg(module):pseudoinvert(interface)]] `(a, pinva [, rtol] [, err])` + +### Arguments + +`a`: Shall be a rank-2, `real` or `complex` array containing the coefficient matrix. +It is an `intent(in)` argument. + +`pinva`: Shall be a rank-2 array of the same kind as `a`, and size equal to that of `transpose(a)`. +On output, it contains the Moore-Penrose pseudo-inverse of `a`. + +`rtol` (optional): Shall be a scalar `real` value specifying the relative tolerance for singular value cutoff. +If not provided, the default threshold is \( \text{max}(m, n) \cdot \epsilon \), where \( \epsilon \) is the +machine precision for the element type of `a`. + +`err` (optional): Shall be a `type(linalg_state_type)` value. It is an `intent(out)` argument. + +### Return value + +Computes the Moore-Penrose pseudo-inverse of the matrix \( A \), \( A^{+} \), and returns it in matrix `pinva`. + +Raises `LINALG_ERROR` if the underlying SVD did not converge. +Raises `LINALG_VALUE_ERROR` if `pinva` and `a` have degenerate or incompatible sizes. +If `err` is not present, exceptions trigger an `error stop`. + +### Example + +```fortran +{!example/linalg/example_pseudoinverse.f90!} +``` + +## `.pinv.` - Moore-Penrose Pseudo-Inverse operator + +### Status + +Experimental + +### Description + +This operator returns the Moore-Penrose pseudo-inverse of a `real` or `complex` matrix \( A \). +The pseudo-inverse \( A^{+} \) is computed using Singular Value Decomposition (SVD), and singular values +below a given threshold are treated as zero. + +This interface is equivalent to the function [[stdlib_linalg(module):pinv(interface)]]. + +### Syntax + +`b = ` [[stdlib_linalg(module):operator(.pinv.)(interface)]] `a` + +### Arguments + +`a`: Shall be a rank-2 array of any `real` or `complex` kinds, with arbitrary dimensions \( m \times n \). It is an `intent(in)` argument. + +### Return value + +Returns a rank-2 array with the same type, kind, and rank as `a`, that contains the Moore-Penrose pseudo-inverse of `a`. + +If an exception occurs, or if the input matrix is degenerate (e.g., rank-deficient), the returned matrix will contain `NaN`s. +For more detailed error handling, it is recommended to use the `subroutine` or `function` interfaces. + +### Example + +```fortran +{!example/linalg/example_pseudoinverse.f90!} +``` + ## `get_norm` - Computes the vector norm of a generic-rank array. ### Status diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 33b973b09..bb3be85b1 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -817,7 +817,7 @@ module stdlib_linalg !! version: experimental !! !! Pseudo-inverse of a matrix - !! ([Specification](../page/specs/stdlib_linalg.html#pinv-pseudo-inverse-of-a-matrix)) + !! ([Specification](../page/specs/stdlib_linalg.html#pinv-moore-penrose-pseudo-inverse-of-a-matrix)) !! !!### Summary !! This interface provides methods for computing the Moore-Penrose pseudo-inverse of a matrix. @@ -862,7 +862,7 @@ module stdlib_linalg !! version: experimental !! !! Computation of the Moore-Penrose pseudo-inverse - !! ([Specification](../page/specs/stdlib_linalg.html#pseudoinvert-computation-of-a-matrix-pseudo-inverse)) + !! ([Specification](../page/specs/stdlib_linalg.html#pseudoinvert-moore-penrose-pseudo-inverse-of-a-matrix)) !! !!### Summary !! This interface provides methods for computing the Moore-Penrose pseudo-inverse of a rectangular @@ -894,7 +894,7 @@ module stdlib_linalg #:for rk,rt,ri in RC_KINDS_TYPES module subroutine stdlib_linalg_pseudoinvert_${ri}$(a,pinva,rtol,err) !> Input matrix a[m,n] - ${rt}$, intent(in) :: a(:,:) + ${rt}$, intent(inout) :: a(:,:) !> Output pseudo-inverse matrix [n,m] ${rt}$, intent(out) :: pinva(:,:) !> [optional] Relative tolerance for singular value cutoff diff --git a/src/stdlib_linalg_pinv.fypp b/src/stdlib_linalg_pinv.fypp index ce201b762..570ce3f72 100644 --- a/src/stdlib_linalg_pinv.fypp +++ b/src/stdlib_linalg_pinv.fypp @@ -7,7 +7,7 @@ submodule(stdlib_linalg) stdlib_linalg_pseudoinverse use stdlib_linalg_lapack use stdlib_linalg_state use stdlib_linalg, only: svd - use iso_fortran_env,only:real32,real64,real128,int8,int16,int32,int64,stderr => error_unit + use ieee_arithmetic, only: ieee_value, ieee_quiet_nan implicit none(type,external) character(*), parameter :: this = 'pseudo-inverse' @@ -21,8 +21,8 @@ submodule(stdlib_linalg) stdlib_linalg_pseudoinverse !> Input matrix a[m,n] ${rt}$, intent(inout) :: a(:,:) !> Output pseudo-inverse matrix - ${rt}$, intent(inout) :: pinva(:,:) - !> [optional] .... + ${rt}$, intent(out) :: pinva(:,:) + !> [optional] Relative tolerance for singular value cutoff real(${rk}$), optional, intent(in) :: rtol !> [optional] state return flag. On error if not requested, the code will stop type(linalg_state_type), optional, intent(out) :: err @@ -89,7 +89,7 @@ submodule(stdlib_linalg) stdlib_linalg_pseudoinverse module function stdlib_linalg_pseudoinverse_${ri}$(a,rtol,err) result(pinva) !> Input matrix a[m,n] ${rt}$, intent(in), target :: a(:,:) - !> [optional] .... + !> [optional] Relative tolerance for singular value cutoff real(${rk}$), optional, intent(in) :: rtol !> [optional] state return flag. On error if not requested, the code will stop type(linalg_state_type), optional, intent(out) :: err @@ -108,7 +108,7 @@ submodule(stdlib_linalg) stdlib_linalg_pseudoinverse module function stdlib_linalg_pinv_${ri}$_operator(a) result(pinva) !> Input matrix a[m,n] ${rt}$, intent(in), target :: a(:,:) - !> Result matrix + !> Result pseudo-inverse matrix ${rt}$ :: pinva(size(a,2,kind=ilp),size(a,1,kind=ilp)) type(linalg_state_type) :: err