Feeds:
Posts
Comments

Archive for the ‘gsl’ Category

gsl provides a random number generator that i use in my work. one has to allocate the random number generator at the start of the program and then deallocate the memory at the end of the program. c++11 provides smart pointers which keep track of the objects they point to for the purpose of memory management. when the smart pointer goes out of scope the resource is automatically deallocated when the associated object is destroyed. that’s one headache less for the programmer!

the following code demonstrates how one uses the gsl random number interface in a c-style code and compares it with the use of smart pointers to automatically manage resource deallocation for you.

#include <iostream>
#include <gsl/gsl_rng.h>
#include <memory> // header file for unique_ptr shared_ptr

// c-style usage
void c() {
  gsl_rng *r = gsl_rng_alloc(gsl_rng_taus); // allocation
  std::cout << gsl_rng_uniform(r) << std::endl;
  gsl_rng_free(r); // de-allocation
}

// c++ style usage with shared_ptr
void cpp_shared() {
  std::shared_ptr<gsl_rng> r(gsl_rng_alloc(gsl_rng_taus), gsl_rng_free);
  std::cout << gsl_rng_uniform(r.get()) << std::endl;
}

// c++ style usage with unique_ptr
void cpp_unique() {
  std::unique_ptr<gsl_rng, void(*)(gsl_rng*)> r(gsl_rng_alloc(gsl_rng_taus), gsl_rng_free);
  std::cout << gsl_rng_uniform(r.get()) << std::endl;
}

int main() {
  c();
  cpp_shared();
  cpp_unique();
  return 0;
}

the functions gsl_rng_alloc() and gsl_rng_free() are the memory allocaters and deallocaters respectively.

std::shared_ptr<gsl_rng> r(gsl_rng_alloc(gsl_rng_taus), gsl_rng_free);

in the code fragment above gsl_rng is the type of object r points to. notice how the smart pointer constructor is provided with two arguments 1) the custom allocator gsl_rng_alloc() and 2) the custom deleter gsl_rng_free(). however gsl_rng_uniform() expects a const gsl_rng* argument for which we use the get() member function of the shared_ptr class i.e. r.get() returns a pointer to the gsl_rng object it is pointing to.

std::unique_ptr<gsl_rng, void(*)(gsl_rng*)> r(gsl_rng_alloc(gsl_rng_taus), gsl_rng_free);

this one is a bit ugly looking. notice that unique_ptr requires the specification of 2 template parameters, the type of object being pointed to and the type of the custom deleter. the function gsl_rng_free() is declared as void gsl_rng_free(gsl_rng*) thus the type of this deleter is specified as void(*)(gsl_rng*) i.e. a function which accepts a gsl_rng* and whose return type is void. here again one must use the get() member function to get a reference to the pointed to object.

leave comments, rants and improvements!

Read Full Post »

to call a LAPACK routine from your C++ code i used to do the following

extern "C"
{
  //http://www.netlib.org/lapack/double/dgesv.f
  void dgesv_(int* n, int *nrhs, double* a, int* lda, int* ipiv, double* b, int* ldb, int* info);
}

int main()
{
  int n = 2;
  int nrhs = 2;
  double a[] = {1,3,2,4}; //colum-major mode 
  int lda = n;
  int* ipiv = new int[n];
  double b[] = {19, 43, 22, 50}; //column-major mode
  int ldb = n;
  int info;
  
  dgesv_(&n, &nrhs, a, &lda, ipiv, b, &ldb, &info);

  delete [] ipiv;
  return 0;
}

the above code solves the system of linear equations \begin{bmatrix} 1 & 2 \\ 3 & 4 \end{bmatrix}\mathbf{X} = \begin{bmatrix} 19 & 22 \\ 43 & 50 \end{bmatrix}. notice that in the code i had to store the matrices in column-major mode that is used by FORTRAN. i generally write helper routines to convert from row-major mode to column major mode before the subroutine call and then convert the results back to row-major mode after returning from the call. also one has to wrap the call to the subroutine in an extern “C” wrapper. both of these were annoyances that more or less kept me away from LAPACK, i could get most of my work done using GSL. while GSL provides a BLAS interface it does not provide any LAPACK interface and i have always wished to be able to use LAPACK.

it turns out that the guys at the intel math kernel library(imkl) and lapack guys collaborated to bring a C language interface to lapack . on archlinux i obtained the lapacke package from the AUR. below is the code using lapacke

#include <lapacke.h> // <--- header file
int main()
{
  double a[] = {1,2,3,4}; //NO need for column-major mode
  double b[] = {19, 22, 43, 50}; //NO need for column-major mode
  int n = 2;
  int nrhs = 2;
  int lda = n;
  int ipiv[n];
  int ldb = n;

  int info = LAPACKE_dgesv(LAPACK_ROW_MAJOR, n, nrhs, a, lda, ipiv, b, ldb);
 
  return 0;
}

compile the above code using g++ main.cpp -llapacke. notice that by using LAPACK_ROW_MAJOR on LAPACK_COL_MAJOR one can specify the storage mode for matrices and therefore there is no need to preprocess the matrices before calling the routine.

i am aware of the CLAPACK project as well which is an f2c’d version of the fortran version but i always wanted a pure C version. this i thought would fill my needs. alas that is NOT the case. i will not be using lapacke because it is not c++0x compatible 😦 and i use a lot of c++0x. i wish GSL would provide a lapack interface as well! my search for a scientific computing library for C++ continues.

Read Full Post »

it is possible to specify your own random number generator to be used with std::random_shuffle() which accepts either of the following:
0) function object a.k.a functor
1) pointer to a function

this article will use the pointer to a function approach.

first we create a unary function that accepts an appropriate type of argument. by appropriate i mean a type which is convertible from std::iterator_traits::difference_type. we will be using the gsl_rng_uniform_int() function provided by the GSL library.

long unsigned myrng(long unsigned n)
{
  return gsl_rng_uniform_int(rng, n);
}

notice how using myrng() we will be able to call a binary function using a unary interface. of course one must ensure that the variable rng is globally defined.

now that we have defined the unary function we need to define a pointer to this function, which we will pass to the
std::random_shuffle() function.

this is how we achieve that

long unsigned (*fp)(long unsigned) = myrng;

all the parts are in place for the final working program which is listed below:

#include <vector>

#include <gsl/gsl_rng.h>
// for gsl_rng_uniform_int

#include <algorithm>
// for random_shuffle

// global definition
gsl_rng * rng;

// unary function 
long unsigned myrng(long unsigned n)
{
  return gsl_rng_uniform_int(rng, n);
}

// pointer to unary function
long unsigned (*fp)(long unsigned) = myrng;

int main()
{
  rng = gsl_rng_alloc(gsl_rng_default);

  // [1.0, 2.0, 3.0, 4.0]
  std::vector <double> a;
  a.push_back(1.0);
  a.push_back(2.0);
  a.push_back(3.0);
  a.push_back(4.0);

  std::random_shuffle(a.begin(), a.end(), fp);

  gsl_rng_free(rng);

  return 0;
  }

i m not too happy with this implementation as this uses a global definition of gsl_rng *. i will also try to write a version using functors at a later date.

Read Full Post »

gsl has helped me liberate myself from the clutches of MATLAB. if you are in academics then i strongly urge you to use open source libraries whenever where ever possible. there are both practical as well as philosophical reasons for this. any ways, end of rant and now onto the topic of this post.

both gsl_rng_uniform() and gsl_rng_uniform_pos() generate random numbers between 0 and 1 while the former includes 0 the latter excludes it and both of them exclude 1. now why is it that it does not generate 1.0 ? probably numbers as high as 0.9999 might be good approximation to 1.0 in which case there is no practical difficulty.

Sampling from a random number generator

Read Full Post »

gsl vector stride

gsl is an excellent excellent library for c/c++ users and if you are in an academic institute and are not using gsl then please hurry up and start using it.

i m on the gsl mailing list where i gained the impression that many new users have problem understanding what stride is. please refer to this page to read about it more.

consider the following array

double array[] = { 1, 2, 3, 4, 5, 6, 7, 8} ;

if you start at 1 and then traverse the array with a stride of 1 you will first go to 2 then to 3 then to 4 then to 5 then to 6 then to 7 and finally to 8.

if you start at 1 and then traverse the array with a stride of 2 you will first go to 3 then to 5 then to 7.

if you start at 1 and then traverse the array with a stride of 3 you will first go to 4 then to 7.

hope that clears it for you !

Read Full Post »