Class for solving linear systems and performing matrix operations. More...

#include <linsys.hpp>

Public Member Functions
	LinSys (const std::vector< std::vector< double >> &mat)
	Constructor for LinSys class. More...

void	display_mtx () const
	Display the augmented matrix. More...

void	display (const std::vector< std::vector< double >> &mat) const

std::vector< double >	solve_gauss ()
	Solve the linear system using Gaussian Elimination. More...

double	determinant () const
	Calculate the determinant of the matrix. More...

void	invert_mtx ()
	Invert the matrix. More...

std::pair< std::vector< std::vector< double > >, std::vector< std::vector< double > > >	lu_decomp ()
	Perform LU decomposition of the matrix. More...

void	solve_lu ()
	Solve the linear system using LU decomposition. More...

void	solve_cholesky ()
	Solve the linear system using Cholesky decomposition. More...

void	solve_jacobi (int maxIterations=100, double tolerance=1e-10)
	Solve the linear system using Jacobi iteration. More...

bool	is_symmetric () const
	Check if the matrix is symmetric. More...

double	frobenius_norm () const
	Calculate the Frobenius norm of the matrix. More...

double	one_norm () const
	Calculate the 1-norm of the matrix. More...

double	inf_norm () const
	Calculate the infinity norm of the matrix. More...

void	gram_schmidt ()
	Perform Gram-Schmidt orthogonalization on the matrix. More...

bool	diagonally_dominant () const
	Check if the matrix is diagonally dominant. More...

bool	is_consistent () const
	Check if the linear system is consistent. More...

bool	is_homogeneous () const
	Check if the linear system is homogeneous. More...

Public Attributes
std::vector< std::vector< double > >	matrix

int	num_rows

int	num_cols

Detailed Description

Class for solving linear systems and performing matrix operations.

Definition at line 48 of file linsys.hpp.

Constructor & Destructor Documentation

◆ LinSys()

gpmp::linalg::LinSys::LinSys ( const std::vector< std::vector< double >> & mat )

Constructor for LinSys class.

Parameters

mat	The augmented matrix [A\|B] representing the linear system

Definition at line 40 of file linsys.cpp.

     : matrix(mat) {
     num_rows = matrix.size();
     num_cols = matrix[0].size();
 }

References matrix, num_cols, and num_rows.

Member Function Documentation

◆ determinant()

double gpmp::linalg::LinSys::determinant ( ) const

Calculate the determinant of the matrix.

Returns: Determinant of the matrix

Definition at line 107 of file linsys.cpp.

                                            {
     // assuming square matrix
     if (num_rows != num_cols) {
         std::cerr << "Error: Determinant is undefined for non-square matrices."
                   << std::endl;
         return 0.0;
     }
  
     std::vector<std::vector<double>> temp_mtx = matrix;
     double det = 1.0;
  
     for (int i = 0; i < num_rows - 1; ++i) {
         for (int k = i + 1; k < num_rows; ++k) {
             double factor = temp_mtx[k][i] / temp_mtx[i][i];
             for (int j = i; j < num_cols; ++j) {
                 temp_mtx[k][j] -= factor * temp_mtx[i][j];
             }
         }
     }
  
     for (int i = 0; i < num_rows; ++i) {
         det *= temp_mtx[i][i];
     }
  
     return det;
 }

Referenced by is_symmetric().

◆ diagonally_dominant()

bool gpmp::linalg::LinSys::diagonally_dominant ( ) const

Check if the matrix is diagonally dominant.

Returns: True if the matrix is diagonally dominant, false otherwise

Definition at line 396 of file linsys.cpp.

                                                  {
     for (int i = 0; i < num_rows; ++i) {
         double diagonalValue = std::abs(matrix[i][i]);
         double rowSum = 0.0;
         for (int j = 0; j < num_cols; ++j) {
             if (j != i) {
                 rowSum += std::abs(matrix[i][j]);
             }
         }
         if (diagonalValue <= rowSum) {
             return false;
         }
     }
     return true;
 }

◆ display()

void gpmp::linalg::LinSys::display ( const std::vector< std::vector< double >> & mat ) const

Definition at line 56 of file linsys.cpp.

                                                    {
     uint32_t rows = mat.size();
     uint32_t cols = mat[0].size();
  
     for (uint32_t i = 0; i < rows; ++i) {
         for (uint32_t j = 0; j < cols; ++j) {
             std::cout << mat[i][j] << " ";
         }
         std::cout << std::endl;
     }
 }

References test_linalg::cols, and test_linalg::rows.

◆ display_mtx()

void gpmp::linalg::LinSys::display_mtx ( ) const

Display the augmented matrix.

Definition at line 47 of file linsys.cpp.

                                          {
     for (int i = 0; i < num_rows; ++i) {
         for (int j = 0; j < num_cols; ++j) {
             std::cout << matrix[i][j] << " ";
         }
         std::cout << std::endl;
     }
 }

◆ frobenius_norm()

double gpmp::linalg::LinSys::frobenius_norm ( ) const

Calculate the Frobenius norm of the matrix.

Returns: The Frobenius norm of the matrix

Definition at line 338 of file linsys.cpp.

                                               {
     double sum = 0.0;
     for (int i = 0; i < num_rows; ++i) {
         for (int j = 0; j < num_cols; ++j) {
             sum += std::pow(matrix[i][j], 2);
         }
     }
     return std::sqrt(sum);
 }

◆ gram_schmidt()

void gpmp::linalg::LinSys::gram_schmidt ( )

Perform Gram-Schmidt orthogonalization on the matrix.

Definition at line 375 of file linsys.cpp.

                                     {
     std::vector<std::vector<double>> orthoBasis(
         num_rows,
         std::vector<double>(num_cols, 0.0));
  
     for (int j = 0; j < num_cols; ++j) {
         for (int i = 0; i < num_rows; ++i) {
             double projection = 0.0;
             for (int k = 0; k < i; ++k) {
                 projection += (matrix[i][j] * orthoBasis[k][j]) /
                               (orthoBasis[k][j] * orthoBasis[k][j]);
             }
             orthoBasis[i][j] = matrix[i][j] - projection;
         }
     }
  
     // Replace the matrix with the orthogonalized basis
     matrix = orthoBasis;
 }

◆ inf_norm()

double gpmp::linalg::LinSys::inf_norm ( ) const

Calculate the infinity norm of the matrix.

Returns: The infinity norm of the matrix

Definition at line 362 of file linsys.cpp.

                                         {
     double maxRowSum = 0.0;
     for (int i = 0; i < num_rows; ++i) {
         double rowSum = 0.0;
         for (int j = 0; j < num_cols; ++j) {
             rowSum += std::abs(matrix[i][j]);
         }
         maxRowSum = std::max(maxRowSum, rowSum);
     }
     return maxRowSum;
 }

◆ invert_mtx()

void gpmp::linalg::LinSys::invert_mtx ( )

Invert the matrix.

Definition at line 135 of file linsys.cpp.

                                   {
     // assuming square matrix
     if (num_rows != num_cols) {
         std::cerr
             << "Error: Matrix inversion is undefined for non-square matrices."
             << std::endl;
         return;
     }
  
     std::vector<std::vector<double>> identity(
         num_rows,
         std::vector<double>(num_cols, 0.0));
  
     // augment the matrix with the identity matrix
     for (int i = 0; i < num_rows; ++i) {
         matrix[i].insert(matrix[i].end(),
                          identity[i].begin(),
                          identity[i].end());
     }
  
     // perform Gaussian elimination
     solve_gauss();
  
     // separate the inverse matrix
     for (int i = 0; i < num_rows; ++i) {
         matrix[i].erase(matrix[i].begin(), matrix[i].begin() + num_rows);
     }
 }

◆ is_consistent()

bool gpmp::linalg::LinSys::is_consistent ( ) const

Check if the linear system is consistent.

Returns: True if the system is consistent, false otherwise

Definition at line 413 of file linsys.cpp.

                                            {
     for (int i = 0; i < num_rows; ++i) {
         bool allZeros = true;
         for (int j = 0; j < num_cols - 1; ++j) {
             if (fabs(matrix[i][j]) > std::numeric_limits<double>::epsilon()) {
  
                 allZeros = false;
                 break;
             }
         }
         if (allZeros && fabs(matrix[i][num_cols - 1]) >
                             std::numeric_limits<double>::epsilon()) {
             return false;
         }
     }
     return true;
 }

◆ is_homogeneous()

bool gpmp::linalg::LinSys::is_homogeneous ( ) const

Check if the linear system is homogeneous.

Returns: True if the system is homogeneous, false otherwise

Definition at line 432 of file linsys.cpp.

                                             {
     for (int i = 0; i < num_rows; ++i) {
         if (fabs(matrix[i][num_cols - 1]) >
             std::numeric_limits<double>::epsilon()) {
             return false;
         }
     }
     return true;
 }

◆ is_symmetric()

bool gpmp::linalg::LinSys::is_symmetric ( ) const

Check if the matrix is symmetric.

Returns: True if the matrix is symmetric, false otherwise

Definition at line 302 of file linsys.cpp.

                                           {
     // assuming square matrix
     if (num_rows != num_cols) {
         return false;
     }
  
     // check symmetry
     for (int i = 0; i < num_rows; ++i) {
         for (int j = 0; j < num_cols; ++j) {
             // if (matrix[i][j] != matrix[j][i]) {
             if (fabs(matrix[i][j] - matrix[j][i]) >
                 std::numeric_limits<double>::epsilon()) {
                 return false;
             }
         }
     }
  
     // check positive definiteness
     for (int i = 0; i < num_rows; ++i) {
         std::vector<std::vector<double>> submatrix(i + 1,
                                                    std::vector<double>(i + 1));
         for (int j = 0; j <= i; ++j) {
             for (int k = 0; k <= i; ++k) {
                 submatrix[j][k] = matrix[j][k];
             }
         }
  
         double det = LinSys(submatrix).determinant();
         if (det <= 0.0) {
             return false;
         }
     }
  
     return true;
 }

References determinant().

◆ lu_decomp()

std::pair< std::vector< std::vector< double > >, std::vector< std::vector< double > > > gpmp::linalg::LinSys::lu_decomp ( )

Perform LU decomposition of the matrix.

Returns: the Lower and Upper matrices of the decomposition

Definition at line 165 of file linsys.cpp.

                             {
     std::vector<std::vector<double>> L(num_rows,
                                        std::vector<double>(num_cols, 0.0));
     std::vector<std::vector<double>> U(num_rows,
                                        std::vector<double>(num_cols, 0.0));
  
     for (int i = 0; i < num_rows; ++i) {
         // Set diagonal elements of L to 1
         L[i][i] = 1.0;
  
         for (int j = i; j < num_cols; ++j) {
             U[i][j] = matrix[i][j];
             for (int k = 0; k < i; ++k) {
                 U[i][j] -= L[i][k] * U[k][j];
             }
         }
  
         for (int j = i + 1; j < num_rows; ++j) {
             L[j][i] = matrix[j][i];
             for (int k = 0; k < i; ++k) {
                 L[j][i] -= L[j][k] * U[k][i];
             }
             L[j][i] /= U[i][i];
         }
     }
  
     return {L, U};
 }

◆ one_norm()

double gpmp::linalg::LinSys::one_norm ( ) const

Calculate the 1-norm of the matrix.

Returns: The 1-norm of the matrix

Definition at line 349 of file linsys.cpp.

                                         {
     double maxColSum = 0.0;
     for (int j = 0; j < num_cols; ++j) {
         double colSum = 0.0;
         for (int i = 0; i < num_rows; ++i) {
             colSum += std::abs(matrix[i][j]);
         }
         maxColSum = std::max(maxColSum, colSum);
     }
     return maxColSum;
 }

◆ solve_cholesky()

void gpmp::linalg::LinSys::solve_cholesky ( )

Solve the linear system using Cholesky decomposition.

Definition at line 214 of file linsys.cpp.

                                       {
     // assuming symmetric and positive-definite matrix
     if (!is_symmetric()) {
         std::cerr << "Error: Cholesky decomposition is only applicable to"
                      "symmetric positive-definite matrices\n";
     }
  
     for (int i = 0; i < num_rows; ++i) {
         for (int j = 0; j <= i; ++j) {
             if (i == j) {
                 double sum = 0.0;
                 for (int k = 0; k < j; ++k) {
                     sum += std::pow(matrix[j][k], 2);
                 }
                 matrix[j][j] = std::sqrt(matrix[j][j] - sum);
             }
  
             else {
                 double sum = 0.0;
                 for (int k = 0; k < j; ++k) {
                     sum += matrix[i][k] * matrix[j][k];
                 }
                 matrix[i][j] = (matrix[i][j] - sum) / matrix[j][j];
             }
         }
     }
  
     // forward substitution
     for (int i = 0; i < num_rows; ++i) {
         double sum = 0.0;
         for (int j = 0; j < i; ++j) {
             sum += matrix[i][j] * matrix[num_cols - 1][j];
         }
         matrix[num_cols - 1][i] =
             (matrix[num_cols - 1][i] - sum) / matrix[i][i];
     }
  
     // back-substitution
     for (int i = num_rows - 1; i >= 0; --i) {
         double sum = 0.0;
         for (int j = i + 1; j < num_cols - 1; ++j) {
             sum += matrix[j][i] * matrix[i][num_cols - 1];
         }
         matrix[i][num_cols - 1] =
             (matrix[i][num_cols - 1] - sum) / matrix[i][i];
     }
 }

◆ solve_gauss()

std::vector< double > gpmp::linalg::LinSys::solve_gauss ( )

Solve the linear system using Gaussian Elimination.

Returns: vector containing solutinos to system

Definition at line 70 of file linsys.cpp.

                                                 {
     // temp matrix because our result is a vector of solutions
     std::vector<std::vector<double>> temp_mtx = matrix;
  
     for (int i = 0; i < num_rows - 1; ++i) {
         // perform partial pivot
         int pivotRow = i;
         for (int k = i + 1; k < num_rows; ++k) {
             if (std::abs(temp_mtx[k][i]) > std::abs(temp_mtx[pivotRow][i]))
                 pivotRow = k;
         }
         std::swap(temp_mtx[i], temp_mtx[pivotRow]);
  
         // eliminate entries below pivot
         for (int k = i + 1; k < num_rows; ++k) {
             double factor = temp_mtx[k][i] / temp_mtx[i][i];
             for (int j = i; j < num_cols; ++j) {
                 temp_mtx[k][j] -= factor * temp_mtx[i][j];
             }
         }
     }
  
     // back-substitution
     std::vector<double> solutions(num_rows, 0);
  
     for (int i = num_rows - 1; i >= 0; --i) {
         double sum = 0.0;
         for (int j = i + 1; j < num_cols - 1; ++j) {
             sum += temp_mtx[i][j] * solutions[j];
         }
         solutions[i] = (temp_mtx[i][num_cols - 1] - sum) / temp_mtx[i][i];
     }
  
     return solutions;
 }

◆ solve_jacobi()

void gpmp::linalg::LinSys::solve_jacobi	(	int	maxIterations = `100`,
		double	tolerance = `1e-10`
	)

Solve the linear system using Jacobi iteration.

Parameters

maxIterations	Maximum number of iterations (default: 100)
tolerance	Convergence tolerance (default: 1e-10)

Definition at line 263 of file linsys.cpp.

                                                                        {
     std::vector<std::vector<double>> x(num_rows, std::vector<double>(1, 0.0));
     std::vector<std::vector<double>> xOld = x;
  
     for (int iter = 0; iter < maxIterations; ++iter) {
         for (int i = 0; i < num_rows; ++i) {
             double sum = 0.0;
             for (int j = 0; j < num_cols - 1; ++j) {
                 if (j != i) {
                     sum += matrix[i][j] * xOld[j][0];
                 }
             }
             x[i][0] = (matrix[i][num_cols - 1] - sum) / matrix[i][i];
         }
  
         // check convergence
         double maxDiff = 0.0;
         for (int i = 0; i < num_rows; ++i) {
             double diff = std::abs(x[i][0] - xOld[i][0]);
             if (diff > maxDiff) {
                 maxDiff = diff;
             }
         }
  
         // converged
         if (maxDiff < tolerance) {
             break;
         }
  
         xOld = x;
     }
  
     // copy the solution back to the augmented matrix
     for (int i = 0; i < num_rows; ++i) {
         matrix[i][num_cols - 1] = x[i][0];
     }
 }

◆ solve_lu()

void gpmp::linalg::LinSys::solve_lu ( )

Solve the linear system using LU decomposition.

Definition at line 195 of file linsys.cpp.

                                 {
     lu_decomp();
  
     // forward substitution
     for (int i = 0; i < num_rows - 1; ++i) {
         for (int j = i + 1; j < num_rows; ++j) {
             matrix[j][num_cols - 1] -= matrix[j][i] * matrix[i][num_cols - 1];
         }
     }
  
     // back-substitution
     for (int i = num_rows - 1; i >= 0; --i) {
         for (int j = i + 1; j < num_cols - 1; ++j) {
             matrix[i][num_cols - 1] -= matrix[i][j] * matrix[j][num_cols - 1];
         }
         matrix[i][num_cols - 1] /= matrix[i][i];
     }
 }

Member Data Documentation

◆ matrix

std::vector<std::vector<double> > gpmp::linalg::LinSys::matrix

Definition at line 50 of file linsys.hpp.

Referenced by LinSys().

◆ num_cols

int gpmp::linalg::LinSys::num_cols

Definition at line 52 of file linsys.hpp.

Referenced by LinSys().

◆ num_rows

int gpmp::linalg::LinSys::num_rows

Definition at line 51 of file linsys.hpp.

Referenced by LinSys().

The documentation for this class was generated from the following files:

include/openGPMP/linalg/linsys.hpp
modules/linalg/linsys.cpp

Public Member Functions

Public Attributes

Detailed Description

Constructor & Destructor Documentation

◆ LinSys()

Member Function Documentation

◆ determinant()

◆ diagonally_dominant()

◆ display()

◆ display_mtx()

◆ frobenius_norm()

◆ gram_schmidt()

◆ inf_norm()

◆ invert_mtx()

◆ is_consistent()

◆ is_homogeneous()

◆ is_symmetric()

◆ lu_decomp()

◆ one_norm()

◆ solve_cholesky()

◆ solve_gauss()

◆ solve_jacobi()

◆ solve_lu()

Member Data Documentation

◆ matrix

◆ num_cols

◆ num_rows