linear.h

Go to the documentation of this file.

#ifndef CHEMMISOL_LINEAR_H
#define CHEMMISOL_LINEAR_H
 
#include <array>
#include <vector>
#include <initializer_list>
#include <iostream>
#include <cmath>
#include <complex>
#include "../logging.h"
 
namespace chemmisol {
    /*
     * Fixed size array linear algebra
     */
 
    template<typename T, std::size_t N>
        using X = std::array<T, N>;
 
    template<typename T, std::size_t N, std::size_t P = N>
        using M = std::array<std::array<T, P>, N>;
 
    template<typename T>
        using VecX = std::vector<T>;
    template<typename T>
        using VecM = std::vector<std::vector<T>>;
 
    template<typename T, std::size_t N, std::size_t P>
        X<T, N> operator*(const M<T, N, P>& m, const X<T, P>& x) {
            X<T, N> x1;
            for(std::size_t i = 0; i < N; i++) {
                x1[i] = 0;
                for(std::size_t j = 0; j < P; j++) {
                    x1[i] += m[i][j] * x[j];
                }
            }
            return x1;
        }
 
    template<typename T, std::size_t N>
        X<T, N> operator-(const X<T, N>& x) {
            X<T, N> x1;
            for(std::size_t i = 0; i < N; i++)
                x1[i] = -x[i];
            return x1;
        }
 
    template<typename T, std::size_t N>
        X<T, N> operator+(const X<T, N>& x1, const X<T, N>& x2) {
            X<T, N> x3;
            for(std::size_t i = 0; i < N; i++)
                x3[i] = x1[i]+x2[i];
            return x3;
        }
 
    template<typename T, std::size_t N>
        X<T, N> operator*(const T& a, const X<T, N>& x) {
            X<T, N> result;
            for(std::size_t i = 0; i < N; i++)
                result[i] = a*x[i];
            return result;
        }
 
    template<typename T, std::size_t N, std::size_t P>
        M<T, N, P> operator*(const T& a, const M<T, N, P>& m) {
            M<T, N, P> result;
            for(std::size_t i = 0; i < N; i++) {
                for(std::size_t j = 0; j < P; j++) {
                    result[i][j] = a * m[i][j];
                }
            }
            return result;
        }
    template<typename T, std::size_t N>
        double norm(const X<T, N>& x) {
            T a = 0;
            for(auto& v : x) {
                a+=std::pow(v, 2);
            }
            return std::sqrt(a);
        }
 
    template<typename T, std::size_t N>
        double norm(const X<std::complex<T>, N>& x) {
            T a = 0;
            for(auto& v : x) {
                a+=std::pow(std::norm(v), 2);
            }
            return std::sqrt(a);
        }
 
    template<typename T, std::size_t N>
        X<T, N> abs(const X<T, N>& x) {
            X<T, N> abs_x;
            for(std::size_t i = 0; i < N; i++)
                abs_x[i] = std::abs(x[i]);
            return abs_x;
        }
 
    template<typename _M>
        M<
            typename _M::value_type::value_type,
            std::tuple_size<_M>::value,
            std::tuple_size<typename _M::value_type>::value+1>
        augment(
                const _M& m,
                const X<typename _M::value_type::value_type, std::tuple_size<_M>::value>& x
                ) {
            M<
                typename _M::value_type::value_type,
                std::tuple_size<_M>::value,
                std::tuple_size<typename _M::value_type>::value+1> a;
                for(std::size_t i = 0; i < std::tuple_size<_M>::value; i++) {
                    for(std::size_t j = 0; j < std::tuple_size<typename _M::value_type>::value; j++) {
                        a[i][j] = m[i][j];
                    }
                    a[i][std::tuple_size<_M>::value] = x[i];
                }
                return a;
        };
 
    template<typename T, std::size_t N>
        std::ostream& operator<<(std::ostream& o, const X<T, N>& x) {
            o << "[";
            for(std::size_t i = 0; i < N-1; i++)
                o << x[i] << ", ";
            if(N > 0)
                o << x[N-1];
            o << "]";
            return o;
        }
 
    /*
     * Vector linear algebra
     */
 
    template<typename T>
        VecX<T> operator*(const VecM<T>& m, const VecX<T>& x) {
            VecX<T> x1(x.size());
            for(std::size_t i = 0; i < x.size(); i++) {
                x1[i] = 0;
                for(std::size_t j = 0; j < m[i].size(); j++) {
                    x1[i] += m[i][j] * x[j];
                }
            }
            return x1;
        }
 
    template<typename T>
        VecX<T> operator-(const VecX<T>& x) {
            VecX<T> x1(x.size());
            for(std::size_t i = 0; i < x.size(); i++)
                x1[i] = -x[i];
            return x1;
        }
 
    template<typename T>
        VecX<T> operator+(const VecX<T>& x1, const VecX<T>& x2) {
            VecX<T> x3(x1.size());
            for(std::size_t i = 0; i < x1.size(); i++)
                x3[i] = x1[i]+x2[i];
            return x3;
        }
 
    template<typename T>
        double norm(const VecX<T>& x) {
            double a = 0;
            for(auto& v : x) {
                a+=std::pow(v, 2);
            }
            return std::sqrt(a);
        }
 
    template<typename T>
        VecX<T> abs(const VecX<T>& x) {
            VecX<T> abs_x(x);
            for(std::size_t i = 0; i < x.size(); i++)
                abs_x[i] = std::abs(x[i]);
            return abs_x;
        }
 
    template<typename T>
        VecM<T>
        augment(const VecM<T>& m, const VecX<T>& x) {
            VecM<T> a(m.size());
            for(std::size_t i = 0; i < m.size(); i++) {
                a[i].resize(m[i].size()+1);
                for(std::size_t j = 0; j < m[i].size(); j++) {
                    a[i][j] = m[i][j];
                }
                a[i][m[i].size()] = x[i];
            }
            return a;
        };
 
    template<typename T>
        std::ostream& operator<<(std::ostream& o, const VecX<T>& x) {
            o << "[";
            for(std::size_t i = 0; i < x.size()-1; i++)
                o << x[i] << ", ";
            if(x.size() > 0)
                o << x[x.size()-1];
            o << "]";
            return o;
        }
 
    /*
     * Scalar values algebra.
     */
 
    template<typename T>
        double norm(const T& x) {
            return std::abs(x);
        }
 
    template<typename T>
        T abs(const T& x) {
            return std::abs(x);
        }
}
#endif