Attachment 173 Details for Bug 275 – Implementation of the Kronecker product for matrices

This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen

Implementation of the Kronecker product for matrices

KroneckerTensorProduct.h (text/x-chdr), 5.30 KB, created by Andreas Platen on 2011-05-20 21:03:28 UTC

(hide)

Description:

Filename:

MIME Type:

Creator: Andreas Platen

Created: 2011-05-20 21:03:28 UTC

Size: 5.30 KB

patch

obsolete

>// 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
>
>
>#include <Eigen/Core>
>
>// Sparse modul is not 100% stable in Eigen 3.0
>#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
>#include <Eigen/Sparse>
>
>
>/*!
> * Computes Kronecker tensor product of two dense matrices
> *
> * \param A   Dense matrix A
> * \param B   Dense matrix B
> * \param AB  Kronecker tensor product of A and B
> */
>template<typename A,typename B,typename C>
>void kroneckerProduct(const Eigen::MatrixBase<A>& a, const Eigen::MatrixBase<B>& b, Eigen::MatrixBase<C>& c) {
>  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 AB  Kronecker tensor product of A and B
> */
>template<typename A,typename B,typename C>
>void kroneckerProduct(const Eigen::MatrixBase<A>& a, const Eigen::SparseMatrixBase<B>& b, Eigen::SparseMatrixBase<C>& c) {
>  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 AB  Kronecker tensor product of A and B
> */
>template<typename A,typename B,typename C>
>void kroneckerProduct(const Eigen::SparseMatrixBase<A>& a, const Eigen::MatrixBase<B>& b, Eigen::SparseMatrixBase<C>& c) {
>  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 AB  Kronecker tensor product of A and B
> */
>template<typename A,typename B,typename C>
>void kroneckerProduct(const Eigen::SparseMatrixBase<A>& a, const Eigen::SparseMatrixBase<B>& b, Eigen::SparseMatrixBase<C>& c) {
>  kroneckerProduct_sparse(a.derived(), b.derived(), c.derived());
>}
>
>
>/*!
> * 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, Eigen::MatrixBase<Derived_AB> const & AB_)
>{
>    const unsigned int Ar = A.rows(),
>                       Ac = A.cols(),
>                       Br = B.rows(),
>                       Bc = B.cols();
>    Eigen::MatrixBase<Derived_AB>& AB = const_cast< Eigen::MatrixBase<Derived_AB>& >(AB_);
>    AB.derived().resize(Ar*Br,Ac*Bc);
>
>    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, Eigen::SparseMatrixBase<Derived_AB> const &AB_)
>{
>    const unsigned int Ar = A.rows(),
>                       Ac = A.cols(),
>                       Br = B.rows(),
>                       Bc = B.cols();
>
>    Eigen::SparseMatrixBase<Derived_AB>& AB = const_cast< Eigen::SparseMatrixBase<Derived_AB>& >(AB_);
>    AB.derived().resize(Ar*Br,Ac*Bc);
>    AB.derived().resizeNonZeros(0);
>    AB.derived().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.derived(),kA); itA; ++itA)
>            {
>                for (typename Derived_B::InnerIterator itB(B.derived(),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.derived().insert(i,j) = itA.value() * itB.value();
>                }
>            }
>        }
>    }
>}
>
>
>#endif // KRONECKER_TENSOR_PRODUCT_H

Actions: View

Attachments on bug 275: 173 | 176 | 177 | 185 | 186