Stan  1.1
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 #include <stan/math/matrix.hpp>
5 #include <stan/math/matrix_error_handling.hpp>
6 
7 #include <stan/agrad/special_functions.hpp>
8 
9 namespace stan {
10  namespace agrad {
11 
12 
13  void fill(var& x, const var& y) {
14  x = y;
15  }
16  template <typename T>
17  void fill(std::vector<T>& x, const var& y) {
18  for (size_t i = 0; i < x.size(); ++i)
19  fill(x[i],y);
20  }
21  template <int R, int C>
22  void fill(Eigen::Matrix<var,R,C>& x, const var& y) {
23  x.fill(y);
24  }
25 
26 
27 
28  class gevv_vvv_vari : public stan::agrad::vari {
29  protected:
30  stan::agrad::vari* alpha_;
31  stan::agrad::vari** v1_;
32  stan::agrad::vari** v2_;
33  double dotval_;
34  size_t length_;
35  inline static double eval_gevv(const stan::agrad::var* alpha,
36  const stan::agrad::var* v1, int stride1,
37  const stan::agrad::var* v2, int stride2,
38  size_t length, double *dotprod) {
39  double result = 0;
40  for (size_t i = 0; i < length; i++)
41  result += v1[i*stride1].vi_->val_ * v2[i*stride2].vi_->val_;
42  *dotprod = result;
43  return alpha->vi_->val_ * result;
44  }
45  public:
46  gevv_vvv_vari(const stan::agrad::var* alpha,
47  const stan::agrad::var* v1, int stride1,
48  const stan::agrad::var* v2, int stride2, size_t length) :
49  vari(eval_gevv(alpha,v1,stride1,v2,stride2,length,&dotval_)), length_(length) {
50  alpha_ = alpha->vi_;
51  v1_ = (stan::agrad::vari**)stan::agrad::memalloc_.alloc(2*length_*sizeof(stan::agrad::vari*));
52  v2_ = v1_ + length_;
53  for (size_t i = 0; i < length_; i++)
54  v1_[i] = v1[i*stride1].vi_;
55  for (size_t i = 0; i < length_; i++)
56  v2_[i] = v2[i*stride2].vi_;
57  }
58  void chain() {
59  const double adj_alpha = adj_ * alpha_->val_;
60  for (size_t i = 0; i < length_; i++) {
61  v1_[i]->adj_ += adj_alpha * v2_[i]->val_;
62  v2_[i]->adj_ += adj_alpha * v1_[i]->val_;
63  }
64  alpha_->adj_ += adj_ * dotval_;
65  }
66  };
67  }
68 }
69 
70 namespace Eigen {
71 
76  template <> struct NumTraits<stan::agrad::var>
77  {
83  typedef stan::agrad::var Real;
84 
90  typedef stan::agrad::var NonInteger;
91 
97  typedef stan::agrad::var Nested;
98 
105  inline static Real epsilon() {
106  return std::numeric_limits<double>::epsilon();
107  }
108 
112  inline static Real dummy_precision() {
113  return 1e-12; // copied from NumTraits.h values for double
114  }
115 
122  inline static Real highest() {
123  return std::numeric_limits<double>::max();
124  }
125 
132  inline static Real lowest() {
133  return -std::numeric_limits<double>::max();
134  }
135 
140  enum {
141  IsInteger = 0,
142  IsSigned = 1,
143  IsComplex = 0,
144  RequireInitialization = 0,
145  ReadCost = 1,
146  AddCost = 1,
147  MulCost = 1,
148  HasFloatingPoint = 1
149  };
150  };
151 
152  namespace internal {
156  template<>
157  struct significant_decimals_default_impl<stan::agrad::var,false>
158  {
159  static inline int run()
160  {
161  using std::ceil;
162  return cast<double,int>(ceil(-log(NumTraits<stan::agrad::var>::epsilon().val())
163  /log(10.0)));
164  }
165  };
166 
171  template <>
172  struct scalar_product_traits<stan::agrad::var,double> {
173  typedef stan::agrad::var ReturnType;
174  };
175 
180  template <>
181  struct scalar_product_traits<double,stan::agrad::var> {
182  typedef stan::agrad::var ReturnType;
183  };
184 
188  template<typename Index, bool ConjugateLhs, bool ConjugateRhs>
189  struct general_matrix_vector_product<Index,stan::agrad::var,ColMajor,ConjugateLhs,stan::agrad::var,ConjugateRhs>
190  {
191  typedef stan::agrad::var LhsScalar;
192  typedef stan::agrad::var RhsScalar;
193  typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar;
194  enum { LhsStorageOrder = ColMajor };
195 
196  EIGEN_DONT_INLINE static void run(
197  Index rows, Index cols,
198  const LhsScalar* lhs, Index lhsStride,
199  const RhsScalar* rhs, Index rhsIncr,
200  ResScalar* res, Index resIncr, const ResScalar &alpha)
201  {
202  for (Index i = 0; i < rows; i++) {
203  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));
204  }
205  }
206  };
207  template<typename Index, bool ConjugateLhs, bool ConjugateRhs>
208  struct general_matrix_vector_product<Index,stan::agrad::var,RowMajor,ConjugateLhs,stan::agrad::var,ConjugateRhs>
209  {
210  typedef stan::agrad::var LhsScalar;
211  typedef stan::agrad::var RhsScalar;
212  typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar;
213  enum { LhsStorageOrder = RowMajor };
214 
215  EIGEN_DONT_INLINE static void run(
216  Index rows, Index cols,
217  const LhsScalar* lhs, Index lhsStride,
218  const RhsScalar* rhs, Index rhsIncr,
219  ResScalar* res, Index resIncr, const RhsScalar &alpha)
220  {
221  for (Index i = 0; i < rows; i++) {
222  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));
223  }
224  }
225  };
226  template<typename Index, int LhsStorageOrder, bool ConjugateLhs, int RhsStorageOrder, bool ConjugateRhs>
227  struct general_matrix_matrix_product<Index,stan::agrad::var,LhsStorageOrder,ConjugateLhs,stan::agrad::var,RhsStorageOrder,ConjugateRhs,ColMajor>
228  {
229  typedef stan::agrad::var LhsScalar;
230  typedef stan::agrad::var RhsScalar;
231  typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar;
232  static void run(Index rows, Index cols, Index depth,
233  const LhsScalar* _lhs, Index lhsStride,
234  const RhsScalar* _rhs, Index rhsStride,
235  ResScalar* res, Index resStride,
236  const ResScalar &alpha,
237  level3_blocking<LhsScalar,RhsScalar>& /* blocking */,
238  GemmParallelInfo<Index>* /* info = 0 */)
239  {
240  for (Index i = 0; i < cols; i++) {
241  general_matrix_vector_product<Index,LhsScalar,LhsStorageOrder,ConjugateLhs,RhsScalar,ConjugateRhs>::run(
242  rows,depth,_lhs,lhsStride,
243  &_rhs[((int)RhsStorageOrder == (int)ColMajor)?(i*rhsStride):(i)],((int)RhsStorageOrder == (int)ColMajor)?(1):(rhsStride),
244  &res[i*resStride],1,alpha);
245  }
246  }
247  };
248  }
249 }
250 
251 namespace stan {
252 
253  namespace agrad {
254 
259  typedef
260  Eigen::Matrix<var,Eigen::Dynamic,Eigen::Dynamic>
261  matrix_v;
262 
267  typedef
268  Eigen::Matrix<var,Eigen::Dynamic,1>
269  vector_v;
270 
275  typedef
276  Eigen::Matrix<var,1,Eigen::Dynamic>
277  row_vector_v;
278 
283  inline void initialize_variable(var& variable, const var& value) {
284  variable = value;
285  }
286 
291  template <int R, int C>
292  inline void initialize_variable(Eigen::Matrix<var,R,C>& matrix, const var& value) {
293  for (int i = 0; i < matrix.size(); ++i)
294  matrix(i) = value;
295  }
296 
300  template <typename T>
301  inline void initialize_variable(std::vector<T>& variables, const var& value) {
302  for (size_t i = 0; i < variables.size(); ++i)
303  initialize_variable(variables[i],value);
304  }
305 
306 
315  inline var to_var(const double& x) {
316  return var(x);
317  }
326  inline var to_var(const var& x) {
327  return x;
328  }
337  inline matrix_v to_var(const stan::math::matrix_d& m) {
338  matrix_v m_v(m.rows(), m.cols());
339  for (int j = 0; j < m.cols(); ++j)
340  for (int i = 0; i < m.rows(); ++i)
341  m_v(i,j) = m(i,j);
342  return m_v;
343  }
352  inline matrix_v to_var(const matrix_v& m) {
353  return m;
354  }
364  inline vector_v to_var(const stan::math::vector_d& v) {
365  vector_v v_v(v.size());
366  for (int i = 0; i < v.size(); ++i)
367  v_v[i] = v[i];
368  return v_v;
369  }
379  inline vector_v to_var(const vector_v& v) {
380  return v;
381  }
391  inline row_vector_v to_var(const stan::math::row_vector_d& rv) {
392  row_vector_v rv_v(rv.size());
393  for (int i = 0; i < rv.size(); ++i)
394  rv_v[i] = rv[i];
395  return rv_v;
396  }
406  inline row_vector_v to_var(const row_vector_v& rv) {
407  return rv;
408  }
409 
410 
411 
412  // scalar returns
413 
414  // /**
415  // * Determinant of the matrix.
416  // *
417  // * Returns the determinant of the specified
418  // * square matrix.
419  // *
420  // * @param m Specified matrix.
421  // * @return Determinant of the matrix.
422  // * @throw std::domain_error if m is not a square matrix
423  // */
424  // var determinant(const Eigen::Matrix<var, Eigen::Dynamic, Eigen::Dynamic>& m);
425 
426  namespace {
427  class dot_self_vari : public vari {
428  protected:
429  vari** v_;
430  size_t size_;
431  public:
432  dot_self_vari(vari** v, size_t size)
433  : vari(var_dot_self(v,size)),
434  v_(v),
435  size_(size) {
436  }
437  template<typename Derived>
438  dot_self_vari(const Eigen::DenseBase<Derived> &v) :
439  vari(var_dot_self(v)), size_(v.size()) {
440  v_ = (vari**)memalloc_.alloc(size_*sizeof(vari*));
441  for (size_t i = 0; i < size_; i++)
442  v_[i] = v[i].vi_;
443  }
444  template <int R, int C>
445  dot_self_vari(const Eigen::Matrix<var,R,C>& v) :
446  vari(var_dot_self(v)), size_(v.size()) {
447  v_ = (vari**) memalloc_.alloc(size_ * sizeof(vari*));
448  for (size_t i = 0; i < size_; ++i)
449  v_[i] = v(i).vi_;
450  }
451  inline static double square(double x) { return x * x; }
452  inline static double var_dot_self(vari** v, size_t size) {
453  double sum = 0.0;
454  for (size_t i = 0; i < size; ++i)
455  sum += square(v[i]->val_);
456  return sum;
457  }
458  template<typename Derived>
459  double var_dot_self(const Eigen::DenseBase<Derived> &v) {
460  double sum = 0.0;
461  for (int i = 0; i < v.size(); ++i)
462  sum += square(v(i).vi_->val_);
463  return sum;
464  }
465  template <int R, int C>
466  inline static double var_dot_self(const Eigen::Matrix<var,R,C> &v) {
467  double sum = 0.0;
468  for (int i = 0; i < v.size(); ++i)
469  sum += square(v(i).vi_->val_);
470  return sum;
471  }
472  virtual void chain() {
473  for (size_t i = 0; i < size_; ++i)
474  v_[i]->adj_ += adj_ * 2.0 * v_[i]->val_;
475  }
476  };
477 
478  class sum_v_vari : public vari{
479  protected:
480  vari** v_;
481  size_t length_;
482  inline static double var_sum(const var *v, size_t length) {
483  double result = 0;
484  for (size_t i = 0; i < length; i++)
485  result += v[i].vi_->val_;
486  return result;
487  }
488  template<typename Derived>
489  inline static double var_sum(const Eigen::DenseBase<Derived> &v) {
490  double result = 0;
491  for (int i = 0; i < v.size(); i++)
492  result += v(i).vi_->val_;
493  return result;
494  }
495  public:
496  template<typename Derived>
497  sum_v_vari(const Eigen::DenseBase<Derived> &v) :
498  vari(var_sum(v)), length_(v.size()) {
499  v_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
500  for (size_t i = 0; i < length_; i++)
501  v_[i] = v(i).vi_;
502  }
503  template<int R1,int C1>
504  sum_v_vari(const Eigen::Matrix<var,R1,C1> &v1) :
505  vari(var_sum(v1)), length_(v1.size()) {
506  v_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
507  for (size_t i = 0; i < length_; i++)
508  v_[i] = v1(i).vi_;
509  }
510  sum_v_vari(const var *v, size_t len) :
511  vari(var_sum(v,len)), length_(len) {
512  v_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
513  for (size_t i = 0; i < length_; i++)
514  v_[i] = v[i].vi_;
515  }
516  virtual void chain() {
517  for (size_t i = 0; i < length_; i++) {
518  v_[i]->adj_ += adj_;
519  }
520  }
521  };
522  class dot_product_vv_vari : public vari {
523  protected:
524  vari** v1_;
525  vari** v2_;
526  size_t length_;
527  inline static double var_dot(const var* v1, const var* v2,
528  size_t length) {
529  double result = 0;
530  for (size_t i = 0; i < length; i++)
531  result += v1[i].vi_->val_ * v2[i].vi_->val_;
532  return result;
533  }
534  template<typename Derived1,typename Derived2>
535  inline static double var_dot(const Eigen::DenseBase<Derived1> &v1,
536  const Eigen::DenseBase<Derived2> &v2) {
537  double result = 0;
538  for (int i = 0; i < v1.size(); i++)
539  result += v1[i].vi_->val_ * v2[i].vi_->val_;
540  return result;
541  }
542  inline static double var_dot(vari** v1, vari** v2, size_t length) {
543  double result = 0;
544  for (size_t i = 0; i < length; ++i)
545  result += v1[i]->val_ * v2[i]->val_;
546  return result;
547  }
548  public:
549  dot_product_vv_vari(vari** v1, vari** v2, size_t length)
550  : vari(var_dot(v1,v2,length)),
551  v1_(v1),
552  v2_(v2),
553  length_(length) {
554 
555  }
556  dot_product_vv_vari(const var* v1, const var* v2, size_t length,
557  dot_product_vv_vari* shared_v1 = NULL,
558  dot_product_vv_vari* shared_v2 = NULL) :
559  vari(var_dot(v1, v2, length)), length_(length) {
560  if (shared_v1 == NULL) {
561  v1_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
562  for (size_t i = 0; i < length_; i++)
563  v1_[i] = v1[i].vi_;
564  }
565  else {
566  v1_ = shared_v1->v1_;
567  }
568  if (shared_v2 == NULL) {
569  v2_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
570  for (size_t i = 0; i < length_; i++)
571  v2_[i] = v2[i].vi_;
572  }
573  else {
574  v2_ = shared_v2->v2_;
575  }
576  }
577  template<typename Derived1,typename Derived2>
578  dot_product_vv_vari(const Eigen::DenseBase<Derived1> &v1,
579  const Eigen::DenseBase<Derived2> &v2,
580  dot_product_vv_vari* shared_v1 = NULL,
581  dot_product_vv_vari* shared_v2 = NULL) :
582  vari(var_dot(v1, v2)), length_(v1.size()) {
583  if (shared_v1 == NULL) {
584  v1_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
585  for (size_t i = 0; i < length_; i++)
586  v1_[i] = v1[i].vi_;
587  }
588  else {
589  v1_ = shared_v1->v1_;
590  }
591  if (shared_v2 == NULL) {
592  v2_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
593  for (size_t i = 0; i < length_; i++)
594  v2_[i] = v2[i].vi_;
595  }
596  else {
597  v2_ = shared_v2->v2_;
598  }
599  }
600  template<int R1,int C1,int R2,int C2>
601  dot_product_vv_vari(const Eigen::Matrix<var,R1,C1> &v1,
602  const Eigen::Matrix<var,R2,C2> &v2,
603  dot_product_vv_vari* shared_v1 = NULL,
604  dot_product_vv_vari* shared_v2 = NULL) :
605  vari(var_dot(v1, v2)), length_(v1.size()) {
606  if (shared_v1 == NULL) {
607  v1_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
608  for (size_t i = 0; i < length_; i++)
609  v1_[i] = v1[i].vi_;
610  }
611  else {
612  v1_ = shared_v1->v1_;
613  }
614  if (shared_v2 == NULL) {
615  v2_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
616  for (size_t i = 0; i < length_; i++)
617  v2_[i] = v2[i].vi_;
618  }
619  else {
620  v2_ = shared_v2->v2_;
621  }
622  }
623  virtual void chain() {
624  for (size_t i = 0; i < length_; i++) {
625  v1_[i]->adj_ += adj_ * v2_[i]->val_;
626  v2_[i]->adj_ += adj_ * v1_[i]->val_;
627  }
628  }
629  };
630 
631  class dot_product_vd_vari : public vari {
632  protected:
633  vari** v1_;
634  double* v2_;
635  size_t length_;
636  inline static double var_dot(const var* v1, const double* v2,
637  size_t length) {
638  double result = 0;
639  for (size_t i = 0; i < length; i++)
640  result += v1[i].vi_->val_ * v2[i];
641  return result;
642  }
643  template<typename Derived1,typename Derived2>
644  inline static double var_dot(const Eigen::DenseBase<Derived1> &v1,
645  const Eigen::DenseBase<Derived2> &v2) {
646  double result = 0;
647  for (int i = 0; i < v1.size(); i++)
648  result += v1[i].vi_->val_ * v2[i];
649  return result;
650  }
651  public:
652  dot_product_vd_vari(const var* v1, const double* v2, size_t length,
653  dot_product_vd_vari *shared_v1 = NULL,
654  dot_product_vd_vari *shared_v2 = NULL) :
655  vari(var_dot(v1, v2, length)), length_(length) {
656  if (shared_v1 == NULL) {
657  v1_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
658  for (size_t i = 0; i < length_; i++)
659  v1_[i] = v1[i].vi_;
660  } else {
661  v1_ = shared_v1->v1_;
662  }
663  if (shared_v2 == NULL) {
664  v2_ = (double*)memalloc_.alloc(length_*sizeof(double));
665  for (size_t i = 0; i < length_; i++)
666  v2_[i] = v2[i];
667  } else {
668  v2_ = shared_v2->v2_;
669  }
670  }
671  template<typename Derived1,typename Derived2>
672  dot_product_vd_vari(const Eigen::DenseBase<Derived1> &v1,
673  const Eigen::DenseBase<Derived2> &v2,
674  dot_product_vd_vari *shared_v1 = NULL,
675  dot_product_vd_vari *shared_v2 = NULL) :
676  vari(var_dot(v1, v2)), length_(v1.size()) {
677  if (shared_v1 == NULL) {
678  v1_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
679  for (size_t i = 0; i < length_; i++)
680  v1_[i] = v1[i].vi_;
681  } else {
682  v1_ = shared_v1->v1_;
683  }
684  if (shared_v2 == NULL) {
685  v2_ = (double*)memalloc_.alloc(length_*sizeof(double));
686  for (size_t i = 0; i < length_; i++)
687  v2_[i] = v2[i];
688  } else {
689  v2_ = shared_v2->v2_;
690  }
691  }
692  template<int R1,int C1,int R2,int C2>
693  dot_product_vd_vari(const Eigen::Matrix<var,R1,C1> &v1,
694  const Eigen::Matrix<double,R2,C2> &v2,
695  dot_product_vd_vari *shared_v1 = NULL,
696  dot_product_vd_vari *shared_v2 = NULL) :
697  vari(var_dot(v1, v2)), length_(v1.size()) {
698  if (shared_v1 == NULL) {
699  v1_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
700  for (size_t i = 0; i < length_; i++)
701  v1_[i] = v1[i].vi_;
702  } else {
703  v1_ = shared_v1->v1_;
704  }
705  if (shared_v2 == NULL) {
706  v2_ = (double*)memalloc_.alloc(length_*sizeof(double));
707  for (size_t i = 0; i < length_; i++)
708  v2_[i] = v2[i];
709  } else {
710  v2_ = shared_v2->v2_;
711  }
712  }
713  virtual void chain() {
714  for (size_t i = 0; i < length_; i++) {
715  v1_[i]->adj_ += adj_ * v2_[i];
716  }
717  }
718  };
719 
720  }
721 
732  template<int R, int C>
733  inline var dot_self(const Eigen::Matrix<var, R, C>& v) {
734  stan::math::validate_vector(v,"dot_self");
735  return var(new dot_self_vari(v));
736  }
737 
743  template<int R,int C>
744  inline Eigen::Matrix<var,1,C>
745  columns_dot_self(const Eigen::Matrix<var,R,C>& x) {
746  Eigen::Matrix<var,1,C> ret(1,x.cols());
747  size_t i;
748  for (i = 0; i < x.cols(); i++) {
749  ret(i) = var(new dot_self_vari(x.col(i)));
750  }
751  return ret;
752  }
753 
754  template<int R1,int C1,int R2, int C2>
755  inline Eigen::Matrix<var, 1, C1>
756  columns_dot_product(const Eigen::Matrix<var, R1, C1>& v1,
757  const Eigen::Matrix<var, R2, C2>& v2) {
758  stan::math::validate_matching_sizes(v1,v2,"columns_dot_product");
759  Eigen::Matrix<var, 1, C1> ret(1,v1.cols());
760  size_t j;
761  for (j = 0; j < v1.cols(); ++j) {
762  ret(j) = var(new dot_product_vv_vari(v1.col(j),v2.col(j)));
763  }
764  return ret;
765  }
766 
767  template<int R1,int C1,int R2, int C2>
768  inline Eigen::Matrix<var, 1, C1>
769  columns_dot_product(const Eigen::Matrix<var, R1, C1>& v1,
770  const Eigen::Matrix<double, R2, C2>& v2) {
771  stan::math::validate_matching_sizes(v1,v2,"columns_dot_product");
772  Eigen::Matrix<var, 1, C1> ret(1,v1.cols());
773  size_t j;
774  for (j = 0; j < v1.cols(); ++j) {
775  ret(j) = var(new dot_product_vd_vari(v1.col(j),v2.col(j)));
776  }
777  return ret;
778  }
779 
780  template<int R1,int C1,int R2, int C2>
781  inline Eigen::Matrix<var, 1, C1>
782  columns_dot_product(const Eigen::Matrix<double, R1, C1>& v1,
783  const Eigen::Matrix<var, R2, C2>& v2) {
784  stan::math::validate_matching_sizes(v1,v2,"columns_dot_product");
785  Eigen::Matrix<var, 1, C1> ret(1,v1.cols());
786  size_t j;
787  for (j = 0; j < v1.cols(); ++j) {
788  ret(j) = var(new dot_product_vd_vari(v2.col(j),v1.col(j)));
789  }
790  return ret;
791  }
792 
801  template<int R1,int C1,int R2, int C2>
802  inline var dot_product(const Eigen::Matrix<var, R1, C1>& v1,
803  const Eigen::Matrix<var, R2, C2>& v2) {
804  stan::math::validate_vector(v1,"dot_product");
805  stan::math::validate_vector(v2,"dot_product");
806  stan::math::validate_matching_sizes(v1,v2,"dot_product");
807  return var(new dot_product_vv_vari(v1,v2));
808  }
818  template<int R1,int C1,int R2, int C2>
819  inline var dot_product(const Eigen::Matrix<var, R1, C1>& v1,
820  const Eigen::Matrix<double, R2, C2>& v2) {
821  stan::math::validate_vector(v1,"dot_product");
822  stan::math::validate_vector(v2,"dot_product");
823  stan::math::validate_matching_sizes(v1,v2,"dot_product");
824  return var(new dot_product_vd_vari(v1,v2));
825  }
835  template<int R1,int C1,int R2, int C2>
836  inline var dot_product(const Eigen::Matrix<double, R1, C1>& v1,
837  const Eigen::Matrix<var, R2, C2>& v2) {
838  stan::math::validate_vector(v1,"dot_product");
839  stan::math::validate_vector(v2,"dot_product");
840  stan::math::validate_matching_sizes(v1,v2,"dot_product");
841  return var(new dot_product_vd_vari(v2,v1));
842  }
851  inline var dot_product(const var* v1, const var* v2, size_t length) {
852  return var(new dot_product_vv_vari(v1, v2, length));
853  }
862  inline var dot_product(const var* v1, const double* v2, size_t length) {
863  return var(new dot_product_vd_vari(v1, v2, length));
864  }
873  inline var dot_product(const double* v1, const var* v2, size_t length) {
874  return var(new dot_product_vd_vari(v2, v1, length));
875  }
884  inline var dot_product(const std::vector<var>& v1,
885  const std::vector<var>& v2) {
886  stan::math::validate_matching_sizes(v1,v2,"dot_product");
887  return var(new dot_product_vv_vari(&v1[0], &v2[0], v1.size()));
888  }
897  inline var dot_product(const std::vector<var>& v1,
898  const std::vector<double>& v2) {
899  stan::math::validate_matching_sizes(v1,v2,"dot_product");
900  return var(new dot_product_vd_vari(&v1[0], &v2[0], v1.size()));
901  }
910  inline var dot_product(const std::vector<double>& v1,
911  const std::vector<var>& v2) {
912  stan::math::validate_matching_sizes(v1,v2,"dot_product");
913  return var(new dot_product_vd_vari(&v2[0], &v1[0], v1.size()));
914  }
915 
922  template <int R, int C>
923  inline var sum(const Eigen::Matrix<var,R,C>& m) {
924  if (m.size() == 0)
925  return 0.0;
926  return var(new sum_v_vari(m));
927  }
928 
929  template <int R1,int C1,int R2,int C2>
930  class mdivide_left_vv_vari : public vari {
931  public:
932  Eigen::Matrix<double,R1,C1> _A;
933  Eigen::Matrix<double,R1,C2> _C;
934 
935  Eigen::Matrix<vari*,R1,C1> _variRefA;
936  Eigen::Matrix<vari*,R2,C2> _variRefB;
937  Eigen::Matrix<vari*,R1,C2> _variRefC;
938 
939  mdivide_left_vv_vari(const Eigen::Matrix<var,R1,C1> &A,
940  const Eigen::Matrix<var,R2,C2> &B)
941  : vari(0.0),
942  _A(A.rows(),A.cols()), _C(B.rows(),B.cols()),
943  _variRefA(A.rows(),A.cols()), _variRefB(B.rows(),B.cols()),
944  _variRefC(B.rows(),B.cols())
945  {
946  size_t i,j;
947 
948  for (j = 0; j < _variRefA.cols(); j++) {
949  for (i = 0; i < _variRefA.rows(); i++) {
950  _variRefA(i,j) = A(i,j).vi_;
951  _A(i,j) = A(i,j).val();
952  }
953  }
954 
955  for (j = 0; j < _variRefB.cols(); j++) {
956  for (i = 0; i < _variRefB.rows(); i++) {
957  _variRefB(i,j) = B(i,j).vi_;
958  _C(i,j) = B(i,j).val();
959  }
960  }
961 
962  _C = _A.colPivHouseholderQr().solve(_C);
963 
964  for (j = 0; j < _variRefC.cols(); j++) {
965  for (i = 0; i < _variRefC.rows(); i++) {
966  _variRefC(i,j) = new vari(_C(i,j),false);
967  }
968  }
969  }
970 
971  virtual void chain() {
972  Eigen::Matrix<double,R1,C1> adjA(_variRefA.rows(),_variRefA.cols());
973  Eigen::Matrix<double,R2,C2> adjB(_variRefB.rows(),_variRefB.cols());
974  Eigen::Matrix<double,R1,C2> adjC(_variRefC.rows(),_variRefC.cols());
975 
976  size_t i,j;
977  for (j = 0; j < adjC.cols(); j++)
978  for (i = 0; i < adjC.rows(); i++)
979  adjC(i,j) = _variRefC(i,j)->adj_;
980 
981 
982  adjB = _A.transpose().colPivHouseholderQr().solve(adjC);
983  adjA.noalias() = -adjB*_C.transpose();
984 
985  for (j = 0; j < adjA.cols(); j++)
986  for (i = 0; i < adjA.rows(); i++)
987  _variRefA(i,j)->adj_ += adjA(i,j);
988 
989  for (j = 0; j < adjB.cols(); j++)
990  for (i = 0; i < adjB.rows(); i++)
991  _variRefB(i,j)->adj_ += adjB(i,j);
992  }
993  };
994 
995  template <int R1,int C1,int R2,int C2>
996  class mdivide_left_dv_vari : public vari {
997  public:
998  Eigen::Matrix<double,R1,C1> _A;
999  Eigen::Matrix<double,R1,C2> _C;
1000 
1001  Eigen::Matrix<vari*,R2,C2> _variRefB;
1002  Eigen::Matrix<vari*,R1,C2> _variRefC;
1003 
1004  mdivide_left_dv_vari(const Eigen::Matrix<double,R1,C1> &A,
1005  const Eigen::Matrix<var,R2,C2> &B)
1006  : vari(0.0),
1007  _A(A), _C(B.rows(),B.cols()),
1008  _variRefB(B.rows(),B.cols()), _variRefC(B.rows(),B.cols())
1009  {
1010  size_t i,j;
1011 
1012  for (j = 0; j < _variRefB.cols(); j++) {
1013  for (i = 0; i < _variRefB.rows(); i++) {
1014  _variRefB(i,j) = B(i,j).vi_;
1015  _C(i,j) = B(i,j).val();
1016  }
1017  }
1018 
1019  _C = _A.colPivHouseholderQr().solve(_C);
1020 
1021  for (j = 0; j < _variRefC.cols(); j++) {
1022  for (i = 0; i < _variRefC.rows(); i++) {
1023  _variRefC(i,j) = new vari(_C(i,j),false);
1024  }
1025  }
1026  }
1027 
1028  virtual void chain() {
1029  Eigen::Matrix<double,R2,C2> adjB(_variRefB.rows(),_variRefB.cols());
1030  Eigen::Matrix<double,R1,C2> adjC(_variRefC.rows(),_variRefC.cols());
1031 
1032  size_t i,j;
1033  for (j = 0; j < adjC.cols(); j++)
1034  for (i = 0; i < adjC.rows(); i++)
1035  adjC(i,j) = _variRefC(i,j)->adj_;
1036 
1037 
1038  adjB = _A.transpose().colPivHouseholderQr().solve(adjC);
1039 
1040  for (j = 0; j < adjB.cols(); j++)
1041  for (i = 0; i < adjB.rows(); i++)
1042  _variRefB(i,j)->adj_ += adjB(i,j);
1043  }
1044  };
1045 
1046  template <int R1,int C1,int R2,int C2>
1047  class mdivide_left_vd_vari : public vari {
1048  public:
1049  Eigen::Matrix<double,R1,C1> _A;
1050  Eigen::Matrix<double,R1,C2> _C;
1051 
1052  Eigen::Matrix<vari*,R1,C1> _variRefA;
1053  Eigen::Matrix<vari*,R1,C2> _variRefC;
1054 
1055  mdivide_left_vd_vari(const Eigen::Matrix<var,R1,C1> &A,
1056  const Eigen::Matrix<double,R2,C2> &B)
1057  : vari(0.0),
1058  _A(A.rows(),A.cols()), _C(B.rows(),B.cols()),
1059  _variRefA(A.rows(),A.cols()), _variRefC(B.rows(),B.cols())
1060  {
1061  size_t i,j;
1062 
1063  for (j = 0; j < _variRefA.cols(); j++) {
1064  for (i = 0; i < _variRefA.rows(); i++) {
1065  _variRefA(i,j) = A(i,j).vi_;
1066  _A(i,j) = A(i,j).val();
1067  }
1068  }
1069 
1070 
1071  _C = _A.colPivHouseholderQr().solve(B);
1072 
1073  for (j = 0; j < _variRefC.cols(); j++) {
1074  for (i = 0; i < _variRefC.rows(); i++) {
1075  _variRefC(i,j) = new vari(_C(i,j),false);
1076  }
1077  }
1078  }
1079 
1080  virtual void chain() {
1081  Eigen::Matrix<double,R1,C1> adjA(_variRefA.rows(),_variRefA.cols());
1082  Eigen::Matrix<double,R1,C2> adjC(_variRefC.rows(),_variRefC.cols());
1083 
1084  size_t i,j;
1085  for (j = 0; j < adjC.cols(); j++)
1086  for (i = 0; i < adjC.rows(); i++)
1087  adjC(i,j) = _variRefC(i,j)->adj_;
1088 
1089  adjA = -_A.transpose().colPivHouseholderQr().solve(adjC*_C.transpose());
1090 
1091  for (j = 0; j < adjA.cols(); j++)
1092  for (i = 0; i < adjA.rows(); i++)
1093  _variRefA(i,j)->adj_ += adjA(i,j);
1094  }
1095  };
1096 
1097  template <int R1,int C1,int R2,int C2>
1098  inline
1099  Eigen::Matrix<var,R1,C2>
1100  mdivide_left(const Eigen::Matrix<var,R1,C1> &A,
1101  const Eigen::Matrix<var,R2,C2> &b) {
1102  Eigen::Matrix<var,R1,C2> res(b.rows(),b.cols());
1103 
1104  stan::math::validate_square(A,"mdivide_left");
1105  stan::math::validate_multiplicable(A,b,"mdivide_left");
1106 
1107  // NOTE: this is not a memory leak, this vari is used in the
1108  // expression graph to evaluate the adjoint, but is not needed
1109  // for the returned matrix. Memory will be cleaned up with the arena allocator.
1110  mdivide_left_vv_vari<R1,C1,R2,C2> *baseVari = new mdivide_left_vv_vari<R1,C1,R2,C2>(A,b);
1111 
1112  size_t i,j;
1113 
1114  for (i = 0; i < res.rows(); i++)
1115  for (j = 0; j < res.cols(); j++)
1116  res(i,j).vi_ = baseVari->_variRefC(i,j);
1117 
1118  return res;
1119  }
1120 
1121  template <int R1,int C1,int R2,int C2>
1122  inline
1123  Eigen::Matrix<var,R1,C2>
1124  mdivide_left(const Eigen::Matrix<var,R1,C1> &A,
1125  const Eigen::Matrix<double,R2,C2> &b) {
1126  Eigen::Matrix<var,R1,C2> res(b.rows(),b.cols());
1127 
1128  stan::math::validate_square(A,"mdivide_left");
1129  stan::math::validate_multiplicable(A,b,"mdivide_left");
1130 
1131  // NOTE: this is not a memory leak, this vari is used in the
1132  // expression graph to evaluate the adjoint, but is not needed
1133  // for the returned matrix. Memory will be cleaned up with the arena allocator.
1134  mdivide_left_vd_vari<R1,C1,R2,C2> *baseVari = new mdivide_left_vd_vari<R1,C1,R2,C2>(A,b);
1135 
1136  size_t i,j;
1137 
1138  for (i = 0; i < res.rows(); i++)
1139  for (j = 0; j < res.cols(); j++)
1140  res(i,j).vi_ = baseVari->_variRefC(i,j);
1141 
1142  return res;
1143  }
1144 
1145  template <int R1,int C1,int R2,int C2>
1146  inline
1147  Eigen::Matrix<var,R1,C2>
1148  mdivide_left(const Eigen::Matrix<double,R1,C1> &A,
1149  const Eigen::Matrix<var,R2,C2> &b) {
1150  Eigen::Matrix<var,R1,C2> res(b.rows(),b.cols());
1151 
1152  stan::math::validate_square(A,"mdivide_left");
1153  stan::math::validate_multiplicable(A,b,"mdivide_left");
1154 
1155  // NOTE: this is not a memory leak, this vari is used in the
1156  // expression graph to evaluate the adjoint, but is not needed
1157  // for the returned matrix. Memory will be cleaned up with the arena allocator.
1158  mdivide_left_dv_vari<R1,C1,R2,C2> *baseVari = new mdivide_left_dv_vari<R1,C1,R2,C2>(A,b);
1159 
1160  size_t i,j;
1161 
1162  for (i = 0; i < res.rows(); i++)
1163  for (j = 0; j < res.cols(); j++)
1164  res(i,j).vi_ = baseVari->_variRefC(i,j);
1165 
1166  return res;
1167  }
1168 
1169  template <int TriView,int R1,int C1,int R2,int C2>
1170  class mdivide_left_tri_vv_vari : public vari {
1171  public:
1172  Eigen::Matrix<double,R1,C1> _A;
1173  Eigen::Matrix<double,R1,C2> _C;
1174 
1175  Eigen::Matrix<vari*,R1,C1> _variRefA;
1176  Eigen::Matrix<vari*,R2,C2> _variRefB;
1177  Eigen::Matrix<vari*,R1,C2> _variRefC;
1178 
1179  mdivide_left_tri_vv_vari(const Eigen::Matrix<var,R1,C1> &A,
1180  const Eigen::Matrix<var,R2,C2> &B)
1181  : vari(0.0),
1182  _A(A.rows(),A.cols()), _C(B.rows(),B.cols()),
1183  _variRefA(A.rows(),A.cols()), _variRefB(B.rows(),B.cols()),
1184  _variRefC(B.rows(),B.cols())
1185  {
1186  size_t i,j;
1187 
1188  for (j = 0; j < _variRefA.cols(); j++) {
1189  for (i = 0; i < _variRefA.rows(); i++) {
1190  _variRefA(i,j) = A(i,j).vi_;
1191  _A(i,j) = A(i,j).val();
1192  }
1193  }
1194 
1195  for (j = 0; j < _variRefB.cols(); j++) {
1196  for (i = 0; i < _variRefB.rows(); i++) {
1197  _variRefB(i,j) = B(i,j).vi_;
1198  _C(i,j) = B(i,j).val();
1199  }
1200  }
1201 
1202  _C = _A.template triangularView<TriView>().solve(_C);
1203 
1204  for (j = 0; j < _variRefC.cols(); j++) {
1205  for (i = 0; i < _variRefC.rows(); i++) {
1206  _variRefC(i,j) = new vari(_C(i,j),false);
1207  }
1208  }
1209  }
1210 
1211  virtual void chain() {
1212  Eigen::Matrix<double,R1,C1> adjA(_variRefA.rows(),_variRefA.cols());
1213  Eigen::Matrix<double,R2,C2> adjB(_variRefB.rows(),_variRefB.cols());
1214  Eigen::Matrix<double,R1,C2> adjC(_variRefC.rows(),_variRefC.cols());
1215 
1216  size_t i,j;
1217  for (j = 0; j < adjC.cols(); j++)
1218  for (i = 0; i < adjC.rows(); i++)
1219  adjC(i,j) = _variRefC(i,j)->adj_;
1220 
1221 
1222  adjB = _A.template triangularView<TriView>().transpose().solve(adjC);
1223  adjA.noalias() = -adjB*_C.transpose();
1224 
1225  if (TriView == Eigen::Lower) {
1226  for (j = 0; j < adjA.cols(); j++)
1227  for (i = j; i < adjA.rows(); i++)
1228  _variRefA(i,j)->adj_ += adjA(i,j);
1229  }
1230  else if (TriView == Eigen::Upper) {
1231  for (j = 0; j < adjA.cols(); j++)
1232  for (i = 0; i < j+1; i++)
1233  _variRefA(i,j)->adj_ += adjA(i,j);
1234  }
1235  else {
1236  assert(TriView == Eigen::Lower || TriView == Eigen::Upper);
1237  }
1238 
1239  for (j = 0; j < adjB.cols(); j++)
1240  for (i = 0; i < adjB.rows(); i++)
1241  _variRefB(i,j)->adj_ += adjB(i,j);
1242  }
1243  };
1244 
1245  template <int TriView,int R1,int C1,int R2,int C2>
1246  class mdivide_left_tri_dv_vari : public vari {
1247  public:
1248  Eigen::Matrix<double,R1,C1> _A;
1249  Eigen::Matrix<double,R1,C2> _C;
1250 
1251  Eigen::Matrix<vari*,R2,C2> _variRefB;
1252  Eigen::Matrix<vari*,R1,C2> _variRefC;
1253 
1254  mdivide_left_tri_dv_vari(const Eigen::Matrix<double,R1,C1> &A,
1255  const Eigen::Matrix<var,R2,C2> &B)
1256  : vari(0.0),
1257  _A(A), _C(B.rows(),B.cols()),
1258  _variRefB(B.rows(),B.cols()), _variRefC(B.rows(),B.cols())
1259  {
1260  size_t i,j;
1261 
1262  for (j = 0; j < _variRefB.cols(); j++) {
1263  for (i = 0; i < _variRefB.rows(); i++) {
1264  _variRefB(i,j) = B(i,j).vi_;
1265  _C(i,j) = B(i,j).val();
1266  }
1267  }
1268 
1269  _C = _A.template triangularView<TriView>().solve(_C);
1270 
1271  for (j = 0; j < _variRefC.cols(); j++) {
1272  for (i = 0; i < _variRefC.rows(); i++) {
1273  _variRefC(i,j) = new vari(_C(i,j),false);
1274  }
1275  }
1276  }
1277 
1278  virtual void chain() {
1279  Eigen::Matrix<double,R2,C2> adjB(_variRefB.rows(),_variRefB.cols());
1280  Eigen::Matrix<double,R1,C2> adjC(_variRefC.rows(),_variRefC.cols());
1281 
1282  size_t i,j;
1283  for (j = 0; j < adjC.cols(); j++)
1284  for (i = 0; i < adjC.rows(); i++)
1285  adjC(i,j) = _variRefC(i,j)->adj_;
1286 
1287 
1288  adjB = _A.template triangularView<TriView>().transpose().solve(adjC);
1289 
1290  for (j = 0; j < adjB.cols(); j++)
1291  for (i = 0; i < adjB.rows(); i++)
1292  _variRefB(i,j)->adj_ += adjB(i,j);
1293  }
1294  };
1295 
1296  template <int TriView,int R1,int C1,int R2,int C2>
1297  class mdivide_left_tri_vd_vari : public vari {
1298  public:
1299  Eigen::Matrix<double,R1,C1> _A;
1300  Eigen::Matrix<double,R1,C2> _C;
1301 
1302  Eigen::Matrix<vari*,R1,C1> _variRefA;
1303  Eigen::Matrix<vari*,R1,C2> _variRefC;
1304 
1305  mdivide_left_tri_vd_vari(const Eigen::Matrix<var,R1,C1> &A,
1306  const Eigen::Matrix<double,R2,C2> &B)
1307  : vari(0.0),
1308  _A(A.rows(),A.cols()), _C(B.rows(),B.cols()),
1309  _variRefA(A.rows(),A.cols()), _variRefC(B.rows(),B.cols())
1310  {
1311  size_t i,j;
1312 
1313  for (j = 0; j < _variRefA.cols(); j++) {
1314  for (i = 0; i < _variRefA.rows(); i++) {
1315  _variRefA(i,j) = A(i,j).vi_;
1316  _A(i,j) = A(i,j).val();
1317  }
1318  }
1319 
1320  _C = _A.template triangularView<TriView>().solve(B);
1321 
1322  for (j = 0; j < _variRefC.cols(); j++) {
1323  for (i = 0; i < _variRefC.rows(); i++) {
1324  _variRefC(i,j) = new vari(_C(i,j),false);
1325  }
1326  }
1327  }
1328 
1329  virtual void chain() {
1330  Eigen::Matrix<double,R1,C1> adjA(_variRefA.rows(),_variRefA.cols());
1331  Eigen::Matrix<double,R1,C2> adjC(_variRefC.rows(),_variRefC.cols());
1332 
1333  size_t i,j;
1334  for (j = 0; j < adjC.cols(); j++)
1335  for (i = 0; i < adjC.rows(); i++)
1336  adjC(i,j) = _variRefC(i,j)->adj_;
1337 
1338 
1339  adjA.noalias() = -_A.template triangularView<TriView>().transpose().solve(adjC*_C.transpose());
1340 
1341  if (TriView == Eigen::Lower) {
1342  for (j = 0; j < adjA.cols(); j++)
1343  for (i = j; i < adjA.rows(); i++)
1344  _variRefA(i,j)->adj_ += adjA(i,j);
1345  }
1346  else if (TriView == Eigen::Upper) {
1347  for (j = 0; j < adjA.cols(); j++)
1348  for (i = 0; i < j+1; i++)
1349  _variRefA(i,j)->adj_ += adjA(i,j);
1350  }
1351  else {
1352  assert(TriView == Eigen::Lower || TriView == Eigen::Upper);
1353  }
1354  }
1355  };
1356 
1357  template <int TriView,int R1,int C1,int R2,int C2>
1358  inline
1359  Eigen::Matrix<var,R1,C2>
1360  mdivide_left_tri(const Eigen::Matrix<var,R1,C1> &A,
1361  const Eigen::Matrix<var,R2,C2> &b) {
1362  Eigen::Matrix<var,R1,C2> res(b.rows(),b.cols());
1363 
1364  stan::math::validate_square(A,"mdivide_left_tri");
1365  stan::math::validate_multiplicable(A,b,"mdivide_left_tri");
1366 
1367  // NOTE: this is not a memory leak, this vari is used in the
1368  // expression graph to evaluate the adjoint, but is not needed
1369  // for the returned matrix. Memory will be cleaned up with the arena allocator.
1370  mdivide_left_tri_vv_vari<TriView,R1,C1,R2,C2> *baseVari = new mdivide_left_tri_vv_vari<TriView,R1,C1,R2,C2>(A,b);
1371 
1372  size_t i,j;
1373 
1374  for (i = 0; i < res.rows(); i++)
1375  for (j = 0; j < res.cols(); j++)
1376  res(i,j).vi_ = baseVari->_variRefC(i,j);
1377 
1378  return res;
1379  }
1380  template <int TriView,int R1,int C1,int R2,int C2>
1381  inline
1382  Eigen::Matrix<var,R1,C2>
1383  mdivide_left_tri(const Eigen::Matrix<double,R1,C1> &A,
1384  const Eigen::Matrix<var,R2,C2> &b) {
1385  Eigen::Matrix<var,R1,C2> res(b.rows(),b.cols());
1386 
1387  stan::math::validate_square(A,"mdivide_left_tri");
1388  stan::math::validate_multiplicable(A,b,"mdivide_left_tri");
1389 
1390  // NOTE: this is not a memory leak, this vari is used in the
1391  // expression graph to evaluate the adjoint, but is not needed
1392  // for the returned matrix. Memory will be cleaned up with the arena allocator.
1393  mdivide_left_tri_dv_vari<TriView,R1,C1,R2,C2> *baseVari = new mdivide_left_tri_dv_vari<TriView,R1,C1,R2,C2>(A,b);
1394 
1395  size_t i,j;
1396 
1397  for (i = 0; i < res.rows(); i++)
1398  for (j = 0; j < res.cols(); j++)
1399  res(i,j).vi_ = baseVari->_variRefC(i,j);
1400 
1401  return res;
1402  }
1403  template <int TriView,int R1,int C1,int R2,int C2>
1404  inline
1405  Eigen::Matrix<var,R1,C2>
1406  mdivide_left_tri(const Eigen::Matrix<var,R1,C1> &A,
1407  const Eigen::Matrix<double,R2,C2> &b) {
1408  Eigen::Matrix<var,R1,C2> res(b.rows(),b.cols());
1409 
1410  stan::math::validate_square(A,"mdivide_left_tri");
1411  stan::math::validate_multiplicable(A,b,"mdivide_left_tri");
1412 
1413  // NOTE: this is not a memory leak, this vari is used in the
1414  // expression graph to evaluate the adjoint, but is not needed
1415  // for the returned matrix. Memory will be cleaned up with the arena allocator.
1416  mdivide_left_tri_vd_vari<TriView,R1,C1,R2,C2> *baseVari = new mdivide_left_tri_vd_vari<TriView,R1,C1,R2,C2>(A,b);
1417 
1418  size_t i,j;
1419 
1420  for (i = 0; i < res.rows(); i++)
1421  for (j = 0; j < res.cols(); j++)
1422  res(i,j).vi_ = baseVari->_variRefC(i,j);
1423 
1424  return res;
1425  }
1426 
1427  template<int R,int C>
1428  class determinant_vari : public vari {
1429  Eigen::Matrix<double,R,C> _A;
1430  Eigen::Matrix<vari*,R,C> _adjARef;
1431  public:
1432  determinant_vari(const Eigen::Matrix<var,R,C> &A)
1433  : vari(determinant_vari_calc(A)), _A(A.rows(),A.cols()), _adjARef(A.rows(),A.cols())
1434  {
1435  size_t i,j;
1436  for (j = 0; j < A.cols(); j++) {
1437  for (i = 0; i < A.rows(); i++) {
1438  _A(i,j) = A(i,j).val();
1439  _adjARef(i,j) = A(i,j).vi_;
1440  }
1441  }
1442  }
1443  static
1444  double determinant_vari_calc(const Eigen::Matrix<var,R,C> &A)
1445  {
1446  Eigen::Matrix<double,R,C> Ad(A.rows(),A.cols());
1447  size_t i,j;
1448  for (j = 0; j < A.cols(); j++)
1449  for (i = 0; i < A.rows(); i++)
1450  Ad(i,j) = A(i,j).val();
1451  return Ad.determinant();
1452  }
1453  virtual void chain() {
1454  Eigen::Matrix<double,R,C> adjA(_A.rows(),_A.cols());
1455  size_t i,j;
1456  adjA = (adj_ * val_) * _A.inverse().transpose();
1457  for (j = 0; j < _adjARef.cols(); j++) {
1458  for (i = 0; i < _adjARef.rows(); i++) {
1459  _adjARef(i,j)->adj_ += adjA(i,j);
1460  }
1461  }
1462  }
1463  };
1464 
1465  template <int R, int C>
1466  inline var determinant(const Eigen::Matrix<var,R,C>& m) {
1467  stan::math::validate_square(m,"determinant");
1468  return var(new determinant_vari<R,C>(m));
1469  }
1470 
1471  template<int R,int C>
1472  class log_determinant_vari : public vari {
1473  Eigen::Matrix<double,R,C> _A;
1474  Eigen::Matrix<vari*,R,C> _adjARef;
1475  public:
1476  log_determinant_vari(const Eigen::Matrix<var,R,C> &A)
1477  : vari(log_determinant_vari_calc(A)), _A(A.rows(),A.cols()), _adjARef(A.rows(),A.cols())
1478  {
1479  size_t i,j;
1480  for (j = 0; j < A.cols(); j++) {
1481  for (i = 0; i < A.rows(); i++) {
1482  _A(i,j) = A(i,j).val();
1483  _adjARef(i,j) = A(i,j).vi_;
1484  }
1485  }
1486  }
1487  static
1488  double log_determinant_vari_calc(const Eigen::Matrix<var,R,C> &A)
1489  {
1490  Eigen::Matrix<double,R,C> Ad(A.rows(),A.cols());
1491  size_t i,j;
1492  for (j = 0; j < A.cols(); j++)
1493  for (i = 0; i < A.rows(); i++)
1494  Ad(i,j) = A(i,j).val();
1495  return Ad.colPivHouseholderQr().logAbsDeterminant();
1496  }
1497  virtual void chain() {
1498  Eigen::Matrix<double,R,C> adjA(_A.rows(),_A.cols());
1499  size_t i,j;
1500  adjA = adj_ * _A.inverse().transpose();
1501  for (j = 0; j < _adjARef.cols(); j++) {
1502  for (i = 0; i < _adjARef.rows(); i++) {
1503  _adjARef(i,j)->adj_ += adjA(i,j);
1504  }
1505  }
1506  }
1507  };
1508 
1509  template <int R, int C>
1510  inline var log_determinant(const Eigen::Matrix<var,R,C>& m) {
1511  stan::math::validate_square(m,"log_determinant");
1512  return var(new log_determinant_vari<R,C>(m));
1513  }
1514 
1522  inline double
1523  divide(double x, double y) {
1524  return x / y;
1525  }
1526  template <typename T1, typename T2>
1527  inline var
1528  divide(const T1& v, const T2& c) {
1529  return to_var(v) / to_var(c);
1530  }
1538  template <typename T1, typename T2, int R, int C>
1539  inline Eigen::Matrix<var,R,C>
1540  divide(const Eigen::Matrix<T1, R,C>& v, const T2& c) {
1541  return to_var(v) / to_var(c);
1542  }
1543 
1544 
1545 
1552  template <typename T1, typename T2>
1553  inline
1554  typename boost::math::tools::promote_args<T1,T2>::type
1555  multiply(const T1& v, const T2& c) {
1556  return v * c;
1557  }
1558 
1565  template<typename T1,typename T2,int R2,int C2>
1566  inline Eigen::Matrix<var,R2,C2> multiply(const T1& c,
1567  const Eigen::Matrix<T2, R2, C2>& m) {
1568  // FIXME: pull out to eliminate overpromotion of one side
1569  // move to matrix.hpp w. promotion?
1570  return to_var(m) * to_var(c);
1571  }
1572 
1579  template<typename T1,int R1,int C1,typename T2>
1580  inline Eigen::Matrix<var,R1,C1> multiply(const Eigen::Matrix<T1, R1, C1>& m,
1581  const T2& c) {
1582  return to_var(m) * to_var(c);
1583  }
1584 
1595  template<int R1,int C1,int R2,int C2>
1596  inline Eigen::Matrix<var,R1,C2> multiply(const Eigen::Matrix<var,R1,C1>& m1,
1597  const Eigen::Matrix<var,R2,C2>& m2) {
1598  stan::math::validate_multiplicable(m1,m2,"multiply");
1599  Eigen::Matrix<var,R1,C2> result(m1.rows(),m2.cols());
1600  for (int i = 0; i < m1.rows(); i++) {
1601  typename Eigen::Matrix<var,R1,C1>::ConstRowXpr crow(m1.row(i));
1602  for (int j = 0; j < m2.cols(); j++) {
1603  typename Eigen::Matrix<var,R2,C2>::ConstColXpr ccol(m2.col(j));
1604  if (j == 0) {
1605  if (i == 0) {
1606  result(i,j) = var(new dot_product_vv_vari(crow,ccol));
1607  }
1608  else {
1609  dot_product_vv_vari *v2 = static_cast<dot_product_vv_vari*>(result(0,j).vi_);
1610  result(i,j) = var(new dot_product_vv_vari(crow,ccol,NULL,v2));
1611  }
1612  }
1613  else {
1614  if (i == 0) {
1615  dot_product_vv_vari *v1 = static_cast<dot_product_vv_vari*>(result(i,0).vi_);
1616  result(i,j) = var(new dot_product_vv_vari(crow,ccol,v1));
1617  }
1618  else /* if (i != 0 && j != 0) */ {
1619  dot_product_vv_vari *v1 = static_cast<dot_product_vv_vari*>(result(i,0).vi_);
1620  dot_product_vv_vari *v2 = static_cast<dot_product_vv_vari*>(result(0,j).vi_);
1621  result(i,j) = var(new dot_product_vv_vari(crow,ccol,v1,v2));
1622  }
1623  }
1624  }
1625  }
1626  return result;
1627  }
1628 
1639  template<int R1,int C1,int R2,int C2>
1640  inline Eigen::Matrix<var,R1,C2> multiply(const Eigen::Matrix<double,R1,C1>& m1,
1641  const Eigen::Matrix<var,R2,C2>& m2) {
1642  stan::math::validate_multiplicable(m1,m2,"multiply");
1643  Eigen::Matrix<var,R1,C2> result(m1.rows(),m2.cols());
1644  for (int i = 0; i < m1.rows(); i++) {
1645  typename Eigen::Matrix<double,R1,C1>::ConstRowXpr crow(m1.row(i));
1646  for (int j = 0; j < m2.cols(); j++) {
1647  typename Eigen::Matrix<var,R2,C2>::ConstColXpr ccol(m2.col(j));
1648 // result(i,j) = dot_product(crow,ccol);
1649  if (j == 0) {
1650  if (i == 0) {
1651  result(i,j) = var(new dot_product_vd_vari(ccol,crow));
1652  }
1653  else {
1654  dot_product_vd_vari *v2 = static_cast<dot_product_vd_vari*>(result(0,j).vi_);
1655  result(i,j) = var(new dot_product_vd_vari(ccol,crow,v2,NULL));
1656  }
1657  }
1658  else {
1659  if (i == 0) {
1660  dot_product_vd_vari *v1 = static_cast<dot_product_vd_vari*>(result(i,0).vi_);
1661  result(i,j) = var(new dot_product_vd_vari(ccol,crow,NULL,v1));
1662  }
1663  else /* if (i != 0 && j != 0) */ {
1664  dot_product_vd_vari *v1 = static_cast<dot_product_vd_vari*>(result(i,0).vi_);
1665  dot_product_vd_vari *v2 = static_cast<dot_product_vd_vari*>(result(0,j).vi_);
1666  result(i,j) = var(new dot_product_vd_vari(ccol,crow,v2,v1));
1667  }
1668  }
1669  }
1670  }
1671  return result;
1672  }
1673 
1684  template<int R1,int C1,int R2,int C2>
1685  inline Eigen::Matrix<var,R1,C2> multiply(const Eigen::Matrix<var,R1,C1>& m1,
1686  const Eigen::Matrix<double,R2,C2>& m2) {
1687  stan::math::validate_multiplicable(m1,m2,"multiply");
1688  Eigen::Matrix<var,R1,C2> result(m1.rows(),m2.cols());
1689  for (int i = 0; i < m1.rows(); i++) {
1690  typename Eigen::Matrix<var,R1,C1>::ConstRowXpr crow(m1.row(i));
1691  for (int j = 0; j < m2.cols(); j++) {
1692  typename Eigen::Matrix<double,R2,C2>::ConstColXpr ccol(m2.col(j));
1693 // result(i,j) = dot_product(crow,ccol);
1694  if (j == 0) {
1695  if (i == 0) {
1696  result(i,j) = var(new dot_product_vd_vari(crow,ccol));
1697  }
1698  else {
1699  dot_product_vd_vari *v2 = static_cast<dot_product_vd_vari*>(result(0,j).vi_);
1700  result(i,j) = var(new dot_product_vd_vari(crow,ccol,NULL,v2));
1701  }
1702  }
1703  else {
1704  if (i == 0) {
1705  dot_product_vd_vari *v1 = static_cast<dot_product_vd_vari*>(result(i,0).vi_);
1706  result(i,j) = var(new dot_product_vd_vari(crow,ccol,v1,NULL));
1707  }
1708  else /* if (i != 0 && j != 0) */ {
1709  dot_product_vd_vari *v1 = static_cast<dot_product_vd_vari*>(result(i,0).vi_);
1710  dot_product_vd_vari *v2 = static_cast<dot_product_vd_vari*>(result(0,j).vi_);
1711  result(i,j) = var(new dot_product_vd_vari(crow,ccol,v1,v2));
1712  }
1713  }
1714  }
1715  }
1716  return result;
1717  }
1718 
1728  template <int C1,int R2>
1729  inline var multiply(const Eigen::Matrix<var, 1, C1>& rv,
1730  const Eigen::Matrix<var, R2, 1>& v) {
1731  if (rv.size() != v.size())
1732  throw std::domain_error("row vector and vector must be same length in multiply");
1733  return dot_product(rv, v);
1734  }
1744  template <int C1,int R2>
1745  inline var multiply(const Eigen::Matrix<double, 1, C1>& rv,
1746  const Eigen::Matrix<var, R2, 1>& v) {
1747  stan::math::validate_multiplicable(rv,v,"multiply");
1748  return dot_product(rv, v);
1749  }
1759  template <int C1,int R2>
1760  inline var multiply(const Eigen::Matrix<var, 1, C1>& rv,
1761  const Eigen::Matrix<double, R2, 1>& v) {
1762  stan::math::validate_multiplicable(rv,v,"multiply");
1763  return dot_product(rv, v);
1764  }
1765 
1766  inline matrix_v
1767  multiply_lower_tri_self_transpose(const matrix_v& L) {
1768 // stan::math::validate_square(L,"multiply_lower_tri_self_transpose");
1769  int K = L.rows();
1770  int J = L.cols();
1771  matrix_v LLt(K,K);
1772  if (K == 0) return LLt;
1773  // if (K == 1) {
1774  // LLt(0,0) = L(0,0) * L(0,0);
1775  // return LLt;
1776  // }
1777  int Knz;
1778  if (K >= J)
1779  Knz = (K-J)*J + (J * (J + 1)) / 2;
1780  else // if (K < J)
1781  Knz = (K * (K + 1)) / 2;
1782  vari** vs = (vari**)memalloc_.alloc( Knz * sizeof(vari*) );
1783  int pos = 0;
1784  for (int m = 0; m < K; ++m)
1785  for (int n = 0; n < ((J < (m+1))?J:(m+1)); ++n) {
1786  vs[pos++] = L(m,n).vi_;
1787  }
1788  for (int m = 0, mpos=0; m < K; ++m, mpos += (J < m)?J:m) {
1789  LLt(m,m) = var(new dot_self_vari(vs + mpos, (J < (m+1))?J:(m+1)));
1790  for (int n = 0, npos = 0; n < m; ++n, npos += (J < n)?J:n) {
1791  LLt(m,n) = LLt(n,m) = var(new dot_product_vv_vari(vs + mpos, vs + npos, (J < (n+1))?J:(n+1)));
1792  }
1793  }
1794  return LLt;
1795  }
1796 
1803  inline matrix_v
1804  tcrossprod(const matrix_v& M) {
1805  if (M.rows() == 0)
1806  return matrix_v(0,0);
1807  if (M.rows() == 1)
1808  return M * M.transpose();
1809 
1810  // WAS JUST THIS
1811  // matrix_v result(M.rows(),M.rows());
1812  // return result.setZero().selfadjointView<Eigen::Upper>().rankUpdate(M);
1813 
1814  matrix_v MMt(M.rows(),M.rows());
1815 
1816  vari** vs
1817  = (vari**)memalloc_.alloc((M.rows() * M.cols() ) * sizeof(vari*));
1818  int pos = 0;
1819  for (int m = 0; m < M.rows(); ++m)
1820  for (int n = 0; n < M.cols(); ++n)
1821  vs[pos++] = M(m,n).vi_;
1822  for (int m = 0; m < M.rows(); ++m)
1823  MMt(m,m) = var(new dot_self_vari(vs + m * M.cols(),M.cols()));
1824  for (int m = 0; m < M.rows(); ++m) {
1825  for (int n = 0; n < m; ++n) {
1826  MMt(m,n) = var(new dot_product_vv_vari(vs + m * M.cols(),
1827  vs + n * M.cols(),
1828  M.cols()));
1829  MMt(n,m) = MMt(m,n);
1830  }
1831  }
1832  return MMt;
1833  }
1834 
1841  inline matrix_v
1842  crossprod(const matrix_v& M) {
1843  return tcrossprod(M.transpose());
1844  }
1845 
1846 
1847  // FIXME: double val?
1848  inline void assign_to_var(stan::agrad::var& var, const double& val) {
1849  var = val;
1850  }
1851  inline void assign_to_var(stan::agrad::var& var, const stan::agrad::var& val) {
1852  var = val;
1853  }
1854  // FIXME: int val?
1855  inline void assign_to_var(int& n_lhs, const int& n_rhs) {
1856  n_lhs = n_rhs; // FIXME: no call -- just filler to instantiate
1857  }
1858  // FIXME: double val?
1859  inline void assign_to_var(double& n_lhs, const double& n_rhs) {
1860  n_lhs = n_rhs; // FIXME: no call -- just filler to instantiate
1861  }
1862 
1863  template <typename LHS, typename RHS>
1864  inline void assign_to_var(std::vector<LHS>& x, const std::vector<RHS>& y) {
1865  stan::math::validate_matching_sizes(x,y,"assign_to_var");
1866  for (size_t i = 0; i < x.size(); ++i)
1867  assign_to_var(x[i],y[i]);
1868  }
1869  template <typename LHS, typename RHS, int R, int C>
1870  inline void assign_to_var(Eigen::Matrix<LHS,R,C>& x,
1871  const Eigen::Matrix<RHS,R,C>& y) {
1872  stan::math::validate_matching_sizes(x,y,"assign_to_var");
1873  typedef typename Eigen::Matrix<LHS,R,C>::size_type size_type;
1874  for (size_type n = 0; n < x.cols(); ++n)
1875  for (size_type m = 0; m < x.rows(); ++m)
1876  assign_to_var(x(m,n),y(m,n));
1877  }
1878 
1879  template <typename LHS, typename RHS>
1880  struct needs_promotion {
1881  enum { value = ( is_constant_struct<RHS>::value
1882  && !is_constant_struct<LHS>::value) };
1883  };
1884 
1885  template <bool PromoteRHS, typename LHS, typename RHS>
1886  struct assigner {
1887  static inline void assign(LHS& /*var*/, const RHS& /*val*/) {
1888  throw std::domain_error("should not call base class of assigner");
1889  }
1890  };
1891 
1892  template <typename LHS, typename RHS>
1893  struct assigner<false,LHS,RHS> {
1894  static inline void assign(LHS& var, const RHS& val) {
1895  var = val; // no promotion of RHS
1896  }
1897  };
1898 
1899  template <typename LHS, typename RHS>
1900  struct assigner<true,LHS,RHS> {
1901  static inline void assign(LHS& var, const RHS& val) {
1902  assign_to_var(var,val); // promote RHS
1903  }
1904  };
1905 
1906 
1907  template <typename LHS, typename RHS>
1908  inline void assign(LHS& var, const RHS& val) {
1909  assigner<needs_promotion<LHS,RHS>::value, LHS, RHS>::assign(var,val);
1910  }
1911 
1912  void stan_print(std::ostream* o, const var& x) {
1913  *o << x.val();
1914  }
1915 
1916  }
1917 }
1918 
1919 
1920 #endif
1921 
     [ Stan Home Page ] © 2011–2013, Stan Development Team.