Runge-Kutta-method/matrix.h at master · user959/Runge-Kutta-method · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
#include <vector>
#include <stdexcept>
#include <iostream>
#include <cmath>

using std::vector;

template <typename T>
vector<T>& operator+=(vector<T>& v1, const vector<T>& v2) {
    if (v1.size() != v2.size())
        throw std::out_of_range{"operator+=: sizes do not match"};

    for (size_t i = 0; i < v1.size(); ++i)
        v1[i] += v2[i];
    return v1;
}

template <typename T>
vector<T> operator+(vector<T> v1, const vector<T>& v2) {
    v1 += v2;
    return v1;
}

template <typename T>
vector<T>& operator-=(vector<T>& v1, const vector<T>& v2) {
    if (v1.size() != v2.size())
        throw std::out_of_range{"operator-=: sizes do not match"};

    for (size_t i = 0; i < v1.size(); ++i)
        v1[i] -= v2[i];
    return v1;
}

template <typename T>
vector<T> operator-(vector<T> v1, const vector<T>& v2) {
    v1 -= v2;
    return v1;
}

template <typename T>
vector<T>& operator*=(vector<T>& v, const T scalar) {
    for (size_t i = 0; i < v.size(); ++i)
        v[i] *= scalar;
    return v;
}

template <typename T>
vector<T>& operator*=(const T scalar, vector<T>& v) {
    return v * scalar;
}

template <typename T>
vector<T> operator*(const T scalar, vector<T> v) {
    v *= scalar;
    return v;
}

// template <typename T>
// vector<vector<T>>& operator*=(vector<vector<T>>& v, vector<vector<T>> v2) {
// }

template <typename T>
vector<vector<T>> operator*(const vector<vector<T>>& lhs, const vector<vector<T>>& rhs) {
    if (lhs.size() && lhs[0].size() != rhs.size())
        throw "operator*: matrices' sizes do not match";
    if (lhs.empty() || lhs[0].empty() || rhs.empty() || rhs[0].empty())
        throw "operator*: degenerate matrix";
    vector<vector<T>> ans(lhs.size(), vector<T>(rhs[0].size(), 0));
    for (size_t i = 0; i < lhs.size(); ++i) {
        for (size_t j = 0; j < rhs[0].size(); ++j) {
            for (size_t k = 0; k < rhs.size(); ++k) {
                ans[i][j] += lhs[i][k] * rhs[k][j];
            }
        }
    }
    return ans;
}

template <typename T>
double euclideanNorm(const vector<T>& v) {
    double ans = 0;
    for (size_t i = 0; i < v.size(); ++i) {
        ans += v[i] * v[i];
    }
    ans /= v.size(); // yes, I know that there isn't such thing in euclidean norm
    return sqrt(ans);
}

template <typename T>
double manhettenNorm(const vector<T>& v) {
    double ans = 0;
    for (size_t i = 0; i < v.size(); ++i) {
        ans += abs(v[i]);
    }
    return ans;
}

template <typename T>
vector<T> operator-(const vector<T>& v) {
    return -1.0 * v;
}

template <typename T>
std::ostream& operator<<(std::ostream& os, const vector<T>& v) {
    os << "{";
    for (size_t i = 0; i < v.size(); ++i) {
        os << v[i];
        if (i + 1 != v.size())
            os << ", ";
    }
    os << "}";
    return os;
}

template <typename T>
std::ostream& operator<<(std::ostream& os, const vector<vector<T>>& v) {
    for (size_t i = 0; i < v.size(); ++i) {
        for (size_t j = 0; j < v[i].size(); ++j) {
            os << v[i][j] << "\t";
        }
        os << "\n";
    }
    return os;
}