GSoC 2012 Proposal: Continued Work on a D Linear Algebra library (SciD - std.linalg)

[Michael Chen]

For the Point 4, I really like to have high order functions like
reduceRow and reduceCol. then the function argument is simply the
reduceRow!foo(0,mat), here the foo is not a function operating on the
whole column but simply a function of two elements (e.g.
reduceRow!("a+b")(0,mat)). Or even better we could have a general
reduce function with row and col as template parameter so that we can
do reduce!(foo,row)(0,mat). I dont know whether we can optimize such
reduce function for different matrix type, but such a function will be
extremely useful from a user perspective.

On Tue, Apr 3, 2012 at 7:20 PM, Cristi Cobzarenco
<cristi.cobzarenco at gmail.com> wrote:
>
>
> On 3 April 2012 02:37, Caligo <iteronvexor at gmail.com> wrote:
>>
>> I've read **Proposed Changes and Additions**, and I would like to
>> comment and ask a few questions if that's okay.  BTW, I've used Eigen
>> a lot and I see some similarities here, but a direct rewrite may not
>> be the best thing because D > C++.
>>
>> 2.  Change the matrix & vector types, adding fixed-sized matrix
>> support in the process.
>>
>> This is a step in the right direction I think, and by that I'm talking
>> about the decision to remove the difference between a Vector and a
>> Matrix.  Also, fixed-size matrices are also a must.  There is
>> compile-time optimization that you won't be able to do for
>> dynamic-size matrices.
>>
>>
>> 3. Add value arrays (or numeric arrays, we can come up with a good name).
>>
>> I really don't see the point for these.  We have the built-in arrays
>> and the one in Phobos (which will get even better soon).
>
>
> The point of these is to have light-weight element wise operation support.
> It's true that in theory the built-in arrays do this. However, this library
> is built on top BLAS/LAPACK, which means operations on large matrices will
> be faster than on D arrays. Also, as far as I know, D doesn't support
> allocating dynamic 2-D arrays (as in not arrays of arrays), not to mention
> 2-D slicing (keeping track of leading dimension).
> Also I'm not sure how a case like this will be compiled, it may or may not
> allocate a temporary:
>
> a[] = b[] * c[] + d[] * 2.0;
>
> The expression templates in SciD mean there will be no temporary allocation
> in this call.
>
>>
>>
>> 4. Add reductions, partial reductions, and broadcasting for matrices and
>> arrays.
>>
>> This one is similar to what we have in Eigen, but I don't understand
>> why the operations are member functions (even in Eigen).  I much
>> rather have something like this:
>>
>> rowwise!sum(mat);
>>
>> Also, that way the users can use their own custom functions with much
>> ease.
>
>
> There is a problem with this design. You want each matrix type (be it
> general, triangular, sparse or even an expression node)  do be able to
> define its own implementation of sum: calling the right BLAS function and
> making whatever specific optimisations they can. Since D doesn't have
> argument dependent look-up (ADL), users can't provide specialisations for
> their own types. The same arguments apply to rowwise() and columnwise()
> which will return proxies specific to the matrix type. You could do
> something like this, in principle:
>
> auto sum( T )( T mat ) {
>   return mat.sum();
> }
>
> And if we want that we can add it, but this will provide no addition in
> extensibility. By the way, you can use std.algorithm with matrices since
> they offer range functionality, but it will be much slower to use
> reduce!mySumFunction(mat) than mat.sum() which uses a BLAS backend.
>
>
>>
>>
>>
>> 6. Add support for interoperation with D built-in arrays (or pointers).
>>
>> So I take that Matrix is not a sub-type?  why? If we have something like
>> this:
>>
>> struct Matrix(Real, size_t row, size_t col) {
>>
>>  Real[row*col] data;
>>  alias data this;
>> }
>>
>> then we wouldn't need any kind of interoperation with built-in arrays,
>> would we?  I think this would save us a lot of headache.
>>
>> That's just me and I could be wrong.
>
>
> Inter-operation referred more to having a matrix object wrapping a pointer
> to an already available piece of memory - maybe allocated through a region
> allocator, maybe resulting from some other library. This means we need to
> take care of different strides and different storage orders which cannot be
> handled by built-in arrays. Right now, matrices wrap ref-counted
> copy-on-write array types (ArrayData in the current code) - we decided last
> year that we don't want to use the garbage collector, because of its current
> issues. Also I would prefer not using the same Matrix type for
> pointer-wrappers and normal matrices because the former must have reference
> semantics while the latter have value semantics. I think it would be
> confusing if some matrices would copy their data and some would share
> memory.
>
>>
>>
>> I've got to tell you though, I'm very excited about this project and
>> I'll be watching it closely.
>>
>> cheers.
>
> ---
> Cristi Cobzarenco
> BSc in Artificial Intelligence and Computer Science
> University of Edinburgh
> Profile: http://www.google.com/profiles/cristi.cobzarenco
>