diff --git a/src/krylov/arnoldi.rs b/src/krylov/arnoldi.rs
new file mode 100644
index 00000000..20cb0098
--- /dev/null
+++ b/src/krylov/arnoldi.rs
@@ -0,0 +1,135 @@
+//! Arnoldi iteration
+
+use super::*;
+use crate::norm::Norm;
+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: Fn(&mut ArrayBase<S, Ix1>),
+    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: Fn(&mut ArrayBase<S, Ix1>),
+    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!(mut v(&mut v) in { *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: Fn(&mut ArrayBase<S, Ix1>),
+    Ortho: Orthogonalizer<Elem = A>,
+{
+    type Item = Array1<A>;
+
+    fn next(&mut self) -> Option<Self::Item> {
+        (self.a)(&mut self.v);
+        let result = self.ortho.div_append(&mut self.v);
+        let norm = self.v.norm_l2();
+        azip!(mut v(&mut self.v) in { *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
+            }
+        }
+    }
+}
+
+/// Interpret a matrix as a linear operator
+pub fn mul_mat<A, S1, S2>(a: ArrayBase<S1, Ix2>) -> impl Fn(&mut ArrayBase<S2, Ix1>)
+where
+    A: Scalar,
+    S1: Data<Elem = A>,
+    S2: DataMut<Elem = A>,
+{
+    let (n, m) = a.dim();
+    assert_eq!(n, m, "Input matrix must be square");
+    move |x| {
+        assert_eq!(m, x.len(), "Input matrix and vector sizes mismatch");
+        let ax = a.dot(x);
+        azip!(mut x(x), ax in { *x = ax });
+    }
+}
+
+/// Utility to execute Arnoldi iteration with Householder reflection
+pub fn arnoldi_householder<A, S1, S2>(a: ArrayBase<S1, Ix2>, v: ArrayBase<S2, Ix1>, tol: A::Real) -> (Q<A>, H<A>)
+where
+    A: Scalar + Lapack,
+    S1: Data<Elem = A>,
+    S2: DataMut<Elem = A>,
+{
+    let householder = Householder::new(v.len(), tol);
+    Arnoldi::new(mul_mat(a), v, householder).complete()
+}
+
+/// Utility to execute Arnoldi iteration with modified Gram-Schmit orthogonalizer
+pub fn arnoldi_mgs<A, S1, S2>(a: ArrayBase<S1, Ix2>, v: ArrayBase<S2, Ix1>, tol: A::Real) -> (Q<A>, H<A>)
+where
+    A: Scalar + Lapack,
+    S1: Data<Elem = A>,
+    S2: DataMut<Elem = A>,
+{
+    let mgs = MGS::new(v.len(), tol);
+    Arnoldi::new(mul_mat(a), v, mgs).complete()
+}
diff --git a/src/krylov/householder.rs b/src/krylov/householder.rs
index 391b18c4..17b7d905 100644
--- a/src/krylov/householder.rs
+++ b/src/krylov/householder.rs
@@ -8,17 +8,17 @@ use crate::{inner::*, norm::*};
 use num_traits::One;
 
 /// Calc a reflactor `w` from a vector `x`
-pub fn calc_reflector<A, S>(x: &mut ArrayBase<S, Ix1>) -> A
+pub fn calc_reflector<A, S>(x: &mut ArrayBase<S, Ix1>)
 where
     A: Scalar + Lapack,
     S: DataMut<Elem = A>,
 {
+    assert!(x.len() > 0);
     let norm = x.norm_l2();
     let alpha = -x[0].mul_real(norm / x[0].abs());
     x[0] -= alpha;
     let inv_rev_norm = A::Real::one() / x.norm_l2();
     azip!(mut a(x) in { *a = a.mul_real(inv_rev_norm)});
-    alpha
 }
 
 /// Take a reflection `P = I - 2ww^T`
@@ -50,12 +50,19 @@ pub struct Householder<A: Scalar> {
     ///
     /// The coefficient is copied into another array, and this does not contain
     v: Vec<Array1<A>>,
+
+    /// Tolerance
+    tol: A::Real,
 }
 
 impl<A: Scalar + Lapack> Householder<A> {
     /// Create a new orthogonalizer
-    pub fn new(dim: usize) -> Self {
-        Householder { dim, v: Vec::new() }
+    pub fn new(dim: usize, tol: A::Real) -> Self {
+        Householder {
+            dim,
+            v: Vec::new(),
+            tol,
+        }
     }
 
     /// Take a Reflection `P = I - 2ww^T`
@@ -92,12 +99,32 @@ impl<A: Scalar + Lapack> Householder<A> {
         }
     }
 
-    fn eval_residual<S>(&self, a: &ArrayBase<S, Ix1>) -> A::Real
+    /// Compose coefficients array using reflected vector
+    fn compose_coefficients<S>(&self, a: &ArrayBase<S, Ix1>) -> Coefficients<A>
     where
         S: Data<Elem = A>,
     {
-        let l = self.v.len();
-        a.slice(s![l..]).norm_l2()
+        let k = self.len();
+        let res = a.slice(s![k..]).norm_l2();
+        let mut c = Array1::zeros(k + 1);
+        azip!(mut c(c.slice_mut(s![..k])), a(a.slice(s![..k])) in { *c = a });
+        if k < a.len() {
+            let ak = a[k];
+            c[k] = -ak.mul_real(res / ak.abs());
+        } else {
+            c[k] = A::from_real(res);
+        }
+        c
+    }
+
+    /// Construct the residual vector from reflected vector
+    fn construct_residual<S>(&self, a: &mut ArrayBase<S, Ix1>)
+    where
+        S: DataMut<Elem = A>,
+    {
+        let k = self.len();
+        azip!(mut a( a.slice_mut(s![..k])) in { *a = A::zero() });
+        self.backward_reflection(a);
     }
 }
 
@@ -112,45 +139,61 @@ impl<A: Scalar + Lapack> Orthogonalizer for Householder<A> {
         self.v.len()
     }
 
+    fn tolerance(&self) -> A::Real {
+        self.tol
+    }
+
+    fn decompose<S>(&self, a: &mut ArrayBase<S, Ix1>) -> Array1<A>
+    where
+        S: DataMut<Elem = A>,
+    {
+        self.forward_reflection(a);
+        let coef = self.compose_coefficients(a);
+        self.construct_residual(a);
+        coef
+    }
+
     fn coeff<S>(&self, a: ArrayBase<S, Ix1>) -> Array1<A>
     where
         S: Data<Elem = A>,
     {
         let mut a = a.into_owned();
         self.forward_reflection(&mut a);
-        let res = self.eval_residual(&a);
-        let k = self.len();
-        let mut c = Array1::zeros(k + 1);
-        azip!(mut c(c.slice_mut(s![..k])), a(a.slice(s![..k])) in { *c = a });
-        c[k] = A::from_real(res);
-        c
+        self.compose_coefficients(&a)
     }
 
-    fn append<S>(&mut self, mut a: ArrayBase<S, Ix1>, rtol: A::Real) -> Result<Array1<A>, Array1<A>>
+    fn div_append<S>(&mut self, a: &mut ArrayBase<S, Ix1>) -> AppendResult<A>
     where
         S: DataMut<Elem = A>,
     {
         assert_eq!(a.len(), self.dim);
         let k = self.len();
-
-        self.forward_reflection(&mut a);
-        let mut coef = Array::zeros(k + 1);
-        for i in 0..k {
-            coef[i] = a[i];
-        }
-        if self.is_full() {
-            return Err(coef); // coef[k] must be zero in this case
+        self.forward_reflection(a);
+        let coef = self.compose_coefficients(a);
+        if coef[k].abs() < self.tol {
+            return AppendResult::Dependent(coef);
         }
+        calc_reflector(&mut a.slice_mut(s![k..]));
+        self.v.push(a.to_owned());
+        self.construct_residual(a);
+        AppendResult::Added(coef)
+    }
 
-        let alpha = calc_reflector(&mut a.slice_mut(s![k..]));
-        coef[k] = alpha;
-
-        if alpha.abs() < rtol {
-            // linearly dependent
-            return Err(coef);
+    fn append<S>(&mut self, a: ArrayBase<S, Ix1>) -> AppendResult<A>
+    where
+        S: Data<Elem = A>,
+    {
+        assert_eq!(a.len(), self.dim);
+        let mut a = a.into_owned();
+        let k = self.len();
+        self.forward_reflection(&mut a);
+        let coef = self.compose_coefficients(&a);
+        if coef[k].abs() < self.tol {
+            return AppendResult::Dependent(coef);
         }
-        self.v.push(a.into_owned());
-        Ok(coef)
+        calc_reflector(&mut a.slice_mut(s![k..]));
+        self.v.push(a.to_owned());
+        AppendResult::Added(coef)
     }
 
     fn get_q(&self) -> Q<A> {
@@ -175,8 +218,8 @@ where
     A: Scalar + Lapack,
     S: Data<Elem = A>,
 {
-    let h = Householder::new(dim);
-    qr(iter, h, rtol, strategy)
+    let h = Householder::new(dim, rtol);
+    qr(iter, h, strategy)
 }
 
 #[cfg(test)]
diff --git a/src/krylov/mgs.rs b/src/krylov/mgs.rs
index 9caff111..b6a8d2bf 100644
--- a/src/krylov/mgs.rs
+++ b/src/krylov/mgs.rs
@@ -5,28 +5,44 @@ use crate::{generate::*, inner::*, norm::Norm};
 
 /// Iterative orthogonalizer using modified Gram-Schmit procedure
 #[derive(Debug, Clone)]
-pub struct MGS<A> {
+pub struct MGS<A: Scalar> {
     /// Dimension of base space
-    dimension: usize,
+    dim: usize,
+
     /// Basis of spanned space
     q: Vec<Array1<A>>,
+
+    /// Tolerance
+    tol: A::Real,
 }
 
 impl<A: Scalar + Lapack> MGS<A> {
     /// Create an empty orthogonalizer
-    pub fn new(dimension: usize) -> Self {
+    pub fn new(dim: usize, tol: A::Real) -> Self {
         Self {
-            dimension,
+            dim,
             q: Vec::new(),
+            tol,
         }
     }
+}
 
-    /// Orthogonalize given vector against to the current basis
-    ///
-    /// - Returned array is coefficients and residual norm
-    /// - `a` will contain the residual vector
-    ///
-    pub fn orthogonalize<S>(&self, a: &mut ArrayBase<S, Ix1>) -> Array1<A>
+impl<A: Scalar + Lapack> Orthogonalizer for MGS<A> {
+    type Elem = A;
+
+    fn dim(&self) -> usize {
+        self.dim
+    }
+
+    fn len(&self) -> usize {
+        self.q.len()
+    }
+
+    fn tolerance(&self) -> A::Real {
+        self.tol
+    }
+
+    fn decompose<S>(&self, a: &mut ArrayBase<S, Ix1>) -> Array1<A>
     where
         S: DataMut<Elem = A>,
     {
@@ -42,18 +58,6 @@ impl<A: Scalar + Lapack> MGS<A> {
         coef[self.len()] = A::from_real(nrm);
         coef
     }
-}
-
-impl<A: Scalar + Lapack> Orthogonalizer for MGS<A> {
-    type Elem = A;
-
-    fn dim(&self) -> usize {
-        self.dimension
-    }
-
-    fn len(&self) -> usize {
-        self.q.len()
-    }
 
     fn coeff<S>(&self, a: ArrayBase<S, Ix1>) -> Array1<A>
     where
@@ -61,24 +65,32 @@ impl<A: Scalar + Lapack> Orthogonalizer for MGS<A> {
         S: Data<Elem = A>,
     {
         let mut a = a.into_owned();
-        self.orthogonalize(&mut a)
+        self.decompose(&mut a)
     }
 
-    fn append<S>(&mut self, a: ArrayBase<S, Ix1>, rtol: A::Real) -> Result<Array1<A>, Array1<A>>
+    fn append<S>(&mut self, a: ArrayBase<S, Ix1>) -> AppendResult<A>
     where
         A: Lapack,
         S: Data<Elem = A>,
     {
         let mut a = a.into_owned();
-        let coef = self.orthogonalize(&mut a);
+        self.div_append(&mut a)
+    }
+
+    fn div_append<S>(&mut self, mut a: &mut ArrayBase<S, Ix1>) -> AppendResult<A>
+    where
+        A: Lapack,
+        S: DataMut<Elem = A>,
+    {
+        let coef = self.decompose(&mut a);
         let nrm = coef[coef.len() - 1].re();
-        if nrm < rtol {
+        if nrm < self.tol {
             // Linearly dependent
-            return Err(coef);
+            return AppendResult::Dependent(coef);
         }
-        azip!(mut a in { *a = *a / A::from_real(nrm) });
-        self.q.push(a);
-        Ok(coef)
+        azip!(mut a(&mut *a) in { *a = *a / A::from_real(nrm) });
+        self.q.push(a.to_owned());
+        AppendResult::Added(coef)
     }
 
     fn get_q(&self) -> Q<A> {
@@ -97,6 +109,6 @@ where
     A: Scalar + Lapack,
     S: Data<Elem = A>,
 {
-    let mgs = MGS::new(dim);
-    qr(iter, mgs, rtol, strategy)
+    let mgs = MGS::new(dim, rtol);
+    qr(iter, mgs, strategy)
 }
diff --git a/src/krylov/mod.rs b/src/krylov/mod.rs
index b53df0c7..14c9f7ef 100644
--- a/src/krylov/mod.rs
+++ b/src/krylov/mod.rs
@@ -1,11 +1,13 @@
-//! Krylov subspace
+//! Krylov subspace methods
 
 use crate::types::*;
 use ndarray::*;
 
+pub mod arnoldi;
 pub mod householder;
 pub mod mgs;
 
+pub use arnoldi::{arnoldi_householder, arnoldi_mgs, Arnoldi};
 pub use householder::{householder, Householder};
 pub use mgs::{mgs, MGS};
 
@@ -23,26 +25,44 @@ pub type Q<A> = Array2<A>;
 ///
 pub type R<A> = Array2<A>;
 
+/// H-matrix
+///
+/// - Maybe **NOT** square
+/// - Hessenberg matrix
+///
+pub type H<A> = Array2<A>;
+
+/// Array type for coefficients to the current basis
+///
+/// - The length must be `self.len() + 1`
+/// - Last component is the residual norm
+///
+pub type Coefficients<A> = Array1<A>;
+
 /// Trait for creating orthogonal basis from iterator of arrays
 ///
+/// Panic
+/// -------
+/// - if the size of the input array mismatches to the dimension
+///
 /// Example
 /// -------
 ///
 /// ```rust
 /// # use ndarray::*;
 /// # use ndarray_linalg::{krylov::*, *};
-/// let mut mgs = MGS::new(3);
-/// let coef = mgs.append(array![0.0, 1.0, 0.0], 1e-9).unwrap();
+/// let mut mgs = MGS::new(3, 1e-9);
+/// let coef = mgs.append(array![0.0, 1.0, 0.0]).into_coeff();
 /// close_l2(&coef, &array![1.0], 1e-9);
 ///
-/// let coef = mgs.append(array![1.0, 1.0, 0.0], 1e-9).unwrap();
+/// let coef = mgs.append(array![1.0, 1.0, 0.0]).into_coeff();
 /// close_l2(&coef, &array![1.0, 1.0], 1e-9);
 ///
 /// // Fail if the vector is linearly dependent
-/// assert!(mgs.append(array![1.0, 2.0, 0.0], 1e-9).is_err());
+/// assert!(mgs.append(array![1.0, 2.0, 0.0]).is_dependent());
 ///
 /// // You can get coefficients of dependent vector
-/// if let Err(coef) = mgs.append(array![1.0, 2.0, 0.0], 1e-9) {
+/// if let AppendResult::Dependent(coef) = mgs.append(array![1.0, 2.0, 0.0]) {
 ///     close_l2(&coef, &array![2.0, 1.0, 0.0], 1e-9);
 /// }
 /// ```
@@ -64,37 +84,35 @@ pub trait Orthogonalizer {
         self.len() == 0
     }
 
-    /// Calculate the coefficient to the given basis and residual norm
+    fn tolerance(&self) -> <Self::Elem as Scalar>::Real;
+
+    /// Decompose given vector into the span of current basis and
+    /// its tangent space
+    ///
+    /// - `a` becomes the tangent vector
+    /// - The Coefficients to the current basis is returned.
     ///
-    /// - The length of the returned array must be `self.len() + 1`
-    /// - Last component is the residual norm
+    fn decompose<S>(&self, a: &mut ArrayBase<S, Ix1>) -> Coefficients<Self::Elem>
+    where
+        S: DataMut<Elem = Self::Elem>;
+
+    /// Calculate the coefficient to the current basis basis
     ///
-    /// Panic
-    /// -------
-    /// - if the size of the input array mismatches to the dimension
+    /// - This will be faster than `decompose` because the construction of the residual vector may
+    ///   requires more Calculation
     ///
-    fn coeff<S>(&self, a: ArrayBase<S, Ix1>) -> Array1<Self::Elem>
+    fn coeff<S>(&self, a: ArrayBase<S, Ix1>) -> Coefficients<Self::Elem>
     where
         S: Data<Elem = Self::Elem>;
 
     /// Add new vector if the residual is larger than relative tolerance
-    ///
-    /// Returns
-    /// --------
-    /// Coefficients to the `i`-th Q-vector
-    ///
-    /// - The size of array must be `self.len() + 1`
-    /// - The last element is the residual norm of input vector
-    ///
-    /// Panic
-    /// -------
-    /// - if the size of the input array mismatches to the dimension
-    ///
-    fn append<S>(
-        &mut self,
-        a: ArrayBase<S, Ix1>,
-        rtol: <Self::Elem as Scalar>::Real,
-    ) -> Result<Array1<Self::Elem>, Array1<Self::Elem>>
+    fn append<S>(&mut self, a: ArrayBase<S, Ix1>) -> AppendResult<Self::Elem>
+    where
+        S: Data<Elem = Self::Elem>;
+
+    /// Add new vector if the residual is larger than relative tolerance,
+    /// and return the residual vector
+    fn div_append<S>(&mut self, a: &mut ArrayBase<S, Ix1>) -> AppendResult<Self::Elem>
     where
         S: DataMut<Elem = Self::Elem>;
 
@@ -102,6 +120,39 @@ pub trait Orthogonalizer {
     fn get_q(&self) -> Q<Self::Elem>;
 }
 
+pub enum AppendResult<A> {
+    Added(Coefficients<A>),
+    Dependent(Coefficients<A>),
+}
+
+impl<A: Scalar> AppendResult<A> {
+    pub fn into_coeff(self) -> Coefficients<A> {
+        match self {
+            AppendResult::Added(c) => c,
+            AppendResult::Dependent(c) => c,
+        }
+    }
+
+    pub fn is_dependent(&self) -> bool {
+        match self {
+            AppendResult::Added(_) => false,
+            AppendResult::Dependent(_) => true,
+        }
+    }
+
+    pub fn coeff(&self) -> &Coefficients<A> {
+        match self {
+            AppendResult::Added(c) => c,
+            AppendResult::Dependent(c) => c,
+        }
+    }
+
+    pub fn residual_norm(&self) -> A::Real {
+        let c = self.coeff();
+        c[c.len() - 1].abs()
+    }
+}
+
 /// Strategy for linearly dependent vectors appearing in iterative QR decomposition
 #[derive(Clone, Copy, Debug, Eq, PartialEq)]
 pub enum Strategy {
@@ -127,7 +178,6 @@ pub enum Strategy {
 pub fn qr<A, S>(
     iter: impl Iterator<Item = ArrayBase<S, Ix1>>,
     mut ortho: impl Orthogonalizer<Elem = A>,
-    rtol: A::Real,
     strategy: Strategy,
 ) -> (Q<A>, R<A>)
 where
@@ -138,9 +188,9 @@ where
 
     let mut coefs = Vec::new();
     for a in iter {
-        match ortho.append(a.into_owned(), rtol) {
-            Ok(coef) => coefs.push(coef),
-            Err(coef) => match strategy {
+        match ortho.append(a.into_owned()) {
+            AppendResult::Added(coef) => coefs.push(coef),
+            AppendResult::Dependent(coef) => match strategy {
                 Strategy::Terminate => break,
                 Strategy::Skip => continue,
                 Strategy::Full => coefs.push(coef),
diff --git a/tests/arnoldi.rs b/tests/arnoldi.rs
new file mode 100644
index 00000000..bbc15553
--- /dev/null
+++ b/tests/arnoldi.rs
@@ -0,0 +1,62 @@
+use ndarray::*;
+use ndarray_linalg::{krylov::*, *};
+
+#[test]
+fn aq_qh_mgs() {
+    let a: Array2<f64> = random((5, 5));
+    let v: Array1<f64> = random(5);
+    let (q, h) = arnoldi_mgs(a.clone(), v, 1e-9);
+    println!("A = \n{:?}", &a);
+    println!("Q = \n{:?}", &q);
+    println!("H = \n{:?}", &h);
+    let aq = a.dot(&q);
+    let qh = q.dot(&h);
+    println!("AQ = \n{:?}", &aq);
+    println!("QH = \n{:?}", &qh);
+    close_l2(&aq, &qh, 1e-9);
+}
+
+#[test]
+fn aq_qh_householder() {
+    let a: Array2<f64> = random((5, 5));
+    let v: Array1<f64> = random(5);
+    let (q, h) = arnoldi_mgs(a.clone(), v, 1e-9);
+    println!("A = \n{:?}", &a);
+    println!("Q = \n{:?}", &q);
+    println!("H = \n{:?}", &h);
+    let aq = a.dot(&q);
+    let qh = q.dot(&h);
+    println!("AQ = \n{:?}", &aq);
+    println!("QH = \n{:?}", &qh);
+    close_l2(&aq, &qh, 1e-9);
+}
+
+#[test]
+fn aq_qh_mgs_complex() {
+    let a: Array2<c64> = random((5, 5));
+    let v: Array1<c64> = random(5);
+    let (q, h) = arnoldi_mgs(a.clone(), v, 1e-9);
+    println!("A = \n{:?}", &a);
+    println!("Q = \n{:?}", &q);
+    println!("H = \n{:?}", &h);
+    let aq = a.dot(&q);
+    let qh = q.dot(&h);
+    println!("AQ = \n{:?}", &aq);
+    println!("QH = \n{:?}", &qh);
+    close_l2(&aq, &qh, 1e-9);
+}
+
+#[test]
+fn aq_qh_householder_complex() {
+    let a: Array2<c64> = random((5, 5));
+    let v: Array1<c64> = random(5);
+    let (q, h) = arnoldi_mgs(a.clone(), v, 1e-9);
+    println!("A = \n{:?}", &a);
+    println!("Q = \n{:?}", &q);
+    println!("H = \n{:?}", &h);
+    let aq = a.dot(&q);
+    let qh = q.dot(&h);
+    println!("AQ = \n{:?}", &aq);
+    println!("QH = \n{:?}", &qh);
+    close_l2(&aq, &qh, 1e-9);
+}