diff options
author | Gael Guennebaud <g.gael@free.fr> | 2011-06-22 14:39:11 +0200 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2011-06-22 14:39:11 +0200 |
commit | 3ecf7e8f6ea95e1dbf9c65ff6b44fd543fcd0de7 (patch) | |
tree | 70e40248fd0465413d73eb98700c21a4f7300993 /unsupported/Eigen | |
parent | 7aabce7c769adc0ad02e45e7ffe7590c5d1412d1 (diff) |
add a KroneckerProduct module (unsupported) from Kolja Brix and Andreas Platen materials.
Diffstat (limited to 'unsupported/Eigen')
-rw-r--r-- | unsupported/Eigen/CMakeLists.txt | 2 | ||||
-rw-r--r-- | unsupported/Eigen/KroneckerProduct | 26 | ||||
-rw-r--r-- | unsupported/Eigen/src/CMakeLists.txt | 1 | ||||
-rw-r--r-- | unsupported/Eigen/src/KroneckerProduct/CMakeLists.txt | 6 | ||||
-rw-r--r-- | unsupported/Eigen/src/KroneckerProduct/KroneckerTensorProduct.h | 168 |
5 files changed, 202 insertions, 1 deletions
diff --git a/unsupported/Eigen/CMakeLists.txt b/unsupported/Eigen/CMakeLists.txt index 219ec8ead..1618525d5 100644 --- a/unsupported/Eigen/CMakeLists.txt +++ b/unsupported/Eigen/CMakeLists.txt @@ -1,6 +1,6 @@ set(Eigen_HEADERS AdolcForward BVH IterativeSolvers MatrixFunctions MoreVectorization AutoDiff AlignedVector3 Polynomials CholmodSupport FFT NonLinearOptimization SparseExtra SuperLUSupport UmfPackSupport IterativeSolvers - NumericalDiff Skyline MPRealSupport OpenGLSupport + NumericalDiff Skyline MPRealSupport OpenGLSupport KroneckerProduct ) install(FILES diff --git a/unsupported/Eigen/KroneckerProduct b/unsupported/Eigen/KroneckerProduct new file mode 100644 index 000000000..291c7fabc --- /dev/null +++ b/unsupported/Eigen/KroneckerProduct @@ -0,0 +1,26 @@ +#ifndef EIGEN_KRONECKER_PRODUCT_MODULE_H +#define EIGEN_KRONECKER_PRODUCT_MODULE_H + +#include "../../Eigen/Core" + +#include "../../Eigen/src/Core/util/DisableStupidWarnings.h" + +namespace Eigen { + +/** \ingroup Unsupported_modules + * \defgroup KroneckerProduct_Module KroneckerProduct module + * + * This module contains an experimental Kronecker product implementation. + * + * \code + * #include <Eigen/KroneckerProduct> + * \endcode + */ + +#include "src/KroneckerProduct/KroneckerTensorProduct.h" + +} // namespace Eigen + +#include "../../Eigen/src/Core/util/ReenableStupidWarnings.h" + +#endif // EIGEN_KRONECKER_PRODUCT_MODULE_H diff --git a/unsupported/Eigen/src/CMakeLists.txt b/unsupported/Eigen/src/CMakeLists.txt index cd442cefa..6fd98b310 100644 --- a/unsupported/Eigen/src/CMakeLists.txt +++ b/unsupported/Eigen/src/CMakeLists.txt @@ -9,3 +9,4 @@ ADD_SUBDIRECTORY(NumericalDiff) ADD_SUBDIRECTORY(Polynomials) ADD_SUBDIRECTORY(Skyline) ADD_SUBDIRECTORY(SparseExtra) +ADD_SUBDIRECTORY(KroneckerProduct)
\ No newline at end of file diff --git a/unsupported/Eigen/src/KroneckerProduct/CMakeLists.txt b/unsupported/Eigen/src/KroneckerProduct/CMakeLists.txt new file mode 100644 index 000000000..4daefebee --- /dev/null +++ b/unsupported/Eigen/src/KroneckerProduct/CMakeLists.txt @@ -0,0 +1,6 @@ +FILE(GLOB Eigen_KroneckerProduct_SRCS "*.h") + +INSTALL(FILES + ${Eigen_KroneckerProduct_SRCS} + DESTINATION ${INCLUDE_INSTALL_DIR}/unsupported/Eigen/src/KroneckerProduct COMPONENT Devel + ) diff --git a/unsupported/Eigen/src/KroneckerProduct/KroneckerTensorProduct.h b/unsupported/Eigen/src/KroneckerProduct/KroneckerTensorProduct.h new file mode 100644 index 000000000..cc5f271f0 --- /dev/null +++ b/unsupported/Eigen/src/KroneckerProduct/KroneckerTensorProduct.h @@ -0,0 +1,168 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2011 Kolja Brix <brix@igpm.rwth-aachen.de> +// Copyright (C) 2011 Andreas Platen <andiplaten@gmx.de> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + + +#ifndef KRONECKER_TENSOR_PRODUCT_H +#define KRONECKER_TENSOR_PRODUCT_H + + +namespace internal { + +/*! + * Kronecker tensor product helper function for dense matrices + * + * \param A Dense matrix A + * \param B Dense matrix B + * \param AB_ Kronecker tensor product of A and B + */ +template<typename Derived_A, typename Derived_B, typename Derived_AB> +void kroneckerProduct_full(const Derived_A& A, const Derived_B& B, Derived_AB & AB) +{ + const unsigned int Ar = A.rows(), + Ac = A.cols(), + Br = B.rows(), + Bc = B.cols(); + for (unsigned int i=0; i<Ar; ++i) + for (unsigned int j=0; j<Ac; ++j) + AB.block(i*Br,j*Bc,Br,Bc) = A(i,j)*B; +} + + +/*! + * Kronecker tensor product helper function for matrices, where at least one is sparse + * + * \param A Matrix A + * \param B Matrix B + * \param AB_ Kronecker tensor product of A and B + */ +template<typename Derived_A, typename Derived_B, typename Derived_AB> +void kroneckerProduct_sparse(const Derived_A &A, const Derived_B &B, Derived_AB &AB) +{ + const unsigned int Ar = A.rows(), + Ac = A.cols(), + Br = B.rows(), + Bc = B.cols(); + AB.resize(Ar*Br,Ac*Bc); + AB.resizeNonZeros(0); + AB.reserve(A.nonZeros()*B.nonZeros()); + + for (int kA=0; kA<A.outerSize(); ++kA) + { + for (int kB=0; kB<B.outerSize(); ++kB) + { + for (typename Derived_A::InnerIterator itA(A,kA); itA; ++itA) + { + for (typename Derived_B::InnerIterator itB(B,kB); itB; ++itB) + { + const unsigned int iA = itA.row(), + jA = itA.col(), + iB = itB.row(), + jB = itB.col(), + i = iA*Br + iB, + j = jA*Bc + jB; + AB.insert(i,j) = itA.value() * itB.value(); + } + } + } + } +} + +} // end namespace internal + + + +/*! + * Computes Kronecker tensor product of two dense matrices + * + * \param a Dense matrix a + * \param b Dense matrix b + * \param c Kronecker tensor product of a and b + */ +template<typename A,typename B,typename CScalar,int CRows,int CCols, int COptions, int CMaxRows, int CMaxCols> +void kroneckerProduct(const MatrixBase<A>& a, const MatrixBase<B>& b, Matrix<CScalar,CRows,CCols,COptions,CMaxRows,CMaxCols>& c) +{ + c.resize(a.rows()*b.rows(),a.cols()*b.cols()); + internal::kroneckerProduct_full(a.derived(), b.derived(), c); +} + +/*! + * Computes Kronecker tensor product of two dense matrices + * + * Remark: this function uses the const cast hack and has been + * implemented to make the function call possible, where the + * output matrix is a submatrix, e.g. + * kroneckerProduct(A,B,AB.block(2,5,6,6)); + * + * \param a Dense matrix a + * \param b Dense matrix b + * \param c Kronecker tensor product of a and b + */ +template<typename A,typename B,typename C> +void kroneckerProduct(const MatrixBase<A>& a, const MatrixBase<B>& b, MatrixBase<C> const & c_) +{ + MatrixBase<C>& c = const_cast<MatrixBase<C>& >(c_); + internal::kroneckerProduct_full(a.derived(), b.derived(), c.derived()); +} + +/*! + * Computes Kronecker tensor product of a dense and a sparse matrix + * + * \param a Dense matrix a + * \param b Sparse matrix b + * \param c Kronecker tensor product of a and b + */ +template<typename A,typename B,typename C> +void kroneckerProduct(const MatrixBase<A>& a, const SparseMatrixBase<B>& b, SparseMatrixBase<C>& c) +{ + internal::kroneckerProduct_sparse(a.derived(), b.derived(), c.derived()); +} + +/*! + * Computes Kronecker tensor product of a sparse and a dense matrix + * + * \param a Sparse matrix a + * \param b Dense matrix b + * \param c Kronecker tensor product of a and b + */ +template<typename A,typename B,typename C> +void kroneckerProduct(const SparseMatrixBase<A>& a, const MatrixBase<B>& b, SparseMatrixBase<C>& c) +{ + internal::kroneckerProduct_sparse(a.derived(), b.derived(), c.derived()); +} + +/*! + * Computes Kronecker tensor product of two sparse matrices + * + * \param a Sparse matrix a + * \param b Sparse matrix b + * \param c Kronecker tensor product of a and b + */ +template<typename A,typename B,typename C> +void kroneckerProduct(const SparseMatrixBase<A>& a, const SparseMatrixBase<B>& b, SparseMatrixBase<C>& c) +{ + internal::kroneckerProduct_sparse(a.derived(), b.derived(), c.derived()); +} + +#endif // KRONECKER_TENSOR_PRODUCT_H |