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