RcdMathLib_doc
Open Source Library for Linear and Non-linear Algebra
matrix.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 
24 #include <stdint.h>
25 #include <string.h>
26 #include <stdio.h>
27 #include <math.h>
28 #include <stdint.h>
29 #include <stdbool.h>
30 #include <float.h>
31 
32 #include "matrix.h"
33 #include "utils.h"
34 #include "svd.h"
35 
36 static int32_t max(int32_t a, int32_t b)
37 {
38  if (a > b) {
39  return a;
40  }
41  else {
42  return b;
43  }
44 }
45 
46 void matrix_clear(uint8_t m, uint8_t n, matrix_t matrix[m][n])
47 {
48  memset(matrix, 0, sizeof(matrix[0][0]) * m * n);
49 }
50 
51 void matrix_init(uint8_t m, uint8_t n, matrix_t matrix[m][n], matrix_t value)
52 {
53  memset(matrix, value, sizeof(matrix[0][0]) * m * n);
54 }
55 
56 void matrix_transpose(uint8_t m, uint8_t n, matrix_t src_matrix[m][n],
57  matrix_t dest_matrix[n][m])
58 {
59  int i, j;
60 
61  for (i = 0; i < m; i++) {
62  for (j = 0; j < n; j++) {
63  dest_matrix[j][i] = src_matrix[i][j];
64  }
65  }
66 }
67 
68 
69 void matrix_in_place_transpose(uint8_t m, matrix_t matrix[][m])
70 {
71  uint8_t i, j;
72  matrix_t tmp;
73 
74  for (i = 0; i < m; i++) {
75  for (j = i + 1; j < m; j++) {
76  tmp = matrix[i][j];
77  matrix[i][j] = matrix[j][i];
78  matrix[j][i] = tmp;
79  }
80  }
81 }
82 
83 void matrix_copy(uint8_t m, uint8_t n, matrix_t src_matrix[m][n],
84  matrix_t dest_matrix[m][n])
85 {
86  memcpy(dest_matrix, src_matrix, sizeof(matrix_t) * m * n);
87 }
88 
89 void matrix_part_copy(uint8_t m, uint8_t n, matrix_t src_matrix[m][n],
90  uint8_t start_row_ind, uint8_t end_row_ind,
91  uint8_t start_col_ind, uint8_t end_col_ind,
92  uint8_t dest_row_num, uint8_t dest_col_num,
93  matrix_t dest_matrix[][dest_col_num])
94 {
95  int16_t i, j;
96 
97  if ((end_row_ind > start_row_ind) && (end_col_ind > start_col_ind)) {
98  for (i = start_row_ind;
99  (i <= end_row_ind)
100  && ((i - start_row_ind)
101  < dest_row_num);
102  i++) {
103  for (j = start_col_ind;
104  (j <= end_col_ind)
105  && ((j - start_col_ind)
106  < dest_col_num);
107  j++) {
108  dest_matrix[i - start_row_ind][j - start_col_ind] =
109  src_matrix[i][j];
110  }
111  }
112  }
113  else {
114  //copy one element
115  if ((start_row_ind == end_row_ind)
116  && (start_col_ind == end_col_ind)) {
117 
118  dest_matrix[start_row_ind][start_row_ind] =
119  src_matrix[start_row_ind][start_row_ind];
120  }
121  //copy a row vector of the matrix
122  else if (start_row_ind == end_row_ind) {
123  i = 0;
124  for (j = start_col_ind; j <= end_col_ind; j++) {
125  dest_matrix[start_row_ind][i++] =
126  src_matrix[start_row_ind][j];
127  }
128  }
129  //copy a column vector of the matrix
130  else if (start_col_ind == end_col_ind) {
131  j = 0;
132  for (i = start_row_ind; i <= end_row_ind; i++) {
133  dest_matrix[j++][start_col_ind] =
134  src_matrix[i][start_col_ind];
135  }
136  }
137  }
138 }
139 
140 
141 void matrix_print(uint8_t m, uint8_t n, matrix_t matrix[m][n])
142 {
143  int i, j;
144 
145  puts("{");
146  for (i = 0; i < m; i++) {
147  printf("{");
148  for (j = 0; j < n; j++) {
149  printf("%7.4f", matrix[i][j]);
150  if (j < n - 1) {
151  printf(", ");
152  }
153  }
154  printf("}");
155  if (i < m - 1) {
156  printf(", ");
157  }
158  puts("");
159  }
160  printf("}");
161 }
162 
163 
164 void matrix_part_print(uint8_t m, uint8_t n, matrix_t matrix[m][n],
165  uint8_t start_row_ind, uint8_t end_row_ind,
166  uint8_t start_col_ind, uint8_t end_col_ind)
167 {
168  int16_t i, j;
169 
170  if ((end_row_ind > start_row_ind) && (end_col_ind > start_col_ind)) {
171  puts("{");
172  for (i = start_row_ind; i <= end_row_ind; i++) {
173  printf("{");
174  for (j = start_col_ind; j <= end_col_ind; j++) {
175  printf("%7.4f", matrix[i][j]);
176  if (j < end_col_ind) {
177  printf(", ");
178  }
179  }
180  printf("}");
181  if (i < end_row_ind) {
182  printf(", ");
183  }
184  puts("");
185  }
186  printf("}");
187  }
188  else {
189  //print only one element of the matrix
190  if ((start_row_ind == end_row_ind)
191  & (start_col_ind == end_col_ind)) {
192 
193  printf("{%3.4f}", matrix[start_row_ind][start_col_ind]);
194  }
195  // print row vector of the matrix
196  else if (start_row_ind == end_row_ind) {
197  printf("{");
198  for (j = start_col_ind; j <= end_col_ind; j++) {
199  printf("%3.4f", matrix[start_row_ind][j]);
200  if (j < end_col_ind) {
201  printf(", ");
202  }
203  }
204  printf("}");
205  }
206  // print column vector of the matrix
207  else if (start_col_ind == end_col_ind) {
208  puts("{");
209  for (i = start_row_ind; i <= end_row_ind; i++) {
210  printf("%3.4f", matrix[i][start_col_ind]);
211  if (i < end_row_ind) {
212  printf(",\n");
213  }
214  }
215  printf("}");
216  }
217  }
218 }
219 
220 void matrix_flex_print(uint8_t m, uint8_t n, matrix_t matrix[m][n],
221  uint8_t before_dot, uint8_t after_dot)
222 {
223  uint16_t i, j;
224  char format_str_buff[13];
225 
226  sprintf(format_str_buff, "%%%u.%uf", before_dot, after_dot);
227 
228  puts("{");
229  for (i = 0; i < m; i++) {
230  printf("%11s", "{");
231  for (j = 0; j < n; j++) {
232  utils_printf(format_str_buff, matrix[i][j]);
233  if (j < n - 1) {
234  printf(", ");
235  }
236  }
237  printf("}");
238  if (i < m - 1) {
239  printf(", ");
240  }
241  puts("");
242  }
243  printf("%7s\n", "};");
244 }
245 
246 
247 void matrix_flex_part_print(uint8_t m, uint8_t n, matrix_t matrix[m][n],
248  uint8_t start_row_ind, uint8_t end_row_ind,
249  uint8_t start_col_ind, uint8_t end_col_ind,
250  uint8_t before_dot, uint8_t after_dot)
251 {
252  int16_t i, j;
253  char format_str_buff[13];
254 
255  sprintf(format_str_buff, "%%%u.%uf", before_dot, after_dot);
256 
257  if ((end_row_ind > start_row_ind) && (end_col_ind > start_col_ind)) {
258  puts("{");
259  for (i = start_row_ind; i <= end_row_ind; i++) {
260  printf("%11s", "{");
261  for (j = start_col_ind; j <= end_col_ind; j++) {
262  utils_printf(format_str_buff, matrix[i][j]);
263  if (j < end_col_ind) {
264  printf(", ");
265  }
266  }
267  printf("}");
268  if (i < end_row_ind) {
269  printf(", ");
270  }
271  puts("");
272  }
273  printf("%11s\n", "};");
274  }
275  else {
276  //print only one element of the matrix
277  if ((start_row_ind == end_row_ind)
278  & (start_col_ind == end_col_ind)) {
279  printf("{");
280  utils_printf(format_str_buff,
281  matrix[start_row_ind][start_col_ind]);
282  printf("}");
283  }
284  // print row vector of the matrix
285  else if (start_row_ind == end_row_ind) {
286  printf("{");
287  for (j = start_col_ind; j <= end_col_ind; j++) {
288  utils_printf(format_str_buff,
289  matrix[start_row_ind][j]);
290  if (j < end_col_ind) {
291  printf(", ");
292  }
293  }
294  printf("}");
295  }
296  // print column vector of the matrix
297  else if (start_col_ind == end_col_ind) {
298  puts("{");
299  for (i = start_row_ind; i <= end_row_ind; i++) {
300  utils_printf(format_str_buff,
301  matrix[i][start_col_ind]);
302  if (i < end_row_ind) {
303  printf(",\n");
304  }
305  }
306  printf("}");
307  }
308  }
309 }
310 
311 
312 uint8_t matrix_get_rank(uint8_t m, uint8_t n, matrix_t singl_values_arr[],
313  uint8_t length)
314 {
315  uint8_t rank = 0;
316  int i;
317  double eps = pow(2.0, -52.0);
318  double tol = max(m, n) * singl_values_arr[0] * eps;
319 
320  for (i = 0; i < length; i++) {
321  if (singl_values_arr[i] > tol) {
322  rank++;
323  }
324  }
325 
326  return rank;
327 }
328 
329 void matrix_sub(uint8_t m, uint8_t n, matrix_t A[m][n], matrix_t B[m][n],
330  matrix_t A_minus_B[m][n])
331 {
332  uint8_t i, j;
333 
334  for (i = 0; i < m; i++) {
335  for (j = 0; j < n; j++) {
336  A_minus_B[i][j] = A[i][j] - B[i][j];
337  }
338  }
339 }
340 
341 void matrix_add(uint8_t m, uint8_t n, matrix_t A[m][n], matrix_t B[m][n],
342  matrix_t A_plus_B[m][n])
343 {
344  uint8_t i, j;
345 
346  for (i = 0; i < m; i++) {
347  for (j = 0; j < n; j++) {
348  A_plus_B[i][j] = A[i][j] + B[i][j];
349  }
350  }
351 }
352 
353 void matrix_add_to_diag(uint8_t n, matrix_t A[][n], uint8_t diag_el_num,
354  matrix_t value)
355 {
356  uint8_t i;
357 
358  for (i = 0; i < diag_el_num; i++) {
359  A[i][i] += value;
360  }
361 }
362 
363 void matrix_mul(uint8_t a_line_num, uint8_t a_col_num,
364  matrix_t a_matrix[a_line_num][a_col_num],
365  uint8_t b_line_num, uint8_t b_col_num,
366  matrix_t b_matrix[b_line_num][b_col_num],
367  matrix_t dest_matrix[a_line_num][b_col_num])
368 {
369 
370  int i, j, k;
371 
372  if (a_col_num != b_line_num) {
373  puts(
374  "the first matrix column number must be equal with the line number"
375  " of the second matrix !!!");
376  }
377  else {
378 
379  for (i = 0; i < a_line_num; i++) {
380  for (j = 0; j < b_col_num; j++) {
381  dest_matrix[i][j] = 0;
382  for (k = 0; k < b_line_num; k++) {
383  dest_matrix[i][j] += a_matrix[i][k]
384  * b_matrix[k][j];
385  }
386  }
387  }
388  }
389 }
390 
391 void matrix_part_mul(uint8_t a_col_num_max, matrix_t a_matrix[][a_col_num_max],
392  uint8_t b_col_num_max, matrix_t b_matrix[][b_col_num_max],
393  uint8_t a_start_row_ind, uint8_t a_end_row_ind,
394  uint8_t a_start_col_ind, uint8_t a_end_col_ind,
395  uint8_t b_start_row_ind, uint8_t b_end_row_ind,
396  uint8_t b_start_col_ind, uint8_t b_end_col_ind,
397  uint8_t dest_col_size,
398  matrix_t dest_matrix[][dest_col_size])
399 {
400  uint8_t a_line_num;
401  uint8_t b_col_num;
402  uint8_t b_line_num;
403 
404  if ((a_end_col_ind - a_start_col_ind)
405  != (b_end_row_ind - b_start_row_ind)) {
406  puts(
407  "the first matrix column number must be equal with the line number"
408  " of the second matrix !!!");
409  }
410  else {
411 
412  a_line_num = a_end_row_ind - a_start_row_ind + 1;
413  b_line_num = b_end_row_ind - b_start_row_ind + 1;
414  b_col_num = b_end_col_ind - b_start_col_ind + 1;
415 
416  for (uint8_t i = 0; i < a_line_num; i++) {
417  for (uint8_t j = 0; j < b_col_num; j++) {
418  dest_matrix[i][j] = 0;
419  for (uint8_t k = 0; k < b_line_num; k++) {
420  dest_matrix[i][j] +=
421  a_matrix[i
422  + a_start_row_ind][k
423  + a_start_col_ind]
424  * b_matrix[k
425  + b_start_row_ind][j
426  + b_start_col_ind];
427  }
428  }
429  }
430  }
431 }
432 
433 
434 void matrix_mul_vec(uint8_t m, uint8_t n, matrix_t matrix[m][n],
435  matrix_t vec[n], matrix_t dst_arr[m])
436 {
437 
438  uint8_t i, j;
439 
440  if (vec != dst_arr) {
441  memset(dst_arr, 0, m * sizeof(matrix_t));
442  }
443 
444  for (i = 0; i < m; i++) {
445  for (j = 0; j < n; j++) {
446  dst_arr[i] += matrix[i][j] * vec[j];
447  }
448  }
449 }
450 
451 
452 void matrix_vec_mul_matr(uint8_t m, uint8_t n, matrix_t vec[m],
453  matrix_t matrix[m][n], matrix_t dst_arr[n])
454 {
455  uint8_t i, j;
456 
457  //Make sense ?
458  if (vec != dst_arr) {
459  memset(dst_arr, 0, n * sizeof(matrix_t));
460  }
461 
462  for (i = 0; i < n; i++) { //column
463  for (j = 0; j < m; j++) { //rows
464  dst_arr[i] += vec[j] * matrix[j][i];
465  }
466  }
467 }
468 
469 
470 void matrix_mul_scalar_vec_matr(uint8_t m, uint8_t n, matrix_t scalar,
471  matrix_t vec[m], matrix_t matrix[m][n],
472  matrix_t dst_arr[n])
473 {
474  uint8_t i, j;
475 
476  //Make sense ?
477  if (vec != dst_arr) {
478  memset(dst_arr, 0, n * sizeof(matrix_t));
479  }
480 
481  for (i = 0; i < n; i++) { //column
482  for (j = 0; j < m; j++) { //rows
483  dst_arr[i] += scalar * vec[j] * matrix[j][i];
484  }
485  }
486 }
487 
488 
489 void matrix_part_mul_scalar_vec_matr(uint8_t max_m, uint8_t max_n,
490  matrix_t scalar, matrix_t vec[max_m],
491  matrix_t matrix[max_m][max_n],
492  uint8_t begin_row, uint8_t begin_column,
493  matrix_t dst_arr[max_n - begin_row])
494 {
495  int16_t i, j;
496 
497  if (vec != dst_arr) {
498  memset(dst_arr, 0, (max_n - begin_row) * sizeof(matrix_t));
499  }
500 
501  for (i = begin_column; i < max_n; i++) { //column
502  for (j = begin_row; j < max_m; j++) { //rows
503  dst_arr[i - begin_column] += scalar * vec[j - begin_row]
504  * matrix[j][i];
505 
506  }
507  }
508 }
509 
510 
511 void matrix_trans_mul_vec(uint8_t m, uint8_t n, matrix_t A[m][n],
512  uint8_t b_size, matrix_t b_vec[m], matrix_t c_vec[n])
513 {
514  uint8_t i, j;
515 
516  if (m != b_size) {
517  puts(
518  "The vector size should be equal the raw number of the matrix !!!");
519  }
520  else {
521  for (j = 0; j < n; j++) {
522  c_vec[j] = 0;
523  for (i = 0; i < m; i++) {
524  c_vec[j] += A[i][j] * b_vec[i];
525  }
526  }
527  }
528 }
529 
530 
531 void matrix_mul_col_vec_row_vec(uint8_t m, matrix_t col_vec[m], uint8_t n,
532  matrix_t row_vec[n], uint8_t max_n,
533  matrix_t res_mat[][max_n])
534 {
535  uint8_t i, j;
536 
537  for (i = 0; i < m; i++) {
538  for (j = 0; j < n; j++) {
539  res_mat[i][j] = col_vec[i] * row_vec[j];
540  }
541  }
542 }
543 
544 
545 void matrix_trans_mul_itself(uint8_t m, uint8_t n, matrix_t A[m][n],
546  matrix_t AT_mul_A[n][n])
547 {
548  uint8_t i, j, k;
549  matrix_t column_vec[m];
550  matrix_t tmp_col_vec[m];
551 
552  for (i = 0; i < n; i++) {
553  matrix_get_column_vec(m, n, A, i, column_vec);
554  for (j = 0; j < n; j++) {
555  matrix_get_column_vec(m, n, A, j, tmp_col_vec);
556  AT_mul_A[i][j] = 0;
557  for (k = 0; k < m; k++) {
558  AT_mul_A[i][j] += column_vec[k]
559  * tmp_col_vec[k];
560  }
561  }
562  }
563 }
564 
565 void matrix_set_diag_elements(uint8_t m, uint8_t n, matrix_t value,
566  matrix_t diag_matrix[m][n])
567 {
568  uint8_t max = 0;
569 
570  if (m < n) {
571  max = m;
572  }
573  else {
574  max = n;
575  }
576  for (uint8_t i = 0; i < max; i++) {
577  diag_matrix[i][i] = value;
578  }
579 }
580 
581 void matrix_get_diag_mat_new(uint8_t m, uint8_t n, matrix_t diag_matrix[m][n],
582  uint8_t length, matrix_t vec[])
583 {
584  uint8_t i;
585 
586  matrix_clear(m, n, diag_matrix);
587 
588  for (i = 0; i < length; i++) {
589  diag_matrix[i][i] = vec[i];
590  }
591 }
592 
593 
594 void matrix_get_diag_mat(uint8_t m, uint8_t n, matrix_t value,
595  matrix_t diag_matrix[m][n])
596 {
597  matrix_clear(m, n, diag_matrix);
598  matrix_set_diag_elements(m, n, value, diag_matrix);
599 }
600 
601 
602 void matrix_mul_scalar(uint8_t m, uint8_t n, matrix_t mat_src[m][n],
603  matrix_t value, matrix_t mat_dest[m][n])
604 {
605  for (int i = 0; i < m; i++) {
606  for (int j = 0; j < n; j++) {
607  mat_dest[i][j] = value * mat_src[i][j];
608  }
609  }
610 }
611 
612 void matrix_get_column_vec(uint8_t m, uint8_t n, matrix_t matrix[m][n],
613  uint8_t col_num, matrix_t col_vec[m])
614 {
615  uint8_t i;
616 
617  for (i = 0; i < m; i++) {
618  col_vec[i] = matrix[i][col_num];
619  }
620 }
621 
622 
623 void matrix_get_part_column_vec(uint8_t max_m, uint8_t max_n,
624  matrix_t matrix[max_m][max_n], uint8_t col_num,
625  uint8_t offset,
626  matrix_t col_vec[max_m - offset])
627 {
628  uint8_t i;
629 
630  for (i = offset; i < max_m; i++) {
631  col_vec[i - offset] = matrix[i][col_num];
632  }
633 }
634 
636  matrix_t matrix[m][n], uint8_t col_num)
637 {
638  matrix_t elem_max;
639  uint8_t i;
640 
641  elem_max = matrix[0][col_num];
642  for (i = 1; i < m; i++) {
643  if (matrix[i][col_num] > elem_max) {
644  elem_max = matrix[i][col_num];
645  }
646  }
647 
648  return elem_max;
649 }
650 
652  matrix_t matrix[m][n],
653  uint8_t col_num)
654 {
655  matrix_t abs_max_elem;
656  uint8_t i;
657 
658  abs_max_elem = fabs(matrix[0][col_num]);
659  for (i = 1; i < m; i++) {
660  if (fabs(matrix[i][col_num]) > abs_max_elem) {
661  abs_max_elem = fabs(matrix[i][col_num]);
662  }
663  }
664 
665  return abs_max_elem;
666 }
667 
668 matrix_t matrix_get_max_elem_in_part_column(uint8_t max_m, uint8_t max_n,
669  matrix_t matrix[max_m][max_n],
670  uint8_t row_num, uint8_t col_num)
671 {
672  matrix_t max_elem;
673  uint8_t i;
674 
675  max_elem = matrix[row_num][col_num];
676  for (i = row_num + 1; i < max_m; i++) {
677  if (matrix[i][col_num] > max_elem) {
678  max_elem = matrix[i][col_num];
679  }
680  }
681 
682  return max_elem;
683 }
684 
686  matrix_t matrix[max_m][max_n],
687  uint8_t row_num,
688  uint8_t col_num)
689 {
690  matrix_t abs_max_elem;
691  uint8_t i;
692 
693  abs_max_elem = fabs(matrix[row_num][col_num]);
694  for (i = row_num + 1; i < max_m; i++) {
695  if (fabs(matrix[i][col_num]) > abs_max_elem) {
696  abs_max_elem = fabs(matrix[i][col_num]);
697  }
698  }
699 
700  return abs_max_elem;
701 }
702 
704  uint8_t max_n,
705  matrix_t matrix[max_m][max_n],
706  uint8_t row_num, uint8_t col_num,
707  uint8_t *index)
708 {
709  matrix_t abs_max_elem;
710  uint8_t i;
711 
712  abs_max_elem = fabs(matrix[row_num][col_num]);
713  *index = row_num;
714  for (i = row_num + 1; i < max_m; i++) {
715  if (fabs(matrix[i][col_num]) > abs_max_elem) {
716  abs_max_elem = fabs(matrix[i][col_num]);
717  *index = i;
718  }
719 
720  }
721 
722  return abs_max_elem;
723 }
724 
725 void matrix_swap_rows(uint8_t n, matrix_t matrix[][n], uint8_t i, uint8_t j)
726 {
727  matrix_t tmp;
728 
729  for (uint8_t k = 0; k < n; k++) {
730  tmp = matrix[i][k];
731  matrix[i][k] = matrix[j][k];
732  matrix[j][k] = tmp;
733  }
734 }
735 
736 void matrix_part_swap_rows(uint8_t n, matrix_t matrix[][n], uint8_t i,
737  uint8_t j, uint8_t col_begin,
738  uint8_t col_end)
739 {
740  matrix_t tmp;
741 
742  if (col_begin <= col_end) {
743  for (uint8_t k = col_begin; k < (col_end + 1); k++) {
744  tmp = matrix[i][k];
745  matrix[i][k] = matrix[j][k];
746  matrix[j][k] = tmp;
747  }
748  }
749 
750 }
751 
752 double matrix_get_two_norm(uint8_t m, uint8_t n, matrix_t A[][n])
753 {
754  double two_norm = 0.0;
755  matrix_dim_t dim;
756  matrix_t copy_A[m][n];
757 
758  svd_get_U_dim(m, n, &dim);
759  matrix_t U[dim.row_num][dim.col_num];
760  matrix_t S[n][n];
761  matrix_t V[n][n];
762  uint8_t s_length = svd_get_single_values_num(m, n);
763  matrix_t s[s_length];
764 
765  matrix_copy(m, n, A, copy_A);
766  svd(m, n, copy_A, dim.row_num, dim.col_num, U, S, V, s_length, s);
767  printf("U1 = ");
768  matrix_print(dim.row_num, dim.col_num, U);
769  puts("");
770  printf("S1 = ");
771  matrix_print(n, n, S);
772  puts("");
773  printf("V1 = ");
774  matrix_print(n, n, V);
775  puts("");
776 
777  two_norm = S[0][0];
778 
779  return two_norm;
780 }
781 
782 double matrix_get_frob_norm(uint8_t m, uint8_t n, matrix_t A[][n])
783 {
784  double frob_norm = 0.0;
785 
786  for (uint8_t i = 0; i < m; i++) {
787  for (uint8_t j = 0; j < n; j++) {
788  frob_norm += pow(A[i][j], 2);
789  }
790  }
791 
792  return sqrt(frob_norm);
793 }
794 
795 void matrix_get_inv_upp_triang(uint8_t m, uint8_t n, matrix_t U[][n],
796  matrix_t inv_U[][m])
797 {
798  //Zero out the inv_U matrix
799  matrix_clear(n, m, inv_U);
800  for (uint8_t i = 0; i < n; i++) {
801  if (U[i][i] != 0) {
802  inv_U[i][i] = 1 / U[i][i];
803  }
804  else {
805  inv_U[i][i] = FLT_MAX;
806  }
807  for (uint8_t j = 0; j <= (i - 1); j++) {
808  matrix_t u = 0.0;
809  for (uint8_t k = j; k <= (i - 1); k++) {
810  u = u + U[k][i] * inv_U[j][k];
811  }
812  inv_U[j][i] = -u * inv_U[i][i];
813  }
814  }
815 }
816 
817 void matrix_get_inv_low_triang(uint8_t m, uint8_t n, matrix_t L[][n],
818  matrix_t inv_L[][m])
819 {
820  //Zero out the inv_L matrix
821  matrix_clear(n, m, inv_L);
822  for (uint8_t i = 0; i < n; i++) {
823  if ((i <= m) & (i <= n)) {
824  if (L[i][i] != 0) {
825  inv_L[i][i] = 1 / L[i][i];
826  }
827  else {
828  inv_L[i][i] = FLT_MAX;
829  }
830  for (uint8_t j = 0; j <= (i - 1); j++) {
831  matrix_t l = 0.0;
832  for (uint8_t k = j; k <= (i - 1); k++) {
833  l = l + L[i][k] * inv_L[k][j];
834  }
835  inv_L[i][j] = -l * inv_L[i][i];
836  }
837  } //if
838  } //for
839 }
840 
841 void matrix_get_upp_triang(uint8_t m, uint8_t n, matrix_t A[][n],
842  matrix_t tr_up_A[][n])
843 {
844  // Same matrix
845  if ((&A[0][0]) == (&tr_up_A[0][0])) {
846  for (uint8_t i = 1; i < m; i++) {
847  for (uint8_t j = 0; (j < i) && (j < n); j++) {
848  A[i][j] = 0;
849  }
850  }
851 
852  }
853  else {
854  for (uint8_t i = 0; i < m; i++) {
855  for (uint8_t j = 0; j < n; j++) {
856  if (i <= j) {
857  tr_up_A[i][j] = A[i][j];
858  }
859  else {
860  tr_up_A[i][j] = 0;
861  }
862  }
863  }
864  }
865 
866 }
867 
868 void matrix_get_low_triang(uint8_t m, uint8_t n, matrix_t A[][n],
869  matrix_t tr_low_A[][n])
870 {
871 
872  if ((&A[0][0]) == (&tr_low_A[0][0])) {
873  // Same matrix
874  for (uint8_t i = 0; i < m; i++) {
875  for (uint8_t j = i + 1; (j > i) && (j < n); j++) {
876  A[i][j] = 0;
877  }
878  }
879 
880  }
881  else {
882  for (uint8_t i = 0; i < m; i++) {
883  for (uint8_t j = 0; j < n; j++) {
884  if (i >= j) {
885  tr_low_A[i][j] = A[i][j];
886  }
887  else {
888  tr_low_A[i][j] = 0;
889  }
890  }
891  }
892  }
893 }
894 
895 matrix_t matrix_read(uint8_t m, uint8_t n, matrix_t matrix[m][n], uint8_t i,
896  uint8_t j)
897 {
898  return matrix[i][j];
899 }
900 
901 void matrix_write(uint8_t m, uint8_t n, matrix_t matrix[m][n], uint8_t i,
902  uint8_t j, matrix_t val)
903 {
904  matrix[i][j] = val;
905 }
matrix_transpose
void matrix_transpose(uint8_t m, uint8_t n, matrix_t src_matrix[m][n], matrix_t dest_matrix[n][m])
Computes the transpose of a matrix.
Definition: matrix.c:56
matrix_mul_col_vec_row_vec
void matrix_mul_col_vec_row_vec(uint8_t m, matrix_t col_vec[m], uint8_t n, matrix_t row_vec[n], uint8_t max_n, matrix_t res_mat[][max_n])
Compute the multiplication of a column and row vector.
Definition: matrix.c:531
matrix_get_abs_max_elem_in_part_column
matrix_t matrix_get_abs_max_elem_in_part_column(uint8_t max_m, uint8_t max_n, matrix_t matrix[max_m][max_n], uint8_t row_num, uint8_t col_num)
Get the maximum absolute value of a column vector in a sub-matrix.
Definition: matrix.c:685
matrix_mul_scalar
void matrix_mul_scalar(uint8_t m, uint8_t n, matrix_t mat_src[m][n], matrix_t value, matrix_t mat_dest[m][n])
Multiply all elements of a matrix with a specified value.
Definition: matrix.c:602
matrix_get_inv_upp_triang
void matrix_get_inv_upp_triang(uint8_t m, uint8_t n, matrix_t U[][n], matrix_t inv_U[][m])
Computes the inverse an upper triangular matrix.
Definition: matrix.c:795
matrix_part_print
void matrix_part_print(uint8_t m, uint8_t n, matrix_t matrix[m][n], uint8_t start_row_ind, uint8_t end_row_ind, uint8_t start_col_ind, uint8_t end_col_ind)
Display the values of the sub-matrix elements.
Definition: matrix.c:164
matrix_print
void matrix_print(uint8_t m, uint8_t n, matrix_t matrix[m][n])
Display the values of the matrix elements.
Definition: matrix.c:141
matrix_in_place_transpose
void matrix_in_place_transpose(uint8_t m, matrix_t matrix[][m])
Computes the in-place transpose of a matrix.
Definition: matrix.c:69
matrix_dim_t::row_num
uint8_t row_num
the row number
Definition: matrix.h:61
matrix_get_two_norm
double matrix_get_two_norm(uint8_t m, uint8_t n, matrix_t A[][n])
Get the 2-norm of a matrix that is equal to the largest singular value.
Definition: matrix.c:752
matrix_mul_scalar_vec_matr
void matrix_mul_scalar_vec_matr(uint8_t m, uint8_t n, matrix_t scalar, matrix_t vec[m], matrix_t matrix[m][n], matrix_t dst_arr[n])
Compute the multiplication of a scalar with row vector and a matrix.
Definition: matrix.c:470
svd
void svd(uint8_t m, uint8_t n, matrix_t A[m][n], uint8_t u_m, uint8_t u_n, matrix_t U[u_m][u_n], matrix_t S[u_n][n], matrix_t V[n][n], uint8_t sing_vec_length, matrix_t singl_values_vec[])
Compute the Singular-Value Decomposition (SVD) of a matrix.
Definition: svd.c:119
matrix_part_copy
void matrix_part_copy(uint8_t m, uint8_t n, matrix_t src_matrix[m][n], uint8_t start_row_ind, uint8_t end_row_ind, uint8_t start_col_ind, uint8_t end_col_ind, uint8_t dest_row_num, uint8_t dest_col_num, matrix_t dest_matrix[][dest_col_num])
Copy a part of a matrix to another matrix or sub-matrix.
Definition: matrix.c:89
matrix_read
matrix_t matrix_read(uint8_t m, uint8_t n, matrix_t matrix[m][n], uint8_t i, uint8_t j)
Get the value of a matrix at the position (i,j).
Definition: matrix.c:895
matrix_init
void matrix_init(uint8_t m, uint8_t n, matrix_t matrix[m][n], matrix_t value)
Initialize all the elements of the matrix with a specified value.
Definition: matrix.c:51
matrix_set_diag_elements
void matrix_set_diag_elements(uint8_t m, uint8_t n, matrix_t value, matrix_t diag_matrix[m][n])
Set all the diagonal elements of a matrix with a specified value.
Definition: matrix.c:565
matrix_flex_part_print
void matrix_flex_part_print(uint8_t m, uint8_t n, matrix_t matrix[m][n], uint8_t start_row_ind, uint8_t end_row_ind, uint8_t start_col_ind, uint8_t end_col_ind, uint8_t before_dot, uint8_t after_dot)
Display the values of the sub-matrix elements.
Definition: matrix.c:247
matrix_swap_rows
void matrix_swap_rows(uint8_t n, matrix_t matrix[][n], uint8_t i, uint8_t j)
Swaps two rows of a matrix.
Definition: matrix.c:725
matrix_mul
void matrix_mul(uint8_t a_line_num, uint8_t a_col_num, matrix_t a_matrix[a_line_num][a_col_num], uint8_t b_line_num, uint8_t b_col_num, matrix_t b_matrix[b_line_num][b_col_num], matrix_t dest_matrix[a_line_num][b_col_num])
Compute the multiplication of two matrices.
Definition: matrix.c:363
matrix_get_abs_max_elem_in_column
matrix_t matrix_get_abs_max_elem_in_column(uint8_t m, uint8_t n, matrix_t matrix[m][n], uint8_t col_num)
Get the maximum absolute value of a column vector in a matrix.
Definition: matrix.c:651
matrix.h
Matrix computations.
matrix_get_column_vec
void matrix_get_column_vec(uint8_t m, uint8_t n, matrix_t matrix[m][n], uint8_t col_num, matrix_t col_vec[m])
Get a column of a matrix.
Definition: matrix.c:612
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
utils.h
Utilities for linear algebra.
matrix_get_part_column_vec
void matrix_get_part_column_vec(uint8_t max_m, uint8_t max_n, matrix_t matrix[max_m][max_n], uint8_t col_num, uint8_t offset, matrix_t col_vec[max_m - offset])
Get a part of a column of a matrix.
Definition: matrix.c:623
matrix_part_mul_scalar_vec_matr
void matrix_part_mul_scalar_vec_matr(uint8_t max_m, uint8_t max_n, matrix_t scalar, matrix_t vec[max_m], matrix_t matrix[max_m][max_n], uint8_t begin_row, uint8_t begin_column, matrix_t dst_arr[max_n - begin_row])
Compute the multiplication of a scalar with row vector and a sub-matrix.
Definition: matrix.c:489
matrix_get_frob_norm
double matrix_get_frob_norm(uint8_t m, uint8_t n, matrix_t A[][n])
Get the Frobenius norm of a matrix.
Definition: matrix.c:782
matrix_get_abs_max_elem_and_index_in_part_column
matrix_t matrix_get_abs_max_elem_and_index_in_part_column(uint8_t max_m, uint8_t max_n, matrix_t matrix[max_m][max_n], uint8_t row_num, uint8_t col_num, uint8_t *index)
Get the maximum absolute value and its position in a column vector in a sub-matrix.
Definition: matrix.c:703
matrix_dim_t
A structure to define the row and column number of a matrix.
Definition: matrix.h:60
matrix_add
void matrix_add(uint8_t m, uint8_t n, matrix_t A[m][n], matrix_t B[m][n], matrix_t A_plus_B[m][n])
Compute the addition of two matrices.
Definition: matrix.c:341
matrix_trans_mul_itself
void matrix_trans_mul_itself(uint8_t m, uint8_t n, matrix_t A[m][n], matrix_t AT_mul_A[n][n])
Compute the multiplication of the transpose of a matrix with itself.
Definition: matrix.c:545
utils_printf
void utils_printf(char *format_str,...)
Print by using variable format string as well as argument lists.
Definition: utils.c:96
matrix_add_to_diag
void matrix_add_to_diag(uint8_t n, matrix_t A[][n], uint8_t diag_el_num, matrix_t value)
Add a number to diagonal elements of a matrix.
Definition: matrix.c:353
matrix_flex_print
void matrix_flex_print(uint8_t m, uint8_t n, matrix_t matrix[m][n], uint8_t before_dot, uint8_t after_dot)
Display the values of the matrix elements.
Definition: matrix.c:220
matrix_trans_mul_vec
void matrix_trans_mul_vec(uint8_t m, uint8_t n, matrix_t A[m][n], uint8_t b_size, matrix_t b_vec[m], matrix_t c_vec[n])
Compute the multiplication of transposed matrix with column vector.
Definition: matrix.c:511
matrix_vec_mul_matr
void matrix_vec_mul_matr(uint8_t m, uint8_t n, matrix_t vec[m], matrix_t matrix[m][n], matrix_t dst_arr[n])
Compute the multiplication of row vector and a matrix.
Definition: matrix.c:452
matrix_part_swap_rows
void matrix_part_swap_rows(uint8_t n, matrix_t matrix[][n], uint8_t i, uint8_t j, uint8_t col_begin, uint8_t col_end)
Swaps two rows of a sub-matrix.
Definition: matrix.c:736
matrix_write
void matrix_write(uint8_t m, uint8_t n, matrix_t matrix[m][n], uint8_t i, uint8_t j, matrix_t val)
Write a value in a matrix at the position (i,j).
Definition: matrix.c:901
svd.h
Algorithm for the Singular Value Decomposition (SVD).
matrix_get_diag_mat_new
void matrix_get_diag_mat_new(uint8_t m, uint8_t n, matrix_t diag_matrix[m][n], uint8_t length, matrix_t vec[])
Set all the diagonal elements of a matrix with values that are saved in a vector.
Definition: matrix.c:581
svd_get_single_values_num
uint8_t svd_get_single_values_num(uint8_t m, uint8_t n)
Calculate the number of the singular values.
Definition: svd.c:114
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
matrix_get_low_triang
void matrix_get_low_triang(uint8_t m, uint8_t n, matrix_t A[][n], matrix_t tr_low_A[][n])
Gets the lower triangular part of a matrix.
Definition: matrix.c:868
matrix_t
#define matrix_t
Define the data type of the matrix elements.
Definition: matrix.h:38
matrix_get_max_elem_in_column
matrix_t matrix_get_max_elem_in_column(uint8_t m, uint8_t n, matrix_t matrix[m][n], uint8_t col_num)
Get the largest element of a column vector in a matrix.
Definition: matrix.c:635
matrix_get_rank
uint8_t matrix_get_rank(uint8_t m, uint8_t n, matrix_t singl_values_arr[], uint8_t length)
Compute the rank of a matrix.
Definition: matrix.c:312
matrix_get_diag_mat
void matrix_get_diag_mat(uint8_t m, uint8_t n, matrix_t value, matrix_t diag_matrix[m][n])
Create a diagonal matrix with a specified value.
Definition: matrix.c:594
matrix_part_mul
void matrix_part_mul(uint8_t a_col_num_max, matrix_t a_matrix[][a_col_num_max], uint8_t b_col_num_max, matrix_t b_matrix[][b_col_num_max], uint8_t a_start_row_ind, uint8_t a_end_row_ind, uint8_t a_start_col_ind, uint8_t a_end_col_ind, uint8_t b_start_row_ind, uint8_t b_end_row_ind, uint8_t b_start_col_ind, uint8_t b_end_col_ind, uint8_t dest_col_size, matrix_t dest_matrix[][dest_col_size])
Compute the partial multiplication of two matrices.
Definition: matrix.c:391
matrix_clear
void matrix_clear(uint8_t m, uint8_t n, matrix_t matrix[m][n])
Clear all the elements of the vector.
Definition: matrix.c:46
matrix_get_inv_low_triang
void matrix_get_inv_low_triang(uint8_t m, uint8_t n, matrix_t L[][n], matrix_t inv_L[][m])
Computes the inverse a lower triangular matrix.
Definition: matrix.c:817
matrix_sub
void matrix_sub(uint8_t m, uint8_t n, matrix_t A[m][n], matrix_t B[m][n], matrix_t A_minus_B[m][n])
Compute the subtraction of two matrices.
Definition: matrix.c:329
matrix_get_upp_triang
void matrix_get_upp_triang(uint8_t m, uint8_t n, matrix_t A[][n], matrix_t tr_up_A[][n])
Gets the upper triangular part of a matrix.
Definition: matrix.c:841
svd_get_U_dim
void svd_get_U_dim(uint8_t m, uint8_t n, matrix_dim_t *u_dim)
Calculate the dimension of the matrix U.
Definition: svd.c:96
matrix_get_max_elem_in_part_column
matrix_t matrix_get_max_elem_in_part_column(uint8_t max_m, uint8_t max_n, matrix_t matrix[max_m][max_n], uint8_t row_num, uint8_t col_num)
Get the largest element of a column vector in a sub-matrix.
Definition: matrix.c:668
matrix_dim_t::col_num
uint8_t col_num
the column number
Definition: matrix.h:62