I am performing these tests using Boost 1.49 for MSVC 2010. Here is what I did to use uBLAS vector and matrix:
#include <boost/numeric/ublas/vector.hpp> // 1 #include <boost/numeric/ublas/matrix.hpp> // 2 #include <boost/numeric/ublas/io.hpp> // 3 namespace ublas = boost::numeric::ublas; // 41. For the uBLAS vector definition.
2. uBLAS matrix.
3. I/O uBLAS support, it makes output to console easier.
4. Synonym definition to the boost numeric ublas namespace.
Error handling
When running in debug mode, uBLAS performs a strict data check. Usually we need to remove all this checks from the production code for efficiency reasons, making the uBLAS code faster, but moving in this way to our code the responsibility of checking all the weirdness that could happen.
If the symbol NDEBUG is defined as preprocessor constant, the uBLAS code generated by the compiler is lean and mean. Otherwise is safer, but bigger and slower.
So, if we plan to keep in production code the safer uBLAS version, we should try/catch for std::exception our uBLAS code, to react adequately in case of anomalies. Otherwise, for NDEBUG generated code, no check is performed, no exception is thrown by uBLAS. Simply if we give garbage in, we get garbage out. It is our responsibility to ensure that this won't happen.
A couple of uBLAS vectors
To begin playing a bit around with uBLAS products, we need a few of vectors:
ublas::vector<double> v1(6); // 1
ublas::vector<double> v2(6);
for(int i = 0; i < 6; ++i)
{
v1[i] = i;
v2[i] = i + 1;
}
ublas::vector<double> vs(3);
for(int i = 0; i < 3; ++i) { vs[i] = i + 1; }
std::cout << v1 << v2 << vs << std::endl; // 2
1. The uBLAS vector behave slightly differently from STL vector. We can specify in the constructor its size, and the number of element passed is allocated without being initialized.2. Dump to standard output console the generated uBLAS vectors:
[6](0,1,2,3,4,5)[6](1,2,3,4,5,6)[3](0,0,0)Inner product
We can calculate the inner product between two vectors through the uBLAS function inner_prod(). It takes in input two vectors with the same size, and returns a scalar value:
std::cout << ublas::inner_prod(v1, v2) << std::endl; // 1
std::cout << std::inner_product(v1.begin(), v1.end(), v2.begin(), 0.0) << std::endl; // 2
try // 3
{
double badValue = ublas::inner_prod(v1, vs); // 4
std::cout << "Bad inner product: " << badValue << std::endl; // 5
}
catch(const std::exception& ex)
{
std::cout << "Can't get inner product: " << ex.what() << std::endl; // 6
}
1. uBLAS inner product, in this case the expected result is 70.2. We could use STL instead.
3. Remember that makes sense try/catching uBLAS functions only if NDEBUG has not been defined
4. The two vectors involved in an inner product should have the same size. Here this is not true, trouble should be expected.
5. This line is going to be executed only if NDEBUG has been defined, and in that case the output will be a meaningless value.
6. If we are not running in NDEBUG mode, the uBLAS check on the vectors size would fail, and an exception is thrown, and here it would be caught.
Outer product
The outer products of two vectors is a matrix where the x,y element is the product of the x element of the first vector by the y element of the second one, as shown by this very basic implementation:
namespace my
{
ublas::matrix<double> outer_prod(const ublas::vector<double>& v1, const ublas::vector<double>& v2)
{
ublas::matrix<double> m(v1.size(), v2.size());
for(unsigned int i = 0; i < v1.size(); ++i)
for(unsigned int j = 0; j < v2.size(); ++j)
m(i, j) = v1[i] * v2[j];
return m;
}
}
There is no constrain on the vector sizes. If they are the same, the result is a square matrix, otherwise we'll get a rectangular one:std::cout << "uBLAS outer product /1: " << ublas::outer_prod(v1, v2) << std::endl; std::cout << "uBLAS outer product /2: " << ublas::outer_prod(v2, vs) << std::endl; ublas::vector<double> v0(0); std::cout << "uBLAS outer product /3: " << ublas::outer_prod(v1, v0) << std::endl; // 11. It is not illegal multiply for a zero-sized vector. In this specific case we get a 6x0 matrix.
Left-right matrix-vector product
The product between matrix and vector is not commutative. The vector could be on the right only if its size matches the matrix size1(), and on the left if it matches size2():
ublas::matrix<double> m1(6,3); // 1
for(unsigned int i = 0; i < m1.size1(); ++i)
for(unsigned int j = 0; j < m1.size2(); ++j)
m1(i, j) = (i + j);
std::cout << ublas::prod(v1, m1) << std::endl; // 2
std::cout << ublas::prod(m1, vs) << std::endl; // 3
try
{
std::cout << "Right " << ublas::prod(vs, m1) << std::endl; // 4
}
catch(const std::exception& ex)
{
std::cout << "Can't get outer product: " << ex.what() << std::endl;
}
1. A 6x3 rectangular matrix.2. Vector v1 size is 6, matching with matrix size1.
3. All right, vector vs size is 3.
4. Size mismatch! If NDEBUG is not defined, uBLAS here throws an exception. Otherwise, we get a meaningless result.
Matrix-matrix product
We can think to the vector-matrix product as a special case of a matrix to matrix product. So, here too we should ensure that first matrix size2 matches to second matrix size1. Otherwise we are bound to get an exception (no NDEBUG) or nonsensical output (if NDEBUG is defined).
ublas::matrix<double> m2(3,6); // 1
for(unsigned int i = 0; i < m2.size1(); ++i)
for(unsigned int j = 0; j < m2.size2(); ++j)
m2(i, j) = (10 - i - j);
std::cout << "Output matrix: " << ublas::prod(m1, m2) << std::endl; // 2
ublas::matrix<double> mq(6,6); // 3
for(unsigned int i = 0; i < mq.size1(); ++i)
for(unsigned int j = 0; j < mq.size2(); ++j)
mq(i, j) = i + j;
try
{
std::cout << "Bad: " << ublas::prod(m1, mq) << std::endl; // 4
}
catch(const std::exception& ex)
{
std::cout << "Can't get matrix-matrix product: " << ex.what() << std::endl;
}
1. A 3x6 rectangular matrix.2. The first matrix is a 6x3 one, so they could be multiplied.
3. A 6x6 square matrix.
4. Here the no NDEBUG version throws an exception, if NDEBUG was defined, we should expect silly values in the resulting matrix.
Element product
The element product takes as operand two vectors, or matrices, having the same sizes, and produces as result a third vector (or matrix) with the same dimensions, and where each element contains a value that is the product of the original elements in the same position.
std::cout << ublas::element_prod(v1, v2) << std::endl; // 1
try
{
std::cout << ublas::element_prod(v1, vs) << std::endl; // 2
}
catch(const std::exception& ex)
{
std::cout << "Bad element product: " << ex.what() << std::endl;
}
1. The vectors v1 and v2 have the same size, we can perform an element_prod() on them.2. Different sizes in the passed vectors result in an exception (if NDEBUG not defined) or in an inconsistent result (NDEBUG defined).
Same for matrices:
ublas::matrix<double> m3(6,3);
for(unsigned int i = 0; i < m3.size1(); ++i)
for(unsigned int j = 0; j < m3.size2(); ++j)
m3(i, j) = (7 - i - j);
std::cout << ublas::element_prod(m1, m3) << std::endl; // 1
try
{
std::cout << ublas::element_prod(m1, m2) << std::endl; // 2
}
catch(const std::exception& ex)
{
std::cout << "Bad element product: " << ex.what() << std::endl;
}
1. This should works fine.2. A mismatch in the matrices' sizes leads to an exception (NDEBUG not defined) or garbage in the output matrix (NDEBUG defined).
Element division
Element division is about the same of element product:
std::cout << ublas::element_div(v1, v2) << std::endl; // 1
std::cout << ublas::element_div(v2, v1) << std::endl; // 2
try
{
std::cout << ublas::element_div(v1, vs) << std::endl; // 3
}
catch(const std::exception& ex)
{
std::cout << "Bad element product: " << ex.what() << std::endl;
}
std::cout << ublas::element_div(m1, m3) << std::endl; // 4
try
{
std::cout << ublas::element_div(m1, m2) << std::endl; // 5
}
catch(const std::exception& ex)
{
std::cout << "Bad element product: " << ex.what() << std::endl;
}
1. This works fine.2. This works fine, too. But given that the first element in the v1 vector is zero, the first element in the resulting vector is not a valid double, represented as 1.#INF in Windows.
3. Vectors with different sizes, their element division leads to an exception (NDEBUG not defined) or to nonsensical result (NDEBUG defined).
4. Matrices with the same sizes, it works alright with the same caveat in (2), given that we have a zero divisor.
5. Like (3), mismatching sizes for matrices lead to troubles.
No comments:
Post a Comment