QZ

.Rprofile

@@ -4,6 +4,24 @@ dllunload <- function(){
   )
 }
 
-makedoc <- function(){
-  roxygen2::roxygenise(load_code = roxygen2::load_installed)
+myinstall <- function() {
+  try(pkgload::unload("EigenR"))
+  if(rstudioapi::isAvailable()) {
+    rstudioapi::restartSession(
+      "devtools::install(quick = TRUE, keep_source = TRUE)"
+    )
+  } else {
+    #try(dllunload())
+    devtools::install(quick = TRUE, keep_source = TRUE)
+  }
+}
+
+mydocument <- function() {
+  if(rstudioapi::isAvailable()) {
+    rstudioapi::restartSession(
+      "roxygen2::roxygenise(load_code = roxygen2::load_installed)"
+    )
+  } else {
+    roxygen2::roxygenise(load_code = roxygen2::load_installed)
+  }
 }

R/EigenR.R

@@ -673,3 +673,76 @@ Eigen_sqrt <- function(M){
   }
   Msqrt
 }
+
+#' Real QZ decomposition
+#' @description Real QZ decomposition of a pair of square matrices.
+#'
+#' @param A,B real square matrices with the same size
+#'
+#' @return A list with the \code{Q}, \code{Z}, \code{S} and \code{T} matrices.
+#' @export
+#'
+#' @examples
+#' library(EigenR)
+#' A <- toeplitz(c(1, 2, 3))
+#' B <- cbind(c(3, 2, 3), c(1, 1, 1), c(5, 0, -2))
+#' qz <- Eigen_realQZ(A, B)
+#' Q <- qz$Q
+#' Z <- qz$Z
+#' S <- qz$S
+#' T <- qz$T
+#' # check decomposition:
+#' A - Q %*% S %*% Z # should be zero matrix
+#' B - Q %*% T %*% Z # should be zero matrix
+#' # check orthogonality of Q and Z:
+#' tcrossprod(Q) # should be identity matrix
+#' tcrossprod(Z) # should be identity matrix
+Eigen_realQZ <- function(A, B) {
+  stopifnot(isSquareMatrix(A), isSquareMatrix(B))
+  stopifnot(isReal(A), isReal(B))
+  if(nrow(A) != nrow(B)) {
+    stop("The matrices `A` and `B` must have the same size.")
+  }
+  EigenR_realQZ(A, B)
+}
+
+#' Real Schur decomposition
+#' @description Real Schur decomposition of a square matrix.
+#'
+#' @param M real square matrix
+#'
+#' @return A list with the \code{T} and \code{U} matrices.
+#' @export
+#'
+#' @examples
+#' library(EigenR)
+#' M <- cbind(c(3, 2, 3), c(1, 1, 1), c(5, 0, -2))
+#' schur <- Eigen_realSchur(M)
+Eigen_realSchur <- function(M) {
+  stopifnot(isSquareMatrix(M), isReal(M))
+  EigenR_realSchur(M)
+}
+
+
+#' Complex Schur decomposition
+#' @description Complex Schur decomposition of a square matrix.
+#'
+#' @param M real or complex square matrix
+#'
+#' @return A list with the \code{T} and \code{U} matrices.
+#' @export
+#'
+#' @examples
+#' library(EigenR)
+#' M <- cbind(c(3, 2i, 1+3i), c(1, 1i, 1), c(5, 0, -2i))
+#' schur <- Eigen_complexSchur(M)
+Eigen_complexSchur <- function(M) {
+  stopifnot(isSquareMatrix(M), isRealOrComplex(M))
+  schur <- EigenR_complexSchur(Re(M), Im(M))
+  T <- schur[["T"]]
+  U <- schur[["U"]]
+  list(
+    "T" = T[["real"]] + 1i * T[["imag"]],
+    "U" = U[["real"]] + 1i * U[["imag"]]
+  )
+}

inst/RealSchur.cpp

@@ -0,0 +1,41 @@
+#include "EigenR.h"
+
+// [[Rcpp::export]]
+Rcpp::List EigenR_realSchur(const Eigen::MatrixXd& M) {
+  Eigen::RealSchur<Eigen::MatrixXd> schur(M); 
+  const Eigen::MatrixXd U = schur.matrixU();
+  const Eigen::MatrixXd T = schur.matrixT();
+  return Rcpp::List::create(Rcpp::Named("U") = U,
+                            Rcpp::Named("T") = T);
+}
+
+
+// [[Rcpp::export]]
+Rcpp::List EigenR_complexSchur(const Eigen::MatrixXd& Re, 
+	                           const Eigen::MatrixXd& Im) {
+  const Eigen::MatrixXcd M = matricesToMatrixXcd(Re, Im);
+  Eigen::ComplexSchur<MatrixXcd> schur(M.rows());
+  schur.compute(M);
+  const Eigen::MatrixXcd U = schur.matrixU();
+  const Eigen::MatrixXcd T = schur.matrixT();
+  return Rcpp::List::create(Rcpp::Named("U") = cplxMatrixToList(U),
+                            Rcpp::Named("T") = cplxMatrixToList(T));
+}
+
+
+std::complex<double> EigenR_det_cplx(const Eigen::MatrixXd& Re,
+                                     const Eigen::MatrixXd& Im) {
+  const Eigen::MatrixXcd M = matricesToMatrixXcd(Re, Im);
+  return determinant<std::complex<double>>(M);
+}
+
+
+MatrixXd A = MatrixXd::Random(6,6);
+cout << "Here is a random 6x6 matrix, A:" << endl << A << endl << endl;
+
+RealSchur<MatrixXd> schur(A);
+cout << "The orthogonal matrix U is:" << endl << schur.matrixU() << endl;
+cout << "The quasi-triangular matrix T is:" << endl << schur.matrixT() << endl << endl;
+
+MatrixXd U = schur.matrixU();
+MatrixXd T = schur.matrixT();

src/RealQZ.cpp

@@ -0,0 +1,15 @@
+#include "EigenR.h"
+
+// [[Rcpp::export]]
+Rcpp::List EigenR_realQZ(const Eigen::MatrixXd& A, const Eigen::MatrixXd& B) {
+  Eigen::RealQZ<Eigen::MatrixXd> qz(A.rows()); // preallocate space 
+  qz.compute(A, B);  // A = Q S Z,  B = Q T Z
+  const Eigen::MatrixXd Q = qz.matrixQ();
+  const Eigen::MatrixXd Z = qz.matrixZ();
+  const Eigen::MatrixXd S = qz.matrixS();
+  const Eigen::MatrixXd T = qz.matrixT();
+  return Rcpp::List::create(Rcpp::Named("Q") = Q,
+                            Rcpp::Named("Z") = Z,
+                            Rcpp::Named("S") = S,
+                            Rcpp::Named("T") = T);
+}