1 Introduction
Structured matrices are two-dimensional arrays which
can be represented by O(n) values, where n is the size of the matrix.
This class provides classical arithmetic operations such
as addition, substraction, multiplication by scalars, vectors, matrices
...The product of two structured matrices is usually not a structured
matrices (at least of the same type) and therefore is not implemented.
The product of a structured matrices by a vector can usually be performed
in almost linear time, ie. in O(n log(n)) ops (or
O(nlog(n)log2(n) according to the field), using for instance FFT.
See (X) and [BP94].
2 Implementation
template< class R> struct MatStruct
Structured matrices.
{
R rep;
typedef typename R::size_type size_type;
typedef typename R::value_type value_type;
typedef MatStruct self_t;
MatStruct() {}
MatStruct(const R & r): rep(r) {}
MatStruct(int m, int n) : rep(n,m) {}
MatStruct(int m, int n, value_type* t):rep(m,n,t){}
MatStruct(const char * str): rep(str) {}
//Copy constructor
MatStruct(const self_t & M): rep(M.rep) {}
template
MatStruct(const VAL & M) {assign(*this,M.rep);}
//Assignement
self_t & operator=(const self_t & M) {rep=M.rep; return *this;}
template
self_t & operator=(const VAL & M) {assign(*this,M.rep); return *this;}
size_type nbrow() const {return rep.nbrow();}
size_type nbcol() const {return rep.nbcol();}
value_type operator() (size_type i, size_type j) const {return rep(i,j);}
VAL > operator()(Range I, Range J);
self_t & operator+=(const self_t & a);
self_t & operator-=(const self_t & a);
self_t & operator*=(const value_type & c);
self_t & transpose() {rep.transpose(); return *this;}
};
No modification operator (i,j) is implemented.
template < class R>
inline VAL< OP< '+',MatStruct< R> ,MatStruct< R> > >
operator+(const MatStruct< R> & a, const MatStruct< R> & b)
Addition operator.
template < class R>
inline VAL< OP< '-',MatStruct< R> ,MatStruct< R> > >
operator-(const MatStruct< R> & a, const MatStruct< R> & b)
Substraction operator.
template < class R>
inline VAL< OP< '-',MatStruct< R> ,MatStruct< R> > >
operator-(const MatStruct< R> & a)
Unary substraction operator.
template < class R,class S> inline
VAL< OP< '*',MatStruct< R> ,VectStd< S> > >
operator*(const MatStruct< R> & a, const VectStd< S> & b)
Matrix vector mutiplication.
template < class R> inline
VAL< OP< '.',MatStruct< R> ,typename MatStruct< R> ::value_type> >
operator*(const MatStruct< R> & a, const typename MatStruct< R> ::value_type & c)
Scalar multiplication.
template < class R> inline
VAL< OP< 's',MatStruct< R> ,Range2d> > Row(const MatStruct< R> & a, const Range & r)
Submatrix from a range of rows.
template < class R> inline
VAL< OP< 's',MatStruct< R> ,Range2d> > Col(const MatStruct< R> & a, const Range & r)
Submatrix from a range of columns.
template< class R>
ostream & operator< < (ostream & os, const MatStruct< R> & p)
Output operator for standard matrices.
template< class R> inline
istream & operator> > (istream & is, MatStruct< R> & M)
Input operator for standard matrices. The input format is of the form
m n c00 c0 ... where m is the number of rows, n is the number of
columns and the coefficients are read by rows.
3 A container for Toeplitz matrices
A matrix T = (ti,j) is a Toeplitz matrix if for
all i, j, the entry ti,j depends only on i-j, that is if ti, j=
ti+1,j+1 for all pairs of (i, j) and (i+1,j+1) for which the
entries ti,j and ti+1,j+1 are defined. Associated with a
Toeplitz matrix, we have the Toeplitz operator which is the
projection of the multiplication by a univariate polynomial in K[x] (see
[MoPa99strmt]). An interesting point of this definition is that
the product of a n × n Toeplitz matrix by a vector
is a subvector of the coefficient vector of the product of two polynomials of
K[x]. This product can be done within O(nlog(n)) (or
O(nlog(n)log2(n) according to the field) operations using FFT.
See (X) and [BP94].
template< class C>
struct toeplitz : public array1d< C>
Containers for Toeplitz matrices, as a subclass of array1d.
Defined by its number of row nbrow_, its number of columns nbcol_.
The elements of the first row and first column are stored, starting
from (0,nbcol) up to (nbrow,0).
The product of a Toeplitz matrix by a vector is performed by FFT
computations (X).
{
unsigned int nbrow_, nbcol_;
typedef C coeff_t;
typedef array1d row_t;
typedef array1d col_t;
toeplitz(): array1d(){}
toeplitz(const toeplitz & v):array1d(v),
nbrow_(v.nbrow_),nbcol_(v.nbcol_){}
toeplitz(unsigned int i, unsigned int j):
array1d(i+j-1),nbrow_(i),nbcol_(j) {}
toeplitz(unsigned int i, unsigned int j, C* t):
array1d(i+j-1,t),nbrow_(i),nbcol_(j){}
void reserve(unsigned int i, unsigned int j);
C operator()(unsigned int i, unsigned int j) const
{return tab_[i-j+nbcol_-1];}
void transpose() ;
size_type nbrow() const {return nbrow_;}
size_type nbcol() const {return nbcol_;}
};
4 A container for Hankel matrices
A matrix H = (hi,j) is a Hankel matrix if its
entries hi,j depends only on i+j, that is, if hi+1,j+1 = hi,j for all pairs (i,j) for which the entries are defined. To this kind
of matrix, we can associate a Hankel operator which is defined as the
projection of the multiplication by a fixed Laurent polynomial (a Laurent
polynomial is a polynomial of K[x,x-1], that is a polynomial in both
the variables x and x-1). So we can compute the product of a
n× n Hankel matrix by a vector as a subvector of the product of a
fixed polynomial h(x) by a polynomial in x-1. This can also be done
within O(nlog(n)) ops (or O(nlog(n)log2(n) according to the
field), using for instance FFT. See (X) and [BP94].
template< class C>
struct hankel : public array1d< C>
Containers for Hankel matrices, as a subclass of array1d.
Defined by its number of row nbrow_, its number of columns nbcol_.
The elements of the first row and last column are stored, starting from
(0,0) up to (nbrow_,nbcol_).
The product of a Hankel matrix by a vector is performed by FFT
computations (X). .
{
unsigned int nbrow_, nbcol_;
typedef C coeff_t;
typedef array1d row_t;
typedef array1d col_t;
hankel(): array1d(){}
hankel(const hankel & v): array1d(v),
nbrow_(v.nbrow_), nbcol_(v.nbcol_){}
hankel(unsigned int i, unsigned int j):
array1d(i+j-1),nbrow_(i),nbcol_(j) {}
hankel(unsigned int i, unsigned int j, C* t):
array1d(i+j-1,t),nbrow_(i),nbcol_(j){}
void reserve(unsigned int i, unsigned int j)
{
this->array1d::reserve(i+j-1); nbrow_ = i; nbcol_=j;
}
C operator()(unsigned int i, unsigned int j) const
{return tab_[i+j];}
size_type nbrow() const {return nbrow_;}
size_type nbcol() const {return nbcol_;}
void transpose() { swap(nbcol_,nbrow_);};
5 How to use them
#include "linalg.H"
#include "linalg/toeplitz.H"
#include "linalg/hankel.H"
typedef VectStd<array1d<double> > VECT;
typedef hankel<double> rep; //Definition of the container.
typedef MatStruct<rep> MAT; //Definition of Hankel matrices.
int main(int argc, char** argv) {
double v[]={1,2,1,1,3,4,2,2,5,5,5,5,-2, 1,2,3,1, -2.8,-2.4,1,.2,5.8};
MAT A(4,4,v), B(A);
B.transpose();
cout<<"A := "<<A<<", B = "<<B<<endl;
| A := |
é
ê
ê
ê
ë |
| 1 |
2 |
1 |
1 |
| 2 |
1 |
1 |
3 |
| 1 |
1 |
3 |
4 |
| 1 |
3 |
4 |
2
|
|
ù
ú
ú
ú
û |
, B = |
é
ê
ê
ê
ë |
| 1 |
2 |
1 |
1 |
| 2 |
1 |
1 |
3 |
| 1 |
1 |
3 |
4 |
| 1 |
3 |
4 |
2
|
|
ù
ú
ú
ú
û |
MAT C;
C= A+B; cout<<C<<endl;
é
ê
ê
ê
ë |
| 2 |
4 |
2 |
2 |
| 4 |
2 |
2 |
6 |
| 2 |
2 |
6 |
8 |
| 2 |
6 |
8 |
4
|
|
ù
ú
ú
ú
û |
C= A-B; cout<<C<<endl;
é
ê
ê
ê
ë |
| 0 |
0 |
0 |
0 |
| 0 |
0 |
0 |
0 |
| 0 |
0 |
0 |
0 |
| 0 |
0 |
0 |
0
|
|
ù
ú
ú
ú
û |
C= A*2; cout<<C<<endl;
é
ê
ê
ê
ë |
| 2 |
4 |
2 |
2 |
| 4 |
2 |
2 |
6 |
| 2 |
2 |
6 |
8 |
| 2 |
6 |
8 |
4
|
|
ù
ú
ú
ú
û |
VECT V(v+1,v+5); cout<<V<<endl;
[2, 1, 1, 3]
//Multiplication using FFT.
VECT W(A*V); cout<<"W(A*V) = "<<W<<endl;
W(A*V) = [8, 15, 18, 15]
//Submatrices
W = Col(A,1); cout<<W<<endl;
[2, 1, 1, 3]
W = Row(A,2); cout<<W<<endl;
[1, 1, 3, 4]
B = A(Range(1,3),Range(1,3)); cout<<"A(1..3,1..3) = "<<B <<endl;
B = Row(A,Range(1,2)); cout<<B<<endl;
B = Col(A,Range(1,2)); cout<<B<<endl;
};