1 Introduction
Dense matrices are two-dimensional arrays in which
all the coefficients are stored (even if they are zero).
They provide classical arithmetic operations such
as addition, subtraction, multiplication by scalars, vectors, matrices
...
2 Implementation
template< class R>
struct MatDense
Interface for dense matrices.
The class R is the container corresponding to the internal
representation. The data are owned by the interface.
The inplace arithmetic operations are members of this class.
The external arithmetic operators +,-,*,/ are defined through
template expression. The differents operations of this class are
implemented in the module MATRIX (X).
{
R rep;
typedef typename R::size_type size_type;
typedef typename R::value_type value_type;
typedef typename R::row_iterator row_iterator;
typedef typename R::col_iterator col_iterator;
typedef VectStd > row_type;
typedef VectStd > col_type;
typedef R rep_type;
typedef MatDense self_t;
MatDense():rep() {}
MatDense(const R & r): rep(r) {}
MatDense(size_type m, size_type n): rep(m,n){}
MatDense(size_type m, size_type n, const value_type & c);
MatDense(size_type m, size_type n, value_type * t, ByRow);
MatDense(size_type m, size_type n, value_type * t, ByCol);
MatDense(size_type m, size_type n, const char * str, ByRow);
MatDense(size_type m, size_type n, const char * str, ByCol);
// Copy constructors.
MatDense(const self_t & M): rep(M.rep) {}
template
MatDense(const VAL & M) {assign(*this,M.rep);}
//Assignement operators
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;}
//Number of rows and columns
size_type nbrow() const {return rep.nbrow();}
size_type nbcol() const {return rep.nbcol();}
value_type & operator() (size_type i, size_type j) {return rep(i,j);}
value_type operator() (size_type i, size_type j) const {return rep(i,j);}
VAL > operator()(Range I, Range J);
//Inplace arithmetic operations
self_t & operator+=(const self_t & a);
self_t & operator-=(const self_t & a);
self_t & operator*=(const value_type & c);
self_t & operator/=(const value_type & c);
self_t & operator*=(const self_t& a);
// Rows and columns operations.
self_t & swapRow(size_type i, size_type j);
self_t & swapCol(size_type, size_type);
self_t & permRow(size_type *pivot);
self_t & permCol(size_type *pivot);
self_t & addRow(size_type i,size_type j, const value_type & c = 1);
self_t & addCol(size_type i,size_type j, const value_type & c = 1);
self_t & multRow(size_type i, const value_type & c = 1);
self_t & multCol(size_type i, const value_type & c = 1);
// Inplace transposition.
void transpose();
};
The value_type is the type of the coefficients. The size_type is the
type of the indices used to access the elements. The row_iterator is the
type of the forward iterators on the elements of a row. The column_iterator
is the type of the forward iterators on the elements of a column.The
row_type is the type of the row views. The col_type is the type of the
column views. The internal representation type R is rep_type.
template < class R>
inline VAL< OP< '+',MatDense< R> ,MatDense< R> > >
operator+(const MatDense< R> & a, const MatDense< R> & b)
Addition operator.
template < class R>
inline VAL< OP< '-',MatDense< R> ,MatDense< R> > >
operator-(const MatDense< R> & a, const MatDense< R> & b)
Substraction operator.
template < class R>
inline VAL< OP< '-',MatDense< R> ,MatDense< R> > >
operator-(const MatDense< R> & a)
Unary substraction operator.
template < class R> inline
VAL< OP< '*',MatDense< R> ,MatDense< R> > >
operator*(const MatDense< R> & a, const MatDense< R> & b)
Matrix multiplication operator.
template < class R,class S> inline
VAL< OP< '*',MatDense< R> ,VectStd< S> > >
operator*(const MatDense< R> & a, const VectStd< S> & b)
Matrix-vector multiplication operator.
template < class R> inline
VAL< OP< '.',MatDense< R> ,typename MatDense< R> ::value_type> >
operator*(const MatDense< R> & a, const typename MatDense< R> ::value_type & c)
Right scalar multiplication operator.
template < class R> inline
VAL< OP< '.',MatDense< R> ,typename MatDense< R> ::value_type> >
operator*(const typename MatDense< R> ::value_type & c, const MatDense< R> & a)
Left scalar multiplication operator.
template < class R> inline
VAL< OP< '/',MatDense< R> ,MatDense< R> > >
operator/(const MatDense< R> & a, const MatDense< R> & b)
Scalar division operator (on the right).
template < class R> inline
VAL< OP< '^',MatDense< R> ,int> > operator^(const MatDense< R> & a, int i)
Power operator.
template < class R> inline
typename MatDense< R> ::row_type Row(const MatDense< R> & a, int i)
Row of index i of a matrix (local view).
template < class R> inline
VAL< OP< 's',MatDense< R> ,Range2d> > Row(const MatDense< R> & a, const Range & r)
Submatrix corresponding to the rows of indices in the range r.
template < class R> inline
typename MatDense< R> ::col_type Col(const MatDense< R> & a, int i)
Column of index i of a matrix (local view).
template < class R> inline
VAL< OP< 's',MatDense< R> ,Range2d> > Col(const MatDense< R> & a, const Range & r)
Submatrix corresponding to the columns of indices in the range r.
template< class R>
ostream & operator< < (ostream & os, const MatDense< R> & p)
Output operator for standard matrices.
template< class R> inline
istream & operator> > (istream & is, MatDense< R> & M)
Input operator for standard matrices. The input format is of the form
m n c00 c01 ... where m is the number of rows, n is the number of
columns and the coefficients are read by rows.
template< class R>
MatDense< R> ::MatDense(size_type m, size_type n, value_type* t, ByRow):rep(m,n)
Construct a m by n matrix from an array of
coefficients. These coefficients are read
row by row.
template< class R>
MatDense< R> ::MatDense(size_type m, size_type n, value_type* t, ByCol):rep(m,n)
Construct a m by n matrix from an array of
coefficients. These coefficients are read
column by column.
template< class R>
MatDense< R> ::MatDense(size_type m, size_type n, const value_type & c): rep(m,n)
Construct a m by n matrix with entries all equal to c.
template< class R>
MatDense< R> ::MatDense(size_type m, size_type n, const char* str, ByRow):
rep(m,n)
Construct a m by n matrix from a string containing the
coefficients of the matrix. These coefficients are read
row by row.
template< class R>
MatDense< R> ::MatDense(size_type m, size_type n, const char* str, ByCol):
rep(m,n)
Same as above but the coefficients are read column by column.
template< class R> inline
VAL< OP< 's',MatDense< R> ,Range2d> > MatDense< R> ::operator()(Range I, Range J)
Submatrix from row and column ranges.
template< class R> inline
MatDense< R> & MatDense< R> ::operator+=(const MatDense< R> & a)
Inplace addition.
template< class R> inline
MatDense< R> & MatDense< R> ::operator-=(const MatDense< R> & a)
Inpace substraction.
template< class R>
MatDense< R> &
MatDense< R> ::operator*=(const typename MatDense< R> ::value_type & c)
Inplace multiplication by a scalar.
template< class R>
MatDense< R> &
MatDense< R> ::operator/=(const typename MatDense< R> ::value_type & c)
Inplace division by a scalar.
template< class R> inline
MatDense< R> & MatDense< R> ::operator*=(const MatDense< R> & a)
Inplace multiplication by a matrix.
template< class R>
MatDense< R> & MatDense< R> ::swapRow(const typename MatDense< R> ::size_type i,
const typename MatDense< R> ::size_type j)
Swap rows i and j.
template< class R>
MatDense< R> & MatDense< R> ::swapCol(const typename MatDense< R> ::size_type i,
const typename MatDense< R> ::size_type j)
Swap columns i and j.
template< class R>
MatDense< R> & MatDense< R> ::addRow(const typename MatDense< R> ::size_type i,
const typename MatDense< R> ::size_type j,
const typename MatDense< R> ::value_type & c)
Add c times row j into row i.
template< class R>
MatDense< R> & MatDense< R> ::addCol(const typename MatDense< R> ::size_type i,
const typename MatDense< R> ::size_type j,
const typename MatDense< R> ::value_type & c)
Add c times column j into column i.
template< class R>
MatDense< R> & MatDense< R> ::multRow(const typename MatDense< R> ::size_type i,
const typename MatDense< R> ::value_type & c)
Multiply by c the row i.
template< class R>
MatDense< R> & MatDense< R> ::multCol(const typename MatDense< R> ::size_type i,
const typename MatDense< R> ::value_type & c)
Multiply by c the column i.
template < class T> inline
void MatDense< T> ::transpose(void)
Inplace transposition.
3 Containers for dense matrices
The container R should provide the following definitions
or methods:
typedef typename R::value_type;
typedef typename R::size_type;
typedef typename R::iterator;
R(size_type i, size_type j, value_type* t);
R(size_type i, size_type j, value_type* t, char c);
R(size_type s, Alloc);
R(size_type s, C *t);
R(iterator b, iterator e);
void R::reserve(size_type m, size_type m);
value_type & operator() (size_type i, size_type j);
const value_type & operator() (size_type i, size_type j) const;
size_type nbrow() const;
size_type nbcol() const;
Possible containers are array2d (X) and
lapack (X).
4 The standard container
The default container is defined here
template < class C>
struct array2d : public array1d< C>
Class of array2d derived from array1d< C> , including the
number of rows nbrow_ and the number of columns nbcol_
of the matrix. The matrix is stored by column in the array array1d< C> .
{
typedef C* col_iterator;
typedef stride_iterator row_iterator;
size_type nbrow_, nbcol_;
array2d(): array1d(){}
array2d(const array2d & v):
array1d(v), nbrow_(v.nbrow_), nbcol_(v.nbcol_) {}
array2d(size_type i,size_type j=1):
array1d(i*j), nbrow_(i),nbcol_(j) {}
array2d(size_type i, size_type j, C* t):
array1d(i*j,t), nbrow_(i), nbcol_(j) {}
array2d(size_type i, size_type j, C* t, char c);
array2d(const char * nm);
void reserve(size_type i,size_type j);
array2d & operator=(const array2d & v);
C & operator()(size_type i, size_type j) {return tab_[i+nbrow_*j];}
const C & operator()(size_type i, size_type j) const {return tab_[i+nbrow_*j];}
size_type nbrow() const {return nbrow_;}
size_type nbcol() const {return nbcol_;}
};
5 The container for the connection with Lapack
template < class C>
struct lapack: public array2d< C>
Container for the connection with the fortran library lapack.
Class derived from array2d< C> , with a new lda_ field.
It indicates the leading dimension array, that is the number
of rows in the internal allocated array. It can differ from the number
of rows returned by nbrow() which is smaller than lda_. This field
is used only in an internal way.
{
unsigned int lda_;
lapack(): array2d() {}
lapack(const lapack & v): array2d(v), lda_(v.lda_){}
lapack(unsigned int i, unsigned int j=1): array2d(i,j), lda_(i){}
lapack(unsigned int i, unsigned int j, C* t): array2d(i,j,t), lda_(i){}
lapack(unsigned int i, unsigned int j, C* t, char c):
array2d(i,j,t,c), lda_(i) {}
lapack(unsigned int i, unsigned int j, unsigned int k, C* t):
array2d(i,j,t), lda_(i) {nbrow_=k;}
lapack(const char * nm):array2d(nm) {lda_ = nbrow_;}
lapack & operator=(const lapack & v);
void reserve(unsigned int i,unsigned int j);
C & operator()(unsigned int i, unsigned int j)
{return tab_[i+lda_*j];}
const C & operator()(unsigned int i, unsigned int j) const
{return tab_[i+lda_*j];}
};
integer LINALG::dcmp_qr(lapack< double> & A)
Compute the QR decomposition of the matrix A, without pivoting.
The operation is done inplace. The output matrix is upper-triangular.
integer LINALG::dcmp_qr(lapack< double> & A, double * tau)
Compute the QR decomposition of the matrix A. The operation is
done inplace, without pivoting. The matrix Q is
represented by the lower part of A and the coefficients in tau.
Use the lapack routine dgeqrf.
integer LINALG::dcmp_qr (lapack< double> & A,
unsigned int * pivot, double * tau)
Compute the QR decomposition of the matrix A. The operation is
done inplace. The pivots are stored in pivot and the matrix Q is
represented by the lower part of A and the coefficients in tau.
Use the lapack routine dgeqpf.
integer LINALG::dcmp_lu(lapack< double> & A)
Compute the LU decomposition of the matrix A. The pivots on the
columns are stored in pivot. Use the lapack routine dgetrf.
integer LINALG::dcmp_lu(lapack< double> & A, lapack< double> & U)
Compute the LU decomposition of the matrix A. The pivots on the
columns are stored in pivot. Use the lapack routine dgetrf.
integer LINALG::dcmp_lu(lapack< double> & A, integer *pivot)
Compute the LU decomposition of the matrix A. The pivots on the
columns are stored in pivot. Use the lapack routine dgetrf.
template < class C>
C LINALG::triang(Gauss m, lapack< C> & A)
Compute a triangulation of A, based on the LU decomposition of the matrix
A. Return the determinant, which is the product of terms on the diagonal.
Triang,Gauss Lapack
void LINALG::solve(LU m, array1d< double> & x,
const lapack< double> & A, const array1d< double> & b)
Solve the system A x = b by LU decomposition of the matrix A.
Use the lapack routine dgesv.
void LINALG::solve(LU m, lapack< double> & X,
const lapack< double> & A)
Solve the system A X = B where B is a matrix by LU decomposition
of the matrix A. Use the lapack routine dgesv.
void LINALG::solve_ls(lapack< double> & A, array1d< double> & b)
Solve the system A x = b by a Least Square method.
Put the result in b.
Use the lapack routine dgelsx.
void LINALG::solve_ls(array1d< double> & x,
const lapack< double> & A, const array1d< double> & b)
Solve the system A x = b by a Least Square method.
A copy of A is made and a copy of b is put in x.
Then the previous function
solve_ls(lapack< double> & A, array1d< double> & b) is used.
void LINALG::svd(array1d< double> & x,const lapack< double> & A)
Compute the Singular values of the matrix A, calling the lapack
subroutine dgesvd. The result, put in x, is the vector of singular values
in decreasing order. The matrix A is copied.
void LINALG::svd(array1d< double> & x,
const lapack< double> & A, lapack< double> & U)
Compute the Singular values of the matrix, calling lapack subroutine
dgesvd. The result, put in x, is the vector of singular values
in decreasing order. The matrix A is copied.
The singular left vectors are stored in U.
void LINALG::svd(array1d< double> & S,
const lapack< double> & A,
lapack< double> & U, lapack< double> & V)
Compute the Singular values of the matrix, calling the lapack
subroutine dgesvd. The result (put in S), is the vector of singular
values in decreasing order. The matrix A is copied in a temporary
matrix.
The matrices U and V are orthogonal matrices, which satisfy
A=U S Vt.
The correct space for the matrix U and V should be allocated
before using this procedure.
void LINALG::svdx(array1d< double> & S,
const lapack< double> & A,
lapack< double> & U, lapack< double> & V)
The same but we have A=U Vt.
void LINALG::eignval(array1d< complex< double> > & V, const lapack< double> & A)
Compute the eigenvalues of A, put them in the complex
vector V. Use the lapack routine dgeev.
The matrix A is copied.
void LINALG::eignval(
array1d< double> & V, const lapack< double> & A, Real)
Compute the real eigenvalues of A and put them in the real
vector V. Use the lapack routine dgeev.
The matrix A is copied.
void LINALG::eignval(array1d< complex< double> > & V, const lapack< double> & A,
array1d< double> & error)
Compute the eigenvalues of A, put them in the complex
vector V and store in error[i] a bound on the error on the
eigenvalue V[i]. Use the lapack routine dgeevx.
The matrix A is copied.
void LINALG::eignval(
array1d< complex< double> > & V,
const lapack< double> & A, const lapack< double> & B)
Compute the eigenvalues of the pencil A - z B and put it in the complex
vector V. Use the lapack routine dgegv.
The matrices A and B are copied.
void LINALG::eignval(
array1d< double> & V,
const lapack< double> & A, const lapack< double> & B, Real real)
Compute the real eigenvalues of the pencil A - z B and put it in
the complex vector V. Use the lapack routine dgegv.
The matrices A and B are copied.
void LINALG::leigvct(lapack< complex< double> > & E,
const lapack< double> & A)
Compute the left eigenvectors of A and put it in the complex
matrix E. Use the lapack routine dgeev. The matrix A
is copied.
void LINALG::leigvct(lapack< complex< double> > & E,
array1d< complex< double> > & V, const lapack< double> & A)
Compute the left eigenvectors and eigenvalues of A.
The matrix of eigenvector is put in the complex matrix E and
the eigenvalues are stored in V. The matrix A
is copied.
void LINALG::reigvct(lapack< complex< double> > & E,
const lapack< double> & A)
Compute the right eigenvectors of the matrix A.
The complex matrix of eigenvector is put in E.
Use the lapack routine dgeev. The matrix A
is copied.
void LINALG::reigvct(lapack< double > & E,
const lapack< double> & A, const Real &)
Compute the real right eigenvectors of the matrix A.
The complex matrix of eigenvector is put in E.
Use the lapack routine dgeev. The matrix A
is copied.
void LINALG::reigvct(lapack< complex< double> > & E,
array1d< complex< double> > & V, const lapack< double> & A)
Compute the left eigenvectors and eigenvalues of A.
The matrix of eigenvector is put in the complex matrix E and
the eigenvalues are stored in V. The lapack routine is
dgeev. The matrix A is copied.
void LINALG::reigvct(lapack< double> & E,
const lapack< double> & A, const lapack< double> & B,
const Real & real,
unsigned int M)
Compute the real eigenvectors of the pencil A - z B and keep only
their m first coordinates. The matrices A and B are copied.
void LINALG::reigvct(lapack< double> & E,
const lapack< double> & A, const lapack< double> & B,
const Real &, array1d< double> & V)
Compute the real eigenvectors of the pencil A - z B. Put the
eigenvalues in V. The matrices A and B are copied.
template< class C>
typename lapack< C> ::size_type LINALG::rank(const lapack< C> & M)
Compute the rank of M, using a QR-decomposition (use the routine
dgeqpf). The matrix M is copied.
template< class C>
typename lapack< C> ::size_type LINALG::rank(const lapack< C> & M,
typename lapack< C> ::size_type * permut_row,
typename lapack< C> ::size_type * permut_col)
Compute the rank of M, using a QR-decomposition (use the routine
dgeqpf. The matrix M is copied. The permutation of rows and
columns needed to get a leading square matrix of maximal rank
are stored in permut_row and permut_col.
integer LINALG::inv_tr(lapack< double> & A, char s, char u)
Inplace inversion of triangular matrices. The value of s should be
'L' for lower or 'U' for upper triangular matrix. The value of u
should be 'U' for unitary (implicitly 1 on the diagonal) or 'N' for
the default. Return an integer which is 0 if the computation
runs correctly.
6 How to use them
Here is an example of dense matrices, with the container lapack.
#include "linalg.H" // General header file for linear algebra.
#include "lapack.H" // We will use lapack containers.
#include "maple.H" // For output in maple format.
typedef VectStd<array1d<double> > VECT;
// Vector definition based on array1d container, with double coefficients.
typedef MatDense<lapack<double> > MAT;
// Marix definition based on lapack container, with double coefficients.
// This will allow us to use the routines of the lapack library.
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,ByRow()),
B(4,4,"1 2 1 1 3 4 2 2 5 5 5 5 -2 1 2 3 1 .2 5.8", ByCol());
cout<<"A="<<A<<", B="<<B<<endl;
| A= |
é
ê
ê
ê
ë |
| 1 |
2 |
1 |
1 |
| 3 |
4 |
2 |
2 |
| 5 |
5 |
5 |
5 |
| -2 |
1 |
2 |
3
|
|
ù
ú
ú
ú
û |
, B= |
é
ê
ê
ê
ë |
| 1 |
3 |
5 |
-2 |
| 2 |
4 |
5 |
1 |
| 1 |
2 |
5 |
2 |
| 1 |
2 |
5 |
3
|
|
ù
ú
ú
ú
û |
// Arithmetic operations.
MAT C;
C= A+B; cout<<C<<endl;
é
ê
ê
ê
ë |
| 2 |
5 |
6 |
-1 |
| 5 |
8 |
7 |
3 |
| 6 |
7 |
10 |
7 |
| -1 |
3 |
7 |
6
|
|
ù
ú
ú
ú
û |
C= A-B; cout<<C<<endl;
é
ê
ê
ê
ë |
| 0 |
-1 |
-4 |
3 |
| 1 |
0 |
-3 |
1 |
| 4 |
3 |
0 |
3 |
| -3 |
-1 |
-3 |
0
|
|
ù
ú
ú
ú
û |
C= A*B; cout<<C<<endl;
é
ê
ê
ê
ë |
| 7 |
15 |
25 |
5 |
| 15 |
33 |
55 |
8 |
| 25 |
55 |
100 |
20 |
| 5 |
8 |
20 |
18
|
|
ù
ú
ú
ú
û |
C= A*B*A; cout<<C<<endl;
é
ê
ê
ê
ë |
| 167 |
204 |
172 |
177 |
| 373 |
445 |
372 |
380 |
| 650 |
790 |
675 |
695 |
| 93 |
160 |
157 |
175
|
|
ù
ú
ú
ú
û |
C= 2*A; cout<<C<<endl;
é
ê
ê
ê
ë |
| 2 |
4 |
2 |
2 |
| 6 |
8 |
4 |
4 |
| 10 |
10 |
10 |
10 |
| -4 |
2 |
4 |
6
|
|
ù
ú
ú
ú
û |
// Construction of a Vector from to iterators.
VECT V(v+1,v+5); cout<<V<<endl;
[2, 1, 1, 3]
VECT W(A*V); cout<<W<<endl;
[8, 18, 35, 8]
//Extractions of sub matrices of vectors:
W = Col(A,1); cout<<"Col(A,1) = "<<W<<endl;
Col(A,1) = [2, 4, 5, 1]
W = Row(A,2); cout<<"Row(A,2) = "<<W<<endl;
Row(A,2) = [5, 5, 5, 5]
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;
//Addition of the row 2 times 1.5 into row 1
A.addRow(1,2,1.5); cout<<A<<endl;
é
ê
ê
ê
ë |
| 1 |
2 |
1 |
1 |
| 10.5 |
11.5 |
9.5 |
9.5 |
| 5 |
5 |
5 |
5 |
| -2 |
1 |
2 |
3
|
|
ù
ú
ú
ú
û |
A.addCol(1,2,-0.5); cout<<A<<endl;
é
ê
ê
ê
ë |
| 1 |
1.5 |
1 |
1 |
| 10.5 |
6.75 |
9.5 |
9.5 |
| 5 |
2.5 |
5 |
5 |
| -2 |
0 |
2 |
3
|
|
ù
ú
ú
ú
û |
// Output to the maple format.
cmpl <<"A := "<<A<<endl;
//A := linalg[matrix](4,4, [[1,1.5,1,1], [10.5,6.75,9.5,9.5], [5,2.5,5,5], [-2,0,2,3]]);
};