Commit 7f4c83d4 authored by snuverink_j's avatar snuverink_j

remove unused (after action removal) classes

parent 96104fa0
set (_SRCS
QRSolver.cpp
)
include_directories (
${CMAKE_CURRENT_SOURCE_DIR}
)
add_opal_sources(${_SRCS})
set (HDRS
QRSolver.h
)
install (FILES ${HDRS} DESTINATION "${CMAKE_INSTALL_PREFIX}/include/Algebra")
# vi: set et ts=4 sw=4 sts=4:
# Local Variables:
# mode: cmake
# cmake-tab-width: 4
# indent-tabs-mode: nil
# require-final-newline: nil
# End:
// ------------------------------------------------------------------------
// $RCSfile: QRSolver.cpp,v $
// ------------------------------------------------------------------------
// $Revision: 1.1.1.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: QRSolver
// This class implements routines for the least-square solution of
// systems of linear equations.
// Given an m by n matrix A, an n by n diagonal matrix D,
// and an m-vector B, two problem can be solved:
//
// 1. Solve the the system
// A*X = B,
// in the least squares sense. The first step to solve this problem is
// QRSolver solver(A, pivot);
// The second step is then
// solver.solveR(X);
//
// 2. Solve the the two systems
// A*X = B, D*X = 0,
// in the least squares sense. The first step to solve this problem is
// QRSolver solver(A, pivot);
// The second step is
// solver.solveS(D, X);
// it can be repeated as many times as required for different diagonal
// matrices D.
//
// In both cases, the method
// solver.getColNorm(C);
// can be called to return the original column norms of A.
//
// ------------------------------------------------------------------------
//
// $Date: 2000/03/27 09:33:35 $
// $Author: Andreas Adelmann $
//
// ------------------------------------------------------------------------
#include "Algebra/QRSolver.h"
#include "Algebra/Matrix.h"
#include "Algebra/Vector.h"
#include <cassert>
#include <cmath>
// This method uses Householder transformations with optional column
// pivoting to compute a QR-factorization of the m-by-n matrix A.
// QRSolver determines an orthogonal matrix Q, a permutation matrix P,
// and an upper trapezoidal matrix R with diagonal elements of
// nonincreasing magnitude, such that A*P = Q*R. The Householder
// transformation for column k, k = 1, 2, ..., min(m,n), is of the form
//
// T
// I - (1/U(k))*U*U
//
// where U has zeros in the first k-1 positions. The form of this
// transformation and the method of pivoting first appeared in the
// corresponding LINPACK subroutine.
//
QRSolver::QRSolver
(const Matrix<double> &A, const Vector<double> &B, bool pivot):
n(A.ncols()), m(A.nrows()), Q(A), R(n, n, 0.0), QtB(B), S(n, n, 0.0),
column_norm(n, 0.0), piv_col(n, 0) {
assert(B.size() >= m);
Vector<double> dot(n);
// Compute the initial and internal column norms.
for(int j = 0; j < n; ++j) {
double temp = 0.0;
for(int i = 0; i < m; ++i) temp += Q(i, j) * Q(i, j);
dot(j) = column_norm(j) = sqrt(temp);
piv_col(j) = j;
}
// Reduce A to R using Householder transformations.
// Up to min(m, n) columns are used.
for(int j = 0; j < std::min(m, n); ++j) {
if(pivot) {
// Bring column of largest norm into pivot position.
int kmax = j;
for(int k = j + 1; k < n; ++k) {
if(dot(k) > dot(kmax)) kmax = k;
}
if(kmax != j) {
std::swap_ranges(Q.col_begin(j), Q.col_end(j), Q.col_begin(kmax));
dot(kmax) = dot(j);
std::swap(piv_col(j), piv_col(kmax));
}
}
// Compute the Householder transformation which reduces the j-th
// column of A to a multiple of the j-th unit vector.
double unorm = 0.0;
for(int k = j; k < m; ++k) unorm += Q(k, j) * Q(k, j);
unorm = (Q(j, j) > 0.0) ? sqrt(unorm) : - sqrt(unorm);
// Build column j of R.
for(int i = 0; i < j; ++i) {
R(i, j) = Q(i, j);
Q(i, j) = 0.0;
}
R(j, j) = - unorm;
if(unorm != 0.0) {
// Modify column j of Q.
for(int i = j; i < m; ++i) {
Q(i, j) /= unorm;
}
Q(j, j) += 1.0;
// Transform remaining columns of Q.
for(int k = j + 1; k < n; ++k) {
// Apply transformation to column k.
double sum = 0.0;
for(int i = j; i < m; ++i) sum += Q(i, j) * Q(i, k);
sum /= Q(j, j);
for(int i = j; i < m; ++i) Q(i, k) -= sum * Q(i, j);
if(pivot) {
// Recompute column norm.
double sum = 0.0;
for(int i = j + 1; i < m; ++i) sum += Q(i, k) * Q(i, k);
dot(k) = sqrt(sum);
}
}
// Transform r.h.s.
if(Q(j, j) != 0.0) {
double sum = 0.0;
for(int i = j; i < m; ++i) sum += Q(i, j) * QtB(i);
sum /= Q(j, j);
for(int i = j; i < m; ++i) QtB(i) -= Q(i, j) * sum;
}
}
}
}
QRSolver::~QRSolver()
{}
void QRSolver::solveR(Vector<double> &X) const {
solve(R, QtB, X);
}
void QRSolver::solveS(const Array1D<double> &D, double p, Vector<double> &X) {
Vector<double> Z(m);
Array1D<double> M(D);
double root = sqrt(p);
for(int i = 0; i < n; ++i) M(i) *= root;
R_to_S(R, M, S, Z);
solve(S, Z, X);
}
void QRSolver::solveRT(Vector<double> &V) const {
solveT(R, V);
}
void QRSolver::solveST(Vector<double> &V) const {
solveT(S, V);
}
void QRSolver::getColNorm(Array1D<double> &norms) const {
norms = column_norm;
}
void QRSolver::R_to_S(const Matrix<double> &R, const Array1D<double> &D,
Matrix<double> &S, Vector<double> &Z) {
// Eliminate the diagonal matrix D using Givens rotations.
assert(D.size() >= n);
S = R;
Z = QtB;
for(int j = 0; j < n; ++j) {
// Prepare the row of D to be eliminated, locating the diagonal
// element using P from the QR-factorization.
int l = piv_col(j);
if(D(l) != 0.0) {
Vector<double> row_D(n, 0.0);
row_D(j) = D(l);
// The transformations to eliminate the row of D; modify only
// a single element of transpose(Q)*B beyond the first n,
// which is initially 0.
double QtBpj = 0.0;
for(int k = j; k < n; ++k) {
// Determine a Givens rotation which eliminates the appropriate
// element in the current row of D.
if(row_D(k) != 0.0) {
double t = sqrt(S(k, k) * S(k, k) + row_D(k) * row_D(k));
double s = row_D(k) / t;
double c = S(k, k) / t;
// Compute the modified diagonal element of S and the modified
// element of ((Q transpose)*B, 0).
S(k, k) = t;
t = c * Z(k) + s * QtBpj;
QtBpj = -s * Z(k) + c * QtBpj;
Z(k) = t;
// Accumulate the transformation in the row of S.
for(int i = k + 1; i < n; ++i) {
t = c * S(k, i) + s * row_D(i);
row_D(i) = -s * S(k, i) + c * row_D(i);
S(k, i) = t;
}
}
}
}
}
}
void QRSolver::solve(const Matrix<double> &R, const Vector<double> &QtB,
Vector<double> &X) const {
// Copy QtB to the solution.
Vector<double> Z(QtB);
// Solve the triangular system for Z.
// If R is singular, then obtain a least squares solution.
int rank = n;
for(int j = 0; j < n; ++j) {
if(R(j, j) == 0.0 && rank == n) rank = j;
if(rank < n) Z(j) = 0.0;
}
for(int j = rank; j-- > 0;) {
double sum = Z(j);
for(int i = j + 1; i < rank; ++i) {
sum -= R(j, i) * Z(i);
}
Z(j) = sum / R(j, j);
}
// Permute the components of Z back to components of X.
X = Vector<double>(n, 0.0);
for(int j = 0; j < n; ++j) {
X(piv_col(j)) = Z(j);
}
}
void QRSolver::solveT(const Matrix<double> &R, Vector<double> &V) const {
assert(V.size() == n);
Vector<double> X(n);
for(int i = 0; i < n; ++i) X(i) = V(piv_col(i));
for(int i = 0; i < n; ++i) {
if(R(i, i) == 0.0) break;
for(int j = 0; j < i; ++j) X(j) -= R(i, j) * V(i);
X(i) /= R(i, i);
}
V = X;
}
#ifndef OPAL_QRSolver_HH
#define OPAL_QRSolver_HH 1
// ------------------------------------------------------------------------
// $RCSfile: QRSolver.h,v $
// ------------------------------------------------------------------------
// $Revision: 1.1.1.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: QRSolver
//
// ------------------------------------------------------------------------
//
// $Date: 2000/03/27 09:33:35 $
// $Author: Andreas Adelmann $
//
// ------------------------------------------------------------------------
#include "Algebra/Matrix.h"
#include "Algebra/Vector.h"
#include <cstddef>
// Class QRSolver
// ------------------------------------------------------------------------
/// Least-square solution of systems of linear equations.
// Given an m by n matrix A, an n by n diagonal matrix D,
// and an m-vector B, two problem can be solved:
// [ol][li]
// Solve the the system $A*X = B$ in the least squares sense.
// The first step to solve this problem is:
// [pre]
// QRSolver solver(A, pivot);
// [/pre]
// The second step is then
// [pre]
// solver.solveR(X);
// [/pre]
// [li]
// Solve the the two systems $A*X = B, D*X = 0$ in the least squares sense.
// The first step to solve this problem is
// [pre]
// QRSolver solver(A, pivot);
// [/pre]
// The second step is
// [pre]
// solver.solveS(D, X);
// [/pre]
// The second step can be repeated as many times as required for different
// diagonal matrices $D$.
// [/ol]
// In both cases, the method
// [pre]
// solver.getColNorm(C);
// [/pre]
// can be called to return the original column norms of $A$.
class QRSolver {
public:
/// Constructor.
// Determine the QR-factorisation $A = Q*R$ of the matrix $A$
// and transform the right-hand side $B$. If [b]pivot[/b] is
// true, then pivot search is done.
// [p]
// This method uses Householder transformations with optional column
// pivoting to compute a QR-factorization of the m-by-n matrix $A$.
// It determines an orthogonal matrix $Q$, a permutation matrix $P$,
// and an upper trapezoidal matrix $R$ with diagonal elements of
// nonincreasing magnitude, such that $A*P = Q*R$. The Householder
// transformation for column $k$ is of the form $I - (1/U(k))*U*U^T$,
// where $U$ has zeros in the first $k-1$ positions.
// The form of this transformation and the method of pivoting first
// appeared in the corresponding LINPACK subroutine.
QRSolver(const Matrix<double> &A, const Vector<double> &B, bool pivot);
~QRSolver();
/// Solution of $A*X = B$ in the least-squares sense.
void solveR(Vector<double> &X) const;
/// Solution of $A*X = B, D*X = 0$ in the least-squares sense.
void solveS(const Array1D<double> &D, double p, Vector<double> &X);
/// Pre-multiply the vector $V$ by $R.transpose()^{-1}$.
void solveRT(Vector<double> &V) const;
/// Pre-multiply the vector $V$ by $S.transpose()^{-1}$.
// Requires prior execution of solveS.
void solveST(Vector<double> &V) const;
/// Return the original column norms of the matrix $A$.
void getColNorm(Array1D<double> &) const;
private:
// Transform R to S, given the diagonal matrix D.
void R_to_S(const Matrix<double> &R,
const Array1D<double> &D,
Matrix<double> &S,
Vector<double> &Z);
// Solution of R*X = QtB.
void solve(const Matrix<double> &R,
const Vector<double> &QtB,
Vector<double> &X) const;
// Solution of R.transpose()**(-1) * V
void solveT(const Matrix<double> &R, Vector<double> &V) const;
// The dimensions of the problem.
int n; // number of columns.
int m; // number of rows.
// The lower trapezoidal part of the matrix Q contains the n vectors
// defining a factored form of the orthogonal matrix Q.
Array2D<double> Q;
// The upper-trapezoidal matrix R contains the matrix R.
Matrix<double> R;
// The transformed r.h.s. Q.transpose() * B.
Vector<double> QtB;
// The upper-trapezoidal matrix S contains the matrix S.
Matrix<double> S;
// The vector column_norm contains the original column norms of M.
Array1D<double> column_norm;
// The vector piv_col defines the permutation matrix P.
// Column j of P is column piv_col(j) of the identity matrix.
Array1D<int> piv_col;
};
#endif // OPAL_QRSolver_HH
......@@ -3,11 +3,9 @@ set (_SRCS
Ctunes.cpp
IndexMap.cpp
Hamiltonian.cpp
LieMapper.cpp
lomb.cpp
NilTracker.cpp
MapAnalyser.cpp
MPSplitIntegrator.cpp
MultiBunchHandler.cpp
OrbitThreader.cpp
ParallelCyclotronTracker.cpp
......@@ -15,7 +13,6 @@ set (_SRCS
StepSizeConfig.cpp
ThickMapper.cpp
ThickTracker.cpp
TransportMapper.cpp
)
include_directories (
......@@ -29,10 +26,8 @@ set (HDRS
Ctunes.h
Hamiltonian.h
IndexMap.h
LieMapper.h
lomb.h
MapAnalyser.h
MPSplitIntegrator.h
MultiBunchHandler.h
NilTracker.h
OrbitThreader.h
......@@ -41,7 +36,6 @@ set (HDRS
StepSizeConfig.h
ThickMapper.h
ThickTracker.h
TransportMapper.h
)
install (FILES ${HDRS} DESTINATION "${CMAKE_INSTALL_PREFIX}/include/Algorithms")
......
This diff is collapsed.
#ifndef OPAL_LieMapper_HH
#define OPAL_LieMapper_HH
// ------------------------------------------------------------------------
// $RCSfile: LieMapper.h,v $
// ------------------------------------------------------------------------
// $Revision: 1.1.1.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: LieMapper
//
// ------------------------------------------------------------------------
//
// $Date: 2000/03/27 09:33:36 $
// $Author: Andreas Adelmann $
//
// ------------------------------------------------------------------------
#include "Algorithms/Mapper.h"
#include "FixedAlgebra/DragtFinnMap.h"
class BMultipoleField;
class PlanarArcGeometry;
// Class LieMapper
// ------------------------------------------------------------------------
/// Build a Lie-algebraic map using a finite-length lens for each elements.
// All maps are factored in ascending Dragt-Finn order.
// [p]
// Phase space coordinates numbering:
// [tab 3 b]
// [row]number [&]name [&]unit [/row]
// [row]0 [&]$x$ [&]metres [/row]
// [row]1 [&]$p_x/p_r$ [&]1 [/row]
// [row]2 [&]$y$ [&]metres [/row]
// [row]3 [&]$p_y/p_r$ [&]1 [/row]
// [row]4 [&]$v*delta_t$ [&]metres [/row]
// [row]5 [&]$delta_p/p_r$ [&]1 [/row]
// [/tab][p]
// Where $p_r$ is the constant reference momentum defining the reference
// frame velocity, $m$ is the rest mass of the particles, and $v$ is the
// instantaneous velocity of the particle.
// [p]
// Other units used:
// [tab 2 b]
// [row]quantity [&]unit [/row]
// [row]reference momentum [&]electron-volts [/row]
// [row]velocity [&]metres/second [/row]
// [row]accelerating voltage [&]volts [/row]
// [row]separator voltage [&]volts [/row]
// [row]frequencies [&]hertz [/row]
// [row]phase lags [&]$2*pi$ [/row]
// [/tab][p]
// All elements are represented by maps for finite-length elements.
// Several important pieces are still MISSING from this algorithm.
class LieMapper: public AbstractMapper {
public:
/// Constructor.
// The beam line to be tracked is "bl".
// The particle reference data are taken from "data".
// If [b]revBeam[/b] is true, the beam runs from s = C to s = 0.
// If [b]revTrack[/b] is true, we track against the beam.
LieMapper(const Beamline &bl, const PartData &data,
bool backBeam, bool backTrack, int order);
virtual ~LieMapper();
/// Return the linear part of the accumulated map.
virtual void getMap(LinearMap<double, 6> &) const;
/// Return the accumulated map.
virtual void getMap(FVps<double, 6> &) const;
/// Return the full map accumulated so far.
virtual void getMap(DragtFinnMap<3> &) const;
/// Reset the linear part of the accumulated map for restart.
virtual void setMap(const LinearMap<double, 6> &);
/// Reset the accumulated map for restart.
virtual void setMap(const FVps<double, 6> &);
/// Reset the full map for restart.
virtual void setMap(const DragtFinnMap<3> &);
/// Apply the algorithm to a BeamBeam.
virtual void visitBeamBeam(const BeamBeam &);
/// Apply the algorithm to a BeamStripping.
virtual void visitBeamStripping(const BeamStripping &);
/// Apply the algorithm to a Collimator.
virtual void visitCCollimator(const CCollimator &);
/// Apply the algorithm to a Component.
virtual void visitComponent(const Component &);
/// Apply the algorithm to a Corrector.
virtual void visitCorrector(const Corrector &);
/// Apply the algorithm to a Degrader
virtual void visitDegrader(const Degrader &);
/// Apply the algorithm to a Diagnostic.
virtual void visitDiagnostic(const Diagnostic &);
/// Apply the algorithm to a Drift.
virtual void visitDrift(const Drift &);
/// Apply the algorithm to a flexible collimator
virtual void visitFlexibleCollimator(const FlexibleCollimator &);
/// Apply the algorithm to a Lambertson.
virtual void visitLambertson(const Lambertson &);
/// Apply the algorithm to a Marker.
virtual void visitMarker(const Marker &);
/// Apply the algorithm to a Monitor.
virtual void visitMonitor(const Monitor &);
/// Apply the algorithm to a Multipole.
virtual void visitMultipole(const Multipole &);
/// Apply the algorithm to a Patch.
virtual void visitPatch(const Patch &);
/// Apply the algorithm to a Probe.
virtual void visitProbe(const Probe &);
/// Apply the algorithm to a RBend.
virtual void visitRBend(const RBend &);
/// Apply the algorithm to a RFCavity.
virtual void visitRFCavity(const RFCavity &);
/// Apply the algorithm to a RFQuadrupole.
virtual void visitRFQuadrupole(const RFQuadrupole &);
/// Apply the algorithm to a SBend.
virtual void visitSBend(const SBend &);
/// Apply the algorithm to a Separator.
virtual void visitSeparator(const Separator &);
/// Apply the algorithm to a Septum.
virtual void visitSeptum(const Septum &);
/// Apply the algorithm to a Solenoid.
virtual void visitSolenoid(const Solenoid &);
/// Apply the algorithm to a ParallelPlate.
virtual void visitParallelPlate(const ParallelPlate &);
private:
// Not implemented.
LieMapper();
LieMapper(const LieMapper &);
void operator=(const LieMapper &);
// Apply drift length.
// Propagate current map through a drift.
void applyDrift(double length);
// Transforms fringing fields.
void applyEntranceFringe(double edge, double curve,
const BMultipoleField &field, double scale);
void applyExitFringe(double edge, double curve,
const BMultipoleField &field, double scale);
/// Apply transform.
// Propagate current map through a geometric transformation.
// Linear approximation for the time being.
virtual void applyTransform(const Euclid3D &, double refLength = 0.0);
// The Lie map being accumulated.
DragtFinnMap<3> itsMap;
// The desired map order.
int itsOrder;
};
#endif // OPAL_LieMapper_HH
\ No newline at end of file
// ------------------------------------------------------------------------
// $RCSfile: MPSplitIntegrator.cpp,v $
// ------------------------------------------------------------------------
// $Revision: 1.1.1.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: MPSplitIntegrator
// A MPSplitIntegrator performs integration through an element using two
// thin lenses of force 1/2, one placed at 1/6 and the other at 5/6 of
// the length respectively.
//
// ------------------------------------------------------------------------
//
// $Date: 2000/03/27 09:33:36 $
// $Author: Andreas Adelmann $