|
3 | 3 |
|
4 | 4 | using namespace cpp4r; |
5 | 5 |
|
6 | | -[[cpp4r::register]] doubles_matrix<> invert_matrix_(doubles_matrix<> A) |
7 | | -{ |
8 | | - // Obtain (X'X)^-1 via Gauss-Jordan |
| 6 | +[[cpp4r::register]] doubles_matrix<> invert_matrix_(doubles_matrix<> A) { |
| 7 | + // Obtain (X'X)^-1 via Gauss-Jordan |
9 | 8 |
|
10 | | - // Check dimensions |
| 9 | + // Check dimensions |
11 | 10 |
|
12 | | - int N = A.nrow(); |
13 | | - int M = A.ncol(); |
| 11 | + int N = A.nrow(); |
| 12 | + int M = A.ncol(); |
14 | 13 |
|
15 | | - if (N != M) |
16 | | - { |
17 | | - stop("X must be a square matrix"); |
18 | | - } |
| 14 | + if (N != M) { |
| 15 | + stop("X must be a square matrix"); |
| 16 | + } |
19 | 17 |
|
20 | | - // Copy the matrix |
| 18 | + // Copy the matrix |
21 | 19 |
|
22 | | - writable::doubles_matrix<> Acopy(N, N); |
| 20 | + writable::doubles_matrix<> Acopy(N, N); |
23 | 21 |
|
24 | | - for (int i = 0; i < N; i++) |
25 | | - { |
26 | | - for (int j = 0; j < N; j++) |
27 | | - { |
28 | | - Acopy(i, j) = A(i, j); |
29 | | - } |
| 22 | + for (int i = 0; i < N; i++) { |
| 23 | + for (int j = 0; j < N; j++) { |
| 24 | + Acopy(i, j) = A(i, j); |
30 | 25 | } |
| 26 | + } |
31 | 27 |
|
32 | | - // Create the identity matrix as a starting point for Gauss-Jordan |
33 | | - |
34 | | - writable::doubles_matrix<> Ainv(N, N); |
35 | | - |
36 | | - for (int i = 0; i < N; i++) |
37 | | - { |
38 | | - for (int j = 0; j < N; j++) |
39 | | - { |
40 | | - if (i == j) |
41 | | - { |
42 | | - Ainv(i, j) = 1.0; |
43 | | - } |
44 | | - else |
45 | | - { |
46 | | - Ainv(i, j) = 0.0; |
47 | | - } |
48 | | - } |
| 28 | + // Create the identity matrix as a starting point for Gauss-Jordan |
| 29 | + |
| 30 | + writable::doubles_matrix<> Ainv(N, N); |
| 31 | + |
| 32 | + for (int i = 0; i < N; i++) { |
| 33 | + for (int j = 0; j < N; j++) { |
| 34 | + if (i == j) { |
| 35 | + Ainv(i, j) = 1.0; |
| 36 | + } else { |
| 37 | + Ainv(i, j) = 0.0; |
| 38 | + } |
49 | 39 | } |
| 40 | + } |
50 | 41 |
|
51 | | - // Overwrite Ainv by steps with the inverse of A |
52 | | - // in other words, find the echelon form of A |
| 42 | + // Overwrite Ainv by steps with the inverse of A |
| 43 | + // in other words, find the echelon form of A |
53 | 44 |
|
54 | | - for (int i = 0; i < M; i++) |
55 | | - { |
56 | | - double a = Acopy(i, i); |
| 45 | + for (int i = 0; i < M; i++) { |
| 46 | + double a = Acopy(i, i); |
57 | 47 |
|
58 | | - for (int j = 0; j < M; j++) |
59 | | - { |
60 | | - Acopy(i, j) /= a; |
61 | | - Ainv(i, j) /= a; |
62 | | - } |
| 48 | + for (int j = 0; j < M; j++) { |
| 49 | + Acopy(i, j) /= a; |
| 50 | + Ainv(i, j) /= a; |
| 51 | + } |
| 52 | + |
| 53 | + for (int j = 0; j < M; j++) { |
| 54 | + if (i != j) { |
| 55 | + a = Acopy(j, i); |
63 | 56 |
|
64 | | - for (int j = 0; j < M; j++) |
65 | | - { |
66 | | - if (i != j) |
67 | | - { |
68 | | - a = Acopy(j, i); |
69 | | - |
70 | | - for (int k = 0; k < M; k++) |
71 | | - { |
72 | | - Acopy(j, k) -= Acopy(i, k) * a; |
73 | | - Ainv(j, k) -= Ainv(i, k) * a; |
74 | | - } |
75 | | - } |
| 57 | + for (int k = 0; k < M; k++) { |
| 58 | + Acopy(j, k) -= Acopy(i, k) * a; |
| 59 | + Ainv(j, k) -= Ainv(i, k) * a; |
76 | 60 | } |
| 61 | + } |
77 | 62 | } |
| 63 | + } |
78 | 64 |
|
79 | | - return Ainv; |
| 65 | + return Ainv; |
80 | 66 | } |
81 | 67 |
|
82 | | -[[cpp4r::register]] doubles_matrix<> solve_system_(doubles_matrix<> A, doubles_matrix<> b) |
83 | | -{ |
84 | | - // Use Gauss-Jordan to solve the system Ax = b |
| 68 | +[[cpp4r::register]] doubles_matrix<> solve_system_(doubles_matrix<> A, |
| 69 | + doubles_matrix<> b) { |
| 70 | + // Use Gauss-Jordan to solve the system Ax = b |
85 | 71 |
|
86 | | - // Check dimensions |
| 72 | + // Check dimensions |
87 | 73 |
|
88 | | - int N1 = A.nrow(); |
89 | | - int M1 = A.ncol(); |
| 74 | + int N1 = A.nrow(); |
| 75 | + int M1 = A.ncol(); |
90 | 76 |
|
91 | | - if (N1 != M1) |
92 | | - { |
93 | | - stop("A must be a square matrix"); |
94 | | - } |
| 77 | + if (N1 != M1) { |
| 78 | + stop("A must be a square matrix"); |
| 79 | + } |
95 | 80 |
|
96 | | - int N2 = b.nrow(); |
97 | | - int M2 = b.ncol(); |
| 81 | + int N2 = b.nrow(); |
| 82 | + int M2 = b.ncol(); |
98 | 83 |
|
99 | | - if (N1 != N2) |
100 | | - { |
101 | | - stop("b must have the same number of rows as A"); |
102 | | - } |
| 84 | + if (N1 != N2) { |
| 85 | + stop("b must have the same number of rows as A"); |
| 86 | + } |
103 | 87 |
|
104 | | - if (M2 != 1) |
105 | | - { |
106 | | - stop("b must be a column vector"); |
107 | | - } |
| 88 | + if (M2 != 1) { |
| 89 | + stop("b must be a column vector"); |
| 90 | + } |
108 | 91 |
|
109 | | - // Obtain the inverse |
| 92 | + // Obtain the inverse |
110 | 93 |
|
111 | | - // writable::doubles_matrix<> Acopy(N1,N1); |
| 94 | + // writable::doubles_matrix<> Acopy(N1,N1); |
112 | 95 |
|
113 | | - // for (int i = 0; i < N1; i++) |
114 | | - // { |
115 | | - // for (int j = 0; j < N1; j++) |
116 | | - // { |
117 | | - // Acopy(i, j) = A(i, j); |
118 | | - // } |
119 | | - // } |
| 96 | + // for (int i = 0; i < N1; i++) |
| 97 | + // { |
| 98 | + // for (int j = 0; j < N1; j++) |
| 99 | + // { |
| 100 | + // Acopy(i, j) = A(i, j); |
| 101 | + // } |
| 102 | + // } |
120 | 103 |
|
121 | | - // doubles_matrix<> Ainv = invert_matrix_(Acopy); |
| 104 | + // doubles_matrix<> Ainv = invert_matrix_(Acopy); |
122 | 105 |
|
123 | | - // I don´t need to create a copy in this case |
| 106 | + // I don´t need to create a copy in this case |
124 | 107 |
|
125 | | - doubles_matrix<> Ainv = invert_matrix_(A); |
| 108 | + doubles_matrix<> Ainv = invert_matrix_(A); |
126 | 109 |
|
127 | | - // Multiply Ainv by b |
| 110 | + // Multiply Ainv by b |
128 | 111 |
|
129 | | - writable::doubles_matrix<> x(N1, 1); |
| 112 | + writable::doubles_matrix<> x(N1, 1); |
130 | 113 |
|
131 | | - for (int i = 0; i < N1; i++) |
132 | | - { |
133 | | - x(i, 0) = 0.0; |
| 114 | + for (int i = 0; i < N1; i++) { |
| 115 | + x(i, 0) = 0.0; |
134 | 116 |
|
135 | | - for (int j = 0; j < N1; j++) |
136 | | - { |
137 | | - x(i, 0) += Ainv(i, j) * b(j, 0); |
138 | | - } |
| 117 | + for (int j = 0; j < N1; j++) { |
| 118 | + x(i, 0) += Ainv(i, j) * b(j, 0); |
139 | 119 | } |
| 120 | + } |
140 | 121 |
|
141 | | - return x; |
| 122 | + return x; |
142 | 123 | } |
0 commit comments