From 6cac0cc36366401b3b4a68544ce1918e42108907 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 3 Jun 2024 09:43:33 +0200 Subject: [PATCH 01/30] add `qr` module --- src/CMakeLists.txt | 1 + src/stdlib_linalg_qr.fypp | 269 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 270 insertions(+) create mode 100644 src/stdlib_linalg_qr.fypp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index ced6e7b43..209883c5b 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -29,6 +29,7 @@ set(fppFiles stdlib_linalg_cross_product.fypp stdlib_linalg_solve.fypp stdlib_linalg_determinant.fypp + stdlib_linalg_qr.fypp stdlib_linalg_state.fypp stdlib_optval.fypp stdlib_selection.fypp diff --git a/src/stdlib_linalg_qr.fypp b/src/stdlib_linalg_qr.fypp new file mode 100644 index 000000000..966aea976 --- /dev/null +++ b/src/stdlib_linalg_qr.fypp @@ -0,0 +1,269 @@ +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +module 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) + private + + !> QR factorization of a matrix + public :: qr + public :: qr_space + + ! Scipy: solve(a, b, lower=False, overwrite_a=False, overwrite_b=False, check_finite=True, assume_a='gen', transposed=False)[source]# + ! IMSL: lu_solve(a, b, transpose=False) + + interface qr + #:for rk,rt,ri in RC_KINDS_TYPES + module procedure stdlib_linalg_${ri}$_qr + #:endfor + end interface qr + + interface qr_space + #:for rk,rt,ri in RC_KINDS_TYPES + module procedure get_qr_${ri}$_workspace + #:endfor + end interface qr_space + + character(*), parameter :: this = 'qr' + + + contains + + 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 subroutine get_qr_${ri}$_workspace(a,lwork,err) + !> Input matrix a[m,n] + ${rt}$, intent(inout), 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) + 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,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 + lwork_ord = -1_ilp + call #{if rt.startswith('complex')}# ungqr #{else}# orgqr #{endif}# & + (m,n,k,a,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 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(inout), 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 sizes + if (m<1 .or. n<1 .or. q1 orthogonal matrix Q + call #{if rt.startswith('complex')}# ungqr #{else}# orgqr #{endif}# & + (m,n,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 module stdlib_linalg_qr From d4b5dbc81c41aec2f9dadc12dd3c1dbea6cf7d72 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 3 Jun 2024 09:46:00 +0200 Subject: [PATCH 02/30] add examples --- example/linalg/CMakeLists.txt | 2 ++ example/linalg/example_qr.f90 | 15 +++++++++++++++ example/linalg/example_qr_space.f90 | 24 ++++++++++++++++++++++++ 3 files changed, 41 insertions(+) create mode 100644 example/linalg/example_qr.f90 create mode 100644 example/linalg/example_qr_space.f90 diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index 2d80f067c..74d554ff4 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -25,3 +25,5 @@ ADD_EXAMPLE(solve2) ADD_EXAMPLE(solve3) ADD_EXAMPLE(determinant) ADD_EXAMPLE(determinant2) +ADD_EXAMPLE(qr) +ADD_EXAMPLE(qr_space) diff --git a/example/linalg/example_qr.f90 b/example/linalg/example_qr.f90 new file mode 100644 index 000000000..df24c7a84 --- /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..c74f71310 --- /dev/null +++ b/example/linalg/example_qr_space.f90 @@ -0,0 +1,24 @@ +! QR example with pre-allocated storage +program example_qr_space + use stdlib_linalg, only: ilp, 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) :: 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 From fe0e243794d80f77e34604707b83dac2d881c177 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 3 Jun 2024 09:53:19 +0200 Subject: [PATCH 03/30] submodule --- src/stdlib_linalg.fypp | 41 +++++++++++++++++++++++++++++++++++++++ src/stdlib_linalg_qr.fypp | 25 ++---------------------- 2 files changed, 43 insertions(+), 23 deletions(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 134bf7af5..c5b4389a4 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -29,6 +29,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 @@ -450,6 +452,45 @@ module stdlib_linalg #:endfor end interface lstsq_space + ! QR factorization of rank-2 array A + interface qr + #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" + 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(inout), 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 + #:endif + #:endfor + end interface qr + + ! Return the working array space required by the QR factorization solver + interface qr_space + #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" + pure module subroutine get_qr_${ri}$_workspace(a,lwork,err) + !> Input matrix a[m,n] + ${rt}$, intent(inout), 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 + #:endif + #:endfor + end interface qr_space + + interface det !! version: experimental !! diff --git a/src/stdlib_linalg_qr.fypp b/src/stdlib_linalg_qr.fypp index 966aea976..beb4f7ef6 100644 --- a/src/stdlib_linalg_qr.fypp +++ b/src/stdlib_linalg_qr.fypp @@ -1,35 +1,14 @@ #:include "common.fypp" #:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES -module stdlib_linalg_qr +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) - private - !> QR factorization of a matrix - public :: qr - public :: qr_space - - ! Scipy: solve(a, b, lower=False, overwrite_a=False, overwrite_b=False, check_finite=True, assume_a='gen', transposed=False)[source]# - ! IMSL: lu_solve(a, b, transpose=False) - - interface qr - #:for rk,rt,ri in RC_KINDS_TYPES - module procedure stdlib_linalg_${ri}$_qr - #:endfor - end interface qr - - interface qr_space - #:for rk,rt,ri in RC_KINDS_TYPES - module procedure get_qr_${ri}$_workspace - #:endfor - end interface qr_space - character(*), parameter :: this = 'qr' - contains elemental subroutine handle_orgqr_info(info,m,n,k,lwork,err) @@ -266,4 +245,4 @@ module stdlib_linalg_qr #:endfor -end module stdlib_linalg_qr +end submodule stdlib_linalg_qr From ec6f68754c0c7838a578c618fc45c5910e3fe265 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 3 Jun 2024 09:54:35 +0200 Subject: [PATCH 04/30] fix the examples --- example/linalg/example_qr.f90 | 2 +- example/linalg/example_qr_space.f90 | 5 +++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/example/linalg/example_qr.f90 b/example/linalg/example_qr.f90 index df24c7a84..4e8d19e3e 100644 --- a/example/linalg/example_qr.f90 +++ b/example/linalg/example_qr.f90 @@ -1,4 +1,4 @@ -program example_qr +program example_qr use stdlib_linalg, only: qr implicit none(type,external) real :: A(104, 32), Q(104,32), R(32,32) diff --git a/example/linalg/example_qr_space.f90 b/example/linalg/example_qr_space.f90 index c74f71310..22235ca12 100644 --- a/example/linalg/example_qr_space.f90 +++ b/example/linalg/example_qr_space.f90 @@ -1,11 +1,12 @@ ! QR example with pre-allocated storage program example_qr_space - use stdlib_linalg, only: ilp, qr, qr_space, linalg_state_type + 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) :: err + type(linalg_state_type) :: err ! Create a random matrix call random_number(A) From a6a6838bfa911065aa3c06a57ef21cf57c82d62d Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 3 Jun 2024 10:17:52 +0200 Subject: [PATCH 05/30] exclude `xdp` --- src/stdlib_linalg_qr.fypp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/stdlib_linalg_qr.fypp b/src/stdlib_linalg_qr.fypp index beb4f7ef6..697e86041 100644 --- a/src/stdlib_linalg_qr.fypp +++ b/src/stdlib_linalg_qr.fypp @@ -58,6 +58,7 @@ submodule (stdlib_linalg) stdlib_linalg_qr end subroutine handle_geqrf_info #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" ! Get workspace size for QR operations pure subroutine get_qr_${ri}$_workspace(a,lwork,err) @@ -243,6 +244,7 @@ submodule (stdlib_linalg) stdlib_linalg_qr end subroutine stdlib_linalg_${ri}$_qr + #:endif #:endfor end submodule stdlib_linalg_qr From feecc14a8d520d2b746403f100198af749c8b521 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 3 Jun 2024 10:29:20 +0200 Subject: [PATCH 06/30] fix `module` `subroutine` --- src/stdlib_linalg_qr.fypp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/stdlib_linalg_qr.fypp b/src/stdlib_linalg_qr.fypp index 697e86041..7b0a0bf1d 100644 --- a/src/stdlib_linalg_qr.fypp +++ b/src/stdlib_linalg_qr.fypp @@ -61,7 +61,7 @@ submodule (stdlib_linalg) stdlib_linalg_qr #:if rk!="xdp" ! Get workspace size for QR operations - pure subroutine get_qr_${ri}$_workspace(a,lwork,err) + pure module subroutine get_qr_${ri}$_workspace(a,lwork,err) !> Input matrix a[m,n] ${rt}$, intent(inout), target :: a(:,:) !> Minimum workspace size for both operations @@ -107,7 +107,7 @@ submodule (stdlib_linalg) stdlib_linalg_qr end subroutine get_qr_${ri}$_workspace ! Compute the solution to a real system of linear equations A * X = B - pure subroutine stdlib_linalg_${ri}$_qr(a,q,r,overwrite_a,storage,err) + 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) From f3d74fff3ffcf7b2803385634cf3a7ff2cc33ba9 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 5 Jun 2024 09:59:22 +0200 Subject: [PATCH 07/30] enforce full/reduced problem sizes --- src/stdlib_linalg_qr.fypp | 46 +++++++++++++++++++++++++++++---------- 1 file changed, 35 insertions(+), 11 deletions(-) diff --git a/src/stdlib_linalg_qr.fypp b/src/stdlib_linalg_qr.fypp index 7b0a0bf1d..b54c3ec14 100644 --- a/src/stdlib_linalg_qr.fypp +++ b/src/stdlib_linalg_qr.fypp @@ -11,6 +11,37 @@ submodule (stdlib_linalg) stdlib_linalg_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 @@ -139,20 +170,13 @@ submodule (stdlib_linalg) stdlib_linalg_qr r1 = size(r,1,kind=ilp) r2 = size(r,2,kind=ilp) - ! Check sizes - if (m<1 .or. n<1 .or. q1 0) then + write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!" + error stop + end if +end program test_qr From 7c9d5b48fe3cb6d70731510daeb2a174dd815706 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 5 Jun 2024 10:20:44 +0200 Subject: [PATCH 09/30] remove `xdp` from tests --- test/linalg/test_linalg_qr.fypp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/linalg/test_linalg_qr.fypp b/test/linalg/test_linalg_qr.fypp index e023889e6..e35661d81 100644 --- a/test/linalg/test_linalg_qr.fypp +++ b/test/linalg/test_linalg_qr.fypp @@ -29,6 +29,7 @@ module test_linalg_qr !> QR factorization of a random matrix #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" subroutine test_qr_random_${ri}$(error) type(error_type), allocatable, intent(out) :: error @@ -107,6 +108,7 @@ module test_linalg_qr end subroutine test_qr_random_${ri}$ + #:endif #:endfor From cda4b0d12774fdb32d741eca51facb0617faebff Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 5 Jun 2024 10:33:17 +0200 Subject: [PATCH 10/30] document `qr_space` --- doc/specs/stdlib_linalg.md | 28 ++++++++++++++++++++++++++++ src/stdlib_linalg.fypp | 13 +++++++++++++ 2 files changed, 41 insertions(+) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index ffc4f7c42..ea869fe7c 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -898,3 +898,31 @@ Exceptions trigger an `error stop`. ```fortran {!example/linalg/example_determinant2.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!} +``` diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index c5b4389a4..0780dc35c 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -476,6 +476,19 @@ module stdlib_linalg ! 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 #:if rk!="xdp" pure module subroutine get_qr_${ri}$_workspace(a,lwork,err) From 8df204f108589bb4d00eaf4b976773cc9d692339 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 5 Jun 2024 10:48:27 +0200 Subject: [PATCH 11/30] document `qr` --- doc/specs/stdlib_linalg.md | 48 ++++++++++++++++++++++++++++++++++++++ src/stdlib_linalg.fypp | 22 +++++++++++++++++ 2 files changed, 70 insertions(+) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index ea869fe7c..9a89fe0c7 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -899,6 +899,54 @@ 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, the user may require to return full matrices (provide `shape(Q)==[m,m]`, `shape(R)==[m,n]`) +or reduced matrices: \( A = Q_1 R_1 \) (provide `shape(Q)==[m,k]`, `shape(R)==[k,n]`). + +### 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 normally an `intent(in)` argument. If `overwrite_a=.true.`, 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 shape 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 shape 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. It minimum size can be determined with a call to [[stdlib_linalg(module):qr_space(interface)]]. It is an `intent(inout)` 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. This is an `intent(in)` argument. + +`err` (optional): Shall be a `type(linalg_state_type)` value. This 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 diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 0780dc35c..25f6a3484 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -454,6 +454,28 @@ module stdlib_linalg ! 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`). + !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``). + !! #:for rk,rt,ri in RC_KINDS_TYPES #:if rk!="xdp" pure module subroutine stdlib_linalg_${ri}$_qr(a,q,r,overwrite_a,storage,err) From b4889c5141a93bc9bc800a2b21872aedac2fd39c Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 5 Jun 2024 10:52:36 +0200 Subject: [PATCH 12/30] relax tol --- test/linalg/test_linalg_qr.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/linalg/test_linalg_qr.fypp b/test/linalg/test_linalg_qr.fypp index e35661d81..65de4873e 100644 --- a/test/linalg/test_linalg_qr.fypp +++ b/test/linalg/test_linalg_qr.fypp @@ -36,7 +36,7 @@ module test_linalg_qr integer(ilp), parameter :: m = 15_ilp integer(ilp), parameter :: n = 4_ilp integer(ilp), parameter :: k = min(m,n) - real(${rk}$), parameter :: tol = 10*sqrt(epsilon(0.0_${rk}$)) + 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 From 9ac97d1ae2973349f6d3a25335b1ca0b3b26c1b1 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 5 Jun 2024 11:00:35 +0200 Subject: [PATCH 13/30] deactivate full tol check --- test/linalg/test_linalg_qr.fypp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/test/linalg/test_linalg_qr.fypp b/test/linalg/test_linalg_qr.fypp index 65de4873e..5c46be831 100644 --- a/test/linalg/test_linalg_qr.fypp +++ b/test/linalg/test_linalg_qr.fypp @@ -60,9 +60,8 @@ module test_linalg_qr if (allocated(error)) return ! Check solution - ! Check solution - call check(error, all(abs(a-matmul(q,r)) Date: Wed, 5 Jun 2024 11:18:00 +0200 Subject: [PATCH 14/30] export error amount --- test/linalg/test_linalg_qr.fypp | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/test/linalg/test_linalg_qr.fypp b/test/linalg/test_linalg_qr.fypp index 5c46be831..40120daf2 100644 --- a/test/linalg/test_linalg_qr.fypp +++ b/test/linalg/test_linalg_qr.fypp @@ -4,7 +4,8 @@ module test_linalg_qr use testdrive, only: error_type, check, new_unittest, unittest_type use stdlib_linalg_constants - use stdlib_linalg, only: qr,qr_space,linalg_state_type + use stdlib_linalg_state, only: LINALG_VALUE_ERROR,linalg_state_type + use stdlib_linalg, only: qr,qr_space implicit none (type,external) @@ -53,6 +54,7 @@ module test_linalg_qr aorig = a ! 1) QR factorization with full matrices + q = 1.0_${rk}$ call qr(a,q,r,err=state) ! Check return code @@ -60,8 +62,13 @@ module test_linalg_qr if (allocated(error)) return ! Check solution -! call check(error, all(abs(a-matmul(q,r)) Date: Wed, 5 Jun 2024 11:39:35 +0200 Subject: [PATCH 15/30] Fix: return full `Q` from `*ORGQR` --- src/stdlib_linalg_qr.fypp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/stdlib_linalg_qr.fypp b/src/stdlib_linalg_qr.fypp index b54c3ec14..806d6fecf 100644 --- a/src/stdlib_linalg_qr.fypp +++ b/src/stdlib_linalg_qr.fypp @@ -121,10 +121,10 @@ submodule (stdlib_linalg) stdlib_linalg_qr endif lwork_qr = ceiling(real(work_dummy(1),kind=${rk}$),kind=ilp) - ! Ordering space + ! Ordering space (for full factorization) lwork_ord = -1_ilp call #{if rt.startswith('complex')}# ungqr #{else}# orgqr #{endif}# & - (m,n,k,a,m,tau_dummy,work_dummy,lwork_ord,info) + (m,m,k,a,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) @@ -241,7 +241,7 @@ submodule (stdlib_linalg) stdlib_linalg_qr ! Convert K elementary reflectors tau(1:k) -> orthogonal matrix Q call #{if rt.startswith('complex')}# ungqr #{else}# orgqr #{endif}# & - (m,n,k,amat,lda,tau,work,lwork,info) + (q1,q2,k,amat,lda,tau,work,lwork,info) call handle_orgqr_info(info,m,n,k,lwork,err0) ! Copy result back to Q From 3db022d47b321ce98a300b75b4ca17a1e24c79bf Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 5 Jun 2024 11:39:47 +0200 Subject: [PATCH 16/30] add test: init NaNs, check correct return --- test/linalg/test_linalg_qr.fypp | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/test/linalg/test_linalg_qr.fypp b/test/linalg/test_linalg_qr.fypp index 40120daf2..d615962a8 100644 --- a/test/linalg/test_linalg_qr.fypp +++ b/test/linalg/test_linalg_qr.fypp @@ -53,8 +53,11 @@ module test_linalg_qr #:endif aorig = a - ! 1) QR factorization with full matrices - q = 1.0_${rk}$ + ! 1) QR factorization with full matrices. Input NaNs to be sure Q and R are OK on return + q = 0.0_${rk}$ + q = 1.0_${rk}$/q + r = 0.0_${rk}$ + r = 1.0_${rk}$/r call qr(a,q,r,err=state) ! Check return code @@ -62,13 +65,8 @@ module test_linalg_qr if (allocated(error)) return ! Check solution - if (.not.all(abs(a-matmul(q,r)) Date: Sat, 6 Jul 2024 15:51:45 +0200 Subject: [PATCH 17/30] lift `xdp` restriction --- src/stdlib_linalg.fypp | 5 ----- src/stdlib_linalg_qr.fypp | 2 -- test/linalg/test_linalg_qr.fypp | 4 ---- 3 files changed, 11 deletions(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 7e848bd74..1b8786c9d 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -466,10 +466,8 @@ module stdlib_linalg !! 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`). - !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``). !! #:for rk,rt,ri in RC_KINDS_TYPES - #:if rk!="xdp" pure module subroutine stdlib_linalg_${ri}$_qr(a,q,r,overwrite_a,storage,err) !> Input matrix a[m,n] ${rt}$, intent(inout), target :: a(:,:) @@ -484,7 +482,6 @@ module stdlib_linalg !> [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 - #:endif #:endfor end interface qr @@ -504,7 +501,6 @@ module stdlib_linalg !! are provided, no internal allocations will take place during the factorization. !! #:for rk,rt,ri in RC_KINDS_TYPES - #:if rk!="xdp" pure module subroutine get_qr_${ri}$_workspace(a,lwork,err) !> Input matrix a[m,n] ${rt}$, intent(inout), target :: a(:,:) @@ -513,7 +509,6 @@ module stdlib_linalg !> State return flag. Returns an error if the query failed type(linalg_state_type), optional, intent(out) :: err end subroutine get_qr_${ri}$_workspace - #:endif #:endfor end interface qr_space diff --git a/src/stdlib_linalg_qr.fypp b/src/stdlib_linalg_qr.fypp index 806d6fecf..095e71c44 100644 --- a/src/stdlib_linalg_qr.fypp +++ b/src/stdlib_linalg_qr.fypp @@ -89,7 +89,6 @@ submodule (stdlib_linalg) stdlib_linalg_qr end subroutine handle_geqrf_info #:for rk,rt,ri in RC_KINDS_TYPES - #:if rk!="xdp" ! Get workspace size for QR operations pure module subroutine get_qr_${ri}$_workspace(a,lwork,err) @@ -268,7 +267,6 @@ submodule (stdlib_linalg) stdlib_linalg_qr end subroutine stdlib_linalg_${ri}$_qr - #:endif #:endfor end submodule stdlib_linalg_qr diff --git a/test/linalg/test_linalg_qr.fypp b/test/linalg/test_linalg_qr.fypp index d615962a8..81d3edb0e 100644 --- a/test/linalg/test_linalg_qr.fypp +++ b/test/linalg/test_linalg_qr.fypp @@ -21,16 +21,13 @@ module test_linalg_qr allocate(tests(0)) #:for rk,rt,ri in RC_KINDS_TYPES - #:if rk!="xdp" tests = [tests,new_unittest("qr_random_${ri}$",test_qr_random_${ri}$)] - #:endif #:endfor end subroutine test_qr_factorization !> QR factorization of a random matrix #:for rk,rt,ri in RC_KINDS_TYPES - #:if rk!="xdp" subroutine test_qr_random_${ri}$(error) type(error_type), allocatable, intent(out) :: error @@ -112,7 +109,6 @@ module test_linalg_qr end subroutine test_qr_random_${ri}$ - #:endif #:endfor From 381d02491a828b92e3ce9c5f8378ddfbca0fe75e Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 8 Jul 2024 08:55:38 +0200 Subject: [PATCH 18/30] Update doc/specs/stdlib_linalg.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 76c879223..8448b71a1 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -901,7 +901,7 @@ Exceptions trigger an `error stop`. {!example/linalg/example_determinant2.f90!} ``` -## `qr` - Compute the QR factorization of a matrix. +## `qr` - Compute the QR factorization of a matrix ### Status From 931f07aec1699565e22a71fe1a288e2e739ca43b Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 8 Jul 2024 08:55:48 +0200 Subject: [PATCH 19/30] Update doc/specs/stdlib_linalg.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 8448b71a1..f746905b5 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -909,7 +909,7 @@ Experimental ### Description -This subroutine computes the QR factorization of a `real` or `complex` matrix: \( A = Q R \). where \( Q \) +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 \). From 5f1591323c9d950b7e6e067c95b5c362faef3bbc Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 8 Jul 2024 08:56:20 +0200 Subject: [PATCH 20/30] Update doc/specs/stdlib_linalg.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index f746905b5..43dcce705 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -923,7 +923,7 @@ or reduced matrices: \( A = Q_1 R_1 \) (provide `shape(Q)==[m,k]`, `shape(R)==[k ### Arguments -`a`: Shall be a rank-2 `real` or `complex` array containing the coefficient matrix of size `[m,n]`. It is normally an `intent(in)` argument. If `overwrite_a=.true.`, it is an `intent(inout)` argument, and is destroyed upon return. +`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 shape either `[m,m]` or `[m,k]` whether the full or the reduced problem is sought for. From 997021c3abdb46af99bae874f7823231d00f9ded Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 8 Jul 2024 08:56:30 +0200 Subject: [PATCH 21/30] Update doc/specs/stdlib_linalg.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 43dcce705..07dc472e9 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -925,7 +925,7 @@ or reduced matrices: \( A = Q_1 R_1 \) (provide `shape(Q)==[m,k]`, `shape(R)==[k `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 shape either `[m,m]` or `[m,k]` whether the full or the reduced problem is sought for. +`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 shape either `[m,n]` or `[k,n]` whether the full or the reduced problem is sought for. From 789e6d0c455e3233d83fedb44705635a5df70411 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 8 Jul 2024 08:56:39 +0200 Subject: [PATCH 22/30] Update doc/specs/stdlib_linalg.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 07dc472e9..e39992dc4 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -927,7 +927,7 @@ or reduced matrices: \( A = Q_1 R_1 \) (provide `shape(Q)==[m,k]`, `shape(R)==[k `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 shape either `[m,n]` or `[k,n]` 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. It minimum size can be determined with a call to [[stdlib_linalg(module):qr_space(interface)]]. It is an `intent(inout)` argument. From e5999f2352246739bcb77495b5439632ae140337 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 8 Jul 2024 08:56:51 +0200 Subject: [PATCH 23/30] Update doc/specs/stdlib_linalg.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index e39992dc4..8360caeac 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -929,7 +929,7 @@ or reduced matrices: \( A = Q_1 R_1 \) (provide `shape(Q)==[m,k]`, `shape(R)==[k `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. It minimum size can be determined with a call to [[stdlib_linalg(module):qr_space(interface)]]. It is an `intent(inout)` argument. +`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(inout)` 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. This is an `intent(in)` argument. From 47f70e68e2e4d59d50026fc9d1c964b04f434584 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 8 Jul 2024 08:57:01 +0200 Subject: [PATCH 24/30] Update doc/specs/stdlib_linalg.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 8360caeac..7f626658f 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -931,7 +931,7 @@ or reduced matrices: \( A = Q_1 R_1 \) (provide `shape(Q)==[m,k]`, `shape(R)==[k `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(inout)` 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. This is an `intent(in)` 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. This is an `intent(out)` argument. From 478f68db2d34c596cceacfb185f827baba80b447 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 8 Jul 2024 08:57:09 +0200 Subject: [PATCH 25/30] Update doc/specs/stdlib_linalg.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 7f626658f..7a6750c11 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -933,7 +933,7 @@ or reduced matrices: \( A = Q_1 R_1 \) (provide `shape(Q)==[m,k]`, `shape(R)==[k `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. This is an `intent(out)` argument. +`err` (optional): Shall be a `type(linalg_state_type)` value. It is an `intent(out)` argument. ### Return value From e3db36ddc5067f6f971d04e9b93644eed8146c10 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 8 Jul 2024 09:20:11 +0200 Subject: [PATCH 26/30] improve qr description --- doc/specs/stdlib_linalg.md | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 7a6750c11..9f46ca087 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -914,8 +914,9 @@ is orthonormal and \( R \) is upper-triangular. Matrix \( A \) has size `[m,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, the user may require to return full matrices (provide `shape(Q)==[m,m]`, `shape(R)==[m,n]`) -or reduced matrices: \( A = Q_1 R_1 \) (provide `shape(Q)==[m,k]`, `shape(R)==[k,n]`). +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 From e69e767e9941045b15473e2129595754265f7370 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 8 Jul 2024 09:23:01 +0200 Subject: [PATCH 27/30] `qr` workspace: make `a` `intent(in)` --- src/stdlib_linalg.fypp | 2 +- src/stdlib_linalg_qr.fypp | 10 +++++----- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 1b8786c9d..bdb00b096 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -503,7 +503,7 @@ module stdlib_linalg #: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(inout), target :: a(:,:) + ${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 diff --git a/src/stdlib_linalg_qr.fypp b/src/stdlib_linalg_qr.fypp index 095e71c44..57bdc59a0 100644 --- a/src/stdlib_linalg_qr.fypp +++ b/src/stdlib_linalg_qr.fypp @@ -93,14 +93,14 @@ submodule (stdlib_linalg) stdlib_linalg_qr ! Get workspace size for QR operations pure module subroutine get_qr_${ri}$_workspace(a,lwork,err) !> Input matrix a[m,n] - ${rt}$, intent(inout), target :: a(:,:) + ${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) + ${rt}$ :: work_dummy(1),tau_dummy(1),a_dummy(1,1) type(linalg_state_type) :: err0 lwork = -1_ilp @@ -108,11 +108,11 @@ submodule (stdlib_linalg) stdlib_linalg_qr !> Problem sizes m = size(a,1,kind=ilp) n = size(a,2,kind=ilp) - k = min(m,n) + k = min(m,n) ! QR space lwork_qr = -1_ilp - call geqrf(m,n,a,m,tau_dummy,work_dummy,lwork_qr,info) + 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) @@ -123,7 +123,7 @@ submodule (stdlib_linalg) stdlib_linalg_qr ! Ordering space (for full factorization) lwork_ord = -1_ilp call #{if rt.startswith('complex')}# ungqr #{else}# orgqr #{endif}# & - (m,m,k,a,m,tau_dummy,work_dummy,lwork_ord,info) + (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) From 8827bf77f61f838656ed2ad1570a5e0f6043a063 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 8 Jul 2024 09:35:03 +0200 Subject: [PATCH 28/30] `qr` temporary storage: make `intent(out)` --- doc/specs/stdlib_linalg.md | 2 +- src/stdlib_linalg.fypp | 2 +- src/stdlib_linalg_qr.fypp | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 9f46ca087..aa3dd8df5 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -930,7 +930,7 @@ the full problem is solved. On reduced matrices (`shape(Q)==[m,k]`, `shape(R)==[ `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(inout)` argument. +`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. diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index bdb00b096..ae0e7ed14 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -478,7 +478,7 @@ module stdlib_linalg !> [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(inout), optional, target :: storage(:) + ${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 diff --git a/src/stdlib_linalg_qr.fypp b/src/stdlib_linalg_qr.fypp index 57bdc59a0..776eabe4c 100644 --- a/src/stdlib_linalg_qr.fypp +++ b/src/stdlib_linalg_qr.fypp @@ -147,7 +147,7 @@ submodule (stdlib_linalg) stdlib_linalg_qr !> [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(inout), optional, target :: storage(:) + ${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 From c6bb171de35fcc4ef0bd412d55cb3a02cafbb979 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 11 Jul 2024 15:11:09 +0200 Subject: [PATCH 29/30] Update src/stdlib_linalg_qr.fypp Co-authored-by: jalvesz <102541118+jalvesz@users.noreply.github.com> --- src/stdlib_linalg_qr.fypp | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/src/stdlib_linalg_qr.fypp b/src/stdlib_linalg_qr.fypp index 776eabe4c..06d772d63 100644 --- a/src/stdlib_linalg_qr.fypp +++ b/src/stdlib_linalg_qr.fypp @@ -181,13 +181,8 @@ submodule (stdlib_linalg) stdlib_linalg_qr use_q_matrix = q1>=m .and. q2>=n ! Can A be overwritten? By default, do not overwrite - if (use_q_matrix) then - overwrite_a_ = .false._lk - elseif (present(overwrite_a)) then - overwrite_a_ = overwrite_a - else - overwrite_a_ = .false._lk - endif + 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 From 6fc0e1d87465723abb499dbdd58b034379a793a2 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 11 Jul 2024 15:16:37 +0200 Subject: [PATCH 30/30] use IEEE NaNs --- test/linalg/test_linalg_qr.fypp | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/test/linalg/test_linalg_qr.fypp b/test/linalg/test_linalg_qr.fypp index 81d3edb0e..8db1d32ad 100644 --- a/test/linalg/test_linalg_qr.fypp +++ b/test/linalg/test_linalg_qr.fypp @@ -6,6 +6,7 @@ module test_linalg_qr 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) @@ -51,10 +52,8 @@ module test_linalg_qr aorig = a ! 1) QR factorization with full matrices. Input NaNs to be sure Q and R are OK on return - q = 0.0_${rk}$ - q = 1.0_${rk}$/q - r = 0.0_${rk}$ - r = 1.0_${rk}$/r + 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