#include <iostream>

#include <boost/timer.hpp>

#include <boost/random.hpp>
#include <boost/random/mersenne_twister.hpp>

#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/vector_of_vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_sparse.hpp>

namespace ublas = boost::numeric::ublas;
using namespace std;

int main(int argc, char** argv)
{
	//const int numNodes = 196;
	//const int numTetrahedrons = 617;
	const int numNodes = 920;
	const int numTetrahedrons = 4196;

	//ublas::generalized_vector_of_vector<float, ublas::row_major, ublas::vector<ublas::mapped_vector<float> > > K(numNodes*3, numNodes*3);
	ublas::matrix<float> ke[numTetrahedrons];//array of element stiffness matrices (12x12)
	ublas::matrix<float> K(numNodes*3, numNodes*3);//global stiffness matrix
	ublas::vector<int> index(numTetrahedrons*12);//mapping from local element index to global index, 12 for each element
	boost::mt19937 mtrand;

	//fill index vector with random integers in [0, numNodes-1] and fill element stiffness matrix with random stuff
	for(int i=0; i<numTetrahedrons; ++i)
	{
		for(int j=0; j<4; ++j)
		{
			int r = mtrand();
			int ni = (int)((r*numNodes)/(double)(mtrand.max()));

			for(int k=0; k<3; ++k)
				index(i*12+j*3+k) = ni*3+k;
		}

		ke[i].resize(12, 12);

		for(int j=0; j<12; ++j)
			for(int k=0; k<12; ++k)
			{
				int r = mtrand();
				float v = (r*100000)/(double)(mtrand.max());
				ke[i](j, k) = v;
			}
	}

	boost::timer t;
	t.restart();

	//assemble K a few times
	const int N = 100;

	for(int a=0; a<N; ++a)
	for(int i=0; i<numTetrahedrons; ++i)
	{
		for(unsigned int j=0; j<12; ++j)
		{
			int jj = index(i*12+j);

			for(unsigned int k=0; k<12; ++k)
			{
				int kk = index(i*12+k);
				K(jj, kk) += ke[i](j, k);
				K(kk, jj) += ke[i](k, j);
			}
		}
	}

	cout << "time elapsed: " << t.elapsed()/N << endl;

	return 0;
}

Timeout