[136] | 1 | // This file is part of Eigen, a lightweight C++ template library
|
---|
| 2 | // for linear algebra.
|
---|
| 3 | //
|
---|
| 4 | // Copyright (C) 2011 Kolja Brix <brix@igpm.rwth-aachen.de>
|
---|
| 5 | // Copyright (C) 2011 Andreas Platen <andiplaten@gmx.de>
|
---|
| 6 | // Copyright (C) 2012 Chen-Pang He <jdh8@ms63.hinet.net>
|
---|
| 7 | //
|
---|
| 8 | // This Source Code Form is subject to the terms of the Mozilla
|
---|
| 9 | // Public License v. 2.0. If a copy of the MPL was not distributed
|
---|
| 10 | // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
|
---|
| 11 |
|
---|
| 12 |
|
---|
| 13 | #include "sparse.h"
|
---|
| 14 | #include <Eigen/SparseExtra>
|
---|
| 15 | #include <Eigen/KroneckerProduct>
|
---|
| 16 |
|
---|
| 17 |
|
---|
| 18 | template<typename MatrixType>
|
---|
| 19 | void check_dimension(const MatrixType& ab, const int rows, const int cols)
|
---|
| 20 | {
|
---|
| 21 | VERIFY_IS_EQUAL(ab.rows(), rows);
|
---|
| 22 | VERIFY_IS_EQUAL(ab.cols(), cols);
|
---|
| 23 | }
|
---|
| 24 |
|
---|
| 25 |
|
---|
| 26 | template<typename MatrixType>
|
---|
| 27 | void check_kronecker_product(const MatrixType& ab)
|
---|
| 28 | {
|
---|
| 29 | VERIFY_IS_EQUAL(ab.rows(), 6);
|
---|
| 30 | VERIFY_IS_EQUAL(ab.cols(), 6);
|
---|
| 31 | VERIFY_IS_EQUAL(ab.nonZeros(), 36);
|
---|
| 32 | VERIFY_IS_APPROX(ab.coeff(0,0), -0.4017367630386106);
|
---|
| 33 | VERIFY_IS_APPROX(ab.coeff(0,1), 0.1056863433932735);
|
---|
| 34 | VERIFY_IS_APPROX(ab.coeff(0,2), -0.7255206194554212);
|
---|
| 35 | VERIFY_IS_APPROX(ab.coeff(0,3), 0.1908653336744706);
|
---|
| 36 | VERIFY_IS_APPROX(ab.coeff(0,4), 0.350864567234111);
|
---|
| 37 | VERIFY_IS_APPROX(ab.coeff(0,5), -0.0923032108308013);
|
---|
| 38 | VERIFY_IS_APPROX(ab.coeff(1,0), 0.415417514804677);
|
---|
| 39 | VERIFY_IS_APPROX(ab.coeff(1,1), -0.2369227701722048);
|
---|
| 40 | VERIFY_IS_APPROX(ab.coeff(1,2), 0.7502275131458511);
|
---|
| 41 | VERIFY_IS_APPROX(ab.coeff(1,3), -0.4278731019742696);
|
---|
| 42 | VERIFY_IS_APPROX(ab.coeff(1,4), -0.3628129162264507);
|
---|
| 43 | VERIFY_IS_APPROX(ab.coeff(1,5), 0.2069210808481275);
|
---|
| 44 | VERIFY_IS_APPROX(ab.coeff(2,0), 0.05465890160863986);
|
---|
| 45 | VERIFY_IS_APPROX(ab.coeff(2,1), -0.2634092511419858);
|
---|
| 46 | VERIFY_IS_APPROX(ab.coeff(2,2), 0.09871180285793758);
|
---|
| 47 | VERIFY_IS_APPROX(ab.coeff(2,3), -0.4757066334017702);
|
---|
| 48 | VERIFY_IS_APPROX(ab.coeff(2,4), -0.04773740823058334);
|
---|
| 49 | VERIFY_IS_APPROX(ab.coeff(2,5), 0.2300535609645254);
|
---|
| 50 | VERIFY_IS_APPROX(ab.coeff(3,0), -0.8172945853260133);
|
---|
| 51 | VERIFY_IS_APPROX(ab.coeff(3,1), 0.2150086428359221);
|
---|
| 52 | VERIFY_IS_APPROX(ab.coeff(3,2), 0.5825113847292743);
|
---|
| 53 | VERIFY_IS_APPROX(ab.coeff(3,3), -0.1532433770097174);
|
---|
| 54 | VERIFY_IS_APPROX(ab.coeff(3,4), -0.329383387282399);
|
---|
| 55 | VERIFY_IS_APPROX(ab.coeff(3,5), 0.08665207912033064);
|
---|
| 56 | VERIFY_IS_APPROX(ab.coeff(4,0), 0.8451267514863225);
|
---|
| 57 | VERIFY_IS_APPROX(ab.coeff(4,1), -0.481996458918977);
|
---|
| 58 | VERIFY_IS_APPROX(ab.coeff(4,2), -0.6023482390791535);
|
---|
| 59 | VERIFY_IS_APPROX(ab.coeff(4,3), 0.3435339347164565);
|
---|
| 60 | VERIFY_IS_APPROX(ab.coeff(4,4), 0.3406002157428891);
|
---|
| 61 | VERIFY_IS_APPROX(ab.coeff(4,5), -0.1942526344200915);
|
---|
| 62 | VERIFY_IS_APPROX(ab.coeff(5,0), 0.1111982482925399);
|
---|
| 63 | VERIFY_IS_APPROX(ab.coeff(5,1), -0.5358806424754169);
|
---|
| 64 | VERIFY_IS_APPROX(ab.coeff(5,2), -0.07925446559335647);
|
---|
| 65 | VERIFY_IS_APPROX(ab.coeff(5,3), 0.3819388757769038);
|
---|
| 66 | VERIFY_IS_APPROX(ab.coeff(5,4), 0.04481475387219876);
|
---|
| 67 | VERIFY_IS_APPROX(ab.coeff(5,5), -0.2159688616158057);
|
---|
| 68 | }
|
---|
| 69 |
|
---|
| 70 |
|
---|
| 71 | template<typename MatrixType>
|
---|
| 72 | void check_sparse_kronecker_product(const MatrixType& ab)
|
---|
| 73 | {
|
---|
| 74 | VERIFY_IS_EQUAL(ab.rows(), 12);
|
---|
| 75 | VERIFY_IS_EQUAL(ab.cols(), 10);
|
---|
| 76 | VERIFY_IS_EQUAL(ab.nonZeros(), 3*2);
|
---|
| 77 | VERIFY_IS_APPROX(ab.coeff(3,0), -0.04);
|
---|
| 78 | VERIFY_IS_APPROX(ab.coeff(5,1), 0.05);
|
---|
| 79 | VERIFY_IS_APPROX(ab.coeff(0,6), -0.08);
|
---|
| 80 | VERIFY_IS_APPROX(ab.coeff(2,7), 0.10);
|
---|
| 81 | VERIFY_IS_APPROX(ab.coeff(6,8), 0.12);
|
---|
| 82 | VERIFY_IS_APPROX(ab.coeff(8,9), -0.15);
|
---|
| 83 | }
|
---|
| 84 |
|
---|
| 85 |
|
---|
| 86 | void test_kronecker_product()
|
---|
| 87 | {
|
---|
| 88 | // DM = dense matrix; SM = sparse matrix
|
---|
| 89 |
|
---|
| 90 | Matrix<double, 2, 3> DM_a;
|
---|
| 91 | SparseMatrix<double> SM_a(2,3);
|
---|
| 92 | SM_a.insert(0,0) = DM_a.coeffRef(0,0) = -0.4461540300782201;
|
---|
| 93 | SM_a.insert(0,1) = DM_a.coeffRef(0,1) = -0.8057364375283049;
|
---|
| 94 | SM_a.insert(0,2) = DM_a.coeffRef(0,2) = 0.3896572459516341;
|
---|
| 95 | SM_a.insert(1,0) = DM_a.coeffRef(1,0) = -0.9076572187376921;
|
---|
| 96 | SM_a.insert(1,1) = DM_a.coeffRef(1,1) = 0.6469156566545853;
|
---|
| 97 | SM_a.insert(1,2) = DM_a.coeffRef(1,2) = -0.3658010398782789;
|
---|
| 98 |
|
---|
| 99 | MatrixXd DM_b(3,2);
|
---|
| 100 | SparseMatrix<double> SM_b(3,2);
|
---|
| 101 | SM_b.insert(0,0) = DM_b.coeffRef(0,0) = 0.9004440976767099;
|
---|
| 102 | SM_b.insert(0,1) = DM_b.coeffRef(0,1) = -0.2368830858139832;
|
---|
| 103 | SM_b.insert(1,0) = DM_b.coeffRef(1,0) = -0.9311078389941825;
|
---|
| 104 | SM_b.insert(1,1) = DM_b.coeffRef(1,1) = 0.5310335762980047;
|
---|
| 105 | SM_b.insert(2,0) = DM_b.coeffRef(2,0) = -0.1225112806872035;
|
---|
| 106 | SM_b.insert(2,1) = DM_b.coeffRef(2,1) = 0.5903998022741264;
|
---|
| 107 |
|
---|
| 108 | SparseMatrix<double,RowMajor> SM_row_a(SM_a), SM_row_b(SM_b);
|
---|
| 109 |
|
---|
| 110 | // test kroneckerProduct(DM_block,DM,DM_fixedSize)
|
---|
| 111 | Matrix<double, 6, 6> DM_fix_ab = kroneckerProduct(DM_a.topLeftCorner<2,3>(),DM_b);
|
---|
| 112 |
|
---|
| 113 | CALL_SUBTEST(check_kronecker_product(DM_fix_ab));
|
---|
| 114 |
|
---|
| 115 | for(int i=0;i<DM_fix_ab.rows();++i)
|
---|
| 116 | for(int j=0;j<DM_fix_ab.cols();++j)
|
---|
| 117 | VERIFY_IS_APPROX(kroneckerProduct(DM_a,DM_b).coeff(i,j), DM_fix_ab(i,j));
|
---|
| 118 |
|
---|
| 119 | // test kroneckerProduct(DM,DM,DM_block)
|
---|
| 120 | MatrixXd DM_block_ab(10,15);
|
---|
| 121 | DM_block_ab.block<6,6>(2,5) = kroneckerProduct(DM_a,DM_b);
|
---|
| 122 | CALL_SUBTEST(check_kronecker_product(DM_block_ab.block<6,6>(2,5)));
|
---|
| 123 |
|
---|
| 124 | // test kroneckerProduct(DM,DM,DM)
|
---|
| 125 | MatrixXd DM_ab = kroneckerProduct(DM_a,DM_b);
|
---|
| 126 | CALL_SUBTEST(check_kronecker_product(DM_ab));
|
---|
| 127 |
|
---|
| 128 | // test kroneckerProduct(SM,DM,SM)
|
---|
| 129 | SparseMatrix<double> SM_ab = kroneckerProduct(SM_a,DM_b);
|
---|
| 130 | CALL_SUBTEST(check_kronecker_product(SM_ab));
|
---|
| 131 | SparseMatrix<double,RowMajor> SM_ab2 = kroneckerProduct(SM_a,DM_b);
|
---|
| 132 | CALL_SUBTEST(check_kronecker_product(SM_ab2));
|
---|
| 133 |
|
---|
| 134 | // test kroneckerProduct(DM,SM,SM)
|
---|
| 135 | SM_ab.setZero();
|
---|
| 136 | SM_ab.insert(0,0)=37.0;
|
---|
| 137 | SM_ab = kroneckerProduct(DM_a,SM_b);
|
---|
| 138 | CALL_SUBTEST(check_kronecker_product(SM_ab));
|
---|
| 139 | SM_ab2.setZero();
|
---|
| 140 | SM_ab2.insert(0,0)=37.0;
|
---|
| 141 | SM_ab2 = kroneckerProduct(DM_a,SM_b);
|
---|
| 142 | CALL_SUBTEST(check_kronecker_product(SM_ab2));
|
---|
| 143 |
|
---|
| 144 | // test kroneckerProduct(SM,SM,SM)
|
---|
| 145 | SM_ab.resize(2,33);
|
---|
| 146 | SM_ab.insert(0,0)=37.0;
|
---|
| 147 | SM_ab = kroneckerProduct(SM_a,SM_b);
|
---|
| 148 | CALL_SUBTEST(check_kronecker_product(SM_ab));
|
---|
| 149 | SM_ab2.resize(5,11);
|
---|
| 150 | SM_ab2.insert(0,0)=37.0;
|
---|
| 151 | SM_ab2 = kroneckerProduct(SM_a,SM_b);
|
---|
| 152 | CALL_SUBTEST(check_kronecker_product(SM_ab2));
|
---|
| 153 |
|
---|
| 154 | // test kroneckerProduct(SM,SM,SM) with sparse pattern
|
---|
| 155 | SM_a.resize(4,5);
|
---|
| 156 | SM_b.resize(3,2);
|
---|
| 157 | SM_a.resizeNonZeros(0);
|
---|
| 158 | SM_b.resizeNonZeros(0);
|
---|
| 159 | SM_a.insert(1,0) = -0.1;
|
---|
| 160 | SM_a.insert(0,3) = -0.2;
|
---|
| 161 | SM_a.insert(2,4) = 0.3;
|
---|
| 162 | SM_a.finalize();
|
---|
| 163 |
|
---|
| 164 | SM_b.insert(0,0) = 0.4;
|
---|
| 165 | SM_b.insert(2,1) = -0.5;
|
---|
| 166 | SM_b.finalize();
|
---|
| 167 | SM_ab.resize(1,1);
|
---|
| 168 | SM_ab.insert(0,0)=37.0;
|
---|
| 169 | SM_ab = kroneckerProduct(SM_a,SM_b);
|
---|
| 170 | CALL_SUBTEST(check_sparse_kronecker_product(SM_ab));
|
---|
| 171 |
|
---|
| 172 | // test dimension of result of kroneckerProduct(DM,DM,DM)
|
---|
| 173 | MatrixXd DM_a2(2,1);
|
---|
| 174 | MatrixXd DM_b2(5,4);
|
---|
| 175 | MatrixXd DM_ab2 = kroneckerProduct(DM_a2,DM_b2);
|
---|
| 176 | CALL_SUBTEST(check_dimension(DM_ab2,2*5,1*4));
|
---|
| 177 | DM_a2.resize(10,9);
|
---|
| 178 | DM_b2.resize(4,8);
|
---|
| 179 | DM_ab2 = kroneckerProduct(DM_a2,DM_b2);
|
---|
| 180 | CALL_SUBTEST(check_dimension(DM_ab2,10*4,9*8));
|
---|
| 181 | }
|
---|