00001 // This matrix class is a C++ wrapper for the GNU Scientific Library 00002 // Copyright (C) ULP-IPB Strasbourg 00003 00004 // This program is free software; you can redistribute it and/or modify 00005 // it under the terms of the GNU General Public License as published by 00006 // the Free Software Foundation; either version 2 of the License, or 00007 // (at your option) any later version. 00008 00009 // This program is distributed in the hope that it will be useful, 00010 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00011 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00012 // GNU General Public License for more details. 00013 00014 // You should have received a copy of the GNU General Public License 00015 // along with this program; if not, write to the Free Software 00016 // Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. 00017 00018 #ifndef _vector_int_h 00019 #define _vector_int_h 00020 00021 #include <iostream.h> 00022 #include <gsl/gsl_math.h> 00023 #include <gsl/gsl_vector_int.h> 00024 #include <gsl/gsl_blas.h> 00025 #include <gslwrap/vector_double.h> 00026 00027 namespace gsl 00028 { 00029 00030 class vector_int_view; 00031 00032 class vector_int 00033 { 00034 protected: 00035 gsl_vector_int *gsldata; 00036 void free(){if(gsldata) gsl_vector_int_free(gsldata);gsldata=NULL;} 00037 void alloc(size_t n) {gsldata=gsl_vector_int_alloc(n);} 00038 void calloc(size_t n){gsldata=gsl_vector_int_calloc(n);} 00039 public: 00040 vector_int() : gsldata(NULL) {;} 00041 vector_int( const vector_int &other ):gsldata(NULL) {copy(other);} 00042 template<class oclass> 00043 vector_int( const oclass &other ):gsldata(NULL) {copy(other);} 00044 ~vector_int(){free();} 00045 vector_int(const size_t& n,bool clear=true) 00046 { 00047 if(clear){this->calloc(n);} 00048 else {this->alloc(n);} 00049 } 00050 vector_int(const int& n,bool clear=true) 00051 { 00052 if(clear){this->calloc(n);} 00053 else {this->alloc(n);} 00054 } 00055 00056 void resize(size_t n); 00057 00058 template <class oclass> 00059 void copy(const oclass &other) 00060 { 00061 resize(other.size()); 00062 for (size_t i=0;i<size();i++) 00063 { 00064 gsl_vector_int_set(gsldata, i, (int)other[i]); 00065 } 00066 } 00067 void copy(const vector_int& other); 00068 // void clone(vector_int& other); 00069 00070 size_t size() const {if (!gsldata) {cout << "vector_int::size vector not initialized" << endl; exit(-1);}return gsldata->size;} 00071 00073 gsl_vector_int *gslobj() {if (!gsldata){cout << "vector_int::gslobj ERROR, data not initialized!! " << endl; exit(-1);}return gsldata;} 00074 const gsl_vector_int *gslobj() const {if (!gsldata){cout << "vector_int::gslobj ERROR, data not initialized!! " << endl; exit(-1);}return gsldata;} 00075 00076 00077 static vector_int_view create_vector_view( const gsl_vector_int_view &other ); 00078 00079 // ********Accessing vector elements 00080 00081 // Unlike FORTRAN compilers, C compilers do not usually provide support for range checking of vectors and matrices (2). However, the functions gsl_vector_int_get and gsl_vector_int_set can perform range checking for you and report an error if you attempt to access elements outside the allowed range. 00082 00083 // The functions for accessing the elements of a vector or matrix are defined in `gsl_vector_int.h' and declared extern inline to eliminate function-call overhead. If necessary you can turn off range checking completely without modifying any source files by recompiling your program with the preprocessor definition GSL_RANGE_CHECK_OFF. Provided your compiler supports inline functions the effect of turning off range checking is to replace calls to gsl_vector_int_get(v,i) by v->data[i*v->stride] and and calls to gsl_vector_int_set(v,i,x) by v->data[i*v->stride]=x. Thus there should be no performance penalty for using the range checking functions when range checking is turned off. 00084 00085 // This function returns the i-th element of a vector v. If i lies outside the allowed range of 0 to n-1 then the error handler is invoked and 0 is returned. 00086 int get(size_t i) const {return gsl_vector_int_get(gsldata,i);} 00087 00088 // This function sets the value of the i-th element of a vector v to x. If i lies outside the allowed range of 0 to n-1 then the error handler is invoked. 00089 void set(size_t i,int x){gsl_vector_int_set(gsldata,i,x);} 00090 00091 // These functions return a pointer to the i-th element of a vector v. If i lies outside the allowed range of 0 to n-1 then the error handler is invoked 00092 int &operator[](size_t i) { return *gsl_vector_int_ptr(gsldata,i);} 00093 const int &operator[](size_t i) const { return *gsl_vector_int_ptr(gsldata,i);} 00094 00095 int &operator()(size_t i) { return *gsl_vector_int_ptr(gsldata,i);} 00096 const int &operator()(size_t i) const { return *gsl_vector_int_ptr(gsldata,i);} 00097 00098 00099 00100 // ***** Initializing vector elements 00101 00102 // This function sets all the elements of the vector v to the value x. 00103 void set_all(int x){gsl_vector_int_set_all (gsldata,x);} 00104 // This function sets all the elements of the vector v to zero. 00105 void set_zero(){gsl_vector_int_set_zero (gsldata);} 00106 00107 // This function makes a basis vector by setting all the elements of the vector v to zero except for the i-th element which is set to one. 00108 int set_basis (size_t i) {return gsl_vector_int_set_basis (gsldata,i);} 00109 00110 // **** Reading and writing vectors 00111 00112 // The library provides functions for reading and writing vectors to a file as binary data or formatted text. 00113 00114 00115 // This function writes the elements of the vector v to the stream stream in binary format. The return value is 0 for success and GSL_EFAILED if there was a problem writing to the file. Since the data is written in the native binary format it may not be portable between different architectures. 00116 int fwrite (FILE * stream) const {return gsl_vector_int_fwrite (stream, gsldata);} 00117 00118 // This function reads into the vector v from the open stream stream in binary format. The vector v must be preallocated with the correct length since the function uses the size of v to determine how many bytes to read. The return value is 0 for success and GSL_EFAILED if there was a problem reading from the file. The data is assumed to have been written in the native binary format on the same architecture. 00119 int fread (FILE * stream) {return gsl_vector_int_fread (stream, gsldata);} 00120 00121 // This function writes the elements of the vector v line-by-line to the stream stream using the format specifier format, which should be one of the %g, %e or %f formats for floating point numbers and %d for integers. The function returns 0 for success and GSL_EFAILED if there was a problem writing to the file. 00122 int fprintf (FILE * stream, const char * format) const {return gsl_vector_int_fprintf (stream, gsldata,format) ;} 00123 00124 // This function reads formatted data from the stream stream into the vector v. The vector v must be preallocated with the correct length since the function uses the size of v to determine how many numbers to read. The function returns 0 for success and GSL_EFAILED if there was a problem reading from the file. 00125 int fscanf (FILE * stream) {return gsl_vector_int_fscanf (stream, gsldata); } 00126 00127 00128 00129 00130 // ******* Vector views 00131 00132 // In addition to creating vectors from slices of blocks it is also possible to slice vectors and create vector views. For example, a subvector of another vector can be described with a view, or two views can be made which provide access to the even and odd elements of a vector. 00133 00134 // A vector view is a temporary object, stored on the stack, which can be used to operate on a subset of vector elements. Vector views can be defined for both constant and non-constant vectors, using separate types that preserve constness. A vector view has the type gsl_vector_int_view and a constant vector view has the type gsl_vector_int_const_view. In both cases the elements of the view can be accessed as a gsl_vector_int using the vector component of the view object. A pointer to a vector of type gsl_vector_int * or const gsl_vector_int * can be obtained by taking the address of this component with the & operator. 00135 00136 // These functions return a vector view of a subvector of another vector v. The start of the new vector is offset by offset elements from the start of the original 00137 // vector. The new vector has n elements. Mathematically, the i-th element of the new vector v' is given by, 00138 00139 // v'(i) = v->data[(offset + i)*v->stride] 00140 00141 // where the index i runs from 0 to n-1. 00142 00143 // The data pointer of the returned vector struct is set to null if the combined parameters (offset,n) overrun the end of the original vector. 00144 00145 // The new vector is only a view of the block underlying the original vector, v. The block containing the elements of v is not owned by the new vector. When the 00146 // new vector goes out of scope the original vector v and its block will continue to exist. The original memory can only be deallocated by freeing the original vector. 00147 // Of course, the original vector should not be deallocated while the new vector is still in use. 00148 00149 // The function gsl_vector_int_const_subvector is equivalent to gsl_vector_int_subvector but can be used for vectors which are declared const. 00150 vector_int_view subvector (size_t offset, size_t n); 00151 const vector_int_view subvector (size_t offset, size_t n) const; 00152 // vector_int_const_view subvector (size_t offset, size_t n) const; 00153 00154 // class view 00155 // { 00156 // gsl_vector_int_view *gsldata; 00157 // public: 00158 // view(); 00159 // }; 00160 // view subvector(size_t offset, size_t n) 00161 // { 00162 // return view(gsl_vector_int_subvector(gsldata,offset,n); 00163 // } 00164 // const view subvector(size_t offset, size_t n) const 00165 // { 00166 // return view(gsl_vector_int_const_subvector(gsldata,offset,n); 00167 // } 00168 00169 00170 00171 // Function: gsl_vector_int gsl_vector_int_subvector_with_stride (gsl_vector_int *v, size_t offset, size_t stride, size_t n) 00172 // Function: gsl_vector_int_const_view gsl_vector_int_const_subvector_with_stride (const gsl_vector_int * v, size_t offset, size_t stride, size_t n) 00173 // These functions return a vector view of a subvector of another vector v with an additional stride argument. The subvector is formed in the same way as for 00174 // gsl_vector_int_subvector but the new vector has n elements with a step-size of stride from one element to the next in the original vector. Mathematically, 00175 // the i-th element of the new vector v' is given by, 00176 00177 // v'(i) = v->data[(offset + i*stride)*v->stride] 00178 00179 // where the index i runs from 0 to n-1. 00180 00181 // Note that subvector views give direct access to the underlying elements of the original vector. For example, the following code will zero the even elements of the 00182 // vector v of length n, while leaving the odd elements untouched, 00183 00184 // gsl_vector_int_view v_even = gsl_vector_int_subvector_with_stride (v, 0, 2, n/2); 00185 // gsl_vector_int_set_zero (&v_even.vector); 00186 00187 // A vector view can be passed to any subroutine which takes a vector argument just as a directly allocated vector would be, using &view.vector. For example, the 00188 // following code computes the norm of odd elements of v using the BLAS routine DNRM2, 00189 00190 // gsl_vector_int_view v_odd = gsl_vector_int_subvector_with_stride (v, 1, 2, n/2); 00191 // double r = gsl_blas_dnrm2 (&v_odd.vector); 00192 00193 // The function gsl_vector_int_const_subvector_with_stride is equivalent to gsl_vector_int_subvector_with_stride but can be used for 00194 // vectors which are declared const. 00195 00196 // Function: gsl_vector_int_view gsl_vector_int_complex_real (gsl_vector_int_complex *v) 00197 // Function: gsl_vector_int_const_view gsl_vector_int_complex_const_real (const gsl_vector_int_complex *v) 00198 // These functions return a vector view of the real parts of the complex vector v. 00199 00200 // The function gsl_vector_int_complex_const_real is equivalent to gsl_vector_int_complex_real but can be used for vectors which are declared 00201 // const. 00202 00203 // Function: gsl_vector_int_view gsl_vector_int_complex_imag (gsl_vector_int_complex *v) 00204 // Function: gsl_vector_int_const_view gsl_vector_int_complex_const_imag (const gsl_vector_int_complex *v) 00205 // These functions return a vector view of the imaginary parts of the complex vector v. 00206 00207 // The function gsl_vector_int_complex_const_imag is equivalent to gsl_vector_int_complex_imag but can be used for vectors which are declared 00208 // const. 00209 00210 // Function: gsl_vector_int_view gsl_vector_int_view_array (double *base, size_t n) 00211 // Function: gsl_vector_int_const_view gsl_vector_int_const_view_array (const double *base, size_t n) 00212 // These functions return a vector view of an array. The start of the new vector is given by base and has n elements. Mathematically, the i-th element of the new 00213 // vector v' is given by, 00214 00215 // v'(i) = base[i] 00216 00217 // where the index i runs from 0 to n-1. 00218 00219 // The array containing the elements of v is not owned by the new vector view. When the view goes out of scope the original array will continue to exist. The 00220 // original memory can only be deallocated by freeing the original pointer base. Of course, the original array should not be deallocated while the view is still in use. 00221 00222 // The function gsl_vector_int_const_view_array is equivalent to gsl_vector_int_view_array but can be used for vectors which are declared const. 00223 00224 // Function: gsl_vector_int_view gsl_vector_int_view_array_with_stride (double * base, size_t stride, size_t n) 00225 // Function: gsl_vector_int_const_view gsl_vector_int_const_view_array_with_stride (const double * base, size_t stride, size_t n) 00226 // These functions return a vector view of an array base with an additional stride argument. The subvector is formed in the same way as for 00227 // gsl_vector_int_view_array but the new vector has n elements with a step-size of stride from one element to the next in the original array. Mathematically, 00228 // the i-th element of the new vector v' is given by, 00229 00230 // v'(i) = base[i*stride] 00231 00232 // where the index i runs from 0 to n-1. 00233 00234 // Note that the view gives direct access to the underlying elements of the original array. A vector view can be passed to any subroutine which takes a vector 00235 // argument just as a directly allocated vector would be, using &view.vector. 00236 00237 // The function gsl_vector_int_const_view_array_with_stride is equivalent to gsl_vector_int_view_array_with_stride but can be used for 00238 // arrays which are declared const. 00239 00240 00241 // ************* Copying vectors 00242 00243 // Common operations on vectors such as addition and multiplication are available in the BLAS part of the library (see section BLAS Support). However, it is useful to have a small number of utility functions which do not require the full BLAS code. The following functions fall into this category. 00244 00245 // This function copies the elements of the vector src into the vector dest. 00246 vector_int& operator=(const vector_int& other){copy(other);return (*this);} 00247 00248 // Function: int gsl_vector_int_swap (gsl_vector_int * v, gsl_vector_int * w) 00249 // This function exchanges the elements of the vectors v and w by copying. The two vectors must have the same length. 00250 00251 // ***** Exchanging elements 00252 00253 // The following function can be used to exchange, or permute, the elements of a vector. 00254 00255 // Function: int gsl_vector_int_swap_elements (gsl_vector_int * v, size_t i, size_t j) 00256 // This function exchanges the i-th and j-th elements of the vector v in-place. 00257 int swap_elements (size_t i, size_t j) {return gsl_vector_int_swap_elements (gsldata, i,j);} 00258 00259 // Function: int gsl_vector_int_reverse (gsl_vector_int * v) 00260 // This function reverses the order of the elements of the vector v. 00261 int reverse () {return gsl_vector_int_reverse (gsldata) ;} 00262 00263 // ******* Vector operations 00264 00265 // The following operations are only defined for real vectors. 00266 00267 // This function adds the elements of vector b to the elements of vector a, a'_i = a_i + b_i. The two vectors must have the same length. 00268 int operator+=(const vector_int &other) {return gsl_vector_int_add (gsldata, other.gsldata);} 00269 00270 // This function subtracts the elements of vector b from the elements of vector a, a'_i = a_i - b_i. The two vectors must have the same length. 00271 int operator-=(const vector_int &other) {return gsl_vector_int_sub (gsldata, other.gsldata);} 00272 00273 // Function: int gsl_vector_int_mul (gsl_vector_int * a, const gsl_vector_int * b) 00274 // This function multiplies the elements of vector a by the elements of vector b, a'_i = a_i * b_i. The two vectors must have the same length. 00275 int operator*=(const vector_int &other) {return gsl_vector_int_mul (gsldata, other.gsldata);} 00276 00277 // This function divides the elements of vector a by the elements of vector b, a'_i = a_i / b_i. The two vectors must have the same length. 00278 int operator/=(const vector_int &other) {return gsl_vector_int_div (gsldata, other.gsldata);} 00279 00280 // This function multiplies the elements of vector a by the constant factor x, a'_i = x a_i. 00281 int operator*=(int x) {return gsl_vector_int_scale (gsldata, x);} 00282 00283 // Function: int gsl_vector_int_add_constant (gsl_vector_int * a, const double x) 00284 // This function adds the constant value x to the elements of the vector a, a'_i = a_i + x. 00285 int operator+=(int x) {return gsl_vector_int_add_constant (gsldata,x);} 00286 00287 // This function multiplies the elements of vector a by the constant factor x, a'_i = x a_i. 00288 int operator/=(int x) {return gsl_vector_int_scale (gsldata, 1/x);} 00289 00290 // bool operators: 00291 bool operator==(const vector_int& other) const; 00292 bool operator!=(const vector_int& other) const { return (!((*this)==other));} 00293 00294 // stream output: 00295 // friend ostream& operator<< ( ostream& os, const vector_int& vect ); 00297 int sum() const; 00298 double norm2() const; 00299 00300 00301 // **** Finding maximum and minimum elements of vectors 00302 00303 // Function: double gsl_vector_int_max (const gsl_vector_int * v) 00304 // This function returns the maximum value in the vector v. 00305 00306 // Function: double gsl_vector_int_min (const gsl_vector_int * v) 00307 // This function returns the minimum value in the vector v. 00308 00309 // Function: void gsl_vector_int_minmax (const gsl_vector_int * v, double * min_out, double * max_out) 00310 // This function returns the minimum and maximum values in the vector v, storing them in min_out and max_out. 00311 00312 // Function: size_t gsl_vector_int_max_index (const gsl_vector_int * v) 00313 // This function returns the index of the maximum value in the vector v. When there are several equal maximum elements then the lowest index is returned. 00314 00315 // Function: size_t gsl_vector_int_min_index (const gsl_vector_int * v) 00316 // This function returns the index of the minimum value in the vector v. When there are several equal minimum elements then the lowest index is returned. 00317 00318 // Function: void gsl_vector_int_minmax_index (const gsl_vector_int * v, size_t * imin, size_t * imax) 00319 // This function returns the indices of the minimum and maximum values in the vector v, storing them in imin and imax. When there are several equal minimum 00320 // or maximum elements then the lowest indices are returned. 00321 00322 // Vector properties 00323 00324 // Function: int gsl_vector_int_isnull (const gsl_vector_int * v) 00325 // This function returns 1 if all the elements of the vector v are zero, and 0 otherwise. }; 00326 00327 }; 00328 00329 // When you add create a view it will stick to its with the view until you call change_view 00330 // ex: 00331 // matrix_float m(5,5); 00332 // vector_float v(5); 00333 // // ... 00334 // m.column(3) = v; //the 3rd column of the matrix m will equal v. 00335 class vector_int_view : public vector_int 00336 { 00337 public: 00338 vector_int_view(const vector_int& other):vector_int(){init(other);/*copy(other);*/} 00339 vector_int_view(const vector_int_view& other):vector_int(){init(other);/*copy(other);*/} 00340 00341 void init(const vector_int& other); 00342 void change_view(const vector_int& other){init(other);} 00343 private: 00344 }; 00345 00346 ostream& operator<< ( ostream& os, const vector_int & vect ); 00347 00348 00349 // vector_type<>::type is a template interface to vector_? 00350 // it is usefull for in templated situations for getting the correct vector type 00351 #define tmp_type_is_int 00352 #ifdef tmp_type_is 00353 typedef vector vector_double; 00354 template<class T> 00355 struct vector_type {typedef vector_double type;}; 00356 #else 00357 template<> struct vector_type<int> {typedef vector_int type;}; 00358 #endif 00359 #undef tmp_type_is_int 00360 00361 } 00362 #endif// _vector_int_h