Accessor functions

In C++, a temporary cannot be passed by non-const reference. If you are slightly inexperienced and (for the sake of encapsulation) writing your accessors like this:

class Object
{
....
MPI::Intracomm GetRowWorld() { return rowWorld; }
MPI::Intracomm GetColWorld() { return colWorld; }
...
}

You'll be in deep trouble when you want to send these return values to a function. Such as this:

void foo(MPI::Intracomm & rowcomm, MPI::Intracomm & colcomm)
{
...
}

Object A;
foo( A.GetRowWorld(), A.GetColWorld());
// Won't compile (and you should be happy about it)


So, what's the fix? You can try to change the signature of foo:
void foo(const MPI::Intracomm & rowcomm, const MPI::Intracomm & colcomm);

You might think this is smarter as it compiles but now you are in even deeper trouble because your accessor functions return values and they go out of scope immediately. Therefore, when foo starts executing, the references are dangling. You might try to change the signature to this:
void foo(MPI::Intracomm rowcomm, MPI::Intracomm colcomm);

This works, but hey it might be terribly inefficient if MPI::Intracomm is a big object because you're creating temporaries everytime you pass by value. Just leave the signature as is and go fix your class definition instead:

class Object
{
....
MPI::Intracomm & GetRowWorld() { return rowWorld; }
MPI::Intracomm & GetColWorld() { return colWorld; }
MPI::Intracomm GetRowWorld() const { return rowWorld; }
MPI::Intracomm GetColWorld() const { return colWorld; }
...
}

static_cast and CRTP

I have been using CRTP (curiously recurring template pattern) for my Sparse Matrix classes for a number of reasons:

1) Performance, as this is a high-performance computing issue.
2) Parameter type checking without the need to inspect RTTI (run-time type information).

Here is the deal:

template < class DER >
class SpMat
{ ... }

template < class DER>
template <typename SR>
void SpMat<DER>::SpGEMM(SpMat<DER> & A, SpMat<DER> & B, ...)
{
// check for conformance, etc.
...
static_cast< DER* >(this)->template PlusEq_AtXBt<SR>(static_cast< DER >(A), static_cast< DER >(B));
...
}

class SpDCCols: public SpMat <SpDCCols >
{
public:
template <typename SR>
int PlusEq_AtXBt(const SpDCCols & A, const SpDCCols & B);
....
}

This won't compile because static_cast< DER >(A) fails as no explicit conversion operator from SpMat< SpDCCols > to SpDCCols is supplied, although the latter class is derived from the former. This also has to do with the semantics of static_cast() operator:
Stroustrup says "For a built-in type T, T(e) is equivalent to static_cast(e)"

So, are we gonna provide a conversion constructor? Of course not. One can convert a pointer/reference of a type A to a pointer/reference of a type B if A is a base class of B.
Therefore, the immediate solution is to cast to a reference rather than a concrete object, which will also avoid unnecessary explicit conversion and make your code faster:
static_cast< DER & >(A)

Template argument deduction

C++ has an excellent mechanism to deduce arguments. However, it is not free of pitfalls. Suppose you are writing a library function to perform a matrix-vector multiplication on a particular semiring (which is defined as a class with static ::multiply and ::add functions):

// Perform y <- y + Ax template < typename T, typename SR>
void MultiplyAdd(Matrix<T> & A, Vector<T> & x, Vector<T> & y){
.... // usual conformance checks
for(int i=0; i< A.m; ++i){
for(int j=0; j< A.n; ++j){
temp = SR::multiply(A[i][j], x[j]);
y[i] = SR::add(temp, y[i]);
}}}

Here,the template argument T can be deduced from the function arguments but SR can not. If you release your library like this, people will need to explicitly call it like:
MultiplyAdd <int, minplus>

which is kind of annoying. Normally, you'd like to list deducible arguments first so that the user can get away with declaring only the non-deducible ones:
MultiplyAdd <minplus>

----------------------

My second point will be on using function pointers in STL library calls. If you are using a templated function pointer, don't expect the compiler to deduce the arguments for you. Example:

template <typename T>
bool compare_second(pair<T,T> t1, pair<T,T> t2)
{
return t1.second < t2.second;
}

const int N = 100;
pair pairarray[N];
.... // somehow populate the array
sort(pairarray, pairarray +N, compare_second); // cryptic compiler error!

The correct way to do it is obviously:
sort(pairarray, pairarray +N, compare_second<int>);

How to get rid of dot/slash when running executables on Linux

... and why you shouldn't be doing this unless it is absolutely necessary.

But anyway, sometimes there is this old/legacy program that tries to execute another program inside the same directory without the dot/slash, and we don't have access to its source code or we don't wanna access its source code (for reasons apparent to a CS person)

Here: http://www.linfo.org/dot_slash.html

Have fun.

Heap Memory in Multithreaded NUMA Programs

Ok, here is the deal: You want to parallelize you application (say matrix multiplication) and you want asynchronous computation. If you have a shared memory NUMA machine, threads seem to be the obvious choice. The caveat is memory management. Calls to "new" and "delete" are serialized because heap memory is always shared among threads. Also, and probably more costly, there are memory contentions, false sharing, etc.
I tried to use hoard library here is what happened:

Without hoard, using plain GNU C++ library (which might be using ptmalloc?)

Loading Matrices
Loaded !
Loading took 14.019800 seconds
Multiplications started
Transposition took 0.423375 seconds
Multiplication took 4.816739 seconds
Retransposition took 0.486832 seconds
Multiplications finished
8.713002 seconds elapsed (including thread creation cost)


Using hoard:

Loading Matrices
Loaded !
Loading took 2.489643 seconds
Multiplications started
Transposition took 0.463726 seconds
Multiplication took 5.337530 seconds
Retransposition took 0.493123 seconds
Multiplications finished
9.304911 seconds elapsed (including thread creation cost)


What's wrong here? Hoard makes loading matrices really fast (2.4 instead of 14 seconds), but it slowed down the multiplication at the same time (5.3 sec instead of 4.8).
Note that the code uses 16 threads on a 16 core machine.

Boost::bind

Boost::bind library is a useful library that I usually use for my multithreaded code, even though it is not really required in any sense. Instead of writing my own function objects to be called by boost::thread, I delegate the job to boost::bind.
However, there is a caveat and I want to share it with you today.

"The arguments that bind takes are copied and held internally by the returned function object".
So,that means the copy constructors are called for your arguments. If you intentionally used reference parameters to avoid data copying during function calls, using boost::bind blindly just nullifies your efforts :(

Two solutions exist as far as I know:

1) Plain and simple, pass pointers instead. Ok, this is too C-like.
2) Force bind to hold a reference instead of calling the copy constructors and keeping copies of function arguments.boost::ref and boost::cref does that for you

Analytical Geometry vs. Programming Languages

I guess anyone with a little bit of programming experience knows the storage of multidimensional arrays. For simplicity, I will just deal with 2D arrays.
In C,C++,etc. memory storage is row-major ordered, meaning that rows are stored consecutively. A 3X4 grid would be stored in the memory as follows:

1 2 3 4
5 6 7 8
9 10 11 12

In Fortran or Matlab, that would be:

1 5 9
2 6 10
3 7 11
4 8 12

So now, my old knowledge from analytical geometry forces me to think about the following axis system:

y
|
|
|
|
|
---------------- x

That certainly doesn't fit to the row-major storage system. I have seen so many people with math background that would access the ith row and the jth column like A[j][i] in C++. Why? Because that is the geometric way to do it. By moving along the x-axis, you change the current column and the x-axis naturally seems to be the first axis to be written in row-major storage.

A[x-axis value][y-axis value]

But that's plain wrong, in fact all accesses are of the form A[y-axis value][x-axis value]. If you want to increase the current column you're in, you change the second dimension like A[i][j++]

Weird, right?