Pivoting QR decomposition #1045

loiseaujc · 2025-10-20T12:50:06Z

Following #930, this PR extends the qr interface to enable the computation of the QR factorization with column pivoting. It is based on the xGEQP3 driver from lapack.

Proposed interfaces

call qr(a, q, r, pivots [, overwrite_a, storage, err])

where a is the matrix to be factorized, q the orthonormal basis for colspan(a), r the upper triangular matrix and pivots an integer array with indices of the pivoted columns.

Progress

Ping: @perazz, @jvdp1, @jalvesz

loiseaujc · 2025-10-22T10:18:01Z

I think this PR is almost ready to be reviewed. I have the specs left to write but this should pretty quick as it mainly is a small modification of the qr specs. I'll try to do that by the end of the week. It seems like there is a small issue with the unbuntu-22.04/cmake/inter-classic 2021.10 setup though. Not sure what it is and I don't have the intel compilers on my laptop so I can't really dig into it at the moment.

Copilot

Pull Request Overview

This PR adds support for QR factorization with column pivoting to the linear algebra library. The implementation leverages the LAPACK geqp3 routine and extends the existing qr interface to support pivot tracking.

Adds pivoting QR factorization via geqp3 LAPACK routine
Extends qr() and qr_space() interfaces to support column pivoting with pivot output
Includes comprehensive test coverage for tall, wide, and rank-deficient matrices

Reviewed Changes

Copilot reviewed 9 out of 9 changed files in this pull request and generated 4 comments.

Show a summary per file

File	Description
test/linalg/test_linalg_pivoting_qr.fypp	Comprehensive test suite for pivoting QR factorization across different matrix types and configurations
test/linalg/CMakeLists.txt	Adds pivoting QR test to build system
src/stdlib_linalg_qr.fypp	Implements pivoting QR factorization routines and workspace calculation
src/stdlib_linalg_lapack.fypp	Adds LAPACK `geqp3` interface for QR with column pivoting
src/stdlib_linalg.fypp	Extends public interface with pivoting QR subroutines
src/lapack/stdlib_linalg_lapack_aux.fypp	Adds error handler for `geqp3` LAPACK routine
example/linalg/example_pivoting_qr_space.f90	Example demonstrating pivoting QR with pre-allocated storage
example/linalg/example_pivoting_qr.f90	Basic example demonstrating pivoting QR factorization
example/linalg/CMakeLists.txt	Adds pivoting QR examples to build system

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

test/linalg/test_linalg_pivoting_qr.fypp

src/stdlib_linalg_lapack.fypp

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>

codecov · 2025-11-04T14:53:12Z

Codecov Report

✅ All modified and coverable lines are covered by tests.
⚠️ Please upload report for BASE (master@9e4c230). Learn more about missing BASE report.

Additional details and impacted files

@@            Coverage Diff            @@
##             master    #1045   +/-   ##
=========================================
  Coverage          ?   25.12%           
=========================================
  Files             ?      570           
  Lines             ?   234201           
  Branches          ?    41275           
=========================================
  Hits              ?    58834           
  Misses            ?   175367           
  Partials          ?        0

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

loiseaujc · 2025-11-05T09:00:59Z

Modulo the issue with the unbuntu-22.04/cmake/intel-classic 2021.10 setup which I haven't been able to investigate yet, this PR is ready for review.

perazz · 2025-11-05T09:18:47Z

doc/specs/stdlib_linalg.md

 ### Syntax

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


Suggested change

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

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

Excellent work @loiseaujc -
Just a preliminary comment - I will have more time this weekend. I would try to not change the previous API. Adding pivots after overwrite_a instead of swapping it still allows running (no keywords)

call qr(a, q, r, pivots)

because storage is kind(A), so, there is no ambiguity.

I made the change in the docs because the actual interface in stdlib_linalg.fypp reads

pure module subroutine stdlib_linalg_${ri}$_qr(a, q, r, overwrite_a, storage, err)

I am not sure why there would be ambiguity though. pivots is an integer array. The call qr(a, q, r, pivots) is exactly the one used in the example.

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

jvdp1

Thank you @loiseaujc . Overall LGTM. I just have some minor suggestions.

jvdp1 · 2025-11-08T13:50:37Z

doc/specs/stdlib_linalg.md

 `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.
+`pivots` (optional): Shall be an `integer` array of size `n`. If provided, QR factorization with column-pivoting is being computed.


Suggested change

`pivots` (optional): Shall be an `integer` array of size `n`. If provided, QR factorization with column-pivoting is being computed.

`pivots` (optional): Shall be an `integer` array of size `n`. If provided, QR factorization with column-pivoting is being computed. It is an `intent(out)` argument.

jvdp1 · 2025-11-08T13:52:05Z

doc/specs/stdlib_linalg.md


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

+`pivoting` (optional): Shall a `logical` flag (default: `.false.`). If `.true.`, on exit `lwork` is the optimal workspace size for the QR factorization with column pivoting. If `.false.`, `lwork` is the optimal workspace size for the standard QR factorization.


Suggested change

`pivoting` (optional): Shall a `logical` flag (default: `.false.`). If `.true.`, on exit `lwork` is the optimal workspace size for the QR factorization with column pivoting. If `.false.`, `lwork` is the optimal workspace size for the standard QR factorization.

`pivoting` (optional): Shall a `logical` flag (default: `.false.`). If `.true.`, on exit `lwork` is the optimal workspace size for the QR factorization with column pivoting. If `.false.`, `lwork` is the optimal workspace size for the standard QR factorization. It is an `intent(in)` argument.

jvdp1 · 2025-11-08T13:57:07Z

src/stdlib_linalg_lapack.fypp

+            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


It seems that all of these interfaces could be defined using a fypp pre-processing approach.

jvdp1 · 2025-11-08T14:02:41Z

test/linalg/test_linalg_pivoting_qr.fypp

+        !-------------------------------
+        !-----     Tall matrix     -----
+        !-------------------------------
+        block


I suggest that each block has its own subroutine. I would facilitate the review and maintainance IMHO.

loiseaujc added 9 commits October 20, 2025 14:17

Function interface for pivoting QR.

98f1708

Added interfaces for xGEQP3.

4b05928

Fix typos in geqp3 interface definition.

0c99ca1

Added handle info function for geqp3.

19cc0f0

Full implementation of pivoting QR.

55bc413

Added test for pivoting QR.

dd9b022

Added test for a tall matrix with rank deficiency.

bd80504

Added example for pivoting QR.

dd785a1

Added pivoting_qr_space example.

e72d2f4

loiseaujc linked an issue Oct 22, 2025 that may be closed by this pull request

Rank-revealing QR decomposition #930

Open

jalvesz requested a review from Copilot November 1, 2025 19:58

Copilot AI reviewed Nov 1, 2025

View reviewed changes

loiseaujc and others added 3 commits November 4, 2025 15:14

Update test/linalg/test_linalg_pivoting_qr.fypp

bf328f6

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>

Update test/linalg/test_linalg_pivoting_qr.fypp

7454ecd

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>

Update test/linalg/test_linalg_pivoting_qr.fypp

07efc4a

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>

loiseaujc added 2 commits November 4, 2025 17:18

Clarify that GEQP3 works both for real and complex matrices.

2bb062a

Added specs.

e2eac02

loiseaujc marked this pull request as ready for review November 5, 2025 09:00

perazz reviewed Nov 5, 2025

View reviewed changes

jvdp1 approved these changes Nov 8, 2025

View reviewed changes

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

	`pivots` (optional): Shall be an `integer` array of size `n`. If provided, QR factorization with column-pivoting is being computed.
	`pivots` (optional): Shall be an `integer` array of size `n`. If provided, QR factorization with column-pivoting is being computed. It is an `intent(out)` 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`.

		`pivoting` (optional): Shall a `logical` flag (default: `.false.`). If `.true.`, on exit `lwork` is the optimal workspace size for the QR factorization with column pivoting. If `.false.`, `lwork` is the optimal workspace size for the standard QR factorization.

Pivoting QR decomposition #1045

Are you sure you want to change the base?

Pivoting QR decomposition #1045

Uh oh!

Conversation

loiseaujc commented Oct 20, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Proposed interfaces

Progress

Uh oh!

loiseaujc commented Oct 22, 2025

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull Request Overview

Reviewed Changes

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

codecov bot commented Nov 4, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

loiseaujc commented Nov 5, 2025

Uh oh!

perazz Nov 5, 2025

Choose a reason for hiding this comment

Uh oh!

loiseaujc Nov 5, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

jvdp1 left a comment

Choose a reason for hiding this comment

Uh oh!

jvdp1 Nov 8, 2025

Choose a reason for hiding this comment

Uh oh!

jvdp1 Nov 8, 2025

Choose a reason for hiding this comment

Uh oh!

jvdp1 Nov 8, 2025

Choose a reason for hiding this comment

Uh oh!

jvdp1 Nov 8, 2025

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

loiseaujc commented Oct 20, 2025 •

edited

Loading

codecov bot commented Nov 4, 2025 •

edited

Loading

loiseaujc Nov 5, 2025 •

edited

Loading