Multidimensional Arrays
Oskar Linde
oskar.lindeREM at OVEgmail.com
Fri Sep 8 07:33:59 PDT 2006
I have been implementing a strided multidimensional array type and am
interested in comments.
It is basically designed after Norbert Nemec's suggestion:
http://homepages.uni-regensburg.de/~nen10015/documents/D-multidimarray.html
+/- some things.
The basic array type looks like this:
struct Array(T,int dims = 1) {
T* ptr;
size_t[dims] range;
ptrdiff_t[dims] stride;
const dimensions = dims;
alias T ElementType;
...
}
Usage:
Initialization
// Create a 3x3 array
auto A = Array!(int,2)(3,3);
//This could use a nicer syntax
A.assign(5,6,7,
3,1,9,
8,8,1);
writefln("A = %s",A.toString);
/*
A = (3x3)
5 6 7
3 1 9
8 8 1
*/
auto B = Array!(int,2)(3,3);
B.assign(0,0,1,
0,2,1,
7,8,9);
Element wise operators:
writefln("A*B = %s",(A*B).toString);
Scalar operators:
writefln("5*B = %s",(5*B).toString);
As in Norberts proposal:
C.transpose returns a transposed array view (all dimensions reversed)
C.transpose(dimA,dimB) swaps dimension dimA and dimB
C.diag returns a N-1 dimensional view of the diagonal of the matrix.
E.g. a 5x5 identity matrix:
// Zeros is just a convenience wrapper to generate a zero initialized array
Auto A = zeros!(int)(5,5);
A.diag[] = 1;
writefln("A = %s",A.toString);
A = (5x5)
1 0 0 0 0
0 1 0 0 0
0 0 1 0 0
0 0 0 1 0
0 0 0 0 1
C.reverse(dim) reverses the given dimension (inverts the stride)
Indexing and slicing:
auto D = Array!(int,3)(7,7,7); // 7x7x7 dimensional array
D[1,2,3] => int
D[1,2,all] => 7 dimensional slice
D[1,all,all] => 7x7 dimensional slice
D[range(1,end-1), all, range(1,end-1)] => 5x7x5 dimensional subarray
D[end-1, range(end-3,end-1), all] => 2x7 dimensional slice
End works in the same way as $ does for regular arrays.
range(a,b) is the replacement for a..b (D doesn't support mixed
indexing/slicing otherwise)
Slices/subarrays are of course views that share the same data:
auto M = zeros!(int)(4,4);
M[range(0,2), range(0,2)][] = 1;
M[range(2,end), range(2,end)][] = 1;
writefln("M = %s", M.toString);
M = (4x4)
1 1 0 0
1 1 0 0
0 0 1 1
0 0 1 1
.dup => Returns a compact copy of the array
In place map over a function:
M.doMap((int y){return 1<<y;});
Functional map:
writefln("cos(M) = %s",M.map(&cos).toString);
cos(M) = (4x4)
0.5403 0.5403 1.0000 1.0000
0.5403 0.5403 1.0000 1.0000
1.0000 1.0000 0.5403 0.5403
1.0000 1.0000 0.5403 0.5403
Foreach iteration:
auto N = zeros!(int)(5,5);
int i = 0;
foreach(inout x; N)
x = ++i;
Of course, custom types are supported. Some examples just for fun: :)
auto Q = ones!(Rational!(long))(5,5) * 3;
writefln("Q/N = %s",(Q/N).toString);
Q/N = (5x5)
3 3/2 1 3/4 3/5
1/2 3/7 3/8 1/3 3/10
3/11 1/4 3/13 3/14 1/5
3/16 3/17 1/6 3/19 3/20
1/7 3/22 3/23 1/8 3/25
auto U = zeros!(Uncertain!(double))(4,4);
double val = -8;
foreach(inout T;U) {
// +- 10% uncertainty + 0.1
T += Uncertain!(double)(val,abs(0.1*val)+0.1);
val += 1;
}
writefln("U = %s",U.toString);
U = (4x4)
(-8.0 ± 0.9) (-7.0 ± 0.8) (-6.0 ± 0.7) (-5.0 ± 0.6)
(-4.0 ± 0.5) (-3.0 ± 0.4) (-2.0 ± 0.3) (-1.0 ± 0.2)
(0. ± 0.1) (1.0 ± 0.2) (2.0 ± 0.3) (3.0 ± 0.4)
(4.0 ± 0.5) (5.0 ± 0.6) (6.0 ± 0.7) (7.0 ± 0.8)
writefln("U/U = %s",(U/U).toString);
U/U = (4x4)
(1.0 ± 0.2) (1.0 ± 0.2) (1.0 ± 0.2) (1.0 ± 0.2)
(1.0 ± 0.3) (1.0 ± 0.3) (1.0 ± 0.3) (1.1 ± 0.4)
(nan ± inf) (1.1 ± 0.4) (1.0 ± 0.3) (1.0 ± 0.3)
(1.0 ± 0.3) (1.0 ± 0.2) (1.0 ± 0.2) (1.0 ± 0.2)
auto V = zeros!(Voltage)(4,4);
val = 0;
foreach(inout v; V) {
v = val * volt;
val += 0.1;
}
auto R = zeros!(Resistance)(4,4);
val = 100;
foreach(inout r; R) {
r = val * ohm;
val *= 2;
}
writefln("V = %s",V.toString);
writefln("R = %s",R.toString);
writefln("V/R = %s",(V/R).toString);
V = (4x4)
0 V 100 mV 200 mV 300 mV
400 mV 500 mV 600 mV 700 mV
800 mV 900 mV 1e+03 mV 1.1 V
1.2 V 1.3 V 1.4 V 1.5 V
R = (4x4)
100 Ω 200 Ω 400 Ω 800 Ω
1.6 kΩ 3.2 kΩ 6.4 kΩ 12.8 kΩ
25.6 kΩ 51.2 kΩ 102 kΩ 205 kΩ
410 kΩ 819 kΩ 1.64 MΩ 3.28 MΩ
V/R = (4x4)
0 A 500 µA 500 µA 375 µA
250 µA 156 µA 93.7 µA 54.7 µA
31.2 µA 17.6 µA 9.77 µA 5.37 µA
2.93 µA 1.59 µA 854 nA 458 nA
Comments:
The only tricky part about all this was getting opIndex working with
compile time variadic arguments. I have not yet gotten DMD to inline the
generated indexing functions, but it should be doable as they are all
trivial. I have not been able to pursue this due to a DMD inlining bug
that in some cases generates "xx is not an lvalue" when compiling with
-inline.
The main issue that remains is op_r operators that have to work a bit
differently from non _r operators(*).
I see this all more as a proof of concept than anything else. If anyone
is interested in the sources, I will find a way to put them up somewhere.
Future extensions would be writing Matrix types, and possibly hook them
up with a BLAS library. Anyone know if there exists BLAS implementations
that supports 80-bit reals?
Other things I want to look at are adding stride support to the range()
type and investigate the possibilities of generating expression templates.
*) Currently, if you have an:
opAdd(T)(T x) {...}
you can't also have an:
opAdd_r(T)(T x) {...}
since this will create ambiguous overloads. I find myself wanting to
define opAdd_r for all cases where there are no regular non _r operator,
but there is (AFAIK) no way to do that.
I would go as far as suggesting a change to D: If both a matching opX
and a matching opX_r overload exists, the opX overload is chosen. This
will demote the opX_r overloads a bit and is consistent with the way
commutative operators work, where if for example no a.opAdd(b) exists,
b.opAdd(a) is called.
/Oskar
More information about the Digitalmars-d-announce
mailing list