Commit ee44770a authored by COULAUD Olivier's avatar COULAUD Olivier

stdComplex is now a class that extends std:complex class to avoid using the...

stdComplex is now a class that extends std:complex class to avoid using the function __muldc3 from libgcc
parent de09ab38
......@@ -277,12 +277,12 @@ public:
{
// Perform entrywise product manually
for (unsigned int j=0; j<opt_rc; ++j){
// FX[j].addMul(newFComplex<FReal>(scale*FC[idx*opt_rc + j].real(),
// scale*FC[idx*opt_rc + j].imag()),
// FY[j]);
FX[j] += stdComplex<FReal>(scale*FC[idx*opt_rc + j].real(),
scale*FC[idx*opt_rc + j].imag())
*FY[j];
FX[j].addMul(stdComplex<FReal>(scale*FC[idx*opt_rc + j].real(),
scale*FC[idx*opt_rc + j].imag()),
FY[j]);
// FX[j] += stdComplex<FReal>(scale*FC[idx*opt_rc + j].real(),
// scale*FC[idx*opt_rc + j].imag())
// *FY[j];
}
}
......@@ -463,8 +463,8 @@ public:
{
// Perform entrywise product manually
for (unsigned int j=0; j<opt_rc; ++j){
// FX[j].addMul(FC[TreeLevel][idx*opt_rc + j],FY[j]);
FX[j] += (FC[TreeLevel][idx*opt_rc + j]*FY[j]);
FX[j].addMul(FC[TreeLevel][idx*opt_rc + j],FY[j]);
// FX[j] += (FC[TreeLevel][idx*opt_rc + j]*FY[j]);
}
}
......
......@@ -17,10 +17,35 @@
#include <complex>
template<typename T>
using stdComplex = std::complex<T> ;
//using stdComplex = std::complex<T> ;
class stdComplex : public std::complex<T> {
public:
using Base = std::complex<T> ;
using Base::Base ;
/** Mul other and another and add the result to current complexe
*/
/*!
* \brief addMul perform z += other*another without using the function __muldc3 from libgcc. ;
*
* "Without -ffast-math or other, complex multiplication yields a call to the function __muldc3 from libgcc. "
* @see https://stackoverflow.com/questions/42659668/stdcomplex-multiplication-is-extremely-slow
*
* \param other
* \param another
*/
void addMul(const stdComplex<T>& other, const stdComplex<T>& another){
// this->complex[0] += (other.complex[0] * another.complex[0]) - (other.complex[1] * another.complex[1]);
// this->complex[1] += (other.complex[0] * another.complex[1]) + (other.complex[1] * another.complex[0]);
T realPart = this->real() + (other.real() * another.real()) - (other.imag() * another.imag());
T imagPart = this->imag() + (other.real() * another.imag()) + (other.imag() * another.real());
this->real(realPart) ;
this->imag(imagPart) ;
}
} ;
#endif //STDCOMPLEXE_HPP
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