# 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).


### 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);


### 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);


### 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);


### 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);


### 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);


### 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);


### 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);


### 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);


### 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);


### 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);


### 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);


### 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);


### 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;


### 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


### 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;


### 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
// etc



### 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);


### 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);