This commit is contained in:
Benoit Jacob 2010-01-06 08:24:42 -05:00
commit e77748ef96
6 changed files with 209 additions and 99 deletions

View File

@ -450,6 +450,12 @@ class SparseMatrix
return *this;
}
template<typename Lhs, typename Rhs>
inline SparseMatrix& operator=(const SparseProduct<Lhs,Rhs>& product)
{
return Base::operator=(product);
}
template<typename OtherDerived>
EIGEN_DONT_INLINE SparseMatrix& operator=(const SparseMatrixBase<OtherDerived>& other)
{

View File

@ -147,13 +147,16 @@ static void ei_sparse_product_impl2(const Lhs& lhs, const Rhs& rhs, ResultType&
ei_assert(lhs.outerSize() == rhs.innerSize());
std::vector<bool> mask(rows,false);
Matrix<Scalar,Dynamic,1> values(rows);
Matrix<int,Dynamic,1> indices(rows);
// estimate the number of non zero entries
float ratioLhs = float(lhs.nonZeros())/(float(lhs.rows())*float(lhs.cols()));
float avgNnzPerRhsColumn = float(rhs.nonZeros())/float(cols);
float ratioRes = std::min(ratioLhs * avgNnzPerRhsColumn, 1.f);
float ratio;
int t200 = rows/(log2(200)*1.39);
int t = (rows*100)/139;
res.resize(rows, cols);
res.reserve(int(ratioRes*rows*cols));
@ -162,6 +165,7 @@ static void ei_sparse_product_impl2(const Lhs& lhs, const Rhs& rhs, ResultType&
{
res.startVec(j);
int nnz = 0;
for (typename Rhs::InnerIterator rhsIt(rhs, j); rhsIt; ++rhsIt)
{
Scalar y = rhsIt.value();
@ -173,42 +177,42 @@ static void ei_sparse_product_impl2(const Lhs& lhs, const Rhs& rhs, ResultType&
if(!mask[i])
{
mask[i] = true;
values[i] = x * y;
res.insertBackNoCheck(j,i);
// values[i] = x * y;
// indices[nnz] = i;
++nnz;
}
else
res._valuePtr()[mask[i]] += x* y;
values[i] += x * y;
}
}
// FIXME reserve nnz non zeros
// FIXME implement fast sort algorithms for very small nnz
// if the result is sparse enough => use a quick sort
// otherwise => loop through the entire vector
SparseInnerVectorSet<ResultType,1> vec(res,j);
int nnz = vec.nonZeros();
if(rows/1.39 > nnz * log2(nnz))
{
std::sort(vec._innerIndexPtr(), vec._innerIndexPtr()+vec.nonZeros());
for (typename ResultType::InnerIterator it(res, j); it; ++it)
{
it.valueRef() = values[it.index()];
mask[it.index()] = false;
}
}
else
{
// dense path
int count = 0;
for(int i=0; i<rows; ++i)
{
if(mask[i])
{
mask[i] = false;
vec._innerIndexPtr()[count] = i;
vec._valuePtr()[count] = i;
++count;
}
}
}
// In order to avoid to perform an expensive log2 when the
// result is clearly very sparse we use a linear bound up to 200.
// if((nnz<200 && nnz<t200) || nnz * log2(nnz) < t)
// {
// if(nnz>1) std::sort(indices.data(),indices.data()+nnz);
// for(int k=0; k<nnz; ++k)
// {
// int i = indices[k];
// res.insertBackNoCheck(j,i) = values[i];
// mask[i] = false;
// }
// }
// else
// {
// // dense path
// for(int i=0; i<rows; ++i)
// {
// if(mask[i])
// {
// mask[i] = false;
// res.insertBackNoCheck(j,i) = values[i];
// }
// }
// }
}
res.finalize();
@ -218,6 +222,8 @@ static void ei_sparse_product_impl2(const Lhs& lhs, const Rhs& rhs, ResultType&
template<typename Lhs, typename Rhs, typename ResultType>
static void ei_sparse_product_impl(const Lhs& lhs, const Rhs& rhs, ResultType& res)
{
// return ei_sparse_product_impl2(lhs,rhs,res);
typedef typename ei_traits<typename ei_cleantype<Lhs>::type>::Scalar Scalar;
// make sure to call innerSize/outerSize since we fake the storage order.
@ -274,6 +280,7 @@ struct ei_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,ColMajor,ColMajor>
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
{
// std::cerr << __LINE__ << "\n";
typename ei_cleantype<ResultType>::type _res(res.rows(), res.cols());
ei_sparse_product_impl<Lhs,Rhs,ResultType>(lhs, rhs, _res);
res.swap(_res);
@ -285,6 +292,7 @@ struct ei_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,ColMajor,RowMajor>
{
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
{
// std::cerr << __LINE__ << "\n";
// we need a col-major matrix to hold the result
typedef SparseMatrix<typename ResultType::Scalar> SparseTemporaryType;
SparseTemporaryType _res(res.rows(), res.cols());
@ -298,6 +306,7 @@ struct ei_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,RowMajor,RowMajor>
{
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
{
// std::cerr << __LINE__ << "\n";
// let's transpose the product to get a column x column product
typename ei_cleantype<ResultType>::type _res(res.rows(), res.cols());
ei_sparse_product_impl<Rhs,Lhs,ResultType>(rhs, lhs, _res);
@ -310,11 +319,20 @@ struct ei_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,RowMajor,ColMajor>
{
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
{
// std::cerr << "here...\n";
typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix;
ColMajorMatrix colLhs(lhs);
ColMajorMatrix colRhs(rhs);
// std::cerr << "more...\n";
ei_sparse_product_impl<ColMajorMatrix,ColMajorMatrix,ResultType>(colLhs, colRhs, res);
// std::cerr << "OK.\n";
// let's transpose the product to get a column x column product
typedef SparseMatrix<typename ResultType::Scalar> SparseTemporaryType;
SparseTemporaryType _res(res.cols(), res.rows());
ei_sparse_product_impl<Rhs,Lhs,SparseTemporaryType>(rhs, lhs, _res);
res = _res.transpose();
// typedef SparseMatrix<typename ResultType::Scalar> SparseTemporaryType;
// SparseTemporaryType _res(res.cols(), res.rows());
// ei_sparse_product_impl<Rhs,Lhs,SparseTemporaryType>(rhs, lhs, _res);
// res = _res.transpose();
}
};
@ -327,6 +345,7 @@ template<typename Derived>
template<typename Lhs, typename Rhs>
inline Derived& SparseMatrixBase<Derived>::operator=(const SparseProduct<Lhs,Rhs>& product)
{
// std::cerr << "there..." << typeid(Lhs).name() << " " << typeid(Lhs).name() << " " << (Derived::Flags&&RowMajorBit) << "\n";
ei_sparse_product_selector<
typename ei_cleantype<Lhs>::type,
typename ei_cleantype<Rhs>::type,
@ -348,7 +367,7 @@ struct ei_sparse_product_selector2<Lhs,Rhs,ResultType,ColMajor,ColMajor,ColMajor
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
{
ei_sparse_product_impl2<Lhs,Rhs,ResultType>(lhs, rhs, res, 0);
ei_sparse_product_impl2<Lhs,Rhs,ResultType>(lhs, rhs, res);
}
};
@ -357,11 +376,11 @@ struct ei_sparse_product_selector2<Lhs,Rhs,ResultType,RowMajor,ColMajor,ColMajor
{
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
{
typedef SparseMatrix<typename ResultType::Scalar,RowMajor> RowMajorMatrix;
RowMajorMatrix rhsRow = rhs;
RowMajorMatrix resRow(res.rows(), res.cols());
ei_sparse_product_impl2<RowMajorMatrix,Lhs,RowMajorMatrix>(rhsRow, lhs, resRow);
res = resRow;
// typedef SparseMatrix<typename ResultType::Scalar,RowMajor> RowMajorMatrix;
// RowMajorMatrix rhsRow = rhs;
// RowMajorMatrix resRow(res.rows(), res.cols());
// ei_sparse_product_impl2<RowMajorMatrix,Lhs,RowMajorMatrix>(rhsRow, lhs, resRow);
// res = resRow;
}
};

View File

@ -0,0 +1,73 @@
namespace Eigen {
/** \page Eigen2ToEigen3 Porting from Eigen2 to Eigen3
The goals of this page is to enumerate the API changes between Eigen2 and Eigen3,
and to help porting an application from Eigen2 to Eigen3.
\b Table \b of \b contents
- \ref summary
- \ref modules
- \ref core
\section CompatibilitySupport Eigen2 compatibility support
In order to ease th switch from Eigen2 to Eigen3, Eigen3 features a compatibility mode which can be enabled by defining the EIGEN2_SUPPORT preprocessor token \b before including any Eigen's header (typically it should be set in your project options).
\section ChangeList List of changes
<table>
<tr><td>Eigen 2</td><td>Eigen 3</td></tr>
<tr><td>vec.start(length)</td><td>vec.head(length)</td></tr>
<tr><td>vec.start<length>()</td><td>vec.head<length>()</td></tr>
<tr><td>vec.end(length)</td><td>vec.tail(length)</td></tr>
<tr><td>vec.end<length>()</td><td>vec.tail<length>()</td></tr>
<tr></tr>
<tr><td>mat.cwise().XXX()</td><td>mat.array().XXX() \redstar</td></tr>
<tr></tr>
<tr><td>\code c = (a * b).lazy();\endcode</td><td>\code c.noalias() = a * b;\endcode</td></tr>
</table>
\redstar See \ref CoefficientWiseOperations for details.
\section CoefficientWiseOperations Coefficient wise operations
In Eigen2, coefficient wise operations which have no proper mathematical definiton (as a coeff wise product)
were achied using the .cwise() prefix, e.g.:
\code a.cwise() * b \code
In Eigen3 this .cwise() prefix has been supersed by a new kind of matrix type called
Array for which all operations are performed coefficient wise. You can easily view a matrix as an array and vice versa using
the MatrixBase::array() and ArrayBase::matrix() functions respectively. Here is an example:
\code
Vector4f a, b, c;
c = a.array() * b.array();
\endcode
Note that the .array() function is not at all a synonym of the deprecated .cwise() prefix.
While the .cwise() prefix changed the behavior of the following operator, the array() function performs
a permanent convertion to the array world. Therefore, for binary operations such as the coeff wise product,
both sides must be converted to an \em array as in the above example. On the other hand, when you
concatenates multiple coeff wise operations you only have to do the conversion once, e.g.:
\code
Vector4f a, b, c;
c = a.array().abs().pow(3) * b.array().abs().sin();
\endcode
With Eigen2 you would have written:
\code
c = (a.cwise().abs().cwise().pow(3)).cwise() * (b.cwise().abs().cwise().sin());
\endcode
\section LazyVsNoalias Lazy evaluation versus noalias
In Eigen all operations are performed in a lazy fashion expected the matrix products which are always evaluated to a temporary by default.
In Eigen2, lazy evaluation could be enforced by tagging a product using the .lazy() function. However, in complex expressions it was not
easy to determine where to put the lazy() function. In Eigen3, the lazy() feature has been supersed by the MatrixBase::noalias() function
which can be used on the left hand side of an assignment when no aliasing can occur. Here is an example:
\code
MatrixXf a, b, c;
...
c.noalias() += 2 * a.transpose() * b;
\endcode
*/
}

View File

@ -39,10 +39,10 @@ A.setIdentity(); // Fill A with the identity.
// Templated size versions are faster. Note that Matlab is 1-based (a size N
// vector is x(1)...x(N)).
// Eigen // Matlab
x.start(n) // x(1:n)
x.start<n>() // x(1:n)
x.tail(n) // N = rows(x); x(N - n: N)
x.tail<n>() // N = rows(x); x(N - n: N)
x.head(n) // x(1:n)
x.head<n>() // x(1:n)
x.tail(n) // N = rows(x); x(N - n: N)
x.tail<n>() // N = rows(x); x(N - n: N)
x.segment(i, n) // x(i+1 : i+n)
x.segment<n>(i) // x(i+1 : i+n)
P.block(i, j, rows, cols) // P(i+1 : i+rows, j+1 : j+cols)
@ -81,30 +81,36 @@ a *= M; R = P + Q; R = P/s;
// Vectorized operations on each element independently
// (most require #include <Eigen/Array>)
// Eigen // Matlab
R = P.cwise() * Q; // R = P .* Q
R = P.cwise() / Q; // R = P ./ Q
R = P.cwise() + s; // R = P + s
R = P.cwise() - s; // R = P - s
R.cwise() += s; // R = R + s
R.cwise() -= s; // R = R - s
R.cwise() *= s; // R = R * s
R.cwise() /= s; // R = R / s
R.cwise() < Q; // R < Q
R.cwise() <= Q; // R <= Q
R.cwise().inverse(); // 1 ./ P
R.cwise().sin() // sin(P)
R.cwise().cos() // cos(P)
R.cwise().pow(s) // P .^ s
R.cwise().square() // P .^ 2
R.cwise().cube() // P .^ 3
R.cwise().sqrt() // sqrt(P)
R.cwise().exp() // exp(P)
R.cwise().log() // log(P)
R.cwise().max(P) // max(R, P)
R.cwise().min(P) // min(R, P)
R.cwise().abs() // abs(P)
R.cwise().abs2() // abs(P.^2)
(R.cwise() < s).select(P,Q); // (R < s ? P : Q)
R = P.cwiseProduct(Q); // R = P .* Q
R = P.array() * s.array();// R = P .* s
R = P.cwiseQuotient(Q); // R = P ./ Q
R = P.array() / Q.array();// R = P ./ Q
R = P.array() + s.array();// R = P + s
R = P.array() - s.array();// R = P - s
R.array() += s; // R = R + s
R.array() -= s; // R = R - s
R.array() < Q.array(); // R < Q
R.array() <= Q.array(); // R <= Q
R.cwiseInverse(); // 1 ./ P
R.array().inverse(); // 1 ./ P
R.array().sin() // sin(P)
R.array().cos() // cos(P)
R.array().pow(s) // P .^ s
R.array().square() // P .^ 2
R.array().cube() // P .^ 3
R.cwiseSqrt() // sqrt(P)
R.array().sqrt() // sqrt(P)
R.array().exp() // exp(P)
R.array().log() // log(P)
R.cwiseMax(P) // max(R, P)
R.array().max(P.array()) // max(R, P)
R.cwiseMin(P) // min(R, P)
R.array().min(P.array()) // min(R, P)
R.cwiseAbs() // abs(P)
R.array().abs() // abs(P)
R.cwiseAbs2() // abs(P.^2)
R.array().abs2() // abs(P.^2)
(R.array() < s).select(P,Q); // (R < s ? P : Q)
// Reductions.
int r, c;

View File

@ -40,7 +40,7 @@ Here are some quick compilation instructions with GCC. To quickly test an exampl
There is no library to link to. For good performance, add the \c -O2 compile-flag. Note however that this makes it impossible to debug inside Eigen code, as many functions get inlined. In some cases, performance can be further improved by disabling Eigen assertions: use \c -DEIGEN_NO_DEBUG or \c -DNDEBUG to disable them.
On the x86 architecture, the SSE2 instruction set is not enabled by default. Use \c -msse2 to enable it, and Eigen will then automatically enable its vectorized paths. On x86-64 and AltiVec-based architectures, vectorization is enabled by default.
On the x86 architecture, the SSE2 instruction set is not enabled by default. Use \c -msse2 to enable it, and Eigen will then automatically enable its vectorized paths. On x86-64 and AltiVec-based architectures, vectorization is enabled by default. To turn off Eigen's vectorization use \c-DEIGEN_DONT_VECTORIZE.
\section TutorialCoreSimpleExampleFixedSize Simple example with fixed-size matrices and vectors
@ -54,7 +54,8 @@ output:
\include Tutorial_simple_example_fixed_size.out
</td></tr></table>
<a href="#" class="top">top</a>\section TutorialCoreSimpleExampleDynamicSize Simple example with dynamic-size matrices and vectors
<a href="#" class="top">top</a>
\section TutorialCoreSimpleExampleDynamicSize Simple example with dynamic-size matrices and vectors
By dynamic-size, we mean that the numbers of rows and columns are not fixed at compile-time. In this case, they are stored as runtime variables and the arrays are dynamically allocated.
@ -74,11 +75,14 @@ output:
This slows compilation down but at least you don't have to worry anymore about including the correct files! There also is the Eigen/Dense header including all dense functionality i.e. leaving out the Sparse module.
<a href="#" class="top">top</a>\section TutorialCoreMatrixTypes Matrix and vector types
<a href="#" class="top">top</a>
\section TutorialCoreMatrixTypes Array, matrix and vector types
In Eigen, all kinds of dense matrices and vectors are represented by the template class Matrix. In most cases, you can simply use one of the \ref matrixtypedefs "convenience typedefs".
Eigen provides two kinds of dense objects: mathematical matrices and vectors which are both represented by the template class Matrix, and 1D and 2D arrays represented by the template class Array. While the former (Matrix) is specialized for the representation of mathematical objects, the latter (Array) represents a collection of scalar values arranged in a 1D or 2D fashion. In particular, all operations performed on arrays are coefficient wise. Conversion between the two worlds can be done using the MatrixBase::array() and ArrayBase::matrix() functions respectively without any overhead. See \ref TutorialCoreArithmeticOperators for further details.
The template class Matrix takes a number of template parameters, but for now it is enough to understand the 3 first ones (and the others can then be left unspecified):
In most cases, you can simply use one of the \ref matrixtypedefs "convenience typedefs".
The template class Matrix, just like the class Array) take a number of template parameters, but for now it is enough to understand the 3 first ones (and the others can then be left unspecified):
\code Matrix<Scalar, RowsAtCompileTime, ColsAtCompileTime> \endcode
@ -368,24 +372,24 @@ most common coefficient-wise operators.
<tr><td>
Add a scalar to all coefficients \redstar</td><td>\code
mat3 = mat1.cwise() + scalar;
mat3.cwise() += scalar;
mat3.cwise() -= scalar;
mat3.array() += scalar;
mat3.array() -= scalar;
\endcode
</td></tr>
<tr><td>
Coefficient wise \link Cwise::operator/() division \endlink \redstar</td><td>\code
mat3 = mat1.cwise() / mat2; \endcode
mat3 = mat1.array() / mat2.array(); \endcode
</td></tr>
<tr><td>
Coefficient wise \link Cwise::inverse() reciprocal \endlink \redstar</td><td>\code
mat3 = mat1.cwise().inverse(); \endcode
mat3 = mat1.array().inverse(); \endcode
</td></tr>
<tr><td>
Coefficient wise comparisons \redstar \n
(support all operators)</td><td>\code
mat3 = mat1.cwise() < mat2;
mat3 = mat1.cwise() <= mat2;
mat3 = mat1.cwise() > mat2;
mat3 = mat1.array() < mat2.array();
mat3 = mat1.array() <= mat2.array();
mat3 = mat1.array() > mat2.array();
etc.
\endcode
</td></tr></table>
@ -394,20 +398,20 @@ etc.
<tr><td>
\b Trigo \redstar: \n
\link Cwise::sin sin \endlink, \link Cwise::cos cos \endlink</td><td>\code
mat3 = mat1.cwise().sin();
mat3 = mat1.array().sin();
etc.
\endcode
</td></tr>
<tr><td>
\b Power \redstar: \n \link Cwise::pow() pow \endlink,
\link Cwise::square square \endlink,
\link Cwise::cube cube \endlink, \n
\link Cwise::sqrt sqrt \endlink,
\link Cwise::exp exp \endlink,
\link Cwise::log log \endlink </td><td>\code
mat3 = mat1.cwise().square();
mat3 = mat1.cwise().pow(5);
mat3 = mat1.cwise().log();
\link ArrayBase::square square \endlink,
\link ArrayBase::cube cube \endlink, \n
\link ArrayBase::sqrt sqrt \endlink,
\link ArrayBase::exp exp \endlink,
\link ArrayBase::log log \endlink </td><td>\code
mat3 = mat1.array().square();
mat3 = mat1.array().pow(5);
mat3 = mat1.array().log();
etc.
\endcode
</td></tr>
@ -415,10 +419,10 @@ etc.
\link Cwise::min min \endlink, \link Cwise::max max \endlink, \n
absolute value (\link Cwise::abs() abs \endlink, \link Cwise::abs2() abs2 \endlink)
</td><td>\code
mat3 = mat1.cwise().min(mat2);
mat3 = mat1.cwise().max(mat2);
mat3 = mat1.cwise().abs();
mat3 = mat1.cwise().abs2();
mat3 = mat1.cwiseMin(mat2);
mat3 = mat1.cwiseMax(mat2);
mat3 = mat1.cwiseAbs();
mat3 = mat1.cwiseAbs2();
\endcode</td></tr>
</table>
</td></tr></table>
@ -492,7 +496,7 @@ Read-write access to sub-vectors:
<td>Optimized versions when the size \n is known at compile time</td></tr>
<td></td>
<tr><td>\code vec1.start(n)\endcode</td><td>\code vec1.start<n>()\endcode</td><td>the first \c n coeffs </td></tr>
<tr><td>\code vec1.head(n)\endcode</td><td>\code vec1.head<n>()\endcode</td><td>the first \c n coeffs </td></tr>
<tr><td>\code vec1.tail(n)\endcode</td><td>\code vec1.tail<n>()\endcode</td><td>the last \c n coeffs </td></tr>
<tr><td>\code vec1.segment(pos,n)\endcode</td><td>\code vec1.segment<n>(pos)\endcode</td>
<td>the \c size coeffs in \n the range [\c pos : \c pos + \c n [</td></tr>

View File

@ -9,15 +9,17 @@ o /** \mainpage Eigen
| \ref TutorialSparse "Sparse matrix"
</div>
This is the API documentation for Eigen.
This is the API documentation for Eigen3.
You come from Eigen2? Here is a \ref Eigen2ToEigen3 guide for porting your application from Eigen2 to Eigen3.
For a first contact with Eigen, the best place is to have a look at the \ref TutorialCore "tutorial". For an even shorter overview, we have an <a href="AsciiQuickReference.txt">ASCII quick reference</a> with Matlab translations.
Most of the API is available as methods in MatrixBase, so this is a good starting point for browsing. Also have a look at Matrix, as a few methods and the matrix constructors are there. Other notable classes for the Eigen API are Cwise, which contains the methods for doing certain coefficient-wise operations, and Part.
Most of the API is available as methods in DenseBase and MatrixBase, so this is a good starting point for browsing. Also have a look at Matrix, as a few methods and the matrix constructors are there. Other notable classes for the Eigen API are ArrayBase, which contains the methods for doing certain coefficient-wise operations, and TriangularView.
In fact, except for advanced use, the only class that you'll have to explicitly name in your program, i.e. of which you'll explicitly contruct objects, is Matrix. For instance, vectors are handled as a special case of Matrix with one column. Typedefs are provided, e.g. Vector2f is a typedef for Matrix<float, 2, 1>.
In fact, except for advanced use, the only classes that you'll have to explicitly name in your program, i.e. of which you'll explicitly contruct objects, is Matrix and Array. For instance, vectors are handled as a special case of Matrix with one column. Typedefs are provided, e.g. Vector2f is a typedef for Matrix<float, 2, 1>.
Most of the other classes are just return types for MatrixBase methods.
Most of the other classes are just return types for DenseBase and MatrixBase methods.
Want more? Checkout the \ref Unsupported_modules "unsupported modules" <a href="unsupported/index.html">documentation</a>.