diff --git a/src/lapack/stdlib_linalg_lapack_aux.fypp b/src/lapack/stdlib_linalg_lapack_aux.fypp
index 475424c91..c63194344 100644
--- a/src/lapack/stdlib_linalg_lapack_aux.fypp
+++ b/src/lapack/stdlib_linalg_lapack_aux.fypp
@@ -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. 
@@ -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
diff --git a/src/linalg/CMakeLists.txt b/src/linalg/CMakeLists.txt
index e0a6c10dd..b216803b2 100644
--- a/src/linalg/CMakeLists.txt
+++ b/src/linalg/CMakeLists.txt
@@ -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
diff --git a/src/linalg/stdlib_linalg.fypp b/src/linalg/stdlib_linalg.fypp
index 06169e071..c81a56ea2 100644
--- a/src/linalg/stdlib_linalg.fypp
+++ b/src/linalg/stdlib_linalg.fypp
@@ -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
@@ -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
 
 
diff --git a/src/linalg/stdlib_linalg_matrix_factorizations.fypp b/src/linalg/stdlib_linalg_matrix_factorizations.fypp
new file mode 100644
index 000000000..b02f4c4b8
--- /dev/null
+++ b/src/linalg/stdlib_linalg_matrix_factorizations.fypp
@@ -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
+
diff --git a/test/linalg/CMakeLists.txt b/test/linalg/CMakeLists.txt
index f27dedd4f..f25ad0870 100644
--- a/test/linalg/CMakeLists.txt
+++ b/test/linalg/CMakeLists.txt
@@ -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
@@ -55,4 +56,5 @@ ADDTEST(linalg_sparse)
 if (STDLIB_SPECIALMATRICES)
   ADDTEST(linalg_specialmatrices)
 endif()
+ADDTEST(linalg_matrix_factorizations)
 ADDTESTPP(blas_lapack)
diff --git a/test/linalg/test_linalg_matrix_factorizations.fypp b/test/linalg/test_linalg_matrix_factorizations.fypp
new file mode 100644
index 000000000..2358bd1e6
--- /dev/null
+++ b/test/linalg/test_linalg_matrix_factorizations.fypp
@@ -0,0 +1,298 @@
+#:include "common.fypp"
+#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
+! Test QR factorization
+module test_linalg_matrix_factorizations
+    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: qrfact, qr, mnorm
+    use stdlib_math, only: all_close
+    use ieee_arithmetic, only: ieee_value,ieee_quiet_nan
+
+    implicit none (type,external)
+
+    public :: test_qr_factorization
+
+    #:for rk, rt, ri in REAL_KINDS_TYPES
+    ${rt}$, parameter :: rel_tol_${rk}$ = sqrt(epsilon(1.0_${rk}$))
+    #:endfor
+
+    contains
+
+    !> QR factorization tests
+    subroutine test_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("qr_random_tall_matrix_${ri}$",test_qr_random_tall_matrix_${ri}$))
+        call add_test(tests,new_unittest("qr_random_wide_matrix_${ri}$",test_qr_random_wide_matrix_${ri}$))
+        call add_test(tests,new_unittest("qr_random_tall_rank_deficient_matrix_${ri}$",test_qr_random_tall_rank_deficient_matrix_${ri}$))
+        #:endfor
+    end subroutine test_qr_factorization
+
+    !> QR factorization of a random matrix
+    #:for rk,rt,ri in RC_KINDS_TYPES
+    subroutine test_qr_random_tall_matrix_${ri}$(error)
+        use stdlib_linalg, only: qr_${rt[0]}$${rk}$_type
+        type(error_type), allocatable, intent(out) :: error
+        integer(ilp), parameter :: m   = 15_ilp
+        integer(ilp), parameter :: n   =  4_ilp
+        integer(ilp), parameter :: k   = min(m,n)
+        ${rt}$ :: a(m,n), q(m, k), r(k, n)
+        real(${rk}$) :: rea(m,n),ima(m,n)
+        type(linalg_state_type) :: state
+        type(qr_${rt[0]}$${rk}$_type) :: F
+
+        ! Random matrix with unit 2-norm.
+        call random_number(rea)
+        #:if rt.startswith('complex')
+        call random_number(ima)
+        a = cmplx(rea,ima,kind=${rk}$)
+        #:else
+        a = rea
+        #:endif
+        a = a / mnorm(a, 2)
+
+        ! Reference QR decomposition.
+        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)
+        call check(error,state%ok(),state%print())
+        if (allocated(error)) return
+
+        ! Instantiate derived-type for the QR.
+        F = qrfact(A)
+
+        ! Check the R matrix with reference.
+        call check(error, all(abs(R - F%R()) < rel_tol_${rk}$))
+        if (allocated(error)) then
+            block
+            integer(ilp) :: i, j
+            ${rt}$, allocatable :: rnew(:, :)
+            print *, "Reference R matrix. shape(R) =", shape(R)
+            do i = 1, k
+                print *, (R(i, j), j = 1, n)
+            enddo
+            print *
+            print *, "New R matrix. shape(R) =", shape(F%R())
+            rnew = F%R()
+            do i = 1, k
+                print *, (rnew(i, j), j = 1, n)
+            enddo
+            print *
+            print *, "Difference between the two."
+            rnew = rnew - r
+            do i = 1, k
+                print *, (abs(rnew(i, j)), j=1, n)
+            enddo
+            end block
+            return
+        endif
+
+        ! Check the Q matrix with reference.
+        call check(error, all(abs(Q - F%Q()) < rel_tol_${rk}$))
+        if (allocated(error)) then
+            block
+            integer(ilp) :: i, j
+            ${rt}$, allocatable :: qnew(:, :)
+            print *, "Reference Q matrix. shape(Q) =", shape(Q)
+            do i = 1, m
+                print *, (Q(i, j), j = 1, k)
+            enddo
+            print *
+            print *, "New Q matrix. shape(Q) =", shape(F%Q())
+            qnew = F%Q()
+            do i = 1, m
+                print *, (qnew(i, j), j = 1, k)
+            enddo
+            print *
+            print *, "Difference between the two."
+            qnew = qnew - q
+            do i = 1, m
+                print *, (abs(qnew(i, j)), j=1, k)
+            enddo
+            end block
+            return
+        endif
+    end subroutine test_qr_random_tall_matrix_${ri}$
+
+    subroutine test_qr_random_wide_matrix_${ri}$(error)
+        use stdlib_linalg, only: qr_${rt[0]}$${rk}$_type
+        type(error_type), allocatable, intent(out) :: error
+        integer(ilp), parameter :: m   = 4_ilp
+        integer(ilp), parameter :: n   = 15_ilp
+        integer(ilp), parameter :: k   = min(m,n)
+        ${rt}$ :: a(m, n), q(m, k), r(k, n)
+        real(${rk}$) :: rea(m, n),ima(m, n)
+        type(linalg_state_type) :: state
+        type(qr_${rt[0]}$${rk}$_type) :: F
+
+        ! Random matrix with unit 2-norm.
+        call random_number(rea)
+        #:if rt.startswith('complex')
+        call random_number(ima)
+        a = cmplx(rea,ima,kind=${rk}$)
+        #:else
+        a = rea
+        #:endif
+        a = a / mnorm(a, 2)
+
+        ! Reference QR decomposition.
+        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)
+        call check(error,state%ok(),state%print())
+        if (allocated(error)) return
+
+        ! Instantiate derived-type for the QR.
+        F = qrfact(A)
+
+        ! Check the R matrix with reference.
+        call check(error, all(abs(R - F%R()) < rel_tol_${rk}$))
+        if (allocated(error)) then
+            block
+            integer(ilp) :: i, j
+            ${rt}$, allocatable :: rnew(:, :)
+            print *, "Reference R matrix. shape(R) =", shape(R)
+            do i = 1, k
+                print *, (R(i, j), j = 1, n)
+            enddo
+            print *
+            print *, "New R matrix. shape(R) =", shape(F%R())
+            rnew = F%R()
+            do i = 1, k
+                print *, (rnew(i, j), j = 1, n)
+            enddo
+            print *
+            print *, "Difference between the two."
+            rnew = rnew - r
+            do i = 1, k
+                print *, (abs(rnew(i, j)), j=1, n)
+            enddo
+            end block
+            return
+        endif
+
+        ! Check the Q matrix with reference.
+        call check(error, all(abs(Q - F%Q()) < rel_tol_${rk}$))
+        if (allocated(error)) then
+            block
+            integer(ilp) :: i, j
+            ${rt}$, allocatable :: qnew(:, :)
+            print *, "Reference Q matrix. shape(Q) =", shape(Q)
+            do i = 1, m
+                print *, (Q(i, j), j = 1, k)
+            enddo
+            print *
+            print *, "New Q matrix. shape(Q) =", shape(F%Q())
+            qnew = F%Q()
+            do i = 1, m
+                print *, (qnew(i, j), j = 1, k)
+            enddo
+            print *
+            print *, "Difference between the two."
+            qnew = qnew - q
+            do i = 1, m
+                print *, (abs(qnew(i, j)), j=1, k)
+            enddo
+            end block
+            return
+        endif
+    end subroutine test_qr_random_wide_matrix_${ri}$
+    #:endfor
+
+    !> QR factorization of a tall random rank-deficient matrix
+    #:for rk,rt,ri in RC_KINDS_TYPES
+    subroutine test_qr_random_tall_rank_deficient_matrix_${ri}$(error)
+        use stdlib_linalg, only: qr_${rt[0]}$${rk}$_type
+        type(error_type), allocatable, intent(out) :: error
+        integer(ilp), parameter :: m   = 15_ilp
+        integer(ilp), parameter :: n   =  4_ilp
+        integer(ilp), parameter :: k   = min(m,n)
+        ${rt}$ :: a(m,n), q(m, k), r(k, n)
+        real(${rk}$) :: rea(m,n),ima(m,n)
+        type(linalg_state_type) :: state
+        type(qr_${rt[0]}$${rk}$_type) :: F
+
+        ! Random matrix with unit 2-norm and redundant column.
+        call random_number(rea)
+        #:if rt.startswith('complex')
+        call random_number(ima)
+        a = cmplx(rea,ima,kind=${rk}$)
+        #:else
+        a = rea
+        #:endif
+        a(:, 3) = a(:, 2)
+        a = a / mnorm(a, 2)
+
+        ! Instantiate derived-type for the QR.
+        F = qrfact(A)
+
+        ! Check the factorized matrix.
+        call check(error, mnorm(A - matmul(F%Q(), F%R()), 2) < rel_tol_${rk}$)
+        if (allocated(error)) then
+            block
+            integer(ilp) :: i, j
+            ${rt}$, allocatable :: anew(:, :)
+            anew = matmul(F%q(), F%r())
+            anew = anew - a
+            do i = 1, m
+                print *, (abs(anew(i, j)), j=1, k)
+            enddo
+            end block
+            return
+        endif
+    end subroutine test_qr_random_tall_rank_deficient_matrix_${ri}$
+    #:endfor
+
+
+    ! gcc-15 bugfix utility
+    subroutine add_test(tests,new_test)
+        type(unittest_type), allocatable, intent(inout) :: tests(:)
+        type(unittest_type), intent(in) :: new_test
+
+        integer :: n
+        type(unittest_type), allocatable :: new_tests(:)
+
+        if (allocated(tests)) then
+            n = size(tests)
+        else
+            n = 0
+        end if
+
+        allocate(new_tests(n+1))
+        if (n>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_matrix_factorizations
+
+program test_matrix_factorizations
+    use, intrinsic :: iso_fortran_env, only : error_unit
+    use testdrive, only : run_testsuite, new_testsuite, testsuite_type
+    use test_linalg_matrix_factorizations, only : test_qr_factorization
+    implicit none
+    integer :: stat, is
+    type(testsuite_type), allocatable :: testsuites(:)
+    character(len=*), parameter :: fmt = '("#", *(1x, a))'
+
+    stat = 0
+
+    testsuites = [ &
+        new_testsuite("linalg_matrix_factorizations: qr", test_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_matrix_factorizations