Here is a C# implementation of a general rank matrix. The data is stored in a 1D array, and the function int GetIndex(int,int,int...) finds the index in the 1D array that corresponds to the n-th rank tensor indexes. The reverse is called int[] GetIndexes(int)
The code also supports both column-major ordering (the default) and row-major ordering.
public enum IndexOrdering
{
ColumnMajor,
RowMajor,
}
public class Matrix<T>
{
readonly T[] _data;
readonly int[] _shape;
public Matrix(IndexOrdering ordering, int[] shape)
{
Ordering = ordering;
Rank = shape.Length;
_shape = shape;
Size = shape.Aggregate(1, (s, l) => s * l);
_data = new T[Size];
}
public Matrix(params int[] shape)
: this(IndexOrdering.ColumnMajor, shape) { }
public Matrix(IndexOrdering ordering, int[] shape, T[] data)
: this(ordering, shape)
{
Array.Copy(data, _data, Size);
}
public int Rank { get; }
public int Size { get; }
public IReadOnlyList<int> Shape { get => _shape; }
internal T[] Data { get => _data; }
public IndexOrdering Ordering { get; set; }
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public int GetIndex(params int[] indexes)
{
switch (Ordering)
{
case IndexOrdering.ColumnMajor:
{
int index = 0;
for (int i = 0; i < Rank; i++)
{
index = _shape[i] * index + indexes[i];
}
return index;
}
case IndexOrdering.RowMajor:
{
int index = 0;
for (int i = Rank - 1; i >= 0; i--)
{
index = _shape[i] * index + indexes[i];
}
return index;
}
default:
throw new NotSupportedException();
}
}
public int[] GetIndexes(int index)
{
int[] indexes = new int[Rank];
switch (Ordering)
{
case IndexOrdering.ColumnMajor:
{
for (int i = Rank - 1; i >= 0; i--)
{
// div = index/shape[i]
// indexes[i] = index - shape[i]*div
// index = div
index = Math.DivRem(index, _shape[i], out indexes[i]);
}
return indexes;
}
case IndexOrdering.RowMajor:
{
for (int i = 0; i < Rank; i++)
{
// div = index/shape[i]
// indexes[i] = index - shape[i]*div
// index = div
index = Math.DivRem(index, _shape[i], out indexes[i]);
}
return indexes;
}
default:
throw new NotSupportedException();
}
}
public T this[params int[] indexes]
{
get => _data[GetIndex(indexes)];
set => _data[GetIndex(indexes)] = value;
}
public override string ToString()
{
return $"[{string.Join(",", _data)}]";
}
}
To test the code with a rank=3 matrix (3D matrix), I used:
static void Main(string[] args)
{
const int L0 = 4;
const int L1 = 3;
const int L2 = 2;
var A = new Matrix<int>(L0, L1, L2);
Console.WriteLine($"rank(A)={A.Rank}");
Console.WriteLine($"size(A)={A.Size}");
Console.WriteLine($"shape(A)=[{string.Join(",", A.Shape)}]");
int index = 0;
for (int i = 0; i < L0; i++)
{
for (int j = 0; j < L1; j++)
{
for (int k = 0; k < L2; k++)
{
A[i, j, k] = index++;
}
}
}
Console.WriteLine($"A=");
for (int i = 0; i < L0; i++)
{
for (int j = 0; j < L1; j++)
{
Console.Write($"[{i},{j},..] = ");
for (int k = 0; k < L2; k++)
{
Console.Write($"{A[i, j, k]} ");
}
Console.WriteLine();
}
Console.WriteLine();
}
for (int idx = 0; idx < A.Size; idx++)
{
var ixs = A.GetIndexes(idx);
Console.WriteLine($"A[{string.Join(",", ixs)}] = {A.Data[idx]}");
}
}
with output
rank(A)=3
size(A)=24
shape(A)=[4,3,2]
A=
[0,0,..] = 0 1
[0,1,..] = 2 3
[0,2,..] = 4 5
[1,0,..] = 6 7
[1,1,..] = 8 9
[1,2,..] = 10 11
[2,0,..] = 12 13
[2,1,..] = 14 15
[2,2,..] = 16 17
[3,0,..] = 18 19
[3,1,..] = 20 21
[3,2,..] = 22 23
A[0,0,0] = 0
A[0,0,1] = 1
A[0,1,0] = 2
A[0,1,1] = 3
A[0,2,0] = 4
A[0,2,1] = 5
A[1,0,0] = 6
A[1,0,1] = 7
...
A[3,1,1] = 21
A[3,2,0] = 22
A[3,2,1] = 23
