| The linear algebra package |
1 Fast Fourier Transform
template< class T> inline void Root_Unity(complex< T> &c, int n)
Root of Unity. Replace c by the primitive root of order n.
inline unsigned int BitReverse(int k, int n)
Functions which permute an array
by Bit Reversal method.
template < class T>
inline void Permut(T *Tab, int N)
Shuflle on N entries.
template < class C, class U>
void FFT(int n, const complex< U> & w, complex< C> *W, int inverse)
Implementation of FFT computation of a complex array.
It is an "in-place" transformation, based on the Butterfly algorithm.
A complex array W, of size n (a power of 2), is replaced by
the FFT of its values. The complex number w is a primitive
root of unity of order n.
2 Triangular decomposition
template< class A, class R>
VAL< R*> Triang(const R & mat)
Triangulation of the matrix m. The first argument of
of template specifies the type of triangulation which is applied.
It should called as follows Triang< Method> (mat);.
The function LINALG::triang(Method(), ...) is called on a copy of mat.
The different methods available are Gauss, QR, Bareiss,... Example:
<B><PRE>
T = Triang$<$ Gauss$>$ (mat)
</PRE></B>
template < class M> inline
VAL< M*> USchur(M & m, int n)
Compute the upper schur complement of size n of the matrix mt, given
by the formula A -B D-1 C where A is the first diagonal submatrix
block of size n.
USchur
3 Solving linear systems
template < class M, class R> inline
VAL< VectStd< R> *> Solve(M & A, const VectStd< R> & b, LU mth)
Resolution by LU decomposition. Calling the function LINALG::solve_lu.
Solve, LU
template < class M, class R> inline
VAL< R*> Solve(M & A, R & B, LU mth)
Resolution of A X = B by LU decomposition.
Calling the function LINALG::solve_lu(X.rep,A.rep),
where X is a copy of B.
Solve, LU
template < class M, class V> inline
VAL< V*> Solve(M & A, V & b, LS mth)
Resolution by LeastSquare method. Calling the function
LINALG::solve_ls(X.rep,A.rep,b.rep) where X is a copy of A.
Solve, LS
4 Determinants
Different methods for computing the determinant of a matrix are discribed
here: Gauss, Bareiss or Berkowitz methods.
- Gauss method uses classical pivoting techniques (ie. LU decomposition)
to compute the determinant.
- Bareiss method applies the following pivoting scheme:
|
ai,j(k) = |
| ai,j(k-1) ak,k(k-1) - ai,k(k-1)
ak,j(k-1) |
|
| ak-1,k-1(k-2) |
|
(1) |
where A(k)=(ai,j(k)) is the transformed matrix at step k and
A(0) is the initial matrix. At the end of this process, after some
permutation of the rows, the matrix is upper triangular. If A is a square
matrix, the last entry An,n(n) on the diagonal at the end of the
process is the determinant of the matrix A. Moreover, if
A has its coefficients in a ring, it is the same for the matrices
A(k), since the division in (X) is exact.
- Berkowitz method is based on the following recursion formula, based on
Cayley-Hamilton identity:
| d0 := Trace(A); A0 := A; Ak+1 := Ak*A- |
|
Trace(Ak) A
|
- A mixed strategy, which expands the determinant with respect
to the first row and Bareiss method is applied to the subdeterminants.
template < class R>
typename R::value_type Det(R & M)
Compute the determinant of M, using Bareiss method.
Example:
<B><PRE>
typedef MatDense$<$ array2d$<$ double$>$ $>$ MAT
MAT M(3,3,"1 2 3 4 5 6 7 8 9", ByRow())
double det
det=Det(M)
</PRE></B>
template < class R>
typename R::value_type Det(const Mixed & mth, R & M)
Compute the determinant of M by mixing two methods:
A pseudo-triangulation is first computed.
If the size of M is less or equal to 4, the determinant is computed by
expansion along the columns. Else a development is done along the first
row with triangulation by Bareiss method of the submatrix corresponding to each column that is removed. Example:
<B><PRE>
typedef MatDense$<$ array2d$<$ double$>$ $>$ MAT
MAT M(3,3, "1 2 3 4 5 6 7 8 9",ByRow())
double det=Det(Mixed(),M)
</PRE></B>
template< class R>
typename R::value_type Det(const Berko & mth, const R & M)
Compute the determinant of M by Berkowitz method.
Example:
<B><PRE>
typedef MatDense$<$ array2d$<$ double$>$ $>$ MAT
MAT M(3,3,"1 2 3 4 5 6 7 8 9")
double det=Det(Berko(), M)
</PRE></B>
5 Rank
template< class R> inline
typename R::size_type Rank(const R & M)
Compute the rank of M, using LINALG::rank.
template< class R,class V>
typename R::size_type Rank(const R & M,
V & I, V & J)
Compute the rank of M, using LINALG::rank. The columns (resp. row)
indices of a permutation which put the non-zero minor in the upper corner
are stored in permut_col and permut_row. They have to be initialized
before, to a size larger than min(M.nbcol(), M.nbrow()).
template< class R,class I>
VAL< R*> Extract(R & M, const I & ipvt, const I & jpvt,
const typename R::size_type & rank)
Extraction of a square submatrix of size rank from the matrix M.
This is performed by keeping the columns whose index is one of the
first rank elements of jpvt,
and the rows whose index is one of the first rank elements of ipvt,
in the same orders.
6 Eigenvectors and eigenvalues
Classical eigenproblem.
template < class M> inline
VAL< VectStd< array1d< AlgClos< typename M::value_type> ::TYPE > > *>
Eigenval(M & A)
Compute the eigenvalues of A.
template < class M> inline
VAL< VectStd< array1d< typename M::value_type> > *> Eigenval(M & A, Real)
Compute the real eigenvalues of A.
template < class M> inline
VAL< VectStd< array1d< AlgClos< typename M::value_type> ::TYPE> > *>
Eigenval(M & A, VectStd< array1d< typename M::value_type> > & Er)
Compute the eigenvalues of A and the error Er on these eigenvalues.
template < class M> inline
VAL< MatDense< lapack< AlgClos< typename M::value_type> ::TYPE> > *>
Eigenvct(M & A)
Compute the eigenvectors of A.
template < class M> inline
VAL< MatDense< lapack< typename M::value_type> > *> Eigenvct(M & A, Real)
Compute the real eigenvectors of A.
Generalized eigenproblem.
template < class M> inline
VAL< VectStd< array1d< AlgClos< typename M::value_type> ::TYPE> > *>
Eigenval(const M & A, const M & B)
Compute the eigenvalues of (A,B).
template < class M> inline
VAL< VectStd< array1d< typename M::value_type> > *>
Eigenval(const M & A, const M & B, Real real)
Compute the real eigenvalues of (A,B).
template < class M,class V> inline
VAL< MatDense< lapack< typename M::value_type> > *>
Eigenvct(MatDense< M> & A,
MatDense< M> & B, Real r,
VectStd< V> & v)
Compute the real eigenvectors of (A,B) and the eigenvalues V.
template < class M> inline
VAL< MatDense< lapack< typename M::value_type> > *>
Eigenvct(MatDense< M> & A,
MatDense< M> & B, Real real,
const typename M::size_type & m)
Compute the real eigenvectors of (A,B).
The first m coordinates of these eigenvectors are returned.
7 Exemple
#include "linalg.H" // General header file for linear algebra.
#include "lapack.H" // We will use lapack containers.
typedef VectStd<array1d<double> > VECT;
// Vector definition.
typedef VectStd<array1d<complex<double> > > VECTC;
// The same with complex coefficients.
typedef MatDense<lapack<double> > MAT;
// Marix definition.
typedef MatDense<lapack<complex<double> > > MATC;
// The same with complex coefficients.
int main(int argc, char** argv) {
double v[]={1,-3,1,2,3,4,2,2,4,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 |
-3 |
1 |
2 |
| 3 |
4 |
2 |
2 |
| 4 |
5 |
5 |
5 |
| 5 |
-2 |
1 |
2
|
|
ù
ú
ú
ú
û |
, B= |
é
ê
ê
ê
ë |
| 1 |
3 |
5 |
-2 |
| 2 |
4 |
5 |
1 |
| 1 |
2 |
5 |
2 |
| 1 |
2 |
5 |
3
|
|
ù
ú
ú
ú
û |
VECT V(v+1,v+5); cout<<V<<endl;
[-3, 1, 2, 3]
VECT W(A*V); cout<<W<<endl;
[2, 5, 18, -9]
W =Svd(A); cout<<"Svd(A) = "<<W<<endl;
Svd(A) = [11.4091, 6.03278, 2.53035, 0.189481]
cout<<"Rank(A) = "<<Rank(A)<<endl;
>> Rank
Rank(A) = 4
//Triangulation by Gauss method.
B = Triang<Gauss>(A); cout<<B<<endl;
é
ê
ê
ê
ë |
| 5 |
-2 |
1 |
2 |
| 0 |
6.6 |
4.2 |
3.4 |
| 0 |
0 |
2.45455 |
2.93939 |
| 0 |
0 |
0 |
0.407407
|
|
ù
ú
ú
ú
û |
//Determinant computed by Bareiss method.
cout<<"Det = "<<Det(A)<<endl;
Bareiss
Det = -33
//Eigenvalue and eigenvector computations.
VECTC L=Eigenval(A); cout<<L<<endl;
[(-1.21002,0), (6.26876,1.12049), (6.26876,-1.12049), (0.672512,0)]
MATC E=Eigenvct(A); cout<<E<<endl;
é
ê
ê
ê
ë |
| 0.377172 |
-0.132449+i*(0.0803325) |
-0.132449+i*(-0.0803325) |
-0.0578894 |
| -0.09358 |
0.457738+i*(-0.0040333) |
0.457738+i*(0.0040333) |
-0.105242 |
| 0.470337 |
0.855002 |
0.855002 |
0.820628 |
| -0.792318 |
-0.134821+i*(0.131372) |
-0.134821+i*(-0.131372) |
-0.558698
|
|
ù
ú
ú
ú
û |
};