#include <Eigen/Dense>

// fixed-size matrix (rows, cols, type)
Eigen::Matrix<float, 3, 3> A;          // 3x3 float matrix
Eigen::Matrix<double, 2, 4> B;         // 2x4 double matrix

// dynamic-size matrix
Eigen::MatrixXd C;                     // double, dynamic size
Eigen::MatrixXf D;                     // float, dynamic size

// vector types (shortcuts)
Eigen::Vector3d v;                     // 3x1 double vector
Eigen::RowVectorXf r;                  // 1xN float row vector

// constructors
Eigen::Matrix3d I = Eigen::Matrix3d::Identity();    // identity
Eigen::MatrixXd Z = Eigen::MatrixXd::Zero(3, 3);    // zero
Eigen::MatrixXd O = Eigen::MatrixXd::Ones(2, 5);    // ones
Eigen::MatrixXd R = Eigen::MatrixXd::Random(4, 4);  // random values

// direct initialization
Eigen::Matrix3d M;
M << 1, 2, 3,
     4, 5, 6,
     7, 8, 9;
     
// map (wrap existing data)
double array[6] = {1, 2, 3, 4, 5, 6};
Eigen::Map<Eigen::Matrix<double, 2, 3>> mat(array); // use array as 2x3 matrix
Eigen::Map<Eigen::MatrixXd> mat(vec.data(), h, w); // convert vector to h*w mat

// resize (dynamic only)
C.resize(5, 5);

// access
M(i, j);        // element at row i, column j (0-based)
M.row(i);       // i-th row (as a vector)
M.col(j);       // j-th column (as a vector)
v(i);           // vector element
v.head(i)       // first i elements (from a vector)
v.segment(1, 3) // 3 elements, starting at index 1 (from a vector)
v.tail(i)       // last i elements (from a vector)

// block operations
M.block(i, j, p, q);  // block of size p x q starting at (i, j)
M.topLeftCorner(p, q);
M.bottomRightCorner(p, q);
M.leftCols(i); // first i cols
M.rightCols(j); // last j cols
M.row(i);
M.col(j);

M.diagonal(-1) // first sub-diagonal
M.diagonal<-1>() // templating is normally faster
M.triangularView<Eigen::Upper>();
M.triangularView<Eigen::StrictlyLower>(); // without diagonal

// operations
A.transpose();        // transpose
A.conjugate();        // complex conjugate
A.adjoint();          // conjugate transpose
A.inverse();          // matrix inverse
A.determinant();      // determinant
A.trace();            // trace

// functions
A.setZero() // fill matrix with zeroes
A.setConstant(c) // fill matrix with c

// reduction
A.size();
A.sum();
A.prod();
A.minCoeff(); // min entry
A.maxCoeff();
A.norm(); // frobenius norm
A.mean();
// lu decomposition
Eigen::FullPivLU<Eigen::MatrixXd> lu(A);

// qr decomposition
Eigen::HouseholderQR<Eigen::MatrixXd> qr(A);

// eigenvalues
Eigen::EigenSolver<Eigen::MatrixXd> es(A);

// singular value decomposition
Eigen::JacobiSVD<Eigen::MatrixXd> svd(A, Eigen::ComputeFullU | Eigen::ComputeFullV);

// apply operation to whole vector
Eigen::VectorXd b = y.array().log().matrix();

X = A.lu().solve(B)

// DFT (discrete fourier transform)
Eigen::FFT<double> fft;
fft.inv(((fft.fwd(u)).cwiseProduct(fft.fwd(x))).eval())

//or, if complex:
(fft.inv(
	fft.fwd(u).cwiseProduct(fft.fwd(x))
)).real()

// periodic convolution
pconvfft(vec1, vec2)
# least squares
Eigen::VectorXd x = Eigen::VectorXd::Zero(n);
x = (A.transpose() * A).lu().solve(A.transpose() * b);
O(mn2multiplication+n3lu+mnmatvec+n2solve)=O(mn2+n3)=O(m)

Normally not recommended for least squares (amplifies conditioning issue)