Line data Source code
1 : /*
2 : * File: TestMatrix.cpp
3 : * Author: Pantelis Sopasakis
4 : *
5 : * Created on Sep 14, 2015, 3:34:41 PM
6 : *
7 : * ForBES is free software: you can redistribute it and/or modify
8 : * it under the terms of the GNU Lesser General Public License as published by
9 : * the Free Software Foundation, either version 3 of the License, or
10 : * (at your option) any later version.
11 : *
12 : * ForBES is distributed in the hope that it will be useful,
13 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 : * GNU Lesser General Public License for more details.
16 : *
17 : * You should have received a copy of the GNU Lesser General Public License
18 : * along with ForBES. If not, see <http://www.gnu.org/licenses/>.
19 : */
20 :
21 : #include "TestMatrix.h"
22 :
23 :
24 1 : CPPUNIT_TEST_SUITE_REGISTRATION(TestMatrix);
25 :
26 : void _testSubmatrixMultiply(Matrix& A, Matrix& B);
27 :
28 100 : TestMatrix::TestMatrix() {
29 100 : }
30 :
31 200 : TestMatrix::~TestMatrix() {
32 200 : }
33 :
34 100 : void TestMatrix::setUp() {
35 100 : }
36 :
37 100 : void TestMatrix::tearDown() {
38 100 : Matrix::destroy_handle();
39 100 : }
40 :
41 1 : void TestMatrix::testToggleDiagonal() {
42 : // Diagonal to vector
43 1 : size_t n = 10;
44 1 : double alpha = 3.0;
45 1 : Matrix X = MatrixFactory::MakeIdentity(n, alpha);
46 1 : _ASSERT_EQ(Matrix::MATRIX_DIAGONAL, X.getType());
47 1 : X.toggle_diagonal();
48 1 : _ASSERT_EQ(Matrix::MATRIX_DENSE, X.getType());
49 1 : _ASSERT(X.isColumnVector());
50 1 : _ASSERT_EQ(n, X.getNrows());
51 11 : for (size_t i = 0; i < n; i++) {
52 10 : _ASSERT_EQ(alpha, X[i]);
53 : }
54 :
55 2 : Matrix x = MatrixFactory::MakeRandomMatrix(n, 1, 0.0, 1.0);
56 1 : _ASSERT_EQ(Matrix::MATRIX_DENSE, x.getType());
57 1 : _ASSERT(x.isColumnVector());
58 1 : x.toggle_diagonal();
59 1 : _ASSERT_EQ(Matrix::MATRIX_DIAGONAL, x.getType());
60 1 : x.toggle_diagonal();
61 1 : _ASSERT(x.isColumnVector());
62 :
63 :
64 2 : Matrix y = MatrixFactory::MakeRandomMatrix(2, n, 0.0, 1.0);
65 1 : _ASSERT_EXCEPTION(y.toggle_diagonal(), std::invalid_argument);
66 :
67 2 : Matrix z = MatrixFactory::MakeRandomMatrix(n, n, 0.0, 1.0, Matrix::MATRIX_LOWERTR);
68 2 : _ASSERT_EXCEPTION(z.toggle_diagonal(), std::invalid_argument);
69 :
70 :
71 1 : }
72 :
73 1 : void TestMatrix::testGetSet() {
74 1 : const size_t n = 10;
75 1 : Matrix f(n, n);
76 11 : for (size_t i = 0; i < n; i++) {
77 110 : for (size_t j = 0; j < 10; j++) {
78 100 : f.set(i, j, static_cast<double> (3 * i + 5 * j + 13));
79 : }
80 : }
81 11 : for (size_t i = 0; i < n; i++) {
82 110 : for (size_t j = 0; j < n; j++) {
83 100 : _ASSERT_EQ(static_cast<double> (3 * i + 5 * j + 13), f.get(i, j));
84 : }
85 : }
86 :
87 2 : Matrix o;
88 1 : _ASSERT(o.isEmpty());
89 2 : _ASSERT_EXCEPTION(o.get(0, 0), std::out_of_range);
90 1 : }
91 :
92 1 : void TestMatrix::testOpplus() {
93 1 : size_t n = 1000;
94 1 : size_t m = 50;
95 1 : Matrix x = MatrixFactory::MakeRandomMatrix(n, m, -1.5, 2.0);
96 2 : Matrix y(x);
97 :
98 1 : x.plusop();
99 :
100 1 : const double tol = 1e-12;
101 1001 : for (size_t i = 0; i < x.getNrows(); i++) {
102 51000 : for (size_t j = 0; j < x.getNcols(); j++) {
103 50000 : if (y.get(i, j) >= 0.0) {
104 12556 : _ASSERT_NUM_EQ(y.get(i, j), x.get(i, j), tol);
105 : } else {
106 37444 : _ASSERT_NUM_EQ(0.0, x.get(i, j), tol);
107 : }
108 : }
109 1 : }
110 1 : }
111 :
112 1 : void TestMatrix::testOpplus2() {
113 1 : size_t n = 10;
114 1 : size_t m = 1000;
115 1 : Matrix x = MatrixFactory::MakeRandomMatrix(n, m, -1.5, 2.0);
116 :
117 2 : Matrix z(n, m);
118 1 : x.plusop(&z);
119 :
120 1 : const double tol = 1e-12;
121 11 : for (size_t i = 0; i < x.getNrows(); i++) {
122 10010 : for (size_t j = 0; j < x.getNcols(); j++) {
123 10000 : if (x.get(i, j) >= 0.0) {
124 2480 : _ASSERT_NUM_EQ(x.get(i, j), z.get(i, j), tol);
125 : } else {
126 7520 : _ASSERT_NUM_EQ(0.0, z.get(i, j), tol);
127 : }
128 : }
129 1 : }
130 1 : }
131 :
132 1 : void TestMatrix::testOpplusSparse() {
133 1 : size_t n = 1000;
134 1 : size_t m = 50;
135 1 : size_t nnz = 200;
136 1 : Matrix x = MatrixFactory::MakeRandomSparse(n, m, nnz, -1.5, 2.0);
137 2 : Matrix y(x);
138 :
139 1 : x.plusop();
140 :
141 1 : const double tol = 1e-12;
142 1001 : for (size_t i = 0; i < x.getNrows(); i++) {
143 51000 : for (size_t j = 0; j < x.getNcols(); j++) {
144 50000 : if (y.get(i, j) >= 0.0) {
145 49847 : _ASSERT_NUM_EQ(y.get(i, j), x.get(i, j), tol);
146 : } else {
147 153 : _ASSERT_NUM_EQ(0.0, x.get(i, j), tol);
148 : }
149 : }
150 1 : }
151 1 : }
152 :
153 1 : void TestMatrix::testQuadratic() {
154 : /* Test quadratic with diagonal matrices */
155 : Matrix *f;
156 8 : for (size_t n = 5; n < 12; n++) {
157 7 : Matrix *x = new Matrix(n, 1);
158 7 : f = new Matrix(n, n);
159 63 : for (size_t i = 0; i < n; i++) {
160 56 : _ASSERT_OK(f -> set(i, i, 1.0));
161 56 : (*x)[i] = i + 1;
162 : }
163 7 : double r = f -> quad(*x);
164 7 : _ASSERT_EQ(static_cast<double> (n * (n + 1)* (2 * n + 1) / 6), 2 * r);
165 7 : _ASSERT_OK(delete f);
166 7 : _ASSERT_OK(delete x);
167 : }
168 1 : Matrix A(5, 6);
169 2 : Matrix y(6, 1);
170 1 : _ASSERT_EXCEPTION(A.quad(y), std::invalid_argument);
171 :
172 2 : Matrix B(5, 5);
173 2 : Matrix z(5, 2);
174 1 : _ASSERT_EXCEPTION(B.quad(z), std::invalid_argument);
175 1 : _ASSERT_EXCEPTION(B.quad(z, z), std::invalid_argument);
176 :
177 2 : Matrix C(5, 5);
178 2 : Matrix s(7, 1);
179 2 : _ASSERT_EXCEPTION(C.quad(s), std::invalid_argument);
180 :
181 1 : }
182 :
183 1 : void TestMatrix::testQuadratic2() {
184 1 : double fdata[9] = {-1.0, 3.0, 1.0, 2.0, -1.0, -1.0, 5.0, 2.0, -5.0};
185 1 : double xdata[3] = {1.0, 2.0, 3.0};
186 :
187 1 : Matrix f(3, 3, fdata);
188 2 : Matrix x(3, 1, xdata);
189 :
190 : double r;
191 1 : _ASSERT_OK(r = f.quad(x));
192 2 : _ASSERT_EQ(-8.0, r);
193 1 : }
194 :
195 1 : void TestMatrix::testQuadratic3() {
196 1 : double fdata[9] = {-1.0, -3.0, 7.5, 2.0, -1.0, -1.0, 5.0, 2.0, -5.0};
197 1 : double xdata[3] = {-1.5, 2.0, 3.0};
198 1 : double qdata[3] = {5.0, -6.0, 13.5};
199 :
200 1 : Matrix f(3, 3, fdata);
201 2 : Matrix x(3, 1, xdata);
202 2 : Matrix q(3, 1, qdata);
203 :
204 :
205 : double r;
206 1 : _ASSERT_OK(r = f.quad(x, q));
207 2 : _ASSERT_EQ(-28.25, r);
208 1 : }
209 :
210 1 : void TestMatrix::testQuadraticDot() {
211 1 : double ydata[4] = {-2.0, 5.0, -6.0, -11.0};
212 1 : double xdata[4] = {10.0, 2.0, 3.0, 4.0};
213 :
214 1 : const size_t n = 4;
215 1 : const size_t m = 1;
216 :
217 1 : Matrix y(n, m, ydata);
218 2 : Matrix x(n, m, xdata);
219 :
220 2 : Matrix r = y*x;
221 :
222 1 : _ASSERT_EQ(m, r.getNcols());
223 1 : _ASSERT_EQ(m, r.getNrows());
224 1 : _ASSERT_EQ(m, r.length());
225 2 : _ASSERT_EQ(-72.0, r[0]);
226 1 : }
227 :
228 1 : void TestMatrix::test_MDD1() {
229 1 : double fdata[9] = {-1.0, 3.0, 1.0, 2.0, -1.0, -1.0, 5.0, 2.0, -5.0};
230 1 : double xdata[3] = {1.0, 2.0, 3.0};
231 :
232 1 : const size_t n = 3;
233 1 : const size_t m = 1;
234 :
235 1 : Matrix f(n, n, fdata);
236 2 : Matrix x(n, m, xdata);
237 :
238 2 : Matrix r;
239 1 : r = f*x;
240 :
241 1 : _ASSERT_EQ(n, r.getNrows());
242 1 : _ASSERT_EQ(m, r.getNcols());
243 1 : _ASSERT_EQ(n, r.length());
244 1 : _ASSERT_EQ(18.0, r[0]);
245 1 : _ASSERT_EQ(7.0, r[1]);
246 1 : _ASSERT_EQ(-16.0, r[2]);
247 :
248 2 : Matrix A(5, 6);
249 2 : Matrix B(3, 5);
250 2 : Matrix C;
251 2 : _ASSERT_EXCEPTION(C = A*B, std::invalid_argument);
252 1 : }
253 :
254 1 : void TestMatrix::testAssignment() {
255 1 : const size_t nRows = 3;
256 1 : const size_t nCols = 2;
257 1 : const size_t size = nRows * nCols;
258 :
259 1 : Matrix f = MatrixFactory::MakeRandomMatrix(nRows, nCols, 0.0, 1.0, Matrix::MATRIX_DENSE);
260 2 : Matrix g;
261 1 : g = f;
262 :
263 1 : _ASSERT_EQ(Matrix::MATRIX_DENSE, g.getType());
264 1 : _ASSERT_EQ(size, g.length());
265 7 : for (size_t i = 0; i < size; i++) {
266 6 : _ASSERT(g[i] >= 0);
267 6 : _ASSERT(g[i] <= 1);
268 : }
269 2 : _ASSERT_OK(f = f);
270 1 : }
271 :
272 1 : void TestMatrix::testAdditionBad() {
273 1 : Matrix A(5, 6);
274 2 : Matrix B(7, 8);
275 2 : Matrix C;
276 1 : CPPUNIT_ASSERT_THROW(C = A + B, std::invalid_argument);
277 1 : CPPUNIT_ASSERT_THROW(C = A - B, std::invalid_argument);
278 1 : CPPUNIT_ASSERT_THROW(A += B, std::invalid_argument);
279 2 : CPPUNIT_ASSERT_THROW(A -= B, std::invalid_argument);
280 1 : }
281 :
282 1 : void TestMatrix::test_ADD1() {
283 1 : const size_t nRows = 3;
284 1 : const size_t nCols = 2;
285 1 : const size_t size = nRows * nCols;
286 : double *a, *b;
287 1 : a = new double[size];
288 1 : b = new double[size];
289 7 : for (size_t i = 0; i < size; i++) {
290 6 : a[i] = i;
291 6 : b[i] = 3 * i + 7;
292 : }
293 :
294 1 : Matrix A(nRows, nCols, a);
295 2 : Matrix B(nRows, nCols, b);
296 :
297 :
298 2 : Matrix C;
299 1 : C = A + B;
300 :
301 7 : for (size_t i = 0; i < size; i++) {
302 6 : _ASSERT_EQ(a[i] + b[i], C[i]);
303 : }
304 :
305 2 : Matrix D = A + B + C;
306 7 : for (size_t i = 0; i < size; i++) {
307 6 : _ASSERT_EQ(2 * C[i], D[i]);
308 : }
309 1 : delete[] a;
310 2 : delete[] b;
311 1 : }
312 :
313 1 : void TestMatrix::testFBMatrix() {
314 : /* Test FBMatrix() - empty constructor */
315 1 : Matrix *fBMatrix = new Matrix();
316 1 : _ASSERT_EQ(static_cast<size_t> (0), fBMatrix->getNcols());
317 1 : _ASSERT_EQ(static_cast<size_t> (0), fBMatrix->getNrows());
318 1 : delete fBMatrix;
319 :
320 : /* Test FBMatrix(size_t, size_t, double*) - Provide data */
321 1 : const size_t n = 5;
322 1 : double *x = new double[n];
323 6 : for (size_t i = 0; i < n; i++) {
324 5 : x[i] = 1 + 7 * i;
325 : }
326 1 : Matrix f(n, 1, x);
327 1 : delete[] x;
328 6 : for (size_t i = 0; i < n; i++) {
329 5 : _ASSERT_EQ(static_cast<double> (1 + 7 * i), f[i]);
330 : }
331 :
332 1 : double s = 0.0;
333 :
334 1 : _ASSERT_EXCEPTION(fBMatrix = new Matrix(3, 4, Matrix::MATRIX_DIAGONAL), std::invalid_argument);
335 1 : _ASSERT_EXCEPTION(fBMatrix = new Matrix(3, 4, Matrix::MATRIX_SYMMETRIC), std::invalid_argument);
336 1 : _ASSERT_EXCEPTION(fBMatrix = new Matrix(3, 4, Matrix::MATRIX_LOWERTR), std::invalid_argument);
337 : // _ASSERT_EXCEPTION(s = f[-1], std::out_of_range);
338 : // _ASSERT_EXCEPTION(s = f[n], std::out_of_range);
339 1 : _ASSERT_NUM_EQ(0.0, s, 1e-10);
340 1 : _ASSERT_OK(Matrix::destroy_handle());
341 1 : _ASSERT_EQ(0, Matrix::destroy_handle());
342 :
343 2 : Matrix E;
344 2 : Matrix T(E);
345 2 : _ASSERT(T.isEmpty());
346 1 : }
347 :
348 1 : void TestMatrix::testMakeRandomFBMatrix() {
349 1 : const size_t nRows = 10;
350 1 : const size_t nCols = 20;
351 1 : const double offset = 0.1;
352 1 : const double scale = 3.5;
353 1 : Matrix f = MatrixFactory::MakeRandomMatrix(nRows, nCols, offset, scale, Matrix::MATRIX_DENSE);
354 :
355 1 : _ASSERT_EQ(nCols, f.getNcols());
356 1 : _ASSERT_EQ(nRows, f.getNrows());
357 201 : for (size_t i = 0; i < nRows * nCols; i++) {
358 200 : if (i > 0) {
359 199 : _ASSERT_NEQ(f[i], f[i - 1]);
360 : }
361 200 : _ASSERT(f[i] >= offset);
362 200 : _ASSERT(f[i] <= offset + scale);
363 1 : }
364 1 : }
365 :
366 1 : void TestMatrix::testGetData() {
367 1 : const size_t n = 100;
368 1 : double *myData = new double[n];
369 101 : for (size_t j = 0; j < n; j++) {
370 100 : myData[j] = j;
371 : }
372 1 : Matrix * mat = new Matrix(n, 1, myData);
373 :
374 : double * retrievedData;
375 1 : retrievedData = mat->getData();
376 :
377 101 : for (size_t j = 0; j < n; j++) {
378 100 : _ASSERT_EQ(j, static_cast<size_t> (retrievedData[j]));
379 : }
380 :
381 1 : retrievedData[0] = 666;
382 1 : _ASSERT_EQ(666, static_cast<int> ((*mat)[0]));
383 :
384 1 : delete(mat);
385 1 : }
386 :
387 1 : void TestMatrix::testGetNcols() {
388 1 : size_t y = static_cast<size_t> (5 + 50 * static_cast<double> (std::rand()) / static_cast<double> (RAND_MAX));
389 1 : Matrix mat(10, y);
390 1 : _ASSERT_EQ(y, mat.getNcols());
391 1 : }
392 :
393 1 : void TestMatrix::testGetNrows() {
394 1 : size_t r = static_cast<size_t> (5 + 50 * static_cast<double> (std::rand()) / static_cast<double> (RAND_MAX));
395 1 : Matrix mat(r, 10);
396 1 : _ASSERT_EQ(r, mat.getNrows());
397 :
398 1 : const size_t n = 5;
399 1 : Matrix *gm = new Matrix(n, n);
400 1 : _ASSERT_NOT(gm->isEmpty());
401 1 : _ASSERT_EQ(n, gm->getNrows());
402 1 : _ASSERT_EQ(n, gm->getNcols());
403 1 : delete gm;
404 1 : }
405 :
406 1 : void TestMatrix::testIsColumnVector() {
407 1 : Matrix rowVec(2345, 1);
408 1 : _ASSERT(rowVec.isColumnVector());
409 1 : }
410 :
411 1 : void TestMatrix::testIsRowVector() {
412 1 : Matrix rowVec(1, 100);
413 1 : _ASSERT(rowVec.isRowVector());
414 1 : }
415 :
416 1 : void TestMatrix::testLength() {
417 1 : size_t nRep = 20;
418 :
419 21 : for (size_t i = 0; i < nRep; i++) {
420 20 : size_t x = static_cast<size_t> (5 + 50 * static_cast<double> (std::rand()) / static_cast<double> (RAND_MAX));
421 20 : size_t y = static_cast<size_t> (5 + 50 * static_cast<double> (std::rand()) / static_cast<double> (RAND_MAX));
422 20 : Matrix * f = new Matrix(x, y);
423 20 : _ASSERT_EQ(x*y, f->length());
424 20 : delete f;
425 : }
426 :
427 1 : }
428 :
429 1 : void TestMatrix::testReshape() {
430 1 : size_t n = 45;
431 1 : size_t m = 56;
432 1 : Matrix f(n, m);
433 1 : size_t k = 5;
434 1 : int status = f.reshape(k, k);
435 1 : _ASSERT_EQ(0, status);
436 1 : _ASSERT_EQ(k, f.getNcols());
437 1 : }
438 :
439 1 : void TestMatrix::testReshapeBad() {
440 1 : size_t n = 45;
441 1 : size_t m = 56;
442 1 : Matrix f(n, m);
443 1 : int status = f.reshape(n, m + 1);
444 1 : _ASSERT_EQ(-2, status);
445 1 : _ASSERT_EQ(n, f.getNrows());
446 1 : _ASSERT_EQ(m, f.getNcols());
447 :
448 1 : status = f.reshape(0, 1);
449 1 : _ASSERT_EQ(-1, status);
450 1 : }
451 :
452 1 : void TestMatrix::testDiagonalGetSet() {
453 1 : const size_t n = 10;
454 1 : double *myData = new double[n];
455 11 : for (size_t j = 0; j < n; j++) {
456 10 : myData[j] = j + 1.0;
457 : }
458 :
459 1 : Matrix *A = new Matrix(n, n, myData, Matrix::MATRIX_DIAGONAL);
460 1 : delete[] myData;
461 11 : for (size_t i = 0; i < n; i++) {
462 110 : for (size_t j = 0; j < n; j++) {
463 100 : _ASSERT_EQ(static_cast<double> (i + 1)*(i == j), A -> get(i, j));
464 : }
465 : }
466 :
467 1 : _ASSERT_EQ(Matrix::MATRIX_DIAGONAL, A -> getType());
468 :
469 1 : double t = -1.234;
470 1 : _ASSERT_OK(A -> set(3, 3, t));
471 1 : _ASSERT_EQ(t, A -> get(3, 3));
472 1 : _ASSERT_EXCEPTION(A -> set(3, 4, 1.0), std::invalid_argument);
473 1 : _ASSERT_EXCEPTION(A -> get(n, n - 1), std::out_of_range);
474 1 : _ASSERT_EXCEPTION(A -> set(n, n, 100.0), std::out_of_range);
475 :
476 1 : delete A;
477 1 : }
478 :
479 1 : void TestMatrix::testDiagonalMultiplication() {
480 : /* Diag * Dense = Dense */
481 1 : const size_t n = 10;
482 1 : const size_t m = 3;
483 1 : double *myData = new double[n];
484 11 : for (size_t j = 0; j < n; j++) {
485 10 : myData[j] = j + 1.0;
486 : }
487 1 : Matrix *A = new Matrix(n, n, myData, Matrix::MATRIX_DIAGONAL);
488 1 : Matrix B(n, m);
489 :
490 11 : for (size_t i = 0; i < n; i++) {
491 40 : for (size_t j = 0; j < m; j++) {
492 30 : B.set(i, j, i + 1.0);
493 : }
494 : }
495 :
496 2 : Matrix C;
497 1 : _ASSERT_OK(C = (*A) * B);
498 1 : _ASSERT_EQ(Matrix::MATRIX_DENSE, C.getType());
499 1 : _ASSERT_EQ(n, C.getNrows());
500 1 : _ASSERT_EQ(m, C.getNcols());
501 1 : _ASSERT_NOT(C.isEmpty());
502 1 : _ASSERT_EQ(n*m, C.length());
503 :
504 11 : for (size_t i = 0; i < n; i++) {
505 40 : for (size_t j = 0; j < m; j++) {
506 30 : _ASSERT_EQ(std::pow(i + 1.0, 2.0), C.get(i, j));
507 : }
508 : }
509 :
510 1 : _ASSERT_OK(delete A);
511 2 : delete[] myData;
512 1 : }
513 :
514 1 : void TestMatrix::testDiagonalMultiplication2() {
515 : /* Diag * Diag = Diag */
516 1 : const size_t n = 10;
517 1 : double *aData = new double[n];
518 1 : double *bData = new double[n];
519 :
520 11 : for (size_t i = 0; i < n; i++) {
521 10 : aData[i] = i + 1.0;
522 10 : bData[i] = n - i;
523 : }
524 1 : Matrix *A = new Matrix(n, n, aData, Matrix::MATRIX_DIAGONAL);
525 1 : Matrix *B = new Matrix(n, n, bData, Matrix::MATRIX_DIAGONAL);
526 :
527 1 : delete[] aData;
528 1 : delete[] bData;
529 :
530 1 : Matrix C = (*A) * (*B);
531 1 : _ASSERT_EQ(Matrix::MATRIX_DIAGONAL, A -> getType());
532 1 : _ASSERT_EQ(Matrix::MATRIX_DIAGONAL, B -> getType());
533 1 : _ASSERT_EQ(Matrix::MATRIX_DIAGONAL, C.getType());
534 :
535 11 : for (size_t i = 0; i < n; i++) {
536 10 : _ASSERT_EQ((i + 1.0)*(n - i), C.get(i, i));
537 1 : }
538 1 : }
539 :
540 1 : void TestMatrix::testDenseTimesDiagonal() {
541 1 : const size_t nRows = 7;
542 1 : const size_t nCols = 3;
543 1 : const size_t size = nRows * nCols;
544 : double *aData;
545 : double *bData;
546 :
547 1 : aData = new double[size];
548 22 : for (size_t i = 0; i < size; i++) {
549 21 : aData[i] = i;
550 : }
551 :
552 1 : bData = new double[nCols];
553 4 : for (size_t i = 0; i < nCols; i++) {
554 3 : bData[i] = 3.0 * (i + 1);
555 : }
556 :
557 :
558 1 : Matrix A(nRows, nCols, aData, Matrix::MATRIX_DENSE);
559 2 : Matrix B(nCols, nCols, bData, Matrix::MATRIX_DIAGONAL);
560 :
561 2 : Matrix C = A * B;
562 :
563 8 : for (size_t i = 0; i < C.getNrows(); i++) {
564 28 : for (size_t j = 0; j < C.getNcols(); j++) {
565 21 : _ASSERT_EQ(A.get(i, j) * B.get(j, j), C.get(i, j));
566 : }
567 : }
568 :
569 1 : delete[] aData;
570 2 : delete[] bData;
571 1 : }
572 :
573 1 : void TestMatrix::testQuadDiagonal() {
574 1 : const size_t n = 10;
575 1 : double *aData = new double[n];
576 1 : double *xData = new double[n];
577 11 : for (size_t i = 0; i < n; i++) {
578 10 : aData[i] = i + 1.0;
579 10 : xData[i] = 3.0 * (i + 1.0);
580 : }
581 1 : Matrix *A = new Matrix(n, n, aData, Matrix::MATRIX_DIAGONAL);
582 1 : Matrix *x = new Matrix(n, 1, xData);
583 1 : double val = A -> quad(*x, *x);
584 :
585 1 : const double correct = 17077.5;
586 1 : _ASSERT_EQ(correct, val);
587 :
588 1 : val = A -> quad(*x);
589 1 : const double correct2 = 13612.5;
590 1 : _ASSERT_EQ(correct2, val);
591 :
592 1 : delete A;
593 1 : delete x;
594 1 : _ASSERT_OK(delete[] aData);
595 1 : _ASSERT_OK(delete[] xData);
596 1 : }
597 :
598 1 : void TestMatrix::testSubtract() {
599 1 : const size_t n = 10;
600 1 : Matrix X = MatrixFactory::MakeRandomMatrix(n, n, 0.0, 1.0, Matrix::MATRIX_DENSE);
601 2 : Matrix Y = X - X;
602 101 : for (size_t i = 0; i < n * n; i++) {
603 100 : _ASSERT_EQ(0.0, Y.get(i));
604 : }
605 :
606 1 : X = MatrixFactory::MakeRandomMatrix(n, n, 0.0, 1.0, Matrix::MATRIX_DENSE);
607 1 : X -= X;
608 101 : for (size_t i = 0; i < n * n; i++) {
609 100 : _ASSERT_EQ(0.0, X.get(i));
610 : }
611 :
612 1 : X = MatrixFactory::MakeRandomMatrix(n, n, 0.0, 1.0, Matrix::MATRIX_DENSE);
613 2 : Matrix B = MatrixFactory::MakeRandomMatrix(n, n, 0.0, 1.0, Matrix::MATRIX_DENSE);
614 2 : Matrix Z = X - B;
615 :
616 11 : for (size_t i = 0; i < n; i++) {
617 110 : for (size_t j = 0; j < n; j++) {
618 100 : _ASSERT_EQ(X.get(i, j) - B.get(i, j), Z.get(i, j));
619 : }
620 1 : }
621 :
622 1 : }
623 :
624 1 : void TestMatrix::testLowerTriangular_getSet() {
625 20 : for (size_t n = 1; n < 20; n++) {
626 19 : Matrix *A = new Matrix(n, n, Matrix::MATRIX_LOWERTR);
627 :
628 209 : for (size_t i = 0; i < n; i++) {
629 1520 : for (size_t j = 0; j <= i; j++) {
630 1330 : A -> set(i, j, 3.2 * i + 7.5 * j + 0.45);
631 : }
632 : }
633 :
634 19 : const double tol = 1e-6;
635 209 : for (size_t i = 0; i < n; i++) {
636 2660 : for (size_t j = 0; j < n; j++) {
637 2470 : if (i >= j) {
638 1330 : _ASSERT_NUM_EQ(3.2 * i + 7.5 * j + 0.45, A->get(i, j), tol);
639 : } else {
640 1140 : _ASSERT_EQ(0.0, A->get(i, j));
641 : }
642 : }
643 : }
644 :
645 19 : size_t exp_length = n * (n + 1) / 2;
646 19 : _ASSERT_EQ(exp_length, A->length());
647 19 : delete A;
648 : }
649 1 : }
650 :
651 1 : void TestMatrix::testLowerTriangularTraspose_getSet() {
652 1 : const size_t n = 10;
653 1 : Matrix A = MatrixFactory::MakeRandomMatrix(n, n, 0.0, 1.0, Matrix::MATRIX_LOWERTR);
654 2 : Matrix AT(A);
655 1 : AT.transpose();
656 11 : for (size_t i = 0; i < n; i++) {
657 110 : for (size_t j = 0; j < n; j++) {
658 100 : _ASSERT_EQ(A.get(i, j), AT.get(j, i));
659 : }
660 : }
661 1 : AT.transpose();
662 2 : _ASSERT_EQ(A, AT);
663 1 : }
664 :
665 1 : void TestMatrix::testSymmetric_getSet() {
666 1 : const size_t n = 4;
667 1 : Matrix *A = new Matrix(n, n, Matrix::MATRIX_SYMMETRIC);
668 1 : _ASSERT(A->isSymmetric());
669 5 : for (size_t i = 0; i < n; i++) {
670 14 : for (size_t j = 0; j <= i; j++) {
671 10 : A -> set(i, j, 3.2 * i + 0.2 * j + 0.45);
672 : }
673 : }
674 1 : const double tol = 1e-6;
675 5 : for (size_t i = 0; i < n; i++) {
676 20 : for (size_t j = 0; j < n; j++) {
677 16 : if (i >= j) {
678 10 : _ASSERT_NUM_EQ(3.2f * i + 0.2f * j + 0.45f, A->get(i, j), tol);
679 : } else {
680 6 : _ASSERT_NUM_EQ(A -> get(j, i), A->get(i, j), tol);
681 : }
682 : }
683 : }
684 1 : size_t exp_length = n * (n + 1) / 2;
685 1 : _ASSERT_EQ(exp_length, A->length());
686 :
687 1 : delete A;
688 1 : }
689 :
690 1 : void TestMatrix::testTranspose() {
691 1 : size_t n = 5;
692 1 : size_t m = 6;
693 1 : Matrix A(n, m);
694 1 : A.transpose();
695 1 : _ASSERT_EQ(n, A.getNcols());
696 1 : _ASSERT_EQ(m, A.getNrows());
697 1 : A.transpose();
698 1 : _ASSERT_EQ(n, A.getNrows());
699 1 : _ASSERT_EQ(m, A.getNcols());
700 :
701 2 : Matrix X(n, m);
702 6 : for (size_t i = 0; i < n; i++) {
703 35 : for (size_t j = 0; j < m; j++) {
704 30 : X.set(i, j, 7.5 * i - 2.8 * j - 1.0);
705 : }
706 : }
707 1 : X.transpose();
708 7 : for (size_t i = 0; i < m; i++) {
709 36 : for (size_t j = 0; j < n; j++) {
710 30 : _ASSERT_OK(X.get(i, j));
711 : }
712 : }
713 :
714 2 : Matrix Y(n, m);
715 6 : for (size_t i = 0; i < n; i++) {
716 35 : for (size_t j = 0; j < m; j++) {
717 30 : Y.set(i, j, 7.5 * i - 2.8 * j - 1.0);
718 : }
719 : }
720 :
721 2 : Matrix YT(m, n); // Construct the transpose of Y as YT
722 7 : for (size_t i = 0; i < m; i++) {
723 36 : for (size_t j = 0; j < n; j++) {
724 30 : YT.set(i, j, Y.get(j, i));
725 : }
726 : }
727 :
728 1 : Y.transpose();
729 7 : for (size_t i = 0; i < m; i++) {
730 36 : for (size_t j = 0; j < n; j++) {
731 30 : _ASSERT_EQ(YT.get(i, j), Y.get(i, j));
732 : }
733 1 : }
734 :
735 1 : }
736 :
737 1 : void TestMatrix::test_MXH() {
738 1 : const size_t n = 10;
739 1 : Matrix D(n, n, Matrix::MATRIX_DIAGONAL);
740 11 : for (size_t i = 0; i < n; i++) {
741 10 : D[i] = i + 1.0;
742 : }
743 :
744 2 : Matrix S(n, n, Matrix::MATRIX_SYMMETRIC);
745 1 : _ASSERT(S.isSymmetric());
746 11 : for (size_t i = 0; i < n; i++) {
747 65 : for (size_t j = i; j < n; j++) {
748 55 : S.set(i, j, -3.1 * i + 3.25 * j + 5.35);
749 : }
750 : }
751 :
752 2 : Matrix DS = D * S;
753 1 : _ASSERT_EQ(Matrix::MATRIX_SYMMETRIC, DS.getType());
754 11 : for (size_t i = 0; i < n; i++) {
755 65 : for (size_t j = i; j < n; j++) {
756 55 : _ASSERT_EQ((i + 1.0)*(-3.1 * i + 3.25 * j + 5.35), DS.get(i, j));
757 : }
758 1 : }
759 1 : }
760 :
761 1 : void TestMatrix::test_MDX() { // dense * diagonal
762 1 : const size_t n = 10;
763 1 : Matrix D = MatrixFactory::MakeRandomMatrix(n, n, 0.0, 2.0, Matrix::MATRIX_DENSE);
764 2 : Matrix X = MatrixFactory::MakeRandomMatrix(n, n, 0.0, 2.0, Matrix::MATRIX_DIAGONAL);
765 2 : Matrix result = D*X;
766 2 : Matrix XX(n, n, Matrix::MATRIX_DENSE);
767 11 : for (size_t i = 0; i < n; i++) {
768 110 : for (size_t j = 0; j < n; j++) {
769 100 : XX.set(i, j, X.get(i, j));
770 : }
771 : }
772 2 : Matrix correct = D*XX;
773 1 : _ASSERT_EQ(correct, result);
774 :
775 1 : result = X*D;
776 1 : correct = XX*D;
777 1 : _ASSERT_EQ(correct, result);
778 :
779 1 : const size_t m = 12;
780 1 : D = MatrixFactory::MakeRandomMatrix(n, m, 0.0, 2.0, Matrix::MATRIX_DENSE);
781 1 : X = MatrixFactory::MakeRandomMatrix(m, m, 0.0, 2.0, Matrix::MATRIX_DIAGONAL);
782 1 : XX = Matrix(m, m, Matrix::MATRIX_DENSE);
783 13 : for (size_t i = 0; i < m; i++) {
784 156 : for (size_t j = 0; j < m; j++) {
785 144 : XX.set(i, j, X.get(i, j));
786 : }
787 : }
788 :
789 1 : correct = D*XX;
790 1 : result = D*X;
791 1 : _ASSERT_EQ(correct, result);
792 :
793 1 : _ASSERT_EQ(n, D.getNrows());
794 1 : _ASSERT_EQ(m, D.getNcols());
795 1 : D.transpose();
796 1 : _ASSERT_EQ(m, D.getNrows());
797 1 : _ASSERT_EQ(n, D.getNcols()); /* Now D is m-by-n */
798 1 : X = MatrixFactory::MakeRandomMatrix(n, n, 0.0, 2.0, Matrix::MATRIX_DIAGONAL); /* X is n-by-n */
799 1 : result = D*X;
800 11 : for (size_t j = 0; j < n; j++) {
801 110 : for (size_t i = 0; i < n; i++) {
802 100 : _ASSERT_NUM_EQ(D.get(i, j) * X.get(j, j), result.get(i, j), 1e-8);
803 : }
804 1 : }
805 1 : }
806 :
807 1 : void TestMatrix::test_MSX() { /* Sparse * Diagonal */
808 1 : const size_t n = 10;
809 1 : Matrix X(n, n, Matrix::MATRIX_DIAGONAL);
810 11 : for (size_t i = 0; i < n; i++) {
811 10 : X.set(i, i, i + 1.0);
812 : }
813 2 : Matrix S = MatrixFactory::MakeRandomSparse(n, n, static_cast<size_t> (2.5 * n), 0.0, 1.0);
814 2 : Matrix R = S*X;
815 1 : _ASSERT_EQ(Matrix::MATRIX_SPARSE, R.getType());
816 11 : for (size_t j = 0; j < n; j++) {
817 110 : for (size_t i = 0; i < n; i++) {
818 100 : _ASSERT_NUM_EQ(S.get(i, j) * X.get(j, j), R.get(i, j), 1e-8);
819 : }
820 1 : }
821 1 : }
822 :
823 1 : void TestMatrix::test_MXL() {
824 1 : size_t n = 6;
825 1 : Matrix D(n, n, Matrix::MATRIX_DIAGONAL);
826 7 : for (size_t i = 0; i < n; i++) {
827 6 : D[i] = i + 1.0;
828 : }
829 :
830 2 : Matrix L(n, n, Matrix::MATRIX_LOWERTR);
831 7 : for (size_t i = 0; i < n; i++) {
832 27 : for (size_t j = 0; j <= i; j++) {
833 21 : L.set(i, j, -0.1 * i + 0.45 * j + 1.01);
834 : }
835 : }
836 :
837 2 : Matrix DL = D * L;
838 1 : _ASSERT_EQ(Matrix::MATRIX_LOWERTR, DL.getType());
839 1 : _ASSERT_EQ(n, DL.getNcols());
840 1 : _ASSERT_EQ(n, DL.getNrows());
841 1 : _ASSERT_NOT(DL.isEmpty());
842 :
843 7 : for (size_t i = 0; i < n; i++) {
844 27 : for (size_t j = 0; j <= i; j++) {
845 21 : _ASSERT_EQ(0.0, DL.get(i, j) - D[i] * L.get(i, j));
846 : }
847 1 : }
848 1 : }
849 :
850 1 : void TestMatrix::test_MDH() {
851 1 : const size_t n = 4;
852 : double a[n * n] = {5, 11, -2, -1,
853 : 6, 3, 7, -1,
854 : -21, 0, 13, -1,
855 1 : -18, -1, -17, 30};
856 1 : Matrix A(n, n, a, Matrix::MATRIX_DENSE);
857 :
858 2 : Matrix S(n, n, Matrix::MATRIX_SYMMETRIC);
859 5 : for (size_t i = 0; i < n; i++) {
860 14 : for (size_t j = 0; j <= i; j++) {
861 10 : S.set(i, j, 3 * i + 2 * j + 4.0);
862 : }
863 : }
864 :
865 2 : Matrix AS = A*S;
866 1 : _ASSERT_NOT(AS.isEmpty());
867 1 : _ASSERT_NOT(AS.isColumnVector());
868 1 : _ASSERT_NOT(AS.isRowVector());
869 1 : _ASSERT_EQ(Matrix::MATRIX_DENSE, AS.getType());
870 1 : _ASSERT_EQ(n, AS.getNcols());
871 1 : _ASSERT_EQ(n, AS.getNrows());
872 :
873 : double asData[n * n] = {-382.0, 52.0, -50.0, 369.0,
874 : -433.0, 89.0, -50.0, 422.0,
875 : -478.0, 129.0, -43.0, 474.0,
876 1 : -544.0, 169.0, -23.0, 525.0};
877 2 : Matrix AS_correct(n, n, asData, Matrix::MATRIX_DENSE);
878 :
879 5 : for (size_t i = 0; i < n; i++) {
880 20 : for (size_t j = 0; j < n; j++) {
881 16 : _ASSERT_EQ(AS_correct.get(i, j), AS.get(i, j));
882 : }
883 1 : }
884 1 : }
885 :
886 1 : void TestMatrix::test_MDL() {
887 1 : const size_t n = 4;
888 : double a[n * n] = {5, 11, -2, -1,
889 : 6, 3, 7, -1,
890 : -21, 0, 13, -1,
891 1 : -18, -1, -17, 30};
892 1 : Matrix A(n, n, a, Matrix::MATRIX_DENSE);
893 :
894 2 : Matrix L(n, n, Matrix::MATRIX_LOWERTR);
895 5 : for (size_t i = 0; i < n; i++) {
896 14 : for (size_t j = 0; j <= i; j++) {
897 10 : L.set(i, j, 3 * i + 2 * j + 4.0);
898 : }
899 : }
900 :
901 2 : Matrix AL = A*L;
902 : double alCorrect[n * n] = {-382.0, 52.0, -50.0, 369.0, -468.0, 12.0, -36.0, 429.0,
903 1 : -600.0, -17.0, -107.0, 496.0, -342.0, -19.0, -323.0, 570.0};
904 2 : Matrix AL_correct(n, n, alCorrect, Matrix::MATRIX_DENSE);
905 :
906 5 : for (size_t i = 0; i < n; i++) {
907 14 : for (size_t j = 0; j <= i; j++) {
908 10 : _ASSERT_EQ(AL.get(i, j), AL_correct.get(i, j));
909 : }
910 1 : }
911 1 : }
912 :
913 1 : void TestMatrix::testQuadSymmetric() {
914 1 : const size_t n = 10;
915 :
916 1 : Matrix A(n, n, Matrix::MATRIX_SYMMETRIC);
917 2 : Matrix x(n, 1);
918 :
919 11 : for (size_t i = 0; i < n; i++) {
920 10 : x.set(i, 0, 3.0 * (i + 1.0));
921 65 : for (size_t j = 0; j <= i; j++) {
922 55 : A.set(i, j, 4.0 * i + 7.0 * j + 1.0);
923 : }
924 : }
925 1 : double correct_val = 855904.5f;
926 1 : double val = A.quad(x);
927 1 : _ASSERT_EQ(correct_val, val);
928 :
929 2 : Matrix q(n, 1);
930 11 : for (size_t i = 0; i < n; i++) {
931 10 : q.set(i, 0, -5.0 * i - 1.0);
932 : }
933 :
934 1 : val = A.quad(x, q);
935 1 : correct_val = 850789.5;
936 2 : _ASSERT_EQ(correct_val, val);
937 1 : }
938 :
939 1 : void TestMatrix::testLeftTransposeMultiply() {
940 1 : size_t n = 10;
941 1 : size_t m = 5;
942 1 : size_t k = 8;
943 :
944 1 : Matrix A = MatrixFactory::MakeRandomMatrix(m, n, 00., 1.0, Matrix::MATRIX_DENSE);
945 2 : Matrix B = MatrixFactory::MakeRandomMatrix(m, k, 00., 1.0, Matrix::MATRIX_DENSE);
946 :
947 1 : A.transpose();
948 1 : _ASSERT_EQ(n, A.getNrows());
949 1 : _ASSERT_EQ(m, A.getNcols());
950 :
951 2 : Matrix C = A * B;
952 :
953 1 : _ASSERT_EQ(n, C.getNrows());
954 1 : _ASSERT_EQ(k, C.getNcols());
955 :
956 1 : A.transpose();
957 2 : Matrix AT(n, m, Matrix::MATRIX_DENSE);
958 11 : for (size_t i = 0; i < n; i++) {
959 60 : for (size_t j = 0; j < m; j++) {
960 50 : AT.set(i, j, A.get(j, i));
961 : }
962 : }
963 2 : Matrix C0 = AT * B;
964 :
965 11 : for (size_t i = 0; i < n; i++) {
966 60 : for (size_t j = 0; j < m; j++) {
967 50 : _ASSERT_EQ(C0.get(i, j), C.get(i, j));
968 : }
969 1 : }
970 :
971 1 : }
972 :
973 1 : void TestMatrix::testRightTransposeMultiply() {
974 1 : const size_t n = 10;
975 1 : const size_t m = 5;
976 1 : const size_t k = 8;
977 :
978 1 : Matrix A = MatrixFactory::MakeRandomMatrix(n, m, 00., 1.0, Matrix::MATRIX_DENSE);
979 2 : Matrix B = MatrixFactory::MakeRandomMatrix(k, m, 00., 1.0, Matrix::MATRIX_DENSE);
980 :
981 1 : B.transpose();
982 2 : Matrix C = A*B;
983 :
984 1 : B.transpose();
985 2 : Matrix BT(m, k);
986 6 : for (size_t i = 0; i < m; i++) {
987 45 : for (size_t j = 0; j < k; j++) {
988 40 : BT.set(i, j, B.get(j, i));
989 : }
990 : }
991 :
992 2 : Matrix C_correct = A*BT;
993 1 : _ASSERT_EQ(C_correct.getNrows(), C.getNrows());
994 1 : _ASSERT_EQ(C_correct.getNcols(), C.getNcols());
995 11 : for (size_t i = 0; i < n; i++) {
996 90 : for (size_t j = 0; j < k; j++) {
997 80 : _ASSERT_EQ(C_correct.get(i, j), C.get(i, j));
998 : }
999 1 : }
1000 1 : }
1001 :
1002 1 : void TestMatrix::testLeftSymmetricMultiply() {
1003 1 : const size_t n = 5;
1004 1 : Matrix S = MatrixFactory::MakeRandomMatrix(n, n, 0.0, 1.0, Matrix::MATRIX_SYMMETRIC);
1005 2 : Matrix A = MatrixFactory::MakeRandomMatrix(n, n, 0.0, 1.0, Matrix::MATRIX_DENSE);
1006 2 : Matrix C = S*A;
1007 :
1008 2 : Matrix S2D(n, n, Matrix::MATRIX_DENSE);
1009 6 : for (size_t i = 0; i < n; i++) {
1010 30 : for (size_t j = 0; j < n; j++) {
1011 25 : S2D.set(i, j, S.get(i, j));
1012 : }
1013 : }
1014 :
1015 2 : Matrix C_correct = S2D*A;
1016 6 : for (size_t i = 0; i < n; i++) {
1017 30 : for (size_t j = 0; j < n; j++) {
1018 25 : _ASSERT_EQ(C_correct.get(i, j), C.get(i, j));
1019 : }
1020 : }
1021 :
1022 2 : Matrix x = MatrixFactory::MakeRandomMatrix(n, 1, 0.0, 1.0, Matrix::MATRIX_DENSE);
1023 1 : C = S*x;
1024 :
1025 1 : C_correct = S2D * x;
1026 :
1027 6 : for (size_t i = 0; i < n; i++) {
1028 5 : _ASSERT_NUM_EQ(C_correct.get(i, 0), C.get(i, 0), 1e-4);
1029 1 : }
1030 1 : }
1031 :
1032 1 : void TestMatrix::testSparseGetSet() {
1033 1 : size_t n = 5;
1034 1 : size_t m = 10;
1035 1 : const size_t max_nnz = 3;
1036 1 : Matrix M = MatrixFactory::MakeSparse(n, m, max_nnz, Matrix::SPARSE_UNSYMMETRIC);
1037 :
1038 1 : double r[3] = {4.576, 3.645, 1.092};
1039 1 : M.set(0, 0, r[0]);
1040 1 : M.set(0, 1, r[1]);
1041 1 : M.set(1, 1, r[2]);
1042 :
1043 1 : _ASSERT_EQ(r[0], M.get(0, 0));
1044 1 : _ASSERT_EQ(r[1], M.get(0, 1));
1045 1 : _ASSERT_EQ(r[2], M.get(1, 1));
1046 1 : _ASSERT_EQ(n, M.getNrows());
1047 1 : _ASSERT_EQ(m, M.getNcols());
1048 1 : _ASSERT_EQ(Matrix::MATRIX_SPARSE, M.getType());
1049 1 : _ASSERT_EQ(0, M.cholmod_handle()->status);
1050 :
1051 1 : _ASSERT_OK(Matrix::destroy_handle());
1052 :
1053 2 : Matrix S = MatrixFactory::MakeSparse(n, m, max_nnz, Matrix::SPARSE_UNSYMMETRIC);
1054 22 : for (size_t i = 0; i <= 20; i++) {
1055 21 : _ASSERT_OK(S.set(0, 0, 0.5 * i));
1056 : }
1057 1 : _ASSERT_EQ(10.0, S.get(0, 0));
1058 :
1059 1 : _ASSERT_OK(S.set(1, 0, 1.0));
1060 1 : _ASSERT_EQ(1.0, S.get(1, 0));
1061 :
1062 1 : _ASSERT_OK(S.set(2, 2, 5.0));
1063 2 : _ASSERT_EQ(5.0, S.get(2, 2));
1064 :
1065 :
1066 :
1067 1 : }
1068 :
1069 1 : void TestMatrix::test_MSD() {
1070 :
1071 1 : size_t n = 3;
1072 1 : size_t m = 3;
1073 1 : size_t max_nnz = 4;
1074 :
1075 1 : Matrix A = MatrixFactory::MakeSparse(n, m, max_nnz, Matrix::SPARSE_UNSYMMETRIC);
1076 1 : A.set(0, 0, 4.0);
1077 1 : A.set(1, 0, 1.0); // A is declared as SPARSE_SYMMETRIC - no need to define A(0,1).
1078 1 : A.set(1, 1, 5.0);
1079 1 : A.set(2, 2, 10.0);
1080 :
1081 :
1082 2 : Matrix B(m, n);
1083 4 : for (size_t i = 0; i < m; i++) {
1084 12 : for (size_t j = 0; j < n; j++) {
1085 9 : B.set(i, j, 3.0 * i + 4.23 * j + 1.10);
1086 : }
1087 : }
1088 :
1089 :
1090 1 : _ASSERT_EQ(m, B.getNrows());
1091 1 : _ASSERT_EQ(n, B.getNcols());
1092 :
1093 2 : Matrix C;
1094 1 : _ASSERT_OK(C = A * B);
1095 :
1096 : double correctData[9] = {
1097 : 4.4000, 21.6000, 71.0000,
1098 : 21.3200, 46.9800, 113.3000,
1099 : 38.2400, 72.3600, 155.6000
1100 1 : };
1101 :
1102 2 : Matrix C_correct(n, n, correctData, Matrix::MATRIX_DENSE);
1103 :
1104 :
1105 1 : _ASSERT_EQ(C_correct, C);
1106 2 : _ASSERT_EQ(0, Matrix::cholmod_handle()->status);
1107 :
1108 1 : }
1109 :
1110 1 : void TestMatrix::test_MSDT() {
1111 1 : const size_t n = 13;
1112 1 : const size_t m = 25;
1113 1 : const size_t nnz = 18;
1114 1 : const size_t rep = 50;
1115 :
1116 51 : for (size_t k = 0; k < rep; k++) {
1117 50 : Matrix S = MatrixFactory::MakeRandomSparse(m, n, nnz, 0.2, 2.0);
1118 100 : Matrix D = MatrixFactory::MakeRandomMatrix(m, n, 0.5, 3.0, Matrix::MATRIX_DENSE);
1119 :
1120 100 : Matrix S2D(m, n);
1121 1300 : for (size_t i = 0; i < S.getNrows(); i++) {
1122 17500 : for (size_t j = 0; j < S.getNcols(); j++) {
1123 16250 : S2D.set(i, j, S.get(i, j));
1124 : }
1125 : }
1126 :
1127 50 : D.transpose();
1128 :
1129 100 : Matrix R = S*D; /* R = S * D' */
1130 100 : Matrix correct = S2D*D;
1131 :
1132 50 : const double tol = 1e-8;
1133 1300 : for (size_t i = 0; i < R.getNrows(); i++) {
1134 32500 : for (size_t j = 0; j < R.getNcols(); j++) {
1135 31250 : _ASSERT_NUM_EQ(correct.get(i, j), R.get(i, j), tol);
1136 : }
1137 : }
1138 :
1139 50 : }
1140 1 : }
1141 :
1142 1 : void TestMatrix::test_MSTDT() {
1143 1 : const size_t n = 5;
1144 1 : const size_t m = 8;
1145 1 : const size_t nnz = 23;
1146 1 : Matrix S = MatrixFactory::MakeRandomSparse(m, n, nnz, 0.2, 2.0);
1147 2 : Matrix D = MatrixFactory::MakeRandomMatrix(n, m, 0.5, 3.0, Matrix::MATRIX_DENSE);
1148 :
1149 2 : Matrix DS = D*S;
1150 1 : DS.transpose();
1151 :
1152 1 : S.transpose(); // n-m
1153 1 : D.transpose(); // m-n
1154 :
1155 2 : Matrix R = S * D;
1156 :
1157 : /*
1158 : * S'*D' = (D*S)'
1159 : */
1160 1 : _ASSERT_EQ(DS.getNrows(), R.getNrows());
1161 1 : _ASSERT_EQ(DS.getNcols(), R.getNcols());
1162 6 : for (size_t i = 0; i < R.getNrows(); i++) {
1163 30 : for (size_t j = 0; j < R.getNcols(); j++) {
1164 25 : _ASSERT_NUM_EQ(DS.get(i, j), R.get(i, j), 1e-8);
1165 : }
1166 1 : }
1167 :
1168 1 : }
1169 :
1170 1 : void TestMatrix::test_MDS() {
1171 1 : const size_t n = 5;
1172 1 : const size_t m = 7;
1173 :
1174 1 : Matrix D = MatrixFactory::MakeRandomMatrix(n, m, 0.5, 3.0, Matrix::MATRIX_DENSE);
1175 2 : Matrix S = MatrixFactory::MakeRandomSparse(m, m, static_cast<size_t> (2.5 * (n + m)), 0.2, 2.0);
1176 :
1177 2 : Matrix R = D*S;
1178 :
1179 1 : D.transpose();
1180 1 : S.transpose();
1181 :
1182 2 : Matrix C = S * D;
1183 :
1184 1 : C.transpose();
1185 :
1186 1 : _ASSERT_EQ(C.getNrows(), R.getNrows());
1187 1 : _ASSERT_EQ(C.getNrows(), n);
1188 1 : _ASSERT_EQ(C.getNcols(), R.getNcols());
1189 1 : _ASSERT_EQ(C.getNcols(), m);
1190 6 : for (size_t i = 0; i < C.getNrows(); i++) {
1191 40 : for (size_t j = 0; j < C.getNcols(); j++) {
1192 35 : _ASSERT_NUM_EQ(C.get(i, j), R.get(i, j), 1e-8);
1193 : }
1194 1 : }
1195 :
1196 1 : }
1197 :
1198 1 : void TestMatrix::test_MSS() {
1199 1 : const size_t n = 10;
1200 1 : const size_t m = 12;
1201 1 : const size_t k = 9;
1202 1 : const size_t nnz_L = 50;
1203 1 : const size_t nnz_R = 40;
1204 :
1205 1 : Matrix L = MatrixFactory::MakeRandomSparse(n, m, nnz_L, 1.0, 2.0);
1206 2 : Matrix R = MatrixFactory::MakeRandomSparse(m, k, nnz_R, 1.0, 2.0);
1207 :
1208 2 : Matrix Y;
1209 1 : _ASSERT_OK(Y = L * R);
1210 1 : _ASSERT_EQ(n, Y.getNrows());
1211 1 : _ASSERT_EQ(k, Y.getNcols());
1212 1 : _ASSERT_EQ(Matrix::MATRIX_SPARSE, Y.getType());
1213 :
1214 2 : Matrix Ld(n, m);
1215 2 : Matrix Rd(m, k);
1216 :
1217 11 : for (size_t i = 0; i < n; i++) {
1218 130 : for (size_t j = 0; j < m; j++) {
1219 120 : Ld.set(i, j, L.get(i, j));
1220 : }
1221 : }
1222 :
1223 13 : for (size_t i = 0; i < m; i++) {
1224 120 : for (size_t j = 0; j < k; j++) {
1225 108 : Rd.set(i, j, R.get(i, j));
1226 : }
1227 : }
1228 :
1229 2 : Matrix Z = Ld*Rd;
1230 11 : for (size_t i = 0; i < n; i++) {
1231 100 : for (size_t j = 0; j < k; j++) {
1232 90 : _ASSERT_NUM_EQ(Z.get(i, j), Y.get(i, j), 1e-6);
1233 : }
1234 1 : }
1235 1 : }
1236 :
1237 1 : void TestMatrix::testSparseAddDense() {
1238 1 : const size_t n = 5;
1239 1 : const size_t m = 7;
1240 1 : const size_t nnz = 6;
1241 1 : Matrix A = MatrixFactory::MakeSparse(n, m, nnz, Matrix::SPARSE_UNSYMMETRIC);
1242 1 : A.set(0, 0, 0.5);
1243 1 : A.set(0, 1, 0.2);
1244 1 : A.set(1, 0, 1.2);
1245 1 : A.set(1, 1, 3.6);
1246 1 : A.set(4, 5, 6.2);
1247 1 : A.set(3, 6, 9.9);
1248 :
1249 2 : Matrix B(n, m, Matrix::MATRIX_DENSE);
1250 6 : for (size_t i = 0; i < n; i++) {
1251 40 : for (size_t j = 0; j < m; j++) {
1252 35 : B.set(i, j, -i - 3 * j - 1.5f);
1253 : }
1254 : }
1255 :
1256 2 : Matrix A_init(A);
1257 1 : A += B;
1258 :
1259 6 : for (size_t i = 0; i < n; i++) {
1260 40 : for (size_t j = 0; j < m; j++) {
1261 35 : _ASSERT_EQ(A.get(i, j), A_init.get(i, j) + B.get(i, j));
1262 : }
1263 1 : }
1264 :
1265 1 : }
1266 :
1267 1 : void TestMatrix::testSparseAddSparse() {
1268 : FILE *fp_A, *fp_B;
1269 1 : fp_A = fopen("matrices/sparse2.mx", "r");
1270 1 : fp_B = fopen("matrices/sparse3.mx", "r");
1271 :
1272 1 : CPPUNIT_ASSERT_MESSAGE("Can't open sparse2.mx", fp_A != NULL);
1273 1 : CPPUNIT_ASSERT_MESSAGE("Can't open sparse3.mx", fp_B != NULL);
1274 :
1275 1 : Matrix A = MatrixFactory::ReadSparse(fp_A);
1276 2 : Matrix B = MatrixFactory::ReadSparse(fp_B);
1277 :
1278 : /* close file handlers */
1279 1 : _ASSERT_EQ(0, fclose(fp_A));
1280 1 : _ASSERT_EQ(0, fclose(fp_B));
1281 :
1282 2 : Matrix A0 = A;
1283 1 : _ASSERT_OK(A += B);
1284 :
1285 3 : for (size_t i = 0; i < A.getNrows(); i++) {
1286 10 : for (size_t j = 0; j < A.getNcols(); j++) {
1287 8 : _ASSERT_EQ(A0.get(i, j) + B.get(i, j), A.get(i, j));
1288 : }
1289 1 : }
1290 :
1291 1 : }
1292 :
1293 1 : void TestMatrix::testSparseAddSparse2() {
1294 1 : const size_t n = 100;
1295 1 : const size_t m = 250;
1296 1 : const size_t nnz = 200;
1297 1 : Matrix A = MatrixFactory::MakeRandomSparse(n, m, nnz, -1.0, 2.0);
1298 2 : Matrix B = MatrixFactory::MakeRandomSparse(n, m, nnz, -1.0, 2.0);
1299 2 : Matrix C;
1300 1 : _ASSERT_OK(C = A + B);
1301 1 : _ASSERT_EQ(Matrix::MATRIX_SPARSE, C.getType());
1302 1 : _ASSERT_EQ(n, C.getNrows());
1303 1 : _ASSERT_EQ(m, C.getNcols());
1304 :
1305 101 : for (size_t i = 0; i < n; i++) {
1306 25100 : for (size_t j = 0; j < m; j++) {
1307 25000 : _ASSERT_EQ(A.get(i, j) + B.get(i, j), C.get(i, j));
1308 : }
1309 : }
1310 :
1311 1 : _ASSERT_OK(C += A);
1312 :
1313 101 : for (size_t i = 0; i < n; i++) {
1314 25100 : for (size_t j = 0; j < m; j++) {
1315 25000 : _ASSERT_EQ(2 * A.get(i, j) + B.get(i, j), C.get(i, j));
1316 : }
1317 : }
1318 :
1319 2 : _ASSERT_OK(Matrix::destroy_handle());
1320 1 : }
1321 :
1322 1 : void TestMatrix::testSparseQuad() {
1323 1 : size_t n = 60;
1324 1 : size_t nnz = std::floor(1.2f * n);
1325 1 : size_t tests = 20;
1326 :
1327 1 : Matrix *A = new Matrix();
1328 1 : Matrix *x = new Matrix();
1329 :
1330 21 : for (size_t k = 0; k < tests; k++) {
1331 20 : _ASSERT_OK(*A = MatrixFactory::MakeRandomSparse(n, n, nnz, 0.0, 10.0));
1332 20 : _ASSERT_OK(*x = MatrixFactory::MakeRandomMatrix(n, 1, 1.0, 2.0, Matrix::MATRIX_DENSE));
1333 :
1334 20 : double r = A->quad(*x);
1335 :
1336 20 : double r_exp = 0.0;
1337 1220 : for (size_t i = 0; i < n; i++) {
1338 73200 : for (size_t j = 0; j < n; j++) {
1339 72000 : r_exp += x->get(i, 0) * x->get(j, 0) * A->get(i, j);
1340 : }
1341 : }
1342 20 : r *= 2;
1343 20 : double tol = 1e-6;
1344 20 : _ASSERT(std::abs(r) > tol);
1345 20 : _ASSERT_NUM_EQ(r_exp, r, tol);
1346 20 : _ASSERT_EQ(0, Matrix::cholmod_handle()->status);
1347 : }
1348 :
1349 1 : delete A;
1350 1 : delete x;
1351 :
1352 1 : _ASSERT_EQ(0, Matrix::cholmod_handle()->status);
1353 1 : _ASSERT_OK(Matrix::destroy_handle());
1354 1 : }
1355 :
1356 1 : void TestMatrix::testSparseQuadSparseX() {
1357 1 : const double tol = 1e-6;
1358 1 : const size_t runs = 10;
1359 11 : for (size_t p = 0; p < runs; p++) {
1360 240 : for (size_t n = 2; n < 70; n += 3) {
1361 230 : size_t nnz_A = static_cast<size_t> (std::ceil(1.2 * n));
1362 230 : size_t nnz_x = std::max(static_cast<size_t> (1), static_cast<size_t> (std::ceil(0.75 * n)));
1363 230 : Matrix A = MatrixFactory::MakeRandomSparse(n, n, nnz_A, 0.0, 10.0);
1364 460 : Matrix x = MatrixFactory::MakeRandomSparse(n, 1, nnz_x, 0.0, 10.0);
1365 230 : double r, r_exp = 0.0;
1366 230 : _ASSERT_OK(r = A.quad(x));
1367 8280 : for (size_t i = 0; i < n; i++) {
1368 380880 : for (size_t j = 0; j < n; j++) {
1369 372830 : r_exp += x.get(i, 0) * x.get(j, 0) * A.get(i, j) / 2;
1370 : }
1371 : }
1372 230 : _ASSERT_NUM_EQ(r_exp, r, tol);
1373 230 : }
1374 : }
1375 1 : }
1376 :
1377 1 : void TestMatrix::testSparseQuad_q() {
1378 1 : size_t n = 6;
1379 1 : size_t nnz_A = 6;
1380 1 : size_t nnz_x = 4;
1381 1 : size_t nnz_q = 1;
1382 :
1383 1 : Matrix A = MatrixFactory::MakeSparse(n, n, nnz_A, Matrix::SPARSE_UNSYMMETRIC);
1384 1 : A.set(0, 0, 4);
1385 1 : A.set(0, 1, 5);
1386 1 : A.set(1, 0, 8);
1387 1 : A.set(1, 2, 10);
1388 1 : A.set(5, 3, 100);
1389 :
1390 2 : Matrix x = MatrixFactory::MakeSparse(n, 1, nnz_x, Matrix::SPARSE_UNSYMMETRIC);
1391 1 : x.set(0, 0, 1);
1392 1 : x.set(1, 0, 2);
1393 1 : x.set(5, 0, 9);
1394 1 : x.set(4, 0, 3);
1395 :
1396 2 : Matrix q = MatrixFactory::MakeSparse(n, 1, nnz_q, Matrix::SPARSE_UNSYMMETRIC);
1397 1 : q.set(4, 0, 1);
1398 :
1399 1 : double r = A.quad(x, q);
1400 1 : const double r_exp = 18.0;
1401 1 : const double tol = 1e-10;
1402 :
1403 1 : _ASSERT_NUM_EQ(r_exp, r, tol);
1404 1 : _ASSERT_EQ(0, Matrix::cholmod_handle()->status);
1405 2 : _ASSERT_OK(Matrix::destroy_handle());
1406 :
1407 1 : }
1408 :
1409 1 : void TestMatrix::testSparseDotProd() {
1410 1 : const size_t n = 10;
1411 1 : Matrix x = MatrixFactory::MakeSparse(n, 1, 0, Matrix::SPARSE_UNSYMMETRIC);
1412 2 : Matrix result;
1413 1 : _ASSERT_OK(result = x * x);
1414 1 : _ASSERT_NOT(result.isEmpty());
1415 1 : _ASSERT_EQ(static_cast<size_t> (1), result.getNcols());
1416 1 : _ASSERT_EQ(static_cast<size_t> (1), result.getNrows());
1417 1 : _ASSERT_EQ(0.0, result.get(0, 0));
1418 :
1419 1 : x = MatrixFactory::MakeSparse(n, 1, 1, Matrix::SPARSE_UNSYMMETRIC);
1420 1 : x.set(0, 0, 2.0);
1421 1 : result = x*x;
1422 2 : _ASSERT_EQ(4.0, result.get(0, 0));
1423 1 : }
1424 :
1425 : /* Tests R = DENSE + (?) */
1426 :
1427 1 : void TestMatrix::test_ADD2() {
1428 1 : size_t n = 80;
1429 1 : size_t m = 5;
1430 1 : const double tol = 1e-10;
1431 1 : Matrix D1 = MatrixFactory::MakeRandomMatrix(n, m, -5.0, 6.0, Matrix::MATRIX_DENSE);
1432 2 : Matrix D2 = MatrixFactory::MakeRandomMatrix(n, m, -5.0, 6.0, Matrix::MATRIX_DENSE);
1433 :
1434 2 : Matrix R;
1435 1 : _ASSERT_OK(R = D1 + D2);
1436 1 : _ASSERT_EQ(Matrix::MATRIX_DENSE, R.getType());
1437 81 : for (size_t i = 0; i < n; i++) {
1438 480 : for (size_t j = 0; j < m; j++) {
1439 400 : _ASSERT_NUM_EQ(D1.get(i, j) + D2.get(i, j), R.get(i, j), tol);
1440 : }
1441 1 : }
1442 1 : }
1443 :
1444 1 : void TestMatrix::test_ADH() {
1445 1 : size_t n = 100;
1446 1 : const double tol = 1e-10;
1447 1 : Matrix D = MatrixFactory::MakeRandomMatrix(n, n, -5.0, 6.0, Matrix::MATRIX_DENSE);
1448 2 : Matrix H = MatrixFactory::MakeRandomMatrix(n, n, -5.0, 6.0, Matrix::MATRIX_SYMMETRIC);
1449 :
1450 2 : Matrix R;
1451 1 : _ASSERT_OK(R = D + H);
1452 1 : _ASSERT_EQ(Matrix::MATRIX_DENSE, R.getType());
1453 101 : for (size_t i = 0; i < n; i++) {
1454 10100 : for (size_t j = 0; j < n; j++) {
1455 10000 : _ASSERT_NUM_EQ(D.get(i, j) + H.get(i, j), R.get(i, j), tol);
1456 : }
1457 1 : }
1458 1 : }
1459 :
1460 1 : void TestMatrix::test_ADL() {
1461 1 : size_t n = 80;
1462 1 : const double tol = 1e-10;
1463 1 : Matrix D = MatrixFactory::MakeRandomMatrix(n, n, -5.0, 6.0, Matrix::MATRIX_DENSE);
1464 2 : Matrix L = MatrixFactory::MakeRandomMatrix(n, n, -5.0, 6.0, Matrix::MATRIX_LOWERTR);
1465 :
1466 2 : Matrix R;
1467 1 : _ASSERT_OK(R = D + L);
1468 1 : _ASSERT_EQ(Matrix::MATRIX_DENSE, R.getType());
1469 81 : for (size_t i = 0; i < n; i++) {
1470 6480 : for (size_t j = 0; j < n; j++) {
1471 6400 : _ASSERT_NUM_EQ(D.get(i, j) + L.get(i, j), R.get(i, j), tol);
1472 6400 : if (i < j) {
1473 3160 : _ASSERT_NUM_EQ(R.get(i, j), D.get(i, j), tol);
1474 : }
1475 : }
1476 1 : }
1477 :
1478 :
1479 1 : }
1480 :
1481 1 : void TestMatrix::test_ADS() {
1482 1 : size_t n = 50;
1483 1 : size_t m = 120;
1484 1 : size_t nnz = 200;
1485 1 : const double tol = 1e-10;
1486 1 : Matrix D = MatrixFactory::MakeRandomMatrix(n, m, -5.0, 6.0, Matrix::MATRIX_DENSE);
1487 2 : Matrix S = MatrixFactory::MakeRandomSparse(n, m, nnz, 0.0, 1.0);
1488 :
1489 2 : Matrix R;
1490 1 : _ASSERT_OK(R = D + S);
1491 1 : _ASSERT_EQ(Matrix::MATRIX_DENSE, R.getType());
1492 51 : for (size_t i = 0; i < n; i++) {
1493 6050 : for (size_t j = 0; j < m; j++) {
1494 6000 : _ASSERT_NUM_EQ(D.get(i, j) + S.get(i, j), R.get(i, j), tol);
1495 : }
1496 1 : }
1497 1 : }
1498 :
1499 1 : void TestMatrix::test_ADW() {
1500 : //std::cerr << "test_ADW:::EMPTY TEST!!!!\n";
1501 1 : }
1502 :
1503 1 : void TestMatrix::test_ADX() {
1504 1 : size_t n = 100;
1505 1 : const double tol = 1e-10;
1506 1 : Matrix D = MatrixFactory::MakeRandomMatrix(n, n, -5.0, 6.0, Matrix::MATRIX_DENSE);
1507 2 : Matrix X = MatrixFactory::MakeRandomMatrix(n, n, 1.0, 10.0, Matrix::MATRIX_DIAGONAL);
1508 2 : Matrix R;
1509 1 : _ASSERT_OK(R = D + X);
1510 1 : _ASSERT_EQ(Matrix::MATRIX_DENSE, R.getType());
1511 101 : for (size_t i = 0; i < n; i++) {
1512 10100 : for (size_t j = 0; j < n; j++) {
1513 10000 : _ASSERT_NUM_EQ(D.get(i, j) + X.get(i, j), R.get(i, j), tol);
1514 : }
1515 1 : }
1516 1 : }
1517 :
1518 1 : void TestMatrix::test_EH() {
1519 1 : size_t n = 50;
1520 1 : Matrix H = MatrixFactory::MakeRandomMatrix(n, n, -5.0, 6.0, Matrix::MATRIX_SYMMETRIC);
1521 2 : Matrix H2(H);
1522 2 : Matrix H3 = H2;
1523 1 : _ASSERT(H2.getNrows() == n);
1524 1 : _ASSERT(H2.getNcols() == n);
1525 1 : _ASSERT(H2.getType() == Matrix::MATRIX_SYMMETRIC);
1526 1 : _ASSERT(H3.getNrows() == n);
1527 1 : _ASSERT(H3.getNcols() == n);
1528 1 : _ASSERT(H3.getType() == Matrix::MATRIX_SYMMETRIC);
1529 1 : _ASSERT(H2.isSymmetric());
1530 1 : _ASSERT(H3.isSymmetric());
1531 51 : for (size_t i = 0; i < n; i++) {
1532 2550 : for (size_t j = 0; j < n; j++) {
1533 2500 : _ASSERT_EQ(H.get(i, j), H2.get(i, j));
1534 2500 : _ASSERT_EQ(H.get(i, j), H3.get(i, j));
1535 : }
1536 1 : }
1537 1 : }
1538 :
1539 1 : void TestMatrix::test_EL() {
1540 1 : size_t n = 50;
1541 1 : Matrix L = MatrixFactory::MakeRandomMatrix(n, n, -5.0, 6.0, Matrix::MATRIX_LOWERTR);
1542 2 : Matrix L2(L);
1543 2 : Matrix L3 = L2;
1544 1 : _ASSERT(L2.getNrows() == n);
1545 1 : _ASSERT(L2.getNcols() == n);
1546 1 : _ASSERT(L2.getType() == Matrix::MATRIX_LOWERTR);
1547 1 : _ASSERT(L3.getNrows() == n);
1548 1 : _ASSERT(L3.getNcols() == n);
1549 1 : _ASSERT(L3.getType() == Matrix::MATRIX_LOWERTR);
1550 51 : for (size_t i = 0; i < n; i++) {
1551 2550 : for (size_t j = 0; j < n; j++) {
1552 2500 : _ASSERT_EQ(L.get(i, j), L2.get(i, j));
1553 2500 : _ASSERT_EQ(L.get(i, j), L3.get(i, j));
1554 : }
1555 1 : }
1556 1 : }
1557 :
1558 1 : void TestMatrix::test_EX() {
1559 1 : size_t n = 50;
1560 1 : Matrix X = MatrixFactory::MakeRandomMatrix(n, n, -5.0, 6.0, Matrix::MATRIX_DIAGONAL);
1561 2 : Matrix X2(X);
1562 2 : Matrix X3 = X2;
1563 1 : _ASSERT(X2.getNrows() == n);
1564 1 : _ASSERT(X2.getNcols() == n);
1565 1 : _ASSERT(X2.getType() == Matrix::MATRIX_DIAGONAL);
1566 1 : _ASSERT(X3.getNrows() == n);
1567 1 : _ASSERT(X3.getNcols() == n);
1568 1 : _ASSERT(X3.getType() == Matrix::MATRIX_DIAGONAL);
1569 51 : for (size_t i = 0; i < n; i++) {
1570 2550 : for (size_t j = 0; j < n; j++) {
1571 2500 : _ASSERT_EQ(X.get(i, j), X2.get(i, j));
1572 2500 : _ASSERT_EQ(X.get(i, j), X3.get(i, j));
1573 2500 : if (i != j) {
1574 2450 : _ASSERT_EQ(0.0, X.get(i, j));
1575 : }
1576 : }
1577 1 : }
1578 1 : }
1579 :
1580 1 : void TestMatrix::test_EHT() {
1581 1 : size_t n = 50;
1582 1 : Matrix H = MatrixFactory::MakeRandomMatrix(n, n, -5.0, 6.0, Matrix::MATRIX_SYMMETRIC);
1583 2 : Matrix H2(H);
1584 1 : H2.transpose();
1585 2 : Matrix H3 = H2;
1586 1 : H3.transpose();
1587 1 : _ASSERT(H2.getNrows() == n);
1588 1 : _ASSERT(H2.getNcols() == n);
1589 1 : _ASSERT(H2.getType() == Matrix::MATRIX_SYMMETRIC);
1590 1 : _ASSERT(H3.getNrows() == n);
1591 1 : _ASSERT(H3.getNcols() == n);
1592 1 : _ASSERT(H3.getType() == Matrix::MATRIX_SYMMETRIC);
1593 1 : _ASSERT(H2.isSymmetric());
1594 1 : _ASSERT(H3.isSymmetric());
1595 51 : for (size_t i = 0; i < n; i++) {
1596 2550 : for (size_t j = 0; j < n; j++) {
1597 2500 : _ASSERT_EQ(H.get(i, j), H2.get(i, j));
1598 2500 : _ASSERT_EQ(H.get(i, j), H2.get(j, i));
1599 2500 : _ASSERT_EQ(H.get(i, j), H3.get(i, j));
1600 2500 : _ASSERT_EQ(H.get(i, j), H3.get(j, i));
1601 : }
1602 1 : }
1603 1 : }
1604 :
1605 1 : void TestMatrix::test_EDT() {
1606 1 : size_t n = 100;
1607 1 : Matrix D = MatrixFactory::MakeRandomMatrix(n, n, -5.0, 6.0, Matrix::MATRIX_DENSE);
1608 2 : Matrix D1 = D;
1609 2 : Matrix D2(D);
1610 1 : D1.transpose();
1611 1 : D2.transpose();
1612 101 : for (size_t i = 0; i < n; i++) {
1613 10100 : for (size_t j = 0; j < n; j++) {
1614 10000 : _ASSERT_EQ(D.get(i, j), D1.get(j, i));
1615 10000 : _ASSERT_EQ(D.get(i, j), D2.get(j, i));
1616 : }
1617 1 : }
1618 1 : }
1619 :
1620 1 : void TestMatrix::test_ADDT() {
1621 1 : size_t n = 10;
1622 1 : const double tol = 1e-10;
1623 1 : Matrix D1 = MatrixFactory::MakeRandomMatrix(n, n, -5.0, 6.0, Matrix::MATRIX_DENSE);
1624 2 : Matrix D2 = MatrixFactory::MakeRandomMatrix(n, n, -5.0, 6.0, Matrix::MATRIX_DENSE);
1625 2 : Matrix D2t(D2); // D2t = D2
1626 1 : D2t.transpose(); // D2t = D2t'
1627 :
1628 11 : for (size_t i = 0; i < n; i++) {
1629 110 : for (size_t j = 0; j < n; j++) {
1630 100 : _ASSERT_NUM_EQ(D2.get(i, j), D2t.get(j, i), tol);
1631 : }
1632 : }
1633 :
1634 2 : Matrix R;
1635 1 : _ASSERT_OK(R = D1 + D2t);
1636 1 : _ASSERT_EQ(Matrix::MATRIX_DENSE, R.getType());
1637 :
1638 11 : for (size_t i = 0; i < n; i++) {
1639 110 : for (size_t j = 0; j < n; j++) {
1640 100 : _ASSERT_NUM_EQ(D1.get(i, j) + D2.get(j, i), R.get(i, j), tol);
1641 : }
1642 1 : }
1643 1 : }
1644 :
1645 1 : void TestMatrix::test_ADHT() {
1646 1 : size_t n = 10;
1647 1 : const double tol = 1e-10;
1648 1 : Matrix D = MatrixFactory::MakeRandomMatrix(n, n, -5.0, 6.0, Matrix::MATRIX_DENSE);
1649 2 : Matrix H = MatrixFactory::MakeRandomMatrix(n, n, -5.0, 6.0, Matrix::MATRIX_SYMMETRIC);
1650 2 : Matrix Ht(H); // D2t = D2
1651 1 : Ht.transpose(); // D2t = D2t'
1652 :
1653 11 : for (size_t i = 0; i < n; i++) {
1654 110 : for (size_t j = 0; j < n; j++) {
1655 100 : _ASSERT_NUM_EQ(H.get(i, j), Ht.get(j, i), tol);
1656 : }
1657 : }
1658 :
1659 2 : Matrix R;
1660 1 : _ASSERT_OK(R = D + Ht);
1661 1 : _ASSERT_EQ(Matrix::MATRIX_DENSE, R.getType());
1662 :
1663 11 : for (size_t i = 0; i < n; i++) {
1664 110 : for (size_t j = 0; j < n; j++) {
1665 100 : _ASSERT_NUM_EQ(D.get(i, j) + H.get(j, i), R.get(i, j), tol);
1666 : }
1667 1 : }
1668 1 : }
1669 :
1670 1 : void TestMatrix::test_ADLT() {
1671 1 : size_t n = 10;
1672 1 : const double tol = 1e-10;
1673 1 : Matrix D = MatrixFactory::MakeRandomMatrix(n, n, -5.0, 6.0, Matrix::MATRIX_DENSE);
1674 2 : Matrix L = MatrixFactory::MakeRandomMatrix(n, n, -5.0, 6.0, Matrix::MATRIX_LOWERTR);
1675 2 : Matrix Lt(L); // D2t = D2
1676 1 : Lt.transpose(); // D2t = D2t'
1677 :
1678 11 : for (size_t i = 0; i < n; i++) {
1679 110 : for (size_t j = 0; j < n; j++) {
1680 100 : _ASSERT_NUM_EQ(L.get(i, j), Lt.get(j, i), tol);
1681 : }
1682 : }
1683 :
1684 2 : Matrix R;
1685 1 : _ASSERT_OK(R = D + Lt);
1686 1 : _ASSERT_EQ(Matrix::MATRIX_DENSE, R.getType());
1687 :
1688 11 : for (size_t i = 0; i < n; i++) {
1689 110 : for (size_t j = 0; j < n; j++) {
1690 100 : _ASSERT_NUM_EQ(D.get(i, j) + L.get(j, i), R.get(i, j), tol);
1691 : }
1692 1 : }
1693 1 : }
1694 :
1695 1 : void TestMatrix::test_ADST() {
1696 1 : size_t n = 50;
1697 1 : size_t nnz = 200;
1698 1 : const double tol = 1e-10;
1699 1 : Matrix D = MatrixFactory::MakeRandomMatrix(n, n, -5.0, 6.0, Matrix::MATRIX_DENSE);
1700 2 : Matrix S = MatrixFactory::MakeRandomSparse(n, n, nnz, 0.0, 10.0);
1701 2 : Matrix St = S; // D2t = D2
1702 1 : St.transpose(); // D2t = D2t'
1703 :
1704 51 : for (size_t i = 0; i < n; i++) {
1705 2550 : for (size_t j = 0; j < n; j++) {
1706 2500 : _ASSERT_NUM_EQ(S.get(i, j), St.get(j, i), tol);
1707 : }
1708 : }
1709 :
1710 2 : Matrix R;
1711 1 : _ASSERT_OK(R = D + St);
1712 1 : _ASSERT_EQ(Matrix::MATRIX_DENSE, R.getType());
1713 :
1714 51 : for (size_t i = 0; i < n; i++) {
1715 2550 : for (size_t j = 0; j < n; j++) {
1716 2500 : _ASSERT_NUM_EQ(D.get(i, j) + S.get(j, i), R.get(i, j), tol);
1717 : }
1718 1 : }
1719 1 : }
1720 :
1721 1 : void TestMatrix::test_ADWT() {
1722 :
1723 1 : }
1724 :
1725 1 : void TestMatrix::test_ADXT() {
1726 1 : size_t n = 10;
1727 1 : const double tol = 1e-10;
1728 1 : Matrix D = MatrixFactory::MakeRandomMatrix(n, n, -5.0, 6.0, Matrix::MATRIX_DENSE);
1729 2 : Matrix X = MatrixFactory::MakeRandomMatrix(n, n, -5.0, 6.0, Matrix::MATRIX_DIAGONAL);
1730 2 : Matrix Xt(X); // D2t = D2
1731 1 : Xt.transpose(); // D2t = D2t'
1732 :
1733 11 : for (size_t i = 0; i < n; i++) {
1734 110 : for (size_t j = 0; j < n; j++) {
1735 100 : _ASSERT_NUM_EQ(X.get(i, j), Xt.get(j, i), tol);
1736 : }
1737 : }
1738 :
1739 2 : Matrix R;
1740 1 : _ASSERT_OK(R = D + Xt);
1741 1 : _ASSERT_EQ(Matrix::MATRIX_DENSE, R.getType());
1742 :
1743 11 : for (size_t i = 0; i < n; i++) {
1744 110 : for (size_t j = 0; j < n; j++) {
1745 100 : _ASSERT_NUM_EQ(D.get(i, j) + X.get(j, i), R.get(i, j), tol);
1746 : }
1747 1 : }
1748 1 : }
1749 :
1750 1 : void TestMatrix::test_AHH() {
1751 1 : size_t n = 100;
1752 1 : const double tol = 1e-10;
1753 1 : Matrix H1 = MatrixFactory::MakeRandomMatrix(n, n, -5.0, 6.0, Matrix::MATRIX_SYMMETRIC);
1754 2 : Matrix H2 = MatrixFactory::MakeRandomMatrix(n, n, -5.0, 6.0, Matrix::MATRIX_SYMMETRIC);
1755 :
1756 2 : Matrix R;
1757 1 : _ASSERT_OK(R = H1 + H2);
1758 1 : _ASSERT_EQ(Matrix::MATRIX_SYMMETRIC, R.getType());
1759 101 : for (size_t i = 0; i < n; i++) {
1760 10100 : for (size_t j = 0; j < n; j++) {
1761 10000 : _ASSERT_NUM_EQ(H1.get(i, j) + H2.get(i, j), R.get(i, j), tol);
1762 : }
1763 1 : }
1764 1 : }
1765 :
1766 1 : void TestMatrix::test_AHX() {
1767 1 : size_t n = 100;
1768 1 : const double tol = 1e-10;
1769 1 : Matrix H = MatrixFactory::MakeRandomMatrix(n, n, -5.0, 6.0, Matrix::MATRIX_SYMMETRIC);
1770 2 : Matrix X = MatrixFactory::MakeRandomMatrix(n, n, -5.0, 6.0, Matrix::MATRIX_DIAGONAL);
1771 :
1772 2 : Matrix R;
1773 1 : _ASSERT_EQ(Matrix::MATRIX_SYMMETRIC, H.getType());
1774 1 : _ASSERT_EQ(Matrix::MATRIX_DIAGONAL, X.getType());
1775 1 : _ASSERT_OK(R = H + X);
1776 1 : _ASSERT_EQ(n, R.getNrows());
1777 1 : _ASSERT_EQ(n, R.getNcols());
1778 1 : _ASSERT(R.isSymmetric());
1779 1 : _ASSERT(H.isSymmetric());
1780 1 : _ASSERT(X.isSymmetric());
1781 1 : _ASSERT_EQ(Matrix::MATRIX_SYMMETRIC, R.getType());
1782 101 : for (size_t i = 0; i < n; i++) {
1783 10100 : for (size_t j = 0; j < n; j++) {
1784 10000 : _ASSERT_NUM_EQ(H.get(i, j) + X.get(i, j), R.get(i, j), tol);
1785 : }
1786 1 : }
1787 1 : }
1788 :
1789 1 : void TestMatrix::test_AHD() {
1790 1 : size_t n = 100;
1791 1 : const double tol = 1e-10;
1792 1 : Matrix H = MatrixFactory::MakeRandomMatrix(n, n, -5.0, 6.0, Matrix::MATRIX_SYMMETRIC);
1793 2 : Matrix D = MatrixFactory::MakeRandomMatrix(n, n, -5.0, 6.0, Matrix::MATRIX_DENSE);
1794 :
1795 2 : Matrix R;
1796 1 : _ASSERT_EQ(Matrix::MATRIX_SYMMETRIC, H.getType());
1797 1 : R = H + D;
1798 1 : _ASSERT_EQ(Matrix::MATRIX_DENSE, R.getType()); /* H + D = D */
1799 1 : _ASSERT_EQ(n, R.getNrows());
1800 1 : _ASSERT_EQ(n, R.getNcols());
1801 1 : _ASSERT_NOT(R.isSymmetric());
1802 101 : for (size_t i = 0; i < n; i++) {
1803 10100 : for (size_t j = 0; j < n; j++) {
1804 10000 : _ASSERT_NUM_EQ(H.get(i, j) + D.get(i, j), R.get(i, j), tol);
1805 : }
1806 1 : }
1807 1 : }
1808 :
1809 1 : void TestMatrix::test_AHL() {
1810 1 : size_t n = 100;
1811 1 : const double tol = 1e-10;
1812 1 : Matrix H = MatrixFactory::MakeRandomMatrix(n, n, -5.0, 6.0, Matrix::MATRIX_SYMMETRIC);
1813 2 : Matrix L = MatrixFactory::MakeRandomMatrix(n, n, -5.0, 6.0, Matrix::MATRIX_LOWERTR);
1814 :
1815 2 : Matrix R;
1816 1 : _ASSERT_EQ(Matrix::MATRIX_SYMMETRIC, H.getType());
1817 1 : _ASSERT_OK(R = H + L);
1818 1 : _ASSERT_EQ(Matrix::MATRIX_DENSE, R.getType()); /* H + L = D */
1819 1 : _ASSERT_EQ(n, R.getNrows());
1820 1 : _ASSERT_EQ(n, R.getNcols());
1821 1 : _ASSERT_NOT(R.isSymmetric());
1822 1 : _ASSERT(H.isSymmetric());
1823 1 : _ASSERT_NOT(L.isSymmetric());
1824 101 : for (size_t i = 0; i < n; i++) {
1825 10100 : for (size_t j = 0; j < n; j++) {
1826 10000 : _ASSERT_NUM_EQ(H.get(i, j) + L.get(i, j), R.get(i, j), tol);
1827 : }
1828 1 : }
1829 1 : }
1830 :
1831 1 : void TestMatrix::test_AHS() {
1832 1 : size_t n = 40;
1833 1 : size_t nnz = 90;
1834 1 : const double tol = 1e-10;
1835 1 : Matrix H = MatrixFactory::MakeRandomMatrix(n, n, -5.0, 6.0, Matrix::MATRIX_SYMMETRIC);
1836 2 : Matrix S = MatrixFactory::MakeRandomSparse(n, n, nnz, 0.0, 1.0);
1837 :
1838 2 : Matrix R;
1839 1 : _ASSERT_EQ(Matrix::MATRIX_SYMMETRIC, H.getType());
1840 1 : _ASSERT_OK(R = H + S);
1841 1 : _ASSERT_EQ(Matrix::MATRIX_DENSE, R.getType()); /* H + S = D */
1842 1 : _ASSERT_EQ(n, R.getNrows());
1843 1 : _ASSERT_EQ(n, R.getNcols());
1844 1 : _ASSERT_NOT(R.isSymmetric());
1845 1 : _ASSERT(H.isSymmetric());
1846 1 : _ASSERT_NOT(S.isSymmetric());
1847 41 : for (size_t i = 0; i < n; i++) {
1848 1640 : for (size_t j = 0; j < n; j++) {
1849 1600 : _ASSERT_NUM_EQ(H.get(i, j) + S.get(i, j), R.get(i, j), tol);
1850 : }
1851 1 : }
1852 1 : }
1853 :
1854 1 : void TestMatrix::test_ASD() { /* Sparse + Diagonal */
1855 1 : size_t n = 120;
1856 1 : size_t nnz = 60;
1857 1 : const double tol = 1e-10;
1858 1 : Matrix S = MatrixFactory::MakeRandomSparse(n, n, nnz, 0.0, 1.0);
1859 2 : Matrix D = MatrixFactory::MakeRandomMatrix(n, n, 2.0, 10.0, Matrix::MATRIX_DENSE);
1860 :
1861 2 : Matrix R;
1862 1 : _ASSERT_OK(R = S + D);
1863 1 : _ASSERT_EQ(Matrix::MATRIX_DENSE, R.getType());
1864 1 : _ASSERT_EQ(n, R.getNrows());
1865 1 : _ASSERT_EQ(n, R.getNcols());
1866 121 : for (size_t i = 0; i < n; i++) {
1867 14520 : for (size_t j = 0; j < n; j++) {
1868 14400 : _ASSERT_NUM_EQ(S.get(i, j) + D.get(i, j), R.get(i, j), tol);
1869 : }
1870 1 : }
1871 1 : }
1872 :
1873 1 : void TestMatrix::test_ASH() {
1874 1 : size_t n = 120;
1875 1 : size_t nnz = 60;
1876 1 : const double tol = 1e-10;
1877 1 : Matrix S = MatrixFactory::MakeRandomSparse(n, n, nnz, 0.0, 1.0);
1878 2 : Matrix H = MatrixFactory::MakeRandomMatrix(n, n, 2.0, 10.0, Matrix::MATRIX_SYMMETRIC);
1879 :
1880 2 : Matrix R;
1881 1 : _ASSERT_OK(R = S + H);
1882 1 : _ASSERT_EQ(Matrix::MATRIX_DENSE, R.getType());
1883 1 : _ASSERT_EQ(n, R.getNrows());
1884 1 : _ASSERT_EQ(n, R.getNcols());
1885 121 : for (size_t i = 0; i < n; i++) {
1886 14520 : for (size_t j = 0; j < n; j++) {
1887 14400 : _ASSERT_NUM_EQ(S.get(i, j) + H.get(i, j), R.get(i, j), tol);
1888 : }
1889 1 : }
1890 1 : }
1891 :
1892 1 : void TestMatrix::test_ASX() {
1893 1 : size_t n = 120;
1894 1 : size_t nnz = 60;
1895 1 : const double tol = 1e-10;
1896 1 : Matrix S = MatrixFactory::MakeRandomSparse(n, n, nnz, 0.0, 1.0);
1897 2 : Matrix X = MatrixFactory::MakeRandomMatrix(n, n, 2.0, 10.0, Matrix::MATRIX_DIAGONAL);
1898 :
1899 2 : Matrix R;
1900 1 : _ASSERT_OK(R = S + X); // SPARSE + DIAGONAL = SPARSE
1901 1 : _ASSERT_EQ(Matrix::MATRIX_SPARSE, R.getType());
1902 1 : _ASSERT_EQ(n, R.getNrows());
1903 1 : _ASSERT_EQ(n, R.getNcols());
1904 121 : for (size_t i = 0; i < n; i++) {
1905 14520 : for (size_t j = 0; j < n; j++) {
1906 14400 : _ASSERT_NUM_EQ(S.get(i, j) + X.get(i, j), R.get(i, j), tol);
1907 : }
1908 1 : }
1909 1 : }
1910 :
1911 1 : void TestMatrix::test_ASL() {
1912 1 : size_t n = 20;
1913 1 : size_t nnz = 60;
1914 1 : const double tol = 1e-10;
1915 1 : Matrix S = MatrixFactory::MakeRandomSparse(n, n, nnz, 0.0, 1.0);
1916 2 : Matrix L = MatrixFactory::MakeRandomMatrix(n, n, 2.0, 10.0, Matrix::MATRIX_LOWERTR);
1917 :
1918 2 : Matrix R;
1919 1 : _ASSERT_OK(R = S + L); // SPARSE + LOWER TRIANGULAR = SPARSE
1920 1 : _ASSERT_EQ(Matrix::MATRIX_SPARSE, R.getType());
1921 1 : _ASSERT_EQ(n, R.getNrows());
1922 1 : _ASSERT_EQ(n, R.getNcols());
1923 21 : for (size_t i = 0; i < n; i++) {
1924 420 : for (size_t j = 0; j < n; j++) {
1925 400 : _ASSERT_NUM_EQ(S.get(i, j) + L.get(i, j), R.get(i, j), tol);
1926 : }
1927 1 : }
1928 1 : }
1929 :
1930 1 : void TestMatrix::test_AXX() {
1931 1 : size_t n = 120;
1932 1 : const double tol = 1e-10;
1933 1 : Matrix X1 = MatrixFactory::MakeRandomMatrix(n, n, 2.0, 10.0, Matrix::MATRIX_DIAGONAL);
1934 2 : Matrix X2 = MatrixFactory::MakeRandomMatrix(n, n, 2.0, 10.0, Matrix::MATRIX_DIAGONAL);
1935 :
1936 2 : Matrix R;
1937 1 : _ASSERT_OK(R = X1 + X2); // DIAGONAL + DIAGONAL = DIAGONAL
1938 1 : _ASSERT_EQ(Matrix::MATRIX_DIAGONAL, R.getType());
1939 1 : _ASSERT_EQ(n, R.getNrows());
1940 1 : _ASSERT_EQ(n, R.getNcols());
1941 121 : for (size_t i = 0; i < n; i++) {
1942 14520 : for (size_t j = 0; j < n; j++) {
1943 14400 : _ASSERT_NUM_EQ(X1.get(i, j) + X2.get(i, j), R.get(i, j), tol);
1944 : }
1945 1 : }
1946 1 : }
1947 :
1948 1 : void TestMatrix::test_ALL() {
1949 1 : size_t n = 120;
1950 1 : const double tol = 1e-10;
1951 1 : Matrix L1 = MatrixFactory::MakeRandomMatrix(n, n, 2.0, 10.0, Matrix::MATRIX_LOWERTR);
1952 2 : Matrix L2 = MatrixFactory::MakeRandomMatrix(n, n, 2.0, 10.0, Matrix::MATRIX_LOWERTR);
1953 2 : Matrix L;
1954 1 : L = L1 + L2;
1955 121 : for (size_t i = 0; i < n; i++) {
1956 14520 : for (size_t j = 0; j < n; j++) {
1957 14400 : _ASSERT_NUM_EQ(L1.get(i, j) + L2.get(i, j), L.get(i, j), tol);
1958 : }
1959 1 : }
1960 1 : }
1961 :
1962 1 : void TestMatrix::test_ALX() {
1963 1 : size_t n = 120;
1964 1 : const double tol = 1e-10;
1965 1 : Matrix L = MatrixFactory::MakeRandomMatrix(n, n, 2.0, 10.0, Matrix::MATRIX_LOWERTR);
1966 2 : Matrix X = MatrixFactory::MakeRandomMatrix(n, n, 2.0, 10.0, Matrix::MATRIX_DIAGONAL);
1967 2 : Matrix R;
1968 1 : R = L + X;
1969 121 : for (size_t i = 0; i < n; i++) {
1970 14520 : for (size_t j = 0; j < n; j++) {
1971 14400 : _ASSERT_NUM_EQ(L.get(i, j) + X.get(i, j), R.get(i, j), tol);
1972 : }
1973 1 : }
1974 1 : }
1975 :
1976 1 : void TestMatrix::test_CD() {
1977 1 : size_t n = 120;
1978 1 : const double tol = 1e-10;
1979 1 : double alpha = 1.55;
1980 1 : Matrix D = MatrixFactory::MakeRandomMatrix(n, n, 2.0, 10.0, Matrix::MATRIX_DENSE);
1981 2 : Matrix DD = D;
1982 1 : D *= alpha;
1983 121 : for (size_t i = 0; i < n; i++) {
1984 14520 : for (size_t j = 0; j < n; j++) {
1985 14400 : _ASSERT_NUM_EQ(alpha * DD.get(i, j), D.get(i, j), tol);
1986 : }
1987 1 : }
1988 :
1989 1 : }
1990 :
1991 1 : void TestMatrix::test_CH() {
1992 1 : size_t n = 120;
1993 1 : const double tol = 1e-10;
1994 1 : double alpha = -2.3456;
1995 1 : Matrix H = MatrixFactory::MakeRandomMatrix(n, n, 2.0, 10.0, Matrix::MATRIX_DENSE);
1996 2 : Matrix HH = H;
1997 1 : H *= alpha;
1998 121 : for (size_t i = 0; i < n; i++) {
1999 14520 : for (size_t j = 0; j < n; j++) {
2000 14400 : _ASSERT_NUM_EQ(alpha * HH.get(i, j), H.get(i, j), tol);
2001 : }
2002 1 : }
2003 1 : }
2004 :
2005 1 : void TestMatrix::test_CS() {
2006 1 : size_t n = 200;
2007 1 : size_t nnz = 60;
2008 1 : const double tol = 1e-10;
2009 1 : double alpha = -2.3456;
2010 1 : Matrix S = MatrixFactory::MakeRandomSparse(n, n, nnz, 0.0, 1.0);
2011 2 : Matrix SS = S;
2012 1 : S *= alpha;
2013 201 : for (size_t i = 0; i < n; i++) {
2014 40200 : for (size_t j = 0; j < n; j++) {
2015 40000 : _ASSERT_NUM_EQ(alpha * SS.get(i, j), S.get(i, j), tol);
2016 : }
2017 1 : }
2018 1 : }
2019 :
2020 1 : void TestMatrix::test_ASST() {
2021 1 : const size_t n = 10;
2022 1 : const size_t m = n + 2;
2023 1 : const size_t nnz = 2 * n + 10;
2024 1 : const double tol = 1e-10;
2025 :
2026 1 : Matrix A = MatrixFactory::MakeRandomSparse(n, n, nnz, 0.0, 1.0);
2027 2 : Matrix At(A);
2028 1 : At.transpose();
2029 :
2030 1 : A += At;
2031 11 : for (size_t i = 0; i < n; i++) {
2032 110 : for (size_t j = 0; j < n; j++) {
2033 100 : _ASSERT_NUM_EQ(A.get(i, j), A.get(j, i), tol);
2034 : }
2035 : }
2036 :
2037 : // and some more testing:
2038 1 : const size_t repetitions = 20;
2039 2 : Matrix B, C, Corig, X;
2040 21 : for (size_t k = 0; k < repetitions; k++) {
2041 20 : B = MatrixFactory::MakeRandomSparse(n, m, nnz, 0.0, 1.0);
2042 20 : C = MatrixFactory::MakeRandomSparse(m, n, nnz, 0.0, 1.0);
2043 20 : Corig = C;
2044 20 : C.transpose(); // C is now n-by-m
2045 20 : _ASSERT_EQ(n, C.getNrows());
2046 20 : _ASSERT_EQ(m, C.getNcols());
2047 :
2048 20 : X = B + C;
2049 :
2050 220 : for (size_t i = 0; i < n; i++) {
2051 2600 : for (size_t j = 0; j < m; j++) {
2052 2400 : _ASSERT_NUM_EQ(X.get(i, j), B.get(i, j) + Corig.get(j, i), tol);
2053 : }
2054 : }
2055 1 : }
2056 :
2057 1 : }
2058 :
2059 1 : void TestMatrix::test_EST() {
2060 1 : const size_t n = 100;
2061 1 : const size_t m = 80;
2062 1 : const size_t nnz = 300;
2063 1 : Matrix S = MatrixFactory::MakeRandomSparse(n, m, nnz, 0.0, 1.0);
2064 2 : Matrix S1 = S;
2065 2 : Matrix S2(S);
2066 1 : S1.transpose();
2067 1 : S2.transpose();
2068 101 : for (size_t i = 0; i < n; i++) {
2069 8100 : for (size_t j = 0; j < m; j++) {
2070 8000 : _ASSERT_EQ(S.get(i, j), S1.get(j, i));
2071 8000 : _ASSERT_EQ(S.get(i, j), S2.get(j, i));
2072 : }
2073 1 : }
2074 1 : }
2075 :
2076 1 : void TestMatrix::test_ES() {
2077 1 : const size_t n = 100;
2078 1 : const size_t m = 80;
2079 1 : const size_t nnz = 300;
2080 1 : Matrix S = MatrixFactory::MakeRandomSparse(n, m, nnz, 0.0, 1.0);
2081 2 : Matrix S1 = S;
2082 2 : Matrix S2(S);
2083 101 : for (size_t i = 0; i < n; i++) {
2084 8100 : for (size_t j = 0; j < m; j++) {
2085 8000 : _ASSERT_EQ(S.get(i, j), S1.get(i, j));
2086 8000 : _ASSERT_EQ(S.get(i, j), S2.get(i, j));
2087 : }
2088 1 : }
2089 1 : }
2090 :
2091 1 : void TestMatrix::test_ASS() {
2092 1 : const double tol = 1e-8;
2093 1 : const size_t n = 100;
2094 1 : const size_t m = 80;
2095 1 : const size_t nnz = 300;
2096 1 : Matrix A = MatrixFactory::MakeRandomSparse(n, m, nnz, 0.0, 1.0);
2097 2 : Matrix B = MatrixFactory::MakeRandomSparse(n, m, nnz, 0.0, 1.0);
2098 2 : Matrix X = A + B;
2099 101 : for (size_t i = 0; i < n; i++) {
2100 8100 : for (size_t j = 0; j < m; j++) {
2101 8000 : _ASSERT_NUM_EQ(A.get(i, j) + B.get(i, j), X.get(i, j), tol);
2102 : }
2103 1 : }
2104 1 : }
2105 :
2106 1 : void TestMatrix::testSubmatrix() {
2107 1 : const double tol = 1e-12;
2108 :
2109 1 : const size_t n = 14;
2110 1 : const size_t m = 11;
2111 :
2112 1 : const size_t row_start = 1;
2113 1 : const size_t row_end = 1;
2114 1 : const size_t col_start = 8;
2115 1 : const size_t col_end = 8;
2116 :
2117 1 : Matrix A = MatrixFactory::MakeRandomMatrix(m, n, 2.0, 10.0, Matrix::MATRIX_DENSE);
2118 2 : Matrix subA = A.submatrixCopy(row_start, row_end, col_start, col_end);
2119 1 : _ASSERT_EQ(row_end - row_start + 1, subA.getNrows());
2120 1 : _ASSERT_EQ(col_end - col_start + 1, subA.getNcols());
2121 :
2122 2 : for (size_t i = 0; i < subA.getNrows(); i++) {
2123 2 : for (size_t j = 0; j < subA.getNcols(); j++) {
2124 1 : _ASSERT_NUM_EQ(A.get(row_start + i, col_start + j), subA.get(i, j), tol);
2125 : }
2126 1 : }
2127 1 : }
2128 :
2129 1 : void TestMatrix::testSubmatrixTranspose() {
2130 1 : const size_t n = 15;
2131 1 : const size_t m = 12;
2132 1 : const double tol = 1e-12;
2133 1 : Matrix A = MatrixFactory::MakeRandomMatrix(m, n, 2.0, 10.0, Matrix::MATRIX_DENSE);
2134 :
2135 1 : A.transpose();
2136 :
2137 : size_t row_start_tr;
2138 : size_t col_start_tr;
2139 : size_t rows_tr;
2140 : size_t cols_tr;
2141 2 : Matrix subA;
2142 :
2143 6 : for (rows_tr = 0; rows_tr <= 4; rows_tr++) {
2144 30 : for (cols_tr = 0; cols_tr <= 4; cols_tr++) {
2145 325 : for (row_start_tr = 0; row_start_tr < n - rows_tr - 1; row_start_tr++) {
2146 3000 : for (col_start_tr = 0; col_start_tr < m - cols_tr - 1; col_start_tr++) {
2147 2700 : subA = A.submatrixCopy(row_start_tr, row_start_tr + rows_tr, col_start_tr, col_start_tr + cols_tr);
2148 10350 : for (size_t i = 0; i < subA.getNrows(); i++) {
2149 28900 : for (size_t j = 0; j < subA.getNcols(); j++) {
2150 21250 : _ASSERT_NUM_EQ(A.get(row_start_tr + i, col_start_tr + j), subA.get(i, j), tol);
2151 : }
2152 : }
2153 : }
2154 : }
2155 : }
2156 1 : }
2157 :
2158 1 : }
2159 :
2160 4 : void _testSubmatrixMultiply(Matrix& A, Matrix& B) {
2161 4 : Matrix Asub;
2162 8 : Matrix Bsub;
2163 8 : Matrix exact;
2164 8 : Matrix result;
2165 :
2166 :
2167 16 : for (size_t k = 0; k < 3; k++) {
2168 48 : for (size_t ds = 0; ds < 3; ds++) {
2169 261 : for (size_t s = 0; s < A.getNcols() - ds - 1; s++) {
2170 1575 : for (size_t dj = 0; dj < 6; dj++) {
2171 8397 : for (size_t j = 0; j < B.getNcols() - dj - 1; j++) {
2172 21141 : for (size_t di = 0; di < 2; di++) {
2173 122499 : for (size_t i = 0; i < A.getNrows() - di - 1; i++) {
2174 108405 : Asub = A.submatrixCopy(i, i + di, s, s + ds);
2175 108405 : Bsub = B.submatrixCopy(k, k + ds, j, j + dj);
2176 108405 : exact = Asub*Bsub;
2177 216810 : result = A.multiplySubmatrix(B,
2178 : i, i + di, s, s + ds,
2179 108405 : k, k + ds, j, j + dj);
2180 108405 : _ASSERT_EQ(exact, result);
2181 108405 : _ASSERT_EQ(di + 1, result.getNrows());
2182 108405 : _ASSERT_EQ(dj + 1, result.getNcols());
2183 : }
2184 : }
2185 : }
2186 : }
2187 : }
2188 : }
2189 4 : }
2190 4 : }
2191 :
2192 1 : void TestMatrix::testSubmatrixMultiply() {
2193 1 : const size_t MA = 11;
2194 1 : const size_t NA = 10;
2195 1 : const size_t MB = 11;
2196 1 : const size_t NB = 9;
2197 1 : Matrix A = MatrixFactory::MakeRandomMatrix(MA, NA, 0.0, 10.0, Matrix::MATRIX_DENSE);
2198 2 : Matrix B = MatrixFactory::MakeRandomMatrix(MB, NB, 0.0, 2.0, Matrix::MATRIX_DENSE);
2199 2 : _testSubmatrixMultiply(A, B);
2200 1 : }
2201 :
2202 1 : void TestMatrix::testSubmatrixMultiplyTr() {
2203 1 : const size_t MA = 7;
2204 1 : const size_t NA = 9;
2205 1 : const size_t MB = 8;
2206 1 : const size_t NB = 10;
2207 1 : Matrix A = MatrixFactory::MakeRandomMatrix(MA, NA, 0.0, 10.0, Matrix::MATRIX_DENSE);
2208 2 : Matrix B = MatrixFactory::MakeRandomMatrix(MB, NB, 0.0, 2.0, Matrix::MATRIX_DENSE);
2209 1 : A.transpose();
2210 1 : _testSubmatrixMultiply(A, B); // A is transposed
2211 1 : A.transpose();
2212 1 : B.transpose();
2213 1 : _testSubmatrixMultiply(A, B); // B is transposed
2214 1 : A.transpose();
2215 2 : _testSubmatrixMultiply(A, B); // both A and B are transposed
2216 1 : }
2217 :
2218 1 : void TestMatrix::testSubmatrixSparse() {
2219 1 : const double tol = 1e-12;
2220 :
2221 1 : const size_t n = 14;
2222 1 : const size_t m = 11;
2223 :
2224 1 : const size_t row_start = 1;
2225 1 : const size_t row_end = 1;
2226 1 : const size_t col_start = 8;
2227 1 : const size_t col_end = 8;
2228 :
2229 1 : Matrix A = MatrixFactory::MakeRandomSparse(m, n, m + n, 0.0, 1.0);
2230 2 : Matrix subA = A.submatrixCopy(row_start, row_end, col_start, col_end);
2231 1 : _ASSERT_EQ(row_end - row_start + 1, subA.getNrows());
2232 1 : _ASSERT_EQ(col_end - col_start + 1, subA.getNcols());
2233 :
2234 2 : for (size_t i = 0; i < subA.getNrows(); i++) {
2235 2 : for (size_t j = 0; j < subA.getNcols(); j++) {
2236 1 : _ASSERT_NUM_EQ(A.get(row_start + i, col_start + j), subA.get(i, j), tol);
2237 : }
2238 1 : }
2239 1 : }
2240 :
2241 1 : void TestMatrix::test_ADTDT() {
2242 1 : size_t n = 15;
2243 1 : size_t m = 25;
2244 :
2245 1 : Matrix D1 = MatrixFactory::MakeRandomMatrix(n, m, 0.0, 1.0, Matrix::MATRIX_DENSE);
2246 2 : Matrix D2 = MatrixFactory::MakeRandomMatrix(n, m, 0.0, 1.0, Matrix::MATRIX_DENSE);
2247 2 : Matrix C1(D1);
2248 2 : Matrix C2(D2);
2249 1 : D1.transpose();
2250 1 : D2.transpose();
2251 :
2252 2 : Matrix R = D1 + D2;
2253 2 : Matrix C = C1 + C2;
2254 1 : C.transpose();
2255 2 : _ASSERT_EQ(R, C);
2256 :
2257 1 : }
2258 :
2259 1 : void TestMatrix::testGetSetTranspose() {
2260 1 : size_t n = 5;
2261 1 : size_t m = 6;
2262 1 : Matrix A(n, m);
2263 1 : A.set(1, 3, 10.4);
2264 :
2265 1 : A.transpose();
2266 1 : _ASSERT_EQ(10.4, A.get(3, 1));
2267 :
2268 1 : A.set(4, 3, 0.5);
2269 1 : _ASSERT_EQ(0.5, A.get(4, 3));
2270 1 : _ASSERT_EQ(0.0, A.get(3, 4));
2271 1 : }
2272 :
2273 1 : void TestMatrix::test_ASTDT() {
2274 1 : size_t n = 5;
2275 1 : size_t m = 6;
2276 :
2277 1 : Matrix D1 = MatrixFactory::MakeRandomSparse(n, m, 30, 0.0, 1.0);
2278 2 : Matrix D2 = MatrixFactory::MakeRandomMatrix(n, m, 0.0, 1.0, Matrix::MATRIX_DENSE);
2279 :
2280 1 : D1.transpose();
2281 1 : D2.transpose();
2282 :
2283 2 : Matrix F = D1 + D2;
2284 :
2285 1 : _ASSERT_EQ(m, F.getNrows());
2286 1 : _ASSERT_EQ(n, F.getNcols());
2287 :
2288 1 : const double tol = 1e-8;
2289 7 : for (size_t i = 0; i < m; i++) {
2290 36 : for (size_t j = 0; j < n; j++) {
2291 30 : _ASSERT_NUM_EQ(D1.get(i, j) + D2.get(i, j), F.get(i, j), tol);
2292 : }
2293 1 : }
2294 4 : }
2295 :
2296 :
2297 :
|
