use pade approx for matrix_exp_2x2 when cosh & sinh ops overflow

Author: Yi Zhang

stan/math/prim/fun/matrix_exp_2x2.hpp

@@ -47,7 +47,12 @@ matrix_exp_2x2(const EigMat& A) {
   B(1, 0) = c * Two_exp_sinh;
   B(1, 1) = exp_half_a_plus_d * (delta_cosh - ad_sinh_half_delta);
 
-  return B / delta;
+  // use pade approximation if cosh & sinh ops overflow to NaN
+  if ((B.array() != B.array()).any()) {
+    return matrix_exp_pade(A);
+  } else {
+    return B / delta;    
+  }
 }
 
 }  // namespace math

test/unit/math/prim/fun/matrix_exp_2x2_test.cpp

@@ -20,3 +20,12 @@ TEST(MathMatrixPrimMat, matrix_exp_2x2_2x2_2) {
 
   EXPECT_MATRIX_FLOAT_EQ(m2, stan::math::matrix_exp_2x2(m1));
 }
+
+TEST(MathMatrixPrimMat, matrix_exp_2x2_2x2_overflow) {
+  // example from issue https://github.com/stan-dev/math/issues/2615
+  Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> m1(2, 2), m2(2, 2);
+  m1 << -3.3228,0.533302,1.2242,-4.04844;
+  double t = 800.0;
+  m2 << 0.0, 0.0, 0.0, 0.0;
+  EXPECT_MATRIX_FLOAT_EQ(m2, stan::math::matrix_exp_2x2(m1 * t));
+}