| // This file is part of Eigen, a lightweight C++ template library | |
| // for linear algebra. | |
| // | |
| // Copyright (C) 2024 The Eigen Authors | |
| // | |
| // This Source Code Form is subject to the terms of the Mozilla | |
| // Public License v. 2.0. If a copy of the MPL was not distributed | |
| // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. | |
| // SPDX-License-Identifier: MPL-2.0 | |
| #include "lapack_common.h" | |
| // Architectural Note: | |
| // Eigen's internal Householder routines (Eigen/src/Householder/) implement a subset of | |
| // operations (primarily Left / Forward / Columnwise reflectors). LAPACK routines (larfg, | |
| // larf, larft, larfb) require the complete Cartesian product of side (Left/Right) x | |
| // trans (NoTrans/ConjTrans) x direct (Forward/Backward) x storev (Columnwise/Rowwise), | |
| // as well as Netlib-compliant scaling (safmin/xlamch) and reflector sign conventions. | |
| // They are implemented natively here with zero dynamic heap allocations. | |
| // LARFG generates an elementary reflector H of order n, such that | |
| // H * [alpha, x]^T = [beta, 0]^T, H^* * H = I. | |
| EIGEN_LAPACK_FUNC(larfg)(int* n, RealScalar* palpha, RealScalar* px, int* incx, RealScalar* ptau) { | |
| Scalar* alpha = reinterpret_cast<Scalar*>(palpha); | |
| Scalar* x = reinterpret_cast<Scalar*>(px); | |
| Scalar* tau = reinterpret_cast<Scalar*>(ptau); | |
| if (*n <= 0) { | |
| *tau = Scalar(0); | |
| return; | |
| } | |
| int m = *n - 1; | |
| RealScalar xnorm = RealScalar(0); | |
| if (m > 0) { | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, 1>, 0, Eigen::InnerStride<Eigen::Dynamic> > xvec( | |
| x, m, Eigen::InnerStride<Eigen::Dynamic>(*incx)); | |
| xnorm = xvec.stableNorm(); | |
| } | |
| #if !ISCOMPLEX | |
| if (m == 0 || xnorm == RealScalar(0)) { | |
| *tau = Scalar(0); | |
| return; | |
| } | |
| RealScalar alphar = *alpha; | |
| RealScalar beta = -Eigen::numext::copysign(Eigen::internal::xlapy2(alphar, xnorm), alphar); | |
| RealScalar safmin = Eigen::internal::xlamch<RealScalar>('S') / Eigen::internal::xlamch<RealScalar>('E'); | |
| int knt = 0; | |
| if (Eigen::numext::abs(beta) < safmin) { | |
| RealScalar rsafmn = RealScalar(1) / safmin; | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, 1>, 0, Eigen::InnerStride<Eigen::Dynamic> > xvec( | |
| x, m, Eigen::InnerStride<Eigen::Dynamic>(*incx)); | |
| do { | |
| ++knt; | |
| xvec *= rsafmn; | |
| beta *= rsafmn; | |
| alphar *= rsafmn; | |
| } while (Eigen::numext::abs(beta) < safmin); | |
| xnorm = xvec.stableNorm(); | |
| beta = -Eigen::numext::copysign(Eigen::internal::xlapy2(alphar, xnorm), alphar); | |
| } | |
| *tau = (beta - alphar) / beta; | |
| { | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, 1>, 0, Eigen::InnerStride<Eigen::Dynamic> > xvec( | |
| x, m, Eigen::InnerStride<Eigen::Dynamic>(*incx)); | |
| xvec *= (RealScalar(1) / (alphar - beta)); | |
| } | |
| for (int j = 0; j < knt; ++j) { | |
| beta *= safmin; | |
| } | |
| *alpha = beta; | |
| #else | |
| RealScalar alphr = Eigen::numext::real(*alpha); | |
| RealScalar alphi = Eigen::numext::imag(*alpha); | |
| if (Eigen::numext::is_exactly_zero(xnorm) && Eigen::numext::is_exactly_zero(alphi)) { | |
| *tau = Scalar(0); | |
| return; | |
| } | |
| RealScalar beta = -Eigen::numext::copysign(Eigen::internal::xlapy3(alphr, alphi, xnorm), alphr); | |
| RealScalar safmin = Eigen::internal::xlamch<RealScalar>('S') / Eigen::internal::xlamch<RealScalar>('E'); | |
| RealScalar rsafmn = RealScalar(1) / safmin; | |
| int knt = 0; | |
| if (Eigen::numext::abs(beta) < safmin) { | |
| do { | |
| ++knt; | |
| if (m > 0) { | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, 1>, 0, Eigen::InnerStride<Eigen::Dynamic> > xvec( | |
| x, m, Eigen::InnerStride<Eigen::Dynamic>(*incx)); | |
| xvec *= rsafmn; | |
| } | |
| beta *= rsafmn; | |
| alphi *= rsafmn; | |
| alphr *= rsafmn; | |
| } while (Eigen::numext::abs(beta) < safmin); | |
| if (m > 0) { | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, 1>, 0, Eigen::InnerStride<Eigen::Dynamic> > xvec( | |
| x, m, Eigen::InnerStride<Eigen::Dynamic>(*incx)); | |
| xnorm = xvec.stableNorm(); | |
| } | |
| *alpha = Scalar(alphr, alphi); | |
| beta = -Eigen::numext::copysign(Eigen::internal::xlapy3(alphr, alphi, xnorm), alphr); | |
| } | |
| *tau = Scalar((beta - alphr) / beta, -alphi / beta); | |
| Scalar denom = *alpha - beta; | |
| // For consistency with LAPACK (CLARFG/ZLARFG), compute 1 / (alpha - beta) via xladiv | |
| // (Smith's algorithm) to avoid intermediate overflow when scaling the reflector vector. | |
| Scalar scale = Eigen::internal::xladiv(Scalar(1), denom); | |
| if (m > 0) { | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, 1>, 0, Eigen::InnerStride<Eigen::Dynamic> > xvec( | |
| x, m, Eigen::InnerStride<Eigen::Dynamic>(*incx)); | |
| xvec *= scale; | |
| } | |
| for (int j = 0; j < knt; ++j) { | |
| beta *= safmin; | |
| } | |
| *alpha = Scalar(beta, RealScalar(0)); | |
| #endif | |
| } | |
| // LARF applies an elementary reflector H = I - tau * v * v^* to matrix C from either left or right. | |
| EIGEN_LAPACK_FUNC(larf) | |
| (const char* side, int* m, int* n, RealScalar* pv, int* incv, RealScalar* ptau, RealScalar* pc, int* ldc, | |
| RealScalar* pwork) { | |
| if (*m <= 0 || *n <= 0) return; | |
| Scalar* v = reinterpret_cast<Scalar*>(pv); | |
| Scalar* tau = reinterpret_cast<Scalar*>(ptau); | |
| Scalar* c = reinterpret_cast<Scalar*>(pc); | |
| Scalar* work = reinterpret_cast<Scalar*>(pwork); | |
| if (*tau == Scalar(0)) return; | |
| bool applyleft = (SIDE(*side) == LEFT); | |
| int lastv = applyleft ? *m : *n; | |
| if (*incv > 0) { | |
| int idx = (lastv - 1) * *incv; | |
| while (lastv > 0 && v[idx] == Scalar(0)) { | |
| --lastv; | |
| idx -= *incv; | |
| } | |
| } else { | |
| int idx = 0; | |
| while (lastv > 0 && v[idx] == Scalar(0)) { | |
| --lastv; | |
| idx -= *incv; | |
| } | |
| } | |
| if (lastv == 0) return; | |
| int lastc = 0; | |
| if (applyleft) { | |
| lastc = Eigen::internal::ilaxlc(lastv, *n, c, *ldc); | |
| } else { | |
| lastc = Eigen::internal::ilaxlr(*m, lastv, c, *ldc); | |
| } | |
| if (lastc == 0) return; | |
| if (applyleft) { | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > C_sub( | |
| c, lastv, lastc, Eigen::OuterStride<>(*ldc)); | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, 1> > w(work, lastc); | |
| if (*incv > 0) { | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, 1>, 0, Eigen::InnerStride<Eigen::Dynamic> > v_sub( | |
| v, lastv, Eigen::InnerStride<Eigen::Dynamic>(*incv)); | |
| w.noalias() = C_sub.adjoint() * v_sub; | |
| C_sub.noalias() -= (*tau) * (v_sub * w.adjoint()); | |
| } else { | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, 1>, 0, Eigen::InnerStride<Eigen::Dynamic> > v_fwd( | |
| v, lastv, Eigen::InnerStride<Eigen::Dynamic>(-*incv)); | |
| auto v_sub = v_fwd.reverse(); | |
| w.noalias() = C_sub.adjoint() * v_sub; | |
| C_sub.noalias() -= (*tau) * (v_sub * w.adjoint()); | |
| } | |
| } else { | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > C_sub( | |
| c, lastc, lastv, Eigen::OuterStride<>(*ldc)); | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, 1> > w(work, lastc); | |
| if (*incv > 0) { | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, 1>, 0, Eigen::InnerStride<Eigen::Dynamic> > v_sub( | |
| v, lastv, Eigen::InnerStride<Eigen::Dynamic>(*incv)); | |
| w.noalias() = C_sub * v_sub; | |
| C_sub.noalias() -= (*tau) * (w * v_sub.adjoint()); | |
| } else { | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, 1>, 0, Eigen::InnerStride<Eigen::Dynamic> > v_fwd( | |
| v, lastv, Eigen::InnerStride<Eigen::Dynamic>(-*incv)); | |
| auto v_sub = v_fwd.reverse(); | |
| w.noalias() = C_sub * v_sub; | |
| C_sub.noalias() -= (*tau) * (w * v_sub.adjoint()); | |
| } | |
| } | |
| } | |
| // LARFT forms the triangular factor T of a block reflector H of order n. | |
| EIGEN_LAPACK_FUNC(larft) | |
| (const char* direct, const char* storev, int* n, int* k, RealScalar* pv, int* ldv, RealScalar* ptau, RealScalar* pt, | |
| int* ldt) { | |
| if (*n == 0 || *k == 0) return; | |
| const Scalar* v = reinterpret_cast<const Scalar*>(pv); | |
| Scalar* tau = reinterpret_cast<Scalar*>(ptau); | |
| Scalar* t = reinterpret_cast<Scalar*>(pt); | |
| bool forward = (*direct == 'F' || *direct == 'f'); | |
| bool colwise = (*storev == 'C' || *storev == 'c'); | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > T_mat( | |
| t, *k, *k, Eigen::OuterStride<>(*ldt)); | |
| if (forward) { | |
| int prevlastv = *n; | |
| for (int i = 0; i < *k; ++i) { | |
| prevlastv = Eigen::numext::maxi(prevlastv, i + 1); | |
| if (tau[i] == Scalar(0)) { | |
| for (int j = 0; j <= i; ++j) { | |
| T_mat(j, i) = Scalar(0); | |
| } | |
| } else { | |
| int lastv = i; | |
| if (colwise) { | |
| for (lastv = *n - 1; lastv >= i + 1; --lastv) { | |
| if (v[lastv + i * (*ldv)] != Scalar(0)) break; | |
| } | |
| int j = Eigen::numext::mini(lastv, prevlastv - 1); | |
| if (i > 0) { | |
| // Reflectors v(i, 0..i-1) have implicit unit diagonal at row i. | |
| // Split dot product: V_prev(i, :)^H * 1 + V_prev(i+1:j, :)^H * V_cur(i+1:j) | |
| // avoiding mutation of v[i + i*ldv]. | |
| Eigen::Map<const Eigen::Matrix<Scalar, 1, Eigen::Dynamic>, 0, Eigen::InnerStride<Eigen::Dynamic> > | |
| V_prev_row(&v[i], i, Eigen::InnerStride<Eigen::Dynamic>(*ldv)); | |
| T_mat.col(i).head(i).noalias() = -tau[i] * V_prev_row.conjugate().transpose(); | |
| if (j > i) { | |
| int len_sub = j - i; | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, | |
| Eigen::OuterStride<> > | |
| V_prev_sub(&v[i + 1], len_sub, i, Eigen::OuterStride<>(*ldv)); | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, 1> > V_cur_sub(&v[(i + 1) + i * (*ldv)], len_sub); | |
| T_mat.col(i).head(i).noalias() += -tau[i] * (V_prev_sub.adjoint() * V_cur_sub); | |
| } | |
| } | |
| } else { | |
| for (lastv = *n - 1; lastv >= i + 1; --lastv) { | |
| if (v[i + lastv * (*ldv)] != Scalar(0)) break; | |
| } | |
| int j = Eigen::numext::mini(lastv, prevlastv - 1); | |
| if (i > 0) { | |
| // Reflectors v(0..i-1, i) have implicit unit diagonal at col i. | |
| // Split dot product: V_prev(:, i) * 1 + V_prev(:, i+1:j) * V_cur(i+1:j)^* | |
| // avoiding mutation of v[i + i*ldv]. | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, 1> > V_prev_col(&v[i * (*ldv)], i); | |
| T_mat.col(i).head(i).noalias() = -tau[i] * V_prev_col; | |
| if (j > i) { | |
| int len_sub = j - i; | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, | |
| Eigen::OuterStride<> > | |
| V_prev_sub(&v[(i + 1) * (*ldv)], i, len_sub, Eigen::OuterStride<>(*ldv)); | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, 1>, 0, Eigen::InnerStride<Eigen::Dynamic> > | |
| V_cur_sub(&v[i + (i + 1) * (*ldv)], len_sub, Eigen::InnerStride<Eigen::Dynamic>(*ldv)); | |
| T_mat.col(i).head(i).noalias() += -tau[i] * (V_prev_sub * V_cur_sub.conjugate()); | |
| } | |
| } | |
| } | |
| if (i > 0) { | |
| // In-place triangular update: T(0:i-1, i) = T(0:i-1, 0:i-1) * T(0:i-1, i) | |
| for (int r = 0; r < i; ++r) { | |
| Scalar sum = Scalar(0); | |
| for (int c = r; c < i; ++c) { | |
| sum += T_mat(r, c) * T_mat(c, i); | |
| } | |
| T_mat(r, i) = sum; | |
| } | |
| } | |
| T_mat(i, i) = tau[i]; | |
| if (i > 0) { | |
| prevlastv = Eigen::numext::maxi(prevlastv, lastv + 1); | |
| } else { | |
| prevlastv = lastv + 1; | |
| } | |
| } | |
| } | |
| } else { | |
| int prevlastv = 0; | |
| for (int i = *k - 1; i >= 0; --i) { | |
| if (tau[i] == Scalar(0)) { | |
| for (int j = i; j < *k; ++j) { | |
| T_mat(j, i) = Scalar(0); | |
| } | |
| } else { | |
| if (i < *k - 1) { | |
| int lastv = 0; | |
| int num_vecs = *k - 1 - i; | |
| if (colwise) { | |
| int diag_row = *n - *k + i; | |
| for (lastv = 0; lastv < i; ++lastv) { | |
| if (v[lastv + i * (*ldv)] != Scalar(0)) break; | |
| } | |
| int j = Eigen::numext::maxi(lastv, prevlastv); | |
| // Reflectors v(j..diag_row, i+1..k-1) have implicit unit diagonal at row diag_row. | |
| // Split dot product: V_rest(diag_row, :)^H * 1 + V_rest(j:diag_row-1, :)^H * V_cur(j:diag_row-1) | |
| // avoiding mutation of v[diag_row + i*ldv]. | |
| Eigen::Map<const Eigen::Matrix<Scalar, 1, Eigen::Dynamic>, 0, Eigen::InnerStride<Eigen::Dynamic> > | |
| V_rest_row(&v[diag_row + (i + 1) * (*ldv)], num_vecs, Eigen::InnerStride<Eigen::Dynamic>(*ldv)); | |
| T_mat.col(i).tail(num_vecs).noalias() = -tau[i] * V_rest_row.conjugate().transpose(); | |
| if (diag_row > j) { | |
| int len_sub = diag_row - j; | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, | |
| Eigen::OuterStride<> > | |
| V_rest_sub(&v[j + (i + 1) * (*ldv)], len_sub, num_vecs, Eigen::OuterStride<>(*ldv)); | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, 1> > V_cur_sub(&v[j + i * (*ldv)], len_sub); | |
| T_mat.col(i).tail(num_vecs).noalias() += -tau[i] * (V_rest_sub.adjoint() * V_cur_sub); | |
| } | |
| } else { | |
| int diag_col = *n - *k + i; | |
| for (lastv = 0; lastv < i; ++lastv) { | |
| if (v[i + lastv * (*ldv)] != Scalar(0)) break; | |
| } | |
| int j = Eigen::numext::maxi(lastv, prevlastv); | |
| // Reflectors v(i+1..k-1, j..diag_col) have implicit unit diagonal at col diag_col. | |
| // Split dot product: V_rest(:, diag_col) * 1 + V_rest(:, j:diag_col-1) * V_cur(j:diag_col-1)^* | |
| // avoiding mutation of v[i + diag_col*ldv]. | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, 1> > V_rest_col(&v[(i + 1) + diag_col * (*ldv)], | |
| num_vecs); | |
| T_mat.col(i).tail(num_vecs).noalias() = -tau[i] * V_rest_col; | |
| if (diag_col > j) { | |
| int len_sub = diag_col - j; | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, | |
| Eigen::OuterStride<> > | |
| V_rest_sub(&v[(i + 1) + j * (*ldv)], num_vecs, len_sub, Eigen::OuterStride<>(*ldv)); | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, 1>, 0, Eigen::InnerStride<Eigen::Dynamic> > | |
| V_cur_sub(&v[i + j * (*ldv)], len_sub, Eigen::InnerStride<Eigen::Dynamic>(*ldv)); | |
| T_mat.col(i).tail(num_vecs).noalias() += -tau[i] * (V_rest_sub * V_cur_sub.conjugate()); | |
| } | |
| } | |
| // In-place triangular update: T(i+1:k-1, i) = T(i+1:k-1, i+1:k-1) * T(i+1:k-1, i) | |
| for (int r = num_vecs - 1; r >= 0; --r) { | |
| int row = (i + 1) + r; | |
| Scalar sum = Scalar(0); | |
| for (int c = 0; c <= r; ++c) { | |
| sum += T_mat(row, (i + 1) + c) * T_mat((i + 1) + c, i); | |
| } | |
| T_mat(row, i) = sum; | |
| } | |
| if (i < *k - 1) { | |
| prevlastv = Eigen::numext::mini(prevlastv, lastv); | |
| } else { | |
| prevlastv = lastv; | |
| } | |
| } | |
| T_mat(i, i) = tau[i]; | |
| } | |
| } | |
| } | |
| } | |
| namespace { | |
| // In-place triangular matrix multiplications: W = W * Tri | |
| // Updates W column-by-column directly within caller-supplied work buffer with zero heap allocations. | |
| template <bool UnitDiag, typename MatrixW, typename MatrixTri> | |
| void apply_triangular_right_upper(MatrixW& W, const MatrixTri& U) { | |
| const Eigen::Index k = W.cols(); | |
| for (Eigen::Index j = k - 1; j >= 1; --j) { | |
| if (!UnitDiag) { | |
| W.col(j) *= U(j, j); | |
| } | |
| W.col(j).noalias() += W.leftCols(j) * U.col(j).head(j); | |
| } | |
| if (!UnitDiag && k > 0) { | |
| W.col(0) *= U(0, 0); | |
| } | |
| } | |
| template <bool UnitDiag, typename MatrixW, typename MatrixTri> | |
| void apply_triangular_right_lower(MatrixW& W, const MatrixTri& L) { | |
| const Eigen::Index k = W.cols(); | |
| for (Eigen::Index j = 0; j < k - 1; ++j) { | |
| const Eigen::Index tail_len = k - 1 - j; | |
| if (!UnitDiag) { | |
| W.col(j) *= L(j, j); | |
| } | |
| W.col(j).noalias() += W.rightCols(tail_len) * L.col(j).tail(tail_len); | |
| } | |
| if (!UnitDiag && k > 0) { | |
| W.col(k - 1) *= L(k - 1, k - 1); | |
| } | |
| } | |
| template <bool UnitDiag, typename MatrixW, typename MatrixTri> | |
| void apply_triangular_right_upper_adjoint(MatrixW& W, const MatrixTri& U) { | |
| // U is Upper Triangular, so U^H is Lower Triangular. | |
| const Eigen::Index k = W.cols(); | |
| for (Eigen::Index j = 0; j < k - 1; ++j) { | |
| const Eigen::Index tail_len = k - 1 - j; | |
| if (!UnitDiag) { | |
| W.col(j) *= Eigen::numext::conj(U(j, j)); | |
| } | |
| W.col(j).noalias() += W.rightCols(tail_len) * U.row(j).tail(tail_len).adjoint(); | |
| } | |
| if (!UnitDiag && k > 0) { | |
| W.col(k - 1) *= Eigen::numext::conj(U(k - 1, k - 1)); | |
| } | |
| } | |
| template <bool UnitDiag, typename MatrixW, typename MatrixTri> | |
| void apply_triangular_right_lower_adjoint(MatrixW& W, const MatrixTri& L) { | |
| // L is Lower Triangular, so L^H is Upper Triangular. | |
| const Eigen::Index k = W.cols(); | |
| for (Eigen::Index j = k - 1; j >= 1; --j) { | |
| if (!UnitDiag) { | |
| W.col(j) *= Eigen::numext::conj(L(j, j)); | |
| } | |
| W.col(j).noalias() += W.leftCols(j) * L.row(j).head(j).adjoint(); | |
| } | |
| if (!UnitDiag && k > 0) { | |
| W.col(0) *= Eigen::numext::conj(L(0, 0)); | |
| } | |
| } | |
| } // namespace | |
| // LARFB applies a block reflector H or its transpose/adjoint to matrix C. | |
| EIGEN_LAPACK_FUNC(larfb) | |
| (const char* side, const char* trans, const char* direct, const char* storev, int* m, int* n, int* k, RealScalar* pv, | |
| int* ldv, RealScalar* pt, int* ldt, RealScalar* pc, int* ldc, RealScalar* pwork, int* ldwork) { | |
| if (*m <= 0 || *n <= 0 || *k <= 0) return; | |
| const Scalar* v = reinterpret_cast<const Scalar*>(pv); | |
| const Scalar* t = reinterpret_cast<const Scalar*>(pt); | |
| Scalar* c = reinterpret_cast<Scalar*>(pc); | |
| Scalar* work = reinterpret_cast<Scalar*>(pwork); | |
| bool applyleft = (SIDE(*side) == LEFT); | |
| bool notran = (OP(*trans) == NOTR); | |
| bool forward = (*direct == 'F' || *direct == 'f'); | |
| bool colwise = (*storev == 'C' || *storev == 'c'); | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > | |
| T_mat(t, *k, *k, Eigen::OuterStride<>(*ldt)); | |
| int lastv = 0; | |
| int lastc = 0; | |
| if (colwise) { | |
| if (forward) { | |
| if (applyleft) { | |
| lastv = Eigen::numext::maxi(*k, Eigen::internal::ilaxlr(*m, *k, v, *ldv)); | |
| lastc = Eigen::internal::ilaxlc(lastv, *n, c, *ldc); | |
| if (lastc == 0) return; | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > W( | |
| work, lastc, *k, Eigen::OuterStride<>(*ldwork)); | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > C1( | |
| c, *k, lastc, Eigen::OuterStride<>(*ldc)); | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, | |
| Eigen::OuterStride<> > | |
| V1(v, *k, *k, Eigen::OuterStride<>(*ldv)); | |
| W.noalias() = C1.adjoint() * V1.template triangularView<Eigen::UnitLower>(); | |
| if (lastv > *k) { | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > | |
| C2(&c[*k], lastv - *k, lastc, Eigen::OuterStride<>(*ldc)); | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, | |
| Eigen::OuterStride<> > | |
| V2(&v[*k], lastv - *k, *k, Eigen::OuterStride<>(*ldv)); | |
| W.noalias() += C2.adjoint() * V2; | |
| } | |
| if (notran) { | |
| apply_triangular_right_upper_adjoint<false>(W, T_mat); | |
| } else { | |
| apply_triangular_right_upper<false>(W, T_mat); | |
| } | |
| if (lastv > *k) { | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > | |
| C2(&c[*k], lastv - *k, lastc, Eigen::OuterStride<>(*ldc)); | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, | |
| Eigen::OuterStride<> > | |
| V2(&v[*k], lastv - *k, *k, Eigen::OuterStride<>(*ldv)); | |
| C2.noalias() -= V2 * W.adjoint(); | |
| } | |
| apply_triangular_right_lower_adjoint<true>(W, V1); | |
| C1.noalias() -= W.adjoint(); | |
| } else { | |
| lastv = Eigen::numext::maxi(*k, Eigen::internal::ilaxlr(*n, *k, v, *ldv)); | |
| lastc = Eigen::internal::ilaxlr(*m, lastv, c, *ldc); | |
| if (lastc == 0) return; | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > W( | |
| work, lastc, *k, Eigen::OuterStride<>(*ldwork)); | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > C1( | |
| c, lastc, *k, Eigen::OuterStride<>(*ldc)); | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, | |
| Eigen::OuterStride<> > | |
| V1(v, *k, *k, Eigen::OuterStride<>(*ldv)); | |
| W.noalias() = C1 * V1.template triangularView<Eigen::UnitLower>(); | |
| if (lastv > *k) { | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > | |
| C2(&c[*k * (*ldc)], lastc, lastv - *k, Eigen::OuterStride<>(*ldc)); | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, | |
| Eigen::OuterStride<> > | |
| V2(&v[*k], lastv - *k, *k, Eigen::OuterStride<>(*ldv)); | |
| W.noalias() += C2 * V2; | |
| } | |
| if (notran) { | |
| apply_triangular_right_upper<false>(W, T_mat); | |
| } else { | |
| apply_triangular_right_upper_adjoint<false>(W, T_mat); | |
| } | |
| if (lastv > *k) { | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > | |
| C2(&c[*k * (*ldc)], lastc, lastv - *k, Eigen::OuterStride<>(*ldc)); | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, | |
| Eigen::OuterStride<> > | |
| V2(&v[*k], lastv - *k, *k, Eigen::OuterStride<>(*ldv)); | |
| C2.noalias() -= W * V2.adjoint(); | |
| } | |
| apply_triangular_right_lower_adjoint<true>(W, V1); | |
| C1.noalias() -= W; | |
| } | |
| } else { | |
| // colwise && !forward | |
| if (applyleft) { | |
| lastv = Eigen::numext::maxi(*k, Eigen::internal::ilaxlr(*m, *k, v, *ldv)); | |
| lastc = Eigen::internal::ilaxlc(lastv, *n, c, *ldc); | |
| if (lastc == 0) return; | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > W( | |
| work, lastc, *k, Eigen::OuterStride<>(*ldwork)); | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > C2( | |
| &c[lastv - *k], *k, lastc, Eigen::OuterStride<>(*ldc)); | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, | |
| Eigen::OuterStride<> > | |
| V2(&v[lastv - *k], *k, *k, Eigen::OuterStride<>(*ldv)); | |
| W.noalias() = C2.adjoint() * V2.template triangularView<Eigen::UnitUpper>(); | |
| if (lastv > *k) { | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > | |
| C1(c, lastv - *k, lastc, Eigen::OuterStride<>(*ldc)); | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, | |
| Eigen::OuterStride<> > | |
| V1(v, lastv - *k, *k, Eigen::OuterStride<>(*ldv)); | |
| W.noalias() += C1.adjoint() * V1; | |
| } | |
| if (notran) { | |
| apply_triangular_right_lower_adjoint<false>(W, T_mat); | |
| } else { | |
| apply_triangular_right_lower<false>(W, T_mat); | |
| } | |
| if (lastv > *k) { | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > | |
| C1(c, lastv - *k, lastc, Eigen::OuterStride<>(*ldc)); | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, | |
| Eigen::OuterStride<> > | |
| V1(v, lastv - *k, *k, Eigen::OuterStride<>(*ldv)); | |
| C1.noalias() -= V1 * W.adjoint(); | |
| } | |
| apply_triangular_right_upper_adjoint<true>(W, V2); | |
| C2.noalias() -= W.adjoint(); | |
| } else { | |
| lastv = Eigen::numext::maxi(*k, Eigen::internal::ilaxlr(*n, *k, v, *ldv)); | |
| lastc = Eigen::internal::ilaxlr(*m, lastv, c, *ldc); | |
| if (lastc == 0) return; | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > W( | |
| work, lastc, *k, Eigen::OuterStride<>(*ldwork)); | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > C2( | |
| &c[(lastv - *k) * (*ldc)], lastc, *k, Eigen::OuterStride<>(*ldc)); | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, | |
| Eigen::OuterStride<> > | |
| V2(&v[lastv - *k], *k, *k, Eigen::OuterStride<>(*ldv)); | |
| W.noalias() = C2 * V2.template triangularView<Eigen::UnitUpper>(); | |
| if (lastv > *k) { | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > | |
| C1(c, lastc, lastv - *k, Eigen::OuterStride<>(*ldc)); | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, | |
| Eigen::OuterStride<> > | |
| V1(v, lastv - *k, *k, Eigen::OuterStride<>(*ldv)); | |
| W.noalias() += C1 * V1; | |
| } | |
| if (notran) { | |
| apply_triangular_right_lower<false>(W, T_mat); | |
| } else { | |
| apply_triangular_right_lower_adjoint<false>(W, T_mat); | |
| } | |
| if (lastv > *k) { | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > | |
| C1(c, lastc, lastv - *k, Eigen::OuterStride<>(*ldc)); | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, | |
| Eigen::OuterStride<> > | |
| V1(v, lastv - *k, *k, Eigen::OuterStride<>(*ldv)); | |
| C1.noalias() -= W * V1.adjoint(); | |
| } | |
| apply_triangular_right_upper_adjoint<true>(W, V2); | |
| C2.noalias() -= W; | |
| } | |
| } | |
| } else { | |
| // storev == 'R' | |
| if (forward) { | |
| if (applyleft) { | |
| lastv = Eigen::numext::maxi(*k, Eigen::internal::ilaxlc(*k, *m, v, *ldv)); | |
| lastc = Eigen::internal::ilaxlc(lastv, *n, c, *ldc); | |
| if (lastc == 0) return; | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > W( | |
| work, lastc, *k, Eigen::OuterStride<>(*ldwork)); | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > C1( | |
| c, *k, lastc, Eigen::OuterStride<>(*ldc)); | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, | |
| Eigen::OuterStride<> > | |
| V1(v, *k, *k, Eigen::OuterStride<>(*ldv)); | |
| W.noalias() = C1.adjoint() * V1.template triangularView<Eigen::UnitUpper>().adjoint(); | |
| if (lastv > *k) { | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > | |
| C2(&c[*k], lastv - *k, lastc, Eigen::OuterStride<>(*ldc)); | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, | |
| Eigen::OuterStride<> > | |
| V2(&v[*k * (*ldv)], *k, lastv - *k, Eigen::OuterStride<>(*ldv)); | |
| W.noalias() += C2.adjoint() * V2.adjoint(); | |
| } | |
| if (notran) { | |
| apply_triangular_right_upper_adjoint<false>(W, T_mat); | |
| } else { | |
| apply_triangular_right_upper<false>(W, T_mat); | |
| } | |
| if (lastv > *k) { | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > | |
| C2(&c[*k], lastv - *k, lastc, Eigen::OuterStride<>(*ldc)); | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, | |
| Eigen::OuterStride<> > | |
| V2(&v[*k * (*ldv)], *k, lastv - *k, Eigen::OuterStride<>(*ldv)); | |
| C2.noalias() -= V2.adjoint() * W.adjoint(); | |
| } | |
| apply_triangular_right_upper<true>(W, V1); | |
| C1.noalias() -= W.adjoint(); | |
| } else { | |
| lastv = Eigen::numext::maxi(*k, Eigen::internal::ilaxlc(*k, *n, v, *ldv)); | |
| lastc = Eigen::internal::ilaxlr(*m, lastv, c, *ldc); | |
| if (lastc == 0) return; | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > W( | |
| work, lastc, *k, Eigen::OuterStride<>(*ldwork)); | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > C1( | |
| c, lastc, *k, Eigen::OuterStride<>(*ldc)); | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, | |
| Eigen::OuterStride<> > | |
| V1(v, *k, *k, Eigen::OuterStride<>(*ldv)); | |
| W.noalias() = C1 * V1.template triangularView<Eigen::UnitUpper>().adjoint(); | |
| if (lastv > *k) { | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > | |
| C2(&c[*k * (*ldc)], lastc, lastv - *k, Eigen::OuterStride<>(*ldc)); | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, | |
| Eigen::OuterStride<> > | |
| V2(&v[*k * (*ldv)], *k, lastv - *k, Eigen::OuterStride<>(*ldv)); | |
| W.noalias() += C2 * V2.adjoint(); | |
| } | |
| if (notran) { | |
| apply_triangular_right_upper<false>(W, T_mat); | |
| } else { | |
| apply_triangular_right_upper_adjoint<false>(W, T_mat); | |
| } | |
| if (lastv > *k) { | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > | |
| C2(&c[*k * (*ldc)], lastc, lastv - *k, Eigen::OuterStride<>(*ldc)); | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, | |
| Eigen::OuterStride<> > | |
| V2(&v[*k * (*ldv)], *k, lastv - *k, Eigen::OuterStride<>(*ldv)); | |
| C2.noalias() -= W * V2; | |
| } | |
| apply_triangular_right_upper<true>(W, V1); | |
| C1.noalias() -= W; | |
| } | |
| } else { | |
| // storev == 'R' && !forward | |
| if (applyleft) { | |
| lastv = Eigen::numext::maxi(*k, Eigen::internal::ilaxlc(*k, *m, v, *ldv)); | |
| lastc = Eigen::internal::ilaxlc(lastv, *n, c, *ldc); | |
| if (lastc == 0) return; | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > W( | |
| work, lastc, *k, Eigen::OuterStride<>(*ldwork)); | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > C2( | |
| &c[lastv - *k], *k, lastc, Eigen::OuterStride<>(*ldc)); | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, | |
| Eigen::OuterStride<> > | |
| V2(&v[(lastv - *k) * (*ldv)], *k, *k, Eigen::OuterStride<>(*ldv)); | |
| W.noalias() = C2.adjoint() * V2.template triangularView<Eigen::UnitLower>().adjoint(); | |
| if (lastv > *k) { | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > | |
| C1(c, lastv - *k, lastc, Eigen::OuterStride<>(*ldc)); | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, | |
| Eigen::OuterStride<> > | |
| V1(v, *k, lastv - *k, Eigen::OuterStride<>(*ldv)); | |
| W.noalias() += C1.adjoint() * V1.adjoint(); | |
| } | |
| if (notran) { | |
| apply_triangular_right_lower_adjoint<false>(W, T_mat); | |
| } else { | |
| apply_triangular_right_lower<false>(W, T_mat); | |
| } | |
| if (lastv > *k) { | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > | |
| C1(c, lastv - *k, lastc, Eigen::OuterStride<>(*ldc)); | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, | |
| Eigen::OuterStride<> > | |
| V1(v, *k, lastv - *k, Eigen::OuterStride<>(*ldv)); | |
| C1.noalias() -= V1.adjoint() * W.adjoint(); | |
| } | |
| apply_triangular_right_lower<true>(W, V2); | |
| C2.noalias() -= W.adjoint(); | |
| } else { | |
| lastv = Eigen::numext::maxi(*k, Eigen::internal::ilaxlc(*k, *n, v, *ldv)); | |
| lastc = Eigen::internal::ilaxlr(*m, lastv, c, *ldc); | |
| if (lastc == 0) return; | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > W( | |
| work, lastc, *k, Eigen::OuterStride<>(*ldwork)); | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > C2( | |
| &c[(lastv - *k) * (*ldc)], lastc, *k, Eigen::OuterStride<>(*ldc)); | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, | |
| Eigen::OuterStride<> > | |
| V2(&v[(lastv - *k) * (*ldv)], *k, *k, Eigen::OuterStride<>(*ldv)); | |
| W.noalias() = C2 * V2.template triangularView<Eigen::UnitLower>().adjoint(); | |
| if (lastv > *k) { | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > | |
| C1(c, lastc, lastv - *k, Eigen::OuterStride<>(*ldc)); | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, | |
| Eigen::OuterStride<> > | |
| V1(v, *k, lastv - *k, Eigen::OuterStride<>(*ldv)); | |
| W.noalias() += C1 * V1.adjoint(); | |
| } | |
| if (notran) { | |
| apply_triangular_right_lower<false>(W, T_mat); | |
| } else { | |
| apply_triangular_right_lower_adjoint<false>(W, T_mat); | |
| } | |
| if (lastv > *k) { | |
| Eigen::Map<Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, Eigen::OuterStride<> > | |
| C1(c, lastc, lastv - *k, Eigen::OuterStride<>(*ldc)); | |
| Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>, 0, | |
| Eigen::OuterStride<> > | |
| V1(v, *k, lastv - *k, Eigen::OuterStride<>(*ldv)); | |
| C1.noalias() -= W * V1; | |
| } | |
| apply_triangular_right_lower<true>(W, V2); | |
| C2.noalias() -= W; | |
| } | |
| } | |
| } | |
| } |