Commit 0a2b939b authored by Claus Kleinwort's avatar Claus Kleinwort
Browse files

Fixes for geometric constraint with external measurement and internal iterations (setZero)

git-svn-id: http://svnsrv.desy.de/public/GeneralBrokenLines/trunk@155 281f6f2b-e318-4fd1-8bce-1a4ba7aab212
parent 79853f6f
......@@ -7,7 +7,7 @@ PROJECT(GBL)
# project version
SET( ${PROJECT_NAME}_VERSION_MAJOR 2 )
SET( ${PROJECT_NAME}_VERSION_MINOR 3 )
SET( ${PROJECT_NAME}_VERSION_PATCH 0 )
SET( ${PROJECT_NAME}_VERSION_PATCH 1 )
# make life easier and simply use the ilcsoft default settings
# load default ilcsoft settings (install prefix, build type, rpath, etc.)
......
......@@ -79,6 +79,7 @@ public:
virtual ~BorderedBandMatrix();
void resize(unsigned int nSize, unsigned int nBorder = 1,
unsigned int nBand = 5);
void setZero();
void solveAndInvertBorderedBand(const VVector &aRightHandSide,
VVector &aSolution);
void addBlockMatrix(double aWeight,
......
......@@ -61,7 +61,7 @@ typedef Eigen::Matrix<double, 5, 5> Matrix5d;
*/
class GblPoint {
public:
GblPoint(const Matrix5d &aJacobian);
GblPoint(const Matrix5d &aJacobian, unsigned int numMeasReserve = 0);
GblPoint(const GblPoint&) = default;
GblPoint& operator=(const GblPoint&) = default;
GblPoint(GblPoint&&) = default;
......
......@@ -46,6 +46,7 @@ public:
VVector(const VVector &aVector);
virtual ~VVector();
void resize(const unsigned int nRows);
void setZero();
VVector getVec(unsigned int len, unsigned int start = 0) const;
void putVec(const VVector &aVector, unsigned int start = 0);
inline double &operator()(unsigned int i);
......@@ -66,6 +67,7 @@ public:
VMatrix(const VMatrix &aMatrix);
virtual ~VMatrix();
void resize(const unsigned int nRows, const unsigned int nCols);
void setZero();
VMatrix transpose() const;
inline double &operator()(unsigned int i, unsigned int j);
inline double operator()(unsigned int i, unsigned int j) const;
......@@ -88,6 +90,7 @@ public:
VSymMatrix(const unsigned int nRows = 0);
virtual ~VSymMatrix();
void resize(const unsigned int nRows);
void setZero();
unsigned int invert();
inline double &operator()(unsigned int i, unsigned int j);
inline double operator()(unsigned int i, unsigned int j) const;
......
......@@ -58,6 +58,13 @@ void BorderedBandMatrix::resize(unsigned int nSize, unsigned int nBorder,
theBand.resize((nBand + 1), numCol);
}
/// Set content to 0.
void BorderedBandMatrix::setZero() {
theBorder.setZero();
theMixed.setZero();
theBand.setZero();
}
/// Add symmetric block matrix.
/**
* Add (extended) block matrix defined by 'aVector * aWeight * aVector.T'
......@@ -68,8 +75,8 @@ void BorderedBandMatrix::resize(unsigned int nSize, unsigned int nBorder,
* \param aVector [in] Vector
*/
void BorderedBandMatrix::addBlockMatrix(double aWeight,
const std::vector<unsigned int>* anIndex,
const std::vector<double>* aVector) {
const std::vector<unsigned int> *anIndex,
const std::vector<double> *aVector) {
int nBorder = numBorder;
for (unsigned int i = 0; i < anIndex->size(); ++i) {
int iIndex = (*anIndex)[i] - 1; // anIndex has to be sorted
......@@ -102,7 +109,7 @@ void BorderedBandMatrix::addBlockMatrix(double aWeight,
* \param aVector [in] Vector
*/
void BorderedBandMatrix::addBlockMatrix(double aWeight, unsigned int aSize,
unsigned int* anIndex, double* aVector) {
unsigned int *anIndex, double *aVector) {
int nBorder = numBorder;
for (unsigned int i = 0; i < aSize; ++i) {
int iIndex = anIndex[i] - 1; // anIndex has to be sorted
......@@ -158,7 +165,7 @@ MatrixXd BorderedBandMatrix::getBlockMatrix(
* \param anIndex [in] Array of rows/colums to be used
*/
MatrixXd BorderedBandMatrix::getBlockMatrix(unsigned int aSize,
unsigned int* anIndex) const {
unsigned int *anIndex) const {
MatrixXd aMatrix(aSize, aSize);
int nBorder = numBorder;
......
......@@ -37,10 +37,12 @@ namespace gbl {
/**
* Create point on (initial) trajectory. Needs transformation jacobian from previous point.
* \param [in] aJacobian Transformation jacobian from previous point
* \param [in] numMeasReserve number of measurements to reserve (space for)
*/
GblPoint::GblPoint(const Matrix5d &aJacobian) :
GblPoint::GblPoint(const Matrix5d &aJacobian, unsigned int numMeasReserve) :
theLabel(0), theOffset(0), p2pJacobian(aJacobian), scatFlag(
false) {
theMeasurements.reserve(numMeasReserve);
}
#ifdef GBL_EIGEN_SUPPORT_ROOT
......
......@@ -1113,7 +1113,9 @@ void GblTrajectory::getResAndErr(unsigned int aData, double &aResidual,
void GblTrajectory::buildLinearEquationSystem() {
unsigned int nBorder = numCurvature + numLocals;
theVector.resize(numParameters);
theVector.setZero();
theMatrix.resize(numParameters, nBorder);
theMatrix.setZero();
double aValue, aWeight;
unsigned int *indLocal;
double *derLocal;
......@@ -1139,8 +1141,9 @@ void GblTrajectory::buildLinearEquationSystem() {
* Generate data (blocks) from measurements, kinks, external seed and measurements.
*
* \exception 10 : inner transformation matrix with invalid number of rows (valid are 5=kinematic or 2=geometric constraint)
* \exception 11 : inner transformation matrix with to few columns (must be >= number of rows)
* \exception 11 : inner transformation matrix with too few columns (must be >= number of rows)
* \exception 12 : inner transformation matrices with varying sizes
* \exception 13 : too many external derivatives (must be <= number of columns of inner transformation matrix)
*/
void GblTrajectory::prepare() {
unsigned int nDim = theDimension.size();
......@@ -1406,10 +1409,18 @@ void GblTrajectory::prepare() {
// external measurements
unsigned int nExt = externalMeasurements.rows();
if (nExt > 0) {
std::vector<unsigned int> index(numCurvature);
std::vector<double> derivatives(numCurvature);
unsigned int nInnerTransCols = innerTransformations[0].cols();
unsigned int nExtDer = externalDerivatives.cols();
if (nExtDer > nInnerTransCols) {
std::cout
<< " GblTrajectory::prepare external measurement with too many derivatives: "
<< nExtDer << ", defined: " << nInnerTransCols << std::endl;
throw 13;
}
std::vector<unsigned int> index(nExtDer);
std::vector<double> derivatives(nExtDer);
for (unsigned int iExt = 0; iExt < nExt; ++iExt) {
for (unsigned int iCol = 0; iCol < numCurvature; ++iCol) {
for (unsigned int iCol = 0; iCol < nExtDer; ++iCol) {
index[iCol] = numLocals + iCol + 1;
derivatives[iCol] = externalDerivatives(iExt, iCol);
}
......
......@@ -58,6 +58,11 @@ void VMatrix::resize(const unsigned int nRows, const unsigned int nCols) {
theVec.resize(nRows * nCols);
}
/// Set content to 0.
void VMatrix::setZero() {
std::fill(theVec.begin(), theVec.end(), 0.);
}
/// Get transposed matrix.
/**
* \return Transposed matrix.
......@@ -145,7 +150,7 @@ VMatrix VMatrix::operator+(const VMatrix &aMatrix) const {
}
/// Assignment Matrix=Matrix.
VMatrix &VMatrix::operator=(const VMatrix &aMatrix) {
VMatrix& VMatrix::operator=(const VMatrix &aMatrix) {
if (this != &aMatrix) { // Gracefully handle self assignment
numRows = aMatrix.getNumRows();
numCols = aMatrix.getNumCols();
......@@ -177,6 +182,11 @@ void VSymMatrix::resize(const unsigned int nRows) {
theVec.resize((nRows * nRows + nRows) / 2);
}
/// Set content to 0.
void VSymMatrix::setZero() {
std::fill(theVec.begin(), theVec.end(), 0.);
}
/// Get number of rows (= number of colums).
/**
* \return Number of rows.
......@@ -264,6 +274,11 @@ void VVector::resize(const unsigned int nRows) {
theVec.resize(nRows);
}
/// Set content to 0.
void VVector::setZero() {
std::fill(theVec.begin(), theVec.end(), 0.);
}
/// Get part of vector.
/**
* \param [in] len Length of part.
......@@ -318,7 +333,7 @@ VVector VVector::operator-(const VVector &aVector) const {
}
/// Assignment Vector=Vector.
VVector &VVector::operator=(const VVector &aVector) {
VVector& VVector::operator=(const VVector &aVector) {
if (this != &aVector) { // Gracefully handle self assignment
numRows = aVector.getNumRows();
theVec.resize(numRows);
......
Supports Markdown
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