igraph/linalg/
eigen.rs

1//! Generic eigensolvers (`igraph_eigen.h`).
2
3use super::{
4    ArpackOptions, Complex, SparseMat,
5    arpack::{matvec_trampoline, with_finite_guard},
6    check_finite,
7    lapack::DgeevxBalance,
8    to_c_int,
9};
10use crate::{
11    error::{Error, Result},
12    ffi::*,
13    igraph_call,
14    matrix::{Matrix, MatrixComplex},
15    vector::{Vector, VectorComplex},
16};
17use std::{ffi::c_void, ops::Range, ptr};
18
19crate::ffi_enum! {
20    /// Which numerical library [`eigen_matrix_symmetric`], [`eigen_matrix`]
21    /// and [`Graph::eigen_adjacency`](crate::Graph::eigen_adjacency) use (`igraph_eigen_algorithm_t`).
22    pub enum EigenAlgorithm: igraph_eigen_algorithm_t {
23        /// Choose automatically. For symmetric matrices: LAPACK when the matrix
24        /// is small (`n < 100`), when as many eigenvalues as the order are
25        /// requested, and for [`EigenWhich::All`], [`EigenWhich::Interval`]
26        /// and [`EigenWhich::Select`] (which only LAPACK supports); ARPACK
27        /// otherwise. ARPACK for [`Graph::eigen_adjacency`](crate::Graph::eigen_adjacency).
28        /// Not implemented by igraph (1.0.0 and 1.0.1) for [`eigen_matrix`].
29        Auto = igraph_eigen_algorithm_t_IGRAPH_EIGEN_AUTO,
30        /// Dense LAPACK routines (the matrix is formed explicitly).
31        Lapack = igraph_eigen_algorithm_t_IGRAPH_EIGEN_LAPACK,
32        /// ARPACK (matrix-free, only matrix-vector products).
33        Arpack = igraph_eigen_algorithm_t_IGRAPH_EIGEN_ARPACK,
34        /// Complex variant of `Auto` (not implemented by igraph 1.0.0 and 1.0.1).
35        CompAuto = igraph_eigen_algorithm_t_IGRAPH_EIGEN_COMP_AUTO,
36        /// Complex variant of `Lapack` (not implemented by igraph 1.0.0 and 1.0.1).
37        CompLapack = igraph_eigen_algorithm_t_IGRAPH_EIGEN_COMP_LAPACK,
38        /// Complex variant of `Arpack` (not implemented by igraph 1.0.0 and 1.0.1).
39        CompArpack = igraph_eigen_algorithm_t_IGRAPH_EIGEN_COMP_ARPACK,
40    }
41}
42
43/// Which eigenvalues to compute (`igraph_eigen_which_t`).
44///
45/// Counts must be between 1 and the order of the matrix. Not every choice is
46/// supported by every solver: the symmetric solvers accept the magnitude,
47/// algebraic, [`BothEnds`](Self::BothEnds), [`All`](Self::All),
48/// [`Interval`](Self::Interval) and [`Select`](Self::Select) choices (the
49/// last three only with LAPACK); the non-symmetric solver accepts the
50/// magnitude, real part, imaginary part, [`All`](Self::All) and
51/// [`Select`](Self::Select) choices.
52#[derive(Debug, Clone, PartialEq)]
53pub enum EigenWhich {
54    /// The given number of eigenvalues of largest magnitude.
55    LargestMagnitude(usize),
56    /// The given number of eigenvalues of smallest magnitude.
57    SmallestMagnitude(usize),
58    /// The given number of largest (algebraic) eigenvalues; symmetric only.
59    LargestAlgebraic(usize),
60    /// The given number of smallest (algebraic) eigenvalues; symmetric only.
61    SmallestAlgebraic(usize),
62    /// The given number of eigenvalues, alternating from the two ends of the
63    /// spectrum (largest first); symmetric only. LAPACK needs at least 2.
64    BothEnds(usize),
65    /// The given number of eigenvalues with largest real part; non-symmetric only.
66    LargestReal(usize),
67    /// The given number of eigenvalues with smallest real part; non-symmetric only.
68    SmallestReal(usize),
69    /// The given number of eigenvalues with largest imaginary part; non-symmetric only.
70    LargestImaginary(usize),
71    /// The given number of eigenvalues with smallest imaginary part; non-symmetric only.
72    SmallestImaginary(usize),
73    /// All the eigenvalues.
74    All,
75    /// All the eigenvalues in the half-open interval `(low, high]`;
76    /// symmetric LAPACK only.
77    Interval {
78        /// Exclusive lower bound.
79        low: f64,
80        /// Inclusive upper bound.
81        high: f64,
82    },
83    /// Eigenvalues by (zero-based) position: for symmetric matrices in
84    /// increasing algebraic order, for non-symmetric ones in increasing
85    /// magnitude. LAPACK only.
86    Select(Range<usize>),
87}
88
89impl EigenWhich {
90    pub(crate) fn to_raw(&self, n: usize) -> Result<igraph_eigen_which_t> {
91        use EigenWhich::*;
92        let mut raw = igraph_eigen_which_t {
93            pos: igraph_eigen_which_position_t_IGRAPH_EIGEN_ALL,
94            howmany: 0,
95            il: 0,
96            iu: 0,
97            vl: 0.0,
98            vu: 0.0,
99            vestimate: 0,
100            balance: DgeevxBalance::None.into(),
101        };
102        let count = |k: usize| -> Result<std::ffi::c_int> {
103            if k == 0 || k > n {
104                return Err(Error::invalid(format!(
105                    "cannot compute {k} eigenvalues of a matrix of order {n}"
106                )));
107            }
108            to_c_int(k, "number of eigenvalues")
109        };
110        let (pos, howmany) = match self {
111            LargestMagnitude(k) => (igraph_eigen_which_position_t_IGRAPH_EIGEN_LM, count(*k)?),
112            SmallestMagnitude(k) => (igraph_eigen_which_position_t_IGRAPH_EIGEN_SM, count(*k)?),
113            LargestAlgebraic(k) => (igraph_eigen_which_position_t_IGRAPH_EIGEN_LA, count(*k)?),
114            SmallestAlgebraic(k) => (igraph_eigen_which_position_t_IGRAPH_EIGEN_SA, count(*k)?),
115            BothEnds(k) => (igraph_eigen_which_position_t_IGRAPH_EIGEN_BE, count(*k)?),
116            LargestReal(k) => (igraph_eigen_which_position_t_IGRAPH_EIGEN_LR, count(*k)?),
117            SmallestReal(k) => (igraph_eigen_which_position_t_IGRAPH_EIGEN_SR, count(*k)?),
118            LargestImaginary(k) => (igraph_eigen_which_position_t_IGRAPH_EIGEN_LI, count(*k)?),
119            SmallestImaginary(k) => (igraph_eigen_which_position_t_IGRAPH_EIGEN_SI, count(*k)?),
120            All => (igraph_eigen_which_position_t_IGRAPH_EIGEN_ALL, 0),
121            Interval { low, high } => {
122                if low.is_nan() || high.is_nan() || low >= high {
123                    return Err(Error::invalid(
124                        "the eigenvalue interval must satisfy low < high",
125                    ));
126                }
127                raw.vl = *low;
128                raw.vu = *high;
129                // An upper bound on the number of eigenvalues found: always correct.
130                raw.vestimate = to_c_int(n, "matrix order")?;
131                (igraph_eigen_which_position_t_IGRAPH_EIGEN_INTERVAL, 0)
132            }
133            Select(r) => {
134                if r.start >= r.end || r.end > n {
135                    return Err(Error::invalid(format!(
136                        "invalid eigenvalue range {r:?} for order {n}"
137                    )));
138                }
139                raw.il = to_c_int(r.start + 1, "range start")?;
140                raw.iu = to_c_int(r.end, "range end")?;
141                (igraph_eigen_which_position_t_IGRAPH_EIGEN_SELECT, 0)
142            }
143        };
144        raw.pos = pos;
145        raw.howmany = howmany;
146        Ok(raw)
147    }
148}
149
150/// Eigenvalues and eigenvectors of a real symmetric matrix.
151#[derive(Debug, Clone, PartialEq)]
152pub struct SymmetricEigen {
153    /// The eigenvalues, in the order implied by the [`EigenWhich`] choice:
154    /// decreasing magnitude for `LargestMagnitude`, increasing magnitude for
155    /// `SmallestMagnitude`, decreasing for `LargestAlgebraic`, increasing for
156    /// `SmallestAlgebraic`, `All`, `Interval` and `Select`, and alternating
157    /// largest / smallest for `BothEnds` — whichever algorithm is used.
158    pub values: Vec<f64>,
159    /// The unit eigenvectors, in the columns (`n` × number of eigenvalues).
160    pub vectors: Matrix,
161}
162
163/// Eigenvalues and eigenvectors of a general real matrix.
164#[derive(Debug, Clone, PartialEq)]
165pub struct ComplexEigen {
166    /// The (possibly complex) eigenvalues.
167    pub values: Vec<Complex>,
168    /// The eigenvectors: `vectors[k]` belongs to `values[k]`.
169    pub vectors: Vec<Vec<Complex>>,
170}
171
172enum Operator<'a> {
173    Dense(&'a Matrix),
174    Sparse(&'a SparseMat),
175    Fn(igraph_arpack_function_t, *mut c_void),
176}
177
178fn symmetric_common(
179    op: Operator<'_>,
180    n: usize,
181    which: &EigenWhich,
182    algorithm: EigenAlgorithm,
183    options: &ArpackOptions,
184) -> Result<SymmetricEigen> {
185    if n == 0 {
186        return Err(Error::invalid(
187            "cannot compute the eigenvalues of an empty matrix",
188        ));
189    }
190    let cn = to_c_int(n, "matrix order")?;
191    // Resolve `Auto` here: igraph (1.0.0 and 1.0.1) compares the (unused,
192    // zero) `howmany` of `All`, `Interval` and `Select` with `n` and so picks
193    // ARPACK for them when `n >= 100`, where they fail. And for `n <= 2`
194    // igraph's "ARPACK" code path is a closed-form shortcut that confuses
195    // magnitude with algebraic order (see `uses_2x2_shortcut`): LAPACK is
196    // exact there.
197    let algorithm = match (algorithm, which) {
198        (
199            EigenAlgorithm::Auto,
200            EigenWhich::All | EigenWhich::Interval { .. } | EigenWhich::Select(_),
201        ) => EigenAlgorithm::Lapack,
202        (
203            EigenAlgorithm::Auto,
204            EigenWhich::LargestMagnitude(k)
205            | EigenWhich::SmallestMagnitude(k)
206            | EigenWhich::LargestAlgebraic(k)
207            | EigenWhich::SmallestAlgebraic(k)
208            | EigenWhich::BothEnds(k),
209        ) if *k == n || n < 100 => EigenAlgorithm::Lapack,
210        (EigenAlgorithm::Auto, _) => EigenAlgorithm::Arpack,
211        (EigenAlgorithm::Arpack, _) if n <= 2 => EigenAlgorithm::Lapack,
212        (other, _) => other,
213    };
214    // igraph's (1.0.0 and 1.0.1) LAPACK code path for "smallest magnitude" walks outside the
215    // eigenvalue array (and the eigenvector matrix) whenever the eigenvalue
216    // of smallest magnitude is at either end of the spectrum, e.g. for every
217    // positive definite matrix. Compute the whole spectrum with LAPACK
218    // instead and pick the eigenpairs here.
219    let sm_with_lapack = match which {
220        EigenWhich::SmallestMagnitude(_) => {
221            which.to_raw(n)?;
222            algorithm == EigenAlgorithm::Lapack
223        }
224        _ => false,
225    };
226    let raw_which = if sm_with_lapack {
227        EigenWhich::All.to_raw(n)?
228    } else {
229        which.to_raw(n)?
230    };
231    let mut raw_opts = options.to_raw_unchecked_which(n)?;
232    let (a, sa, fun, extra) = match op {
233        Operator::Dense(m) => (m as *const Matrix, ptr::null(), None, ptr::null_mut()),
234        Operator::Sparse(s) => (ptr::null(), s as *const SparseMat, None, ptr::null_mut()),
235        Operator::Fn(f, e) => (ptr::null(), ptr::null(), f, e),
236    };
237    let mut values = Vector::new();
238    let mut vectors = Matrix::new();
239    // Only the ARPACK code path uses ARPACK's non-re-entrant state.
240    let _arpack = if algorithm == EigenAlgorithm::Lapack {
241        None
242    } else {
243        Some(super::arpack::ArpackGuard::enter()?)
244    };
245    igraph_call!(igraph_eigen_matrix_symmetric(
246        a,
247        sa,
248        fun,
249        cn,
250        extra,
251        algorithm.into(),
252        &raw_which,
253        &mut raw_opts,
254        ptr::null_mut(),
255        &mut values,
256        &mut vectors
257    ))?;
258    SymmetricEigen {
259        values: values.into(),
260        vectors,
261    }
262    .ordered(which)
263}
264
265impl SymmetricEigen {
266    /// Puts the eigenpairs in the documented order for `which` (igraph's
267    /// LAPACK and ARPACK code paths do not agree on it), keeping only the
268    /// requested number of them.
269    fn ordered(self, which: &EigenWhich) -> Result<Self> {
270        let key: fn(f64, f64) -> std::cmp::Ordering = match which {
271            EigenWhich::LargestMagnitude(_) => |a, b| b.abs().total_cmp(&a.abs()),
272            EigenWhich::SmallestMagnitude(_) => |a, b| a.abs().total_cmp(&b.abs()),
273            EigenWhich::LargestAlgebraic(_) => |a, b| b.total_cmp(&a),
274            EigenWhich::SmallestAlgebraic(_) => |a, b| a.total_cmp(&b),
275            _ => return Ok(self),
276        };
277        let k = match which {
278            EigenWhich::LargestMagnitude(k)
279            | EigenWhich::SmallestMagnitude(k)
280            | EigenWhich::LargestAlgebraic(k)
281            | EigenWhich::SmallestAlgebraic(k) => (*k).min(self.values.len()),
282            _ => self.values.len(),
283        };
284        let mut order: Vec<usize> = (0..self.values.len()).collect();
285        // Stable: ties keep igraph's order.
286        order.sort_by(|&i, &j| key(self.values[i], self.values[j]));
287        order.truncate(k);
288        let n = self.vectors.nrow();
289        let has_vectors = self.vectors.ncol() == self.values.len();
290        let mut data = Vec::with_capacity(n * k);
291        let mut values = Vec::with_capacity(k);
292        for &i in &order {
293            values.push(self.values[i]);
294            if has_vectors {
295                let v = self.vectors.column(i);
296                let norm = v.iter().map(|x| x * x).sum::<f64>().sqrt();
297                if n <= 2 && norm > 0.0 {
298                    // igraph's closed-form 2 x 2 "ARPACK" shortcut does not
299                    // normalize the eigenvectors.
300                    data.extend(v.iter().map(|x| x / norm));
301                } else {
302                    data.extend_from_slice(v);
303                }
304            }
305        }
306        let vectors = if has_vectors {
307            Matrix::from_column_major(n, k, &data)?
308        } else {
309            self.vectors
310        };
311        Ok(Self { values, vectors })
312    }
313}
314
315/// Selected eigenvalues and eigenvectors of a real **symmetric** dense
316/// matrix (`igraph_eigen_matrix_symmetric`), with LAPACK or ARPACK.
317///
318/// Only the upper triangle is used by LAPACK. `options` tunes ARPACK: only
319/// its `tol` and `mxiter` are used (`which`, `nev` and `ncv` are derived
320/// from `which`, and ARPACK starts from a random vector of the calling
321/// thread's RNG). Matrices of order at most 2 always use LAPACK.
322///
323/// Binds `igraph_eigen_matrix_symmetric` (see the
324/// [linear algebra chapter](https://igraph.org/c/html/latest/igraph-Linalg.html)).
325///
326/// # Errors
327/// If the matrix is not square or empty, the choice is invalid for the
328/// matrix order, or not supported by the algorithm (e.g. `Select` with
329/// ARPACK gives [`ErrorKind::Unimplemented`](crate::ErrorKind::Unimplemented)).
330///
331/// # Examples
332///
333/// ```
334/// use igraph::{linalg::*, prelude::*};
335/// // Eigenvalues of [[2, 1], [1, 2]] are 1 and 3.
336/// let a = Matrix::from_rows(&[[2.0, 1.0], [1.0, 2.0]]).unwrap();
337/// let e = eigen_matrix_symmetric(&a, &EigenWhich::All, EigenAlgorithm::Lapack, &ArpackOptions::default()).unwrap();
338/// let mut values = e.values.clone();
339/// values.sort_by(f64::total_cmp);
340/// assert!((values[0] - 1.0).abs() < 1e-12 && (values[1] - 3.0).abs() < 1e-12);
341/// ```
342pub fn eigen_matrix_symmetric(
343    a: &Matrix,
344    which: &EigenWhich,
345    algorithm: EigenAlgorithm,
346    options: &ArpackOptions,
347) -> Result<SymmetricEigen> {
348    if a.nrow() != a.ncol() {
349        return Err(Error::invalid("eigenvalues need a square matrix"));
350    }
351    check_finite(a.as_slice(), "the matrix")?;
352    symmetric_common(Operator::Dense(a), a.nrow(), which, algorithm, options)
353}
354
355/// Selected eigenvalues and eigenvectors of a real **symmetric** operator
356/// given by a matrix-vector product closure (`igraph_eigen_matrix_symmetric`
357/// with a callback), see [`arpack_rssolve`](super::arpack_rssolve) for the
358/// closure contract. With LAPACK the matrix is first formed by applying the
359/// closure to the unit vectors.
360///
361/// Binds `igraph_eigen_matrix_symmetric` (see the
362/// [linear algebra chapter](https://igraph.org/c/html/latest/igraph-Linalg.html)).
363///
364/// # Errors
365/// As [`eigen_matrix_symmetric`], plus non-finite values produced by the
366/// closure.
367///
368/// # Examples
369///
370/// The operator `x -> (sum x) 1` (the all-ones matrix `J_4`) has eigenvalue
371/// `4` once and `0` three times.
372///
373/// ```
374/// use igraph::linalg::*;
375/// let ones = |x: &[f64], y: &mut [f64]| {
376///     let s: f64 = x.iter().sum();
377///     y.iter_mut().for_each(|yi| *yi = s);
378/// };
379/// let e = eigen_symmetric_fn(4, ones, &EigenWhich::All, EigenAlgorithm::Lapack, &ArpackOptions::default())
380///     .unwrap();
381/// let mut values = e.values.clone();
382/// values.sort_by(f64::total_cmp);
383/// assert!(values[..3].iter().all(|x| x.abs() < 1e-12) && (values[3] - 4.0).abs() < 1e-12);
384/// ```
385pub fn eigen_symmetric_fn<F: FnMut(&[f64], &mut [f64])>(
386    n: usize,
387    matvec: F,
388    which: &EigenWhich,
389    algorithm: EigenAlgorithm,
390    options: &ArpackOptions,
391) -> Result<SymmetricEigen> {
392    with_finite_guard(matvec, |mut f| {
393        let op = Operator::Fn(
394            Some(matvec_trampoline::<&mut dyn FnMut(&[f64], &mut [f64])>),
395            &mut f as *mut &mut dyn FnMut(&[f64], &mut [f64]) as *mut c_void,
396        );
397        symmetric_common(op, n, which, algorithm, options)
398    })
399}
400
401fn complex_common(
402    op: Operator<'_>,
403    n: usize,
404    which: &EigenWhich,
405    algorithm: EigenAlgorithm,
406) -> Result<ComplexEigen> {
407    if n == 0 {
408        return Err(Error::invalid(
409            "cannot compute the eigenvalues of an empty matrix",
410        ));
411    }
412    let cn = to_c_int(n, "matrix order")?;
413    let raw_which = which.to_raw(n)?;
414    let (a, sa) = match op {
415        Operator::Dense(m) => (m as *const Matrix, ptr::null()),
416        Operator::Sparse(s) => (ptr::null(), s as *const SparseMat),
417        Operator::Fn(..) => unreachable!("closures are not supported by igraph_eigen_matrix"),
418    };
419    let mut raw_opts = ArpackOptions::default().to_raw_unchecked_which(n)?;
420    let mut values = VectorComplex::new();
421    let mut vectors = MatrixComplex::new();
422    igraph_call!(igraph_eigen_matrix(
423        a,
424        sa,
425        None,
426        cn,
427        ptr::null_mut(),
428        algorithm.into(),
429        &raw_which,
430        &mut raw_opts,
431        ptr::null_mut(),
432        &mut values,
433        &mut vectors
434    ))?;
435    Ok(ComplexEigen {
436        values: values.to_vec(),
437        vectors: vectors.columns().map(<[Complex]>::to_vec).collect(),
438    })
439}
440
441/// Selected eigenvalues and eigenvectors of a general real dense matrix
442/// (`igraph_eigen_matrix`); only [`EigenAlgorithm::Lapack`] is implemented
443/// by igraph 1.0.0 and 1.0.1.
444///
445/// The eigenvalues are ordered according to `which`; [`EigenWhich::All`]
446/// and [`EigenWhich::Select`] use increasing magnitude. Among eigenvalues of
447/// equal magnitude (resp. real or imaginary part), real ones come first for
448/// the "largest" choices and complex ones come first for the "smallest"
449/// choices, `All` and `Select`, so that each "smallest" order is the exact
450/// reverse of the corresponding "largest" one.
451///
452/// Binds `igraph_eigen_matrix` (see the
453/// [linear algebra chapter](https://igraph.org/c/html/latest/igraph-Linalg.html)).
454/// The callback form of the C function is not exposed: igraph 1.0.0 and 1.0.1
455/// dereferences a null matrix in that code path.
456///
457/// # Examples
458///
459/// A rotation by 90 degrees has eigenvalues `±i`.
460///
461/// ```
462/// use igraph::{linalg::*, prelude::*};
463/// let r = Matrix::from_rows(&[[0.0, -1.0], [1.0, 0.0]]).unwrap();
464/// let e = eigen_matrix(&r, &EigenWhich::All, EigenAlgorithm::Lapack).unwrap();
465/// assert_eq!(e.values.len(), 2);
466/// for v in &e.values {
467///     assert!(v.re().abs() < 1e-12 && (v.im().abs() - 1.0).abs() < 1e-12);
468/// }
469/// ```
470pub fn eigen_matrix(
471    a: &Matrix,
472    which: &EigenWhich,
473    algorithm: EigenAlgorithm,
474) -> Result<ComplexEigen> {
475    if a.nrow() != a.ncol() {
476        return Err(Error::invalid("eigenvalues need a square matrix"));
477    }
478    check_finite(a.as_slice(), "the matrix")?;
479    complex_common(Operator::Dense(a), a.nrow(), which, algorithm)
480}
481
482impl igraph_sparsemat_t {
483    /// Selected eigenvalues and eigenvectors of a symmetric sparse matrix
484    /// (`igraph_eigen_matrix_symmetric` with a sparse input); see
485    /// [`eigen_matrix_symmetric`].
486    pub fn eigen_symmetric(
487        &self,
488        which: &EigenWhich,
489        algorithm: EigenAlgorithm,
490        options: &ArpackOptions,
491    ) -> Result<SymmetricEigen> {
492        if self.nrow() != self.ncol() {
493            return Err(Error::invalid("eigenvalues need a square matrix"));
494        }
495        let cc = self.compress()?;
496        check_finite(&cc.getelements()?.x, "the matrix")?;
497        symmetric_common(
498            Operator::Sparse(&cc),
499            self.nrow(),
500            which,
501            algorithm,
502            options,
503        )
504    }
505
506    /// Selected eigenvalues and eigenvectors of a general sparse matrix
507    /// (`igraph_eigen_matrix` with a sparse input); see [`eigen_matrix`].
508    pub fn eigen(&self, which: &EigenWhich, algorithm: EigenAlgorithm) -> Result<ComplexEigen> {
509        if self.nrow() != self.ncol() {
510            return Err(Error::invalid("eigenvalues need a square matrix"));
511        }
512        let cc = self.compress()?;
513        check_finite(&cc.getelements()?.x, "the matrix")?;
514        complex_common(Operator::Sparse(&cc), self.nrow(), which, algorithm)
515    }
516}
517
518impl igraph_t {
519    /// Eigenvalues and eigenvectors of the adjacency matrix of an undirected
520    /// graph (`igraph_eigen_adjacency`), computed with ARPACK (the only
521    /// algorithm implemented by igraph 1.0.0 and 1.0.1, also chosen by
522    /// [`EigenAlgorithm::Auto`]). Self-loops count once on the diagonal,
523    /// multi-edges add up.
524    ///
525    /// Supported choices: largest/smallest magnitude, largest/smallest
526    /// algebraic.
527    ///
528    /// Binds `igraph_eigen_adjacency` (see the
529    /// [linear algebra chapter](https://igraph.org/c/html/latest/igraph-Linalg.html)).
530    ///
531    /// ARPACK starts from a random vector drawn from the calling thread's
532    /// igraph RNG ([`rng::seed`](crate::rng::seed) makes runs reproducible).
533    /// On small graphs whose wanted eigenvalue is highly repeated (e.g. the
534    /// eigenvalue `-1` of a complete graph) it can stop with "maximum number
535    /// of iterations reached" for some start vectors; LAPACK on the dense
536    /// matrix is the robust choice there.
537    ///
538    /// See also [`Graph::eigenvector_centrality`](crate::Graph::eigenvector_centrality)
539    /// (the scaled leading eigenvector, also for directed and weighted
540    /// graphs) and [`Graph::get_adjacency`](crate::Graph::get_adjacency) to
541    /// form the dense matrix for [`eigen_matrix_symmetric`] or
542    /// [`lapack_dsyevr`](super::lapack_dsyevr).
543    ///
544    /// # Errors
545    /// For directed graphs and unsupported choices
546    /// ([`ErrorKind::Unimplemented`](crate::ErrorKind::Unimplemented)).
547    ///
548    /// # Examples
549    ///
550    /// The complete graph `K_n` has adjacency eigenvalues `n - 1` (once) and
551    /// `-1`.
552    ///
553    /// ```
554    /// use igraph::{linalg::*, prelude::*};
555    /// let k6 = Graph::full(6, false, false).unwrap();
556    /// let opts = ArpackOptions::default();
557    /// let e = k6.eigen_adjacency(&EigenWhich::LargestAlgebraic(1), EigenAlgorithm::Auto, &opts).unwrap();
558    /// assert!((e.values[0] - 5.0).abs() < 1e-10);
559    /// ```
560    pub fn eigen_adjacency(
561        &self,
562        which: &EigenWhich,
563        algorithm: EigenAlgorithm,
564        options: &ArpackOptions,
565    ) -> Result<SymmetricEigen> {
566        let n = self.vcount();
567        if n == 0 {
568            return Err(Error::invalid("the graph has no vertices"));
569        }
570        let raw_which = which.to_raw(n)?;
571        let mut raw_opts = options.to_raw_unchecked_which(n)?;
572        let mut values = Vector::new();
573        let mut vectors = Matrix::new();
574        let mut cvalues = VectorComplex::new();
575        let mut cvectors = MatrixComplex::new();
576        let _arpack = super::arpack::ArpackGuard::enter()?;
577        igraph_call!(igraph_eigen_adjacency(
578            self,
579            algorithm.into(),
580            &raw_which,
581            &mut raw_opts,
582            ptr::null_mut(),
583            &mut values,
584            &mut vectors,
585            &mut cvalues,
586            &mut cvectors
587        ))?;
588        SymmetricEigen {
589            values: values.into(),
590            vectors,
591        }
592        .ordered(which)
593    }
594}
igraph/linalg/eigen.rs

igraph/linalg/
eigen.rs
