-
Notifications
You must be signed in to change notification settings - Fork 116
6. Advanced topics
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 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). We 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 similar to
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 Stack Overflow
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
CPP auto mat = A * B + C;
statement.
Whenever testing your application, we 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 we know, this extra copy operation
is not elided. If your version of Eigen 3
supports move semantics, the additional assignment operator is for “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.
Most Quantum++ functions 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();
// EXCEPTION CHECKS
// check zero-size
if (!internal::check_nonzero_size(rA))
throw Exception::ZeroSize("qpp::transpose()", "A");
// END EXCEPTION CHECKS
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, we use this reference binding trick all over
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. We 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. We hope you find Quantum++ useful and we wish you a happy usage!
Copyright (c) 2017 - 2024 softwareQ Inc. All rights reserved.