diff --git a/include/irlba/irlba.hpp b/include/irlba/irlba.hpp
index 228c7fa..763aa9d 100644
--- a/include/irlba/irlba.hpp
+++ b/include/irlba/irlba.hpp
@@ -154,7 +154,7 @@ class Irlba {
     /** 
      * Run IRLBA on an input matrix to perform an approximate SVD, with arbitrary centering and scaling operations.
      *
-     * @tparam M Matrix class, typically from the **Eigen** matrix manipulation library.
+     * @tparam Input Matrix class, typically from the **Eigen** matrix manipulation library.
      * However, other classes are also supported, see the other `run()` methods for details.
      * @tparam Engine A (pseudo-)random number generator class, returning a randomly sampled value when called as a functor with no arguments.
      *
@@ -182,19 +182,22 @@ class Irlba {
      * Note that `scale=true` requires `center=true` to guarantee unit variance along each column. 
      * No scaling is performed when the variance of a column is zero, so as to avoid divide-by-zero errors. 
      */
-    template<class M, class Engine = std::mt19937_64>
+    template<class Input, class Engine = std::mt19937_64>
     std::pair<bool, int> run(
-        const M& mat, 
+        const Input& mat, 
         bool center, 
         bool scale, 
-        Eigen::MatrixXd& outU, 
-        Eigen::MatrixXd& outV, 
-        Eigen::VectorXd& outD, 
+        MatrixOf<typename Input::Scalar>& outU, 
+        MatrixOf<typename Input::Scalar>& outV, 
+        RealVectorOf<typename Input::Scalar>& outD, 
         Engine* eng = null_rng(),
-        Eigen::VectorXd* init = NULL) 
+        VectorOf<typename Input::Scalar>* init = NULL) 
     const {
+        typedef VectorOf<typename Input::Scalar> Vector;
+        typedef typename Input::Scalar Scalar;
+
         if (scale || center) {
-            Eigen::VectorXd center0, scale0;
+            Vector center0, scale0;
 
             if (center) {
                 if (mat.rows() < 1) {
@@ -211,14 +214,14 @@ class Irlba {
             }
 
             for (Eigen::Index i = 0; i < mat.cols(); ++i) {
-                double mean = 0;
+                Scalar mean = 0;
                 if (center) {
                     mean = mat.col(i).sum() / mat.rows();
                     center0[i] = mean;
                 }
                 if (scale) {
-                    Eigen::VectorXd current = mat.col(i); // force it to be a VectorXd, even if it's a sparse matrix.
-                    double var = 0;
+                    Vector current = mat.col(i); // force the column to be a Vector, even if it's coming from a sparse matrix.
+                    Scalar var = 0;
                     for (auto x : current) {
                         var += (x - mean)*(x - mean);
                     }
@@ -232,7 +235,7 @@ class Irlba {
             }
 
             if (center) {
-                Centered<M> centered(&mat, &center0);
+                Centered<Input> centered(&mat, &center0);
                 if (scale) {
                     Scaled<decltype(centered)> centered_scaled(&centered, &scale0);
                     return run(centered_scaled, outU, outV, outD, eng, init);
@@ -240,7 +243,7 @@ class Irlba {
                     return run(centered, outU, outV, outD, eng, init);
                 }
             } else {
-                Scaled<M> scaled(&mat, &scale0);
+                Scaled<Input> scaled(&mat, &scale0);
                 return run(scaled, outU, outV, outD, eng, init);
             }
         } else {
@@ -251,8 +254,8 @@ class Irlba {
     /** 
      * Run IRLBA on an input matrix to perform an approximate SVD.
      *
-     * @tparam M Matrix class, most typically from the **Eigen** matrix manipulation library.
-     * However, custom classes are also supported, see below for details.
+     * @tparam Input Matrix class, typically from the **Eigen** matrix manipulation library.
+     * However, other classes are also supported, see below for details.
      * @tparam Engine A (pseudo-)random number generator class, returning a randomly sampled value when called as a functor with no arguments.
      *
      * @param[in] mat Input matrix.
@@ -277,31 +280,31 @@ class Irlba {
      * - A `rows()` method that returns the number of rows.
      * - A `cols()` method that returns the number of columns.
      * - One of the following for matrix-vector multiplication:
-     *   - `multiply(rhs, out)`, which should compute the product of the matrix with `rhs`, a `Eigen::VectorXd`-equivalent of length equal to the number of columns;
-     *     and stores the result in `out`, an `Eigen::VectorXd` of length equal to the number of rows.
-     *   - A `*` method where the right-hand side is an `Eigen::VectorXd` (or equivalent expression) of length equal to the number of columsn,
-     *     and returns an `Eigen::VectorXd`-equivalent of length equal to the number of rows.
+     *   - `multiply(rhs, out)`, which should compute the product of the matrix with `rhs`, a `Vector`-equivalent of length equal to the number of columns;
+     *     and stores the result in `out`, an `Vector` of length equal to the number of rows.
+     *   - A `*` method where the right-hand side is an `Vector` (or equivalent expression) of length equal to the number of columsn,
+     *     and returns an `Vector`-equivalent of length equal to the number of rows.
      * - One of the following for matrix transpose-vector multiplication:
-     *   - `adjoint_multiply(rhs, out)`, which should compute the product of the matrix transpose with `rhs`, a `Eigen::VectorXd`-equivalent of length equal to the number of rows;
-     *     and stores the result in `out`, an `Eigen::VectorXd` of length equal to the number of columns.
+     *   - `adjoint_multiply(rhs, out)`, which should compute the product of the matrix transpose with `rhs`, a `Vector`-equivalent of length equal to the number of rows;
+     *     and stores the result in `out`, an `Vector` of length equal to the number of columns.
      *   - An `adjoint()` method that returns an instance of any class that has a `*` method for matrix-vector multiplication.
-     *     The method should accept an `Eigen::VectorXd`-equivalent of length equal to the number of rows,
-     *     and return an `Eigen::VectorXd`-equvialent of length equal to the number of columns.
-     * - A `realize()` method that returns an `Eigen::MatrixXd` object representing the modified matrix.
-     *   This can be omitted if an `Eigen::MatrixXd` can be copy-constructed from the class.
+     *     The method should accept an `Vector`-equivalent of length equal to the number of rows,
+     *     and return an `Vector`-equvialent of length equal to the number of columns.
+     * - A `realize()` method that returns an `MatrixOf<typename Input::Scalar>` object representing the modified matrix.
+     *   This can be omitted if an `MatrixOf<typename Input::Scalar>` can be copy-constructed from the class.
      *
      * See the `Centered` and `Scaled` classes for more details.
      *
      * If the smallest dimension of `mat` is below 6, this method falls back to performing an exact SVD.
      */
-    template<class Matrix, class Engine = std::mt19937_64>
+    template<class Input, class Engine = std::mt19937_64>
     std::pair<bool, int> run(
-        const Matrix& mat, 
-        Eigen::MatrixXd& outU, 
-        Eigen::MatrixXd& outV, 
-        Eigen::VectorXd& outD, 
+        const Input& mat, 
+        MatrixOf<typename Input::Scalar>& outU, 
+        MatrixOf<typename Input::Scalar>& outV, 
+        RealVectorOf<typename Input::Scalar>& outD, 
         Engine* eng = null_rng(),
-        Eigen::VectorXd* init = NULL) 
+        VectorOf<typename Input::Scalar>* init = NULL) 
     const {
         if (eng == NULL) {
             std::mt19937_64 rng(seed);
@@ -312,15 +315,19 @@ class Irlba {
     }
 
 private:
-    template<class M, class Engine>
+    template<class Input, class Engine>
     std::pair<bool, int> run_internal(
-        const M& mat, 
+        const Input& mat, 
         Engine& eng, 
-        Eigen::MatrixXd& outU, 
-        Eigen::MatrixXd& outV, 
-        Eigen::VectorXd& outD, 
-        Eigen::VectorXd* init)
+        MatrixOf<typename Input::Scalar>& outU, 
+        MatrixOf<typename Input::Scalar>& outV, 
+        RealVectorOf<typename Input::Scalar>& outD, 
+        VectorOf<typename Input::Scalar>* init)
     const {
+        typedef MatrixOf<typename Input::Scalar> Matrix;
+        typedef VectorOf<typename Input::Scalar> Vector;
+        typedef RealVectorOf<typename Input::Scalar> RealVector;
+
         const int smaller = std::min(mat.rows(), mat.cols());
         if (number >= smaller) {
             throw std::runtime_error("requested number of singular values must be less than smaller matrix dimension");
@@ -334,7 +341,7 @@ class Irlba {
 
         const int work = std::min(number + extra_work, smaller);
 
-        Eigen::MatrixXd V(mat.cols(), work);
+        Matrix V(mat.cols(), work);
         if (init) {
             if (init->size() != V.rows()) {
                 throw std::runtime_error("initialization vector does not have expected number of rows");
@@ -347,20 +354,20 @@ class Irlba {
 
         bool converged = false;
         int iter = 0, k =0;
-        Eigen::JacobiSVD<Eigen::MatrixXd> svd(work, work, Eigen::ComputeThinU | Eigen::ComputeThinV);
+        Eigen::JacobiSVD<Matrix> svd(work, work, Eigen::ComputeThinU | Eigen::ComputeThinV);
 
         auto lptmp = lp.initialize(mat);
 
-        Eigen::MatrixXd W(mat.rows(), work);
-        Eigen::MatrixXd Wtmp(mat.rows(), work);
-        Eigen::MatrixXd Vtmp(mat.cols(), work);
+        Matrix W(mat.rows(), work);
+        Matrix Wtmp(mat.rows(), work);
+        Matrix Vtmp(mat.cols(), work);
 
-        Eigen::MatrixXd B(work, work);
+        Matrix B(work, work);
         B.setZero(work, work);
-        Eigen::VectorXd res(work);
-        Eigen::VectorXd F(mat.cols());
+        Vector res(work);
+        Vector F(mat.cols());
 
-        Eigen::VectorXd prevS(work);
+        RealVector prevS(work);
 
         for (; iter < maxit; ++iter) {
             // Technically, this is only a 'true' Lanczos bidiagonalization
@@ -380,7 +387,9 @@ class Irlba {
             const auto& BV = svd.matrixV();
 
             // Checking for convergence.
-            if (B(work-1, work-1) == 0) { // a.k.a. the final value of 'S' from the Lanczos iterations.
+            // Bottom-left element of B is the final value of 'S' from the Lanczos iterations.
+            // Note that this still works for complex numbers, as 0+0i == 0.0.
+            if (B(work-1, work-1) == 0.0) { 
                 converged = true;
                 break;
             }
@@ -455,20 +464,19 @@ class Irlba {
     }
 
 private:
-    template<class M>
-    void exact(const M& mat, Eigen::MatrixXd& outU, Eigen::MatrixXd& outV, Eigen::VectorXd& outD) const {
-        Eigen::BDCSVD<Eigen::MatrixXd> svd(mat.rows(), mat.cols(), Eigen::ComputeThinU | Eigen::ComputeThinV);
-
-        if constexpr(std::is_same<M, Eigen::MatrixXd>::value) {
+    template<class Input>
+    void exact(const Input& mat, MatrixOf<typename Input::Scalar>& outU, MatrixOf<typename Input::Scalar>& outV, RealVectorOf<typename Input::Scalar>& outD) const {
+        typedef MatrixOf<typename Input::Scalar> Matrix;
+        Eigen::BDCSVD<Matrix> svd(mat.rows(), mat.cols(), Eigen::ComputeThinU | Eigen::ComputeThinV);
+
+        if constexpr(has_realize_method<Input>::value) {
+            Matrix adjusted = mat.realize();
+            svd.compute(adjusted);
+        } else if constexpr(std::is_same<Matrix, Input>::value) {
             svd.compute(mat);
         } else {
-            if constexpr(has_realize_method<M>::value) {
-                Eigen::MatrixXd adjusted = mat.realize();
-                svd.compute(adjusted);
-            } else {
-                Eigen::MatrixXd adjusted(mat);
-                svd.compute(adjusted);
-            }
+            Matrix adjusted(mat);
+            svd.compute(adjusted);
         }
 
         outD.resize(number);
@@ -486,27 +494,30 @@ class Irlba {
 public:
     /**
      * Result of the IRLBA-based decomposition.
+     *
+     * @tparam Scalar Scalar type for the left/right singular vectors.
      */
+    template<class Scalar>
     struct Results {
         /**
          * The left singular vectors, stored as columns of `U`.
          * The number of rows in `U` is equal to the number of rows in the input matrix,
          * and the number of columns is equal to the number of requested vectors.
          */
-        Eigen::MatrixXd U;
+        MatrixOf<Scalar> U;
 
         /**
          * The right singular vectors, stored as columns of `U`.
          * The number of rows in `U` is equal to the number of columns in the input matrix,
          * and the number of columns is equal to the number of requested vectors.
          */
-        Eigen::MatrixXd V;
+        MatrixOf<Scalar> V;
 
         /**
          * The requested number of singular values, ordered by decreasing value.
          * These correspond to the vectors in `U` and `V`.
          */
-        Eigen::VectorXd D;
+        RealVectorOf<Scalar> D;
 
         /**
          * Whether the algorithm converged.
@@ -522,7 +533,7 @@ class Irlba {
     /** 
      * Run IRLBA on an input matrix to perform an approximate SVD with centering and scaling.
      * 
-     * @tparam M Matrix class, most typically from the **Eigen** matrix manipulation library.
+     * @tparam Input Matrix class, typically from the **Eigen** matrix manipulation library.
      * However, other classes are also supported, see the other `run()` methods for details.
      * @tparam Engine A (pseudo-)random number generator class, returning a randomly sampled value when called as a functor with no arguments.
      *
@@ -536,9 +547,9 @@ class Irlba {
      *
      * @return A `Results` object containing the singular vectors and values, as well as some statistics on convergence.
      */
-    template<class M, class Engine = std::mt19937_64>
-    Results run(const M& mat, bool center, bool scale, Engine* eng = null_rng(), Eigen::VectorXd* init = NULL) const {
-        Results output;
+    template<class Input, class Engine = std::mt19937_64>
+    Results<typename Input::Scalar> run(const Input& mat, bool center, bool scale, Engine* eng = null_rng(), VectorOf<typename Input::Scalar>* init = NULL) {
+        Results<typename Input::Scalar> output;
         auto stats = run(mat, center, scale, output.U, output.V, output.D, eng, init);
         output.converged = stats.first;
         output.iterations = stats.second;
@@ -548,7 +559,7 @@ class Irlba {
     /** 
      * Run IRLBA on an input matrix to perform an approximate SVD, see the `run()` method for more details.
      *
-     * @tparam M Matrix class,  most typically from the **Eigen** matrix manipulation library.
+     * @tparam Input Matrix class, typically from the **Eigen** matrix manipulation library.
      * However, other classes are also supported, see the other `run()` methods for details.
      * @tparam Engine A (pseudo-)random number generator class, returning a randomly sampled value when called as a functor with no arguments.
      *
@@ -560,9 +571,9 @@ class Irlba {
      *
      * @return A `Results` object containing the singular vectors and values, as well as some statistics on convergence.
      */
-    template<class M, class Engine = std::mt19937_64>
-    Results run(const M& mat, Engine* eng = null_rng(), Eigen::VectorXd* init = NULL) const {
-        Results output;
+    template<class Input, class Engine = std::mt19937_64>
+    Results<typename Input::Scalar> run(const Input& mat, Engine* eng = null_rng(), VectorOf<typename Input::Scalar>* init = NULL) {
+        Results<typename Input::Scalar> output;
         auto stats = run(mat, output.U, output.V, output.D, eng, init);
         output.converged = stats.first;
         output.iterations = stats.second;
diff --git a/include/irlba/lanczos.hpp b/include/irlba/lanczos.hpp
index abe9b29..e4ce40a 100644
--- a/include/irlba/lanczos.hpp
+++ b/include/irlba/lanczos.hpp
@@ -45,39 +45,43 @@ class LanczosBidiagonalization {
 public:
     /**
      * @brief Intermediate data structures to avoid repeated allocations.
+     *
+     * @tparam Scalar Scalar type for the intermediate results.
      */
+    template<class Scalar>
     struct Intermediates {
-        /**
-         * @tparam M Matrix class, most typically from the **Eigen** library.
-         *
-         * @param mat Instance of a matrix class `M`.
-         */
-        template<class M>
-        Intermediates(const M& mat) : F(mat.cols()), W_next(mat.rows()), orthog_tmp(mat.cols()) {}
-
         /**
          * Obtain the residual vector, see algorithm 2.1 of Baglama and Reichel (2005).
          *
          * @return Vector of residuals of length equal to the number of columns of `mat` in `run()`.
          */
-        const Eigen::VectorXd& residuals() const {
+        const VectorOf<Scalar>& residuals() const {
             return F;
         }
 
         /**
          * @cond
          */
-        Eigen::VectorXd F; 
-        Eigen::VectorXd W_next;
-        Eigen::VectorXd orthog_tmp;
+        Intermediates(size_t nrow, size_t ncol) : F(ncol), W_next(nrow), orthog_tmp(ncol) {}
+
+        VectorOf<Scalar> F; 
+        VectorOf<Scalar> W_next;
+        VectorOf<Scalar> orthog_tmp;
         /**
          * @endcond
          */
     };
 
-    template<class M>
-    Intermediates initialize(const M& mat) const {
-        return Intermediates(mat);
+    /**
+     * @param mat Input matrix.
+     *
+     * @return An `Intermediates` instance with appropriately allocated vectors.
+     *
+     * @tparam Input Input matrix class.
+     */
+    template<class Input>
+    Intermediates<typename Input::Scalar> initialize(const Input& mat) const {
+        return Intermediates<typename Input::Scalar>(mat.rows(), mat.cols());
     }
 
 public:
@@ -87,7 +91,7 @@ class LanczosBidiagonalization {
      * Support is provided for centering and scaling without modifying `mat`.
      * Protection against invariant subspaces is also implemented.
      *
-     * @tparam M Matrix class, most typically from the **Eigen** library.
+     * @tparam Input Input matrix class, most typically from the **Eigen** library.
      * See the `Irlba` documentation for a detailed description of the expected methods.
      * @tparam Engine A functor that, when called with no arguments, returns a random integer from a discrete uniform distribution.
      *
@@ -107,16 +111,18 @@ class LanczosBidiagonalization {
      * `W` is filled with orthonormal vectors, as is `V`.
      * `B` is filled with upper diagonal entries.
      */
-    template<class M, class Engine>
+    template<class Input, class Engine>
     void run(
-        const M& mat, 
-        Eigen::MatrixXd& W, 
-        Eigen::MatrixXd& V, 
-        Eigen::MatrixXd& B, 
+        const Input& mat, 
+        MatrixOf<typename Input::Scalar>& W, 
+        MatrixOf<typename Input::Scalar>& V, 
+        MatrixOf<typename Input::Scalar>& B, 
         Engine& eng, 
-        Intermediates& inter, 
+        Intermediates<typename Input::Scalar>& inter, 
         int start = 0) 
     const {
+        typedef MatrixOf<typename Input::Scalar> Matrix;
+        typedef VectorOf<typename Input::Scalar> Vector;
         const double eps = (epsilon < 0 ? std::pow(std::numeric_limits<double>::epsilon(), 0.8) : epsilon);
 
         int work = W.cols();
@@ -125,7 +131,7 @@ class LanczosBidiagonalization {
         auto& otmp = inter.orthog_tmp;
 
         F = V.col(start);
-        if constexpr(has_multiply_method<M>::value) {
+        if constexpr(has_multiply_method<Input, Vector>::value) {
             W_next.noalias() = mat * F;
         } else {
             mat.multiply(F, W_next);
@@ -145,7 +151,7 @@ class LanczosBidiagonalization {
 
         // The Lanczos iterations themselves.
         for (int j = start; j < work; ++j) {
-            if constexpr(has_adjoint_multiply_method<M>::value) {
+            if constexpr(has_adjoint_multiply_method<Input, Matrix>::value) {
                 F.noalias() = mat.adjoint() * W.col(j);
             } else {
                 mat.adjoint_multiply(W.col(j), F);
@@ -172,7 +178,7 @@ class LanczosBidiagonalization {
                 B(j, j) = S;
                 B(j, j + 1) = R_F;
 
-                if constexpr(has_multiply_method<M>::value) {
+                if constexpr(has_multiply_method<Input, Vector>::value) {
                     W_next.noalias() = mat * F;
                 } else {
                     mat.multiply(F, W_next);
diff --git a/include/irlba/utils.hpp b/include/irlba/utils.hpp
index d4132d4..4b3ece1 100644
--- a/include/irlba/utils.hpp
+++ b/include/irlba/utils.hpp
@@ -14,6 +14,30 @@
 
 namespace irlba {
 
+/**
+ * **Eigen** `Vector` type containing the real component of the specified scalar.
+ *
+ * @tparam Scalar Scalar type, typically a float or a complex number.
+ */
+template<class Scalar>
+using RealVectorOf = Eigen::Vector<typename Eigen::NumTraits<Scalar>::Real, Eigen::Dynamic>;
+
+/**
+ * **Eigen** `Vector` type containing the specified scalar.
+ *
+ * @tparam Scalar Scalar type, typically a float or a complex number.
+ */
+template<class Scalar>
+using VectorOf = Eigen::Vector<Scalar, Eigen::Dynamic>;
+
+/**
+ * Dense **Eigen** `Matrix` type containing the specified scalar.
+ *
+ * @tparam Scalar Scalar type, typically a float or a complex number.
+ */
+template<class Scalar>
+using MatrixOf = Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>;
+
 /**
  * Orthogonalize a vector against a set of orthonormal column vectors. 
  *
@@ -24,8 +48,11 @@ namespace irlba {
  *
  * @return `vec` is modified to contain `vec - mat0 * t(mat0) * vec`, where `mat0` is defined as the first `ncols` columns of `mat`.
  * This ensures that it is orthogonal to each column of `mat0`.
+ *
+ * @tparam Matrix Eigen `Matrix` type.
  */
-inline void orthogonalize_vector(const Eigen::MatrixXd& mat, Eigen::VectorXd& vec, size_t ncols, Eigen::VectorXd& tmp) {
+template<typename Matrix>
+inline void orthogonalize_vector(const Matrix& mat, VectorOf<typename Matrix::Scalar>& vec, size_t ncols, VectorOf<typename Matrix::Scalar>& tmp) {
     tmp.head(ncols).noalias() = mat.leftCols(ncols).adjoint() * vec;
     vec.noalias() -= mat.leftCols(ncols) * tmp.head(ncols);
     return;
@@ -34,16 +61,16 @@ inline void orthogonalize_vector(const Eigen::MatrixXd& mat, Eigen::VectorXd& ve
 /** 
  * Fill an **Eigen** vector with random normals via **aarand**.
  *
- * @param Vec Any **Eigen** vector class or equivalent proxy object.
+ * @param Vector Any **Eigen** vector class or equivalent proxy object.
  * @param Engine A (pseudo-)random number generator class that returns a random number when called with no arguments.
  *
- * @param vec Instance of a `Vec` class.
+ * @param vec Instance of a `Vector` class.
  * @param eng Instance of an `Engine` class.
  *
  * @return `vec` is filled with random draws from a standard normal distribution.
  */
-template<class Vec, class Engine>
-void fill_with_random_normals(Vec& vec, Engine& eng) {
+template<class Vector, class Engine>
+void fill_with_random_normals(Vector& vec, Engine& eng) {
     Eigen::Index i = 0;
     while (i < vec.size() - 1) {
         auto paired = aarand::standard_normal(eng);
@@ -56,6 +83,7 @@ void fill_with_random_normals(Vec& vec, Engine& eng) {
         auto paired = aarand::standard_normal(eng);
         vec[i] = paired.first;
     }
+
     return;
 }
 
@@ -77,7 +105,7 @@ struct ColumnVectorProxy {
 /** 
  * Fill a column of an **Eigen** matrix with random normals via **aarand**.
  *
- * @param Matrix Any **Eigen** matrix class or equivalent proxy object.
+ * @param Matrix Any **Eigen** matrix class.
  * @param Engine A (pseudo-)random number generator class that returns a random number when called with no arguments.
  *
  * @param mat Instance of a `Matrix` class.
@@ -145,10 +173,15 @@ struct ConvergenceTest {
      *
      * @param sv Vector of singular values.
      * @param residuals Vector of residuals for each singular value/vector.
+     * @param last Vector of singular values from the previous iteration.
      *
      * @return The number of singular values/vectors that have achieved convergence.
+     *
+     * @tparam RealVector An **Eigen** `Vector` class containing the singular values.
+     * @tparam Vector An **Eigen** `Vector` class for the residuals.
      */
-    int run(const Eigen::VectorXd& sv, const Eigen::VectorXd& residuals, const Eigen::VectorXd& last) const {
+    template<typename RealVector, typename Vector>
+    int run(const RealVector& sv, const Vector& residuals, const RealVector& last) const {
         int counter = 0;
         double Smax = *std::max_element(sv.begin(), sv.end());
         double svtol_actual = (svtol >= 0 ? svtol : tol);
@@ -176,23 +209,23 @@ constexpr std::mt19937_64* null_rng() { return NULL; }
 /**
  * @cond
  */
-template<class M, typename = int>
+template<typename M, typename V, typename = int>
 struct has_multiply_method {
     static constexpr bool value = false;
 };
 
-template<class M>
-struct has_multiply_method<M, decltype((void) (std::declval<M>() * std::declval<Eigen::VectorXd>()), 0)> {
+template<typename M, typename V>
+struct has_multiply_method<M, V, decltype((void) (std::declval<M>() * std::declval<V>()), 0)> {
     static constexpr bool value = true;
 };
 
-template<class M, typename = int>
+template<typename M, typename V, typename = int>
 struct has_adjoint_multiply_method {
     static constexpr bool value = false;
 };
 
-template<class M>
-struct has_adjoint_multiply_method<M, decltype((void) (std::declval<M>().adjoint() * std::declval<Eigen::VectorXd>()), 0)> {
+template<typename M, typename V>
+struct has_adjoint_multiply_method<M, V, decltype((void) (std::declval<M>().adjoint() * std::declval<V>()), 0)> {
     static constexpr bool value = true;
 };
 
@@ -203,7 +236,10 @@ struct has_realize_method {
 
 template<class M>
 struct has_realize_method<M, decltype((void) std::declval<M>().realize(), 0)> {
-    static constexpr bool value = std::is_same<decltype(std::declval<M>().realize()), Eigen::MatrixXd>::value;
+    // Only is 'true' if the realize() method actually returns a realization as
+    // an Eigen::Matrix. This protects against future problems if Eigen itself
+    // defines a realize() method on its classes.
+    static constexpr bool value = std::is_same<decltype(std::declval<M>().realize()), Eigen::Matrix<typename M::Scalar, Eigen::Dynamic, Eigen::Dynamic> >::value;
 };
 /**
  * @endcond
diff --git a/include/irlba/wrappers.hpp b/include/irlba/wrappers.hpp
index f0eb717..d76bf19 100644
--- a/include/irlba/wrappers.hpp
+++ b/include/irlba/wrappers.hpp
@@ -17,11 +17,12 @@
  *
  * - `mat.rows()`, returning the number of rows.
  * - `mat.cols()`, returning the number of columns.
+ * - a `Scalar` typedef, returning the type of scalar in the modified matrix.
  * - `mat.multiply(rhs, out)`, which computes the matrix product `mat * rhs` and stores it in `out`.
- * `rhs` should be an `Eigen::VectorXd` (or an expression equivalent) while `out` should be a `Eigen::VectorXd`.
+ * `rhs` and `out` should be `Eigen::Vector<Scalar, Eigen::Dynamic>` instances; alternatively, `rhs` may also be an expression equivalent.
  * - `mat.adjoint_multiply(rhs, out)`, which computes the matrix product `mat.adjoint() * rhs` and stores it in `out`.
- * `rhs` should be an `Eigen::VectorXd` (or an expression equivalent) while `out` should be a `Eigen::VectorXd`.
- * - `mat.realize()`, which returns an `Eigen::MatrixXd` containing the matrix with all modifications applied.
+ * `rhs` and `out` should be `Eigen::Vector<Scalar, Eigen::Dynamic>` instances; alternatively, `rhs` may also be an expression equivalent.
+ * - `mat.realize()`, which returns an `Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>` object containing the matrix with all modifications applied.
  */
 
 namespace irlba {
@@ -29,19 +30,24 @@ namespace irlba {
 /**
  * @brief Wrapper for a centered matrix.
  *
- * @tparam Matrix An **Eigen** matrix class - or alternatively, a wrapper class around such a class.
+ * @tparam Input An **Eigen** `Matrix` class, or a wrapper around such a class.
  * 
  * This modification involves centering all columns, i.e., subtracting the mean of each column from the values of that column.
  * Naively doing such an operation would involve loss of sparsity, which we avoid by deferring the subtraction into the subspace defined by `rhs`.
  */
-template<class Matrix>
+template<typename Input>
 struct Centered {
+    /**
+     * Type of scalar in the centered matrix.
+     */
+    typedef typename Input::Scalar Scalar;
+
     /**
      * @param m Underlying matrix to be column-centered.
      * @param c Vector of length equal to the number of columns of `m`,
      * containing the value to subtract from each column.
      */
-    Centered(const Matrix* m, const Eigen::VectorXd* c) : mat(m), center(c) {}
+    Centered(const Input* m, const Eigen::Vector<Scalar, Eigen::Dynamic>* c) : mat(m), center(c) {}
 
     /**
      * @return Number of rows in the matrix.
@@ -54,7 +60,7 @@ struct Centered {
     auto cols() const { return mat->cols(); }
 
     /**
-     * @tparam Right An `Eigen::VectorXd` or equivalent expression.
+     * @tparam Right An Eigen `Vector` class or equivalent expression.
      *
      * @param[in] rhs The right-hand side of the matrix product.
      * This should be a vector or have only one column.
@@ -63,8 +69,8 @@ struct Centered {
      * @return `out` is filled with the product of this matrix and `rhs`.
      */
     template<class Right>
-    void multiply(const Right& rhs, Eigen::VectorXd& out) const {
-        if constexpr(has_multiply_method<Matrix>::value) {
+    void multiply(const Right& rhs, Eigen::Vector<Scalar, Eigen::Dynamic>& out) const {
+        if constexpr(has_multiply_method<Input, Right>::value) {
             out.noalias() = *mat * rhs;
         } else {
             mat->multiply(rhs, out);
@@ -78,7 +84,7 @@ struct Centered {
     }
 
     /**
-     * @tparam Right An `Eigen::VectorXd` or equivalent expression.
+     * @tparam Right An Eigen `Vector` class or equivalent expression.
      *
      * @param[in] rhs The right-hand side of the matrix product.
      * This should be a vector or have only one column.
@@ -87,8 +93,8 @@ struct Centered {
      * @return `out` is filled with the product of the transpose of this matrix and `rhs`.
      */
     template<class Right>
-    void adjoint_multiply(const Right& rhs, Eigen::VectorXd& out) const {
-        if constexpr(has_adjoint_multiply_method<Matrix>::value) {
+    void adjoint_multiply(const Right& rhs, Eigen::Vector<Scalar, Eigen::Dynamic>& out) const {
+        if constexpr(has_adjoint_multiply_method<Input, Right>::value) {
             out.noalias() = mat->adjoint() * rhs;
         } else {
             mat->adjoint_multiply(rhs, out);
@@ -103,8 +109,8 @@ struct Centered {
      * @return A realized version of the centered matrix,
      * where the centering has been explicitly applied.
      */
-    Eigen::MatrixXd realize() const {
-        auto subtractor = [&](Eigen::MatrixXd& m) -> void {
+    Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> realize() const {
+        auto subtractor = [&](auto& m) -> void {
             for (Eigen::Index c = 0; c < m.cols(); ++c) {
                 for (Eigen::Index r = 0; r < m.rows(); ++r) {
                     m(r, c) -= (*center)[c];
@@ -112,38 +118,43 @@ struct Centered {
             }
         };
 
-        if constexpr(has_realize_method<Matrix>::value) {
-            Eigen::MatrixXd output = mat->realize();
+        if constexpr(has_realize_method<Input>::value) {
+            auto output = mat->realize(); // Scalars are the same, so the output should also be the same type of Eigen matrix.
             subtractor(output);
             return output;
         } else {
-            Eigen::MatrixXd output(*mat);
+            Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> output(*mat);
             subtractor(output);
             return output;
         }
     }
 
 private:
-    const Matrix* mat;
-    const Eigen::VectorXd* center;
+    const Input* mat;
+    const Eigen::Vector<Scalar, Eigen::Dynamic>* center;
 };
 
 /**
  * @brief Wrapper for a scaled matrix.
  *
- * @tparam Matrix An **Eigen** matrix class - or alternatively, a wrapper class around such a class.
+ * @tparam Input An **Eigen** matrix class - or alternatively, a wrapper class around such a class.
  * 
  * This modification involves scaling all columns, i.e., dividing the values of each column by the standard deviation of that column to achieve unit variance.
  * Naively doing such an operation would involve a copy of the matrix, which we avoid by deferring the scaling into the subspace defined by `rhs`.
  */
-template<class Matrix>
+template<typename Input>
 struct Scaled {
+    /**
+     * Type of scalar in the centered matrix.
+     */
+    typedef typename Input::Scalar Scalar;
+
     /**
      * @param m Underlying matrix to be column-scaled.
      * @param s Vector of length equal to the number of columns of `m`,
      * containing the value to scale each column.
      */
-    Scaled(const Matrix* m, const Eigen::VectorXd* s) : mat(m), scale(s) {}
+    Scaled(const Input* m, const Eigen::Vector<Scalar, Eigen::Dynamic>* s) : mat(m), scale(s) {}
 
     /**
      * @return Number of rows in the matrix.
@@ -156,7 +167,7 @@ struct Scaled {
     auto cols() const { return mat->cols(); }
 
     /**
-     * @tparam Right An `Eigen::VectorXd` or equivalent expression.
+     * @tparam Right An Eigen `Vector` or equivalent expression.
      *
      * @param[in] rhs The right-hand side of the matrix product.
      * This should be a vector or have only one column.
@@ -165,8 +176,8 @@ struct Scaled {
      * @return `out` is filled with the product of this matrix and `rhs`.
      */
     template<class Right>
-    void multiply(const Right& rhs, Eigen::VectorXd& out) const {
-        if constexpr(has_multiply_method<Matrix>::value) {
+    void multiply(const Right& rhs, Eigen::Vector<Scalar, Eigen::Dynamic>& out) const {
+        if constexpr(has_multiply_method<Input, Right>::value) {
             out.noalias() = *mat * rhs.cwiseQuotient(*scale);
         } else {
             mat->multiply(rhs.cwiseQuotient(*scale), out);
@@ -175,7 +186,7 @@ struct Scaled {
     }
 
     /**
-     * @tparam Right An `Eigen::VectorXd` or equivalent expression.
+     * @tparam Right An Eigen `Vector` or equivalent expression.
      *
      * @param[in] rhs The right-hand side of the matrix product.
      * This should be a vector or have only one column.
@@ -184,8 +195,8 @@ struct Scaled {
      * @return `out` is filled with the product of the transpose of this matrix and `rhs`.
      */
     template<class Right>
-    void adjoint_multiply(const Right& rhs, Eigen::VectorXd& out) const {
-        if constexpr(has_adjoint_multiply_method<Matrix>::value) {
+    void adjoint_multiply(const Right& rhs, Eigen::Vector<Scalar, Eigen::Dynamic>& out) const {
+        if constexpr(has_adjoint_multiply_method<Input, Right>::value) {
             out.noalias() = mat->adjoint() * rhs;
         } else {
             mat->adjoint_multiply(rhs, out);
@@ -198,8 +209,8 @@ struct Scaled {
      * @return A realized version of the scaled matrix,
      * where the scaling has been explicitly applied.
      */
-    Eigen::MatrixXd realize() const {
-        auto scaler = [&](Eigen::MatrixXd& m) -> void {
+    Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> realize() const {
+        auto scaler = [&](auto& m) -> void {
             for (Eigen::Index c = 0; c < m.cols(); ++c) {
                 for (Eigen::Index r = 0; r < m.rows(); ++r) {
                     m(r, c) /= (*scale)[c];
@@ -207,20 +218,20 @@ struct Scaled {
             }
         };
 
-        if constexpr(has_realize_method<Matrix>::value) {
-            Eigen::MatrixXd output = mat->realize();
+        if constexpr(has_realize_method<Input>::value) {
+            auto output = mat->realize(); // Scalars are the same, so the output should also be the same type of Eigen matrix.
             scaler(output);
             return output;
         } else {
-            Eigen::MatrixXd output(*mat);
+            Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> output(*mat);
             scaler(output);
             return output;
         }
     }
 
 private:
-    const Matrix* mat;
-    const Eigen::VectorXd* scale;
+    const Input* mat;
+    const Eigen::Vector<Scalar, Eigen::Dynamic>* scale;
 };
 
 }
diff --git a/tests/src/NormalSampler.h b/tests/src/NormalSampler.h
index 8b1785a..8a5a83e 100644
--- a/tests/src/NormalSampler.h
+++ b/tests/src/NormalSampler.h
@@ -22,6 +22,17 @@ inline Eigen::MatrixXd create_random_matrix(size_t nr, size_t nc, int seed = 42)
     return A;
 }
 
+inline Eigen::MatrixXcd create_random_complex_matrix(size_t nr, size_t nc, int seed = 42) {
+    Eigen::MatrixXcd A(nr, nc);
+    NormalSampler norm(seed); 
+    for (size_t i = 0; i < nc; ++i) {
+        for (size_t j = 0; j < nr; ++j) {
+            A(j, i) = std::complex<double>(norm(), norm());
+        }
+    }
+    return A;
+}
+
 inline Eigen::VectorXd create_random_vector(size_t n, int seed = 50) {
     NormalSampler norm(seed);
     Eigen::VectorXd output(n);
diff --git a/tests/src/compare.h b/tests/src/compare.h
index 6436c38..fa79f77 100644
--- a/tests/src/compare.h
+++ b/tests/src/compare.h
@@ -4,12 +4,13 @@
 #include <gtest/gtest.h>
 #include "Eigen/Dense"
 
-inline bool same_same(double left, double right, double tol) {
+template<typename T>
+inline bool same_same(T left, T right, double tol) {
     return std::abs(left - right) <= (std::abs(left) + std::abs(right)) * tol;
 }
 
-template<bool centered=false>
-void expect_equal_column_vectors(const Eigen::MatrixXd& left, const Eigen::MatrixXd& right, double tol=1e-8) {
+template<bool centered=false, typename LeftType, typename RightType>
+void expect_equal_column_vectors(const LeftType& left, const RightType& right, double tol=1e-8) {
     ASSERT_EQ(left.cols(), right.cols());
     ASSERT_EQ(left.rows(), right.rows());
 
@@ -30,7 +31,8 @@ void expect_equal_column_vectors(const Eigen::MatrixXd& left, const Eigen::Matri
     return;
 }
 
-inline void expect_equal_matrix(const Eigen::MatrixXd& left, const Eigen::MatrixXd& right, double tol=1e-8) {
+template<typename LeftType, typename RightType>
+inline void expect_equal_matrix(const LeftType& left, const RightType& right, double tol=1e-8) {
     ASSERT_EQ(left.cols(), right.cols());
     ASSERT_EQ(left.rows(), right.rows());
     for (size_t i = 0; i < left.cols(); ++i) {
@@ -40,7 +42,8 @@ inline void expect_equal_matrix(const Eigen::MatrixXd& left, const Eigen::Matrix
     }
 }
 
-inline void expect_equal_vectors(const Eigen::VectorXd& left, const Eigen::VectorXd& right, double tol=1e-8) {
+template<typename LeftType, typename RightType>
+inline void expect_equal_vectors(const LeftType& left, const RightType& right, double tol=1e-8) {
     ASSERT_EQ(left.size(), right.size());
     for (size_t i = 0; i < left.size(); ++i) {
         EXPECT_TRUE(same_same(left[i], right[i], tol));
diff --git a/tests/src/irlba.cpp b/tests/src/irlba.cpp
index af25943..ffc4fa0 100644
--- a/tests/src/irlba.cpp
+++ b/tests/src/irlba.cpp
@@ -25,6 +25,22 @@ TEST(IrlbaTest, Exact) {
     expect_equal_column_vectors(res.V, svd.matrixV().leftCols(rank), 1e-8);
 }
 
+TEST(IrlbaComplexTest, Exact) {
+    // For the test, the key is that rank + workspace > min(nr, nc), in which
+    // case we can be pretty confident of getting a near-exact match of the
+    // true SVD. Otherwise it's more approximate and the test is weaker.
+    int rank = 5;
+    auto A = create_random_complex_matrix(20, 10);
+
+    irlba::Irlba irb;
+    auto res = irb.set_number(rank).run(A);
+
+    Eigen::BDCSVD svd(A, Eigen::ComputeThinU | Eigen::ComputeThinV);
+    expect_equal_vectors(res.D, svd.singularValues().head(rank), 1e-8);
+    expect_equal_column_vectors(res.U, svd.matrixU().leftCols(rank), 1e-8);
+    expect_equal_column_vectors(res.V, svd.matrixV().leftCols(rank), 1e-8);
+}
+
 class IrlbaTester : public ::testing::TestWithParam<std::tuple<int, int, int> > {
 protected:
     template<class Param>
diff --git a/tests/src/utils.cpp b/tests/src/utils.cpp
index 39095b0..6f287a5 100644
--- a/tests/src/utils.cpp
+++ b/tests/src/utils.cpp
@@ -1,5 +1,6 @@
 #include <gtest/gtest.h>
 #include "irlba/utils.hpp"
+#include "irlba/wrappers.hpp"
 #include "Eigen/Dense"
 #include "Eigen/Sparse"
 #include "NormalSampler.h"
@@ -95,12 +96,18 @@ TEST(UtilsTest, FillNormals) {
 
 struct dummy_class {
     Eigen::MatrixXd realize() { return Eigen::MatrixXd(); }
+    typedef double Scalar;
 };
 
 TEST(UtilsTest, ConvertibleCheck) {
-    EXPECT_FALSE(irlba::has_realize_method<Eigen::MatrixXd>::value);
-    EXPECT_FALSE(irlba::has_realize_method<Eigen::SparseMatrix<double> >::value);
-    EXPECT_TRUE(irlba::has_realize_method<dummy_class>::value);
+    bool const realize_xd_xd = irlba::has_realize_method<Eigen::MatrixXd>::value;
+    EXPECT_FALSE(realize_xd_xd);
+    bool const realize_sp_sp = irlba::has_realize_method<Eigen::SparseMatrix<double> >::value;
+    EXPECT_FALSE(realize_sp_sp);
+    bool const realize_dummy = irlba::has_realize_method<dummy_class>::value;
+    EXPECT_TRUE(realize_dummy);
+    bool const realize_scaled = irlba::has_realize_method<irlba::Scaled<Eigen::MatrixXd> >::value;
+    EXPECT_TRUE(realize_scaled);
 }
 
 TEST(UtilsTest, NormalSampler) {
diff --git a/tests/src/wrappers.cpp b/tests/src/wrappers.cpp
index 9744aca..b420fe3 100644
--- a/tests/src/wrappers.cpp
+++ b/tests/src/wrappers.cpp
@@ -17,14 +17,14 @@ TEST(WrapperTest, Centering) {
         auto C = create_random_vector(10, 1234);
         Eigen::VectorXd output(20);
         centered.multiply(C, output);
-        expect_equal_matrix(ref * C, output);
+        expect_equal_matrix<Eigen::MatrixXd>(ref * C, output);
     }
 
     {
         auto C = create_random_vector(20, 1234);
         Eigen::VectorXd output(10);
         centered.adjoint_multiply(C, output);
-        expect_equal_matrix(ref.adjoint() * C, output);
+        expect_equal_matrix<Eigen::MatrixXd>(ref.adjoint() * C, output);
     }
 }
 
@@ -41,13 +41,13 @@ TEST(WrapperTest, Scaling) {
         auto C = create_random_vector(10, 1234);
         Eigen::VectorXd output(20);
         centered.multiply(C, output);
-        expect_equal_matrix(ref * C, output);
+        expect_equal_matrix<Eigen::MatrixXd>(ref * C, output);
     }
 
     {
         auto C = create_random_vector(20, 1234);
         Eigen::VectorXd output(10);
         centered.adjoint_multiply(C, output);
-        expect_equal_matrix(ref.adjoint() * C, output);
+        expect_equal_matrix<Eigen::MatrixXd>(ref.adjoint() * C, output);
     }
 }