RcdMathLib_doc
Open Source Library for Linear and Non-linear Algebra
qr_givens.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 
25 #include <stdio.h>
26 #include <math.h>
27 #include <float.h>
28 
29 #include "matrix.h"
30 #include "qr_givens.h"
31 
32 /*R will be stored in A */
33 int8_t qr_givens_decomp(uint8_t m, uint8_t n, matrix_t A[][n],
34  uint8_t q_col_num,
35  matrix_t Q[][q_col_num], bool reduced)
36 {
37  if (m < n) {
38  puts("A has more columns than rows");
39  return -1;
40  }
41 
42  for (uint8_t j = 0; j < n; j++) {
43  for (uint8_t i = j + 1; i < m; i++) {
44  matrix_t c_s_t_r_vec[4];
45  qr_givens_get_params(A[j][j], A[i][j], c_s_t_r_vec);
46  A[i][j] = c_s_t_r_vec[2]; //A[i][j] = t
47  A[j][j] = c_s_t_r_vec[3]; //A[j][j] = r
48  if (j < n) {
49  for (uint8_t k = j + 1; k < n; k++) {
50  matrix_t c, s;
51  c = c_s_t_r_vec[0];
52  s = c_s_t_r_vec[1];
53  matrix_t tmp = A[j][k];
54  A[j][k] = c * A[j][k] + s * A[i][k];
55  A[i][k] = -s * tmp + c * A[i][k];
56  }
57  }
58 
59  }
60  }
61 
62  uint8_t max_m;
63  if (reduced) {
64  max_m = n;
65  }
66  else {
67  max_m = m;
68  }
69 
70  if (!reduced) {
71  matrix_clear(m, max_m, Q);
72  }
73 
74  matrix_part_copy(m, n, A, 0, m - 1, 0, n - 1, m, max_m, Q);
75 
76  // Build R
77  matrix_get_upp_triang(m, n, A, A);
78 
79  // Build Q
80  for (int16_t j = max_m - 1; j >= 0; j--) {
81  // Zero out column j from row 0 up to row j
82  for (uint8_t k = 0; k <= j; k++) {
83  Q[k][j] = 0;
84  }
85  Q[j][j] = 1;
86 
87  for (int16_t i = m - 1; i >= j + 1; i--) {
88  matrix_t t = Q[i][j];
89  Q[i][j] = 0;
90  matrix_t c = 1 / sqrt(1 + t * t);
91  matrix_t s;
92  if (c != 0) {
93  s = c * t;
94  }
95  else {
96  s = 1;
97  }
98 
99  for (uint8_t k = j; k < max_m; k++) {
100  matrix_t tmp = Q[j][k];
101  Q[j][k] = c * Q[j][k] - s * Q[i][k];
102  Q[i][k] = s * tmp + c * Q[i][k];
103  }
104  } //for
105  }
106 
107  return 1;
108 }
109 
110 void qr_givens_get_params(matrix_t xjj, matrix_t xij, matrix_t c_s_t_r_vec[])
111 {
112  if (xjj != 0) {
113  matrix_t u;
114 
115  c_s_t_r_vec[2] = xij / xjj; //t = xij/xjj
116  u = sqrt(1 + c_s_t_r_vec[2] * c_s_t_r_vec[2]); //u = sqrt(1 + t*t)
117  c_s_t_r_vec[0] = 1 / u; //c = 1/u
118  c_s_t_r_vec[1] = c_s_t_r_vec[0] * c_s_t_r_vec[2]; //s = c*t
119  c_s_t_r_vec[3] = u * xjj; //r = u*xjj
120  }
121  else {
122  c_s_t_r_vec[0] = 0; //c = 0
123  c_s_t_r_vec[1] = 1; //s = 1
124  c_s_t_r_vec[2] = FLT_MAX; //t = infinity
125  c_s_t_r_vec[3] = xij; //r = xij
126  }
127 }
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
qr_givens_get_params
void qr_givens_get_params(matrix_t xjj, matrix_t xij, matrix_t c_s_t_r_vec[])
Compute the Givens parameters.
Definition: qr_givens.c:110
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
matrix.h
Matrix computations.
qr_givens_decomp
int8_t qr_givens_decomp(uint8_t m, uint8_t n, matrix_t A[][n], uint8_t q_col_num, matrix_t Q[][q_col_num], bool reduced)
Computes the QR decomposition of the matrix A by using the Givens algorithm.
Definition: qr_givens.c:33
matrix_t
#define matrix_t
Define the data type of the matrix elements.
Definition: matrix.h:38
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
qr_givens.h
Givens algorithm for the QR-decomposition. Provide necessary methods to construct Q- and R- matrices ...