GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-14 12:46:00 +08:00
Code Issues Packages Projects Releases Wiki Activity

bug #1544: Generate correct Q matrix in complex case. Original patch was by Jeff Trull in PR-386.

Browse Source
This commit is contained in:
Christoph Hertzberg 2018-05-17 19:17:01 +02:00
parent 927d023cea
commit 39125654ce
1 changed files with 9 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
Eigen/src/SparseQR/SparseQR.h
View File

@ -52,7 +52,7 @@ namespace internal {
* rank-revealing permutations. Use colsPermutation() to get it. * rank-revealing permutations. Use colsPermutation() to get it.
* *
* Q is the orthogonal matrix represented as products of Householder reflectors. * Q is the orthogonal matrix represented as products of Householder reflectors.
* Use matrixQ() to get an expression and matrixQ().transpose() to get the transpose. * Use matrixQ() to get an expression and matrixQ().adjoint() to get the adjoint.
* You can then apply it to a vector. * You can then apply it to a vector.
* *
* R is the sparse triangular or trapezoidal matrix. The later occurs when A is rank-deficient. * R is the sparse triangular or trapezoidal matrix. The later occurs when A is rank-deficient.
@ -65,6 +65,7 @@ namespace internal {
* \implsparsesolverconcept * \implsparsesolverconcept
* *
* \warning The input sparse matrix A must be in compressed mode (see SparseMatrix::makeCompressed()). * \warning The input sparse matrix A must be in compressed mode (see SparseMatrix::makeCompressed()).
* \warning For complex matrices matrixQ().transpose() will actually return the adjoint matrix.
* *
*/ */
template<typename _MatrixType, typename _OrderingType> template<typename _MatrixType, typename _OrderingType>
@ -196,9 +197,9 @@ class SparseQR : public SparseSolverBase<SparseQR<_MatrixType,_OrderingType> >
Index rank = this->rank(); Index rank = this->rank();
// Compute Q^T * b; // Compute Q^* * b;
typename Dest::PlainObject y, b; typename Dest::PlainObject y, b;
y = this->matrixQ().transpose() * B; y = this->matrixQ().adjoint() * B;
b = y; b = y;
// Solve with the triangular matrix R // Solve with the triangular matrix R
@ -641,7 +642,7 @@ struct SparseQR_QProduct : ReturnByValue<SparseQR_QProduct<SparseQRType, Derived
Scalar tau = Scalar(0); Scalar tau = Scalar(0);
tau = m_qr.m_Q.col(k).dot(res.col(j)); tau = m_qr.m_Q.col(k).dot(res.col(j));
if(tau==Scalar(0)) continue; if(tau==Scalar(0)) continue;
tau = tau * m_qr.m_hcoeffs(k); tau = tau * numext::conj(m_qr.m_hcoeffs(k));
res.col(j) -= tau * m_qr.m_Q.col(k); res.col(j) -= tau * m_qr.m_Q.col(k);
} }
} }
@ -650,7 +651,7 @@ struct SparseQR_QProduct : ReturnByValue<SparseQR_QProduct<SparseQRType, Derived
const SparseQRType& m_qr; const SparseQRType& m_qr;
const Derived& m_other; const Derived& m_other;
bool m_transpose; bool m_transpose; // TODO this actually means adjoint
}; };
template<typename SparseQRType> template<typename SparseQRType>
@ -668,13 +669,14 @@ struct SparseQRMatrixQReturnType : public EigenBase<SparseQRMatrixQReturnType<Sp
{ {
return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(),false); return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(),false);
} }
// To use for operations with the adjoint of Q
SparseQRMatrixQTransposeReturnType<SparseQRType> adjoint() const SparseQRMatrixQTransposeReturnType<SparseQRType> adjoint() const
{ {
return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr); return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
} }
inline Index rows() const { return m_qr.rows(); } inline Index rows() const { return m_qr.rows(); }
inline Index cols() const { return (std::min)(m_qr.rows(),m_qr.cols()); } inline Index cols() const { return (std::min)(m_qr.rows(),m_qr.cols()); }
// To use for operations with the transpose of Q // To use for operations with the transpose of Q FIXME this is the same as adjoint at the moment
SparseQRMatrixQTransposeReturnType<SparseQRType> transpose() const SparseQRMatrixQTransposeReturnType<SparseQRType> transpose() const
{ {
return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr); return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
@ -682,6 +684,7 @@ struct SparseQRMatrixQReturnType : public EigenBase<SparseQRMatrixQReturnType<Sp
const SparseQRType& m_qr; const SparseQRType& m_qr;
}; };
// TODO this actually represents the adjoint of Q
template<typename SparseQRType> template<typename SparseQRType>
struct SparseQRMatrixQTransposeReturnType struct SparseQRMatrixQTransposeReturnType
{ {