Subject: Re: [boost] different matrix library?
From: joel (joel.falcou_at_[hidden])
Date: 2009-08-17 02:45:13
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