Re: GSoC 2012 Proposal: Continued Work on a D Linear Algebra library (SciD - std.linalg)
Cristi Cobzarenco wrote: Unfortunately, my proposal was not picked up this year. I might try to work on these ideas this summer anyway so I would still be very much interested in ideas and feedback, but I will probably have much less time if I'll be working somewhere else. I'm less familiar with your SciD code base but I have used Eigen regularly. Maybe you can answer my questions right away: 1. How are expression evaluated? Eigen uses no BLAS backend. All code is generated by themselves. Do you plan for such an option? 2. What is your goal for SciD? Do you want to have an Eigen in D? Are there places where a D solution may improved over the Eigen? 3. I'm not so sure about the array() stuff. I never liked in Eigen. Being able to call std.algorithm on the raw memory may be sufficient for the time being. 4. How about a sparse storage backend for sparse matrices? I'm missing sparse matrices in Eigen even though the situation is improving but they're not fully integrated yet. I'd like to support your efforts for including your work in Phobos. How about you clone Phobos and gradually move your work into a work in progress branch? Really little steps. This allows me to follow you closely. To familiarize myself with the code and I would then fill in unittest, documentation and benchmarking code where missing as we go. Jens PS It's great to hear that you plan for including your work in Phobos.
Re: GSoC 2012 Proposal: Continued Work on a D Linear Algebra library (SciD - std.linalg)
On 24 April 2012 08:32, Jens Mueller jens.k.muel...@gmx.de wrote: Cristi Cobzarenco wrote: Unfortunately, my proposal was not picked up this year. I might try to work on these ideas this summer anyway so I would still be very much interested in ideas and feedback, but I will probably have much less time if I'll be working somewhere else. I'm less familiar with your SciD code base but I have used Eigen regularly. Maybe you can answer my questions right away: 1. How are expression evaluated? Eigen uses no BLAS backend. All code is generated by themselves. Do you plan for such an option? The expressions are built in a similar way to how they are in Eigen (any many other linear algebra libraries), using expression templates that build, essentially, abstract syntax trees which get evaluated on assignment (or with the eval() function). The back-end for evaluation can be specified using a version flag. Specifying version=nodeps results in code being generated by the library with no external dependencies - this is much slower, as it doesn't use SIMD operations or anything of the sort. In the current state of the library, this is done by essentially providing naive implementations of the BLAS functions. In the revamped version, I plan this to be done in slightly different way (more similar to the way Eigen works) which removes the need for temporaries in some cases where using BLAS/LAPACK makes temporary allocation inevitable. In the immediate future I do not plan to include SIMD operations for the version=nodeps version and the library will be the most efficient when using BLAS LAPACK. 2. What is your goal for SciD? Do you want to have an Eigen in D? Are there places where a D solution may improved over the Eigen? These are really two questions: 2.1. What is your goal for SciD? Do you want to have an Eigen in D? Short anwer: no. We already have a working library, it's mostly the interface that's getting a face-lift. Some of the design decisions we made early on (such as support for views) and compiler bugs (such as a clash between templated opAssign and postblit constructor) forced us to make some unfortunate interface decisions. In the long run it became obvious that views are pretty much useless and a lot of the compiler bugs that were hindering progress were fixed - so I thought I would improve the library's interface. The more I worked on it, the more obvious it became I was converging onto Eigen's interface, with bits from MATLAB. 2.2. Are there places where a D solution may improved over the Eigen? I'm not entirely sure. There's the obvious stuff like more natural slicing syntax which D allows us to do. Easily swappable backends is another thing - I'm writing everything to allow custom backends for the matrices (the backend will be a template parameter in the new version). This would in theory, allow GPU matrices to be easily written if anyone wanted (particularly if someone wrapped something like CUBLAS). I am actually very interested in hearing from other people if they think there's stuff which we could do better thanks to D. 3. I'm not so sure about the array() stuff. I never liked in Eigen. Being able to call std.algorithm on the raw memory may be sufficient for the time being. I am starting to lean against array myself. I wanted mostly to allow for arbitrary typed, higher-dimensional matrices with element-wise operations. But having worked a bit on the changes already, I realised that it's very easy to disallow certain matrix operations on types which do not support multiplication for example. Also I think I found an elegant way of supporting higher-dimensional matrices using the same type, but I have to do a bit more research into that. In terms of running std.algorihtm on raw memory, that is and will be possible. Matrices provide two properties .data and .cdata, which provide a pointer to the raw memory (.cdata provides a constant one which may point to memory shared by multiple matrices, while .data makes sure the memory is unique). 4. How about a sparse storage backend for sparse matrices? I'm missing sparse matrices in Eigen even though the situation is improving but they're not fully integrated yet. Definitely. Sparse storage is a long term goal for the library, but not in the first iteration. I'd like to support your efforts for including your work in Phobos. How about you clone Phobos and gradually move your work into a work in progress branch? Really little steps. This allows me to follow you closely. To familiarize myself with the code and I would then fill in unittest, documentation and benchmarking code where missing as we go. That is actually a better idea than what I had in mind. I will do that. I am really glad you're interested in the project, but unfortunately I won't be able to start work on the library until the beginning of June. Also the proposal for the project was meant as my full-time job this summer,
Re: GSoC 2012 Proposal: Continued Work on a D Linear Algebra library (SciD - std.linalg)
Cristi Cobzarenco wrote: On 24 April 2012 08:32, Jens Mueller jens.k.muel...@gmx.de wrote: Cristi Cobzarenco wrote: Unfortunately, my proposal was not picked up this year. I might try to work on these ideas this summer anyway so I would still be very much interested in ideas and feedback, but I will probably have much less time if I'll be working somewhere else. I'm less familiar with your SciD code base but I have used Eigen regularly. Maybe you can answer my questions right away: 1. How are expression evaluated? Eigen uses no BLAS backend. All code is generated by themselves. Do you plan for such an option? The expressions are built in a similar way to how they are in Eigen (any many other linear algebra libraries), using expression templates that build, essentially, abstract syntax trees which get evaluated on assignment (or with the eval() function). The back-end for evaluation can be specified using a version flag. Specifying version=nodeps results in code being generated by the library with no external dependencies - this is much slower, as it doesn't use SIMD operations or anything of the sort. In the current state of the library, this is done by essentially providing naive implementations of the BLAS functions. In the revamped version, I plan this to be done in slightly different way (more similar to the way Eigen works) which removes the need for temporaries in some cases where using BLAS/LAPACK makes temporary allocation inevitable. In the immediate future I do not plan to include SIMD operations for the version=nodeps version and the library will be the most efficient when using BLAS LAPACK. I see. That means you have already an expression evaluator. I still wonder why a[] = b[] + c[] shouldn't be handled completely by the compiler. I thought that this is sufficient for the compiler for vectorization. snip 2.2. Are there places where a D solution may improved over the Eigen? I'm not entirely sure. There's the obvious stuff like more natural slicing syntax which D allows us to do. Easily swappable backends is another thing - I'm writing everything to allow custom backends for the matrices (the backend will be a template parameter in the new version). This would in theory, allow GPU matrices to be easily written if anyone wanted (particularly if someone wrapped something like CUBLAS). Sounds interesting. I am actually very interested in hearing from other people if they think there's stuff which we could do better thanks to D. 3. I'm not so sure about the array() stuff. I never liked in Eigen. Being able to call std.algorithm on the raw memory may be sufficient for the time being. I am starting to lean against array myself. I wanted mostly to allow for arbitrary typed, higher-dimensional matrices with element-wise operations. But having worked a bit on the changes already, I realised that it's very easy to disallow certain matrix operations on types which do not support multiplication for example. Also I think I found an elegant way of supporting higher-dimensional matrices using the same type, but I have to do a bit more research into that. Nice. 4. How about a sparse storage backend for sparse matrices? I'm missing sparse matrices in Eigen even though the situation is improving but they're not fully integrated yet. Definitely. Sparse storage is a long term goal for the library, but not in the first iteration. Don't want it in the first iteration. Just want to make sure that current design allows a straightforward extension in that direction. I'd like to support your efforts for including your work in Phobos. How about you clone Phobos and gradually move your work into a work in progress branch? Really little steps. This allows me to follow you closely. To familiarize myself with the code and I would then fill in unittest, documentation and benchmarking code where missing as we go. That is actually a better idea than what I had in mind. I will do that. I am really glad you're interested in the project, but unfortunately I won't be able to start work on the library until the beginning of June. Also the proposal for the project was meant as my full-time job this summer, but since I'll probably be working somewhere else now you should take into account that I probably won't be able to stick to the schedule I gave in the proposal. Take your time. I need to familiarize myself with the code anyway. Jens
Re: GSoC 2012 Proposal: Continued Work on a D Linear Algebra library (SciD - std.linalg)
Unfortunately, my proposal was not picked up this year. I might try to work
on these ideas this summer anyway so I would still be very much interested
in ideas and feedback, but I will probably have much less time if I'll be
working somewhere else.
---
Cristi Cobzarenco
BSc in Artificial Intelligence and Computer Science
University of Edinburgh
Profile: http://www.google.com/profiles/cristi.cobzarenco
On 10 April 2012 12:31, Cristi Cobzarenco cristi.cobzare...@gmail.comwrote:
Timon is, of course, right. I got a bit confused when trying to simplify
in a hurry. What I meant was actually something like this:
ops.d:
import std.stdio;
int sum( T )( T mat ){
writeln(ops.sum);
// return reduce!a+b( 0, mat );
return 1;
}
int numelems( T )( T mat ) {
// return mat.rows * mat.columns;
return 1;
}
int mean( T )( T mat ) {
return sum( mat ) / numelems( mat );
}
diagonal.d:
import ops;
import std.stdio;
struct DiagonalMatrix {
}
int sum( T : DiagonalMatrix )( T mat ) {
writeln( diagonal.sum );
// return reduce!a+b( 0, mat.diagonal() );
return 1;
}
main.d:
import ops;
import diagonal;
void main() {
DiagonalMatrix mat;
sum( mat ); // this will not compile, but if removed mean
mean( mat ); // will use ops.sum, not diagonal.sum anyway
}
The first problem could, in theory, be fixed using a member flag,
something like enum specializedSum = true; and ops.sum could be redefined
as: int sum( T )( T mat ) if( !T.specializedSum ) {}. But in this case
mean() would still try to use ops.sum which will fail. All of this can be
easily fixed by providing member functions, and having the free functions
call the member ones when they exist.
---
Cristi Cobzarenco
BSc in Artificial Intelligence and Computer Science
University of Edinburgh
Profile: http://www.google.com/profiles/cristi.cobzarenco
On 10 April 2012 10:35, Timon Gehr timon.g...@gmx.ch wrote:
On 04/10/2012 03:24 AM, Cristi Cobzarenco wrote:
Thanks for the suggestions!
I don't think UFCS would help us. Our problem is that we can't do this:
triangular.d:
struct TriangularMatrix {
}
void sum( T )( T x ) if( is( T : TriangularMatrix ) ) {
}
diagonal.d:
struct DiagonalMatrix {
}
void sum( T )( T x ) if( is( T : DiagonalMatrix ) ) {
}
main.d:
import diagonal;
import triangular;
void bar() {
TriangularMatrix a;
Diagonal b;
sum( a ); // this does not compile because sum() is ambiguous
sum( b ); // nor does this
}
There are no ambiguities in that example, and if ambiguities occur, they
can always be fixed manually.
This, AFAIK, is deliberate to avoid name hijacking - ADL in C++ had its
share of criticism. I doubt we will ever get this behaviour in D and
that is perhaps a good thing. I may have misunderstood UFCS though - or
what you meant by making non-member function calls look nicer - please
correct me if that's the case.
Don't worry about long names, t() is already the way transposition is
defined SciD. Moreover, it's a property so you can actually do a.t * a
- convenience galore. I'm also considering of creating a submodule like
std.linalg.short which defines aliases with short names for types and
free functions - this will allow particularly numerics-heavy functions
to be written more compactly. I'm not entirely sure it would be a good
idea though as it may sacrifice readability where it's most needed.
---
Cristi Cobzarenco
BSc in Artificial Intelligence and Computer Science
University of Edinburgh
Profile:
http://www.google.com/**profiles/cristi.cobzarencohttp://www.google.com/profiles/cristi.cobzarenco
If you change 'Diagonal' to 'DiagonalMatrix', this compiles fine.
Re: GSoC 2012 Proposal: Continued Work on a D Linear Algebra library (SciD - std.linalg)
On 04/10/2012 03:24 AM, Cristi Cobzarenco wrote:
Thanks for the suggestions!
I don't think UFCS would help us. Our problem is that we can't do this:
triangular.d:
struct TriangularMatrix {
}
void sum( T )( T x ) if( is( T : TriangularMatrix ) ) {
}
diagonal.d:
struct DiagonalMatrix {
}
void sum( T )( T x ) if( is( T : DiagonalMatrix ) ) {
}
main.d:
import diagonal;
import triangular;
void bar() {
TriangularMatrix a;
Diagonal b;
sum( a ); // this does not compile because sum() is ambiguous
sum( b ); // nor does this
}
There are no ambiguities in that example, and if ambiguities occur, they
can always be fixed manually.
This, AFAIK, is deliberate to avoid name hijacking - ADL in C++ had its
share of criticism. I doubt we will ever get this behaviour in D and
that is perhaps a good thing. I may have misunderstood UFCS though - or
what you meant by making non-member function calls look nicer - please
correct me if that's the case.
Don't worry about long names, t() is already the way transposition is
defined SciD. Moreover, it's a property so you can actually do a.t * a
- convenience galore. I'm also considering of creating a submodule like
std.linalg.short which defines aliases with short names for types and
free functions - this will allow particularly numerics-heavy functions
to be written more compactly. I'm not entirely sure it would be a good
idea though as it may sacrifice readability where it's most needed.
---
Cristi Cobzarenco
BSc in Artificial Intelligence and Computer Science
University of Edinburgh
Profile: http://www.google.com/profiles/cristi.cobzarenco
If you change 'Diagonal' to 'DiagonalMatrix', this compiles fine.
Re: GSoC 2012 Proposal: Continued Work on a D Linear Algebra library (SciD - std.linalg)
Timon is, of course, right. I got a bit confused when trying to simplify in
a hurry. What I meant was actually something like this:
ops.d:
import std.stdio;
int sum( T )( T mat ){
writeln(ops.sum);
// return reduce!a+b( 0, mat );
return 1;
}
int numelems( T )( T mat ) {
// return mat.rows * mat.columns;
return 1;
}
int mean( T )( T mat ) {
return sum( mat ) / numelems( mat );
}
diagonal.d:
import ops;
import std.stdio;
struct DiagonalMatrix {
}
int sum( T : DiagonalMatrix )( T mat ) {
writeln( diagonal.sum );
// return reduce!a+b( 0, mat.diagonal() );
return 1;
}
main.d:
import ops;
import diagonal;
void main() {
DiagonalMatrix mat;
sum( mat ); // this will not compile, but if removed mean
mean( mat ); // will use ops.sum, not diagonal.sum anyway
}
The first problem could, in theory, be fixed using a member flag, something
like enum specializedSum = true; and ops.sum could be redefined as: int
sum( T )( T mat ) if( !T.specializedSum ) {}. But in this case mean() would
still try to use ops.sum which will fail. All of this can be easily fixed
by providing member functions, and having the free functions call the
member ones when they exist.
---
Cristi Cobzarenco
BSc in Artificial Intelligence and Computer Science
University of Edinburgh
Profile: http://www.google.com/profiles/cristi.cobzarenco
On 10 April 2012 10:35, Timon Gehr timon.g...@gmx.ch wrote:
On 04/10/2012 03:24 AM, Cristi Cobzarenco wrote:
Thanks for the suggestions!
I don't think UFCS would help us. Our problem is that we can't do this:
triangular.d:
struct TriangularMatrix {
}
void sum( T )( T x ) if( is( T : TriangularMatrix ) ) {
}
diagonal.d:
struct DiagonalMatrix {
}
void sum( T )( T x ) if( is( T : DiagonalMatrix ) ) {
}
main.d:
import diagonal;
import triangular;
void bar() {
TriangularMatrix a;
Diagonal b;
sum( a ); // this does not compile because sum() is ambiguous
sum( b ); // nor does this
}
There are no ambiguities in that example, and if ambiguities occur, they
can always be fixed manually.
This, AFAIK, is deliberate to avoid name hijacking - ADL in C++ had its
share of criticism. I doubt we will ever get this behaviour in D and
that is perhaps a good thing. I may have misunderstood UFCS though - or
what you meant by making non-member function calls look nicer - please
correct me if that's the case.
Don't worry about long names, t() is already the way transposition is
defined SciD. Moreover, it's a property so you can actually do a.t * a
- convenience galore. I'm also considering of creating a submodule like
std.linalg.short which defines aliases with short names for types and
free functions - this will allow particularly numerics-heavy functions
to be written more compactly. I'm not entirely sure it would be a good
idea though as it may sacrifice readability where it's most needed.
---
Cristi Cobzarenco
BSc in Artificial Intelligence and Computer Science
University of Edinburgh
Profile:
http://www.google.com/**profiles/cristi.cobzarencohttp://www.google.com/profiles/cristi.cobzarenco
If you change 'Diagonal' to 'DiagonalMatrix', this compiles fine.
Re: GSoC 2012 Proposal: Continued Work on a D Linear Algebra library (SciD - std.linalg)
Hi, Cristi,
From change log of D 2.059. I saw that uniform function call syntax
was implemented. I hope you can leverage this feature to make non
member function calls look nicer. Another suggestion, please use
shorter function name for example M.t() instead of M.transpose() so
that long expression is easy to read.
Best,
Mo
On Mon, Apr 9, 2012 at 6:52 AM, Cristi Cobzarenco
cristi.cobzare...@gmail.com wrote:
On 8 April 2012 18:59, Caligo iteronve...@gmail.com wrote:
On Tue, Apr 3, 2012 at 6:20 AM, Cristi Cobzarenco
cristi.cobzare...@gmail.com wrote:
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.
I can't agree with building it on top of LAPACK or any other BLAS
implementation, but perhaps I shouldn't complain because I'm not the
one who's going to be implementing it. I like the approach Eigen has
taken where it offers its own BLAS implementation and, iirc, other
BLAS libraries can be used as optional back-ends.
Yes, I agree with you. I already built naive implementations for BLAS
LAPACK functions as part of my last year project, using external libraries
is optional (building with version=nodeps ensures absolutely no dependencies
are needed). The argument still stands, though. If you _do_ use external
BLAS libraries then using value arrays will _still_ be faster.
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).
I fail to see why there is any need for 2D arrays. We need to make
sure multidimensional arrays (matrices) have data in very good
arrangements. This is called tiling and it requires 1D arrays: 2D
arrays are stored as 1D arrays together with an indexing mechanism.
That is precisely what I meant. We need wrappers around D arrays, because
they, by themselves, do not support 2-D indexing. By providing a wrapper we
allow the nice syntax of matrix[ a, b ]. This kind of wrapping is already in
SciD by CowMatrix. I meant we shouldn't use D built-in arrays directly, not
not at all. Also, as mentioned before, we can't use new for allocation,
because we want the library to be GC-independent.
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.
Why are expression templates used?
As H. S. Teoh rightfully pointed out, it is important not to allocate
temporaries in matrix operations. You want to evaluate A = B * 3.0 + D * 2.0
(where .* is element-wise multiplication) as (in BLAS terms):
copy( B, A );
scal( 3.0, A );
axpy( D, 2.0, A );
Or with naive expression template evaluation:
for( i = 0 ; i n ; ++ i ) {
A[ i ] = B[ i ] * 3.0 + D * 2.0;
}
The D compiler would instead evaluate this as:
// allocate temporaries
allocate( tmp1, A.length );
allocate( tmp2, A.length );
allocate( tmp3, A.length );
// compute tmp1 = B * 3.0
copy( B, tmp1 );
scal( 3.0, tmp1 );
// compute tmp2 = D * 2.0;
copy( D, tmp2 );
scal( 2.0, tmp2 );
// compute tmp3 = tmp1 + tmp2;
copy( tmp1, tmp3 );
axpy( tmp2, 1.0, tmp1 );
// copy tmp3 into A
copy( tmp3, A );
Plenty of useless computation. Note this is not a fault of the compiler, it
has no way of knowing which temporaries can be optimised away for user types
(i.e. our matrices).
Are you pretty much rewriting Eigen in D?
No. It is just the interface - and even that only to a certain extent - that
will mimic Eigen. The core the library would be very similar to what I
implemented for SciD last year, which, you will find, is very D-like and not
at all like Eigen.
Re: GSoC 2012 Proposal: Continued Work on a D Linear Algebra library (SciD - std.linalg)
Thanks for the suggestions!
I don't think UFCS would help us. Our problem is that we can't do this:
triangular.d:
struct TriangularMatrix {
}
void sum( T )( T x ) if( is( T : TriangularMatrix ) ) {
}
diagonal.d:
struct DiagonalMatrix {
}
void sum( T )( T x ) if( is( T : DiagonalMatrix ) ) {
}
main.d:
import diagonal;
import triangular;
void bar() {
TriangularMatrix a;
Diagonal b;
sum( a ); // this does not compile because sum() is ambiguous
sum( b ); // nor does this
}
This, AFAIK, is deliberate to avoid name hijacking - ADL in C++ had its
share of criticism. I doubt we will ever get this behaviour in D and that
is perhaps a good thing. I may have misunderstood UFCS though - or what you
meant by making non-member function calls look nicer - please correct me if
that's the case.
Don't worry about long names, t() is already the way transposition is
defined SciD. Moreover, it's a property so you can actually do a.t * a -
convenience galore. I'm also considering of creating a submodule like
std.linalg.short which defines aliases with short names for types and free
functions - this will allow particularly numerics-heavy functions to be
written more compactly. I'm not entirely sure it would be a good idea
though as it may sacrifice readability where it's most needed.
---
Cristi Cobzarenco
BSc in Artificial Intelligence and Computer Science
University of Edinburgh
Profile: http://www.google.com/profiles/cristi.cobzarenco
On 10 April 2012 01:36, Michael Chen sth4...@gmail.com wrote:
Hi, Cristi,
From change log of D 2.059. I saw that uniform function call syntax
was implemented. I hope you can leverage this feature to make non
member function calls look nicer. Another suggestion, please use
shorter function name for example M.t() instead of M.transpose() so
that long expression is easy to read.
Best,
Mo
On Mon, Apr 9, 2012 at 6:52 AM, Cristi Cobzarenco
cristi.cobzare...@gmail.com wrote:
On 8 April 2012 18:59, Caligo iteronve...@gmail.com wrote:
On Tue, Apr 3, 2012 at 6:20 AM, Cristi Cobzarenco
cristi.cobzare...@gmail.com wrote:
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.
I can't agree with building it on top of LAPACK or any other BLAS
implementation, but perhaps I shouldn't complain because I'm not the
one who's going to be implementing it. I like the approach Eigen has
taken where it offers its own BLAS implementation and, iirc, other
BLAS libraries can be used as optional back-ends.
Yes, I agree with you. I already built naive implementations for BLAS
LAPACK functions as part of my last year project, using external
libraries
is optional (building with version=nodeps ensures absolutely no
dependencies
are needed). The argument still stands, though. If you _do_ use external
BLAS libraries then using value arrays will _still_ be faster.
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).
I fail to see why there is any need for 2D arrays. We need to make
sure multidimensional arrays (matrices) have data in very good
arrangements. This is called tiling and it requires 1D arrays: 2D
arrays are stored as 1D arrays together with an indexing mechanism.
That is precisely what I meant. We need wrappers around D arrays, because
they, by themselves, do not support 2-D indexing. By providing a wrapper
we
allow the nice syntax of matrix[ a, b ]. This kind of wrapping is
already in
SciD by CowMatrix. I meant we shouldn't use D built-in arrays directly,
not
not at all. Also, as mentioned before, we can't use new for allocation,
because we want the library to be GC-independent.
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.
Why are expression templates used?
As H. S. Teoh rightfully pointed out, it is important not to allocate
temporaries in matrix operations. You want to evaluate A = B * 3.0 + D *
2.0
(where .* is element-wise multiplication) as (in BLAS terms):
copy( B, A );
scal( 3.0, A );
axpy( D, 2.0, A );
Or with naive expression template evaluation:
for( i = 0 ; i n ; ++ i ) {
A[ i ] = B[ i ] * 3.0 + D * 2.0;
}
The D compiler would instead evaluate this as:
// allocate temporaries
allocate( tmp1, A.length );
allocate( tmp2, A.length );
allocate( tmp3, A.length );
// compute tmp1 = B * 3.0
copy( B, tmp1 );
scal( 3.0, tmp1 );
// compute tmp2 = D *
Re: GSoC 2012 Proposal: Continued Work on a D Linear Algebra library (SciD - std.linalg)
On Tue, Apr 3, 2012 at 6:20 AM, Cristi Cobzarenco cristi.cobzare...@gmail.com wrote: 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. I can't agree with building it on top of LAPACK or any other BLAS implementation, but perhaps I shouldn't complain because I'm not the one who's going to be implementing it. I like the approach Eigen has taken where it offers its own BLAS implementation and, iirc, other BLAS libraries can be used as optional back-ends. 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). I fail to see why there is any need for 2D arrays. We need to make sure multidimensional arrays (matrices) have data in very good arrangements. This is called tiling and it requires 1D arrays: 2D arrays are stored as 1D arrays together with an indexing mechanism. 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. Why are expression templates used? Are you pretty much rewriting Eigen in D ? I don't see why expression templates would be needed with D.
Re: GSoC 2012 Proposal: Continued Work on a D Linear Algebra library (SciD - std.linalg)
On Sun, Apr 08, 2012 at 12:59:08PM -0500, Caligo wrote: On Tue, Apr 3, 2012 at 6:20 AM, Cristi Cobzarenco cristi.cobzare...@gmail.com wrote: [...] 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. Why are expression templates used? Are you pretty much rewriting Eigen in D ? I don't see why expression templates would be needed with D. Expression templates are necessary to eliminate temporaries and optimize loops. The reason libraries like BLAS are so efficient is because they directly implement common operations like A = k*B + C. Without expression templates, a D implementation would create a temporary for k*B and then k*B + C and then assign the result, whereas expression templates allow you to directly compute the elements of A in a single loop over the matrix elements. T -- Change is inevitable, except from a vending machine.
Re: GSoC 2012 Proposal: Continued Work on a D Linear Algebra library (SciD - std.linalg)
On 8 April 2012 18:59, Caligo iteronve...@gmail.com wrote:
On Tue, Apr 3, 2012 at 6:20 AM, Cristi Cobzarenco
cristi.cobzare...@gmail.com wrote:
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.
I can't agree with building it on top of LAPACK or any other BLAS
implementation, but perhaps I shouldn't complain because I'm not the
one who's going to be implementing it. I like the approach Eigen has
taken where it offers its own BLAS implementation and, iirc, other
BLAS libraries can be used as optional back-ends.
Yes, I agree with you. I already built naive implementations for BLAS
LAPACK functions as part of my last year project, using external libraries
is optional (building with version=nodeps ensures absolutely no
dependencies are needed). The argument still stands, though. If you _do_
use external BLAS libraries then using value arrays will _still_ be faster.
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).
I fail to see why there is any need for 2D arrays. We need to make
sure multidimensional arrays (matrices) have data in very good
arrangements. This is called tiling and it requires 1D arrays: 2D
arrays are stored as 1D arrays together with an indexing mechanism.
That is precisely what I meant. We need wrappers around D arrays, because
they, by themselves, do not support 2-D indexing. By providing a wrapper we
allow the nice syntax of matrix[ a, b ]. This kind of wrapping is already
in SciD by CowMatrix. I meant we shouldn't use D built-in arrays directly,
not not at all. Also, as mentioned before, we can't use new for allocation,
because we want the library to be GC-independent.
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.
Why are expression templates used?
As H. S. Teoh rightfully pointed out, it is important not to allocate
temporaries in matrix operations. You want to evaluate A = B * 3.0 + D *
2.0 (where .* is element-wise multiplication) as (in BLAS terms):
copy( B, A );
scal( 3.0, A );
axpy( D, 2.0, A );
Or with naive expression template evaluation:
for( i = 0 ; i n ; ++ i ) {
A[ i ] = B[ i ] * 3.0 + D * 2.0;
}
The D compiler would instead evaluate this as:
// allocate temporaries
allocate( tmp1, A.length );
allocate( tmp2, A.length );
allocate( tmp3, A.length );
// compute tmp1 = B * 3.0
copy( B, tmp1 );
scal( 3.0, tmp1 );
// compute tmp2 = D * 2.0;
copy( D, tmp2 );
scal( 2.0, tmp2 );
// compute tmp3 = tmp1 + tmp2;
copy( tmp1, tmp3 );
axpy( tmp2, 1.0, tmp1 );
// copy tmp3 into A
copy( tmp3, A );
Plenty of useless computation. Note this is not a fault of the compiler, it
has no way of knowing which temporaries can be optimised away for user
types (i.e. our matrices).
Are you pretty much rewriting Eigen in D?
No. It is just the interface - and even that only to a certain extent -
that will mimic Eigen. The core the library would be very similar to what I
implemented for SciD last year, which, you will find, is very D-like and
not at all like Eigen.
Re: GSoC 2012 Proposal: Continued Work on a D Linear Algebra library (SciD - std.linalg)
Thanks for the explanation, now I get it. In case you are interested,
there is excellent article about monad style c++ template meta
programming by Bartosz Milewski which might be helpful for compile
time optimization for evaluation order.
Really looking forward to official release of the SciD. D is really
suitable for scientific computation. It will be great to have an
efficient and easy-to-use linear algebra library.
On Thu, Apr 5, 2012 at 7:42 AM, Cristi Cobzarenco
cristi.cobzare...@gmail.com wrote:
Thanks for the feedback!
On 4 April 2012 10:21, Michael Chen sth4...@gmail.com wrote:
another btw, there is also another great c++ linear algebra library
besides Eigen: Amardillo which has very simple matlab like interface
and performs on a par with Eigen.
I'll look into it, thanks.
On Wed, Apr 4, 2012 at 5:14 PM, Michael Chen sth4...@gmail.com wrote:
Btw, I really don't like the matrix api to be member functions. It is
hard for user to extend the library in an unified way. It is also ugly
when you want to chain the function calls.
On Wed, Apr 4, 2012 at 5:09 PM, Michael Chen sth4...@gmail.com wrote:
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
Well, as I said before, there's nothing stopping us from providing the free
functions that call the member functionsl. However there is something
stopping us from not providing a member function alternative. D's lack of
ADL means we would not allow easy extensibility of possible operations. Say
Alice invents a new kind of matrix type, the DiagonalMatrix type which
stores its elements in a 1D array (this isn't exactly how it would work with
the design we have, she would in fact have to define a storage type, but
bear with me). If she wants to implement sum on her matrix, she can't simply
add a specialisation to the sum( T ) function because of the lack of ADL. If
instead we implemented sum(T) as:
auto sum( T )( T matrix ) {
static if( is( typeof( T.init.sum() ) ) )
return matrix.sum;
else
return reduce!a+b( 0, matrix.elements() );
}
then Alice could simply define a DiagonalMatrix.sum() method that wouldn't
need to go through all the zero elements, for example. The idea of rowReduce
and columnReduce still works - we can provide this for when a user wants to
use her own operation for reducing. We could also have a way of
automatically optimising rowReduce!a+b( mat ) by calling
mat.rowwise().sum() if the operation is available - but that wouldn't be
exactly high priority.
.
On Tue, Apr 3, 2012 at 7:20 PM, Cristi Cobzarenco
cristi.cobzare...@gmail.com wrote:
On 3 April 2012 02:37, Caligo iteronve...@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
Re: GSoC 2012 Proposal: Continued Work on a D Linear Algebra library (SciD - std.linalg)
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.cobzare...@gmail.com wrote:
On 3 April 2012 02:37, Caligo iteronve...@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
Re: GSoC 2012 Proposal: Continued Work on a D Linear Algebra library (SciD - std.linalg)
Thanks for the feedback!
On 4 April 2012 10:21, Michael Chen sth4...@gmail.com wrote:
another btw, there is also another great c++ linear algebra library
besides Eigen: Amardillo which has very simple matlab like interface
and performs on a par with Eigen.
I'll look into it, thanks.
On Wed, Apr 4, 2012 at 5:14 PM, Michael Chen sth4...@gmail.com wrote:
Btw, I really don't like the matrix api to be member functions. It is
hard for user to extend the library in an unified way. It is also ugly
when you want to chain the function calls.
On Wed, Apr 4, 2012 at 5:09 PM, Michael Chen sth4...@gmail.com wrote:
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
Well, as I said before, there's nothing stopping us from providing the free
functions that call the member functionsl. However there is something
stopping us from not providing a member function alternative. D's lack of
ADL means we would not allow easy extensibility of possible operations. Say
Alice invents a new kind of matrix type, the DiagonalMatrix type which
stores its elements in a 1D array (this isn't exactly how it would work
with the design we have, she would in fact have to define a storage type,
but bear with me). If she wants to implement sum on her matrix, she can't
simply add a specialisation to the sum( T ) function because of the lack of
ADL. If instead we implemented sum(T) as:
auto sum( T )( T matrix ) {
static if( is( typeof( T.init.sum() ) ) )
return matrix.sum;
else
return reduce!a+b( 0, matrix.elements() );
}
then Alice could simply define a DiagonalMatrix.sum() method that wouldn't
need to go through all the zero elements, for example. The idea of
rowReduce and columnReduce still works - we can provide this for when a
user wants to use her own operation for reducing. We could also have a way
of automatically optimising rowReduce!a+b( mat ) by calling
mat.rowwise().sum() if the operation is available - but that wouldn't be
exactly high priority.
.
On Tue, Apr 3, 2012 at 7:20 PM, Cristi Cobzarenco
cristi.cobzare...@gmail.com wrote:
On 3 April 2012 02:37, Caligo iteronve...@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
Re: GSoC 2012 Proposal: Continued Work on a D Linear Algebra library (SciD - std.linalg)
Sorry for the follow-up, I just wanted to mention that I also uploaded the proposal on Melange at http://www.google-melange.com/gsoc/proposal/review/google/gsoc2012/cristicbz/41002#, but I still recommend reading it on github. Thank you, Cristian. --- Cristi Cobzarenco BSc in Artificial Intelligence and Computer Science University of Edinburgh Profile: http://www.google.com/profiles/cristi.cobzarenco On 27 March 2012 11:44, Cristi Cobzarenco cristi.cobzare...@gmail.com wrote: Hello everyone, My name is Cristian Cobzarenco and last year, mentored by David Simcha and co-mentored by Fawzi Mohamed and Andrei Alexandrescu, I worked on a fork of SciD as part of Google's Summer of Code. While we've got a lot done last summer there's still plenty to do and my proposal this year is to get the library into a finished state where it can be reviewed for inclusion in Phobos as std.linalg. The easiest place to read my proposal is at the library's wiki: https://github.com/cristicbz/scid/wiki/GSoC-2012-Proposal, I will convert the formatting and add it to Melange soon. Feedback would be most welcome, I'm interested particularly in if you think the priorities I set are correct or if there is something you would (or wouldn't) rather see in the library. Let me know if anyone wants me to attach the content to an e-mail (it would be a pretty long e-mail), to allow in-line commenting. Yours, Cristian. --- Cristi Cobzarenco BSc in Artificial Intelligence and Computer Science University of Edinburgh Profile: http://www.google.com/profiles/cristi.cobzarenco
