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 //
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>
check_dimension(const MatrixType & ab,const unsigned int rows,const unsigned int cols)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>
check_kronecker_product(const MatrixType & ab)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>
check_sparse_kronecker_product(const MatrixType & ab)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
test_kronecker_product()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