This file provides a short description of algorithms implemented in file "capd/matrixAlgorithms/floatMatrixAlgorithms.hpp" from CAPD package - for more details see the source file.
The following function solves system of linear equations by using Gauss elimination algorithm
template<typename MatrixType, typename VectorType> VectorType gauss(MatrixType a, VectorType b);
Since both arguments are modified in this function, they are copied to the actual arguments a,b
Example of usage:
// we create a matrix double data[] = {1,2,1,2,2,2,2,0,1,97,3,4,1,19,1,15}; Matrix<double,0,0> P(4,4,data); Vector<double,0> v = P.row(0); // here we call the function gauss Vector<double,0> solution = gauss(P,v); std::cout << "gauss(P,v)=" << solution << std::endl; // this should be close to v std::cout << "verification of the result:\nP*gauss(P,v)=" << P*solution << std::endl;
The following function orthonormalizes columns of a nonsingular matrix. The function throws an exception when the matrix is singular.
template<typename MatrixType> void orthonormalize(MatrixType& Q);
The above function performs Gramm-Schmidt orthonormalization of columns of matrix Q sorted by decreasing norm. The argument Q is modified by the function and after calling orthonormalize(Q), the matrix Q contains result of this function.
Example of usage:
double data[] = {1,2,1,2,2,2,2,0,1,97,3,4,1,19,1,15}; Matrix<double,0,0> Q(4,4,data); std::cout << "orthonormalization of Q=" << Q << std::endl; orthonormalize(Q); std::cout << "is equal to " << Q << std::endl;
The following function computes QR-decomposition of a square nonsingular matrix.
template<typename MatrixType> void QR_decompose(const MatrixType& A, MatrixType& Q, MatrixType& R);
As a result we obtain orthogonal and upper triangular matrix. First argument of the function is a matrix to be decomposed. Matrices Q and R are modified by the function and contain the result.
Example of usage:
double data[] = {1,2,1,2,2,2,2,0,1,97,3,4,1,19,1,15}; Matrix<double,0,0> P(4,4,data), Q(4,4), R(4,4); QR_decompose(P,Q,R); // P should be equal to Q*R std::cout << "QR decomposition of P:" << std::endl; std::cout << "Q=" << Q << std::endl; std::cout << "R=" << R << std::endl; std::cout << "verification if P-Q*R is close to zero\n"; std::cout << "P-Q*R=" << P-Q*R << std::endl;
The function throws an exception when the matrix is singular.
The following function performs Jacobi rotation algorithm for computing diagonalization of a square symmetric matrix.
template<typename MatrixType> int symMatrixDiagonalize(const MatrixType& A, MatrixType& D, typename MatrixType::ScalarType diagonalizingRelTolerance = typename MatrixType::ScalarType(1e-5));
The argument D contains the result of the function. The last argument is a diagonalizing tolerance we want to obtain, namely if by 'd' and 's' we denote
the algorithm stops when s < d*diagonalizingRelTolerance or the next iteration would worsen the estimation. It should be noted that the procedure does not guarantee that such tolerance is achieved. This procedure is the first step in functions which compute spectral radius and maximal eigenvalue of a symmetric matrix.
Example of usage:
double data[] = {1,2,1,2,2,2,2,0,1,97,3,4,1,19,1,15}; Matrix<double,0,0> P(4,4,data), D(4,4); Matrix<double,0,0> S = 0.5*(P+Transpose(P)); // diagonalization of S=0.5*(P+P^T) symMatrixDiagonalize(S,D); std::cout << "Diagonalization of 0.5*(P+P^T)=" << D << std::endl;
The following function computes an upper bound for spectral radius of a symmetric matrix.
template<typename MatrixType> typename MatrixType::ScalarType spectralRadiusOfSymMatrix(const MatrixType &A, typename MatrixType::ScalarType diagonalizingRelTolerance = typename MatrixType::ScalarType(1e-5));
First, the function computes matrix which has the same eigenvalues and which is close to diagonal, next upper bound is computed from Gerschgorin theorem
Example of usage:
double data[] = {1,2,1,2,2,2,2,0,1,97,3,4,1,19,1,15}; Matrix<double,0,0> P(4,4,data); Matrix<double,0,0> S = 0.5*(P+Transpose(P)); std::cout << "spectralRadiusOfSymMatrix(S)=" << spectralRadiusOfSymMatrix(S) << std::endl;
The following function computes an upper bound for maximal eigenvalue of a symmetric matrix.
template<typename MatrixType> typename MatrixType::ScalarType maxEigenValueOfSymMatrix(const MatrixType &A, typename MatrixType::ScalarType diagonalizingRelTolerance = typename MatrixType::ScalarType(1e-5));
First, the function computes matrix which has the same eigenvalues and which is close to diagonal, next upper bound is computed from Gerschgorin theorem
Example of usage:
double data[] = {1,2,1,2,2,2,2,0,1,97,3,4,1,19,1,15}; Matrix<double,0,0> P(4,4,data); Matrix<double,0,0> S = 0.5*(P+Transpose(P)); std::cout << "maxEigenValueOfSymMatrix(S)=" << maxEigenValueOfSymMatrix(S) << std::endl;
The matrixAlgorithm module implements two functions for computing the inverse of a matrix. First function decomposes a matrix A=Q*R, where Q-orthogonal, R-upper diagonal and compute A-1 = R-1*QT
template<typename MatrixType> MatrixType inverseMatrix(const MatrixType &A);
The second function (recommended) computes the inverse matrix using Gauss elimination algorithm
template<typename MatrixType> MatrixType gaussInverseMatrix(const MatrixType &A);
Example of usage:
double data[] = {1,2,1,2,2,2,2,0,1,97,3,4,1,19,1,15}; Matrix<double,0,0> P(4,4,data); Matrix<double,0,0> invP = inverseMatrix(P); std::cout << "Inverse of P by using QR-decomposition:" << invP << std::endl; std::cout << "verification: P*inverseMatrix(P)=" << P*invP << std::endl; invP = gaussInverseMatrix(P); std::cout << "Inverse of P by using Gauss elimination:" << invP << std::endl; std::cout << "verification: P*gaussInverseMatrix(P)=" << P*invP << std::endl;
Author: Daniel Wilczak, last modified on July 7, 2008.