Vectors and matrices#
Caution
If you are not familiar with C++ syntax, we recommend you to see C++ syntax before reading this page.
The Vector class#
The purpose of the Vector class is mainly to deal with complex or real vector. In
particular, this class is used in the definition of the user functions (see Functions and kernels).
It is a templated class, inheriting from std::Vector<T>, mainly used as a real or complex vector, through its aliases Reals and Complexes:
Reals u; // u=[0.]
Reals v(3); // v=[0. 0. 0.]
Reals w(3, 2.5); // w=[2.5 2.5 2.5]
Reals e{1.,2.,3.}; // from initialization list, e=[1. 2. 3.]
Reals f={1.,2.,3.}; // alternate construction
Complexes cu; // cu=[(0.,0.)]
Complexes cv(3); // cv=[(0.,0.) [(0.,0.) (0.,0.)]
Complex c(2.,1.); // the complex 2+i
Complexes cw(3, c); // cv=[(2.,1.) [(2.,1.) (2.,1.)]
Complex ce(e); // ce=[(1.,0.) (2.,0.) (3.,0.)], Reals to Complexes
Complex cf = f; // cf=[(1.,0.) (2.,0.) (3.,0.)], Reals to Complexes
Of course, you can define by yourself a vector of anything, such as Vector<Number> (that is different from Numbers), but we will focus in the following on vectors of reals or complexes.
It is also possible to deal with vector of vectors, with, for instance, its aliases RealVectors or ComplexVectors for instance:
Reals ones(3, 1); // ones=[1. 1. 1.]
RealVectors U(3, ones); // U=[[1. 1. 1.] [1. 1. 1.] [1. 1. 1.]]
RealVectors V={{11.,12.,13.}, {21.,22.,23.}}; // V=[[11.,12.,13.] [21.,22.,23.]]
To access to a vector component (both read and write access) use the operator () with index from 1 to the vector length:
Reals v(3); // v=[0. 0. 0.]
v(1)=1.;v(2)=2.;v(3)=3.; // v=[1. 2. 3.]
Complexes cv(3); // cv=[(0.,0.) [(0.,0.) (0.,0.)]
cv(2)=Complex(1,1); // cv=[(0.,0.) [(1.,1.) (0.,0.)]
Number s = v.size(); // vector size
Note that access using operator [] with index from 0 to the vector length -1, is also possible.
Advanced users can use member functions begin and end returning respectively iterators (or const iterators) to the beginning and the end of the vector.
It is also possible to extract some vector components in a new vector or to set some vector components by specifying a set of indices either given by lower and upper indices or given by a vector of indices:
Reals v(5); // v =[0. 0. 0. 0. 0.]
for (Number i=1; i<=5; i++)
v(i)=i*i; // v =[1. 5. 9. 16. 25.]
Reals w=v(3, 5); // w =[9. 16. 25.]
Vector<Number> is(3);
is(1)=1; is(2)=3; is(3)=5; // is=[1 3 5]
w=v(is); // w =[1. 9. 25.]
w=v.get(is); // w =[1. 9. 25.]
v.set(1, 3, w); // v =[1. 9. 25. 16. 25.]
Reals z(3, 0.); // w =[0. 0. 0.]
v.set(is, z); // v =[0. 9. 0. 16. 0.]
Standard algebraic operations ( +=, -=, *=, /=, +, -, *, /) are supported by the Vector. Some shortcuts are also possible, for instance a vector plus a scalar, a scalar plus a vector, … Here are a few examples:
Reals u(3, 1); // u=[1. 1. 1.]
Reals v(3); // v=[0. 0. 0.]
v=2.; // v=[2. 2. 2.]
Reals w=u+v; // v=[3. 3. 3.]
w=2.*v; // v=[4. 4. 4.]
w-=2.; // v=[2. 2. 2.]
Complex i(0,1); // complex number
Complexes cv(3,i); // cv=[(0.,1.) [(0.,1.) (0.,1.)]
cv=cv*i; // cv=[(-1.,0.) [(-1.,0.) (-1.,0.)]
cv/=2.; // cv=[(-0.5,0.) [(-0.5.,0.) (-0.5.,0.)]
For algebraic operations involving two vectors, the compatibility of the size of vectors is checked. All the algebraic operations involving a real vector (resp. a complex vector) and a real scalar (resp. a complex scalar) are supported. Be cautious, as an integer value is not always certainly cast to a real value, some operations may fail during the compiling process. For instance, the addition between a real vector and an integer does not work, cast explicitly to a real!
Reals u(3,1); // u=[1. 1. 1.]
Reals v(3); // v=[0. 0. 0.]
v=u+2; // DOES NOT WORK (2 is an int)
v=u+2.; // v=[3. 3. 3.]
Automatic cast from real vector to complex vector is supported. For instance, the following instructions are legal:
Complex i(0,1); // complex number i
Reals u(3,1); // u=[1. 1. 1.]
Complexes cv(3); // cv=[(0.,0.) [(0.,0.) (0.,0.)]
cv=u+2.*i; // cv=[(1.,2.) [(1.,2.) (1.,2.)]
cv*=3.; // cv=[(3.,6.) [(3.,6.) (3.,6.)]
Caution
Automatic cast is not supported for vector of vectors.
The class also provides some various functionalities:
Reals u(3, 1); // u=[1. 1. 1.]
u.norminfty(); // the sup norm of u
u.norm2squared(); // squared l2 norm of u
u.norm2(); // l2 norm of u
u.norm1(); // l1 norm of u
cout<<u; // output the vector u: [1 1 1]
u(1)=1.e-8; // u=[1.e-8 1. 1.]
Real eps=1.e-6;
u.roundToZero(eps); // u=[0. 1. 1.] replace by 0 all |ui| < eps
u.round(eps); // round to nearest eps : round(ui/eps)*eps
u.normalize(); // u/u.norm2(u)
u.normalize(0); // u/u.norminfty(u)
Complexes cv(3, i); // cv=[(0.,1.) [(0.,1.) (0.,1.)]
conj(u); // conjugate of cv
cv=cmplx(u); // transform a real vector in a complex one
u=real(cv); // take the real parts
u=imag(cv); // take the imaginary parts
Reals s = sqrt(u); // sqrt(ui), ui>=0 , works in complex
s = log(u); // log(ui), ui>=0 , works in complex
s = power(u,2.5); // ui^2.5, ui>=0 , works in complex
dotRC(u,s); // inner product returning a complex
dotC(u,cv); // hermitian product returning a complex
dot(u,s); // inner product returning a real or a complex
crossProduct(u,s); // cross product restricted to 3d vector
crossProduct2D(u,s); // cross product restricted to 2d vector
norm1(u); norm2(u), norminfty(u); // norm of vector
trivialNumbering(2.,5.); // return [2. 3. 4. 5.]
lineSpace(0.,1.,5); // return [0. 0.25 0.5 0.75 1.]
Contrary to the Point class, the Vector class offers no comparison function.
Note also that there is no direct link between these two classes except that a Point may be automatically constructed from a Vector:
Reals u(3, 1);
Point P=u;
The Matrix class#
The purpose of the Matrix class is mainly to deal with complex or real dense
matrices. In particular, this class is used in the definition of the user functions (see Functions and kernels). This template class inherits from std::vector<T>
and is compliant with the Vector class. Although, it can deal with matrices of anything, it is only fully functional for real or complex matrices.
Matrix<Real> and Matrix<Complex> are aliases by RealMatrix and ComplexMatrix:
RealMatrix rA; // an empty matrix
RealMatrix rB(3, 2); // a 3x2 zeros matrix
RealMatrix rC(3, 2, 1); // a 3x2 ones matrix
RealMatrix rE = {11.,12},{21.,22}; // a 2x2 matrix from explicit
ComplexMatrix cA(3, 2, i_); // a 3x2 i matrix
It is possible to construct diagonal matrix from a Vector or a matrix from a Vector of Vector, to load (and save) a matrix from a file and to construct particular matrices (zeroMatrix, onesMatrix, idMatrix, hilbertMatrix):
Reals u(3, 2.); // vector [2. 2. 2.]
RealMatrix rA(u); // 3x3 matrix with u as diagonal
Vector<Reals> us(2,u); // Vector of Vector
RealMatrix rU(us); // 2x3 matrix from Vector of Vector
RealMatrix rB("mat.dat"); // matrix loaded from "mat.dat" file
RealMatrix rO(3, zeroMatrix); // a 3x3 zeros matrix
RealMatrix r1(3, onesMatrix); // a 3x3 ones matrix
RealMatrix rI(3, idMatrix); // a 3x3 identity matrix
RealMatrix rH(3, hilbertMatrix); // the 3x3 Hilbert matrix
Tip
Matrix constructors work also with std::vector and std::vector<std::vector>.
Construction of complex matrix from real data are allowed (automatic cast). But the contrary is not.
There are some member functions to access to the matrix properties:
vsize()/numberOfRows(), hsize()/numberOfColumns(), dims() // returning a dimPair
isSymmetric(), isSkewSymmetric(), isSelfAdjoint(), isSkewAdjoint()
and some utilities to access to a coefficient, a row or a column or the diagonal of the matrix:
RealMatrix A(2, 2, 1); // a 2x2 ones matrix
A(1, 1)=2.; // change the coefficient A11
Reals r=A.row(1); // first row of A
r=A.column(2); // second column of A
r=A.diag(); // diagonal of A
A.column(1,r); // assign a vector to the first column
A.row(2,r); // assign a vector to the second row
A.diag(r); // assign a vector to the diagonal
All these functions support automatic cast from real to complex but not the contrary.
Advanced users can use member functions begin() and end() returning respectively iterators (or const iterators) to the beginning and the end of the Matrix.
The data values of the matrix are stored according to the C convention, i.e. row-wise.
There are also generalized access tools either to extract submatrix (get() or operator()) or to set submatrix of matrix (set()):
RealMatrix M(3, 3);
for (Number i=1; i<=3; i++)
for(Number j=1; j<=3; j++)
M(i,j )=i+j; // M=[2 3 4; 3 4 5; 4 5 6]
RealMatrix N= M.get(2, 3, 2, 3); // N=[4 5; 5 6]
// N=M(2, 3, 2, 3) or N=M.subMatrix(2, 3, 2, 3) give the same
Vector<Number> is(2);
is(1)=1; is(2)=3;
N= M.get(is, is); // N=[2 4; 4 6]
// N=M(is, is) gives the same
RealMatrix Z(2, 2, 0); // Z=[0 0; 0 0]
M.set(1, 2, 1, 2, Z); // M=[0 0 4; 0 0 5; 4 5 6]
RealMatrix U(2, 2, 1); // U=[1 1; 1 1]
M.set(is, is, Z); // M=[1 0 1; 0 0 5; 1 5 1]
RealMatrix L=M.lower(); // M=[1 0 0; 0 0 0; 1 5 1]
L=M.lower(2.); // M=[2 0 0; 0 2 0; 1 5 2]
RealMatrix U=M.upper(); // M=[1 0 1; 0 0 5; 0 0 1]
U=M.upper(2.); // M=[2 0 1; 0 2 0; 0 0 2]
Besides, the Matrix class proposes some transformations either as internal functions
or external functions:
RealMatrix A(2, 2, 1), B;
ComplexMatrix C(2, 2, i), D;
A.roundToZero(1.e-6); // replace Aij by 0 when |Aij| < 1e.-6
A.transpose(); // self transposition of A
B=transpose(A); // transposition of A, A not changed
C.adjoint(); // self transposition and conjugate C
D=adjoint(A); // transpose and conjugate, C not changed
B=diag(A); // from diagonal of A to a diagonal matrix
A=real(C); // real part of C
B=imag(C); // imaginary part of C
D=conj(C); // conjugate of C
D=cmplx(A); // forced casting from real to complex
Real n2=norm2(A); // Frobenius norm
Real ninf=norminfty(A); // infinite norm
Standard algebraic operations ( += , -= , *= , /= , + , - , * , / ) are supported by the Matrix class.
Some shortcuts are also possible, for instance a matrix plus a scalar, a scalar plus a matrix,
… Automatic cast from real to complex is supported. There is no comparison operator.
The Matrix proposes some solvers:
RealMatrix A, M;
Reals B;
RealMatrix Q , R, U, V, LU;
Complexes lambdas; // eigenvalues
Vector<Complexes> Es; // eigenvectors
Reals S; // singular values
std::vector<Dimen> P; // permutation vector
...
gaussSolver(A, B, piv, row); // solve AX=B using Gauss reduction
gaussMultipleSolver(A, B, nbrhs, piv, row); // solveAXs=Bxs using Gauss reduction
RealMatrix invA=inverse(A); // inverse of a square matrix
lu(A, LU); // LU factorization, A preserved
A=lu(A); // LU factorization, A not preserved
lu(A, LU, P); // LU factorization with pivoting
qr(A, Q, R); // QR factorization
eigs(A,lambdas,Es); // eigen elements AX=lambda X
eigs(A,lambdas,M, Es); // eigen elements AX=lambda MX (real!)
Number r = 10; // max number of singular values
Real eps=1e-6; // singular value threshold
svd(A, U, S, V, r, eps); // restricted SVD factorization A= U S V*
svd(A, U, S, V); // full SVD factorization A= U S V*
rsvd(A, U, S, V, r, eps); // restricted random SVD factorization
rsvd(A, U, S, V); // full random SVD factorization
r3svd(A, U, S, V); // Improved random SVD factorization
Warning
QR, EIGS and SVD are available only if Eigen library is set on.
Tip
It is also possible to deal with matrix of matrices, (aliases RealMatrices and ComplexMatrices) for instance:
RealMatrix ones(2, 2, 1); // a 2x2 ones matrix
RealMatrices A(2, 2, ones); // a 2x2 matrix of 2x2 ones matrix
but all operations are not supported for such matrices!
The SparseMatrix class#
The SparseMatrix class is only provided for user own business. It is not used by XLiFE++ stuff. It handles sparse matrix: to each stored coefficient is associated its location index pair \((i,j)\) (indices start from 1).
It extends automatically its storage, when new index pair is handled in writing mode. SparseMatrix<T> is a template class with T a real or a complex type.
Standard algebraic operations ( += , -= , *= , /= , + , - , * , / ) are supported by the SparseMatrix class.
SparseMatrix<Real> s0; // null matrix
SparseMatrix<Real> s42(4,2); // void 4x2 matrix
Vector<Real> d3(3,1);
SparseMatrix<Real> s3(d3); // constructor from diagonal
SparseMatrix<Real> sd5(5,_idMatrix); // Id constructor
s42(3,2)=32; s42(2,1)=21; // alllocate new coefficients
s42(5,1)=51; // dim increases!
theCout << s42; // output on stream
Vector<Real> r1=s42.getRow(1); // acces to row 1
Vector<Real> c2=s42.getCol(2); // acces to col 2
bool sym = s3.isSymmetric(); // check symmetry
s42+=sd5; s42-=sd5; // += -= operators
s42*=10; s42/=10; // *= /= operators
theCout << (s42*4.-2.*s42)/2.; // linear combination
theCout << s42.real() << s42.imag(); // real and imag part
s42.saveToFile("s42.dat"); // save to file, format i j aij
SparseMatrix<Real> s42f;
s42f.loadFromFile("s42.dat"); // load from file
theCout <<abs(s42); // sparse matrix |aij|
bool same =(s42==s42f); // comparison
Vector<Real> r5(5), v2(2);
theCout << s42*v2; // matrix x vector
theCout << v5*s42; // vector x matrix
theCout << transpose(s42); // transpose
theCout << adjoint(s42); // adjoint sparse matrix
theCout << conjugate(s42); // conjugate sparse matrix
theCout << dimsOf(s42); // access to dimension (pair)
Matrix<Real> m42=s42.toMatrix(); // convert to dense matrix
SparseMatrix<Real> sm42(m42); // sparse matrix from dense
// matrix (zero's not included)