ARTS 2.5.4 (git: bcd8c674)
test_complex.cc
Go to the documentation of this file.
1/* Copyright (C) 2002-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
26#include "matpack_complex.h"
27#include <iostream>
28
29using std::cout;
30
31void test01() {
32 Complex a;
33 Complex b(3., 0.);
34
35 //cout << "a = " << a << "\n";
36 cout << "b = " << b << "\n";
37
38 a = b;
39 cout << "a = " << a << "\n";
40
41 a = Complex(1., 1.);
42 cout << "a = " << a << "\n";
43 cout << "a.abs() = " << abs(a) << "\n";
44 cout << "a.arg()= " << arg(a) * 57.2957914331333 << " degree"
45 << "\n";
46
47 Complex c;
48 c = a + b;
49 cout << "c = " << c << "\n";
50
51 cout << "a += b: " << (a += b) << "\n";
52
53 cout << "c + 3 = " << (c + 3.) << "\n";
54
55 cout << "c + 3i = c + Complex (0,3) = " << (c + Complex(0., 3.)) << "\n";
56
57 Complex d;
58 a = Complex(1., 1.);
59 b = Complex(-1., -1.);
60 d = a * b;
61 cout << "d = " << d << "\n";
62 Complex e;
63 e = a / b;
64 cout << "e = " << e << "\n";
65 Complex f;
66 f = pow(a, 0) + pow(b, 0);
67 cout << "f = " << f << "\n";
68 Complex g;
69 g = pow(a, 1) + pow(b, 1);
70 cout << "g = " << g << "\n";
71 Complex h;
72 h = pow(a, 2) + pow(b, 2) + pow(a, 3) + pow(b, 3);
73 cout << "h = " << h << "\n";
74
75 ComplexMatrix A(4, 4);
76 ComplexMatrix B(4, 4);
77 ComplexMatrix X1(4, 1);
78 ComplexVector X2(4);
79 ComplexMatrix C(4, 4);
80 for (Index i = 0; i < 4; i++) {
81 X1(i, 0) = Complex((Numeric)i * 5.0 + 2.0, (Numeric)i + 1.0);
82 X2[i] = Complex((Numeric)i * 5.0 + 2.0, (Numeric)i + 1.0);
83 for (Index j = 0; j < 4; j++) {
84 A(i, j) = Complex((Numeric)i + (Numeric)j, 0.0) + Complex(1.0, 1.0);
85 B(i, j) =
86 Complex(2.0 * (Numeric)i + 4.0 * (Numeric)j, 0.0) + Complex(3.0, 3.0);
87 }
88 }
89
90 std::cout << "\n";
91 mult(C, A, B);
92 std::cout << MapToEigen(A) << "\nx\n"
93 << MapToEigen(B) << "\n=\n"
94 << MapToEigen(C) << "\n\n";
95 mult(C, C, C);
96 mult(C, C, C);
97 std::cout
98 << "Same matrix can be both input and output, here is the above to the power of 4:\n"
99 << MapToEigen(C) << "\n\n";
100 mult(C(joker, 1), A, X2);
101 std::cout << MapToEigen(A) << "\nx\n"
102 << MapToEigen(X2) << "\n=\n"
103 << MapToEigen(C(joker, 1)) << "\n\n";
104 mult(C(Range(1, 3), Range(0, 3)),
105 A(Range(1, 3), Range(1, 2)),
106 B(Range(1, 2), Range(1, 3)));
107 std::cout << MapToEigen(A(Range(1, 3), Range(1, 2))) << "\nx\n"
108 << MapToEigen(B(Range(1, 2), Range(1, 3))) << "\n=\n"
109 << MapToEigen(C(Range(1, 3), Range(0, 3))) << "\n\n";
110}
111
112void test02() {
113 constexpr Numeric nul = 0.f;
114 constexpr Numeric one = 1.f;
115 constexpr Numeric two = 2.f;
116 constexpr Numeric tre = 3.f;
117 constexpr Numeric fyr = 4.f;
118
119 constexpr Complex r(one, nul);
120 constexpr Complex i(nul, one);
121
122 constexpr Complex a(one, two);
123 constexpr Complex b(tre, two);
124 constexpr Complex c(tre, fyr);
125
126 // Test that every operator is treated as Numeric regardless of basic types
127 static_assert(tre + a == 3 + a, "Bad operator+ int Complex");
128 static_assert(tre - a == 3 - a, "Bad operator- int Complex");
129 static_assert(tre * a == 3 * a, "Bad operator* int Complex");
130 static_assert(tre / a == 3 / a, "Bad operator/ int Complex");
131 static_assert(a + fyr == a + 4, "Bad operator+ Complex int");
132 static_assert(a - fyr == a - 4, "Bad operator- Complex int");
133 static_assert(a * fyr == a * 4, "Bad operator* Complex int");
134 static_assert(a / fyr == a / 4, "Bad operator/ Complex int");
135 static_assert(tre + a == 3.f + a, "Bad operator+ float Complex");
136 static_assert(tre - a == 3.f - a, "Bad operator- float Complex");
137 static_assert(tre * a == 3.f * a, "Bad operator* float Complex");
138 static_assert(tre / a == 3.f / a, "Bad operator/ float Complex");
139 static_assert(a + fyr == a + 4.f, "Bad operator+ Complex float");
140 static_assert(a - fyr == a - 4.f, "Bad operator- Complex float");
141 static_assert(a * fyr == a * 4.f, "Bad operator* Complex float");
142 static_assert(a / fyr == a / 4.f, "Bad operator/ Complex float");
143
144 // Test the most basic test of (1, 0)*z and (0, 1)*z
145 static_assert(r * a == Complex(one * one, one * two),
146 "Bad Complex(1, 0) Complex ");
147 static_assert(i * a == Complex(-one * two, one * one),
148 "Bad Complex(0, 1) Complex ");
149
150 // Test key helper function
151 static_assert(abs2(c) == tre * tre + fyr * fyr, "Bad Numeric abs2(Complex)");
152
153 // Test Complex-Complex operators
154 static_assert(a + b == Complex(one + tre, two + two),
155 "Bad operator+ Complex Complex");
156 static_assert(a - b == Complex(one - tre, two - two),
157 "Bad operator- Complex Complex");
158 static_assert(a * b == Complex(one * tre - two * two, one * two + two * tre),
159 "Bad operator* Complex Complex");
160 static_assert(
161 a / b == Complex(one * tre + two * two, two * tre - one * two) / abs2(b),
162 "Bad operator/ Complex Complex");
163}
164
165void test03()
166{
167 constexpr Index n=4;
168 constexpr Index col=2;
169 constexpr Index row=2;
170 static_assert(n>=4, "Req for size of val to follow");
171 static_assert(col<n, "Req for size of val to follow");
172 static_assert(row<n, "Req for size of val to follow");
173
174 constexpr Range r0 = Range(joker);
175 constexpr Range r1 = Range(1, 2);
176 constexpr Range r2 = Range(1, 2, 2);
177 constexpr Range r3 = Range(n-1, 2, -1);
178 constexpr Range r4 = Range(n-1, 2, -2);
179 constexpr Range r5 = Range(n-1, n, -1);
180 const auto ranges={r0, r1, r2, r3, r4, r5};
181
182 ComplexVector V(n, n);
183 for (Index i=0; i<n; i++)
184 V[i] = Complex(Numeric(i), Numeric(i+2*n));
185 const ComplexVector cV=V;
186
187 ComplexMatrix M(n, n);
188 for (Index i=0; i<n; i++)
189 for (Index j=0; j<n; j++)
190 M(i, j) = Complex(Numeric(i), Numeric(j+2*n));
191 const ComplexMatrix cM=M;
192
193 std::cout<<"Vector: ";
194 std::cout<<MapToEigen(V).transpose()<<'\n';
195 std::cout<<"----------------------------------------------------------------------------------------------------\n";
196
197 std::cout<<"Matrix:\n";
198 std::cout<<MapToEigen(M)<<"\n";
199 std::cout<<"----------------------------------------------------------------------------------------------------\n";
200 std::cout<<"Vector\n";
201 std::cout<<"----------------------------------------------------------------------------------------------------\n";
202 for (auto r: ranges) {
203 std::cout<<"Vector["<<r<<"]: "<<MapToEigen(V[r]).transpose()<<'\n';
204 std::cout<<"Real: "<<MapToEigen(V[r].real()).transpose()<<'\n';
205 std::cout<<"Imag: "<<MapToEigen(V[r].imag()).transpose()<<'\n';
206 std::cout<<"---------------------------------------------------\n";
207 }
208 std::cout<<"----------------------------------------------------------------------------------------------------\n";
209 std::cout<<"const Vector\n";
210 std::cout<<"----------------------------------------------------------------------------------------------------\n";
211 for (auto r: ranges) {
212 std::cout<<"const Vector["<<r<<"]: "<<MapToEigen(cV[r]).transpose()<<'\n';
213 std::cout<<"Real: "<<MapToEigen(V[r].real()).transpose()<<'\n';
214 std::cout<<"Imag: "<<MapToEigen(V[r].imag()).transpose()<<'\n';
215 std::cout<<"---------------------------------------------------\n";
216 }
217 std::cout<<"----------------------------------------------------------------------------------------------------\n";
218 std::cout<<"RowView Matrix\n";
219 std::cout<<"----------------------------------------------------------------------------------------------------\n";
220 for (auto r: ranges) {
221 std::cout<<"Matrix("<<row<<", "<<r<<"): "<<MapToEigen(M(row,r)).transpose()<<'\n';
222 std::cout<<"Real: "<<MapToEigen(M(row,r).real()).transpose()<<'\n';
223 std::cout<<"Imag: "<<MapToEigen(M(row,r).imag()).transpose()<<'\n';
224 std::cout<<"---------------------------------------------------\n";
225 }
226 std::cout<<"----------------------------------------------------------------------------------------------------\n";
227 std::cout<<"ColView Matrix\n";
228 std::cout<<"----------------------------------------------------------------------------------------------------\n";
229 for (auto r: ranges) {
230 std::cout<<"Matrix("<<r<<", "<<col<<"): "<<MapToEigen(M(r, col)).transpose()<<'\n';
231 std::cout<<"Real: "<<MapToEigen(M(r, col).real()).transpose()<<'\n';
232 std::cout<<"Imag: "<<MapToEigen(M(r, col).imag()).transpose()<<'\n';
233 std::cout<<"---------------------------------------------------\n";
234 }
235 std::cout<<"----------------------------------------------------------------------------------------------------\n";
236 std::cout<<"RowView Matrix-transpose\n";
237 std::cout<<"----------------------------------------------------------------------------------------------------\n";
238 for (auto r: ranges) {
239 std::cout<<"transpose(Matrix)("<<row<<", "<<r<<"): "<<MapToEigen(transpose(M)(row,r)).transpose()<<'\n';
240 std::cout<<"Real: "<<MapToEigen(transpose(M)(row,r).real()).transpose()<<'\n';
241 std::cout<<"Imag: "<<MapToEigen(transpose(M)(row,r).imag()).transpose()<<'\n';
242 std::cout<<"---------------------------------------------------\n";
243 }
244 std::cout<<"----------------------------------------------------------------------------------------------------\n";
245 std::cout<<"ColView Matrix transpose\n";
246 std::cout<<"----------------------------------------------------------------------------------------------------\n";
247 for (auto r: ranges) {
248 std::cout<<"transpose(Matrix)("<<r<<", "<<col<<"): "<<MapToEigen(transpose(M)(r, col)).transpose()<<'\n';
249 std::cout<<"Real: "<<MapToEigen(transpose(M)(r, col).real()).transpose()<<'\n';
250 std::cout<<"Imag: "<<MapToEigen(transpose(M)(r, col).imag()).transpose()<<'\n';
251 std::cout<<"---------------------------------------------------\n";
252 }
253 std::cout<<"----------------------------------------------------------------------------------------------------\n";
254 std::cout<<"Sub-Matrix\n";
255 std::cout<<"----------------------------------------------------------------------------------------------------\n";
256 for (auto rr: ranges) {
257 for (auto rc: ranges) {
258 std::cout<<"Matrix("<<rr<<", "<<rc<<"):\n"<<MapToEigen(M(rr,rc)).transpose()<<'\n';
259 std::cout<<"Real:\n"<<MapToEigen(M(rr,rc).real()).transpose()<<'\n';
260 std::cout<<"Imag:\n"<<MapToEigen(M(rr,rc).imag()).transpose()<<'\n';
261 std::cout<<"---------------------------------------------------\n";
262 }
263 }
264 std::cout<<"----------------------------------------------------------------------------------------------------\n";
265 std::cout<<"Sub-Matrix transpose\n";
266 std::cout<<"----------------------------------------------------------------------------------------------------\n";
267 for (auto rr: ranges) {
268 for (auto rc: ranges) {
269 std::cout<<"transpose(Matrix)("<<rr<<", "<<rc<<"):\n"<<MapToEigen(transpose(M)(rr,rc)).transpose()<<'\n';
270 std::cout<<"Real:\n"<<MapToEigen(transpose(M)(rr,rc).real()).transpose()<<'\n';
271 std::cout<<"Imag:\n"<<MapToEigen(transpose(M)(rr,rc).imag()).transpose()<<'\n';
272 std::cout<<"---------------------------------------------------\n";
273 }
274 }
275 std::cout<<"----------------------------------------------------------------------------------------------------\n";
276 std::cout<<"RowView Matrix const\n";
277 std::cout<<"----------------------------------------------------------------------------------------------------\n";
278 for (auto r: ranges) {
279 std::cout<<"const Matrix("<<row<<", "<<r<<"): "<<MapToEigen(cM(row,r)).transpose()<<'\n';
280 std::cout<<"Real: "<<MapToEigen(cM(row,r).real()).transpose()<<'\n';
281 std::cout<<"Imag: "<<MapToEigen(cM(row,r).imag()).transpose()<<'\n';
282 std::cout<<"---------------------------------------------------\n";
283 }
284 std::cout<<"----------------------------------------------------------------------------------------------------\n";
285 std::cout<<"ColView Matrix const\n";
286 std::cout<<"----------------------------------------------------------------------------------------------------\n";
287 for (auto r: ranges) {
288 std::cout<<"const Matrix("<<r<<", "<<col<<"): "<<MapToEigen(cM(r, col)).transpose()<<'\n';
289 std::cout<<"Real: "<<MapToEigen(cM(r, col).real()).transpose()<<'\n';
290 std::cout<<"Imag: "<<MapToEigen(cM(r, col).imag()).transpose()<<'\n';
291 std::cout<<"---------------------------------------------------\n";
292 }
293 std::cout<<"----------------------------------------------------------------------------------------------------\n";
294 std::cout<<"RowView Matrix-transpose const\n";
295 std::cout<<"----------------------------------------------------------------------------------------------------\n";
296 for (auto r: ranges) {
297 std::cout<<"transpose(const Matrix)("<<row<<", "<<r<<"): "<<MapToEigen(transpose(cM)(row,r)).transpose()<<'\n';
298 std::cout<<"Real: "<<MapToEigen(transpose(cM)(row,r).real()).transpose()<<'\n';
299 std::cout<<"Imag: "<<MapToEigen(transpose(cM)(row,r).imag()).transpose()<<'\n';
300 std::cout<<"---------------------------------------------------\n";
301 }
302 std::cout<<"----------------------------------------------------------------------------------------------------\n";
303 std::cout<<"ColView Matrix-transpose const\n";
304 std::cout<<"----------------------------------------------------------------------------------------------------\n";
305 for (auto r: ranges) {
306 std::cout<<"transpose(const Matrix)("<<r<<", "<<col<<"): "<<MapToEigen(transpose(cM)(r, col)).transpose()<<'\n';
307 std::cout<<"Real: "<<MapToEigen(transpose(cM)(r, col).real()).transpose()<<'\n';
308 std::cout<<"Imag: "<<MapToEigen(transpose(cM)(r, col).imag()).transpose()<<'\n';
309 std::cout<<"---------------------------------------------------\n";
310 }
311 std::cout<<"----------------------------------------------------------------------------------------------------\n";
312 std::cout<<"Sub-Matrix const\n";
313 std::cout<<"----------------------------------------------------------------------------------------------------\n";
314 for (auto rr: ranges) {
315 for (auto rc: ranges) {
316 std::cout<<"Matrix("<<rr<<", "<<rc<<"):\n"<<MapToEigen(M(rr,rc)).transpose()<<'\n';
317 std::cout<<"Real:\n"<<MapToEigen(M(rr,rc).real()).transpose()<<'\n';
318 std::cout<<"Imag:\n"<<MapToEigen(M(rr,rc).imag()).transpose()<<'\n';
319 std::cout<<"---------------------------------------------------\n";
320 }
321 }
322 std::cout<<"----------------------------------------------------------------------------------------------------\n";
323 std::cout<<"Sub-Matrix transpose const\n";
324 std::cout<<"----------------------------------------------------------------------------------------------------\n";
325 for (auto rr: ranges) {
326 for (auto rc: ranges) {
327 std::cout<<"transpose(const Matrix)("<<rr<<", "<<rc<<"):\n"<<MapToEigen(transpose(cM)(rr,rc)).transpose()<<'\n';
328 std::cout<<"Real:\n"<<MapToEigen(transpose(cM)(rr,rc).real()).transpose()<<'\n';
329 std::cout<<"Imag:\n"<<MapToEigen(transpose(cM)(rr,rc).imag()).transpose()<<'\n';
330 std::cout<<"---------------------------------------------------\n";
331 }
332 }
333}
334
335int main() {
336// test01();
337 test03();
338
339 Complex a(0, 0);
340 real_val(a) += 1;
341 imag_val(a) -= 2;
342 imag_val(a) *= 2;
343 std::cout << a << '\n';
344 return 0;
345}
The ComplexMatrix class.
The ComplexVector class.
The range class.
Definition: matpackI.h:177
void mult(MatrixView C, ConstMatrixView A, const Block &B)
#define abs(x)
MatrixViewMap MapToEigen(Tensor3View &A, const Index &i)
Definition: lin_alg.cc:684
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.
Numeric & real_val(Complex &c) noexcept
Return a reference to the real value of c.
constexpr Numeric imag(Complex c) noexcept
imag
const Joker joker
Numeric & imag_val(Complex &c) noexcept
Return a reference to the imaginary value of c.
constexpr Numeric abs2(Complex c) noexcept
squared magnitude of c
std::complex< Numeric > Complex
constexpr Numeric real(Complex c) noexcept
real
constexpr Numeric r0
The reference radius in IGRF13.
Definition: igrf13.cc:203
#define d
#define a
#define c
#define b
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:736
#define M
Definition: rng.cc:165
void test03()
void test02()
void test01()
Definition: test_complex.cc:31
int main()