Hessenberg

DESCRIPTION

@@ -6,10 +6,10 @@ Author: Stéphane Laurent
 Maintainer: Stéphane Laurent <laurent_step@outlook.fr>
 Description: Matrix algebra using the 'Eigen' C++ library: determinant,
     rank, inverse, pseudo-inverse, kernel and image, QR decomposition,
-    Cholesky decomposition, Schur decomposition, QZ decomposition, linear
-    least-squares problems. Also provides matrix functions such as
-    exponential, logarithm, power, sine and cosine. Complex matrices are
-    supported, except for the QZ decomposition.
+    Cholesky decomposition, Schur decomposition, QZ decomposition, 
+    Hessenberg decomposition, linear least-squares problems. Also provides 
+    matrix functions such as exponential, logarithm, power, sine and cosine. 
+    Complex matrices are supported, except for the QZ decomposition.
 License: GPL-3
 URL: https://github.com/stla/EigenR
 BugReports: https://github.com/stla/EigenR/issues

NAMESPACE

@@ -1,6 +1,7 @@
 # Generated by roxygen2: do not edit by hand
 
 S3method(print,SparseMatrix)
+export(Eigen_Hessenberg)
 export(Eigen_QR)
 export(Eigen_UtDU)
 export(Eigen_absdet)

NEWS.md

@@ -1,6 +1,7 @@
 # EigenR 1.3.0
 
-New features: Schur decomposition, QZ decomposition (only for real matrices).
+New features: Schur decomposition, QZ decomposition (only for real matrices), 
+Hessenberg decomposition.
 
 
 # EigenR 1.2.3

R/EigenR.R

@@ -755,3 +755,31 @@ Eigen_complexSchur <- function(M) {
     "U" = U[["real"]] + 1i * U[["imag"]]
   )
 }
+
+#' Hessenberg decomposition
+#' @description Hessenberg decomposition of a square matrix.
+#'
+#' @param M real or complex square matrix
+#'
+#' @return A list with the \code{H} and \code{Q} matrices.
+#' @export
+#' @details See \href{https://eigen.tuxfamily.org/dox/classEigen_1_1HessenbergDecomposition.html}{Eigen::HessenbergDecomposition}.
+#'
+#' @examples
+#' library(EigenR)
+#' M <- cbind(c(3, 2i, 1+3i), c(1, 1i, 1), c(5, 0, -2i))
+#' Eigen_Hessenberg(M)
+Eigen_Hessenberg <- function(M) {
+  stopifnot(isSquareMatrix(M), isRealOrComplex(M))
+  if(isReal(M)) {
+    EigenR_Hessenberg_real(M)
+  } else {
+    hd <- EigenR_Hessenberg_cplx(Re(M), Im(M))
+    H <- hd[["H"]]
+    Q <- hd[["Q"]]
+    list(
+      "H" = H[["real"]] + 1i * H[["imag"]],
+      "Q" = Q[["real"]] + 1i * Q[["imag"]]
+    )
+  }
+}

R/RcppExports.R

@@ -49,6 +49,14 @@ EigenR_logabsdet <- function(M) {
     .Call('_EigenR_EigenR_logabsdet', PACKAGE = 'EigenR', M)
 }
 
+EigenR_Hessenberg_real <- function(M) {
+    .Call('_EigenR_EigenR_Hessenberg_real', PACKAGE = 'EigenR', M)
+}
+
+EigenR_Hessenberg_cplx <- function(Re, Im) {
+    .Call('_EigenR_EigenR_Hessenberg_cplx', PACKAGE = 'EigenR', Re, Im)
+}
+
 EigenR_image_LU_real <- function(M) {
     .Call('_EigenR_EigenR_image_LU_real', PACKAGE = 'EigenR', M)
 }

src/Hessenberg.cpp

@@ -0,0 +1,23 @@
+#include "EigenR.h"
+
+// [[Rcpp::export]]
+Rcpp::List EigenR_Hessenberg_real(const Eigen::MatrixXd& M) {
+  Eigen::HessenbergDecomposition<Eigen::MatrixXd> hd(M.rows());
+  hd.compute(M);
+  const Eigen::MatrixXd H = hd.matrixH();
+  const Eigen::MatrixXd Q = hd.matrixQ();
+  return Rcpp::List::create(Rcpp::Named("H") = H,
+                            Rcpp::Named("Q") = Q);
+}
+
+// [[Rcpp::export]]
+Rcpp::List EigenR_Hessenberg_cplx(const Eigen::MatrixXd& Re, 
+                                  const Eigen::MatrixXd& Im) {
+  const Eigen::MatrixXcd M = matricesToMatrixXcd(Re, Im);
+  Eigen::HessenbergDecomposition<Eigen::MatrixXcd> hd(M.rows());
+  hd.compute(M);
+  const Eigen::MatrixXcd H = hd.matrixH();
+  const Eigen::MatrixXcd Q = hd.matrixQ();
+  return Rcpp::List::create(Rcpp::Named("H") = cplxMatrixToList(H),
+                            Rcpp::Named("Q") = cplxMatrixToList(Q));
+}

src/RcppExports.cpp

@@ -162,6 +162,29 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// EigenR_Hessenberg_real
+Rcpp::List EigenR_Hessenberg_real(const Eigen::MatrixXd& M);
+RcppExport SEXP _EigenR_EigenR_Hessenberg_real(SEXP MSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< const Eigen::MatrixXd& >::type M(MSEXP);
+    rcpp_result_gen = Rcpp::wrap(EigenR_Hessenberg_real(M));
+    return rcpp_result_gen;
+END_RCPP
+}
+// EigenR_Hessenberg_cplx
+Rcpp::List EigenR_Hessenberg_cplx(const Eigen::MatrixXd& Re, const Eigen::MatrixXd& Im);
+RcppExport SEXP _EigenR_EigenR_Hessenberg_cplx(SEXP ReSEXP, SEXP ImSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< const Eigen::MatrixXd& >::type Re(ReSEXP);
+    Rcpp::traits::input_parameter< const Eigen::MatrixXd& >::type Im(ImSEXP);
+    rcpp_result_gen = Rcpp::wrap(EigenR_Hessenberg_cplx(Re, Im));
+    return rcpp_result_gen;
+END_RCPP
+}
 // EigenR_image_LU_real
 Eigen::MatrixXd EigenR_image_LU_real(const Eigen::MatrixXd& M);
 RcppExport SEXP _EigenR_EigenR_image_LU_real(SEXP MSEXP) {
@@ -748,6 +771,8 @@ static const R_CallMethodDef CallEntries[] = {
     {"_EigenR_EigenR_det_sparse_cplx", (DL_FUNC) &_EigenR_EigenR_det_sparse_cplx, 5},
     {"_EigenR_EigenR_absdet", (DL_FUNC) &_EigenR_EigenR_absdet, 1},
     {"_EigenR_EigenR_logabsdet", (DL_FUNC) &_EigenR_EigenR_logabsdet, 1},
+    {"_EigenR_EigenR_Hessenberg_real", (DL_FUNC) &_EigenR_EigenR_Hessenberg_real, 1},
+    {"_EigenR_EigenR_Hessenberg_cplx", (DL_FUNC) &_EigenR_EigenR_Hessenberg_cplx, 2},
     {"_EigenR_EigenR_image_LU_real", (DL_FUNC) &_EigenR_EigenR_image_LU_real, 1},
     {"_EigenR_EigenR_image_LU_cplx", (DL_FUNC) &_EigenR_EigenR_image_LU_cplx, 2},
     {"_EigenR_EigenR_image_QR_real", (DL_FUNC) &_EigenR_EigenR_image_QR_real, 1},