1 //===================================================== 2 // File : ublas_interface.hh 3 // Author : L. Plagne <laurent.plagne (at) edf.fr)> 4 // Copyright (C) EDF R&D, lun sep 30 14:23:27 CEST 2002 5 //===================================================== 6 // 7 // This program is free software; you can redistribute it and/or 8 // modify it under the terms of the GNU General Public License 9 // as published by the Free Software Foundation; either version 2 10 // of the License, or (at your option) any later version. 11 // 12 // This program is distributed in the hope that it will be useful, 13 // but WITHOUT ANY WARRANTY; without even the implied warranty of 14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15 // GNU General Public License for more details. 16 // You should have received a copy of the GNU General Public License 17 // along with this program; if not, write to the Free Software 18 // Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. 19 // 20 #ifndef UBLAS_INTERFACE_HH 21 #define UBLAS_INTERFACE_HH 22 23 #include <boost/numeric/ublas/vector.hpp> 24 #include <boost/numeric/ublas/matrix.hpp> 25 #include <boost/numeric/ublas/io.hpp> 26 #include <boost/numeric/ublas/triangular.hpp> 27 28 using namespace boost::numeric; 29 30 template <class real> 31 class ublas_interface{ 32 33 public : 34 35 typedef real real_type ; 36 37 typedef std::vector<real> stl_vector; 38 typedef std::vector<stl_vector> stl_matrix; 39 40 typedef typename boost::numeric::ublas::matrix<real,boost::numeric::ublas::column_major> gene_matrix; 41 typedef typename boost::numeric::ublas::vector<real> gene_vector; 42 43 static inline std::string name( void ) { return "ublas"; } 44 45 static void free_matrix(gene_matrix & A, int N) {} 46 47 static void free_vector(gene_vector & B) {} 48 49 static inline void matrix_from_stl(gene_matrix & A, stl_matrix & A_stl){ 50 A.resize(A_stl.size(),A_stl[0].size()); 51 for (int j=0; j<A_stl.size() ; j++) 52 for (int i=0; i<A_stl[j].size() ; i++) 53 A(i,j)=A_stl[j][i]; 54 } 55 56 static inline void vector_from_stl(gene_vector & B, stl_vector & B_stl){ 57 B.resize(B_stl.size()); 58 for (int i=0; i<B_stl.size() ; i++) 59 B(i)=B_stl[i]; 60 } 61 62 static inline void vector_to_stl(gene_vector & B, stl_vector & B_stl){ 63 for (int i=0; i<B_stl.size() ; i++) 64 B_stl[i]=B(i); 65 } 66 67 static inline void matrix_to_stl(gene_matrix & A, stl_matrix & A_stl){ 68 int N=A_stl.size(); 69 for (int j=0;j<N;j++) 70 { 71 A_stl[j].resize(N); 72 for (int i=0;i<N;i++) 73 A_stl[j][i]=A(i,j); 74 } 75 } 76 77 static inline void copy_vector(const gene_vector & source, gene_vector & cible, int N){ 78 for (int i=0;i<N;i++){ 79 cible(i) = source(i); 80 } 81 } 82 83 static inline void copy_matrix(const gene_matrix & source, gene_matrix & cible, int N){ 84 for (int i=0;i<N;i++){ 85 for (int j=0;j<N;j++){ 86 cible(i,j) = source(i,j); 87 } 88 } 89 } 90 91 static inline void matrix_vector_product_slow(gene_matrix & A, gene_vector & B, gene_vector & X, int N){ 92 X = prod(A,B); 93 } 94 95 static inline void matrix_matrix_product_slow(gene_matrix & A, gene_matrix & B, gene_matrix & X, int N){ 96 X = prod(A,B); 97 } 98 99 static inline void axpy_slow(const real coef, const gene_vector & X, gene_vector & Y, int N){ 100 Y+=coef*X; 101 } 102 103 // alias free assignements 104 105 static inline void matrix_vector_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N){ 106 X.assign(prod(A,B)); 107 } 108 109 static inline void atv_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N){ 110 X.assign(prod(trans(A),B)); 111 } 112 113 static inline void matrix_matrix_product(gene_matrix & A, gene_matrix & B, gene_matrix & X, int N){ 114 X.assign(prod(A,B)); 115 } 116 117 static inline void axpy(const real coef, const gene_vector & X, gene_vector & Y, int N){ 118 Y.plus_assign(coef*X); 119 } 120 121 static inline void axpby(real a, const gene_vector & X, real b, gene_vector & Y, int N){ 122 Y = a*X + b*Y; 123 } 124 125 static inline void ata_product(gene_matrix & A, gene_matrix & X, int N){ 126 // X = prod(trans(A),A); 127 X.assign(prod(trans(A),A)); 128 } 129 130 static inline void aat_product(gene_matrix & A, gene_matrix & X, int N){ 131 // X = prod(A,trans(A)); 132 X.assign(prod(A,trans(A))); 133 } 134 135 static inline void trisolve_lower(const gene_matrix & L, const gene_vector& B, gene_vector & X, int N){ 136 X = solve(L, B, ublas::lower_tag ()); 137 } 138 139 }; 140 141 #endif 142
