Skip to content

Commit 8a030ea

Browse files
perazzperazz
committed
document inverse
1 parent 26d2a1e commit 8a030ea

File tree

3 files changed

+186
-53
lines changed

3 files changed

+186
-53
lines changed

‎README.md

Lines changed: 92 additions & 37 deletions
Original file line numberDiff line numberDiff line change
@@ -100,73 +100,128 @@ The function returns a `real` scalar value representing the determinant of the i
100100
- The determinant of the matrix is computed using the LAPACK [getrf](@ref la_lapack::getrf) backend.
101101
- If `overwrite_a` is enabled, the input matrix `a` will be destroyed during the computation process.
102102

103-
## `inv(A)`
104-
**Type**: Function
105-
**Description**: Inverse of a scalar or square matrix.
106-
**Optional arguments**:
107-
- `err`: Return state handler.
108103

109-
## `invert(A)`
110-
**Type**: Subroutine
111-
**Description**: In-place inverse of a scalar or square matrix.
112-
**Optional arguments**:
113-
- `err`: Return state handler.
114104

115-
**Usage**: `call invert(A, err=err)` where `A` is replaced with $A^{-1}$.
105+
## [inv](@ref la_inverse::inv) - Inverse of a square matrix.
106+
107+
### Syntax
108+
109+
`inv_a = inv(a [, err])`
116110

117-
## `.inv.A`
118-
**Type**: Operator
119-
**Description**: Inverse of a scalar or square matrix.
111+
### Description
112+
113+
This function computes the inverse \f$ A^{-1} \f$ of a real or complex square matrix \f$ A \f$, provided that \f$ A \f$ is non-singular.
114+
The inverse of a matrix is defined as:
115+
116+
\f[
117+
A A^{-1} = A^{-1} A = I
118+
\f]
119+
120+
where \f$ I \f$ is the identity matrix of the same size as \f$ A \f$. The inverse exists only if \f$ A \f$ is square and has full rank (i.e., all its singular values are nonzero).
120121

121-
**Effect**: `A` is replaced with $A^{-1}$.
122+
The computation is performed using LU decomposition.
122123

123-
## [pseudoinvert](@ref la_pseudoinverse::pseudoinvert) - Moore-Penrose pseudo-inverse of a matrix.
124+
### Arguments
125+
126+
- `a`: A `real` or `complex` square matrix of size \f$ [n,n] \f$, representing the matrix to be inverted.
127+
- `err` (optional): A [type(la_state)](@ref la_state_type::la_state) variable that returns the error state. If not provided, the function will stop execution on error.
128+
129+
### Return value
130+
131+
- `inv_a`: A `real` or `complex` square matrix of size \f$ [n,n] \f$, representing the inverse of `a`.
132+
133+
### Errors
134+
135+
- Raises [LINALG_VALUE_ERROR](@ref la_state_type::linalg_value_error) if `a` is singular or has invalid size.
136+
- If `err` is not provided, exceptions will trigger an `error stop`.
137+
138+
### Notes
139+
140+
- This function computes the inverse using LAPACK's LU decomposition routine [GETRF](@ref la_lapack::getrf) followed by [GETRI](@ref la_lapack::getri).
141+
- The inverse should be used with caution in numerical computations. For solving linear systems, using [solve](@ref la_solve::solve) is usually more stable and efficient than explicitly computing the inverse.
142+
143+
## [invert](@ref la_inverse::invert) - In-place matrix inversion
124144

125145
### Syntax
126146

127-
`call pseudoinvert(a, pinva [, rtol] [, err])`
147+
`call invert(a [, err])`
128148

129149
### Description
130150

131-
This subroutine computes the Moore-Penrose pseudo-inverse \f$ A^+ \f$ of a real or complex matrix \f$ A \f$ using Singular Value Decomposition (SVD). The pseudo-inverse is a generalization of the matrix inverse that can be computed for non-square and singular matrices. It is particularly useful for solving least-squares problems and underdetermined systems.
151+
This subroutine computes the inverse \\( A^{-1} \\) of a real or complex square matrix \\( A \\) **in-place**, modifying `a` directly. It uses the LU decomposition method via LAPACK's [GETRF](@ref la_lapack::getrf) and [GETRI](@ref la_lapack::getri) routines.
132152

133-
The computation is based on the singular value decomposition (SVD):
153+
Given a square matrix \\( A \\), the LU decomposition factorizes it as:
134154

135155
\f[
136-
A = U \Sigma V^T
156+
A = P L U
137157
\f]
138158

139-
where \f$ U \f$ and \f$ V \f$ are orthogonal matrices, and \f$ \Sigma \f$ is a diagonal matrix containing the singular values. The pseudo-inverse is computed as:
159+
where:
160+
- \\( P \\) is a permutation matrix,
161+
- \\( L \\) is a lower triangular matrix with unit diagonal,
162+
- \\( U \\) is an upper triangular matrix.
163+
164+
The inverse is then obtained by solving \\( A X = I \\) using the LU factors.
165+
166+
### Arguments
167+
168+
- `a`: A `real` or `complex` square matrix of size \\( [n,n] \\). On output, it is replaced with its inverse \\( A^{-1} \\).
169+
- `err` (optional): A [type(la_state)](@ref la_state_type::la_state) variable that returns the error state. If not provided, the function will stop execution on error.
170+
171+
### Errors
172+
173+
- Raises [LINALG_ERROR](@ref la_state_type::linalg_error) if the matrix is singular.
174+
- Raises [LINALG_VALUE_ERROR](@ref la_state_type::linalg_value_error) if `a` has invalid size.
175+
- If `err` is not provided, exceptions will trigger an `error stop`.
176+
177+
### Notes
178+
179+
- This subroutine modifies `a` in-place. If the original matrix needs to be preserved, use [inv](@ref la_inverse::inv) instead.
180+
- The determinant of `a` can be computed before inversion using [det](@ref la_determinant::det) to check for singularity.
181+
- The computational complexity is \\( O(n^3) \\), making it expensive for large matrices.
182+
- It is recommended to use matrix factorizations (e.g., LU or QR) for solving linear systems instead of computing the inverse explicitly, as it is numerically more stable and efficient.
183+
184+
185+
## [operator(.inv.)](@ref la_inverse::operator(.inv.)) - Compute the inverse of a square matrix.
186+
187+
### Syntax
188+
189+
```fortran
190+
invA = .inv. A
191+
```
192+
193+
### Description
194+
195+
This operator computes the inverse \f$ A^{-1} \f$ of a square, non-singular real or complex matrix \f$ A \f$ using an LU decomposition. The inversion satisfies:
140196

141197
\f[
142-
A^+ = V \Sigma^+ U^T
198+
A A^{-1} = I
143199
\f]
144200

145-
where \f$ \Sigma^+ \f$ is obtained by inverting the nonzero singular values.
201+
where \f$ I \f$ is the identity matrix of appropriate size.
202+
203+
This operator is functionally equivalent to [inv](@ref la_inverse::inv) but provides a more convenient syntax. It supports operator chaining, allowing multiple inversions within expressions:
146204

147205
### Arguments
148206

149-
- `a`: A `real` or `complex` matrix of size \f$ [m,n] \f$, representing the input matrix to be inverted.
150-
- `pinva`: A `real` or `complex` matrix of size \f$ [n,m] \f$, representing the pseudo-inverse of `a`. This is an output argument.
151-
- `rtol` (optional): A real scalar specifying the relative tolerance for singular value truncation. Singular values smaller than `rtol * max(singular_values(A))` are set to zero. If not provided, a default machine-precision-based tolerance is used.
152-
- `err` (optional): A [type(la_state)](@ref la_state_type::la_state) variable that returns the error state. If not provided, the function will stop execution on error.
207+
- `A`: A `real` or `complex` square matrix of size \f$ [n,n] \f$, representing the input matrix to be inverted.
153208

154209
### Return value
155210

156-
The pseudo-inverse of the input matrix `a` is returned in `pinva`.
211+
- `invA`: A `real` or `complex` square matrix of size \f$ [n,n] \f$, and same kind as `A` representing its inverse.
212+
- If `A` is singular or the inversion fails, an **empty matrix** (size \f$ [0,0] \f$) is returned instead of raising an error.
157213

158214
### Errors
159215

160-
- Raises [LINALG_VALUE_ERROR](@ref la_state_type::linalg_value_error) if the dimensions of `pinva` do not match the expected output size.
161-
- Raises [LINALG_ERROR](@ref la_state_type::linalg_error) if the SVD decomposition fails.
162-
- If `err` is not provided, exceptions will trigger an `error stop`.
216+
- Unlike [inv](@ref la_inverse::inv), this operator **does not provide explicit error handling**.
217+
- If `A` is singular or an error occurs during inversion, the function **returns an empty matrix** (size \f$ [0,0] \f$) instead of raising an exception.
218+
- The caller should check the size of the returned matrix to determine if inversion was successful.
163219

164220
### Notes
165221

166-
- This subroutine computes the pseudo-inverse using LAPACK’s SVD decomposition routine [`*GESVD`](@ref la_lapack::gesvd).
167-
- The choice of `rtol` affects numerical stability and rank estimation: setting it too high may result in an inaccurate inverse, while setting it too low may amplify numerical noise.
168-
- This version requires `pinva` to be pre-allocated. To obtain the required size before allocation, use [`pseudoinvert_space`](@ref la_pseudoinvert::pseudoinvert_space).
169-
222+
- This operator internally calls LAPACK's LU decomposition routine [GETRF](@ref la_lapack::getrf) followed by [GETRI](@ref la_lapack::getri).
223+
- The chaining property allows for concise expressions but requires caution: if any intermediate inversion fails, subsequent operations may propagate errors due to empty matrix results.
224+
- If strict error handling is required, use [inv](@ref la_inverse::inv) instead.
170225

171226
## [pinv](@ref la_pseudoinverse::pinv) - Moore-Penrose pseudo-inverse of a matrix (function).
172227

@@ -210,7 +265,7 @@ where \f$ \Sigma^+ \f$ is obtained by inverting the nonzero singular values.
210265

211266
### Notes
212267

213-
- This function computes the pseudo-inverse using LAPACK’s SVD decomposition routine [`*GESVD`](@ref la_lapack::gesvd).
268+
- This function computes the pseudo-inverse using LAPACK's SVD decomposition routine [`*GESVD`](@ref la_lapack::gesvd).
214269
- The choice of `rtol` affects numerical stability and rank estimation: setting it too high may result in an inaccurate inverse, while setting it too low may amplify numerical noise.
215270
- This function returns a newly allocated matrix. For an in-place version, use [`pseudoinvert`](@ref la_pseudoinverse::pseudoinvert).
216271

@@ -255,7 +310,7 @@ where \f$ \Sigma^+ \f$ is the inverse of the nonzero singular values in \f$ \Sig
255310
### Notes
256311

257312
- This operator internally calls [pinv](@ref la_pseudoinverse::pinv) and behaves identically.
258-
- The pseudo-inverse is computed using LAPACK’s SVD decomposition routine [GESVD](@ref la_lapack::gesvd).
313+
- The pseudo-inverse is computed using LAPACK's SVD decomposition routine [GESVD](@ref la_lapack::gesvd).
259314
- This operator is a convenient shorthand for calling the functional interface `pinv(a)`.
260315

261316

‎fypp/src/la_inverse.fypp

Lines changed: 46 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
#:include "common.fypp"
2-
! Compute matrix inverse
2+
!> Compute matrix inverse
33
module la_inverse
44
use la_constants
55
use la_blas
@@ -9,16 +9,54 @@ module la_inverse
99
implicit none(type,external)
1010
private
1111

12-
!> Function interface return the matrix inverse
12+
!> @brief Compute the inverse of a square matrix.
13+
!!
14+
!! This function computes the inverse of a real or complex square matrix \f$ A \f$.
15+
!! The inverse is computed using an LU decomposition with partial pivoting.
16+
!!
17+
!! @param[in] A The input square matrix of size \f$ [n, n] \f$.
18+
!! @param[out] err (Optional) A state return flag. If an error occurs and `err` is not provided,
19+
!! the function will stop execution.
20+
!!
21+
!! @return The inverse matrix \f$ A^{-1} \f$ of size \f$ [n, n] \f$.
22+
!!
23+
!! @note This function relies on LAPACK's LU decomposition routines ([GETRF](@ref la_lapack::getrf)
24+
!! and [GETRI](@ref la_lapack::getri)).
25+
!! @warning The matrix \f$ A \f$ must be non-singular. If it is singular or nearly singular,
26+
!! the function will fail.
27+
!!
1328
public :: inv
14-
!> Subroutine interface: invert matrix inplace
29+
30+
!> @brief Compute the inverse of a square matrix in-place.
31+
!!
32+
!! This subroutine computes the inverse of a real or complex square matrix \f$ A \f$ in-place.
33+
!! The inverse is computed using an LU decomposition with partial pivoting.
34+
!!
35+
!! @param[in,out] A The input square matrix of size \f$ [n, n] \f$. It is replaced by its inverse \f$ A^{-1} \f$.
36+
!! @param[out] err (Optional) A state return flag. If an error occurs and `err` is not provided,
37+
!! the function will stop execution.
38+
!!
39+
!! @note This subroutine is useful when memory efficiency is a priority, as it avoids additional allocations.
40+
!! @warning The matrix \f$ A \f$ must be non-singular. If it is singular or nearly singular,
41+
!! the computation will fail.
42+
!!
1543
public :: invert
16-
!> Operator interface: .inv.A returns the matrix inverse of A
17-
public :: operator(.inv.)
1844

19-
! Numpy: inv(a)
20-
! Scipy: inv(a, overwrite_a=False, check_finite=True)
21-
! IMSL: .i.a
45+
!> @brief Compute the inverse of a square matrix using the `.inv.` operator.
46+
!!
47+
!! This operator computes the inverse of a real or complex square matrix \f$ A \f$ using
48+
!! an LU decomposition with partial pivoting.
49+
!!
50+
!! @param[in] A The input square matrix of size \f$ [n, n] \f$.
51+
!!
52+
!! @return The inverse matrix \f$ A^{-1} \f$ of size \f$ [n, n] \f$.
53+
!!
54+
!! @note This operator is a shorthand for calling `inv(A)`, allowing expressions such as:
55+
!! \f$ X = .inv.A \f$
56+
!! @warning The matrix \f$ A \f$ must be non-singular. If it is singular or nearly singular,
57+
!! the computation will fail.
58+
!!
59+
public :: operator(.inv.)
2260

2361
! Function interface
2462
interface inv

‎src/la_inverse.f90

Lines changed: 48 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
! Compute matrix inverse
1+
!> Inverse of a square matrix
22
module la_inverse
33
use la_constants
44
use la_blas
@@ -8,16 +8,56 @@ module la_inverse
88
implicit none(type,external)
99
private
1010

11-
!> Function interface return the matrix inverse
11+
!> @brief Compute the inverse of a square matrix.
12+
!!
13+
!! This function computes the inverse of a real or complex square matrix \f$ A \f$.
14+
!! The inverse is computed using an LU decomposition with partial pivoting.
15+
!!
16+
!! @param[in] A The input square matrix of size \f$ [n, n] \f$.
17+
!! @param[out] err (Optional) A state return flag. If an error occurs and `err` is not provided,
18+
!! the function will stop execution.
19+
!!
20+
!! @return The inverse matrix \f$ A^{-1} \f$ of size \f$ [n, n] \f$.
21+
!!
22+
!! @note This function relies on LAPACK's LU decomposition routines ([GETRF](@ref la_lapack::getrf)
23+
!! and [GETRI](@ref la_lapack::getri)).
24+
!! @warning The matrix \f$ A \f$ must be non-singular. If it is singular or nearly singular,
25+
!! the function will fail.
26+
!!
1227
public :: inv
13-
!> Subroutine interface: invert matrix inplace
28+
29+
30+
!> @brief Compute the inverse of a square matrix in-place.
31+
!!
32+
!! This subroutine computes the inverse of a real or complex square matrix \f$ A \f$ in-place.
33+
!! The inverse is computed using an LU decomposition with partial pivoting.
34+
!!
35+
!! @param[in,out] A The input square matrix of size \f$ [n, n] \f$. It is replaced by its inverse \f$ A^{-1} \f$.
36+
!! @param[out] err (Optional) A state return flag. If an error occurs and `err` is not provided,
37+
!! the function will stop execution.
38+
!!
39+
!! @note This subroutine is useful when memory efficiency is a priority, as it avoids additional allocations.
40+
!! @warning The matrix \f$ A \f$ must be non-singular. If it is singular or nearly singular,
41+
!! the computation will fail.
42+
!!
1443
public :: invert
15-
!> Operator interface: .inv.A returns the matrix inverse of A
16-
public :: operator(.inv.)
1744

18-
! Numpy: inv(a)
19-
! Scipy: inv(a, overwrite_a=False, check_finite=True)
20-
! IMSL: .i.a
45+
46+
!> @brief Compute the inverse of a square matrix using the `.inv.` operator.
47+
!!
48+
!! This operator computes the inverse of a real or complex square matrix \f$ A \f$ using
49+
!! an LU decomposition with partial pivoting.
50+
!!
51+
!! @param[in] A The input square matrix of size \f$ [n, n] \f$.
52+
!!
53+
!! @return The inverse matrix \f$ A^{-1} \f$ of size \f$ [n, n] \f$.
54+
!!
55+
!! @note This operator is a shorthand for calling `inv(A)`, allowing expressions such as:
56+
!! \f$ X = .inv.A \f$
57+
!! @warning The matrix \f$ A \f$ must be non-singular. If it is singular or nearly singular,
58+
!! the computation will fail.
59+
!!
60+
public :: operator(.inv.)
2161

2262
! Function interface
2363
interface inv

0 commit comments

Comments
 (0)