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