Commit 9b8c6a9c authored by kraus's avatar kraus

fixing yves patch

parent cea68850
/*
/*
* Copyright (c) 2015, Chris Rogers
* All rights reserved.
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
* 1. Redistributions of source code must retain the above copyright notice,
* this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
* 3. Neither the name of STFC nor the names of its contributors may be used to
* endorse or promote products derived from this software without specific
* 3. Neither the name of STFC nor the names of its contributors may be used to
* endorse or promote products derived from this software without specific
* prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
* POSSIBILITY OF SUCH DAMAGE.
*/
......@@ -61,23 +61,23 @@ void MMatrix<m_complex>::delete_matrix()
template <>
MMatrix<double>& MMatrix<double>::operator= (const MMatrix<double>& mm)
{
{
if (&mm == this) return *this;
delete_matrix();
if(!mm._matrix) { _matrix = NULL; return *this; }
_matrix = gsl_matrix_alloc( mm.num_row(), mm.num_col() );
gsl_matrix_memcpy((gsl_matrix*)_matrix, (const gsl_matrix*)mm._matrix);
_matrix = gsl_matrix_alloc( mm.num_row(), mm.num_col() );
gsl_matrix_memcpy((gsl_matrix*)_matrix, (const gsl_matrix*)mm._matrix);
return *this;
}
template <>
MMatrix<m_complex>& MMatrix<m_complex>::operator= (const MMatrix<m_complex>& mm)
{
{
if (&mm == this) return *this;
delete_matrix();
if(!mm._matrix) { _matrix = NULL; return *this; }
_matrix = gsl_matrix_complex_alloc( mm.num_row(), mm.num_col() );
gsl_matrix_complex_memcpy((gsl_matrix_complex*)_matrix, (const gsl_matrix_complex*)mm._matrix);
_matrix = gsl_matrix_complex_alloc( mm.num_row(), mm.num_col() );
gsl_matrix_complex_memcpy((gsl_matrix_complex*)_matrix, (const gsl_matrix_complex*)mm._matrix);
return *this;
}
......@@ -147,7 +147,7 @@ template MMatrix<m_complex> MMatrix<m_complex>::Diagonal(size_t i, m_complex dia
template <class Tmplt>
MMatrix<Tmplt>::~MMatrix()
{
delete_matrix();
delete_matrix();
}
template MMatrix<double>::~MMatrix();
template MMatrix<m_complex>::~MMatrix();
......@@ -223,7 +223,7 @@ template MMatrix<m_complex> MMatrix<m_complex>::inverse() const;
template <>
void MMatrix<m_complex>::invert()
void MMatrix<m_complex>::invert()
{
if(num_row() != num_col()) throw(GeneralClassicException("MMatrix::invert()", "Attempt to get inverse of non-square matrix"));
gsl_permutation* perm = gsl_permutation_alloc(num_row());
......@@ -232,7 +232,7 @@ void MMatrix<m_complex>::invert()
gsl_linalg_complex_LU_decomp ( (gsl_matrix_complex*)lu._matrix, perm, &s);
gsl_linalg_complex_LU_invert ( (gsl_matrix_complex*)lu._matrix, perm, (gsl_matrix_complex*)_matrix);
gsl_permutation_free( perm );
for(size_t i=1; i<=num_row(); i++)
for(size_t i=1; i<=num_row(); i++)
for(size_t j=1; j<=num_row(); j++)
if(operator()(i,j) != operator()(i,j)) throw(GeneralClassicException("MMatrix::invert()", "Failed to invert matrix - singular?"));
}
......@@ -247,7 +247,7 @@ void MMatrix<double>::invert()
gsl_linalg_LU_decomp ( (gsl_matrix*)lu._matrix, perm, &s);
gsl_linalg_LU_invert ( (gsl_matrix*)lu._matrix, perm, (gsl_matrix*)_matrix);
gsl_permutation_free( perm );
for(size_t i=1; i<=num_row(); i++)
for(size_t i=1; i<=num_row(); i++)
for(size_t j=1; j<=num_row(); j++)
if(operator()(i,j) != operator()(i,j)) throw(GeneralClassicException("MMatrix::invert()", "Failed to invert matrix - singular?"));
}
......@@ -288,7 +288,7 @@ Tmplt MMatrix<Tmplt>::trace() const
for(size_t i=2; i<=num_row() && i<=num_col(); i++) t+= operator()(i,i);
return t;
}
template double MMatrix<double>::trace() const;
template double MMatrix<double>::trace() const;
template m_complex MMatrix<m_complex>::trace() const;
......@@ -325,19 +325,19 @@ template std::pair<MVector<m_complex>, MMatrix<m_complex> > MMatrix<double>::eig
///////////// OPERATORS
MMatrix<double>& operator *= (MMatrix<double>& m1, MMatrix<double> m2)
{
MMatrix<double> out(m1.num_row(), m2.num_col());
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1., (gsl_matrix*)m1._matrix, (gsl_matrix*)m2._matrix, 0., (gsl_matrix*)out._matrix);
{
MMatrix<double> out(m1.num_row(), m2.num_col());
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1., (gsl_matrix*)m1._matrix, (gsl_matrix*)m2._matrix, 0., (gsl_matrix*)out._matrix);
m1 = out;
return m1;
}
MMatrix<m_complex> & operator *= (MMatrix<m_complex>& m1, MMatrix<m_complex> m2)
{
MMatrix<m_complex> out(m1.num_row(), m2.num_col());
gsl_blas_zgemm( CblasNoTrans, CblasNoTrans, m_complex_build(1.), (gsl_matrix_complex*)m1._matrix,
(gsl_matrix_complex*)m2._matrix, m_complex_build(0.), (gsl_matrix_complex*)out._matrix);
m1 = out;
{
MMatrix<m_complex> out(m1.num_row(), m2.num_col());
gsl_blas_zgemm( CblasNoTrans, CblasNoTrans, m_complex_build(1.), (gsl_matrix_complex*)m1._matrix,
(gsl_matrix_complex*)m2._matrix, m_complex_build(0.), (gsl_matrix_complex*)out._matrix);
m1 = out;
return m1;
}
......@@ -371,13 +371,13 @@ template std::istream& operator>>(std::istream& in, MMatrix<m_complex>& mat);
///////////////// INTERFACES
const gsl_matrix* MMatrix_to_gsl(const MMatrix<double>& m)
const gsl_matrix* MMatrix_to_gsl(const MMatrix<double>& m)
{
if(m._matrix == NULL) throw(GeneralClassicException("MMatrix_to_gsl", "Attempt to reference uninitialised matrix"));
return (gsl_matrix*)m._matrix;
}
const gsl_matrix_complex* MMatrix_to_gsl(const MMatrix<m_complex>& m)
const gsl_matrix_complex* MMatrix_to_gsl(const MMatrix<m_complex>& m)
{
if(m._matrix == NULL) throw(GeneralClassicException("MMatrix_to_gsl", "Attempt to reference uninitialised matrix"));
return (gsl_matrix_complex*)m._matrix;
......@@ -412,7 +412,7 @@ MMatrix<m_complex> complex(MMatrix<double> real)
MMatrix<m_complex> complex(MMatrix<double> real, MMatrix<double> imaginary)
{
if(real.num_row() != imaginary.num_row() || real.num_col() != imaginary.num_col())
if(real.num_row() != imaginary.num_row() || real.num_col() != imaginary.num_col())
throw(GeneralClassicException("MMatrix<m_complex>::complex", "Attempt to build complex matrix when real and imaginary matrix don't have the same size"));
MMatrix<m_complex> mc(real.num_row(), real.num_col());
for(size_t i=1; i<=mc.num_row(); i++)
......@@ -431,4 +431,4 @@ MVector<Tmplt> MMatrix<Tmplt>::get_mvector(size_t column) const
template MVector<double> MMatrix<double>::get_mvector(size_t column) const;
template MVector<m_complex> MMatrix<m_complex>::get_mvector(size_t column) const;
}
}
\ No newline at end of file
/*
/*
* Copyright (c) 2015, Chris Rogers
* All rights reserved.
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
* 1. Redistributions of source code must retain the above copyright notice,
* this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
* 3. Neither the name of STFC nor the names of its contributors may be used to
* endorse or promote products derived from this software without specific
* 3. Neither the name of STFC nor the names of its contributors may be used to
* endorse or promote products derived from this software without specific
* prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
* POSSIBILITY OF SUCH DAMAGE.
*/
......@@ -31,16 +31,16 @@
namespace interpolation {
///////////////////////// m_complex ////////////////////////////
std::ostream& operator<<(std::ostream& out, m_complex c)
{
out << re(c) << " r " << im(c) << " i";
return out;
std::ostream& operator<<(std::ostream& out, m_complex c)
{
out << re(c) << " r " << im(c) << " i";
return out;
}
std::istream& operator>>(std::istream& in, m_complex& c)
{
std::string dummy;
in >> re(c) >> dummy >> im(c) >> dummy;
std::istream& operator>>(std::istream& in, m_complex& c)
{
std::string dummy;
in >> re(c) >> dummy >> im(c) >> dummy;
return in;
}
......@@ -50,14 +50,14 @@ std::istream& operator>>(std::istream& in, m_complex& c)
template <typename Tmplt>
MVector<Tmplt>::MVector( size_t i ) : _vector(NULL)
{
build_vector(i);
{
build_vector(i);
}
template MVector<double> ::MVector(size_t i);
template MVector<m_complex>::MVector(size_t i);
template <typename Tmplt>
MVector<Tmplt>::MVector( const MVector<Tmplt>& mv) : _vector(NULL)
MVector<Tmplt>::MVector( const MVector<Tmplt>& mv) : _vector(NULL)
{ *this = mv; }
template MVector<double> ::MVector(const MVector<double>&);
template MVector<m_complex>::MVector(const MVector<m_complex>&);
......@@ -130,10 +130,10 @@ template MMatrix<m_complex> MVector<m_complex>::T() const;
template <class Tmplt> bool operator==(const MVector<Tmplt>& c1, const MVector<Tmplt>& c2)
{
if(c1.num_row() != c2.num_row())
return false;
for(size_t i=0; i<c1.num_row(); i++)
if(c1(i+1) != c2(i+1)) return false;
if(c1.num_row() != c2.num_row())
return false;
for(size_t i=0; i<c1.num_row(); i++)
if(c1(i+1) != c2(i+1)) return false;
return true;
}
template bool operator==(const MVector<double>& c1, const MVector<double>& c2);
......@@ -141,23 +141,23 @@ template bool operator==(const MVector<m_complex>& c1, const MVector<m_complex>&
template <>
MVector<double>& MVector<double>::operator= (const MVector<double>& mv)
{
{
if (&mv == this) return *this;
delete_vector();
if(!mv._vector) { _vector = NULL; return *this; }
_vector = gsl_vector_alloc( mv.num_row() );
gsl_vector_memcpy((gsl_vector*)_vector, (const gsl_vector*)mv._vector);
_vector = gsl_vector_alloc( mv.num_row() );
gsl_vector_memcpy((gsl_vector*)_vector, (const gsl_vector*)mv._vector);
return *this;
}
template <>
MVector<m_complex>& MVector<m_complex>::operator= (const MVector<m_complex>& mv)
{
{
if (&mv == this) return *this;
delete_vector();
if(!mv._vector) { _vector = NULL; return *this; }
_vector = gsl_vector_complex_alloc( mv.num_row() );
gsl_vector_complex_memcpy((gsl_vector_complex*)_vector, (const gsl_vector_complex*)mv._vector);
_vector = gsl_vector_complex_alloc( mv.num_row() );
gsl_vector_complex_memcpy((gsl_vector_complex*)_vector, (const gsl_vector_complex*)mv._vector);
return *this;
}
......@@ -165,7 +165,7 @@ template <class Tmplt> std::ostream& operator<<(std::ostream& out, MVector<Tmplt
{
out << v.num_row() << "\n";
for(size_t i=0; i<v.num_row(); i++) out << " " << v(i+1) << "\n";
return out;
return out;
}
template std::ostream& operator<<(std::ostream& out, MVector<double> v);
template std::ostream& operator<<(std::ostream& out, MVector<m_complex> v);
......@@ -176,14 +176,14 @@ template <class Tmplt> std::istream& operator>>(std::istream& in, MVector<Tmplt>
in >> n;
v = MVector<Tmplt>(n);
for(size_t i=1; i<=v.num_row(); i++) in >> v(i);
return in;
return in;
}
template std::istream& operator>>(std::istream& out, MVector<double>& v);
template std::istream& operator>>(std::istream& out, MVector<m_complex>& v);
const gsl_vector* MVector_to_gsl(const MVector<double>& vd)
const gsl_vector* MVector_to_gsl(const MVector<double>& vd)
{return vd.get_vector(vd);}
const gsl_vector_complex* MVector_to_gsl(const MVector<gsl_complex>& vc)
const gsl_vector_complex* MVector_to_gsl(const MVector<gsl_complex>& vc)
{return vc.get_vector(vc);}
template <class Tmplt>
......@@ -224,5 +224,4 @@ MVector<double> im (MVector<m_complex> c)
return d;
}
}
}
\ No newline at end of file
/*
/*
* Copyright (c) 2015, Chris Rogers
* All rights reserved.
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
* 1. Redistributions of source code must retain the above copyright notice,
* this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
* 3. Neither the name of STFC nor the names of its contributors may be used to
* endorse or promote products derived from this software without specific
* 3. Neither the name of STFC nor the names of its contributors may be used to
* endorse or promote products derived from this software without specific
* prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
* POSSIBILITY OF SUCH DAMAGE.
*/
......@@ -43,11 +43,11 @@ bool SquarePolynomialVector::_printHeaders=true;
std::vector< std::vector< std::vector<int> > > SquarePolynomialVector::_polyKeyByPower;
std::vector< std::vector< std::vector<int> > > SquarePolynomialVector::_polyKeyByVector;
SquarePolynomialVector::SquarePolynomialVector()
SquarePolynomialVector::SquarePolynomialVector()
: _pointDim(0), _polyCoeffs() {
}
SquarePolynomialVector::SquarePolynomialVector (const SquarePolynomialVector& pv)
SquarePolynomialVector::SquarePolynomialVector (const SquarePolynomialVector& pv)
: _pointDim(pv._pointDim),
_polyCoeffs(pv._polyCoeffs.num_row(), pv._polyCoeffs.num_col(), 0.) {
SetCoefficients(pv._polyCoeffs);
......@@ -114,7 +114,7 @@ void SquarePolynomialVector::F(const double* point, double* value)
void SquarePolynomialVector::F(const MVector<double>& point, MVector<double>& value) const
{
MVector<double> polyVector(_polyCoeffs.num_col(), 1);
MVector<double> polyVector(_polyCoeffs.num_col(), 1);
MakePolyVector(point, polyVector);
MVector<double> answer = _polyCoeffs*polyVector;
for(unsigned int i=0; i<ValueDimension(); i++) value(i+1) = answer(i+1);
......@@ -227,7 +227,7 @@ void SquarePolynomialVector::PrintContainer(std::ostream& out, const Container&
std::stringstream strstr1("");
std::stringstream strstr2("");
class Container::const_iterator it1 = container.begin();
class Container::const_iterator it2 = it1;
class Container::const_iterator it2 = it1;
while(it1!=container.end())
{
it2++;
......@@ -245,6 +245,4 @@ void SquarePolynomialVector::PrintContainer(std::ostream& out, const Container&
if(pad_at_start) out << strstr2.str();
}
}
}
\ No newline at end of file
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