RcdMathLib_doc
Open Source Library for Linear and Non-linear Algebra
svd.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 
27 #include <stdbool.h>
28 #include <inttypes.h>
29 #include <math.h>
30 #include <string.h>
31 #include <stdio.h>
32 
33 #include "svd.h"
34 #include "matrix.h"
35 #include "vector.h"
36 
37 static void svd_bidiagonal_trans_and_store_diag_elem(uint8_t m, uint8_t n,
38  matrix_t A[m][n],
39  uint8_t u_m, uint8_t u_n,
40  matrix_t U0[u_m][u_n],
41  matrix_t V0[n][n],
42  matrix_t sup_diag_elem_vec[],
43  matrix_t diag_elem_vec[]
44  );
45 
46 static void svd_setup_bi_diagonal_matrix(uint8_t m, uint8_t n, matrix_t A[m][n],
47  uint8_t u_n, matrix_t U1[m][u_n],
48  matrix_t V1[n][n],
49  matrix_t sup_diag_elem_vec[],
50  matrix_t diag_elem_vec[]
51  );
52 
53 static void svd_compute_singular_values_U_V(uint8_t m, uint8_t n, uint8_t u_n,
54  matrix_t U[m][u_n],
55  matrix_t V[n][n],
56  matrix_t sup_diag_elem_vec[],
57  matrix_t diag_elem_vec[]);
58 
59 static void svd_generate_S(uint8_t m, uint8_t n, matrix_t S[n][n],
60  uint8_t sing_vec_length, matrix_t singl_values_vec[]);
61 
62 static int32_t min(int32_t a, int32_t b)
63 {
64  if (a < b) {
65  return a;
66  }
67  else {
68  return b;
69  }
70 }
71 
72 static int32_t max(int32_t a, int32_t b)
73 {
74  if (a > b) {
75  return a;
76  }
77  else {
78  return b;
79  }
80 }
81 
82 static matrix_t get_abs_max(matrix_t arr[], uint8_t length)
83 {
84  uint8_t i;
85  matrix_t max = fabs(arr[0]);
86 
87  for (i = 1; i < length; i++) {
88  if (fabs(arr[i]) > max) {
89  max = fabs(arr[i]);
90  }
91  }
92 
93  return max;
94 }
95 
96 void svd_get_U_dim(uint8_t m, uint8_t n, matrix_dim_t *u_dim)
97 {
98  u_dim->row_num = m;
99  u_dim->col_num = min(m, n);
100 }
101 
102 void svd_get_S_dim(uint8_t m, uint8_t n, matrix_dim_t *s_dim)
103 {
104  s_dim->row_num = min(m, n);
105  s_dim->col_num = n;
106 }
107 
108 void svd_get_V_dim(uint8_t m, uint8_t n, matrix_dim_t *v_dim)
109 {
110  v_dim->row_num = m;
111  v_dim->col_num = n;
112 }
113 
114 uint8_t svd_get_single_values_num(uint8_t m, uint8_t n)
115 {
116  return min(m + 1, n);
117 }
118 
119 void svd(uint8_t m, uint8_t n, matrix_t A[m][n],
120  uint8_t u_m, uint8_t u_n, matrix_t U[u_m][u_n],
121  matrix_t S[u_n][n], matrix_t V[n][n],
122  uint8_t sing_vec_length, matrix_t singl_values_vec[])
123 {
124 
125  int8_t u_n_min = min(m, n);
126  matrix_t sup_diag_elem_vec[n];
127 
128  matrix_clear(m, u_n_min, U);
129  matrix_clear(n, n, S);
130  matrix_clear(n, n, V);
131  memset(singl_values_vec, 0, sizeof(matrix_t) * sing_vec_length);
132 
133  svd_bidiagonal_trans_and_store_diag_elem(m, n, A, u_m, u_n, U, V,
134  sup_diag_elem_vec, singl_values_vec);
135 
136  svd_setup_bi_diagonal_matrix(m, n, A, u_n, U, V, sup_diag_elem_vec,
137  singl_values_vec);
138 
139  svd_compute_singular_values_U_V(m, n, u_n, U, V, sup_diag_elem_vec,
140  singl_values_vec);
141 
142  svd_generate_S(u_n, n, S, sing_vec_length, singl_values_vec);
143 }
144 
166 static void svd_bidiagonal_trans_and_store_diag_elem(uint8_t m, uint8_t n,
167  matrix_t A[m][n],
168  uint8_t u_m, uint8_t u_n,
169  matrix_t U0[u_m][u_n],
170  matrix_t V0[n][n],
171  matrix_t sup_diag_elem_vec[],
172  matrix_t diag_elem_vec[])
173 {
174  int col_trans_num = min(m - 1, n);
175  int row_trans_num = max(0, min(n - 2, m));
176  int upper_lim = max(col_trans_num, row_trans_num);
177 
178  for (int index = 0; index < upper_lim; index++) {
179  if (index < col_trans_num) {
180  diag_elem_vec[index] = 0;
181  for (int row_num = index; row_num < m; row_num++) {
182  diag_elem_vec[index] = hypot(
183  diag_elem_vec[index],
184  A[row_num][index]);
185  }
186  if (diag_elem_vec[index] != 0.0) {
187  if (A[index][index] < 0.0) {
188  diag_elem_vec[index] =
189  -diag_elem_vec[index];
190  }
191  for (int row_num = index; row_num < m;
192  row_num++) {
193  A[row_num][index] /=
194  diag_elem_vec[index];
195  }
196  A[index][index] += 1.0;
197  }
198  diag_elem_vec[index] = -diag_elem_vec[index];
199  }
200  for (int col_num = index + 1; col_num < n; col_num++) {
201  if ((index < col_trans_num)
202  & (diag_elem_vec[index] != 0.0)) {
203 
204  /* Transformation */
205  matrix_t transf_val = 0;
206  for (int i = index; i < m; i++) {
207  transf_val += A[i][index]
208  * A[i][col_num];
209  }
210  transf_val = -transf_val / A[index][index];
211  for (int i = index; i < m; i++) {
212  A[i][col_num] += transf_val
213  * A[i][index];
214  }
215  }
216  sup_diag_elem_vec[col_num] = A[index][col_num];
217  }
218  if (index < col_trans_num) {
219 
220  /* Place a part of the bi-diagonal transformation in U0 */
221  for (int row_num = index; row_num < m; row_num++) {
222  U0[row_num][index] = A[row_num][index];
223  }
224  }
225  if (index < row_trans_num) {
226  sup_diag_elem_vec[index] = 0;
227  for (int i = index + 1; i < n; i++) {
228  sup_diag_elem_vec[index] = hypot(
229  sup_diag_elem_vec[index],
230  sup_diag_elem_vec[i]);
231  }
232  if (sup_diag_elem_vec[index] != 0.0) {
233  if (sup_diag_elem_vec[index + 1] < 0.0) {
234  sup_diag_elem_vec[index] =
235  -sup_diag_elem_vec[index];
236  }
237  for (int i = index + 1; i < n; i++) {
238  sup_diag_elem_vec[i] /=
239  sup_diag_elem_vec[index];
240  }
241  sup_diag_elem_vec[index + 1] += 1.0;
242  }
243  sup_diag_elem_vec[index] = -sup_diag_elem_vec[index];
244  if ((index + 1 < m)
245  & (sup_diag_elem_vec[index] != 0.0)) {
246  matrix_t transf_vec[m];
247 
248  /* Compute the transformation vector */
249  for (int i = index + 1; i < m; i++) {
250  transf_vec[i] = 0.0;
251  }
252  for (int j = index + 1; j < n; j++) {
253  for (int i = index + 1; i < m; i++) {
254  transf_vec[i] +=
255  sup_diag_elem_vec[j]
256  * A[i][j];
257  }
258  }
259 
260  /* Reduce the matrix */
261  for (int col_num = index + 1; col_num < n;
262  col_num++) {
263  matrix_t t =
264  -sup_diag_elem_vec[col_num]
265  / sup_diag_elem_vec[index
266  + 1];
267  for (int row_num = index + 1;
268  row_num < m;
269  row_num++) {
270  A[row_num][col_num] +=
271  t
272  * transf_vec[row_num];
273  }
274  }
275  }
276 
277  /* Place a part of the bi-diagonal transformation in V0 */
278  for (int row_num = index + 1; row_num < n; row_num++) {
279  V0[row_num][index] = sup_diag_elem_vec[row_num];
280  }
281  }
282  }
283 }
284 
304 static void svd_setup_bi_diagonal_matrix(uint8_t m, uint8_t n, matrix_t A[m][n],
305  uint8_t u_n, matrix_t U1[m][u_n],
306  matrix_t V1[n][n],
307  matrix_t sup_diag_elem_vec[],
308  matrix_t diag_elem_vec[])
309 {
310 
311  int col_trans_num = min(m - 1, n);
312  int row_trans_num = max(0, min(n - 2, m));
313  int8_t upp_col_num = min(m, n);
314  int p_index = min(n, m + 1);
315 
316  if (col_trans_num < n) {
317  diag_elem_vec[col_trans_num] = A[col_trans_num][col_trans_num];
318  }
319  if (m < p_index) {
320  diag_elem_vec[p_index - 1] = 0.0;
321  }
322  if (row_trans_num + 1 < p_index) {
323  sup_diag_elem_vec[row_trans_num] =
324  A[row_trans_num][p_index - 1];
325  }
326  sup_diag_elem_vec[p_index - 1] = 0.0;
327 
328  /* Generate the matrix U1 */
329  for (int j = col_trans_num; j < upp_col_num; j++) {
330  for (int i = 0; i < m; i++) {
331  U1[i][j] = 0.0;
332  }
333  U1[j][j] = 1.0;
334  }
335  for (int k = col_trans_num - 1; k >= 0; k--) {
336  if (diag_elem_vec[k] != 0.0) {
337  for (int j = k + 1; j < upp_col_num; j++) {
338  matrix_t factor_t = 0;
339  for (int i = k; i < m; i++) {
340  factor_t += U1[i][k] * U1[i][j];
341  }
342  factor_t = -factor_t / U1[k][k];
343  for (int i = k; i < m; i++) {
344  U1[i][j] += factor_t * U1[i][k];
345  }
346  }
347  for (int i = k; i < m; i++) {
348  U1[i][k] = -U1[i][k];
349  }
350  U1[k][k] += 1.0;
351  for (int i = 0; i < k - 1; i++) {
352  U1[i][k] = 0.0;
353  }
354  }
355  else {
356  for (int i = 0; i < m; i++) {
357  U1[i][k] = 0.0;
358  }
359  U1[k][k] = 1.0;
360  }
361  }
362 
363  /* Generate the matrix V1 */
364  for (int k = n - 1; k >= 0; k--) {
365  if ((k < row_trans_num) & (sup_diag_elem_vec[k] != 0.0)) {
366  for (int j = k + 1; j < n; j++) {
367  matrix_t t = 0;
368  for (int i = k + 1; i < n; i++) {
369  t += V1[i][k] * V1[i][j];
370  }
371  t = -t / V1[k + 1][k];
372  for (int i = k + 1; i < n; i++) {
373  V1[i][j] += t * V1[i][k];
374  }
375  }
376  }
377  for (int i = 0; i < n; i++) {
378  V1[i][k] = 0.0;
379  }
380  V1[k][k] = 1.0;
381  }
382 }
383 
406 static void svd_compute_singular_values_U_V(uint8_t m, uint8_t n, uint8_t u_n,
407  matrix_t U[m][u_n],
408  matrix_t V[n][n],
409  matrix_t sup_diag_elem_vec[],
410  matrix_t diag_elem_vec[])
411 {
412  int p_index = min(n, m + 1);
413  int max_k = p_index - 1;
414  int iter_num = 0;
415  matrix_t eps = pow(2.0, -52.0);
416  matrix_t tol = pow(2.0, -966.0);
417 
418  while (p_index > 0) {
419  int k_index, case_num;
420 
421  for (k_index = p_index - 2; k_index >= -1; k_index--) {
422  if (k_index == -1) {
423  break;
424  }
425  if (fabs(sup_diag_elem_vec[k_index])
426  <=
427  tol
428  + eps
429  * (fabs(
430  diag_elem_vec[k_index])
431  + fabs(
432  diag_elem_vec[k_index
433  + 1]))) {
434  sup_diag_elem_vec[k_index] = 0.0;
435  break;
436  }
437  }
438  if (k_index == p_index - 2) {
439  case_num = SVD_ORDER_ABSOLUTE_SING_VALUES; //4
440  }
441  else {
442  int ks_index;
443  for (ks_index = p_index - 1; ks_index >= k_index;
444  ks_index--) {
445  if (ks_index == k_index) {
446  break;
447  }
448  matrix_t t =
449  (ks_index != p_index ?
450  fabs(
451  sup_diag_elem_vec[ks_index]) :
452  0.)
453  +
454  (ks_index
455  != k_index
456  + 1 ?
457  fabs(
458  sup_diag_elem_vec[ks_index
459  - 1]) :
460  0.);
461  if (fabs(diag_elem_vec[ks_index])
462  <= tol + eps * t) {
463  diag_elem_vec[ks_index] = 0.0;
464  break;
465  }
466  }
467  if (ks_index == k_index) {
468  case_num = SVD_QR_STEP; //3
469  }
470  else if (ks_index == p_index - 1) {
471  case_num = SVD_COMPUTE_NEGLIGIBLE_VALUES; //1
472  }
473  else {
474  case_num = SVD_SPLIT_AT_NEGLIGIBLE_VALUES; //2
475  k_index = ks_index;
476  }
477  }
478  k_index++;
479 
480  switch (case_num) {
481 
482  /* Compute negligible values */
484  matrix_t f_val = sup_diag_elem_vec[p_index - 2];
485  sup_diag_elem_vec[p_index - 2] = 0.0;
486  for (int j = p_index - 2; j >= k_index; j--) {
487  matrix_t t_val = hypot(diag_elem_vec[j], f_val);
488  matrix_t split_cs_fac = diag_elem_vec[j]
489  / t_val;
490  matrix_t split_sn_fac = f_val / t_val;
491  diag_elem_vec[j] = t_val;
492  if (j != k_index) {
493  f_val = -split_sn_fac
494  * sup_diag_elem_vec[j
495  - 1];
496  sup_diag_elem_vec[j - 1] = split_cs_fac
497  * sup_diag_elem_vec[j
498  - 1];
499  }
500  /* Compute the V matrix */
501  for (int i = 0; i < n; i++) {
502  t_val =
503  split_cs_fac * V[i][j]
504  + split_sn_fac
505  * V[i][p_index
506  - 1];
507  V[i][p_index - 1] =
508  -split_sn_fac * V[i][j]
509  + split_cs_fac
510  * V[i][p_index
511  - 1];
512  V[i][j] = t_val;
513  }
514  }
515  }
516  break;
517 
518  /* Split at negligible values sup_diag_elem_vec(k) */
520  matrix_t f_val = sup_diag_elem_vec[k_index - 1];
521  sup_diag_elem_vec[k_index - 1] = 0.0;
522  for (int j = k_index; j < p_index; j++) {
523  matrix_t t_val = hypot(diag_elem_vec[j], f_val);
524  matrix_t split_cs_fac = diag_elem_vec[j]
525  / t_val;
526  matrix_t split_sn_fac = f_val / t_val;
527  diag_elem_vec[j] = t_val;
528  f_val = -split_sn_fac * sup_diag_elem_vec[j];
529  sup_diag_elem_vec[j] = split_cs_fac
530  * sup_diag_elem_vec[j];
531 
532  /* Compute the U matrix */
533  for (int i = 0; i < m; i++) {
534  t_val =
535  split_cs_fac * U[i][j]
536  + split_sn_fac
537  * U[i][k_index
538  - 1];
539  U[i][k_index - 1] =
540  -split_sn_fac * U[i][j]
541  + split_cs_fac
542  * U[i][k_index
543  - 1];
544  U[i][j] = t_val;
545  }
546  }
547  }
548  break;
549 
550  /* The QR Step */
551  case SVD_QR_STEP: {
552 
553  matrix_t buff[5] = { diag_elem_vec[p_index - 1],
554  diag_elem_vec[p_index - 2],
555  sup_diag_elem_vec[p_index - 2],
556  diag_elem_vec[k_index],
557  sup_diag_elem_vec[k_index] };
558  matrix_t scale = get_abs_max(buff, 5);
559 
560  matrix_t p_diag_elem = diag_elem_vec[p_index - 1]
561  / scale;
562  matrix_t p_minus_diag_elem = diag_elem_vec[p_index - 2]
563  / scale;
564  matrix_t p_minus_sup_diag_elem =
565  sup_diag_elem_vec[p_index - 2] / scale;
566  matrix_t k_diag_elem = diag_elem_vec[k_index] / scale;
567  matrix_t k_sup_diag_elem = sup_diag_elem_vec[k_index]
568  / scale;
569  matrix_t d_val = ((p_minus_diag_elem + p_diag_elem)
570  * (p_minus_diag_elem - p_diag_elem)
571  + p_minus_sup_diag_elem
572  * p_minus_sup_diag_elem)
573  / 2.0;
574  matrix_t e_val = (p_diag_elem * p_minus_sup_diag_elem)
575  * (p_diag_elem * p_minus_sup_diag_elem);
576  matrix_t wilkinson_shift = 0.0;
577  if ((d_val != 0.0) | (e_val != 0.0)) {
578  wilkinson_shift = sqrt(d_val * d_val + e_val);
579  if (d_val < 0.0) {
580  wilkinson_shift = -wilkinson_shift;
581  }
582  wilkinson_shift = e_val
583  / (d_val + wilkinson_shift);
584  }
585  matrix_t f_val = (k_diag_elem + p_diag_elem)
586  * (k_diag_elem - p_diag_elem)
587  + wilkinson_shift;
588  matrix_t g_val = k_diag_elem * k_sup_diag_elem;
589 
590  for (int j = k_index; j < p_index - 1; j++) {
591  matrix_t t = hypot(f_val, g_val);
592  matrix_t split_cs_fac = f_val / t;
593  matrix_t split_sn_fac = g_val / t;
594  if (j != k_index) {
595  sup_diag_elem_vec[j - 1] = t;
596  }
597  f_val =
598  split_cs_fac * diag_elem_vec[j]
599  + split_sn_fac
600  * sup_diag_elem_vec[j];
601  sup_diag_elem_vec[j] =
602  split_cs_fac
603  * sup_diag_elem_vec[j]
604  - split_sn_fac
605  * diag_elem_vec[j];
606  g_val = split_sn_fac * diag_elem_vec[j + 1];
607  diag_elem_vec[j + 1] = split_cs_fac
608  * diag_elem_vec[j + 1];
609 
610  /* Compute the V matrix */
611  for (int i = 0; i < n; i++) {
612  t =
613  split_cs_fac * V[i][j]
614  + split_sn_fac
615  * V[i][j
616  + 1];
617  V[i][j + 1] =
618  -split_sn_fac * V[i][j]
619  + split_cs_fac
620  * V[i][j
621  + 1];
622  V[i][j] = t;
623  }
624 
625  t = hypot(f_val, g_val);
626  split_cs_fac = f_val / t;
627  split_sn_fac = g_val / t;
628  diag_elem_vec[j] = t;
629  f_val =
630  split_cs_fac
631  * sup_diag_elem_vec[j]
632  + split_sn_fac
633  * diag_elem_vec[j
634  + 1];
635  diag_elem_vec[j + 1] =
636  -split_sn_fac
637  * sup_diag_elem_vec[j]
638  + split_cs_fac
639  * diag_elem_vec[j
640  + 1];
641  g_val = split_sn_fac * sup_diag_elem_vec[j + 1];
642  sup_diag_elem_vec[j + 1] = split_cs_fac
643  * sup_diag_elem_vec[j + 1];
644  if (j < m - 1) {
645  for (int i = 0; i < m; i++) {
646  t =
647  split_cs_fac
648  * U[i][j]
649  + split_sn_fac
650  * U[i][j
651  + 1];
652  U[i][j + 1] =
653  -split_sn_fac
654  * U[i][j]
655  + split_cs_fac
656  * U[i][j
657  + 1];
658  U[i][j] = t;
659  }
660  }
661  }
662  sup_diag_elem_vec[p_index - 2] = f_val;
663  iter_num++;
664  }
665  break;
666 
667  /* Order the absolute singular values */
669 
670  if (diag_elem_vec[k_index] <= 0.0) {
671  diag_elem_vec[k_index] =
672  (diag_elem_vec[k_index] < 0.0 ?
673  -diag_elem_vec[k_index] :
674  0.0);
675 
676  for (int i = 0; i < n; i++) {
677  V[i][k_index] = -V[i][k_index];
678  }
679  }
680 
681  while (k_index < max_k) {
682  if (diag_elem_vec[k_index]
683  >= diag_elem_vec[k_index + 1]) {
684  break;
685  }
686  matrix_t k_diag_elem = diag_elem_vec[k_index];
687  diag_elem_vec[k_index] = diag_elem_vec[k_index
688  + 1];
689  diag_elem_vec[k_index + 1] = k_diag_elem;
690  if (k_index < n - 1) {
691  for (int i = 0; i < n; i++) {
692  k_diag_elem = V[i][k_index + 1];
693  V[i][k_index + 1] =
694  V[i][k_index];
695  V[i][k_index] = k_diag_elem;
696  }
697  }
698  if (k_index < m - 1) {
699  for (int i = 0; i < m; i++) {
700  k_diag_elem = U[i][k_index + 1];
701  U[i][k_index + 1] =
702  U[i][k_index];
703  U[i][k_index] = k_diag_elem;
704  }
705  }
706  k_index++;
707  }
708  iter_num = 0;
709  p_index--;
710  }
711  break;
712  }
713  }
714 }
715 
732 static void svd_generate_S(uint8_t m, uint8_t n, matrix_t S[n][n],
733  uint8_t sing_vec_length, matrix_t singl_values_vec[])
734 {
735  uint8_t i;
736 
737  for (i = 0; i < m; i++) {
738  for (uint8_t j = 0; j < n; j++) {
739  S[i][j] = 0.0;
740  }
741  if (i < sing_vec_length) {
742  S[i][i] = singl_values_vec[i];
743  }
744  }
745 }
746 
747 void svd_get_reciproc_singular_values(uint8_t m, uint8_t n, uint8_t length,
748  matrix_t singl_values_arr[],
749  matrix_t recip_singl_values_arr[]
750  )
751 {
752  uint8_t i;
753  matrix_t eps = pow(2.0, -52.0);
754  matrix_t tol = max(m, n) * singl_values_arr[0] * eps;
755 
756  for (i = 0; i < length; i++) {
757  if (fabs(singl_values_arr[i]) >= tol) {
758  recip_singl_values_arr[i] = 1.0 / singl_values_arr[i];
759  }
760  else {
761  recip_singl_values_arr[i] = 0;
762  }
763  }
764 }
765 
766 void svd_compute_print_U_S_V_s(uint8_t m, uint8_t n, matrix_t matrix_arr[m][n],
767  uint8_t i)
768 {
769  printf(
770  "########################## Test %d ##########################\n",
771  i);
772  matrix_dim_t u_dim;
773  svd_get_U_dim(m, n, &u_dim);
774  matrix_t U[u_dim.row_num][u_dim.col_num];
775  matrix_t S[u_dim.col_num][n];
776  matrix_t V[n][n];
777  uint8_t s_length = svd_get_single_values_num(m, n);
778 
779  printf("matrix%d =\n", i);
780  matrix_print(m, n, matrix_arr);
781  puts("");
782  matrix_t s[s_length];
783  svd(m, n, matrix_arr, u_dim.row_num, u_dim.col_num, U, S, V, s_length,
784  s);
785 
786  printf("U%d =\n", i);
787  matrix_print(u_dim.row_num, u_dim.col_num, U);
788  puts("");
789 
790  printf("S%d =\n", i);
791  matrix_print(u_dim.col_num, n, S);
792  puts("");
793 
794  printf("V%d =\n", i);
795  matrix_print(n, n, V);
796  puts("");
797 
798  printf("s%d = ", i);
799  printf("{");
800  for (uint8_t i = 0; i < s_length; i++) {
801  printf("%5.4f ", s[i]);
802  if (i < s_length - 1) {
803  printf(", ");
804  }
805  }
806  printf("}");
807 
808  puts("");
809 }
svd_compute_print_U_S_V_s
void svd_compute_print_U_S_V_s(uint8_t m, uint8_t n, matrix_t matrix_arr[m][n], uint8_t i)
Compute and print the SVD of a matrix.
Definition: svd.c:766
svd_get_reciproc_singular_values
void svd_get_reciproc_singular_values(uint8_t m, uint8_t n, uint8_t length, matrix_t singl_values_arr[], matrix_t recip_singl_values_arr[])
Compute the reciprocal singular values.
Definition: svd.c:747
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_dim_t::row_num
uint8_t row_num
the row number
Definition: matrix.h:61
SVD_SPLIT_AT_NEGLIGIBLE_VALUES
#define SVD_SPLIT_AT_NEGLIGIBLE_VALUES
The case of splitting at negligible values.
Definition: svd.h:43
SVD_QR_STEP
#define SVD_QR_STEP
The case of the QR-step.
Definition: svd.h:48
SVD_COMPUTE_NEGLIGIBLE_VALUES
#define SVD_COMPUTE_NEGLIGIBLE_VALUES
The case of computing negligible values.
Definition: svd.h:38
matrix.h
Matrix computations.
matrix_dim_t
A structure to define the row and column number of a matrix.
Definition: matrix.h:60
svd_get_V_dim
void svd_get_V_dim(uint8_t m, uint8_t n, matrix_dim_t *v_dim)
Calculate the dimension of the matrix V.
Definition: svd.c:108
svd.h
Algorithm for the Singular Value Decomposition (SVD).
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_t
#define matrix_t
Define the data type of the matrix elements.
Definition: matrix.h:38
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
vector.h
Vector computations.
svd_get_S_dim
void svd_get_S_dim(uint8_t m, uint8_t n, matrix_dim_t *s_dim)
Calculate the dimension of the matrix S.
Definition: svd.c:102
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_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
SVD_ORDER_ABSOLUTE_SING_VALUES
#define SVD_ORDER_ABSOLUTE_SING_VALUES
The case of the order of absolute singular values.
Definition: svd.h:53
matrix_dim_t::col_num
uint8_t col_num
the column number
Definition: matrix.h:62