Mirrors/eigen
0
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-01-18 17:31:19 +01:00
Code Packages Releases Activity

Merge branch eigen:master into master

Browse Source
This commit is contained in:
Rasmus Munk Larsen
2025-10-13 21:39:39 +00:00
parent 0068623c7e 99f8512985
commit b0ffc9cf7e
29 changed files with 1580 additions and 202 deletions
Download Patch File Download Diff File Expand all files Collapse all files

54
.gitlab/issue_templates/Bug Report.md
View File

@@ -1,42 +1,37 @@
<!--
Please read this!
Thank you for submitting an issue!
Before opening a new issue, make sure to search for keywords in the issues
filtered by "bug::confirmed" or "bug::unconfirmed" and "bugzilla" label:
- https://gitlab.com/libeigen/eigen/-/issues?scope=all&utf8=%E2%9C%93&state=opened&label_name[]=bug%3A%3Aconfirmed
- https://gitlab.com/libeigen/eigen/-/issues?scope=all&utf8=%E2%9C%93&state=opened&label_name[]=bug%3A%3Aunconfirmed
- https://gitlab.com/libeigen/eigen/-/issues?scope=all&utf8=%E2%9C%93&state=opened&label_name[]=bugzilla
and verify the issue you're about to submit isn't a duplicate. -->
Before opening a new issue, please search for keywords in the existing [list of issues](https://gitlab.com/libeigen/eigen/-/issues?state=opened) to verify it isn't a duplicate.
-->
### Summary
<!-- Summarize the bug encountered concisely. -->
### Environment
<!-- Please provide your development environment here -->
<!-- Please provide your development environment. -->
- **Operating System** : Windows/Linux
- **Architecture** : x64/Arm64/PowerPC ...
- **Eigen Version** : 3.3.9
- **Compiler Version** : Gcc7.0
- **Eigen Version** : 5.0.0
- **Compiler Version** : gcc-12.0
- **Compile Flags** : -O3 -march=native
- **Vector Extension** : SSE/AVX/NEON ...
### Minimal Example
<!-- If possible, please create a minimal example here that exhibits the problematic behavior.
You can also link to [godbolt](https://godbolt.org). But please note that you need to click
the "Share" button in the top right-hand corner of the godbolt page where you reproduce the sample
code to get the share link instead of in your browser address bar.
<!--
Please create a minimal reproducing example here that exhibits the problematic behavior.
The example should be complete, in that it can fully build and run. See the [the guidelines on stackoverflow](https://stackoverflow.com/help/minimal-reproducible-example) for how to create a good minimal example.
You can read [the guidelines on stackoverflow](https://stackoverflow.com/help/minimal-reproducible-example)
on how to create a good minimal example. -->
You can also link to [godbolt](https://godbolt.org). Note that you need to click
the "Share" button in the top right-hand corner of the godbolt page to get the share link
instead of the URL in your browser address bar.
-->
```cpp
//show your code here
// Insert your code here.
```
### Steps to reproduce
<!-- Describe how one can reproduce the issue - this is very important. Please use an ordered list. -->
### Steps to reproduce the issue
<!-- Describe the necessary steps to reproduce the issue. -->
1. first step
2. second step
@@ -49,21 +44,16 @@ on how to create a good minimal example. -->
<!-- Describe what you should see instead. -->
### Relevant logs
<!-- Add relevant code snippets or program output within blocks marked by " ``` " -->
<!-- Add relevant build logs or program output within blocks marked by " ``` " -->
<!-- OPTIONAL: remove this section if you are not reporting a compilation warning issue.-->
### Warning Messages
<!-- Show us the warning messages you got! -->
<!-- OPTIONAL: remove this section if you are not reporting a performance issue. -->
### Benchmark scripts and results
### [Optional] Benchmark scripts and results
<!-- Please share any benchmark scripts - either standalone, or using [Google Benchmark](https://github.com/google/benchmark). -->
### Anything else that might help
<!-- It will be better to provide us more information to help narrow down the cause.
<!--
It will be better to provide us more information to help narrow down the cause.
Including but not limited to the following:
- lines of code that might help us diagnose the problem.
- potential ways to address the issue.
- last known working/first broken version (release number or commit hash). -->
- [ ] Have a plan to fix this issue.
- last known working/first broken version (release number or commit hash).
-->

9
.gitlab/issue_templates/Feature Request.md
View File

@@ -1,6 +1,13 @@
<!--
Thank you for submitting a Feature Request!
If you want to run ideas by the maintainers and the Eigen community first,
you can chat about them on the [Eigen Discord server](https://discord.gg/2SkEJGqZjR).
-->
### Describe the feature you would like to be implemented.
### Would such a feature be useful for other users? Why?
### Why Would such a feature be useful for other users?
### Any hints on how to implement the requested feature?

26
.gitlab/merge_request_templates/Merge Request Template.md
View File

@@ -1,26 +0,0 @@
<!--
Thanks for contributing a merge request! Please name and fully describe your MR as you would for a commit message.
If the MR fixes an issue, please include "Fixes #issue" in the commit message and the MR description.
In addition, we recommend that first-time contributors read our [contribution guidelines](https://eigen.tuxfamily.org/index.php?title=Contributing_to_Eigen) and [git page](https://eigen.tuxfamily.org/index.php?title=Git), which will help you submit a more standardized MR.
Before submitting the MR, you also need to complete the following checks:
- Make one PR per feature/bugfix (don't mix multiple changes into one PR). Avoid committing unrelated changes.
- Rebase before committing
- For code changes, run the test suite (at least the tests that are likely affected by the change).
See our [test guidelines](https://eigen.tuxfamily.org/index.php?title=Tests).
- If possible, add a test (both for bug-fixes as well as new features)
- Make sure new features are documented
Note that we are a team of volunteers; we appreciate your patience during the review process.
Again, thanks for contributing! -->
### Reference issue
<!-- You can link to a specific issue using the gitlab syntax #<issue number> -->
### What does this implement/fix?
<!--Please explain your changes.-->
### Additional information
<!--Any additional information you think is important.-->

30
.gitlab/merge_request_templates/Merge Request.md Normal file
View File

@@ -0,0 +1,30 @@
<!--
Thanks for contributing a merge request!
We recommend that first-time contributors read our [contribution guidelines](https://eigen.tuxfamily.org/index.php?title=Contributing_to_Eigen).
Before submitting the MR, please complete the following checks:
- Create one PR per feature or bugfix,
- Run the test suite to verify your changes.
See our [test guidelines](https://eigen.tuxfamily.org/index.php?title=Tests).
- Add tests to cover the bug addressed or any new feature.
- Document new features. If it is a substantial change, add it to the [Changelog](https://gitlab.com/libeigen/eigen/-/blob/master/CHANGELOG.md).
- Leave the following box checked when submitting: `Allow commits from members who can merge to the target branch`.
This allows us to rebase and merge your change.
Note that we are a team of volunteers; we appreciate your patience during the review process.
-->
### Description
<!--Please explain your changes.-->
%{first_multiline_commit}
### Reference issue
<!--
You can link to a specific issue using the gitlab syntax #<issue number>.
If the MR fixes an issue, write "Fixes #<issue number>" to have the issue automatically closed on merge.
-->
### Additional information
<!--Any additional information you think is important.-->

1
Eigen/Eigenvalues
View File

@@ -44,6 +44,7 @@
#include "src/Eigenvalues/ComplexSchur.h"
#include "src/Eigenvalues/ComplexEigenSolver.h"
#include "src/Eigenvalues/RealQZ.h"
#include "src/Eigenvalues/ComplexQZ.h"
#include "src/Eigenvalues/GeneralizedEigenSolver.h"
#include "src/Eigenvalues/MatrixBaseEigenvalues.h"
#ifdef EIGEN_USE_LAPACKE

18
Eigen/src/Core/arch/Default/ConjHelper.h
View File

@@ -17,6 +17,9 @@
EIGEN_STRONG_INLINE PACKET_CPLX pmadd(const PACKET_REAL& x, const PACKET_CPLX& y, const PACKET_CPLX& c) const { \
return padd(c, this->pmul(x, y)); \
} \
EIGEN_STRONG_INLINE PACKET_CPLX pmsub(const PACKET_REAL& x, const PACKET_CPLX& y, const PACKET_CPLX& c) const { \
return psub(this->pmul(x, y), c); \
} \
EIGEN_STRONG_INLINE PACKET_CPLX pmul(const PACKET_REAL& x, const PACKET_CPLX& y) const { \
return PACKET_CPLX(Eigen::internal::pmul<PACKET_REAL>(x, y.v)); \
} \
@@ -27,6 +30,9 @@
EIGEN_STRONG_INLINE PACKET_CPLX pmadd(const PACKET_CPLX& x, const PACKET_REAL& y, const PACKET_CPLX& c) const { \
return padd(c, this->pmul(x, y)); \
} \
EIGEN_STRONG_INLINE PACKET_CPLX pmsub(const PACKET_CPLX& x, const PACKET_REAL& y, const PACKET_CPLX& c) const { \
return psub(this->pmul(x, y), c); \
} \
EIGEN_STRONG_INLINE PACKET_CPLX pmul(const PACKET_CPLX& x, const PACKET_REAL& y) const { \
return PACKET_CPLX(Eigen::internal::pmul<PACKET_REAL>(x.v, y)); \
} \
@@ -76,6 +82,11 @@ struct conj_helper {
return this->pmul(x, y) + c;
}
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE ResultType pmsub(const LhsType& x, const RhsType& y,
const ResultType& c) const {
return this->pmul(x, y) - c;
}
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE ResultType pmul(const LhsType& x, const RhsType& y) const {
return conj_if<ConjLhs>()(x) * conj_if<ConjRhs>()(y);
}
@@ -104,6 +115,10 @@ struct conj_helper<Packet, Packet, ConjLhs, ConjRhs> {
return Eigen::internal::pmadd(conj_if<ConjLhs>().pconj(x), conj_if<ConjRhs>().pconj(y), c);
}
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet pmsub(const Packet& x, const Packet& y, const Packet& c) const {
return Eigen::internal::pmsub(conj_if<ConjLhs>().pconj(x), conj_if<ConjRhs>().pconj(y), c);
}
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet pmul(const Packet& x, const Packet& y) const {
return Eigen::internal::pmul(conj_if<ConjLhs>().pconj(x), conj_if<ConjRhs>().pconj(y));
}
@@ -116,6 +131,9 @@ struct conj_helper<Packet, Packet, true, true> {
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet pmadd(const Packet& x, const Packet& y, const Packet& c) const {
return Eigen::internal::pmadd(pconj(x), pconj(y), c);
}
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet pmsub(const Packet& x, const Packet& y, const Packet& c) const {
return Eigen::internal::pmsub(pconj(x), pconj(y), c);
}
// We save a conjuation by using the identity conj(a)*conj(b) = conj(a*b).
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet pmul(const Packet& x, const Packet& y) const {
return pconj(Eigen::internal::pmul(x, y));

661
Eigen/src/Eigenvalues/ComplexQZ.h Normal file
View File

@@ -0,0 +1,661 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2012 Alexey Korepanov
// Copyright (C) 2025 Ludwig Striet <ludwig.striet@mathematik.uni-freiburg.de>
//
// 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
// https://mozilla.org/MPL/2.0/.
//
// Derived from: Eigen/src/Eigenvalues/RealQZ.h
#include "Eigen/Core"
#include "Eigen/Jacobi"
#ifndef EIGEN_COMPLEX_QZ_H_
#define EIGEN_COMPLEX_QZ_H_
#include "../../SparseQR"
// IWYU pragma: private
#include "./InternalHeaderCheck.h"
/** \eigenvalues_module \ingroup Eigenvalues_Module
*
*
* \class ComplexQZ
*
* \brief Performs a QZ decomposition of a pair of matrices A, B
*
* \tparam MatrixType_ the type input type of the matrix.
*
* Given to complex square matrices A and B, this class computes the QZ decomposition
* \f$ A = Q S Z \f$, \f$ B = Q T Z\f$ where Q and Z are unitary matrices and
* S and T a re upper-triangular matrices. More precisely, Q and Z fulfill
* \f$ Q Q* = Id\f$ and \f$ Z Z* = Id\f$. The generalized Eigenvalues are then
* obtained as ratios of corresponding diagonal entries, lambda(i) = S(i,i) / T(i, i).
*
* The QZ algorithm was introduced in the seminal work "An Algorithm for
* Generalized Matrix Eigenvalue Problems" by Moler & Stewart in 1973. The matrix
* pair S = A, T = B is first transformed to Hessenberg-Triangular form where S is an
* upper Hessenberg matrix and T is an upper Triangular matrix.
*
* This pair is subsequently reduced to the desired form using implicit QZ shifts as
* described in the original paper. The algorithms to find small entries on the
* diagonals and subdiagonals are based on the variants in the implementation
* for Real matrices in the RealQZ class.
*
* \sa class RealQZ
*/
namespace Eigen {
template <typename MatrixType_>
class ComplexQZ {
public:
using MatrixType = MatrixType_;
using Scalar = typename MatrixType_::Scalar;
using RealScalar = typename MatrixType_::RealScalar;
enum {
RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
Options = internal::traits<MatrixType>::Options,
MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
};
using Vec = Matrix<Scalar, Dynamic, 1>;
using Vec2 = Matrix<Scalar, 2, 1>;
using Vec3 = Matrix<Scalar, 3, 1>;
using Row2 = Matrix<Scalar, 1, 2>;
using Mat2 = Matrix<Scalar, 2, 2>;
/** \brief Returns matrix Q in the QZ decomposition.
*
* \returns A const reference to the matrix Q.
*/
const MatrixType& matrixQ() const {
eigen_assert(m_isInitialized && "ComplexQZ is not initialized.");
eigen_assert(m_computeQZ && "The matrices Q and Z have not been computed during the QZ decomposition.");
return m_Q;
}
/** \brief Returns matrix Z in the QZ decomposition.
*
* \returns A const reference to the matrix Z.
*/
const MatrixType& matrixZ() const {
eigen_assert(m_isInitialized && "ComplexQZ is not initialized.");
eigen_assert(m_computeQZ && "The matrices Q and Z have not been computed during the QZ decomposition.");
return m_Z;
}
/** \brief Returns matrix S in the QZ decomposition.
*
* \returns A const reference to the matrix S.
*/
const MatrixType& matrixS() const {
eigen_assert(m_isInitialized && "ComplexQZ is not initialized.");
return m_S;
}
/** \brief Returns matrix S in the QZ decomposition.
*
* \returns A const reference to the matrix S.
*/
const MatrixType& matrixT() const {
eigen_assert(m_isInitialized && "ComplexQZ is not initialized.");
return m_T;
}
/** \brief Constructor
*
* \param[in] n size of the matrices whose QZ decomposition we compute
*
* This constructor is used when we use the compute(...) method later,
* especially when we aim to compute the decomposition of two sparse
* matrices.
*/
ComplexQZ(Index n, bool computeQZ = true, unsigned int maxIters = 400)
: m_n(n),
m_S(n, n),
m_T(n, n),
m_Q(computeQZ ? n : (MatrixType::RowsAtCompileTime == Eigen::Dynamic ? 0 : MatrixType::RowsAtCompileTime),
computeQZ ? n : (MatrixType::ColsAtCompileTime == Eigen::Dynamic ? 0 : MatrixType::ColsAtCompileTime)),
m_Z(computeQZ ? n : (MatrixType::RowsAtCompileTime == Eigen::Dynamic ? 0 : MatrixType::RowsAtCompileTime),
computeQZ ? n : (MatrixType::ColsAtCompileTime == Eigen::Dynamic ? 0 : MatrixType::ColsAtCompileTime)),
m_ws(2 * n),
m_computeQZ(computeQZ),
m_maxIters(maxIters){
};
/** \brief Constructor. computes the QZ decomposition of given matrices
* upon creation
*
* \param[in] A input matrix A
* \param[in] B input matrix B
* \param[in] computeQZ If false, the matrices Q and Z are not computed
*
* This constructor calls the compute() method to compute the QZ decomposition.
* If input matrices are sparse, call the constructor that uses only the
* size as input the computeSparse(...) method.
*/
ComplexQZ(const MatrixType& A, const MatrixType& B, bool computeQZ = true, unsigned int maxIters = 400)
: m_n(A.rows()),
m_maxIters(maxIters),
m_computeQZ(computeQZ),
m_S(A.rows(), A.cols()),
m_T(A.rows(), A.cols()),
m_Q(computeQZ ? m_n : (MatrixType::RowsAtCompileTime == Eigen::Dynamic ? 0 : MatrixType::RowsAtCompileTime),
computeQZ ? m_n : (MatrixType::ColsAtCompileTime == Eigen::Dynamic ? 0 : MatrixType::ColsAtCompileTime)),
m_Z(computeQZ ? m_n : (MatrixType::RowsAtCompileTime == Eigen::Dynamic ? 0 : MatrixType::RowsAtCompileTime),
computeQZ ? m_n : (MatrixType::ColsAtCompileTime == Eigen::Dynamic ? 0 : MatrixType::ColsAtCompileTime)),
m_ws(2 * m_n) {
compute(A, B, computeQZ);
}
/** \brief Compute the QZ decomposition of complex input matrices
*
* \param[in] A Matrix A.
* \param[in] B Matrix B.
* \param[in] computeQZ If false, the matrices Q and Z are not computed.
*/
void compute(const MatrixType& A, const MatrixType& B, bool computeQZ = true);
/** \brief Compute the decomposition of sparse complex input matrices.
* Main difference to the compute(...) method is that it computes a
* SparseQR decomposition of B
*
* \param[in] A Matrix A.
* \param[in] B Matrix B.
* \param[in] computeQZ If false, the matrices Q and Z are not computed.
*/
template <typename SparseMatrixType_>
void computeSparse(const SparseMatrixType_& A, const SparseMatrixType_& B, bool computeQZ = true);
/** \brief Reports whether the last computation was successfull.
*
* \returns \c Success if computation was successfull, \c NoConvergence otherwise.
*/
ComputationInfo info() const { return m_info; };
/** \brief number of performed QZ steps
*/
unsigned int iterations() const {
eigen_assert(m_isInitialized && "ComplexQZ is not initialized.");
return m_global_iter;
};
private:
Index m_n;
const unsigned int m_maxIters;
unsigned int m_global_iter;
bool m_isInitialized;
bool m_computeQZ;
ComputationInfo m_info;
MatrixType m_S, m_T, m_Q, m_Z;
RealScalar m_normOfT, m_normOfS;
Vec m_ws;
// Test if a Scalar is 0 up to a certain tolerance
static bool is_negligible(const Scalar x, const RealScalar tol = NumTraits<RealScalar>::epsilon()) {
return numext::abs(x) <= tol;
}
void do_QZ_step(Index p, Index q);
inline Mat2 computeZk2(const Row2& b);
// This is basically taken from from Eigen3::RealQZ
void hessenbergTriangular(const MatrixType& A, const MatrixType& B);
// This function can be called when m_Q and m_Z are initialized and m_S, m_T
// are in hessenberg-triangular form
void reduceHessenbergTriangular();
// Sparse variant of the above method.
template <typename SparseMatrixType_>
void hessenbergTriangularSparse(const SparseMatrixType_& A, const SparseMatrixType_& B);
void computeNorms();
Index findSmallSubdiagEntry(Index l);
Index findSmallDiagEntry(Index f, Index l);
void push_down_zero_ST(Index k, Index l);
void reduceDiagonal2x2block(Index i);
};
template <typename MatrixType_>
void ComplexQZ<MatrixType_>::compute(const MatrixType& A, const MatrixType& B, bool computeQZ) {
m_computeQZ = computeQZ;
m_n = A.rows();
eigen_assert(m_n == A.cols() && "A is not a square matrix");
eigen_assert(m_n == B.rows() && m_n == B.cols() && "B is not a square matrix or B is not of the same size as A");
m_isInitialized = true;
m_global_iter = 0;
// This will initialize m_Q and m_Z and bring m_S, m_T to hessenberg-triangular form
hessenbergTriangular(A, B);
// We assume that we already have that S is upper-Hessenberg and T is
// upper-triangular. This is what the hessenbergTriangular(...) method does
reduceHessenbergTriangular();
}
// This is basically taken from from Eigen3::RealQZ
template <typename MatrixType_>
void ComplexQZ<MatrixType_>::hessenbergTriangular(const MatrixType& A, const MatrixType& B) {
// Copy A and B, these will be the matrices on which we operate later
m_S = A;
m_T = B;
// Perform QR decomposition of the matrix Q
HouseholderQR<MatrixType> qr(m_T);
m_T = qr.matrixQR();
m_T.template triangularView<StrictlyLower>().setZero();
if (m_computeQZ) m_Q = qr.householderQ();
// overwrite S with Q* x S
m_S.applyOnTheLeft(qr.householderQ().adjoint());
if (m_computeQZ) m_Z = MatrixType::Identity(m_n, m_n);
int steps = 0;
// reduce S to upper Hessenberg with Givens rotations
for (Index j = 0; j <= m_n - 3; j++) {
for (Index i = m_n - 1; i >= j + 2; i--) {
JacobiRotation<Scalar> G;
// delete S(i,j)
if (!numext::is_exactly_zero(m_S.coeff(i, j))) {
G.makeGivens(m_S.coeff(i - 1, j), m_S.coeff(i, j), &m_S.coeffRef(i - 1, j));
m_S.coeffRef(i, j) = Scalar(0);
m_T.rightCols(m_n - i + 1).applyOnTheLeft(i - 1, i, G.adjoint());
m_S.rightCols(m_n - j - 1).applyOnTheLeft(i - 1, i, G.adjoint());
// This is what we want to achieve
if (!is_negligible(m_S(i, j)))
m_info = ComputationInfo::NumericalIssue;
else
m_S(i, j) = Scalar(0);
// update Q
if (m_computeQZ) m_Q.applyOnTheRight(i - 1, i, G);
}
if (!numext::is_exactly_zero(m_T.coeff(i, i - 1))) {
// Compute rotation and update matrix T
G.makeGivens(m_T.coeff(i, i), m_T.coeff(i, i - 1), &m_T.coeffRef(i, i));
m_T.topRows(i).applyOnTheRight(i - 1, i, G.adjoint());
m_T.coeffRef(i, i - 1) = Scalar(0);
// Update matrix S
m_S.applyOnTheRight(i - 1, i, G.adjoint());
// update Z
if (m_computeQZ) m_Z.applyOnTheLeft(i - 1, i, G);
}
steps++;
}
}
}
template <typename MatrixType>
template <typename SparseMatrixType_>
void ComplexQZ<MatrixType>::hessenbergTriangularSparse(const SparseMatrixType_& A, const SparseMatrixType_& B) {
m_S = A.toDense();
SparseQR<SparseMatrix<Scalar, ColMajor>, NaturalOrdering<Index>> sparseQR;
eigen_assert(B.isCompressed() &&
"SparseQR requires a sparse matrix in compressed mode."
"Call .makeCompressed() before passing it to SparseQR");
// Computing QR decomposition of T...
sparseQR.setPivotThreshold(RealScalar(0)); // This prevends algorithm from doing pivoting
sparseQR.compute(B);
// perform QR decomposition of T, overwrite T with R, save Q
// HouseholderQR<Mat> qrT(m_T);
m_T = sparseQR.matrixR();
m_T.template triangularView<StrictlyLower>().setZero();
if (m_computeQZ) m_Q = sparseQR.matrixQ();
// overwrite S with Q* S
m_S = sparseQR.matrixQ().adjoint() * m_S;
if (m_computeQZ) m_Z = MatrixType::Identity(m_n, m_n);
unsigned int steps = 0;
// reduce S to upper Hessenberg with Givens rotations
for (Index j = 0; j <= m_n - 3; j++) {
for (Index i = m_n - 1; i >= j + 2; i--) {
JacobiRotation<Scalar> G;
// kill S(i,j)
// if(!numext::is_exactly_zero(_S.coeff(i, j)))
if (m_S.coeff(i, j) != Scalar(0)) {
// This is the adapted code
G.makeGivens(m_S.coeff(i - 1, j), m_S.coeff(i, j), &m_S.coeffRef(i - 1, j));
m_S.coeffRef(i, j) = Scalar(0);
m_T.rightCols(m_n - i + 1).applyOnTheLeft(i - 1, i, G.adjoint());
m_S.rightCols(m_n - j - 1).applyOnTheLeft(i - 1, i, G.adjoint());
// This is what we want to achieve
if (!is_negligible(m_S(i, j))) {
m_info = ComputationInfo::NumericalIssue;
}
m_S(i, j) = Scalar(0);
// update Q
if (m_computeQZ) m_Q.applyOnTheRight(i - 1, i, G);
}
if (!numext::is_exactly_zero(m_T.coeff(i, i - 1))) {
// Compute rotation and update matrix T
G.makeGivens(m_T.coeff(i, i), m_T.coeff(i, i - 1), &m_T.coeffRef(i, i));
m_T.topRows(i).applyOnTheRight(i - 1, i, G.adjoint());
m_T.coeffRef(i, i - 1) = Scalar(0);
// Update matrix S
m_S.applyOnTheRight(i - 1, i, G.adjoint());
// update Z
if (m_computeQZ) m_Z.applyOnTheLeft(i - 1, i, G);
}
steps++;
}
}
}
template <typename MatrixType>
template <typename SparseMatrixType_>
void ComplexQZ<MatrixType>::computeSparse(const SparseMatrixType_& A, const SparseMatrixType_& B, bool computeQZ) {
m_computeQZ = computeQZ;
m_n = A.rows();
eigen_assert(m_n == A.cols() && "A is not a square matrix");
eigen_assert(m_n == B.rows() && m_n == B.cols() && "B is not a square matrix or B is not of the same size as A");
m_isInitialized = true;
m_global_iter = 0;
hessenbergTriangularSparse(A, B);
// We assume that we already have that A is upper-Hessenberg and B is
// upper-triangular. This is what the hessenbergTriangular(...) method does
reduceHessenbergTriangular();
}
template <typename MatrixType_>
void ComplexQZ<MatrixType_>::reduceHessenbergTriangular() {
Index l = m_n - 1, f;
unsigned int local_iter = 0;
computeNorms();
while (l > 0 && local_iter < m_maxIters) {
f = findSmallSubdiagEntry(l);
// Subdiag entry is small -> can be safely set to 0
if (f > 0) {
m_S.coeffRef(f, f - 1) = Scalar(0);
}
if (f == l) { // One root found
l--;
local_iter = 0;
} else if (f == l - 1) { // Two roots found
// We found an undesired non-zero at (f+1,f) in S and eliminate it immediately
reduceDiagonal2x2block(f);
l -= 2;
local_iter = 0;
} else {
Index z = findSmallDiagEntry(f, l);
if (z >= f) {
push_down_zero_ST(z, l);
} else {
do_QZ_step(f, m_n - l - 1);
local_iter++;
m_global_iter++;
}
}
}
m_info = (local_iter < m_maxIters) ? Success : NoConvergence;
}
template <typename MatrixType_>
inline typename ComplexQZ<MatrixType_>::Mat2 ComplexQZ<MatrixType_>::computeZk2(const Row2& b) {
Mat2 S;
S << Scalar(0), Scalar(1), Scalar(1), Scalar(0);
Vec2 bprime = S * b.adjoint();
JacobiRotation<Scalar> J;
J.makeGivens(bprime(0), bprime(1));
Mat2 Z = S;
Z.applyOnTheLeft(0, 1, J);
Z = S * Z;
return Z;
}
template <typename MatrixType_>
void ComplexQZ<MatrixType_>::do_QZ_step(Index p, Index q) {
// This is certainly not the most efficient way of doing this,
// but a readable one.
const auto a = [p, this](Index i, Index j) { return m_S(p + i - 1, p + j - 1); };
const auto b = [p, this](Index i, Index j) { return m_T(p + i - 1, p + j - 1); };
const Index m = m_n - p - q; // Size of the inner block
Scalar x, y, z;
// We could introduce doing exceptional shifts from time to time.
Scalar W1 = a(m - 1, m - 1) / b(m - 1, m - 1) - a(1, 1) / b(1, 1), W2 = a(m, m) / b(m, m) - a(1, 1) / b(1, 1),
W3 = a(m, m - 1) / b(m - 1, m - 1);
x = (W1 * W2 - a(m - 1, m) / b(m, m) * W3 + W3 * b(m - 1, m) / b(m, m) * a(1, 1) / b(1, 1)) * b(1, 1) / a(2, 1) +
a(1, 2) / b(2, 2) - a(1, 1) / b(1, 1) * b(1, 2) / b(2, 2);
y = (a(2, 2) / b(2, 2) - a(1, 1) / b(1, 1)) - a(2, 1) / b(1, 1) * b(1, 2) / b(2, 2) - W1 - W2 +
W3 * (b(m - 1, m) / b(m, m));
z = a(3, 2) / b(2, 2);
Vec3 X;
const PermutationMatrix<3, 3, int> S3(Vector3i(2, 0, 1));
for (Index k = p; k < p + m - 2; k++) {
X << x, y, z;
Vec2 ess;
Scalar tau;
RealScalar beta;
X.makeHouseholder(ess, tau, beta);
// The permutations are needed because the makeHouseHolder-method computes
// the householder transformation in a way that the vector is reflected to
// (1 0 ... 0) instead of (0 ... 0 1)
m_S.template middleRows<3>(k)
.rightCols((std::min)(m_n, m_n - k + 1))
.applyHouseholderOnTheLeft(ess, tau, m_ws.data());
m_T.template middleRows<3>(k).rightCols(m_n - k).applyHouseholderOnTheLeft(ess, tau, m_ws.data());
if (m_computeQZ) m_Q.template middleCols<3>(k).applyHouseholderOnTheRight(ess, std::conj(tau), m_ws.data());
// Compute Matrix Zk1 s.t. (b(k+2,k) ... b(k+2, k+2)) Zk1 = (0,0,*)
Vec3 bprime = (m_T.template block<1, 3>(k + 2, k) * S3).adjoint();
bprime.makeHouseholder(ess, tau, beta);
m_S.template middleCols<3>(k).topRows((std::min)(k + 4, m_n)).applyOnTheRight(S3);
m_S.template middleCols<3>(k)
.topRows((std::min)(k + 4, m_n))
.applyHouseholderOnTheRight(ess, std::conj(tau), m_ws.data());
m_S.template middleCols<3>(k).topRows((std::min)(k + 4, m_n)).applyOnTheRight(S3.transpose());
m_T.template middleCols<3>(k).topRows((std::min)(k + 3, m_n)).applyOnTheRight(S3);
m_T.template middleCols<3>(k)
.topRows((std::min)(k + 3, m_n))
.applyHouseholderOnTheRight(ess, std::conj(tau), m_ws.data());
m_T.template middleCols<3>(k).topRows((std::min)(k + 3, m_n)).applyOnTheRight(S3.transpose());
if (m_computeQZ) {
m_Z.template middleRows<3>(k).applyOnTheLeft(S3.transpose());
m_Z.template middleRows<3>(k).applyHouseholderOnTheLeft(ess, tau, m_ws.data());
m_Z.template middleRows<3>(k).applyOnTheLeft(S3);
}
Mat2 Zk2 = computeZk2(m_T.template block<1, 2>(k + 1, k));
m_S.template middleCols<2>(k).topRows((std::min)(k + 4, m_n)).applyOnTheRight(Zk2);
m_T.template middleCols<2>(k).topRows((std::min)(k + 3, m_n)).applyOnTheRight(Zk2);
if (m_computeQZ) m_Z.template middleRows<2>(k).applyOnTheLeft(Zk2.adjoint());
x = m_S(k + 1, k);
y = m_S(k + 2, k);
if (k < p + m - 3) {
z = m_S(k + 3, k);
}
};
// Find a Householdermartirx Qn1 s.t. Qn1 (x y)^T = (* 0)
JacobiRotation<Scalar> J;
J.makeGivens(x, y);
m_S.template middleRows<2>(p + m - 2).applyOnTheLeft(0, 1, J.adjoint());
m_T.template middleRows<2>(p + m - 2).applyOnTheLeft(0, 1, J.adjoint());
if (m_computeQZ) m_Q.template middleCols<2>(p + m - 2).applyOnTheRight(0, 1, J);
// Find a Householdermatrix Zn1 s.t. (b(n,n-1) b(n,n)) * Zn1 = (0 *)
Mat2 Zn1 = computeZk2(m_T.template block<1, 2>(p + m - 1, p + m - 2));
m_S.template middleCols<2>(p + m - 2).applyOnTheRight(Zn1);
m_T.template middleCols<2>(p + m - 2).applyOnTheRight(Zn1);
if (m_computeQZ) m_Z.template middleRows<2>(p + m - 2).applyOnTheLeft(Zn1.adjoint());
}
/** \internal we found an undesired non-zero at (i+1,i) on the subdiagonal of S and reduce the block */
template <typename MatrixType_>
void ComplexQZ<MatrixType_>::reduceDiagonal2x2block(Index i) {
// We have found a non-zero on the subdiagonal and want to eliminate it
Mat2 Si = m_S.template block<2, 2>(i, i), Ti = m_T.template block<2, 2>(i, i);
if (is_negligible(Ti(0, 0)) && !is_negligible(Ti(1, 1))) {
Eigen::JacobiRotation<Scalar> G;
G.makeGivens(m_S(i, i), m_S(i + 1, i));
m_S.applyOnTheLeft(i, i + 1, G.adjoint());
m_T.applyOnTheLeft(i, i + 1, G.adjoint());
if (m_computeQZ) m_Q.applyOnTheRight(i, i + 1, G);
} else if (!is_negligible(Ti(0, 0)) && is_negligible(Ti(1, 1))) {
Eigen::JacobiRotation<Scalar> G;
G.makeGivens(m_S(i + 1, i + 1), m_S(i + 1, i));
m_S.applyOnTheRight(i, i + 1, G.adjoint());
m_T.applyOnTheRight(i, i + 1, G.adjoint());
if (m_computeQZ) m_Z.applyOnTheLeft(i, i + 1, G);
} else if (!is_negligible(Ti(0, 0)) && !is_negligible((Ti(1, 1)))) {
Scalar mu = Si(0, 0) / Ti(0, 0);
Scalar a12_bar = Si(0, 1) - mu * Ti(0, 1);
Scalar a22_bar = Si(1, 1) - mu * Ti(1, 1);
Scalar p = Scalar(0.5) * (a22_bar / Ti(1, 1) - Ti(0, 1) * Si(1, 0) / (Ti(0, 0) * Ti(1, 1)));
RealScalar sgn_p = p.real() >= RealScalar(0) ? RealScalar(1) : RealScalar(-1);
Scalar q = Si(1, 0) * a12_bar / (Ti(0, 0) * Ti(1, 1));
Scalar r = p * p + q;
Scalar lambda = mu + p + sgn_p * numext::sqrt(r);
Mat2 E = Si - lambda * Ti;
Index l;
E.rowwise().norm().maxCoeff(&l);
JacobiRotation<Scalar> G;
G.makeGivens(E(l, 1), E(l, 0));
m_S.applyOnTheRight(i, i + 1, G.adjoint());
m_T.applyOnTheRight(i, i + 1, G.adjoint());
if (m_computeQZ) m_Z.applyOnTheLeft(i, i + 1, G);
Mat2 tildeSi = m_S.template block<2, 2>(i, i), tildeTi = m_T.template block<2, 2>(i, i);
Mat2 C = tildeSi.norm() < (lambda * tildeTi).norm() ? tildeSi : lambda * tildeTi;
G.makeGivens(C(0, 0), C(1, 0));
m_S.applyOnTheLeft(i, i + 1, G.adjoint());
m_T.applyOnTheLeft(i, i + 1, G.adjoint());
if (m_computeQZ) m_Q.applyOnTheRight(i, i + 1, G);
}
if (!is_negligible(m_S(i + 1, i), m_normOfS * NumTraits<RealScalar>::epsilon())) {
m_info = ComputationInfo::NumericalIssue;
} else {
m_S(i + 1, i) = Scalar(0);
}
}
/** \internal We found a zero at T(k,k) and want to "push it down" to T(l,l) */
template <typename MatrixType_>
void ComplexQZ<MatrixType_>::push_down_zero_ST(Index k, Index l) {
// Test Preconditions
JacobiRotation<Scalar> J;
for (Index j = k + 1; j <= l; j++) {
// Create a 0 at _T(j, j)
J.makeGivens(m_T(j - 1, j), m_T(j, j), &m_T.coeffRef(j - 1, j));
m_T.rightCols(m_n - j - 1).applyOnTheLeft(j - 1, j, J.adjoint());
m_T.coeffRef(j, j) = Scalar(0);
m_S.applyOnTheLeft(j - 1, j, J.adjoint());
if (m_computeQZ) m_Q.applyOnTheRight(j - 1, j, J);
// Delete the non-desired non-zero at _S(j, j-2)
if (j > 1) {
J.makeGivens(std::conj(m_S(j, j - 1)), std::conj(m_S(j, j - 2)));
m_S.applyOnTheRight(j - 1, j - 2, J);
m_S(j, j - 2) = Scalar(0);
m_T.applyOnTheRight(j - 1, j - 2, J);
if (m_computeQZ) m_Z.applyOnTheLeft(j - 1, j - 2, J.adjoint());
}
}
// Assume we have the desired structure now, up to the non-zero entry at
// _S(l, l-1) which we will delete through a last right-jacobi-rotation
J.makeGivens(std::conj(m_S(l, l)), std::conj(m_S(l, l - 1)));
m_S.topRows(l + 1).applyOnTheRight(l, l - 1, J);
if (!is_negligible(m_S(l, l - 1), m_normOfS * NumTraits<Scalar>::epsilon())) {
m_info = ComputationInfo::NumericalIssue;
} else {
m_S(l, l - 1) = Scalar(0);
}
m_T.topRows(l + 1).applyOnTheRight(l, l - 1, J);
if (m_computeQZ) m_Z.applyOnTheLeft(l, l - 1, J.adjoint());
// Ensure postconditions
if (!is_negligible(m_T(l, l)) || !is_negligible(m_S(l, l - 1))) {
m_info = ComputationInfo::NumericalIssue;
} else {
m_T(l, l) = Scalar(0);
m_S(l, l - 1) = Scalar(0);
}
};
/** \internal Computes vector L1 norms of S and T when in Hessenberg-Triangular form already */
template <typename MatrixType_>
void ComplexQZ<MatrixType_>::computeNorms() {
const Index size = m_S.cols();
m_normOfS = RealScalar(0);
m_normOfT = RealScalar(0);
for (Index j = 0; j < size; ++j) {
m_normOfS += m_S.col(j).segment(0, (std::min)(size, j + 2)).cwiseAbs().sum();
m_normOfT += m_T.row(j).segment(j, size - j).cwiseAbs().sum();
}
};
/** \internal Look for single small sub-diagonal element S(res, res-1) and return res (or 0). Copied from Eigen3 RealQZ
* implementation */
template <typename MatrixType_>
inline Index ComplexQZ<MatrixType_>::findSmallSubdiagEntry(Index iu) {
Index res = iu;
while (res > 0) {
RealScalar s = numext::abs(m_S.coeff(res - 1, res - 1)) + numext::abs(m_S.coeff(res, res));
if (s == Scalar(0)) s = m_normOfS;
if (numext::abs(m_S.coeff(res, res - 1)) < NumTraits<RealScalar>::epsilon() * s) break;
res--;
}
return res;
}
//
/** \internal Look for single small diagonal element T(res, res) for res between f and l, and return res (or f-1).
* Copied from Eigen3 RealQZ implementation. */
template <typename MatrixType_>
inline Index ComplexQZ<MatrixType_>::findSmallDiagEntry(Index f, Index l) {
Index res = l;
while (res >= f) {
if (numext::abs(m_T.coeff(res, res)) <= NumTraits<RealScalar>::epsilon() * m_normOfT) break;
res--;
}
return res;
}
} // namespace Eigen
#endif // _COMPLEX_QZ_H_

20
Eigen/src/Jacobi/Jacobi.h
View File

@@ -335,8 +335,8 @@ struct apply_rotation_in_the_plane_selector<Scalar, OtherScalar, SizeAtCompileTi
for (Index i = alignedStart; i < alignedEnd; i += PacketSize) {
Packet xi = pload<Packet>(px);
Packet yi = pload<Packet>(py);
pstore(px, padd(pm.pmul(pc, xi), pcj.pmul(ps, yi)));
pstore(py, psub(pcj.pmul(pc, yi), pm.pmul(ps, xi)));
pstore(px, pm.pmadd(pc, xi, pcj.pmul(ps, yi)));
pstore(py, pcj.pmsub(pc, yi, pm.pmul(ps, xi)));
px += PacketSize;
py += PacketSize;
}
@@ -347,18 +347,18 @@ struct apply_rotation_in_the_plane_selector<Scalar, OtherScalar, SizeAtCompileTi
Packet xi1 = ploadu<Packet>(px + PacketSize);
Packet yi = pload<Packet>(py);
Packet yi1 = pload<Packet>(py + PacketSize);
pstoreu(px, padd(pm.pmul(pc, xi), pcj.pmul(ps, yi)));
pstoreu(px + PacketSize, padd(pm.pmul(pc, xi1), pcj.pmul(ps, yi1)));
pstore(py, psub(pcj.pmul(pc, yi), pm.pmul(ps, xi)));
pstore(py + PacketSize, psub(pcj.pmul(pc, yi1), pm.pmul(ps, xi1)));
pstoreu(px, pm.pmadd(pc, xi, pcj.pmul(ps, yi)));
pstoreu(px + PacketSize, pm.pmadd(pc, xi1, pcj.pmul(ps, yi1)));
pstore(py, pcj.pmsub(pc, yi, pm.pmul(ps, xi)));
pstore(py + PacketSize, pcj.pmsub(pc, yi1, pm.pmul(ps, xi1)));
px += Peeling * PacketSize;
py += Peeling * PacketSize;
}
if (alignedEnd != peelingEnd) {
Packet xi = ploadu<Packet>(x + peelingEnd);
Packet yi = pload<Packet>(y + peelingEnd);
pstoreu(x + peelingEnd, padd(pm.pmul(pc, xi), pcj.pmul(ps, yi)));
pstore(y + peelingEnd, psub(pcj.pmul(pc, yi), pm.pmul(ps, xi)));
pstoreu(x + peelingEnd, pm.pmadd(pc, xi, pcj.pmul(ps, yi)));
pstore(y + peelingEnd, pcj.pmsub(pc, yi, pm.pmul(ps, xi)));
}
}
@@ -381,8 +381,8 @@ struct apply_rotation_in_the_plane_selector<Scalar, OtherScalar, SizeAtCompileTi
for (Index i = 0; i < size; i += PacketSize) {
Packet xi = pload<Packet>(px);
Packet yi = pload<Packet>(py);
pstore(px, padd(pm.pmul(pc, xi), pcj.pmul(ps, yi)));
pstore(py, psub(pcj.pmul(pc, yi), pm.pmul(ps, xi)));
pstore(px, pm.pmadd(pc, xi, pcj.pmul(ps, yi)));
pstore(py, pcj.pmsub(pc, yi, pm.pmul(ps, xi)));
px += PacketSize;
py += PacketSize;
}

2
Eigen/src/KLUSupport/KLUSupport.h
View File

@@ -182,7 +182,7 @@ class KLU : public SparseSolverBase<KLU<MatrixType_> > {
/** Performs a numeric decomposition of \a matrix
*
* The given matrix must has the same sparsity than the matrix on which the pattern anylysis has been performed.
* The given matrix must have the same sparsity than the matrix on which the pattern anylysis has been performed.
*
* \sa analyzePattern(), compute()
*/

3
Eigen/src/PardisoSupport/PardisoSupport.h
View File

@@ -157,7 +157,8 @@ class PardisoImpl : public SparseSolverBase<Derived> {
/** Performs a numeric decomposition of \a matrix
*
* The given matrix must has the same sparsity than the matrix on which the symbolic decomposition has been performed.
* The given matrix must have the same sparsity than the matrix on which the symbolic decomposition has been
* performed.
*
* \sa analyzePattern()
*/

15
Eigen/src/SparseCholesky/SimplicialCholesky.h
View File

@@ -416,7 +416,8 @@ class SimplicialLLT : public SimplicialCholeskyBase<SimplicialLLT<MatrixType_, U
/** Performs a numeric decomposition of \a matrix
*
* The given matrix must has the same sparsity than the matrix on which the symbolic decomposition has been performed.
* The given matrix must have the same sparsity than the matrix on which the symbolic decomposition has been
* performed.
*
* \sa analyzePattern()
*/
@@ -504,7 +505,8 @@ class SimplicialLDLT : public SimplicialCholeskyBase<SimplicialLDLT<MatrixType_,
/** Performs a numeric decomposition of \a matrix
*
* The given matrix must has the same sparsity than the matrix on which the symbolic decomposition has been performed.
* The given matrix must have the same sparsity than the matrix on which the symbolic decomposition has been
* performed.
*
* \sa analyzePattern()
*/
@@ -585,7 +587,8 @@ class SimplicialNonHermitianLLT
/** Performs a numeric decomposition of \a matrix
*
* The given matrix must has the same sparsity than the matrix on which the symbolic decomposition has been performed.
* The given matrix must have the same sparsity than the matrix on which the symbolic decomposition has been
* performed.
*
* \sa analyzePattern()
*/
@@ -674,7 +677,8 @@ class SimplicialNonHermitianLDLT
/** Performs a numeric decomposition of \a matrix
*
* The given matrix must has the same sparsity than the matrix on which the symbolic decomposition has been performed.
* The given matrix must have the same sparsity than the matrix on which the symbolic decomposition has been
* performed.
*
* \sa analyzePattern()
*/
@@ -757,7 +761,8 @@ class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<Matr
/** Performs a numeric decomposition of \a matrix
*
* The given matrix must has the same sparsity than the matrix on which the symbolic decomposition has been performed.
* The given matrix must have the same sparsity than the matrix on which the symbolic decomposition has been
* performed.
*
* \sa analyzePattern()
*/

6
Eigen/src/SuperLUSupport/SuperLUSupport.h
View File

@@ -487,7 +487,8 @@ class SuperLU : public SuperLUBase<MatrixType_, SuperLU<MatrixType_> > {
/** Performs a numeric decomposition of \a matrix
*
* The given matrix must has the same sparsity than the matrix on which the symbolic decomposition has been performed.
* The given matrix must have the same sparsity than the matrix on which the symbolic decomposition has been
* performed.
*
* \sa analyzePattern()
*/
@@ -791,7 +792,8 @@ class SuperILU : public SuperLUBase<MatrixType_, SuperILU<MatrixType_> > {
/** Performs a numeric decomposition of \a matrix
*
* The given matrix must has the same sparsity than the matrix on which the symbolic decomposition has been performed.
* The given matrix must have the same sparsity than the matrix on which the symbolic decomposition has been
* performed.
*
* \sa analyzePattern()
*/

2
Eigen/src/UmfPackSupport/UmfPackSupport.h
View File

@@ -425,7 +425,7 @@ class UmfPackLU : public SparseSolverBase<UmfPackLU<MatrixType_> > {
/** Performs a numeric decomposition of \a matrix
*
* The given matrix must has the same sparsity than the matrix on which the pattern anylysis has been performed.
* The given matrix must have the same sparsity than the matrix on which the pattern anylysis has been performed.
*
* \sa analyzePattern(), compute()
*/

5
blas/CMakeLists.txt
View File

@@ -19,9 +19,10 @@ add_library(eigen_blas_static STATIC ${EigenBlas_SRCS})
list(APPEND EIGEN_BLAS_TARGETS eigen_blas_static)
if (EIGEN_BUILD_SHARED_LIBS)
add_library(eigen_blas SHARED ${EigenBlas_SRCS})
add_library(eigen_blas SHARED ${EigenBlas_SRCS} "eigen_blas.def")
target_compile_definitions(eigen_blas PUBLIC "EIGEN_BLAS_BUILD_DLL")
set_target_properties(eigen_blas PROPERTIES CXX_VISIBILITY_PRESET hidden)
list(APPEND EIGEN_BLAS_TARGETS eigen_blas)
target_compile_definitions(eigen_blas PUBLIC "EIGEN_BUILD_DLL")
endif()
foreach(target IN LISTS EIGEN_BLAS_TARGETS)

12
blas/blas.h
View File

@@ -2,13 +2,15 @@
#define BLAS_H
#if defined(_WIN32)
#if defined(EIGEN_BUILD_DLL)
#if defined(EIGEN_BLAS_BUILD_DLL)
#define EIGEN_BLAS_API __declspec(dllexport)
#elif defined(EIGEN_LINK_DLL)
#elif defined(EIGEN_BLAS_LINK_DLL)
#define EIGEN_BLAS_API __declspec(dllimport)
#else
#define EIGEN_BLAS_API
#endif
#elif ((defined(__GNUC__) && __GNUC__ >= 4) || defined(__clang__)) && defined(EIGEN_BLAS_BUILD_DLL)
#define EIGEN_BLAS_API __attribute__((visibility("default")))
#else
#define EIGEN_BLAS_API
#endif
@@ -27,6 +29,7 @@ typedef long BLASLONG;
typedef unsigned long BLASULONG;
#endif
EIGEN_BLAS_API int BLASFUNC(lsame)(const char *, const char *);
EIGEN_BLAS_API void BLASFUNC(xerbla)(const char *, int *info);
EIGEN_BLAS_API float BLASFUNC(sdot)(int *, float *, int *, float *, int *);
@@ -36,6 +39,11 @@ EIGEN_BLAS_API double BLASFUNC(dsdot)(int *, float *, int *, float *, int *);
EIGEN_BLAS_API double BLASFUNC(ddot)(int *, double *, int *, double *, int *);
EIGEN_BLAS_API double BLASFUNC(qdot)(int *, double *, int *, double *, int *);
EIGEN_BLAS_API void BLASFUNC(cdotu)(int *, float *, int *, float *, int *);
EIGEN_BLAS_API void BLASFUNC(cdotc)(int *, float *, int *, float *, int *);
EIGEN_BLAS_API void BLASFUNC(zdotu)(int *, double *, int *, double *, int *);
EIGEN_BLAS_API void BLASFUNC(zdotc)(int *, double *, int *, double *, int *);
EIGEN_BLAS_API void BLASFUNC(cdotuw)(int *, float *, int *, float *, int *, float *);
EIGEN_BLAS_API void BLASFUNC(cdotcw)(int *, float *, int *, float *, int *, float *);
EIGEN_BLAS_API void BLASFUNC(zdotuw)(int *, double *, int *, double *, int *, double *);

222
blas/eigen_blas.def Normal file
View File

@@ -0,0 +1,222 @@
; Definition file for eigen_blas.dll.
LIBRARY eigen_blas
EXPORTS
; Utilities
lsame_
xerbla_
; Level 1
saxpy_
daxpy_
caxpy_
zaxpy_
; caxpyc_
; zaxpyc_
scopy_
dcopy_
ccopy_
zcopy_
sdot_
sdsdot_
dsdot_
ddot_
cdotc_
zdotc_
cdotu_
zdotu_
cdotcw_
zdotcw_
cdotuw_
zdotuw_
snrm2_
dnrm2_
scnrm2_
dznrm2_
srot_
drot_
csrot_
zdrot_
srotg_
drotg_
crotg_
zrotg_
srotm_
drotm_
srotmg_
drotmg_
sscal_
dscal_
cscal_
zscal_
csscal_
zdscal_
sswap_
dswap_
cswap_
zswap_
sasum_
scasum_
dasum_
dzasum_
; ismax_
; idmax_
; icmax_
; izmax_
isamax_
idamax_
icamax_
izamax_
isamin_
idamin_
icamin_
izamin_
; ismin_
; idmin_
; icmin_
; izmin_
; samax_
; damax_
; scamax_
; dzamax_
; samin_
; damin_
; scamin_
; dzamin_
; smax_
; dmax_
; cmax_
; zmax_
; smin_
; dmin_
; cmin_
; zmin_
; Level 2
sgemv_
dgemv_
cgemv_
zgemv_
sger_
dger_
cgerc_
zgerc_
cgeru_
zgeru_
ssymv_
dsymv_
ssyr_
dsyr_
ssyr2_
dsyr2_
; csyr2_
; zsyr2_
strmv_
dtrmv_
ctrmv_
ztrmv_
strsv_
dtrsv_
ctrsv_
ztrsv_
stpsv_
dtpsv_
ctpsv_
ztpsv_
stpmv_
dtpmv_
ctpmv_
ztpmv_
stbmv_
dtbmv_
ctbmv_
ztbmv_
stbsv_
dtbsv_
ctbsv_
ztbsv_
sspmv_
dspmv_
sspr_
dspr_
; cspr_
; zspr_
sspr2_
dspr2_
; cspr2_
; zspr2_
cher_
zher_
chpr_
zhpr_
cher2_
zher2_
chpr2_
zhpr2_
chemv_
zhemv_
chpmv_
zhpmv_
; snorm_
; dnorm_
; cnorm_
; znorm_
sgbmv_
dgbmv_
cgbmv_
zgbmv_
ssbmv_
dsbmv_
; csbmv_
; zsbmv_
chbmv_
zhbmv_
; Level 3 BLAS
sgemm_
dgemm_
cgemm_
zgemm_
; cgemm3m_
; zgemm3m_
; sge2mm_
; dge2mm_
; cge2mm_
; zge2mm_
ssymm_
dsymm_
csymm_
zsymm_
; csymm3m_
; zsymm3m_
ssyrk_
dsyrk_
csyrk_
zsyrk_
ssyr2k_
dsyr2k_
csyr2k_
zsyr2k_
strmm_
dtrmm_
ctrmm_
ztrmm_
strsm_
dtrsm_
ctrsm_
ztrsm_
chemm_
zhemm_
; chemm3m_
; zhemm3m_
cherk_
zherk_
cher2k_
zher2k_
; cher2m_
; zher2m_
sgemmtr_
dgemmtr_
cgemmtr_
zgemmtr_

2
blas/xerbla.cpp
View File

@@ -1,6 +1,8 @@
#include <stdio.h>
#include "blas.h"
#if (defined __GNUC__) && (!defined __MINGW32__) && (!defined __CYGWIN__)
#define EIGEN_WEAK_LINKING __attribute__((weak))
#else

15
lapack/CMakeLists.txt
View File

@@ -22,9 +22,11 @@ add_custom_target(lapack)
include_directories(../blas)
set(EigenLapack_SRCS
dsecnd_INT_CPU_TIME.cpp second_INT_CPU_TIME.cpp single.cpp double.cpp complex_single.cpp complex_double.cpp ../blas/xerbla.cpp
dsecnd_INT_CPU_TIME.cpp second_INT_CPU_TIME.cpp single.cpp double.cpp complex_single.cpp complex_double.cpp
)
set(EIGEN_LAPACK_DEF "eigen_lapack_cpp.def")
if(EIGEN_Fortran_COMPILER_WORKS)
set(EigenLapack_SRCS ${EigenLapack_SRCS}
@@ -40,6 +42,8 @@ set(EigenLapack_SRCS ${EigenLapack_SRCS}
slamch.f dlamch.f
)
set(EIGEN_LAPACK_DEF "eigen_lapack.def")
option(EIGEN_ENABLE_LAPACK_TESTS OFF "Enable the Lapack unit tests")
if(EIGEN_ENABLE_LAPACK_TESTS)
@@ -98,13 +102,16 @@ endif()
set(EIGEN_LAPACK_TARGETS "")
add_library(eigen_lapack_static STATIC ${EigenLapack_SRCS} ${ReferenceLapack_SRCS})
target_link_libraries(eigen_lapack_static eigen_blas_static)
list(APPEND EIGEN_LAPACK_TARGETS eigen_lapack_static)
if (EIGEN_BUILD_SHARED_LIBS)
add_library(eigen_lapack SHARED ${EigenLapack_SRCS})
target_compile_definitions(eigen_lapack PUBLIC "EIGEN_BUILD_DLL")
list(APPEND EIGEN_LAPACK_TARGETS eigen_lapack)
add_library(eigen_lapack SHARED ${EigenLapack_SRCS} ${EIGEN_LAPACK_DEF})
# Build LAPACK but link BLAS.
target_compile_definitions(eigen_lapack PUBLIC "EIGEN_BLAS_LINK_DLL" "EIGEN_LAPACK_BUILD_DLL")
target_link_libraries(eigen_lapack eigen_blas)
set_target_properties(eigen_lapack PROPERTIES CXX_VISIBILITY_PRESET hidden)
list(APPEND EIGEN_LAPACK_TARGETS eigen_lapack)
endif()
foreach(target IN LISTS EIGEN_LAPACK_TARGETS)

5
lapack/cholesky.inc
View File

@@ -11,7 +11,7 @@
#include <Eigen/Cholesky>
// POTRF computes the Cholesky factorization of a real symmetric positive definite matrix A.
EIGEN_LAPACK_FUNC(potrf)(char *uplo, int *n, RealScalar *pa, int *lda, int *info) {
EIGEN_LAPACK_FUNC(potrf)(const char *uplo, int *n, RealScalar *pa, int *lda, int *info) {
*info = 0;
if (UPLO(*uplo) == INVALID)
*info = -1;
@@ -38,7 +38,8 @@ EIGEN_LAPACK_FUNC(potrf)(char *uplo, int *n, RealScalar *pa, int *lda, int *info
// POTRS solves a system of linear equations A*X = B with a symmetric
// positive definite matrix A using the Cholesky factorization
// A = U**T*U or A = L*L**T computed by DPOTRF.
EIGEN_LAPACK_FUNC(potrs)(char *uplo, int *n, int *nrhs, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, int *info) {
EIGEN_LAPACK_FUNC(potrs)
(const char *uplo, int *n, int *nrhs, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, int *info) {
*info = 0;
if (UPLO(*uplo) == INVALID)
*info = -1;

2
lapack/dsecnd_INT_CPU_TIME.cpp
View File

@@ -15,6 +15,8 @@
#include <ctime>
#endif
#include "lapack.h"
extern "C" {
double dsecnd_();
}

143
lapack/eigen_lapack.def Normal file
View File

@@ -0,0 +1,143 @@
; Definition file for eigen_lapack.dll.
LIBRARY eigen_lapack
EXPORTS
; Eigen C/C++ implementations
; Utilities
xerbla_
; Eigenvalues
ssyev_
dsyev_
; LU
sgetrf_
sgetrs_
dgetrf_
dgetrs_
cgetrf_
cgetrs_
zgetrf_
zgetrs_
; QR
spotrf_
spotrs_
dpotrf_
dpotrs_
cpotrf_
cpotrs_
zpotrf_
zpotrs_
; SVD
sgesdd_
sgesvd_
dgesdd_
dgesvd_
cgesdd_
cgesvd_
zgesdd_
zgesvd_
; Time
second_
dsecnd_
; Fortran implementations
clacgv_
zlacgv_
sladiv_
dladiv_
cladiv_
zladiv_
slamch_
dlamch_
slamc3_
dlamc3_
slapy2_
dlapy2_
slapy3_
dlapy3_
slarf_
dlarf_
clarf_
zlarf_
slarfb_
dlarfb_
clarfb_
zlarfb_
slarfg_
dlarfg_
clarfg_
zlarfg_
slarft_
dlarft_
clarft_
zlarft_
ilaclc_
ilaclr_
iladlc_
iladlr_
ilaslc_
ilaslr_
ilazlc_
ilazlr_
; Missing
; csymv_
; zsymv_
; cspmv_
; zspmv_
; csyr_
; zsyr_
; cspr_
; zspr_
; sgemt_
; dgemt_
; cgemt_
; zgemt_
; sgema_
; dgema_
; cgema_
; zgema_
; sgems_
; dgems_
; cgems_
; zgems_
; sgetf2_
; dgetf2_
; cgetf2_
; zgetf2_
; slaswp_
; dlaswp_
; claswp_
; zlaswp_
; sgesv_
; dgesv_
; cgesv_
; zgesv_
; spotf2_
; dpotf2_
; cpotf2_
; zpotf2_
; slauu2_
; dlauu2_
; clauu2_
; zlauu2_
; slauum_
; dlauum_
; clauum_
; zlauum_
; strti2_
; dtrti2_
; ctrti2_
; ztrti2_
; strtri_
; dtrtri_
; ctrtri_
; ztrtri_
; spotri_
; dpotri_
; cpotri_
; zpotri_

143
lapack/eigen_lapack_cpp.def Normal file
View File

@@ -0,0 +1,143 @@
; Definition file for eigen_lapack.dll, containing only the C++ implementations.
LIBRARY eigen_lapack
EXPORTS
; Eigen C/C++ implementations
; Utilities
xerbla_
; Eigenvalues
ssyev_
dsyev_
; LU
sgetrf_
sgetrs_
dgetrf_
dgetrs_
cgetrf_
cgetrs_
zgetrf_
zgetrs_
; QR
spotrf_
spotrs_
dpotrf_
dpotrs_
cpotrf_
cpotrs_
zpotrf_
zpotrs_
; SVD
sgesdd_
sgesvd_
dgesdd_
dgesvd_
cgesdd_
cgesvd_
zgesdd_
zgesvd_
; Time
second_
dsecnd_
; Fortran implementations
; clacgv_
; zlacgv_
; sladiv_
; dladiv_
; cladiv_
; zladiv_
; slamch_
; dlamch_
; slamc3_
; dlamc3_
; slapy2_
; dlapy2_
; slapy3_
; dlapy3_
; slarf_
; dlarf_
; clarf_
; zlarf_
; slarfb_
; dlarfb_
; clarfb_
; zlarfb_
; slarfg_
; dlarfg_
; clarfg_
; zlarfg_
; slarft_
; dlarft_
; clarft_
; zlarft_
; ilaclc_
; ilaclr_
; iladlc_
; iladlr_
; ilaslc_
; ilaslr_
; ilazlc_
; ilazlr_
; Missing
; csymv_
; zsymv_
; cspmv_
; zspmv_
; csyr_
; zsyr_
; cspr_
; zspr_
; sgemt_
; dgemt_
; cgemt_
; zgemt_
; sgema_
; dgema_
; cgema_
; zgema_
; sgems_
; dgems_
; cgems_
; zgems_
; sgetf2_
; dgetf2_
; cgetf2_
; zgetf2_
; slaswp_
; dlaswp_
; claswp_
; zlaswp_
; sgesv_
; dgesv_
; cgesv_
; zgesv_
; spotf2_
; dpotf2_
; cpotf2_
; zpotf2_
; slauu2_
; dlauu2_
; clauu2_
; zlauu2_
; slauum_
; dlauum_
; clauum_
; zlauum_
; strti2_
; dtrti2_
; ctrti2_
; ztrti2_
; strtri_
; dtrtri_
; ctrtri_
; ztrtri_
; spotri_
; dpotri_
; cpotri_
; zpotri_

6
lapack/eigenvalues.inc
View File

@@ -12,7 +12,8 @@
// computes eigen values and vectors of a general N-by-N matrix A
EIGEN_LAPACK_FUNC(syev)
(char* jobz, char* uplo, int* n, Scalar* a, int* lda, Scalar* w, Scalar* /*work*/, int* lwork, int* info) {
(const char* jobz, const char* uplo, int* n, RealScalar* ra, int* lda, RealScalar* rw, RealScalar* /*work*/, int* lwork,
int* info) {
// TODO exploit the work buffer
bool query_size = *lwork == -1;
@@ -40,6 +41,9 @@ EIGEN_LAPACK_FUNC(syev)
if (*n == 0) return;
Scalar* a = reinterpret_cast<Scalar*>(ra);
Scalar* w = reinterpret_cast<Scalar*>(rw);
PlainMatrixType mat(*n, *n);
if (UPLO(*uplo) == UP)
mat = matrix(a, *n, *n, *lda).adjoint();

271
lapack/lapack.h
View File

@@ -3,135 +3,192 @@
#include "../blas/blas.h"
#if defined(_WIN32)
#if defined(EIGEN_LAPACK_BUILD_DLL)
#define EIGEN_LAPACK_API __declspec(dllexport)
#elif defined(EIGEN_LAPACK_LINK_DLL)
#define EIGEN_LAPACK_API __declspec(dllimport)
#else
#define EIGEN_LAPACK_API
#endif
#elif ((defined(__GNUC__) && __GNUC__ >= 4) || defined(__clang__)) && defined(EIGEN_LAPACK_BUILD_DLL)
#define EIGEN_LAPACK_API __attribute__((visibility("default")))
#else
#define EIGEN_LAPACK_API
#endif
#ifdef __cplusplus
extern "C" {
#endif
EIGEN_BLAS_API void BLASFUNC(csymv)(const char *, const int *, const float *, const float *, const int *, const float *,
const int *, const float *, float *, const int *);
EIGEN_BLAS_API void BLASFUNC(zsymv)(const char *, const int *, const double *, const double *, const int *,
const double *, const int *, const double *, double *, const int *);
EIGEN_BLAS_API void BLASFUNC(xsymv)(const char *, const int *, const double *, const double *, const int *,
const double *, const int *, const double *, double *, const int *);
EIGEN_LAPACK_API void BLASFUNC(csymv)(const char *, const int *, const float *, const float *, const int *,
const float *, const int *, const float *, float *, const int *);
EIGEN_LAPACK_API void BLASFUNC(zsymv)(const char *, const int *, const double *, const double *, const int *,
const double *, const int *, const double *, double *, const int *);
EIGEN_LAPACK_API void BLASFUNC(xsymv)(const char *, const int *, const double *, const double *, const int *,
const double *, const int *, const double *, double *, const int *);
EIGEN_BLAS_API void BLASFUNC(cspmv)(char *, int *, float *, float *, float *, int *, float *, float *, int *);
EIGEN_BLAS_API void BLASFUNC(zspmv)(char *, int *, double *, double *, double *, int *, double *, double *, int *);
EIGEN_BLAS_API void BLASFUNC(xspmv)(char *, int *, double *, double *, double *, int *, double *, double *, int *);
EIGEN_LAPACK_API void BLASFUNC(cspmv)(char *, int *, float *, float *, float *, int *, float *, float *, int *);
EIGEN_LAPACK_API void BLASFUNC(zspmv)(char *, int *, double *, double *, double *, int *, double *, double *, int *);
EIGEN_LAPACK_API void BLASFUNC(xspmv)(char *, int *, double *, double *, double *, int *, double *, double *, int *);
EIGEN_BLAS_API void BLASFUNC(csyr)(char *, int *, float *, float *, int *, float *, int *);
EIGEN_BLAS_API void BLASFUNC(zsyr)(char *, int *, double *, double *, int *, double *, int *);
EIGEN_BLAS_API void BLASFUNC(xsyr)(char *, int *, double *, double *, int *, double *, int *);
EIGEN_LAPACK_API void BLASFUNC(csyr)(char *, int *, float *, float *, int *, float *, int *);
EIGEN_LAPACK_API void BLASFUNC(zsyr)(char *, int *, double *, double *, int *, double *, int *);
EIGEN_LAPACK_API void BLASFUNC(xsyr)(char *, int *, double *, double *, int *, double *, int *);
EIGEN_BLAS_API void BLASFUNC(cspr)(char *, int *, float *, float *, int *, float *);
EIGEN_BLAS_API void BLASFUNC(zspr)(char *, int *, double *, double *, int *, double *);
EIGEN_BLAS_API void BLASFUNC(xspr)(char *, int *, double *, double *, int *, double *);
EIGEN_LAPACK_API void BLASFUNC(cspr)(char *, int *, float *, float *, int *, float *);
EIGEN_LAPACK_API void BLASFUNC(zspr)(char *, int *, double *, double *, int *, double *);
EIGEN_LAPACK_API void BLASFUNC(xspr)(char *, int *, double *, double *, int *, double *);
EIGEN_BLAS_API void BLASFUNC(sgemt)(char *, int *, int *, float *, float *, int *, float *, int *);
EIGEN_BLAS_API void BLASFUNC(dgemt)(char *, int *, int *, double *, double *, int *, double *, int *);
EIGEN_BLAS_API void BLASFUNC(cgemt)(char *, int *, int *, float *, float *, int *, float *, int *);
EIGEN_BLAS_API void BLASFUNC(zgemt)(char *, int *, int *, double *, double *, int *, double *, int *);
EIGEN_LAPACK_API void BLASFUNC(sgemt)(char *, int *, int *, float *, float *, int *, float *, int *);
EIGEN_LAPACK_API void BLASFUNC(dgemt)(char *, int *, int *, double *, double *, int *, double *, int *);
EIGEN_LAPACK_API void BLASFUNC(cgemt)(char *, int *, int *, float *, float *, int *, float *, int *);
EIGEN_LAPACK_API void BLASFUNC(zgemt)(char *, int *, int *, double *, double *, int *, double *, int *);
EIGEN_BLAS_API void BLASFUNC(sgema)(char *, char *, int *, int *, float *, float *, int *, float *, float *, int *,
float *, int *);
EIGEN_BLAS_API void BLASFUNC(dgema)(char *, char *, int *, int *, double *, double *, int *, double *, double *, int *,
double *, int *);
EIGEN_BLAS_API void BLASFUNC(cgema)(char *, char *, int *, int *, float *, float *, int *, float *, float *, int *,
float *, int *);
EIGEN_BLAS_API void BLASFUNC(zgema)(char *, char *, int *, int *, double *, double *, int *, double *, double *, int *,
double *, int *);
EIGEN_LAPACK_API void BLASFUNC(sgema)(char *, char *, int *, int *, float *, float *, int *, float *, float *, int *,
float *, int *);
EIGEN_LAPACK_API void BLASFUNC(dgema)(char *, char *, int *, int *, double *, double *, int *, double *, double *,
int *, double *, int *);
EIGEN_LAPACK_API void BLASFUNC(cgema)(char *, char *, int *, int *, float *, float *, int *, float *, float *, int *,
float *, int *);
EIGEN_LAPACK_API void BLASFUNC(zgema)(char *, char *, int *, int *, double *, double *, int *, double *, double *,
int *, double *, int *);
EIGEN_BLAS_API void BLASFUNC(sgems)(char *, char *, int *, int *, float *, float *, int *, float *, float *, int *,
float *, int *);
EIGEN_BLAS_API void BLASFUNC(dgems)(char *, char *, int *, int *, double *, double *, int *, double *, double *, int *,
double *, int *);
EIGEN_BLAS_API void BLASFUNC(cgems)(char *, char *, int *, int *, float *, float *, int *, float *, float *, int *,
float *, int *);
EIGEN_BLAS_API void BLASFUNC(zgems)(char *, char *, int *, int *, double *, double *, int *, double *, double *, int *,
double *, int *);
EIGEN_LAPACK_API void BLASFUNC(sgems)(char *, char *, int *, int *, float *, float *, int *, float *, float *, int *,
float *, int *);
EIGEN_LAPACK_API void BLASFUNC(dgems)(char *, char *, int *, int *, double *, double *, int *, double *, double *,
int *, double *, int *);
EIGEN_LAPACK_API void BLASFUNC(cgems)(char *, char *, int *, int *, float *, float *, int *, float *, float *, int *,
float *, int *);
EIGEN_LAPACK_API void BLASFUNC(zgems)(char *, char *, int *, int *, double *, double *, int *, double *, double *,
int *, double *, int *);
EIGEN_BLAS_API void BLASFUNC(sgetf2)(int *, int *, float *, int *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(dgetf2)(int *, int *, double *, int *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(qgetf2)(int *, int *, double *, int *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(cgetf2)(int *, int *, float *, int *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(zgetf2)(int *, int *, double *, int *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(xgetf2)(int *, int *, double *, int *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(sgetf2)(int *, int *, float *, int *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(dgetf2)(int *, int *, double *, int *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(qgetf2)(int *, int *, double *, int *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(cgetf2)(int *, int *, float *, int *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(zgetf2)(int *, int *, double *, int *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(xgetf2)(int *, int *, double *, int *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(sgetrf)(int *, int *, float *, int *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(dgetrf)(int *, int *, double *, int *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(qgetrf)(int *, int *, double *, int *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(cgetrf)(int *, int *, float *, int *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(zgetrf)(int *, int *, double *, int *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(xgetrf)(int *, int *, double *, int *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(qgetrf)(int *, int *, double *, int *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(xgetrf)(int *, int *, double *, int *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(slaswp)(int *, float *, int *, int *, int *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(dlaswp)(int *, double *, int *, int *, int *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(qlaswp)(int *, double *, int *, int *, int *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(claswp)(int *, float *, int *, int *, int *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(zlaswp)(int *, double *, int *, int *, int *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(xlaswp)(int *, double *, int *, int *, int *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(slaswp)(int *, float *, int *, int *, int *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(dlaswp)(int *, double *, int *, int *, int *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(qlaswp)(int *, double *, int *, int *, int *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(claswp)(int *, float *, int *, int *, int *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(zlaswp)(int *, double *, int *, int *, int *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(xlaswp)(int *, double *, int *, int *, int *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(sgetrs)(char *, int *, int *, float *, int *, int *, float *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(dgetrs)(char *, int *, int *, double *, int *, int *, double *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(qgetrs)(char *, int *, int *, double *, int *, int *, double *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(cgetrs)(char *, int *, int *, float *, int *, int *, float *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(zgetrs)(char *, int *, int *, double *, int *, int *, double *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(xgetrs)(char *, int *, int *, double *, int *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(qgetrs)(char *, int *, int *, double *, int *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(xgetrs)(char *, int *, int *, double *, int *, int *, double *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(sgesv)(int *, int *, float *, int *, int *, float *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(dgesv)(int *, int *, double *, int *, int *, double *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(qgesv)(int *, int *, double *, int *, int *, double *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(cgesv)(int *, int *, float *, int *, int *, float *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(zgesv)(int *, int *, double *, int *, int *, double *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(xgesv)(int *, int *, double *, int *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(sgesv)(int *, int *, float *, int *, int *, float *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(dgesv)(int *, int *, double *, int *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(qgesv)(int *, int *, double *, int *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(cgesv)(int *, int *, float *, int *, int *, float *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(zgesv)(int *, int *, double *, int *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(xgesv)(int *, int *, double *, int *, int *, double *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(spotf2)(char *, int *, float *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(dpotf2)(char *, int *, double *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(qpotf2)(char *, int *, double *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(cpotf2)(char *, int *, float *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(zpotf2)(char *, int *, double *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(xpotf2)(char *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(spotf2)(char *, int *, float *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(dpotf2)(char *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(qpotf2)(char *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(cpotf2)(char *, int *, float *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(zpotf2)(char *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(xpotf2)(char *, int *, double *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(spotrf)(char *, int *, float *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(dpotrf)(char *, int *, double *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(qpotrf)(char *, int *, double *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(cpotrf)(char *, int *, float *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(zpotrf)(char *, int *, double *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(xpotrf)(char *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(qpotrf)(char *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(xpotrf)(char *, int *, double *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(slauu2)(char *, int *, float *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(dlauu2)(char *, int *, double *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(qlauu2)(char *, int *, double *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(clauu2)(char *, int *, float *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(zlauu2)(char *, int *, double *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(xlauu2)(char *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(slauu2)(char *, int *, float *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(dlauu2)(char *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(qlauu2)(char *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(clauu2)(char *, int *, float *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(zlauu2)(char *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(xlauu2)(char *, int *, double *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(slauum)(char *, int *, float *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(dlauum)(char *, int *, double *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(qlauum)(char *, int *, double *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(clauum)(char *, int *, float *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(zlauum)(char *, int *, double *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(xlauum)(char *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(slauum)(char *, int *, float *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(dlauum)(char *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(qlauum)(char *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(clauum)(char *, int *, float *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(zlauum)(char *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(xlauum)(char *, int *, double *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(strti2)(char *, char *, int *, float *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(dtrti2)(char *, char *, int *, double *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(qtrti2)(char *, char *, int *, double *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(ctrti2)(char *, char *, int *, float *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(ztrti2)(char *, char *, int *, double *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(xtrti2)(char *, char *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(strti2)(char *, char *, int *, float *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(dtrti2)(char *, char *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(qtrti2)(char *, char *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(ctrti2)(char *, char *, int *, float *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(ztrti2)(char *, char *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(xtrti2)(char *, char *, int *, double *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(strtri)(char *, char *, int *, float *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(dtrtri)(char *, char *, int *, double *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(qtrtri)(char *, char *, int *, double *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(ctrtri)(char *, char *, int *, float *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(ztrtri)(char *, char *, int *, double *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(xtrtri)(char *, char *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(strtri)(char *, char *, int *, float *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(dtrtri)(char *, char *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(qtrtri)(char *, char *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(ctrtri)(char *, char *, int *, float *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(ztrtri)(char *, char *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(xtrtri)(char *, char *, int *, double *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(spotri)(char *, int *, float *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(dpotri)(char *, int *, double *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(qpotri)(char *, int *, double *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(cpotri)(char *, int *, float *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(zpotri)(char *, int *, double *, int *, int *);
EIGEN_BLAS_API void BLASFUNC(xpotri)(char *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(spotri)(char *, int *, float *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(dpotri)(char *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(qpotri)(char *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(cpotri)(char *, int *, float *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(zpotri)(char *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(xpotri)(char *, int *, double *, int *, int *);
//-----------------------------------------------------------------------------
// Eigen C++ implementations.
//-----------------------------------------------------------------------------
// Cholesky.
EIGEN_LAPACK_API void BLASFUNC(spotrf)(const char *, int *, float *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(dpotrf)(const char *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(cpotrf)(const char *, int *, float *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(zpotrf)(const char *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(spotrs)(const char *, int *, int *, float *, int *, float *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(dpotrs)(const char *, int *, int *, double *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(cpotrs)(const char *, int *, int *, float *, int *, float *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(zpotrs)(const char *, int *, int *, double *, int *, double *, int *, int *);
// Eigenvalues.
EIGEN_LAPACK_API void BLASFUNC(ssyev)(const char *, const char *, int *, float *, int *, float *, float *, int *,
int *);
EIGEN_LAPACK_API void BLASFUNC(dsyev)(const char *, const char *, int *, double *, int *, double *, double *, int *,
int *);
// LU.
EIGEN_LAPACK_API void BLASFUNC(sgetrf)(int *, int *, float *, int *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(dgetrf)(int *, int *, double *, int *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(cgetrf)(int *, int *, float *, int *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(zgetrf)(int *, int *, double *, int *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(sgetrs)(const char *, int *, int *, float *, int *, int *, float *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(dgetrs)(const char *, int *, int *, double *, int *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(cgetrs)(const char *, int *, int *, float *, int *, int *, float *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(zgetrs)(const char *, int *, int *, double *, int *, int *, double *, int *, int *);
// SVD.
EIGEN_LAPACK_API void BLASFUNC(sgesdd)(const char *, int *, int *, float *, int *, float *, float *, int *, float *,
int *, float *, int *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(dgesdd)(const char *, int *, int *, double *, int *, double *, double *, int *, double *,
int *, double *, int *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(cgesdd)(const char *, int *, int *, float *, int *, float *, float *, int *, float *,
int *, float *, int *, float *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(zgesdd)(const char *, int *, int *, double *, int *, double *, double *, int *, double *,
int *, double *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(sgesvd)(const char *, const char *, int *, int *, float *, int *, float *, float *,
int *, float *, int *, float *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(dgesvd)(const char *, const char *, int *, int *, double *, int *, double *, double *,
int *, double *, int *, double *, int *, int *);
EIGEN_LAPACK_API void BLASFUNC(cgesvd)(const char *, const char *, int *, int *, float *, int *, float *, float *,
int *, float *, int *, float *, int *, float *, int *);
EIGEN_LAPACK_API void BLASFUNC(zgesvd)(const char *, const char *, int *, int *, double *, int *, double *, double *,
int *, double *, int *, double *, int *, double *, int *);
// Time.
EIGEN_LAPACK_API float BLASFUNC(second)();
EIGEN_LAPACK_API double BLASFUNC(dsecnd)();
#ifdef __cplusplus
}

2
lapack/lu.inc
View File

@@ -40,7 +40,7 @@ EIGEN_LAPACK_FUNC(getrf)(int *m, int *n, RealScalar *pa, int *lda, int *ipiv, in
// A * X = B or A' * X = B
// with a general N-by-N matrix A using the LU factorization computed by GETRF
EIGEN_LAPACK_FUNC(getrs)
(char *trans, int *n, int *nrhs, RealScalar *pa, int *lda, int *ipiv, RealScalar *pb, int *ldb, int *info) {
(const char *trans, int *n, int *nrhs, RealScalar *pa, int *lda, int *ipiv, RealScalar *pb, int *ldb, int *info) {
*info = 0;
if (OP(*trans) == INVALID)
*info = -1;

2
lapack/second_INT_CPU_TIME.cpp
View File

@@ -15,6 +15,8 @@
#include <ctime>
#endif
#include "lapack.h"
extern "C" {
float second_();
}

18
lapack/svd.inc
View File

@@ -18,8 +18,9 @@
// computes the singular values/vectors a general M-by-N matrix A using divide-and-conquer
EIGEN_LAPACK_FUNC(gesdd)
(char *jobz, int *m, int *n, Scalar *a, int *lda, RealScalar *s, Scalar *u, int *ldu, Scalar *vt, int *ldvt,
Scalar * /*work*/, int *lwork, EIGEN_LAPACK_ARG_IF_COMPLEX(RealScalar * /*rwork*/) int * /*iwork*/, int *info) {
(const char *jobz, int *m, int *n, RealScalar *ra, int *lda, RealScalar *s, RealScalar *ru, int *ldu, RealScalar *rvt,
int *ldvt, RealScalar * /*work*/, int *lwork, EIGEN_LAPACK_ARG_IF_COMPLEX(RealScalar * /*rwork*/) int * /*iwork*/,
int *info) {
// TODO exploit the work buffer
bool query_size = *lwork == -1;
int diag_size = (std::min)(*m, *n);
@@ -53,6 +54,10 @@ EIGEN_LAPACK_FUNC(gesdd)
if (*n == 0 || *m == 0) return;
Scalar *a = reinterpret_cast<Scalar *>(ra);
Scalar *u = reinterpret_cast<Scalar *>(ru);
Scalar *vt = reinterpret_cast<Scalar *>(rvt);
PlainMatrixType mat(*m, *n);
mat = matrix(a, *m, *n, *lda);
@@ -84,8 +89,9 @@ EIGEN_LAPACK_FUNC(gesdd)
// computes the singular values/vectors a general M-by-N matrix A using two sided jacobi algorithm
EIGEN_LAPACK_FUNC(gesvd)
(char *jobu, char *jobv, int *m, int *n, Scalar *a, int *lda, RealScalar *s, Scalar *u, int *ldu, Scalar *vt, int *ldvt,
Scalar * /*work*/, int *lwork, EIGEN_LAPACK_ARG_IF_COMPLEX(RealScalar * /*rwork*/) int *info) {
(const char *jobu, const char *jobv, int *m, int *n, RealScalar *ra, int *lda, RealScalar *s, RealScalar *ru, int *ldu,
RealScalar *rvt, int *ldvt, RealScalar * /*work*/, int *lwork,
EIGEN_LAPACK_ARG_IF_COMPLEX(RealScalar * /*rwork*/) int *info) {
// TODO exploit the work buffer
bool query_size = *lwork == -1;
int diag_size = (std::min)(*m, *n);
@@ -118,6 +124,10 @@ EIGEN_LAPACK_FUNC(gesvd)
if (*n == 0 || *m == 0) return;
Scalar *a = reinterpret_cast<Scalar *>(ra);
Scalar *u = reinterpret_cast<Scalar *>(ru);
Scalar *vt = reinterpret_cast<Scalar *>(rvt);
PlainMatrixType mat(*m, *n);
mat = matrix(a, *m, *n, *lda);

1
test/CMakeLists.txt
View File

@@ -257,6 +257,7 @@ ei_add_test(eigensolver_selfadjoint)
ei_add_test(eigensolver_generic)
ei_add_test(eigensolver_complex)
ei_add_test(real_qz)
ei_add_test(complex_qz)
ei_add_test(eigensolver_generalized_real)
ei_add_test(jacobi)
ei_add_test(jacobisvd)

86
test/complex_qz.cpp Normal file
View File

@@ -0,0 +1,86 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2012 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/.
#define EIGEN_RUNTIME_NO_MALLOC
#include "main.h"
#include <Eigen/Eigenvalues>
/* this test covers the following files:
ComplexQZ.h
*/
template <typename MatrixType>
void generate_random_matrix_pair(const Index dim, MatrixType& A, MatrixType& B) {
A.resize(dim, dim);
B.resize(dim, dim);
A.setRandom();
B.setRandom();
// Set each row of B with a probability of 10% to 0
for (int i = 0; i < dim; i++) {
if (internal::random<int>(0, 10) == 0) B.row(i).setZero();
}
}
template <typename MatrixType>
void complex_qz(const MatrixType& A, const MatrixType& B) {
using std::abs;
const Index dim = A.rows();
ComplexQZ<MatrixType> qz(A, B);
VERIFY_IS_EQUAL(qz.info(), Success);
auto T = qz.matrixT(), S = qz.matrixS();
bool is_all_zero_T = true, is_all_zero_S = true;
using RealScalar = typename MatrixType::RealScalar;
RealScalar tol = dim * 10 * NumTraits<RealScalar>::epsilon();
for (Index j = 0; j < dim; j++) {
for (Index i = j + 1; i < dim; i++) {
if (std::abs(T(i, j)) > tol) {
std::cerr << std::abs(T(i, j)) << std::endl;
is_all_zero_T = false;
}
if (std::abs(S(i, j)) > tol) {
std::cerr << std::abs(S(i, j)) << std::endl;
is_all_zero_S = false;
}
}
}
VERIFY_IS_EQUAL(is_all_zero_T, true);
VERIFY_IS_EQUAL(is_all_zero_S, true);
VERIFY_IS_APPROX(qz.matrixQ() * qz.matrixS() * qz.matrixZ(), A);
VERIFY_IS_APPROX(qz.matrixQ() * qz.matrixT() * qz.matrixZ(), B);
VERIFY_IS_APPROX(qz.matrixQ() * qz.matrixQ().adjoint(), MatrixType::Identity(dim, dim));
VERIFY_IS_APPROX(qz.matrixZ() * qz.matrixZ().adjoint(), MatrixType::Identity(dim, dim));
}
EIGEN_DECLARE_TEST(complex_qz) {
// const Index dim1 = 15;
// const Index dim2 = 80;
for (int i = 0; i < g_repeat; i++) {
// Check for very small, fixed-sized double- and float complex matrices
Eigen::Matrix2cd A_2x2, B_2x2;
A_2x2.setRandom();
B_2x2.setRandom();
B_2x2.row(1).setZero();
Eigen::Matrix3cf A_3x3, B_3x3;
A_3x3.setRandom();
B_3x3.setRandom();
B_3x3.col(i % 3).setRandom();
// Test for small float complex matrices
Eigen::MatrixXcf A_float, B_float;
const Index dim1 = internal::random<Index>(15, 80), dim2 = internal::random<Index>(15, 80);
generate_random_matrix_pair(dim1, A_float, B_float);
// Test for a bit larger double complex matrices
Eigen::MatrixXcd A_double, B_double;
generate_random_matrix_pair(dim2, A_double, B_double);
CALL_SUBTEST_1(complex_qz(A_2x2, B_2x2));
CALL_SUBTEST_2(complex_qz(A_3x3, B_3x3));
CALL_SUBTEST_3(complex_qz(A_float, B_float));
CALL_SUBTEST_4(complex_qz(A_double, B_double));
}
}