27 #ifndef TOMOGRAPHER_MATHTOOLS_SOLVECLYAP_H 28 #define TOMOGRAPHER_MATHTOOLS_SOLVECLYAP_H 40 #include <Eigen/Eigenvalues> 45 extern "C" void ztrsyl_(
char *TRANA,
char *TRANB,
int *ISGN,
int *M,
int *N,
46 double *A,
int *lda,
double *B,
int *ldb,
47 double *C,
int *ldc,
double *scale,
int *info);
51 namespace SolveCLyap {
65 virtual const char * what()
const throw() {
70 namespace tomo_internal {
73 "TRANA",
"TRANB",
"ISGN",
"M",
"N",
"A",
"lda",
"B",
"ldb",
"C",
"ldc",
"scale",
"info" 77 template<
bool debug_perform_check>
78 struct solve_check_helper {
79 template<
typename XT,
typename AT,
typename CT,
typename LoggerType>
80 static inline void check(
const XT&,
const AT&,
const CT&, LoggerType&)
87 struct solve_check_helper<true> {
88 template<
typename XT,
typename AT,
typename CT,
typename LoggerType>
89 static inline void check(
const XT& X,
const AT& A,
const CT& C, LoggerType & logger)
93 double n1 = (A.adjoint() * X + X * A - C).norm();
94 double n2 = A.norm() + C.norm();
96 logger.debug(
"LyapCSolve::solve/check",
97 "(A.adjoint() * X + X * A - C).norm()/(A.norm()+C.norm()) == %g/%g == %g ; norm(C)=%g",
98 n1, n2, ok, (
double)C.norm());
100 logger.warning(
"LyapCSolve::solve/check",
"Bad solution quality! rel norm error = %g", (
double)ok);
127 template<
bool debug_perform_check = false,
typename LoggerType>
130 LoggerType & logger,
const double tol = 1e-8)
134 const int d = A.rows();
136 tomographer_assert(X.rows() == X.cols() && X.rows() == d);
137 tomographer_assert(A.rows() == A.cols() && A.rows() == d);
138 tomographer_assert(C.rows() == C.cols() && C.rows() == d);
145 throw SolveError(
"Can't diagonalize matrix A: No Convergence");
152 for (
int k = 0; k < d; ++k) {
153 if (eigvalsA(k) > tol) {
160 MatType D = MatType::Zero(M, M);
161 MatType W = MatType::Zero(d, M);
166 for (
int k = 0; k < d; ++k) {
167 if (eigvalsA(k) < tol) {
171 W.col(counter) = eigU.col(k);
172 D(counter, counter) = eigvalsA(k);
188 Z = W.adjoint() * C * W;
195 int lda = D.outerStride();
197 int ldc = Z.outerStride();
199 ztrsyl_((
char*)
"C", (
char*)
"N", &plusone, &M, &M, (
double*)D.data(), &lda, (
double*)D.data(), &ldb,
200 (
double*)Z.data(), &ldc, &scale, &info);
204 logger.warning(
"SolveCLyap::solve()",
"Warning: A and B have common or very close eigenvalues; " 205 "perturbed values were used to solve the equation");
208 X.noalias() = W * Z * W.adjoint() / scale;
211 tomo_internal::solve_check_helper<debug_perform_check>::check(X, A, C, logger);
216 throw SolveError(
"Argument " + tomo_internal::ztrsyl_argnames[(-info)-1] +
" to ztrsyl_ was invalid.");
ComputationInfo info() const
Base namespace for the Tomographer project.
const RealVectorType & eigenvalues() const
Some C++ utilities, with a tad of C++11 tricks.
Matrix< std::complex< double >, Dynamic, Dynamic > MatrixXcd
internal::plain_col_type< MatrixType, RealScalar >::type RealVectorType
const MatrixType & eigenvectors() const