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