diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index b03fdd192..47ca5102e 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -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 diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index c8c676043..487e2197b 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -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) diff --git a/example/linalg/example_qr.f90 b/example/linalg/example_qr.f90 new file mode 100644 index 000000000..4e8d19e3e --- /dev/null +++ b/example/linalg/example_qr.f90 @@ -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 diff --git a/example/linalg/example_qr_space.f90 b/example/linalg/example_qr_space.f90 new file mode 100644 index 000000000..22235ca12 --- /dev/null +++ b/example/linalg/example_qr_space.f90 @@ -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 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index ef11b642e..82cb2c450 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -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 diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 4574568e5..cd4de131c 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -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 !! diff --git a/src/stdlib_linalg_qr.fypp b/src/stdlib_linalg_qr.fypp new file mode 100644 index 000000000..06d772d63 --- /dev/null +++ b/src/stdlib_linalg_qr.fypp @@ -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 diff --git a/test/linalg/CMakeLists.txt b/test/linalg/CMakeLists.txt index 0cdbbaa3c..adf8b5db2 100644 --- a/test/linalg/CMakeLists.txt +++ b/test/linalg/CMakeLists.txt @@ -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) diff --git a/test/linalg/test_linalg_qr.fypp b/test/linalg/test_linalg_qr.fypp new file mode 100644 index 000000000..8db1d32ad --- /dev/null +++ b/test/linalg/test_linalg_qr.fypp @@ -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)) 0) then + write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!" + error stop + end if +end program test_qr