blitz
blitz copied to clipboard
Wrong assignment between Array, and TinyMatrix for a certain layout
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 ]
}
}