6. Advanced topics

Aliasing

Aliasing occurs whenever the same Eigen 3 matrix/vector appears on both sides of the assignment operator, and happens because of Eigen 3’s lazy evaluation system. Examples that exhibit aliasing:

mat = 2 * mat;

or

mat = mat.transpose();

Aliasing does not occur in statements like

mat = f(mat);

where f() returns by value. Aliasing produces in general unexpected results, and should be avoided at all costs.

Whereas the first line produces aliasing, it is not dangerous, since the assignment is done in a one-to-one manner, i.e. each element $(i,j)$ on the left hand side of the assignment operator is solely a function of the the same $(i,j)$ element on the right hand side, i.e. $mat(i,j) = f(mat(i,j))$, $\forall i,j$. The problem appears whenever coefficients are being combined and overlap, such as in the second example, where $mat(i,j) = mat(j,i)$, $\forall i,j$. To avoid aliasing, use the Eigen member function eval() to force the evaluation of the right hand side into another (temporary) object, such as

mat = 2 * mat.eval();

In general, aliasing can not be detected at compile time, but can be detected at runtime whenever the macro EIGEN_NO_DEBUG is not set (e.g. in debug mode). I highly recommend first to compile your program in debug mode to detect aliasing run-time assertions, as well as other possible issues that may have escaped you, such as assigning to a matrix another matrix of different dimension etc.

For more details about aliasing, see the official Eigen 3 aliasing documentation.

Type deduction via auto

Avoid the usage of auto when working with Eigen 3 expressions, e.g. avoid writing code like

auto mat = A * B + C;

but write instead

cmat mat = A * B + C;

or

auto mat = (A * B + C).eval();

to force evaluation, as otherwise you may get unexpected results. The “problem" lies in the Eigen 3 lazy evaluation system and reference binding, see e.g. this question on StackOverflow for more details. In short, the reference to the internal data represented by the expression A * B + C is dangling at the end of the auto mat = A * B + C; statement.

Optimizations

Whenever testing your application, I recommend compiling in debug mode, as Eigen 3 run-time assertions can provide extremely helpful feedback on potential issues. Whenever the code is production-ready, you should always compile with optimization flags turned on, such as -O3 (for g++) and -DEIGEN_NO_DEBUG. You should also turn on the OpenMP (if available) multi-processing flag (-fopenmp for g++), as it enables multi-core/multi-processing with shared memory. Eigen 3 uses multi-processing when available, e.g. in matrix multiplication. Quantum++ also uses multi-processing in computationally-intensive functions.

Since most Quantum++ functions return by value, in assignments of the form

mat = f(another_mat);

there is an additional copy assignment operator when assigning the temporary returned by f() back to mat. As far as I know, this extra copy operation is not elided. Unfortunately, Eigen 3 does not yet support move semantics, which would have got rid of this additional assignment via the corresponding move assignment operator. If in the future Eigen 3 will support move semantics, the additional assignment operator will be “free", and you won’t have to modify any existing code to enable the optimization; the Eigen 3 move assignment operator should take care of it for you.

Note that in a line of the form

cmat mat = f(another_mat);

most compilers perform copy elision, i.e. the temporary on the right hand side is constructed directly inside the object mat, the copy constructor being elided.

Extending Quantum++

Most Quantum++ operate on Eigen 3 matrices/vectors, and return either a matrix or a scalar. In principle, you may be tempted to write a new function such as

cmat f(const cmat& A){...}

The problem with the approach above is that Eigen 3 uses expression templates as the type of each expression, i.e. different expressions have in general different types, see the official Eigen 3 documentation for more details. The correct way to write a generic function that is guaranteed to work with any matrix expression is to make your function template and declare the input parameter as Eigen::MatrixBase<Derived>, where Derived is the template parameter. For example, the Quantum++ qpp::transpose() function is defined as

template<typename Derived> 
dyn_mat<typename Derived::Scalar> 
transpose(const Eigen::MatrixBase<Derived>& A)
{
    const dyn_mat<typename Derived::Scalar>& rA = A.derived();

    // check zero-size
    if (!internal::check_nonzero_size(rA))
        throw Exception("qpp::transpose()", Exception::Type::ZERO_SIZE);

    return rA.transpose();
}

It takes an Eigen 3 matrix expression and returns a dynamic matrix over the scalar field of the expression. In the line 

const dyn_mat<typename Derived::Scalar>& rA = A.derived();

we implicitly convert the input expression A to a dynamic matrix rA over the same scalar field as the expression, via binding to a const reference, therefore paying no copying cost. We then use rA instead of the original expression A in the rest of the function. Note that most of the time it is OK to use the original expression, however there are some cases where you may get a compile time error if the expression is not explicitly casted to a matrix. For consistency, I use this reference binding trick in the code of all Quantum++ functions.

As you may have already seen, Quantum++ consists mainly of a collection of functions and few classes. There is no complicated class hierarchy, and you can regard the Quantum++ API as a low/medium-level API. You may extend it to incorporate graphical input, e.g. use a graphical library such as Qt, or build a more sophisticated library on top of it. I recommend to read the source code and make yourself familiar with the library before deciding to extend it. You should also check the official API documentation for an extensive documentation of all functions and classes. I hope you find Quantum++ useful and I wish you a happy usage!