ARTS  2.4.0(git:4fb77825)
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 
38 using std::cout;
39 using std::endl;
40 using std::max;
41 using std::setw;
42 
43 void 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 
183 void 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 
211 void 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 
239 void 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 
254 void 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 
278 void 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 
296 void 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 << "...";
306  xml_read_from_file(a, A, Verbosity());
307  cout << "done.\n Reading " << b << "...";
308  xml_read_from_file(b, B, Verbosity());
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 
329 void 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 << "...";
337  xml_read_from_file(a, A, Verbosity());
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 
358 void 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 
369 void 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 
384 void 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 
411 Numeric 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);
426  xml_read_from_file(a, B, Verbosity());
427  xml_write_to_file("B.xml", B, FILE_TYPE_ASCII, 0, Verbosity());
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 
457 Numeric 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 
519 Numeric 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 
570 Numeric 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 
1028 Numeric 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 
1111 int 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 }
Matrix
The Matrix class.
Definition: matpackI.h:1193
main
int main()
Definition: test_sparse.cc:1111
FILE_TYPE_ASCII
@ FILE_TYPE_ASCII
Definition: xml_io.h:38
MatrixView
The MatrixView class.
Definition: matpackI.h:1093
Sparse::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackII.cc:69
test_utils.h
Utility functions for testing.
id_mat
void id_mat(MatrixView I)
Identity Matrix.
Definition: lin_alg.cc:2240
test_sparse_multiplication
Numeric test_sparse_multiplication(Index m, Index n, Index ntests, bool verbose)
Test sparse multiplication.
Definition: test_sparse.cc:932
joker
const Joker joker
abs
#define abs(x)
Definition: legacy_continua.cc:20626
mult
void mult(ComplexVectorView y, const ConstComplexMatrixView &M, const ConstComplexVectorView &x)
Matrix-Vector Multiplication.
Definition: complex.cc:1579
Sparse
The Sparse class.
Definition: matpackII.h:60
random_fill_vector
void random_fill_vector(VectorView v, Numeric range, bool positive)
Fill vector with random values.
Definition: test_utils.cc:230
ARTS::Var::y
Vector y(Workspace &ws) noexcept
Definition: autoarts.h:7401
Vector::resize
void resize(Index n)
Resize function.
Definition: matpackI.cc:404
test_identity
Numeric test_identity(Index ntests, bool verbose)
Test sparse identity matrix.
Definition: test_sparse.cc:519
test43
void test43()
Definition: test_sparse.cc:254
test_sparse_unary_operations
Numeric test_sparse_unary_operations(Index m, Index n, Index ntests, bool verbose)
Test unary operation on sparse matrices.
Definition: test_sparse.cc:613
Sparse::insert_row
void insert_row(Index r, Vector v)
Insert row function.
Definition: matpackII.cc:314
test40
void test40()
Definition: test_sparse.cc:183
test49
void test49()
Definition: test_sparse.cc:384
test_sparse_dense_multiplication
Numeric test_sparse_dense_multiplication(Index m, Index n, Index ntests, bool verbose)
Test sparse-dense multiplication.
Definition: test_sparse.cc:805
matpackI.h
Implementation of Matrix, Vector, and such stuff.
test47
void test47()
Definition: test_sparse.cc:358
test3
void test3()
Definition: test_sparse.cc:43
min
#define min(a, b)
Definition: legacy_continua.cc:20628
random_fill_matrix
void random_fill_matrix(MatrixView A, Numeric range, bool positive)
Fill matrix with random values.
Definition: test_utils.cc:59
xml_read_from_file
void xml_read_from_file(const String &filename, T &type, const Verbosity &verbosity)
Reads data from XML file.
Definition: xml_io.cc:901
my_basic_string< char >
sub
void sub(Sparse &A, const Sparse &B, const Sparse &C)
Sparse - Sparse subtraction.
Definition: matpackII.cc:667
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
lin_alg.h
Linear algebra functions.
test46
void test46()
Definition: test_sparse.cc:329
test_sparse_construction
Numeric test_sparse_construction(Index m, Index n, Index ntests, bool verbose)
Test sparse matrix construction.
Definition: test_sparse.cc:570
max
#define max(a, b)
Definition: legacy_continua.cc:20629
test41
void test41()
Definition: test_sparse.cc:211
Matrix::resize
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1056
test44
void test44()
Definition: test_sparse.cc:278
test_xml_io
Numeric test_xml_io(Index ntests, bool verbose)
Definition: test_sparse.cc:411
Sparse::resize
void resize(Index r, Index c)
Resize function.
Definition: matpackII.cc:361
Range
The range class.
Definition: matpackI.h:160
test48
void test48()
Definition: test_sparse.cc:369
xml_write_to_file
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.cc:972
transpose_mult
void transpose_mult(VectorView y, const Sparse &M, ConstVectorView x)
Sparse matrix - Vector multiplication.
Definition: matpackII.cc:452
Sparse::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackII.cc:66
Sparse::rw
Numeric & rw(Index r, Index c)
Read and write index operator.
Definition: matpackII.cc:91
matpackII.h
Header file for sparse matrices.
get_maximum_error
Numeric get_maximum_error(ConstVectorView v1, ConstVectorView v2, bool relative)
Maximum element-wise error of two vectors.
Definition: test_utils.cc:297
test_sparse_arithmetic
Numeric test_sparse_arithmetic(Index m, Index n, Index ntests, bool verbose)
Test sparse matrix arithmetic.
Definition: test_sparse.cc:1028
ARTS::Var::x
Vector x(Workspace &ws) noexcept
Definition: autoarts.h:7346
test_insert_row
Numeric test_insert_row(Index ntests, bool verbose)
Test sparse insert_row function.
Definition: test_sparse.cc:457
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
transpose
ConstComplexMatrixView transpose(ConstComplexMatrixView m)
Const version of transpose.
Definition: complex.cc:1509
M
#define M
Definition: rng.cc:165
test42
void test42()
Definition: test_sparse.cc:239
Vector
The Vector class.
Definition: matpackI.h:860
test_dense_sparse_multiplication
Numeric test_dense_sparse_multiplication(Index m, Index n, Index ntests, bool verbose)
Test dense-sparse multiplication.
Definition: test_sparse.cc:685
ARTS::Group::Verbosity
Verbosity Verbosity
Definition: autoarts.h:114
test45
void test45()
Definition: test_sparse.cc:296
xml_io.h
This file contains basic functions to handle XML data files.
add
void add(Sparse &A, const Sparse &B, const Sparse &C)
Sparse - Sparse addition.
Definition: matpackII.cc:630