<commit_before>#ifndef __STAN__AGRAD__MATRIX_HPP__
#define __STAN__AGRAD__MATRIX_HPP__

#include <stan/math/matrix.hpp>
#include <stan/math/matrix_error_handling.hpp>

#include <stan/agrad/special_functions.hpp>

namespace stan {
  namespace agrad {
    typedef Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic>::size_type size_type;

    void fill(var& x, const var& y) {
      x = y;
    }
    template <typename T>
    void fill(std::vector<T>& x, const var& y) {
      for (size_t i = 0; i < x.size(); ++i)
        fill(x[i],y);
    }
    template <int R, int C>
    void fill(Eigen::Matrix<var,R,C>& x, const var& y) {
      x.fill(y);
    }



    class gevv_vvv_vari : public stan::agrad::vari {
    protected:
      stan::agrad::vari* alpha_;
      stan::agrad::vari** v1_;
      stan::agrad::vari** v2_;
      double dotval_;
      size_t length_;
      inline static double eval_gevv(const stan::agrad::var* alpha,
                                     const stan::agrad::var* v1, int stride1,
                                     const stan::agrad::var* v2, int stride2,
                                     size_t length, double *dotprod) {
        double result = 0;
        for (size_t i = 0; i < length; i++)
          result += v1[i*stride1].vi_->val_ * v2[i*stride2].vi_->val_;
        *dotprod = result;
        return alpha->vi_->val_ * result;
      }
    public:
      gevv_vvv_vari(const stan::agrad::var* alpha, 
                    const stan::agrad::var* v1, int stride1, 
                    const stan::agrad::var* v2, int stride2, size_t length) : 
        vari(eval_gevv(alpha,v1,stride1,v2,stride2,length,&dotval_)), length_(length) {
        alpha_ = alpha->vi_;
        v1_ = (stan::agrad::vari**)stan::agrad::memalloc_.alloc(2*length_*sizeof(stan::agrad::vari*));
        v2_ = v1_ + length_;
        for (size_t i = 0; i < length_; i++)
          v1_[i] = v1[i*stride1].vi_;
        for (size_t i = 0; i < length_; i++)
          v2_[i] = v2[i*stride2].vi_;
      }
      void chain() {
        const double adj_alpha = adj_ * alpha_->val_;
        for (size_t i = 0; i < length_; i++) {
          v1_[i]->adj_ += adj_alpha * v2_[i]->val_;
          v2_[i]->adj_ += adj_alpha * v1_[i]->val_;
        }
        alpha_->adj_ += adj_ * dotval_;
      }
    };
  }
}

namespace Eigen {

  /**
   * Numerical traits template override for Eigen for automatic
   * gradient variables.
   */
  template <> struct NumTraits<stan::agrad::var>
  {
    /**
     * Real-valued variables.
     *
     * Required for numerical traits.
     */
    typedef stan::agrad::var Real;

    /**
     * Non-integer valued variables.
     *
     * Required for numerical traits.
     */
    typedef stan::agrad::var NonInteger;

    /**
     * Nested variables.
     *
     * Required for numerical traits.
     */
    typedef stan::agrad::var Nested;

    /**
     * Return standard library's epsilon for double-precision floating
     * point, <code>std::numeric_limits&lt;double&gt;::epsilon()</code>.
     *
     * @return Same epsilon as a <code>double</code>.
     */
    inline static Real epsilon() { 
      return std::numeric_limits<double>::epsilon(); 
    }

    /**
     * Return dummy precision
     */
    inline static Real dummy_precision() {
      return 1e-12; // copied from NumTraits.h values for double
    }

    /**
     * Return standard library's highest for double-precision floating
     * point, <code>std::numeric_limits&lt;double&gt;::max()</code>.
     *
     * @return Same highest value as a <code>double</code>.
     */
    inline static Real highest() {
      return std::numeric_limits<double>::max();
    }

    /**
     * Return standard library's lowest for double-precision floating
     * point, <code>&#45;std::numeric_limits&lt;double&gt;::max()</code>.
     *
     * @return Same lowest value as a <code>double</code>.
     */    
    inline static Real lowest() {
      return -std::numeric_limits<double>::max();
    }

    /**
     * Properties for automatic differentiation variables
     * read by Eigen matrix library.
     */
    enum {
      IsInteger = 0,
      IsSigned = 1,
      IsComplex = 0,
      RequireInitialization = 0,
      ReadCost = 1,
      AddCost = 1,
      MulCost = 1,
      HasFloatingPoint = 1
    };
  };

  namespace internal {
    /**
     * Implemented this for printing to stream.
     */
    template<>
    struct significant_decimals_default_impl<stan::agrad::var,false>
    {
      static inline int run()
      {
        using std::ceil;
        return cast<double,int>(ceil(-log(NumTraits<stan::agrad::var>::epsilon().val())
                                     /log(10.0)));
      }
    };

    /**
     * Scalar product traits override for Eigen for automatic
     * gradient variables.
     */
    template <>  
    struct scalar_product_traits<stan::agrad::var,double> {
      typedef stan::agrad::var ReturnType;
    };

    /**
     * Scalar product traits override for Eigen for automatic
     * gradient variables.
     */
    template <>  
    struct scalar_product_traits<double,stan::agrad::var> {
      typedef stan::agrad::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::agrad::var,ColMajor,ConjugateLhs,stan::agrad::var,ConjugateRhs>
    {
      typedef stan::agrad::var LhsScalar;
      typedef stan::agrad::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::agrad::var(new stan::agrad::gevv_vvv_vari(&alpha,((int)LhsStorageOrder == (int)ColMajor)?(&lhs[i]):(&lhs[i*lhsStride]),((int)LhsStorageOrder == (int)ColMajor)?(lhsStride):(1),rhs,rhsIncr,cols));
        }
      }
    };
    template<typename Index, bool ConjugateLhs, bool ConjugateRhs>
    struct general_matrix_vector_product<Index,stan::agrad::var,RowMajor,ConjugateLhs,stan::agrad::var,ConjugateRhs>
    {
      typedef stan::agrad::var LhsScalar;
      typedef stan::agrad::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::agrad::var(new stan::agrad::gevv_vvv_vari(&alpha,((int)LhsStorageOrder == (int)ColMajor)?(&lhs[i]):(&lhs[i*lhsStride]),((int)LhsStorageOrder == (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::agrad::var,LhsStorageOrder,ConjugateLhs,stan::agrad::var,RhsStorageOrder,ConjugateRhs,ColMajor>
    {
      typedef stan::agrad::var LhsScalar;
      typedef stan::agrad::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[((int)RhsStorageOrder == (int)ColMajor)?(i*rhsStride):(i)],((int)RhsStorageOrder == (int)ColMajor)?(1):(rhsStride),
              &res[i*resStride],1,alpha);
        }
      }
    };
  }
}

namespace stan {

  namespace agrad {

    /**
     * The type of a matrix holding <code>stan::agrad::var</code>
     * values.
     */
    typedef 
    Eigen::Matrix<var,Eigen::Dynamic,Eigen::Dynamic>
    matrix_v;

    /**
     * The type of a (column) vector holding <code>stan::agrad::var</code>
     * values.
     */
    typedef 
    Eigen::Matrix<var,Eigen::Dynamic,1>
    vector_v;

    /**
     * The type of a row vector holding <code>stan::agrad::var</code>
     * values.
     */
    typedef 
    Eigen::Matrix<var,1,Eigen::Dynamic>
    row_vector_v;

    /**
     * Initialize variable to value.  (Function may look pointless, but
     * its needed to bottom out recursion.)
     */
    inline void initialize_variable(var& variable, const var& value) {
      variable = value;
    }

    /**
     * Initialize every cell in the matrix to the specified value.
     * 
     */
    template <int R, int C>
    inline void initialize_variable(Eigen::Matrix<var,R,C>& matrix, const var& value) {
      for (int i = 0; i < matrix.size(); ++i)
        matrix(i) = value;
    }

    /**
     * Initialize the variables in the standard vector recursively.
     */
    template <typename T>
    inline void initialize_variable(std::vector<T>& variables, const var& value) {
      for (size_t i = 0; i < variables.size(); ++i)
        initialize_variable(variables[i],value);
    }


    /**
     * Converts argument to an automatic differentiation variable.
     *
     * Returns a stan::agrad::var variable with the input value.
     *
     * @param[in] x A scalar value
     * @return An automatic differentiation variable with the input value.
     */
    inline var to_var(const double& x) {
      return var(x);
    }
    /**
     * Converts argument to an automatic differentiation variable.
     *
     * Returns a stan::agrad::var variable with the input value.
     *
     * @param[in] x An automatic differentiation variable.
     * @return An automatic differentiation variable with the input value.
     */    
    inline var to_var(const var& x) {
      return x;
    }
    /**
     * Converts argument to an automatic differentiation variable.
     *
     * Returns a stan::agrad::var variable with the input value.
     *
     * @param[in] m A Matrix with scalars
     * @return A Matrix with automatic differentiation variables
     */
    inline matrix_v to_var(const stan::math::matrix_d& m) {
      matrix_v m_v(m.rows(), m.cols());
      for (int j = 0; j < m.cols(); ++j)
        for (int i = 0; i < m.rows(); ++i)
          m_v(i,j) = m(i,j);
      return m_v;
    }
    /**
     * Converts argument to an automatic differentiation variable.
     *
     * Returns a stan::agrad::var variable with the input value.
     * 
     * @param[in] m A Matrix with automatic differentiation variables.
     * @return A Matrix with automatic differentiation variables.
     */
    inline matrix_v to_var(const matrix_v& m) {
      return m;
    }
    /**
     * Converts argument to an automatic differentiation variable.
     *
     * Returns a stan::agrad::var variable with the input value.     
     *
     * @param[in] v A Vector of scalars
     * @return A Vector of automatic differentiation variables with
     *   values of v
     */
    inline vector_v to_var(const stan::math::vector_d& v) {
      vector_v v_v(v.size());
      for (int i = 0; i < v.size(); ++i)
        v_v[i] = v[i];
      return v_v;
    }
    /**
     * Converts argument to an automatic differentiation variable.
     *
     * Returns a stan::agrad::var variable with the input value.     
     *
     * @param[in] v A Vector of automatic differentiation variables
     * @return A Vector of automatic differentiation variables with
     *   values of v
     */
    inline vector_v to_var(const vector_v& v) {
      return v;
    }
    /**
     * Converts argument to an automatic differentiation variable.
     *
     * Returns a stan::agrad::var variable with the input value.     
     *
     * @param[in] rv A row vector of scalars
     * @return A row vector of automatic differentation variables with 
     *   values of rv.
     */
    inline row_vector_v to_var(const stan::math::row_vector_d& rv) {
      row_vector_v rv_v(rv.size());
      for (int i = 0; i < rv.size(); ++i)
        rv_v[i] = rv[i];
      return rv_v;
    }
    /**
     * Converts argument to an automatic differentiation variable.
     *
     * Returns a stan::agrad::var variable with the input value.     
     *
     * @param[in] rv A row vector with automatic differentiation variables
     * @return A row vector with automatic differentiation variables
     *    with values of rv.
     */
    inline row_vector_v to_var(const row_vector_v& rv) {
      return rv;
    }



    // scalar returns

    // /**
    //  * Determinant of the matrix.
    //  *
    //  * Returns the determinant of the specified
    //  * square matrix.
    //  *
    //  * @param m Specified matrix.
    //  * @return Determinant of the matrix.
    //  * @throw std::domain_error if m is not a square matrix
    //  */
    // var determinant(const Eigen::Matrix<var, Eigen::Dynamic, Eigen::Dynamic>& m);
    
    namespace {
      class dot_self_vari : public vari {
      protected:
        vari** v_;
        size_t size_;
      public:
        dot_self_vari(vari** v, size_t size) 
          : vari(var_dot_self(v,size)), 
            v_(v),
            size_(size) {
        }
        template<typename Derived>
        dot_self_vari(const Eigen::DenseBase<Derived> &v) : 
          vari(var_dot_self(v)), size_(v.size()) {
          v_ = (vari**)memalloc_.alloc(size_*sizeof(vari*));
          for (size_t i = 0; i < size_; i++)
            v_[i] = v[i].vi_;
        }
        template <int R, int C>
        dot_self_vari(const Eigen::Matrix<var,R,C>& v) :
          vari(var_dot_self(v)), size_(v.size()) {
          v_ = (vari**) memalloc_.alloc(size_ * sizeof(vari*));
          for (size_t i = 0; i < size_; ++i)
            v_[i] = v(i).vi_;
        }
        inline static double square(double x) { return x * x; }
        inline static double var_dot_self(vari** v, size_t size) {
          double sum = 0.0;
          for (size_t i = 0; i < size; ++i)
            sum += square(v[i]->val_);
          return sum;
        }
        template<typename Derived>
        double var_dot_self(const Eigen::DenseBase<Derived> &v) {
          double sum = 0.0;
          for (int i = 0; i < v.size(); ++i)
            sum += square(v(i).vi_->val_);
          return sum;
        }
        template <int R, int C>
        inline static double var_dot_self(const Eigen::Matrix<var,R,C> &v) {
          double sum = 0.0;
          for (int i = 0; i < v.size(); ++i)
            sum += square(v(i).vi_->val_);
          return sum;
        }
        virtual void chain() {
          for (size_t i = 0; i < size_; ++i) 
            v_[i]->adj_ += adj_ * 2.0 * v_[i]->val_;
        }
      };

      class sum_v_vari : public vari{
      protected:
        vari** v_;
        size_t length_;
        inline static double var_sum(const var *v, size_t length) {
          double result = 0;
          for (size_t i = 0; i < length; i++)
            result += v[i].vi_->val_;
          return result;
        } 
        template<typename Derived>
        inline static double var_sum(const Eigen::DenseBase<Derived> &v) {
          double result = 0;
          for (int i = 0; i < v.size(); i++)
            result += v(i).vi_->val_;
          return result;
        } 
      public:
        template<typename Derived>
        sum_v_vari(const Eigen::DenseBase<Derived> &v) :
          vari(var_sum(v)), length_(v.size()) {
          v_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
          for (size_t i = 0; i < length_; i++)
            v_[i] = v(i).vi_;
        }
        template<int R1,int C1>
        sum_v_vari(const Eigen::Matrix<var,R1,C1> &v1) :
          vari(var_sum(v1)), length_(v1.size()) {
          v_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
          for (size_t i = 0; i < length_; i++)
            v_[i] = v1(i).vi_;
        }
        sum_v_vari(const var *v, size_t len) :
          vari(var_sum(v,len)), length_(len) {
          v_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
          for (size_t i = 0; i < length_; i++)
            v_[i] = v[i].vi_;
        }
        virtual void chain() {
          for (size_t i = 0; i < length_; i++) {
            v_[i]->adj_ += adj_;
          }
        }
      };
      class dot_product_vv_vari : public vari {
      protected:
        vari** v1_;
        vari** v2_;
        size_t length_;
        inline static double var_dot(const var* v1, const var* v2,
                                     size_t length) {
          double result = 0;
          for (size_t i = 0; i < length; i++)
            result += v1[i].vi_->val_ * v2[i].vi_->val_;
          return result;
        }
        template<typename Derived1,typename Derived2>
        inline static double var_dot(const Eigen::DenseBase<Derived1> &v1,
                                     const Eigen::DenseBase<Derived2> &v2) {
          double result = 0;
          for (int i = 0; i < v1.size(); i++)
            result += v1[i].vi_->val_ * v2[i].vi_->val_;
          return result;
        }
        inline static double var_dot(vari** v1, vari** v2, size_t length) {
          double result = 0;
          for (size_t i = 0; i < length; ++i)
            result += v1[i]->val_ * v2[i]->val_;
          return result;
        }
      public:
        dot_product_vv_vari(vari** v1, vari** v2, size_t length)
          : vari(var_dot(v1,v2,length)),
            v1_(v1), 
            v2_(v2), 
            length_(length) {

        }
        dot_product_vv_vari(const var* v1, const var* v2, size_t length,
                            dot_product_vv_vari* shared_v1 = NULL,
                            dot_product_vv_vari* shared_v2 = NULL) : 
          vari(var_dot(v1, v2, length)), length_(length) {
          if (shared_v1 == NULL) {
            v1_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
            for (size_t i = 0; i < length_; i++)
              v1_[i] = v1[i].vi_;
          }
          else {
            v1_ = shared_v1->v1_;
          }
          if (shared_v2 == NULL) {
            v2_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
            for (size_t i = 0; i < length_; i++)
              v2_[i] = v2[i].vi_;
          }
          else {
            v2_ = shared_v2->v2_;
          }
        }
        template<typename Derived1,typename Derived2>
        dot_product_vv_vari(const Eigen::DenseBase<Derived1> &v1,
                            const Eigen::DenseBase<Derived2> &v2,
                            dot_product_vv_vari* shared_v1 = NULL,
                            dot_product_vv_vari* shared_v2 = NULL) : 
          vari(var_dot(v1, v2)), length_(v1.size()) {
          if (shared_v1 == NULL) {
            v1_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
            for (size_t i = 0; i < length_; i++)
              v1_[i] = v1[i].vi_;
          }
          else {
            v1_ = shared_v1->v1_;
          }
          if (shared_v2 == NULL) {
            v2_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
            for (size_t i = 0; i < length_; i++)
              v2_[i] = v2[i].vi_;
          }
          else {
            v2_ = shared_v2->v2_;
          }
        }
        template<int R1,int C1,int R2,int C2>
        dot_product_vv_vari(const Eigen::Matrix<var,R1,C1> &v1,
                            const Eigen::Matrix<var,R2,C2> &v2,
                            dot_product_vv_vari* shared_v1 = NULL,
                            dot_product_vv_vari* shared_v2 = NULL) : 
          vari(var_dot(v1, v2)), length_(v1.size()) {
          if (shared_v1 == NULL) {
            v1_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
            for (size_t i = 0; i < length_; i++)
              v1_[i] = v1[i].vi_;
          }
          else {
            v1_ = shared_v1->v1_;
          }
          if (shared_v2 == NULL) {
            v2_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
            for (size_t i = 0; i < length_; i++)
              v2_[i] = v2[i].vi_;
          }
          else {
            v2_ = shared_v2->v2_;
          }
        }
        virtual void chain() {
          for (size_t i = 0; i < length_; i++) {
            v1_[i]->adj_ += adj_ * v2_[i]->val_;
            v2_[i]->adj_ += adj_ * v1_[i]->val_;
          }
        }
      };

      class dot_product_vd_vari : public vari {
      protected:
        vari** v1_;
        double* v2_;
        size_t length_;
        inline static double var_dot(const var* v1, const double* v2,
                                     size_t length) {
          double result = 0;
          for (size_t i = 0; i < length; i++)
            result += v1[i].vi_->val_ * v2[i];
          return result;
        }
        template<typename Derived1,typename Derived2>
        inline static double var_dot(const Eigen::DenseBase<Derived1> &v1,
                                     const Eigen::DenseBase<Derived2> &v2) {
          double result = 0;
          for (int i = 0; i < v1.size(); i++)
            result += v1[i].vi_->val_ * v2[i];
          return result;
        }
      public:
        dot_product_vd_vari(const var* v1, const double* v2, size_t length,
                            dot_product_vd_vari *shared_v1 = NULL,
                            dot_product_vd_vari *shared_v2 = NULL) : 
          vari(var_dot(v1, v2, length)), length_(length) {
          if (shared_v1 == NULL) {
            v1_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
            for (size_t i = 0; i < length_; i++)
              v1_[i] = v1[i].vi_;
          } else {
            v1_ = shared_v1->v1_;
          }
          if (shared_v2 == NULL) {
            v2_ = (double*)memalloc_.alloc(length_*sizeof(double));
            for (size_t i = 0; i < length_; i++)
              v2_[i] = v2[i];
          } else {
            v2_ = shared_v2->v2_;
          }
        }
        template<typename Derived1,typename Derived2>
        dot_product_vd_vari(const Eigen::DenseBase<Derived1> &v1,
                            const Eigen::DenseBase<Derived2> &v2,
                            dot_product_vd_vari *shared_v1 = NULL,
                            dot_product_vd_vari *shared_v2 = NULL) : 
          vari(var_dot(v1, v2)), length_(v1.size()) {
          if (shared_v1 == NULL) {
            v1_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
            for (size_t i = 0; i < length_; i++)
              v1_[i] = v1[i].vi_;
          } else {
            v1_ = shared_v1->v1_;
          }
          if (shared_v2 == NULL) {
            v2_ = (double*)memalloc_.alloc(length_*sizeof(double));
            for (size_t i = 0; i < length_; i++)
              v2_[i] = v2[i];
          } else {
            v2_ = shared_v2->v2_;
          }
        }
        template<int R1,int C1,int R2,int C2>
        dot_product_vd_vari(const Eigen::Matrix<var,R1,C1> &v1,
                            const Eigen::Matrix<double,R2,C2> &v2,
                            dot_product_vd_vari *shared_v1 = NULL,
                            dot_product_vd_vari *shared_v2 = NULL) : 
          vari(var_dot(v1, v2)), length_(v1.size()) {
          if (shared_v1 == NULL) {
            v1_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
            for (size_t i = 0; i < length_; i++)
              v1_[i] = v1[i].vi_;
          } else {
            v1_ = shared_v1->v1_;
          }
          if (shared_v2 == NULL) {
            v2_ = (double*)memalloc_.alloc(length_*sizeof(double));
            for (size_t i = 0; i < length_; i++)
              v2_[i] = v2[i];
          } else {
            v2_ = shared_v2->v2_;
          }
        }
        virtual void chain() {
          for (size_t i = 0; i < length_; i++) {
            v1_[i]->adj_ += adj_ * v2_[i];
          }
        }
      };

    }

    /**
     * Returns the dot product of a vector with itself.
     *
     * @param[in] v Vector.
     * @return Dot product of the vector with itself.
     * @tparam R number of rows or <code>Eigen::Dynamic</code> for
     * dynamic; one of R or C must be 1
     * @tparam C number of rows or <code>Eigen::Dyanmic</code> for
     * dynamic; one of R or C must be 1
     */
    template<int R, int C>
    inline var dot_self(const Eigen::Matrix<var, R, C>& v) {
      stan::math::validate_vector(v,"dot_self");
      return var(new dot_self_vari(v));
    }
    
    /**
     * Returns the dot product of each column of a matrix with itself.
     * @param x Matrix.
     * @tparam T scalar type
     */
    template<int R,int C>
    inline Eigen::Matrix<var,1,C> 
    columns_dot_self(const Eigen::Matrix<var,R,C>& x) {
      Eigen::Matrix<var,1,C> ret(1,x.cols());
      for (size_type i = 0; i < x.cols(); i++) {
        ret(i) = var(new dot_self_vari(x.col(i)));
      }
      return ret;
    }
    
    template<int R1,int C1,int R2, int C2>
    inline Eigen::Matrix<var, 1, C1>
    columns_dot_product(const Eigen::Matrix<var, R1, C1>& v1, 
                        const Eigen::Matrix<var, R2, C2>& v2) {
      stan::math::validate_matching_sizes(v1,v2,"columns_dot_product");
      Eigen::Matrix<var, 1, C1> ret(1,v1.cols());
      for (size_type j = 0; j < v1.cols(); ++j) {
        ret(j) = var(new dot_product_vv_vari(v1.col(j),v2.col(j)));
      }
      return ret;
    }
    
    template<int R1,int C1,int R2, int C2>
    inline Eigen::Matrix<var, 1, C1>
    columns_dot_product(const Eigen::Matrix<var, R1, C1>& v1, 
                        const Eigen::Matrix<double, R2, C2>& v2) {
      stan::math::validate_matching_sizes(v1,v2,"columns_dot_product");
      Eigen::Matrix<var, 1, C1> ret(1,v1.cols());
      size_t j;
      for (j = 0; j < v1.cols(); ++j) {
        ret(j) = var(new dot_product_vd_vari(v1.col(j),v2.col(j)));
      }
      return ret;
    }

    template<int R1,int C1,int R2, int C2>
    inline Eigen::Matrix<var, 1, C1>
    columns_dot_product(const Eigen::Matrix<double, R1, C1>& v1, 
                        const Eigen::Matrix<var, R2, C2>& v2) {
      stan::math::validate_matching_sizes(v1,v2,"columns_dot_product");
      Eigen::Matrix<var, 1, C1> ret(1,v1.cols());
      for (size_type j = 0; j < v1.cols(); ++j) {
        ret(j) = var(new dot_product_vd_vari(v2.col(j),v1.col(j)));
      }
      return ret;
    }
    
    /**
     * Returns the dot product.
     *
     * @param[in] v1 First column vector.
     * @param[in] v2 Second column vector.
     * @return Dot product of the vectors.
     * @throw std::domain_error if length of v1 is not equal to length of v2.
     */
    template<int R1,int C1,int R2, int C2>
    inline var dot_product(const Eigen::Matrix<var, R1, C1>& v1, 
                           const Eigen::Matrix<var, R2, C2>& v2) {
      stan::math::validate_vector(v1,"dot_product");
      stan::math::validate_vector(v2,"dot_product");
      stan::math::validate_matching_sizes(v1,v2,"dot_product");
      return var(new dot_product_vv_vari(v1,v2));
    }
    /**
     * Returns the dot product.
     *
     * @param[in] v1 First column vector.
     * @param[in] v2 Second column vector.
     * @return Dot product of the vectors.
     * @throw std::domain_error if length of v1 is not equal to length of v2
     * or either v1 or v2 are not vectors.
     */
    template<int R1,int C1,int R2, int C2>
    inline var dot_product(const Eigen::Matrix<var, R1, C1>& v1, 
                           const Eigen::Matrix<double, R2, C2>& v2) {
      stan::math::validate_vector(v1,"dot_product");
      stan::math::validate_vector(v2,"dot_product");
      stan::math::validate_matching_sizes(v1,v2,"dot_product");
      return var(new dot_product_vd_vari(v1,v2));
    }
    /**
     * Returns the dot product.
     *
     * @param[in] v1 First column vector.
     * @param[in] v2 Second column vector.
     * @return Dot product of the vectors.
     * @throw std::domain_error if length of v1 is not equal to length of v2
     * or either v1 or v2 are not vectors.
     */
    template<int R1,int C1,int R2, int C2>
    inline var dot_product(const Eigen::Matrix<double, R1, C1>& v1, 
                           const Eigen::Matrix<var, R2, C2>& v2) {
      stan::math::validate_vector(v1,"dot_product");
      stan::math::validate_vector(v2,"dot_product");
      stan::math::validate_matching_sizes(v1,v2,"dot_product");
      return var(new dot_product_vd_vari(v2,v1));
    }
    /**
     * Returns the dot product.
     *
     * @param[in] v1 First array.
     * @param[in] v2 Second array.
     * @param[in] length Length of both arrays.
     * @return Dot product of the arrays.
     */
    inline var dot_product(const var* v1, const var* v2, size_t length) {
      return var(new dot_product_vv_vari(v1, v2, length));
    }
    /**
     * Returns the dot product.
     *
     * @param[in] v1 First array.
     * @param[in] v2 Second array.
     * @param[in] length Length of both arrays.
     * @return Dot product of the arrays.
     */
    inline var dot_product(const var* v1, const double* v2, size_t length) {
      return var(new dot_product_vd_vari(v1, v2, length));
    }
    /**
     * Returns the dot product.
     *
     * @param[in] v1 First array.
     * @param[in] v2 Second array.
     * @param[in] length Length of both arrays.
     * @return Dot product of the arrays.
     */
    inline var dot_product(const double* v1, const var* v2, size_t length) {
      return var(new dot_product_vd_vari(v2, v1, length));
    }
    /**
     * Returns the dot product.
     *
     * @param[in] v1 First vector.
     * @param[in] v2 Second vector.
     * @return Dot product of the vectors.
     * @throw std::domain_error if sizes of v1 and v2 do not match.
     */
    inline var dot_product(const std::vector<var>& v1,
                           const std::vector<var>& v2) {
      stan::math::validate_matching_sizes(v1,v2,"dot_product");
      return var(new dot_product_vv_vari(&v1[0], &v2[0], v1.size()));
    }
    /**
     * Returns the dot product.
     *
     * @param[in] v1 First vector.
     * @param[in] v2 Second vector.
     * @return Dot product of the vectors.
     * @throw std::domain_error if sizes of v1 and v2 do not match.
     */
    inline var dot_product(const std::vector<var>& v1,
                           const std::vector<double>& v2) {
      stan::math::validate_matching_sizes(v1,v2,"dot_product");
      return var(new dot_product_vd_vari(&v1[0], &v2[0], v1.size()));
    }
    /**
     * Returns the dot product.
     *
     * @param[in] v1 First vector.
     * @param[in] v2 Second vector.
     * @return Dot product of the vectors.
     * @throw std::domain_error if sizes of v1 and v2 do not match.
     */
    inline var dot_product(const std::vector<double>& v1,
                           const std::vector<var>& v2) {
      stan::math::validate_matching_sizes(v1,v2,"dot_product");
      return var(new dot_product_vd_vari(&v2[0], &v1[0], v1.size()));
    }

    /**
     * Returns the sum of the coefficients of the specified
     * matrix, column vector or row vector.
     * @param m Specified matrix or vector.
     * @return Sum of coefficients of matrix.
     */
    template <int R, int C>
    inline var sum(const Eigen::Matrix<var,R,C>& m) {
      if (m.size() == 0)
        return 0.0;
      return var(new sum_v_vari(m));
    }

    template <int R1,int C1,int R2,int C2>
    class mdivide_left_vv_vari : public vari {
    public:
      int _M; // A.rows() = A.cols() = B.rows()
      int _N; // B.cols()
      double* _A;
      double* _C;
      vari** _variRefA;
      vari** _variRefB;
      vari** _variRefC;

      mdivide_left_vv_vari(const Eigen::Matrix<var,R1,C1> &A,
                           const Eigen::Matrix<var,R2,C2> &B)
      : vari(0.0),
 	_M(A.rows()),
	_N(B.cols()),
	_A((double*)stan::agrad::memalloc_.alloc(sizeof(double) 
                                                 * A.rows() * A.cols())),
	_C((double*)stan::agrad::memalloc_.alloc(sizeof(double) 
                                                 * B.rows() * B.cols())),
        _variRefA((vari**)stan::agrad::memalloc_.alloc(sizeof(vari*) 
						       * A.rows() * A.cols())),
        _variRefB((vari**)stan::agrad::memalloc_.alloc(sizeof(vari*) 
						       * B.rows() * B.cols())),
        _variRefC((vari**)stan::agrad::memalloc_.alloc(sizeof(vari*) 
                                                      * B.rows() * B.cols()))
      {
	using Eigen::Matrix;
        using Eigen::Map;

	size_t pos = 0;
	for (size_type j = 0; j < _M; j++) {
          for (size_type i = 0; i < _M; i++) {
            _variRefA[pos] = A(i,j).vi_;
	    _A[pos++] = A(i,j).val();
          }
        }
	
	pos = 0;
	for (size_type j = 0; j < _N; j++) {
          for (size_type i = 0; i < _M; i++) {
            _variRefB[pos] = B(i,j).vi_;
            _C[pos++] = B(i,j).val();
          }
        }
        
	Matrix<double,R1,C2> C(_M,_N);
	C = Map<Matrix<double,R1,C2> >(_C,_M,_N);

	C = Map<Matrix<double,R1,C1> >(_A,_M,_M)
	  .colPivHouseholderQr().solve(C);

	pos = 0;
        for (size_type j = 0; j < _N; j++) {
          for (size_type i = 0; i < _M; i++) {
	    _C[pos] = C(i,j);
            _variRefC[pos] = new vari(_C[pos],false);
	    pos++;
          }
        }
      }
      
      virtual void chain() {
	using Eigen::Matrix;
        using Eigen::Map;
	Eigen::Matrix<double,R1,C1> adjA(_M,_M);
        Eigen::Matrix<double,R2,C2> adjB(_M,_N);
        Eigen::Matrix<double,R1,C2> adjC(_M,_N);

	size_t pos = 0;
        for (size_type j = 0; j < adjC.cols(); j++)
          for (size_type i = 0; i < adjC.rows(); i++)
            adjC(i,j) = _variRefC[pos++]->adj_;
        
        
	adjB = Map<Matrix<double,R1,C1> >(_A,_M,_M)
	  .transpose().colPivHouseholderQr().solve(adjC);
	adjA.noalias() = -adjB
	  * Map<Matrix<double,R1,C2> >(_C,_M,_N).transpose();
        
	pos = 0;
        for (size_type j = 0; j < adjA.cols(); j++)
          for (size_type i = 0; i < adjA.rows(); i++)
            _variRefA[pos++]->adj_ += adjA(i,j);
        
	pos = 0;
        for (size_type j = 0; j < adjB.cols(); j++)
          for (size_type i = 0; i < adjB.rows(); i++)
            _variRefB[pos++]->adj_ += adjB(i,j);
      }
    };
    
    template <int R1,int C1,int R2,int C2>
    class mdivide_left_dv_vari : public vari {
    public:
      int _M; // A.rows() = A.cols() = B.rows()
      int _N; // B.cols()
      double* _A;
      double* _C;
      vari** _variRefB;
      vari** _variRefC;
      
      mdivide_left_dv_vari(const Eigen::Matrix<double,R1,C1> &A,
                           const Eigen::Matrix<var,R2,C2> &B)
      : vari(0.0),
 	_M(A.rows()),
	_N(B.cols()),
	_A((double*)stan::agrad::memalloc_.alloc(sizeof(double) 
                                                 * A.rows() * A.cols())),
	_C((double*)stan::agrad::memalloc_.alloc(sizeof(double) 
                                                 * B.rows() * B.cols())),
        _variRefB((vari**)stan::agrad::memalloc_.alloc(sizeof(vari*) 
						       * B.rows() * B.cols())),
        _variRefC((vari**)stan::agrad::memalloc_.alloc(sizeof(vari*) 
						       * B.rows() * B.cols()))
      {
	using Eigen::Matrix;
	using Eigen::Map;
	
	size_t pos = 0;
	for (size_type j = 0; j < _M; j++) {
	  for (size_type i = 0; i < _M; i++) {
	    _A[pos++] = A(i,j);
	  }
	}
	
	pos = 0;
	for (size_type j = 0; j < _N; j++) {
          for (size_type i = 0; i < _M; i++) {
            _variRefB[pos] = B(i,j).vi_;
            _C[pos++] = B(i,j).val();
          }
        }
                
	Matrix<double,R1,C2> C(_M,_N);
	C = Map<Matrix<double,R1,C2> >(_C,_M,_N);

	C = Map<Matrix<double,R1,C1> >(_A,_M,_M)
	  .colPivHouseholderQr().solve(C);
	
	pos = 0;
        for (size_type j = 0; j < _N; j++) {
          for (size_type i = 0; i < _M; i++) {
	    _C[pos] = C(i,j);
            _variRefC[pos] = new vari(_C[pos],false);
	    pos++;
          }
        }
      }
      
      virtual void chain() {
	using Eigen::Matrix;
        using Eigen::Map;
        Eigen::Matrix<double,R2,C2> adjB(_M,_N);
        Eigen::Matrix<double,R1,C2> adjC(_M,_N);

	size_t pos = 0;
        for (size_type j = 0; j < adjC.cols(); j++)
          for (size_type i = 0; i < adjC.rows(); i++)
            adjC(i,j) = _variRefC[pos++]->adj_;

        adjB = Map<Matrix<double,R1,C1> >(_A,_M,_M)
	  .transpose().colPivHouseholderQr().solve(adjC);

	pos = 0;
        for (size_type j = 0; j < adjB.cols(); j++)
          for (size_type i = 0; i < adjB.rows(); i++)
            _variRefB[pos++]->adj_ += adjB(i,j);
      }
    };
    
    template <int R1,int C1,int R2,int C2>
    class mdivide_left_vd_vari : public vari {
    public:
      Eigen::Matrix<double,R1,C1> _A;
      Eigen::Matrix<double,R1,C2> _C;
      
      Eigen::Matrix<vari*,R1,C1> _variRefA;
      Eigen::Matrix<vari*,R1,C2> _variRefC;
      
      mdivide_left_vd_vari(const Eigen::Matrix<var,R1,C1> &A,
                           const Eigen::Matrix<double,R2,C2> &B)
      : vari(0.0),
      _A(A.rows(),A.cols()), _C(B.rows(),B.cols()),
      _variRefA(A.rows(),A.cols()), _variRefC(B.rows(),B.cols())
      {
        for (size_type j = 0; j < _variRefA.cols(); j++) {
          for (size_type i = 0; i < _variRefA.rows(); i++) {
            _variRefA(i,j) = A(i,j).vi_;
            _A(i,j) = A(i,j).val();
          }
        }
        
        
        _C = _A.colPivHouseholderQr().solve(B);
        
        for (size_type j = 0; j < _variRefC.cols(); j++) {
          for (size_type i = 0; i < _variRefC.rows(); i++) {
            _variRefC(i,j) = new vari(_C(i,j),false);
          }
        }
      }
      
      virtual void chain() {
        Eigen::Matrix<double,R1,C1> adjA(_variRefA.rows(),_variRefA.cols());
        Eigen::Matrix<double,R1,C2> adjC(_variRefC.rows(),_variRefC.cols());
        
        for (size_type j = 0; j < adjC.cols(); j++)
          for (size_type i = 0; i < adjC.rows(); i++)
            adjC(i,j) = _variRefC(i,j)->adj_;
        
        adjA = -_A.transpose().colPivHouseholderQr().solve(adjC*_C.transpose());
        
        for (size_type j = 0; j < adjA.cols(); j++)
          for (size_type i = 0; i < adjA.rows(); i++)
            _variRefA(i,j)->adj_ += adjA(i,j);
      }
    };

    template <int R1,int C1,int R2,int C2>
    inline 
    Eigen::Matrix<var,R1,C2>
    mdivide_left(const Eigen::Matrix<var,R1,C1> &A,
                 const Eigen::Matrix<var,R2,C2> &b) {
      Eigen::Matrix<var,R1,C2> res(b.rows(),b.cols());
      
      stan::math::validate_square(A,"mdivide_left");
      stan::math::validate_multiplicable(A,b,"mdivide_left");
      
      // NOTE: this is not a memory leak, this vari is used in the 
      // expression graph to evaluate the adjoint, but is not needed
      // for the returned matrix.  Memory will be cleaned up with the arena allocator.
      mdivide_left_vv_vari<R1,C1,R2,C2> *baseVari = new mdivide_left_vv_vari<R1,C1,R2,C2>(A,b);
      
      size_t pos = 0;
      for (size_type j = 0; j < res.cols(); j++)
	for (size_type i = 0; i < res.rows(); i++)
          res(i,j).vi_ = baseVari->_variRefC[pos++];
      
      return res;
    }

    template <int R1,int C1,int R2,int C2>
    inline 
    Eigen::Matrix<var,R1,C2>
    mdivide_left(const Eigen::Matrix<var,R1,C1> &A,
                 const Eigen::Matrix<double,R2,C2> &b) {
      Eigen::Matrix<var,R1,C2> res(b.rows(),b.cols());
      
      stan::math::validate_square(A,"mdivide_left");
      stan::math::validate_multiplicable(A,b,"mdivide_left");
      
      // NOTE: this is not a memory leak, this vari is used in the 
      // expression graph to evaluate the adjoint, but is not needed
      // for the returned matrix.  Memory will be cleaned up with the arena allocator.
      mdivide_left_vd_vari<R1,C1,R2,C2> *baseVari = new mdivide_left_vd_vari<R1,C1,R2,C2>(A,b);
      
      for (size_type i = 0; i < res.rows(); i++)
        for (size_type j = 0; j < res.cols(); j++)
          res(i,j).vi_ = baseVari->_variRefC(i,j);
      
      return res;
    }
    
    template <int R1,int C1,int R2,int C2>
    inline 
    Eigen::Matrix<var,R1,C2>
    mdivide_left(const Eigen::Matrix<double,R1,C1> &A,
                 const Eigen::Matrix<var,R2,C2> &b) {
      Eigen::Matrix<var,R1,C2> res(b.rows(),b.cols());
      
      stan::math::validate_square(A,"mdivide_left");
      stan::math::validate_multiplicable(A,b,"mdivide_left");
      
      // NOTE: this is not a memory leak, this vari is used in the 
      // expression graph to evaluate the adjoint, but is not needed
      // for the returned matrix.  Memory will be cleaned up with the arena allocator.<commit_msg>fixed memory leak in mdivide_left_vd_vari<commit_after>#ifndef __STAN__AGRAD__MATRIX_HPP__
#define __STAN__AGRAD__MATRIX_HPP__

#include <stan/math/matrix.hpp>
#include <stan/math/matrix_error_handling.hpp>

#include <stan/agrad/special_functions.hpp>

namespace stan {
  namespace agrad {
    typedef Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic>::size_type size_type;

    void fill(var& x, const var& y) {
      x = y;
    }
    template <typename T>
    void fill(std::vector<T>& x, const var& y) {
      for (size_t i = 0; i < x.size(); ++i)
        fill(x[i],y);
    }
    template <int R, int C>
    void fill(Eigen::Matrix<var,R,C>& x, const var& y) {
      x.fill(y);
    }



    class gevv_vvv_vari : public stan::agrad::vari {
    protected:
      stan::agrad::vari* alpha_;
      stan::agrad::vari** v1_;
      stan::agrad::vari** v2_;
      double dotval_;
      size_t length_;
      inline static double eval_gevv(const stan::agrad::var* alpha,
                                     const stan::agrad::var* v1, int stride1,
                                     const stan::agrad::var* v2, int stride2,
                                     size_t length, double *dotprod) {
        double result = 0;
        for (size_t i = 0; i < length; i++)
          result += v1[i*stride1].vi_->val_ * v2[i*stride2].vi_->val_;
        *dotprod = result;
        return alpha->vi_->val_ * result;
      }
    public:
      gevv_vvv_vari(const stan::agrad::var* alpha, 
                    const stan::agrad::var* v1, int stride1, 
                    const stan::agrad::var* v2, int stride2, size_t length) : 
        vari(eval_gevv(alpha,v1,stride1,v2,stride2,length,&dotval_)), length_(length) {
        alpha_ = alpha->vi_;
        v1_ = (stan::agrad::vari**)stan::agrad::memalloc_.alloc(2*length_*sizeof(stan::agrad::vari*));
        v2_ = v1_ + length_;
        for (size_t i = 0; i < length_; i++)
          v1_[i] = v1[i*stride1].vi_;
        for (size_t i = 0; i < length_; i++)
          v2_[i] = v2[i*stride2].vi_;
      }
      void chain() {
        const double adj_alpha = adj_ * alpha_->val_;
        for (size_t i = 0; i < length_; i++) {
          v1_[i]->adj_ += adj_alpha * v2_[i]->val_;
          v2_[i]->adj_ += adj_alpha * v1_[i]->val_;
        }
        alpha_->adj_ += adj_ * dotval_;
      }
    };
  }
}

namespace Eigen {

  /**
   * Numerical traits template override for Eigen for automatic
   * gradient variables.
   */
  template <> struct NumTraits<stan::agrad::var>
  {
    /**
     * Real-valued variables.
     *
     * Required for numerical traits.
     */
    typedef stan::agrad::var Real;

    /**
     * Non-integer valued variables.
     *
     * Required for numerical traits.
     */
    typedef stan::agrad::var NonInteger;

    /**
     * Nested variables.
     *
     * Required for numerical traits.
     */
    typedef stan::agrad::var Nested;

    /**
     * Return standard library's epsilon for double-precision floating
     * point, <code>std::numeric_limits&lt;double&gt;::epsilon()</code>.
     *
     * @return Same epsilon as a <code>double</code>.
     */
    inline static Real epsilon() { 
      return std::numeric_limits<double>::epsilon(); 
    }

    /**
     * Return dummy precision
     */
    inline static Real dummy_precision() {
      return 1e-12; // copied from NumTraits.h values for double
    }

    /**
     * Return standard library's highest for double-precision floating
     * point, <code>std::numeric_limits&lt;double&gt;::max()</code>.
     *
     * @return Same highest value as a <code>double</code>.
     */
    inline static Real highest() {
      return std::numeric_limits<double>::max();
    }

    /**
     * Return standard library's lowest for double-precision floating
     * point, <code>&#45;std::numeric_limits&lt;double&gt;::max()</code>.
     *
     * @return Same lowest value as a <code>double</code>.
     */    
    inline static Real lowest() {
      return -std::numeric_limits<double>::max();
    }

    /**
     * Properties for automatic differentiation variables
     * read by Eigen matrix library.
     */
    enum {
      IsInteger = 0,
      IsSigned = 1,
      IsComplex = 0,
      RequireInitialization = 0,
      ReadCost = 1,
      AddCost = 1,
      MulCost = 1,
      HasFloatingPoint = 1
    };
  };

  namespace internal {
    /**
     * Implemented this for printing to stream.
     */
    template<>
    struct significant_decimals_default_impl<stan::agrad::var,false>
    {
      static inline int run()
      {
        using std::ceil;
        return cast<double,int>(ceil(-log(NumTraits<stan::agrad::var>::epsilon().val())
                                     /log(10.0)));
      }
    };

    /**
     * Scalar product traits override for Eigen for automatic
     * gradient variables.
     */
    template <>  
    struct scalar_product_traits<stan::agrad::var,double> {
      typedef stan::agrad::var ReturnType;
    };

    /**
     * Scalar product traits override for Eigen for automatic
     * gradient variables.
     */
    template <>  
    struct scalar_product_traits<double,stan::agrad::var> {
      typedef stan::agrad::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::agrad::var,ColMajor,ConjugateLhs,stan::agrad::var,ConjugateRhs>
    {
      typedef stan::agrad::var LhsScalar;
      typedef stan::agrad::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::agrad::var(new stan::agrad::gevv_vvv_vari(&alpha,((int)LhsStorageOrder == (int)ColMajor)?(&lhs[i]):(&lhs[i*lhsStride]),((int)LhsStorageOrder == (int)ColMajor)?(lhsStride):(1),rhs,rhsIncr,cols));
        }
      }
    };
    template<typename Index, bool ConjugateLhs, bool ConjugateRhs>
    struct general_matrix_vector_product<Index,stan::agrad::var,RowMajor,ConjugateLhs,stan::agrad::var,ConjugateRhs>
    {
      typedef stan::agrad::var LhsScalar;
      typedef stan::agrad::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::agrad::var(new stan::agrad::gevv_vvv_vari(&alpha,((int)LhsStorageOrder == (int)ColMajor)?(&lhs[i]):(&lhs[i*lhsStride]),((int)LhsStorageOrder == (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::agrad::var,LhsStorageOrder,ConjugateLhs,stan::agrad::var,RhsStorageOrder,ConjugateRhs,ColMajor>
    {
      typedef stan::agrad::var LhsScalar;
      typedef stan::agrad::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[((int)RhsStorageOrder == (int)ColMajor)?(i*rhsStride):(i)],((int)RhsStorageOrder == (int)ColMajor)?(1):(rhsStride),
              &res[i*resStride],1,alpha);
        }
      }
    };
  }
}

namespace stan {

  namespace agrad {

    /**
     * The type of a matrix holding <code>stan::agrad::var</code>
     * values.
     */
    typedef 
    Eigen::Matrix<var,Eigen::Dynamic,Eigen::Dynamic>
    matrix_v;

    /**
     * The type of a (column) vector holding <code>stan::agrad::var</code>
     * values.
     */
    typedef 
    Eigen::Matrix<var,Eigen::Dynamic,1>
    vector_v;

    /**
     * The type of a row vector holding <code>stan::agrad::var</code>
     * values.
     */
    typedef 
    Eigen::Matrix<var,1,Eigen::Dynamic>
    row_vector_v;

    /**
     * Initialize variable to value.  (Function may look pointless, but
     * its needed to bottom out recursion.)
     */
    inline void initialize_variable(var& variable, const var& value) {
      variable = value;
    }

    /**
     * Initialize every cell in the matrix to the specified value.
     * 
     */
    template <int R, int C>
    inline void initialize_variable(Eigen::Matrix<var,R,C>& matrix, const var& value) {
      for (int i = 0; i < matrix.size(); ++i)
        matrix(i) = value;
    }

    /**
     * Initialize the variables in the standard vector recursively.
     */
    template <typename T>
    inline void initialize_variable(std::vector<T>& variables, const var& value) {
      for (size_t i = 0; i < variables.size(); ++i)
        initialize_variable(variables[i],value);
    }


    /**
     * Converts argument to an automatic differentiation variable.
     *
     * Returns a stan::agrad::var variable with the input value.
     *
     * @param[in] x A scalar value
     * @return An automatic differentiation variable with the input value.
     */
    inline var to_var(const double& x) {
      return var(x);
    }
    /**
     * Converts argument to an automatic differentiation variable.
     *
     * Returns a stan::agrad::var variable with the input value.
     *
     * @param[in] x An automatic differentiation variable.
     * @return An automatic differentiation variable with the input value.
     */    
    inline var to_var(const var& x) {
      return x;
    }
    /**
     * Converts argument to an automatic differentiation variable.
     *
     * Returns a stan::agrad::var variable with the input value.
     *
     * @param[in] m A Matrix with scalars
     * @return A Matrix with automatic differentiation variables
     */
    inline matrix_v to_var(const stan::math::matrix_d& m) {
      matrix_v m_v(m.rows(), m.cols());
      for (int j = 0; j < m.cols(); ++j)
        for (int i = 0; i < m.rows(); ++i)
          m_v(i,j) = m(i,j);
      return m_v;
    }
    /**
     * Converts argument to an automatic differentiation variable.
     *
     * Returns a stan::agrad::var variable with the input value.
     * 
     * @param[in] m A Matrix with automatic differentiation variables.
     * @return A Matrix with automatic differentiation variables.
     */
    inline matrix_v to_var(const matrix_v& m) {
      return m;
    }
    /**
     * Converts argument to an automatic differentiation variable.
     *
     * Returns a stan::agrad::var variable with the input value.     
     *
     * @param[in] v A Vector of scalars
     * @return A Vector of automatic differentiation variables with
     *   values of v
     */
    inline vector_v to_var(const stan::math::vector_d& v) {
      vector_v v_v(v.size());
      for (int i = 0; i < v.size(); ++i)
        v_v[i] = v[i];
      return v_v;
    }
    /**
     * Converts argument to an automatic differentiation variable.
     *
     * Returns a stan::agrad::var variable with the input value.     
     *
     * @param[in] v A Vector of automatic differentiation variables
     * @return A Vector of automatic differentiation variables with
     *   values of v
     */
    inline vector_v to_var(const vector_v& v) {
      return v;
    }
    /**
     * Converts argument to an automatic differentiation variable.
     *
     * Returns a stan::agrad::var variable with the input value.     
     *
     * @param[in] rv A row vector of scalars
     * @return A row vector of automatic differentation variables with 
     *   values of rv.
     */
    inline row_vector_v to_var(const stan::math::row_vector_d& rv) {
      row_vector_v rv_v(rv.size());
      for (int i = 0; i < rv.size(); ++i)
        rv_v[i] = rv[i];
      return rv_v;
    }
    /**
     * Converts argument to an automatic differentiation variable.
     *
     * Returns a stan::agrad::var variable with the input value.     
     *
     * @param[in] rv A row vector with automatic differentiation variables
     * @return A row vector with automatic differentiation variables
     *    with values of rv.
     */
    inline row_vector_v to_var(const row_vector_v& rv) {
      return rv;
    }



    // scalar returns

    // /**
    //  * Determinant of the matrix.
    //  *
    //  * Returns the determinant of the specified
    //  * square matrix.
    //  *
    //  * @param m Specified matrix.
    //  * @return Determinant of the matrix.
    //  * @throw std::domain_error if m is not a square matrix
    //  */
    // var determinant(const Eigen::Matrix<var, Eigen::Dynamic, Eigen::Dynamic>& m);
    
    namespace {
      class dot_self_vari : public vari {
      protected:
        vari** v_;
        size_t size_;
      public:
        dot_self_vari(vari** v, size_t size) 
          : vari(var_dot_self(v,size)), 
            v_(v),
            size_(size) {
        }
        template<typename Derived>
        dot_self_vari(const Eigen::DenseBase<Derived> &v) : 
          vari(var_dot_self(v)), size_(v.size()) {
          v_ = (vari**)memalloc_.alloc(size_*sizeof(vari*));
          for (size_t i = 0; i < size_; i++)
            v_[i] = v[i].vi_;
        }
        template <int R, int C>
        dot_self_vari(const Eigen::Matrix<var,R,C>& v) :
          vari(var_dot_self(v)), size_(v.size()) {
          v_ = (vari**) memalloc_.alloc(size_ * sizeof(vari*));
          for (size_t i = 0; i < size_; ++i)
            v_[i] = v(i).vi_;
        }
        inline static double square(double x) { return x * x; }
        inline static double var_dot_self(vari** v, size_t size) {
          double sum = 0.0;
          for (size_t i = 0; i < size; ++i)
            sum += square(v[i]->val_);
          return sum;
        }
        template<typename Derived>
        double var_dot_self(const Eigen::DenseBase<Derived> &v) {
          double sum = 0.0;
          for (int i = 0; i < v.size(); ++i)
            sum += square(v(i).vi_->val_);
          return sum;
        }
        template <int R, int C>
        inline static double var_dot_self(const Eigen::Matrix<var,R,C> &v) {
          double sum = 0.0;
          for (int i = 0; i < v.size(); ++i)
            sum += square(v(i).vi_->val_);
          return sum;
        }
        virtual void chain() {
          for (size_t i = 0; i < size_; ++i) 
            v_[i]->adj_ += adj_ * 2.0 * v_[i]->val_;
        }
      };

      class sum_v_vari : public vari{
      protected:
        vari** v_;
        size_t length_;
        inline static double var_sum(const var *v, size_t length) {
          double result = 0;
          for (size_t i = 0; i < length; i++)
            result += v[i].vi_->val_;
          return result;
        } 
        template<typename Derived>
        inline static double var_sum(const Eigen::DenseBase<Derived> &v) {
          double result = 0;
          for (int i = 0; i < v.size(); i++)
            result += v(i).vi_->val_;
          return result;
        } 
      public:
        template<typename Derived>
        sum_v_vari(const Eigen::DenseBase<Derived> &v) :
          vari(var_sum(v)), length_(v.size()) {
          v_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
          for (size_t i = 0; i < length_; i++)
            v_[i] = v(i).vi_;
        }
        template<int R1,int C1>
        sum_v_vari(const Eigen::Matrix<var,R1,C1> &v1) :
          vari(var_sum(v1)), length_(v1.size()) {
          v_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
          for (size_t i = 0; i < length_; i++)
            v_[i] = v1(i).vi_;
        }
        sum_v_vari(const var *v, size_t len) :
          vari(var_sum(v,len)), length_(len) {
          v_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
          for (size_t i = 0; i < length_; i++)
            v_[i] = v[i].vi_;
        }
        virtual void chain() {
          for (size_t i = 0; i < length_; i++) {
            v_[i]->adj_ += adj_;
          }
        }
      };
      class dot_product_vv_vari : public vari {
      protected:
        vari** v1_;
        vari** v2_;
        size_t length_;
        inline static double var_dot(const var* v1, const var* v2,
                                     size_t length) {
          double result = 0;
          for (size_t i = 0; i < length; i++)
            result += v1[i].vi_->val_ * v2[i].vi_->val_;
          return result;
        }
        template<typename Derived1,typename Derived2>
        inline static double var_dot(const Eigen::DenseBase<Derived1> &v1,
                                     const Eigen::DenseBase<Derived2> &v2) {
          double result = 0;
          for (int i = 0; i < v1.size(); i++)
            result += v1[i].vi_->val_ * v2[i].vi_->val_;
          return result;
        }
        inline static double var_dot(vari** v1, vari** v2, size_t length) {
          double result = 0;
          for (size_t i = 0; i < length; ++i)
            result += v1[i]->val_ * v2[i]->val_;
          return result;
        }
      public:
        dot_product_vv_vari(vari** v1, vari** v2, size_t length)
          : vari(var_dot(v1,v2,length)),
            v1_(v1), 
            v2_(v2), 
            length_(length) {

        }
        dot_product_vv_vari(const var* v1, const var* v2, size_t length,
                            dot_product_vv_vari* shared_v1 = NULL,
                            dot_product_vv_vari* shared_v2 = NULL) : 
          vari(var_dot(v1, v2, length)), length_(length) {
          if (shared_v1 == NULL) {
            v1_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
            for (size_t i = 0; i < length_; i++)
              v1_[i] = v1[i].vi_;
          }
          else {
            v1_ = shared_v1->v1_;
          }
          if (shared_v2 == NULL) {
            v2_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
            for (size_t i = 0; i < length_; i++)
              v2_[i] = v2[i].vi_;
          }
          else {
            v2_ = shared_v2->v2_;
          }
        }
        template<typename Derived1,typename Derived2>
        dot_product_vv_vari(const Eigen::DenseBase<Derived1> &v1,
                            const Eigen::DenseBase<Derived2> &v2,
                            dot_product_vv_vari* shared_v1 = NULL,
                            dot_product_vv_vari* shared_v2 = NULL) : 
          vari(var_dot(v1, v2)), length_(v1.size()) {
          if (shared_v1 == NULL) {
            v1_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
            for (size_t i = 0; i < length_; i++)
              v1_[i] = v1[i].vi_;
          }
          else {
            v1_ = shared_v1->v1_;
          }
          if (shared_v2 == NULL) {
            v2_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
            for (size_t i = 0; i < length_; i++)
              v2_[i] = v2[i].vi_;
          }
          else {
            v2_ = shared_v2->v2_;
          }
        }
        template<int R1,int C1,int R2,int C2>
        dot_product_vv_vari(const Eigen::Matrix<var,R1,C1> &v1,
                            const Eigen::Matrix<var,R2,C2> &v2,
                            dot_product_vv_vari* shared_v1 = NULL,
                            dot_product_vv_vari* shared_v2 = NULL) : 
          vari(var_dot(v1, v2)), length_(v1.size()) {
          if (shared_v1 == NULL) {
            v1_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
            for (size_t i = 0; i < length_; i++)
              v1_[i] = v1[i].vi_;
          }
          else {
            v1_ = shared_v1->v1_;
          }
          if (shared_v2 == NULL) {
            v2_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
            for (size_t i = 0; i < length_; i++)
              v2_[i] = v2[i].vi_;
          }
          else {
            v2_ = shared_v2->v2_;
          }
        }
        virtual void chain() {
          for (size_t i = 0; i < length_; i++) {
            v1_[i]->adj_ += adj_ * v2_[i]->val_;
            v2_[i]->adj_ += adj_ * v1_[i]->val_;
          }
        }
      };

      class dot_product_vd_vari : public vari {
      protected:
        vari** v1_;
        double* v2_;
        size_t length_;
        inline static double var_dot(const var* v1, const double* v2,
                                     size_t length) {
          double result = 0;
          for (size_t i = 0; i < length; i++)
            result += v1[i].vi_->val_ * v2[i];
          return result;
        }
        template<typename Derived1,typename Derived2>
        inline static double var_dot(const Eigen::DenseBase<Derived1> &v1,
                                     const Eigen::DenseBase<Derived2> &v2) {
          double result = 0;
          for (int i = 0; i < v1.size(); i++)
            result += v1[i].vi_->val_ * v2[i];
          return result;
        }
      public:
        dot_product_vd_vari(const var* v1, const double* v2, size_t length,
                            dot_product_vd_vari *shared_v1 = NULL,
                            dot_product_vd_vari *shared_v2 = NULL) : 
          vari(var_dot(v1, v2, length)), length_(length) {
          if (shared_v1 == NULL) {
            v1_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
            for (size_t i = 0; i < length_; i++)
              v1_[i] = v1[i].vi_;
          } else {
            v1_ = shared_v1->v1_;
          }
          if (shared_v2 == NULL) {
            v2_ = (double*)memalloc_.alloc(length_*sizeof(double));
            for (size_t i = 0; i < length_; i++)
              v2_[i] = v2[i];
          } else {
            v2_ = shared_v2->v2_;
          }
        }
        template<typename Derived1,typename Derived2>
        dot_product_vd_vari(const Eigen::DenseBase<Derived1> &v1,
                            const Eigen::DenseBase<Derived2> &v2,
                            dot_product_vd_vari *shared_v1 = NULL,
                            dot_product_vd_vari *shared_v2 = NULL) : 
          vari(var_dot(v1, v2)), length_(v1.size()) {
          if (shared_v1 == NULL) {
            v1_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
            for (size_t i = 0; i < length_; i++)
              v1_[i] = v1[i].vi_;
          } else {
            v1_ = shared_v1->v1_;
          }
          if (shared_v2 == NULL) {
            v2_ = (double*)memalloc_.alloc(length_*sizeof(double));
            for (size_t i = 0; i < length_; i++)
              v2_[i] = v2[i];
          } else {
            v2_ = shared_v2->v2_;
          }
        }
        template<int R1,int C1,int R2,int C2>
        dot_product_vd_vari(const Eigen::Matrix<var,R1,C1> &v1,
                            const Eigen::Matrix<double,R2,C2> &v2,
                            dot_product_vd_vari *shared_v1 = NULL,
                            dot_product_vd_vari *shared_v2 = NULL) : 
          vari(var_dot(v1, v2)), length_(v1.size()) {
          if (shared_v1 == NULL) {
            v1_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
            for (size_t i = 0; i < length_; i++)
              v1_[i] = v1[i].vi_;
          } else {
            v1_ = shared_v1->v1_;
          }
          if (shared_v2 == NULL) {
            v2_ = (double*)memalloc_.alloc(length_*sizeof(double));
            for (size_t i = 0; i < length_; i++)
              v2_[i] = v2[i];
          } else {
            v2_ = shared_v2->v2_;
          }
        }
        virtual void chain() {
          for (size_t i = 0; i < length_; i++) {
            v1_[i]->adj_ += adj_ * v2_[i];
          }
        }
      };

    }

    /**
     * Returns the dot product of a vector with itself.
     *
     * @param[in] v Vector.
     * @return Dot product of the vector with itself.
     * @tparam R number of rows or <code>Eigen::Dynamic</code> for
     * dynamic; one of R or C must be 1
     * @tparam C number of rows or <code>Eigen::Dyanmic</code> for
     * dynamic; one of R or C must be 1
     */
    template<int R, int C>
    inline var dot_self(const Eigen::Matrix<var, R, C>& v) {
      stan::math::validate_vector(v,"dot_self");
      return var(new dot_self_vari(v));
    }
    
    /**
     * Returns the dot product of each column of a matrix with itself.
     * @param x Matrix.
     * @tparam T scalar type
     */
    template<int R,int C>
    inline Eigen::Matrix<var,1,C> 
    columns_dot_self(const Eigen::Matrix<var,R,C>& x) {
      Eigen::Matrix<var,1,C> ret(1,x.cols());
      for (size_type i = 0; i < x.cols(); i++) {
        ret(i) = var(new dot_self_vari(x.col(i)));
      }
      return ret;
    }
    
    template<int R1,int C1,int R2, int C2>
    inline Eigen::Matrix<var, 1, C1>
    columns_dot_product(const Eigen::Matrix<var, R1, C1>& v1, 
                        const Eigen::Matrix<var, R2, C2>& v2) {
      stan::math::validate_matching_sizes(v1,v2,"columns_dot_product");
      Eigen::Matrix<var, 1, C1> ret(1,v1.cols());
      for (size_type j = 0; j < v1.cols(); ++j) {
        ret(j) = var(new dot_product_vv_vari(v1.col(j),v2.col(j)));
      }
      return ret;
    }
    
    template<int R1,int C1,int R2, int C2>
    inline Eigen::Matrix<var, 1, C1>
    columns_dot_product(const Eigen::Matrix<var, R1, C1>& v1, 
                        const Eigen::Matrix<double, R2, C2>& v2) {
      stan::math::validate_matching_sizes(v1,v2,"columns_dot_product");
      Eigen::Matrix<var, 1, C1> ret(1,v1.cols());
      size_t j;
      for (j = 0; j < v1.cols(); ++j) {
        ret(j) = var(new dot_product_vd_vari(v1.col(j),v2.col(j)));
      }
      return ret;
    }

    template<int R1,int C1,int R2, int C2>
    inline Eigen::Matrix<var, 1, C1>
    columns_dot_product(const Eigen::Matrix<double, R1, C1>& v1, 
                        const Eigen::Matrix<var, R2, C2>& v2) {
      stan::math::validate_matching_sizes(v1,v2,"columns_dot_product");
      Eigen::Matrix<var, 1, C1> ret(1,v1.cols());
      for (size_type j = 0; j < v1.cols(); ++j) {
        ret(j) = var(new dot_product_vd_vari(v2.col(j),v1.col(j)));
      }
      return ret;
    }
    
    /**
     * Returns the dot product.
     *
     * @param[in] v1 First column vector.
     * @param[in] v2 Second column vector.
     * @return Dot product of the vectors.
     * @throw std::domain_error if length of v1 is not equal to length of v2.
     */
    template<int R1,int C1,int R2, int C2>
    inline var dot_product(const Eigen::Matrix<var, R1, C1>& v1, 
                           const Eigen::Matrix<var, R2, C2>& v2) {
      stan::math::validate_vector(v1,"dot_product");
      stan::math::validate_vector(v2,"dot_product");
      stan::math::validate_matching_sizes(v1,v2,"dot_product");
      return var(new dot_product_vv_vari(v1,v2));
    }
    /**
     * Returns the dot product.
     *
     * @param[in] v1 First column vector.
     * @param[in] v2 Second column vector.
     * @return Dot product of the vectors.
     * @throw std::domain_error if length of v1 is not equal to length of v2
     * or either v1 or v2 are not vectors.
     */
    template<int R1,int C1,int R2, int C2>
    inline var dot_product(const Eigen::Matrix<var, R1, C1>& v1, 
                           const Eigen::Matrix<double, R2, C2>& v2) {
      stan::math::validate_vector(v1,"dot_product");
      stan::math::validate_vector(v2,"dot_product");
      stan::math::validate_matching_sizes(v1,v2,"dot_product");
      return var(new dot_product_vd_vari(v1,v2));
    }
    /**
     * Returns the dot product.
     *
     * @param[in] v1 First column vector.
     * @param[in] v2 Second column vector.
     * @return Dot product of the vectors.
     * @throw std::domain_error if length of v1 is not equal to length of v2
     * or either v1 or v2 are not vectors.
     */
    template<int R1,int C1,int R2, int C2>
    inline var dot_product(const Eigen::Matrix<double, R1, C1>& v1, 
                           const Eigen::Matrix<var, R2, C2>& v2) {
      stan::math::validate_vector(v1,"dot_product");
      stan::math::validate_vector(v2,"dot_product");
      stan::math::validate_matching_sizes(v1,v2,"dot_product");
      return var(new dot_product_vd_vari(v2,v1));
    }
    /**
     * Returns the dot product.
     *
     * @param[in] v1 First array.
     * @param[in] v2 Second array.
     * @param[in] length Length of both arrays.
     * @return Dot product of the arrays.
     */
    inline var dot_product(const var* v1, const var* v2, size_t length) {
      return var(new dot_product_vv_vari(v1, v2, length));
    }
    /**
     * Returns the dot product.
     *
     * @param[in] v1 First array.
     * @param[in] v2 Second array.
     * @param[in] length Length of both arrays.
     * @return Dot product of the arrays.
     */
    inline var dot_product(const var* v1, const double* v2, size_t length) {
      return var(new dot_product_vd_vari(v1, v2, length));
    }
    /**
     * Returns the dot product.
     *
     * @param[in] v1 First array.
     * @param[in] v2 Second array.
     * @param[in] length Length of both arrays.
     * @return Dot product of the arrays.
     */
    inline var dot_product(const double* v1, const var* v2, size_t length) {
      return var(new dot_product_vd_vari(v2, v1, length));
    }
    /**
     * Returns the dot product.
     *
     * @param[in] v1 First vector.
     * @param[in] v2 Second vector.
     * @return Dot product of the vectors.
     * @throw std::domain_error if sizes of v1 and v2 do not match.
     */
    inline var dot_product(const std::vector<var>& v1,
                           const std::vector<var>& v2) {
      stan::math::validate_matching_sizes(v1,v2,"dot_product");
      return var(new dot_product_vv_vari(&v1[0], &v2[0], v1.size()));
    }
    /**
     * Returns the dot product.
     *
     * @param[in] v1 First vector.
     * @param[in] v2 Second vector.
     * @return Dot product of the vectors.
     * @throw std::domain_error if sizes of v1 and v2 do not match.
     */
    inline var dot_product(const std::vector<var>& v1,
                           const std::vector<double>& v2) {
      stan::math::validate_matching_sizes(v1,v2,"dot_product");
      return var(new dot_product_vd_vari(&v1[0], &v2[0], v1.size()));
    }
    /**
     * Returns the dot product.
     *
     * @param[in] v1 First vector.
     * @param[in] v2 Second vector.
     * @return Dot product of the vectors.
     * @throw std::domain_error if sizes of v1 and v2 do not match.
     */
    inline var dot_product(const std::vector<double>& v1,
                           const std::vector<var>& v2) {
      stan::math::validate_matching_sizes(v1,v2,"dot_product");
      return var(new dot_product_vd_vari(&v2[0], &v1[0], v1.size()));
    }

    /**
     * Returns the sum of the coefficients of the specified
     * matrix, column vector or row vector.
     * @param m Specified matrix or vector.
     * @return Sum of coefficients of matrix.
     */
    template <int R, int C>
    inline var sum(const Eigen::Matrix<var,R,C>& m) {
      if (m.size() == 0)
        return 0.0;
      return var(new sum_v_vari(m));
    }

    template <int R1,int C1,int R2,int C2>
    class mdivide_left_vv_vari : public vari {
    public:
      int _M; // A.rows() = A.cols() = B.rows()
      int _N; // B.cols()
      double* _A;
      double* _C;
      vari** _variRefA;
      vari** _variRefB;
      vari** _variRefC;

      mdivide_left_vv_vari(const Eigen::Matrix<var,R1,C1> &A,
                           const Eigen::Matrix<var,R2,C2> &B)
      : vari(0.0),
 	_M(A.rows()),
	_N(B.cols()),
	_A((double*)stan::agrad::memalloc_.alloc(sizeof(double) 
                                                 * A.rows() * A.cols())),
	_C((double*)stan::agrad::memalloc_.alloc(sizeof(double) 
                                                 * B.rows() * B.cols())),
        _variRefA((vari**)stan::agrad::memalloc_.alloc(sizeof(vari*) 
						       * A.rows() * A.cols())),
        _variRefB((vari**)stan::agrad::memalloc_.alloc(sizeof(vari*) 
						       * B.rows() * B.cols())),
        _variRefC((vari**)stan::agrad::memalloc_.alloc(sizeof(vari*) 
                                                      * B.rows() * B.cols()))
      {
	using Eigen::Matrix;
        using Eigen::Map;

	size_t pos = 0;
	for (size_type j = 0; j < _M; j++) {
          for (size_type i = 0; i < _M; i++) {
            _variRefA[pos] = A(i,j).vi_;
	    _A[pos++] = A(i,j).val();
          }
        }
	
	pos = 0;
	for (size_type j = 0; j < _N; j++) {
          for (size_type i = 0; i < _M; i++) {
            _variRefB[pos] = B(i,j).vi_;
            _C[pos++] = B(i,j).val();
          }
        }
        
	Matrix<double,R1,C2> C(_M,_N);
	C = Map<Matrix<double,R1,C2> >(_C,_M,_N);

	C = Map<Matrix<double,R1,C1> >(_A,_M,_M)
	  .colPivHouseholderQr().solve(C);

	pos = 0;
        for (size_type j = 0; j < _N; j++) {
          for (size_type i = 0; i < _M; i++) {
	    _C[pos] = C(i,j);
            _variRefC[pos] = new vari(_C[pos],false);
	    pos++;
          }
        }
      }
      
      virtual void chain() {
	using Eigen::Matrix;
        using Eigen::Map;
	Eigen::Matrix<double,R1,C1> adjA(_M,_M);
        Eigen::Matrix<double,R2,C2> adjB(_M,_N);
        Eigen::Matrix<double,R1,C2> adjC(_M,_N);

	size_t pos = 0;
        for (size_type j = 0; j < adjC.cols(); j++)
          for (size_type i = 0; i < adjC.rows(); i++)
            adjC(i,j) = _variRefC[pos++]->adj_;
        
        
	adjB = Map<Matrix<double,R1,C1> >(_A,_M,_M)
	  .transpose().colPivHouseholderQr().solve(adjC);
	adjA.noalias() = -adjB
	  * Map<Matrix<double,R1,C2> >(_C,_M,_N).transpose();
        
	pos = 0;
        for (size_type j = 0; j < adjA.cols(); j++)
          for (size_type i = 0; i < adjA.rows(); i++)
            _variRefA[pos++]->adj_ += adjA(i,j);
        
	pos = 0;
        for (size_type j = 0; j < adjB.cols(); j++)
          for (size_type i = 0; i < adjB.rows(); i++)
            _variRefB[pos++]->adj_ += adjB(i,j);
      }
    };
    
    template <int R1,int C1,int R2,int C2>
    class mdivide_left_dv_vari : public vari {
    public:
      int _M; // A.rows() = A.cols() = B.rows()
      int _N; // B.cols()
      double* _A;
      double* _C;
      vari** _variRefB;
      vari** _variRefC;
      
      mdivide_left_dv_vari(const Eigen::Matrix<double,R1,C1> &A,
                           const Eigen::Matrix<var,R2,C2> &B)
      : vari(0.0),
 	_M(A.rows()),
	_N(B.cols()),
	_A((double*)stan::agrad::memalloc_.alloc(sizeof(double) 
                                                 * A.rows() * A.cols())),
	_C((double*)stan::agrad::memalloc_.alloc(sizeof(double) 
                                                 * B.rows() * B.cols())),
        _variRefB((vari**)stan::agrad::memalloc_.alloc(sizeof(vari*) 
						       * B.rows() * B.cols())),
        _variRefC((vari**)stan::agrad::memalloc_.alloc(sizeof(vari*) 
						       * B.rows() * B.cols()))
      {
	using Eigen::Matrix;
	using Eigen::Map;
	
	size_t pos = 0;
	for (size_type j = 0; j < _M; j++) {
	  for (size_type i = 0; i < _M; i++) {
	    _A[pos++] = A(i,j);
	  }
	}
	
	pos = 0;
	for (size_type j = 0; j < _N; j++) {
          for (size_type i = 0; i < _M; i++) {
            _variRefB[pos] = B(i,j).vi_;
            _C[pos++] = B(i,j).val();
          }
        }
                
	Matrix<double,R1,C2> C(_M,_N);
	C = Map<Matrix<double,R1,C2> >(_C,_M,_N);

	C = Map<Matrix<double,R1,C1> >(_A,_M,_M)
	  .colPivHouseholderQr().solve(C);
	
	pos = 0;
        for (size_type j = 0; j < _N; j++) {
          for (size_type i = 0; i < _M; i++) {
	    _C[pos] = C(i,j);
            _variRefC[pos] = new vari(_C[pos],false);
	    pos++;
          }
        }
      }
      
      virtual void chain() {
	using Eigen::Matrix;
        using Eigen::Map;
        Eigen::Matrix<double,R2,C2> adjB(_M,_N);
        Eigen::Matrix<double,R1,C2> adjC(_M,_N);

	size_t pos = 0;
        for (size_type j = 0; j < adjC.cols(); j++)
          for (size_type i = 0; i < adjC.rows(); i++)
            adjC(i,j) = _variRefC[pos++]->adj_;

        adjB = Map<Matrix<double,R1,C1> >(_A,_M,_M)
	  .transpose().colPivHouseholderQr().solve(adjC);

	pos = 0;
        for (size_type j = 0; j < adjB.cols(); j++)
          for (size_type i = 0; i < adjB.rows(); i++)
            _variRefB[pos++]->adj_ += adjB(i,j);
      }
    };
    
    template <int R1,int C1,int R2,int C2>
    class mdivide_left_vd_vari : public vari {
    public:
      int _M; // A.rows() = A.cols() = B.rows()
      int _N; // B.cols()
      double* _A;
      double* _C;
      vari** _variRefA;
      vari** _variRefC;
      
      mdivide_left_vd_vari(const Eigen::Matrix<var,R1,C1> &A,
                           const Eigen::Matrix<double,R2,C2> &B)
      : vari(0.0),
 	_M(A.rows()),
	_N(B.cols()),
	_A((double*)stan::agrad::memalloc_.alloc(sizeof(double) 
                                                 * A.rows() * A.cols())),
	_C((double*)stan::agrad::memalloc_.alloc(sizeof(double) 
                                                 * B.rows() * B.cols())),
        _variRefA((vari**)stan::agrad::memalloc_.alloc(sizeof(vari*) 
						       * A.rows() * A.cols())),
        _variRefC((vari**)stan::agrad::memalloc_.alloc(sizeof(vari*) 
                                                      * B.rows() * B.cols()))
      {
	using Eigen::Matrix;
        using Eigen::Map;

	size_t pos = 0;
	for (size_type j = 0; j < _M; j++) {
          for (size_type i = 0; i < _M; i++) {
            _variRefA[pos] = A(i,j).vi_;
	    _A[pos++] = A(i,j).val();
          }
        }
	
        Matrix<double,R1,C2> C(_M,_N);
	C = Map<Matrix<double,R1,C1> >(_A,_M,_M)
	  .colPivHouseholderQr().solve(B);

	pos = 0;
        for (size_type j = 0; j < _N; j++) {
          for (size_type i = 0; i < _M; i++) {
	    _C[pos] = C(i,j);
            _variRefC[pos] = new vari(_C[pos],false);
	    pos++;
          }
        }
      }
      
      virtual void chain() {
	using Eigen::Matrix;
        using Eigen::Map;
	Eigen::Matrix<double,R1,C1> adjA(_M,_M);
        Eigen::Matrix<double,R1,C2> adjC(_M,_N);

	size_t pos = 0;
        for (size_type j = 0; j < adjC.cols(); j++)
          for (size_type i = 0; i < adjC.rows(); i++)
            adjC(i,j) = _variRefC[pos++]->adj_;
        
	adjA = -Map<Matrix<double,R1,C1> >(_A,_M,_M)
	  .transpose()
	  .colPivHouseholderQr()
	  .solve(adjC*Map<Matrix<double,R1,C2> >(_C,_M,_N).transpose());

	pos = 0;
        for (size_type j = 0; j < adjA.cols(); j++)
          for (size_type i = 0; i < adjA.rows(); i++)
            _variRefA[pos++]->adj_ += adjA(i,j);
      }
    };

    template <int R1,int C1,int R2,int C2>
    inline 
    Eigen::Matrix<var,R1,C2>
    mdivide_left(const Eigen::Matrix<var,R1,C1> &A,
                 const Eigen::Matrix<var,R2,C2> &b) {
      Eigen::Matrix<var,R1,C2> res(b.rows(),b.cols());
      
      stan::math::validate_square(A,"mdivide_left");
      stan::math::validate_multiplicable(A,b,"mdivide_left");
      
      // NOTE: this is not a memory leak, this vari is used in the 
      // expression graph to evaluate the adjoint, but is not needed
      // for the returned matrix.  Memory will be cleaned up with the arena allocator.
      mdivide_left_vv_vari<R1,C1,R2,C2> *baseVari = new mdivide_left_vv_vari<R1,C1,R2,C2>(A,b);
      
      size_t pos = 0;
      for (size_type j = 0; j < res.cols(); j++)
	for (size_type i = 0; i < res.rows(); i++)
          res(i,j).vi_ = baseVari->_variRefC[pos++];
      
      return res;
    }

    template <int R1,int C1,int R2,int C2>
    inline 
    Eigen::Matrix<var,R1,C2>
    mdivide_left(const Eigen::Matrix<var,R1,C1> &A,
                 const Eigen::Matrix<double,R2,C2> &b) {
      Eigen::Matrix<var,R1,C2> res(b.rows(),b.cols());
      
      stan::math::validate_square(A,"mdivide_left");
      stan::math::validate_multiplicable(A,b,"mdivide_left");
      
      // NOTE: this is not a memory leak, this vari is used in the 
      // expression graph to evaluate the adjoint, but is not needed
      // for the returned matrix.  Memory will be cleaned up with the arena allocator.
      mdivide_left_vd_vari<R1,C1,R2,C2> *baseVari = new mdivide_left_vd_vari<R1,C1,R2,C2>(A,b);
      
      size_t pos = 0;
      for (size_type j = 0; j < res.cols(); j++)
	for (size_type i = 0; i < res.rows(); i++)
          res(i,j).vi_ = baseVari->_variRefC[pos++];
      
      return res;
    }
    
    template <int R1,int C1,int R2,int C2>
    inline 
    Eigen::Matrix<var,R1,C2>
    mdivide_left(const Eigen::Matrix<double,R1,C1> &A,
                 const Eigen::Matrix<var,R2,C2> &b) {
      Eigen::Matrix<var,R1,C2> res(b.rows(),b.cols());
      
      stan::math::validate_square(A,"mdivide_left");
      stan::math::validate_multiplicable(A,b,"mdivide_left");
      
      // NOTE: this is not a memory leak, this vari is used in the 
      // expression graph to evaluate the adjoint, but is not needed
      // for the returned matrix.  Memory will be cleaned up with the arena allocator.
{'fraction_non_alphanumeric': 0.141948310139165, 'fraction_numerical': 0.01656726308813784, 'mean_word_length': 2.0255176436278797, 'pattern_counts': {'":': 0, '<': 548, '<?xml version=': 0, '>': 488, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 42, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4860918', 'n_tokens_mistral': 29015, 'n_tokens_neox': 28106, 'n_words': 6554}