From 98f1708834e094b3bfe58a77c7a7604d95965cac Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Mon, 20 Oct 2025 14:17:35 +0200 Subject: [PATCH 1/9] Function interface for pivoting QR. --- src/stdlib_linalg.fypp | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 98fdfa172..93e197174 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -614,6 +614,23 @@ 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 + + pure module subroutine stdlib_linalg_${ri}$_pivoting_qr(a, q, r, pivots, 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(:, :) + !> Pivots. + integer(ilp), intent(out) :: pivots(:) + !> [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 pivoting_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}$_pivoting_qr #:endfor end interface qr From 4b05928456994dce0c982d1c697a614ea0d3feb2 Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Wed, 22 Oct 2025 10:21:56 +0200 Subject: [PATCH 2/9] Added interfaces for xGEQP3. --- src/stdlib_linalg_lapack.fypp | 71 +++++++++++++++++++++++++++++++++++ 1 file changed, 71 insertions(+) diff --git a/src/stdlib_linalg_lapack.fypp b/src/stdlib_linalg_lapack.fypp index d906e9df0..4ea07402e 100644 --- a/src/stdlib_linalg_lapack.fypp +++ b/src/stdlib_linalg_lapack.fypp @@ -3230,6 +3230,77 @@ module stdlib_linalg_lapack #:endfor end interface geqrt3 + interface geqp3 + !! GEQP3 computes a QR factorization with column pivoting of a complex + !! M-by-N matrix A: + !! + !! A * P = Q * R, + !! + !! where: + !! Q is an M-by-min(M, N) orthogonal matrix + !! R is an min(M, N)-by-N upper triangular matrix; +#:for ik, it, ii in LINALG_INT_KINDS_TYPES +#ifdef STDLIB_EXTERNAL_LAPACK${ii}$ + pure subroutine sgeqp3(m, n, a, lda, jpvt, tau, work, lwork, info) + import sp, dp, qp, ${ik}$, lk + implicit none + integer(${ik}$), intent(in) :: m, n, lda, lwork + integer(${ik}$), intent(out) :: info + integer(${ik}$), intent(inout) :: jpvt(*) + real(sp), intent(inout) :: a(lda, *) + real(sp), intent(out) :: tau(*), work(*) + end subroutine sgeqp3 + + pure subroutine dgeqp3(m, n, a, lda, jpvt, tau, work, lwork, info) + import sp, dp, qp, ${ik}$, lk + implicit none + integer(${ik}$), intent(in) :: m, n, lda, lwork + integer(${ik}$), intent(out) :: info + integer(${ik}$), intent(inout) :: jpvt(*) + real(dp), intent(inout) :: a(lda, *) + real(dp), intent(out) :: tau(*), work(*) + end subroutine dgeqp3 + + pure subroutine cgeqp3(m, n, a, lda, jpvt, tau, work, lwork, rwork, info) + import sp, dp, qp, ${ik}$, lk + implicit none + integer(${ik}$), intent(in) :: m, n, lda, lwork + integer(${ik}$), intent(out) :: info + integer(${ik}$), intent(inout) :: jpvt(*) + complex(sp), intent(inout) :: a(lda, *) + complex(sp), intent(out) :: tau(*), work(*) + real(sp), intent(out) :: rwork(*) + end subroutine cgeqp3 + + pure subroutine zgeqp3(m, n, a, lda, jpvt, tau, work, lwork, rwork, info) + import sp, dp, qp, ${ik}$, lk + implicit none + integer(${ik}$), intent(in) :: m, n, lda, lwork + integer(${ik}$), intent(out) :: info + integer(${ik}$), intent(inout) :: jpvt(*) + complex(dp), intent(inout) :: a(lda, *) + complex(dp), intent(out) :: tau(*), work(*) + real(dp), intent(out) :: rwork(*) + end subroutine zgeqp3 +#else + module procedure stdlib${ii}$_sgeqp3 + module procedure stdlib${ii}$_dgeqp3 + module procedure stdlib${ii}$_cgeqp3 + module procedure stdlib${ii}$_zgeqp3 +#endif +#:for rk, rt, ri in REAL_KINDS_TYPES +#:if not rk in ["sp", "dp"] + module procedure stdlib${ii}$_${ri}$geqrf +#:endif +#:endfor +#:for rk, rt, ri in CMPLX_KINDS_TYPES +#:if not rk in ["sp", "dp"] + module procedure stdlib${ii}$_${ri}$geqrf +#:endif +#:endfor +#:endfor + end interface geqp3 + interface gerfs !! GERFS improves the computed solution to a system of linear !! equations and provides error bounds and backward error estimates for From 0c99ca1981f1bbb58a2b1b53b53b45079280ca1a Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Wed, 22 Oct 2025 10:38:08 +0200 Subject: [PATCH 3/9] Fix typos in geqp3 interface definition. --- src/stdlib_linalg_lapack.fypp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/stdlib_linalg_lapack.fypp b/src/stdlib_linalg_lapack.fypp index 4ea07402e..a2c891bd1 100644 --- a/src/stdlib_linalg_lapack.fypp +++ b/src/stdlib_linalg_lapack.fypp @@ -3290,12 +3290,12 @@ module stdlib_linalg_lapack #endif #:for rk, rt, ri in REAL_KINDS_TYPES #:if not rk in ["sp", "dp"] - module procedure stdlib${ii}$_${ri}$geqrf + module procedure stdlib${ii}$_${ri}$geqp3 #:endif #:endfor #:for rk, rt, ri in CMPLX_KINDS_TYPES #:if not rk in ["sp", "dp"] - module procedure stdlib${ii}$_${ri}$geqrf + module procedure stdlib${ii}$_${ri}$geqp3 #:endif #:endfor #:endfor From 19cc0f0b9bb65ea729593bf6bd2de39b9cdb8029 Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Wed, 22 Oct 2025 10:46:14 +0200 Subject: [PATCH 4/9] Added handle info function for geqp3. --- src/lapack/stdlib_linalg_lapack_aux.fypp | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/src/lapack/stdlib_linalg_lapack_aux.fypp b/src/lapack/stdlib_linalg_lapack_aux.fypp index 442999ca3..640a2025b 100644 --- a/src/lapack/stdlib_linalg_lapack_aux.fypp +++ b/src/lapack/stdlib_linalg_lapack_aux.fypp @@ -46,6 +46,7 @@ module stdlib_linalg_lapack_aux public :: handle_gesv_info public :: handle_gees_info public :: handle_geqrf_info + public :: handle_geqp3_info public :: handle_orgqr_info public :: handle_gelsd_info public :: handle_geev_info @@ -1462,6 +1463,27 @@ module stdlib_linalg_lapack_aux end subroutine handle_geqrf_info + elemental subroutine handle_geqp3_info(this, info, m, n, lwork, err) + character(len=*), intent(in) :: this + 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_VALUE_ERROR, 'invalid input for lwork=', lwork) + case default + err = linalg_state_type(this, LINALG_VALUE_ERROR, 'catastrophic error') + end select + end subroutine handle_geqp3_info + elemental subroutine handle_orgqr_info(this,info,m,n,k,lwork,err) character(len=*), intent(in) :: this integer(ilp), intent(in) :: info,m,n,k,lwork From 55bc413de77395ffa69f734d6d526fc6f5ef4c6f Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Wed, 22 Oct 2025 11:42:59 +0200 Subject: [PATCH 5/9] Full implementation of pivoting QR. --- src/stdlib_linalg.fypp | 11 +++ src/stdlib_linalg_qr.fypp | 201 +++++++++++++++++++++++++++++++++++++- 2 files changed, 210 insertions(+), 2 deletions(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 93e197174..784836d14 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -658,6 +658,17 @@ 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 + + pure module subroutine get_pivoting_qr_${ri}$_workspace(a, lwork, pivoting, err) + !> Input matrix a[m, n] + ${rt}$, intent(in), target :: a(:, :) + !> Minimum workspace size for both operations. + integer(ilp), intent(out) :: lwork + !> Pivoting flag. + logical(lk), intent(in) :: pivoting + !> State return flag. Returns an error if the query failed. + type(linalg_state_type), optional, intent(out) :: err + end subroutine get_pivoting_qr_${ri}$_workspace #:endfor end interface qr_space diff --git a/src/stdlib_linalg_qr.fypp b/src/stdlib_linalg_qr.fypp index 2212c39c3..eb1d603e0 100644 --- a/src/stdlib_linalg_qr.fypp +++ b/src/stdlib_linalg_qr.fypp @@ -2,8 +2,8 @@ #: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_lapack_aux, only: handle_geqrf_info, handle_orgqr_info + use stdlib_linalg_lapack, only: geqrf, geqp3, orgqr, ungqr + use stdlib_linalg_lapack_aux, only: handle_geqrf_info, handle_orgqr_info, handle_geqp3_info use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, & LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR implicit none @@ -220,4 +220,201 @@ submodule (stdlib_linalg) stdlib_linalg_qr #:endfor + !--------------------------------------------------------- + !----- QR decomposition with column pivoting ----- + !--------------------------------------------------------- + + #:for rk, rt, ri in RC_KINDS_TYPES + ! Get workspace size for QR operations + pure module subroutine get_pivoting_qr_${ri}$_workspace(a,lwork,pivoting,err) + !> Input matrix a[m,n] + ${rt}$, intent(in), target :: a(:,:) + !> Minimum workspace size for both operations + integer(ilp), intent(out) :: lwork + !> Pivoting flag. + logical(lk), intent(in) :: pivoting + !> 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) + integer(ilp) :: jpvt_dummy(1) + real(${rk}$) :: rwork_dummy(1) + type(linalg_state_type) :: err0 + + if (pivoting) then + 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 + #:if rt.startswith('complex') + call geqp3(m, n, a_dummy, m, jpvt_dummy, tau_dummy, work_dummy, lwork_qr, rwork_dummy, info) + #:else + call geqp3(m, n, a_dummy, m, jpvt_dummy, tau_dummy, work_dummy, lwork_qr, info) + #:endif + call handle_geqp3_info(this, 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(this,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) + else + call qr_space(a, lwork, err) + endif + + end subroutine get_pivoting_qr_${ri}$_workspace + + pure module subroutine stdlib_linalg_${ri}$_pivoting_qr(a, q, r, pivots, 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(:, :) + !> Pivots. + integer(ilp), intent(out) :: pivots(:) + !> [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 pivoting_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(:) + #:if rt.startswith('complex') + real(${rk}$) :: rwork(2*size(a, 2, kind=ilp)) + #:endif + + !> 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) + pivots = 0_ilp + + !> Full or thin QR factorization ? + call check_problem_size(m, n, q1, q2, r1, r2, err0, reduced) + if (err0%error()) then + call linalg_error_handling(err0, err) + return + endif + + !> Can Q be used as temporary storage for A, + ! to be destroyed by *GEQP3. + use_q_matrix = q1 >= m .and. q2 >= n + + !> Can A be overwritten ? (By default, no). + overwrite_a_ = .false._lk + if (present(overwrite_a) .and. .not. use_q_matrix) overwrite_a_ = overwrite_a + + !> Initialize a temporary matrix or reuse available storage if possible. + if (use_q_matrix) then + amat => q + q(:m, :n) = a + else if (overwrite_a_) then + amat => a + else + allocate(amat(m, n), source=a) + endif + lda = size(amat, 1, kind=ilp) + + !> Store the elementary reflectors. + 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 + ! r11 element. Use the full column and store it in a dummy variable. + tau(1:k) => r(1:k, 1) + endif + + ! Retrieve workspace size. + call get_pivoting_qr_${ri}$_workspace(a, lwork, .true., 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. + #:if rt.startswith('complex') + call geqp3(m, n, amat, m, pivots, tau, work, lwork, rwork, info) + #:else + call geqp3(m, n, amat, m, pivots, tau, work, lwork, info) + #:endif + call handle_geqp3_info(this, 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 used by `tau` + r11 = amat(1, 1) + do concurrent(i=1:min(r1, m), j=2:n) + r(i, j) = merge(amat(i, j), zero, i <= j) + enddo + + ! 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(this,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}$_pivoting_qr + #:endfor + end submodule stdlib_linalg_qr From dd9b0225f078ce7bc22948d990f2356cf14ee688 Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Wed, 22 Oct 2025 11:43:05 +0200 Subject: [PATCH 6/9] Added test for pivoting QR. --- test/linalg/CMakeLists.txt | 2 + test/linalg/test_linalg_pivoting_qr.fypp | 245 +++++++++++++++++++++++ 2 files changed, 247 insertions(+) create mode 100644 test/linalg/test_linalg_pivoting_qr.fypp diff --git a/test/linalg/CMakeLists.txt b/test/linalg/CMakeLists.txt index a03c031ec..ae0ca0b34 100644 --- a/test/linalg/CMakeLists.txt +++ b/test/linalg/CMakeLists.txt @@ -10,6 +10,7 @@ set( "test_linalg_mnorm.fypp" "test_linalg_determinant.fypp" "test_linalg_qr.fypp" + "test_linalg_pivoting_qr.fypp" "test_linalg_schur.fypp" "test_linalg_solve_iterative.fypp" "test_linalg_svd.fypp" @@ -42,6 +43,7 @@ ADDTEST(linalg_mnorm) ADDTEST(linalg_solve) ADDTEST(linalg_lstsq) ADDTEST(linalg_qr) +ADDTEST(linalg_pivoting_qr) ADDTEST(linalg_schur) ADDTEST(linalg_solve_iterative) ADDTEST(linalg_svd) diff --git a/test/linalg/test_linalg_pivoting_qr.fypp b/test/linalg/test_linalg_pivoting_qr.fypp new file mode 100644 index 000000000..db1fc6ded --- /dev/null +++ b/test/linalg/test_linalg_pivoting_qr.fypp @@ -0,0 +1,245 @@ +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +! Test QR factorization +module test_linalg_pivoting_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_pivoting_qr_factorization + + contains + + !> QR factorization tests + subroutine test_pivoting_qr_factorization(tests) + !> Collection of tests + type(unittest_type), allocatable, intent(out) :: tests(:) + + allocate(tests(0)) + + #:for rk,rt,ri in RC_KINDS_TYPES + call add_test(tests,new_unittest("pivoting_qr_random_${ri}$",test_pivoting_qr_random_${ri}$)) + #:endfor + + end subroutine test_pivoting_qr_factorization + + !> QR factorization of a random matrix + #:for rk,rt,ri in RC_KINDS_TYPES + subroutine test_pivoting_qr_random_${ri}$(error) + type(error_type), allocatable, intent(out) :: error + + !------------------------------- + !----- Tall matrix ----- + !------------------------------- + block + 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) :: pivots(n), i, j + 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, pivots, err=state) + + ! Check return code + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + ! Check solution + call check(error, all(abs(a(:, pivots)-matmul(q,r))0) new_tests(1:n) = tests(1:n) + new_tests(1+n) = new_test + call move_alloc(from=new_tests,to=tests) + + end subroutine add_test + +end module test_linalg_pivoting_qr + +program test_pivoting_qr + use, intrinsic :: iso_fortran_env, only : error_unit + use testdrive, only : run_testsuite, new_testsuite, testsuite_type + use test_linalg_pivoting_qr, only : test_pivoting_qr_factorization + implicit none + integer :: stat, is + type(testsuite_type), allocatable :: testsuites(:) + character(len=*), parameter :: fmt = '("#", *(1x, a))' + + stat = 0 + + testsuites = [ & + new_testsuite("linalg_pivoting_qr", test_pivoting_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_pivoting_qr From bd805041433e3cf751e9f924bb7a108268fc7e4d Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Wed, 22 Oct 2025 11:46:02 +0200 Subject: [PATCH 7/9] Added test for a tall matrix with rank deficiency. --- test/linalg/test_linalg_pivoting_qr.fypp | 81 ++++++++++++++++++++++++ 1 file changed, 81 insertions(+) diff --git a/test/linalg/test_linalg_pivoting_qr.fypp b/test/linalg/test_linalg_pivoting_qr.fypp index db1fc6ded..3f19b6160 100644 --- a/test/linalg/test_linalg_pivoting_qr.fypp +++ b/test/linalg/test_linalg_pivoting_qr.fypp @@ -112,6 +112,87 @@ module test_linalg_pivoting_qr if (allocated(error)) return end block + !---------------------------------------------------- + !----- Tall matrix with rank deficiency ----- + !---------------------------------------------------- + block + 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) :: pivots(n), i, j + 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 + a(:, 3) = 0.0_${rk}$ ! Zero-out column to test rank-deficiency. + 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, pivots, err=state) + + ! Check return code + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + ! Check solution + call check(error, all(abs(a(:, pivots)-matmul(q,r)) Date: Wed, 22 Oct 2025 12:00:32 +0200 Subject: [PATCH 8/9] Added example for pivoting QR. --- example/linalg/CMakeLists.txt | 1 + example/linalg/example_pivoting_qr.f90 | 16 ++++++++++++++++ 2 files changed, 17 insertions(+) create mode 100644 example/linalg/example_pivoting_qr.f90 diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index b00cfe1d4..bb467a888 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -54,6 +54,7 @@ ADD_EXAMPLE(svdvals) ADD_EXAMPLE(determinant) ADD_EXAMPLE(determinant2) ADD_EXAMPLE(qr) +ADD_EXAMPLE(pivoting_qr) ADD_EXAMPLE(qr_space) ADD_EXAMPLE(cholesky) ADD_EXAMPLE(chol) diff --git a/example/linalg/example_pivoting_qr.f90 b/example/linalg/example_pivoting_qr.f90 new file mode 100644 index 000000000..b2f8e1f23 --- /dev/null +++ b/example/linalg/example_pivoting_qr.f90 @@ -0,0 +1,16 @@ +program example_pivoting_qr + use stdlib_linalg, only: qr + implicit none + real :: A(104, 32), Q(104, 32), R(32, 32) + integer :: pivots(32) + + ! Create a random matrix + call random_number(A) + + ! Compute its QR factorization (reduced) + call qr(A, Q, R, pivots) + + ! Test factorization: Q*R = A + print *, maxval(abs(matmul(Q, R) - A(:, pivots))) + +end program example_pivoting_qr From e72d2f4e82026a6f9399f58c2e4f2900121ba3b7 Mon Sep 17 00:00:00 2001 From: Jean-Christophe Date: Wed, 22 Oct 2025 12:04:42 +0200 Subject: [PATCH 9/9] Added pivoting_qr_space example. --- example/linalg/CMakeLists.txt | 1 + example/linalg/example_pivoting_qr_space.f90 | 25 ++++++++++++++++++++ 2 files changed, 26 insertions(+) create mode 100644 example/linalg/example_pivoting_qr_space.f90 diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index bb467a888..fe243ecd3 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -56,6 +56,7 @@ ADD_EXAMPLE(determinant2) ADD_EXAMPLE(qr) ADD_EXAMPLE(pivoting_qr) ADD_EXAMPLE(qr_space) +ADD_EXAMPLE(pivoting_qr_space) ADD_EXAMPLE(cholesky) ADD_EXAMPLE(chol) ADD_EXAMPLE(expm) diff --git a/example/linalg/example_pivoting_qr_space.f90 b/example/linalg/example_pivoting_qr_space.f90 new file mode 100644 index 000000000..920a91d88 --- /dev/null +++ b/example/linalg/example_pivoting_qr_space.f90 @@ -0,0 +1,25 @@ +! QR example with pre-allocated storage +program example_pivoting_qr_space + use stdlib_linalg_constants, only: ilp + use stdlib_linalg, only: qr, qr_space, linalg_state_type + implicit none + real :: A(104, 32), Q(104, 32), R(32, 32) + real, allocatable :: work(:) + integer(ilp) :: lwork, pivots(32) + type(linalg_state_type) :: err + + ! Create a random matrix + call random_number(A) + + ! Prepare QR workspace + call qr_space(A, lwork, pivoting=.true.) + allocate (work(lwork)) + + ! Compute its QR factorization (reduced) + call qr(A, Q, R, pivots, storage=work, err=err) + + ! Test factorization: Q*R = A + print *, maxval(abs(matmul(Q, R) - A(:, pivots))) + print *, err%print() + +end program example_pivoting_qr_space