Additional Seldon-like functions for linear algebra
Additional functions related to linear algebra are available in Montjoie.
| GetCol | extracts several columns of a sparse matrix |
| EraseRow | erases several rows of a sparse matrix |
| EraseCol | erases several columns of a sparse matrix |
| GetRowSum | sums absolute values of non-zero entries by row |
| GetColSum | sums absolute values of non-zero entries by column |
| GetRowColSum | sums absolute values of non-zero entries by row and by column |
| GetSubMatrix | extracts a sub-matrix from a sparse matrix |
| CopyReal | extracts the real part of a complex matrix |
| ExtractDistributedInfo | constructs distributed integer arrays from a subset of rows |
| CompressMatrix | compresses a matrix by keeping a subset of rows |
| ExtractSubMatrix | extracts a sub-matrix from a sparse matrix |
| GetAbsoluteValue | replaces a matrix by its absolute value |
| GetSquareRoot | replaces a matrix by its square root |
| GetExponential | replaces a matrix by its exponential |
| GetFunctionMatrix | replaces a complex matrix A by f(A) for any function f |
| GetFunctionMatrixReal | replaces a real matrix A by f(A) for any function f |
| GetHouseholderNormale | Computes the normale appearing in Householder transformation |
| GetReorthogonalization | Orthogonalize a vector with previous stored Householder transformations |
| ApplyHouseholderTransformation | Applies Householder transformation to a single vector |
GetCol
Syntax
void GetCol(const Matrix& A, const IVect& col_num, Vector<VectorSparse>& V);
This function extracts the asked columns of A in a sparse vector of sparse vectors (denoted V). The column numbers to extract are provided in the array col_num. These numbers should be provided in ascending order, otherwise the columns of V will not be sorted.
Example :
// for a sparse matrix Matrix<double, General, ArrayRowSparse> A; // if you want to extract the column 2 and 5 IVect col_num(2); col_num(0) = 2; col_num(1) = 5; // columns are stored in V (which is sparse) Vector<Vector<double, VectSparse>, VectSparse> V; GetCol(A, col_num, V); // the columns are stored in V(2) and V(5), other components of V are empty (sparse structure).
Location :
FunctionMatrixExtraction.hxx
FunctionMatrixExtraction.cxx
EraseRow
Syntax
void EraseRow(IVect& num, Matrix& A);
This function erases some rows of the matrix A. The numbers of the rows to erase are provided in the array num.
Example :
// filling a sparse Matrix<double, General, ArrayRowSparse> A; // then you can erase some rows IVect num(2); num(0) = 1; num(1) = 4; // for example the row 1 and 4 EraseRow(num, A);
Location :
FunctionMatrixExtraction.hxx
FunctionMatrixExtraction.cxx
EraseCol
Syntax
void EraseCol(IVect& num, Matrix& A);
This function erases some columns of the matrix A. The numbers of the columns to erase are provided in the array num.
Example :
// filling a sparse Matrix<double, General, ArrayRowSparse> A; // then you can erase some columns IVect num(2); num(0) = 1; num(1) = 4; // for example the column 1 and 4 EraseCol(num, A);
Location :
FunctionMatrixExtraction.hxx
FunctionMatrixExtraction.cxx
GetColSum
Syntax
void GetColSum(Vector& V, const Matrix& A);
This function computes the sum of absolute values of elements of each column.
Example :
// for a sparse matrix
Matrix<double, General, ArrayRowSparse> A;
Vector<double> V;
// computing V_j = \sum_i |a_{i,j}|
GetColSum(V, A);
Location :
FunctionMatrixExtraction.hxx
FunctionMatrixExtraction.cxx
GetRowSum
Syntax
void GetRowSum(Vector& V, const Matrix& A);
This function computes the sum of absolute values of elements of each row.
Example :
// for a sparse matrix
Matrix<double, General, ArrayRowSparse> A;
Vector<double> V;
// computing V_i = \sum_j |a_{i,j}|
GetRowSum(V, A);
Location :
FunctionMatrixExtraction.hxx
FunctionMatrixExtraction.cxx
GetRowColSum
Syntax
void GetRowColSum(Vector& Vrow, Vector& Vcol, const Matrix& A);
This function computes the sum of absolute values of elements of each row and column.
Example :
// for a sparse matrix
Matrix<double, General, ArrayRowSparse> A;
Vector<double> Vrow, Vcol;
// computing Vrow_i = \sum_j |a_{i,j}|
// computing Vcol_j = \sum_i |a_{i,j}|
GetRowSum(Vrow, Vcol, A);
Location :
FunctionMatrixExtraction.hxx
FunctionMatrixExtraction.cxx
GetSubMatrix
Syntax
void GetSubMatrix(const Matrix& A, int m, int n, Matrix& B); void GetSubMatrix(const Matrix& A, int m1, int m2, int n1, int n2, Matrix& B);
This function extracts a sub-matrix B from the global matrix A.
Example :
// for a sparse matrix Matrix<double, General, ArrayRowSparse> A; // extracting the first 10 rows and columns of A // B = A[0:10, 0:10] (Python notation) GetSubMatrix(A, 0, 10, B); // B = A[m1:m2, n1:n2] (Python notation) GetSubMatrix(A, m1, m2, n1, n2, B);
Location :
FunctionMatrixExtraction.hxx
FunctionMatrixExtraction.cxx
CopyReal
Syntax
void CopyReal(const Matrix& A, Matrix& B);
This function copies the real part of A in matrix B.
Example :
// for a sparse matrix Matrix<complex<double>, General, ArrayRowSparse> A; Matrix<double, General, ArrayRowSparse> B; // only the real part of A is copied into B CopyReal(A, B);
Location :
FunctionMatrixExtraction.hxx
FunctionMatrixExtraction.cxx
GetAbsoluteValue
Syntax
void GetAbsoluteValue(Matrix& A);
his function replaces the matrix A by its absolute value. The absolute value of a matrix is computed by diagonalizing the matrix A, and replacing the diagonal D by its absolute value. This function is implemented for dense general or dense symmetric matrices, it is not implemented for hermitian, triangular or sparse matrices.
Example :
// for example, we fill a symmetric matrix
Matrix<double, Symmetric, RowSymPacked> A(n, n);
// if A = P D P^{-1} where D is a diagonal matrix
// A is replaced by P |D| P^{-1}
GetAbsoluteValue(A);
// if the second argument is false, instead of the modulus
// we take f(x) = -x if real(x) < 0, x otherwise
Matrix<complex<double>, General, RowMajor> B(n, n);
GetAbsoluteValue(B, false);
Location :
GetSquareRoot
Syntax
void GetSquareRoot(Matrix& A);
This function replaces the matrix A by its square root. The square root of a matrix is computed by diagonalizing the matrix A, and replacing the diagonal D by its square root. For a real matrix, this function returns the result as a real matrix even if the result is complex. This function is implemented for dense general or dense symmetric matrices, it is not implemented for hermitian, triangular or sparse matrices.
Example :
// for example, we fill a symmetric matrix
Matrix<double, Symmetric, RowSymPacked> A(n, n);
// if A = P D P^{-1} where D is a diagonal matrix
// A is replaced by P sqrt(D) P^{-1}
// here the result should be correct if all eigenvalues are positive
GetSquareRoot(A);
Location :
GetExponential
Syntax
void GetExponential(Matrix& A);
This function replaces the matrix A by its exponential. The exponential of a matrix is computed by diagonalizing the matrix A, and replacing the diagonal D by its exponential. This function is implemented for dense general or dense symmetric matrices, it is not implemented for hermitian, triangular or sparse matrices.
Example :
// for example, we fill a complex symmetric matrix
Matrix<complex<double>, Symmetric, RowSymPacked> A(n, n);
// if A = P D P^{-1} where D is a diagonal matrix
// A is replaced by P exp(D) P^{-1}
GetExponential(A);
Location :
GetFunctionMatrix
Syntax
void GetFunctionMatrix(Matrix& A, f);
This function replaces the complex matrix A by f(A). This result matrix is computed by diagonalizing the matrix A, and replacing the diagonal D by f(D). The function f is a complex-valued function (it returns a complex and takes a complex as argument). GetFunctionMatrix is implemented for dense general or dense symmetric matrices, it is not implemented for hermitian, triangular or sparse matrices.
Example :
// for example, we fill a complex matrix
Matrix<complex<double>, General, ColMajor> A(n, n);
// if A = P D P^{-1} where D is a diagonal matrix
// A is replaced by P f(D) P^{-1}
// the function f is given as the second argument (e.g. the sine function)
GetFunctionMatrix(A, sin);
Location :
GetFunctionMatrixReal
Syntax
void GetFunctionMatrixReal(Matrix& A; f);
This function replaces the real matrix A by real(f(A)). This result matrix is computed by diagonalizing the matrix A, and replacing the diagonal D by f(D). The function f is a complex-valued function (it returns a complex and takes a complex as argument). The resulting matrix is always real because we overwrite A with the real part of f(A). This function is implemented for dense general or dense symmetric matrices, it is not implemented for hermitian, triangular or sparse matrices. Here, the result is assumed to be real.
Example :
// for example, we fill a real matrix
Matrix<double, General, ColMajor> A(n, n);
// if A = P D P^{-1} where D is a diagonal matrix
// A is replaced by P f(D) P^{-1}
// the function f is given as the second argument (e.g. the cosine function)
// this function is defined for complex arguments
GetFunctionMatrixReal(A, cos);
Location :
GetHouseholderNormale
Syntax
T GetHouseholderNormale(Vector<T>& U);
This function replaces the vector U by the normale N of the householder transformation that transforms the vector into the vector V = (alpha, 0, ..., 0). The transformation is given as H = I - 2 N N^T. The function returns the first component of vector V (coefficient alpha).
Example :
// vector U Vector<double> U(n), N; U.FillRand(); // we compute Householder transformation trough normale N N = U; double alpha = GetHouseholderNormale(N); // then you can apply Householder transformation for any vector X Vector<double> X(n); X.FillRand(); ApplyHouseholderTransformation(N, 0, X); // if you apply it to the original vector U, you should obtain // the vector V = (alpha, 0, 0, ..., 0) ApplyHouseholderTransformation(N, 0, U); cout << "V = " << U << endl;
Location :
FactorisationLU.hxx
FactorisationLU.cxx
GetReorthogonalization
Syntax
T GetReorthogonalization(Vector<Vector<T> >& H, int k, Vector<T>& Q);
This function orthogonalizes Q with respect to householder transformations stored in H and appends a new Householder transformation defined from Q in H.
Example :
// to be done
Location :
FactorisationLU.hxx
FactorisationLU.cxx
ApplyHouseholderTransformation
Syntax
void ApplyHouseholderTransformation(const Vector<T>& N, int k, Vector<T>& X);
This function applies the Householder transformation H = I - 2 N N^T given trough the normale N. The transformation is applied to vector X(k:end).
Example :
// vector U Vector<double> U(n), N; U.FillRand(); // we compute Householder transformation trough normale N N = U; double alpha = GetHouseholderNormale(N); // then you can apply Householder transformation for any vector X Vector<double> X(n); X.FillRand(); ApplyHouseholderTransformation(N, 0, X); // if you apply it to the original vector U, you should obtain // the vector V = (alpha, 0, 0, ..., 0) ApplyHouseholderTransformation(N, 0, U); cout << "V = " << U << endl;
Location :
FactorisationLU.hxx
FactorisationLU.cxx
ExtractDistributedInfo
Syntax
void ExtractDistributedInfo(const DistributedMatrixIntegerArray& A,
const Vector<int>& num, const Vector<int>& index,
DistributedMatrixIntegerArray& B)
This function constructs the object B by selecting a subset of rows of the object A. The local row numbers that are kept are given in the second argument (array num). The third argument is the reciprocal array (i.e. index(num(i)) is equal to i, index(j) is equal to -1 if j is not contained in num). This function is useful, if you want to construct a submatrix of a distributed matrix.
Example :
DistributedMatrix<double, General, ArrayRowSparse> A;
DistributedMatrixIntegerArray Ainfo;
// let's say that we start from a matrix with 8 rows
// proc 0 owns rows 0, 2, 4, 5, 7
// proc 1 owns rows 1, 3, 4, 6, 7
// rows 4 and 7 are shared
Vector<int> glob_num(5);
if (MPI::COMM_WORLD.Get_rank() == 0)
{
glob_num(0) = 0; glob_num(0) = 2; glob_num(0) = 4; glob_num(0) = 5; glob_num(0) = 7;
}
else
{
glob_num(0) = 1; glob_num(0) = 3; glob_num(0) = 4; glob_num(0) = 6; glob_num(0) = 7;
}
// we construct the object Ainfo
A.Init(glob_num, MPI::COMM_WORLD, Ainfo);
// then we want to construct a sub-matrix B with only rows 0, 1, 4, 6
// proc 0 will own rows 0, 4
// proc 1 will own rows 1, 4, 6
DistributedMatrix<double, General, ArrayRowSparse> B;
DistributedMatrixIntegerArray Binfo;
// here num will be the local numbers that are kept
Vector<int> num, index(glob_num.GetM());
index.Fill(-1); // we fill with -1 to detect rows that are not kept in B
if (MPI::COMM_WORLD.Get_rank() == 0)
{
num.Reallocate(2);
num(0) = 0; num(1) = 2; // 0, 2 are local numbers of global rows 0, 4
index(0) = 0; index(2) = 1; // reciprocal array index(num(i)) = i
}
else
{
num.Reallocate(3);
num(0) = 0; num(1) = 2; num(2) = 3; // 0, 2, 3 are local numbers of global rows 1, 4, 6
index(0) = 0; index(2) = 1; index(3) = 2; // reciprocal array index(num(i)) = i
}
// then instead of computing Binfo by calling B.Init, you can call ExtractDistributedInfo
// it should be faster (less communication)
ExtractDistributedInfo(Ainfo, num, index, Binfo);
// and then you can initialize B with Binfo
B.Reallocate(num.GetM(), num.GetM());
B.Init(Binfo);
// and then fill B
B.AddInteraction(0, 0, 2.5);
// etc
Location :
FunctionMatrixExtraction.hxx
FunctionMatrixExtraction.cxx
CompressMatrix
Syntax
void CompressMatrix(Matrix& A, Vector<int>& index);
This function compresses the current matrix into a smaller one. Only rows such that index(i) ≥ 0 are conserved. When index(i) is greater or equal to 0, it represents the new row number of the compressed matrix.
Example :
Matrix<double, General, ArrayRowSparse> A; A.Reallocate(n, n); A.Fill(); // if you want to keep rows 1, 2, 4, 6, 7 Vector<int> row_to_keep(5); row_to_keep(0) = 1; row_to_keep(1) = 2; row_to_keep(2) = 4; row_to_keep(3) = 6; row_to_keep(4) = 7; // you need to construct the index array (reciprocal array with -1 for rows to remove) Vector<int> index(n); index.Fill(-1); for (int i = 0; i < row_to_keep.GetM(); i++) index(row_to_keep(i)) = i; // you can compress the matrix // only rows 1, 2, 4, 6, 7 are kept and their new numbers are 0, 1, 2, 3, 4 CompressMatrix(A, index);
Location :
FunctionMatrixExtraction.hxx
FunctionMatrixExtraction.cxx
ExtractSubMatrix
Syntax
void ExtractSubMatrix(const Matrix& A, const Vector<int>& row_num, const Vector<int> index_row, const Vector<int>& col_num, const Vector<int> index_col, Matrix& B);
This function extracts a sub-matrix B from a larger matrix A. Only rows (respectively columns) such that index_row(i) ≥ 0 (resp. index_col(i)) are conserved. When index_row(i) is greater or equal to 0, it represents the new row number of the sub-matrix. The array row_num (resp. col_num) is the reciprocal array of index_row (resp. index_col), and stores the row numbers of the extracted rows. If you want to overwrite the matrix A by the sub-matrix B, and that row_num is equal to col_num, then you should call CompressMatrix instead.
Example :
Matrix<double, General, ArrayRowSparse> A; A.Reallocate(n, n); A.Fill(); // if you want to extract rows 1, 2, 4, 6, 7 and columns 2, 3, 5, 7 Vector<int> row_num(5), col_num(4); row_num(0) = 1; row_num(1) = 2; row_num(2) = 4; row_num(3) = 6; row_num(4) = 7; col_num(0) = 2; col_num(1) = 3; col_num(2) = 5; col_num(3) = 7; // you need to construct the index array (reciprocal array) Vector<int> index_row(n), index_col(n); index_row.Fill(-1); index_col.Fill(-1); for (int i = 0; i < row_num.GetM(); i++) index_row(row_num(i)) = i; for (int i = 0; i < col_num.GetM(); i++) index_col(col_num(i)) = i; // you can extract the sub-matrix, result is B Matrix<double, General, ArrayRowSparse> B; ExtractSubMatrix(A, row_num, index_row, col_num, index_col, B);