Skip to content
GitLab
Projects
Groups
Snippets
Help
Loading...
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Sign in
Toggle navigation
Open sidebar
kraus
src
Commits
a49eca79
Commit
a49eca79
authored
Mar 28, 2019
by
snuverink_j
Browse files
remove loading std functions into global namespace in headers
parent
c46e953c
Changes
4
Hide whitespace changes
Inline
Side-by-side
Showing
4 changed files
with
112 additions
and
135 deletions
+112
-135
src/Classic/FixedAlgebra/DragtFinnMap.h
src/Classic/FixedAlgebra/DragtFinnMap.h
+20
-19
src/Classic/FixedAlgebra/DragtFinnNormalForm.h
src/Classic/FixedAlgebra/DragtFinnNormalForm.h
+26
-34
src/Classic/FixedAlgebra/FDoubleEigen.h
src/Classic/FixedAlgebra/FDoubleEigen.h
+38
-46
src/Classic/FixedAlgebra/FNormalForm.h
src/Classic/FixedAlgebra/FNormalForm.h
+28
-36
No files found.
src/Classic/FixedAlgebra/DragtFinnMap.h
View file @
a49eca79
...
...
@@ -23,6 +23,7 @@
#include "Algebra/Array1D.h"
#include "FixedAlgebra/FMatrix.h"
#include "FixedAlgebra/FTps.h"
#include <complex>
template
<
class
T
,
int
N
>
class
FLieGenerator
;
...
...
@@ -432,7 +433,7 @@ factorBerzForestIrwin(const FTps<double, 2 * N> &HH) {
FMatrix
<
double
,
2
*
N
,
2
*
N
>
M
=
makeMatrix
(
H_2
);
// Eigenvalues and eigenvectors of M.
FVector
<
complex
<
double
>
,
2
*
N
>
mu
;
// eigenvalues
FVector
<
std
::
complex
<
double
>
,
2
*
N
>
mu
;
// eigenvalues
FMatrix
<
double
,
2
*
N
,
2
*
N
>
V
;
// eigenvectors
FMatrix
<
double
,
2
*
N
,
2
*
N
>
V_inv
;
// inverse of V
{
...
...
@@ -473,10 +474,10 @@ factorBerzForestIrwin(const FTps<double, 2 * N> &HH) {
R_inv
(
i
+
1
,
i
)
=
+
1.0
;
R_inv
(
i
+
1
,
i
+
1
)
=
-
1.0
;
if
(
std
::
abs
(
imag
(
mu
[
i
]))
>
tol
)
{
if
(
std
::
abs
(
std
::
imag
(
mu
[
i
]))
>
tol
)
{
// "stable" eigenvalue pair.
double
c
=
cos
(
imag
(
mu
[
i
]));
double
s
=
s
in
(
imag
(
mu
[
i
]));
double
c
=
std
::
cos
(
std
::
imag
(
mu
[
i
]));
double
s
=
s
td
::
sin
(
std
::
imag
(
mu
[
i
]));
Rot
(
i
,
i
)
=
+
c
;
Rot
(
i
,
i
+
1
)
=
+
s
;
Rot
(
i
+
1
,
i
)
=
-
s
;
...
...
@@ -485,8 +486,8 @@ factorBerzForestIrwin(const FTps<double, 2 * N> &HH) {
I_dir
(
i
+
1
,
i
)
=
1.0
;
}
else
{
// "unstable" eigenvalue pair.
double
ch
=
cosh
(
real
(
mu
[
i
]));
double
sh
=
s
inh
(
real
(
mu
[
i
]));
double
ch
=
std
::
cosh
(
std
::
real
(
mu
[
i
]));
double
sh
=
s
td
::
sinh
(
std
::
real
(
mu
[
i
]));
Rot
(
i
,
i
)
=
+
ch
;
Rot
(
i
,
i
+
1
)
=
-
sh
;
Rot
(
i
+
1
,
i
)
=
-
sh
;
...
...
@@ -515,22 +516,22 @@ factorBerzForestIrwin(const FTps<double, 2 * N> &HH) {
for
(
int
m
=
f
.
getBottomIndex
();
m
<
f
.
getTopIndex
();
m
++
)
{
const
FMonomial
<
2
*
N
>
&
index
=
FTpsData
<
2
*
N
>::
getExponents
(
m
);
complex
<
double
>
factor
=
0.0
;
std
::
complex
<
double
>
factor
=
0.0
;
int
count
=
0
;
for
(
int
j
=
0
;
j
<
2
*
modes
;
j
+=
2
)
{
if
(
std
::
abs
(
imag
(
mu
[
j
]))
>
tol
)
count
+=
index
[
j
+
1
];
if
(
std
::
abs
(
std
::
imag
(
mu
[
j
]))
>
tol
)
count
+=
index
[
j
+
1
];
factor
+=
double
(
index
[
j
])
*
mu
[
j
]
+
double
(
index
[
j
+
1
])
*
mu
[
j
+
1
];
}
if
(
std
::
abs
(
factor
)
>
tol
)
{
// Term can be removed.
factor
=
1.0
/
factor
;
a
[
m
]
=
real
(
factor
);
b
[
m
]
=
imag
(
factor
);
a
[
m
]
=
std
::
real
(
factor
);
b
[
m
]
=
std
::
imag
(
factor
);
}
pi
[
m
]
=
pow
(
-
1.0
,
(
count
+
1
)
/
2
);
pi
[
m
]
=
std
::
pow
(
-
1.0
,
(
count
+
1
)
/
2
);
}
// Compute cal_T^(-1) * f.
...
...
@@ -1084,7 +1085,7 @@ catenateZero(const DragtFinnMap<N> &g) const {
h
.
assign
(
g_mat
*
f_mat
);
// Extract the generators.
int
order
=
max
(
getOrder
(),
g
.
getOrder
());
int
order
=
std
::
max
(
getOrder
(),
g
.
getOrder
());
FMatrix
<
double
,
2
*
N
,
2
*
N
>
g_inv
=
FLUMatrix
<
double
,
2
*
N
>
(
g_mat
).
inverse
();
FTps
<
double
,
2
*
N
>
f_gen
=
itsGenerators
.
filter
(
2
,
itsGenerators
.
getMaxOrder
()).
substitute
(
g_inv
)
...
...
@@ -1355,12 +1356,12 @@ move_g_1(DragtFinnMap<N> &f, DragtFinnMap<N> &g) {
template
<
int
N
>
int
DragtFinnMap
<
N
>::
orderModes
(
FMatrix
<
double
,
2
*
N
,
2
*
N
>
&
V
,
FVector
<
complex
<
double
>
,
2
*
N
>
&
mu
)
{
orderModes
(
FMatrix
<
double
,
2
*
N
,
2
*
N
>
&
V
,
FVector
<
std
::
complex
<
double
>
,
2
*
N
>
&
mu
)
{
// Static constant.
static
const
double
tol
=
1.0e-12
;
FMatrix
<
double
,
2
*
N
,
2
*
N
>
tmat
(
V
);
FVector
<
complex
<
double
>
,
2
*
N
>
tmu
(
mu
);
FVector
<
std
::
complex
<
double
>
,
2
*
N
>
tmu
(
mu
);
int
nDim
=
2
*
N
;
int
n_c
=
0
;
int
n_r
=
0
;
...
...
@@ -1373,7 +1374,7 @@ orderModes(FMatrix<double, 2 * N, 2 * N> &V, FVector<complex<double>, 2 * N> &mu
for
(
int
j
=
0
;
j
<
2
*
N
;
++
j
)
V
(
j
,
nDim
)
=
0.0
;
V
(
nDim
,
nDim
)
=
1.0
;
i
++
;
}
else
if
(
std
::
abs
(
imag
(
tmu
[
i
]))
<
tol
)
{
}
else
if
(
std
::
abs
(
std
::
imag
(
tmu
[
i
]))
<
tol
)
{
// Collect "unstable" modes in lower indices of tmat.
if
(
n_r
!=
i
)
{
tmu
[
n_r
]
=
tmu
[
i
];
...
...
@@ -1417,14 +1418,14 @@ orderModes(FMatrix<double, 2 * N, 2 * N> &V, FVector<complex<double>, 2 * N> &mu
// Swap values to make pair.
if
(
m
!=
i
+
1
)
{
swap
(
tmu
[
m
],
tmu
[
i
+
1
]);
std
::
swap
(
tmu
[
m
],
tmu
[
i
+
1
]);
tmat
.
swapColumns
(
m
,
i
+
1
);
}
// Take positive eigenvalue first.
int
i1
=
i
;
int
i2
=
i
;
if
(
real
(
tmu
[
i
])
>
0.0
)
{
if
(
std
::
real
(
tmu
[
i
])
>
0.0
)
{
++
i2
;
}
else
{
++
i1
;
...
...
@@ -1466,8 +1467,8 @@ orderModes(FMatrix<double, 2 * N, 2 * N> &V, FVector<complex<double>, 2 * N> &mu
if
(
k
!=
i
)
{
// Move eigenvector pair to its place.
swap
(
mu
[
i
],
mu
[
k
]);
swap
(
mu
[
i
+
1
],
mu
[
k
+
1
]);
std
::
swap
(
mu
[
i
],
mu
[
k
]);
std
::
swap
(
mu
[
i
+
1
],
mu
[
k
+
1
]);
V
.
swapColumns
(
i
,
k
);
V
.
swapColumns
(
i
+
1
,
k
+
1
);
}
...
...
src/Classic/FixedAlgebra/DragtFinnNormalForm.h
View file @
a49eca79
...
...
@@ -24,14 +24,6 @@
#include <algorithm>
#include <complex>
using
std
::
complex
;
using
std
::
abs
;
using
std
::
arg
;
using
std
::
imag
;
using
std
::
pow
;
using
std
::
real
;
using
std
::
swap
;
template
<
int
>
class
DragtFinnMap
;
template
<
class
T
,
int
>
class
FLieGenerator
;
template
<
class
T
,
int
,
int
>
class
FMatrix
;
...
...
@@ -77,7 +69,7 @@ public:
const
DragtFinnMap
<
N
>
&
normalisingMap
()
const
;
/// Get eigenvalues of the linear part as a complex vector.
const
FVector
<
complex
<
double
>
,
2
*
N
>
&
eigenValues
()
const
;
const
FVector
<
std
::
complex
<
double
>
,
2
*
N
>
&
eigenValues
()
const
;
/// Get eigenvectors of the linear part in packed form.
const
FMatrix
<
double
,
2
*
N
,
2
*
N
>
&
eigenVectors
()
const
;
...
...
@@ -91,7 +83,7 @@ public:
protected:
// Order the modes of the map and associate them to the planes.
void
orderModes
(
FVector
<
complex
<
double
>
,
2
*
N
>
,
FMatrix
<
double
,
2
*
N
,
2
*
N
>
);
void
orderModes
(
FVector
<
std
::
complex
<
double
>
,
2
*
N
>
,
FMatrix
<
double
,
2
*
N
,
2
*
N
>
);
private:
...
...
@@ -110,7 +102,7 @@ private:
DragtFinnMap
<
N
>
N_scr
;
// The vector of eigenvalues.
FVector
<
complex
<
double
>
,
2
*
N
>
lambda
;
FVector
<
std
::
complex
<
double
>
,
2
*
N
>
lambda
;
// The matrix of eigenvectors.
FMatrix
<
double
,
2
*
N
,
2
*
N
>
V
;
...
...
@@ -175,18 +167,18 @@ DragtFinnNormalForm<N>::DragtFinnNormalForm(const DragtFinnMap<N> &map):
R_inv
(
i
,
i
)
=
R_inv
(
i
,
i
+
1
)
=
R_inv
(
i
+
1
,
i
)
=
1.0
;
R_inv
(
i
+
1
,
i
+
1
)
=
-
1.0
;
if
(
std
::
abs
(
imag
(
lambda
[
i
]))
>
tol
)
{
if
(
std
::
abs
(
std
::
imag
(
lambda
[
i
]))
>
tol
)
{
// Complex eigenvalue pair.
Rot
(
i
,
i
)
=
Rot
(
i
+
1
,
i
+
1
)
=
real
(
lambda
[
i
]);
Rot
(
i
,
i
+
1
)
=
imag
(
lambda
[
i
]);
Rot
(
i
+
1
,
i
)
=
-
imag
(
lambda
[
i
]);
Rot
(
i
,
i
)
=
Rot
(
i
+
1
,
i
+
1
)
=
std
::
real
(
lambda
[
i
]);
Rot
(
i
,
i
+
1
)
=
std
::
imag
(
lambda
[
i
]);
Rot
(
i
+
1
,
i
)
=
-
std
::
imag
(
lambda
[
i
]);
I_dir
(
i
,
i
)
=
I_dir
(
i
+
1
,
i
+
1
)
=
0.0
;
I_dir
(
i
,
i
+
1
)
=
I_dir
(
i
+
1
,
i
)
=
1.0
;
}
else
{
// Real eigenvalue pair.
Rot
(
i
,
i
)
=
Rot
(
i
+
1
,
i
+
1
)
=
real
(
lambda
[
i
]
+
lambda
[
i
+
1
])
*
0.5
;
Rot
(
i
,
i
+
1
)
=
Rot
(
i
+
1
,
i
)
=
real
(
lambda
[
i
]
-
lambda
[
i
+
1
])
*
0.5
;
Rot
(
i
,
i
)
=
Rot
(
i
+
1
,
i
+
1
)
=
std
::
real
(
lambda
[
i
]
+
lambda
[
i
+
1
])
*
0.5
;
Rot
(
i
,
i
+
1
)
=
Rot
(
i
+
1
,
i
)
=
std
::
real
(
lambda
[
i
]
-
lambda
[
i
+
1
])
*
0.5
;
}
}
...
...
@@ -210,12 +202,12 @@ DragtFinnNormalForm<N>::DragtFinnNormalForm(const DragtFinnMap<N> &map):
for
(
int
m
=
f
.
getBottomIndex
();
m
<
f
.
getTopIndex
();
++
m
)
{
const
FMonomial
<
2
*
N
>
&
index
=
FTpsData
<
2
*
N
>::
getExponents
(
m
);
complex
<
double
>
factor
=
1.0
;
std
::
complex
<
double
>
factor
=
1.0
;
int
count
=
0
;
for
(
int
j
=
0
;
j
<
2
*
freedom
;
j
+=
2
)
{
if
(
std
::
abs
(
imag
(
lambda
[
j
]))
>
tol
)
count
+=
index
[
j
+
1
];
factor
*=
(
pow
(
lambda
[
j
],
index
[
j
])
*
pow
(
lambda
[
j
+
1
],
index
[
j
+
1
]));
if
(
std
::
abs
(
std
::
imag
(
lambda
[
j
]))
>
tol
)
count
+=
index
[
j
+
1
];
factor
*=
(
std
::
pow
(
lambda
[
j
],
index
[
j
])
*
std
::
pow
(
lambda
[
j
+
1
],
index
[
j
+
1
]));
}
if
(
std
::
abs
(
1.0
-
factor
)
<
tol
)
{
...
...
@@ -224,11 +216,11 @@ DragtFinnNormalForm<N>::DragtFinnNormalForm(const DragtFinnMap<N> &map):
}
else
{
// Term can be removed.
factor
=
1.0
/
(
1.0
-
factor
);
a
[
m
]
=
real
(
factor
);
b
[
m
]
=
imag
(
factor
);
a
[
m
]
=
std
::
real
(
factor
);
b
[
m
]
=
std
::
imag
(
factor
);
}
pi
[
m
]
=
pow
(
-
1.0
,
int
(
count
+
1
)
/
2
);
pi
[
m
]
=
std
::
pow
(
-
1.0
,
int
(
count
+
1
)
/
2
);
}
// Compute cal_T^(-1) * f and T_omega.
...
...
@@ -299,7 +291,7 @@ FLieGenerator<double, N> DragtFinnNormalForm<N>::invariant(int mode) const {
b
[
j
+
2
]
=
V
(
j
,
2
*
mode
+
1
);
}
if
(
std
::
abs
(
imag
(
lambda
[
2
*
mode
]))
>
tol
)
{
if
(
std
::
abs
(
std
::
imag
(
lambda
[
2
*
mode
]))
>
tol
)
{
return
(
a
*
a
+
b
*
b
);
}
else
{
return
(
a
*
a
-
b
*
b
);
...
...
@@ -320,7 +312,7 @@ const DragtFinnMap<N> &DragtFinnNormalForm<N>::normalisingMap() const {
template
<
int
N
>
const
FVector
<
complex
<
double
>
,
2
*
N
>
&
DragtFinnNormalForm
<
N
>::
eigenValues
()
const
{
const
FVector
<
std
::
complex
<
double
>
,
2
*
N
>
&
DragtFinnNormalForm
<
N
>::
eigenValues
()
const
{
return
lambda
;
}
...
...
@@ -333,7 +325,7 @@ const FMatrix<double, 2 * N, 2 * N> &DragtFinnNormalForm<N>::eigenVectors() cons
template
<
int
N
>
void
DragtFinnNormalForm
<
N
>::
orderModes
(
FVector
<
complex
<
double
>
,
2
*
N
>
tlam
,
FMatrix
<
double
,
2
*
N
,
2
*
N
>
tmat
)
{
(
FVector
<
std
::
complex
<
double
>
,
2
*
N
>
tlam
,
FMatrix
<
double
,
2
*
N
,
2
*
N
>
tmat
)
{
// Static constant.
int
nDim
=
2
*
N
;
int
n_c
=
0
;
...
...
@@ -347,7 +339,7 @@ void DragtFinnNormalForm<N>::orderModes
std
::
fill
(
V
.
col_begin
(
nDim
),
V
.
col_end
(
nDim
),
0.0
);
V
(
nDim
,
nDim
)
=
1.0
;
i
++
;
}
else
if
(
std
::
abs
(
imag
(
tlam
[
i
]))
<
tol
)
{
}
else
if
(
std
::
abs
(
std
::
imag
(
tlam
[
i
]))
<
tol
)
{
// Collect "unstable" modes in lower indices of tmat.
if
(
n_r
!=
i
)
{
tlam
[
n_r
]
=
tlam
[
i
];
...
...
@@ -399,7 +391,7 @@ void DragtFinnNormalForm<N>::orderModes
// Swap values to make pair.
if
(
m
!=
i
+
1
)
{
swap
(
tlam
[
m
],
tlam
[
i
+
1
]);
std
::
swap
(
tlam
[
m
],
tlam
[
i
+
1
]);
tmat
.
swapColumns
(
m
,
i
+
1
);
}
...
...
@@ -453,21 +445,21 @@ void DragtFinnNormalForm<N>::orderModes
if
(
k
!=
i
)
{
// Move eigenvector pair to its place.
swap
(
lambda
[
i
],
lambda
[
k
]);
swap
(
lambda
[
i
+
1
],
lambda
[
k
+
1
]);
std
::
swap
(
lambda
[
i
],
lambda
[
k
]);
std
::
swap
(
lambda
[
i
+
1
],
lambda
[
k
+
1
]);
V
.
swapColumns
(
i
,
k
);
V
.
swapColumns
(
i
+
1
,
k
+
1
);
}
if
(
std
::
abs
(
imag
(
lambda
[
i
]))
>
tol
)
{
if
(
std
::
abs
(
std
::
imag
(
lambda
[
i
]))
>
tol
)
{
// Rotate complex eigenvectors to make their main components real.
double
re
=
V
(
i
,
i
)
/
sqrt
(
V
(
i
,
i
)
*
V
(
i
,
i
)
+
V
(
i
,
i
+
1
)
*
V
(
i
,
i
+
1
));
double
re
=
V
(
i
,
i
)
/
sqrt
(
V
(
i
,
i
)
*
V
(
i
,
i
)
+
V
(
i
,
i
+
1
)
*
V
(
i
,
i
+
1
));
double
im
=
V
(
i
,
i
+
1
)
/
sqrt
(
V
(
i
,
i
)
*
V
(
i
,
i
)
+
V
(
i
,
i
+
1
)
*
V
(
i
,
i
+
1
));
for
(
int
j
=
0
;
j
<
2
*
N
;
j
++
)
{
double
real_part
=
re
*
V
(
j
,
i
)
+
im
*
V
(
j
,
i
+
1
);
double
real_part
=
re
*
V
(
j
,
i
)
+
im
*
V
(
j
,
i
+
1
);
double
imag_part
=
re
*
V
(
j
,
i
+
1
)
-
im
*
V
(
j
,
i
);
V
(
j
,
i
)
=
real_part
;
V
(
j
,
i
)
=
real_part
;
V
(
j
,
i
+
1
)
=
imag_part
;
}
}
...
...
src/Classic/FixedAlgebra/FDoubleEigen.h
View file @
a49eca79
...
...
@@ -24,14 +24,6 @@
#include "FixedAlgebra/FVector.h"
#include <complex>
using
std
::
complex
;
using
std
::
abs
;
using
std
::
imag
;
using
std
::
max
;
using
std
::
real
;
using
std
::
swap
;
// Class FDoubleEigen
// ------------------------------------------------------------------------
/// Eigenvalues and eigenvectors for a real general matrix.
...
...
@@ -55,11 +47,11 @@ public:
/// Get eigenvalues.
// Return eigenvalues as a complex vector.
FVector
<
complex
<
double
>
,
N
>
eigenValues
()
const
;
FVector
<
std
::
complex
<
double
>
,
N
>
eigenValues
()
const
;
/// Get eigenvectors.
// Return eigenvectors as a complex matrix.
FMatrix
<
complex
<
double
>
,
N
,
N
>
eigenVectors
()
const
;
FMatrix
<
std
::
complex
<
double
>
,
N
,
N
>
eigenVectors
()
const
;
/// Get eigenvectors.
// Return eigenvectors packed in a real matrix.
...
...
@@ -95,7 +87,7 @@ private:
void
balbak
(
int
low
,
int
high
,
double
scale
[
N
]);
// Representation of the eigenvalues and eigenvectors.
FVector
<
complex
<
double
>
,
N
>
lambda
;
FVector
<
std
::
complex
<
double
>
,
N
>
lambda
;
FMatrix
<
double
,
N
,
N
>
vectors
;
};
...
...
@@ -188,7 +180,7 @@ FDoubleEigen<N>::~FDoubleEigen()
template
<
int
N
>
FVector
<
complex
<
double
>
,
N
>
FDoubleEigen
<
N
>::
eigenValues
()
const
FVector
<
std
::
complex
<
double
>
,
N
>
FDoubleEigen
<
N
>::
eigenValues
()
const
// Return Eigenvalues as complex vector.
{
return
lambda
;
...
...
@@ -196,22 +188,22 @@ FVector<complex<double>, N> FDoubleEigen<N>::eigenValues() const
template
<
int
N
>
FMatrix
<
complex
<
double
>
,
N
,
N
>
FDoubleEigen
<
N
>::
eigenVectors
()
const
FMatrix
<
std
::
complex
<
double
>
,
N
,
N
>
FDoubleEigen
<
N
>::
eigenVectors
()
const
// Return eigenvectors as a complex matrix.
{
FMatrix
<
complex
<
double
>
,
N
,
N
>
R
;
FMatrix
<
std
::
complex
<
double
>
,
N
,
N
>
R
;
for
(
int
i
=
0
;
i
<
N
;
i
++
)
{
if
(
imag
(
lambda
[
i
])
==
0.0
)
{
if
(
std
::
imag
(
lambda
[
i
])
==
0.0
)
{
// One real eigenvector.
for
(
int
j
=
0
;
j
<
N
;
j
++
)
{
R
[
j
][
i
]
=
complex
<
double
>
(
vectors
[
j
][
i
]);
R
[
j
][
i
]
=
std
::
complex
<
double
>
(
vectors
[
j
][
i
]);
}
}
else
{
// Two complex eigenvectors.
for
(
int
j
=
0
;
j
<
N
;
j
++
)
{
R
[
j
][
i
]
=
complex
<
double
>
(
vectors
[
j
][
i
],
+
vectors
[
j
][
i
+
1
]);
R
[
j
][
i
+
1
]
=
complex
<
double
>
(
vectors
[
j
][
i
],
-
vectors
[
j
][
i
+
1
]);
R
[
j
][
i
]
=
std
::
complex
<
double
>
(
vectors
[
j
][
i
],
+
vectors
[
j
][
i
+
1
]);
R
[
j
][
i
+
1
]
=
std
::
complex
<
double
>
(
vectors
[
j
][
i
],
-
vectors
[
j
][
i
+
1
]);
}
i
++
;
}
...
...
@@ -439,8 +431,8 @@ void FDoubleEigen<N>::elmhes(FMatrix<double, N, N> ©, int low,
index
[
m
]
=
i
;
if
(
i
!=
m
)
{
for
(
int
j
=
m
-
1
;
j
<
N
;
j
++
)
swap
(
copy
[
i
][
j
],
copy
[
m
][
j
]);
for
(
int
j
=
0
;
j
<=
upp
;
j
++
)
swap
(
copy
[
j
][
i
],
copy
[
j
][
m
]);
for
(
int
j
=
m
-
1
;
j
<
N
;
j
++
)
std
::
swap
(
copy
[
i
][
j
],
copy
[
m
][
j
]);
for
(
int
j
=
0
;
j
<=
upp
;
j
++
)
std
::
swap
(
copy
[
j
][
i
],
copy
[
j
][
m
]);
}
if
(
x
!=
0.0
)
{
...
...
@@ -508,8 +500,8 @@ template <int N>
void
FDoubleEigen
<
N
>::
exchange
(
FMatrix
<
double
,
N
,
N
>
&
copy
,
int
j
,
int
m
,
int
low
,
int
upp
)
{
if
(
j
!=
m
)
{
for
(
int
i
=
0
;
i
<=
upp
;
i
++
)
swap
(
copy
[
i
][
j
],
copy
[
i
][
m
]);
for
(
int
i
=
low
;
i
<
N
;
i
++
)
swap
(
copy
[
j
][
i
],
copy
[
m
][
i
]);
for
(
int
i
=
0
;
i
<=
upp
;
i
++
)
std
::
swap
(
copy
[
i
][
j
],
copy
[
i
][
m
]);
for
(
int
i
=
low
;
i
<
N
;
i
++
)
std
::
swap
(
copy
[
j
][
i
],
copy
[
m
][
i
]);
}
}
...
...
@@ -550,8 +542,8 @@ int FDoubleEigen<N>::hqr(FMatrix<double, N, N> &h, int low, int upp)
// while the j-th eigenvalue is being sought.
{
// Store roots isolated by "balance".
for
(
int
i
=
0
;
i
<
low
;
i
++
)
lambda
[
i
]
=
complex
<
double
>
(
h
[
i
][
i
]);
for
(
int
i
=
upp
+
1
;
i
<
N
;
i
++
)
lambda
[
i
]
=
complex
<
double
>
(
h
[
i
][
i
]);
for
(
int
i
=
0
;
i
<
low
;
i
++
)
lambda
[
i
]
=
std
::
complex
<
double
>
(
h
[
i
][
i
]);
for
(
int
i
=
upp
+
1
;
i
<
N
;
i
++
)
lambda
[
i
]
=
std
::
complex
<
double
>
(
h
[
i
][
i
]);
// Compute matrix norm.
double
norm
=
0.0
;
...
...
@@ -680,7 +672,7 @@ int FDoubleEigen<N>::hqr(FMatrix<double, N, N> &h, int low, int upp)
}
one_root:
lambda
[
en
]
=
complex
<
double
>
(
x
+
t
);
lambda
[
en
]
=
std
::
complex
<
double
>
(
x
+
t
);
continue
;
two_roots:
...
...
@@ -692,12 +684,12 @@ two_roots:
if
(
q
>=
0.0
)
{
// Real pair.
z
=
(
p
>
0.0
)
?
(
p
+
z
)
:
(
p
-
z
);
lambda
[
en
-
1
]
=
complex
<
double
>
(
x
+
z
);
lambda
[
en
]
=
complex
<
double
>
((
z
!=
0.0
)
?
(
x
-
w
/
z
)
:
x
+
z
);
lambda
[
en
-
1
]
=
std
::
complex
<
double
>
(
x
+
z
);
lambda
[
en
]
=
std
::
complex
<
double
>
((
z
!=
0.0
)
?
(
x
-
w
/
z
)
:
x
+
z
);
}
else
{
// Complex pair.
lambda
[
en
-
1
]
=
complex
<
double
>
(
x
+
p
,
+
z
);
lambda
[
en
]
=
complex
<
double
>
(
x
+
p
,
-
z
);
lambda
[
en
-
1
]
=
std
::
complex
<
double
>
(
x
+
p
,
+
z
);
lambda
[
en
]
=
std
::
complex
<
double
>
(
x
+
p
,
-
z
);
}
en
--
;
...
...
@@ -759,8 +751,8 @@ int FDoubleEigen<N>::hqr2(FMatrix<double, N, N> &h, int low, int upp)
double
tst1
,
tst2
;
// Store roots isolated by "balance".
for
(
int
i
=
0
;
i
<
low
;
i
++
)
lambda
[
i
]
=
complex
<
double
>
(
h
[
i
][
i
]);
for
(
int
i
=
upp
+
1
;
i
<
N
;
i
++
)
lambda
[
i
]
=
complex
<
double
>
(
h
[
i
][
i
]);
for
(
int
i
=
0
;
i
<
low
;
i
++
)
lambda
[
i
]
=
std
::
complex
<
double
>
(
h
[
i
][
i
]);
for
(
int
i
=
upp
+
1
;
i
<
N
;
i
++
)
lambda
[
i
]
=
std
::
complex
<
double
>
(
h
[
i
][
i
]);
// Compute matrix norm.
double
norm
=
0.0
;
...
...
@@ -900,7 +892,7 @@ int FDoubleEigen<N>::hqr2(FMatrix<double, N, N> &h, int low, int upp)
}
one_root:
lambda
[
en
]
=
complex
<
double
>
(
h
[
en
][
en
]
=
x
+
t
);
lambda
[
en
]
=
std
::
complex
<
double
>
(
h
[
en
][
en
]
=
x
+
t
);
continue
;
two_roots:
...
...
@@ -945,8 +937,8 @@ two_roots:
}
}
else
{
// Complex pair.
lambda
[
en
-
1
]
=
complex
<
double
>
(
x
+
p
,
+
z
);
lambda
[
en
]
=
complex
<
double
>
(
x
+
p
,
-
z
);
lambda
[
en
-
1
]
=
std
::
complex
<
double
>
(
x
+
p
,
+
z
);
lambda
[
en
]
=
std
::
complex
<
double
>
(
x
+
p
,
-
z
);
}
en
--
;
...
...
@@ -957,8 +949,8 @@ two_roots:
if
(
norm
==
0.0
)
return
0
;
for
(
int
en
=
N
;
en
--
>
0
;)
{
p
=
real
(
lambda
[
en
]);
q
=
imag
(
lambda
[
en
]);
p
=
std
::
real
(
lambda
[
en
]);
q
=
std
::
imag
(
lambda
[
en
]);
if
(
q
<
0.0
)
{
// Complex vector.
...
...
@@ -988,21 +980,21 @@ two_roots:
sa
+=
h
[
i
][
j
]
*
h
[
j
][
en
];
}
if
(
imag
(
lambda
[
i
])
<
0.0
)
{
if
(
std
::
imag
(
lambda
[
i
])
<
0.0
)
{
z
=
w
;
r
=
ra
;
s
=
sa
;
}
else
{
m
=
i
;
if
(
imag
(
lambda
[
i
])
==
0.0
)
{
if
(
std
::
imag
(
lambda
[
i
])
==
0.0
)
{
cdiv
(
-
ra
,
-
sa
,
w
,
q
,
h
[
i
][
en
-
1
],
h
[
i
][
en
]);
}
else
{
// Solve complex equations.
x
=
h
[
i
][
i
+
1
];
y
=
h
[
i
+
1
][
i
];
double
vr
=
(
real
(
lambda
[
i
])
-
p
)
*
(
real
(
lambda
[
i
])
-
p
)
+
imag
(
lambda
[
i
])
*
imag
(
lambda
[
i
])
-
q
*
q
;
double
vi
=
real
(
lambda
[
i
]
-
p
)
*
2.0
*
q
;
double
vr
=
(
std
::
real
(
lambda
[
i
])
-
p
)
*
(
std
::
real
(
lambda
[
i
])
-
p
)
+
std
::
imag
(
lambda
[
i
])
*
std
::
imag
(
lambda
[
i
])
-
q
*
q
;
double
vi
=
std
::
real
(
lambda
[
i
]
-
p
)
*
2.0
*
q
;
if
(
vr
==
0.0
&&
vi
==
0.0
)
{
tst1
=
...
...
@@ -1028,7 +1020,7 @@ two_roots:
}
// Overflow control.
t
=
max
(
std
::
abs
(
h
[
i
][
en
-
1
]),
std
::
abs
(
h
[
i
][
en
]));
t
=
std
::
max
(
std
::
abs
(
h
[
i
][
en
-
1
]),
std
::
abs
(
h
[
i
][
en
]));
if
(
t
!=
0.0
&&
(
t
+
1.0
/
t
)
<=
t
)
{
for
(
int
j
=
i
;
j
<=
en
;
j
++
)
{
...
...
@@ -1051,12 +1043,12 @@ two_roots:
r
=
0.0
;
for
(
int
j
=
m
;
j
<=
en
;
j
++
)
r
+=
h
[
i
][
j
]
*
h
[
j
][
en
];
if
(
imag
(
lambda
[
i
])
<
0.0
)
{
if
(
std
::
imag
(
lambda
[
i
])
<
0.0
)
{
z
=
w
;
s
=
r
;
}
else
{
m
=
i
;
if
(
imag
(
lambda
[
i
])
==
0.0
)
{
if
(
std
::
imag
(
lambda
[
i
])
==
0.0
)
{
t
=
w
;
if
(
t
==
0.0
)
{
t
=
tst1
=
norm
;
...
...
@@ -1072,8 +1064,8 @@ two_roots:
// Solve real equations.
x
=
h
[
i
][
i
+
1
];
y
=
h
[
i
+
1
][
i
];
q
=
(
real
(
lambda
[
i
])
-
p
)
*
(
real
(
lambda
[
i
])
-
p
)
+
imag
(
lambda
[
i
])
*
imag
(
lambda
[
i
]);
q
=
(
std
::
real
(
lambda
[
i
])
-
p
)
*
(
std
::
real
(
lambda
[
i
])
-
p
)
+
std
::
imag
(
lambda
[
i
])
*
std
::
imag
(
lambda
[
i
]);
t
=
(
x
*
s
-
z
*
r
)
/
q
;
h
[
i
][
en
]
=
t
;
h
[
i
+
1
][
en
]
=
(
std
::
abs
(
x
)
>
std
::
abs
(
z
))
?
...
...
src/Classic/FixedAlgebra/FNormalForm.h
View file @
a49eca79
...
...
@@ -25,14 +25,6 @@
#include <algorithm>
#include <complex>
using
std
::
complex
;
using
std
::
abs
;
using
std
::
arg
;
using
std
::
imag
;
using
std
::
pow
;
using
std
::
real
;
using
std
::
swap
;
template
<
class
T
,
int
,
int
>
class
FMatrix
;
template
<
class
T
,
int
>
class
FTps
;
template
<
class
T
,
int
>
class
FVps
;
...
...
@@ -78,7 +70,7 @@ public:
const
FTps
<
double
,
N
>
&
normalisingMap
()
const
;
/// Get eigenvalues of the linear part as a complex vector.
const
FVector
<
complex
<
double
>
,
N
>
&
eigenValues
()
const
;
const
FVector
<
std
::
complex
<
double
>
,
N
>
&
eigenValues
()
const
;
/// Get eigenvectors of the linear part in packed form.
const
FMatrix
<
double
,
N
,
N
>
&
eigenVectors
()
const
;
...
...
@@ -92,7 +84,7 @@ public:
protected:
// Order the modes of the map and associate them to the planes.
void
orderModes
(
FVector
<
complex
<
double
>
,
N
>
,
FMatrix
<
double
,
N
,
N
>
);
void
orderModes
(
FVector
<
std
::
complex
<
double
>
,
N
>
,
FMatrix
<
double
,
N
,
N
>
);
private:
...
...
@@ -111,7 +103,7 @@ private:
FTps
<
double
,
N
>
N_Lie
;
// The vector of eigenvalues.
FVector
<
complex
<
double
>
,
N
>
lambda
;
FVector
<
std
::
complex
<
double
>
,
N
>
lambda
;
// The matrix of eigenvectors.
FMatrix
<
double
,
N
,
N
>
V
;
...
...
@@ -180,8 +172,8 @@ FNormalForm<N>::FNormalForm(const FVps<double, N> &M_scr):
// Store linear part of normal form.
double
c1
,
c2
;
if
(
std
::
abs
(
imag
(
lambda
[
i
]))
>
tol
)
{
c1
=
c2
=
arg
(
lambda
[
i
])
/
2.0
;
if
(
std
::
abs
(
std
::
imag
(
lambda
[
i
]))
>
tol
)
{
c1
=
c2
=
std
::
arg
(
lambda
[
i
])
/
2.0
;
}
else
{
c1
=
log
(
std
::
abs
(
lambda
[
i
]))
/
2.0
;