mfem: Potential bug when manipulating Vectors?

Hi,

I noticed this behavior when trying to write a new Riemann solver. I’m using GCC 9.1.0 with MFEM 4.4 installed through spack. Basically, if you add two vectors and assign to a third, you get a compiler error (fair enough). But if you subtract two vectors, you get no compiler error but end up getting a wrong result.

Vector A(5), B(5), C(5); 
// Set some values in A and B 
for (int i = 0; i < 5; i++){
        A(i) = i;
        B(i) = i*i;
    }

// error: no match for ‘operator+’ (operand types are ‘mfem::Vector’ and ‘mfem::Vector’)
C = A + B; 

// No compiler error, but incorrect result 
C = A - B;
C.Print(); 
// Prints out -6 -6 -6 -6 -6

// No compiler error and correct result with expected behavior (i.e. no change in A)
C = A; 
C -= B; 
A.Print();
// 0 1 2 3 4
C.Print(); 
// 0 0 -2 -6 -12

About this issue

Original URL
State: closed
Created 2 years ago
Comments: 18 (17 by maintainers)

Most upvoted comments

As this is a very technical problem, let me quickly explain what the issue is. The vector class overloaded the indirection operator (i.e. operator*double()). Since the subtraction operator is not overloaded, the expression Vector - Vector is resolved by the compiler to the pointer arithmetic expression double* - double*. Note that expressions of type double* + double* are not valid pointer arithmetic expressions (as explained in §8.5.6 of the C++ standard). This pointer difference the triggers the overloaded assignment operator on Vector.

An MWE exhibiting this issue is given by:

#include <iostream>

struct A {
    double data;
    std::string me;

    inline double operator=(double val)
    {
        data = val;
        std::cout << "Assigning to " << me << "(" << this << ")" << std::endl;
        return data;
    }

    inline operator double *()
    {
        std::cout << "Dereferencing to " << me << "(" << this << ")" << std::endl;
        return &data;
    }
};

int main()
{
    A a1, a2, a3;
    a1.me = "1";
    a2.me = "2";
    a3.me = "3";

    a1.data = 1.0;
    a2.data = 2.0;

    a3 = a2 - a1;

    return a3.data;
}

Edit: Will was faster. 😃

termi-official on May 20, 2022

Hi @varchasgopalaswamy,

In general, I wouldn’t recommend using expressions like C = A - B; with Vectors in MFEM. That this even compiles is actually because A and B are implicitly cast to double *, and so A - B becomes a difference of pointers (integer type), which is then promoted to double, and C is given that constant value.

It’s probably better to work with explicit for loops.

pazner on May 20, 2022

This indicates that we should overload operator- with a static_assert so this fails at compile time, right?

We could, an implementation would look like:

template <typename...>
struct always_false { static constexpr bool value = false; };

template <typename... Ts>
Vector &operator-(const Vector &v)
{ static_assert(always_false<Ts...>::value, "operator- not supported for Vector."); }

I don’t know if there are other such “dangerous” operators.

pazner on May 20, 2022

2B would be great for serac. If mfem provides those operator overloads, we would adopt them immediately and probably remove the expression template implementation.

samuelpmishLLNL on Aug 31, 2022

This indicates that we should overload operator- with a static_assert so this fails at compile time, right?

jandrej on May 20, 2022

There are also helper functions like

 /// Set v = v1 + alpha * v2.
void add(const Vector &v1, double alpha, const Vector &v2, Vector &v);

tzanio on May 20, 2022