blitzpp / blitz

Blitz++ Multi-Dimensional Array Library for C++
https://github.com/blitzpp/blitz/wiki
Other
405 stars 84 forks source link

Wrong assignment between Array, and TinyMatrix for a certain layout #67

Open slayoo opened 5 years ago

slayoo commented 5 years ago

Migrated from: https://sourceforge.net/p/blitz/bugs/66/

As reported by @aherrmann (?) or @AndreasHerrmann (?):

I believe that I found a bug in the assignment of a TinyMatrix into a sub-array of non-standard storage order.

The attached code demonstrates it in detail.

In short: Take a rank-3 Array and assign a TinyMatrix into a subarray of appropriate shape. This will work if the Array has C-, or Fortran-order. However, it will fail if the Array has a different storage-order.

However, it will work, if you replace the TinyMatrix by another Array.

The issue was observed in Blitz-0.10 with both GCC and Clang.

Best, Andreas

#include <iostream>
#include <blitz/blitz.h>
#include <blitz/array.h>
#include <blitz/tinymat2.h>

int main() {
    using std::cout;
    using namespace blitz;
    auto all = Range::all();

    { // The broken case:
        TinyMatrix<int, 3, 3> m;
        m = 1, 2, 3,
            4, 5, 6,
            7, 8, 9;
        cout << "m:\n" << m << "\n";
        // Output:
        //   m:
        //   (1,2,3; 4,5,6; 7,8,9)

        GeneralArrayStorage<3> storage;
        // storage.ordering() = thirdDim, secondDim, firstDim;  // works
        // storage.ordering() = firstDim, secondDim, thirdDim;  // works
        storage.ordering() = thirdDim, firstDim, secondDim;  // error!
        Array<int, 3> A(2, 3, 3, storage);
        A = 10;
        A(0, all, all) = m;
        cout << "A:\n" << A << "\n";
        // Expected output:
        //   A:
        //   (0,1) x (0,2) x (0,2)
        //   [ 1 2 3
        //     4 5 6
        //     7 8 9
        //     10 10 10
        //     10 10 10
        //     10 10 10 ]
        //
        // Actual output:
        //   A:
        //   (0,1) x (0,2) x (0,2)
        //   [ 1 2 3       <
        //     4 5 3       < wrong!
        //     7 8 3       <
        //     10 10 10
        //     10 10 10
        //     10 10 10 ]
    }

    { // For reference, the working case:
        Array<int, 2> m(3, 3);
        m = 1, 2, 3,
            4, 5, 6,
            7, 8, 9;
        cout << "m:\n" << m << "\n";
        //   Output
        //   (0,2) x (0,2)
        //   [ 1 2 3
        //     4 5 6
        //     7 8 9 ]
        //

        GeneralArrayStorage<3> storage;
        // storage.ordering() = thirdDim, secondDim, firstDim;
        // storage.ordering() = firstDim, secondDim, thirdDim;
        storage.ordering() = thirdDim, firstDim, secondDim;
        Array<int, 3> A(2, 3, 3, storage);
        A = 10;
        A(0, all, all) = m;
        cout << "A:\n" << A << "\n";
        // Correct output:
        //   A:
        //   (0,1) x (0,2) x (0,2)
        //   [ 1 2 3       <
        //     4 5 6       < correct.
        //     7 8 9       <
        //     10 10 10
        //     10 10 10
        //     10 10 10 ]
    }
}