1 // This is core/vnl/tests/test_fastops.cxx
2 #include "testlib/testlib_test.h"
3 //:
4 // \file
5 // \author Peter Vanroose
6 // \date 17 June 2004
7 #include "vnl/vnl_fastops.h"
8 
9 void
test_fastops()10 test_fastops()
11 {
12   // The data to work with
13   vnl_matrix<double> result_m, X;
14   vnl_vector<double> result_v, Y;
15   vnl_matrix<double> id1x1(1, 1);
16   id1x1.set_identity();
17   vnl_matrix<double> id2x2(2, 2);
18   id2x2.set_identity();
19   vnl_matrix<double> id3x3(3, 3);
20   id3x3.set_identity();
21   vnl_matrix<double> id10x10(10, 10);
22   id10x10.set_identity();
23   vnl_matrix<double> id99x99(99, 99);
24   id99x99.set_identity();
25   double data[10] = { -7.5, 1.5, 0., -107.25, 511., -509.75, 1.25, -1., 0., 1. };
26   vnl_matrix<double> m1x1(data, 1, 1); // the number -7.5
27   vnl_matrix<double> m2x2(data, 2, 2);
28   vnl_matrix<double> m2x2t = m2x2.transpose();
29   vnl_matrix<double> m3x3(data, 3, 3);
30   vnl_matrix<double> m3x3t = m3x3.transpose();
31   vnl_vector<double> v1(data, 1); // the number -7.5
32   vnl_vector<double> v2(data, 2);
33   vnl_vector<double> v3(data, 3);
34   vnl_vector<double> v10(data, 10);
35   vnl_matrix<double> m10x2(10, 2), m2x10(2, 10);
36   for (unsigned int i = 0; i < 10; ++i)
37     for (unsigned int j = 0; j < 2; ++j)
38       m10x2[i][j] = m2x10[j][i] = (i + 1) * 1.5 + (j + 1) * (j + i);
39 
40   // First test: $I \times I$
41   result_m.set_size(1, 1);
42   vnl_fastops::AtA(result_m, id1x1);
43   TEST("vnl_fastops::AtA(id1x1)", result_m, id1x1);
44   vnl_fastops::AB(result_m, id1x1, id1x1);
45   TEST("vnl_fastops::AB(id1x1,id1x1)", result_m, id1x1);
46   vnl_fastops::AtB(result_m, id1x1, id1x1);
47   TEST("vnl_fastops::AtB(id1x1,id1x1)", result_m, id1x1);
48   vnl_fastops::ABt(result_m, id1x1, id1x1);
49   TEST("vnl_fastops::ABt(id1x1,id1x1)", result_m, id1x1);
50   result_m.set_size(2, 2);
51   vnl_fastops::AtA(result_m, id2x2);
52   TEST("vnl_fastops::AtA(id2x2)", result_m, id2x2);
53   vnl_fastops::AB(result_m, id2x2, id2x2);
54   TEST("vnl_fastops::AB(id2x2,id2x2)", result_m, id2x2);
55   vnl_fastops::AtB(result_m, id2x2, id2x2);
56   TEST("vnl_fastops::AtB(id2x2,id2x2)", result_m, id2x2);
57   vnl_fastops::ABt(result_m, id2x2, id2x2);
58   TEST("vnl_fastops::ABt(id2x2,id2x2)", result_m, id2x2);
59   result_m.set_size(99, 99);
60   vnl_fastops::AtA(result_m, id99x99);
61   TEST("vnl_fastops::AtA(id99x99)", result_m, id99x99);
62   vnl_fastops::AB(result_m, id99x99, id99x99);
63   TEST("vnl_fastops::AB(id99x99,id99x99)", result_m, id99x99);
64   vnl_fastops::AtB(result_m, id99x99, id99x99);
65   TEST("vnl_fastops::AtB(id99x99,id99x99)", result_m, id99x99);
66   vnl_fastops::ABt(result_m, id99x99, id99x99);
67   TEST("vnl_fastops::ABt(id99x99,id99x99)", result_m, id99x99);
68 
69   // Second test: $I \times M$ and $M \times I$ and $M^\top \times I$
70   result_m.set_size(1, 1);
71   result_v.set_size(1);
72   vnl_fastops::AB(result_m, id1x1, m1x1);
73   TEST("vnl_fastops::AB(id1x1,m1x1)", result_m, m1x1);
74   vnl_fastops::AB(result_m, m1x1, id1x1);
75   TEST("vnl_fastops::AB(m1x1,id1x1)", result_m, m1x1);
76   vnl_fastops::AtB(result_v, id1x1, v1);
77   TEST("vnl_fastops::AtB(id1x1,v1)", result_v, v1);
78   vnl_fastops::AtB(result_m, id1x1, m1x1);
79   TEST("vnl_fastops::AtB(id1x1,m1x1)", result_m, m1x1);
80   vnl_fastops::ABt(result_m, m1x1, id1x1);
81   TEST("vnl_fastops::ABt(m1x1,id1x1)", result_m, m1x1);
82   vnl_fastops::AtB(result_m, m1x1, id1x1);
83   TEST("vnl_fastops::AtB(m1x1,id1x1)", result_m, m1x1);
84   vnl_fastops::ABt(result_m, id1x1, m1x1);
85   TEST("vnl_fastops::ABt(id1x1,m1x1)", result_m, m1x1);
86   result_m.set_size(2, 2);
87   result_v.set_size(2);
88   vnl_fastops::AB(result_m, id2x2, m2x2);
89   TEST("vnl_fastops::AB(id2x2,m2x2)", result_m, m2x2);
90   vnl_fastops::AB(result_m, m2x2, id2x2);
91   TEST("vnl_fastops::AB(m2x2,id2x2)", result_m, m2x2);
92   vnl_fastops::AtB(result_v, id2x2, v2);
93   TEST("vnl_fastops::AtB(id2x2,v2)", result_v, v2);
94   vnl_fastops::AtB(result_m, id2x2, m2x2);
95   TEST("vnl_fastops::AtB(id2x2,m2x2)", result_m, m2x2);
96   vnl_fastops::ABt(result_m, m2x2, id2x2);
97   TEST("vnl_fastops::ABt(m2x2,id2x2)", result_m, m2x2);
98   vnl_fastops::AtB(result_m, m2x2, id2x2);
99   TEST("vnl_fastops::AtB(m2x2,id2x2)", result_m, m2x2t);
100   vnl_fastops::ABt(result_m, id2x2, m2x2);
101   TEST("vnl_fastops::ABt(id2x2,m2x2)", result_m, m2x2t);
102   result_m.set_size(3, 3);
103   result_v.set_size(3);
104   vnl_fastops::AB(result_m, id3x3, m3x3);
105   TEST("vnl_fastops::AB(id3x3,m3x3)", result_m, m3x3);
106   vnl_fastops::AB(result_m, m3x3, id3x3);
107   TEST("vnl_fastops::AB(m3x3,id3x3)", result_m, m3x3);
108   vnl_fastops::AtB(result_v, id3x3, v3);
109   TEST("vnl_fastops::AtB(id3x3,v3)", result_v, v3);
110   vnl_fastops::AtB(result_m, id3x3, m3x3);
111   TEST("vnl_fastops::AtB(id3x3,m3x3)", result_m, m3x3);
112   vnl_fastops::ABt(result_m, m3x3, id3x3);
113   TEST("vnl_fastops::ABt(m3x3,id3x3)", result_m, m3x3);
114   vnl_fastops::AtB(result_m, m3x3, id3x3);
115   TEST("vnl_fastops::AtB(m3x3,id3x3)", result_m, m3x3t);
116   vnl_fastops::ABt(result_m, id3x3, m3x3);
117   TEST("vnl_fastops::ABt(id3x3,m3x3)", result_m, m3x3t);
118   result_v.set_size(10);
119   vnl_fastops::AtB(result_v, id10x10, v10);
120   TEST("vnl_fastops::AtB(id10x10,v10)", result_v, v10);
121 
122   // Third test: $M \times M$ and $M^\top \times M$ and $M \times M^\top$
123   result_m.set_size(1, 1);
124   result_v.set_size(1);
125   vnl_fastops::AtA(result_m, m1x1);
126   TEST("vnl_fastops::AtA(m1x1)", result_m, m1x1 * m1x1);
127   vnl_fastops::AB(result_m, m1x1, m1x1);
128   TEST("vnl_fastops::AB(m1x1,m1x1)", result_m, m1x1 * m1x1);
129   vnl_fastops::AtB(result_v, m1x1, v1);
130   TEST("vnl_fastops::AtB(m1x1,v1)", result_v, m1x1 * v1);
131   vnl_fastops::AtB(result_m, m1x1, m1x1);
132   TEST("vnl_fastops::AtB(m1x1,m1x1)", result_m, m1x1 * m1x1);
133   vnl_fastops::ABt(result_m, m1x1, m1x1);
134   TEST("vnl_fastops::ABt(m1x1,m1x1)", result_m, m1x1 * m1x1);
135   result_m.set_size(2, 2);
136   result_v.set_size(2);
137   vnl_fastops::AtA(result_m, m2x2);
138   TEST("vnl_fastops::AtA(m2x2)", result_m, m2x2t * m2x2);
139   vnl_fastops::AB(result_m, m2x2, m2x2);
140   TEST("vnl_fastops::AB(m2x2,m2x2)", result_m, m2x2 * m2x2);
141   vnl_fastops::AtB(result_v, m2x2, v2);
142   TEST("vnl_fastops::AtB(m2x2,v2)", result_v, m2x2t * v2);
143   vnl_fastops::AtB(result_m, m2x2, m2x2);
144   TEST("vnl_fastops::AtB(m2x2,m2x2)", result_m, m2x2t * m2x2);
145   vnl_fastops::ABt(result_m, m2x2, m2x2);
146   TEST("vnl_fastops::ABt(m2x2,m2x2)", result_m, m2x2 * m2x2t);
147   result_m.set_size(3, 3);
148   result_v.set_size(3);
149   vnl_fastops::AtA(result_m, m3x3);
150   TEST("vnl_fastops::AtA(m3x3)", result_m, m3x3t * m3x3);
151   vnl_fastops::AB(result_m, m3x3, m3x3);
152   TEST("vnl_fastops::AB(m3x3,m3x3)", result_m, m3x3 * m3x3);
153   vnl_fastops::AtB(result_v, m3x3, v3);
154   TEST("vnl_fastops::AtB(m3x3,v3)", result_v, m3x3t * v3);
155   vnl_fastops::AtB(result_m, m3x3, m3x3);
156   TEST("vnl_fastops::AtB(m3x3,m3x3)", result_m, m3x3t * m3x3);
157   vnl_fastops::ABt(result_m, m3x3, m3x3);
158   TEST("vnl_fastops::ABt(m3x3,m3x3)", result_m, m3x3 * m3x3t);
159   result_m.set_size(2, 2);
160   vnl_fastops::AB(result_m, m2x10, m10x2);
161   TEST("vnl_fastops::AB(m2x10,m10x2)", result_m, m2x10 * m10x2);
162   vnl_fastops::AtB(result_m, m10x2, m10x2);
163   TEST("vnl_fastops::AtB(m10x2,m10x2)", result_m, m2x10 * m10x2);
164   vnl_fastops::ABt(result_m, m2x10, m2x10);
165   TEST("vnl_fastops::ABt(m2x10,m2x10)", result_m, m2x10 * m10x2);
166   result_m.set_size(10, 10);
167   vnl_fastops::AB(result_m, m10x2, m2x10);
168   TEST("vnl_fastops::AB(m10x2,m2x10)", result_m, m10x2 * m2x10);
169   vnl_fastops::AtB(result_m, m2x10, m2x10);
170   TEST("vnl_fastops::AtB(m2x10,m2x10)", result_m, m10x2 * m2x10);
171   vnl_fastops::ABt(result_m, m10x2, m10x2);
172   TEST("vnl_fastops::ABt(m10x2,m10x2)", result_m, m10x2 * m2x10);
173   result_v.set_size(2);
174   vnl_fastops::AtB(result_v, m10x2, v10);
175   TEST("vnl_fastops::AtB(m10x2,v10)", result_v, m2x10 * v10);
176   result_v.set_size(10);
177   vnl_fastops::AtB(result_v, m2x10, v2);
178   TEST("vnl_fastops::AtB(m2x10,v2)", result_v, m10x2 * v2);
179 
180   // Fourth test: increments and decrements
181   X = m1x1;
182   Y = v1;
183   vnl_fastops::inc_X_by_AtA(X, m1x1);
184   TEST("vnl_fastops::inc_X_by_AtA(X, m1x1)", X, m1x1 * m1x1 + m1x1);
185   vnl_fastops::dec_X_by_AtA(X, m1x1);
186   TEST("vnl_fastops::dec_X_by_AtA(X, m1x1)", X, m1x1);
187   vnl_fastops::inc_X_by_AB(X, m1x1, m1x1);
188   TEST("vnl_fastops::inc_X_by_AB(X, m1x1,m1x1)", X, m1x1 * m1x1 + m1x1);
189   vnl_fastops::dec_X_by_AB(X, m1x1, m1x1);
190   TEST("vnl_fastops::dec_X_by_AB(X, m1x1,m1x1)", X, m1x1);
191   vnl_fastops::inc_X_by_AtB(Y, m1x1, v1);
192   TEST("vnl_fastops::inc_X_by_AtB(X, m1x1,v1)", Y, m1x1 * v1 + v1);
193   vnl_fastops::dec_X_by_AtB(Y, m1x1, v1);
194   TEST("vnl_fastops::dec_X_by_AtB(X, m1x1,v1)", Y, v1);
195   vnl_fastops::inc_X_by_AtB(X, m1x1, m1x1);
196   TEST("vnl_fastops::inc_X_by_AtB(X, m1x1,m1x1)", X, m1x1 * m1x1 + m1x1);
197   vnl_fastops::dec_X_by_AtB(X, m1x1, m1x1);
198   TEST("vnl_fastops::dec_X_by_AtB(X, m1x1,m1x1)", X, m1x1);
199   vnl_fastops::inc_X_by_ABt(X, m1x1, m1x1);
200   TEST("vnl_fastops::inc_X_by_ABt(X, m1x1,m1x1)", X, m1x1 * m1x1 + m1x1);
201   vnl_fastops::dec_X_by_ABt(X, m1x1, m1x1);
202   TEST("vnl_fastops::dec_X_by_ABt(X, m1x1,m1x1)", X, m1x1);
203   X = m2x2;
204   Y = v2;
205   vnl_fastops::inc_X_by_AtA(X, m2x2);
206   TEST("vnl_fastops::inc_X_by_AtA(X, m2x2)", X, m2x2t * m2x2 + m2x2);
207   vnl_fastops::dec_X_by_AtA(X, m2x2);
208   TEST("vnl_fastops::dec_X_by_AtA(X, m2x2)", X, m2x2);
209   vnl_fastops::inc_X_by_AB(X, m2x2, m2x2);
210   TEST("vnl_fastops::inc_X_by_AB(X, m2x2,m2x2)", X, m2x2 * m2x2 + m2x2);
211   vnl_fastops::dec_X_by_AB(X, m2x2, m2x2);
212   TEST("vnl_fastops::dec_X_by_AB(X, m2x2,m2x2)", X, m2x2);
213   vnl_fastops::inc_X_by_AtB(Y, m2x2, v2);
214   TEST("vnl_fastops::inc_X_by_AtB(X, m2x2,v2)", Y, m2x2t * v2 + v2);
215   vnl_fastops::dec_X_by_AtB(Y, m2x2, v2);
216   TEST("vnl_fastops::dec_X_by_AtB(X, m2x2,v2)", Y, v2);
217   vnl_fastops::inc_X_by_AtB(X, m2x2, m2x2);
218   TEST("vnl_fastops::inc_X_by_AtB(X, m2x2,m2x2)", X, m2x2t * m2x2 + m2x2);
219   vnl_fastops::dec_X_by_AtB(X, m2x2, m2x2);
220   TEST("vnl_fastops::dec_X_by_AtB(X, m2x2,m2x2)", X, m2x2);
221   vnl_fastops::inc_X_by_ABt(X, m2x2, m2x2);
222   TEST("vnl_fastops::inc_X_by_ABt(X, m2x2,m2x2)", X, m2x2 * m2x2t + m2x2);
223   vnl_fastops::dec_X_by_ABt(X, m2x2, m2x2);
224   TEST("vnl_fastops::dec_X_by_ABt(X, m2x2,m2x2)", X, m2x2);
225   X = m3x3;
226   Y = v3;
227   vnl_fastops::inc_X_by_AtA(X, m3x3);
228   TEST("vnl_fastops::inc_X_by_AtA(X, m3x3)", X, m3x3t * m3x3 + m3x3);
229   vnl_fastops::dec_X_by_AtA(X, m3x3);
230   TEST("vnl_fastops::dec_X_by_AtA(X, m3x3)", X, m3x3);
231   vnl_fastops::inc_X_by_AB(X, m3x3, m3x3);
232   TEST("vnl_fastops::inc_X_by_AB(X, m3x3,m3x3)", X, m3x3 * m3x3 + m3x3);
233   vnl_fastops::dec_X_by_AB(X, m3x3, m3x3);
234   TEST("vnl_fastops::dec_X_by_AB(X, m3x3,m3x3)", X, m3x3);
235   vnl_fastops::inc_X_by_AtB(Y, m3x3, v3);
236   TEST("vnl_fastops::inc_X_by_AtB(X, m3x3,v3)", Y, m3x3t * v3 + v3);
237   vnl_fastops::dec_X_by_AtB(Y, m3x3, v3);
238   TEST("vnl_fastops::dec_X_by_AtB(X, m3x3,v3)", Y, v3);
239   vnl_fastops::inc_X_by_AtB(X, m3x3, m3x3);
240   TEST("vnl_fastops::inc_X_by_AtB(X, m3x3,m3x3)", X, m3x3t * m3x3 + m3x3);
241   vnl_fastops::dec_X_by_AtB(X, m3x3, m3x3);
242   TEST("vnl_fastops::dec_X_by_AtB(X, m3x3,m3x3)", X, m3x3);
243   vnl_fastops::inc_X_by_ABt(X, m3x3, m3x3);
244   TEST("vnl_fastops::inc_X_by_ABt(X, m3x3,m3x3)", X, m3x3 * m3x3t + m3x3);
245   vnl_fastops::dec_X_by_ABt(X, m3x3, m3x3);
246   TEST("vnl_fastops::dec_X_by_ABt(X, m3x3,m3x3)", X, m3x3);
247   X = m2x2;
248   vnl_fastops::inc_X_by_AB(X, m2x10, m10x2);
249   TEST("vnl_fastops::inc_X_by_AB(X, m2x10,m10x2)", X, m2x10 * m10x2 + m2x2);
250   vnl_fastops::dec_X_by_AB(X, m2x10, m10x2);
251   TEST("vnl_fastops::dec_X_by_AB(X, m2x10,m10x2)", X, m2x2);
252   vnl_fastops::inc_X_by_AtB(X, m10x2, m10x2);
253   TEST("vnl_fastops::inc_X_by_AtB(X, m10x2,m10x2)", X, m2x10 * m10x2 + m2x2);
254   vnl_fastops::dec_X_by_AtB(X, m10x2, m10x2);
255   TEST("vnl_fastops::dec_X_by_AtB(X, m10x2,m10x2)", X, m2x2);
256   vnl_fastops::inc_X_by_ABt(X, m2x10, m2x10);
257   TEST("vnl_fastops::inc_X_by_ABt(X, m2x10,m2x10)", X, m2x10 * m10x2 + m2x2);
258   vnl_fastops::dec_X_by_ABt(X, m2x10, m2x10);
259   TEST("vnl_fastops::dec_X_by_ABt(X, m2x10,m2x10)", X, m2x2);
260   X = m10x2 * m2x10;
261   vnl_fastops::inc_X_by_AB(X, m10x2, m2x10);
262   TEST("vnl_fastops::inc_X_by_AB(X, m10x2,m2x10)", X, m10x2 * m2x10 * 2);
263   vnl_fastops::dec_X_by_AB(X, m10x2, m2x10);
264   TEST("vnl_fastops::dec_X_by_AB(X, m10x2,m2x10)", X, m10x2 * m2x10);
265   vnl_fastops::inc_X_by_AtB(X, m2x10, m2x10);
266   TEST("vnl_fastops::inc_X_by_AtB(X, m2x10,m2x10)", X, m10x2 * m2x10 * 2);
267   vnl_fastops::dec_X_by_AtB(X, m2x10, m2x10);
268   TEST("vnl_fastops::dec_X_by_AtB(X, m2x10,m2x10)", X, m10x2 * m2x10);
269   vnl_fastops::inc_X_by_ABt(X, m10x2, m10x2);
270   TEST("vnl_fastops::inc_X_by_ABt(X, m10x2,m10x2)", X, m10x2 * m2x10 * 2);
271   vnl_fastops::dec_X_by_ABt(X, m10x2, m10x2);
272   TEST("vnl_fastops::dec_X_by_ABt(X, m10x2,m10x2)", X, m10x2 * m2x10);
273   Y = v2;
274   vnl_fastops::inc_X_by_AtB(Y, m10x2, v10);
275   TEST("vnl_fastops::inc_X_by_AtB(X, m10x2,v10)", Y, m2x10 * v10 + v2);
276   vnl_fastops::dec_X_by_AtB(Y, m10x2, v10);
277   TEST("vnl_fastops::dec_X_by_AtB(X, m10x2,v10)", Y, v2);
278   Y = v10;
279   vnl_fastops::inc_X_by_AtB(Y, m2x10, v2);
280   TEST("vnl_fastops::inc_X_by_AtB(X, m2x10,v2)", Y, m10x2 * v2 + v10);
281   vnl_fastops::dec_X_by_AtB(Y, m2x10, v2);
282   TEST("vnl_fastops::dec_X_by_AtB(X, m2x10,v2)", Y, v10);
283 }
284 
285 TESTMAIN(test_fastops);
286