Feeds:
Posts
Comments

Archive for the ‘c++’ 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 »

the following code will result in a compilation error with gcc 4.7

#include <iostream>
int main() {
  std::cout << "filename::"__FILE__ << std::endl; // error
  return 0;
}

the error message being

error: unable to find string literal operator ‘operator"" __FILE__’

fortunately the solution is provided on this page. look under User-defined literals and whitespace.

all you got to do is to introduce space between the string literal “filename::” and the macro __FILE__ to read

std::cout << "filename::" __FILE__ << std::endl;

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 »

what a weird topic to talk about. i don’t even know if this piece
of code is any good. it will probably serve to remind me of the
usage of certain functions etc at a later date. may be it will
scare somebody off and may be some kind stranger will correct my
mistakes and point me to enlightened realms. all this does is to
compute the column sums of a matrix. notice that we maintain our
distance from the for and while loops. i have been admonished
myriads of time for not using the standard library functions so i
have avoided them here. what can i say i simply love the lambda
functions
over functors. although the compiler ultimately uses
functors behind the scenes but i love how clean and separated the
code looks.

#include <iostream>
#include <vector>
#include <algorithm>
#include <iterator>
#include <cstdlib>

int main(){
  // matrix
  std::vector<std::vector<int> > a = {{0,0,0}, {0,0,0}, {0,0,0}};
// c++0x allows initialization lists

  // randomly initialize using lambda expressions
  std::for_each(a.begin(), a.end(), [](std::vector<int>& i){
      std::generate(i.begin(), i.end(), [](){return rand()%10;});
    });

  // display using lamda expressions
  std::for_each(a.begin(), a.end(), [](std::vector<int>& i){
      std::copy(i.begin(), i.end(), std::ostream_iterator<int>(std::cout, "\t"));
      std::cout << std::endl;
    });

  // zero vector
  std::vector<int> zeros = {0,0,0};
  
  // column sums
  std::vector<int> csums = std::accumulate(a.begin(), a.end(), zeros, 
[](std::vector<int>& i, std::vector<int>& j)->std::vector<int>{
std::transform(i.begin(), i.end(), j.begin(), i.begin(), std::plus<int>());
return i;
});

  // display the column sums using ostream_iterators
  std::copy(csums.begin(), csums.end(), std::ostream_iterator<int>(std::cout, "\t"));

  return 0;
}

i now explain bits and parts of the mess

  // randomly initialize
  std::for_each(a.begin(), a.end(), [](std::vector<int>& i){
      std::generate(i.begin(), i.end(), [](){return rand()%10;});
    });

the above initializes the matrix a. since each element of a is itself a vector i use a lambda function that initializes the individual vectors. the lambda is as follows:

[](std::vector<int>& i){std::generate(i.begin(), i.end(), [](){return rand()%10;});}

the actual adding of columns is done using the std::accumulate() library function

  // column sums
  std::vector<int> csums = std::accumulate(a.begin(), a.end(), zeros, 
[](std::vector<int>& i, std::vector<int>& j)->std::vector<int>{
std::transform(i.begin(), i.end(), j.begin(), i.begin(), std::plus<int>());
return i;
});

that uses the lamda

[](std::vector<int>& i, std::vector<int>& j)->std::vector<int>{
std::transform(i.begin(), i.end(), j.begin(), i.begin(), std::plus<int>());
return i;
}

Read Full Post »

whenever possible one should use the standard library tools provided. std::mismatch() returns a pair of iterators to the first mismatch. examples of this are abundant on the WWW but i wanted to find all the mismatches and could not find an example online, hence this post. i don’t think the code is elegant so please let me know better versions of my code.

#include <iostream>
#include <vector>
#include <algorithm>
 
int main()
{
  std::vector<unsigned> a = {1,2,3,4,5,6};
  std::vector<unsigned> b = {1,2,33,4,55,66};
 
  // find the first mismatch
  auto p = std::mismatch(a.begin(), a.end(), b.begin());

  while(p.first != a.end())
    {
      std::cout << *p.first << std::endl;
      std::cout << *p.second << std::endl;
      p = std::mismatch(++p.first, a.end(), ++p.second);
    }
 
return 0;
}

havent’ tried the above on list container but i believe list iterators support the increment operator.

don’t forget to compile the above using -std=c++0x

Read Full Post »

guys when i run a code i get a lot of text dump on the console so i m in the habbit of issuing

“tmux clearhist && ./a.out”

so that when i scroll back i only see dump pertaining to only the current execution of a.out . turns out that this interacts weirdly with the standard input buffer. for a simple code such as the follows:

#include <iostream>
#include <string>
int main()
{
  std::string s;
  while(std::getline(std::cin,s))
    {
      std::cout << s << std::endl;
    }
return 0;
}

does not work at all with “tmux clearhist && ./a.out”.

by “does not work” i mean the program does not wait for any input and simply proceeds to completion. at first i thought that probably the standard input buffer has whitespaces or newlines sitting which i needed to clear before doing an input, but that did not work either.

so currently i have resorted to first clearing up history and then running my executable instead of chaining them together as one code. or one can use

tmux clearhist >/dev/null 2>&1 </dev/null && ./a.out

Read Full Post »

i have always defined matrices of two dimensions as a vector of a vector container as follows

vector< vector <double> > v; // note the distance between the two right angle brakets just before v

however i recently learned that c++0x allows one to get rid of the extra space between the right angle brackets as in the following

vector< vector <double> > v_old_style; 
vector< vector <double>> v_new_style; 

however emacs 23 breaks indentation when faced with this format

i also use 3 and 4 dimensional tensor like objects for which i would much like the following:

vector < vector < vector < vector <double> > > > v1; // don't like it
vector<vector<vector<vector<double>>>> v2; // like this a lot 

for me at least the saving in horizontal space matters when i am comparing code side by side … also its quite easy to figure out the dimensionality of the tensor at one glance.

i don’t know how to fix this 😦 help me please

Read Full Post »

Older Posts »