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 // Copyright (C) 2012 Chen-Pang He <jdh8 (at) 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 #ifdef EIGEN_TEST_PART_1 13 14 #include "sparse.h" 15 #include <Eigen/SparseExtra> 16 #include <Eigen/KroneckerProduct> 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 DM_fixedSize = kroneckerProduct(DM_block,DM) 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 CALL_SUBTEST(check_kronecker_product(kroneckerProduct(DM_a.topLeftCorner<2,3>(),DM_b))); 115 116 for(int i=0;i<DM_fix_ab.rows();++i) 117 for(int j=0;j<DM_fix_ab.cols();++j) 118 VERIFY_IS_APPROX(kroneckerProduct(DM_a,DM_b).coeff(i,j), DM_fix_ab(i,j)); 119 120 // test DM_block = kroneckerProduct(DM,DM) 121 MatrixXd DM_block_ab(10,15); 122 DM_block_ab.block<6,6>(2,5) = kroneckerProduct(DM_a,DM_b); 123 CALL_SUBTEST(check_kronecker_product(DM_block_ab.block<6,6>(2,5))); 124 125 // test DM = kroneckerProduct(DM,DM) 126 MatrixXd DM_ab = kroneckerProduct(DM_a,DM_b); 127 CALL_SUBTEST(check_kronecker_product(DM_ab)); 128 CALL_SUBTEST(check_kronecker_product(kroneckerProduct(DM_a,DM_b))); 129 130 // test SM = kroneckerProduct(SM,DM) 131 SparseMatrix<double> SM_ab = kroneckerProduct(SM_a,DM_b); 132 CALL_SUBTEST(check_kronecker_product(SM_ab)); 133 SparseMatrix<double,RowMajor> SM_ab2 = kroneckerProduct(SM_a,DM_b); 134 CALL_SUBTEST(check_kronecker_product(SM_ab2)); 135 CALL_SUBTEST(check_kronecker_product(kroneckerProduct(SM_a,DM_b))); 136 137 // test SM = kroneckerProduct(DM,SM) 138 SM_ab.setZero(); 139 SM_ab.insert(0,0)=37.0; 140 SM_ab = kroneckerProduct(DM_a,SM_b); 141 CALL_SUBTEST(check_kronecker_product(SM_ab)); 142 SM_ab2.setZero(); 143 SM_ab2.insert(0,0)=37.0; 144 SM_ab2 = kroneckerProduct(DM_a,SM_b); 145 CALL_SUBTEST(check_kronecker_product(SM_ab2)); 146 CALL_SUBTEST(check_kronecker_product(kroneckerProduct(DM_a,SM_b))); 147 148 // test SM = kroneckerProduct(SM,SM) 149 SM_ab.resize(2,33); 150 SM_ab.insert(0,0)=37.0; 151 SM_ab = kroneckerProduct(SM_a,SM_b); 152 CALL_SUBTEST(check_kronecker_product(SM_ab)); 153 SM_ab2.resize(5,11); 154 SM_ab2.insert(0,0)=37.0; 155 SM_ab2 = kroneckerProduct(SM_a,SM_b); 156 CALL_SUBTEST(check_kronecker_product(SM_ab2)); 157 CALL_SUBTEST(check_kronecker_product(kroneckerProduct(SM_a,SM_b))); 158 159 // test SM = kroneckerProduct(SM,SM) with sparse pattern 160 SM_a.resize(4,5); 161 SM_b.resize(3,2); 162 SM_a.resizeNonZeros(0); 163 SM_b.resizeNonZeros(0); 164 SM_a.insert(1,0) = -0.1; 165 SM_a.insert(0,3) = -0.2; 166 SM_a.insert(2,4) = 0.3; 167 SM_a.finalize(); 168 169 SM_b.insert(0,0) = 0.4; 170 SM_b.insert(2,1) = -0.5; 171 SM_b.finalize(); 172 SM_ab.resize(1,1); 173 SM_ab.insert(0,0)=37.0; 174 SM_ab = kroneckerProduct(SM_a,SM_b); 175 CALL_SUBTEST(check_sparse_kronecker_product(SM_ab)); 176 177 // test dimension of result of DM = kroneckerProduct(DM,DM) 178 MatrixXd DM_a2(2,1); 179 MatrixXd DM_b2(5,4); 180 MatrixXd DM_ab2 = kroneckerProduct(DM_a2,DM_b2); 181 CALL_SUBTEST(check_dimension(DM_ab2,2*5,1*4)); 182 DM_a2.resize(10,9); 183 DM_b2.resize(4,8); 184 DM_ab2 = kroneckerProduct(DM_a2,DM_b2); 185 CALL_SUBTEST(check_dimension(DM_ab2,10*4,9*8)); 186 187 for(int i = 0; i < g_repeat; i++) 188 { 189 double density = Eigen::internal::random<double>(0.01,0.5); 190 int ra = Eigen::internal::random<int>(1,50); 191 int ca = Eigen::internal::random<int>(1,50); 192 int rb = Eigen::internal::random<int>(1,50); 193 int cb = Eigen::internal::random<int>(1,50); 194 SparseMatrix<float,ColMajor> sA(ra,ca), sB(rb,cb), sC; 195 SparseMatrix<float,RowMajor> sC2; 196 MatrixXf dA(ra,ca), dB(rb,cb), dC; 197 initSparse(density, dA, sA); 198 initSparse(density, dB, sB); 199 200 sC = kroneckerProduct(sA,sB); 201 dC = kroneckerProduct(dA,dB); 202 VERIFY_IS_APPROX(MatrixXf(sC),dC); 203 204 sC = kroneckerProduct(sA.transpose(),sB); 205 dC = kroneckerProduct(dA.transpose(),dB); 206 VERIFY_IS_APPROX(MatrixXf(sC),dC); 207 208 sC = kroneckerProduct(sA.transpose(),sB.transpose()); 209 dC = kroneckerProduct(dA.transpose(),dB.transpose()); 210 VERIFY_IS_APPROX(MatrixXf(sC),dC); 211 212 sC = kroneckerProduct(sA,sB.transpose()); 213 dC = kroneckerProduct(dA,dB.transpose()); 214 VERIFY_IS_APPROX(MatrixXf(sC),dC); 215 216 sC2 = kroneckerProduct(sA,sB); 217 dC = kroneckerProduct(dA,dB); 218 VERIFY_IS_APPROX(MatrixXf(sC2),dC); 219 220 sC2 = kroneckerProduct(dA,sB); 221 dC = kroneckerProduct(dA,dB); 222 VERIFY_IS_APPROX(MatrixXf(sC2),dC); 223 224 sC2 = kroneckerProduct(sA,dB); 225 dC = kroneckerProduct(dA,dB); 226 VERIFY_IS_APPROX(MatrixXf(sC2),dC); 227 228 sC2 = kroneckerProduct(2*sA,sB); 229 dC = kroneckerProduct(2*dA,dB); 230 VERIFY_IS_APPROX(MatrixXf(sC2),dC); 231 } 232 } 233 234 #endif 235 236 #ifdef EIGEN_TEST_PART_2 237 238 // simply check that for a dense kronecker product, sparse module is not needed 239 240 #include "main.h" 241 #include <Eigen/KroneckerProduct> 242 243 void test_kronecker_product() 244 { 245 MatrixXd a(2,2), b(3,3), c; 246 a.setRandom(); 247 b.setRandom(); 248 c = kroneckerProduct(a,b); 249 VERIFY_IS_APPROX(c.block(3,3,3,3), a(1,1)*b); 250 } 251 252 #endif 253
