ARTS 2.5.9 (git: 825fa5f2)
test_sparse.cc
Go to the documentation of this file.
1/* Copyright (C) 2003-2012
2 Stefan Buehler <sbuehler@ltu.se>
3 Mattias Ekstroem <ekstrom@rss.chalmers.se>
4
5 This program is free software; you can redistribute it and/or modify it
6 under the terms of the GNU General Public License as published by the
7 Free Software Foundation; either version 2, or (at your option) any
8 later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with this program; if not, write to the Free Software
17 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
18 USA. */
19
30#include <iostream>
31#include <stdexcept>
32#include "lin_alg.h"
33#include "matpackI.h"
34#include "matpackII.h"
35#include "test_utils.h"
36#include "xml_io.h"
37
38using std::cout;
39using std::endl;
40using std::max;
41using std::setw;
42
43void test3() {
44 Sparse M(10, 15);
45
46 /*
47 cout << "M.nrows(), M.ncols() = "
48 << M.nrows() << ", " << M.ncols() << "\n";
49 */
50 for (Index i = 3; i < 10; ++i) M.rw(i, i) = (Numeric)i + 1;
51 M.rw(0, 0) = 1;
52 M.rw(0, 1) = 2;
53 M.rw(0, 2) = 3;
54
55 cout << "\nM = \n" << M;
56
57 /*
58 // Test Sparse matrix-Matrix multiplication
59 Matrix A(10,5);
60 Matrix C(15,5,2.0);
61 // C = 2;
62
63 mult(A, M, C(Range(joker), Range(joker)));
64 cout << "\nA = \n" << A << "\n";
65
66 */
67
68 // Test Sparse-Sparse multiplication
69 // Sparse A(10,5);
70 // Sparse C(15,5);
71 // for (Index i=0; i<5; i++) {
72 // C.rw(i*3,i) = i*3+1;
73 // C.rw(i*3+1,i) = i*3+2;
74 // C.rw(i*3+2,i) = i*3+3;
75 // }
76
77 // mult(A,M,C);
78
79 // cout << "\nA = \n" << A;
80
81 /*
82 // Test transpose
83 Sparse B(15,10);
84
85 transpose(B,M);
86 cout << "\nM' = \n" << B;
87 */
88
89 /*
90 // Test rw-operator
91 Sparse S(M);
92 S.rw(2,0) = 5;
93
94 cout << "\nS(2,0) = " << S.ro(2,0) << "\n";
95
96 cout << "\nS = \n" << S;
97 */
98
99 /*
100 // Test vector multiplication
101 Vector y(20, 0.0);
102 Vector x(1,30,1);
103
104 mult(y[Range(1,10,2)], S, x[Range(1,15,2)]);
105
106 cout << "\ny = \n" << y << "\n";
107 */
108}
109
110// void test38()
111// {
112// cout << "Test sparse matrix - sparse matrix multiplication\n";
113
114// Sparse B(757,2271);
115// Sparse C(B.ncols(),B.ncols());
116// Sparse A(B.nrows(),B.ncols());
117
118// Index i=0;
119// for (Index j=0; j<B.ncols(); j++) {
120// B.rw(i,j) = (j+1.0);
121// if( i<B.nrows()-1 )
122// i++;
123// else
124// i=0;
125// }
126
127// for (Index i=0; i<C.nrows(); i++)
128// C.rw(i,i) = 1;
129
130// mult(A,B,C);
131
132// cout << "\n(Sparse) A = \n" << A;
133// cout << "\n(Sparse) B = \n" << B;
134// cout << "\n(Sparse) C = \n" << C;
135
136// Matrix a(5,15), b(5,15), c(15,15);
137
138// i=0;
139// for (Index j=0; j<15; j++) {
140// b(i,j) = j+1;
141// if( i<4 )
142// i++;
143// else
144// i=0;
145// }
146
147// for (Index i=0; i<15; i++)
148// c(i,i) = 1;
149
150// mult(a,b,c);
151
152// //cout << "\n(Full) a = \n" << a << "\n";
153
154// // cout << "\n(Sparse) B = \n" << B << "\n";
155// // cout << "\n(Full) b = \n" << b << "\n";
156// // cout << "\n(Sparse) C = \n" << C << "\n";
157// // cout << "\n(Full) c = \n" << c << "\n";
158
159// }
160
161// void test39()
162// {
163// //Test sparse transpose function
164// Sparse B(1000,2000);
165// Sparse Bt(B.ncols(), B.nrows());
166
167// Index i=0;
168// for (Index j=0; j<B.ncols(); j=j+2) {
169// B.rw(i,j) = j+1;
170// if( i<B.nrows()-2 )
171// i += 2;
172// else
173// i=0;
174// }
175
176// cout << "\nB = \n" << B;
177
178// transpose(Bt,B);
179
180// cout << "\ntranspose(B) = \n" << Bt << "\n";
181// }
182
183void test40() {
184 cout << "Testing the new simplified Sparse matrices:\n";
185
186 Sparse A(3, 3);
187 cout << "Empty A: " << A << "\n";
188
189 A.rw(0, 0) = 11;
190 A.rw(1, 1) = 22;
191 A.rw(2, 2) = 33;
192 cout << "Diagonal A:\n" << A << "\n";
193
194 Vector b(1, 3, 1), c(3);
195 cout << "b:\n" << b << "\n";
196
197 mult(c, A, b);
198 cout << "c = A*b (should be [11,44,99]):\n" << c << "\n";
199
200 Matrix D(3, 2);
201 D(joker, 0) = b;
202 D(joker, 1) = b;
203 D(joker, 1) *= 2;
204 cout << "D:\n" << D << "\n";
205
206 Matrix E(3, 2);
207 mult(E, A, D);
208 cout << "E = A*D (should be [11,22],[44,88],[99,198]):\n" << E << "\n";
209}
210
211void test41() {
212 cout << "Testing transpose for the new simplified sparse matrices:\n";
213
214 Sparse B(4, 5);
215 Index r[] = {0, 1, 1, 2, 2, 2, 3, 1, 3};
216 Index c[] = {0, 0, 1, 1, 2, 3, 3, 4, 4};
217 for (Index i = 0; i < 9; i++) B.rw(r[i], c[i]) = (Numeric)(i + 1);
218
219 cout << "B:\n" << B << "\n";
220
221 Sparse A(5, 4);
222
223 transpose(A, B);
224
225 cout << "A:\n" << A << "\n";
226
227 cout << "Testing with a fully occupied matrix:\n";
228
229 for (Index ri = 0; ri < 4; ri++)
230 for (Index ci = 0; ci < 5; ci++) {
231 B.rw(ri, ci) = (Numeric)(ri * 10 + ci);
232 }
233
234 cout << "B:\n" << B << "\n";
235 transpose(A, B);
236 cout << "A:\n" << A << "\n";
237}
238
239void test42() {
240 cout << "Testing sparse-sparse matrix multiplication:\n";
241
242 Sparse B(4, 5);
243 Index r[] = {0, 1, 1, 2, 2, 2, 3, 1, 3};
244 Index c[] = {0, 0, 1, 1, 2, 3, 3, 4, 4};
245 for (Index i = 0; i < 9; i++) B.rw(r[i], c[i]) = (Numeric)(i + 1);
246
247 Sparse A(4, 4), Bt(5, 4);
248 transpose(Bt, B);
249 mult(A, B, Bt);
250
251 cout << "A:\n" << A << "\n";
252}
253
254void test43() {
255 cout << "Testing sparse copying:\n";
256
257 Sparse B(4, 5);
258 Index r[] = {0, 1, 1, 2, 2, 2, 3, 1, 3};
259 Index c[] = {0, 0, 1, 1, 2, 3, 3, 4, 4};
260 for (Index i = 0; i < 9; i++) B.rw(r[i], c[i]) = (Numeric)(i + 1);
261
262 cout << "B:\n" << B << "\n";
263
264 Sparse A;
265
266 A = B;
267
268 cout << "A:\n" << A << "\n";
269
270 for (Index i = 0; i < 100; ++i) {
271 B.rw(0, 0) += 1;
272 A = B;
273 }
274
275 cout << "A now:\n" << A << "\n";
276}
277
278void test44() {
279 cout << "Test to insert row in sparse:\n";
280
281 Vector v(5, 10);
282
283 Sparse B(4, 5);
284 Index r[] = {0, 1, 1, 2, 2, 2, 3, 1, 3};
285 Index c[] = {0, 0, 1, 1, 2, 3, 3, 4, 4};
286 for (Index i = 0; i < 9; i++) B.rw(r[i], c[i]) = (Numeric)(i + 1);
287
288 cout << "B[" << B.nrows() << "," << B.ncols() << "]:\n" << B << "\n";
289 cout << "v:\n" << v << "\n";
290
291 B.insert_row(3, v);
292
293 cout << "B (after insertion):\n" << B << "\n";
294}
295
296void test45() {
297 cout << "Test Sparse-Sparse multiplication reading matrices from xml "
298 "files:\n";
299
300 Sparse A, B;
301 String a = "antenna.xml";
302 String b = "backend.xml";
303
304 try {
305 cout << " Reading " << a << "...";
307 cout << "done.\n Reading " << b << "...";
309 cout << "done.\n";
310 } catch (const std::runtime_error &e) {
311 cerr << e.what() << endl;
312 }
313
314 Sparse C(B.nrows(), A.ncols());
315 cout << " Performing multiplication...";
316 mult(C, B, A);
317 cout << "done.\n";
318
319 //cout << "C=A*B:\n" << A << "\n";
320 try {
321 cout << " Writing product to file: test45.xml...";
322 xml_write_to_file("test45.xml", C, FILE_TYPE_ASCII, 0, Verbosity());
323 cout << "done.\n";
324 } catch (const std::runtime_error &e) {
325 cerr << e.what() << endl;
326 }
327}
328
329void test46() {
330 cout << "Test transpose with large matrix read from xml file:\n";
331
332 Sparse A;
333 String a = "backend.xml";
334
335 try {
336 cout << " Reading " << a << "...";
338 cout << "done.\n";
339 } catch (const std::runtime_error &e) {
340 cerr << e.what() << endl;
341 }
342
343 //cout << "A:\n" << A << endl;
344
345 Sparse B(A.ncols(), A.nrows());
346 transpose(B, A);
347
348 try {
349 cout << " Writing transpose(A) to file test46.xml" << endl;
350 xml_write_to_file("test46.xml", B, FILE_TYPE_ASCII, 0, Verbosity());
351 } catch (const std::runtime_error &e) {
352 cerr << e.what() << endl;
353 }
354
355 //cout << "transpose(A):\n" << B << endl;
356}
357
358void test47() {
359 cout << "Test make Identity matrix:\n";
360
361 Sparse A;
362
363 A.resize(5, 5);
364 id_mat(A);
365
366 cout << "A:\n" << A << endl;
367}
368
369void test48() {
370 cout << "Test absolute values of sparse matrix:\n";
371
372 Sparse B(4, 5);
373 Index r[] = {0, 1, 1, 2, 2, 2, 3, 1, 3};
374 Index c[] = {0, 0, 1, 1, 2, 3, 3, 4, 4};
375 for (Index i = 0; i < 9; i++) B.rw(r[i], c[i]) = -(Numeric)i * 0.5;
376 cout << "B:\n" << B << endl;
377
378 Sparse A(B);
379 abs(A, B);
380
381 cout << "abs(B):\n" << A << endl;
382}
383
384void test49() {
385 cout << "Testing sparse adding:\n";
386
387 Sparse B(4, 5);
388 Index rb[] = {1, 3};
389 Index cb[] = {1, 3};
390 for (Index i = 0; i < 2; i++) B.rw(rb[i], cb[i]) = (Numeric)(i + 1);
391
392 Sparse C(4, 5);
393 Index rc[] = {0, 1, 2};
394 Index cc[] = {0, 1, 2};
395 for (Index i = 0; i < 3; i++) {
396 C.rw(rc[i], cc[i]) = (Numeric)(i + 1);
397 cout << C.rw(rc[i], cc[i]) << endl;
398 }
399
400 cout << "B:\n" << B << "\n";
401 cout << "C:\n" << C << "\n";
402
403 Sparse A;
404 add(A, B, C);
405 cout << "A=B+C:\n" << A << "\n";
406 Sparse D;
407 sub(D, B, C);
408 cout << "D=B-C:\n" << D << "\n";
409}
410
411Numeric test_xml_io(Index ntests, bool verbose) {
412 if (verbose) cout << "Testing xml IO." << endl << endl;
413
414 Index n, m;
415
416 Numeric err_max = 0.0;
417
418 for (int i = 0; i < ntests; i++) {
419 m = rand() % 1000 + 1;
420 n = rand() % 1000 + 1;
421
422 Sparse A(m, n), B;
423 String a("A.xml");
424 random_fill_matrix(A, 10, false);
428
429 Numeric err =
430 get_maximum_error(static_cast<Matrix>(B), static_cast<Matrix>(A), true);
431 if (err > err_max) err_max = err;
432
433 if (verbose)
434 cout << "Test " << i << ": "
435 << "Max. Rel. Error = " << err << endl;
436 }
437
438 return err_max;
439}
440
442
457Numeric test_insert_row(Index ntests, bool verbose) {
458 Numeric err_max = 0.0;
459
460 Vector v;
461 Matrix A;
462 Sparse A_sparse;
463
464 if (verbose) cout << endl << "Testing insert_row:" << endl << endl;
465 ;
466
467 for (Index i = 0; i < ntests; i++) {
468 Index m = (std::rand() % 10) + 1;
469 Index n = (std::rand() % 10) + 1;
470 Index r = std::rand() % m;
471
472 v.resize(n);
473 random_fill_vector(v, 10, false);
474
475 A.resize(m, n);
476
477 // Test resizing as well.
478 A_sparse.resize(m, n);
479 if ((A_sparse.nrows() != m) || (A_sparse.ncols() != n)) {
480 if (verbose) {
481 cout << "FAILED: Resize failed." << endl;
482 return 1.0;
483 }
484 }
485
486 A_sparse.insert_row(r, v);
487
488 A = static_cast<Matrix>(A_sparse);
489
490 Numeric err = get_maximum_error(A(r, joker), v, true);
491 if (err > err_max) err_max = err;
492
493 if (verbose) {
494 cout << endl;
495 cout << "Maximum relative error: " << err << endl;
496 }
497 }
498
499 return err_max;
500}
501
503
519Numeric test_identity(Index ntests, bool verbose) {
520 Numeric err_max = 0.0;
521
522 Matrix A;
523 Sparse A_sparse;
524
525 if (verbose) cout << endl << "Testing sparse identity:" << endl << endl;
526
527 for (Index i = 0; i < ntests; i++) {
528 Index n = (std::rand() % 1000) + 1;
529
530 A.resize(n, n);
531 A_sparse.resize(n, n);
532
533 if ((A_sparse.nrows() != n) || (A_sparse.ncols() != n)) {
534 if (verbose) {
535 cout << endl << "FAILED: Resize failed." << endl;
536 return 1.0;
537 }
538 }
539
540 id_mat(A);
541 id_mat(A_sparse);
542
543 Numeric err = get_maximum_error(static_cast<Matrix>(A_sparse), A, true);
544 if (err > err_max) err_max = err;
545
546 if (verbose) {
547 cout << "\t Maximum relative error: " << err << endl;
548 }
549 }
550
551 return err_max;
552}
553
555
570Numeric test_sparse_construction(Index m, Index n, Index ntests, bool verbose) {
571 Numeric err_max = 0.0;
572
573 if (verbose) {
574 cout << endl;
575 cout << "Testing sparse matrix construction and conversion:" << endl;
576 }
577
578 for (Index i = 0; i < ntests; i++) {
579 Matrix A;
580 Sparse A_sparse(m, n);
581
582 random_fill_matrix(A_sparse, 10, false);
583 A = static_cast<Matrix>(A_sparse);
584
585 Numeric err = get_maximum_error(static_cast<Matrix>(A_sparse), A, true);
586 if (err > err_max) err_max = err;
587
588 if (verbose) {
589 cout << "Test No.: " << setw(5) << i;
590 cout << " Maximum rel. error: " << err_max << endl;
591 }
592 }
593 return err_max;
594}
595
597
614 Index n,
615 Index ntests,
616 bool verbose) {
617 Numeric err_max = 0.0;
618
619 if (verbose) {
620 cout << endl;
621 cout << "Testing sparse unary operations:" << endl << endl;
622 cout << setw(5) << "Test " << setw(15) << "Construction";
623 cout << setw(15) << "abs" << setw(15) << "transpose" << endl;
624 cout << std::string(55, '-') << endl;
625 }
626
627 for (Index i = 0; i < ntests; i++) {
628 Matrix A(m, n), B(m, n);
629 Sparse A_sparse(m, n), B_sparse(m, n), C_sparse(n, m);
630
631 A = 0.0;
632 random_fill_matrix(A, A_sparse, 10, false);
633
634 // Construction
635
636 Numeric err = get_maximum_error(static_cast<Matrix>(A_sparse), A, true);
637 if (err > err_max) err_max = err;
638
639 if (verbose) cout << setw(5) << i << setw(15) << err;
640
641 // abs
642
643 for (Index j = 0; j < m; j++) {
644 for (Index k = 0; k < n; k++) {
645 B(j, k) = abs(A(j, k));
646 }
647 }
648
649 abs(B_sparse, A_sparse);
650
651 err = get_maximum_error(static_cast<Matrix>(B_sparse), B, true);
652 if (err > err_max) err_max = err;
653
654 if (verbose) cout << setw(15) << err;
655
656 // transpose
657
658 transpose(C_sparse, A_sparse);
659
660 err = get_maximum_error(static_cast<Matrix>(C_sparse), transpose(A), true);
661 if (err > err_max) err_max = err;
662
663 if (verbose) cout << setw(15) << err << endl;
664 }
665 return err_max;
666}
667
669
686 Index n,
687 Index ntests,
688 bool verbose) {
689 Numeric err_max = 0.0;
690
691 if (verbose) {
692 cout << endl;
693 cout << "Testing sparse multiplication:" << endl << endl;
694 cout << setw(5) << "Test " << setw(5) << "m1" << setw(5) << "m2";
695 cout << setw(10) << "c_stride" << setw(5) << "n1" << setw(5) << "n2";
696 cout << setw(10) << setw(10) << "r_stride" << setw(15) << "A = B x C";
697 cout << setw(18) << "A^T = C^T x B^T" << endl;
698 cout << std::string(80, '-') << endl;
699 }
700
701 Matrix A(m, n);
702 Matrix A_ref(m, n);
703 Matrix B(m, n);
704 random_fill_matrix(B, 10.0, true);
705
706 for (Index i = 0; i < ntests; i++) {
707 //
708 // Generate random sub-matrix.
709 //
710
711 Index m1, dm, n1, dn, c_stride, r_stride;
712
713 m1 = rand() % (m - 1);
714 n1 = rand() % (n - 1);
715 dm = (rand() % (m - m1)) + 1;
716 dn = (rand() % (n - n1)) + 1;
717 c_stride = (rand() % 10) + 1;
718 r_stride = (rand() % 10) + 1;
719
720 c_stride = min(dn, c_stride);
721 r_stride = min(dm, r_stride);
722 dn = dn / c_stride;
723 dm = dm / r_stride;
724
725 if (verbose) {
726 cout << setw(5) << i << setw(5) << m1 << setw(5) << dm;
727 cout << setw(10) << r_stride << setw(5) << n1 << setw(5) << dn;
728 cout << setw(10) << c_stride;
729 }
730
731 MatrixView A_view = A(Range(m1, dm, r_stride), joker);
732 MatrixView A_ref_view = A_ref(Range(m1, dm, r_stride), joker);
733 Matrix C(dn, n);
734 MatrixView B_mul = B(Range(m1, dm, r_stride), Range(n1, dn, c_stride));
735 Sparse C_sparse(dn, n), C_sparse_transpose(dn, n);
736
737 //
738 // Test standard dense-sparse multiplication.
739 //
740
741 random_fill_matrix(C_sparse, 10, false);
742 C = static_cast<Matrix>(C_sparse);
743
744 mult(A_ref_view, B_mul, C);
745 mult(A_view, B_mul, C_sparse);
746
747 // Compare results
748 Numeric err = get_maximum_error(A_view, A_ref_view, true);
749 if (err > err_max) err_max = err;
750
751 if (verbose) {
752 cout << setw(15) << err;
753 }
754
755 //
756 // Test transposed dense-sparse multiplication.
757 //
758
759 A.resize(n, m);
760 A_ref.resize(n, m);
761 MatrixView A_view_transp = A(joker, Range(m1, dm, r_stride));
762 MatrixView A_ref_view_transp = A_ref(joker, Range(m1, dm, r_stride));
763
764 C_sparse.resize(n, dn);
765 C.resize(n, dn);
766 random_fill_matrix(C_sparse, 10, false);
767 transpose(C_sparse_transpose, C_sparse);
768 C = static_cast<Matrix>(C_sparse);
769
770 MatrixView B_mul_transp =
771 B(Range(n1, dn, c_stride), Range(m1, dm, r_stride));
772
773 mult(transpose(A_ref_view_transp), transpose(B_mul_transp), transpose(C));
774 mult(transpose(A_view_transp), transpose(B_mul_transp), C_sparse_transpose);
775
776 // Compare results
777
778 err = get_maximum_error(A_view, A_ref_view, true);
779 if (err > err_max) err_max = err;
780
781 if (verbose) {
782 cout << setw(15) << err << endl;
783 }
784 }
785 return err_max;
786}
787
789
806 Index n,
807 Index ntests,
808 bool verbose) {
809 Numeric err_max = 0.0;
810
811 if (verbose) {
812 cout << endl;
813 cout << "Testing sparse multiplication:" << endl << endl;
814 cout << setw(5) << "Test " << setw(5) << "m1" << setw(5) << "m2";
815 cout << setw(10) << "c_stride" << setw(5) << "n1" << setw(5) << "n2";
816 cout << setw(10) << setw(10) << "r_stride" << setw(15) << "A = B x C";
817 cout << setw(18) << "A^T = C^T x B^T" << endl;
818 cout << std::string(80, '-') << endl;
819 }
820
821 Matrix C(m, n);
822 Matrix A(m, n);
823 Matrix A_ref(m, n);
824 random_fill_matrix(C, 10.0, true);
825
826 for (Index i = 0; i < ntests; i++) {
827 //
828 // Generate random sub-matrix.
829 //
830
831 Index m1, dm, n1, dn, c_stride, r_stride;
832
833 m1 = rand() % (m - 1);
834 n1 = rand() % (n - 1);
835 dm = (rand() % (m - m1)) + 1;
836 dn = (rand() % (n - n1)) + 1;
837 c_stride = (rand() % 10) + 1;
838 r_stride = (rand() % 10) + 1;
839
840 if (verbose) {
841 cout << setw(5) << i << setw(5) << m1 << setw(5) << dm;
842 cout << setw(10) << c_stride << setw(5) << n1 << setw(5) << dn;
843 cout << setw(10) << r_stride;
844 }
845
846 c_stride = min(dn, c_stride);
847 r_stride = min(dm, r_stride);
848 dn = dn / c_stride;
849 dm = dm / r_stride;
850
851 Matrix B(m, dm);
852 MatrixView A_view = A(joker, Range(n1, dn, c_stride));
853 MatrixView A_ref_view = A_ref(joker, Range(n1, dn, c_stride));
854 MatrixView C_mul = C(Range(m1, dm, r_stride), Range(n1, dn, c_stride));
855 Sparse B_sparse(m, dm), B_sparse_transpose(m, dm);
856
857 //
858 // Test standard multiplication.
859 //
860
861 random_fill_matrix(static_cast<Matrix>(B_sparse), 10, false);
862 B = static_cast<Matrix>(B_sparse);
863
864 // Sparse-sparse multiplication
865 mult(A_ref_view, B, C_mul);
866 mult(A_view, B_sparse, C_mul);
867
868 // Compare results
869 Numeric err = get_maximum_error(A_view, A_ref_view, true);
870 if (err > err_max) err_max = err;
871
872 if (verbose) {
873 cout << setw(15) << err;
874 }
875
876 //
877 // Test transposed matrix views.
878 //
879
880 A.resize(n, m);
881 A_ref.resize(n, m);
882 MatrixView A_view_transp = A(Range(n1, dn, c_stride), joker);
883 MatrixView A_ref_view_transp = A_ref(Range(n1, dn, c_stride), joker);
884
885 B_sparse.resize(dm, m);
886 B.resize(dm, m);
887 random_fill_matrix(B_sparse, 10, false);
888 transpose(B_sparse_transpose, B_sparse);
889 B = static_cast<Matrix>(B_sparse);
890
891 MatrixView C_mul2 = C(Range(n1, dn, c_stride), Range(m1, dm, r_stride));
892
893 // Transposed sparse-dense multiplication
894 mult(transpose(A_ref_view_transp), transpose(B), transpose(C_mul2));
895 mult(transpose(A_view_transp), B_sparse_transpose, transpose(C_mul2));
896
897 // Compare results
898 err = get_maximum_error(A_view, A_ref_view, true);
899 if (err > err_max) err_max = err;
900
901 if (verbose) {
902 cout << setw(15) << err << endl;
903 }
904 }
905 return err_max;
906}
907
909
933 Index n,
934 Index ntests,
935 bool verbose) {
936 Numeric err_max = 0.0;
937
938 if (verbose) {
939 cout << endl;
940 cout << "Testing sparse multiplication:" << endl << endl;
941 cout << setw(5) << "Test " << setw(15) << "sparse-sparse";
942 cout << setw(15) << "sparse-dense" << setw(15) << "matrix-vector";
943 cout << "transposed matrix-vector";
944 cout << endl << std::string(80, '-') << endl;
945 }
946
947 for (Index i = 0; i < ntests; i++) {
948 Index k;
949
950 k = (rand() % 1000) + 1;
951
952 Matrix A(m, n), A2(m, n), B, C;
953 Vector y(m), y2(m), yt(k), yt2(k), x(k), xt(m);
954 Sparse A_sparse(m, n), B_sparse(m, k), C_sparse(k, n);
955
956 random_fill_matrix(B_sparse, 10, false);
957 random_fill_matrix(C_sparse, 10, false);
958 random_fill_vector(x, 10, false);
959 random_fill_vector(xt, 10, false);
960 B = static_cast<Matrix>(B_sparse);
961 C = static_cast<Matrix>(C_sparse);
962
963 // Sparse-sparse
964 mult(A, B, C);
965 mult(A_sparse, B_sparse, C_sparse);
966
967 Numeric err = get_maximum_error(static_cast<Matrix>(A_sparse), A, true);
968 if (err > err_max) err_max = err;
969
970 if (verbose) {
971 cout << setw(5) << i << setw(15) << err;
972 }
973
974 // Sparse-dense
975 mult(A, B, C);
976 mult(A2, B_sparse, C);
977
978 err = get_maximum_error(A2, A, true);
979 if (err > err_max) err_max = err;
980
981 if (verbose) {
982 cout << setw(15) << err;
983 }
984
985 // Matrix-vector
986 mult(y, B, x);
987 mult(y2, B_sparse, x);
988
989 err = get_maximum_error(y2, y, true);
990 if (err > err_max) err_max = err;
991
992 if (verbose) {
993 cout << setw(15) << err;
994 }
995
996 // transposed Matrix-vector
997 mult(yt, transpose(B), xt);
998 transpose_mult(yt2, B_sparse, xt);
999
1000 err = get_maximum_error(yt2, yt, true);
1001 if (err > err_max) err_max = err;
1002
1003 if (verbose) {
1004 cout << setw(15) << err << endl;
1005 }
1006 }
1007 return err_max;
1008}
1009
1011
1028Numeric test_sparse_arithmetic(Index m, Index n, Index ntests, bool verbose) {
1029 Numeric err_max = 0.0;
1030
1031 if (verbose) {
1032 cout << endl;
1033 cout << "Testing sparse arithmetic:" << endl << endl;
1034 cout << setw(5) << "Test " << setw(15) << "Addition";
1035 cout << setw(15) << "Multiplication" << setw(15) << "Subtraction";
1036 cout << setw(15) << "Scaling" << endl;
1037 cout << std::string(70, '-') << endl;
1038 }
1039
1040 for (Index i = 0; i < ntests; i++) {
1041 Matrix A, B, C(m, n), D;
1042 Sparse A_sparse(m, m), B_sparse(m, n), C_sparse(m, n), D_sparse(m, n);
1043
1044 random_fill_matrix(A_sparse, 10, false);
1045 random_fill_matrix(B_sparse, 10, false);
1046 random_fill_matrix(D_sparse, 10, false);
1047 A = static_cast<Matrix>(A_sparse);
1048 B = static_cast<Matrix>(B_sparse);
1049 D = static_cast<Matrix>(D_sparse);
1050
1051 // Multiplication
1052
1053 mult(C, A, B);
1054 mult(C_sparse, A_sparse, B_sparse);
1055
1056 Numeric err = get_maximum_error(static_cast<Matrix>(C_sparse), C, true);
1057 if (err > err_max) err_max = err;
1058
1059 if (verbose) cout << setw(5) << i << setw(15) << err;
1060
1061 // Addition
1062
1063 C += D;
1064 add(B_sparse, C_sparse, D_sparse);
1065
1066 err = get_maximum_error(static_cast<Matrix>(B_sparse), C, true);
1067 if (err > err_max) err_max = err;
1068
1069 C_sparse += D_sparse;
1070
1071 err = get_maximum_error(static_cast<Matrix>(C_sparse), C, true);
1072 if (err > err_max) err_max = err;
1073
1074 if (verbose) cout << setw(15) << err;
1075
1076 // Subtraction
1077
1078 C = static_cast<Matrix>(C_sparse);
1079 C -= D;
1080 sub(B_sparse, C_sparse, D_sparse);
1081
1082 err = get_maximum_error(static_cast<Matrix>(B_sparse), C, true);
1083 if (err > err_max) err_max = err;
1084
1085 C_sparse -= D_sparse;
1086
1087 err = get_maximum_error(static_cast<Matrix>(C_sparse), C, true);
1088 if (err > err_max) err_max = err;
1089
1090 if (verbose) cout << setw(15) << err;
1091
1092 // Scaling
1093
1094 C_sparse *= 1.2;
1095 C *= 1.2;
1096
1097 err = get_maximum_error(static_cast<Matrix>(C_sparse), C, true);
1098 if (err > err_max) err_max = err;
1099
1100 C_sparse /= 1.2;
1101 C /= 1.2;
1102
1103 err = get_maximum_error(static_cast<Matrix>(C_sparse), C, true);
1104 if (err > err_max) err_max = err;
1105
1106 if (verbose) cout << setw(15) << err << endl;
1107 }
1108 return err_max;
1109}
1110
1111int main() {
1112 // test3();
1113 // test38();
1114 // test39();
1115 // test40();
1116 // test41();
1117 // test42();
1118 // test43();
1119 // test44();
1120 // test45();
1121 // test46();
1122 // test47();
1123 // test48();
1124
1125 Numeric err;
1126
1127 cout << "Testing sparse arithmetic: ";
1128 err = test_sparse_arithmetic(1000, 1000, 100, false);
1129 if (err < 1e-11)
1130 cout << "PASSED" << endl;
1131 else
1132 cout << "FAILED (Error: " << err << ")" << endl;
1133
1134 cout << "Testing xml IO: ";
1135 err = test_xml_io(100, false);
1136 if (err < 1e-11)
1137 cout << "PASSED" << endl;
1138 else
1139 cout << "FAILED (Error: " << err << ")" << endl;
1140
1141 cout << "Testing identity matrix: ";
1142 err = test_identity(1000, false);
1143 if (err < 1e-11)
1144 cout << "PASSED" << endl;
1145 else
1146 cout << "FAILED (Error: " << err << ")" << endl;
1147
1148 cout << "Testing inserting of rows: ";
1149 err = test_insert_row(1000, false);
1150 if (err < 1e-11)
1151 cout << "PASSED" << endl;
1152 else
1153 cout << "FAILED (Error: " << err << ")" << endl;
1154
1155 cout << "Testing abs(...) and transpose(...): ";
1156 err = test_sparse_unary_operations(1000, 1000, 1000, false);
1157 if (err < 1e-11)
1158 cout << "PASSED" << endl;
1159 else
1160 cout << "FAILED (Error: " << err << ")" << endl;
1161
1162 cout << "Testing multiplication: ";
1163 err = test_sparse_multiplication(1000, 1000, 1000, false);
1164 if (err < 1e-11)
1165 cout << "PASSED" << endl;
1166 else
1167 cout << "FAILED (Error: " << err << ")" << endl;
1168
1169 cout << "Testing sparse-dense multiplication: ";
1170 err = test_sparse_dense_multiplication(1000, 1000, 1000, false);
1171 if (err < 1e-11)
1172 cout << "PASSED" << endl;
1173 else
1174 cout << "FAILED (Error: " << err << ")" << endl;
1175
1176 cout << "Testing dense-sparse multiplication: ";
1177 err = test_dense_sparse_multiplication(1000, 1000, 1000, false);
1178 if (err < 1e-11)
1179 cout << "PASSED" << endl;
1180 else
1181 cout << "FAILED (Error: " << err << ")" << endl;
1182
1183 return 0;
1184}
base min(const Array< base > &x)
Min function.
Definition: array.h:161
The MatrixView class.
Definition: matpackI.h:1188
The Matrix class.
Definition: matpackI.h:1285
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1011
The range class.
Definition: matpackI.h:160
The Vector class.
Definition: matpackI.h:910
void mult(MatrixView C, ConstMatrixView A, const Block &B)
void id_mat(MatrixView I)
Identity Matrix.
Definition: lin_alg.cc:488
Linear algebra functions.
void abs(Sparse &A, const Sparse &B)
Absolute value of sparse matrix elements.
Definition: matpackII.cc:394
void transpose_mult(VectorView y, const Sparse &M, ConstVectorView x)
Sparse matrix - Vector multiplication.
Definition: matpackII.cc:452
void sub(Sparse &A, const Sparse &B, const Sparse &C)
Sparse - Sparse subtraction.
Definition: matpackII.cc:667
void add(Sparse &A, const Sparse &B, const Sparse &C)
Sparse - Sparse addition.
Definition: matpackII.cc:630
Header file for sparse matrices.
Implementation of Matrix, Vector, and such stuff.
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
ConstComplexMatrixView transpose(ConstComplexMatrixView m)
Const version of transpose.
const Joker joker
#define M
Definition: rng.cc:165
The Sparse class.
Definition: matpackII.h:75
Numeric & rw(Index r, Index c)
Read and write index operator.
Definition: matpackII.cc:87
void resize(Index r, Index c)
Resize function.
Definition: matpackII.cc:361
void insert_row(Index r, Vector v)
Insert row function.
Definition: matpackII.cc:314
Index nrows() const
Returns the number of rows.
Definition: matpackII.cc:62
Index ncols() const
Returns the number of columns.
Definition: matpackII.cc:65
Numeric test_sparse_unary_operations(Index m, Index n, Index ntests, bool verbose)
Test unary operation on sparse matrices.
Definition: test_sparse.cc:613
void test46()
Definition: test_sparse.cc:329
Numeric test_sparse_dense_multiplication(Index m, Index n, Index ntests, bool verbose)
Test sparse-dense multiplication.
Definition: test_sparse.cc:805
void test42()
Definition: test_sparse.cc:239
void test48()
Definition: test_sparse.cc:369
void test3()
Definition: test_sparse.cc:43
Numeric test_sparse_arithmetic(Index m, Index n, Index ntests, bool verbose)
Test sparse matrix arithmetic.
void test44()
Definition: test_sparse.cc:278
Numeric test_insert_row(Index ntests, bool verbose)
Test sparse insert_row function.
Definition: test_sparse.cc:457
void test47()
Definition: test_sparse.cc:358
Numeric test_dense_sparse_multiplication(Index m, Index n, Index ntests, bool verbose)
Test dense-sparse multiplication.
Definition: test_sparse.cc:685
Numeric test_xml_io(Index ntests, bool verbose)
Definition: test_sparse.cc:411
void test49()
Definition: test_sparse.cc:384
void test45()
Definition: test_sparse.cc:296
Numeric test_identity(Index ntests, bool verbose)
Test sparse identity matrix.
Definition: test_sparse.cc:519
void test40()
Definition: test_sparse.cc:183
Numeric test_sparse_multiplication(Index m, Index n, Index ntests, bool verbose)
Test sparse multiplication.
Definition: test_sparse.cc:932
int main()
void test43()
Definition: test_sparse.cc:254
void test41()
Definition: test_sparse.cc:211
Numeric test_sparse_construction(Index m, Index n, Index ntests, bool verbose)
Test sparse matrix construction.
Definition: test_sparse.cc:570
void random_fill_matrix(MatrixView A, Numeric range, bool positive)
Fill matrix with random values.
Definition: test_utils.cc:59
Numeric get_maximum_error(ConstVectorView v1, ConstVectorView v2, bool relative)
Maximum element-wise error of two vectors.
Definition: test_utils.cc:297
void random_fill_vector(VectorView v, Numeric range, bool positive)
Fill vector with random values.
Definition: test_utils.cc:230
Utility functions for testing.
#define v
#define a
#define c
#define b
This file contains basic functions to handle XML data files.
void xml_write_to_file(const String &filename, const T &type, const FileType ftype, const Index no_clobber, const Verbosity &verbosity)
Write data to XML file.
Definition: xml_io.h:172
void xml_read_from_file(const String &filename, T &type, const Verbosity &verbosity)
Reads data from XML file.
Definition: xml_io.h:151
@ FILE_TYPE_ASCII
Definition: xml_io_base.h:43