gauss.h

Go to the documentation of this file.

#ifndef CHEMMISOL_GAUSS_H
#define CHEMMISOL_GAUSS_H
 
#include <iostream>
#include <array>
#include <complex>
#include "../logging.h"
 
namespace chemmisol { namespace gauss {
 
    template<typename M, typename X>
        X solve(const M& m, const X& y) {
            CHEM_LOG(TRACE) << "[GAUSS START]";
            auto _m = augment(m, y);
            std::size_t n = _m.size();
            CHEM_LOG(TRACE) << "Step 0:";
            for(std::size_t i = 0; i < n; i++) {
                CHEM_LOG(TRACE) << _m[i];
            }
 
            for(std::size_t i = 0; i < n; i++) {
                for(std::size_t j = i+1; j < n; j++) {
                    auto ratio = _m[j][i]/_m[i][i];
                    for(std::size_t k = 0; k < _m[i].size(); k++) {
                        _m[j][k] = _m[j][k] - ratio*_m[i][k];
                    }
                }
            }
 
            CHEM_LOG(TRACE) << "Step 1:";
            for(std::size_t i = 0; i < n; i++) {
                CHEM_LOG(TRACE) << _m[i];
            }
            
            X x = y;
            x[n-1] = _m[n-1][n]/_m[n-1][n-1];
 
            for(int i = n-2; i>=0; i--) {
                x[i] = _m[i][n];
                for(std::size_t j = i+1; j < n; j++) {
                    x[i] = x[i] - _m[i][j]*x[j];
                }
                x[i] = x[i]/_m[i][i];
            }
            CHEM_LOG(TRACE) << "Step 2:";
            for(std::size_t i = 0; i < x.size(); i++)
                CHEM_LOG(TRACE) << "x[" << i << "]=" << x[i];
            CHEM_LOG(TRACE) << "[GAUSS END]";
            return x;
        }
 
    template<>
        double solve<double, double>(const double& f, const double& y);
 
    template<>
        std::complex<double> solve<std::complex<double>, std::complex<double>>(
                const std::complex<double>& f, const std::complex<double>& y);
 
}}
#endif