Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

linalg: QR factorization #832

Merged
merged 35 commits into from
Sep 18, 2024
Merged
Changes from all commits
Commits
Show all changes
35 commits
Select commit Hold shift + click to select a range
6cac0cc
add `qr` module
perazz Jun 3, 2024
d4b5dbc
add examples
perazz Jun 3, 2024
fe0e243
submodule
perazz Jun 3, 2024
ec6f687
fix the examples
perazz Jun 3, 2024
a6a6838
exclude `xdp`
perazz Jun 3, 2024
feecc14
fix `module` `subroutine`
perazz Jun 3, 2024
f3d74ff
enforce full/reduced problem sizes
perazz Jun 5, 2024
5735b05
add tests
perazz Jun 5, 2024
7c9d5b4
remove `xdp` from tests
perazz Jun 5, 2024
cda4b0d
document `qr_space`
perazz Jun 5, 2024
8df204f
document `qr`
perazz Jun 5, 2024
b4889c5
relax tol
perazz Jun 5, 2024
9ac97d1
deactivate full tol check
perazz Jun 5, 2024
364b9ca
export error amount
perazz Jun 5, 2024
8ed1fe4
Fix: return full `Q` from `*ORGQR`
perazz Jun 5, 2024
3db022d
add test: init NaNs, check correct return
perazz Jun 5, 2024
fcb5306
Merge branch 'master' into qr
perazz Jun 21, 2024
68ece36
Merge branch 'master' into qr
perazz Jul 6, 2024
86d15bc
lift `xdp` restriction
perazz Jul 6, 2024
381d024
Update doc/specs/stdlib_linalg.md
perazz Jul 8, 2024
931f07a
Update doc/specs/stdlib_linalg.md
perazz Jul 8, 2024
5f15913
Update doc/specs/stdlib_linalg.md
perazz Jul 8, 2024
997021c
Update doc/specs/stdlib_linalg.md
perazz Jul 8, 2024
789e6d0
Update doc/specs/stdlib_linalg.md
perazz Jul 8, 2024
e5999f2
Update doc/specs/stdlib_linalg.md
perazz Jul 8, 2024
47f70e6
Update doc/specs/stdlib_linalg.md
perazz Jul 8, 2024
478f68d
Update doc/specs/stdlib_linalg.md
perazz Jul 8, 2024
e3db36d
improve qr description
perazz Jul 8, 2024
e69e767
`qr` workspace: make `a` `intent(in)`
perazz Jul 8, 2024
8827bf7
`qr` temporary storage: make `intent(out)`
perazz Jul 8, 2024
7393a18
Merge branch 'master' into qr
perazz Jul 8, 2024
c6bb171
Update src/stdlib_linalg_qr.fypp
perazz Jul 11, 2024
6fc0e1d
use IEEE NaNs
perazz Jul 11, 2024
df4016a
Merge branch 'master' into qr
perazz Aug 19, 2024
d3ae5a7
Merge branch 'fortran-lang:master' into qr
perazz Aug 23, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
79 changes: 78 additions & 1 deletion doc/specs/stdlib_linalg.md
Original file line number Diff line number Diff line change
@@ -902,6 +902,83 @@ Exceptions trigger an `error stop`.
{!example/linalg/example_determinant2.f90!}
```

## `qr` - Compute the QR factorization of a matrix

### Status

Experimental

### Description

This subroutine computes the QR factorization of a `real` or `complex` matrix: \( A = Q R \) where \( Q \)
is orthonormal and \( R \) is upper-triangular. Matrix \( A \) has size `[m,n]`, with \( m \ge n \).

The results are returned in output matrices \( Q \) and \(R \), that have the same type and kind as \( A \).
Given `k = min(m,n)`, one can write \( A = \( Q_1 Q_2 \) \cdot \( \frac{R_1}{0}\) \).
Because the lower rows of \( R \) are zeros, a reduced problem \( A = Q_1 R_1 \) may be solved. The size of
the input arguments determines what problem is solved: on full matrices (`shape(Q)==[m,m]`, `shape(R)==[m,n]`),
the full problem is solved. On reduced matrices (`shape(Q)==[m,k]`, `shape(R)==[k,n]`), the reduced problem is solved.

### Syntax

`call ` [[stdlib_linalg(module):qr(interface)]] `(a, q, r, [, storage] [, overwrite_a] [, err])`

### Arguments

`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.

`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.

`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.

`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.

`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.

`err` (optional): Shall be a `type(linalg_state_type)` value. It is an `intent(out)` argument.

### Return value

Returns the QR factorization matrices into the \( Q \) and \( R \) arguments.

Raises `LINALG_VALUE_ERROR` if any of the matrices has invalid or unsuitable size for the full/reduced problem.
Raises `LINALG_ERROR` on insufficient user storage space.
If the state argument `err` is not present, exceptions trigger an `error stop`.

### Example

```fortran
{!example/linalg/example_qr.f90!}
```

## `qr_space` - Compute internal working space requirements for the QR factorization.

### Status

Experimental

### Description

This subroutine computes the internal working space requirements for the QR factorization, [[stdlib_linalg(module):qr(interface)]] .

### Syntax

`call ` [[stdlib_linalg(module):qr_space(interface)]] `(a, lwork, [, err])`

### Arguments

`a`: Shall be a rank-2 `real` or `complex` array containing the coefficient matrix. It is an `intent(in)` argument.

`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`.

`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument.

### Example

```fortran
{!example/linalg/example_qr_space.f90!}
```

## `eig` - Eigenvalues and Eigenvectors of a Square Matrix

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


### Example

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

### Class

Subroutine

### Arguments
2 changes: 2 additions & 0 deletions example/linalg/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -35,5 +35,7 @@ ADD_EXAMPLE(svd)
ADD_EXAMPLE(svdvals)
ADD_EXAMPLE(determinant)
ADD_EXAMPLE(determinant2)
ADD_EXAMPLE(qr)
ADD_EXAMPLE(qr_space)
ADD_EXAMPLE(cholesky)
ADD_EXAMPLE(chol)
15 changes: 15 additions & 0 deletions example/linalg/example_qr.f90
Original file line number Diff line number Diff line change
@@ -0,0 +1,15 @@
program example_qr
use stdlib_linalg, only: qr
implicit none(type,external)
real :: A(104, 32), Q(104,32), R(32,32)

! Create a random matrix
call random_number(A)

! Compute its QR factorization (reduced)
call qr(A,Q,R)

! Test factorization: Q*R = A
print *, maxval(abs(matmul(Q,R)-A))

end program example_qr
25 changes: 25 additions & 0 deletions example/linalg/example_qr_space.f90
Original file line number Diff line number Diff line change
@@ -0,0 +1,25 @@
! QR example with pre-allocated storage
program example_qr_space
use stdlib_linalg_constants, only: ilp
use stdlib_linalg, only: qr, qr_space, linalg_state_type
implicit none(type,external)
real :: A(104, 32), Q(104,32), R(32,32)
real, allocatable :: work(:)
integer(ilp) :: lwork
type(linalg_state_type) :: err

! Create a random matrix
call random_number(A)

! Prepare QR workspace
call qr_space(A,lwork)
allocate(work(lwork))

! Compute its QR factorization (reduced)
call qr(A,Q,R,storage=work,err=err)

! Test factorization: Q*R = A
print *, maxval(abs(matmul(Q,R)-A))
print *, err%print()

end program example_qr_space
1 change: 1 addition & 0 deletions src/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -30,6 +30,7 @@ set(fppFiles
stdlib_linalg_eigenvalues.fypp
stdlib_linalg_solve.fypp
stdlib_linalg_determinant.fypp
stdlib_linalg_qr.fypp
stdlib_linalg_inverse.fypp
stdlib_linalg_state.fypp
stdlib_linalg_svd.fypp
71 changes: 71 additions & 0 deletions src/stdlib_linalg.fypp
Original file line number Diff line number Diff line change
@@ -40,6 +40,8 @@ module stdlib_linalg
public :: outer_product
public :: kronecker_product
public :: cross_product
public :: qr
public :: qr_space
public :: is_square
public :: is_diagonal
public :: is_symmetric
@@ -526,6 +528,75 @@ module stdlib_linalg
#:endfor
end interface lstsq_space

! QR factorization of rank-2 array A
interface qr
!! version: experimental
!!
!! Computes the QR factorization of matrix \( A = Q R \).
!! ([Specification](../page/specs/stdlib_linalg.html#qr-compute-the-qr-factorization-of-a-matrix))
!!
!!### Summary
!! Compute the QR factorization of a `real` or `complex` matrix: \( A = Q R \), where \( Q \) is orthonormal
!! and \( R \) is upper-triangular. Matrix \( A \) has size `[m,n]`, with \( m\ge n \).
!!
!!### Description
!!
!! This interface provides methods for computing the QR factorization of a matrix.
!! Supported data types include `real` and `complex`. If a pre-allocated work space
!! is provided, no internal memory allocations take place when using this interface.
!!
!! Given `k = min(m,n)`, one can write \( A = \( Q_1 Q_2 \) \cdot \( \frac{R_1}{0}\) \).
!! The user may want the full problem (provide `shape(Q)==[m,m]`, `shape(R)==[m,n]`) or the reduced
!! problem only: \( A = Q_1 R_1 \) (provide `shape(Q)==[m,k]`, `shape(R)==[k,n]`).
!!
!!@note The solution is based on LAPACK's QR factorization (`*GEQRF`) and ordered matrix output (`*ORGQR`, `*UNGQR`).
!!
#:for rk,rt,ri in RC_KINDS_TYPES
pure module subroutine stdlib_linalg_${ri}$_qr(a,q,r,overwrite_a,storage,err)
!> Input matrix a[m,n]
${rt}$, intent(inout), target :: a(:,:)
!> Orthogonal matrix Q ([m,m], or [m,k] if reduced)
${rt}$, intent(out), contiguous, target :: q(:,:)
!> Upper triangular matrix R ([m,n], or [k,n] if reduced)
${rt}$, intent(out), contiguous, target :: r(:,:)
!> [optional] Can A data be overwritten and destroyed?
logical(lk), optional, intent(in) :: overwrite_a
!> [optional] Provide pre-allocated workspace, size to be checked with qr_space
${rt}$, intent(out), optional, target :: storage(:)
!> [optional] state return flag. On error if not requested, the code will stop
type(linalg_state_type), optional, intent(out) :: err
end subroutine stdlib_linalg_${ri}$_qr
#:endfor
end interface qr

! Return the working array space required by the QR factorization solver
interface qr_space
!! version: experimental
!!
!! Computes the working array space required by the QR factorization solver
!! ([Specification](../page/specs/stdlib_linalg.html#qr-space-compute-internal-working-space-requirements-for-the-qr-factorization))
!!
!!### Description
!!
!! This interface returns the size of the `real` or `complex` working storage required by the
!! QR factorization solver. The working size only depends on the kind (`real` or `complex`) and size of
!! the matrix being factorized. Storage size can be used to pre-allocate a working array in case several
!! repeated QR factorizations to a same-size matrix are sought. If pre-allocated working arrays
!! are provided, no internal allocations will take place during the factorization.
!!
#:for rk,rt,ri in RC_KINDS_TYPES
pure module subroutine get_qr_${ri}$_workspace(a,lwork,err)
!> Input matrix a[m,n]
${rt}$, intent(in), target :: a(:,:)
!> Minimum workspace size for both operations
integer(ilp), intent(out) :: lwork
!> State return flag. Returns an error if the query failed
type(linalg_state_type), optional, intent(out) :: err
end subroutine get_qr_${ri}$_workspace
#:endfor
end interface qr_space


interface det
!! version: experimental
!!
267 changes: 267 additions & 0 deletions src/stdlib_linalg_qr.fypp
Original file line number Diff line number Diff line change
@@ -0,0 +1,267 @@
#:include "common.fypp"
#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
submodule (stdlib_linalg) stdlib_linalg_qr
use stdlib_linalg_constants
use stdlib_linalg_lapack, only: geqrf, orgqr, ungqr
use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, &
LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR
implicit none(type,external)

character(*), parameter :: this = 'qr'

contains

! Check problem size and evaluate whether full/reduced problem is requested
pure subroutine check_problem_size(m,n,q1,q2,r1,r2,err,reduced)
integer(ilp), intent(in) :: m,n,q1,q2,r1,r2
type(linalg_state_type), intent(out) :: err
logical, intent(out) :: reduced

integer(ilp) :: k

k = min(m,n)
reduced = .false.

! Check sizes
if (m<1 .or. n<1) then
err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size: a(m,n)=',[m,n])
else

! Check if we should operate on reduced full QR
! - Reduced: shape(Q)==[m,k], shape(R)==[k,n]
! - Full : shape(Q)==[m,m], shape(R)==[m,n]
if (all([q1,q2]==[m,k] .and. [r1,r2]==[k,n])) then
reduced = .true.
elseif (all([q1,q2]==[m,m] .and. [r1,r2]==[m,n])) then
reduced = .false.
else
err = linalg_state_type(this,LINALG_VALUE_ERROR,'with a=',[m,n],'q=',[q1,q2],'r=',[r1,r2], &
'problem is neither full (q=',[m,m],'r=',[m,n],') nor reduced (q=',[m,m],'r=',[m,n],')')
endif
end if

end subroutine check_problem_size

elemental subroutine handle_orgqr_info(info,m,n,k,lwork,err)
integer(ilp), intent(in) :: info,m,n,k,lwork
type(linalg_state_type), intent(out) :: err

! Process output
select case (info)
case (0)
! Success
case (-1)
err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size m=',m)
case (-2)
err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size n=',n)
case (-4)
err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid k=min(m,n)=',k)
case (-5)
err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size a=',[m,n])
case (-8)
err = linalg_state_type(this,LINALG_ERROR,'invalid input for lwork=',lwork)
case default
err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error')
end select

end subroutine handle_orgqr_info

elemental subroutine handle_geqrf_info(info,m,n,lwork,err)
integer(ilp), intent(in) :: info,m,n,lwork
type(linalg_state_type), intent(out) :: err

! Process output
select case (info)
case (0)
! Success
case (-1)
err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size m=',m)
case (-2)
err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size n=',n)
case (-4)
err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size a=',[m,n])
case (-7)
err = linalg_state_type(this,LINALG_ERROR,'invalid input for lwork=',lwork)
case default
err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error')
end select

end subroutine handle_geqrf_info

#:for rk,rt,ri in RC_KINDS_TYPES

! Get workspace size for QR operations
pure module subroutine get_qr_${ri}$_workspace(a,lwork,err)
!> Input matrix a[m,n]
${rt}$, intent(in), target :: a(:,:)
!> Minimum workspace size for both operations
integer(ilp), intent(out) :: lwork
!> State return flag. Returns an error if the query failed
type(linalg_state_type), optional, intent(out) :: err

integer(ilp) :: m,n,k,info,lwork_qr,lwork_ord
${rt}$ :: work_dummy(1),tau_dummy(1),a_dummy(1,1)
type(linalg_state_type) :: err0

lwork = -1_ilp

!> Problem sizes
m = size(a,1,kind=ilp)
n = size(a,2,kind=ilp)
k = min(m,n)

! QR space
lwork_qr = -1_ilp
call geqrf(m,n,a_dummy,m,tau_dummy,work_dummy,lwork_qr,info)
call handle_geqrf_info(info,m,n,lwork_qr,err0)
if (err0%error()) then
call linalg_error_handling(err0,err)
return
endif
lwork_qr = ceiling(real(work_dummy(1),kind=${rk}$),kind=ilp)

! Ordering space (for full factorization)
lwork_ord = -1_ilp
call #{if rt.startswith('complex')}# ungqr #{else}# orgqr #{endif}# &
(m,m,k,a_dummy,m,tau_dummy,work_dummy,lwork_ord,info)
call handle_orgqr_info(info,m,n,k,lwork_ord,err0)
if (err0%error()) then
call linalg_error_handling(err0,err)
return
endif
lwork_ord = ceiling(real(work_dummy(1),kind=${rk}$),kind=ilp)

! Pick the largest size, so two operations can be performed with the same allocation
lwork = max(lwork_qr, lwork_ord)

end subroutine get_qr_${ri}$_workspace

! Compute the solution to a real system of linear equations A * X = B
pure module subroutine stdlib_linalg_${ri}$_qr(a,q,r,overwrite_a,storage,err)
!> Input matrix a[m,n]
${rt}$, intent(inout), target :: a(:,:)
!> Orthogonal matrix Q ([m,m], or [m,k] if reduced)
${rt}$, intent(out), contiguous, target :: q(:,:)
!> Upper triangular matrix R ([m,n], or [k,n] if reduced)
${rt}$, intent(out), contiguous, target :: r(:,:)
!> [optional] Can A data be overwritten and destroyed?
logical(lk), optional, intent(in) :: overwrite_a
!> [optional] Provide pre-allocated workspace, size to be checked with qr_space
${rt}$, intent(out), optional, target :: storage(:)
!> [optional] state return flag. On error if not requested, the code will stop
type(linalg_state_type), optional, intent(out) :: err

!> Local variables
type(linalg_state_type) :: err0
integer(ilp) :: i,j,m,n,k,q1,q2,r1,r2,lda,lwork,info
logical(lk) :: overwrite_a_,use_q_matrix,reduced
${rt}$ :: r11
${rt}$, parameter :: zero = 0.0_${rk}$

${rt}$, pointer :: amat(:,:),tau(:),work(:)

!> Problem sizes
m = size(a,1,kind=ilp)
n = size(a,2,kind=ilp)
k = min(m,n)
q1 = size(q,1,kind=ilp)
q2 = size(q,2,kind=ilp)
r1 = size(r,1,kind=ilp)
r2 = size(r,2,kind=ilp)

! Check if we should operate on reduced full QR
call check_problem_size(m,n,q1,q2,r1,r2,err0,reduced)
if (err0%error()) then
call linalg_error_handling(err0,err)
return
end if

! Check if Q can be used as temporary storage for A,
! to be destroyed by *GEQRF
use_q_matrix = q1>=m .and. q2>=n

! Can A be overwritten? By default, do not overwrite
overwrite_a_ = .false._lk
if (present(overwrite_a) .and. .not.use_q_matrix) overwrite_a_ = overwrite_a

! Initialize a matrix temporary, or reuse available
! storage if possible
if (use_q_matrix) then
amat => q
q(:m,:n) = a
elseif (overwrite_a_) then
amat => a
else
allocate(amat(m,n),source=a)
endif
lda = size(amat,1,kind=ilp)

! To store the elementary reflectors, we need a [1:k] column.
if (.not.use_q_matrix) then
! Q is not being used as the storage matrix
tau(1:k) => q(1:k,1)
else
! R has unused contiguous storage in the 1st column, except for the
! diagonal element. So, use the full column and store it in a dummy variable
tau(1:k) => r(1:k,1)
endif

! Retrieve workspace size
call get_qr_${ri}$_workspace(a,lwork,err0)

if (err0%ok()) then

if (present(storage)) then
work => storage
else
allocate(work(lwork))
endif
if (.not.size(work,kind=ilp)>=lwork) then
err0 = linalg_state_type(this,LINALG_ERROR,'insufficient workspace: should be at least ',lwork)
call linalg_error_handling(err0,err)
return
endif

! Compute factorization.
call geqrf(m,n,amat,m,tau,work,lwork,info)
call handle_geqrf_info(info,m,n,lwork,err0)

if (err0%ok()) then

! Get R matrix out before overwritten.
! Do not copy the first column at this stage: it may be being used by `tau`
r11 = amat(1,1)
forall(i=1:min(r1,m),j=2:n) r(i,j) = merge(amat(i,j),zero,i<=j)

! Convert K elementary reflectors tau(1:k) -> orthogonal matrix Q
call #{if rt.startswith('complex')}# ungqr #{else}# orgqr #{endif}# &
(q1,q2,k,amat,lda,tau,work,lwork,info)
call handle_orgqr_info(info,m,n,k,lwork,err0)

! Copy result back to Q
if (.not.use_q_matrix) q = amat(:q1,:q2)

! Copy first column of R
r(1,1) = r11
r(2:r1,1) = zero

! Ensure last m-n rows of R are zeros,
! if full matrices were provided
if (.not.reduced) r(k+1:m,1:n) = zero

endif

if (.not.present(storage)) deallocate(work)

endif

if (.not.(use_q_matrix.or.overwrite_a_)) deallocate(amat)

! Process output and return
call linalg_error_handling(err0,err)

end subroutine stdlib_linalg_${ri}$_qr

#:endfor

end submodule stdlib_linalg_qr
2 changes: 2 additions & 0 deletions test/linalg/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -8,6 +8,7 @@ set(
"test_linalg_inverse.fypp"
"test_linalg_lstsq.fypp"
"test_linalg_determinant.fypp"
"test_linalg_qr.fypp"
"test_linalg_svd.fypp"
"test_linalg_matrix_property_checks.fypp"
)
@@ -21,5 +22,6 @@ ADDTEST(linalg_matrix_property_checks)
ADDTEST(linalg_inverse)
ADDTEST(linalg_solve)
ADDTEST(linalg_lstsq)
ADDTEST(linalg_qr)
ADDTEST(linalg_svd)
ADDTEST(blas_lapack)
140 changes: 140 additions & 0 deletions test/linalg/test_linalg_qr.fypp
Original file line number Diff line number Diff line change
@@ -0,0 +1,140 @@
#:include "common.fypp"
#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
! Test QR factorization
module test_linalg_qr
use testdrive, only: error_type, check, new_unittest, unittest_type
use stdlib_linalg_constants
use stdlib_linalg_state, only: LINALG_VALUE_ERROR,linalg_state_type
use stdlib_linalg, only: qr,qr_space
use ieee_arithmetic, only: ieee_value,ieee_quiet_nan

implicit none (type,external)

public :: test_qr_factorization

contains

!> QR factorization tests
subroutine test_qr_factorization(tests)
!> Collection of tests
type(unittest_type), allocatable, intent(out) :: tests(:)

allocate(tests(0))

#:for rk,rt,ri in RC_KINDS_TYPES
tests = [tests,new_unittest("qr_random_${ri}$",test_qr_random_${ri}$)]
#:endfor

end subroutine test_qr_factorization

!> QR factorization of a random matrix
#:for rk,rt,ri in RC_KINDS_TYPES
subroutine test_qr_random_${ri}$(error)
type(error_type), allocatable, intent(out) :: error

integer(ilp), parameter :: m = 15_ilp
integer(ilp), parameter :: n = 4_ilp
integer(ilp), parameter :: k = min(m,n)
real(${rk}$), parameter :: tol = 100*sqrt(epsilon(0.0_${rk}$))
${rt}$ :: a(m,n),aorig(m,n),q(m,m),r(m,n),qred(m,k),rred(k,n),qerr(m,6),rerr(6,n)
real(${rk}$) :: rea(m,n),ima(m,n)
integer(ilp) :: lwork
${rt}$, allocatable :: work(:)
type(linalg_state_type) :: state

call random_number(rea)
#:if rt.startswith('complex')
call random_number(ima)
a = cmplx(rea,ima,kind=${rk}$)
#:else
a = rea
#:endif
aorig = a

! 1) QR factorization with full matrices. Input NaNs to be sure Q and R are OK on return
q = ieee_value(0.0_${rk}$,ieee_quiet_nan)
r = ieee_value(0.0_${rk}$,ieee_quiet_nan)
call qr(a,q,r,err=state)

! Check return code
call check(error,state%ok(),state%print())
if (allocated(error)) return

! Check solution
call check(error, all(abs(a-matmul(q,r))<tol), 'converged solution (fulle)')
if (allocated(error)) return

! 2) QR factorization with reduced matrices
call qr(a,qred,rred,err=state)

! Check return code
call check(error,state%ok(),state%print())
if (allocated(error)) return

! Check solution
call check(error, all(abs(a-matmul(qred,rred))<tol), 'converged solution (reduced)')
if (allocated(error)) return

! 3) overwrite A
call qr(a,qred,rred,overwrite_a=.true.,err=state)

! Check return code
call check(error,state%ok(),state%print())
if (allocated(error)) return

! Check solution
call check(error, all(abs(aorig-matmul(qred,rred))<tol), 'converged solution (overwrite A)')
if (allocated(error)) return

! 4) External storage option
a = aorig
call qr_space(a,lwork)
allocate(work(lwork))
call qr(a,q,r,storage=work,err=state)

! Check return code
call check(error,state%ok(),state%print())
if (allocated(error)) return

! Check solution
call check(error, all(abs(a-matmul(q,r))<tol), 'converged solution (external storage)')
if (allocated(error)) return

! Check that an invalid problem size returns an error
a = aorig
call qr(a,qerr,rerr,err=state)
call check(error,state%error(),'invalid matrix sizes')
if (allocated(error)) return

end subroutine test_qr_random_${ri}$

#:endfor


end module test_linalg_qr

program test_qr
use, intrinsic :: iso_fortran_env, only : error_unit
use testdrive, only : run_testsuite, new_testsuite, testsuite_type
use test_linalg_qr, only : test_qr_factorization
implicit none
integer :: stat, is
type(testsuite_type), allocatable :: testsuites(:)
character(len=*), parameter :: fmt = '("#", *(1x, a))'

stat = 0

testsuites = [ &
new_testsuite("linalg_qr", test_qr_factorization) &
]

do is = 1, size(testsuites)
write(error_unit, fmt) "Testing:", testsuites(is)%name
call run_testsuite(testsuites(is)%collect, error_unit, stat)
end do

if (stat > 0) then
write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!"
error stop
end if
end program test_qr