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