Merge commit '3dc38852' into HEAD

stan/math/fwd/mat/fun/Eigen_NumTraits.hpp

@@ -52,14 +52,12 @@ struct NumTraits<stan::math::fvar<T>> : GenericNumTraits<stan::math::fvar<T>> {
   static int digits10() { return std::numeric_limits<double>::digits10; }
 };
 
-namespace internal {
-
 /**
  * Scalar product traits specialization for Eigen for forward-mode
  * autodiff variables.
  */
-template <typename T>
-struct scalar_product_traits<stan::math::fvar<T>, double> {
+template <typename T, typename BinaryOp>
+struct ScalarBinaryOpTraits<stan::math::fvar<T>, double, BinaryOp> {
   typedef stan::math::fvar<T> ReturnType;
 };
 
@@ -67,11 +65,10 @@ struct scalar_product_traits<stan::math::fvar<T>, double> {
  * Scalar product traits specialization for Eigen for forward-mode
  * autodiff variables.
  */
-template <typename T>
-struct scalar_product_traits<double, stan::math::fvar<T>> {
+template <typename T, typename BinaryOp>
+struct ScalarBinaryOpTraits<double, stan::math::fvar<T>, BinaryOp> {
   typedef stan::math::fvar<T> ReturnType;
 };
-}  // namespace internal
 
 }  // namespace Eigen
 #endif

stan/math/fwd/scal/fun/log_diff_exp.hpp

@@ -12,9 +12,11 @@ namespace math {
 
 template <typename T>
 inline fvar<T> log_diff_exp(const fvar<T>& x1, const fvar<T>& x2) {
-  using std::exp;
-  if (x1.val_ <= x2.val_)
+  if (x1.val_ <= x2.val_) {
+    if (x1.val_ < INFTY && x1.val_ == x2.val_)
+      return fvar<T>(NEGATIVE_INFTY, NOT_A_NUMBER);
     return fvar<T>(NOT_A_NUMBER, NOT_A_NUMBER);
+  }
   return fvar<T>(
       log_diff_exp(x1.val_, x2.val_),
       -(x1.d_ / expm1(x2.val_ - x1.val_) + x2.d_ / expm1(x1.val_ - x2.val_)));
@@ -22,17 +24,24 @@ inline fvar<T> log_diff_exp(const fvar<T>& x1, const fvar<T>& x2) {
 
 template <typename T1, typename T2>
 inline fvar<T2> log_diff_exp(const T1& x1, const fvar<T2>& x2) {
-  using std::exp;
-  if (x1 <= x2.val_)
+  if (x1 <= x2.val_) {
+    if (x1 < INFTY && x1 == x2.val_)
+      return fvar<T2>(NEGATIVE_INFTY, x2.d_ * NEGATIVE_INFTY);
     return fvar<T2>(NOT_A_NUMBER, NOT_A_NUMBER);
+  }
   return fvar<T2>(log_diff_exp(x1, x2.val_), -x2.d_ / expm1(x1 - x2.val_));
 }
 
 template <typename T1, typename T2>
 inline fvar<T1> log_diff_exp(const fvar<T1>& x1, const T2& x2) {
-  using std::exp;
-  if (x1.val_ <= x2)
+  if (x1.val_ <= x2) {
+    if (x1.val_ < INFTY && x1.val_ == x2) {
+      if (x2 == NEGATIVE_INFTY)
+        return fvar<T1>(NEGATIVE_INFTY, x1.d_);
+      return fvar<T1>(NEGATIVE_INFTY, x1.d_ * INFTY);
+    }
     return fvar<T1>(NOT_A_NUMBER, NOT_A_NUMBER);
+  }
   return fvar<T1>(log_diff_exp(x1.val_, x2), -x1.d_ / expm1(x2 - x1.val_));
 }
 }  // namespace math

stan/math/fwd/scal/fun/log_sum_exp.hpp

@@ -3,6 +3,7 @@
 
 #include <stan/math/fwd/meta.hpp>
 #include <stan/math/fwd/core.hpp>
+#include <stan/math/prim/scal/fun/constants.hpp>
 #include <stan/math/prim/scal/fun/log_sum_exp.hpp>
 
 namespace stan {
@@ -19,12 +20,16 @@ inline fvar<T> log_sum_exp(const fvar<T>& x1, const fvar<T>& x2) {
 template <typename T>
 inline fvar<T> log_sum_exp(double x1, const fvar<T>& x2) {
   using std::exp;
+  if (x1 == NEGATIVE_INFTY)
+    return fvar<T>(x2.val_, x2.d_);
   return fvar<T>(log_sum_exp(x1, x2.val_), x2.d_ / (exp(x1 - x2.val_) + 1));
 }
 
 template <typename T>
 inline fvar<T> log_sum_exp(const fvar<T>& x1, double x2) {
   using std::exp;
+  if (x2 == NEGATIVE_INFTY)
+    return fvar<T>(x1.val_, x1.d_);
   return fvar<T>(log_sum_exp(x1.val_, x2), x1.d_ / (1 + exp(x2 - x1.val_)));
 }
 

stan/math/opencl/matrix_cl_view.hpp

@@ -70,7 +70,9 @@ inline const matrix_cl_view transpose(const matrix_cl_view view) {
  */
 inline const matrix_cl_view invert(const matrix_cl_view view) {
   typedef typename std::underlying_type<matrix_cl_view>::type underlying;
-  return static_cast<matrix_cl_view>(~static_cast<underlying>(view));
+  return static_cast<matrix_cl_view>(
+      static_cast<underlying>(matrix_cl_view::Entire)
+      & ~static_cast<underlying>(view));
 }
 
 /**

stan/math/prim/mat/fun/LDLT_factor.hpp

@@ -99,19 +99,11 @@ class LDLT_factor {
     ldltP_->solveInPlace(invA);
   }
 
-#if EIGEN_VERSION_AT_LEAST(3, 3, 0)
   template <typename Rhs>
   inline const Eigen::Solve<ldlt_t, Rhs> solve(
       const Eigen::MatrixBase<Rhs>& b) const {
     return ldltP_->solve(b);
   }
-#else
-  template <typename Rhs>
-  inline const Eigen::internal::solve_retval<ldlt_t, Rhs> solve(
-      const Eigen::MatrixBase<Rhs>& b) const {
-    return ldltP_->solve(b);
-  }
-#endif
 
   inline matrix_t solveRight(const matrix_t& B) const {
     return ldltP_->solve(B.transpose()).transpose();

stan/math/prim/mat/fun/cholesky_decompose.hpp

@@ -6,10 +6,7 @@
 #include <stan/math/prim/mat/err/check_square.hpp>
 #include <stan/math/prim/mat/err/check_symmetric.hpp>
 #ifdef STAN_OPENCL
-#include <stan/math/opencl/opencl_context.hpp>
-#include <stan/math/opencl/err/check_symmetric.hpp>
-#include <stan/math/opencl/cholesky_decompose.hpp>
-#include <stan/math/opencl/copy.hpp>
+#include <stan/math/opencl/opencl.hpp>
 #endif
 
 #include <cmath>

stan/math/prim/mat/fun/log_sum_exp.hpp

@@ -25,6 +25,8 @@ namespace math {
 template <int R, int C>
 double log_sum_exp(const Eigen::Matrix<double, R, C>& x) {
   const double max = x.maxCoeff();
+  if (!std::isfinite(max))
+    return max;
   return max + std::log((x.array() - max).exp().sum());
 }
 

stan/math/prim/mat/fun/multiply.hpp

@@ -5,6 +5,9 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/arr/err/check_matching_sizes.hpp>
 #include <stan/math/prim/mat/err/check_multiplicable.hpp>
+#ifdef STAN_OPENCL
+#include <stan/math/opencl/opencl.hpp>
+#endif
 #include <type_traits>
 
 namespace stan {
@@ -55,7 +58,19 @@ template <int R1, int C1, int R2, int C2, typename T1, typename T2,
 inline Eigen::Matrix<return_type_t<T1, T2>, R1, C2> multiply(
     const Eigen::Matrix<T1, R1, C1>& m1, const Eigen::Matrix<T2, R2, C2>& m2) {
   check_multiplicable("multiply", "m1", m1, "m2", m2);
+#ifdef STAN_OPENCL
+  if (m1.rows() * m1.cols() * m2.cols()
+      > opencl_context.tuning_opts().multiply_dim_prod_worth_transfer) {
+    matrix_cl<double> m1_cl(m1);
+    matrix_cl<double> m2_cl(m2);
+    matrix_cl<double> m3_cl = m1_cl * m2_cl;
+    return from_matrix_cl(m3_cl);
+  } else {
+    return m1 * m2;
+  }
+#else
   return m1 * m2;
+#endif
 }
 
 /**

stan/math/prim/scal/fun/log_diff_exp.hpp

@@ -11,16 +11,16 @@ namespace math {
 
 /**
  * The natural logarithm of the difference of the natural exponentiation
- * of x1 and the natural exponentiation of x2
+ * of x and the natural exponentiation of y
  *
- * This function is only defined for x<0
+ * This function is only defined for x >= y
  *
  *
    \f[
    \mbox{log\_diff\_exp}(x, y) =
    \begin{cases}
-     \textrm{NaN} & \mbox{if } x \leq y\\
-     \ln(\exp(x)-\exp(y)) & \mbox{if } x > y \\[6pt]
+     \textrm{NaN} & \mbox{if } x < y\\
+     \ln(\exp(x)-\exp(y)) & \mbox{if } x \geq y \\[6pt]
      \textrm{NaN} & \mbox{if } x = \textrm{NaN or } y = \textrm{NaN}
    \end{cases}
    \f]
@@ -47,7 +47,7 @@ namespace math {
 template <typename T1, typename T2>
 inline return_type_t<T1, T2> log_diff_exp(const T1 x, const T2 y) {
   if (x <= y)
-    return NOT_A_NUMBER;
+    return (x < INFTY && x == y) ? NEGATIVE_INFTY : NOT_A_NUMBER;
   return x + log1m_exp(y - x);
 }
 

stan/math/prim/scal/fun/log_sum_exp.hpp

@@ -3,8 +3,8 @@
 
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/scal/fun/log1p_exp.hpp>
+#include <stan/math/prim/scal/fun/constants.hpp>
 #include <boost/math/tools/promotion.hpp>
-#include <limits>
 
 namespace stan {
 namespace math {
@@ -46,7 +46,10 @@ namespace math {
  */
 template <typename T1, typename T2>
 inline return_type_t<T1, T2> log_sum_exp(const T2& a, const T1& b) {
-  using std::exp;
+  if (a == NEGATIVE_INFTY)
+    return b;
+  if (a == INFTY && b == INFTY)
+    return INFTY;
   if (a > b)
     return a + log1p_exp(b - a);
   return b + log1p_exp(a - b);

stan/math/rev/mat/fun/Eigen_NumTraits.hpp

@@ -79,33 +79,32 @@ struct NumTraits<stan::math::var> : GenericNumTraits<stan::math::var> {
   static int digits10() { return std::numeric_limits<double>::digits10; }
 };
 
-namespace internal {
 /**
- * Partial specialization of Eigen's remove_all struct to stop
- * Eigen removing pointer from vari* variables
+ * Scalar product traits specialization for Eigen for reverse-mode
+ * autodiff variables.
  */
-template <>
-struct remove_all<stan::math::vari*> {
-  typedef stan::math::vari* type;
+template <typename BinaryOp>
+struct ScalarBinaryOpTraits<stan::math::var, double, BinaryOp> {
+  typedef stan::math::var ReturnType;
 };
 
-#if EIGEN_VERSION_AT_LEAST(3, 3, 0)
 /**
  * Scalar product traits specialization for Eigen for reverse-mode
  * autodiff variables.
  */
-template <>
-struct scalar_product_traits<stan::math::var, double> {
+template <typename BinaryOp>
+struct ScalarBinaryOpTraits<double, stan::math::var, BinaryOp> {
   typedef stan::math::var ReturnType;
 };
 
+namespace internal {
 /**
- * Scalar product traits specialization for Eigen for reverse-mode
- * autodiff variables.
+ * Partial specialization of Eigen's remove_all struct to stop
+ * Eigen removing pointer from vari* variables
  */
 template <>
-struct scalar_product_traits<double, stan::math::var> {
-  typedef stan::math::var ReturnType;
+struct remove_all<stan::math::vari*> {
+  typedef stan::math::vari* type;
 };
 
 /**
@@ -248,129 +247,6 @@ struct general_matrix_matrix_product<Index, stan::math::var, LhsStorageOrder,
         alpha, blocking, info);
   }
 };
-#else
-/**
- * Implemented this for printing to stream.
- */
-template <>
-struct significant_decimals_default_impl<stan::math::var, false> {
-  static inline int run() {
-    using std::ceil;
-    using std::log;
-    return cast<double, int>(
-        ceil(-log(std::numeric_limits<double>::epsilon()) / log(10.0)));
-  }
-};
-
-/**
- * Scalar product traits override for Eigen for automatic
- * gradient variables.
- */
-template <>
-struct scalar_product_traits<stan::math::var, double> {
-  typedef stan::math::var ReturnType;
-};
-
-/**
- * Scalar product traits override for Eigen for automatic
- * gradient variables.
- */
-template <>
-struct scalar_product_traits<double, stan::math::var> {
-  typedef stan::math::var ReturnType;
-};
-
-/**
- * Override matrix-vector and matrix-matrix products to use more efficient
- * implementation.
- */
-template <typename Index, bool ConjugateLhs, bool ConjugateRhs>
-struct general_matrix_vector_product<Index, stan::math::var, ColMajor,
-                                     ConjugateLhs, stan::math::var,
-                                     ConjugateRhs> {
-  typedef stan::math::var LhsScalar;
-  typedef stan::math::var RhsScalar;
-  typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType
-      ResScalar;
-  enum { LhsStorageOrder = ColMajor };
-
-  EIGEN_DONT_INLINE static void run(Index rows, Index cols,
-                                    const LhsScalar* lhs, Index lhsStride,
-                                    const RhsScalar* rhs, Index rhsIncr,
-                                    ResScalar* res, Index resIncr,
-                                    const ResScalar& alpha) {
-    for (Index i = 0; i < rows; i++) {
-      res[i * resIncr] += stan::math::var(new stan::math::gevv_vvv_vari(
-          &alpha,
-          (static_cast<int>(LhsStorageOrder) == static_cast<int>(ColMajor))
-              ? (&lhs[i])
-              : (&lhs[i * lhsStride]),
-          (static_cast<int>(LhsStorageOrder) == static_cast<int>(ColMajor))
-              ? (lhsStride)
-              : (1),
-          rhs, rhsIncr, cols));
-    }
-  }
-};
-template <typename Index, bool ConjugateLhs, bool ConjugateRhs>
-struct general_matrix_vector_product<Index, stan::math::var, RowMajor,
-                                     ConjugateLhs, stan::math::var,
-                                     ConjugateRhs> {
-  typedef stan::math::var LhsScalar;
-  typedef stan::math::var RhsScalar;
-  typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType
-      ResScalar;
-  enum { LhsStorageOrder = RowMajor };
-
-  EIGEN_DONT_INLINE static void run(Index rows, Index cols,
-                                    const LhsScalar* lhs, Index lhsStride,
-                                    const RhsScalar* rhs, Index rhsIncr,
-                                    ResScalar* res, Index resIncr,
-                                    const RhsScalar& alpha) {
-    for (Index i = 0; i < rows; i++) {
-      res[i * resIncr] += stan::math::var(new stan::math::gevv_vvv_vari(
-          &alpha,
-          (static_cast<int>(LhsStorageOrder) == static_cast<int>(ColMajor))
-              ? (&lhs[i])
-              : (&lhs[i * lhsStride]),
-          (static_cast<int>(LhsStorageOrder) == static_cast<int>(ColMajor))
-              ? (lhsStride)
-              : (1),
-          rhs, rhsIncr, cols));
-    }
-  }
-};
-template <typename Index, int LhsStorageOrder, bool ConjugateLhs,
-          int RhsStorageOrder, bool ConjugateRhs>
-struct general_matrix_matrix_product<Index, stan::math::var, LhsStorageOrder,
-                                     ConjugateLhs, stan::math::var,
-                                     RhsStorageOrder, ConjugateRhs, ColMajor> {
-  typedef stan::math::var LhsScalar;
-  typedef stan::math::var RhsScalar;
-  typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType
-      ResScalar;
-  static void run(Index rows, Index cols, Index depth, const LhsScalar* lhs,
-                  Index lhsStride, const RhsScalar* rhs, Index rhsStride,
-                  ResScalar* res, Index resStride, const ResScalar& alpha,
-                  level3_blocking<LhsScalar, RhsScalar>& /* blocking */,
-                  GemmParallelInfo<Index>* /* info = 0 */) {
-    for (Index i = 0; i < cols; i++) {
-      general_matrix_vector_product<
-          Index, LhsScalar, LhsStorageOrder, ConjugateLhs, RhsScalar,
-          ConjugateRhs>::run(rows, depth, lhs, lhsStride,
-                             &rhs[(static_cast<int>(RhsStorageOrder)
-                                   == static_cast<int>(ColMajor))
-                                      ? (i * rhsStride)
-                                      : (i)],
-                             (static_cast<int>(RhsStorageOrder)
-                              == static_cast<int>(ColMajor))
-                                 ? (1)
-                                 : (rhsStride),
-                             &res[i * resStride], 1, alpha);
-    }
-  }
-};
-#endif
 }  // namespace internal
 }  // namespace Eigen
 #endif

stan/math/rev/mat/fun/log_sum_exp.hpp

@@ -6,6 +6,7 @@
 #include <stan/math/rev/mat/fun/typedefs.hpp>
 #include <stan/math/prim/scal/fun/log_sum_exp.hpp>
 #include <stan/math/prim/mat/fun/Eigen.hpp>
+#include <cmath>
 #include <limits>
 
 namespace stan {
@@ -19,6 +20,8 @@ namespace internal {
 template <int R, int C>
 inline double log_sum_exp_as_double(const Eigen::Matrix<var, R, C>& x) {
   const double max = x.val().maxCoeff();
+  if (!std::isfinite(max))
+    return max;
   return max + std::log((x.val().array() - max).exp().sum());
 }
 

stan/math/rev/scal/fun/log_diff_exp.hpp

@@ -4,6 +4,7 @@
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/core.hpp>
 #include <stan/math/rev/scal/fun/calculate_chain.hpp>
+#include <stan/math/prim/scal/fun/constants.hpp>
 #include <stan/math/prim/scal/fun/expm1.hpp>
 #include <stan/math/prim/scal/fun/log_diff_exp.hpp>
 
@@ -24,13 +25,23 @@ class log_diff_exp_vd_vari : public op_vd_vari {
  public:
   log_diff_exp_vd_vari(vari* avi, double b)
       : op_vd_vari(log_diff_exp(avi->val_, b), avi, b) {}
-  void chain() { avi_->adj_ += adj_ * calculate_chain(avi_->val_, val_); }
+  void chain() {
+    if (val_ == NEGATIVE_INFTY)
+      avi_->adj_ += (bd_ == NEGATIVE_INFTY) ? adj_ : adj_ * INFTY;
+    else
+      avi_->adj_ += adj_ * calculate_chain(avi_->val_, val_);
+  }
 };
 class log_diff_exp_dv_vari : public op_dv_vari {
  public:
   log_diff_exp_dv_vari(double a, vari* bvi)
       : op_dv_vari(log_diff_exp(a, bvi->val_), a, bvi) {}
-  void chain() { bvi_->adj_ -= adj_ / expm1(ad_ - bvi_->val_); }
+  void chain() {
+    if (val_ == NEGATIVE_INFTY)
+      bvi_->adj_ -= adj_ * INFTY;
+    else
+      bvi_->adj_ -= adj_ / expm1(ad_ - bvi_->val_);
+  }
 };
 }  // namespace internal
 
@@ -39,7 +50,7 @@ class log_diff_exp_dv_vari : public op_dv_vari {
  *
  * @param[in] a First argument.
  * @param[in] b Second argument.
- * @return Log difference of the expnoentiated arguments.
+ * @return Log difference of the exponentiated arguments.
  */
 inline var log_diff_exp(const var& a, const var& b) {
   return var(new internal::log_diff_exp_vv_vari(a.vi_, b.vi_));
@@ -50,7 +61,7 @@ inline var log_diff_exp(const var& a, const var& b) {
  *
  * @param[in] a First argument.
  * @param[in] b Second argument.
- * @return Log difference of the expnoentiated arguments.
+ * @return Log difference of the exponentiated arguments.
  */
 inline var log_diff_exp(const var& a, double b) {
   return var(new internal::log_diff_exp_vd_vari(a.vi_, b));
@@ -61,7 +72,7 @@ inline var log_diff_exp(const var& a, double b) {
  *
  * @param[in] a First argument.
  * @param[in] b Second argument.
- * @return Log difference of the expnoentiated arguments.
+ * @return Log difference of the exponentiated arguments.
  */
 inline var log_diff_exp(double a, const var& b) {
   return var(new internal::log_diff_exp_dv_vari(a, b.vi_));

stan/math/rev/scal/fun/log_sum_exp.hpp

@@ -4,6 +4,7 @@
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/core.hpp>
 #include <stan/math/rev/scal/fun/calculate_chain.hpp>
+#include <stan/math/prim/scal/fun/constants.hpp>
 #include <stan/math/prim/scal/fun/log_sum_exp.hpp>
 
 namespace stan {
@@ -24,13 +25,23 @@ class log_sum_exp_vd_vari : public op_vd_vari {
  public:
   log_sum_exp_vd_vari(vari* avi, double b)
       : op_vd_vari(log_sum_exp(avi->val_, b), avi, b) {}
-  void chain() { avi_->adj_ += adj_ * calculate_chain(avi_->val_, val_); }
+  void chain() {
+    if (val_ == NEGATIVE_INFTY)
+      avi_->adj_ += adj_;
+    else
+      avi_->adj_ += adj_ * calculate_chain(avi_->val_, val_);
+  }
 };
 class log_sum_exp_dv_vari : public op_dv_vari {
  public:
   log_sum_exp_dv_vari(double a, vari* bvi)
       : op_dv_vari(log_sum_exp(a, bvi->val_), a, bvi) {}
-  void chain() { bvi_->adj_ += adj_ * calculate_chain(bvi_->val_, val_); }
+  void chain() {
+    if (val_ == NEGATIVE_INFTY)
+      bvi_->adj_ += adj_;
+    else
+      bvi_->adj_ += adj_ * calculate_chain(bvi_->val_, val_);
+  }
 };
 
 }  // namespace internal

test/unit/math/fwd/scal/fun/log_diff_exp_test.cpp

 #include <stan/math/fwd/scal.hpp>
 #include <gtest/gtest.h>
 #include <test/unit/math/fwd/scal/fun/nan_util.hpp>
+#include <limits>
 
 TEST(AgradFwdLogDiffExp, Fvar) {
   using stan::math::fvar;
@@ -25,6 +26,26 @@ TEST(AgradFwdLogDiffExp, Fvar) {
   fvar<double> c = log_diff_exp(z, y);
   EXPECT_FLOAT_EQ(log_diff_exp(1.1, 0.5), c.val_);
   EXPECT_FLOAT_EQ(2 / (1 - exp(1.1 - 0.5)), c.d_);
+
+  double ninf = -std::numeric_limits<double>::infinity();
+  fvar<double> d = log_diff_exp(x, ninf);
+  EXPECT_FLOAT_EQ(x.val_, d.val_);
+  EXPECT_FLOAT_EQ(1.0, d.d_);
+
+  fvar<double> e = log_diff_exp(x, 1.2);
+  EXPECT_FLOAT_EQ(-std::numeric_limits<double>::infinity(), e.val_);
+  EXPECT_FLOAT_EQ(+std::numeric_limits<double>::infinity(), e.d_);
+
+  fvar<double> w(-std::numeric_limits<double>::infinity());
+  w.d_ = 1.0;
+
+  fvar<double> g = log_diff_exp(w, ninf);
+  EXPECT_FLOAT_EQ(w.val_, g.val_);
+  EXPECT_FLOAT_EQ(1.0, g.d_);
+
+  fvar<double> f = log_diff_exp(0.5, y);
+  EXPECT_FLOAT_EQ(-std::numeric_limits<double>::infinity(), f.val_);
+  EXPECT_FLOAT_EQ(-std::numeric_limits<double>::infinity(), f.d_);
 }
 
 TEST(AgradFwdLogDiffExp, AgradFvar_exception) {

test/unit/math/fwd/scal/fun/log_sum_exp_test.cpp

 #include <stan/math/fwd/scal.hpp>
 #include <gtest/gtest.h>
 #include <test/unit/math/fwd/scal/fun/nan_util.hpp>
+#include <limits>
 
 TEST(AgradFwdLogSumExp, Fvar) {
   using stan::math::fvar;
@@ -23,6 +24,17 @@ TEST(AgradFwdLogSumExp, Fvar) {
   fvar<double> c = log_sum_exp(z, x);
   EXPECT_FLOAT_EQ(log_sum_exp(1.4, 0.5), c.val_);
   EXPECT_FLOAT_EQ(1.0 * exp(0.5) / (exp(0.5) + exp(1.4)), c.d_);
+
+  double ninf = -std::numeric_limits<double>::infinity();
+  fvar<double> w(-std::numeric_limits<double>::infinity(), 1.5);
+
+  fvar<double> d = log_sum_exp(w, ninf);
+  EXPECT_FLOAT_EQ(-std::numeric_limits<double>::infinity(), d.val_);
+  EXPECT_FLOAT_EQ(1.5, d.d_);
+
+  fvar<double> e = log_sum_exp(ninf, w);
+  EXPECT_FLOAT_EQ(-std::numeric_limits<double>::infinity(), e.val_);
+  EXPECT_FLOAT_EQ(1.5, e.d_);
 }
 
 TEST(AgradFwdLogSumExp, FvarFvarDouble) {

test/unit/math/opencl/matrix_cl_view.cpp → test/unit/math/opencl/matrix_cl_view_test.cpp

@@ -93,15 +93,13 @@ TEST(matrix_cl_view, invert) {
   EXPECT_EQ(matrix_cl_view::Diagonal, invert(matrix_cl_view::Entire));
 }
 
-TEST(matrix_cl_view, from_eigen_triangular_type) {
-  using stan::math::from_eigen_triangular_type;
+TEST(matrix_cl_view, from_eigen_uplo_type) {
+  using stan::math::from_eigen_uplo_type;
   using stan::math::matrix_cl_view;
-  EXPECT_EQ(matrix_cl_view::Lower, from_eigen_triangular_type(Eigen::Lower));
-  EXPECT_EQ(matrix_cl_view::Upper, from_eigen_triangular_type(Eigen::Upper));
-  EXPECT_EQ(matrix_cl_view::Entire,
-            from_eigen_triangular_type(Eigen::SelfAdjoint));
-  EXPECT_EQ(matrix_cl_view::Entire,
-            from_eigen_triangular_type(Eigen::UnitDiag));
+  EXPECT_EQ(matrix_cl_view::Lower, from_eigen_uplo_type(Eigen::Lower));
+  EXPECT_EQ(matrix_cl_view::Upper, from_eigen_uplo_type(Eigen::Upper));
+  EXPECT_EQ(matrix_cl_view::Entire, from_eigen_uplo_type(Eigen::SelfAdjoint));
+  EXPECT_EQ(matrix_cl_view::Entire, from_eigen_uplo_type(Eigen::UnitDiag));
 }
 
 #endif

test/unit/math/prim/mat/fun/log_sum_exp_test.cpp

 #include <stan/math/prim/mat.hpp>
 #include <gtest/gtest.h>
+#include <limits>
 
 template <int R, int C>
 void test_log_sum_exp(const Eigen::Matrix<double, R, C>& as) {
@@ -25,11 +26,19 @@ TEST(MathFunctions, log_sum_exp) {
   v << 1, 2, 3;
   test_log_sum_exp(v);
 
-  Matrix<double, Dynamic, 1> rv(3);
+  Matrix<double, 1, Dynamic> rv(3);
   rv << 1, 2, 3;
   test_log_sum_exp(rv);
 
   Matrix<double, Dynamic, Dynamic> m_trivial(1, 1);
   m_trivial << 2;
   EXPECT_FLOAT_EQ(2, log_sum_exp(m_trivial));
+
+  Matrix<double, Dynamic, 1> i(3);
+  i << 1, 2, -std::numeric_limits<double>::infinity();
+  test_log_sum_exp(i);
+
+  Matrix<double, Dynamic, 1> ii(1);
+  ii << -std::numeric_limits<double>::infinity();
+  test_log_sum_exp(ii);
 }

test/unit/math/prim/mat/fun/mdivide_left_tri_test.cpp

@@ -2,10 +2,7 @@
 #include <gtest/gtest.h>
 
 #ifdef STAN_OPENCL
-#include <stan/math/opencl/opencl_context.hpp>
-#include <stan/math/opencl/multiply.hpp>
-#include <stan/math/opencl/copy.hpp>
-#include <stan/math/opencl/tri_inverse.hpp>
+#include <stan/math/opencl/opencl.hpp>
 #include <boost/random/mersenne_twister.hpp>
 #endif
 

test/unit/math/prim/mat/fun/multiply_test.cpp

@@ -133,3 +133,31 @@ TEST(AgradRevMatrix, multiply_vector_int) {
   EXPECT_EQ(4.0, prod_vec[1]);
   EXPECT_EQ(6.0, prod_vec[2]);
 }
+
+#ifdef STAN_OPENCL
+#define EXPECT_MATRIX_NEAR(A, B, DELTA) \
+  for (int i = 0; i < A.size(); i++)    \
+    EXPECT_NEAR(A(i), B(i), DELTA);
+
+TEST(MathMatrix, multiply_opencl) {
+  int multiply_dim_prod_worth_transfer
+      = stan::math::opencl_context.tuning_opts()
+            .multiply_dim_prod_worth_transfer;
+  stan::math::opencl_context.tuning_opts().multiply_dim_prod_worth_transfer = 0;
+  using stan::math::multiply;
+  int size = 400;
+  stan::math::matrix_d m1 = stan::math::matrix_d::Random(size, size).eval();
+  stan::math::matrix_d m2 = stan::math::matrix_d::Random(size, size).eval();
+
+  stan::math::matrix_d m3_cl = multiply(m1, m2);
+  // to make sure we dont use OpenCL
+  stan::math::opencl_context.tuning_opts().multiply_dim_prod_worth_transfer
+      = INT_MAX;
+
+  stan::math::matrix_d m3 = multiply(m1, m2);
+
+  EXPECT_MATRIX_NEAR(m3_cl, m3, 1e-10);
+  stan::math::opencl_context.tuning_opts().multiply_dim_prod_worth_transfer
+      = multiply_dim_prod_worth_transfer;
+}
+#endif