C and C++

Here are some tips for those who wish to use the compiled programming language C or C++ for scientific work instead of some interpreted programming language (Matlab, ...). I assume, that you have some understanding of the fundamentals of the syntax of C and C++ but you are looking for high level functionality such as plotting, random number generation and handling of matrices.

I recommend two styles:

  1. If you are using plain C, then I suggest to use the GSL (Gnu Scientific Library).
  2. If you are using C++, then I suggest to use the boost library.

I explain my ideas by giving examples.

Visualization

Since most scientific writing is done in LaTeX I suggest to visualize data as follows. In the program, the data is written to a text file. Then we use the LaTeX package pgfplots to plot this data. (See my examples in the Graphics in LaTeX with PGF section.)

Consider the following three files:

  1. The C style file viz.c:

    #include <stdio.h>
#include <math.h>

int main(int argc, char **argv) {
  unsigned int N = 30, n;
  float Ts = 0.05, f = 1.7, t, x, y;
  FILE *fp;
  char fname[] = "viz.dat";
  fp = fopen(fname,"w");
  if (fp == NULL) {
    fprintf(stderr, "Error: Failed to open file "
	    "\"%s\" in write mode.\n", fname);
    return -1;
  }
  fprintf(fp, "      t\t     sin\t     exp\n");
  for(n = 0; n < N; n = n++) {
    t = n*Ts;
    x = cos(t*f*2*M_PI);
    y = exp(-t);
    fprintf(fp, "%1.5f\t%+1.5f\t%+1.5f\n", t, x, y);
  }
  fclose(fp);
  printf("done\n");
  return 0;
}

  2. The C++ style file viz.cpp:

    #include <iostream>
#include <fstream>
#include <cmath>
 
int main(int argc, char **argv) {
  unsigned int N = 30, n;
  float Ts = 0.05, f = 1.7, t, x, y;
  char fname[] = "viz.dat";
  std::ofstream file(fname);
  if (! file.is_open()) {
    std::cerr << "Error: Failed to open file \"" 
	      << fname << "\" in write mode." << std::endl;
    return -1;
  }
  file << std::scientific;
  file << "t\t" << "sin\t" << "exp" << std::endl;
  for(n = 0; n < N; n = n++) {
    t = n*Ts;
    x = cos(t*f*2*M_PI);
    y = exp(-t);
    file << t << "\t" << x << "\t" << y << std::endl;
  }
  file.close();
  std::cout << "done" << std::endl;
  return 0;
}

  3. Finally the LaTeX file viz.tex:

    \documentclass{article}
\usepackage{pgfplots}

\begin{document}

\begin{tikzpicture}
  \begin{axis}
    \addplot table[x=t,y=sin]{viz.dat};
    \addplot table[x=t,y=exp]{viz.dat};
    \legend{Sinusoid,Exponential}
  \end{axis}
\end{tikzpicture}

\end{document}

We can compile the C file or the C++ file with the GNU compiler. The resulting program produces our data file viz.dat which we want to plot. Finally, we use LaTeX to produce a pdf file containing the plot.

The sequence of commands for the C style file is:

gcc viz.c -lm
./a.out
pdflatex viz.tex

The sequence of commands for the C++ style file is:

g++ viz.cpp
./a.out
pdflatex viz.tex

The plot should look as follows:

\begin{axis} \addplot table[x=t,y=sin] {$wd/comp/prog/c-and-c++/viz.dat}; \addplot table[x=t,y=exp] {$wd/comp/prog/c-and-c++/viz.dat}; \legend{Sinusoid,Exponential} \end{axis}

Random Numbers

In plain C, GSL provides good random number generators. The following example uses the default generator: (Compile with gcc random.c -lgsl -lgslcblas.)

#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

const gsl_rng_type *T;
gsl_rng *r;

int main (void) {
  /* Here, we initialize the random number generator. */
  gsl_rng_env_setup();
  T = gsl_rng_default;
  r = gsl_rng_alloc(T);

  unsigned int N = 50, n;
  float Ts = 0.01, f = 5.0, t, x, y;
  FILE *fp;
  char fname[] = "random.dat";
  fp = fopen(fname, "w");
  if (fp == NULL) {
    fprintf(stderr, "Error: Failed to open file "
	    "\"%s\" in write mode.\n", fname);
    return -1;
  }
  fprintf(fp, "      t\t     sin\t   noisy\n");
  for(n = 0; n < N; n = n++) {
    t = n*Ts;
    x = cos(t*f*2*M_PI);
    y = x + gsl_ran_gaussian(r, 0.1);
    fprintf(fp, "%1.5f\t%+1.5f\t%+1.5f\n", t, x, y);
  }
  fclose(fp);
  gsl_rng_free(r);
  printf("done\n");
  return 0;
}

In C++, boost provides good random number generators. The following example uses a Mersenne Twister generator. (Compile with g++ random.cpp.)

#include <iostream>
#include <fstream>
#include <boost/random.hpp>

namespace bst = boost;

// Here, we construct and initialize the random number generator.
static bst::mt19937 rng;
static bst::normal_distribution<double> norm_dist(0, 0.1);
static bst::variate_generator<bst::mt19937&,
  bst::normal_distribution<double> > randn(rng, norm_dist);

 
int main(int argc, char **argv) {
  unsigned int N = 50;
  float Ts = 0.01, f = 5.0, t, x, y;
  char fname[] = "random.dat";
  std::ofstream file(fname);
  if (! file.is_open()) {
    std::cerr << "Error: Failed to open file \"" << fname
	      << "\" in write mode." << std::endl;
    return -1;
  }
  file << std::scientific;
  file << "t\t" << "sin\t" << "noisy" << std::endl;
  for(int n = 0; n < N; n = n++) {
    t = n*Ts;
    x = cos(t*f*2*M_PI);
    y = x + randn();
    file << t << "\t" << x << "\t" << y << std::endl;
  }
  file.close();
  std::cout << "done" << std::endl;
  return 0;
}

Note

Starting from boost version 1.47, the random number part of the boost library has been moved to the namespace boost/random and header files are split up in several files. The first lines of this example have to be changed to

#include <iostream>
#include <fstream>
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/variate_generator.hpp>
#include <boost/random/normal_distribution.hpp>

using namespace std;
using namespace boost::random;

The resulting data should look as follows:

\begin{axis} \addplot table[x=t,y=sin]{$wd/comp/prog/c-and-c++/random.dat}; \addplot table[x=t,y=noisy]{$wd/comp/prog/c-and-c++/random.dat}; \legend{Sinusoid, Noisy} \end{axis}

Matrix Inversion

In plain C, GSL provides structures for matrices and vectors. A general non-singular matrix can be inverted using the LU-factorization as shown in the following example matrix-inverse.c. (Compile with gcc matrix-inverse.c -lgsl -lgslcblas.)

Example output

With high probability the matrix A is non-singular and the output is, e.g.:

A =
0.0133919
-0.0088101
0.167441
0.0733641
0.0997525
-0.12775
-0.239672
-0.067928
-0.00390913

A*Ai=
1
0
-1.38778e-17
0
1
5.55112e-17
-2.08167e-17
1.03216e-16
1

Ai*A=
1
-5.55112e-17
-7.97973e-17
0
1
-3.46945e-18
0
6.93889e-18
1

done

#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>

const gsl_rng_type *T;
gsl_rng *r;

/* matrix inversion routing using LU-factorization */
int invert_matrix(const gsl_matrix *A, gsl_matrix *Ai)
{
  /* check matrix sizes */
  if(! A->size1 == A->size2
     && Ai->size1 == Ai->size2
     && A->size1 == Ai->size1)
    return (-1);
  /* create a permutation matrix for the LU-factorization */
  gsl_permutation * P = gsl_permutation_alloc(A->size1);
  int s; 			/* a dummy variable */
  /* allocate a buffer matrix */
  gsl_matrix *AA = gsl_matrix_alloc(A->size1, A->size2);
  gsl_matrix_memcpy(AA, A);	/* copy A to AA */
  /* perform LU-factorization */
  gsl_linalg_LU_decomp(AA, P, &s);
  /* backsubstitute to get the inverse */
  gsl_linalg_LU_invert(AA, P, Ai);
  return 0;
}

int main (void) 
{
  /* initialize the random number generator */
  gsl_rng_env_setup();
  T = gsl_rng_default;
  r = gsl_rng_alloc(T);
  
  /* generate a random matrix */
  gsl_matrix * A = gsl_matrix_alloc(3,3),
    * B = gsl_matrix_alloc(3,3),
    * Ai = gsl_matrix_alloc(3,3);
  int n, m;
  for (n = 0; n < 3; n++)
    for (m = 0; m < 3; m++)
      gsl_matrix_set(A,n,m,gsl_ran_gaussian(r, 0.1));
  printf("A =\n");
  gsl_matrix_fprintf(stdout,A,"%g");
  invert_matrix(A,Ai);
  /* compute B = A*Ai */
  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, A, Ai, 0.0, B);
  printf("\nA*Ai=\n");
  gsl_matrix_fprintf(stdout,B,"%g");
  /* compute B = Ai*A */
  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, Ai, A, 0.0, B); 
  printf("\nAi*A=\n");
  gsl_matrix_fprintf(stdout,B,"%g");

  /* free all the allocated objects */
  gsl_matrix_free(A);
  gsl_matrix_free(B);
  gsl_matrix_free(Ai);
  gsl_rng_free(r);
  printf("\ndone\n");
  return 0;
}

In C++, we use again boost and boost::ublas to perform the same matrix inversion based on the LU-factorization as shown in the following example matrix-inverse.cpp. (Compile with g++ matrix-inverse.cpp.)

#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/random.hpp>

namespace bst = boost;
namespace ubl = boost::numeric::ublas;

static bst::mt19937 rng;
static bst::normal_distribution<double> norm_dist(0,1);
static bst::variate_generator<bst::mt19937&,
  bst::normal_distribution<double> > randn(rng, norm_dist);

// matrix inversion routine using lu_factorize and lu_substitute
template<class T>
bool invert_matrix(ubl::matrix<T> A, ubl::matrix<T>& Ai)
{
  // create a permutation matrix for the LU-factorization
  ubl::permutation_matrix<std::size_t> P(A.size1());
  // perform LU-factorization, this will change A
  if(ubl::lu_factorize(A, P))
    return false;
  // create identity matrix
  Ai.assign(ubl::identity_matrix<T> (A.size1()));
  // backsubstitute to get the inverse
  ubl::lu_substitute(A, P, Ai);
  return true;
}

int main(int argc, char **argv)
{
  ubl::matrix<double> A(3, 3), Ai(3, 3);
  // generate a random matrix
  for(int i=0; i<3; i++)
    for(int j=0; j<3; j++)
      A(i,j) = randn();
  std::cout << "A    = " << A << std::endl;
  if(invert_matrix(A, Ai))
    {
      std::cout << "Ai   = " << Ai << std::endl;
      std::cout << "Ai*A = " << prod(Ai,A) << std::endl;
      std::cout << "A*Ai = " << prod(A,Ai) << std::endl;
    }
  else
    std::cout << "A is singular" << std::endl;
  return 0;
}

An example output is, e.g.:

A    = [3,3]((0.357628,0.14061,1.39772),(0.0460402,-1.46791,-0.770404),(0.274781,-0.892802,1.63291))
Ai   = [3,3]((4.84101,2.31866,-3.04983),(0.450192,-0.313718,-0.533363),(-0.568484,-0.561703,0.833999))
Ai*A = [3,3]((1,8.88178e-16,-1.77636e-15),(0,1,-1.11022e-16),(-2.77556e-17,-1.11022e-16,1))
A*Ai = [3,3]((1,0,0),(0,1,-1.11022e-16),(-3.33067e-16,0,1))