-
Notifications
You must be signed in to change notification settings - Fork 198
Pivoting QR decomposition #1045
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
base: master
Are you sure you want to change the base?
Conversation
|
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 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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
geqp3LAPACK routine - Extends
qr()andqr_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.
Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## master #1045 +/- ##
=========================================
Coverage ? 25.12%
=========================================
Files ? 570
Lines ? 234201
Branches ? 41275
=========================================
Hits ? 58845
Misses ? 175356
Partials ? 0 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
|
Modulo the issue with the |
jvdp1
left a comment
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thank you @loiseaujc . Overall LGTM. I just have some minor suggestions.
| 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 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It seems that all of these interfaces could be defined using a fypp pre-processing approach.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yeah sure. I copied what looked like the standard way to do it in the stdlib_linalg_lapack.fypp file, eg.
#:for ik,it,ii in LINALG_INT_KINDS_TYPES
#ifdef STDLIB_EXTERNAL_LAPACK${ii}$
pure recursive subroutine cgeqrt3( m, n, a, lda, t, ldt, info )
import sp,dp,qp,${ik}$,lk
implicit none
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda,m,n,ldt
complex(sp), intent(inout) :: a(lda,*)
complex(sp), intent(out) :: t(ldt,*)
end subroutine cgeqrt3
#else
module procedure stdlib${ii}$_cgeqrt3
#endif
#ifdef STDLIB_EXTERNAL_LAPACK${ii}$
pure recursive subroutine dgeqrt3( m, n, a, lda, t, ldt, info )
import sp,dp,qp,${ik}$,lk
implicit none
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda,m,n,ldt
real(dp), intent(inout) :: a(lda,*)
real(dp), intent(out) :: t(ldt,*)
end subroutine dgeqrt3
#else
module procedure stdlib${ii}$_dgeqrt3
#endif
...
but simply lumped everything into a single #ifdef block. Are there any technical reasons why you guys went this way initially rather than fypp-ified the whole prefix thing ?
| !------------------------------- | ||
| !----- Tall matrix ----- | ||
| !------------------------------- | ||
| block |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I suggest that each block has its own subroutine. I would facilitate the review and maintainance IMHO.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Will do right away.
perazz
left a comment
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thank you for this implementation @loiseaujc, great work.
I stand by @jvdp1's comments, and I added a few edits/typos.
Overall, looks good to me with the above.
Co-authored-by: Federico Perini <federico.perini@gmail.com>
Co-authored-by: Federico Perini <federico.perini@gmail.com>
Co-authored-by: Federico Perini <federico.perini@gmail.com>
Co-authored-by: Jeremie Vandenplas <jeremie.vandenplas@gmail.com>
Co-authored-by: Jeremie Vandenplas <jeremie.vandenplas@gmail.com>
Co-authored-by: Federico Perini <federico.perini@gmail.com>
Co-authored-by: Federico Perini <federico.perini@gmail.com>
Following #930, this PR extends the
qrinterface to enable the computation of the QR factorization with column pivoting. It is based on thexGEQP3driver fromlapack.Proposed interfaces
call qr(a, q, r, pivots [, overwrite_a, storage, err])where
ais the matrix to be factorized,qthe orthonormal basis forcolspan(a),rthe upper triangular matrix andpivotsan integer array with indices of the pivoted columns.Progress
Ping: @perazz, @jvdp1, @jalvesz