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 #ifdef EIGEN_TEST_PART_1
13 
14 #include "sparse.h"
15 #include <Eigen/SparseExtra>
16 #include <Eigen/KroneckerProduct>
17 
18 template<typename MatrixType>
check_dimension(const MatrixType & ab,const int rows,const int cols)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>
check_kronecker_product(const MatrixType & ab)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>
check_sparse_kronecker_product(const MatrixType & ab)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 
test_kronecker_product()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 
test_kronecker_product()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