Files
adler32
aho_corasick
alga
approx
ascii
atty
backtrace
backtrace_sys
base64
bitflags
blas_src
block_buffer
block_padding
brotli2
brotli_sys
buf_redux
byte_tools
byteorder
cauchy
cblas_sys
cfg_if
chrono
chunked_transfer
colored
crc32fast
crossbeam
crossbeam_channel
crossbeam_deque
crossbeam_epoch
crossbeam_queue
crossbeam_utils
ctrlc
deflate
digest
dirs
error_chain
filetime
futures
generic_array
getrandom
gzip_header
hex
httparse
hyper
idna
itoa
language_tags
lapack_src
lapacke
lapacke_sys
lazy_static
libc
libm
linked_hash_map
log
matches
matrixmultiply
maybe_uninit
md5
memchr
memoffset
mime
mime_guess
multipart
nalgebra
base
geometry
linalg
ndarray
ndarray_linalg
net2
netlib_src
nix
num_complex
num_cpus
num_integer
num_rational
num_traits
opaque_debug
percent_encoding
phf
phf_shared
ppv_lite86
proc_macro2
quick_error
quote
rand
rand_chacha
rand_core
rand_distr
rawpointer
regex
regex_syntax
remove_dir_all
rosrust
rosrust_codegen
rosrust_msg
rouille
rustc_demangle
rustros_tf
ryu
safemem
scopeguard
serde
serde_bytes
serde_derive
serde_json
serde_xml_rs
sha1
siphasher
smallvec
syn
tempdir
term
thread_local
threadpool
time
tiny_http
traitobject
twoway
typeable
typenum
ucd_util
unicase
unicode_bidi
unicode_normalization
unicode_xid
url
utf8_ranges
void
xml
xml_rpc
yaml_rust
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
//! Arnoldi iteration

use super::*;
use crate::{norm::Norm, operator::LinearOperator};
use num_traits::One;
use std::iter::*;

/// Execute Arnoldi iteration as Rust iterator
///
/// - [Arnoldi iteration - Wikipedia](https://en.wikipedia.org/wiki/Arnoldi_iteration)
///
pub struct Arnoldi<A, S, F, Ortho>
where
    A: Scalar,
    S: DataMut<Elem = A>,
    F: LinearOperator<Elem = A>,
    Ortho: Orthogonalizer<Elem = A>,
{
    a: F,
    /// Next vector (normalized `|v|=1`)
    v: ArrayBase<S, Ix1>,
    /// Orthogonalizer
    ortho: Ortho,
    /// Coefficients to be composed into H-matrix
    h: Vec<Array1<A>>,
}

impl<A, S, F, Ortho> Arnoldi<A, S, F, Ortho>
where
    A: Scalar + Lapack,
    S: DataMut<Elem = A>,
    F: LinearOperator<Elem = A>,
    Ortho: Orthogonalizer<Elem = A>,
{
    /// Create an Arnoldi iterator from any linear operator `a`
    pub fn new(a: F, mut v: ArrayBase<S, Ix1>, mut ortho: Ortho) -> Self {
        assert_eq!(ortho.len(), 0);
        assert!(ortho.tolerance() < One::one());
        // normalize before append because |v| may be smaller than ortho.tolerance()
        let norm = v.norm_l2();
        azip!((v in &mut v)  *v = v.div_real(norm));
        ortho.append(v.view());
        Arnoldi {
            a,
            v,
            ortho,
            h: Vec::new(),
        }
    }

    /// Dimension of Krylov subspace
    pub fn dim(&self) -> usize {
        self.ortho.len()
    }

    /// Iterate until convergent
    pub fn complete(mut self) -> (Q<A>, H<A>) {
        for _ in &mut self {} // execute iteration until convergent
        let q = self.ortho.get_q();
        let n = self.h.len();
        let mut h = Array2::zeros((n, n).f());
        for (i, hc) in self.h.iter().enumerate() {
            let m = std::cmp::min(n, i + 2);
            for j in 0..m {
                h[(j, i)] = hc[j];
            }
        }
        (q, h)
    }
}

impl<A, S, F, Ortho> Iterator for Arnoldi<A, S, F, Ortho>
where
    A: Scalar + Lapack,
    S: DataMut<Elem = A>,
    F: LinearOperator<Elem = A>,
    Ortho: Orthogonalizer<Elem = A>,
{
    type Item = Array1<A>;

    fn next(&mut self) -> Option<Self::Item> {
        self.a.apply_mut(&mut self.v);
        let result = self.ortho.div_append(&mut self.v);
        let norm = self.v.norm_l2();
        azip!((v in &mut self.v) *v = v.div_real(norm));
        match result {
            AppendResult::Added(coef) => {
                self.h.push(coef.clone());
                Some(coef)
            }
            AppendResult::Dependent(coef) => {
                self.h.push(coef);
                None
            }
        }
    }
}

/// Utility to execute Arnoldi iteration with Householder reflection
pub fn arnoldi_householder<A, S>(a: impl LinearOperator<Elem = A>, v: ArrayBase<S, Ix1>, tol: A::Real) -> (Q<A>, H<A>)
where
    A: Scalar + Lapack,
    S: DataMut<Elem = A>,
{
    let householder = Householder::new(v.len(), tol);
    Arnoldi::new(a, v, householder).complete()
}

/// Utility to execute Arnoldi iteration with modified Gram-Schmit orthogonalizer
pub fn arnoldi_mgs<A, S>(a: impl LinearOperator<Elem = A>, v: ArrayBase<S, Ix1>, tol: A::Real) -> (Q<A>, H<A>)
where
    A: Scalar + Lapack,
    S: DataMut<Elem = A>,
{
    let mgs = MGS::new(v.len(), tol);
    Arnoldi::new(a, v, mgs).complete()
}