Re: [eigen] Incorrect result of multiplication with scalar (mingw gcc 4.5 x64)

[Eamon Nerbonne <eamon.nerbonne@xxxxxxxxx>]

On Thu, Aug 26, 2010 at 15:22, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> wrote:

Hi,

Since we're discussing a possible compiler bug here, can you please
make a compilable testcase and give us your complete g++ command line?

Benoit

2010/8/26 Eamon Nerbonne <eamon.nerbonne@xxxxxxxxx>:

> Hi all,
>
> I stumbled across a compiler-specific bug (Current Hg tip version 'fix bad
> "using typename"'):
>
>     static MatrixXd Covariance(Eigen::MatrixBase<MatrixXd> const & points,
> VectorXd const & mean) {
>         return (points.colwise() - mean) * (points.colwise() -
> mean).transpose() * (1.0/(points.cols()-1.0));
>     }
>
>     static MatrixXd CovarianceB(Eigen::MatrixBase<MatrixXd>const & points,
> VectorXd const & mean) {
>         VectorXd diff = VectorXd::Zero(points.rows());
>         MatrixXd cov = MatrixXd::Zero(points.rows(),points.rows());
>         for(int i=0;i<points.cols();++i) {
>             diff.noalias() = points.col(i) - mean;
>             cov.noalias() += diff * diff.transpose();
>         }
>         return cov * (1.0/(points.cols()-1.0));
>     }
>
> These two methods return approximately equal results in MSC but not in GCC;
> turns out the multiplication by the scalar "1.0/(points.cols()-1.0" is
> ignored by GCC in the first example.
>
> Rearranging the order of the multiplication (i.e. putting the scalar factor
> in front or in the middle) provides a workaround at no performance cost, as
> does placing the result into an intermediate variable with the assignment
> operator rather than returning it directly (with perf cost).  On the other
> hand, this is also broken:
>
>     static MatrixXd Covariance(Eigen::MatrixBase<MatrixXd> const & points,
> VectorXd const & mean) {
>         return  (1.0/(points.cols()-1.0)) *  ((points.colwise() - mean) *
> (points.colwise() - mean).transpose());
>     }
>
> I'm guessing somewhere some specialization of GeneralProduct is wrong, but
> as to why this is compiler specific or what exactly the problem is, I don't
> know.
>
> --eamon@xxxxxxxxxxxx - Tel#:+31-6-15142163
>

//A quick test-case
//Compile with g++ -DSTANDALONE -I %EIGEN3_DIR% -march=native CovarianceStandalone.cpp && a
//or appropriately adapted...
#include <bench/BenchTimer.h>
#include <Eigen/Core>
using namespace Eigen;

static MatrixXd Covariance(Eigen::MatrixBase<MatrixXd> const & points, VectorXd const & mean) {
    return (points.colwise() - mean) * (points.colwise() - mean).transpose() * (1.0/(points.cols()-1.0))  ;
}

static MatrixXd CovarianceB(Eigen::MatrixBase<MatrixXd>const & points, VectorXd const & mean) {
    VectorXd diff = VectorXd::Zero(points.rows());
    MatrixXd cov = MatrixXd::Zero(points.rows(),points.rows());
    for(int i=0;i<points.cols();++i) {
        diff.noalias() = points.col(i) - mean;
        cov.noalias() += diff * diff.transpose();
    }
    return cov * (1.0/(points..cols()-1.0));
}

static MatrixXd CovarianceC(Eigen::MatrixBase<MatrixXd> const & points, VectorXd const & mean) {
    return (points.colwise() - mean) * (1.0/(points.cols()-1..0)) * (points.colwise() - mean).transpose()   ;
}


bool testIt() {
    MatrixXd points = Eigen::MatrixXd::Random(7,1001);
    VectorXd mean = points.rowwise().sum() * (1.0/points.cols());
    return Covariance(points,mean).isApprox(CovarianceB(points,mean));
}

bool testIt2() {
    MatrixXd points = Eigen::MatrixXd::Random(7,1001);
    VectorXd mean = points.rowwise().sum() * (1.0/points.cols());
    return Covariance(points,mean).isApprox(CovarianceC(points,mean));
}

bool testIt3() {
    MatrixXd points = Eigen::MatrixXd::Random(7,1001);
    VectorXd mean = points.rowwise().sum() * (1.0/points.cols());
    return CovarianceB(points,mean).isApprox(CovarianceC(points,mean));
}



#ifndef STANDALONE

#define BOOST_TEST_INCLUDED
#include <boost/test/unit_test.hpp>

BOOST_AUTO_TEST_CASE( covariance_standalone_test ){ 
    BOOST_CHECK(testIt());
    BOOST_CHECK(testIt2());
    BOOST_CHECK(testIt3());
}

#else
#include <iostream>
int main(int argc, char *argv[] ) {
    bool isOk =  testIt();
    bool isOk2 =  testIt2();
    bool isOk3 =  testIt3();
    std::cout<<"testIt(): "<< isOk<<"\n";
    std::cout<<"testIt2():"<< isOk2<<"\n";
    std::cout<<"testIt3():"<< isOk3<<"\n";
    if(!isOk || !isOk2 || !isOk3) {
        std::cout << "Bad output!\n";
        return 1;
    }

    return 0;
}

#endif