41 Matrix<T, M, N>* mat_r, T
const& epsilon = T(1e-12));
57 T* givens_c, T* givens_s, T
const& epsilon)
66 if (std::abs(beta) > std::abs(alpha))
68 T tao = -alpha / beta;
69 *givens_s = T(1) / std::sqrt(T(1) + tao * tao);
70 *givens_c = *givens_s * tao;
74 T tao = -beta / alpha;
75 *givens_c = T(1) / std::sqrt(T(1) + tao * tao);
76 *givens_s = *givens_c * tao;
126 T* mat_q, T* mat_r, T
const& epsilon)
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);
134 std::vector<T> matrix_subset(2 * (cols + 1), T(0));
135 for (
int j = 0; j < cols; ++j)
137 for (
int i = rows - 1; i > j; --i)
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);
143 int submatrix_width = cols - j;
144 for (
int k = 0; k < submatrix_width; ++k)
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)];
153 for (
int k = 0; k < (cols - j); ++k)
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];
164 internal::matrix_apply_givens_column(mat_q, rows, rows, i - 1, i,
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 ...
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^...
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...
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.