ARTS 2.5.4 (git: 4c0d3b4d)
test_matpack.cc
Go to the documentation of this file.
1/* Copyright (C) 2001-2012 Stefan Buehler <sbuehler@ltu.se>
2
3 This program is free software; you can redistribute it and/or modify it
4 under the terms of the GNU General Public License as published by the
5 Free Software Foundation; either version 2, or (at your option) any
6 later version.
7
8 This program is distributed in the hope that it will be useful,
9 but WITHOUT ANY WARRANTY; without even the implied warranty of
10 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 GNU General Public License for more details.
12
13 You should have received a copy of the GNU General Public License
14 along with this program; if not, write to the Free Software
15 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
16 USA. */
17
18#include <algorithm>
19#include <cmath>
20#include <cstdlib>
21#include <exception>
22#include <iostream>
23#include "array.h"
24#include "describe.h"
25#include "exceptions.h"
26#include "logic.h"
27#include "math_funcs.h"
28#include "matpackII.h"
29#include "matpackVII.h"
30#include "mystring.h"
31#include "rational.h"
32#include "test_utils.h"
33#include "time.h"
34#include "wigner_functions.h"
35
36using std::cout;
37using std::endl;
38using std::runtime_error;
39
40Numeric by_reference(const Numeric& x) { return x + 1; }
41
42Numeric by_value(Numeric x) { return x + 1; }
43
44void fill_with_junk(VectorView x) { x = 999; }
45
46void fill_with_junk(MatrixView x) { x = 888; }
47
48int test1() {
49 Vector v(20);
50
51 cout << "v.nelem() = " << v.nelem() << "\n";
52
53 for (Index i = 0; i < v.nelem(); ++i) v[i] = (Numeric)i;
54
55 cout << "v.begin() = " << *v.begin() << "\n";
56
57 cout << "v = \n" << v << "\n";
58
59 fill_with_junk(v[Range(1, 8, 2)][Range(2, joker)]);
60 // fill_with_junk(v);
61
62 Vector v2 = v[Range(2, 4)];
63
64 cout << "v2 = \n" << v2 << "\n";
65
66 for (Index i = 0; i < 1000; ++i) {
67 Vector v3(1000);
68 v3 = (Numeric)i;
69 }
70
71 v2[Range(joker)] = 88;
72
73 v2[Range(0, 2)] = 77;
74
75 cout << "v = \n" << v << "\n";
76 cout << "v2 = \n" << v2 << "\n";
77 cout << "v2.nelem() = \n" << v2.nelem() << "\n";
78
79 Vector v3;
80 v3.resize(v2.nelem());
81 v3 = v2;
82
83 cout << "\nv3 = \n" << v3 << "\n";
85 cout << "\nv3 after junking v2 = \n" << v3 << "\n";
86 v3 *= 2;
87 cout << "\nv3 after *2 = \n" << v3 << "\n";
88
89 Matrix M(10, 15);
90 {
91 Numeric n = 0;
92 for (Index i = 0; i < M.nrows(); ++i)
93 for (Index j = 0; j < M.ncols(); ++j) M(i, j) = ++n;
94 }
95
96 cout << "\nM =\n" << M << "\n";
97
98 cout << "\nM(Range(2,4),Range(2,4)) =\n"
99 << M(Range(2, 4), Range(2, 4)) << "\n";
100
101 cout << "\nM(Range(2,4),Range(2,4))(Range(1,2),Range(1,2)) =\n"
102 << M(Range(2, 4), Range(2, 4))(Range(1, 2), Range(1, 2)) << "\n";
103
104 cout << "\nM(1,Range(joker)) =\n" << M(1, Range(joker)) << "\n";
105
106 cout << "\nFilling M(1,Range(1,2)) with junk.\n";
107 fill_with_junk(M(1, Range(1, 2)));
108
109 cout << "\nM(Range(0,4),Range(0,4)) =\n"
110 << M(Range(0, 4), Range(0, 4)) << "\n";
111
112 cout << "\nFilling M(Range(4,2,2),Range(6,3)) with junk.\n";
113
114 MatrixView s = M(Range(4, 2, 2), Range(6, 3));
116
117 cout << "\nM =\n" << M << "\n";
118
119 const Matrix C = M;
120
121 cout << "\nC(Range(3,4,2),Range(2,3,3)) =\n"
122 << C(Range(3, 4, 2), Range(2, 3, 3)) << "\n";
123
124 cout << "\nC(Range(3,4,2),Range(2,3,3)).transpose() =\n"
125 << transpose(C(Range(3, 4, 2), Range(2, 3, 3))) << "\n";
126
127 return 0;
128}
129
130void test2() {
131 Vector v(50000000);
132
133 cout << "v.nelem() = " << v.nelem() << "\n";
134
135 cout << "Filling\n";
136 // for (Index i=0; i<v.nelem(); ++i )
137 // v[i] = sqrt(i);
138 v = 1.;
139 cout << "Done\n";
140}
141
142void test4() {
143 Vector a(10);
144 Vector b(a.nelem());
145
146 for (Index i = 0; i < a.nelem(); ++i) {
147 a[i] = (Numeric)(i + 1);
148 b[i] = (Numeric)(a.nelem() - i);
149 }
150
151 cout << "a = \n" << a << "\n";
152 cout << "b = \n" << b << "\n";
153 cout << "a*b \n= " << a * b << "\n";
154
155 Matrix A(11, 6);
156 Matrix B(10, 20);
157 Matrix C(20, 5);
158
159 B = 2;
160 C = 3;
161 mult(A(Range(1, joker), Range(1, joker)), B, C);
162
163 // cout << "\nB =\n" << B << "\n";
164 // cout << "\nC =\n" << C << "\n";
165 cout << "\nB*C =\n" << A << "\n";
166}
167
168void test5() {
169 Vector a(10);
170 Vector b(20);
171 Matrix M(10, 20);
172
173 // Fill b and M with a constant number:
174 b = 1;
175 M = 2;
176
177 cout << "b = \n" << b << "\n";
178 cout << "M =\n" << M << "\n";
179
180 mult(a, M, b); // a = M*b
181 cout << "\na = M*b = \n" << a << "\n";
182
183 mult(transpose((MatrixView)b), transpose((MatrixView)a), M); // b^t = a^t * M
184 cout << "\nb^t = a^t * M = \n" << transpose((MatrixView)b) << "\n";
185}
186
187void test6() {
188 Index n = 5000;
189 Vector x(1, n, 1), y(n);
190 Matrix M(n, n);
191 M = 1;
192 // cout << "x = \n" << x << "\n";
193
194 cout << "Transforming.\n";
195 // transform(x,sin,x);
196 // transform(transpose(y),sin,transpose(x));
197 // cout << "sin(x) =\n" << y << "\n";
198 for (Index i = 0; i < 1000; ++i) {
199 // mult(y,M,x);
200 transform(y, sin, static_cast<MatrixView>(x));
201 x += 1;
202 }
203 // cout << "y =\n" << y << "\n";
204
205 cout << "Done.\n";
206}
207
208void test7() {
209 Vector x(1, 20000000, 1);
210 Vector y(x.nelem());
211 transform(y, sin, x);
212 cout << "min(sin(x)), max(sin(x)) = " << min(y) << ", " << max(y) << "\n";
213}
214
215void test8() {
216 Vector x(80000000);
217 for (Index i = 0; i < x.nelem(); ++i) x[i] = (Numeric)i;
218 cout << "Done."
219 << "\n";
220}
221
222void test9() {
223 // Initialization of Matrix with view of other Matrix:
224 Matrix A(4, 8);
225 Matrix B(A(Range(joker), Range(0, 3)));
226 cout << "B = " << B << "\n";
227}
228
229void test10() {
230 // Initialization of Matrix with a vector (giving a 1 column Matrix).
231
232 // At the moment doing this with a non-const Vector will result in a
233 // warning message.
234 Vector v(1, 8, 1);
235 Matrix M((const Vector)v);
236 cout << "M = " << M << "\n";
237}
238
239void test11() {
240 // Assignment between Vector and Matrix:
241
242 // At the moment doing this with a non-const Vector will result in a
243 // warning message.
244 Vector v(1, 8, 1);
245 Matrix M(v.nelem(), 1);
246 M = v;
247 cout << "M = " << M << "\n";
248}
249
250void test12() {
251 // Copying of Arrays
252
253 Array<String> sa(3);
254 sa[0] = "It's ";
255 sa[1] = "a ";
256 sa[2] = "test.";
257
258 Array<String> sb(sa), sc(sa.nelem());
259
260 cout << "sb = \n" << sb << "\n";
261
262 sc = sa;
263
264 cout << "sc = \n" << sc << "\n";
265}
266
267void test13() {
268 // Mix vector and one-column matrix in += operator.
269 const Vector v(1, 8, 1); // The const is necessary here to
270 // avoid compiler warnings about
271 // different conversion paths.
272 Matrix M(v);
273 M += v;
274 cout << "M = \n" << M << "\n";
275}
276
277void test14() {
278 // Test explicit Array constructors.
279 Array<String> a{"Test"};
280 Array<Index> b{1, 2};
281 Array<Numeric> c{1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0};
282 cout << "a = \n" << a << "\n";
283 cout << "b = \n" << b << "\n";
284 cout << "c = \n" << c << "\n";
285}
286
287void test15() {
288 // Test String.
289 String a = "Nur ein Test.";
290 cout << "a = " << a << "\n";
291 String b(a, 5, -1);
292 cout << "b = " << b << "\n";
293}
294
295/*void test16()
296{
297 // Test interaction between Array<Numeric> and Vector.
298 Vector a;
299 Array<Numeric> b;
300 b.push_back(1);
301 b.push_back(2);
302 b.push_back(3);
303 a.resize(b.nelem());
304 a = b;
305 cout << "b =\n" << b << "\n";
306 cout << "a =\n" << a << "\n";
307}*/
308
309void test17() {
310 // Test Sum.
311 Vector a(1, 10, 1);
312 cout << "a.sum() = " << a.sum() << "\n";
313}
314
315void test18() {
316 // Test elementvise square of a vector:
317 Vector a(1, 10, 1);
318 a *= a;
319 cout << "a *= a =\n" << a << "\n";
320}
321
322void test19() {
323 // There exists no explicit filling constructor of the form
324 // Vector a(3,1.7).
325 // But you can use the more general filling constructor with 3 arguments.
326
327 Vector a(1, 10, 1);
328 Vector b(5.3, 10, 0);
329 cout << "a =\n" << a << "\n";
330 cout << "b =\n" << b << "\n";
331}
332
333void test20() {
334 // Test initialization list constructor:
335 Vector a{1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
336 cout << "a =\n" << a << "\n";
337}
338
339void test21() {
340 Numeric s = 0;
341 // Test speed of call by reference:
342 cout << "By reference:\n";
343 for (Index i = 0; i < (Index)1e8; ++i) {
344 s += by_reference(s);
345 s -= by_reference(s);
346 }
347 cout << "s = " << s << "\n";
348}
349
350void test22() {
351 Numeric s = 0;
352 // Test speed of call by value:
353 cout << "By value:\n";
354 for (Index i = 0; i < (Index)1e8; ++i) {
355 s += by_value(s);
356 s -= by_value(s);
357 }
358 cout << "s = " << s << "\n";
359}
360
361void test23() {
362 // Test constructors that fills with constant:
363 Vector a(10, 3.5);
364 cout << "a =\n" << a << "\n";
365 Matrix b(10, 10, 4.5);
366 cout << "b =\n" << b << "\n";
367}
368
369void test24() {
370 // Try element-vise multiplication of Matrix and Vector:
371 Matrix a(5, 1, 2.5);
372 Vector b(1, 5, 1);
373 a *= b;
374 cout << "a*=b =\n" << a << "\n";
375 a /= b;
376 cout << "a/=b =\n" << a << "\n";
377 a += b;
378 cout << "a+=b =\n" << a << "\n";
379 a -= b;
380 cout << "a-=b =\n" << a << "\n";
381}
382
383void test25() {
384 // Test min and max for Array:
385 Array<Index> a{1, 2, 3, 4, 5, 6, 5, 4, 3, 2, 1};
386 cout << "min/max of a = " << min(a) << "/" << max(a) << "\n";
387}
388
389void test26() {
390 cout << "Test filling constructor for Array:\n";
391 Array<String> a(4, "Hello");
392 cout << "a =\n" << a << "\n";
393}
394
395void test27() {
396 cout << "Test Arrays of Vectors:\n";
398 a.push_back({1.0, 2.0});
399 a.push_back(Vector(1.0, 10, 1.0));
400 cout << "a =\n" << a << "\n";
401}
402
403void test28() {
404 cout << "Test default constructor for Matrix:\n";
405 Matrix a;
406 Matrix b(a);
407 cout << "b =\n" << b << "\n";
408}
409
410void test29() {
411 cout << "Test Arrays of Matrix:\n";
413 Matrix b;
414
415 b.resize(2, 2);
416 b(0, 0) = 1;
417 b(0, 1) = 2;
418 b(1, 0) = 3;
419 b(1, 1) = 4;
420 a.push_back(b);
421 b *= 2;
422 a.push_back(b);
423
424 a[0].resize(2, 3);
425 a[0] = 4;
426
427 a.resize(3);
428 a[2].resize(4, 5);
429 a[2] = 5;
430
431 cout << "a =\n" << a << "\n";
432}
433
434void test30() {
435 cout << "Test Matrices of size 0:\n";
436 Matrix a(0, 0);
437 // cout << "a(0,0) =\n" << a(0,0) << "\n";
438 a.resize(2, 2);
439 a = 1;
440 cout << "a =\n" << a << "\n";
441
442 Matrix b(3, 0);
443 // cout << "b(0,0) =\n" << b(0,0) << "\n";
444 b.resize(b.nrows(), b.ncols() + 3);
445 b = 2;
446 cout << "b =\n" << b << "\n";
447
448 Matrix c(0, 3);
449 // cout << "c(0,0) =\n" << c(0,0) << "\n";
450 c.resize(c.nrows() + 3, c.ncols());
451 c = 3;
452 cout << "c =\n" << c << "\n";
453}
454
455void test31() {
456 cout << "Test Tensor3:\n";
457
458 Tensor3 a(2, 3, 4, 1.0);
459
460 Index fill = 0;
461
462 // Fill with some numbers
463 for (Index i = 0; i < a.npages(); ++i)
464 for (Index j = 0; j < a.nrows(); ++j)
465 for (Index k = 0; k < a.ncols(); ++k) a(i, j, k) = (Numeric)(++fill);
466
467 cout << "a =\n" << a << "\n";
468
469 cout << "Taking out first row of first page:\n"
470 << a(0, 0, Range(joker)) << "\n";
471
472 cout << "Taking out last column of second page:\n"
473 << a(1, Range(joker), a.ncols() - 1) << "\n";
474
475 cout << "Taking out the first letter on every page:\n"
476 << a(Range(joker), 0, 0) << "\n";
477
478 cout << "Taking out first page:\n"
479 << a(0, Range(joker), Range(joker)) << "\n";
480
481 cout << "Taking out last row of all pages:\n"
482 << a(Range(joker), a.nrows() - 1, Range(joker)) << "\n";
483
484 cout << "Taking out second column of all pages:\n"
485 << a(Range(joker), Range(joker), 1) << "\n";
486
487 a *= 2;
488
489 cout << "After element-vise multiplication with 2:\n" << a << "\n";
490
491 transform(a, sqrt, a);
492
493 cout << "After taking the square-root:\n" << a << "\n";
494
495 Index s = 200;
496 cout << "Let's allocate a large tensor, "
497 << (Numeric)(s * s * s * 8) / 1024. / 1024. << " MB...\n";
498
499 a.resize(s, s, s);
500
501 cout << "Set it to 1...\n";
502
503 a = 1;
504
505 cout << "a(900,900,900) = " << a(90, 90, 90) << "\n";
506
507 fill = 0;
508
509 cout << "Fill with running numbers, using for loops...\n";
510 for (Index i = 0; i < a.npages(); ++i)
511 for (Index j = 0; j < a.nrows(); ++j)
512 for (Index k = 0; k < a.ncols(); ++k) a(i, j, k) = (Numeric)(++fill);
513
514 cout << "Max(a) = ...\n";
515
516 cout << max(a) << "\n";
517}
518
519void test32() {
520 cout << "Test von X = A*X:\n";
521 Matrix X(3, 3), A(3, 3), B(3, 3);
522
523 for (Index j = 0; j < A.nrows(); ++j)
524 for (Index k = 0; k < A.ncols(); ++k) {
525 X(j, k) = 1;
526 A(j, k) = (Numeric)(j + k);
527 }
528 cout << "A:\n" << A << "\n";
529 cout << "X:\n" << X << "\n";
530
531 mult(B, A, X);
532 cout << "B = A*X:\n" << B << "\n";
533 mult(X, A, X);
534 cout << "X = A*X:\n" << X << "\n";
535
536 cout << "This is not the same, and should not be, because you\n"
537 << "are not allowed to use the same matrix as input and output\n"
538 << "for mult!\n";
539}
540
541void test33() {
542 cout << "Making things look bigger than they are...\n";
543
544 {
545 cout << "1. Make a scalar look like a vector:\n";
546 Numeric a = 3.1415; // Just any number here.
547 VectorView av(a);
548 cout << "a, viewed as a vector: " << av << "\n";
549 cout << "Describe a: " << describe(a) << "\n";
550 av[0] += 1;
551 cout << "a, after the first element\n"
552 << "of the vector has been increased by 1: " << a << "\n";
553 }
554
555 {
556 cout
557 << "\n2. Make a vector look like a matrix:\n"
558 << "This is an exception, because the new dimension is added at the end.\n";
559 Vector a{1, 2, 3, 4, 5};
560 MatrixView am = a;
561 cout << "a, viewed as a matrix:\n" << am << "\n";
562 cout << "Trasnpose view:\n" << transpose(am) << "\n";
563 cout << "Describe transpose(am): " << describe(transpose(am)) << "\n";
564
565 cout << "\n3. Make a matrix look like a tensor3:\n";
566 Tensor3View at3 = am;
567 cout << "at3 = \n" << at3 << "\n";
568 cout << "Describe at3: " << describe(at3) << "\n";
569 at3(0, 2, 0) += 1;
570 cout << "a after Increasing element at3(0,2,0) by 1: \n" << a << "\n\n";
571
572 Tensor4View at4 = at3;
573 cout << "at4 = \n" << at4 << "\n";
574 cout << "Describe at4: " << describe(at4) << "\n";
575
576 Tensor5View at5 = at4;
577 cout << "at5 = \n" << at5 << "\n";
578 cout << "Describe at5: " << describe(at5) << "\n";
579
580 Tensor6View at6 = at5;
581 cout << "at6 = \n" << at6 << "\n";
582 cout << "Describe at6: " << describe(at6) << "\n";
583
584 Tensor7View at7 = at6;
585 cout << "at7 = \n" << at7 << "\n";
586 cout << "Describe at7: " << describe(at7) << "\n";
587
588 at7(0, 0, 0, 0, 0, 2, 0) -= 1;
589
590 cout << "After subtracting 1 from at7(0,0,0,0,0,2,0)\n"
591 << "a = " << a << "\n";
592
593 cout << "\nAll in one go:\n";
594 Numeric b = 3.1415; // Just any number here.
597 cout << "bt7:\n" << bt7 << "\n";
598 cout << "Describe bt7: " << describe(bt7) << "\n";
599 }
600}
601
602void junk4(Tensor4View a) { cout << "Describe a: " << describe(a) << "\n"; }
603
604void junk2(ConstVectorView a) { cout << "Describe a: " << describe(a) << "\n"; }
605
606void test34() {
607 cout << "Test, if dimension expansion works implicitly.\n";
608
609 Tensor3 t3(2, 3, 4);
610 junk4(t3);
611
612 Numeric x;
614}
615
616void test35() {
617 cout << "Test the new copy semantics.\n";
618 Vector a(1, 4, 1);
619 Vector b;
620
621 b = a;
622 cout << "b = " << b << "\n";
623
624 Vector aa(1, 5, 1);
625 ConstVectorView c = aa;
626 b = c;
627 cout << "b = " << b << "\n";
628
629 Vector aaa(1, 6, 1);
630 VectorView d = aaa;
631 b = d;
632 cout << "b = " << b << "\n";
633}
634
635void test36() {
636 cout << "Test using naked joker on Vector.\n";
637 Vector a(1, 4, 1);
638 VectorView b = a[joker];
639 cout << "a = " << a << "\n";
640 cout << "b = " << b << "\n";
641}
642
643void test37(const Index& i) {
644 Vector v1(5e-15, 10, 0.42e-15 / 11);
645 Vector v2 = v1;
646 // Index i = 10000000;
647
648 v1 /= (Numeric)i;
649 v2 /= (Numeric)i;
650 cout.precision(12);
651 // cout.setf(ios_base::scientific,ios_base::floatfield);
652 v1 *= v1;
653 v2 *= v2;
654 cout << v1 << endl;
655 cout << v2 << endl;
656}
657
658void test38() {
659 Vector v(5, 0.);
660 Numeric* const a = v.get_c_array();
661
662 a[4] = 5.;
663
664 cout << v << endl;
665 cout << endl << "========================" << endl << endl;
666
667 Matrix m(5, 5, 0.);
668 Numeric* const b = m.get_c_array();
669
670 b[4] = 5.;
671
672 cout << m << endl;
673 cout << endl << "========================" << endl << endl;
674
675 Tensor3 t3(5, 6, 7, 0.);
676 Numeric* const c = t3.get_c_array();
677
678 c[6] = 5.;
679
680 cout << t3 << endl;
681}
682
683void test39() {
684 Vector v1(1, 5, 1), v2(5);
685
686 v2 = v1 * 2;
687 // Unfortunately, this thing compiles, but at least it gives an
688 // assertion failure at runtime. I think what happens is that it
689 // implicitly creates a one element vector out of the "2", then
690 // tries to do a scalar product with v1.
691
692 cout << v2 << endl;
693}
694
695void test40() {
696 Vector v(100);
697
698 v = 5;
699
700 cout << v << endl;
701}
702
703void test41() {
704 const Vector v1(10, 1);
705
706 ConstVectorView vv = v1[Range(0, 5)];
707
708 cout << "Vector: " << v1 << endl;
709 cout << "VectorView: " << vv << endl;
710
711 try {
712 vv = 2;
713 } catch (const std::runtime_error& e) {
714 std::cerr << e.what() << endl;
715 exit(EXIT_FAILURE);
716 }
717
718 cout << "VectorView: " << vv << endl;
719 cout << "Vector: " << v1 << endl;
720}
721
722// Test behaviour of VectorView::operator MatrixView, for which I have fixed
723// a bug today.
724// SAB 2013-01-18
725void test42() {
726 cout << "test42\n";
727 Vector x{1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
728 cout << "x: " << x << endl;
729
730 VectorView y = x[Range(2, 4, 2)];
731 cout << "y: " << y << endl;
732
734 cout << "z:\n" << z << endl;
735
736 cout << "Every other z:\n" << z(Range(1, 2, 2), joker) << endl;
737
738 // Try write access also:
739 MatrixView zz(y);
740 zz(Range(1, 2, 2), joker) = 0;
741 cout << "zz:\n" << zz << endl;
742 cout << "New x: " << x << endl;
743}
744
745
746// Test function for internal use, must return 2 as last n2r
747constexpr Rational test_numeric2rational(const Index i, const Index maxi, const Rational r=0, const Rational n2r=0) {
748 if (i > maxi)
749 return n2r;
750 else {
751 return
752 (r == n2r) ?
753 test_numeric2rational(i+1, maxi, Rational(maxi + i, maxi), numeric2rational(1 + Numeric(i)/Numeric(maxi), 12)) :
754 throw std::logic_error("Fail to initialize");;
755 }
756}
757
758
759void test43() {
760 // Simple construction compile-time test
761 constexpr Rational r=3_2; // should be 3/2
762 static_assert(r.Nom() == 3, "Rational fail to initialize properly");
763 static_assert(r.Denom() == 2, "Rational fail to initialize properly");
764
765 // Simple expression compile-time test
766 constexpr Rational r2 = 1 + r; // should be 5/2
767 static_assert(r2.Nom() == 5, "Rational fail to initialize properly");
768 static_assert(r2.Denom() == 2, "Rational fail to initialize properly");
769
770 // Div-zero by index generates an undefined rational: 0/0
771 constexpr Rational r3 = r / 0; // should be undefined
772 static_assert(r3.Nom() == 0, "Rational fail to initialize properly");
773 static_assert(r3.Denom() == 0, "Rational fail to initialize properly");
774 static_assert(r3 not_eq r3, "Undefined rational does not equal self");
775 static_assert(r3.isUndefined(), "Rational fail to initialize properly");
776
777 // Mul-zero gives 0/1
778 constexpr Rational r4 = r * 0; // should be undefined
779 static_assert(r4.Nom() == 0, "Rational fail to initialize properly");
780 static_assert(r4.Denom() == 1, "Rational fail to initialize properly");
781
782 // Complicated expression compile-time test
783 constexpr Rational r5 = (Rational(1, 9) * r) % 3; // should be 1/6
784 static_assert(r5.Nom() == 1, "Rational fail to initialize properly");
785 static_assert(r5.Denom() == 6, "Rational fail to initialize properly");
786
787 // The simplify operation simplifies expressions
788 constexpr Rational r6 = 10 % r; // should be 1
789 static_assert(r6.Nom() == 1, "Rational fail to initialize properly");
790 static_assert(r6.Denom() == 1, "Rational fail to initialize properly");
791 static_assert(r6.toInt() == 1, "Rational fail to initialize properly");
792 static_assert(r6.toIndex() == 1, "Rational fail to initialize properly");
793 static_assert(r6.toNumeric() == 1e0, "Rational fail to initialize properly");
794
795 constexpr Index rattest=1<<8;
796 constexpr Rational r7 = test_numeric2rational(0, rattest);
797 static_assert(r7 == 2, "Rational fail to initialize properly");
798}
799
800void test44() {
801#define docheck(fn, val, expect) \
802 cout << #fn << "(" << val << ") = " << fn(x) << " (expected " << #expect \
803 << ")" << std::endl;
804
805 Vector x{1, 2, 3};
807 docheck(is_sorted, x, 1)
808
809 x = {3, 2, 1};
811 docheck(is_sorted, x, 0)
812
813 x = {1, 2, 2};
815 docheck(is_sorted, x, 1)
816
817 x = {2, 2, 1};
819 docheck(is_sorted, x, 0)
820
821 x = {1, NAN, 2};
823 docheck(is_sorted, x, 0)
824
825 x = {2, NAN, 1};
827 docheck(is_sorted, x, 0)
828
829 x = {NAN, NAN, NAN};
831 docheck(is_sorted, x, 0)
832
833#undef docheck
834}
835
836void test46() {
837 Vector v(5, 0.);
838 nlinspace(v, 1, 10, 10);
839 VectorView v1 = v;
840 auto compare_func = [](Numeric n) { return n != 0; };
841 cout << std::any_of(v1.begin(), v1.end(), compare_func) << endl;
842 v1 = 0.;
843 cout << std::any_of(v1.begin(), v1.end(), compare_func) << endl;
844}
845
847
855bool test_diagonal(Index ntests) {
856 Rand<Index> rand(1, 100);
857
858 Index n_diag;
859
860 bool pass = true;
861
862 for (Index i = 0; i < ntests; i++) {
863 Matrix A(rand(), rand());
864 n_diag = std::min(A.ncols(), A.nrows());
865
867 ConstMatrixView B = A;
868 if (n_diag > 3) {
869 B = A(Range(1, Joker(), 2), Range(1, Joker(), 2));
870 }
871
873 ConstVectorView vt(AT.diagonal());
875
876 random_fill_matrix(A, 10, true);
877
878 for (Index j = 0; j < n_diag; j++) {
879 pass = pass && (v[j] == A(j, j));
880 pass = pass && (vt[j] == AT(j, j));
881
882 if (j < vb.nelem()) pass = pass && (vb[j] == B(j, j));
883 }
884 cout << endl;
885 }
886
887 if (pass)
888 cout << "test diagonal: Passed all tests." << endl;
889 else
890 cout << "test diagonal: Failed." << endl;
891
892 return pass;
893}
894
895Numeric matrix_vector_mult(Index m, Index n, Index ntests, bool verbose) {
896 Numeric max_err = 0;
897 Matrix A(m, n);
898 Vector x(n), xt, x_ref(n), y(m), y_ref(m);
899
900 Rand<Index> random_row_stride(1, n / 4);
901 Rand<Index> random_column_stride(1, m / 4);
902
903 for (Index i = 0; i < ntests; i++) {
904 // A * x
905
906 random_fill_matrix(A, 1000, false);
907 random_fill_vector(x, 1000, false);
908
909 mult(y, A, x);
910 mult_general(y_ref, A, x);
911
912 Numeric err_mul = get_maximum_error(y, y_ref, true);
913 if (err_mul > max_err) max_err = err_mul;
914
915 if (verbose) {
916 cout << "\t A * x: max. rel. error = " << err_mul << endl;
917 }
918
919 // A^T * y
920
921 mult(x, transpose(A), y);
922 mult_general(x_ref, transpose(A), y);
923
924 err_mul = get_maximum_error(x, x_ref, true);
925 if (err_mul > max_err) max_err = err_mul;
926
927 if (verbose) {
928 cout << "\t A^T * x: max. rel. error = " << err_mul << endl;
929 }
930
931 // Random stride
932 Index column_stride(random_column_stride()),
933 row_stride(random_row_stride());
934
935 Index m_sub, n_sub;
936 m_sub = (m - 1) / row_stride + 1;
937 n_sub = (n - 1) / column_stride + 1;
938
939 mult(y[Range(0, joker, row_stride)],
940 A(Range(0, m_sub), Range(0, n_sub)),
941 x[Range(0, joker, column_stride)]);
942 mult_general(y_ref[Range(0, joker, row_stride)],
943 A(Range(0, m_sub), Range(0, n_sub)),
944 x[Range(0, joker, column_stride)]);
945
946 err_mul = get_maximum_error(y[Range(0, joker, row_stride)],
947 y_ref[Range(0, joker, row_stride)],
948 true);
949 if (err_mul > max_err) max_err = err_mul;
950
951 if (verbose) {
952 cout << "\t Random stride: max. rel. error = " << err_mul << endl << endl;
953 }
954
955 // Random offset
956 if ((m > 1) && (n > 1)) {
957 Index y_offset = rand() % (m - 1);
958 Index x_offset = rand() % (n - 1);
959
960 m_sub = m - y_offset - 1;
961 n_sub = n - x_offset - 1;
962
963 mult(y[Range(y_offset, m_sub)],
964 A(Range(y_offset, m_sub), Range(x_offset, n_sub)),
965 x[Range(x_offset, n_sub)]);
966 mult_general(y_ref[Range(y_offset, m_sub)],
967 A(Range(y_offset, m_sub), Range(x_offset, n_sub)),
968 x[Range(x_offset, n_sub)]);
969
970 err_mul = get_maximum_error(
971 y[Range(y_offset, m_sub)], y_ref[Range(y_offset, m_sub)], true);
972 if (err_mul > max_err) max_err = err_mul;
973
974 if (verbose) {
975 cout << "\t Random offset: max. rel. error = " << err_mul << endl
976 << endl;
977 }
978 }
979 }
980
981 return max_err;
982}
983
985 Numeric err, max_err;
986
987 if (verbose)
988 cout << "Matrix-Vector Multiplication: n = m = 100, ntests = 100" << endl;
989
990 max_err = matrix_vector_mult(100, 100, 100, verbose);
991 if (verbose) {
992 cout << endl;
993 cout << "Matrix-Vector Multiplication: n = 100, m = 20, ntests = 100"
994 << endl;
995 }
996
997 err = matrix_vector_mult(100, 20, 100, verbose);
998 if (err > max_err) max_err = err;
999 if (verbose) {
1000 cout << endl;
1001 cout << "Matrix-Vector Multiplication: n = 20, m = 100, ntests = 100"
1002 << endl;
1003 }
1004
1005 err = matrix_vector_mult(20, 100, 100, verbose);
1006 if (err > max_err) max_err = err;
1007 if (verbose) {
1008 if (max_err < 1e-9)
1009 cout << endl << "Matrix Vector Multiplication: PASSED" << endl;
1010 else
1011 cout << endl << "Matrix Vector Multiplication: FAILED" << endl;
1012 }
1013
1014 err = matrix_vector_mult(100, 1, 100, verbose);
1015 if (err > max_err) max_err = err;
1016 if (verbose) {
1017 cout << endl;
1018 cout << "Matrix-Vector Multiplication: n = 20, m = 100, ntests = 100"
1019 << endl;
1020 }
1021
1022 err = matrix_vector_mult(1, 100, 100, verbose);
1023 if (err > max_err) max_err = err;
1024 if (verbose) {
1025 cout << endl;
1026 cout << "Matrix-Vector Multiplication: n = 20, m = 100, ntests = 100"
1027 << endl;
1028 }
1029
1030 if (verbose) {
1031 if (max_err < 1e-9)
1032 cout << endl << "Matrix Vector Multiplication: PASSED" << endl;
1033 else
1034 cout << endl << "Matrix Vector Multiplication: FAILED" << endl;
1035 }
1036 return max_err;
1037}
1038
1040
1060 Index k, Index m, Index n, Index ntests, Index nsubtests, bool verbose) {
1061 Numeric max_err = 0;
1062 Matrix A(m, k), B(k, n), C(m, n), C_ref(m, n);
1063
1064 for (Index i = 0; i < ntests; i++) {
1065 random_fill_matrix(A, 1000, false);
1066 random_fill_matrix(B, 1000, false);
1067
1068 if (verbose) {
1069 cout.precision(15);
1070 cout << "MATRIX MULTIPLICATION: m = " << m << ", k = " << k;
1071 cout << ", n = " << n << endl;
1072 }
1073
1074 // A * B
1075
1076 random_fill_matrix(C, 10, false);
1077 random_fill_matrix(C_ref, 10, false);
1078
1079 mult(C, A, B);
1080 mult_general(C_ref, A, B);
1081
1082 Numeric err_mul = get_maximum_error(C, C_ref, true);
1083 if (err_mul > max_err) max_err = err_mul;
1084
1085 if (verbose) {
1086 cout << "\t"
1087 << "A * B: max. rel. error = " << err_mul << endl;
1088 }
1089
1090 // A^T * B
1091
1092 Matrix AT(k, m);
1093 random_fill_matrix(AT, 100.0, false);
1094
1095 mult(C, transpose(AT), B);
1096 mult_general(C_ref, transpose(AT), B);
1097 Numeric err_trans1 = get_maximum_error(C, C_ref, true);
1098 if (err_trans1 > max_err) max_err = err_trans1;
1099
1100 if (verbose) {
1101 cout << "\t"
1102 << "A^T * B: max. rel. err = " << err_trans1 << endl;
1103 }
1104
1105 // A * B^T
1106
1107 Matrix BT(n, k);
1108 random_fill_matrix(BT, 100.0, false);
1109
1110 mult(C, A, transpose(BT));
1111 mult_general(C_ref, A, transpose(BT));
1112 Numeric err_trans2 = get_maximum_error(C, C_ref, true);
1113 if (err_trans2 > max_err) max_err = err_trans2;
1114
1115 if (verbose) {
1116 cout << "\t"
1117 << "A * B^T: max. rel. err = " << err_trans2 << endl;
1118 }
1119
1120 // A^T * B^T
1121
1122 mult(C, transpose(AT), transpose(BT));
1123 mult_general(C_ref, transpose(AT), transpose(BT));
1124 Numeric err_trans3 = get_maximum_error(C, C_ref, true);
1125 if (err_trans3 > max_err) max_err = err_trans3;
1126
1127 if (verbose) {
1128 cout << "\t"
1129 << "A^T * B^T: max. rel. err = " << err_trans3 << endl;
1130 cout << endl;
1131 }
1132 }
1133
1134 // Multiplication of submatrices.
1135 Rand<Index> k_rand(1, k - 1), m_rand(1, m - 1), n_rand(1, n - 1);
1136
1137 for (Index i = 0; i < nsubtests - 1; i++) {
1138 Index k1(k_rand()), m1(m_rand()), n1(n_rand());
1139
1140 Range r1(random_range(m1));
1141 Range r2(random_range(k1));
1142 Range r3(random_range(n1));
1143
1144 MatrixView C_sub(C(r1, r3));
1145 MatrixView C_sub_ref(C_ref(r1, r3));
1146
1147 ConstMatrixView A_sub(A(r1, r2));
1148 ConstMatrixView B_sub(B(r2, r3));
1149
1150 mult(C_sub, A_sub, B_sub);
1151 mult_general(C_sub_ref, A_sub, B_sub);
1152
1153 Numeric err = get_maximum_error(C_sub, C_sub_ref, true);
1154 if (err > max_err) max_err = err;
1155
1156 if (verbose) {
1157 cout << "\t"
1158 << "Submatrix multiplication: max. rel. err = " << err;
1159 cout << endl;
1160 }
1161 }
1162
1163 if (verbose) {
1164 cout << endl;
1165 }
1166
1167 return max_err;
1168}
1169
1171
1183 Numeric max_err, err;
1184
1185 // Trivial case k = m = n = 0.
1186 max_err = matrix_mult(0, 0, 0, 10, 0, verbose);
1187
1188 // k = 1, m = 1, n = 1.
1189 err = matrix_mult(1, 1, 1, 10, 0, verbose);
1190 if (err > max_err) max_err = err;
1191
1192 // k = 10, m = 1, n = 10.
1193 err = matrix_mult(10, 1, 10, 20, 20, verbose);
1194 if (err > max_err) max_err = err;
1195
1196 // k = 10, m = 1, n = 10.
1197 err = matrix_mult(10, 1, 10, 20, 20, verbose);
1198 if (err > max_err) max_err = err;
1199
1200 // k = 10, m = 1, n = 1.
1201 err = matrix_mult(10, 1, 1, 20, 20, verbose);
1202 if (err > max_err) max_err = err;
1203
1204 // k = 100, m = 100, n = 100.
1205 err = matrix_mult(200, 200, 200, 20, 20, verbose);
1206 if (err > max_err) max_err = err;
1207
1208 // k = 10, m = 100, n = 10.
1209 err = matrix_mult(10, 100, 10, 20, 20, verbose);
1210 if (err > max_err) max_err = err;
1211
1212 // k = 10, m = 100, n = 100.
1213 err = matrix_mult(10, 100, 100, 20, 20, verbose);
1214 if (err > max_err) max_err = err;
1215
1216 // k = 100, m = 100, n = 100, 100 submatrix multiplications.
1217 err = matrix_mult(100, 100, 100, 0, 100, verbose);
1218 if (err > max_err) max_err = err;
1219
1220 return max_err;
1221}
1222
1225 {
1226 cout << "Array" << endl;
1228 v.resize(0);
1229 cout << v.empty() << endl;
1230 v.resize(1);
1231 cout << v.empty() << endl;
1232 }
1233 {
1234 cout << "Vector" << endl;
1235 Vector v;
1236 v.resize(0);
1237 cout << v.empty() << endl;
1238 v.resize(1);
1239 cout << v.empty() << endl;
1240 }
1241 {
1242 cout << "Matrix" << endl;
1243 Matrix v;
1244 v.resize(0, 1);
1245 cout << v.empty() << endl;
1246 v.resize(1, 0);
1247 cout << v.empty() << endl;
1248 v.resize(1, 1);
1249 cout << v.empty() << endl;
1250 }
1251 {
1252 cout << "Tensor3" << endl;
1253 Tensor3 v;
1254 v.resize(0, 1, 1);
1255 cout << v.empty() << endl;
1256 v.resize(1, 0, 1);
1257 cout << v.empty() << endl;
1258 v.resize(1, 1, 0);
1259 cout << v.empty() << endl;
1260 v.resize(1, 1, 1);
1261 cout << v.empty() << endl;
1262 }
1263 {
1264 cout << "Tensor4" << endl;
1265 Tensor4 v;
1266 v.resize(0, 1, 1, 1);
1267 cout << v.empty() << endl;
1268 v.resize(1, 0, 1, 1);
1269 cout << v.empty() << endl;
1270 v.resize(1, 1, 0, 1);
1271 cout << v.empty() << endl;
1272 v.resize(1, 1, 1, 0);
1273 cout << v.empty() << endl;
1274 v.resize(1, 1, 1, 1);
1275 cout << v.empty() << endl;
1276 }
1277 {
1278 cout << "Tensor5" << endl;
1279 Tensor5 v;
1280 v.resize(0, 1, 1, 1, 1);
1281 cout << v.empty() << endl;
1282 v.resize(1, 0, 1, 1, 1);
1283 cout << v.empty() << endl;
1284 v.resize(1, 1, 0, 1, 1);
1285 cout << v.empty() << endl;
1286 v.resize(1, 1, 1, 0, 1);
1287 cout << v.empty() << endl;
1288 v.resize(1, 1, 1, 1, 0);
1289 cout << v.empty() << endl;
1290 v.resize(1, 1, 1, 1, 1);
1291 cout << v.empty() << endl;
1292 }
1293 {
1294 cout << "Tensor6" << endl;
1295 Tensor6 v;
1296 v.resize(0, 1, 1, 1, 1, 1);
1297 cout << v.empty() << endl;
1298 v.resize(1, 0, 1, 1, 1, 1);
1299 cout << v.empty() << endl;
1300 v.resize(1, 1, 0, 1, 1, 1);
1301 cout << v.empty() << endl;
1302 v.resize(1, 1, 1, 0, 1, 1);
1303 cout << v.empty() << endl;
1304 v.resize(1, 1, 1, 1, 0, 1);
1305 cout << v.empty() << endl;
1306 v.resize(1, 1, 1, 1, 1, 0);
1307 cout << v.empty() << endl;
1308 v.resize(1, 1, 1, 1, 1, 1);
1309 cout << v.empty() << endl;
1310 }
1311 {
1312 cout << "Tensor7" << endl;
1313 Tensor7 v;
1314 v.resize(0, 1, 1, 1, 1, 1, 1);
1315 cout << v.empty() << endl;
1316 v.resize(1, 0, 1, 1, 1, 1, 1);
1317 cout << v.empty() << endl;
1318 v.resize(1, 1, 0, 1, 1, 1, 1);
1319 cout << v.empty() << endl;
1320 v.resize(1, 1, 1, 0, 1, 1, 1);
1321 cout << v.empty() << endl;
1322 v.resize(1, 1, 1, 1, 0, 1, 1);
1323 cout << v.empty() << endl;
1324 v.resize(1, 1, 1, 1, 1, 0, 1);
1325 cout << v.empty() << endl;
1326 v.resize(1, 1, 1, 1, 1, 1, 0);
1327 cout << v.empty() << endl;
1328 v.resize(1, 1, 1, 1, 1, 1, 1);
1329 cout << v.empty() << endl;
1330 }
1331}
1332
1334 const Numeric start,
1335 const Numeric stop,
1336 const Index n) {
1337 ARTS_ASSERT(1 < n); // Number of points must be greater 1.
1338 x.resize(n);
1339 Numeric step = (stop - start) / ((double)n - 1);
1340 for (Index i = 0; i < n - 1; i++) x[i] = start + (double)i * step;
1341 x[n - 1] = stop;
1342}
1343
1344void test47() {
1345 // Selecting empty matpack dimensions with Joker shouldn't fail
1346 Vector v;
1347 Matrix m;
1348
1349 VectorView vv = v[joker];
1350 MatrixView mv = m(joker, joker);
1351
1352 std::cout << "vv.nelem: " << vv.nelem() << std::endl;
1353 std::cout << "mv.nrows: " << mv.nrows() << " ";
1354 std::cout << "mv.ncols: " << mv.ncols() << std::endl;
1355
1356 Matrix m2(0, 5);
1357 MatrixView mv2 = m2(joker, 3);
1358 std::cout << "mv2.nrows: " << mv2.nrows() << " ";
1359 std::cout << "mv2.ncols: " << mv2.ncols() << std::endl;
1360}
1361
1362void test48() {
1363
1364 // Test each simple reduction
1365 std::cout << Tensor7(2, 2, 2, 2, 2, 2, 1, 8).reduce_rank<0, 1, 2, 3, 4, 5>() << '\n';
1366 std::cout << Tensor7(2, 2, 2, 2, 2, 1, 1, 9).reduce_rank<0, 1, 2, 3, 4>() << '\n';
1367 std::cout << Tensor7(2, 2, 2, 2, 1, 1, 1, 10).reduce_rank<0, 1, 2, 3>() << '\n';
1368 std::cout << Tensor7(2, 2, 2, 1, 1, 1, 1, 11).reduce_rank<0, 1, 2>() << '\n';
1369 std::cout << Tensor7(2, 2, 1, 1, 1, 1, 1, 12).reduce_rank<0, 1>() << '\n';
1370 std::cout << Tensor7(2, 1, 1, 1, 1, 1, 1, 13).reduce_rank<0>() << '\n';
1371 std::cout << Tensor7(1, 1, 1, 1, 1, 1, 1, 14) << '\n';
1372 std::cout << Tensor6(2, 2, 2, 2, 2, 1, 15).reduce_rank<0, 1, 2, 3, 4>() << '\n';
1373 std::cout << Tensor6(2, 2, 2, 2, 1, 1, 16).reduce_rank<0, 1, 2, 3>() << '\n';
1374 std::cout << Tensor6(2, 2, 2, 1, 1, 1, 17).reduce_rank<0, 1, 2>() << '\n';
1375 std::cout << Tensor6(2, 2, 1, 1, 1, 1, 18).reduce_rank<0, 1>() << '\n';
1376 std::cout << Tensor6(2, 1, 1, 1, 1, 1, 19).reduce_rank<0>() << '\n';
1377 std::cout << Tensor6(1, 1, 1, 1, 1, 1, 20) << '\n';
1378 std::cout << Tensor5(2, 2, 2, 2, 1, 21).reduce_rank<0, 1, 2, 3>() << '\n';
1379 std::cout << Tensor5(2, 2, 2, 1, 1, 22).reduce_rank<0, 1, 2>() << '\n';
1380 std::cout << Tensor5(2, 2, 1, 1, 1, 23).reduce_rank<0, 1>() << '\n';
1381 std::cout << Tensor5(2, 1, 1, 1, 1, 24).reduce_rank<0>() << '\n';
1382 std::cout << Tensor5(1, 1, 1, 1, 1, 25) << '\n';
1383 std::cout << Tensor4(2, 2, 2, 1, 26).reduce_rank<0, 1, 2>() << '\n';
1384 std::cout << Tensor4(2, 2, 1, 1, 27).reduce_rank<0, 1>() << '\n';
1385 std::cout << Tensor4(2, 1, 1, 1, 28).reduce_rank<0>() << '\n';
1386 std::cout << Tensor4(1, 1, 1, 1, 29) << '\n';
1387 std::cout << Tensor3(2, 2, 1, 30).reduce_rank<0, 1>() << '\n';
1388 std::cout << Tensor3(2, 1, 1, 31).reduce_rank<0>() << '\n';
1389 std::cout << Tensor3(1, 1, 1, 32) << '\n';
1390 std::cout << Matrix(2, 1, 33).reduce_rank<0>() << '\n';
1391 std::cout << Matrix(1, 1, 34) << '\n';
1392
1393 // Test that the reductions work along different axis for Vector
1394 std::cout << MapToEigen(Tensor7(2, 1, 1, 1, 1, 1, 1, 35).reduce_rank<0, 1>()).transpose() << '\n'
1395 << MapToEigen(Tensor7(2, 1, 1, 1, 1, 1, 1, 35).reduce_rank<0, 2>()).transpose() << '\n'
1396 << MapToEigen(Tensor7(2, 1, 1, 1, 1, 1, 1, 35).reduce_rank<0, 3>()).transpose() << '\n'
1397 << MapToEigen(Tensor7(2, 1, 1, 1, 1, 1, 1, 35).reduce_rank<0, 4>()).transpose() << '\n'
1398 << MapToEigen(Tensor7(2, 1, 1, 1, 1, 1, 1, 35).reduce_rank<0, 5>()).transpose() << '\n'
1399 << MapToEigen(Tensor7(2, 1, 1, 1, 1, 1, 1, 35).reduce_rank<0, 6>()).transpose() << '\n';
1400
1401 // Test that the reductions work along different axis for Matrix
1402 std::cout << MapToEigen(Tensor7(2, 2, 1, 1, 1, 1, 1, 36).reduce_rank<0, 1>()) << '\n'
1403 << MapToEigen(Tensor7(2, 1, 2, 1, 1, 1, 1, 36).reduce_rank<0, 2>()) << '\n'
1404 << MapToEigen(Tensor7(2, 1, 1, 2, 1, 1, 1, 36).reduce_rank<0, 3>()) << '\n'
1405 << MapToEigen(Tensor7(2, 1, 1, 1, 2, 1, 1, 36).reduce_rank<0, 4>()) << '\n'
1406 << MapToEigen(Tensor7(2, 1, 1, 1, 1, 2, 1, 36).reduce_rank<0, 5>()) << '\n'
1407 << MapToEigen(Tensor7(2, 1, 1, 1, 1, 1, 2, 36).reduce_rank<0, 6>()) << '\n';
1408
1409 // Test that the right elements are accessed
1410 {
1411 Index val=1;
1412 Tensor7 testvar(9, 2, 4, 3, 5, 7, 11);for (Index i=0; i<9; i++)
1413 for (Index j=0; j<2; j++) {
1414 for (Index k=0; k<4; k++) {
1415 for (Index l=0; l<3; l++) {
1416 for (Index m=0; m<5; m++) {
1417 for (Index n=0; n<7; n++) {
1418 for (Index o=0; o<11; o++) {
1419 testvar(i, j, k, l, m, n, o) = Numeric(val++);
1420 }
1421 }
1422 }
1423 }
1424 }
1425 }
1426
1427 // Vector test
1428 for (Index j=0; j<2; j++) {
1429 for (Index k=0; k<4; k++) {
1430 for (Index l=0; l<3; l++) {
1431 for (Index m=0; m<5; m++) {
1432 for (Index n=0; n<7; n++) {
1433 for (Index o=0; o<11; o++) {
1434 std::cout << Tensor7(testvar(joker, Range(j, 1), Range(k, 1), Range(l, 1), Range(m, 1), Range(n, 1), Range(o, 1))).reduce_rank<0>() << '\n';
1435 std::cout << testvar(joker, j, k, l, m, n, o) << '\n';
1436 std::cout << '\n';
1437 }
1438 }
1439 }
1440 }
1441 }
1442 }
1443 for (Index i=0; i<9; i++) {
1444 for (Index k=0; k<4; k++) {
1445 for (Index l=0; l<3; l++) {
1446 for (Index m=0; m<5; m++) {
1447 for (Index n=0; n<7; n++) {
1448 for (Index o=0; o<11; o++) {
1449 std::cout << Tensor7(testvar(Range(i, 1), joker, Range(k, 1), Range(l, 1), Range(m, 1), Range(n, 1), Range(o, 1))).reduce_rank<1>() << '\n';
1450 std::cout << testvar(i, joker, k, l, m, n, o) << '\n';
1451 std::cout << '\n';
1452 }
1453 }
1454 }
1455 }
1456 }
1457 }
1458 for (Index i=0; i<9; i++) {
1459 for (Index j=0; j<2; j++) {
1460 for (Index l=0; l<3; l++) {
1461 for (Index m=0; m<5; m++) {
1462 for (Index n=0; n<7; n++) {
1463 for (Index o=0; o<11; o++) {
1464 std::cout << Tensor7(testvar(Range(i, 1), Range(j, 1), joker, Range(l, 1), Range(m, 1), Range(n, 1), Range(o, 1))).reduce_rank<2>() << '\n';
1465 std::cout << testvar(i, j, joker, l, m, n, o) << '\n';
1466 std::cout << '\n';
1467 }
1468 }
1469 }
1470 }
1471 }
1472 }
1473 for (Index i=0; i<9; i++) {
1474 for (Index j=0; j<2; j++) {
1475 for (Index k=0; k<4; k++) {
1476 for (Index m=0; m<5; m++) {
1477 for (Index n=0; n<7; n++) {
1478 for (Index o=0; o<11; o++) {
1479 std::cout << Tensor7(testvar(Range(i, 1), Range(j, 1), Range(k, 1), joker, Range(m, 1), Range(n, 1), Range(o, 1))).reduce_rank<3>() << '\n';
1480 std::cout << testvar(i, j, k, joker, m, n, o) << '\n';
1481 std::cout << '\n';
1482 }
1483 }
1484 }
1485 }
1486 }
1487 }
1488 for (Index i=0; i<9; i++) {
1489 for (Index j=0; j<2; j++) {
1490 for (Index k=0; k<4; k++) {
1491 for (Index l=0; l<3; l++) {
1492 for (Index n=0; n<7; n++) {
1493 for (Index o=0; o<11; o++) {
1494 std::cout << Tensor7(testvar(Range(i, 1), Range(j, 1), Range(k, 1), Range(l, 1), joker, Range(n, 1), Range(o, 1))).reduce_rank<4>() << '\n';
1495 std::cout << testvar(i, j, k, l, joker, n, o) << '\n';
1496 std::cout << '\n';
1497 }
1498 }
1499 }
1500 }
1501 }
1502 }
1503 for (Index i=0; i<9; i++) {
1504 for (Index j=0; j<2; j++) {
1505 for (Index k=0; k<4; k++) {
1506 for (Index l=0; l<3; l++) {
1507 for (Index m=0; m<5; m++) {
1508 for (Index o=0; o<11; o++) {
1509 std::cout << Tensor7(testvar(Range(i, 1), Range(j, 1), Range(k, 1), Range(l, 1), Range(m, 1), joker, Range(o, 1))).reduce_rank<5>() << '\n';
1510 std::cout << testvar(i, j, k, l, m, joker, o) << '\n';
1511 std::cout << '\n';
1512 }
1513 }
1514 }
1515 }
1516 }
1517 }
1518 for (Index i=0; i<9; i++) {
1519 for (Index j=0; j<2; j++) {
1520 for (Index k=0; k<4; k++) {
1521 for (Index l=0; l<3; l++) {
1522 for (Index m=0; m<5; m++) {
1523 for (Index n=0; n<7; n++) {
1524 std::cout << Tensor7(testvar(Range(i, 1), Range(j, 1), Range(k, 1), Range(l, 1), Range(m, 1), Range(n, 1), joker)).reduce_rank<6>() << '\n';
1525 std::cout << testvar(i, j, k, l, m, n, joker) << '\n';
1526 std::cout << '\n';
1527 }
1528 }
1529 }
1530 }
1531 }
1532 }
1533
1534 // Matrix test
1535 for (Index k=0; k<4; k++) {
1536 for (Index l=0; l<3; l++) {
1537 for (Index m=0; m<5; m++) {
1538 for (Index n=0; n<7; n++) {
1539 for (Index o=0; o<11; o++) {
1540 std::cout << Tensor7(testvar(joker, joker, Range(k, 1), Range(l, 1), Range(m, 1), Range(n, 1), Range(o, 1))).reduce_rank<0, 1>() << '\n';
1541 std::cout << testvar(joker, joker, k, l, m, n, o) << '\n';
1542 std::cout << '\n';
1543 }
1544 }
1545 }
1546 }
1547 }
1548 for (Index j=0; j<2; j++) {
1549 for (Index l=0; l<3; l++) {
1550 for (Index m=0; m<5; m++) {
1551 for (Index n=0; n<7; n++) {
1552 for (Index o=0; o<11; o++) {
1553 std::cout << Tensor7(testvar(joker, Range(j, 1), joker, Range(l, 1), Range(m, 1), Range(n, 1), Range(o, 1))).reduce_rank<0, 2>() << '\n';
1554 std::cout << testvar(joker, j, joker, l, m, n, o) << '\n';
1555 std::cout << '\n';
1556 }
1557 }
1558 }
1559 }
1560 }
1561 for (Index j=0; j<2; j++) {
1562 for (Index k=0; k<4; k++) {
1563 for (Index m=0; m<5; m++) {
1564 for (Index n=0; n<7; n++) {
1565 for (Index o=0; o<11; o++) {
1566 std::cout << Tensor7(testvar(joker, Range(j, 1), Range(k, 1), joker, Range(m, 1), Range(n, 1), Range(o, 1))).reduce_rank<0, 3>() << '\n';
1567 std::cout << testvar(joker, j, k, joker, m, n, o) << '\n';
1568 std::cout << '\n';
1569 }
1570 }
1571 }
1572 }
1573 }
1574 for (Index j=0; j<2; j++) {
1575 for (Index k=0; k<4; k++) {
1576 for (Index l=0; l<3; l++) {
1577 for (Index n=0; n<7; n++) {
1578 for (Index o=0; o<11; o++) {
1579 std::cout << Tensor7(testvar(joker, Range(j, 1), Range(k, 1), Range(l, 1), joker, Range(n, 1), Range(o, 1))).reduce_rank<0, 4>() << '\n';
1580 std::cout << testvar(joker, j, k, l, joker, n, o) << '\n';
1581 std::cout << '\n';
1582 }
1583 }
1584 }
1585 }
1586 }
1587 for (Index j=0; j<2; j++) {
1588 for (Index k=0; k<4; k++) {
1589 for (Index l=0; l<3; l++) {
1590 for (Index m=0; m<5; m++) {
1591 for (Index o=0; o<11; o++) {
1592 std::cout << Tensor7(testvar(joker, Range(j, 1), Range(k, 1), Range(l, 1), Range(m, 1), joker, Range(o, 1))).reduce_rank<0, 5>() << '\n';
1593 std::cout << testvar(joker, j, k, l, m, joker, o) << '\n';
1594 std::cout << '\n';
1595 }
1596 }
1597 }
1598 }
1599 }
1600 for (Index j=0; j<2; j++) {
1601 for (Index k=0; k<4; k++) {
1602 for (Index l=0; l<3; l++) {
1603 for (Index m=0; m<5; m++) {
1604 for (Index n=0; n<7; n++) {
1605 std::cout << Tensor7(testvar(joker, Range(j, 1), Range(k, 1), Range(l, 1), Range(m, 1), Range(n, 1), joker)).reduce_rank<0, 6>() << '\n';
1606 std::cout << testvar(joker, j, k, l, m, n, joker) << '\n';
1607 std::cout << '\n';
1608 }
1609 }
1610 }
1611 }
1612 }
1613 for (Index i=0; i<9; i++) {
1614 for (Index l=0; l<3; l++) {
1615 for (Index m=0; m<5; m++) {
1616 for (Index n=0; n<7; n++) {
1617 for (Index o=0; o<11; o++) {
1618 std::cout << Tensor7(testvar(Range(i, 1), joker, joker, Range(l, 1), Range(m, 1), Range(n, 1), Range(o, 1))).reduce_rank<1, 2>() << '\n';
1619 std::cout << testvar(i, joker, joker, l, m, n, o) << '\n';
1620 std::cout << '\n';
1621 }
1622 }
1623 }
1624 }
1625 }
1626 for (Index i=0; i<9; i++) {
1627 for (Index k=0; k<4; k++) {
1628 for (Index m=0; m<5; m++) {
1629 for (Index n=0; n<7; n++) {
1630 for (Index o=0; o<11; o++) {
1631 std::cout << Tensor7(testvar(Range(i, 1), joker, Range(k, 1), joker, Range(m, 1), Range(n, 1), Range(o, 1))).reduce_rank<1, 3>() << '\n';
1632 std::cout << testvar(i, joker, k, joker, m, n, o) << '\n';
1633 std::cout << '\n';
1634 }
1635 }
1636 }
1637 }
1638 }
1639 for (Index i=0; i<9; i++) {
1640 for (Index k=0; k<4; k++) {
1641 for (Index l=0; l<3; l++) {
1642 for (Index n=0; n<7; n++) {
1643 for (Index o=0; o<11; o++) {
1644 std::cout << Tensor7(testvar(Range(i, 1), joker, Range(k, 1), Range(l, 1), joker, Range(n, 1), Range(o, 1))).reduce_rank<1, 4>() << '\n';
1645 std::cout << testvar(i, joker, k, l, joker, n, o) << '\n';
1646 std::cout << '\n';
1647 }
1648 }
1649 }
1650 }
1651 }
1652 for (Index i=0; i<9; i++) {
1653 for (Index k=0; k<4; k++) {
1654 for (Index l=0; l<3; l++) {
1655 for (Index m=0; m<5; m++) {
1656 for (Index o=0; o<11; o++) {
1657 std::cout << Tensor7(testvar(Range(i, 1), joker, Range(k, 1), Range(l, 1), Range(m, 1), joker, Range(o, 1))).reduce_rank<1, 5>() << '\n';
1658 std::cout << testvar(i, joker, k, l, m, joker, o) << '\n';
1659 std::cout << '\n';
1660 }
1661 }
1662 }
1663 }
1664 }
1665 for (Index i=0; i<9; i++) {
1666 for (Index k=0; k<4; k++) {
1667 for (Index l=0; l<3; l++) {
1668 for (Index m=0; m<5; m++) {
1669 for (Index n=0; n<7; n++) {
1670 std::cout << Tensor7(testvar(Range(i, 1), joker, Range(k, 1), Range(l, 1), Range(m, 1), Range(n, 1), joker)).reduce_rank<1, 6>() << '\n';
1671 std::cout << testvar(i, joker, k, l, m, n, joker) << '\n';
1672 std::cout << '\n';
1673 }
1674 }
1675 }
1676 }
1677 }
1678 for (Index i=0; i<9; i++) {
1679 for (Index j=0; j<2; j++) {
1680 for (Index m=0; m<5; m++) {
1681 for (Index n=0; n<7; n++) {
1682 for (Index o=0; o<11; o++) {
1683 std::cout << Tensor7(testvar(Range(i, 1), Range(j, 1), joker, joker, Range(m, 1), Range(n, 1), Range(o, 1))).reduce_rank<2, 3>() << '\n';
1684 std::cout << testvar(i, j, joker, joker, m, n, o) << '\n';
1685 std::cout << '\n';
1686 }
1687 }
1688 }
1689 }
1690 }
1691 for (Index i=0; i<9; i++) {
1692 for (Index j=0; j<2; j++) {
1693 for (Index l=0; l<3; l++) {
1694 for (Index n=0; n<7; n++) {
1695 for (Index o=0; o<11; o++) {
1696 std::cout << Tensor7(testvar(Range(i, 1), Range(j, 1), joker, Range(l, 1), joker, Range(n, 1), Range(o, 1))).reduce_rank<2, 4>() << '\n';
1697 std::cout << testvar(i, j, joker, l, joker, n, o) << '\n';
1698 std::cout << '\n';
1699 }
1700 }
1701 }
1702 }
1703 }
1704 for (Index i=0; i<9; i++) {
1705 for (Index j=0; j<2; j++) {
1706 for (Index l=0; l<3; l++) {
1707 for (Index m=0; m<5; m++) {
1708 for (Index o=0; o<11; o++) {
1709 std::cout << Tensor7(testvar(Range(i, 1), Range(j, 1), joker, Range(l, 1), Range(m, 1), joker, Range(o, 1))).reduce_rank<2, 5>() << '\n';
1710 std::cout << testvar(i, j, joker, l, m, joker, o) << '\n';
1711 std::cout << '\n';
1712 }
1713 }
1714 }
1715 }
1716 }
1717 for (Index i=0; i<9; i++) {
1718 for (Index j=0; j<2; j++) {
1719 for (Index l=0; l<3; l++) {
1720 for (Index m=0; m<5; m++) {
1721 for (Index n=0; n<7; n++) {
1722 std::cout << Tensor7(testvar(Range(i, 1), Range(j, 1), joker, Range(l, 1), Range(m, 1), Range(n, 1), joker)).reduce_rank<2, 6>() << '\n';
1723 std::cout << testvar(i, j, joker, l, m, n, joker) << '\n';
1724 std::cout << '\n';
1725 }
1726 }
1727 }
1728 }
1729 }
1730 for (Index i=0; i<9; i++) {
1731 for (Index j=0; j<2; j++) {
1732 for (Index k=0; k<4; k++) {
1733 for (Index n=0; n<7; n++) {
1734 for (Index o=0; o<11; o++) {
1735 std::cout << Tensor7(testvar(Range(i, 1), Range(j, 1), Range(k, 1), joker, joker, Range(n, 1), Range(o, 1))).reduce_rank<3, 4>() << '\n';
1736 std::cout << testvar(i, j, k, joker, joker, n, o) << '\n';
1737 std::cout << '\n';
1738 }
1739 }
1740 }
1741 }
1742 }
1743 for (Index i=0; i<9; i++) {
1744 for (Index j=0; j<2; j++) {
1745 for (Index k=0; k<4; k++) {
1746 for (Index m=0; m<5; m++) {
1747 for (Index o=0; o<11; o++) {
1748 std::cout << Tensor7(testvar(Range(i, 1), Range(j, 1), Range(k, 1), joker, Range(m, 1), joker, Range(o, 1))).reduce_rank<3, 5>() << '\n';
1749 std::cout << testvar(i, j, k, joker, m, joker, o) << '\n';
1750 std::cout << '\n';
1751 }
1752 }
1753 }
1754 }
1755 }
1756 for (Index i=0; i<9; i++) {
1757 for (Index j=0; j<2; j++) {
1758 for (Index k=0; k<4; k++) {
1759 for (Index m=0; m<5; m++) {
1760 for (Index n=0; n<7; n++) {
1761 std::cout << Tensor7(testvar(Range(i, 1), Range(j, 1), Range(k, 1), joker, Range(m, 1), Range(n, 1), joker)).reduce_rank<3, 6>() << '\n';
1762 std::cout << testvar(i, j, k, joker, m, n, joker) << '\n';
1763 std::cout << '\n';
1764 }
1765 }
1766 }
1767 }
1768 }
1769 for (Index i=0; i<9; i++) {
1770 for (Index j=0; j<2; j++) {
1771 for (Index k=0; k<4; k++) {
1772 for (Index l=0; l<3; l++) {
1773 for (Index o=0; o<11; o++) {
1774 std::cout << Tensor7(testvar(Range(i, 1), Range(j, 1), Range(k, 1), Range(l, 1), joker, joker, Range(o, 1))).reduce_rank<4, 5>() << '\n';
1775 std::cout << testvar(i, j, k, l, joker, joker, o) << '\n';
1776 std::cout << '\n';
1777 }
1778 }
1779 }
1780 }
1781 }
1782 for (Index i=0; i<9; i++) {
1783 for (Index j=0; j<2; j++) {
1784 for (Index k=0; k<4; k++) {
1785 for (Index l=0; l<3; l++) {
1786 for (Index n=0; n<7; n++) {
1787 std::cout << Tensor7(testvar(Range(i, 1), Range(j, 1), Range(k, 1), Range(l, 1), joker, Range(n, 1), joker)).reduce_rank<4, 6>() << '\n';
1788 std::cout << testvar(i, j, k, l, joker, n, joker) << '\n';
1789 std::cout << '\n';
1790 }
1791 }
1792 }
1793 }
1794 }
1795 for (Index i=0; i<9; i++) {
1796 for (Index j=0; j<2; j++) {
1797 for (Index k=0; k<4; k++) {
1798 for (Index l=0; l<3; l++) {
1799 for (Index m=0; m<5; m++) {
1800 std::cout << Tensor7(testvar(Range(i, 1), Range(j, 1), Range(k, 1), Range(l, 1), Range(m, 1), joker, joker)).reduce_rank<5, 6>() << '\n';
1801 std::cout << testvar(i, j, k, l, m, joker, joker) << '\n';
1802 std::cout << '\n';
1803 }
1804 }
1805 }
1806 }
1807 }
1808 }
1809}
1810
1812 try {
1813 wigner3j(1, 0, 1, 0, 0, 0);
1814 } catch(std::exception& e) {
1815 std::cerr << e.what() << '\n';
1816 }
1817
1818 make_wigner_ready(250, 20000000, 3);
1819 wigner3j(1, 0, 1, 0, 0, 0);
1820}
1821
1822int main() {
1823 // test1();
1824 // test2();
1825 // test3();
1826 // test4();
1827 // test5();
1828 // test6();
1829 // test7();
1830 // test8();
1831 // test9();
1832 // test10();
1833 // test11();
1834 // test12();
1835 // test13();
1836 // test14();
1837 // test15();
1838 // test16();
1839 // test17();
1840 // test18();
1841 // test19();
1842 // test20();
1843 // test21();
1844 // test22();
1845 // test23();
1846 // test24();
1847 // test25();
1848 // test26();
1849 // test27();
1850 // test28();
1851 // test29();
1852 // test30();
1853 // test31();
1854 // test32();
1855 // test33();
1856 // test34();
1857 // test35();
1858 // test36();
1859 // Index i = 10000000;
1860 // test37(i);
1861 // test38();
1862 // test39();
1863 // test40();
1864 // test41();
1865 // test42();
1866 // test43();
1867 // test44();
1868 // test45();
1869 // test46();
1870 // test47();
1871 //test48();
1873
1874 // const double tolerance = 1e-9;
1875 // double error;
1876 //
1877 // // Matrix Vector Multiplication.
1878 // error = test_matrix_vector_multiplication(false);
1879 // cout << "Matrix Vector Multiplication: ";
1880 // if (error > tolerance)
1881 // cout << "FAILED, maximum error: " << error << endl;
1882 // else
1883 // cout << "PASSED." << endl;
1884 //
1885 // // Matrix Matrix Multiplication.
1886 // error = test_matrix_multiplication(false);
1887 //
1888 // cout << "Matrix Matrix Multiplication: ";
1889 // if (error > tolerance)
1890 // cout << "FAILED, maximum error: " << error << endl;
1891 // else
1892 // cout << "PASSED." << endl;
1893
1894 //test_diagonal( 100 );
1895 //test_empty();
1896
1897 return 1;
1898}
This file contains the definition of Array.
Index nelem() const ARTS_NOEXCEPT
Number of elements.
Definition: array.h:197
A constant view of a Matrix.
Definition: matpackI.h:1050
Index nrows() const noexcept
Definition: matpackI.h:1061
Index ncols() const noexcept
Definition: matpackI.h:1062
ConstVectorView diagonal() const ARTS_NOEXCEPT
Matrix diagonal as vector.
Definition: matpackI.cc:489
A constant view of a Vector.
Definition: matpackI.h:517
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:541
The MatrixView class.
Definition: matpackI.h:1169
const Numeric * get_c_array() const ARTS_NOEXCEPT
Conversion to plain C-array.
Definition: matpackI.cc:738
The Matrix class.
Definition: matpackI.h:1270
The range class.
Definition: matpackI.h:166
Implements rational numbers to work with other ARTS types.
Definition: rational.h:52
constexpr bool isUndefined() const noexcept
Is the object not defined.
Definition: rational.h:108
constexpr Index Denom() const noexcept
Denominator.
Definition: rational.h:86
constexpr Numeric toNumeric() const noexcept
Converts this to a Numeric.
Definition: rational.h:142
constexpr Index toIndex(int n=1) const noexcept
Converts the value to index by n-scaled division.
Definition: rational.h:134
constexpr Index Nom() const noexcept
Nominator.
Definition: rational.h:83
constexpr int toInt(int n=1) const noexcept
Converts the value to int by n-scaled division in Index form.
Definition: rational.h:153
The Tensor3View class.
Definition: matpackIII.h:244
const Numeric * get_c_array() const ARTS_NOEXCEPT
Conversion to plain C-array.
Definition: matpackIII.cc:315
The Tensor3 class.
Definition: matpackIII.h:344
Vector reduce_rank() &&ARTS_NOEXCEPT
Definition: matpackIII.h:389
The Tensor4View class.
Definition: matpackIV.h:290
The Tensor4 class.
Definition: matpackIV.h:427
Vector reduce_rank() &&ARTS_NOEXCEPT
Definition: matpackIV.h:477
The Tensor5View class.
Definition: matpackV.h:341
The Tensor5 class.
Definition: matpackV.h:514
Vector reduce_rank() &&ARTS_NOEXCEPT
Definition: matpackV.h:567
The Tensor6View class.
Definition: matpackVI.h:630
The Tensor6 class.
Definition: matpackVI.h:1097
Vector reduce_rank() &&ARTS_NOEXCEPT
Definition: matpackVI.h:1154
The Tensor7View class.
Definition: matpackVII.h:1301
The Tensor7 class.
Definition: matpackVII.h:2397
Vector reduce_rank() &&ARTS_NOEXCEPT
Definition: matpackVII.h:2432
The VectorView class.
Definition: matpackI.h:658
Iterator1D begin() ARTS_NOEXCEPT
Return iterator to first element.
Definition: matpackI.cc:144
Iterator1D end() ARTS_NOEXCEPT
Return iterator behind last element.
Definition: matpackI.cc:148
The Vector class.
Definition: matpackI.h:908
void resize(Index n)
Resize function.
Definition: matpackI.cc:411
void mult(MatrixView C, ConstMatrixView A, const Block &B)
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
String describe(ConstTensor7View x)
Describe Tensor7.
Definition: describe.cc:38
Header file for describe.cc.
The declarations of all the exception classes.
#define min(a, b)
#define max(a, b)
MatrixViewMap MapToEigen(Tensor3View &A, const Index &i)
Definition: lin_alg.cc:745
bool is_increasing(ConstVectorView x)
Checks if a vector is sorted and strictly increasing.
Definition: logic.cc:215
bool is_sorted(ConstVectorView x)
Checks if a vector is sorted in ascending order.
Definition: logic.cc:199
bool is_decreasing(ConstVectorView x)
Checks if a vector is sorted in reversed order and is strictly decreasing.
Definition: logic.cc:252
Header file for logic.cc.
Header file for sparse matrices.
void mult_general(VectorView y, const ConstMatrixView &M, const ConstVectorView &x) ARTS_NOEXCEPT
Matrix Vector multiplication.
Definition: matpackI.cc:1189
void transform(VectorView y, double(&my_func)(double), ConstVectorView x)
A generic transform function for vectors, which can be used to implement mathematical functions opera...
Definition: matpackI.cc:1484
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
This file contains the definition of String, the ARTS string class.
constexpr Numeric l(const Index p0, const Index n, const Numeric x, const SortedVectorType &xi, const Index j, const std::pair< Numeric, Numeric > cycle={ -180, 180}) noexcept
constexpr std::string_view Joker
Definition: isotopologues.h:13
constexpr Rational start(Rational Ju, Rational Jl, Polarization type) noexcept
Gives the lowest M for a polarization type of this transition.
Definition: zeemandata.h:80
invlib::Matrix< ArtsMatrix > Matrix
invlib wrapper type for ARTS matrices.
Definition: oem.h:33
invlib::Vector< ArtsVector > Vector
invlib wrapper type for ARTS vectors.
Definition: oem.h:31
#define d
#define v
#define a
#define c
#define b
Contains the rational class definition.
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:739
constexpr Rational numeric2rational(Numeric x, size_t maxdec=4) noexcept
Rational from Numeric.
Definition: rational.h:361
#define M
Definition: rng.cc:165
void test26()
void test2()
void test29()
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
void test30()
void test46()
void junk2(ConstVectorView a)
void test18()
void test_wigner_error()
#define docheck(fn, val, expect)
void test8()
void test4()
Numeric matrix_vector_mult(Index m, Index n, Index ntests, bool verbose)
void junk4(Tensor4View a)
void test9()
void test42()
void test20()
void test27()
Numeric test_matrix_vector_multiplication(bool verbose)
void test15()
void test48()
int test1()
Definition: test_matpack.cc:48
Numeric matrix_mult(Index k, Index m, Index n, Index ntests, Index nsubtests, bool verbose)
Perform matrix multiplication test.
void test25()
Numeric test_matrix_multiplication(bool verbose)
Perform matrix multiplication tests.
void test31()
Numeric by_value(Numeric x)
Definition: test_matpack.cc:42
void test44()
void test39()
void test23()
void fill_with_junk(VectorView x)
Definition: test_matpack.cc:44
void test5()
void test7()
void test6()
Numeric by_reference(const Numeric &x)
Definition: test_matpack.cc:40
void test33()
void test10()
void test24()
void test47()
void test12()
void test19()
constexpr Rational test_numeric2rational(const Index i, const Index maxi, const Rational r=0, const Rational n2r=0)
void test17()
void test13()
void test38()
void test35()
void test40()
void test14()
bool test_diagonal(Index ntests)
Test diagonal vector.
void test37(const Index &i)
void test11()
int main()
void test28()
void test43()
void test41()
void test22()
void test34()
void test21()
void test_empty()
Check if the empty function is working correctly.
void test32()
void test36()
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
Range random_range(Index n)
Generate random sub-range of the range [0, n-1].
Definition: test_utils.cc:273
Utility functions for testing.
Index make_wigner_ready(int largest, int fastest, int size)
Ready Wigner.
Numeric wigner3j(const Rational j1, const Rational j2, const Rational j3, const Rational m1, const Rational m2, const Rational m3)
Wigner 3J symbol.
Wigner symbol interactions.