ARTS  2.4.0(git:4fb77825)
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 <iostream>
27 #include "complex.h"
28 
29 using std::cout;
30 
31 void 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 
112 void 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 
165 void 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 
335 int main() {
336 // test01();
337  test03();
338  return 0;
339 }
test02
void test02()
Definition: test_complex.cc:112
joker
const Joker joker
ComplexVector
The ComplexVector class.
Definition: complex.h:573
abs
#define abs(x)
Definition: legacy_continua.cc:20626
mult
void mult(ComplexVectorView y, const ConstComplexMatrixView &M, const ConstComplexVectorView &x)
Matrix-Vector Multiplication.
Definition: complex.cc:1579
Complex
std::complex< Numeric > Complex
Definition: complex.h:33
complex.h
A class implementing complex numbers for ARTS.
pow
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:628
MapToEigen
ComplexConstMatrixViewMap MapToEigen(const ConstComplexMatrixView &A)
Definition: complex.cc:1655
test03
void test03()
Definition: test_complex.cc:165
abs2
constexpr Numeric abs2(Complex c)
squared magnitude of c
Definition: complex.h:62
main
int main()
Definition: test_complex.cc:335
test01
void test01()
Definition: test_complex.cc:31
ComplexMatrix
The ComplexMatrix class.
Definition: complex.h:858
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
Range
The range class.
Definition: matpackI.h:160
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
transpose
ConstComplexMatrixView transpose(ConstComplexMatrixView m)
Const version of transpose.
Definition: complex.cc:1509
M
#define M
Definition: rng.cc:165