Commit ff347732 authored by snuverink_j's avatar snuverink_j
Browse files

remove more unused classes

parent 7f4c83d4
......@@ -15,7 +15,6 @@ add_opal_sources (${_SRCS})
set (HDRS
Array1D.h
Array2D.h
LieMap.h
LUMatrix.h
Matrix.h
SliceIterator.h
......@@ -28,7 +27,6 @@ set (HDRS
Vector.h
Vps.h
Vps.hpp
VpsInvMap.h
VpsMap.h
)
......
#ifndef CLASSIC_LieMap_HH
#define CLASSIC_LieMap_HH
// ------------------------------------------------------------------------
// $RCSfile: LieMap.h,v $
// ------------------------------------------------------------------------
// $Revision: 1.1.1.1.2.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Template class: LieMap
//
// ------------------------------------------------------------------------
// Class category: Algebra
// ------------------------------------------------------------------------
//
// $Date: 2004/11/18 22:18:05 $
// $Author: jsberg $
//
// ------------------------------------------------------------------------
#include "Algebra/Matrix.h"
#include "Algebra/Vector.h"
#include "Algebra/Tps.h"
#include "Algebra/TpsData.h"
#include "Algebra/VpsInvMap.h"
#include "Utilities/ConvergenceError.h"
#include "Utilities/LogicalError.h"
#include <iosfwd>
// Template class LieMap<T>
// ------------------------------------------------------------------------
/// Lie algebraic map.
// LieMap<T> is a truncated Taylor series map with coefficients of type
// [b]T[/b]. It implements operations available only for symplectic maps.
template <class T>
class LieMap: public VpsInvMap<T> {
public:
/// Constructor.
// Construct identity with [b]nDim[/b] variables.
LieMap(int nDim);
/// Convert.
// Construct a map from the matrix [b]M[/b], which should be symplectic.
LieMap(const Matrix<T> &M);
/// Convert.
// Construct an identity map with a displacement vector [b]V[/b].
LieMap(const Vector<T> &V);
/// Conversion.
// Convert possibly non-symplectic map.
LieMap(const Vps<T> &rhs);
LieMap();
LieMap(const LieMap<T> &rhs);
~LieMap();
LieMap<T> &operator=(const LieMap<T> &);
/// Check consistency.
void check() const;
/// Lie series.
// Build the exponential series exp(:H:)Z from the Tps<T> H,
// acting on the identity map [b]Z[/b].
static LieMap<T> ExpMap(const Tps<T> &H);
/// Lie series.
// Build the exponential series exp(:H:)M from the Tps<T> H,
// acting on the map M.
static LieMap<T> ExpMap(const Tps<T> &H, const LieMap<T> &M);
private:
// The iteration limit for Lie series convergence.
static const int MAXITER = 100;
};
// Operations on Vps<T>.
// ------------------------------------------------------------------------
/// Poisson bracket.
template <class T>
Tps<T> PoissonBracket(const Tps<T> &x, const Tps<T> &y);
/// Poisson bracket.
template <class T>
Vps<T> PoissonBracket(const Tps<T> &x, const Vps<T> &y);
/// Extract LieMap<T> from stream.
template <class T>
std::istream &operator>>(std::istream &, LieMap<T> &x);
/// Insert LieMap<T> to stream.
template <class T>
std::ostream &operator<<(std::ostream &, const LieMap<T> &x);
// Implementation of global functions for class LieMap<T>.
// ------------------------------------------------------------------------
template <class T>
LieMap<T>::LieMap(): VpsInvMap<T>()
{}
template <class T>
LieMap<T>::LieMap(int nDim): VpsInvMap<T>(nDim) {
check();
}
template <class T>
LieMap<T>::LieMap(const LieMap<T> &rhs): VpsInvMap<T>(rhs)
{}
template <class T>
LieMap<T>::LieMap(const Matrix<T> &x): VpsInvMap<T>(x)
{}
template <class T>
LieMap<T>::LieMap(const Vector<T> &x): VpsInvMap<T>(x)
{}
template <class T>
LieMap<T>::LieMap(const Vps<T> &rhs):
VpsInvMap<T>(rhs) {
check();
}
template <class T> inline
LieMap<T>::~LieMap()
{}
template <class T> inline
LieMap<T> &LieMap<T>::operator=(const LieMap<T> &rhs) {
VpsInvMap<T>::operator=(rhs);
return *this;
}
template <class T>
Tps<T> PoissonBracket(const Tps<T> &x, const Tps<T> &y) {
Tps<T> z;
int nDim = x.getDimension();
for(int q = 0; q < nDim; q += 2) {
int p = q + 1;
z += x.derivative(q) * y.derivative(p) - x.derivative(p) * y.derivative(q);
}
return z;
}
template <class T>
Vps<T> PoissonBracket(const Tps<T> &x, const Vps<T> &y) {
Vps<T> z;
int nDim = x.getDimension();
for(int q = 0; q < nDim; q += 2) {
int p = q + 1;
z += x.derivative(q) * y.derivative(p) - x.derivative(p) * y.derivative(q);
}
return z;
}
template <class T>
LieMap<T> LieMap<T>::ExpMap(const Tps<T> &H) {
return ExpMap(H, VpsInvMap<T>::identity(H.getVariables()));
}
template <class T>
LieMap<T> LieMap<T>::ExpMap(const Tps<T> &H, const LieMap<T> &M) {
int nDim = H.getVariables();
if(nDim % 2 != 0) {
throw LogicalError("LieMap::ExpMap()", "Tps dimension is zero or odd.");
}
if(nDim == 0) return M;
if(nDim != M.getVariables()) {
throw LogicalError
("LieMap::ExpMap()",
"Lie map and Hamiltonian have inconsistent dimensions.");
}
LieMap<T> temp(nDim);
Vps<T> dH(nDim, nDim);
for(int i = 0; i < nDim; i += 2) {
dH[i] = - H.derivative(i + 1);
dH[i+1] = H.derivative(i);
}
for(int var = 0; var < nDim; var++) {
Tps<T> old = 0.0;
Tps<T> u = temp[var] = M[var];
for(int count = 1; temp[var] != old; count++) {
if(count >= MAXITER) {
throw ConvergenceError("LieMap::ExpMap()",
"No convergence in ExpMap()");
}
Tps<T> w = 0.0;
for(int v = 0; v < nDim; v++) w += dH[v] * u.derivative(v);
u = w / T(count);
old = temp[var];
temp[var] += u;
}
}
return temp;
}
template <class T>
std::istream &operator>>(std::istream &is, LieMap<T> &x)
{ return x.get(is); }
template <class T>
std::ostream &operator<<(std::ostream &os, const LieMap<T> &x)
{ return x.put(os); }
template <class T>
void LieMap<T>::check() const {
VpsInvMap<T>::check();
if(this->getDimension() % 2 != 0) {
throw LogicalError("LieMap::check()", "LieMap representation corrupted.");
}
}
#endif // CLASSIC_LieMap_HH
#ifndef CLASSIC_VpsInvMap_HH
#define CLASSIC_VpsInvMap_HH
// ------------------------------------------------------------------------
// $RCSfile: VpsInvMap.h,v $
// ------------------------------------------------------------------------
// $Revision: 1.1.1.1.2.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: VpsInvMap
//
// ------------------------------------------------------------------------
// Class category: Algebra
// ------------------------------------------------------------------------
//
// $Date: 2004/11/18 22:18:06 $
// $Author: jsberg $
//
// ------------------------------------------------------------------------
#include "Algebra/LUMatrix.h"
#include "Algebra/Matrix.h"
#include "Algebra/Vector.h"
#include "Algebra/VpsInvMap.h"
#include "Algebra/Tps.h"
#include "Algebra/VpsMap.h"
#include "Utilities/LogicalError.h"
#include <iosfwd>
// Class VpsInvMap
// ------------------------------------------------------------------------
/// Invertible power series map
// An invertible truncated power series map.
// The number of variables and the dimension are equal.
template <class T>
class VpsInvMap: public VpsMap<T> {
public:
/// Constructor.
// Construct identity in [b]nDim[/b] variables.
VpsInvMap(int nDim);
/// Convert.
// Throw SizeError if [b]rhs[/b] is not R**n --> R**n.
VpsInvMap(const VpsMap<T> &rhs);
/// Convert.
// Throw SizeError if [b]rhs[/b] is not R**n --> R**n.
VpsInvMap(const Vps<T> &rhs);
/// Constructor.
// Construct map of dimension [b]nDim[/b] and fill from [b]rhs[/b].
VpsInvMap(int nDim, const Tps<T> rhs[]);
/// Convert.
// The constant terms are zero.
// The linear terms are taken from [b]M[/b].
// Throw SizeError if [b]M[/b] is not square.
VpsInvMap(const Matrix<T> &M);
/// Convert.
// The constant terms are taken from [b]V[/b].
// The linear terms are the identity map.
VpsInvMap(const Vector<T> &);
VpsInvMap();
VpsInvMap(const VpsInvMap<T> &rhs);
~VpsInvMap();
VpsInvMap<T> &operator=(const VpsInvMap<T> &y);
/// Check consistency.
void check() const;
/// Inverse.
VpsInvMap<T> inverse() const;
/// Set to identity.
static VpsInvMap<T> identity(int nVar);
};
// ------------------------------------------------------------------------
///{ Global Operators on VpsInvMap<T>
/// Multiply.
template <class T>
VpsInvMap<T> operator*(const VpsInvMap<T> &x, const Tps<T> &y);
/// Multiply.
template <class T>
VpsInvMap<T> operator*(const Tps<T> &x, const VpsInvMap<T> &y);
/// Multiply.
template <class T>
VpsInvMap<T> operator*(const VpsInvMap<T> &x, const T &y);
/// Multiply.
template <class T>
VpsInvMap<T> operator*(const T &x, const VpsInvMap<T> &y);
/// Divide.
// Throw DivideError, if constant part of [b]y[/b] is zero.
template <class T>
VpsInvMap<T> operator/(const VpsInvMap<T> &x, const Tps<T> &y);
/// Divide.
// Throw DivideError, if [b]y[/b] is zero.
template <class T>
VpsInvMap<T> operator/(const VpsInvMap<T> &x, const T &y);
/// Add.
template <class T>
VpsInvMap<T> operator+(const VpsInvMap<T> &x, const VpsInvMap<T> &y);
/// Subtract.
template <class T>
VpsInvMap<T> operator-(const VpsInvMap<T> &x, const VpsInvMap<T> &y);
/// Add.
template <class T>
VpsInvMap<T> operator+(const VpsInvMap<T> &x, const Vector<T> &y);
/// Subtract.
template <class T>
VpsInvMap<T> operator-(const VpsInvMap<T> &x, const Vector<T> &y);
/// Add.
template <class T>
VpsInvMap<T> operator+(const Vector<T> &x, const VpsInvMap<T> &y);
/// Subtract.
template <class T>
VpsInvMap<T> operator-(const Vector<T> &x, const VpsInvMap<T> &y);
/// Extract from stream.
template <class T>
std::istream &operator>>(std::istream &, VpsInvMap<T> &x);
/// Insert to stream.
template <class T>
std::ostream &operator<<(std::ostream &, const VpsInvMap<T> &x);
/// Multiply.
template <class T> VpsInvMap<T>
operator*(const Matrix<T> &x, const VpsInvMap<T> &y);
// Implementation of global functions for class Vps<T>.
// ------------------------------------------------------------------------
template <class T> inline
VpsInvMap<T> operator*(const VpsInvMap<T> &x, const Tps<T> &y)
{ VpsInvMap<T> z(x); return z *= y; }
template <class T> inline
VpsInvMap<T> operator*(const Tps<T> &x, const VpsInvMap<T> &y)
{ VpsInvMap<T> z(x); return z *= y; }
template <class T> inline
VpsInvMap<T> operator*(const VpsInvMap<T> &x, const T &y)
{ VpsInvMap<T> z(x); return z *= y; }
template <class T> inline
VpsInvMap<T> operator*(const T &x, const VpsInvMap<T> &y)
{ VpsInvMap<T> z(x); return z *= y; }
template <class T> inline
VpsInvMap<T> operator/(const VpsInvMap<T> &x, const Tps<T> &y)
{ VpsInvMap<T> z(x); return z /= y; }
template <class T> inline
VpsInvMap<T> operator/(const VpsInvMap<T> &x, const T &y)
{ VpsInvMap<T> z(x); return z /= y; }
template <class T> inline
VpsInvMap<T> operator+(const VpsInvMap<T> &x, const VpsInvMap<T> &y)
{ VpsInvMap<T> z(x); return z += y; }
template <class T> inline
VpsInvMap<T> operator-(const VpsInvMap<T> &x, const VpsInvMap<T> &y)
{ VpsInvMap<T> z(x); return z -= y; }
template <class T> inline
VpsInvMap<T> operator+(const VpsInvMap<T> &x, const Vector<T> &y)
{ VpsInvMap<T> z(x); return z += y; }
template <class T> inline
VpsInvMap<T> operator-(const VpsInvMap<T> &x, const Vector<T> &y)
{ VpsInvMap<T> z(x); return z -= y; }
template <class T> inline
VpsInvMap<T> operator+(const Vector<T> &x, const VpsInvMap<T> &y)
{ VpsInvMap<T> z(y); return z += x; }
template <class T> inline
VpsInvMap<T> operator-(const Vector<T> &x, const VpsInvMap<T> &y)
{ VpsInvMap<T> z(- y); return z += x; }
template <class T> inline
VpsInvMap<T> operator*(const Matrix<T> &x, const VpsInvMap<T> &y)
{ return y.substituteInto(x); }
template <class T> inline
std::istream &operator>>(std::istream &is, VpsInvMap<T> &x)
{ return x.get(is); }
template <class T> inline
std::ostream &operator<<(std::ostream &os, const VpsInvMap<T> &x)
{ return x.put(os); }
// Class VpsInvMap, implementation.
// ------------------------------------------------------------------------
template <class T>
VpsInvMap<T>::VpsInvMap():
VpsMap<T>()
{}
template <class T>
VpsInvMap<T>::VpsInvMap(int nDim):
VpsMap<T>(nDim, nDim) {
for(int i = 0; i < nDim; i++) this->data[i] = Tps<T>::makeVariable(nDim, i);
}
template <class T>
VpsInvMap<T>::VpsInvMap(const VpsInvMap<T> &rhs):
VpsMap<T>(rhs)
{}
template <class T>
VpsInvMap<T>::VpsInvMap(const VpsMap<T> &rhs):
VpsMap<T>(rhs)
{}
template <class T>
VpsInvMap<T>::VpsInvMap(const Vps<T> &rhs):
VpsMap<T>(rhs)
{}
template <class T>
VpsInvMap<T>::VpsInvMap(int nDim, const Tps<T> rhs[]):
VpsMap<T>(nDim, rhs) {
check();
}
template <class T>
VpsInvMap<T>::VpsInvMap(const Matrix<T> &x):
VpsMap<T>(x) {
if(x.nrows() != x.ncols()) {
throw LogicalError("VpsInvMap::VpsInvMap()", "Matrix is not n by n.");
}
}
template <class T>
VpsInvMap<T>::VpsInvMap(const Vector<T> &x):
VpsMap<T>(x) {
}
template <class T>
VpsInvMap<T>::~VpsInvMap() {
}
template <class T>
VpsInvMap<T> &VpsInvMap<T>::operator=(const VpsInvMap &y) {
VpsMap<T>::operator=(y);
return *this;
}
template <class T>
void VpsInvMap<T>::check() const {
Vps<T>::check();
if(this->getDimension() != this->variables) {
throw LogicalError("VpsInvMap::VpsInvMap()", "Map is not n->n.");
}
}
template <class T>
VpsInvMap<T> VpsInvMap<T>::inverse() const {
// Separate out linear terms and change sign of all other terms:
// M := A(1), the linear part of map A,
// map1 := - (A - A(1)).
check();
Vector<T> vec = this->constantTerm();
LUMatrix<T> lu(this->linearTerms());
Matrix<T> mat(lu.inverse());
VpsInvMap<T> map1 = - this->filter(2, Tps<T>::EXACT);
// Initialize map B = M**(-1) (with truncation at maxOrder).
VpsInvMap<T> B(mat);
B -= mat * vec;
// Compute inverse order by order.
for(int maxOrder = 2; maxOrder <= this->getTruncOrder(); maxOrder++) {
VpsInvMap<T> map2 = map1.substitute(B, maxOrder);
B = mat * (map2 + identity(this->variables) - vec);
}
return VpsInvMap<T>(B);
}
template <class T>
VpsInvMap<T> VpsInvMap<T>::identity(int nDim) {
return VpsInvMap<T>(nDim);
}
#endif // CLASSIC_VpsInvMap_HH
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment