MVE - Multi-View Environment mve-devel
Loading...
Searching...
No Matches
matrix_qr.h
Go to the documentation of this file.
1/*
2 * Copyright (C) 2015, Daniel Thuerck, Simon Fuhrmann
3 * TU Darmstadt - Graphics, Capture and Massively Parallel Computing
4 * All rights reserved.
5 *
6 * This software may be modified and distributed under the terms
7 * of the BSD 3-Clause license. See the LICENSE.txt file for details.
8 */
9
10#ifndef MATH_MATRIX_QR_HEADER
11#define MATH_MATRIX_QR_HEADER
12
13#include <vector>
14#include <algorithm>
15#include <iostream>
16
17#include "math/defines.h"
18#include "math/matrix.h"
19
21
29template <typename T>
30void
31matrix_qr (T const* mat_a, int rows, int cols,
32 T* mat_q, T* mat_r, T const& epsilon = T(1e-12));
33
38template <typename T, int M, int N>
39void
40matrix_qr (Matrix<T, M, N> const& mat_a, Matrix<T, M, M>* mat_q,
41 Matrix<T, M, N>* mat_r, T const& epsilon = T(1e-12));
42
44
45/* ------------------------- QR Internals ------------------------- */
46
49
54template <typename T>
55void
56matrix_givens_rotation (T const& alpha, T const& beta,
57 T* givens_c, T* givens_s, T const& epsilon)
58{
59 if (MATH_EPSILON_EQ(beta, T(0), epsilon))
60 {
61 *givens_c = T(1);
62 *givens_s = T(0);
63 return;
64 }
65
66 if (std::abs(beta) > std::abs(alpha))
67 {
68 T tao = -alpha / beta;
69 *givens_s = T(1) / std::sqrt(T(1) + tao * tao);
70 *givens_c = *givens_s * tao;
71 }
72 else
73 {
74 T tao = -beta / alpha;
75 *givens_c = T(1) / std::sqrt(T(1) + tao * tao);
76 *givens_s = *givens_c * tao;
77 }
78}
79
84template <typename T>
85void
86matrix_apply_givens_column (T* mat, int rows, int cols, int givens_i,
87 int givens_k, T const& givens_c, T const& givens_s)
88{
89 for (int j = 0; j < rows; ++j)
90 {
91 T const tao1 = mat[j * cols + givens_i];
92 T const tao2 = mat[j * cols + givens_k];
93 mat[j * cols + givens_i] = givens_c * tao1 - givens_s * tao2;
94 mat[j * cols + givens_k] = givens_s * tao1 + givens_c * tao2;
95 }
96}
97
102template <typename T>
103void
104matrix_apply_givens_row (T* mat, int /*rows*/, int cols, int givens_i,
105 int givens_k, T const& givens_c, T const& givens_s)
106{
107 for (int j = 0; j < cols; ++j)
108 {
109 T const tao1 = mat[givens_i * cols + j];
110 T const tao2 = mat[givens_k * cols + j];
111 mat[givens_i * cols + j] = givens_c * tao1 - givens_s * tao2;
112 mat[givens_k * cols + j] = givens_s * tao1 + givens_c * tao2;
113 }
114}
115
118
119/* ------------------------ Implementation ------------------------ */
120
122
123template <typename T>
124void
125matrix_qr (T const* mat_a, int rows, int cols,
126 T* mat_q, T* mat_r, T const& epsilon)
127{
128 /* Prepare Q and R: Copy A to R, init Q's diagonal. */
129 std::copy(mat_a, mat_a + rows * cols, mat_r);
130 std::fill(mat_q, mat_q + rows * rows, T(0));
131 for (int i = 0; i < rows; ++i)
132 mat_q[i * rows + i] = T(1);
133
134 std::vector<T> matrix_subset(2 * (cols + 1), T(0));
135 for (int j = 0; j < cols; ++j)
136 {
137 for (int i = rows - 1; i > j; --i)
138 {
139 T givens_c, givens_s;
140 internal::matrix_givens_rotation(mat_r[(i - 1) * cols + j],
141 mat_r[i * cols + j], &givens_c, &givens_s, epsilon);
142
143 int submatrix_width = cols - j;
144 for (int k = 0; k < submatrix_width; ++k)
145 {
146 matrix_subset[0 * submatrix_width + k] =
147 mat_r[(i - 1) * cols + (j + k)];
148 matrix_subset[1 * submatrix_width + k] =
149 mat_r[i * cols + (j + k)];
150 }
151
152 /* Only apply Givens rotation on a subblock of the matrix. */
153 for (int k = 0; k < (cols - j); ++k)
154 {
155 mat_r[(i - 1) * cols + (j + k)] =
156 givens_c * matrix_subset[0 * submatrix_width + k] -
157 givens_s * matrix_subset[1 * submatrix_width + k];
158 mat_r[i * cols + (j + k)] =
159 givens_s * matrix_subset[0 * submatrix_width + k] +
160 givens_c * matrix_subset[1 * submatrix_width + k];
161 }
162
163 /* Construct Q matrix by repeated Givens rotations. */
164 internal::matrix_apply_givens_column(mat_q, rows, rows, i - 1, i,
165 givens_c, givens_s);
166 }
167 }
168}
169
170template <typename T, int M, int N>
171void
173 Matrix<T, M, N>* mat_r, T const& epsilon)
174{
175 matrix_qr(mat_a.begin(), M, N, mat_q->begin(), mat_r->begin(), epsilon);
176}
177
179
180#endif /* MATH_MATRIX_QR_HEADER */
Matrix class for arbitrary dimensions and types.
Definition matrix.h:54
T * begin(void)
Definition matrix.h:506
#define MATH_INTERNAL_NAMESPACE_BEGIN
Definition defines.h:24
#define MATH_EPSILON_EQ(x, v, eps)
Definition defines.h:96
#define MATH_NAMESPACE_BEGIN
Definition defines.h:15
#define MATH_NAMESPACE_END
Definition defines.h:16
#define MATH_INTERNAL_NAMESPACE_END
Definition defines.h:25
void matrix_apply_givens_row(T *mat, int, int cols, int givens_i, int givens_k, T const &givens_c, T const &givens_s)
Applies a transposed Givens rotation for rows (givens_i, givens_k) by only rotating the required set ...
Definition matrix_qr.h:104
void matrix_givens_rotation(T const &alpha, T const &beta, T *givens_c, T *givens_s, T const &epsilon)
Calculates the Givens rotation coefficients c and s by solving [alpha beta] [c s;-c s] = [sqrt(alpha^...
Definition matrix_qr.h:56
void matrix_apply_givens_column(T *mat, int rows, int cols, int givens_i, int givens_k, T const &givens_c, T const &givens_s)
Applies a Givens rotation for columns (givens_i, givens_k) by only rotating the required set of colum...
Definition matrix_qr.h:86
void matrix_qr(T const *mat_a, int rows, int cols, T *mat_q, T *mat_r, T const &epsilon=T(1e-12))
Calculates a QR decomposition for a given matrix A.
Definition matrix_qr.h:125