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