1 //===================================================== 2 // Copyright (C) 2008 Gael Guennebaud <g.gael (at) free.fr> 3 //===================================================== 4 // 5 // This program is free software; you can redistribute it and/or 6 // modify it under the terms of the GNU General Public License 7 // as published by the Free Software Foundation; either version 2 8 // of the License, or (at your option) any later version. 9 // 10 // This program is distributed in the hope that it will be useful, 11 // but WITHOUT ANY WARRANTY; without even the implied warranty of 12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13 // GNU General Public License for more details. 14 // You should have received a copy of the GNU General Public License 15 // along with this program; if not, write to the Free Software 16 // Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. 17 // 18 #ifndef EIGEN2_INTERFACE_HH 19 #define EIGEN2_INTERFACE_HH 20 // #include <cblas.h> 21 #include <Eigen/Core> 22 #include <Eigen/Cholesky> 23 #include <Eigen/LU> 24 #include <Eigen/QR> 25 #include <vector> 26 #include "btl.hh" 27 28 using namespace Eigen; 29 30 template<class real, int SIZE=Dynamic> 31 class eigen2_interface 32 { 33 34 public : 35 36 enum {IsFixedSize = (SIZE!=Dynamic)}; 37 38 typedef real real_type; 39 40 typedef std::vector<real> stl_vector; 41 typedef std::vector<stl_vector> stl_matrix; 42 43 typedef Eigen::Matrix<real,SIZE,SIZE> gene_matrix; 44 typedef Eigen::Matrix<real,SIZE,1> gene_vector; 45 46 static inline std::string name( void ) 47 { 48 #if defined(EIGEN_VECTORIZE_SSE) 49 if (SIZE==Dynamic) return "eigen2"; else return "tiny_eigen2"; 50 #elif defined(EIGEN_VECTORIZE_ALTIVEC) || defined(EIGEN_VECTORIZE_VSX) 51 if (SIZE==Dynamic) return "eigen2"; else return "tiny_eigen2"; 52 #else 53 if (SIZE==Dynamic) return "eigen2_novec"; else return "tiny_eigen2_novec"; 54 #endif 55 } 56 57 static void free_matrix(gene_matrix & A, int N) {} 58 59 static void free_vector(gene_vector & B) {} 60 61 static BTL_DONT_INLINE void matrix_from_stl(gene_matrix & A, stl_matrix & A_stl){ 62 A.resize(A_stl[0].size(), A_stl.size()); 63 64 for (int j=0; j<A_stl.size() ; j++){ 65 for (int i=0; i<A_stl[j].size() ; i++){ 66 A.coeffRef(i,j) = A_stl[j][i]; 67 } 68 } 69 } 70 71 static BTL_DONT_INLINE void vector_from_stl(gene_vector & B, stl_vector & B_stl){ 72 B.resize(B_stl.size(),1); 73 74 for (int i=0; i<B_stl.size() ; i++){ 75 B.coeffRef(i) = B_stl[i]; 76 } 77 } 78 79 static BTL_DONT_INLINE void vector_to_stl(gene_vector & B, stl_vector & B_stl){ 80 for (int i=0; i<B_stl.size() ; i++){ 81 B_stl[i] = B.coeff(i); 82 } 83 } 84 85 static BTL_DONT_INLINE void matrix_to_stl(gene_matrix & A, stl_matrix & A_stl){ 86 int N=A_stl.size(); 87 88 for (int j=0;j<N;j++){ 89 A_stl[j].resize(N); 90 for (int i=0;i<N;i++){ 91 A_stl[j][i] = A.coeff(i,j); 92 } 93 } 94 } 95 96 static inline void matrix_matrix_product(const gene_matrix & A, const gene_matrix & B, gene_matrix & X, int N){ 97 X = (A*B).lazy(); 98 } 99 100 static inline void transposed_matrix_matrix_product(const gene_matrix & A, const gene_matrix & B, gene_matrix & X, int N){ 101 X = (A.transpose()*B.transpose()).lazy(); 102 } 103 104 static inline void ata_product(const gene_matrix & A, gene_matrix & X, int N){ 105 X = (A.transpose()*A).lazy(); 106 } 107 108 static inline void aat_product(const gene_matrix & A, gene_matrix & X, int N){ 109 X = (A*A.transpose()).lazy(); 110 } 111 112 static inline void matrix_vector_product(const gene_matrix & A, const gene_vector & B, gene_vector & X, int N){ 113 X = (A*B)/*.lazy()*/; 114 } 115 116 static inline void atv_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N){ 117 X = (A.transpose()*B)/*.lazy()*/; 118 } 119 120 static inline void axpy(real coef, const gene_vector & X, gene_vector & Y, int N){ 121 Y += coef * X; 122 } 123 124 static inline void axpby(real a, const gene_vector & X, real b, gene_vector & Y, int N){ 125 Y = a*X + b*Y; 126 } 127 128 static inline void copy_matrix(const gene_matrix & source, gene_matrix & cible, int N){ 129 cible = source; 130 } 131 132 static inline void copy_vector(const gene_vector & source, gene_vector & cible, int N){ 133 cible = source; 134 } 135 136 static inline void trisolve_lower(const gene_matrix & L, const gene_vector& B, gene_vector& X, int N){ 137 X = L.template marked<LowerTriangular>().solveTriangular(B); 138 } 139 140 static inline void trisolve_lower_matrix(const gene_matrix & L, const gene_matrix& B, gene_matrix& X, int N){ 141 X = L.template marked<LowerTriangular>().solveTriangular(B); 142 } 143 144 static inline void cholesky(const gene_matrix & X, gene_matrix & C, int N){ 145 C = X.llt().matrixL(); 146 // C = X; 147 // Cholesky<gene_matrix>::computeInPlace(C); 148 // Cholesky<gene_matrix>::computeInPlaceBlock(C); 149 } 150 151 static inline void lu_decomp(const gene_matrix & X, gene_matrix & C, int N){ 152 C = X.lu().matrixLU(); 153 // C = X.inverse(); 154 } 155 156 static inline void tridiagonalization(const gene_matrix & X, gene_matrix & C, int N){ 157 C = Tridiagonalization<gene_matrix>(X).packedMatrix(); 158 } 159 160 static inline void hessenberg(const gene_matrix & X, gene_matrix & C, int N){ 161 C = HessenbergDecomposition<gene_matrix>(X).packedMatrix(); 162 } 163 164 165 166 }; 167 168 #endif 169
