1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath//===================================================== 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.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 GMM_INTERFACE_HH 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define GMM_INTERFACE_HH 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include <gmm/gmm.h> 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include <vector> 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathusing namespace gmm; 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<class real> 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass gmm_interface { 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathpublic : 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef real real_type ; 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef std::vector<real> stl_vector; 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef std::vector<stl_vector > stl_matrix; 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef gmm::dense_matrix<real> gene_matrix; 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef stl_vector gene_vector; 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static inline std::string name( void ) 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return "gmm"; 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static void free_matrix(gene_matrix & A, int N){ 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return ; 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static void free_vector(gene_vector & B){ 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return ; 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static inline void matrix_from_stl(gene_matrix & A, stl_matrix & A_stl){ 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath A.resize(A_stl[0].size(),A_stl.size()); 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int j=0; j<A_stl.size() ; j++){ 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int i=0; i<A_stl[j].size() ; i++){ 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath A(i,j) = A_stl[j][i]; 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static inline void vector_from_stl(gene_vector & B, stl_vector & B_stl){ 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath B = B_stl; 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static inline void vector_to_stl(gene_vector & B, stl_vector & B_stl){ 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath B_stl = B; 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static inline void matrix_to_stl(gene_matrix & A, stl_matrix & A_stl){ 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int N=A_stl.size(); 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int j=0;j<N;j++){ 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath A_stl[j].resize(N); 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int i=0;i<N;i++){ 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath A_stl[j][i] = A(i,j); 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static inline void matrix_matrix_product(const gene_matrix & A, const gene_matrix & B, gene_matrix & X, int N){ 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath gmm::mult(A,B, X); 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static inline void transposed_matrix_matrix_product(const gene_matrix & A, const gene_matrix & B, gene_matrix & X, int N){ 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath gmm::mult(gmm::transposed(A),gmm::transposed(B), X); 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static inline void ata_product(const gene_matrix & A, gene_matrix & X, int N){ 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath gmm::mult(gmm::transposed(A),A, X); 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static inline void aat_product(const gene_matrix & A, gene_matrix & X, int N){ 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath gmm::mult(A,gmm::transposed(A), X); 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static inline void matrix_vector_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N){ 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath gmm::mult(A,B,X); 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static inline void atv_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N){ 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath gmm::mult(gmm::transposed(A),B,X); 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static inline void axpy(const real coef, const gene_vector & X, gene_vector & Y, int N){ 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath gmm::add(gmm::scaled(X,coef), Y); 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static inline void axpby(real a, const gene_vector & X, real b, gene_vector & Y, int N){ 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath gmm::add(gmm::scaled(X,a), gmm::scaled(Y,b), Y); 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static inline void copy_matrix(const gene_matrix & source, gene_matrix & cible, int N){ 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath gmm::copy(source,cible); 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static inline void copy_vector(const gene_vector & source, gene_vector & cible, int N){ 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath gmm::copy(source,cible); 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static inline void trisolve_lower(const gene_matrix & L, const gene_vector& B, gene_vector & X, int N){ 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath gmm::copy(B,X); 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath gmm::lower_tri_solve(L, X, false); 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static inline void partial_lu_decomp(const gene_matrix & X, gene_matrix & R, int N){ 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath gmm::copy(X,R); 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::vector<int> ipvt(N); 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath gmm::lu_factor(R, ipvt); 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static inline void hessenberg(const gene_matrix & X, gene_matrix & R, int N){ 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath gmm::copy(X,R); 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath gmm::Hessenberg_reduction(R,X,false); 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static inline void tridiagonalization(const gene_matrix & X, gene_matrix & R, int N){ 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath gmm::copy(X,R); 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath gmm::Householder_tridiagonalization(R,X,false); 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif 145