HMSBEAGLE  1.0.0
EigenDecomposition.h
1 /*
2  * EigenDecomposition.h
3  *
4  * Created on: Sep 24, 2009
5  * Author: msuchard
6  */
7 
8 #ifndef EIGENDECOMPOSITION_H_
9 #define EIGENDECOMPOSITION_H_
10 
11 #include <cstdio>
12 #include <cstdlib>
13 #include <iostream>
14 #include <cstring>
15 #include <cmath>
16 #include <cassert>
17 #include <vector>
18 
19 #define BEAGLE_CPU_EIGEN_GENERIC REALTYPE, T_PAD
20 #define BEAGLE_CPU_EIGEN_TEMPLATE template <typename REALTYPE, int T_PAD>
21 
22 namespace beagle {
23 namespace cpu {
24 
25 BEAGLE_CPU_EIGEN_TEMPLATE
27 
28 protected:
29  REALTYPE** gEigenValues;
30  int kStateCount;
31  int kEigenDecompCount;
32  int kCategoryCount;
33  long kFlags;
34  REALTYPE* matrixTmp;
35  REALTYPE* firstDerivTmp;
36  REALTYPE* secondDerivTmp;
37 
38 public:
39  EigenDecomposition(int decompositionCount,
40  int stateCount,
41  int categoryCount,
42  long flags)
43  {
44 
45  kEigenDecompCount = decompositionCount;
46  kStateCount = stateCount;
47  kCategoryCount = categoryCount;
48  kFlags = flags;
49  };
50 
51  virtual ~EigenDecomposition() {};
52 
53  // sets the Eigen decomposition for a given matrix
54  //
55  // matrixIndex the matrix index to update
56  // eigenVectors an array containing the Eigen Vectors
57  // inverseEigenVectors an array containing the inverse Eigen Vectors
58  // eigenValues an array containing the Eigen Values
59  virtual void setEigenDecomposition(int eigenIndex,
60  const double* inEigenVectors,
61  const double* inInverseEigenVectors,
62  const double* inEigenValues) = 0;
63 
64  // calculate a transition probability matrices for a given list of node. This will
65  // calculate for all categories (and all matrices if more than one is being used).
66  //
67  // nodeIndices an array of node indices that require transition probability matrices
68  // edgeLengths an array of expected lengths in substitutions per site
69  // count the number of elements in the above arrays
70  virtual void updateTransitionMatrices(int eigenIndex,
71  const int* probabilityIndices,
72  const int* firstDerivativeIndices,
73  const int* secondDerivativeIndices,
74  const double* edgeLengths,
75  const double* categoryRates,
76  REALTYPE** transitionMatrices,
77  int count) = 0;
78 
79 };
80 
81 }
82 }
83 
84 #endif /* EIGENDECOMPOSITION_H_ */