Skip to main content

mdarray_linalg_lapack/qr/
context.rs

1//! QR Decomposition:
2//!     A = Q * R
3//! where:
4//!     - A is m × n (input matrix)
5//!     - Q is m × k (orthogonal matrix)
6//!     - R is k × n (upper triangular matrix)
7//!
8//! `qr()` returns the reduced factorization where k = min(m, n):
9//! Q is m × k and R is k × n.  `qr_write()` infers the requested mode
10//! from the shapes of the provided output matrices, so callers can request
11//! either reduced or complete QR by choosing output shapes.
12//!
13//! The implementation wraps two LAPACK routines:
14//! - `geqrf`: Computes the QR factorization in a compact (implicit) form using
15//!   blocked Householder reflectors.
16//! - `orgqr` / `ungqr`: Explicitly generates the orthogonal matrix Q from the
17//!   reflectors returned by `geqrf`, producing either the thin or square Q
18//!   depending on the output shape.
19
20use mdarray::{Array, Dim, Layout, Shape, Slice};
21use mdarray_linalg::qr::QR;
22use num_complex::ComplexFloat;
23
24use super::{
25    scalar::{LapackScalar, NeedsRwork},
26    simple::geqrf,
27};
28use crate::{Lapack, QRConfig};
29
30impl<T, D0: Dim, D1: Dim> QR<T, D0, D1> for Lapack
31where
32    T: ComplexFloat + Default + LapackScalar + NeedsRwork,
33    T::Real: Into<T>,
34{
35    fn qr_write<D2: Dim, L: Layout, Lq: Layout, Lr: Layout>(
36        &self,
37        a: &mut Slice<T, (D0, D1), L>,
38        q: &mut Slice<T, (D0, D2), Lq>,
39        r: &mut Slice<T, (D2, D1), Lr>,
40    ) {
41        geqrf(a, q, r, self.qr_config)
42    }
43
44    fn qr<L: Layout>(
45        &self,
46        a: &mut Slice<T, (D0, D1), L>,
47    ) -> (Array<T, (D0, usize)>, Array<T, (usize, D1)>) {
48        let ash = *a.shape();
49        let m = ash.dim(0);
50        let n = ash.dim(1);
51        let min_mn = m.min(n);
52
53        let (q_rows, q_cols, r_rows, r_cols) = match self.qr_config {
54            QRConfig::Reduced => (m, min_mn, min_mn, n),
55            QRConfig::Complete => (m, m, m, n),
56        };
57
58        let q_shape = <(D0, usize) as Shape>::from_dims(&[q_rows, q_cols]);
59        let r_shape = <(usize, D1) as Shape>::from_dims(&[r_rows, r_cols]);
60
61        let mut q: Array<T, (D0, usize)> = Array::from_elem(q_shape, T::default());
62        let mut r: Array<T, (usize, D1)> = Array::from_elem(r_shape, T::default());
63
64        geqrf(a, &mut q, &mut r, self.qr_config);
65
66        (q, r)
67    }
68}