Attached are two files (updates of the same ones from before) with some more
hacking on N-D arrays using expression templates.
The 2D matrices seem to be about 10% slower than static 2D matrices. I
haven't timed the higher dimensions. Supported operations include indexing
and slicing and setting the slice step. For example:
  NDArray!(int,3) A;
  ...
  // index an element
  A[2][3][4];

  // create a sub-array
  A[1..4][3][0..4]

  // create a sub-array with a stepsize of 2 in the first dimension
  A[1..4].by(2)[3][4]

Anyway, I'll probably continue poking around with this stuff and trying to get more of Norbert's ideas prototyped. There are probably lots of performance issues with what I wrote since it relies on the compiler to inline lots of templated code but at least the 2D case seems pretty good for now.

-Ben

Looks very intersting but:
if i try to compile the example without -release i get linker error.
And when i add -release i get compiler error:

ndarray.d(249): 'make(cast(NDArray *)(#this.data),#inds)' is not an lvalue

Why do i get this?


"Ben Hinkle" <bhinkle4@juno.com> wrote in message news:c7v52k$1g3d$1@digitaldaemon.com...
> Attached are two files (updates of the same ones from before) with some
more
> hacking on N-D arrays using expression templates.
> The 2D matrices seem to be about 10% slower than static 2D matrices. I
> haven't timed the higher dimensions. Supported operations include indexing
> and slicing and setting the slice step. For example:
>   NDArray!(int,3) A;
>   ...
>   // index an element
>   A[2][3][4];
>
>   // create a sub-array
>   A[1..4][3][0..4]
>
>   // create a sub-array with a stepsize of 2 in the first dimension
>   A[1..4].by(2)[3][4]
>
> Anyway, I'll probably continue poking around with this stuff and trying to get more of Norbert's ideas prototyped. There are probably lots of performance issues with what I wrote since it relies on the compiler to inline lots of templated code but at least the 2D case seems pretty good for now.
>
> -Ben

Good work. I just hope you don't expect to bring this to production quality. It is a great case study demonstrating some details, but I would not want to use that in production code. :-)

Especially the stepwise indexing/slicing: it is nice to see that this is possible, but if this kind of coding should make it into a standard library, we would be better of sticking with good old C++...

One remark, by the way: If opCall is supposed to create a C-aligned array (which I would think the reasonable default) the *last* stride should be 1. What your code is generating, is a Fortran-aligned array.

Ivan Senji wrote:

> Looks very intersting but:
> if i try to compile the example without -release i get linker error.

What is the link error? I don't get any errors with dmd. I just tried gdc and I had to make the private nd_elems public but otherwise everything looked fine.

> And when i add -release i get compiler error:
> 
> ndarray.d(249): 'make(cast(NDArray *)(#this.data),#inds)' is not an lvalue
> 
> Why do i get this?

I don't know. I'll check it out. I compiled the 2D stuff with -release just
fine before writing the n-d case.

-Ben

"Ben Hinkle" <bhinkle4@juno.com> wrote in message news:c7vqus$2h2o$1@digitaldaemon.com...
>
> Ivan Senji wrote:
>
> > Looks very intersting but:
> > if i try to compile the example without -release i get linker error.
>
> What is the link error? I don't get any errors with dmd. I just tried gdc and I had to make the private nd_elems public but otherwise everything looked fine.

With only this in main:
  NDArray!(double,3) x2 = NDArray!(double,3)(null,3,4,5);
and without -release i get:

E:\DLANGU~1\learn_d\TESTBE~1\testBensMDArray.obj(testBensMDArray) Error 42: Symbol Undefined _array_7ndarray

> > And when i add -release i get compiler error:
> >
> > ndarray.d(249): 'make(cast(NDArray *)(#this.data),#inds)' is not an
lvalue
> >
> > Why do i get this?
>
> I don't know. I'll check it out. I compiled the 2D stuff with -release
just
> fine before writing the n-d case.
>
> -Ben
>

On Thu, 13 May 2004 16:39:06 +0200, "Ivan Senji" <ivan.senji@public.srce.hr> wrote:

>"Ben Hinkle" <bhinkle4@juno.com> wrote in message news:c7vqus$2h2o$1@digitaldaemon.com...
>>
>> Ivan Senji wrote:
>>
>> > Looks very intersting but:
>> > if i try to compile the example without -release i get linker error.
>>
>> What is the link error? I don't get any errors with dmd. I just tried gdc and I had to make the private nd_elems public but otherwise everything looked fine.
>
>With only this in main:
>  NDArray!(double,3) x2 = NDArray!(double,3)(null,3,4,5);
>and without -release i get:
>
>E:\DLANGU~1\learn_d\TESTBE~1\testBensMDArray.obj(testBensMDArray) Error 42: Symbol Undefined _array_7ndarray
>
>> > And when i add -release i get compiler error:
>> >
>> > ndarray.d(249): 'make(cast(NDArray *)(#this.data),#inds)' is not an
>lvalue
>> >
>> > Why do i get this?
>>
>> I don't know. I'll check it out. I compiled the 2D stuff with -release
>just
>> fine before writing the n-d case.
>>
>> -Ben
>>
>

I re-organized some return statements and got rid of the opCall syntax
and that seemed to make the linking happier (at least on Windows).
I couldn't reproduce the "_array_7ndarray undefined" error - did you
include both ndarray.d and nd_example.d on the same call to dmd?
I've attached modified versions of ndarray.d and nd_example.d with the
tweaks for -inline and -release.

Norbert Nemec wrote:

> Good work. I just hope you don't expect to bring this to production quality. It is a great case study demonstrating some details, but I would not want to use that in production code. :-)

The 2D matrices are fast enough to be useful, I think. I haven't done any
profiling of the nd stuff so I have no idea how much could be inlined but it
all uses structs and direct calls so theoretically it should be fast. Of
course building in nd-arrays would most likely give the fastest code and
best error messages, etc.

> Especially the stepwise indexing/slicing: it is nice to see that this is possible, but if this kind of coding should make it into a standard library, we would be better of sticking with good old C++...

Hmm. and I went and added comments and everything :-)
I guess code is always simpler when you write it.

> One remark, by the way: If opCall is supposed to create a C-aligned array (which I would think the reasonable default) the *last* stride should be 1. What your code is generating, is a Fortran-aligned array.

I want Fortran-aligned for 2 reasons: Fortran and MATLAB. For numerical stuff I prefer interoperating with those two over C compatibility. I suppose since I work for the makers of MATLAB I have a bias but I've been trying to keep that in check.

It works now and the slicing works great, i love it :)
But i am having some more problems but this looks like a compiler
bug to me:
I wrote a function to print the array:

void printArray(NDArray!(double,3) y)
{
 int k=0;
 for(int n2=0; n2 < y.length[2]; n2++)
 {
      for(int n1=0; n1 < y.length[1]; n1++)
      {
           for(int n0=0; n0 < y.length[0]; n0++,k++)
           {
                printf("%d",cast(int)y[n0][n1][n2]);
           }
           printf("\n");
      }
      printf("\n");
  }
}

But compiler says:
E:\D language\learn_d\testBensMDArray\ndarray.d(16): template NDArray(T,int
N) is used as a type

To fix it i have to write:
alias NDArray!(double,3) array3;

and change everywhere in the code "NDArray!(double,3)"  to "array3" and it
works then.

On Thu, 13 May 2004 12:36:00 -0400, Ben Hinkle wrote:

> 
> import ndarray;
> 
> 
> int main() {
> 
>   // create a 3D matrix using generic template
>   NDArray!(double,3) x2 = NDArray!(double,3).make(null,3,4,5);
> 
>   for (int n2=0; n2<5; n2++) {
>     for (int n1=0; n1<4; n1++) {
>       for (int n0=0; n0<3; n0++) {
> 	x2[n0][n1][n2] = 100*n2 + 10*n1 + n0;
>       }
>     }
>   }

I really don't like this notation. The operation is parallel in nature but the for loop make it serial. I do not know how easy it is for the compiler to recognize and parallelize code like this but I would prefer a notation which don't specify a order of execution.

Here I use the notation I suggested a while back in a post name Einstein summation

index(3) i0;
index(4) i1;
index(5) i2;

x2[i0][i1][i2]=i2*100+i1*10+i0;

With this notation the loops are implicit so the compiler can parallelize away;

You can think of i0,i1,i2 as named slices.
the notation might be
slice i0=0..2,i1=0..3,i2=0..4;
instead.


Knud