linear_algebra.hh

#ifndef VTESTBED_LINALG_LINEAR_ALGEBRA_HH
#define VTESTBED_LINALG_LINEAR_ALGEBRA_HH

#include <stdexcept>

#include <vtestbed/base/blitz.hh>
#include <vtestbed/linalg/types.hh>

namespace vtb {

    namespace linalg {

        class Invalid_matrix: public std::runtime_error {

        public:

            Invalid_matrix() = default;

            inline explicit
            Invalid_matrix(const char* msg): std::runtime_error(msg) {}

        };

        inline float conj(float x) { return x; }
        inline float conj(double x) { return x; }

        template <class T>
        inline std::complex<T>
        conj(const std::complex<T>& z) {
            return std::conj(z);
        }

        template <class T>
        class Vector: public blitz::Array<T,1> {

        private:
            typedef blitz::Array<T,1> base_type;

        public:
            using base_type::base_type;
            using base_type::operator=;
            using base_type::operator();

            inline
            Vector(base_type rhs):
            base_type(rhs) {}

        };

        template <class T>
        class Matrix: public blitz::Array<T,2> {

        private:
            typedef blitz::Array<T,2> base_type;

        public:
            using base_type::base_type;
            using base_type::operator=;

            inline
            Matrix(base_type rhs):
            base_type(rhs) {}

            inline bool
            is_square() const {
                return this->rows() == this->cols();
            }

            bool
            is_symmetric() const;

            bool
            is_symmetric(T eps) const;

            bool
            is_positive_definite() const;

            bool
            is_toeplitz() const;

            bool
            is_toeplitz(T eps) const;

            inline Square_matrix<T>
            inverse() const;

            inline Square_matrix<T>&
            inverse_self();

            inline Matrix<T>
            transpose() const {
                Matrix<T> tmp(*this);
                tmp.base_type::transpose(1, 0);
                return tmp;
            }

            inline void
            transpose_self() {
                this->base_type::transposeSelf(1, 0);
            }

        };

        template <class T>
        class Square_matrix: public Matrix<T> {

        public:
            using Matrix<T>::Matrix;
            using Matrix<T>::operator=;
            using Matrix<T>::operator();

        };

        template <class T>
        class Lower_triangular_matrix: public Square_matrix<T> {

        private:
            typedef Square_matrix<T> base_type;

        public:
            using base_type::base_type;
            using base_type::operator=;
            using base_type::operator();

            Vector<T>
            forward_substitution(const Vector<T>& b) const;

            Vector<T>
            backward_substitution(const Vector<T>& y) const;

            inline Vector<T>
            solve(const Vector<T>& b) const {
                return backward_substitution(forward_substitution(b));
            }

        };

        template <class T>
        class Lower_triangular_matrix_ldlt: public Lower_triangular_matrix<T> {

        private:
            typedef Lower_triangular_matrix<T> base_type;

            friend class Symmetric_matrix<T>;

        private:
            Vector<T> _diagonal;

        public:
            using base_type::base_type;
            using base_type::operator=;
            using base_type::operator();

            inline Vector<T>
            solve(const Vector<T>& b) const {
                Vector<T> y(this->forward_substitution(b) / this->_diagonal);
                return this->backward_substitution(y);
            }

        };

        template <class T>
        class Symmetric_matrix: public Square_matrix<T> {

        private:
            typedef Square_matrix<T> base_type;

        public:
            using base_type::base_type;
            using base_type::operator=;
            using base_type::operator();

            Lower_triangular_matrix_ldlt<T>
            cholesky_indefinite() const;

            Vector<T>
            eigen_values(T eps, int nsweeps) const;

        };

        template <class T>
        class Positive_definite_matrix: public Symmetric_matrix<T> {

        private:
            typedef Symmetric_matrix<T> base_type;

        public:
            using base_type::base_type;
            using base_type::operator=;
            using base_type::operator();

            Lower_triangular_matrix<T>
            cholesky() const;

        };

        template <class T>
        class Lower_upper_triangular_matrix: public Square_matrix<T> {

        private:
            typedef Square_matrix<T> base_type;
            typedef Vector<int> int_array_type;

        private:
            int_array_type _permutations;
            int _npermutations = 0;

        public:
            using base_type::base_type;
            using base_type::operator=;
            using base_type::operator();

            inline const int_array_type&
            permutations() const noexcept {
                return this->_permutations;
            }

            Square_matrix<T>
            inverse() const;

            inline Square_matrix<T>&
            inverse_self() {
                Square_matrix<T> m = this->inverse();
                this->reference(m);
                return *this;
            }

            Vector<T>
            solve(const Vector<T>& b) const;

            Lower_upper_triangular_matrix<T>&
            lower_upper_triangular_self();

            T
            determinant() const;

        };

        template <class T>
        inline Lower_upper_triangular_matrix<T>
        lower_upper_triangular(const Matrix<T>& rhs) {
            Lower_upper_triangular_matrix<T> tmp(rhs.shape());
            tmp = rhs;
            tmp.lower_upper_triangular_self();
            return tmp;
        }

        template <class NewMatrix, class OldMatrix>
        inline const NewMatrix&
        matrix_cast(const OldMatrix& rhs) {
            return static_cast<const NewMatrix&>(rhs);
        }

        template <class NewMatrix, class OldMatrix>
        inline NewMatrix&
        matrix_cast(OldMatrix& rhs) {
            return static_cast<NewMatrix&>(rhs);
        }

        template <class T>
        Vector<T>
        product(const Matrix<T>& lhs, const Vector<T>& rhs);

        template <class T>
        inline Vector<T>
        operator*(const Matrix<T>& lhs, const Vector<T>& rhs) {
            return product<T>(lhs, rhs);
        }

        template <class T>
        Matrix<T>
        product(const Matrix<T>& lhs, const Matrix<T>& rhs);

        template <class T>
        inline Matrix<T>
        operator*(const Matrix<T>& lhs, const Matrix<T>& rhs) {
            return product<T>(lhs, rhs);
        }

        template <class T>
        inline T
        det(const Lower_upper_triangular_matrix<T>& rhs) {
            return rhs.determinant();
        }

        template <class T>
        inline T
        det(const Square_matrix<T>& rhs) {
            return det(lower_upper_triangular(rhs));
        }

        template <class T>
        inline T
        inverse(const Lower_upper_triangular_matrix<T>& rhs) {
            return rhs.inverse();
        }

        template <class T>
        inline T
        inverse(const Square_matrix<T>& rhs) {
            return inverse(lower_upper_triangular(rhs));
        }

        template <class T>
        inline T
        trace(const Square_matrix<T>& rhs) {
            T sum{0};
            const int n = rhs.nrows();
            for (int i=0; i<n; ++i) {
                sum += rhs(i,i);
            }
            return sum;
        }

        template <class T>
        inline Square_matrix<T>
        Matrix<T>::inverse() const {
            return lower_upper_triangular(*this).inverse();
        }

        template <class T>
        inline Square_matrix<T>&
        Matrix<T>::inverse_self() {
            this->reference(lower_upper_triangular(*this).inverse());
            return matrix_cast<Square_matrix<T>>(*this);
        }

        template <class T>
        inline Vector<T>
        solve(
            const Lower_upper_triangular_matrix<T>& lhs,
            const Vector<T>& rhs
        ) {
            return lhs.solve(rhs);
        }

        template <class T>
        inline Vector<T>
        solve(const Square_matrix<T>& lhs, const Vector<T>& rhs) {
            return solve(lower_upper_triangular(lhs), rhs);
        }

        template <class T>
        inline Vector<T>
        solve(const Symmetric_matrix<T>& lhs, const Vector<T>& rhs) {
            return lhs.cholesky_indefinite().solve(rhs);
        }

        template <class T>
        inline Vector<T>
        solve(const Positive_definite_matrix<T>& lhs, const Vector<T>& rhs) {
            return lhs.cholesky().solve(rhs);
        }

        template <class T>
        inline Vector<T>
        eigen_values(const Symmetric_matrix<T>& rhs, T eps, int nsweeps) {
            return rhs.eigen_values(eps, nsweeps);
        }

    }

}

#endif // vim:filetype=cpp