@@ -76,9 +76,42 @@ where the first argument \c v1 is a vector and the second argument \c 2*v2 is an
These examples are just intended to give the reader a first impression of how functions can be written which take a plain and constant Matrix or Array argument. They are also intended to give the reader an idea about the most common base classes being the optimal candidates for functions. In the next section we will look in more detail at an example and the different ways it can be implemented, while discussing each implementation's problems and advantages. For the discussion below, Matrix and Array as well as MatrixBase and ArrayBase can be exchanged and all arguments still hold.
+
+\section TopicUsingRefClass How to write generic, but non-templated function?
+
+In all the previous examples, the functions had to be template functions. This approach allows to write very generic code, but it is often desirable to write non templated function and still keep some level of genericity to avoid stupid copies of the arguments. The typical example is to write functions accepting both a MatrixXf or a block of a MatrixXf. This exactly the purpose of the Ref class. Here is a simple example:
+
+
+| Example: | Output: |
|---|
+|
+\include function_taking_ref.cpp
+ |
+
+\verbinclude function_taking_ref.out
+ |
+In the first two calls to inv_cond, no copy occur because the memory layout of the arguments matches the memory layout accepted by Ref. However, in the last call, we have a generic expression that will be automatically evaluated into a temporary MatrixXf by the Ref<> object.
+
+A Ref object can also be writable. Here is an example of a function computing the covariance matrix of two input matrices where each row is an observation:
+\code
+void cov(const Ref x, const Ref y, Ref C)
+{
+ const float num_observations = static_cast(x.rows());
+ const RowVectorXf x_mean = x.colwise().sum() / num_observations;
+ const RowVectorXf y_mean = y.colwise().sum() / num_observations;
+ C = (x.rowwise() - x_mean).transpose() * (y.rowwise() - y_mean) / num_observations;
+}
+\endcode
+and here are two examples calling cov without any copy:
+\code
+MatrixXf m1, m2, m3
+cov(m1, m2, m3);
+cov(m1.leftCols<3>(), m2.leftCols<3>(), m3.topLeftCorner<3,3>());
+\endcode
+The Ref<> class has two other optional template arguments allowing to control the kind of memory layout that can be accepted without any copy. See the class Ref documentation for the details.
+
\section TopicPlainFunctionsWorking In which cases do functions taking plain Matrix or Array arguments work?
-Let's assume one wants to write a function computing the covariance matrix of two input matrices where each row is an observation. The implementation of this function might look like this
+Without using template functions, and without the Ref class, a naive implementation of the previous cov function might look like this
\code
MatrixXf cov(const MatrixXf& x, const MatrixXf& y)
{
@@ -88,7 +121,7 @@ MatrixXf cov(const MatrixXf& x, const MatrixXf& y)
return (x.rowwise() - x_mean).transpose() * (y.rowwise() - y_mean) / num_observations;
}
\endcode
-and contrary to what one might think at first, this implementation is fine unless you require a genric implementation that works with double matrices too and unless you do not care about temporary objects. Why is that the case? Where are temporaries involved? How can code as given below compile?
+and contrary to what one might think at first, this implementation is fine unless you require a generic implementation that works with double matrices too and unless you do not care about temporary objects. Why is that the case? Where are temporaries involved? How can code as given below compile?
\code
MatrixXf x,y,z;
MatrixXf C = cov(x,y+z);
@@ -97,6 +130,7 @@ In this special case, the example is fine and will be working because both param
\b Note: Functions taking \e const references to Matrix (or Array) can process expressions at the cost of temporaries.
+
\section TopicPlainFunctionsFailing In which cases do functions taking a plain Matrix or Array argument fail?
Here, we consider a slightly modified version of the function given above. This time, we do not want to return the result but pass an additional non-const paramter which allows us to store the result. A first naive implementation might look as follows.
@@ -149,7 +183,7 @@ MatrixXf y = MatrixXf::Random(100,3);
MatrixXf C;
cov(x, y, C);
\endcode
-This is not the case anymore, when we are using an implementation taking MatrixBase as a parameter. In general, Eigen supports automatic resizing but it is not possible to do so on expressions. Why should resizing of a matrix Block be allowed? It is a reference to a sub-matrix and we definitely don't want to resize that. So how can we incorporate resizing if we cannot resize on MatrixBase? The solution is to resize the derived object as in this implementation.
+This is not the case anymore, when we are using an implementation taking MatrixBase as a parameter. In general, %Eigen supports automatic resizing but it is not possible to do so on expressions. Why should resizing of a matrix Block be allowed? It is a reference to a sub-matrix and we definitely don't want to resize that. So how can we incorporate resizing if we cannot resize on MatrixBase? The solution is to resize the derived object as in this implementation.
\code
template
void cov(const MatrixBase& x, const MatrixBase& y, MatrixBase const & C_)
diff --git a/doc/examples/function_taking_ref.cpp b/doc/examples/function_taking_ref.cpp
new file mode 100644
index 000000000..162a202e4
--- /dev/null
+++ b/doc/examples/function_taking_ref.cpp
@@ -0,0 +1,19 @@
+#include
+#include
+using namespace Eigen;
+using namespace std;
+
+float inv_cond(const Ref& a)
+{
+ const VectorXf sing_vals = a.jacobiSvd().singularValues();
+ return sing_vals(sing_vals.size()-1) / sing_vals(0);
+}
+
+int main()
+{
+ Matrix4f m = Matrix4f::Random();
+ cout << "matrix m:" << endl << m << endl << endl;
+ cout << "inv_cond(m): " << inv_cond(m) << endl;
+ cout << "inv_cond(m(1:3,1:3)): " << inv_cond(m.topLeftCorner(3,3)) << endl;
+ cout << "inv_cond(m+I): " << inv_cond(m+Matrix4f::Identity()) << endl;
+}
|