Bug 12367

Hello,

I'm using the gcc:

*Version: 3.2, as well as 3.3
*System type: linux
*options given when GCC was configured/built:
Configured with:
../configure --prefix=/usr --libdir=/usr/lib --with-slibdir=/lib --mandir=/u
sr/share/man --infodir=/usr/share/info --enable-shared --enable-threads=posi
x --disable-checking --enable-long-long --enable-__cxa_atexit --enable-langu
ages=c,c++,ada,f77,objc,java --host=i586-mandrake-linux-gnu --with-system-zl
ib
Thread model: posix
gcc version 3.2 (Mandrake Linux 9.0 3.2-1mdk)

* compilerflags were
g++ -O2 
g++ -O2 -W -Wall -Winline -finline-limit=5000 -ftemplate-depth-200 -fstrict-alia
sing -malign-double

I've got no error or warnings compiling the attached file. The program's
output is:

K = Matrix<d, 2, 2> = [
  [6.25, -5001.25],
  [-3.75, 2506.25]
]

Where K has the wrong result. It should be:

K = Matrix<d, 2, 2> = [
  [6.25, -3.75],
  [-3.75, 6.25]
]

The problem is related to the temporary temp_lhs used by the prod function. As 
used here in the example, using RVO results in wrong results. Parts of the 
result are correct, other not. The algorithm/code used is correct (checked with 
a lot of regression tests). A simple product of matrizes (using prod(const 
Matrix<>&, const Matrix<>&) not shown here) got the expected result. 

The attached file uses (depends on define) dynamic allocated memory or static 
arrays. valgrind shows some "Invalid read of size 8" (see comments on the 
lines). The interesting behavior is by defined VMET_RVO_BUG_WO_STATIC_TEMP, 
where static temps are used for evaluating temporaries inside the prod() 
function - there are no  errors found by valgrind! These looks for me as wrong 
code produced by gcc. Unfortunally I'm not able to understand the asm code 
produced.

It doesn't seems to be a problem with the FPU - the same with long types too.

Regards
Olaf

----8<----
#include <iostream>
#include <iomanip>

#define TVMET_DYNAMIC_DATA
#define TVMET_RVO_BUG_WO_STATIC_TEMP

#ifndef _tvmet_restrict
#define _tvmet_restrict  __restrict__
#endif


namespace tvmet {

template<class T, std::size_t Rows, std::size_t Cols> class Matrix;

template<class T1, class T2>
struct PromoteTraits { };

template<>
struct PromoteTraits<double, double> {
    typedef double value_type;
};

class XprNull
{
  XprNull& operator=(const XprNull&);

public:
  explicit XprNull() { }
};

template< class T >
static inline
T operator+(const T& lhs, XprNull) { return lhs; } // valgrind: Invalid read of 
size 8


namespace meta {

template<std::size_t Rows, std::size_t Cols,
	 std::size_t RowStride=0, std::size_t ColStride=0>
class Matrix
{
private:
  enum {
    doRows = (RowStride < Rows - 1) ? 1 : 0,
    doCols = (ColStride < Cols - 1) ? 1 : 0
  };

public:
  template<class Mtrx, class E, class Fcnl>
  static inline
  void assign2(Mtrx& mat, const E& expr, const Fcnl& fn) {
    fn.applyOn(mat(RowStride, ColStride), expr(RowStride, ColStride));
    Matrix<Rows * doCols, Cols * doCols,
           RowStride * doCols, (ColStride+1) * doCols>::assign2(mat, expr, fn);
  }

  template<class Mtrx, class E, class Fcnl>
  static inline
  void assign(Mtrx& mat, const E& expr, const Fcnl& fn) {
    Matrix<Rows, Cols,
           RowStride, 0>::assign2(mat, expr, fn);
    Matrix<Rows * doRows, Cols * doRows,
          (RowStride+1) * doRows, 0>::assign(mat, expr, fn);
  }
};

template<>
class Matrix<0, 0, 0, 0>
{
public:
  template<class Mtrx, class E, class Fcnl>
  static inline void assign2(Mtrx&, const E&, const Fcnl&) { }

  template<class Mtrx, class E, class Fcnl>
  static inline void assign(Mtrx&, const E&, const Fcnl&) { }
};

template<std::size_t Rows1, std::size_t Cols1,
	 std::size_t Cols2,
	 std::size_t RowStride1, std::size_t ColStride1,
	 std::size_t RowStride2, std::size_t ColStride2,
	 std::size_t K>
class gemm
{
private:
  enum {
    doIt = (K != Cols1 - 1) 
  };

public:
  template<class T1, class T2>
  static inline
  typename PromoteTraits<T1, T2>::value_type
  prod(const T1* _tvmet_restrict lhs, const T2* _tvmet_restrict rhs, 
std::size_t i, std::size_t j) { // valgrind: Invalid read of size 8
    return lhs[i * RowStride1 + K * ColStride1] * rhs[K * RowStride2 + j * 
ColStride2]
      + gemm<Rows1 * doIt, Cols1 * doIt,
               Cols2 * doIt, RowStride1 * doIt, ColStride1 * doIt,
               RowStride2 * doIt, ColStride2 * doIt, (K+1) * doIt>::prod(lhs, 
rhs, i, j);
  }
};

template<>
class gemm<0,0,0,0,0,0,0,0>
{
public:
  static inline
  XprNull prod(const void*, const void*, std::size_t, std::size_t) {
    return XprNull();
  }
};


} // namespace meta

template <class T1, class T2>
struct fcnl_Assign {
  typedef void						return_type;
  static inline return_type applyOn(T1& _tvmet_restrict lhs, T2 rhs) {
    lhs = static_cast<T1>(rhs);
  }
};

template<class T1, std::size_t Rows1, std::size_t Cols1,
	 class T2, std::size_t Cols2,
	 std::size_t RowStride1, std::size_t ColStride1,
	 std::size_t RowStride2, std::size_t ColStride2>
class XprMMProduct
{
public:
  typedef typename PromoteTraits<T1, T2>::value_type	value_type;

public:
  explicit XprMMProduct(const T1* _tvmet_restrict lhs, const T2* 
_tvmet_restrict rhs)
    : m_lhs(lhs), m_rhs(rhs)
  { }

  value_type operator()(std::size_t i, std::size_t j) const {
    return meta::gemm<Rows1, Cols1,
                      Cols2,
                      RowStride1, ColStride1,
                      RowStride2, ColStride2, 0>::prod(m_lhs, m_rhs, i, j);
  }

private:
  const T1* _tvmet_restrict 				m_lhs;
  const T2* _tvmet_restrict 				m_rhs;
};

template<class E>
class XprMatrixTranspose
{
public:
  typedef E						expr_type;
  typedef typename expr_type::value_type		value_type;

public: 
  explicit XprMatrixTranspose(const expr_type& e)
    : m_expr(e)
  { }

  value_type operator()(std::size_t i, std::size_t j) const { return m_expr(j, 
i); }

private:
  const expr_type& _tvmet_restrict			m_expr;
};

template<class E, std::size_t Rows, std::size_t Cols>
class XprMatrix
{
public:
  typedef E						expr_type;
  typedef typename expr_type::value_type		value_type;

public:
  explicit XprMatrix(const expr_type& e)
    : m_expr(e)
  { }

  value_type operator()(std::size_t i, std::size_t j) const {
    return m_expr(i, j);
  }

private:
  const expr_type& _tvmet_restrict			m_expr;
};

template<class T,
	 std::size_t Rows, std::size_t Cols,
	 std::size_t RowStride=Cols, std::size_t ColStride=1>
class MatrixConstReference
{
public:
  typedef T						value_type;

public:
  explicit MatrixConstReference(const Matrix<T, Rows, Cols>& rhs)
    : m_data(rhs.data())
  { }

  value_type operator()(std::size_t i, std::size_t j) const {
    return m_data[i * RowStride + j * ColStride];
  }

private:
  const value_type* _tvmet_restrict 			m_data;
};

template<class T, std::size_t Rows, std::size_t Cols>
class Matrix
{
public:
  typedef T						value_type;
  typedef Matrix<T, Rows, Cols>				this_type;

public:
  explicit Matrix() {
#if defined(TVMET_DYNAMIC_DATA)
    m_data = new double [Rows*Cols];
#endif
  }

  template<class E>
  explicit Matrix(const XprMatrix<E, Rows, Cols>& expr) {
#if defined(TVMET_DYNAMIC_DATA)
    m_data = new double [Rows*Cols];
#endif
    this->assign(expr, fcnl_Assign<value_type, typename E::value_type>());
  }

  ~Matrix() {
#if defined(TVMET_DYNAMIC_DATA)
    delete [] m_data;
#endif
  }

public: 
  value_type* _tvmet_restrict data() { return m_data; }
  const value_type* _tvmet_restrict data() const { return m_data; }
  value_type& _tvmet_restrict operator()(std::size_t i, std::size_t j) {
    return m_data[i * Cols + j];
  }
  value_type operator()(std::size_t i, std::size_t j) const {
    return m_data[i * Cols + j];
  }

public:
  typedef MatrixConstReference<T, Rows, Cols>   	ConstReference;
  ConstReference const_ref() const { return ConstReference(*this); }

private: 
  template<class E, class Fcnl>
  void assign(const E& expr, const Fcnl& fn) {
    meta::Matrix<Rows, Cols>::assign(*this, expr, fn);
  }
  
public: 
  template <class E> this_type& operator=(const XprMatrix<E, Rows, Cols>& rhs) {
    this->assign(rhs, fcnl_Assign<value_type, typename E::value_type>());
    return *this;
  }

public:
  std::ostream& print_on(std::ostream& os) const;

private:
#if defined(TVMET_DYNAMIC_DATA)
  value_type* 						m_data;
#else
  value_type						m_data[Rows*Cols];
#endif
};

template<class T, std::size_t Rows, std::size_t Cols>
inline
std::ostream& Matrix<T, Rows, Cols>::print_on(std::ostream& os) const
{
  std::streamsize w = os.width();


    os << std::setw(0) << "Matrix<" << typeid(T).name() << ", "
       << Rows << ", " << Cols << "> = [\n";
    for(std::size_t i = 0; i < Rows; ++i) {
      os << " [";
      for(std::size_t j = 0; j < (Cols - 1); ++j) {
	os << std::setw(w) << this->operator()(i, j) << ", ";
      }
      os << std::setw(w) << this->operator()(i, Cols - 1)
	 << (i != (Rows-1) ? "],\n" : "]\n");
    }
    os << "]";

  return os;
}

template<class T, std::size_t Rows, std::size_t Cols>
inline
std::ostream& operator<<(std::ostream& os, const Matrix<T, Rows, Cols>& rhs) {
  return rhs.print_on(os);
}


template<class E1, std::size_t Rows1, std::size_t Cols1,
	 class T2, std::size_t Cols2>
inline
XprMatrix<
  XprMMProduct<
    typename E1::value_type, Rows1, Cols1,
    T2, Cols2,
    Cols1, 1,
    Cols2, 1
  >,
  Rows1, Cols2
>
prod(const XprMatrix<E1, Rows1, Cols1>& lhs, const Matrix<T2, Cols1, Cols2>& 
rhs) {
  typedef Matrix<typename E1::value_type, Rows1, Cols1>	temp_matrix_type;
  typedef XprMMProduct<
    typename E1::value_type, Rows1, Cols1,
    T2, Cols2,
    Cols1, 1,
    Cols2, 1
  >							expr_type;
#if defined(TVMET_RVO_BUG_WO_STATIC_TEMP)
  static temp_matrix_type 				temp_lhs(lhs);
  return XprMatrix<expr_type, Rows1, Cols2>(expr_type(temp_lhs.data(), rhs.data
()));
#else
  return XprMatrix<expr_type, Rows1, Cols2>(expr_type(temp_matrix_type(lhs).data
(), rhs.data()));
#endif
}

template<class T, std::size_t Rows, std::size_t Cols>
inline
XprMatrix<
  XprMatrixTranspose<
    MatrixConstReference<T, Rows, Cols>
  >,
  Cols, Rows
>
trans(const Matrix<T, Rows, Cols>& rhs) {
  typedef XprMatrixTranspose<
    MatrixConstReference<T, Rows, Cols>
  >							expr_type;
  return XprMatrix<expr_type, Cols, Rows>(expr_type(rhs.const_ref()));
}

} // namespace tvmet


/**
 * Test driver
 */
using namespace std;

int main()
{
  tvmet::Matrix<double,3,2>	B;
  tvmet::Matrix<double,3,3>	D;
  tvmet::Matrix<double,2,2>	K;

  B(0,0) = -0.05;	B(0,1) =  0;
  B(1,0) =  0;		B(1,1) =  0.05;
  B(2,0) =  0.05;	B(2,1) = -0.05;

  D(0,0) = 2000;	D(0,1) = 1000;		D(0,2) = 0;
  D(1,0) = 1000;	D(1,1) = 2000;		D(1,2) = 0;
  D(2,0) = 0;		D(2,1) = 0;		D(2,2) = 500;

  K = prod(prod(trans(B), D), B);
  cout << "K = " << K << endl;
}
---->8----