/*
https://libeigen.gitlab.io/eigen/docs-5.0/GettingStarted.html

ADD to CMakeLists.txt:

target_include_directories(Lecture10 PRIVATE /home/zqex/tmp/QT/3rdparty/Eigen)
or
target_include_directories(Lecture10 PRIVATE ${PROJECT_SOURCE_DIR}/../3rdparty/Eigen)

Optionally:
add_compile_options(-std=c++11) #optional to use c++11

Class:

Matrix<number type, rows, columns>

Number type = int, float, double, complex ..
rows = number or Dynamic
cols = number or Dynamic

Number: then matrix size is known at compile time. Dynamic it is
determined during the execution of the program.

Useful typedefs:
typedef Matrix<double, Dynamic, Dynamic> MatrixXd;
typedef Matrix<double, 3, 1> Vector3d;
typedef Matrix<double, 1, 3> RowVector3d;
*/

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

using namespace std;

int main() {
// Define Eigen matrix:
    Eigen::MatrixXd A(2, 2);   // 2x2 matrix   0..1 x 0..1  of doubles

// Initialize matrix element by element:
    A(0, 0) =  3.0;
    A(1, 0) =  2.5;
    A(0, 1) = -1.0;
    A(1, 1) =  2.0;

// A second matrix
    Eigen::MatrixXd B(2, 2);
    B << 1.0 , 2.0 , 3.0, 4.0;     // Special syntax for initializing.

// Print them out:
    cout << "Matrix: A=" << endl
         << A << endl << endl;

    cout << "Matrix: B=" << endl
         << B << endl << endl;

// Vectors
    Eigen::Vector2d v(1.0, 2.0);
    cout << "Matrix*Vector = " << endl
         << B*v << endl << endl;

// Matrix Aritmetrics:
    cout << "Matrix: A*B=" << endl
         << A*B << endl << endl;

    Eigen::MatrixXd C=A*B+B+2*A;
    cout << "Matrix: A*B+B+2A=" << endl
         << C << endl << endl;

// More advanced determinant and inverse:
    cout << "det(C)=" << C.determinant() << endl;
    cout << "The inverse of C is:" << endl
         << C.inverse() << endl << endl;

    cout << "Proof: C*C^-1 =" << endl
         << C*C.inverse() << endl << endl;
    cout << "Proof: C^-1*C =" << endl
         << C.inverse()*C << endl << endl;


// Obtaining the eigenvalues and eigenvectors
    Eigen::EigenSolver<Eigen::MatrixXd> solv(C);
    if (solv.info() != Eigen::Success) abort();
    Eigen::MatrixXcd EVec=solv.eigenvectors();
    Eigen::VectorXcd EVal=solv.eigenvalues();

    cout << "The eigenvalues of A are:\n" << EVal << endl << endl;
    cout << "The eigen vectors are: " << endl
              << EVec << endl << endl;
}