Skip to content

Derived-types for matrix factorizations : QR #1092

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.

Sign up for GitHub

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

Draft
loiseaujc wants to merge 17 commits into
base: master
Choose a base branch
from
+504 −0
Draft

Derived-types for matrix factorizations : QR #1092

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
17 commits
Select commit Hold shift + click to select a range
3d9c025
Merge pull request #2 from loiseaujc/master
loiseaujc Jan 20, 2026
1241e5c
Working implementation for qrfact.
loiseaujc Jan 21, 2026
1459197
Fix error in defining the size of the R matrix.
loiseaujc Jan 21, 2026
6800485
Debug precision error for failing test.
loiseaujc Jan 21, 2026
59ece85
Added test for wide matrices.
loiseaujc Jan 21, 2026
4e9c80f
Re-use stdlib_linalg_eye instead of hand-coded equivalent.
loiseaujc Jan 21, 2026
fdb62ef
Pre-allocate workspace as part of the derived-type.
loiseaujc Jan 21, 2026
c642a21
Add specific tolerance to the test.
loiseaujc Jan 21, 2026
befbd98
resolved conflicts with stdlib v0.8.1
loiseaujc Jan 27, 2026
c251bef
Update CMakeLists to include stdlib_linalg_matrix_factorizations.fypp
loiseaujc Jan 27, 2026
119e08d
Merge pull request #4 from loiseaujc/update_fork
loiseaujc Jan 27, 2026
0fa5fd4
Debug prints.
loiseaujc Jan 27, 2026
d497bf4
More debug prints.
loiseaujc Jan 27, 2026
9b68721
Handling function for gemqrt.
loiseaujc Jan 27, 2026
dbf9fb2
Fix bounds issues in the debug print.
loiseaujc Jan 27, 2026
c5df6d6
Changed error metric.
loiseaujc Jan 27, 2026
cdc4a9c
qrmv and qrmm interfaces.
loiseaujc Jan 28, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
56 changes: 56 additions & 0 deletions src/lapack/stdlib_linalg_lapack_aux.fypp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -53,6 +53,8 @@ module stdlib_linalg_lapack_aux
public :: handle_ggev_info
public :: handle_heev_info
public :: handle_gglse_info
public :: handle_geqrt_info
public :: handle_gemqrt_info

! SELCTG is a LOGICAL FUNCTION of three DOUBLE PRECISION arguments
! used to select eigenvalues to sort to the top left of the Schur form.
Expand Down Expand Up @@ -1654,4 +1656,58 @@ module stdlib_linalg_lapack_aux
end select
end subroutine handle_gglse_info

elemental subroutine handle_geqrt_info(this, info, m, n, nb, ldt, err)
character(len=*), intent(in) :: this
integer(ilp), intent(in) :: info, m, n, nb, ldt
type(linalg_state_type), intent(out) :: err
select case (info)
case(0)
err%state = LINALG_SUCCESS
case(-1, -2)
err = linalg_state_type(this, LINALG_VALUE_ERROR, 'invalid matrix size a=', [m, n])
case(-3)
err = linalg_state_type(this, LINALG_VALUE_ERROR, 'invalid block size nb=', nb)
case(-5)
err = linalg_state_type(this, LINALG_VALUE_ERROR, 'invalid lda=', m)
case(-7)
err = linalg_state_type(this, LINALG_VALUE_ERROR, 'invalid ldt=', ldt)
case default
err = linalg_state_type(this, LINALG_INTERNAL_ERROR, 'catastrophic error.')
end select
end subroutine handle_geqrt_info

elemental subroutine handle_gemqrt_info(this, info, side, trans, m, n, k, nb, ldv, ldt, ldc, err)
character(len=*), intent(in) :: this, side, trans
integer(ilp), intent(in) :: info, m, n, k, nb, ldv, ldt, ldc
type(linalg_state_type), intent(out) :: err
select case (info)
case(0)
err%state = LINALG_SUCCESS
case(-1)
err = linalg_state_type(this, LINALG_VALUE_ERROR, "invalid side argument side=", side)
case(-2)
err = linalg_state_type(this, LINALG_VALUE_ERROR, "invalid operation trans=", trans)
case(-3, -4)
err = linalg_state_type(this, LINALG_VALUE_ERROR, "invalid matrix size c=", [m, n])
case(-5)
err = linalg_state_type(this, LINALG_VALUE_ERROR, "invalid number of elementary reflector k=", k)
case(-6)
err = linalg_state_type(this, LINALG_VALUE_ERROR, "invalid number of blocks nb=", nb)
case(-7)
err = linalg_state_type(this, LINALG_VALUE_ERROR, "invalid matrix size v=", [ldv, k])
case(-8)
err = linalg_state_type(this, LINALG_VALUE_ERROR, "invalid leading dimension for v, ldv=", ldv)
case(-9)
err = linalg_state_type(this, LINALG_VALUE_ERROR, "invalid matrix size t=", [ldt, k])
case(-10)
err = linalg_state_type(this, LINALG_VALUE_ERROR, "invalid leading dimension for t, ldt=", ldc)
case(-11)
err = linalg_state_type(this, LINALG_VALUE_ERROR, "invalid matrix c=", [ldc, n])
case(-12)
err = linalg_state_type(this, LINALG_VALUE_ERROR, "invalid leading dimension for c, ldc=", ldc)
case default
err = linalg_state_type(this, LINALG_INTERNAL_ERROR, "catastrophic error.")
end select
end subroutine handle_gemqrt_info

end module stdlib_linalg_lapack_aux
1 change: 1 addition & 0 deletions src/linalg/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@ set(linalg_fppFiles
stdlib_linalg_kronecker.fypp
stdlib_linalg_least_squares.fypp
stdlib_linalg_matrix_functions.fypp
stdlib_linalg_matrix_factorizations.fypp
stdlib_linalg_norms.fypp
stdlib_linalg_outer_product.fypp
stdlib_linalg_pinv.fypp
Expand Down
98 changes: 98 additions & 0 deletions src/linalg/stdlib_linalg.fypp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,6 +65,8 @@ module stdlib_linalg
public :: is_hermitian
public :: is_triangular
public :: is_hessenberg
! Derived-types for matrix factorizations.
public :: qrfact

! Export linalg error handling
public :: linalg_state_type, linalg_error_handling
Expand Down Expand Up @@ -1911,6 +1913,102 @@ module stdlib_linalg
end subroutine stdlib_linalg_${ri}$_expm
#:endfor
end interface matrix_exp

!----------------------------------------------
!----- DISCIPLINED LINEAR ALGEBRA -----
!----------------------------------------------

#:for rk, rt, ri in RC_KINDS_TYPES
type, public :: qr_${rt[0]}$${rk}$_type
private
integer(ilp) :: m, n
!> Matrix dimensions.
${rt}$, allocatable :: data(:, :)
!> Compact WY representation of the QR factorization.
integer(ilp) :: nb, ldt
!> Block size and leading dimension of the array T.
${rt}$, allocatable :: t(:, :)
!> Upper triangular block reflectors stored in compact format.
${rt}$, allocatable :: work(:)
!> Allocatable workspace.
contains
private
procedure, pass(self) :: get_qfactor_${rt[0]}$${rk}$
generic, public :: Q => get_qfactor_${rt[0]}$${rk}$
!> Construct the Q matrix from the compact representation.
procedure, pass(self) :: get_rfactor_${rt[0]}$${rk}$
generic, public :: R => get_rfactor_${rt[0]}$${rk}$
!> Extract upper triangular R matrix from the compact representation.
end type

interface
module function get_qfactor_${rt[0]}$${rk}$(self) result(Q)
class(qr_${rt[0]}$${rk}$_type), intent(inout) :: self
${rt}$, allocatable :: Q(:, :)
end function get_qfactor_${rt[0]}$${rk}$

pure module function get_rfactor_${rt[0]}$${rk}$(self) result(R)
class(qr_${rt[0]}$${rk}$_type), intent(in) :: self
${rt}$, allocatable :: R(:, :)
end function get_rfactor_${rt[0]}$${rk}$
end interface

#:endfor

interface qrfact
#:for rk, rt, ri in RC_KINDS_TYPES
pure module function stdlib_qrfact_${ri}$(A) result(F)
${rt}$, intent(in) :: A(:, :)
!> Matrix to be factorized.
type(qr_${rt[0]}$${rk}$_type) :: F
!> Derived-type for the QR factorization of A.
end function stdlib_qrfact_${ri}$
#:endfor
end interface

interface qrmv
#:for rk, rt, ri in RC_KINDS_TYPES
module subroutine stdlib_qrmv_${ri}$(A, x, y, side, op, err)
type(qr_${rt[0]}$${rk}$_type), intent(in) :: A
!> QR-factorized matrix.
${rt}$, intent(inout) :: x(:)
!> Vector to be multiplied with (either from left or right).
!> If y is not passed, the matrix-vector product is performed in-place.
${rt}$, optional, intent(out) :: y(:, :)
!> [optional] Resulting vector for out-of-place computations.
character(len=*), intent(in) :: side
!> From which side to multiply (side = "L" or "R", "L" being the default).
character(len=*), intent(in) :: op
!> Whether we multiply by A or A.T (A.H for complex matrices).
!> Choose from op = "N", op = "T" (real matrices) or op = "H" (complex matrices).
!> Default is op = "N".
type(linalg_state_type), optional, intent(out) :: err
!>
end subroutine stdlib_qrmv_${ri}$
#:endfor
end interface qrmv

interface qrmm
#:for rk, rt, ri in RC_KINDS_TYPES
module subroutine stdlib_qrmm_${ri}$(A, X, Y, side, op, err)
type(qr_${rt[0]}$${rk}$_type), intent(in) :: A
!> QR-factorized matrix.
${rt}$, intent(inout) :: X(:, :)
!> Matrix to be multiplied with (either from left or right).
!> If Y is not passed, the matrix-matrix product is performed in-place.
${rt}$, optional, intent(out) :: Y(:, :)
!> [optional] Resulting matrix for out-of-place computations.
character(len=*), intent(in) :: side
!> From which side to multiply (side = "L" or "R", "L" being the default).
character(len=*), intent(in) :: op
!> Whether we multiply by A or A.T (A.H for complex matrices).
!> Choose from op = "N", op = "T" (real matrices) or op = "H" (complex matrices).
!> Default is op = "N".
type(linalg_state_type), optional, intent(out) :: err
end subroutine stdlib_qrmm_${ri}$
#:endfor
end interface qrmm

contains


Expand Down
96 changes: 96 additions & 0 deletions src/linalg/stdlib_linalg_matrix_factorizations.fypp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,96 @@
#:include "common.fypp"
#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
! Cholesky factorization of a matrix, based on LAPACK *POTRF functions
submodule (stdlib_linalg) stdlib_linalg_matrix_factorizations
use stdlib_linalg_lapack, only: geqrt, gemqrt
use stdlib_linalg_lapack_aux, only: handle_geqrt_info, handle_gemqrt_info
implicit none

character(len=*), parameter :: this = 'matrix_factorizations'

contains

!-----------------------------------------------------
!----- DERIVED-TYPE FOR QR FACTORIZATION -----
!-----------------------------------------------------

#:for rk, rt, ri in RC_KINDS_TYPES
pure module function stdlib_qrfact_${ri}$(A) result(F)
${rt}$, intent(in) :: A(:, :)
!> Matrix to be factorized.
type(qr_${rt[0]}$${rk}$_type) :: F

!----- Internal variables -----
integer(ilp), parameter :: min_nb = 36 ! Minimum block size.
integer(ilp) :: m, n, nb, info
${rt}$, allocatable :: work(:)
type(linalg_state_type) :: err

!> Problem dimensions.
m = size(A, 1) ! Number of rows.
n = size(A, 2) ! Number of columns.
nb = min(m, n, min_nb) ! Block size.

!> Copy matrix.
F%data = A
F%ldt = nb

!> Allocate data.
allocate(F%t(F%ldt, min(m, n)), F%work(nb*max(m, n)))

!> QR factorization.
call geqrt(m, n, nb, F%data, m, F%t, F%ldt, F%work, info)
call handle_geqrt_info(this, info, m, n, nb, F%ldt, err)
call linalg_error_handling(err)
end function stdlib_qrfact_${ri}$

module function get_qfactor_${rt[0]}$${rk}$(self) result(Q)
class(qr_${rt[0]}$${rk}$_type), intent(inout) :: self
!> Compact representation of the QR factorization.
${rt}$, allocatable :: Q(:, :)
!> Orthogonal matrix.

!----- Internal variables -----
integer(ilp) :: i, j, k, info
character(len=1), parameter :: side='L', trans='N'
${rt}$, parameter :: zero = #{if rt.startswith("r")}#0.0_${rk}$#{else}#cmplx(0.0_${rk}$, 0.0_${rk}$, kind=${rk}$)#{endif}#
type(linalg_state_type) :: err

associate( m => size(self%data, 1), n => size(self%data, 2), &
nb => self%ldt, ldv => size(self%data, 1), ldt => self%ldt, ldc => size(self%data, 1))
!> Compute the Q matrix.
k = min(m, n)
Q = eye(m, k, mold=zero)
call gemqrt(side, trans, m, k, k, &
nb, self%data(:, :k), ldv, self%t, ldt, Q, ldc, self%work, info)
call handle_gemqrt_info(this, info, side, trans, m, n, k, nb, ldv, ldt, ldc, err)
call linalg_error_handling(err)
end associate
end function get_qfactor_${rt[0]}$${rk}$

pure module function get_rfactor_${rt[0]}$${rk}$(self) result(R)
class(qr_${rt[0]}$${rk}$_type), intent(in) :: self
!> Compact representation of the QR factorization.
${rt}$, allocatable :: R(:, :)
!> Upper-triangular matrix.

!----- Internal variables -----
integer(ilp) :: i, j, k, m, n
${rt}$, parameter :: zero = #{if rt.startswith("r")}# 0.0_${rk}$#{else}#cmplx(0.0_${rk}$, 0.0_${rk}$, kind=${rk}$) #{endif}#

!> Matrix size.
m = size(self%data, 1)
n = size(self%data, 2)
k = min(m, n)

!> Allocate R matrix.
allocate(R(k, n), source=zero)

!> Fill-in matrix.
do concurrent(i=1:k, j=1:n, i <= j)
R(i, j) = self%data(i, j)
enddo
end function get_rfactor_${rt[0]}$${rk}$
#:endfor
end submodule stdlib_linalg_matrix_factorizations

2 changes: 2 additions & 0 deletions test/linalg/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@ set(
"test_linalg_specialmatrices.fypp"
"test_linalg_cholesky.fypp"
"test_linalg_expm.fypp"
"test_linalg_matrix_factorizations.fypp"
)

# Preprocessed files to contain preprocessor directives -> .F90
Expand Down Expand Up @@ -55,4 +56,5 @@ ADDTEST(linalg_sparse)
if (STDLIB_SPECIALMATRICES)
ADDTEST(linalg_specialmatrices)
endif()
ADDTEST(linalg_matrix_factorizations)
ADDTESTPP(blas_lapack)
Loading
Loading