From f350f34560d1cc67a8c220d973b003226ff58892 Mon Sep 17 00:00:00 2001 From: Desire NUENTSA Date: Wed, 20 Mar 2013 18:38:22 +0100 Subject: [PATCH] Add complex support to dgmres and the unit test --- .../Eigen/src/IterativeSolvers/DGMRES.h | 78 +++++++++++-------- unsupported/test/dgmres.cpp | 31 ++++++++ 2 files changed, 77 insertions(+), 32 deletions(-) create mode 100644 unsupported/test/dgmres.cpp diff --git a/unsupported/Eigen/src/IterativeSolvers/DGMRES.h b/unsupported/Eigen/src/IterativeSolvers/DGMRES.h index 2ea0eccf1..7b5b5a91b 100644 --- a/unsupported/Eigen/src/IterativeSolvers/DGMRES.h +++ b/unsupported/Eigen/src/IterativeSolvers/DGMRES.h @@ -41,7 +41,6 @@ void sortWithPermutation (VectorType& vec, IndexType& perm, typename IndexType:: eigen_assert(vec.size() == perm.size()); typedef typename IndexType::Scalar Index; typedef typename VectorType::Scalar Scalar; - Index n = vec.size(); bool flag; for (Index k = 0; k < ncut; k++) { @@ -115,8 +114,10 @@ class DGMRES : public IterativeSolverBase > typedef typename MatrixType::RealScalar RealScalar; typedef _Preconditioner Preconditioner; typedef Matrix DenseMatrix; - typedef Matrix DenseVector; - typedef std::complex ComplexScalar; + typedef Matrix DenseRealMatrix; + typedef Matrix DenseVector; + typedef Matrix DenseRealVector; + typedef Matrix, Dynamic, 1> ComplexVector; /** Default constructor. */ @@ -220,6 +221,8 @@ class DGMRES : public IterativeSolverBase > // Apply deflation to a vector template int dgmresApplyDeflation(const RhsType& In, DestType& Out) const; + ComplexVector schurValues(const ComplexSchur& schurofH) const; + ComplexVector schurValues(const RealSchur& schurofH) const; // Init data for deflation void dgmresInitDeflation(Index& rows) const; mutable DenseMatrix m_V; // Krylov basis vectors @@ -307,9 +310,7 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresCycle(const MatrixType& mat, con int n = mat.rows(); DenseVector tv1(n), tv2(n); //Temporary vectors while (m_info == NoConvergence && it < m_restart && nbIts < m_iterations) - { - int n = m_V.rows(); - + { // Apply preconditioner(s) at right if (m_isDeflInitialized ) { @@ -323,7 +324,7 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresCycle(const MatrixType& mat, con tv1 = mat * tv2; // Orthogonalize it with the previous basis in the basis using modified Gram-Schmidt - RealScalar coef; + Scalar coef; for (int i = 0; i <= it; ++i) { coef = tv1.dot(m_V.col(i)); @@ -398,58 +399,71 @@ void DGMRES<_MatrixType, _Preconditioner>::dgmresInitDeflation(Index& rows) cons } template< typename _MatrixType, typename _Preconditioner> -int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, Index& neig) const +inline typename DGMRES<_MatrixType, _Preconditioner>::ComplexVector DGMRES<_MatrixType, _Preconditioner>::schurValues(const ComplexSchur& schurofH) const { - // First, find the Schur form of the Hessenberg matrix H - RealSchur schurofH; - bool computeU = true; - DenseMatrix matrixQ(it,it); - matrixQ.setIdentity(); - schurofH.computeHessenberg(m_Hes.topLeftCorner(it,it), matrixQ, computeU); + return schurofH.matrixT().diagonal(); +} + +template< typename _MatrixType, typename _Preconditioner> +inline typename DGMRES<_MatrixType, _Preconditioner>::ComplexVector DGMRES<_MatrixType, _Preconditioner>::schurValues(const RealSchur& schurofH) const +{ + typedef typename MatrixType::Index Index; const DenseMatrix& T = schurofH.matrixT(); - - // Extract the schur values from the diagonal of T; - Matrix eig(it); - Matrixperm(it); - int j = 0; + Index it = T.rows(); + ComplexVector eig(it); + Index j = 0; while (j < it-1) { if (T(j+1,j) ==Scalar(0)) { - eig(j) = ComplexScalar(T(j,j),Scalar(0)); + eig(j) = std::complex(T(j,j),RealScalar(0)); j++; } else { - eig(j) = ComplexScalar(T(j,j),T(j+1,j)); - eig(j+1) = ComplexScalar(T(j,j+1),T(j+1,j+1)); + eig(j) = std::complex(T(j,j),T(j+1,j)); + eig(j+1) = std::complex(T(j,j+1),T(j+1,j+1)); j++; } } - if (j < it) eig(j) = ComplexScalar(T(j,j),Scalar(0)); + if (j < it-1) eig(j) = std::complex(T(j,j),RealScalar(0)); + return eig; +} + +template< typename _MatrixType, typename _Preconditioner> +int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, Index& neig) const +{ + // First, find the Schur form of the Hessenberg matrix H + typename internal::conditional::IsComplex, ComplexSchur, RealSchur >::type schurofH; + bool computeU = true; + DenseMatrix matrixQ(it,it); + matrixQ.setIdentity(); + schurofH.computeFromHessenberg(m_Hes.topLeftCorner(it,it), matrixQ, computeU); + + ComplexVector eig(it); + Matrixperm(it); + eig = this->schurValues(schurofH); // Reorder the absolute values of Schur values - DenseVector modulEig(it); + DenseRealVector modulEig(it); for (int j=0; j m_lambdaN) - m_lambdaN = modulEig(i); + m_lambdaN = (std::max)(modulEig.maxCoeff(), m_lambdaN); } //Count the real number of extracted eigenvalues (with complex conjugates) int nbrEig = 0; while (nbrEig < neig) { - if(eig(perm(it-nbrEig-1)).imag() == Scalar(0)) nbrEig++; + if(eig(perm(it-nbrEig-1)).imag() == RealScalar(0)) nbrEig++; else nbrEig += 2; } - // Extract the smallest Schur vectors + // Extract the Schur vectors corresponding to the smallest Ritz values DenseMatrix Sr(it, nbrEig); + Sr.setZero(); for (int j = 0; j < nbrEig; j++) { Sr.col(j) = schurofH.matrixU().col(perm(it-j-1)); @@ -478,7 +492,7 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const Matri MX.col(j) = precond.solve(tv1); } - //Update T = [U'MU U'MX; X'MU X'MX] + //Update m_T = [U'MU U'MX; X'MU X'MX] m_T.block(m_r, m_r, nbrEig, nbrEig) = X.transpose() * MX; if(m_r) { @@ -525,4 +539,4 @@ struct solve_retval, Rhs> } // end namespace internal } // end namespace Eigen -#endif +#endif \ No newline at end of file diff --git a/unsupported/test/dgmres.cpp b/unsupported/test/dgmres.cpp new file mode 100644 index 000000000..2b11807c8 --- /dev/null +++ b/unsupported/test/dgmres.cpp @@ -0,0 +1,31 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2011 Gael Guennebaud +// Copyright (C) 2012 desire Nuentsa + +template void test_dgmres_T() +{ + DGMRES, DiagonalPreconditioner > dgmres_colmajor_diag; + DGMRES, IdentityPreconditioner > dgmres_colmajor_I; + DGMRES, IncompleteLUT > dgmres_colmajor_ilut; + //GMRES, SSORPreconditioner > dgmres_colmajor_ssor; + + CALL_SUBTEST( check_sparse_square_solving(dgmres_colmajor_diag) ); +// CALL_SUBTEST( check_sparse_square_solving(dgmres_colmajor_I) ); + CALL_SUBTEST( check_sparse_square_solving(dgmres_colmajor_ilut) ); + //CALL_SUBTEST( check_sparse_square_solving(dgmres_colmajor_ssor) ); +} + +void test_dgmres() +{ + CALL_SUBTEST_1(test_dgmres_T()); + CALL_SUBTEST_2(test_dgmres_T >()); +}