In that page, we detail functions that are not related to Lapack
| Transpose | replaces a matrix by its transpose |
| TransposeConj | replaces a matrix by its conjugate transpose |
| SetRow | modifies a row of the matrix |
| GetRow | extracts a row from the matrix |
| SetCol | modifies a column of the matrix |
| GetCol | extracts a column from the matrix |
| MaxAbs | returns highest absolute value of A |
| Norm1 | returns 1-norm of A |
| NormInf | returns infinity-norm of A |
| ApplyPermutation | permutes row and column numbers of a matrix |
| ApplyInversePermutation | permutes row and column numbers of a matrix |
| ScaleMatrix | multiplies rows and columns of a matrix by coefficients |
| ScaleLeftMatrix | multiplies rows of a matrix by coefficients |
| ScaleRightMatrix | multiplies columns of a matrix by coefficients |
| Copy | copies a sparse matrix into another one (conversion of format if needed) |
| ConvertMatrix_to_Coordinates | conversion of a sparse matrix into coordinates format |
| ConvertMatrix_from_Coordinates | conversion of a matrix given as a triplet (i, j, val) to a sparse matrix |
| ConvertToCSC | converts a sparse matrix to CSC (Compressed Sparse Column) format |
| ConvertToCSR | converts a sparse matrix to CSR (Compressed Sparse Row) format |
| GetSymmetricPattern | computes the sparsity pattern of A + A' |
| ConvertToSparse | converts dense matrices to sparse matrices by specifying a threshold. |
| Gauss | basic Gauss pivoting for dense matrices. |
| GaussSeidel | performs a Gauss-Seidel iteration for dense matrices. |
| SOR | applies successive over-relaxations to matrix |
| Cg, Gmres, BiCgSTAB, etc | solves iteratively a linear system |
| IsComplexMatrix | returns true if the matrix is complex |
| IsSymmetricMatrix | returns true if the matrix is symmetric |
| GetLowerTriangular | extracts lower triangular part of a matrix |
| GetUpperTriangular | extracts upper triangular part of a matrix |
| ReadCoordinateMatrix | reads a matrix in coordinate format (as in Matlab) |
| WriteCoordinateMatrix | writes a matrix in coordinate format (as in Matlab) |
| ReadHarwellBoeing | reads a matrix in Harwell-Boeing format |
| WriteHarwellBoeing | writes a matrix in Harwell-Boeing format |
| ReadMatrixMarket | reads a matrix in Matrix Market format |
| WriteMatrixMarket | writes a matrix in Matrix Market format |
| 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 |
| CopySubMatrix | extracts a sub-matrix from a sparse matrix |
void Transpose(Matrix&); void Transpose(const Matrix&, Matrix&); void TransposeConj(Matrix&); void TransposeConj(const Matrix&, Matrix&);
Transpose overwrites a matrix by its transpose, while TransposeConj overwrites a matrix by its conjugate transpose. You can also compute the tranpose (or conjugate transpose) and copy this transpose into another matrix.
Matrix<double> A(5, 5); A.Fill(); Transpose(A); Matrix<complex<double> > B(5, 5); // you fill B as you want, then overwrites it by conj(transpose(B)) TransposeConj(B); // You can keep the original matrix Matrix<double, General, ArrayRowSparse> As, Atrans; // the original matrix As is constructed then Transpose(As, Atrans); // Atrans contains now the transpose of As Matrix<complex<double>, General, RowSparse> Bs Btrans; // the original matrix Bs is constructed then TransposeConj(As, Atrans); // Btrans contains now the conjugate transpose of Bs
void SetRow(const Vector&, int, Matrix&);
This function modifies a row in the provided matrix. For sparse matrices, the expected vector is sparse.
Matrix<double> A(5, 5); A.Fill(); // now you construct a row Vector<double> row(5); row.FillRand(); // and you put it in A int irow = 1; SetRow(row, irow, A); // For sparse matrices the row is sparse Matrix<double, General, ArrayRowSparse> B(10, 10); Vector<double, VectSparse> row_sparse(2); row_sparse.Index(0) = 2; row_sparse.Value(0) = 1.2; row_sparse.Index(1) = 7; row_sparse.Value(1) = 2.4; // setting a row of B SetRow(row_sparse, irow, B);
void GetRow(const Matrix&, int, Vector&);
This function extracts a row from the provided matrix. For sparse matrices, the extracted vector is a sparse vector.
Matrix<double> A(5, 5); A.Fill(); // now you extract the first row Vector<double> row; GetRow(A, 0, row); // For sparse matrices, declare a sparse vector Matrix<double, General, RowSparse> B(10, 10); Vector<double, VectSparse> row_sparse; GetRow(B, 0, row_sparse);
void SetCol(const Vector&, int, Matrix&);
This function modifies a column in the provided matrix. For sparse matrices, the expected vector is sparse.
Matrix<double> A(5, 5); A.Fill(); // now you construct a column Vector<double> col(5); col.FillRand(); // and you put it in A int icol = 1; SetCol(col, icol, A); // For sparse matrices the column is sparse Matrix<double, General, ArrayRowSparse> B(10, 10); Vector<double, VectSparse> col_sparse(2); col_sparse.Index(0) = 2; col_sparse.Value(0) = 1.2; col_sparse.Index(1) = 7; col_sparse.Value(1) = 2.4; // setting a column of B SetCol(col_sparse, icol, B);
void GetCol(const Matrix&, int, Vector&);
This function extracts a column from the provided matrix. For sparse matrices, the extracted vector is a sparse vector.
Matrix<double> A(5, 5); A.Fill(); // now you extract the first column Vector<double> col; GetCol(A, 0, col); // For sparse matrices, declare a sparse vector Matrix<double, General, RowSparse> B(10, 10); Vector<double, VectSparse> col_sparse; GetCol(B, 0, col_sparse);
T MaxAbs(const Matrix<T>&);
This function returns the highest absolute value of a matrix.
Matrix<complex<double> > A(5, 5); A.Fill(); double module_max = MaxAbs(A);
T Norm1(const Matrix<T>&);
This function returns the 1-norm of a matrix.
Matrix<complex<double> > A(5, 5); A.Fill(); double anorm_one = Norm1(A);
T NormInf(const Matrix<T>&);
This function returns the infinity-norm of a matrix.
Matrix<complex<double> > A(5, 5); A.Fill(); double anorm_inf = NormInf(A);
void ApplyPermutation(Matrix&, const Vector<int>&, const Vector<int>&); void ApplyInversePermutation(Matrix&, const Vector<int>&, const Vector<int>&);
This function permutes a given matrix with the provided new row numbers and column numbers. ApplyInversePermutation(A, row, col) does the same operation as A(row, col) = A in Matlab, whereas ApplyPermutation is similar to A = A(row, col).
// you fill A as you wish Matrix<double, Symmetric, ArrayRowSymSparse> A(5, 5); // then new row and column numbers IVect row(5); // for symmetric matrix, // only second argument is actually used for the permutation // of both rows and columns in order to preserve symmetry ApplyInversePermutation(A, row, row); // for unsymmetric matrices, you can specify different permutations IVect col(5); Matrix<double, General, ArrayRowSparse> B(5, 5); ApplyInversePermutation(B, row, col); // reciprocal operation ApplyPermutation(B, row, col);
void ScaleMatrix(Matrix&, const Vector& L, const Vector& R);
This function multiplies each row and column with a coefficient, i.e. A is replaced by L*A*R where L and R are diagonal matrices and you provide the diagonal when you call ScaleMatrix.
// you fill A as you wish Matrix<double, Symmetric, ArrayRowSymSparse> A(5, 5); // then scaling vectors Vector<double> scale(5); // for symmetric matrix, row and column scaling must be the same // if you want to keep symmetry ScaleMatrix(A, scale, scale); // for unsymmetric matrices, you can specify different scalings Vector<double> scale_right(5); Matrix<double, General, ArrayRowSparse> B(5, 5); ScaleMatrix(B, scale, scale_right);
void ScaleLeftMatrix(Matrix&, const Vector& L);
This function multiplies each row with a coefficient, i.e. A is replaced by L*A where L is diagonal and you provide the diagonal when you call ScaleLeftMatrix. This function is not available for symmetric matrices since this operation would break the symmetry.
// you fill A as you wish Matrix<double, General, ArrayRowSparse> A(5, 5); // then scaling vector Vector<double> scale(5); ScaleLeftMatrix(A, scale);
void ScaleRightMatrix(Matrix&, const Vector& R);
This function multiplies each column with a coefficient, i.e. A is replaced by A*R where R is diagonal and you provide the diagonal when you call ScaleRightMatrix. This function is not available for symmetric matrices since this operation would break the symmetry.
// you fill A as you wish Matrix<double, General, ArrayRowSparse> A(5, 5); // then scaling vector Vector<double> scale(5); ScaleRightMatrix(A, scale);
void Copy(const Matrix&, Matrix2&);
This function copies a sparse matrix into another one. If the types of the matrices differ, a conversion is performed. However, not all the conversion have been implemented, so you may have a compilation error when copying some matrices.
// you fill A as you wish Matrix<double, General, ArrayRowSparse> A(5, 5); // then you can copy it to another form Matrix<double, General, RowSparse> B; // B does not need to be allocated Copy(A, B); // For other types that don't compile // you can use ConvertMatrix_to_Coordinates : Matrix<double, Symmetric, ColSymSparse> C; // conversion to triplet form (i, j, value) IVect IndRow, IndCol; Vector<double> Val; ConvertMatrix_to_Coordinates(A, IndRow, IndCol, Val, 0, true); // then C is filled ConvertMatrix_from_Coordinates(IndRow, IndCol, Val, C, 0);
void ConvertMatrix_to_Coordinates(const Matrix& A, Vector<int>& IndRow, Vector<int>& IndCol, Vector<T>& Val, int index, bool sym);
This function converts a sparse matrix to the triplet form (i, j, val) (coordinate format). The row and column numbers will start at index, therefore you can switch between 1-based indices and 0-based indices. If sym is true and the storage is symmetric, the lower and upper part of the matrix are transformed in coordinate format. If sym is false and the storage is symmetric, only the upper part of the matrix is converted. If the storage is unsymmetrix, sym is not used. The default values of index and sym are 0 and false.
// you fill A as you wish
Matrix<double, General, ArrayRowSparse> A(5, 5);
// conversion to triplet form (i, j, value)
IVect IndRow, IndCol;
Vector<double> Val;
ConvertMatrix_to_Coordinates(A, IndRow, IndCol, Val, 0, true);
// number of non-zero entries :
int nnz = IndRow.GetM();
for (int i = 0; i < nnz; i++)
{
cout << "Row index : " << IndRow(i) << endl;
cout << "Column index : " << IndCol(i) << endl;
cout << "value : " << Val(i) << endl;
}
void ConvertMatrix_from_Coordinates(Vector<int>& IndRow, Vector<int>& IndCol, Vector<T>& Val, Matrix& A, int index);
This function converts a triplet form (i, j, val) (coordinate format) to a sparse matrix. The row and column numbers are assumed to start at index, therefore you can switch between 1-based indices and 0-based indices. The default value of index is 0.
// creating a sparse matrix // A = | 1 0 0 2| // | -1 2 3 0| // | 0 1.2 0 2.5| int nnz = 7; IVect IndRow(nnz), IndCol(nnz); Vector<double> Val(nnz); IndRow(0) = 0; IndCol(0) = 0; Val(0) = 1.0; IndRow(1) = 0; IndCol(1) = 3; Val(1) = 2.0; IndRow(2) = 1; IndCol(2) = 0; Val(2) = -1.0; IndRow(3) = 1; IndCol(3) = 1; Val(3) = 2.0; IndRow(4) = 1; IndCol(4) = 2; Val(4) = 3.0; IndRow(5) = 2; IndCol(5) = 1; Val(5) = 1.2; IndRow(6) = 2; IndCol(6) = 3; Val(6) = 2.5; // conversion to a Seldon structure Matrix<double, General, RowSparse> ConvertMatrix_from_Coordinates(IndRow, IndCol, Val, A, 0);
void ConvertToCSC(const Matrix& A, Property& sym, Vector<int>& Ptr, Vector<int>& Ind, Vector<T>& Val, bool sym_pat);
This function converts a matrix to Compressed Sparse Column (CSC) format. Val stores the values of non-zero entries, Ind the row indexes of the non-zero entries, and Ptr the locations in Val of non-zero entries starting a column. This is the storage represented by ColSparse in Seldon. If Property is Symmetric, only upper part of the matrix will be converted (ColSymSparse storage). If sym_pat is true, the sparsity pattern is symmetrized, that is to say that if a(i, j) exists, then a(j, i) also exists. Default value of sym_pat is false. This feature is used for exemple in the interface of Pastix solver, since this solver requires that the sparsity pattern is symmetric (values may be not symmetric).
// you fill A as you wish Matrix<double, General, ArrayRowSparse> A(5, 5); // then you can retrieve Ptr, Ind, Val arrays of CSC format General prop; IVect Ptr, Ind; Vector<double> Val; ConvertToCSC(A, prop, Ptr, Ind, Val);
void ConvertToCSR(const Matrix& A, Property& sym, Vector<int>& Ptr, Vector<int>& Ind, Vector<T>& Val);
This function converts a matrix to Compressed Sparse Row (CSR) format. Val stores the values of non-zero entries, Ind the column indexes of the non-zero entries, and Ptr the locations in Val of non-zero entries starting a row. This is the storage represented by RowSparse in Seldon. If Property is Symmetric, only upper part of the matrix will be converted (RowSymSparse storage).
// you fill A as you wish Matrix<double, General, ArrayColSparse> A(5, 5); // then you can retrieve Ptr, Ind, Val arrays of CSR format General prop; IVect Ptr, Ind; Vector<double> Val; ConvertToCSR(A, prop, Ptr, Ind, Val);
void GetSymmetricPattern(const Matrix& A, Vector<int>& Ptr, Vector<int>& Ind);
This function extracts the profile of a given matrix. The pattern is stored in Compressed Sparse Row (CSR) format. Ind stores the column indexes of the non-zero entries, and Ptr the index of starting entries of a row.
// you fill A as you wish Matrix<double, General, ArrayColSparse> A(5, 5); // then you can retrieve Ptr, Ind arrays of CSR format // profile is symmetrized if needed IVect Ptr, Ind; GetSymmetricPattern(A, Ptr, Ind);
void ConvertToSparse(const Matrix& A, Matrix& B, const T& threshold);
This function converts a dense matrix to a sparse matrix. All values whose modulus is below or equal to threshold are skipped.
// you fill A as you wish Matrix<double, General, RowMajor> A(5, 5); // then you can convert it to a sparse matrix Matrix<double, General, RowSparse> B; ConvertToSparse(A, B, 1e-12); // and retrieve the number of non-zero entries int nnz = B.GetDataSize();
void Gauss(A, b);
This function overwrites b with the solution of A x = b by performing Gauss algorithm. This basic method is only implemented for dense matrices and without partial pivoting. It is probably to use GetLU or GetAndSolveLU.
// you fill A as you wish Matrix<double, General, RowMajor> A(5, 5); A.FillRand(); // forming the right hand side Vector<double> b(5); b.FillRand(); // then computation of the solution of A x = b // A is modified during the operation Vector<double> x(b); Gauss(A, x);
void GaussSeidel(A, x, b, nb_iterations); void GaussSeidel(A, x, b, nb_iterations, int type_algo);
This function overwrites b with an approximate solution of A x = b by performing nb_iterations Gauss-Seidel steps. This basic method is only implemented for dense matrices. For sparse matrices, you can use SOR function by choosing a relaxation parameter equal to one. If type_algo is equal to 0, a forward sweep is followed by a backward sweep (so that the associated operator is symmetric) for each iteration. If type_algo is equal to 2, only forward sweeps are applied, and if type_algo is equal to 3, only backward sweeps are applied. The default value of type_algo is equal to 2.
// you fill A as you wish Matrix<double, General, RowMajor> A(5, 5); A.FillRand(); // forming the right hand side Vector<double> b(5); b.FillRand(); // then computation of the solution of A x = b // with Gauss-Seidel algorithm Vector<double> x(5); x.Zero(); int nb_iterations = 100; GaussSeidel(A, x, b, nb_iterations, 0);
void IsComplexMatrix(const Matrix&)
This function returns true if the matrix is complex. It does not check if all the values are real, but merely returns true if the value type is complex (e.g. T = complex<double>).
// complex matrix
Matrix<complex<double>, General, RowMajor> A;
// IsComplexMatrix should return true
if (IsComplexMatrix(A))
{
cout << "A is complex" << endl;
}void IsSymmetricMatrix(const Matrix&)
This function returns true if the matrix is symmetric. It does not check if a(i,j) = a(j,i) for all i and j, but merely returns true if the property of the matrix is set to symmetric.
// complex matrix
Matrix<complex<double>, Symmetric, RowSymPacked> A;
// IsSymmetricMatrix should return true
if (IsSymmetricMatrix(A))
{
cout << "A is symmetric" << endl;
}void GetLowerTriangular(const Matrix&, Matrix&)
This function extracts the lower part of a matrix, this function is implemented only for dense storages.
// dense matrix int n = 6; Matrix<complex<double>, General, RowMajor> A(n, n); A.FillRand(); Matrix<complex<double>, General, RowMajor> L; // equivalent Matlab function L = tril(A) GetLowerTriangular(A, L);
void GetUpperTriangular(const Matrix&, Matrix&)
This function extracts the upper part of a matrix, this function is implemented only for dense storages.
// dense matrix int n = 6; Matrix<complex<double>, General, RowMajor> A(n, n); A.FillRand(); Matrix<complex<double>, General, RowMajor> U; // equivalent Matlab function U = triu(A) GetUpperTriangular(A, U);
void ReadCoordinateMatrix(istream& file_stream, Vector<int>& rows, Vector<int>& cols, Vector<T>& values, bool cplx)
This function reads a sparse matrix from a file (indexes start from 1). The file is expected to look like :
1 1 0.545432 1 3 -0.4349 1 7 33.42343 2 1 -8.43321 2 4 0.987987
The sparse matrix is stored in coordinate format, with the three arrays rows, col, vals, respctively the row indexes, column indexes and the values, that will be equal for this file :
rows = [1, 1, 1, 2, 2] cols = [1, 3, 7, 1, 4] values = [0.545432, -0.4349, 33.42343, -8.43321, 0.987987]
If you know already the number of non-zero entries, it is better to allocate these arrays before calling ReadCoordinateMatrix. If cplx is equal to true, a complex matrix is expected to be stored in the following file
1 1 0.545432 0.0 1 3 -0.4349 -1.23 1 7 33.42343 34.211 2 1 -8.43321 0.0113 2 4 0.987987 0.0
// sparse matrix in coordinate format
Vector<int> rows, cols;
Vector<double> values;
// initializing this matrix by reading it in a file
// first line of the file contains m, n, nnz
int m, n, nnz;
ifstream file_in("mat.dat");
if (!file_in.good())
cout << "File does not exist" << endk;
file_in >> m >> n >> nnz;
rows.Reallocate(nnz);
cols.Reallocate(nnz);
values.Reallocate(nnz);
ReadCoordinateMatrix(file_in, rows, col, values);
file_in.close();void WriteCoordinateMatrix(ostream& file_stream, const Vector<int>& rows, const Vector<int>& cols, const Vector<T>& values, bool cplx)
This function writes a sparse matrix in a file (indexes start from 1). If we have taken the following arguments
rows = [1, 1, 1, 2, 2] cols = [1, 3, 7, 1, 4] values = [0.545432, -0.4349, 33.42343, -8.43321, 0.987987]
The file will look like :
1 1 0.545432 1 3 -0.4349 1 7 33.42343 2 1 -8.43321 2 4 0.987987
If cplx is equal to true, the values of a complex matrix will be stored with two separated columns for the real and imaginary part (no parenthesis).
// sparse matrix in coordinate format
Vector<int> rows, cols;
Vector<double> values;
// initializaing those arrays
int m = 6, n, 8, nnz = 20;
rows.Reallocate(nnz); cols.Reallocate(nnz);
values.Reallocate(nnz);
ofstream file_out("mat.dat");
// writing m, n, nnz
file_out << m << " " << n << " " << nnz << endl;
// then writing triplets i, j, val
WriteCoordinateMatrix(file_out, rows, col, values);
file_out.close();void ReadHarwellBoeing(string file_name, Matrix& A)
This function reads a sparse matrix from a file (Harwell-Boeing format).
Matrix<double, General, ArrayRowSparse> A;
ReadHarwellBoeing("mat.dat", A);void WriteHarwellBoeing(const Matrix& A, const string& file_name)
This function writes a sparse matrix in a file (Harwell-Boeing format).
Matrix<double, General, ArrayRowSparse> A; WriteHarwellBoeing(A, "mat.dat");
void ReadMatrixMarket(string file_name, Matrix& A)
This function reads a sparse matrix from a file (Matrix Market format).
Matrix<double, General, ArrayRowSparse> A;
ReadMatrixMarket("mat.mtx", A);void WriteMatrixMarket(const Matrix& A, const string& file_name)
This function writes a sparse matrix in a file (Matrix Market format).
Matrix<double, General, ArrayRowSparse> A; WriteMatrixMarket(A, "mat.mtx");
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.
// 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);
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.
// 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);
void GetColSum(Vector& V, const Matrix& A);
This function computes the sum of absolute values of elements of each column.
// for a sparse matrix
Matrix<double, General, ArrayRowSparse> A;
Vector<double> V;
// computing V_j = \sum_i |a_{i,j}|
GetColSum(V, A);
void GetRowSum(Vector& V, const Matrix& A);
This function computes the sum of absolute values of elements of each row.
// for a sparse matrix
Matrix<double, General, ArrayRowSparse> A;
Vector<double> V;
// computing V_i = \sum_j |a_{i,j}|
GetRowSum(V, A);
void GetRowColSum(Vector& Vrow, Vector& Vcol, const Matrix& A);
This function computes the sum of absolute values of elements of each row and column.
// 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);
void CopySubMatrix(const Matrix& A, const IVect& row, const IVect& col, Matrix& B);
This function extracts a sub-matrix B from the global matrix A.
// 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) // you need to fill row numbers and col numbers to keep IVect row(10), col(10); row.Fill(); col.Fill(); CopySubMatrix(A, row, col, B); // B = A[m1:m2, n1:n2] (Python notation) row.Reallocate(m2-m1); col.Reallocate(n2-n1); for (int i = m1; i < m2; i++) row(i-m1) = i; for (int i = n1; i < n2; i++) col(i-n1) = i; CopySubMatrix(A, row, col, B);