Skip to content

Commit ccdba91

Browse files
perazzperazz
authored
linalg: QR factorization (#832)
2 parents 2b4d8b2 + d3ae5a7 commit ccdba91

9 files changed

+601
-1
lines changed

‎doc/specs/stdlib_linalg.md

+78-1
Original file line numberDiff line numberDiff line change
@@ -902,6 +902,83 @@ Exceptions trigger an `error stop`.
902902
{!example/linalg/example_determinant2.f90!}
903903
```
904904

905+
## `qr` - Compute the QR factorization of a matrix
906+
907+
### Status
908+
909+
Experimental
910+
911+
### Description
912+
913+
This subroutine computes the QR factorization of a `real` or `complex` matrix: \( A = Q R \) where \( Q \)
914+
is orthonormal and \( R \) is upper-triangular. Matrix \( A \) has size `[m,n]`, with \( m \ge n \).
915+
916+
The results are returned in output matrices \( Q \) and \(R \), that have the same type and kind as \( A \).
917+
Given `k = min(m,n)`, one can write \( A = \( Q_1 Q_2 \) \cdot \( \frac{R_1}{0}\) \).
918+
Because the lower rows of \( R \) are zeros, a reduced problem \( A = Q_1 R_1 \) may be solved. The size of
919+
the input arguments determines what problem is solved: on full matrices (`shape(Q)==[m,m]`, `shape(R)==[m,n]`),
920+
the full problem is solved. On reduced matrices (`shape(Q)==[m,k]`, `shape(R)==[k,n]`), the reduced problem is solved.
921+
922+
### Syntax
923+
924+
`call ` [[stdlib_linalg(module):qr(interface)]] `(a, q, r, [, storage] [, overwrite_a] [, err])`
925+
926+
### Arguments
927+
928+
`a`: Shall be a rank-2 `real` or `complex` array containing the coefficient matrix of size `[m,n]`. It is an `intent(in)` argument, if `overwrite_a=.false.`. Otherwise, it is an `intent(inout)` argument, and is destroyed upon return.
929+
930+
`q`: Shall be a rank-2 array of the same kind as `a`, containing the orthonormal matrix `q`. It is an `intent(out)` argument. It should have a shape equal to either `[m,m]` or `[m,k]`, whether the full or the reduced problem is sought for.
931+
932+
`r`: Shall be a rank-2 array of the same kind as `a`, containing the upper triangular matrix `r`. It is an `intent(out)` argument. It should have a shape equal to either `[m,n]` or `[k,n]`, whether the full or the reduced problem is sought for.
933+
934+
`storage` (optional): Shall be a rank-1 array of the same type and kind as `a`, providing working storage for the solver. Its minimum size can be determined with a call to [[stdlib_linalg(module):qr_space(interface)]]. It is an `intent(out)` argument.
935+
936+
`overwrite_a` (optional): Shall be an input `logical` flag (default: `.false.`). If `.true.`, input matrix `a` will be used as temporary storage and overwritten. This avoids internal data allocation. It is an `intent(in)` argument.
937+
938+
`err` (optional): Shall be a `type(linalg_state_type)` value. It is an `intent(out)` argument.
939+
940+
### Return value
941+
942+
Returns the QR factorization matrices into the \( Q \) and \( R \) arguments.
943+
944+
Raises `LINALG_VALUE_ERROR` if any of the matrices has invalid or unsuitable size for the full/reduced problem.
945+
Raises `LINALG_ERROR` on insufficient user storage space.
946+
If the state argument `err` is not present, exceptions trigger an `error stop`.
947+
948+
### Example
949+
950+
```fortran
951+
{!example/linalg/example_qr.f90!}
952+
```
953+
954+
## `qr_space` - Compute internal working space requirements for the QR factorization.
955+
956+
### Status
957+
958+
Experimental
959+
960+
### Description
961+
962+
This subroutine computes the internal working space requirements for the QR factorization, [[stdlib_linalg(module):qr(interface)]] .
963+
964+
### Syntax
965+
966+
`call ` [[stdlib_linalg(module):qr_space(interface)]] `(a, lwork, [, err])`
967+
968+
### Arguments
969+
970+
`a`: Shall be a rank-2 `real` or `complex` array containing the coefficient matrix. It is an `intent(in)` argument.
971+
972+
`lwork`: Shall be an `integer` scalar, that returns the minimum array size required for the working storage in [[stdlib_linalg(module):qr(interface)]] to factorize `a`.
973+
974+
`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument.
975+
976+
### Example
977+
978+
```fortran
979+
{!example/linalg/example_qr_space.f90!}
980+
```
981+
905982
## `eig` - Eigenvalues and Eigenvectors of a Square Matrix
906983

907984
### Status
@@ -1028,7 +1105,6 @@ Raises `LINALG_ERROR` if the calculation did not converge.
10281105
Raises `LINALG_VALUE_ERROR` if any matrix or arrays have invalid/incompatible sizes.
10291106
If `err` is not present, exceptions trigger an `error stop`.
10301107

1031-
10321108
### Example
10331109

10341110
```fortran
@@ -1096,6 +1172,7 @@ If requested, `vt` contains the right singular vectors, as rows of \( V^T \).
10961172
`call ` [[stdlib_linalg(module):svd(interface)]] `(a, s, [, u, vt, overwrite_a, full_matrices, err])`
10971173

10981174
### Class
1175+
10991176
Subroutine
11001177

11011178
### Arguments

‎example/linalg/CMakeLists.txt

+2
Original file line numberDiff line numberDiff line change
@@ -35,5 +35,7 @@ ADD_EXAMPLE(svd)
3535
ADD_EXAMPLE(svdvals)
3636
ADD_EXAMPLE(determinant)
3737
ADD_EXAMPLE(determinant2)
38+
ADD_EXAMPLE(qr)
39+
ADD_EXAMPLE(qr_space)
3840
ADD_EXAMPLE(cholesky)
3941
ADD_EXAMPLE(chol)

‎example/linalg/example_qr.f90

+15
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,15 @@
1+
program example_qr
2+
use stdlib_linalg, only: qr
3+
implicit none(type,external)
4+
real :: A(104, 32), Q(104,32), R(32,32)
5+
6+
! Create a random matrix
7+
call random_number(A)
8+
9+
! Compute its QR factorization (reduced)
10+
call qr(A,Q,R)
11+
12+
! Test factorization: Q*R = A
13+
print *, maxval(abs(matmul(Q,R)-A))
14+
15+
end program example_qr

‎example/linalg/example_qr_space.f90

+25
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,25 @@
1+
! QR example with pre-allocated storage
2+
program example_qr_space
3+
use stdlib_linalg_constants, only: ilp
4+
use stdlib_linalg, only: qr, qr_space, linalg_state_type
5+
implicit none(type,external)
6+
real :: A(104, 32), Q(104,32), R(32,32)
7+
real, allocatable :: work(:)
8+
integer(ilp) :: lwork
9+
type(linalg_state_type) :: err
10+
11+
! Create a random matrix
12+
call random_number(A)
13+
14+
! Prepare QR workspace
15+
call qr_space(A,lwork)
16+
allocate(work(lwork))
17+
18+
! Compute its QR factorization (reduced)
19+
call qr(A,Q,R,storage=work,err=err)
20+
21+
! Test factorization: Q*R = A
22+
print *, maxval(abs(matmul(Q,R)-A))
23+
print *, err%print()
24+
25+
end program example_qr_space

‎src/CMakeLists.txt

+1
Original file line numberDiff line numberDiff line change
@@ -30,6 +30,7 @@ set(fppFiles
3030
stdlib_linalg_eigenvalues.fypp
3131
stdlib_linalg_solve.fypp
3232
stdlib_linalg_determinant.fypp
33+
stdlib_linalg_qr.fypp
3334
stdlib_linalg_inverse.fypp
3435
stdlib_linalg_state.fypp
3536
stdlib_linalg_svd.fypp

‎src/stdlib_linalg.fypp

+71
Original file line numberDiff line numberDiff line change
@@ -40,6 +40,8 @@ module stdlib_linalg
4040
public :: outer_product
4141
public :: kronecker_product
4242
public :: cross_product
43+
public :: qr
44+
public :: qr_space
4345
public :: is_square
4446
public :: is_diagonal
4547
public :: is_symmetric
@@ -526,6 +528,75 @@ module stdlib_linalg
526528
#:endfor
527529
end interface lstsq_space
528530

531+
! QR factorization of rank-2 array A
532+
interface qr
533+
!! version: experimental
534+
!!
535+
!! Computes the QR factorization of matrix \( A = Q R \).
536+
!! ([Specification](../page/specs/stdlib_linalg.html#qr-compute-the-qr-factorization-of-a-matrix))
537+
!!
538+
!!### Summary
539+
!! Compute the QR factorization of a `real` or `complex` matrix: \( A = Q R \), where \( Q \) is orthonormal
540+
!! and \( R \) is upper-triangular. Matrix \( A \) has size `[m,n]`, with \( m\ge n \).
541+
!!
542+
!!### Description
543+
!!
544+
!! This interface provides methods for computing the QR factorization of a matrix.
545+
!! Supported data types include `real` and `complex`. If a pre-allocated work space
546+
!! is provided, no internal memory allocations take place when using this interface.
547+
!!
548+
!! Given `k = min(m,n)`, one can write \( A = \( Q_1 Q_2 \) \cdot \( \frac{R_1}{0}\) \).
549+
!! The user may want the full problem (provide `shape(Q)==[m,m]`, `shape(R)==[m,n]`) or the reduced
550+
!! problem only: \( A = Q_1 R_1 \) (provide `shape(Q)==[m,k]`, `shape(R)==[k,n]`).
551+
!!
552+
!!@note The solution is based on LAPACK's QR factorization (`*GEQRF`) and ordered matrix output (`*ORGQR`, `*UNGQR`).
553+
!!
554+
#:for rk,rt,ri in RC_KINDS_TYPES
555+
pure module subroutine stdlib_linalg_${ri}$_qr(a,q,r,overwrite_a,storage,err)
556+
!> Input matrix a[m,n]
557+
${rt}$, intent(inout), target :: a(:,:)
558+
!> Orthogonal matrix Q ([m,m], or [m,k] if reduced)
559+
${rt}$, intent(out), contiguous, target :: q(:,:)
560+
!> Upper triangular matrix R ([m,n], or [k,n] if reduced)
561+
${rt}$, intent(out), contiguous, target :: r(:,:)
562+
!> [optional] Can A data be overwritten and destroyed?
563+
logical(lk), optional, intent(in) :: overwrite_a
564+
!> [optional] Provide pre-allocated workspace, size to be checked with qr_space
565+
${rt}$, intent(out), optional, target :: storage(:)
566+
!> [optional] state return flag. On error if not requested, the code will stop
567+
type(linalg_state_type), optional, intent(out) :: err
568+
end subroutine stdlib_linalg_${ri}$_qr
569+
#:endfor
570+
end interface qr
571+
572+
! Return the working array space required by the QR factorization solver
573+
interface qr_space
574+
!! version: experimental
575+
!!
576+
!! Computes the working array space required by the QR factorization solver
577+
!! ([Specification](../page/specs/stdlib_linalg.html#qr-space-compute-internal-working-space-requirements-for-the-qr-factorization))
578+
!!
579+
!!### Description
580+
!!
581+
!! This interface returns the size of the `real` or `complex` working storage required by the
582+
!! QR factorization solver. The working size only depends on the kind (`real` or `complex`) and size of
583+
!! the matrix being factorized. Storage size can be used to pre-allocate a working array in case several
584+
!! repeated QR factorizations to a same-size matrix are sought. If pre-allocated working arrays
585+
!! are provided, no internal allocations will take place during the factorization.
586+
!!
587+
#:for rk,rt,ri in RC_KINDS_TYPES
588+
pure module subroutine get_qr_${ri}$_workspace(a,lwork,err)
589+
!> Input matrix a[m,n]
590+
${rt}$, intent(in), target :: a(:,:)
591+
!> Minimum workspace size for both operations
592+
integer(ilp), intent(out) :: lwork
593+
!> State return flag. Returns an error if the query failed
594+
type(linalg_state_type), optional, intent(out) :: err
595+
end subroutine get_qr_${ri}$_workspace
596+
#:endfor
597+
end interface qr_space
598+
599+
529600
interface det
530601
!! version: experimental
531602
!!

0 commit comments

Comments
 (0)