It is easy to write a pure C function calculating the Frobenius norm:

double frobeniusNorm(double* matrix, int size1, int size2) { double result = 0.0; for(int i = 0; i < size1; ++i) { for(int j = 0; j < size2; ++j) { double value = *(matrix + (i*size2) + j); result += value * value; } } return sqrt(result); }This piece of code would look a bit counterintuitive to someone not used to how multidimensional arrays are managed in C, but once you get the point that they are flatted and seen just like a one dimensional array, the things start looking clearer.

We could have written the same function in a more friendly way, using the array notation, ending up in something like this:

double frobeniusNorm(double matrix[][MAT_SIZE_2], int size1) { double result = 0.0; for(int i = 0; i < size1; ++i) { for(int j = 0; j < MAT_SIZE_2; ++j) { double value = matrix[i][j]; result += value * value; } } return sqrt(result); }But the compiler has to know one matrix dimension, so that it could converts the matrix notation in the actual address of the memory cell containing its value, so we are forced to pass it (here is used an int constant, MAT_SIZE_2, whose definition is showed below).

In C, a 3x3 matrix would be typically defined in one of this two ways:

double dm[] = { // 1 1, 2, 3, 4, 5, 6, 7, 8, 9 }; const int MAT_SIZE_2 = 3; double dm2[][MAT_SIZE_2] = { // 2 {1, 2, 3}, {4, 5, 6}, {7, 8, 9} };1. In this way we keep the representation close to what is the actual way the data are stored in memory, but the code could be a bit confusing: is this a matrix or a vector?

2. Here we are stating clearly that we are dealing with a two dimensional vector, but this is what we have to do to call our generic Frobenius norm calculator:

frobeniusNorm(&dm2[0][0], 3, 3)We can't pass a bidimesional array when a bare pointer is expected, so we have to explicitly say to the compiler that it should get the address to the first element in the array.

We could save all this trouble using the Boost uBLAS matrix class:

// ... #include <boost/numeric/ublas/matrix.hpp> #include <boost/numeric/ublas/io.hpp> namespace ublas = boost::numeric::ublas; // ... ublas::matrix<double> matrix(3, 3); int k = 0; for(unsigned int i = 0; i < matrix.size1(); ++i) for(unsigned int j = 0; j < matrix.size2(); ++j) matrix(i, j) = ++k; // ... double frobeniusNorm(const ublas::matrix<double>& matrix) { double result = 0.0; for(unsigned int i = 0; i < matrix.size1(); ++i) { for(unsigned int j = 0; j < matrix.size2(); ++j) { double value = matrix(i, j); result += value * value; } } return sqrt(result); }Naturally, the above implementation of a Frobenius norm calculator is useful only to show how to work with an uBLAS matrix. In real code, it does not make much sense using uBLAS data structure and not its functions:

std::cout << ublas::norm_frobenius(matrix) << std::endl;

What if the value^2 i.e. Aij^2 resulted in in an overflow or an underflow i.e. inf or 0? How would you avoid this error?

ReplyDeleteHi Tyler. I don't think there is an easy answer for your problem. I would suggest you to have a look to the machine dependent way for enabling a check on floating point exception.

Delete