text: 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.
  -- ___________________________________________