RcdMathLib_doc
Open Source Library for Linear and Non-linear Algebra
qr_householder.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 <stdbool.h>
27 #include <math.h>
28 
29 #include "matrix.h"
30 #include "utils.h"
31 
32 /*R will be stored in A */
33 int8_t qr_householder_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  matrix_t diag_R_vec[n];
42 
43  for (uint8_t k = 0; k < n; k++) {
44  //compute the two-norm of the k-th column of the matrix A
45  double norm_2 = 0.0;
46  for (uint8_t i = k; i < m; i++) {
47  norm_2 = utils_get_save_square_root(norm_2, A[i][k]);
48 
49  }
50 
51  if (norm_2 != 0.0) {
52  // Form k-th Householder vector.
53  if (A[k][k] < 0) {
54  norm_2 = -norm_2;
55  }
56  for (int i = k; i < m; i++) {
57  A[i][k] /= norm_2;
58  }
59 
60  A[k][k] += 1.0;
61 
62  // Apply the Householder transformation to the rest columns of the matrix.
63  for (uint8_t j = k + 1; j < n; j++) {
64  double sigma = 0.0;
65  for (uint8_t i = k; i < m; i++) {
66  sigma += A[i][k] * A[i][j];
67  }
68  sigma = -sigma / A[k][k];
69  for (uint8_t i = k; i < m; i++) {
70  A[i][j] += sigma * A[i][k];
71  }
72  }
73 
74  } // if
75 
76  diag_R_vec[k] = -norm_2;
77 
78  } //for
79 
80  // Build Q
81  uint8_t max_n, max_m;
82  if (!reduced) {
83  max_n = m;
84  max_m = m;
85  }
86  else {
87  max_n = n;
88  max_m = n;
89  }
90 
91  for (int k = max_n - 1; k >= 0; k--) {
92  for (int i = 0; i < m; i++) {
93  Q[i][k] = 0.0;
94  } // clear matrix Q
95  Q[k][k] = 1.0;
96  for (int j = k; j < max_n; j++) {
97  if (k < n) {
98  if (A[k][k] != 0) {
99  double s = 0.0;
100  for (int i = k; i < m; i++) {
101  s += A[i][k] * Q[i][j];
102  }
103  s = -s / A[k][k];
104  for (int i = k; i < m; i++) {
105  Q[i][j] += s * A[i][k];
106  }
107  } //if
108  } //if
109  }
110  }
111 
112  // Build R
113  for (int i = 0; i < max_m; i++) {
114  for (int j = 0; j < n; j++) {
115  if (i == j) {
116  A[i][j] = diag_R_vec[i];
117  }
118  else if (i > j) {
119  A[i][j] = 0;
120  }
121  }
122  }
123 
124  return 1;
125 }
matrix.h
Matrix computations.
utils.h
Utilities for linear algebra.
qr_householder_decomp
int8_t qr_householder_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 Householder algorithm.
Definition: qr_householder.c:33
matrix_t
#define matrix_t
Define the data type of the matrix elements.
Definition: matrix.h:38
utils_get_save_square_root
double utils_get_save_square_root(double x, double y)
Compute the square root without under/overflow.
Definition: utils.c:168