Skip to content

Commit 395057c

Browse files
perazzperazz
authored
More docs: eye, diag, schur (#74)
2 parents 8a030ea + 29519c9 commit 395057c

File tree

5 files changed

+309
-126
lines changed

5 files changed

+309
-126
lines changed

‎README.md

Lines changed: 145 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -331,13 +331,73 @@ where \f$ \Sigma^+ \f$ is the inverse of the nonzero singular values in \f$ \Sig
331331
**Description**: Singular values $S$ from $A = U S V^t$.
332332
**Usage**: `s = svdvals(A)` where `s` is a real array with the same precision as `A`.
333333

334-
## `eye(m)`
335-
**Type**: Function
336-
**Description**: Identity matrix of size `m`.
337-
**Optional arguments**:
338-
- `n`: Optional column size.
339-
- `mold`: Optional datatype (default: real64).
340-
- `err`: Error handler.
334+
## [diag](@ref la_eye::diag) - Diagonal matrix.
335+
336+
### Syntax
337+
338+
`d = diag(n, source [, err])` for scalar input
339+
`d = diag(source(:) [, err])` for array input
340+
341+
### Description
342+
343+
This function generates a square diagonal matrix where the diagonal elements are populated either by a scalar value or an array of values. The size of the matrix is determined by the input parameter \f$n\f$ or the size of the input array.
344+
If a scalar is provided, the diagonal elements are all set to the same value. If an array is provided, its length determines the size of the matrix, and its elements are placed along the diagonal.
345+
346+
### Arguments
347+
348+
- `n`: The size of the square matrix (only used if a scalar is provided for the diagonal).
349+
- `source`:
350+
- If a scalar, this value is used to populate all the diagonal elements of the matrix.
351+
- If an array, the elements of the array are used to populate the diagonal of the matrix. The size of the array determines the matrix size.
352+
- `err` (optional): A state return flag of [type(la_state)](@ref la_state_type::la_state). If an error occurs and `err` is not provided, the function will stop execution.
353+
354+
### Return value
355+
356+
The function returns a matrix of size \f$n \times n\f$, where the diagonal elements are either all equal to the scalar `source` or populated by the values from the input array.
357+
358+
### Errors
359+
360+
- Raises [LINALG_VALUE_ERROR](@ref la_state_type::linalg_value_error) if the dimensions of the matrix are invalid or if the array size does not match the expected matrix size.
361+
- If `err` is not provided, the function will stop execution on errors.
362+
363+
### Notes
364+
365+
- The diagonal elements are set to the specified scalar or the array values in the order they appear in the input.
366+
- If the `err` parameter is provided, the error state of the function will be returned.
367+
368+
369+
## [eye](@ref la_eye::eye) - Identity matrix.
370+
371+
### Syntax
372+
373+
`eye = eye(m [, n] [, mold] [, err])`
374+
375+
### Description
376+
377+
This function constructs an identity matrix of size \f$m \times n\f$, where the diagonal elements are set to 1 and all off-diagonal elements are set to 0. If only the number of rows \f$m\f$ is provided, a square matrix of size \f$m \times m\f$ is returned. The matrix is populated with a real data type, by default `real(real64)`, or a type specified by the user.
378+
379+
### Arguments
380+
381+
- `m`: The number of rows of the identity matrix.
382+
- `n` (optional): The number of columns of the identity matrix. If omitted, the matrix is square (\f$m \times m\f$).
383+
- `mold` (optional): The data type to define the return type. Defaults to `real(real64)`.
384+
- `err` (optional): A state return flag of [type(la_state)](@ref la_state_type::la_state). If an error occurs and `err` is not provided, the function will stop execution.
385+
386+
### Return value
387+
388+
The function returns a matrix of size \f$m \times n\f$ (or \f$m \times m\f$ if \f$n\f$ is omitted) with diagonal elements set to 1 and all other elements set to 0.
389+
390+
### Errors
391+
392+
- Raises [LINALG_VALUE_ERROR](@ref la_state_type::linalg_value_error) if the dimensions of the matrix are invalid (e.g., negative values).
393+
- If `err` is not provided, the function will stop execution on errors.
394+
395+
### Notes
396+
397+
- The identity matrix is constructed with the specified data type, which defaults to `real(real64)` if no type is specified.
398+
- The `mold` scalar is used to provide a function return type.
399+
- If the `err` parameter is provided, the error state of the function will be returned.
400+
341401

342402
## `eigvals(A)`
343403
**Type**: Function
@@ -370,19 +430,9 @@ where \f$ \Sigma^+ \f$ is the inverse of the nonzero singular values in \f$ \Sig
370430
- `overwrite_a`: Option to let A be destroyed.
371431
- `err`: Return state handler.
372432

373-
## `diag(n, source)`
374-
**Type**: Function
375-
**Description**: Diagonal matrix from scalar input value.
376-
**Optional arguments**:
377-
- `err`: Error handler.
378433

379-
## `diag(source)`
380-
**Type**: Function
381-
**Description**: Diagonal matrix from array input values.
382-
**Optional arguments**:
383-
- `err`: Error handler.
384434

385-
## [qr](@ref la_qr::qr) - Compute the QR factorization of a matrix.
435+
## [qr](@ref la_qr::qr) - QR factorization of a matrix.
386436

387437
### Syntax
388438

@@ -459,6 +509,83 @@ The workspace size \f$ lwork \f$ that should be allocated before calling the QR
459509
- This subroutine is useful for preallocating memory for QR factorization in large systems.
460510
- It is important to ensure that the workspace size is correctly allocated before proceeding with QR factorization to avoid memory issues.
461511

512+
## [schur](@ref la_schur::schur) - Schur decomposition of a matrix.
513+
514+
### Syntax
515+
516+
`call schur(a, t, z [, eigvals] [, overwrite_a] [, storage] [, err])`
517+
518+
### Description
519+
520+
This subroutine computes the Schur decomposition of a `real` or `complex` matrix \f$ A = Z T Z^H \f$, where \f$ Z \f$ is an orthonormal/unitary matrix, and \f$ T \f$ is an upper-triangular or quasi-upper-triangular matrix. The matrix \f$ A \f$ has size \f$ [m,m] \f$.
521+
522+
The decomposition produces:
523+
- \f$ T \f$, which is upper-triangular for `complex` matrices and quasi-upper-triangular for `real` matrices (with possible \f$ 2 \times 2 \f$ blocks on the diagonal).
524+
- \f$ Z \f$, the transformation matrix, which is optional.
525+
- Optionally, the eigenvalues corresponding to the diagonal elements of \f$ T \f$.
526+
527+
If a pre-allocated workspace is provided, no internal memory allocations take place.
528+
529+
### Arguments
530+
531+
- `a`: A `real` or `complex` matrix of size \f$ [m,m] \f$. If `overwrite_a = .false.`, this is an input argument. If `overwrite_a = .true.`, it is an `inout` argument and is overwritten upon return.
532+
- `t`: A rank-2 array of the same type and kind as `a`, representing the Schur form of `a`. This is an output argument with shape \f$ [m,m] \f$.
533+
- `z` (optional): A rank-2 array of the same type and kind as `a`, representing the unitary/orthonormal transformation matrix \f$ Z \f$. This is an output argument with shape \f$ [m,m] \f$.
534+
- `eigvals` (optional): A complex array of size \f$ [m] \f$, representing the eigenvalues that appear on the diagonal of \f$ T \f$. This is an output argument.
535+
- `storage` (optional): A rank-1 array of the same type and kind as `a`, providing working storage for the solver. Its minimum size can be determined by a call to [schur_space](@ref la_schur::schur_space). This is an input argument.
536+
- `overwrite_a` (optional, default = `.false.`): A logical flag that determines whether the input matrix `a` can be overwritten. If `.true.`, the matrix `a` is used as temporary storage and overwritten to avoid internal memory allocation. This is an input argument.
537+
- `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.
538+
539+
### Return value
540+
541+
The Schur decomposition matrices \f$ T \f$ and optionally \f$ Z \f$ are returned in the corresponding arguments.
542+
543+
### Errors
544+
545+
- Raises [LINALG_VALUE_ERROR](@ref la_state_type::linalg_value_error) if the sizes of the matrices are incompatible.
546+
- Raises [LINALG_ERROR](@ref la_state_type::linalg_error) if the algorithm did not converge.
547+
- If `err` is not provided, exceptions will trigger an `error stop`.
548+
549+
### Notes
550+
551+
- This subroutine computes the Schur decomposition using LAPACK's Schur decomposition routines ([GEES](@ref la_lapack:gees)).
552+
- Sorting options for eigenvalues can be requested, utilizing LAPACK's eigenvalue sorting mechanism.
553+
- If `overwrite_a` is enabled, the input matrix `a` will be modified during computation.
554+
555+
556+
## [schur_space](@ref la_schur::schur_space) - Workspace size for Schur decomposition.
557+
558+
### Syntax
559+
560+
`call schur_space(a, lwork [, err])`
561+
562+
### Description
563+
564+
This subroutine computes the minimum workspace size required for performing Schur decomposition. The size of the workspace array needed is determined based on the input matrix \f$ A \f$.
565+
566+
The input matrix \f$ A \f$ has size \f$ [m,m] \f$, and the output value \f$ lwork \f$ represents the minimum size of the workspace array that should be allocated for Schur decomposition operations.
567+
568+
### Arguments
569+
570+
- `a`: A `real` or `complex` matrix of size \f$ [m,m] \f$, representing the input matrix used to determine the required workspace size.
571+
- `lwork`: An integer variable that will return the minimum workspace size required for Schur decomposition.
572+
- `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.
573+
574+
### Return value
575+
576+
The workspace size \f$ lwork \f$ that should be allocated before calling the Schur decomposition routine is returned.
577+
578+
### Errors
579+
580+
- Raises [LINALG_ERROR](@ref la_state_type::linalg_error) if there is an issue determining the required workspace size.
581+
- If `err` is not provided, exceptions will trigger an `error stop`.
582+
583+
### Notes
584+
585+
- This subroutine is useful for preallocating memory for Schur decomposition in large systems.
586+
- It is important to ensure that the workspace size is correctly allocated before proceeding with Schur decomposition to avoid memory issues.
587+
588+
462589

463590
# BLAS, LAPACK
464591
Modern Fortran modules with full explicit typing features are available as modules `la_blas` and `la_lapack`.

‎fypp/src/la_eye.fypp

Lines changed: 39 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
#:include "common.fypp"
2-
! Return a 2-D matrix with ones on the diagonal and zeros everywhere else
2+
!> Identity and diagonal matrices.
33
module la_eye
44
use la_constants
55
use la_blas
@@ -9,21 +9,45 @@ module la_eye
99
implicit none(type,external)
1010
private
1111

12-
!> Constructs the identity matrix.
13-
!! This interface provides procedures to generate an identity matrix of a given size.
14-
!! The resulting matrix has 1s on the diagonal and 0s elsewhere.
15-
public :: eye
16-
17-
public :: diag
18-
12+
!> @brief Construct an identity matrix of size \f$m \times n\f$.
13+
!!
14+
!! This function returns a diagonal identity matrix of size \f$m \times n\f$, where all diagonal elements
15+
!! are set to 1 and all off-diagonal elements are set to 0. The number of rows and columns can be specified.
16+
!! If only one parameter is provided, a square matrix of size \f$m \times m\f$ is returned.
17+
!!
18+
!! @param[in] m The number of rows of the identity matrix.
19+
!! @param[in] n (Optional) The number of columns of the identity matrix. If omitted, the matrix is square (\f$m \times m\f$).
20+
!! @param[in] mold (Optional) Data type to define the return type. Defaults to `real(real64)`.
21+
!!
22+
!! @return The identity matrix with size \f$m \times n\f$.
23+
!!
24+
!! @note If the `mold` parameter is omitted, the default type is `real(real64)`. If specified, the return type
25+
!! matches the given type.
26+
!!
27+
!! @warning Ensure that the matrix dimensions are valid and consistent with the type definition.
28+
public :: eye
29+
30+
!> @brief Return a square diagonal matrix with diagonal values.
31+
!!
32+
!! This function generates a square diagonal matrix where the diagonal elements are either
33+
!! equal to the specified scalar value or populated by the input array.
34+
!! The size of the matrix is determined by the input parameter \f$n\f$ or the size of the input array.
35+
!!
36+
!! @param[in] n The size of the square matrix (only used if a scalar is provided for the diagonal).
37+
!! @param[in] source If a scalar, this value is used to populate the diagonal of the matrix.
38+
!! If an array, the elements of the array are used for the diagonal.
39+
!! @param[out] err (Optional) State return flag. If not provided, the function will stop on error.
40+
!! @return The diagonal matrix with size \f$n \times n\f$, where the diagonal elements are populated by \f$source\f$.
41+
!!
42+
!! @note If a scalar value is passed, the diagonal elements of the matrix will all be equal to \f$source\f$.
43+
!! If an array is passed, its length determines the size of the matrix, and the array elements are placed
44+
!! along the diagonal. The `err` parameter is optional. If not requested, the code will stop on error.
45+
!! Otherwise, it returns the error state of the function.
46+
!!
47+
public :: diag
48+
49+
! Identity matrix interface
1950
interface eye
20-
!> Procedure for creating an identity matrix.
21-
!! - **Inputs:**
22-
!! - `m` (integer): Number of rows.
23-
!! - `n` (integer, optional): Number of columns. Defaults to `m` if not provided.
24-
!! - `mold` (datatype, optional): Data type used to define the matrix elements. Defaults to `real(real64)`.
25-
!! - **Outputs:**
26-
!! - Identity matrix of specified size and type.
2751
#:for epref,esuf,earg,epresent in FUNCTION_INTERFACES
2852
#:for rk,rt,ri in ALL_KINDS_TYPES
2953
module procedure la_eye_${ri}$${esuf}$

‎fypp/src/la_schur.fypp

Lines changed: 55 additions & 46 deletions
Original file line numberDiff line numberDiff line change
@@ -11,7 +11,7 @@ module la_schur
1111

1212
character(*), parameter :: this = 'schur'
1313

14-
!> List of internal GEES tasks:
14+
! List of internal GEES tasks:
1515

1616
!> No task request
1717
character, parameter :: GEES_NOT = 'N'
@@ -22,58 +22,67 @@ module la_schur
2222
!> Request Schur vectors to be sorted
2323
character, parameter :: GEES_SORTED_VECTORS = 'S'
2424

25-
! Schur decomposition of rank-2 array A
26-
interface schur
27-
!! version: experimental
25+
!> @brief Compute the Schur decomposition of a matrix.
2826
!!
29-
!! Computes the Schur decomposition of matrix \( A = Z T Z^H \).
30-
!! ([Specification](../page/specs/la_linalg.html#schur-compute-the-schur-decomposition-of-a-matrix))
27+
!! This function computes the Schur decomposition of a real or complex matrix \f$ A \f$:
3128
!!
32-
!!### Summary
33-
!! Compute the Schur decomposition of a `real` or `complex` matrix: \( A = Z T Z^H \), where \( Z \) is
34-
!! orthonormal/unitary and \( T \) is upper-triangular or quasi-upper-triangular. Matrix \( A \) has size `[m,m]`.
29+
!! \f$ A = Z T Z^H \f$
3530
!!
36-
!!### Description
37-
!!
38-
!! This interface provides methods for computing the Schur decomposition of a matrix.
39-
!! Supported data types include `real` and `complex`. If a pre-allocated workspace is provided, no internal
40-
!! memory allocations take place when using this interface.
31+
!! where \f$ Z \f$ is an orthonormal/unitary matrix and \f$ T \f$ is upper-triangular or quasi-upper-triangular.
32+
!! The input matrix \f$ A \f$ has size \f$ [m, m] \f$.
4133
!!
42-
!! The output matrix \( T \) is upper-triangular for `complex` input matrices and quasi-upper-triangular
43-
!! for `real` input matrices (with possible \( 2 \times 2 \) blocks on the diagonal).
44-
!! The user can optionally request sorting of eigenvalues based on conditions, or a custom sorting function.
34+
!! The decomposition is available for both real and complex matrices:
35+
!! - For real matrices, the Schur form \f$ T \f$ may contain 2x2 blocks for complex eigenvalues.
36+
!! - For complex matrices, \f$ T \f$ is always upper-triangular.
4537
!!
46-
!!@note The solution is based on LAPACK's Schur decomposition routines (`*GEES`). Sorting options
47-
!! are implemented using LAPACK's eigenvalue sorting mechanism.
48-
!!
49-
#:for rk,rt,ri in ALL_KINDS_TYPES
50-
module procedure la_${ri}$_schur
51-
module procedure la_real_eig_${ri}$_schur
52-
#:endfor
53-
end interface schur
38+
!! @param[in,out] A The input matrix of size \f$ [m, m] \f$. Can be overwritten if `overwrite_a` is set.
39+
!! @param[out] T The upper-triangular or quasi-upper-triangular matrix of size \f$ [m, m] \f$.
40+
!! @param[out] Z (Optional) The unitary/orthonormal transformation matrix of size \f$ [m, m] \f$.
41+
!! @param[out] eigvals (Optional) The eigenvalues that appear on the diagonal of \f$ T \f$.
42+
!! For real matrices, this is a real-valued array.
43+
!! For complex matrices, this is a complex-valued array.
44+
!! @param[in,out] storage (Optional) Pre-allocated workspace array. If provided, no internal allocations occur.
45+
!! @param[in] overwrite_a (Optional) Logical flag indicating whether \f$ A \f$ can be overwritten.
46+
!! @param[out] err (Optional) State return flag. If not provided, the function will stop on error.
47+
!!
48+
!! @note The computation is based on LAPACK's Schur decomposition routines ([GEES](@ref la_lapack::gees)).
49+
!! @warning Ensure that `overwrite_a` is set correctly to avoid unintended modification of \f$ A \f$.
50+
!!
51+
interface schur
52+
#:for rk,rt,ri in ALL_KINDS_TYPES
53+
module procedure la_${ri}$_schur
54+
module procedure la_real_eig_${ri}$_schur
55+
#:endfor
56+
end interface schur
5457

55-
! Return the working array space required by the Schur decomposition solver
56-
interface schur_space
57-
!! version: experimental
58+
59+
interface schur
60+
#:for rk,rt,ri in ALL_KINDS_TYPES
61+
module procedure la_${ri}$_schur
62+
module procedure la_real_eig_${ri}$_schur
63+
#:endfor
64+
end interface schur
65+
66+
!> @brief Compute the required workspace size for Schur decomposition.
5867
!!
59-
!! Computes the working array space required by the Schur decomposition solver
60-
!! ([Specification](../page/specs/la_linalg.html#schur-space-compute-internal-working-space-requirements-for-the-schur-decomposition))
68+
!! This subroutine determines the minimum workspace size required for the Schur decomposition of a matrix.
69+
!! The required size depends on the matrix dimensions and type (`real` or `complex`).
6170
!!
62-
!!### Description
63-
!!
64-
!! This interface returns the size of the `real` or `complex` working storage required by the
65-
!! Schur decomposition solver. The working size only depends on the kind (`real` or `complex`) and size of
66-
!! the matrix being decomposed. Storage size can be used to pre-allocate a working array in case several
67-
!! repeated Schur decompositions of same-size matrices are sought. If pre-allocated working arrays
68-
!! are provided, no internal allocations will take place during the decomposition.
69-
!!
70-
#:for rk,rt,ri in ALL_KINDS_TYPES
71-
module procedure get_schur_${ri}$_workspace
72-
#:endfor
73-
end interface schur_space
74-
75-
76-
71+
!! @param[in] A The input matrix of size \f$ [m, m] \f$.
72+
!! @param[out] lwork The minimum workspace size required for the decomposition.
73+
!! @param[out] err (Optional) State return flag. If provided, it will return an error status in case of failure.
74+
!!
75+
!! @note This routine is useful for pre-allocating workspace when performing multiple decompositions on matrices
76+
!! of the same size, minimizing dynamic memory allocation overhead.
77+
!!
78+
interface schur_space
79+
80+
interface schur_space
81+
#:for rk,rt,ri in ALL_KINDS_TYPES
82+
module procedure get_schur_${ri}$_workspace
83+
#:endfor
84+
end interface schur_space
85+
7786
contains
7887

7988
!> Wrapper function for Schur vectors request
@@ -363,7 +372,7 @@ module la_schur
363372
end subroutine la_${ri}$_schur
364373

365374
! Schur decomposition subroutine: real eigenvalue interface
366-
module subroutine la_real_eig_${ri}$_schur(a,t,z,eigvals,overwrite_a,storage,err)
375+
subroutine la_real_eig_${ri}$_schur(a,t,z,eigvals,overwrite_a,storage,err)
367376
!> Input matrix a[m,m]
368377
${rt}$, intent(inout), target :: a(:,:)
369378
!> Schur form of A: upper-triangular or quasi-upper-triangular matrix T

0 commit comments

Comments
 (0)