## digitalmars.D - std.linalg

- FreeSlave (60/60) Oct 11 2013 There is "Matrices and linear algebra" module in wish list. Let's
- ponce (22/51) Oct 11 2013 Good idea since there is so much implementations of "fixed-width
- FreeSlave (19/68) Oct 11 2013 Well, pointer cast would work for one-dimensional arrays, but not
- H. S. Teoh (191/246) Oct 11 2013 I think we need to differentiate between multidimensional arrays (as a
- FreeSlave (68/383) Oct 11 2013 The problem is that algorithms still should know matrix template
- H. S. Teoh (123/246) Oct 11 2013 I seem to have lost the link to Denis' code, but you can see the docs
- deadalnix (3/3) Oct 11 2013 First thing. I see std.linarg I have no clue whatsoever what it
- H. S. Teoh (9/12) Oct 11 2013 I picked it up immediately. It's an obvious abbreviation of "linear
- Justin Whear (3/17) Oct 11 2013 Don't be silly. It'd obviously have to be
- Andrei Alexandrescu (3/9) Oct 11 2013 Honestly I first thought it was "linear-time algorithms".
- H. S. Teoh (12/22) Oct 11 2013 [...]
- Andrej Mitrovic (4/6) Oct 11 2013 I thought it was some kind of algae species.
- Brian Rogoff (12/19) Oct 11 2013 Because, amongst people who know, numerical linear algebra
- John Colvin (3/6) Oct 12 2013 linalg is a pretty standard contraction of linear algebra. It's
- SomeDude (19/25) Oct 11 2013 Firstly, my opinion is, it shouldn't be in the std package.
- FreeSlave (10/38) Oct 11 2013 For these cases we may let users to choose low-level backend if
- John Colvin (11/71) Oct 12 2013 I'd like to echo some of the points made by CJS.

There is "Matrices and linear algebra" module in wish list. Let's discuss its design. D is complicated language so it's difficult to choose the right way here. We need to find compromise between efficiency and convenient interface. I'm going to make some suggestions how this module should look like. First of all, it should provide two templates for matrices. Let's call them StaticMatrix and DynamicMatrix. The first one has "templated" size and therefore may use static arrays and compile-time checks. It can be useful when the size is determined by our needs, for example, in graphics. DynamicMatrix has variable size, i.e. it should be created in heap. It can be useful in all other math areas. Both templates should support all floating point types and moreover user-defined (for example wrappers for GMP library and others). For efficiency in both cases matrices should use one-dimensional array for inner representation. But actually I'm not sure if matrices should support other container types besides standard D arrays. The good thing about one-dimensional arrays is that they can be easily exposed to foreign functions, for example, to C libraries and OpenGL. So we should take care about memory layout - at least row-major and column-major. I think it can be templated too. But another question arises - which "majority" should we use in interface? Interface should not depend on inner representation. All functions need unambiguity to avoid complication and repetition of design. Well, actually we can deal with different majority in interface - we can provide something like "asTransposed" adapter, that will be applied by functions if needed, but then we will force user to check majority of matrix interface, it's not very good approach. Sometimes user takes data from some other source and wants to avoid copying in Matrix construction, but she also wants to get matrix functionality. So we should provide "arrayAsMatrix" adapter, that can adopt one-dimensional and two-dimensional arrays making them feel like matrices. It definitely should not make copy of dynamic array, but I'm not sure about static. About operation overloading. It's quite clear about 'add' and 'subtruct' operations, but what's about product? Here I think all 'op'-functions should be 'element by element' operations. So we can use all other operations too without ambiguity. For actual matrix multiplication it can provide 'multiply' or 'product' function. It's similar to Maxima approach, besides Maxima uses dot notation for these needs. Transposition. I've already mentioned "asTransposed" adapter. It should be useful to make matrix feel like transposed without its copying. We also can implement 'transpose' and 'transposed' functions. The first one transposes matrix in place. It's actually not allowed for non-square StaticMatrix since we can't change the size of this type of matrices at runtime. The second one returns copy so it's applicable in all cases. Actually I'm not sure should these functions be member functions or not. Invertible matrix. It must not be allowed for square StaticMatrix. DynamicMatrix may make checks at runtime. Same as for transposition we can implement 'invert' to invert in place and 'inverted' to make copy. There is issue here - what should we do when determinant is 0? I believe the best approach here is to throw exception since if user needs invertible matrix it's definitely exception when it can't be calculated. Please, make your suggestions too.

Oct 11 2013

Good idea since there is so much implementations of "fixed-width vectors/matrices" since the beginning of times. There has been such efforts in the past to makes it more standardized. There will be On Friday, 11 October 2013 at 16:10:21 UTC, FreeSlave wrote:Both templates should support all floating point types and moreover user-defined (for example wrappers for GMP library and others).And integers for eg. video processing. Thankfully it's not that much harder.For efficiency in both cases matrices should use one-dimensional array for inner representation. But actually I'm not sure if matrices should support other container types besides standard D arrays. The good thing about one-dimensional arrays is that they can be easily exposed to foreign functions, for example, to C libraries and OpenGL. So we should take care about memory layout - at least row-major and column-major. I think it can be templated too.About 3D graphics: I don't think choosing one or the other is that big a deal. I tend to prefer row major order but eg. gl3n uses the reverse. It's up to you. If you templated against order, you will be able to implement you "asTransposed" idea.Sometimes user takes data from some other source and wants to avoid copying in Matrix construction, but she also wants to get matrix functionality. So we should provide "arrayAsMatrix" adapter, that can adopt one-dimensional and two-dimensional arrays making them feel like matrices. It definitely should not make copy of dynamic array, but I'm not sure about static.Possible with pointer cast? (strict aliasing rule aside)About operation overloading. It's quite clear about 'add' and 'subtruct' operations, but what's about product? Here I think all 'op'-functions should be 'element by element' operations.The product to me is a special case. Element-wise matrix multiplication is seldom used, it's a lot more common to want normal matrix multiplication. I and others use "*" for that purpose. https://github.com/Dav1dde/gl3n/blob/master/gl3n/linalg.d#L1734 https://github.com/p0nce/gfm/blob/master/gfm/math/matrix.d#L150 https://github.com/gecko0307/dlib/blob/master/dlib/math/matrix.d#L171So we can use all other operations too without ambiguity. For actual matrix multiplication it can provide 'multiply' or 'product' function. It's similar to Maxima approach, besides Maxima uses dot notation for these needs.But it's dissimilar to shader languages and existing D libraries.Invertible matrix. It must not be allowed for square StaticMatrix. DynamicMatrix may make checks at runtime. Same as for transposition we can implement 'invert' to invert in place and 'inverted' to make copy. There is issue here - what should we do when determinant is 0?Use enforce to crash. Dividing by zero is not recoverable.

Oct 11 2013

On Friday, 11 October 2013 at 17:08:26 UTC, ponce wrote:Good idea since there is so much implementations of "fixed-width vectors/matrices" since the beginning of times. There has been such efforts in the past to makes it more standardized. There will be On Friday, 11 October 2013 at 16:10:21 UTC, FreeSlave wrote:Well, pointer cast would work for one-dimensional arrays, but not for dynamic two-dimensional matrices. Anyway it's just implementation specific, not design.Both templates should support all floating point types and moreover user-defined (for example wrappers for GMP library and others).And integers for eg. video processing. Thankfully it's not that much harder.For efficiency in both cases matrices should use one-dimensional array for inner representation. But actually I'm not sure if matrices should support other container types besides standard D arrays. The good thing about one-dimensional arrays is that they can be easily exposed to foreign functions, for example, to C libraries and OpenGL. So we should take care about memory layout - at least row-major and column-major. I think it can be templated too.About 3D graphics: I don't think choosing one or the other is that big a deal. I tend to prefer row major order but eg. gl3n uses the reverse. It's up to you. If you templated against order, you will be able to implement you "asTransposed" idea.Sometimes user takes data from some other source and wants to avoid copying in Matrix construction, but she also wants to get matrix functionality. So we should provide "arrayAsMatrix" adapter, that can adopt one-dimensional and two-dimensional arrays making them feel like matrices. It definitely should not make copy of dynamic array, but I'm not sure about static.Possible with pointer cast? (strict aliasing rule aside)That's actual a common issue to use operator overloading for multiplication or not. Some people does not like it, others prefer this way. That's why I want to discuss it. "*" is more practical, but user always expects "/" as opposite. Matrix has no division operation (except of dividing by scalar), only multiplication by invertible matrix.About operation overloading. It's quite clear about 'add' and 'subtruct' operations, but what's about product? Here I think all 'op'-functions should be 'element by element' operations.The product to me is a special case. Element-wise matrix multiplication is seldom used, it's a lot more common to want normal matrix multiplication. I and others use "*" for that purpose. https://github.com/Dav1dde/gl3n/blob/master/gl3n/linalg.d#L1734 https://github.com/p0nce/gfm/blob/master/gfm/math/matrix.d#L150 https://github.com/gecko0307/dlib/blob/master/dlib/math/matrix.d#L171So we can use all other operations too without ambiguity. For actual matrix multiplication it can provide 'multiply' or 'product' function. It's similar to Maxima approach, besides Maxima uses dot notation for these needs.But it's dissimilar to shader languages and existing D libraries.Invertible matrix. It must not be allowed for square StaticMatrix.I meant non-square of course. Sorry for this mistake. I don't know how to edit my posts here. The another issue is getting of minors. It can use some operator overload magic and look like matrix(start, middle1, middle2, end)(start, middle1, middle2, end) where the first () stands for rows and the second () stands for columns. Gap is between middle1 and middle2. There is efficiency issue - should it make copy or work like some kind of two-dimensional range?

Oct 11 2013

On Fri, Oct 11, 2013 at 06:10:19PM +0200, FreeSlave wrote:There is "Matrices and linear algebra" module in wish list. Let's discuss its design. D is complicated language so it's difficult to choose the right way here. We need to find compromise between efficiency and convenient interface. I'm going to make some suggestions how this module should look like.I think we need to differentiate between multidimensional arrays (as a data storage type) and linear algebra (operations performed on 2D arrays). These two intersect, but they also have areas that are not compatible with each other (e.g. matrix product vs. element-by-element product). Ideally, we should support both in a clean way. As far as the former is concerned, Denis has implemented a multidimensional array library, and I've independently done the same, with a slightly different interface. I think one or two others have implemented similar libraries as well. It would be good if we standardize the API so that our code can become interoperable. As far as the latter is concerned, I've been meaning to implement a double-description convex hull algorithm, but have been too busy to actually work on it. This particular algorithm is interesting, because it stress-tests (1) performance of D algorithms, and (2) challenges the design of matrix APIs because while the input vertices (resp. hyperplanes) can be interpreted as a matrix, the algorithm itself also needs to permute rows, which means it is most efficient when given an array-of-pointers representation, contrary to the usual flattened representations (as proposed below). I think there's a place for both, which is why we need to distinguish between data representation and the algorithms that work on them.First of all, it should provide two templates for matrices. Let's call them StaticMatrix and DynamicMatrix. The first one has "templated" size and therefore may use static arrays and compile-time checks. It can be useful when the size is determined by our needs, for example, in graphics. DynamicMatrix has variable size, i.e. it should be created in heap. It can be useful in all other math areas.I like this idea. Ideally, we should have many possible representations, but all conforming to a single API understood by all algorithms, so that you only have to write algorithms once, and they will work with any data structure. That's one key advantage of D, and we should make good use of it. I do not want to see a repetition of the C++ situation where there are so many different matrix/multidimensional array libraries, and all of them use incompatible representations and you cannot freely pass data from one to algorithms in the other. Then when no library meets precisely what you need, you're forced to reinvent yet another matrix class, which is a waste of time.Both templates should support all floating point types and moreover user-defined (for example wrappers for GMP library and others).Definitely, yes.For efficiency in both cases matrices should use one-dimensional array for inner representation. But actually I'm not sure if matrices should support other container types besides standard D arrays. The good thing about one-dimensional arrays is that they can be easily exposed to foreign functions, for example, to C libraries and OpenGL. So we should take care about memory layout - at least row-major and column-major. I think it can be templated too.We should not tie algorithms to specific data representations (concrete types). One key advantage of D is that you can write algorithms generically, such that they can work with *any* type as long as it conforms to a standard API. One excellent example is the range API: *anything* that conforms to the range API can be used with std.algorithm, not just a specific representation. In fact, std.range provides a whole bunch of different ranges and range wrappers, and all of them automatically can be used with std.algorithm, because the code in std.algorithm uses only the range API and never (at least in theory :P) depends on concrete types. We should take advantage of this feature. It would be good, of course, to provide some standard, commonly-used representations, for example row-major (or column-major) matrix classes / structs, etc.. But the algorithms should not directly depend on these concrete types. An algorithm that works with a matrix stored as a 1D array should also work with a matrix stored as a nested array of arrays, as well as a sparse matrix representation that uses some other kind of storage mechanism. As long as a type conforms to some standard matrix API, it should Just Work(tm) with any std.linalg algorithm.But another question arises - which "majority" should we use in interface? Interface should not depend on inner representation. All functions need unambiguity to avoid complication and repetition of design. Well, actually we can deal with different majority in interface - we can provide something like "asTransposed" adapter, that will be applied by functions if needed, but then we will force user to check majority of matrix interface, it's not very good approach.Algorithms shouldn't even care what majority the data representation is in. It should only access data via the standardized matrix API (whatever it is we decide on). The input type should be templated so that *any* type that conforms to this API will work. Of course, for performance-sensitive code, the user should be aware of which representations are best-performing, and make sure to pass in the appropriate type of representations; but we should not prematurely optimize here. Any linear algebra algorithms should be able to work with *any* type that conforms to a standard matrix API.Sometimes user takes data from some other source and wants to avoid copying in Matrix construction, but she also wants to get matrix functionality. So we should provide "arrayAsMatrix" adapter, that can adopt one-dimensional and two-dimensional arrays making them feel like matrices. It definitely should not make copy of dynamic array, but I'm not sure about static.If a function expects a 1xN matrix, we should be able to pass in an array and it should Just Work. Manually using adapters should not be needed. Of course, standard concrete matrix types provided by the library should have ctors / factory methods for initializing a matrix object that uses some input array as initial data -- if we design this correctly, it should be a cheap operation (the matrix type itself should just be a thin wrapper over the array to provide methods that conform to the standard matrix API). Then if some function F requires a matrix object, we should be able to just create a Matrix instance with our input array as initial data, and pass it to F.About operation overloading. It's quite clear about 'add' and 'subtruct' operations, but what's about product? Here I think all 'op'-functions should be 'element by element' operations. So we can use all other operations too without ambiguity. For actual matrix multiplication it can provide 'multiply' or 'product' function. It's similar to Maxima approach, besides Maxima uses dot notation for these needs.Here is where we see the advantage of separating representation from algorithm. Technically, a matrix is not the same thing as a 2D array, because a matrix has a specific interpretation in linear algebra, whereas a 2D array is just a 2D container of some elements. My suggestion would be to write a Matrix struct that wraps around a 2D array, and provides / overrides the overloaded operators to have a linear algebra interpretation. So, a 2D array type should have per-element operations, but once wrapped in a Matrix struct, it will acquire special matrix algebra operations like matrix products, inversion, etc.. In the most general case, a 2D array should be a specific instance of a multidimensional array, and a Matrix struct should be able to use any underlying representation that conforms to a 2D array API. For example: // Example of a generic multidimensional array type struct Array(int dimension, ElemType) { ... Array opBinary(string op)(Array x) { // implement per-element operations here } } // A matrix wrapper around a 2D array type. struct Matrix(T) if (is2DArray!T) { T representation; Matrix opBinary(string op)(Matrix x) if (op == "*") { // implement matrix multiplication here } Matrix opBinary(string op)(Matrix x) if (op != "*") { // forward to representation.opBinary to default // to per-element operations } // Provide operations specific to matrices that don't // exist in general multidimensional arrays. Matrix invert() { ... } } Array!(2,float) myArray, myOtherArray; auto arrayProd = myArray * myOtherArray; // per-element multiplication auto A = Matrix(myArray); // wrap array in Matrix wrapper auto B = Matrix(myOtherArray); auto C = A * B; // matrix product The idea of the Matrix struct here is that the user should be free to choose any underlying matrix representation: a 1D array in row-major or column-major representation, or a nested array of arrays, or a sparse array with some other kind of representation. As long as they provide a standard way of accessing array elements, Matrix should be able to accept them, and provide matrix algebra semantics for them.Transposition. I've already mentioned "asTransposed" adapter. It should be useful to make matrix feel like transposed without its copying. We also can implement 'transpose' and 'transposed' functions. The first one transposes matrix in place. It's actually not allowed for non-square StaticMatrix since we can't change the size of this type of matrices at runtime. The second one returns copy so it's applicable in all cases. Actually I'm not sure should these functions be member functions or not.The most generic approach to transposition is simply a reordering of indices. This difference is important once you get to 3D arrays and beyond, because then there is no unique transpose, but any permutation of array indices should be permissible. Denis' multidimensional arrays have a method that does O(1) reordering of array indices: basically, you create a "view" of the original array that has its indices swapped around. So there is no data copying; it's just a different "view" into the same underlying data. This approach of using "views" rather than copying data allows for O(1) submatrix extraction: if you have a 50x50 matrix, then you can take arbitrary 10x10 submatrices of it without needing to copy any of the data, which would be very expensive. Avoiding unnecessary copying becomes very important when the dimension of the array increases; if you have a 3D or 5D array, copying subarrays become extremely expensive very quickly. A .dup method should be provided in the cases where you actually *want* to copy the data, of course. Basically, subarrays / transpositions / index reordering should be regarded as generalizations of D's array slices. No data should be copied until necessary.Invertible matrix. It must not be allowed for square StaticMatrix.You mean for non-square StaticMatrix?DynamicMatrix may make checks at runtime. Same as for transposition we can implement 'invert' to invert in place and 'inverted' to make copy. There is issue here - what should we do when determinant is 0? I believe the best approach here is to throw exception since if user needs invertible matrix it's definitely exception when it can't be calculated.[...] Here I'd just like to say that inversion is specific to linear algebra, and doesn't apply to multidimensional arrays in general. So it should be implemented in the Matrix wrapper, not tied to any specific data representation. As for trying to invert non-invertible matrices: it's just another form of "division by zero". I'd say throwing an exception is the right thing to do here. Of course, you might want to provide an alternative API that checks for invertibility when the caller is not sure whether it can be inverted: since this check pretty much already computes the inverse if it exists, it's probably more efficient to have a single method of doing both. Since it's probably cheaper to do this inversion / check in-place, you could just have a method that takes a Matrix by ref: struct Matrix(T) { ... // Returns true if m is invertible, in which case m will // be replaced with its inverse; otherwise returns // false, and the contents of m will be scrambled. Make // a copy of m if you wish the preserve its original // value. static bool invert(Matrix m) { ... } ... } This is a sort of low-level operation for performance-sensitive code, so the "nicer" invert operation can be implemented in terms of this one: struct Matrix(T) { ... static bool invert(Matrix m) { ... } ... property Matrix inverse() { auto result = this.dup; // make a copy first if (!invert(result)) throw new Exception("Tried to invert non-invertible matrix"); return result; } } Matrix!(Array!(2,float)) A = ...; auto B = A.inverse; // this returns a separate inverse matrix A.invert(A); // this modifies A in-place to become its inverse Of course, these method names are tentative; I'm sure you can think of better names for them. T -- Computers shouldn't beep through the keyhole.

Oct 11 2013

On Friday, 11 October 2013 at 17:49:32 UTC, H. S. Teoh wrote:On Fri, Oct 11, 2013 at 06:10:19PM +0200, FreeSlave wrote:Can you please give links to both libraries?There is "Matrices and linear algebra" module in wish list. Let's discuss its design. D is complicated language so it's difficult to choose the right way here. We need to find compromise between efficiency and convenient interface. I'm going to make some suggestions how this module should look like.I think we need to differentiate between multidimensional arrays (as a data storage type) and linear algebra (operations performed on 2D arrays). These two intersect, but they also have areas that are not compatible with each other (e.g. matrix product vs. element-by-element product). Ideally, we should support both in a clean way. As far as the former is concerned, Denis has implemented a multidimensional array library, and I've independently done the same, with a slightly different interface. I think one or two others have implemented similar libraries as well. It would be good if we standardize the API so that our code can become interoperable.As far as the latter is concerned, I've been meaning to implement a double-description convex hull algorithm, but have been too busy to actually work on it. This particular algorithm is interesting, because it stress-tests (1) performance of D algorithms, and (2) challenges the design of matrix APIs because while the input vertices (resp. hyperplanes) can be interpreted as a matrix, the algorithm itself also needs to permute rows, which means it is most efficient when given an array-of-pointers representation, contrary to the usual flattened representations (as proposed below). I think there's a place for both, which is why we need to distinguish between data representation and the algorithms that work on them.The problem is that algorithms still should know matrix template to provide compile-time checks if possible or throw exceptions at runtime if something gone wrong.First of all, it should provide two templates for matrices. Let's call them StaticMatrix and DynamicMatrix. The first one has "templated" size and therefore may use static arrays and compile-time checks. It can be useful when the size is determined by our needs, for example, in graphics. DynamicMatrix has variable size, i.e. it should be created in heap. It can be useful in all other math areas.I like this idea. Ideally, we should have many possible representations, but all conforming to a single API understood by all algorithms, so that you only have to write algorithms once, and they will work with any data structure. That's one key advantage of D, and we should make good use of it.I do not want to see a repetition of the C++ situation where there are so many different matrix/multidimensional array libraries, and all of them use incompatible representations and you cannot freely pass data from one to algorithms in the other. Then when no library meets precisely what you need, you're forced to reinvent yet another matrix class, which is a waste of time.I'm not sure if you understand idea of differences between inner implementation majority and interface majority. I agree that inner majority should be defined by inner type. Interface majority is just choice between matrix[rowIndex, columnIndex] and matrix[columnIndex, rowIndex] In case of interface majority we just must choose the appropriate one and use it all over the library. It does not relate to performance.But another question arises - which "majority" should we use in interface? Interface should not depend on inner representation. All functions need unambiguity to avoid complication and repetition of design. Well, actually we can deal with different majority in interface - we can provide something like "asTransposed" adapter, that will be applied by functions if needed, but then we will force user to check majority of matrix interface, it's not very good approach.Algorithms shouldn't even care what majority the data representation is in. It should only access data via the standardized matrix API (whatever it is we decide on). The input type should be templated so that *any* type that conforms to this API will work. Of course, for performance-sensitive code, the user should be aware of which representations are best-performing, and make sure to pass in the appropriate type of representations; but we should not prematurely optimize here. Any linear algebra algorithms should be able to work with *any* type that conforms to a standard matrix API.Yes, non-square. My bad. Well, ok. We want to abstract from inner representation to provide freedom for users. We fall in metaprogramming and generic programming here, so we need to define concepts just like Boost/STL/std.range do. The good thing is that in D types with different interfaces and syntax constraints can satisfy same concept that would be impossible or very difficult in C++. Thanks to static if and is(typeof()). For example inner representation type can provide [][] operator or [,] operator and Matrix type will understand both cases. Suppose: template canBeMatrixRepresentation(T) { enum bool canBeMatrixRepresentation = is(typeof( { T t; //default constructable const(T) ct; alias ElementType!T E; //has ElementType E e; //element type is default constructable static if (/*has [,] operator*/) { t[0,0] = e; //can be assigned e = ct[0,0]; //can retrive element value from const(T) } else static if (/*has [][] operator*/) { t[0][0] = e; //can be assigned e = ct[0][0]; //can retrive element value from const(T) } else { static assert(false); } size_t rows = ct.rowNum; //has row number size_t cols = ct.columnNum; //has column number t.rowNum = size_t.init; t.columnNum = size_t.init; })); } We see that two-dimensional D array does not satisfy this concept because it has no rowNum and columnNum so it should be handled separately. This concept is not ideal since not all types may provide variable rowNum and columnNum. Also concept should expose information whether is "static" type or not, so algorithms will know can they use compile-time checks or not. Types also can provide copy constructor. If they do then Matrix will use it, if they don't then Matrix will do element-by-element copy. It also can try .dup property. It's just example how it should work, but I hope the point is clear. We also need Matrix concept (or separate concepts for StaticMatrix and DynamicMatrix).About operation overloading. It's quite clear about 'add' and 'subtruct' operations, but what's about product? Here I think all 'op'-functions should be 'element by element' operations. So we can use all other operations too without ambiguity. For actual matrix multiplication it can provide 'multiply' or 'product' function. It's similar to Maxima approach, besides Maxima uses dot notation for these needs.Here is where we see the advantage of separating representation from algorithm. Technically, a matrix is not the same thing as a 2D array, because a matrix has a specific interpretation in linear algebra, whereas a 2D array is just a 2D container of some elements. My suggestion would be to write a Matrix struct that wraps around a 2D array, and provides / overrides the overloaded operators to have a linear algebra interpretation. So, a 2D array type should have per-element operations, but once wrapped in a Matrix struct, it will acquire special matrix algebra operations like matrix products, inversion, etc.. In the most general case, a 2D array should be a specific instance of a multidimensional array, and a Matrix struct should be able to use any underlying representation that conforms to a 2D array API. For example: // Example of a generic multidimensional array type struct Array(int dimension, ElemType) { ... Array opBinary(string op)(Array x) { // implement per-element operations here } } // A matrix wrapper around a 2D array type. struct Matrix(T) if (is2DArray!T) { T representation; Matrix opBinary(string op)(Matrix x) if (op == "*") { // implement matrix multiplication here } Matrix opBinary(string op)(Matrix x) if (op != "*") { // forward to representation.opBinary to default // to per-element operations } // Provide operations specific to matrices that don't // exist in general multidimensional arrays. Matrix invert() { ... } } Array!(2,float) myArray, myOtherArray; auto arrayProd = myArray * myOtherArray; // per-element multiplication auto A = Matrix(myArray); // wrap array in Matrix wrapper auto B = Matrix(myOtherArray); auto C = A * B; // matrix product The idea of the Matrix struct here is that the user should be free to choose any underlying matrix representation: a 1D array in row-major or column-major representation, or a nested array of arrays, or a sparse array with some other kind of representation. As long as they provide a standard way of accessing array elements, Matrix should be able to accept them, and provide matrix algebra semantics for them.Transposition. I've already mentioned "asTransposed" adapter. It should be useful to make matrix feel like transposed without its copying. We also can implement 'transpose' and 'transposed' functions. The first one transposes matrix in place. It's actually not allowed for non-square StaticMatrix since we can't change the size of this type of matrices at runtime. The second one returns copy so it's applicable in all cases. Actually I'm not sure should these functions be member functions or not.The most generic approach to transposition is simply a reordering of indices. This difference is important once you get to 3D arrays and beyond, because then there is no unique transpose, but any permutation of array indices should be permissible. Denis' multidimensional arrays have a method that does O(1) reordering of array indices: basically, you create a "view" of the original array that has its indices swapped around. So there is no data copying; it's just a different "view" into the same underlying data. This approach of using "views" rather than copying data allows for O(1) submatrix extraction: if you have a 50x50 matrix, then you can take arbitrary 10x10 submatrices of it without needing to copy any of the data, which would be very expensive. Avoiding unnecessary copying becomes very important when the dimension of the array increases; if you have a 3D or 5D array, copying subarrays become extremely expensive very quickly. A .dup method should be provided in the cases where you actually *want* to copy the data, of course. Basically, subarrays / transpositions / index reordering should be regarded as generalizations of D's array slices. No data should be copied until necessary.Invertible matrix. It must not be allowed for square StaticMatrix.You mean for non-square StaticMatrix?

Oct 11 2013

On Fri, Oct 11, 2013 at 09:49:22PM +0200, FreeSlave wrote:On Friday, 11 October 2013 at 17:49:32 UTC, H. S. Teoh wrote:I seem to have lost the link to Denis' code, but you can see the docs here: http://denis-sh.bitbucket.org/unstandard/unstd.multidimensionalarray.html I haven't put my code up in any public repo (yet), because it's not quite ready for public consumption. But you can see the docs (somewhat older version) here: http://eusebeia.dyndns.org/~hsteoh/tmp/doc/mda I'm not 100% satisfied with its API design (spanningIndices probably introduces an unacceptable bottleneck in many cases), though it does employ the principle of generalizing D array slices, so all instances of Arrays are "views" into some underlying data. Assignment, subarrays, etc., are therefore cheap and have reference semantics. Subarrays can be employed to allow block-initialization of arrays, etc.. A .dup method is provided to make actual copies of data when needed. There are some ugly hacks (IndexRange) in order to work around language limitations until Kenji's pull for implementing multidimensional slicing is committed. My code is also not as flexible as Denis' when it comes to index reordering, but it does work in the 2D case. Some of the unittests successfully demonstrate that the code is generic enough to handle arrays of non-numeric types, such as arrays of strings, in which case per-element unary/binary operations work as expected (e.g., A ~ "a" appends "a" to all elements, when A is an array of strings). [...]On Fri, Oct 11, 2013 at 06:10:19PM +0200, FreeSlave wrote:Can you please give links to both libraries?There is "Matrices and linear algebra" module in wish list. Let's discuss its design. D is complicated language so it's difficult to choose the right way here. We need to find compromise between efficiency and convenient interface. I'm going to make some suggestions how this module should look like.I think we need to differentiate between multidimensional arrays (as a data storage type) and linear algebra (operations performed on 2D arrays). These two intersect, but they also have areas that are not compatible with each other (e.g. matrix product vs. element-by-element product). Ideally, we should support both in a clean way. As far as the former is concerned, Denis has implemented a multidimensional array library, and I've independently done the same, with a slightly different interface. I think one or two others have implemented similar libraries as well. It would be good if we standardize the API so that our code can become interoperable.We should simply have standard APIs for checking for compile-time available properties, similar to std.range.hasLength. Then algorithms that care about the distinction can make use of it, and other algorithms that don't care will automatically work with both variants. [...]The problem is that algorithms still should know matrix template to provide compile-time checks if possible or throw exceptions at runtime if something gone wrong.First of all, it should provide two templates for matrices. Let's call them StaticMatrix and DynamicMatrix. The first one has "templated" size and therefore may use static arrays and compile-time checks. It can be useful when the size is determined by our needs, for example, in graphics. DynamicMatrix has variable size, i.e. it should be created in heap. It can be useful in all other math areas.I like this idea. Ideally, we should have many possible representations, but all conforming to a single API understood by all algorithms, so that you only have to write algorithms once, and they will work with any data structure. That's one key advantage of D, and we should make good use of it.Sorry, I misunderstood what you meant. I thought you were talking about storage order, but apparently you were talking about index order in the API, which is a different issue. My personal preference is columnIndex first, then rowIndex, but that is opposite of standard matrix notation in math (which is row-first then column). I think we should just choose one way or the other, and just stick with it consistently. The important thing is to be consistent to prevent bugs; exactly which order is the "right" way is just bikeshedding, I think. There's another related point here, related to declaring array dimensions. In my multi-dim array code, I chose to have the index order match the order of the specified array dimensions; e.g., if you declare: auto A = Array!(2,int)([2, 5]); then you'd index A from A[0,0] to A[1,4]. However, this is opposite of what nested arrays in D does: int[2][5] B; Here, the index range for B goes from B[0][0] to B[4][1], because B is a 5-element array of 2-element subarrays. Again, I'm not sure what the "correct" way is; I chose my way because it's easier to remember (the order you specify dimensions is the same order the indices appear in, you don't have to keep remembering when to reverse the order). But either way, a single standard approach should be chosen, and we should stick with it consistently. [...]I'm not sure if you understand idea of differences between inner implementation majority and interface majority. I agree that inner majority should be defined by inner type. Interface majority is just choice between matrix[rowIndex, columnIndex] and matrix[columnIndex, rowIndex] In case of interface majority we just must choose the appropriate one and use it all over the library. It does not relate to performance.But another question arises - which "majority" should we use in interface? Interface should not depend on inner representation. All functions need unambiguity to avoid complication and repetition of design. Well, actually we can deal with different majority in interface - we can provide something like "asTransposed" adapter, that will be applied by functions if needed, but then we will force user to check majority of matrix interface, it's not very good approach.Algorithms shouldn't even care what majority the data representation is in. It should only access data via the standardized matrix API (whatever it is we decide on). The input type should be templated so that *any* type that conforms to this API will work. Of course, for performance-sensitive code, the user should be aware of which representations are best-performing, and make sure to pass in the appropriate type of representations; but we should not prematurely optimize here. Any linear algebra algorithms should be able to work with *any* type that conforms to a standard matrix API.Well, ok. We want to abstract from inner representation to provide freedom for users. We fall in metaprogramming and generic programming here, so we need to define concepts just like Boost/STL/std.range do. The good thing is that in D types with different interfaces and syntax constraints can satisfy same concept that would be impossible or very difficult in C++. Thanks to static if and is(typeof()). For example inner representation type can provide [][] operator or [,] operator and Matrix type will understand both cases.Suppose: template canBeMatrixRepresentation(T)<bikeshed> What about a nicer, more concise name, like is2DArray? :) </bikeshed>{ enum bool canBeMatrixRepresentation = is(typeof( { T t; //default constructable const(T) ct; alias ElementType!T E; //has ElementType E e; //element type is default constructable static if (/*has [,] operator*/) { t[0,0] = e; //can be assigned e = ct[0,0]; //can retrive element value from const(T) } else static if (/*has [][] operator*/) { t[0][0] = e; //can be assigned e = ct[0][0]; //can retrive element value from const(T) } else { static assert(false); } size_t rows = ct.rowNum; //has row number size_t cols = ct.columnNum; //has column numberI'm uneasy about this one. This is too specific to 2D arrays, and would not work with general multidimensional arrays (that just happen to be 2D). I think it's better to use opDollar: size_t rows = ct.opDollar!0; // assume row-major indexing size_t cols = ct.opDollar!1; Requiring user types to implement opDollar has the advantage that then we can provide nice $ notation inside the slicing operator: auto topLeftElem = myMatrix[0,0]; auto bottomRightElem = myMatrix[$-1, $-1]; [...]We see that two-dimensional D array does not satisfy this concept because it has no rowNum and columnNum so it should be handled separately. This concept is not ideal since not all types may provide variable rowNum and columnNum. Also concept should expose information whether is "static" type or not, so algorithms will know can they use compile-time checks or not.I think we should use opDollar instead of arbitrarily-named custom fields like rowNum or columnNum. Now, opDollar doesn't quite support built-in nested arrays, so I'm not sure how to fix that; perhaps provide some wrappers in std.array? It's kinda annoying, though, because nested arrays may be jagged, in which case it cannot be treated as a matrix. But rejecting all nested arrays outright seems a bit too limited. Maybe provide a wrapper to convert nested arrays into "real" 2D arrays (with runtime checking against jaggedness)? That also lets us check if a particular dimension is compile-time known or not: template isStaticDim(Array, int dim) { enum isStaticDim = is(typeof({ // If it's possible to put opDollar!dim into an // enum, then it's statically-known, otherwise // it's dynamic. enum length = Array.init.opDollar!dim; })); } StaticMatrix!(4,4) sm; // both dimensions known at compile-time DynamicMatrix dm(4,4); // both dimensions specified at runtime RowAppendableMatrix!4 am; // one dimension specified at runtime assert(isStaticDim!(typeof(sm), 0)); assert(isStaticDim!(typeof(sm), 1)); assert(!isStaticDim!(typeof(dm), 0)); assert(!isStaticDim!(typeof(dm), 1)); assert(isStaticDim!(typeof(am), 0)); assert(!isStaticDim!(typeof(am), 1));Types also can provide copy constructor. If they do then Matrix will use it, if they don't then Matrix will do element-by-element copy. It also can try .dup property.I think it should be OK to standardize on ".dup", by analogy with built-in arrays. If no such method is provided, then do element-by-element copy.It's just example how it should work, but I hope the point is clear. We also need Matrix concept (or separate concepts for StaticMatrix and DynamicMatrix).Looks like a hierarchy of concepts: Matrix / \ StaticMatrix DynamicMatrix The Matrix concept provides the common core APIs for accessing matrices; the more specific concepts StaticMatrix and DynamicMatrix build on top of Matrix by adding more API methods/properties specific to static or dynamic matrices. If an algorithm doesn't care which kind of matrix is passed to it, it can simply take Matrix (via an isMatrix signature constraint, ala std.algorithm.* specifying isInputRange, etc.). If it wants a specific kind of matrix, it specifies that instead (e.g. isStaticMatrix). If it can handle both, it can use static if to treat each case separately: auto myAlgo(M)(M matrix) if (isMatrix!M) // accept most general form of matrix { ... // code static if (isStaticMatrix!M) ... // stuff specific to static matrices else static if (isDynamicMatrix!M) ... // stuff specific to dynamic matrices else static assert(0); // oops ... // more code } T -- VI = Visual Irritation

Oct 11 2013

First thing. I see std.linarg I have no clue whatsoever what it even may be about. Do we really want to have a standard lib with names that look like base64 encoded ?

Oct 11 2013

On Sat, Oct 12, 2013 at 12:36:08AM +0200, deadalnix wrote:First thing. I see std.linarg I have no clue whatsoever what it even may be about. Do we really want to have a standard lib with names that look like base64 encoded ?I picked it up immediately. It's an obvious abbreviation of "linear algebra". I'm not sure I understand this resentment against abbreviated names. By your argument, we should rename std.xml to std.extensibleMarkupLanguage instead. Bring on the rainbow, the bikeshed is up for painting. T -- You are only young once, but you can stay immature indefinitely. -- azephrahel

Oct 11 2013

On Fri, 11 Oct 2013 15:42:17 -0700, H. S. Teoh wrote:On Sat, Oct 12, 2013 at 12:36:08AM +0200, deadalnix wrote:Don't be silly. It'd obviously have to be "standardLibrary.extensibleMarkupLanguage"First thing. I see std.linarg I have no clue whatsoever what it even may be about. Do we really want to have a standard lib with names that look like base64 encoded ?I picked it up immediately. It's an obvious abbreviation of "linear algebra". I'm not sure I understand this resentment against abbreviated names. By your argument, we should rename std.xml to std.extensibleMarkupLanguage instead. Bring on the rainbow, the bikeshed is up for painting. T

Oct 11 2013

On 10/11/13 3:42 PM, H. S. Teoh wrote:On Sat, Oct 12, 2013 at 12:36:08AM +0200, deadalnix wrote:Honestly I first thought it was "linear-time algorithms". Andrei

Oct 11 2013

On Fri, Oct 11, 2013 at 04:13:41PM -0700, Andrei Alexandrescu wrote:On 10/11/13 3:42 PM, H. S. Teoh wrote:[...] But I thought that was the province of std.algorithm? BTW, are we going to split up std.algorithm anytime soon, now that we have package.d? Jonathan is already working on splitting up std.datetime into more manageable pieces, and it would be nice if std.algorithm can be broken up into saner pieces too. (The last time I brought it up a few people supported it, but I haven't had the time to do anything about it so far.) T -- The easy way is the wrong way, and the hard way is the stupid way. Pick one.On Sat, Oct 12, 2013 at 12:36:08AM +0200, deadalnix wrote:Honestly I first thought it was "linear-time algorithms".

Oct 11 2013

On Friday, 11 October 2013 at 23:13:34 UTC, Andrei Alexandrescu wrote:Honestly I first thought it was "linear-time algorithms". AndreiI thought it was some kind of algae species. Why not just call it std.algebra?

Oct 11 2013

On Friday, 11 October 2013 at 23:24:00 UTC, Andrej Mitrovic wrote:On Friday, 11 October 2013 at 23:13:34 UTC, Andrei Alexandrescu wrote:Because, amongst people who know, numerical linear algebra algorithms form their own little world of computing apart from what may be called 'algebra'. 'Algebra' could mean any number of things, like symbolic algebra, but the entire world of MATLAB/R programmers know exactly what linalg might be. Numerical linear algebra. See the reference text by Golub and Van Loan. I would try to fit this module under some larger numeric or scientific computing hierarchy, but I don't see one in D yet. The Julia language is the competitor here, and it is already well ahead of D, but who knows how the race will end? -- BrianHonestly I first thought it was "linear-time algorithms". AndreiI thought it was some kind of algae species. Why not just call it std.algebra?

Oct 11 2013

On Friday, 11 October 2013 at 22:36:09 UTC, deadalnix wrote:

Oct 12 2013

On Friday, 11 October 2013 at 16:10:21 UTC, FreeSlave wrote:Please, make your suggestions too.Firstly, my opinion is, it shouldn't be in the std package. Secondly, if it's not super highly performant, noone will use it. The highest performing linalg libraries implement the low level BLAS API or are vendor specific, like the Intel MKL. These are super optimized libraries where the asm code is hand optimized by specialists (GotoBLAS and the derived OpenBLAS for instance), and these will be *very* hard to beat. So much so that if BLAS isn't used as the engine, it's very unlikely that your library will be used for serious number crunching. Built above BLAS, you can find a number of libraries that make using linear algebra more pleasant. One such library, currently the best IMHO is Armadillo C++. It has the very big advantage that it mimics Matlab functions. This feature is extremely useful, because scientists very often model their algorithms in Matlab before porting them to C++. This is what is done in the company I work for. So the best I could suggest is to study the Armadillo C++ library and port it to D.

Oct 11 2013

On Saturday, 12 October 2013 at 05:20:11 UTC, SomeDude wrote:On Friday, 11 October 2013 at 16:10:21 UTC, FreeSlave wrote:For these cases we may let users to choose low-level backend if they need. High-level interface and default implementation are needed anyway. I called it std.linalg because there is website http://www.linalg.org/ about C++ library for exact computational linear algebra. Also SciD has module scid.linalg. We can use std.linearalgebra or something else. Names are not really important now. Ok, things are more clear now. Today I look what I can do.Please, make your suggestions too.Firstly, my opinion is, it shouldn't be in the std package. Secondly, if it's not super highly performant, noone will use it. The highest performing linalg libraries implement the low level BLAS API or are vendor specific, like the Intel MKL. These are super optimized libraries where the asm code is hand optimized by specialists (GotoBLAS and the derived OpenBLAS for instance), and these will be *very* hard to beat. So much so that if BLAS isn't used as the engine, it's very unlikely that your library will be used for serious number crunching. Built above BLAS, you can find a number of libraries that make using linear algebra more pleasant. One such library, currently the best IMHO is Armadillo C++. It has the very big advantage that it mimics Matlab functions. This feature is extremely useful, because scientists very often model their algorithms in Matlab before porting them to C++. This is what is done in the company I work for. So the best I could suggest is to study the Armadillo C++ library and port it to D.

Oct 11 2013

On Saturday, 12 October 2013 at 06:24:58 UTC, FreeSlave wrote:For these cases we may let users to choose low-level backend if they need. High-level interface and default implementation are needed anyway. I called it std.linalg because there is website http://www.linalg.org/ about C++ library for exact computational linear algebra. Also SciD has module scid.linalg. We can use std.linearalgebra or something else. Names are not really important now. Ok, things are more clear now. Today I look what I can do.There are litterally dozens of linear algebra packages: Eigen, Armadillo, Blitz++, IT++, etc. I was not complaining about the linalg name, but about the fact that you want to make it a std subpackage. I contend that if you want to make it a std package, it must be nearly perfect, i.e better performing than ALL the other alternatives, even the C++ ones, and that it's really good as an API. Else it will be deprecated because someone will have made a better alternative. Given the number of past tries, I consider this project is very likely doomed to failure. So no std please.

Oct 12 2013

On Saturday, 12 October 2013 at 08:47:33 UTC, SomeDude wrote:On Saturday, 12 October 2013 at 06:24:58 UTC, FreeSlave wrote:It's not my idea to include this kind of module into std. I found it in wish list here http://wiki.dlang.org/Wish_list so I supposed community needs it and it should be discussed and designed, and then implemented by someone, not necessarily by me, though I can and will try to do it by myself. You're right, we should to learn existing solutions firstly and then make the best in D.For these cases we may let users to choose low-level backend if they need. High-level interface and default implementation are needed anyway. I called it std.linalg because there is website http://www.linalg.org/ about C++ library for exact computational linear algebra. Also SciD has module scid.linalg. We can use std.linearalgebra or something else. Names are not really important now. Ok, things are more clear now. Today I look what I can do.There are litterally dozens of linear algebra packages: Eigen, Armadillo, Blitz++, IT++, etc. I was not complaining about the linalg name, but about the fact that you want to make it a std subpackage. I contend that if you want to make it a std package, it must be nearly perfect, i.e better performing than ALL the other alternatives, even the C++ ones, and that it's really good as an API. Else it will be deprecated because someone will have made a better alternative. Given the number of past tries, I consider this project is very likely doomed to failure. So no std please.

Oct 12 2013

I'm with SomeDude: This project is doomed to failure from the get go, and under no circumstances should be put into std. First, some background so this doesn't seem too harsh. Numerical linear algebra is easily the single oldest and most well-developed area of computing. Period. There's a huge history there, and if you don't know it then whatever you make won't be used and will be universally derided. Second, beyond the background of implementations is the theoretical background of what they're used for and how that influences the API. You also have understand that linear algebra algorithms and data structures come in 3 major varieties: Dense, Sparse Direct, and Iterative methods. For just the theoretical background on these and associated floating point issues there is an absolute minimum of the material similar to what's covered in these two standard references: Matrix Computations by Golub and Van Loan: http://www.amazon.com/Computations-Hopkins-Mathematical-Sciences-ebook/dp/B00BD2DVIC/ref=sr_1_6?ie=UTF8&qid=1381614419&sr=8-6&keywords=nick+higham Accuracy and Stability of Numerical Algorithms by Nick Higham: http://www.amazon.com/Accuracy-Stability-Numerical-Algorithms-Nicholas/dp/0898715210/ref=sr_1_7?ie=UTF8&qid=1381614419&sr=8-7&keywords=nick+higham The old, well-tested, and venerable BLAS and LAPACK APIs/packages are just for dense matrices. There are other variants for different kinds of parallelism. And notice there are differences in numerical linear algebra between instruction level parallelism (like SSE and AVX), shared memory, distributed memory, and coprocessors (like OpenCL or CUDA). In terms of popular existing implementations there are more recent packages like Magma (http://icl.utk.edu/magma/), PLASMA(http://icl.cs.utk.edu/plasma/), PETSc (http://www.mcs.anl.gov/petsc/), Trillinos (http://trilinos.sandia.gov/), or Elemental (http://libelemental.org/). The most popular sparse direct package, as far as I know, is SuiteSparse (http://www.cise.ufl.edu/research/sparse/SuiteSparse/), which is wrapped by Eigen (http://eigen.tuxfamily.org/index.php?title=Main_Page). Eigen is so popular in part because it does a good job of balancing expression template magic for dense matrices and wrapping/calling BLAS/LAPACK and SuiteSparse. The very popular NumPy and SciPy packages are also basically wrappers around BLAS and LAPACK. There is also a scipy.sparse package, but its API differences are an ongoing headache. All of which is to say this is a very deep area, the differences in algorithms and storage are sufficiently fundamental to make a unified API very difficult, and in any case the D community simply cannot commit to developing and maintaining such a library or commit to rolling in future developments from ongoing research in the area. This is exactly the sort of thing that should be a third-party extension, like in every other language, rather than in the standard library. On Saturday, 12 October 2013 at 08:47:33 UTC, SomeDude wrote:On Saturday, 12 October 2013 at 06:24:58 UTC, FreeSlave wrote:For these cases we may let users to choose low-level backend if they need. High-level interface and default implementation are needed anyway. I called it std.linalg because there is website http://www.linalg.org/ about C++ library for exact computational linear algebra. Also SciD has module scid.linalg. We can use std.linearalgebra or something else. Names are not really important now. Ok, things are more clear now. Today I look what I can do.There are litterally dozens of linear algebra packages: Eigen, Armadillo, Blitz++, IT++, etc. I was not complaining about the linalg name, but about the fact that you want to make it a std subpackage. I contend that if you want to make it a std package, it must be nearly perfect, i.e better performing than ALL the other alternatives, even the C++ ones, and that it's really good as an API. Else it will be deprecated because someone will have made a better alternative. Given the number of past tries, I consider this project is very likely doomed to failure. So no std please.

Oct 12 2013

On Friday, 11 October 2013 at 16:10:21 UTC, FreeSlave wrote:There is "Matrices and linear algebra" module in wish list. Let's discuss its design. D is complicated language so it's difficult to choose the right way here. We need to find compromise between efficiency and convenient interface. I'm going to make some suggestions how this module should look like. First of all, it should provide two templates for matrices. Let's call them StaticMatrix and DynamicMatrix. The first one has "templated" size and therefore may use static arrays and compile-time checks. It can be useful when the size is determined by our needs, for example, in graphics. DynamicMatrix has variable size, i.e. it should be created in heap. It can be useful in all other math areas. Both templates should support all floating point types and moreover user-defined (for example wrappers for GMP library and others). For efficiency in both cases matrices should use one-dimensional array for inner representation. But actually I'm not sure if matrices should support other container types besides standard D arrays. The good thing about one-dimensional arrays is that they can be easily exposed to foreign functions, for example, to C libraries and OpenGL. So we should take care about memory layout - at least row-major and column-major. I think it can be templated too. But another question arises - which "majority" should we use in interface? Interface should not depend on inner representation. All functions need unambiguity to avoid complication and repetition of design. Well, actually we can deal with different majority in interface - we can provide something like "asTransposed" adapter, that will be applied by functions if needed, but then we will force user to check majority of matrix interface, it's not very good approach. Sometimes user takes data from some other source and wants to avoid copying in Matrix construction, but she also wants to get matrix functionality. So we should provide "arrayAsMatrix" adapter, that can adopt one-dimensional and two-dimensional arrays making them feel like matrices. It definitely should not make copy of dynamic array, but I'm not sure about static. About operation overloading. It's quite clear about 'add' and 'subtruct' operations, but what's about product? Here I think all 'op'-functions should be 'element by element' operations. So we can use all other operations too without ambiguity. For actual matrix multiplication it can provide 'multiply' or 'product' function. It's similar to Maxima approach, besides Maxima uses dot notation for these needs. Transposition. I've already mentioned "asTransposed" adapter. It should be useful to make matrix feel like transposed without its copying. We also can implement 'transpose' and 'transposed' functions. The first one transposes matrix in place. It's actually not allowed for non-square StaticMatrix since we can't change the size of this type of matrices at runtime. The second one returns copy so it's applicable in all cases. Actually I'm not sure should these functions be member functions or not. Invertible matrix. It must not be allowed for square StaticMatrix. DynamicMatrix may make checks at runtime. Same as for transposition we can implement 'invert' to invert in place and 'inverted' to make copy. There is issue here - what should we do when determinant is 0? I believe the best approach here is to throw exception since if user needs invertible matrix it's definitely exception when it can't be calculated. Please, make your suggestions too.I'd like to echo some of the points made by CJS. Rolling your own linear algebra (for anything other than very small matrices) is an exercise in futility unless you are prepared to become a world-class expert in what is an incredibly mature and detailed field. Clever wrapping of industry standard libraries is the way to go. I would love to see this https://github.com/cristicbz/scid taken further, which was a fork of scid that used clever magic to reduce matrix expressions to an semi-optimal set of blas calls. It doesn't appear to have got much love in recent times.

Oct 12 2013