Larger support of structured matrices

[Benoît Lacelle]

Hi,

With Eric Schost, we are developing a library for linear algebra on
structured matrices. Actually, it works only on Toeplitz-like matrices
for the moment. We believe that the LinBox handle only some exact
structures (Toeplitz at least). So we thought it could be interesting
to integrate our developments in the LinBox.

Do you think it is actually a good idea ?
If yes, do you advice a particular way to integrate our library
knowing it is quite big (~10k lines of code) ?

Benoit Lacelle

[Dave Saunders]

Hi, Benoit,

This sounds interesting.

As you say, LinBox has just a few examples of structured matrices, no
displacement rank manipulations, etc. Given a fast blackbox matrix
(having a softO(n) matrix vector product), LinBox provides fast
(softO(n^2)) blackbox algorithms for many standard matrix problems
(rank, det, solve, minpoly, etc) which can be compared to the superfast
(softO(n) for a given displacement rank), algorithms that I suppose you
are implementing. It is likely that there will be a threshold of matrix
size below which an existing LinBox method is good and beyond which a
superfast method is best. Thus linking your library could serve
multiple purposes: algorithm comparisons, hybrid algorithms, and
presenting a common interface for dense, sparse, structured matrices.

Will you be at ISSAC '09? We could discuss it there. In any case, can
you give me a pointer to your library. I'd be happy to consider how we
might integrate.

Best, -dave

[Benoît Lacelle]

Hi Dave,

Thanks for your prompt reply. I won't be at ISSAC.
You can download current version of our library there:
http://benoitlacelle.fr/files/SLA_090722.zip
To be used, it needs the last version of NTL.

We focused on the resolution of systems and the inversion of matrices
in the fields K=Z/pZ and K[X].

here is a quick overview of current architecture:
About file:
files like "c*.hpp" define classes
files like "tools_*.hpp" defines helping functions
most algorithms are in "invers*.hpp"
files like "NTL_*.hpp" and "proxy_NTL_*.hpp" make the bridge with
the NTL

Main classes:
IStructured.hpp : Abstract class for structured matrices (defined
by generators)
cStructured.hpp : templated class for structured matrices where
generators are defined by two matrices
cKaltofen.hpp : templated class for structured matrices where
generators are defined by two vectors of polynomials
cPolynom.hpp : defines the templated class cRichPolynom: represent
a polynomial by both a dense representation and a FFT representation

Main algorithms:
inverse_MBA.hpp: Implementation of the MBA algorithm, works on
cKaltofen
inverse_quadratic.hpp: Quadratic algorithms for the resolution and
inversion of matrices, works on cStructured
inverse_polynomial.hpp: Implementation of a Structured Newton's
Iteration. Works on cStructured
reduction.hpp: reduce the number of generators defining a
structured matrix, Works on IStructured
tools_generators: convert generators between displacement operators
tools_multiply.hpp: all multiplications for the case K = Z/pZ and
the displacement operator is phi+ or phi-
mul_matrix_vector.hpp: some multiplications for the case K = Z/pZ
and K[X] and the displacement operator is a sylvestor operator.

We tried to put some genericity in our library through templates. But
some features are really dependant on the way NTL works. On the whole,
we expect objects with the same interface as NTL ones. We represent
the hierarchy of math objects (Field (K) -> PolynomialField (K[X]))
through templated typedefs (td_polynom in cPolynom.hpp, td_matrix in
cMatrix.hpp).

About a possible integration, (very big) changes will have to be done
about the classes representing structured matrices. We will have time
to discuss about that but current architecture is nothing but "good
enough for current state of affairs". Main points being "how should I
have access to the generators?" and "how should I have the knowledge
about the displacement operator?". Of course, I am open to remarks and
advice, mainly about the general architecture.

Benoit