Numerical libraries in Matrix Problems

11. Numerical libraries in Matrix Problems#

11.1. GSL : Gnu Scientific Library#

https://www.gnu.org/software/gsl/ Documentation: https://www.gnu.org/software/gsl/doc/html/index.html

11.1.1. Special functions example : Bessel function#

See: https://www.gnu.org/software/gsl/doc/html/specfunc.html

// Adapted to c++ from : https://www.gnu.org/software/gsl/doc/html/specfunc.html 
// compile as: g++ -std=c++17 name.cpp -lgsl -lgslcblas
#include <cstdio>
#include <cmath>
#include <gsl/gsl_sf_bessel.h>

int main (void)
{
    double x = 5.0;
    double expected = -0.17759677131433830434739701;
    double fromcmath = std::cyl_bessel_j(0, x);

    double y = gsl_sf_bessel_J0 (x);

    std::printf ("J0(5.0)     = %.18f\n", y);
    std::printf ("exact       = %.18f\n", expected);
    std::printf ("fromcmath   = %.18f\n", fromcmath);
    return 0;
}

Exercise: Write a program that prints the Bessel function of order zero for \(x \in [0, 10]\) in steps \(\Delta x = 0.01\). Plot it.

Exercise: Print and plot some of the Hypergeometric functions.
Exercise: What about linear algebra? check the examples.

11.2. Standard library for special functions#

Check the standard library for implementatios of many special functions. Must compile with -std=c++17

11.3. Eigen : Linear Algebra#

http://eigen.tuxfamily.org/index.php?title=Main_Page Start by: http://eigen.tuxfamily.org/dox/GettingStarted.html

Basic matrix

#include <iostream>
#include <eigen3/Eigen/Dense>

int main(void)
{
    Eigen::MatrixXd m(2, 2);
    m(0, 1) = -9.877;
    std::cout << m << std::endl;

    return 0;
}


0	-9.877
0	0

Exercise: Create a 8x8 matrix and print it to the screen.
random matrix

#include <iostream>
#include <cstdlib>
#include <eigen3/Eigen/Dense>

int main(int argc, char **argv)
{
    //const int SEED = std::atoi(argv[1]);
    //srand(SEED); // control seed
    Eigen::MatrixXd m = Eigen::MatrixXd::Random(5,5);
    std::cout << m << std::endl;

    return 0;
}


0.680375	-0.604897	-0.0452059	0.83239	-0.967399
-0.211234	-0.329554	0.257742	0.271423	-0.514226
0.566198	0.536459	-0.270431	0.434594	-0.725537
0.59688	-0.444451	0.0268018	-0.716795	0.608354
0.823295	0.10794	0.904459	0.213938	-0.686642

Exercise: Create a 8x8 random matrix and print it to the screen. Control the seed, print it twice to check that the elements are the same.
First example from the tutorial:

// First example on how to use eigen
// compile with: g++ -std=c++11 -I /usr/include/eigen3
#include <iostream>
#include <eigen3/Eigen/Dense>

int main()
{
  // X = Dynamic size
  // d = double
  Eigen::MatrixXd m(2,2);
  m(0,0) = 3;
  m(1,0) = 2.5;
  m(0,1) = -1;
  m(1,1) = m(1,0) + m(0,1);
  std::cout << m << std::endl;
}


3	-1
2.5	1.5

Second example from the tutorial: Matrix-vector multiplication

#include <iostream>
#include <eigen3/Eigen/Dense>

int main()
{
    const int N = 3;
    srand(10);
    Eigen::MatrixXd m = Eigen::MatrixXd::Random(N, N);
    m = (m + Eigen::MatrixXd::Constant(N, N, 1.2)) * 50;
    std::cout << "m =" << std::endl << m << std::endl;
    Eigen::VectorXd v(N);
    v << 1, 2, 3;
    std::cout << "m * v =" << std::endl << m * v << std::endl;
}

11.3.1. Linear Algebra#

In this example we are going to use eigen to solve the system \(\boldsymbol{A}\vec x = \vec{x}\) using HouseHolder QR decomposition as done in http://eigen.tuxfamily.org/dox/group__TutorialLinearAlgebra.html .

#include <iostream>
#include <eigen3/Eigen/Dense>

int main()
{
    Eigen::Matrix3f A;
    Eigen::Vector3f b;
    A << 1,2,3,  4,5,6,  7,8,10;
    b << 3, 3, 4;
    std::cout << "Here is the matrix A:\n" << A << std::endl;
    std::cout << "Here is the vector b:\n" << b << std::endl;
    Eigen::Vector3f x = A.colPivHouseholderQr().solve(b);
    std::cout << "The solution is:\n" << x << std::endl;
    std::cout << "Verification:\n" << (A*x - b).norm() << std::endl;
}


Here	is	the	matrix	A:
1	2	3
4	5	6
7	8	10
Here	is	the	vector	b:
3
3
4
The	solution	is:
-2
0.999997
1
Verification:
1.43051e-06

With double :

#include <iostream>
#include <eigen3/Eigen/Dense>

int main()
{
  Eigen::Matrix3d A;
  Eigen::Vector3d b;
  A << 1,2,3,  4,5,6,  7,8,10;
  b << 3, 3, 4;
  std::cout << "Here is the matrix A:\n" << A << std::endl;
  std::cout << "Here is the vector b:\n" << b << std::endl;
  Eigen::Vector3d x = A.colPivHouseholderQr().solve(b);
  std::cout << "The solution is:\n" << x << std::endl;
  std::cout << "Verification:\n" << (A*x - b).norm() << std::endl;
}


Here	is	the	matrix	A:
1	2	3
4	5	6
7	8	10
Here	is	the	vector	b:
3
3
4
The	solution	is:
-2
1
1
Verification:
3.58036e-15

Exercise: Create a random 40x40 matrix \(A\), a random vector \(b\) of size 40, and solve the system \(Ax = b\) using FullPivHouseholderQR. Compare the execution time with a 100x100 problem size.
Exercise: Compile the previous example with profiling, run under valgrind cachegrind, and check the output. Can you change or optimize any of that?

Now use eigen to compute the eigenvalues of a given matrix:

#include <iostream>
#include <eigen3/Eigen/Dense>

int main()
{
    Eigen::Matrix2d A;
    A << 1, 2, 2, 3;
    std::cout << "Here is the matrix A:\n" << A << std::endl;
    Eigen::SelfAdjointEigenSolver<Eigen::Matrix2d> eigensolver(A);
    std::cout << "The eigenvalues of A are:\n" << eigensolver.eigenvalues() << std::endl;
    std::cout << "Here's a matrix whose columns are eigenvectors of A \n"
              << "corresponding to these eigenvalues:\n"
              << eigensolver.eigenvectors() << std::endl;
}

Exercise : Compute the eigen-values of a random matrix of size 20x20.