igraph/linalg/
lapack.rs

1//! The LAPACK interface (`igraph_lapack.h`).
2//!
3//! Invalid arguments make the LAPACK bundled with igraph terminate the
4//! process, so all the dimensions and ranges are validated here first.
5
6use super::{Complex, check_finite, to_c_int, unpack_lapack_vectors};
7use crate::{
8    error::{Error, ErrorKind, Result},
9    ffi::*,
10    igraph_call,
11    matrix::Matrix,
12    vector::{Vector, VectorInt},
13};
14use std::{ffi::c_int, ops::Range, ptr};
15
16crate::ffi_enum! {
17    /// Balancing performed by [`lapack_dgeevx`] (`igraph_lapack_dgeevx_balance_t`).
18    pub enum DgeevxBalance: igraph_lapack_dgeevx_balance_t {
19        /// Neither permute nor scale.
20        None = igraph_lapack_dgeevx_balance_t_IGRAPH_LAPACK_DGEEVX_BALANCE_NONE,
21        /// Permute rows and columns to make the matrix more nearly upper
22        /// triangular; do not scale.
23        Perm = igraph_lapack_dgeevx_balance_t_IGRAPH_LAPACK_DGEEVX_BALANCE_PERM,
24        /// Diagonally scale (`D A D^-1`) to make rows and columns closer in
25        /// norm; do not permute.
26        Scale = igraph_lapack_dgeevx_balance_t_IGRAPH_LAPACK_DGEEVX_BALANCE_SCALE,
27        /// Both permute and scale.
28        Both = igraph_lapack_dgeevx_balance_t_IGRAPH_LAPACK_DGEEVX_BALANCE_BOTH,
29    }
30}
31
32fn square(a: &Matrix, what: &str) -> Result<c_int> {
33    if a.nrow() != a.ncol() {
34        return Err(Error::invalid(format!(
35            "{what} needs a square matrix, got {:?}",
36            a.shape()
37        )));
38    }
39    if a.nrow() == 0 {
40        return Err(Error::invalid(format!("{what} needs a non-empty matrix")));
41    }
42    check_finite(a.as_slice(), "the matrix")?;
43    to_c_int(a.nrow(), "matrix order")
44}
45
46/// An LU factorization `A = P L U` computed by [`lapack_dgetrf`].
47#[derive(Debug, Clone, PartialEq)]
48pub struct LuFactors {
49    /// `L` (below the diagonal, unit diagonal not stored) and `U` (on and
50    /// above the diagonal), packed in one matrix of the shape of `A`.
51    pub lu: Matrix,
52    /// Pivot indices, **one-based** as in LAPACK: row `i` was interchanged
53    /// with row `ipiv[i] - 1`.
54    pub ipiv: Vec<i64>,
55    /// LAPACK's `info`: 0 on success, `i > 0` if `U(i-1, i-1)` is exactly
56    /// zero (the factorization is complete, but `U` is singular).
57    pub info: i32,
58}
59
60impl LuFactors {
61    /// Whether `U` is exactly singular.
62    pub fn is_singular(&self) -> bool {
63        self.info > 0
64    }
65
66    /// Solves `A X = B` (or `A' X = B` if `transpose`) with these factors,
67    /// see [`lapack_dgetrs`].
68    pub fn solve(&self, transpose: bool, b: &Matrix) -> Result<Matrix> {
69        lapack_dgetrs(transpose, &self.lu, &self.ipiv, b)
70    }
71}
72
73/// LU factorization with partial pivoting of a general `m` × `n` matrix,
74/// `A = P L U` (`igraph_lapack_dgetrf`), where `L` is lower triangular
75/// (trapezoidal if `m > n`) with unit diagonal and `U` is upper triangular
76/// (trapezoidal if `m < n`).
77///
78/// A singular matrix is not an error: check [`LuFactors::info`].
79///
80/// Binds [`igraph_lapack_dgetrf`](https://igraph.org/c/html/latest/igraph-Linalg.html#igraph_lapack_dgetrf).
81pub fn lapack_dgetrf(a: &Matrix) -> Result<LuFactors> {
82    to_c_int(a.nrow(), "number of rows")?;
83    to_c_int(a.ncol(), "number of columns")?;
84    let mut lu = a.clone();
85    let mut ipiv = VectorInt::new();
86    let mut info: c_int = 0;
87    igraph_call!(igraph_lapack_dgetrf(&mut lu, &mut ipiv, &mut info))?;
88    Ok(LuFactors {
89        lu,
90        ipiv: ipiv.into(),
91        info,
92    })
93}
94
95/// Solves `A X = B` or `A' X = B` (`transpose`) using the LU factors of a
96/// square `A` computed by [`lapack_dgetrf`] (`igraph_lapack_dgetrs`). No
97/// check is made for singularity.
98///
99/// Binds [`igraph_lapack_dgetrs`](https://igraph.org/c/html/latest/igraph-Linalg.html#igraph_lapack_dgetrs).
100///
101/// # Errors
102/// If `lu` is not square, `b` has the wrong number of rows, or the pivots
103/// are out of range.
104pub fn lapack_dgetrs(transpose: bool, lu: &Matrix, ipiv: &[i64], b: &Matrix) -> Result<Matrix> {
105    if lu.nrow() != lu.ncol() {
106        return Err(Error::invalid("dgetrs needs a square LU matrix"));
107    }
108    let n = lu.nrow();
109    to_c_int(n, "matrix order")?;
110    to_c_int(b.ncol(), "number of right hand sides")?;
111    if b.nrow() != n {
112        return Err(Error::invalid(
113            "the right hand side has the wrong number of rows",
114        ));
115    }
116    if ipiv.len() != n || ipiv.iter().any(|&p| p < 1 || p as usize > n) {
117        return Err(Error::invalid("invalid pivot vector"));
118    }
119    let piv = VectorInt::view(ipiv);
120    let mut x = b.clone();
121    igraph_call!(igraph_lapack_dgetrs(transpose, lu, piv.as_ptr(), &mut x))?;
122    Ok(x)
123}
124
125/// The solution of a linear system computed by [`lapack_dgesv`].
126#[derive(Debug, Clone, PartialEq)]
127pub struct DgesvResult {
128    /// The solution `X` of `A X = B` (same shape as `B`).
129    pub solution: Matrix,
130    /// The LU factors of `A`, see [`LuFactors`].
131    pub factors: LuFactors,
132}
133
134/// Solves the linear system `A X = B` for a square `A` and one or more
135/// right hand sides (the columns of `B`), by LU decomposition with partial
136/// pivoting (`igraph_lapack_dgesv`).
137///
138/// Binds [`igraph_lapack_dgesv`](https://igraph.org/c/html/latest/igraph-Linalg.html#igraph_lapack_dgesv).
139///
140/// # Errors
141/// If the dimensions do not match, or `A` is exactly singular
142/// ([`ErrorKind::Failure`](crate::ErrorKind::Failure)).
143///
144/// # Examples
145///
146/// ```
147/// use igraph::{linalg::lapack_dgesv, prelude::*};
148/// // 2x + y = 5, x + 3y = 10  =>  x = 1, y = 3
149/// let a = Matrix::from_rows(&[[2.0, 1.0], [1.0, 3.0]]).unwrap();
150/// let b = Matrix::from_rows(&[[5.0], [10.0]]).unwrap();
151/// let x = lapack_dgesv(&a, &b).unwrap().solution;
152/// assert!((x[(0, 0)] - 1.0).abs() < 1e-12 && (x[(1, 0)] - 3.0).abs() < 1e-12);
153/// ```
154pub fn lapack_dgesv(a: &Matrix, b: &Matrix) -> Result<DgesvResult> {
155    if a.nrow() != a.ncol() {
156        return Err(Error::invalid("dgesv needs a square coefficient matrix"));
157    }
158    to_c_int(a.nrow(), "matrix order")?;
159    to_c_int(b.ncol(), "number of right hand sides")?;
160    if b.nrow() != a.nrow() {
161        return Err(Error::invalid(
162            "the right hand side has the wrong number of rows",
163        ));
164    }
165    let mut lu = a.clone();
166    let mut x = b.clone();
167    let mut ipiv = VectorInt::new();
168    let mut info: c_int = 0;
169    igraph_call!(igraph_lapack_dgesv(&mut lu, &mut ipiv, &mut x, &mut info))?;
170    if info > 0 {
171        return Err(Error::new(
172            ErrorKind::Failure,
173            format!("the matrix is exactly singular (U({0},{0}) = 0)", info - 1),
174        ));
175    }
176    Ok(DgesvResult {
177        solution: x,
178        factors: LuFactors {
179            lu,
180            ipiv: ipiv.into(),
181            info,
182        },
183    })
184}
185
186/// Solves `A x = b` for a single right hand side (via [`lapack_dgesv`]).
187///
188/// ```
189/// use igraph::{linalg::solve, prelude::*};
190/// let a = Matrix::from_rows(&[[4.0, -2.0], [1.0, 1.0]]).unwrap();
191/// let x = solve(&a, &[2.0, 3.0]).unwrap();
192/// assert!((x[0] - 4.0 / 3.0).abs() < 1e-12 && (x[1] - 5.0 / 3.0).abs() < 1e-12);
193/// ```
194pub fn solve(a: &Matrix, b: &[f64]) -> Result<Vec<f64>> {
195    let bm = Matrix::from_column_major(b.len(), 1, b)?;
196    Ok(lapack_dgesv(a, &bm)?.solution.as_slice().to_vec())
197}
198
199/// Which eigenvalues [`lapack_dsyevr`] computes (`igraph_lapack_dsyev_which_t`
200/// plus its parameters).
201#[derive(Debug, Clone, PartialEq)]
202pub enum SymmetricRange {
203    /// All the eigenvalues (`IGRAPH_LAPACK_DSYEV_ALL`).
204    All,
205    /// The eigenvalues in the half-open interval `(low, high]`
206    /// (`IGRAPH_LAPACK_DSYEV_INTERVAL`).
207    Interval {
208        /// Exclusive lower bound.
209        low: f64,
210        /// Inclusive upper bound.
211        high: f64,
212    },
213    /// The eigenvalues with the given zero-based positions in increasing
214    /// order (`IGRAPH_LAPACK_DSYEV_SELECT` with one-based `il = start + 1`,
215    /// `iu = end`).
216    Select(Range<usize>),
217}
218
219/// Eigenvalues and eigenvectors computed by [`lapack_dsyevr`].
220#[derive(Debug, Clone, PartialEq)]
221pub struct DsyevrResult {
222    /// The selected eigenvalues, in increasing order.
223    pub values: Vec<f64>,
224    /// The orthonormal eigenvectors, in the columns.
225    pub vectors: Matrix,
226}
227
228/// Selected eigenvalues and eigenvectors of a real symmetric matrix, with
229/// LAPACK's relatively robust representations algorithm
230/// (`igraph_lapack_dsyevr`). Only the upper triangle of `a` is used.
231///
232/// `abstol` is the absolute error tolerance for the eigenvalues: an
233/// approximate eigenvalue is accepted when it lies in an interval `[a, b]`
234/// of width at most `abstol + eps * max(|a|, |b|)`.
235///
236/// Binds [`igraph_lapack_dsyevr`](https://igraph.org/c/html/latest/igraph-Linalg.html#igraph_lapack_dsyevr).
237/// The eigenvector support output of the C function is not exposed.
238///
239/// # Examples
240///
241/// This is igraph's `igraph_lapack_dsyevr.c` example: the matrix
242/// `[[2, -1], [-1, 3]]` has eigenvalues `(5 ± √5) / 2`.
243///
244/// ```
245/// use igraph::{linalg::*, prelude::*};
246/// let a = Matrix::from_rows(&[[2.0, -1.0], [-1.0, 3.0]]).unwrap();
247/// let low = lapack_dsyevr(&a, &SymmetricRange::Select(0..1), 1e-10).unwrap();
248/// assert!((low.values[0] - 1.381966).abs() < 1e-6);
249/// let high = lapack_dsyevr(&a, &SymmetricRange::Interval { low: 3.0, high: 4.0 }, 1e-10).unwrap();
250/// assert!((high.values[0] - 3.618034).abs() < 1e-6);
251/// assert_eq!(high.vectors.shape(), (2, 1));
252/// ```
253pub fn lapack_dsyevr(a: &Matrix, range: &SymmetricRange, abstol: f64) -> Result<DsyevrResult> {
254    let n = square(a, "dsyevr")?;
255    let (which, vl, vu, vestimate, il, iu) = match range {
256        SymmetricRange::All => (
257            igraph_lapack_dsyev_which_t_IGRAPH_LAPACK_DSYEV_ALL,
258            0.0,
259            0.0,
260            0,
261            0,
262            0,
263        ),
264        SymmetricRange::Interval { low, high } => {
265            if low.is_nan() || high.is_nan() || low >= high {
266                return Err(Error::invalid(
267                    "the eigenvalue interval must satisfy low < high",
268                ));
269            }
270            // `n` is always a correct upper bound on the number of eigenvalues.
271            (
272                igraph_lapack_dsyev_which_t_IGRAPH_LAPACK_DSYEV_INTERVAL,
273                *low,
274                *high,
275                n,
276                0,
277                0,
278            )
279        }
280        SymmetricRange::Select(r) => {
281            if r.start >= r.end || r.end > n as usize {
282                return Err(Error::invalid(format!(
283                    "invalid eigenvalue range {r:?} for order {n}"
284                )));
285            }
286            (
287                igraph_lapack_dsyev_which_t_IGRAPH_LAPACK_DSYEV_SELECT,
288                0.0,
289                0.0,
290                0,
291                r.start as c_int + 1,
292                r.end as c_int,
293            )
294        }
295    };
296    let mut values = Vector::new();
297    let mut vectors = Matrix::new();
298    igraph_call!(igraph_lapack_dsyevr(
299        a,
300        which,
301        vl,
302        vu,
303        vestimate,
304        il,
305        iu,
306        abstol,
307        &mut values,
308        &mut vectors,
309        ptr::null_mut()
310    ))?;
311    Ok(DsyevrResult {
312        values: values.into(),
313        vectors,
314    })
315}
316
317/// Eigenvalues and eigenvectors of a general real matrix computed by
318/// [`lapack_dgeev`], in LAPACK's packed real format.
319#[derive(Debug, Clone, PartialEq)]
320pub struct DgeevResult {
321    /// Real parts of the eigenvalues.
322    pub values_real: Vec<f64>,
323    /// Imaginary parts of the eigenvalues; complex conjugate pairs appear
324    /// consecutively, the one with positive imaginary part first.
325    pub values_imag: Vec<f64>,
326    /// Left eigenvectors `u` (`u^H A = lambda u^H`), if requested, packed:
327    /// a real eigenvalue `j` owns column `j`; a complex pair `(j, j+1)` has
328    /// its eigenvectors `u_j = col_j + i col_{j+1}` and `u_{j+1} = conj(u_j)`.
329    pub vectors_left: Option<Matrix>,
330    /// Right eigenvectors `v` (`A v = lambda v`), if requested, packed like
331    /// [`vectors_left`](Self::vectors_left).
332    pub vectors_right: Option<Matrix>,
333}
334
335impl DgeevResult {
336    /// The eigenvalues as complex numbers.
337    pub fn values(&self) -> Vec<Complex> {
338        self.values_real
339            .iter()
340            .zip(&self.values_imag)
341            .map(|(&re, &im)| Complex::new(re, im))
342            .collect()
343    }
344
345    /// The unpacked right eigenvectors, `result[k]` belonging to eigenvalue `k`.
346    pub fn right_eigenvectors(&self) -> Option<Vec<Vec<Complex>>> {
347        self.vectors_right
348            .as_ref()
349            .map(|m| unpack_lapack_vectors(&self.values_imag, m))
350    }
351
352    /// The unpacked left eigenvectors, `result[k]` belonging to eigenvalue `k`.
353    pub fn left_eigenvectors(&self) -> Option<Vec<Vec<Complex>>> {
354        self.vectors_left
355            .as_ref()
356            .map(|m| unpack_lapack_vectors(&self.values_imag, m))
357    }
358}
359
360/// Eigenvalues and, optionally, left and/or right eigenvectors of a general
361/// real square matrix (`igraph_lapack_dgeev`). The eigenvectors are
362/// normalized to unit Euclidean norm with largest component real.
363///
364/// Binds [`igraph_lapack_dgeev`](https://igraph.org/c/html/latest/igraph-Linalg.html#igraph_lapack_dgeev).
365///
366/// # Errors
367/// If the matrix is not square or empty, or the QR algorithm fails.
368///
369/// # Examples
370///
371/// igraph's `igraph_lapack_dgeev.c` example: `[[1, 1], [-1, 1]]` has
372/// eigenvalues `1 ± i`.
373///
374/// ```
375/// use igraph::{linalg::lapack_dgeev, prelude::*};
376/// let a = Matrix::from_rows(&[[1.0, 1.0], [-1.0, 1.0]]).unwrap();
377/// let e = lapack_dgeev(&a, true, true).unwrap();
378/// let values = e.values();
379/// assert!((values[0].re() - 1.0).abs() < 1e-12 && (values[0].im() - 1.0).abs() < 1e-12);
380/// assert!((values[1].re() - 1.0).abs() < 1e-12 && (values[1].im() + 1.0).abs() < 1e-12);
381/// // Check A v = (1 + i) v on the first component: (A v)_0 = v_0 + v_1.
382/// let v = &e.right_eigenvectors().unwrap()[0];
383/// let (lhs_re, lhs_im) = (v[0].re() + v[1].re(), v[0].im() + v[1].im());
384/// let (rhs_re, rhs_im) = (v[0].re() - v[0].im(), v[0].re() + v[0].im());
385/// assert!((lhs_re - rhs_re).abs() < 1e-12 && (lhs_im - rhs_im).abs() < 1e-12);
386/// ```
387pub fn lapack_dgeev(a: &Matrix, left: bool, right: bool) -> Result<DgeevResult> {
388    square(a, "dgeev")?;
389    let (mut re, mut im) = (Vector::new(), Vector::new());
390    let mut vl = Matrix::new();
391    let mut vr = Matrix::new();
392    let vlp = if left {
393        &mut vl as *mut Matrix
394    } else {
395        ptr::null_mut()
396    };
397    let vrp = if right {
398        &mut vr as *mut Matrix
399    } else {
400        ptr::null_mut()
401    };
402    // Non-zero on entry: report failures of the QR algorithm as errors.
403    let mut info: c_int = 1;
404    igraph_call!(igraph_lapack_dgeev(
405        a, &mut re, &mut im, vlp, vrp, &mut info
406    ))?;
407    Ok(DgeevResult {
408        values_real: re.into(),
409        values_imag: im.into(),
410        vectors_left: left.then_some(vl),
411        vectors_right: right.then_some(vr),
412    })
413}
414
415/// Output of [`lapack_dgeevx`].
416#[derive(Debug, Clone, PartialEq)]
417pub struct DgeevxResult {
418    /// Eigenvalues and (packed) left and right eigenvectors, see [`DgeevResult`].
419    pub eigen: DgeevResult,
420    /// `ilo` of the balancing (one-based): the balanced `A(i, j) = 0` if
421    /// `i > j` and `j < ilo` or `i > ihi`.
422    pub ilo: i32,
423    /// `ihi` of the balancing (one-based).
424    pub ihi: i32,
425    /// Details of the permutations and scaling factors of the balancing.
426    pub scale: Vec<f64>,
427    /// The one-norm of the balanced matrix (maximum absolute column sum).
428    pub abnrm: f64,
429    /// Reciprocal condition numbers of the eigenvalues.
430    pub rconde: Vec<f64>,
431}
432
433/// Eigenvalues and left and right eigenvectors of a general real matrix,
434/// "expert" version (`igraph_lapack_dgeevx`): it also balances the matrix
435/// (see [`DgeevxBalance`]) and computes reciprocal condition numbers of the
436/// eigenvalues.
437///
438/// Binds [`igraph_lapack_dgeevx`](https://igraph.org/c/html/latest/igraph-Linalg.html#igraph_lapack_dgeevx).
439/// The reciprocal condition numbers of the eigenvectors (`rcondv`) are not
440/// computed: igraph 1.0.0 and 1.0.1 allocate too small an integer workspace for them.
441///
442/// ```
443/// use igraph::{linalg::*, prelude::*};
444/// let a = Matrix::from_rows(&[[1.0, 1e4], [1e-4, 1.0]]).unwrap();
445/// let r = lapack_dgeevx(DgeevxBalance::Both, &a).unwrap();
446/// let mut re = r.eigen.values_real.clone();
447/// re.sort_by(f64::total_cmp);
448/// assert!((re[0] - 0.0).abs() < 1e-9 && (re[1] - 2.0).abs() < 1e-9);
449/// assert!(r.rconde.iter().all(|&c| c > 0.0 && c <= 1.0 + 1e-12));
450/// ```
451pub fn lapack_dgeevx(balance: DgeevxBalance, a: &Matrix) -> Result<DgeevxResult> {
452    let n = square(a, "dgeevx")? as usize;
453    let (mut re, mut im) = (Vector::new(), Vector::new());
454    let (mut vl, mut vr) = (Matrix::new(), Matrix::new());
455    let (mut ilo, mut ihi): (c_int, c_int) = (0, 0);
456    let mut scale = Vector::new();
457    let mut abnrm = 0.0;
458    // igraph does not resize `rconde` before LAPACK writes `n` values into it.
459    let mut rconde = Vector::zeros(n);
460    let mut info: c_int = 1;
461    igraph_call!(igraph_lapack_dgeevx(
462        balance.into(),
463        a,
464        &mut re,
465        &mut im,
466        &mut vl,
467        &mut vr,
468        &mut ilo,
469        &mut ihi,
470        &mut scale,
471        &mut abnrm,
472        &mut rconde,
473        ptr::null_mut(),
474        &mut info
475    ))?;
476    Ok(DgeevxResult {
477        eigen: DgeevResult {
478            values_real: re.into(),
479            values_imag: im.into(),
480            vectors_left: Some(vl),
481            vectors_right: Some(vr),
482        },
483        ilo,
484        ihi,
485        scale: scale.into(),
486        abnrm,
487        rconde: rconde.into(),
488    })
489}
490
491/// Reduces a general square matrix to upper Hessenberg form by an
492/// orthogonal similarity transformation (`igraph_lapack_dgehrd`): the
493/// result `H = Q' A Q` has zeros below the first subdiagonal and the same
494/// eigenvalues as `A`.
495///
496/// `ilo` and `ihi` are one-based, `1 <= ilo <= ihi <= n`; use `1` and `n`
497/// unless `A` is already upper triangular in rows and columns outside
498/// `ilo..=ihi` (e.g. after balancing with [`lapack_dgeevx`]).
499///
500/// Binds [`igraph_lapack_dgehrd`](https://igraph.org/c/html/latest/igraph-Linalg.html#igraph_lapack_dgehrd).
501///
502/// ```
503/// use igraph::{linalg::lapack_dgehrd, prelude::*};
504/// let a = Matrix::from_rows(&[[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 10.0]]).unwrap();
505/// let h = lapack_dgehrd(&a, 1, 3).unwrap();
506/// assert_eq!(h[(2, 0)], 0.0);
507/// // The trace is preserved by similarity transformations.
508/// assert!((h[(0, 0)] + h[(1, 1)] + h[(2, 2)] - 16.0).abs() < 1e-10);
509/// ```
510pub fn lapack_dgehrd(a: &Matrix, ilo: usize, ihi: usize) -> Result<Matrix> {
511    let n = square(a, "dgehrd")? as usize;
512    if ilo < 1 || ihi > n || ilo > ihi {
513        return Err(Error::invalid(format!(
514            "invalid ilo = {ilo} and ihi = {ihi} for order {n}"
515        )));
516    }
517    let mut res = Matrix::new();
518    igraph_call!(igraph_lapack_dgehrd(
519        a,
520        ilo as c_int,
521        ihi as c_int,
522        &mut res
523    ))?;
524    Ok(res)
525}
igraph/linalg/lapack.rs

igraph/linalg/
lapack.rs
