Re: [boost] different matrix library?

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 : [phone]