RcdMathLib_doc
Open Source Library for Linear and Non-linear Algebra
solve_test.c
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
3  * 2020 Freie Universität Berlin
4  *
5  * This file is subject to the terms and conditions of the GNU Lesser General
6  * Public License v2.1. See the file LICENSE in the top level directory for more
7  * details.
8  */
9 
23 #include <stdio.h>
24 
25 #include "solve.h"
26 #include "matrix.h"
27 #include "vector.h"
28 
29 void solve_test(void)
30 {
31  int8_t error_no;
32 
33  puts("############ Solve Linear Algebra Module ############");
34 
35  puts("******* Determined system *******");
36  uint8_t m = 3;
37  uint8_t n = m;
38  enum ALGORITHM algo = Moore_Penrose;
39  matrix_t A1[3][3] = { { 2, 1, 1 },
40  { -1, 1, -1 },
41  { 1, 2, 3 } };
42  matrix_t copy_A1[m][n];
43  matrix_t b1[3] = { 2, 3, -10 };
44  matrix_t x[3];
45 
46  printf("A1 = ");
47  matrix_flex_print(m, n, A1, 3, 1);
48 
49  puts("Moore Penrose: ");
50  matrix_copy(m, n, A1, copy_A1);
51  error_no = solve(m, n, copy_A1, b1, x, algo);
52  puts("");
53 
54  if (error_no > 0) {
55  printf("x1 = ");
56  vector_flex_print(n, x, 3, 1);
57  puts("");
58  }
59 
60  puts("Householder: ");
61  algo = Householder;
62  matrix_copy(m, n, A1, copy_A1);
63  error_no = solve(m, n, copy_A1, b1, x, algo);
64 
65  if (error_no > 0) {
66  puts("");
67  printf("x1 = ");
68  vector_flex_print(n, x, 3, 1);
69  puts("");
70  }
71 
72  puts("Givens: ");
73  algo = Givens;
74  matrix_copy(m, n, A1, copy_A1);
75  error_no = solve(m, n, copy_A1, b1, x, algo);
76 
77  if (error_no > 0) {
78  puts("");
79  printf("x1 = ");
80  vector_flex_print(n, x, 3, 1);
81  puts("");
82  }
83 
84  puts("Gauss: ");
85  algo = Gauss;
86  matrix_copy(m, n, A1, copy_A1);
87  error_no = solve(m, n, copy_A1, b1, x, algo);
88 
89  if (error_no > 0) {
90  puts("");
91  printf("x1 = ");
92  vector_flex_print(n, x, 3, 1);
93  puts("");
94  }
95 
96  puts("******* Overdetermined system *******");
97  m = 4;
98  n = 3;
99  matrix_t A2[4][3] = { { 25, 5, 1 },
100  { 64, 8, 1 },
101  { 100, 10, 1 },
102  { 144, 12, 1 }, };
103  matrix_t copy_A2[m][n];
104  matrix_t b2[4] = { 106.8, 177.2, 232.6, 279.2 };
105  printf("A2 = ");
106  matrix_flex_print(m, n, A2, 3, 1);
107 
108  puts("Moore Penrose: ");
109  algo = Moore_Penrose;
110  matrix_copy(m, n, A2, copy_A2);
111  error_no = solve(m, n, copy_A2, b2, x, algo);
112 
113  if (error_no > 0) {
114  printf("x2 = ");
115  vector_flex_print(n, x, 3, 4);
116  puts("");
117  }
118 
119  puts("Householder: ");
120  algo = Householder;
121  matrix_copy(m, n, A2, copy_A2);
122  error_no = solve(m, n, copy_A2, b2, x, algo);
123 
124  if (error_no > 0) {
125  printf("x2 = ");
126  vector_flex_print(n, x, 3, 4);
127  puts("");
128  }
129 
130  puts("Givens: ");
131  algo = Givens;
132  matrix_copy(m, n, A2, copy_A2);
133  error_no = solve(m, n, copy_A2, b2, x, algo);
134 
135  if (error_no > 0) {
136  printf("x2 = ");
137  vector_flex_print(n, x, 3, 4);
138  puts("");
139  }
140 
141  puts("******* Underdetermined system *******");
142  m = 3;
143  n = 4;
144  matrix_t A3[3][4] = { { 2, -4, 2, -14 },
145  { -1, 2, -2, 11 },
146  { -1, 2, -1, 7 } };
147  matrix_t copy_A3[m][n];
148  matrix_t b3[3] = { 10, -6, -5 };
149  printf("A3 = ");
150  matrix_flex_print(m, n, A3, 3, 1);
151 
152  puts("Moore Penrose: ");
153  algo = Moore_Penrose;
154  matrix_copy(m, n, A3, copy_A3);
155  error_no = solve(m, n, A3, b3, x, algo);
156 
157  if (error_no > 0) {
158  printf("x3 = ");
159  vector_flex_print(n, x, 3, 4);
160  puts("");
161  }
162 
163  puts("Householder: ");
164  algo = Householder;
165  matrix_copy(m, n, A3, copy_A3);
166  error_no = solve(m, n, copy_A3, b3, x, algo);
167 
168  if (error_no > 0) {
169  printf("x3 = ");
170  vector_flex_print(n, x, 3, 4);
171  puts("");
172  }
173 
174  puts("Givens: ");
175  algo = Givens;
176  matrix_copy(m, n, A3, copy_A3);
177  error_no = solve(m, n, copy_A3, b3, x, algo);
178 
179  if (error_no > 0) {
180  printf("x3 = ");
181  vector_flex_print(n, x, 3, 4);
182  puts("");
183  }
184 
185  puts("*** Test: Solve with LU decomposition ***");
186  algo = Gauss;
187  matrix_t A4[5][5] = {
188  { 0.1712, 0.0971, 0.0344, 0.1869, 0.7547 },
189  { 0.7060, 0.8235, 0.4387, 0.4898, 0.2760 },
190  { 0.0318, 0.6948, 0.3816, 0.4456, 0.6797 },
191  { 0.2769, 0.3171, 0.7655, 0.6463, 0.6551 },
192  { 0.0462, 0.9502, 0.7952, 0.7094, 0.1626 }
193  };
194  matrix_t x4[5];
195  matrix_t b4[5] = { 0.1190, 0.4984, 0.9597, 0.3404, 0.5853 };
196 
197  m = 5;
198  n = 5;
199  matrix_t copy_A4[m][m];
200  matrix_copy(m, n, A4, copy_A4);
201  error_no = solve(m, n, copy_A4, b4, x4, algo);
202 
203  if (error_no > 0) {
204  printf("x4 = ");
205  vector_flex_print(n, x4, 3, 4);
206  puts("");
207  }
208 
209  matrix_t dot_b4[5];
210  matrix_mul_vec(m, m, A4, x4, dot_b4);
211  printf("dot_b4 = ");
212  vector_flex_print(n, dot_b4, 3, 4);
213  puts("");
214 }
215 
217 {
218  matrix_t A[10][5] = { { 0.8147, 0.1576, 0.6557, 0.7060, 0.4387 },
219  { 0.9058, 0.9706, 0.0357, 0.0318, 0.3816 },
220  { 0.1270, 0.9572, 0.8491, 0.2769, 0.7655 },
221  { 0.9134, 0.4854, 0.9340, 0.0462, 0.7952 },
222  { 0.6324, 0.8003, 0.6787, 0.0971, 0.1869 },
223  { 0.0975, 0.1419, 0.7577, 0.8235, 0.4898 },
224  { 0.2785, 0.4218, 0.7431, 0.6948, 0.4456 },
225  { 0.5469, 0.9157, 0.3922, 0.3171, 0.6463 },
226  { 0.9575, 0.7922, 0.6555, 0.9502, 0.7094 },
227  { 0.9649, 0.9595, 0.1712, 0.0344, 0.7547 } };
228  uint8_t m = 10, n = 5;
229  int8_t error_no;
230  enum ALGORITHM algo;
231  matrix_t copy_A[m][n];
232  matrix_t b[10] = { 1.2902, 0.8819, 0.9721, 1.2347, 0.9185, 0.9844, 1.0627,
233  1.0280, 1.7283, 1.0618 };
234  matrix_t x[n];
235 
236  printf("A = ");
237  matrix_flex_print(m, n, A, 3, 4);
238 
239  puts("Moore Penrose: ");
240  algo = Moore_Penrose;
241  matrix_copy(m, n, A, copy_A);
242  error_no = solve(m, n, copy_A, b, x, algo);
243  puts("");
244 
245  if (error_no > 0) {
246  printf("x = ");
247  vector_flex_print(n, x, 3, 4);
248  puts("");
249  }
250 
251  puts("Householder: ");
252  algo = Householder;
253  matrix_copy(m, n, A, copy_A);
254  error_no = solve(m, n, copy_A, b, x, algo);
255 
256  if (error_no > 0) {
257  puts("");
258  printf("x = ");
259  vector_flex_print(n, x, 3, 4);
260  puts("");
261  }
262 
263  puts("Givens: ");
264  algo = Givens;
265  matrix_copy(m, n, A, copy_A);
266  error_no = solve(m, n, copy_A, b, x, algo);
267 
268  if (error_no > 0) {
269  puts("");
270  printf("x = ");
271  vector_flex_print(n, x, 3, 4);
272  puts("");
273  }
274 }
matrix_mul_vec
void matrix_mul_vec(uint8_t m, uint8_t n, matrix_t matrix[m][n], matrix_t vec[n], matrix_t dst_arr[m])
Compute the multiplication of a matrix with a column vector.
Definition: matrix.c:434
Givens
Givens algorithm.
Definition: pseudo_inverse.h:33
solve_big_matrix_test
void solve_big_matrix_test(void)
Example of solving an (10,5) linear equation system.
Definition: solve_test.c:216
vector_flex_print
void vector_flex_print(uint32_t length, vector_t arr[], uint8_t before_dot, uint8_t after_dot)
Display the values of the vector's elements.
Definition: vector.c:284
Moore_Penrose
Moore–Penrose algorithm.
Definition: pseudo_inverse.h:31
solve.h
Enables to solve systems of linear equations Ax = b for x.
matrix.h
Matrix computations.
matrix_flex_print
void matrix_flex_print(uint8_t m, uint8_t n, matrix_t matrix[m][n], uint8_t before_dec, uint8_t after_dec)
Display the values of the matrix elements.
Definition: matrix.c:220
Householder
Householder algorithm.
Definition: pseudo_inverse.h:32
ALGORITHM
ALGORITHM
Possible algorithms to compute the pseudo-inverse matrix.
Definition: pseudo_inverse.h:30
solve
int8_t solve(uint8_t m, uint8_t n, matrix_t A[][n], matrix_t b[m], matrix_t x_sol[n], enum ALGORITHM algo)
Solve an (m n) linear system Ax = b by using the Moore–Penrose, Householder, or the Givens algorithm...
Definition: solve.c:35
Gauss
Gaussian elimination with pivoting algorithm.
Definition: pseudo_inverse.h:34
matrix_copy
void matrix_copy(uint8_t m, uint8_t n, matrix_t src_matrix[m][n], matrix_t dest_matrix[m][n])
Copy the elements of a matrix to another matrix.
Definition: matrix.c:83
solve_test
void solve_test(void)
Examples of solving linear equation systems.
Definition: solve_test.c:29
matrix_t
#define matrix_t
Define the data type of the matrix elements.
Definition: matrix.h:38
vector.h
Vector computations.