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