Explicit return types, improved testing

Author: WardBrian

stan/math/prim/fun/eigenvalues.hpp

@@ -17,7 +17,8 @@ namespace math {
  */
 template <typename EigMat, require_eigen_matrix_dynamic_t<EigMat>* = nullptr,
           require_not_vt_complex<EigMat>* = nullptr>
-inline auto eigenvalues(const EigMat& m) {
+inline Eigen::Matrix<complex_return_t<value_type_t<EigMat>>, 1, -1> eigenvalues(
+    const EigMat& m) {
   using PlainMat = plain_type_t<EigMat>;
   const PlainMat& m_eval = m;
   check_nonzero_size("eigenvalues", "m", m_eval);
@@ -30,15 +31,17 @@ inline auto eigenvalues(const EigMat& m) {
 /**
  * Return the eigenvalues of a (complex-valued) matrix
  *
- * @tparam EigMat type of complex matrix argument
+ * @tparam EigCplxMat type of complex matrix argument
  * @param[in] m matrix to find the eigenvectors of. Must be square and have a
  * non-zero size.
  * @return a complex-valued column vector with entries the eigenvectors of `m`
  */
-template <typename EigMat, require_eigen_matrix_dynamic_t<EigMat>* = nullptr,
-          require_vt_complex<EigMat>* = nullptr>
-inline auto eigenvalues(const EigMat& m) {
-  using PlainMat = Eigen::Matrix<scalar_type_t<EigMat>, -1, -1>;
+template <typename EigCplxMat,
+          require_eigen_matrix_dynamic_t<EigCplxMat>* = nullptr,
+          require_vt_complex<EigCplxMat>* = nullptr>
+inline Eigen::Matrix<complex_return_t<value_type_t<EigCplxMat>>, 1, -1>
+eigenvalues(const EigCplxMat& m) {
+  using PlainMat = Eigen::Matrix<scalar_type_t<EigCplxMat>, -1, -1>;
   const PlainMat& m_eval = m;
   check_nonzero_size("eigenvalues", "m", m_eval);
   check_square("eigenvalues", "m", m_eval);

stan/math/prim/fun/eigenvectors.hpp

@@ -18,7 +18,8 @@ namespace math {
  */
 template <typename EigMat, require_eigen_matrix_dynamic_t<EigMat>* = nullptr,
           require_not_vt_complex<EigMat>* = nullptr>
-inline auto eigenvectors(const EigMat& m) {
+inline Eigen::Matrix<complex_return_t<value_type_t<EigMat>>, -1, -1>
+eigenvectors(const EigMat& m) {
   using PlainMat = plain_type_t<EigMat>;
   const PlainMat& m_eval = m;
   check_nonzero_size("eigenvectors", "m", m_eval);
@@ -31,16 +32,18 @@ inline auto eigenvectors(const EigMat& m) {
 /**
  * Return the eigenvectors of a (complex-valued) matrix
  *
- * @tparam EigMat type of complex matrix argument
+ * @tparam EigCplxMat type of complex matrix argument
  * @param[in] m matrix to find the eigenvectors of. Must be square and have a
  * non-zero size.
  * @return a complex-valued matrix with columns representing the eigenvectors of
  * `m`
  */
-template <typename EigMat, require_eigen_matrix_dynamic_t<EigMat>* = nullptr,
-          require_vt_complex<EigMat>* = nullptr>
-inline auto eigenvectors(const EigMat& m) {
-  using PlainMat = Eigen::Matrix<scalar_type_t<EigMat>, -1, -1>;
+template <typename EigCplxMat,
+          require_eigen_matrix_dynamic_t<EigCplxMat>* = nullptr,
+          require_vt_complex<EigCplxMat>* = nullptr>
+inline Eigen::Matrix<complex_return_t<value_type_t<EigCplxMat>>, -1, -1>
+eigenvectors(const EigCplxMat& m) {
+  using PlainMat = Eigen::Matrix<scalar_type_t<EigCplxMat>, -1, -1>;
   const PlainMat& m_eval = m;
   check_nonzero_size("eigenvectors", "m", m_eval);
   check_square("eigenvectors", "m", m_eval);

test/unit/math/mix/fun/eigenvectors_test.cpp

@@ -31,17 +31,17 @@ TEST(mathMixFun, eigenvectorsComplex) {
 
 template <typename T>
 void expect_identity_matrix(const T& x) {
-  EXPECT_EQUAL(x.rows(), x.cols());
+  EXPECT_EQ(x.rows(), x.cols());
   for (int j = 0; j < x.cols(); ++j) {
     for (int i = 0; i < x.rows(); ++i) {
-      EXPECT_NEAR(i == j ? 1 : 0, x(i, j), 1e-6);
+      EXPECT_NEAR(i == j ? 1 : 0, stan::math::value_of_rec(x(i, j)), 1e-6);
     }
   }
 }
 
 template <typename T>
 void expectEigenvectorsId() {
-  for (const auto& m_d : stan::test::square_test_matrices(0, 2)) {
+  for (const auto& m_d : stan::test::square_test_matrices(1, 2)) {
     Eigen::Matrix<T, -1, -1> m(m_d);
     auto vecs = eigenvectors(m).eval();
     auto vals = eigenvalues(m).eval();
@@ -50,7 +50,21 @@ void expectEigenvectorsId() {
   }
 }
 
-// THESE TESTS USED TO WORK STANDALONE
+template <typename T>
+void expectComplexEigenvectorsId() {
+  Eigen::Matrix<std::complex<T>, -1, -1> c22(2, 2);
+  c22 << stan::math::to_complex(T(0), T(-1)),
+      stan::math::to_complex(T(0), T(0)), stan::math::to_complex(T(2), T(0)),
+      stan::math::to_complex(T(4), T(0));
+  auto eigenvalues = stan::math::eigenvalues(c22);
+  auto eigenvectors = stan::math::eigenvectors(c22);
+
+  auto I = (eigenvectors.inverse() * c22 * eigenvectors
+            * eigenvalues.asDiagonal().inverse())
+               .real();
+
+  expect_identity_matrix(I);
+}
 
 TEST(mathMixFun, eigenvectorsId) {
   using d_t = double;
@@ -60,9 +74,15 @@ TEST(mathMixFun, eigenvectorsId) {
   using fv_t = stan::math::fvar<stan::math::var>;
   using ffv_t = stan::math::fvar<fv_t>;
 
-  // expectEigenvectorsId<v_t>();
-  // expectEigenvectorsId<fd_t>();
-  // expectEigenvectorsId<ffd_t>();
-  // expectEigenvectorsId<fv_t>();
-  // expectEigenvectorsId<ffv_t>();
+  expectEigenvectorsId<v_t>();
+  expectEigenvectorsId<fd_t>();
+  expectEigenvectorsId<ffd_t>();
+  expectEigenvectorsId<fv_t>();
+  expectEigenvectorsId<ffv_t>();
+
+  expectComplexEigenvectorsId<v_t>();
+  expectComplexEigenvectorsId<fd_t>();
+  expectComplexEigenvectorsId<ffd_t>();
+  expectComplexEigenvectorsId<fv_t>();
+  expectComplexEigenvectorsId<ffv_t>();
 }