// This file is part of Eigen, a lightweight C++ template library // for linear algebra. // // Copyright (C) 2011 Kolja Brix // Copyright (C) 2011 Andreas Platen // Copyright (C) 2012 Chen-Pang He // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. #ifdef EIGEN_TEST_PART_1 #include "sparse.h" #include #include template void check_dimension(const MatrixType& ab, const int rows, const int cols) { VERIFY_IS_EQUAL(ab.rows(), rows); VERIFY_IS_EQUAL(ab.cols(), cols); } template void check_kronecker_product(const MatrixType& ab) { VERIFY_IS_EQUAL(ab.rows(), 6); VERIFY_IS_EQUAL(ab.cols(), 6); VERIFY_IS_EQUAL(ab.size(), 36); VERIFY_IS_APPROX(ab.coeff(0, 0), -0.4017367630386106); VERIFY_IS_APPROX(ab.coeff(0, 1), 0.1056863433932735); VERIFY_IS_APPROX(ab.coeff(0, 2), -0.7255206194554212); VERIFY_IS_APPROX(ab.coeff(0, 3), 0.1908653336744706); VERIFY_IS_APPROX(ab.coeff(0, 4), 0.350864567234111); VERIFY_IS_APPROX(ab.coeff(0, 5), -0.0923032108308013); VERIFY_IS_APPROX(ab.coeff(1, 0), 0.415417514804677); VERIFY_IS_APPROX(ab.coeff(1, 1), -0.2369227701722048); VERIFY_IS_APPROX(ab.coeff(1, 2), 0.7502275131458511); VERIFY_IS_APPROX(ab.coeff(1, 3), -0.4278731019742696); VERIFY_IS_APPROX(ab.coeff(1, 4), -0.3628129162264507); VERIFY_IS_APPROX(ab.coeff(1, 5), 0.2069210808481275); VERIFY_IS_APPROX(ab.coeff(2, 0), 0.05465890160863986); VERIFY_IS_APPROX(ab.coeff(2, 1), -0.2634092511419858); VERIFY_IS_APPROX(ab.coeff(2, 2), 0.09871180285793758); VERIFY_IS_APPROX(ab.coeff(2, 3), -0.4757066334017702); VERIFY_IS_APPROX(ab.coeff(2, 4), -0.04773740823058334); VERIFY_IS_APPROX(ab.coeff(2, 5), 0.2300535609645254); VERIFY_IS_APPROX(ab.coeff(3, 0), -0.8172945853260133); VERIFY_IS_APPROX(ab.coeff(3, 1), 0.2150086428359221); VERIFY_IS_APPROX(ab.coeff(3, 2), 0.5825113847292743); VERIFY_IS_APPROX(ab.coeff(3, 3), -0.1532433770097174); VERIFY_IS_APPROX(ab.coeff(3, 4), -0.329383387282399); VERIFY_IS_APPROX(ab.coeff(3, 5), 0.08665207912033064); VERIFY_IS_APPROX(ab.coeff(4, 0), 0.8451267514863225); VERIFY_IS_APPROX(ab.coeff(4, 1), -0.481996458918977); VERIFY_IS_APPROX(ab.coeff(4, 2), -0.6023482390791535); VERIFY_IS_APPROX(ab.coeff(4, 3), 0.3435339347164565); VERIFY_IS_APPROX(ab.coeff(4, 4), 0.3406002157428891); VERIFY_IS_APPROX(ab.coeff(4, 5), -0.1942526344200915); VERIFY_IS_APPROX(ab.coeff(5, 0), 0.1111982482925399); VERIFY_IS_APPROX(ab.coeff(5, 1), -0.5358806424754169); VERIFY_IS_APPROX(ab.coeff(5, 2), -0.07925446559335647); VERIFY_IS_APPROX(ab.coeff(5, 3), 0.3819388757769038); VERIFY_IS_APPROX(ab.coeff(5, 4), 0.04481475387219876); VERIFY_IS_APPROX(ab.coeff(5, 5), -0.2159688616158057); } template void check_sparse_kronecker_product(const MatrixType& ab) { VERIFY_IS_EQUAL(ab.rows(), 12); VERIFY_IS_EQUAL(ab.cols(), 10); VERIFY_IS_EQUAL(ab.nonZeros(), 3 * 2); VERIFY_IS_APPROX(ab.coeff(3, 0), -0.04); VERIFY_IS_APPROX(ab.coeff(5, 1), 0.05); VERIFY_IS_APPROX(ab.coeff(0, 6), -0.08); VERIFY_IS_APPROX(ab.coeff(2, 7), 0.10); VERIFY_IS_APPROX(ab.coeff(6, 8), 0.12); VERIFY_IS_APPROX(ab.coeff(8, 9), -0.15); } EIGEN_DECLARE_TEST(kronecker_product) { // DM = dense matrix; SM = sparse matrix Matrix DM_a; SparseMatrix SM_a(2, 3); SM_a.insert(0, 0) = DM_a.coeffRef(0, 0) = -0.4461540300782201; SM_a.insert(0, 1) = DM_a.coeffRef(0, 1) = -0.8057364375283049; SM_a.insert(0, 2) = DM_a.coeffRef(0, 2) = 0.3896572459516341; SM_a.insert(1, 0) = DM_a.coeffRef(1, 0) = -0.9076572187376921; SM_a.insert(1, 1) = DM_a.coeffRef(1, 1) = 0.6469156566545853; SM_a.insert(1, 2) = DM_a.coeffRef(1, 2) = -0.3658010398782789; MatrixXd DM_b(3, 2); SparseMatrix SM_b(3, 2); SM_b.insert(0, 0) = DM_b.coeffRef(0, 0) = 0.9004440976767099; SM_b.insert(0, 1) = DM_b.coeffRef(0, 1) = -0.2368830858139832; SM_b.insert(1, 0) = DM_b.coeffRef(1, 0) = -0.9311078389941825; SM_b.insert(1, 1) = DM_b.coeffRef(1, 1) = 0.5310335762980047; SM_b.insert(2, 0) = DM_b.coeffRef(2, 0) = -0.1225112806872035; SM_b.insert(2, 1) = DM_b.coeffRef(2, 1) = 0.5903998022741264; SparseMatrix SM_row_a(SM_a), SM_row_b(SM_b); // test DM_fixedSize = kroneckerProduct(DM_block,DM) Matrix DM_fix_ab = kroneckerProduct(DM_a.topLeftCorner<2, 3>(), DM_b); CALL_SUBTEST(check_kronecker_product(DM_fix_ab)); CALL_SUBTEST(check_kronecker_product(kroneckerProduct(DM_a.topLeftCorner<2, 3>(), DM_b))); for (int i = 0; i < DM_fix_ab.rows(); ++i) for (int j = 0; j < DM_fix_ab.cols(); ++j) VERIFY_IS_APPROX(kroneckerProduct(DM_a, DM_b).coeff(i, j), DM_fix_ab(i, j)); // test DM_block = kroneckerProduct(DM,DM) MatrixXd DM_block_ab(10, 15); DM_block_ab.block<6, 6>(2, 5) = kroneckerProduct(DM_a, DM_b); CALL_SUBTEST(check_kronecker_product(DM_block_ab.block<6, 6>(2, 5))); // test DM = kroneckerProduct(DM,DM) MatrixXd DM_ab = kroneckerProduct(DM_a, DM_b); CALL_SUBTEST(check_kronecker_product(DM_ab)); CALL_SUBTEST(check_kronecker_product(kroneckerProduct(DM_a, DM_b))); // test SM = kroneckerProduct(SM,DM) SparseMatrix SM_ab = kroneckerProduct(SM_a, DM_b); CALL_SUBTEST(check_kronecker_product(SM_ab)); SparseMatrix SM_ab2 = kroneckerProduct(SM_a, DM_b); CALL_SUBTEST(check_kronecker_product(SM_ab2)); CALL_SUBTEST(check_kronecker_product(kroneckerProduct(SM_a, DM_b))); // test SM = kroneckerProduct(DM,SM) SM_ab.setZero(); SM_ab.insert(0, 0) = 37.0; SM_ab = kroneckerProduct(DM_a, SM_b); CALL_SUBTEST(check_kronecker_product(SM_ab)); SM_ab2.setZero(); SM_ab2.insert(0, 0) = 37.0; SM_ab2 = kroneckerProduct(DM_a, SM_b); CALL_SUBTEST(check_kronecker_product(SM_ab2)); CALL_SUBTEST(check_kronecker_product(kroneckerProduct(DM_a, SM_b))); // test SM = kroneckerProduct(SM,SM) SM_ab.resize(2, 33); SM_ab.insert(0, 0) = 37.0; SM_ab = kroneckerProduct(SM_a, SM_b); CALL_SUBTEST(check_kronecker_product(SM_ab)); SM_ab2.resize(5, 11); SM_ab2.insert(0, 0) = 37.0; SM_ab2 = kroneckerProduct(SM_a, SM_b); CALL_SUBTEST(check_kronecker_product(SM_ab2)); CALL_SUBTEST(check_kronecker_product(kroneckerProduct(SM_a, SM_b))); // test SM = kroneckerProduct(SM,SM) with sparse pattern SM_a.resize(4, 5); SM_b.resize(3, 2); SM_a.resizeNonZeros(0); SM_b.resizeNonZeros(0); SM_a.insert(1, 0) = -0.1; SM_a.insert(0, 3) = -0.2; SM_a.insert(2, 4) = 0.3; SM_a.finalize(); SM_b.insert(0, 0) = 0.4; SM_b.insert(2, 1) = -0.5; SM_b.finalize(); SM_ab.resize(1, 1); SM_ab.insert(0, 0) = 37.0; SM_ab = kroneckerProduct(SM_a, SM_b); CALL_SUBTEST(check_sparse_kronecker_product(SM_ab)); // test dimension of result of DM = kroneckerProduct(DM,DM) MatrixXd DM_a2(2, 1); MatrixXd DM_b2(5, 4); MatrixXd DM_ab2 = kroneckerProduct(DM_a2, DM_b2); CALL_SUBTEST(check_dimension(DM_ab2, 2 * 5, 1 * 4)); DM_a2.resize(10, 9); DM_b2.resize(4, 8); DM_ab2 = kroneckerProduct(DM_a2, DM_b2); CALL_SUBTEST(check_dimension(DM_ab2, 10 * 4, 9 * 8)); for (int i = 0; i < g_repeat; i++) { double density = Eigen::internal::random(0.01, 0.5); int ra = Eigen::internal::random(1, 50); int ca = Eigen::internal::random(1, 50); int rb = Eigen::internal::random(1, 50); int cb = Eigen::internal::random(1, 50); SparseMatrix sA(ra, ca), sB(rb, cb), sC; SparseMatrix sC2; MatrixXf dA(ra, ca), dB(rb, cb), dC; initSparse(density, dA, sA); initSparse(density, dB, sB); sC = kroneckerProduct(sA, sB); dC = kroneckerProduct(dA, dB); VERIFY_IS_APPROX(MatrixXf(sC), dC); sC = kroneckerProduct(sA.transpose(), sB); dC = kroneckerProduct(dA.transpose(), dB); VERIFY_IS_APPROX(MatrixXf(sC), dC); sC = kroneckerProduct(sA.transpose(), sB.transpose()); dC = kroneckerProduct(dA.transpose(), dB.transpose()); VERIFY_IS_APPROX(MatrixXf(sC), dC); sC = kroneckerProduct(sA, sB.transpose()); dC = kroneckerProduct(dA, dB.transpose()); VERIFY_IS_APPROX(MatrixXf(sC), dC); sC2 = kroneckerProduct(sA, sB); dC = kroneckerProduct(dA, dB); VERIFY_IS_APPROX(MatrixXf(sC2), dC); sC2 = kroneckerProduct(dA, sB); dC = kroneckerProduct(dA, dB); VERIFY_IS_APPROX(MatrixXf(sC2), dC); sC2 = kroneckerProduct(sA, dB); dC = kroneckerProduct(dA, dB); VERIFY_IS_APPROX(MatrixXf(sC2), dC); sC2 = kroneckerProduct(2 * sA, sB); dC = kroneckerProduct(2 * dA, dB); VERIFY_IS_APPROX(MatrixXf(sC2), dC); } } #endif #ifdef EIGEN_TEST_PART_2 // simply check that for a dense kronecker product, sparse module is not needed #include "main.h" #include EIGEN_DECLARE_TEST(kronecker_product) { MatrixXd a(2, 2), b(3, 3), c; a.setRandom(); b.setRandom(); c = kroneckerProduct(a, b); VERIFY_IS_APPROX(c.block(3, 3, 3, 3), a(1, 1) * b); } #endif