Rutger ter Borg wrote:
With that I meant that every linear algebra library tends to provide its own implementation of vector and matrix types. Sometimes with different access semantics. One could use std::vector for dense vectors, (unordered) map for sparse vectors/matrices. Perhaps what's missing in the standard containers are matrices? In other words, I would prefer using a traits system to access any container.
Suppose if only standard containers were to be used, perhaps
let(C) = a*A*B + b*C
will be enough to inject template expressions?
NT2 as a as_matrix<T> function std::vector<float> u(100*100); // fill u somehow; matrix<float> r; r = cos(as_matrix(u)) * as_matrix(u); No copy of course of the original vector
I've also noted that matrix libraries are usually not exhaustive with their matrix classes. There's all kinds of stuff, like symmetric, triangular, diagonal, banded, bidiagonal, tridiagonal, .... Adaptors / math property wrapper come in handy too, if you have another container which happens to be a symmetric matrix or a symmetric matrix which happens to be positive definite.
E.g., suppose A is dense,
x = solve( A, b ) x = solve( positive_definite( upper( A ) ), b )
in the second case you say to use the upper triangular part, and make that pos def. in my work for writing a python translator from LAPACK's Fortran source base to C++ bindings code, I've come across the specializations available in LAPACK, that might give a clue for an exhaustive list and write down a good set of orthogonal concepts (e.g., storage traits, structure traits, math properties, etc.).
Policies & lazy-typer. x = solve( as_matrix<settings(positive_definite)>(A), b );
Although this sounds good, I think in some cases there will be multiple matches of available algorithms. Then the specialization for the expression with the highest numerical complexity should "win" (assuming the gain is highest there). Is this trivial, too?
This would be really neat, I think an entire area of research is writing algorithms for special cases.
Introspection on the AST properties of matrix , checking for pattern OR settings usage somewhere. NT2 solve for example is a mapping over the 20+ LAPACK solver and use CT/RT checks to select the best one.
Given the previous point, it's trivial to plug-in BLAS and LAPACK, so plus minus some handling of temporaries, we might already be almost there.
Well, don't forget SIMDifed C is faster than FORTRAN ;)
What I forgot to mention was error bounds and error propagation. For many users this may be more important than speed. Many algorithms are specialized for this purpose only (same speed, but higher precision).
Just drop them in various namespaces. NT2 has plan to have a set of toolbox in which the speed/precision trade-off avry while the base namespace has a middle ground implementation. There is nothing ahrd per se in those, it's just tedious considering the large number of variation. -- ___________________________________________ Joel Falcou - Assistant Professor PARALL Team - LRI - Universite Paris Sud XI Tel : (+33)1 69 15 66 35