From e3d329ff4f7ce700e7a92815a07a67778c4e34e6 Mon Sep 17 00:00:00 2001 From: stla Date: Thu, 7 Mar 2024 13:46:06 +0100 Subject: [PATCH] QZ --- .Rprofile | 22 ++++++++++++-- R/EigenR.R | 73 +++++++++++++++++++++++++++++++++++++++++++++ inst/RealSchur.cpp | 41 +++++++++++++++++++++++++ man/Eigen_kernel.Rd | 3 ++ man/Eigen_realQZ.Rd | 33 ++++++++++++++++++++ src/RealQZ.cpp | 15 ++++++++++ 6 files changed, 185 insertions(+), 2 deletions(-) create mode 100644 inst/RealSchur.cpp create mode 100644 man/Eigen_realQZ.Rd create mode 100644 src/RealQZ.cpp diff --git a/.Rprofile b/.Rprofile index 49fb520..38d52f5 100644 --- a/.Rprofile +++ b/.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) + } } diff --git a/R/EigenR.R b/R/EigenR.R index 8c2d638..96ddeef 100644 --- a/R/EigenR.R +++ b/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"]] + ) +} diff --git a/inst/RealSchur.cpp b/inst/RealSchur.cpp new file mode 100644 index 0000000..102901c --- /dev/null +++ b/inst/RealSchur.cpp @@ -0,0 +1,41 @@ +#include "EigenR.h" + +// [[Rcpp::export]] +Rcpp::List EigenR_realSchur(const Eigen::MatrixXd& M) { + Eigen::RealSchur 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 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 EigenR_det_cplx(const Eigen::MatrixXd& Re, + const Eigen::MatrixXd& Im) { + const Eigen::MatrixXcd M = matricesToMatrixXcd(Re, Im); + return determinant>(M); +} + + +MatrixXd A = MatrixXd::Random(6,6); +cout << "Here is a random 6x6 matrix, A:" << endl << A << endl << endl; + +RealSchur 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(); \ No newline at end of file diff --git a/man/Eigen_kernel.Rd b/man/Eigen_kernel.Rd index 823c86a..6c881e1 100644 --- a/man/Eigen_kernel.Rd +++ b/man/Eigen_kernel.Rd @@ -28,3 +28,6 @@ Eigen_kernel(M, method = "LU") # orthonormal basis of the kernel of `M`: Eigen_kernel(M, method = "COD") } +\seealso{ +\code{\link{Eigen_kernelDimension}}. +} diff --git a/man/Eigen_realQZ.Rd b/man/Eigen_realQZ.Rd new file mode 100644 index 0000000..3deb5d7 --- /dev/null +++ b/man/Eigen_realQZ.Rd @@ -0,0 +1,33 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/EigenR.R +\name{Eigen_realQZ} +\alias{Eigen_realQZ} +\title{Real QZ decomposition} +\usage{ +Eigen_realQZ(A, B) +} +\arguments{ +\item{A, B}{real square matrices with the same size} +} +\value{ +A list with the \code{Q}, \code{Z}, \code{S} and \code{T} matrices. +} +\description{ +Real QZ decomposition of a pair of square matrices. +} +\examples{ +library(EigenR) +A <- toeplitz(c(1, 2, 3)) +B <- toeplitz(c(3, 2, 1)) +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 S and T: +tcrossprod(S) # should be identity matrix +tcrossprod(T) # should be identity matrix +} diff --git a/src/RealQZ.cpp b/src/RealQZ.cpp new file mode 100644 index 0000000..fca6826 --- /dev/null +++ b/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 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); +}