linalg/matrix.cpp at main · sirohikartik/linalg · 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
#include "matrix.h"

using namespace std;


void test_determinant(){
    auto approx = [](double a, double b){ return abs(a-b) < 1e-9; };

    vector<vector<double>> v1 = {{1,2},{3,4}};
    M m1(v1);
    assert(approx(m1.det(0,1,0,1), -2.0));

    vector<vector<double>> v2 = {{1,0,0},{0,1,0},{0,0,1}};
    M m2(v2);
    assert(approx(m2.det(0,2,0,2), 1.0));

    vector<vector<double>> v3 = {{6,1,1},{4,-2,5},{2,8,7}};
    M m3(v3);
    assert(approx(m3.det(0,2,0,2), -306.0));
    assert(approx(m3.det(0,1,0,1), -16.0));

    vector<vector<double>> v4 = {{0,1},{1,0}};

    M m4(v4);
    assert(approx(m4.det(),-1.0));

    cout << "\033[1;32m [5/5] passed \033[0m" << endl;
    //cout << " [4/4] passed" << endl;
}


void test_operators(){
    auto approx = [](double a, double b){ return abs(a-b) < 1e-9; };
    auto mat_approx = [&](M a, M b, int rows, int cols){
        for(int i=0;i<rows;i++)
            for(int j=0;j<cols;j++)
                if(!approx(a.det(j,j,i,i), b.det(j,j,i,i))) return false;
        return true;
    };

    vector<vector<double>> v1 = {{1,2},{3,4}};
    vector<vector<double>> v2 = {{5,6},{7,8}};
    M m1(v1), m2(v2);

    // test +
    vector<vector<double>> exp_add = {{6,8},{10,12}};
    M m_add = m1 + m2;
    assert(approx(m_add.det(), M(exp_add).det()));

    // test -
    vector<vector<double>> exp_sub = {{-4,-4},{-4,-4}};
    M m_sub = m1 - m2;
    assert(approx(m_sub.det(), M(exp_sub).det()));

    // test * scalar (M * lambda)
    vector<vector<double>> exp_scale = {{2,4},{6,8}};
    M m_scale = m1 * 2.0;
    assert(approx(m_scale.det(), M(exp_scale).det()));

    // test * scalar (lambda * M)
    M m_scale2 = 2.0 * m1;
    assert(approx(m_scale2.det(), M(exp_scale).det()));

    // test * matrix
    // [1,2][3,4] * [5,6][7,8] = [19,22][43,50]
    vector<vector<double>> exp_mul = {{19,22},{43,50}};
    M m_mul = m1 * m2;
    assert(approx(m_mul.det(), M(exp_mul).det()));

    // test size mismatch throws
    vector<vector<double>> v3 = {{1,2,3},{4,5,6}};
    M m3(v3);
    bool threw = false;
    try { M bad = m1 + m3; }
    catch(runtime_error&) { threw = true; }
    assert(threw);

    cout << "\033[1;32m [6/6] passed \033[0m" << endl;
}

void test_inverse(){
    auto approx = [](double a, double b){ return abs(a-b) < 1e-9; };

    // A * A_inv should equal identity
    vector<vector<double>> v = {{1,2},{3,4}};
    M m(v);
    M m_inv = m.inv();
    M should_be_I = m * m_inv;

    assert(approx(should_be_I.det(), 1.0));  // det of identity = 1

    // singular matrix should throw
    bool threw = false;
    vector<vector<double>> singular = {{1,2},{2,4}};
    M ms(singular);
    try { ms.inv(); }
    catch(runtime_error&) { threw = true; }
    assert(threw);

    cout << "\033[1;32m [2/2] inverse passed \033[0m" << endl;
}
int main() {
	test_determinant();
	test_operators();
	test_inverse();
	return 0;
}