Stan  1.0
probability, sampling & optimization
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros Pages
matrix.hpp
Go to the documentation of this file.
1 #ifndef __STAN__AGRAD__MATRIX_HPP__
2 #define __STAN__AGRAD__MATRIX_HPP__
3 
4 // global include
5 #include <sstream>
6 
7 #include <stan/agrad/agrad.hpp>
8 #include <stan/agrad/special_functions.hpp>
9 #include <stan/math/boost_error_handling.hpp>
10 #include <stan/math/matrix.hpp>
11 #include <stan/math/matrix_error_handling.hpp>
12 
13 namespace stan {
14  namespace agrad {
15  class gevv_vvv_vari : public stan::agrad::vari {
16  protected:
17  stan::agrad::vari* alpha_;
18  stan::agrad::vari** v1_;
19  stan::agrad::vari** v2_;
20  double dotval_;
21  size_t length_;
22  inline static double eval_gevv(const stan::agrad::var* alpha,
23  const stan::agrad::var* v1, int stride1,
24  const stan::agrad::var* v2, int stride2,
25  size_t length, double *dotprod) {
26  double result = 0;
27  for (size_t i = 0; i < length; i++)
28  result += v1[i*stride1].vi_->val_ * v2[i*stride2].vi_->val_;
29  *dotprod = result;
30  return alpha->vi_->val_ * result;
31  }
32  public:
33  gevv_vvv_vari(const stan::agrad::var* alpha,
34  const stan::agrad::var* v1, int stride1,
35  const stan::agrad::var* v2, int stride2, size_t length) :
36  vari(eval_gevv(alpha,v1,stride1,v2,stride2,length,&dotval_)), length_(length) {
37  alpha_ = alpha->vi_;
38  v1_ = (stan::agrad::vari**)stan::agrad::memalloc_.alloc(2*length_*sizeof(stan::agrad::vari*));
39  v2_ = v1_ + length_;
40  for (size_t i = 0; i < length_; i++)
41  v1_[i] = v1[i*stride1].vi_;
42  for (size_t i = 0; i < length_; i++)
43  v2_[i] = v2[i*stride2].vi_;
44  }
45  void chain() {
46  const double adj_alpha = adj_ * alpha_->val_;
47  for (size_t i = 0; i < length_; i++) {
48  v1_[i]->adj_ += adj_alpha * v2_[i]->val_;
49  v2_[i]->adj_ += adj_alpha * v1_[i]->val_;
50  }
51  alpha_->adj_ += adj_ * dotval_;
52  }
53  };
54  }
55 }
56 
57 namespace Eigen {
58 
63  template <> struct NumTraits<stan::agrad::var>
64  {
70  typedef stan::agrad::var Real;
71 
77  typedef stan::agrad::var NonInteger;
78 
84  typedef stan::agrad::var Nested;
85 
92  inline static Real epsilon() {
93  return std::numeric_limits<double>::epsilon();
94  }
95 
99  inline static Real dummy_precision() {
100  return 1e-12; // copied from NumTraits.h values for double
101  }
102 
109  inline static Real highest() {
110  return std::numeric_limits<double>::max();
111  }
112 
119  inline static Real lowest() {
120  return -std::numeric_limits<double>::max();
121  }
122 
127  enum {
128  IsInteger = 0,
129  IsSigned = 1,
130  IsComplex = 0,
131  RequireInitialization = 0,
132  ReadCost = 1,
133  AddCost = 1,
134  MulCost = 1,
135  HasFloatingPoint = 1,
136  };
137  };
138 
139  namespace internal {
143  template<>
144  struct significant_decimals_default_impl<stan::agrad::var,false>
145  {
146  static inline int run()
147  {
148  using std::ceil;
149  return cast<double,int>(ceil(-log(NumTraits<stan::agrad::var>::epsilon().val())
150  /log(10.0)));
151  }
152  };
153 
158  template <>
159  struct scalar_product_traits<stan::agrad::var,double> {
160  typedef stan::agrad::var ReturnType;
161  };
162 
167  template <>
168  struct scalar_product_traits<double,stan::agrad::var> {
169  typedef stan::agrad::var ReturnType;
170  };
171 
175  template<typename Index, bool ConjugateLhs, bool ConjugateRhs>
176  struct general_matrix_vector_product<Index,stan::agrad::var,ColMajor,ConjugateLhs,stan::agrad::var,ConjugateRhs>
177  {
178  typedef stan::agrad::var LhsScalar;
179  typedef stan::agrad::var RhsScalar;
180  typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar;
181  enum { LhsStorageOrder = ColMajor };
182 
183  EIGEN_DONT_INLINE static void run(
184  Index rows, Index cols,
185  const LhsScalar* lhs, Index lhsStride,
186  const RhsScalar* rhs, Index rhsIncr,
187  ResScalar* res, Index resIncr, const ResScalar &alpha)
188  {
189  for (Index i = 0; i < rows; i++) {
190  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));
191  }
192  }
193  };
194  template<typename Index, bool ConjugateLhs, bool ConjugateRhs>
195  struct general_matrix_vector_product<Index,stan::agrad::var,RowMajor,ConjugateLhs,stan::agrad::var,ConjugateRhs>
196  {
197  typedef stan::agrad::var LhsScalar;
198  typedef stan::agrad::var RhsScalar;
199  typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar;
200  enum { LhsStorageOrder = RowMajor };
201 
202  EIGEN_DONT_INLINE static void run(
203  Index rows, Index cols,
204  const LhsScalar* lhs, Index lhsStride,
205  const RhsScalar* rhs, Index rhsIncr,
206  ResScalar* res, Index resIncr, const RhsScalar &alpha)
207  {
208  for (Index i = 0; i < rows; i++) {
209  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));
210  }
211  }
212  };
213  template<typename Index, int LhsStorageOrder, bool ConjugateLhs, int RhsStorageOrder, bool ConjugateRhs>
214  struct general_matrix_matrix_product<Index,stan::agrad::var,LhsStorageOrder,ConjugateLhs,stan::agrad::var,RhsStorageOrder,ConjugateRhs,ColMajor>
215  {
216  typedef stan::agrad::var LhsScalar;
217  typedef stan::agrad::var RhsScalar;
218  typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar;
219  static void run(Index rows, Index cols, Index depth,
220  const LhsScalar* _lhs, Index lhsStride,
221  const RhsScalar* _rhs, Index rhsStride,
222  ResScalar* res, Index resStride,
223  const ResScalar &alpha,
224  level3_blocking<LhsScalar,RhsScalar>& blocking,
225  GemmParallelInfo<Index>* info = 0)
226  {
227  for (Index i = 0; i < cols; i++) {
228  general_matrix_vector_product<Index,LhsScalar,LhsStorageOrder,ConjugateLhs,RhsScalar,ConjugateRhs>::run(
229  rows,depth,_lhs,lhsStride,
230  &_rhs[((int)RhsStorageOrder == (int)ColMajor)?(i*rhsStride):(i)],((int)RhsStorageOrder == (int)ColMajor)?(1):(rhsStride),
231  &res[i*resStride],1,alpha);
232  }
233  }
234  };
235  }
236 }
237 
238 namespace stan {
239 
240  namespace agrad {
241 
246  typedef stan::math::EigenType<var>::matrix matrix_v;
247 
252  typedef stan::math::EigenType<var>::vector vector_v;
253 
258  typedef stan::math::EigenType<var>::row_vector row_vector_v;
259 
264  inline void initialize_variable(var& variable, const var& value) {
265  variable = value;
266  }
267 
272  template <int R, int C>
273  inline void initialize_variable(Eigen::Matrix<var,R,C>& matrix, const var& value) {
274  for (int i = 0; i < matrix.size(); ++i)
275  matrix(i) = value;
276  }
277 
281  template <typename T>
282  inline void initialize_variable(std::vector<T>& variables, const var& value) {
283  for (size_t i = 0; i < variables.size(); ++i)
284  initialize_variable(variables[i],value);
285  }
286 
287 
296  inline var to_var(const double& x) {
297  return var(x);
298  }
308  inline void to_var(const double& x, var& var_x) {
309  var_x = x;
310  }
319  inline var to_var(const var& x) {
320  return x;
321  }
331  inline void to_var(const var& var_in, var& var_out) {
332  var_out = var_in;
333  }
342  inline matrix_v to_var(const stan::math::matrix_d& m) {
343  matrix_v m_v(m.rows(), m.cols());
344  for (int j = 0; j < m.cols(); ++j)
345  for (int i = 0; i < m.rows(); ++i)
346  m_v(i,j) = m(i,j);
347  return m_v;
348  }
349 
359  inline void to_var (const stan::math::matrix_d& m, matrix_v& m_v) {
360  m_v.resize(m.rows(), m.cols());
361  for (int j = 0; j < m.cols(); ++j)
362  for (int i = 0; i < m.rows(); ++i)
363  m_v(i,j) = m(i,j);
364  }
373  inline matrix_v to_var(const matrix_v& m) {
374  return m;
375  }
385  inline void to_var(const matrix_v& m_in,
386  matrix_v& m_out) {
387  m_out = m_in;
388  }
398  inline vector_v to_var(const stan::math::vector_d& v) {
399  vector_v v_v(v.size());
400  for (int i = 0; i < v.size(); ++i)
401  v_v[i] = v[i];
402  return v_v;
403  }
413  inline void to_var(const stan::math::vector_d& v,
414  vector_v& v_v) {
415  v_v.resize(v.size());
416  for (int i = 0; i < v.size(); ++i)
417  v_v[i] = v[i];
418  }
428  inline vector_v to_var(const vector_v& v) {
429  return v;
430  }
440  inline void to_var(const vector_v& v_in,
441  vector_v& v_out) {
442  v_out = v_in;
443  }
453  inline row_vector_v to_var(const stan::math::row_vector_d& rv) {
454  row_vector_v rv_v(rv.size());
455  for (int i = 0; i < rv.size(); ++i)
456  rv_v[i] = rv[i];
457  return rv_v;
458  }
468  inline void to_var(const stan::math::row_vector_d& rv,
469  row_vector_v& rv_v) {
470  rv_v.resize(rv.size());
471  for (int i = 0; i < rv.size(); ++i)
472  rv_v[i] = rv[i];
473  }
483  inline row_vector_v to_var(const row_vector_v& rv) {
484  return rv;
485  }
495  inline void to_var(const row_vector_v& rv_in,
496  row_vector_v& rv_out) {
497  rv_out = rv_in;
498  }
499 
509  inline size_t rows(const vector_v& v) {
510  return v.size();
511  }
521  inline size_t rows(const row_vector_v& rv) {
522  return 1;
523  }
532  inline size_t rows(const matrix_v& m) {
533  return m.rows();
534  }
535 
545  inline size_t cols(const vector_v& v) {
546  return 1;
547  }
557  inline size_t cols(const row_vector_v& rv) {
558  return rv.size();
559  }
568  inline size_t cols(const matrix_v& m) {
569  return m.cols();
570  }
571 
572  // scalar returns
573 
584  var determinant(const Eigen::Matrix<var, Eigen::Dynamic, Eigen::Dynamic>& m);
585 
586  namespace {
587  class dot_self_vari : public vari {
588  protected:
589  vari** v_;
590  size_t size_;
591  public:
592  dot_self_vari(vari** v, size_t size)
593  : vari(var_dot_self(v,size)),
594  v_(v),
595  size_(size) {
596  }
597  template <int R, int C>
598  dot_self_vari(const Eigen::Matrix<var,R,C>& v) :
599  vari(var_dot_self(v)), size_(v.size()) {
600  v_ = (vari**) memalloc_.alloc(size_ * sizeof(vari*));
601  for (size_t i = 0; i < size_; ++i)
602  v_[i] = v(i).vi_;
603  }
604  inline static double square(double x) { return x * x; }
605  inline static double var_dot_self(vari** v, size_t size) {
606  double sum = 0.0;
607  for (size_t i = 0; i < size; ++i)
608  sum += square(v[i]->val_);
609  return sum;
610  }
611  template <int R, int C>
612  inline static double var_dot_self(const Eigen::Matrix<var,R,C> &v) {
613  double sum = 0.0;
614  for (int i = 0; i < v.size(); ++i)
615  sum += square(v(i).vi_->val_);
616  return sum;
617  }
618  void chain() {
619  for (size_t i = 0; i < size_; ++i)
620  v_[i]->adj_ += adj_ * 2.0 * v_[i]->val_;
621  }
622  };
623 
624  class sum_v_vari : public vari{
625  protected:
626  vari** v_;
627  size_t length_;
628  inline static double var_sum(const var *v, size_t length) {
629  double result = 0;
630  for (size_t i = 0; i < length; i++)
631  result += v[i].vi_->val_;
632  return result;
633  }
634  template<typename Derived>
635  inline static double var_sum(const Eigen::DenseBase<Derived> &v) {
636  double result = 0;
637  for (int i = 0; i < v.size(); i++)
638  result += v(i).vi_->val_;
639  return result;
640  }
641  public:
642  template<typename Derived>
643  sum_v_vari(const Eigen::DenseBase<Derived> &v) :
644  vari(var_sum(v)), length_(v.size()) {
645  v_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
646  for (size_t i = 0; i < length_; i++)
647  v_[i] = v(i).vi_;
648  }
649  template<int R1,int C1>
650  sum_v_vari(const Eigen::Matrix<var,R1,C1> &v1) :
651  vari(var_sum(v1)), length_(v1.size()) {
652  v_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
653  for (size_t i = 0; i < length_; i++)
654  v_[i] = v1(i).vi_;
655  }
656  sum_v_vari(const var *v, size_t len) :
657  vari(var_sum(v,len)), length_(len) {
658  v_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
659  for (size_t i = 0; i < length_; i++)
660  v_[i] = v[i].vi_;
661  }
662  void chain() {
663  for (size_t i = 0; i < length_; i++) {
664  v_[i]->adj_ += adj_;
665  }
666  }
667  };
668  class dot_product_vv_vari : public vari {
669  protected:
670  vari** v1_;
671  vari** v2_;
672  size_t length_;
673  inline static double var_dot(const var* v1, const var* v2,
674  size_t length) {
675  double result = 0;
676  for (size_t i = 0; i < length; i++)
677  result += v1[i].vi_->val_ * v2[i].vi_->val_;
678  return result;
679  }
680  template<typename Derived1,typename Derived2>
681  inline static double var_dot(const Eigen::DenseBase<Derived1> &v1,
682  const Eigen::DenseBase<Derived2> &v2) {
683  double result = 0;
684  for (int i = 0; i < v1.size(); i++)
685  result += v1[i].vi_->val_ * v2[i].vi_->val_;
686  return result;
687  }
688  inline static double var_dot(vari** v1, vari** v2, size_t length) {
689  double result = 0;
690  for (size_t i = 0; i < length; ++i)
691  result += v1[i]->val_ * v2[i]->val_;
692  return result;
693  }
694  public:
695  dot_product_vv_vari(vari** v1, vari** v2, size_t length)
696  : vari(var_dot(v1,v2,length)),
697  v1_(v1),
698  v2_(v2),
699  length_(length) {
700 
701  }
702  dot_product_vv_vari(const var* v1, const var* v2, size_t length,
703  dot_product_vv_vari* shared_v1 = NULL,
704  dot_product_vv_vari* shared_v2 = NULL) :
705  vari(var_dot(v1, v2, length)), length_(length) {
706  if (shared_v1 == NULL) {
707  v1_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
708  for (size_t i = 0; i < length_; i++)
709  v1_[i] = v1[i].vi_;
710  }
711  else {
712  v1_ = shared_v1->v1_;
713  }
714  if (shared_v2 == NULL) {
715  v2_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
716  for (size_t i = 0; i < length_; i++)
717  v2_[i] = v2[i].vi_;
718  }
719  else {
720  v2_ = shared_v2->v2_;
721  }
722  }
723  template<typename Derived1,typename Derived2>
724  dot_product_vv_vari(const Eigen::DenseBase<Derived1> &v1,
725  const Eigen::DenseBase<Derived2> &v2,
726  dot_product_vv_vari* shared_v1 = NULL,
727  dot_product_vv_vari* shared_v2 = NULL) :
728  vari(var_dot(v1, v2)), length_(v1.size()) {
729  if (shared_v1 == NULL) {
730  v1_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
731  for (size_t i = 0; i < length_; i++)
732  v1_[i] = v1[i].vi_;
733  }
734  else {
735  v1_ = shared_v1->v1_;
736  }
737  if (shared_v2 == NULL) {
738  v2_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
739  for (size_t i = 0; i < length_; i++)
740  v2_[i] = v2[i].vi_;
741  }
742  else {
743  v2_ = shared_v2->v2_;
744  }
745  }
746  template<int R1,int C1,int R2,int C2>
747  dot_product_vv_vari(const Eigen::Matrix<var,R1,C1> &v1,
748  const Eigen::Matrix<var,R2,C2> &v2,
749  dot_product_vv_vari* shared_v1 = NULL,
750  dot_product_vv_vari* shared_v2 = NULL) :
751  vari(var_dot(v1, v2)), length_(v1.size()) {
752  if (shared_v1 == NULL) {
753  v1_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
754  for (size_t i = 0; i < length_; i++)
755  v1_[i] = v1[i].vi_;
756  }
757  else {
758  v1_ = shared_v1->v1_;
759  }
760  if (shared_v2 == NULL) {
761  v2_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
762  for (size_t i = 0; i < length_; i++)
763  v2_[i] = v2[i].vi_;
764  }
765  else {
766  v2_ = shared_v2->v2_;
767  }
768  }
769  void chain() {
770  for (size_t i = 0; i < length_; i++) {
771  v1_[i]->adj_ += adj_ * v2_[i]->val_;
772  v2_[i]->adj_ += adj_ * v1_[i]->val_;
773  }
774  }
775  };
776 
777  class dot_product_vd_vari : public vari {
778  protected:
779  vari** v1_;
780  double* v2_;
781  size_t length_;
782  inline static double var_dot(const var* v1, const double* v2,
783  size_t length) {
784  double result = 0;
785  for (size_t i = 0; i < length; i++)
786  result += v1[i].vi_->val_ * v2[i];
787  return result;
788  }
789  template<typename Derived1,typename Derived2>
790  inline static double var_dot(const Eigen::DenseBase<Derived1> &v1,
791  const Eigen::DenseBase<Derived2> &v2) {
792  double result = 0;
793  for (int i = 0; i < v1.size(); i++)
794  result += v1[i].vi_->val_ * v2[i];
795  return result;
796  }
797  public:
798  dot_product_vd_vari(const var* v1, const double* v2, size_t length,
799  dot_product_vd_vari *shared_v1 = NULL,
800  dot_product_vd_vari *shared_v2 = NULL) :
801  vari(var_dot(v1, v2, length)), length_(length) {
802  if (shared_v1 == NULL) {
803  v1_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
804  for (size_t i = 0; i < length_; i++)
805  v1_[i] = v1[i].vi_;
806  } else {
807  v1_ = shared_v1->v1_;
808  }
809  if (shared_v2 == NULL) {
810  v2_ = (double*)memalloc_.alloc(length_*sizeof(double));
811  for (size_t i = 0; i < length_; i++)
812  v2_[i] = v2[i];
813  } else {
814  v2_ = shared_v2->v2_;
815  }
816  }
817  template<typename Derived1,typename Derived2>
818  dot_product_vd_vari(const Eigen::DenseBase<Derived1> &v1,
819  const Eigen::DenseBase<Derived2> &v2,
820  dot_product_vd_vari *shared_v1 = NULL,
821  dot_product_vd_vari *shared_v2 = NULL) :
822  vari(var_dot(v1, v2)), length_(v1.size()) {
823  if (shared_v1 == NULL) {
824  v1_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
825  for (size_t i = 0; i < length_; i++)
826  v1_[i] = v1[i].vi_;
827  } else {
828  v1_ = shared_v1->v1_;
829  }
830  if (shared_v2 == NULL) {
831  v2_ = (double*)memalloc_.alloc(length_*sizeof(double));
832  for (size_t i = 0; i < length_; i++)
833  v2_[i] = v2[i];
834  } else {
835  v2_ = shared_v2->v2_;
836  }
837  }
838  template<int R1,int C1,int R2,int C2>
839  dot_product_vd_vari(const Eigen::Matrix<var,R1,C1> &v1,
840  const Eigen::Matrix<double,R2,C2> &v2,
841  dot_product_vd_vari *shared_v1 = NULL,
842  dot_product_vd_vari *shared_v2 = NULL) :
843  vari(var_dot(v1, v2)), length_(v1.size()) {
844  if (shared_v1 == NULL) {
845  v1_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
846  for (size_t i = 0; i < length_; i++)
847  v1_[i] = v1[i].vi_;
848  } else {
849  v1_ = shared_v1->v1_;
850  }
851  if (shared_v2 == NULL) {
852  v2_ = (double*)memalloc_.alloc(length_*sizeof(double));
853  for (size_t i = 0; i < length_; i++)
854  v2_[i] = v2[i];
855  } else {
856  v2_ = shared_v2->v2_;
857  }
858  }
859  void chain() {
860  for (size_t i = 0; i < length_; i++) {
861  v1_[i]->adj_ += adj_ * v2_[i];
862  }
863  }
864  };
865 
866  }
867 
878  template<int R, int C>
879  inline var dot_self(const Eigen::Matrix<var, R, C>& v) {
880  stan::math::validate_vector(v,"dot_self");
881  return var(new dot_self_vari(v));
882  }
883 
892  template<int R1,int C1,int R2, int C2>
893  inline var dot_product(const Eigen::Matrix<var, R1, C1>& v1,
894  const Eigen::Matrix<var, R2, C2>& v2) {
895  stan::math::validate_vector(v1,"dot_product");
896  stan::math::validate_vector(v2,"dot_product");
897  stan::math::validate_matching_sizes(v1,v2,"dot_product");
898  return var(new dot_product_vv_vari(v1,v2));
899  }
909  template<int R1,int C1,int R2, int C2>
910  inline var dot_product(const Eigen::Matrix<var, R1, C1>& v1,
911  const Eigen::Matrix<double, R2, C2>& v2) {
912  stan::math::validate_vector(v1,"dot_product");
913  stan::math::validate_vector(v2,"dot_product");
914  stan::math::validate_matching_sizes(v1,v2,"dot_product");
915  return var(new dot_product_vd_vari(v1,v2));
916  }
926  template<int R1,int C1,int R2, int C2>
927  inline var dot_product(const Eigen::Matrix<double, R1, C1>& v1,
928  const Eigen::Matrix<var, R2, C2>& v2) {
929  stan::math::validate_vector(v1,"dot_product");
930  stan::math::validate_vector(v2,"dot_product");
931  stan::math::validate_matching_sizes(v1,v2,"dot_product");
932  return var(new dot_product_vd_vari(v2,v1));
933  }
942  inline var dot_product(const var* v1, const var* v2, size_t length) {
943  return var(new dot_product_vv_vari(v1, v2, length));
944  }
953  inline var dot_product(const var* v1, const double* v2, size_t length) {
954  return var(new dot_product_vd_vari(v1, v2, length));
955  }
964  inline var dot_product(const double* v1, const var* v2, size_t length) {
965  return var(new dot_product_vd_vari(v2, v1, length));
966  }
975  inline var dot_product(const std::vector<var>& v1,
976  const std::vector<var>& v2) {
977  stan::math::validate_matching_sizes(v1,v2,"dot_product");
978  return var(new dot_product_vv_vari(&v1[0], &v2[0], v1.size()));
979  }
988  inline var dot_product(const std::vector<var>& v1,
989  const std::vector<double>& v2) {
990  stan::math::validate_matching_sizes(v1,v2,"dot_product");
991  return var(new dot_product_vd_vari(&v1[0], &v2[0], v1.size()));
992  }
1001  inline var dot_product(const std::vector<double>& v1,
1002  const std::vector<var>& v2) {
1003  stan::math::validate_matching_sizes(v1,v2,"dot_product");
1004  return var(new dot_product_vd_vari(&v2[0], &v1[0], v1.size()));
1005  }
1006 
1013  template <int R, int C>
1014  inline var sum(const Eigen::Matrix<var,R,C>& m) {
1015  if (m.size() == 0)
1016  return 0.0;
1017  return var(new sum_v_vari(m));
1018  }
1019 
1020 
1021 
1031  template <typename T1, typename T2, int R, int C>
1032  inline Eigen::Matrix<typename boost::math::tools::promote_args<T1,T2>::type, R, C>
1033  add(const Eigen::Matrix<T1, R, C>& v1,
1034  const Eigen::Matrix<T2, R, C>& v2) {
1035  stan::math::validate_matching_dims(v1,v2,"add");
1036  return to_var(v1) + to_var(v2);
1037  }
1044  template <typename T1, typename T2, int R, int C>
1045  inline Eigen::Matrix<typename boost::math::tools::promote_args<T1,T2>::type, R, C>
1046  add(const Eigen::Matrix<T1, R, C>& m,
1047  const T2& c) {
1048  return (to_var(m).array() + to_var(c)).matrix();
1049  }
1056  template <typename T1, typename T2, int R, int C>
1057  inline Eigen::Matrix<typename boost::math::tools::promote_args<T1,T2>::type, R, C>
1058  add(const T1& c,
1059  const Eigen::Matrix<T2, R, C>& m) {
1060  return (to_var(c) + to_var(m).array()).matrix();
1061  }
1062 
1063 
1072  template <typename T1, typename T2, int R, int C>
1073  inline Eigen::Matrix<typename boost::math::tools::promote_args<T1,T2>::type, R, C>
1074  subtract(const Eigen::Matrix<T1, R, C>& v1,
1075  const Eigen::Matrix<T2, R, C>& v2) {
1076  stan::math::validate_matching_dims(v1,v2,"subtract");
1077  return to_var(v1) - to_var(v2);
1078  }
1085  template <typename T1, typename T2, int R, int C>
1086  inline Eigen::Matrix<typename boost::math::tools::promote_args<T1,T2>::type, R, C>
1087  subtract(const Eigen::Matrix<T1, R, C>& m,
1088  const T2& c) {
1089  return (to_var(m).array() - to_var(c)).matrix();
1090  }
1097  template <typename T1, typename T2, int R, int C>
1098  inline Eigen::Matrix<typename boost::math::tools::promote_args<T1,T2>::type, R, C>
1099  subtract(const T1& c,
1100  const Eigen::Matrix<T2, R, C>& m) {
1101  return (to_var(c) - to_var(m).array()).matrix();
1102  }
1103 
1110  template <typename T>
1111  inline T minus(const T& v) {
1112  return -v;
1113  }
1120  template <typename T, int R, int C>
1121  inline Eigen::Matrix<T,R,C> minus(const Eigen::Matrix<T, R, C>& v) {
1122  return -v;
1123  }
1124 
1132  inline double
1133  divide(double x, double y) {
1134  return x / y;
1135  }
1136  template <typename T1, typename T2>
1137  inline var
1138  divide(const T1& v, const T2& c) {
1139  return to_var(v) / to_var(c);
1140  }
1148  template <typename T1, typename T2, int R, int C>
1149  inline Eigen::Matrix<var,R,C>
1150  divide(const Eigen::Matrix<T1, R,C>& v, const T2& c) {
1151  return to_var(v) / to_var(c);
1152  }
1153 
1160  template <typename T1, typename T2>
1161  inline Eigen::Matrix<stan::agrad::var,Eigen::Dynamic,Eigen::Dynamic>
1162  elt_multiply(const Eigen::Matrix<T1,Eigen::Dynamic,Eigen::Dynamic>& m1,
1163  const Eigen::Matrix<T2,Eigen::Dynamic,Eigen::Dynamic>& m2) {
1164  stan::math::validate_matching_dims(m1,m2,"elt_multiply");
1165  Eigen::Matrix<stan::agrad::var,Eigen::Dynamic, Eigen::Dynamic>
1166  result(m1.rows(),m1.cols());
1167  for (int j = 0; j < m1.cols(); ++j)
1168  for (int i = 0; i < m1.rows(); ++i)
1169  result(i,j) = m1(i,j) * m2(i,j);
1170  return result;
1171  }
1178  template <typename T1, typename T2>
1179  inline Eigen::Matrix<stan::agrad::var,Eigen::Dynamic,1>
1180  elt_multiply(const Eigen::Matrix<T1,Eigen::Dynamic,1>& v1,
1181  const Eigen::Matrix<T2,Eigen::Dynamic,1>& v2) {
1182  stan::math::validate_matching_dims(v1,v2,"elt_multiply");
1183  Eigen::Matrix<stan::agrad::var,Eigen::Dynamic,1> result(v1.size());
1184  for (int i = 0; i < v1.size(); ++i)
1185  result(i) = v1(i) * v2(i);
1186  return result;
1187  }
1194  template <typename T1, typename T2>
1195  inline Eigen::Matrix<stan::agrad::var,1,Eigen::Dynamic>
1196  elt_multiply(const Eigen::Matrix<T1,1,Eigen::Dynamic>& v1,
1197  const Eigen::Matrix<T2,1,Eigen::Dynamic>& v2) {
1198  stan::math::validate_matching_dims(v1,v2,"elt_multiply");
1199  Eigen::Matrix<stan::agrad::var,1,Eigen::Dynamic> result(v1.size());
1200  for (int i = 0; i < v1.size(); ++i)
1201  result(i) = v1(i) * v2(i);
1202  return result;
1203  }
1204 
1205 
1212  template <typename T1, typename T2>
1213  inline Eigen::Matrix<stan::agrad::var,Eigen::Dynamic,Eigen::Dynamic>
1214  elt_divide(const Eigen::Matrix<T1,Eigen::Dynamic,Eigen::Dynamic>& m1,
1215  const Eigen::Matrix<T2,Eigen::Dynamic,Eigen::Dynamic>& m2) {
1216  stan::math::validate_matching_dims(m1,m2,"elt_divide");
1217  Eigen::Matrix<stan::agrad::var,Eigen::Dynamic, Eigen::Dynamic>
1218  result(m1.rows(),m1.cols());
1219  for (int j = 0; j < m1.cols(); ++j)
1220  for (int i = 0; i < m1.rows(); ++i)
1221  result(i,j) = m1(i,j) / m2(i,j);
1222  return result;
1223  }
1230  template <typename T1, typename T2>
1231  inline Eigen::Matrix<stan::agrad::var,Eigen::Dynamic,1>
1232  elt_divide(const Eigen::Matrix<T1,Eigen::Dynamic,1>& v1,
1233  const Eigen::Matrix<T2,Eigen::Dynamic,1>& v2) {
1234  stan::math::validate_matching_dims(v1,v2,"elt_divide");
1235  Eigen::Matrix<stan::agrad::var,Eigen::Dynamic,1> result(v1.size());
1236  for (int i = 0; i < v1.size(); ++i)
1237  result(i) = v1(i) / v2(i);
1238  return result;
1239  }
1246  template <typename T1, typename T2>
1247  inline Eigen::Matrix<stan::agrad::var,1,Eigen::Dynamic>
1248  elt_divide(const Eigen::Matrix<T1,1,Eigen::Dynamic>& v1,
1249  const Eigen::Matrix<T2,1,Eigen::Dynamic>& v2) {
1250  stan::math::validate_matching_dims(v1,v2,"elt_divide");
1251  Eigen::Matrix<stan::agrad::var,1,Eigen::Dynamic> result(v1.size());
1252  for (int i = 0; i < v1.size(); ++i)
1253  result(i) = v1(i) / v2(i);
1254  return result;
1255  }
1256 
1263  template <typename T1, typename T2>
1264  inline
1265  typename boost::math::tools::promote_args<T1,T2>::type
1266  multiply(const T1& v, const T2& c) {
1267  return v * c;
1268  }
1269 
1276  template<typename T1,typename T2,int R2,int C2>
1277  inline Eigen::Matrix<var,R2,C2> multiply(const T1& c,
1278  const Eigen::Matrix<T2, R2, C2>& m) {
1279  // FIXME: pull out to eliminate overpromotion of one side
1280  // move to matrix.hpp w. promotion?
1281  return to_var(m) * to_var(c);
1282  }
1283 
1290  template<typename T1,int R1,int C1,typename T2>
1291  inline Eigen::Matrix<var,R1,C1> multiply(const Eigen::Matrix<T1, R1, C1>& m,
1292  const T2& c) {
1293  return to_var(m) * to_var(c);
1294  }
1295 
1306  template<int R1,int C1,int R2,int C2>
1307  inline Eigen::Matrix<var,R1,C2> multiply(const Eigen::Matrix<var,R1,C1>& m1,
1308  const Eigen::Matrix<var,R2,C2>& m2) {
1309  stan::math::validate_multiplicable(m1,m2,"multiply");
1310  Eigen::Matrix<var,R1,C2> result(m1.rows(),m2.cols());
1311  for (int i = 0; i < m1.rows(); i++) {
1312  typename Eigen::Matrix<var,R1,C1>::ConstRowXpr crow(m1.row(i));
1313  for (int j = 0; j < m2.cols(); j++) {
1314  typename Eigen::Matrix<var,R2,C2>::ConstColXpr ccol(m2.col(j));
1315  if (j == 0) {
1316  if (i == 0) {
1317  result(i,j) = var(new dot_product_vv_vari(crow,ccol));
1318  }
1319  else {
1320  dot_product_vv_vari *v2 = static_cast<dot_product_vv_vari*>(result(0,j).vi_);
1321  result(i,j) = var(new dot_product_vv_vari(crow,ccol,NULL,v2));
1322  }
1323  }
1324  else {
1325  if (i == 0) {
1326  dot_product_vv_vari *v1 = static_cast<dot_product_vv_vari*>(result(i,0).vi_);
1327  result(i,j) = var(new dot_product_vv_vari(crow,ccol,v1));
1328  }
1329  else /* if (i != 0 && j != 0) */ {
1330  dot_product_vv_vari *v1 = static_cast<dot_product_vv_vari*>(result(i,0).vi_);
1331  dot_product_vv_vari *v2 = static_cast<dot_product_vv_vari*>(result(0,j).vi_);
1332  result(i,j) = var(new dot_product_vv_vari(crow,ccol,v1,v2));
1333  }
1334  }
1335  }
1336  }
1337  return result;
1338  }
1339 
1350  template<int R1,int C1,int R2,int C2>
1351  inline Eigen::Matrix<var,R1,C2> multiply(const Eigen::Matrix<double,R1,C1>& m1,
1352  const Eigen::Matrix<var,R2,C2>& m2) {
1353  stan::math::validate_multiplicable(m1,m2,"multiply");
1354  Eigen::Matrix<var,R1,C2> result(m1.rows(),m2.cols());
1355  for (int i = 0; i < m1.rows(); i++) {
1356  typename Eigen::Matrix<double,R1,C1>::ConstRowXpr crow(m1.row(i));
1357  for (int j = 0; j < m2.cols(); j++) {
1358  typename Eigen::Matrix<var,R2,C2>::ConstColXpr ccol(m2.col(j));
1359 // result(i,j) = dot_product(crow,ccol);
1360  if (j == 0) {
1361  if (i == 0) {
1362  result(i,j) = var(new dot_product_vd_vari(ccol,crow));
1363  }
1364  else {
1365  dot_product_vd_vari *v2 = static_cast<dot_product_vd_vari*>(result(0,j).vi_);
1366  result(i,j) = var(new dot_product_vd_vari(ccol,crow,v2,NULL));
1367  }
1368  }
1369  else {
1370  if (i == 0) {
1371  dot_product_vd_vari *v1 = static_cast<dot_product_vd_vari*>(result(i,0).vi_);
1372  result(i,j) = var(new dot_product_vd_vari(ccol,crow,NULL,v1));
1373  }
1374  else /* if (i != 0 && j != 0) */ {
1375  dot_product_vd_vari *v1 = static_cast<dot_product_vd_vari*>(result(i,0).vi_);
1376  dot_product_vd_vari *v2 = static_cast<dot_product_vd_vari*>(result(0,j).vi_);
1377  result(i,j) = var(new dot_product_vd_vari(ccol,crow,v2,v1));
1378  }
1379  }
1380  }
1381  }
1382  return result;
1383  }
1384 
1395  template<int R1,int C1,int R2,int C2>
1396  inline Eigen::Matrix<var,R1,C2> multiply(const Eigen::Matrix<var,R1,C1>& m1,
1397  const Eigen::Matrix<double,R2,C2>& m2) {
1398  stan::math::validate_multiplicable(m1,m2,"multiply");
1399  Eigen::Matrix<var,R1,C2> result(m1.rows(),m2.cols());
1400  for (int i = 0; i < m1.rows(); i++) {
1401  typename Eigen::Matrix<var,R1,C1>::ConstRowXpr crow(m1.row(i));
1402  for (int j = 0; j < m2.cols(); j++) {
1403  typename Eigen::Matrix<double,R2,C2>::ConstColXpr ccol(m2.col(j));
1404 // result(i,j) = dot_product(crow,ccol);
1405  if (j == 0) {
1406  if (i == 0) {
1407  result(i,j) = var(new dot_product_vd_vari(crow,ccol));
1408  }
1409  else {
1410  dot_product_vd_vari *v2 = static_cast<dot_product_vd_vari*>(result(0,j).vi_);
1411  result(i,j) = var(new dot_product_vd_vari(crow,ccol,NULL,v2));
1412  }
1413  }
1414  else {
1415  if (i == 0) {
1416  dot_product_vd_vari *v1 = static_cast<dot_product_vd_vari*>(result(i,0).vi_);
1417  result(i,j) = var(new dot_product_vd_vari(crow,ccol,v1,NULL));
1418  }
1419  else /* if (i != 0 && j != 0) */ {
1420  dot_product_vd_vari *v1 = static_cast<dot_product_vd_vari*>(result(i,0).vi_);
1421  dot_product_vd_vari *v2 = static_cast<dot_product_vd_vari*>(result(0,j).vi_);
1422  result(i,j) = var(new dot_product_vd_vari(crow,ccol,v1,v2));
1423  }
1424  }
1425  }
1426  }
1427  return result;
1428  }
1429 
1439  template <int C1,int R2>
1440  inline var multiply(const Eigen::Matrix<var, 1, C1>& rv,
1441  const Eigen::Matrix<var, R2, 1>& v) {
1442  if (rv.size() != v.size())
1443  throw std::domain_error("row vector and vector must be same length in multiply");
1444  return dot_product(rv, v);
1445  }
1455  template <int C1,int R2>
1456  inline var multiply(const Eigen::Matrix<double, 1, C1>& rv,
1457  const Eigen::Matrix<var, R2, 1>& v) {
1458  stan::math::validate_multiplicable(rv,v,"multiply");
1459  return dot_product(rv, v);
1460  }
1470  template <int C1,int R2>
1471  inline var multiply(const Eigen::Matrix<var, 1, C1>& rv,
1472  const Eigen::Matrix<double, R2, 1>& v) {
1473  stan::math::validate_multiplicable(rv,v,"multiply");
1474  return dot_product(rv, v);
1475  }
1476 
1477  inline matrix_v
1478  multiply_lower_tri_self_transpose(const matrix_v& L) {
1479  stan::math::validate_square(L,"multiply_lower_tri_self_transpose");
1480  int K = L.rows();
1481  matrix_v LLt(K,K);
1482  if (K == 0) return LLt;
1483  // if (K == 1) {
1484  // LLt(0,0) = L(0,0) * L(0,0);
1485  // return LLt;
1486  // }
1487  int Knz = (K * (K + 1)) / 2; // nonzero: (K choose 2) below
1488  // diag + K on diag
1489  vari** vs = (vari**)memalloc_.alloc( Knz * sizeof(vari*) );
1490  int pos = 0;
1491  for (int m = 0; m < K; ++m)
1492  for (int n = 0; n <= m; ++n)
1493  vs[pos++] = L(m,n).vi_;
1494  for (int m = 0, mpos=0; m < K; ++m, mpos += m) {
1495  // FIXME: replace with dot_self
1496  LLt(m,m) = var(new dot_self_vari(vs + mpos, m + 1));
1497  for (int n = 0, npos = 0; n < m; ++n, npos += n)
1498  LLt(m,n) = LLt(n,m) = var(new dot_product_vv_vari(vs + mpos, vs + npos, n + 1));
1499  }
1500  return LLt;
1501  }
1502 
1509  inline matrix_v
1510  tcrossprod(const matrix_v& M) {
1511  if(M.rows() == 0)
1512  return matrix_v(0,0);
1513  if(M.rows() == 1) {
1514  return M * M.transpose(); // FIXME: replace with dot_self
1515  }
1516  matrix_v result(M.rows(),M.rows());
1517  return result.setZero().selfadjointView<Eigen::Upper>().rankUpdate(M);
1518  }
1519 
1526  inline matrix_v
1527  crossprod(const matrix_v& M) {
1528  return tcrossprod(M.transpose());
1529  }
1530 
1542  row_vector_v row(const matrix_v& m, size_t i);
1543 
1555  vector_v col(const matrix_v& m, size_t j);
1556 
1563  vector_v diagonal(const matrix_v& m);
1564 
1571  matrix_v diag_matrix(const vector_v& v);
1572 
1579  row_vector_v transpose(const vector_v& v);
1586  vector_v transpose(const row_vector_v& rv);
1592  matrix_v transpose(const matrix_v& m);
1593 
1599  matrix_v inverse(const matrix_v& m);
1600 
1606  vector_v softmax(const vector_v& v);
1607 
1608  template<int R1,int C1,int R2,int C2>
1609  inline Eigen::Matrix<var,R1,C2> mdivide_left_tri_low(const Eigen::Matrix<var,R1,C1> &A,
1610  const Eigen::Matrix<var,R2,C2> &b) {
1611  stan::math::validate_square(A,"mdivide_left_tri_low");
1612  stan::math::validate_multiplicable(A,b,"mdivide_left_tri_low");
1613  return A.template triangularView<Eigen::Lower>().solve(b);
1614  }
1615  // FIXME: fold next two into above by templating out scalar types
1616  // & emplying to_var to both variables
1617  template<int R1,int C1,int R2,int C2>
1618  inline Eigen::Matrix<var,R1,C2> mdivide_left_tri_low(const Eigen::Matrix<var,R1,C1> &A,
1619  const Eigen::Matrix<double,R2,C2> &b) {
1620  stan::math::validate_square(A,"mdivide_left_tri_low");
1621  stan::math::validate_multiplicable(A,b,"mdivide_left_tri_low");
1622  return A.template triangularView<Eigen::Lower>().solve(to_var(b));
1623  }
1624  template<int R1,int C1,int R2,int C2>
1625  inline Eigen::Matrix<var,R1,C2> mdivide_left_tri_low(const Eigen::Matrix<double,R1,C1> &A,
1626  const Eigen::Matrix<var,R2,C2> &b) {
1627  stan::math::validate_square(A,"mdivide_left_tri_low");
1628  stan::math::validate_multiplicable(A,b,"mdivide_left_tri_low");
1629  return to_var(A).template triangularView<Eigen::Lower>().solve(b);
1630  }
1631 
1632  inline Eigen::Matrix<var,Eigen::Dynamic,Eigen::Dynamic>
1633  mdivide_left_tri_low(const Eigen::Matrix<var,Eigen::Dynamic,Eigen::Dynamic> &A) {
1634  stan::math::validate_square(A,"mdivide_left_tri_low");
1635  int n = A.rows();
1636  Eigen::Matrix<var,Eigen::Dynamic,Eigen::Dynamic> b;
1637  b.setIdentity(n,n);
1638  A.triangularView<Eigen::Lower>().solveInPlace(b);
1639  return(b);
1640  }
1641 
1651  template<int TriView,int R1,int C1,int R2,int C2>
1652  inline Eigen::Matrix<var,R1,C2> mdivide_left_tri(const Eigen::Matrix<var,R1,C1> &A,
1653  const Eigen::Matrix<var,R2,C2> &b) {
1654  stan::math::validate_square(A,"mdivide_left_tri");
1655  stan::math::validate_multiplicable(A,b,"mdivide_left_tri");
1656  return A.template triangularView<TriView>().solve(b);
1657  }
1667  template<int TriView,int R1,int C1,int R2,int C2>
1668  inline Eigen::Matrix<var,R1,C2> mdivide_left_tri(const Eigen::Matrix<var,R1,C1> &A,
1669  const Eigen::Matrix<double,R2,C2> &b) {
1670  stan::math::validate_square(A,"mdivide_left_tri");
1671  stan::math::validate_multiplicable(A,b,"mdivide_left_tri");
1672  return A.template triangularView<TriView>().solve(to_var(b));
1673  }
1683  template<int TriView,int R1,int C1,int R2,int C2>
1684  inline Eigen::Matrix<var,R1,C2> mdivide_left_tri(const Eigen::Matrix<double,R1,C1> &A,
1685  const Eigen::Matrix<var,R2,C2> &b) {
1686  stan::math::validate_square(A,"mdivide_left_tri");
1687  stan::math::validate_multiplicable(A,b,"mdivide_left_tri");
1688  return to_var(A).template triangularView<TriView>().solve(b);
1689  }
1690 
1698  template<int TriView>
1699  inline Eigen::Matrix<var,Eigen::Dynamic,Eigen::Dynamic>
1700  mdivide_left_tri(const Eigen::Matrix<var,Eigen::Dynamic,Eigen::Dynamic> &A) {
1701  stan::math::validate_square(A,"mdivide_left_tri_low");
1702  int n = A.rows();
1703  Eigen::Matrix<var,Eigen::Dynamic,Eigen::Dynamic> b;
1704  b.setIdentity(n,n);
1705  A.triangularView<TriView>().solveInPlace(b);
1706  return b;
1707  }
1708 
1717  template<int R1,int C1,int R2,int C2>
1718  inline Eigen::Matrix<var,R1,C2> mdivide_left(const Eigen::Matrix<var,R1,C1> &A,
1719  const Eigen::Matrix<var,R2,C2> &b) {
1720  stan::math::validate_square(A,"mdivide_left");
1721  stan::math::validate_multiplicable(A,b,"mdivide_left");
1722  return A.lu().solve(b);
1723  }
1732  template<int R1, int C1, int R2, int C2>
1733  inline Eigen::Matrix<var,R1,C2> mdivide_left(const Eigen::Matrix<double,R1,C1> &A,
1734  const Eigen::Matrix<var,R2,C2> &b) {
1735  stan::math::validate_square(A,"mdivide_left");
1736  stan::math::validate_multiplicable(A,b,"mdivide_left");
1737  // FIXME: it would be much faster to do LU, then convert to var
1738  return to_var(A).lu().solve(b);
1739  }
1748  template<int R1, int C1, int R2, int C2>
1749  inline Eigen::Matrix<var,R1,C2> mdivide_left(const Eigen::Matrix<var,R1,C1> &A,
1750  const Eigen::Matrix<double,R2,C2> &b) {
1751  stan::math::validate_square(A,"mdivide_left");
1752  stan::math::validate_multiplicable(A,b,"mdivide_left");
1753  return A.lu().solve(to_var(b));
1754  }
1755 
1764  template<int R1,int C1,int R2,int C2>
1765  inline Eigen::Matrix<var,R1,C2> mdivide_right(const Eigen::Matrix<var,R1,C1> &b,
1766  const Eigen::Matrix<var,R2,C2> &A) {
1767  stan::math::validate_square(A,"mdivide_right");
1768  stan::math::validate_multiplicable(b,A,"mdivide_right");
1769  return A.transpose().lu().solve(b.transpose()).transpose();
1770  }
1779  template<int R1,int C1,int R2,int C2>
1780  inline Eigen::Matrix<var,R1,C2> mdivide_right(const Eigen::Matrix<double,R1,C1> &b,
1781  const Eigen::Matrix<var,R2,C2> &A) {
1782  stan::math::validate_square(A,"mdivide_right");
1783  stan::math::validate_multiplicable(b,A,"mdivide_right");
1784  return A.transpose().lu().solve(to_var(b).transpose()).transpose();
1785  }
1794  template<int R1,int C1,int R2,int C2>
1795  inline Eigen::Matrix<var,R1,C2> mdivide_right(const Eigen::Matrix<var,R1,C1> &b,
1796  const Eigen::Matrix<double,R2,C2> &A) {
1797  stan::math::validate_square(A,"mdivide_right");
1798  stan::math::validate_multiplicable(b,A,"mdivide_right");
1799  return to_var(A).transpose().lu().solve(b.transpose()).transpose();
1800  }
1801 
1811  template<int TriView,typename T1,typename T2,int R1,int C1,int R2,int C2>
1812  inline Eigen::Matrix<typename boost::math::tools::promote_args<T1,T2>::type,R1,C2>
1813  mdivide_right_tri(const Eigen::Matrix<T1,R1,C1> &b,
1814  const Eigen::Matrix<T2,R2,C2> &A) {
1815  stan::math::validate_square(A,"mdivide_right_tri_low");
1816  stan::math::validate_multiplicable(b,A,"mdivide_right_tri_low");
1817  return to_var(A)
1818  .template triangularView<TriView>()
1819  .transpose()
1820  .solve(to_var(b).transpose())
1821  .transpose();
1822  }
1833  template<typename T1,typename T2,int R1,int C1,int R2,int C2>
1834  inline Eigen::Matrix<typename boost::math::tools::promote_args<T1,T2>::type,R1,C2>
1835  mdivide_right_tri_low(const Eigen::Matrix<T1,R1,C1> &b,
1836  const Eigen::Matrix<T2,R2,C2> &A) {
1837  return mdivide_right_tri<Eigen::Lower,T1,T2,R1,C1,R2,C2>(b,A);
1838  }
1839 
1847  vector_v eigenvalues(const matrix_v& m);
1848 
1856  matrix_v eigenvectors(const matrix_v& m);
1874  void eigen_decompose(const matrix_v& m,
1875  vector_v& eigenvalues,
1876  matrix_v& eigenvectors);
1877 
1887  vector_v eigenvalues_sym(const matrix_v& m);
1897  matrix_v eigenvectors_sym(const matrix_v& m);
1910  void eigen_decompose_sym(const matrix_v& m,
1911  vector_v& eigenvalues,
1912  matrix_v& eigenvectors);
1913 
1924  matrix_v cholesky_decompose(const matrix_v& m);
1925 
1934  vector_v singular_values(const matrix_v& m);
1935 
1957  void svd(const matrix_v& m,
1958  matrix_v& u,
1959  matrix_v& v,
1960  vector_v& s);
1961 
1962 
1963  // FIXME: double val?
1964  inline void assign_to_var(stan::agrad::var& var, const double& val) {
1965  var = val;
1966  }
1967  inline void assign_to_var(stan::agrad::var& var, const stan::agrad::var& val) {
1968  var = val;
1969  }
1970  // FIXME: int val?
1971  inline void assign_to_var(int& n_lhs, const int& n_rhs) {
1972  n_lhs = n_rhs; // FIXME: no call -- just filler to instantiate
1973  }
1974  // FIXME: double val?
1975  inline void assign_to_var(double& n_lhs, const double& n_rhs) {
1976  n_lhs = n_rhs; // FIXME: no call -- just filler to instantiate
1977  }
1978 
1979  template <typename T_result, class Policy>
1980  inline bool check_pos_definite(const char* function,
1981  const Eigen::Matrix<var,Eigen::Dynamic,Eigen::Dynamic>& y,
1982  const char* name,
1983  T_result* result,
1984  const Policy&) {
1985  typedef
1986  typename Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic>::size_type
1987  size_type;
1988  Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> y_d(y.rows(),y.cols());
1989  for (size_type i = 0; i < y_d.rows(); i++)
1990  for (size_type j = 0; j < y_d.cols(); j++)
1991  y_d(i,j) = y(i,j).val();
1992  return stan::math::check_pos_definite(function,y_d,name,result,Policy());
1993  }
1994 
1995  template <typename LHS, typename RHS>
1996  inline void assign_to_var(std::vector<LHS>& x, const std::vector<RHS>& y) {
1997  for (size_t i = 0; i < x.size(); ++i)
1998  assign_to_var(x[i],y[i]);
1999  }
2000  template <typename LHS, typename RHS>
2001  inline void assign_to_var(Eigen::Matrix<LHS,Eigen::Dynamic,1>& x,
2002  const Eigen::Matrix<RHS,Eigen::Dynamic,1>& y) {
2003  typedef
2004  typename Eigen::Matrix<LHS,Eigen::Dynamic,1>::size_type
2005  size_type;
2006  for (size_type i = 0; i < x.size(); ++i)
2007  assign_to_var(x(i),y(i));
2008  }
2009  template <typename LHS, typename RHS>
2010  inline void assign_to_var(Eigen::Matrix<LHS,1,Eigen::Dynamic>& x,
2011  const Eigen::Matrix<RHS,1,Eigen::Dynamic>& y) {
2012  typedef
2013  typename Eigen::Matrix<LHS,1,Eigen::Dynamic>::size_type
2014  size_type;
2015  for (size_type i = 0; i < x.size(); ++i)
2016  assign_to_var(x(i),y(i));
2017  }
2018  template <typename LHS, typename RHS>
2019  inline void assign_to_var(Eigen::Matrix<LHS,Eigen::Dynamic,Eigen::Dynamic>& x,
2020  const Eigen::Matrix<RHS,Eigen::Dynamic,Eigen::Dynamic>& y) {
2021  typedef
2022  typename Eigen::Matrix<LHS,Eigen::Dynamic,Eigen::Dynamic>::size_type
2023  size_type1;
2024  typedef
2025  typename Eigen::Matrix<RHS,Eigen::Dynamic,Eigen::Dynamic>::size_type
2026  size_type2;
2027  for (size_type1 n = 0; n < x.cols(); ++n)
2028  for (size_type2 m = 0; m < x.rows(); ++m)
2029  assign_to_var(x(m,n),y(m,n));
2030  }
2031 
2032 
2033  template <typename LHS, typename RHS>
2034  struct needs_promotion {
2035  enum { value = ( is_constant_struct<RHS>::value
2036  && !is_constant_struct<LHS>::value) };
2037  };
2038 
2039  template <bool PromoteRHS, typename LHS, typename RHS>
2040  struct assigner {
2041  static inline void assign(LHS& var, const RHS& val) {
2042  throw std::domain_error("should not call base class of assigner");
2043  }
2044  };
2045 
2046  template <typename LHS, typename RHS>
2047  struct assigner<false,LHS,RHS> {
2048  static inline void assign(LHS& var, const RHS& val) {
2049  var = val; // no promotion of RHS
2050  }
2051  };
2052 
2053  template <typename LHS, typename RHS>
2054  struct assigner<true,LHS,RHS> {
2055  static inline void assign(LHS& var, const RHS& val) {
2056  assign_to_var(var,val); // promote RHS
2057  }
2058  };
2059 
2060 
2061  template <typename LHS, typename RHS>
2062  inline void assign(LHS& var, const RHS& val) {
2063  assigner<needs_promotion<LHS,RHS>::value, LHS, RHS>::assign(var,val);
2064  }
2065 
2066  void stan_print(std::ostream* o, const var& x);
2067 
2068  }
2069 }
2070 
2071 
2072 #endif
2073 
     [ Stan Home Page ] © 2011–2012, Stan Development Team.