From 73c658c711138a9f0e1627f44731683444f69e41 Mon Sep 17 00:00:00 2001 From: Jack Poulson Date: Thu, 9 Mar 2017 07:58:01 -0800 Subject: [PATCH] Fixing issue #222 by avoiding the QR algorithm calls within D&C from only computing subsets --- src/lapack_like/spectral/BidiagSVD.cpp | 2 +- src/lapack_like/spectral/HermitianEig.cpp | 58 +++--- .../spectral/HermitianTridiagEig.cpp | 193 +++++++++++------- .../HermitianTridiagEig/DivideAndConquer.hpp | 24 +-- 4 files changed, 161 insertions(+), 116 deletions(-) diff --git a/src/lapack_like/spectral/BidiagSVD.cpp b/src/lapack_like/spectral/BidiagSVD.cpp index 5d40414ddd..62fb8e2de0 100644 --- a/src/lapack_like/spectral/BidiagSVD.cpp +++ b/src/lapack_like/spectral/BidiagSVD.cpp @@ -1,5 +1,5 @@ /* - Copyright (c) 2009-2016, Jack Poulson + Copyright (c) 2009-2017, Jack Poulson All rights reserved. This file is part of Elemental and is under the BSD 2-Clause License, diff --git a/src/lapack_like/spectral/HermitianEig.cpp b/src/lapack_like/spectral/HermitianEig.cpp index 5ab31344d9..556b72a595 100644 --- a/src/lapack_like/spectral/HermitianEig.cpp +++ b/src/lapack_like/spectral/HermitianEig.cpp @@ -1,9 +1,9 @@ /* - Copyright (c) 2009-2016, Jack Poulson + Copyright (c) 2009-2017, Jack Poulson All rights reserved. - This file is part of Elemental and is under the BSD 2-Clause License, - which can be found in the LICENSE file in the root directory, or at + This file is part of Elemental and is under the BSD 2-Clause License, + which can be found in the LICENSE file in the root directory, or at http://opensource.org/licenses/BSD-2-Clause */ #include @@ -78,8 +78,8 @@ MRRRPostEstimate namespace herm_eig { -// We create specialized redistribution routines for redistributing the -// real eigenvectors of the symmetric tridiagonal matrix at the core of our +// We create specialized redistribution routines for redistributing the +// real eigenvectors of the symmetric tridiagonal matrix at the core of our // eigensolver in order to minimize the temporary memory usage. template void InPlaceRedist( DistMatrix& Q, Int rowAlign, const Base* readBuffer ) @@ -101,7 +101,7 @@ void InPlaceRedist( DistMatrix& Q, Int rowAlign, const Base* readBuffer ) const Int maxHeight = MaxLength(height,r); const Int maxWidth = MaxLength(width,p); const Int portionSize = mpi::Pad( maxHeight*maxWidth ); - + // Allocate our send/recv buffers vector buffer(2*r*portionSize); Real* sendBuffer = &buffer[0]; @@ -119,7 +119,7 @@ void InPlaceRedist( DistMatrix& Q, Int rowAlign, const Base* readBuffer ) EL_INNER_PARALLEL_FOR_COLLAPSE2 for( Int j=0; j HermitianEigInfo BlackBox ( UpperOrLower uplo, - Matrix& A, + Matrix& A, Matrix>& w, const HermitianEigCtrl& ctrl ) { @@ -205,13 +205,13 @@ BlackBox const Int n = A.Height(); Int numValid = n; if( subset.indexSubset ) - { + { numValid = subset.upperIndex-subset.lowerIndex+1; } else if( subset.rangeSubset ) { if( subset.lowerBound >= Real(0) || subset.upperBound < Real(0) ) - { + { numValid = 0; } } @@ -256,7 +256,7 @@ template HermitianEigInfo HermitianEig ( UpperOrLower uplo, - Matrix& A, + Matrix& A, Matrix>& w, const HermitianEigCtrl& ctrl ) { @@ -279,7 +279,7 @@ template HermitianEigInfo SequentialHelper ( UpperOrLower uplo, - AbstractDistMatrix& A, + AbstractDistMatrix& A, AbstractDistMatrix>& w, const HermitianEigCtrl& ctrl ) { @@ -426,7 +426,7 @@ BlackBox else if( subset.rangeSubset ) { if( subset.lowerBound >= Real(0) || subset.upperBound < Real(0) ) - { + { numValid = 0; } } @@ -451,7 +451,7 @@ BlackBox if( A.Grid().Rank() == 0 ) timer.Start(); } - + // Tridiagonalize A herm_tridiag::ExplicitCondensed( uplo, A, ctrl.tridiagCtrl ); @@ -540,7 +540,7 @@ BlackBox ( UpperOrLower uplo, Matrix& A, Matrix>& w, - Matrix& Q, + Matrix& Q, const HermitianEigCtrl& ctrl ) { EL_DEBUG_CSE @@ -566,14 +566,14 @@ BlackBox ( UpperOrLower uplo, AbstractDistMatrix& APre, AbstractDistMatrix>& w, - AbstractDistMatrix& QPre, + AbstractDistMatrix& QPre, const HermitianEigCtrl& ctrl ) { EL_DEBUG_CSE const Grid& g = APre.Grid(); HermitianEigInfo info; - DistMatrixReadProxy AProx( APre ); + DistMatrixReadProxy AProx( APre ); auto& A = AProx.Get(); // TODO(poulson): Extend interface to support ctrl.tridiagCtrl @@ -613,7 +613,7 @@ HermitianEig ( UpperOrLower uplo, Matrix& A, Matrix>& w, - Matrix& Q, + Matrix& Q, const HermitianEigCtrl& ctrl ) { EL_DEBUG_CSE @@ -648,13 +648,13 @@ HermitianEig const Int n = A.Height(); Int numValid = n; if( subset.indexSubset ) - { + { numValid = subset.upperIndex-subset.lowerIndex+1; } else if( subset.rangeSubset ) { if( subset.lowerBound >= Real(0) || subset.upperBound < Real(0) ) - { + { numValid = 0; } } @@ -701,9 +701,9 @@ template HermitianEigInfo SequentialHelper ( UpperOrLower uplo, - AbstractDistMatrix& A, + AbstractDistMatrix& A, AbstractDistMatrix>& w, - AbstractDistMatrix& Q, + AbstractDistMatrix& Q, const HermitianEigCtrl& ctrl ) { EL_DEBUG_CSE @@ -728,7 +728,7 @@ SequentialHelper Matrix> wProx; Matrix QProx; wProx.Resize( n, 1 ); - QProx.Resize( n, n ); + QProx.Resize( n, n ); info = HermitianEig( uplo, A.Matrix(), wProx, QProx, ctrl ); @@ -890,18 +890,18 @@ MRRR else kEst = n; - // We will use the same buffer for Q in the vector distribution used by - // PMRRR as for the matrix distribution used by Elemental. In order to + // We will use the same buffer for Q in the vector distribution used by + // PMRRR as for the matrix distribution used by Elemental. In order to // do so, we must pad Q's dimensions slightly. const Int N = MaxLength(n,g.Height())*g.Height(); - const Int K = MaxLength(kEst,g.Size())*g.Size(); + const Int K = MaxLength(kEst,g.Size())*g.Size(); ElementalProxyCtrl proxCtrl; proxCtrl.colConstrain = true; proxCtrl.rowConstrain = true; proxCtrl.colAlign = 0; proxCtrl.rowAlign = 0; - + DistMatrixWriteProxy QProx( QPre, proxCtrl ); auto& Q = QProx.Get(); @@ -941,14 +941,14 @@ MRRR const Int k = w.Height(); { - // Redistribute Q piece-by-piece in place. This is to keep the + // Redistribute Q piece-by-piece in place. This is to keep the // send/recv buffer memory usage low. const Int p = g.Size(); const Int numEqualPanels = K/p; const Int numPanelsPerComm = (numEqualPanels / TARGET_CHUNKS) + 1; const Int nbProp = numPanelsPerComm*p; - // Manually maintain information about the implicit Q[* ,VR] stored + // Manually maintain information about the implicit Q[* ,VR] stored // at the end of the Q[MC,MR] buffers. Int alignment = 0; const Real* readBuffer = Q_STAR_VR.LockedBuffer(); diff --git a/src/lapack_like/spectral/HermitianTridiagEig.cpp b/src/lapack_like/spectral/HermitianTridiagEig.cpp index dc8600b370..97d73d1cc4 100644 --- a/src/lapack_like/spectral/HermitianTridiagEig.cpp +++ b/src/lapack_like/spectral/HermitianTridiagEig.cpp @@ -1,9 +1,9 @@ /* - Copyright (c) 2009-2016, Jack Poulson + Copyright (c) 2009-2017, Jack Poulson All rights reserved. - This file is part of Elemental and is under the BSD 2-Clause License, - which can be found in the LICENSE file in the root directory, or at + This file is part of Elemental and is under the BSD 2-Clause License, + which can be found in the LICENSE file in the root directory, or at http://opensource.org/licenses/BSD-2-Clause */ #include @@ -315,7 +315,7 @@ QRHelper { EL_DEBUG_CSE HermitianTridiagEigInfo info; - w = d; + w = d; info.qrInfo = QRAlg( w, dSub, ctrl ); herm_eig::SortAndFilter( w, ctrl ); return info; @@ -332,6 +332,8 @@ DCHelper EL_DEBUG_CSE HermitianTridiagEigInfo info; auto ctrlMod( ctrl ); + ctrlMod.subset.indexSubset = false; + ctrlMod.subset.rangeSubset = false; ctrlMod.wantEigVecs = false; Matrix Q; info.dcInfo = DivideAndConquer( d, dSub, w, Q, ctrlMod ); @@ -339,7 +341,8 @@ DCHelper return info; } -template>> +template>> HermitianTridiagEigInfo LAPACKHelper ( Matrix& d, @@ -355,7 +358,7 @@ LAPACKHelper if( ctrl.subset.rangeSubset ) { const Int k = lapack::SymmetricTridiagEig - ( BlasInt(n), d.Buffer(), dSub.Buffer(), w.Buffer(), + ( BlasInt(n), d.Buffer(), dSub.Buffer(), w.Buffer(), ctrl.subset.lowerBound, ctrl.subset.upperBound ); w.Resize( k, 1 ); } @@ -363,7 +366,7 @@ LAPACKHelper { const Int numEig = ctrl.subset.upperIndex-ctrl.subset.lowerIndex+1; lapack::SymmetricTridiagEig - ( BlasInt(n), d.Buffer(), dSub.Buffer(), w.Buffer(), + ( BlasInt(n), d.Buffer(), dSub.Buffer(), w.Buffer(), BlasInt(ctrl.subset.lowerIndex), BlasInt(ctrl.subset.upperIndex) ); w.Resize( numEig, 1 ); @@ -376,7 +379,8 @@ LAPACKHelper return info; } -template>> +template>> HermitianTridiagEigInfo Helper ( const Matrix& d, @@ -403,7 +407,9 @@ Helper return LAPACKHelper( dMod, dSubMod, w, ctrl ); } -template>,typename=void> +template>, + typename=void> HermitianTridiagEigInfo Helper ( const Matrix& d, @@ -462,7 +468,7 @@ HermitianTridiagEigInfo QRHelper ( const AbstractDistMatrix& d, const AbstractDistMatrix& dSub, - AbstractDistMatrix& w, + AbstractDistMatrix& w, const HermitianTridiagEigCtrl& ctrl ) { EL_DEBUG_CSE @@ -515,6 +521,8 @@ DCHelper auto& w = wProx.Get(); auto ctrlMod( ctrl ); + ctrlMod.subset.indexSubset = false; + ctrlMod.subset.rangeSubset = false; ctrlMod.wantEigVecs = false; DistMatrix Q(w.Grid()); info.dcInfo = @@ -546,6 +554,8 @@ DCHelper auto& w = wProx.Get(); auto ctrlMod( ctrl ); + ctrlMod.subset.indexSubset = false; + ctrlMod.subset.rangeSubset = false; ctrlMod.wantEigVecs = false; DistMatrix Q(w.Grid()); info.dcInfo = @@ -555,7 +565,8 @@ DCHelper return info; } -template>> +template>> HermitianTridiagEigInfo MRRRHelper ( const AbstractDistMatrix& d, @@ -586,17 +597,17 @@ MRRRHelper herm_tridiag_eig::Info rangeInfo; if( ctrl.subset.rangeSubset ) rangeInfo = herm_tridiag_eig::Eig - ( int(n), d_STAR_STAR.Buffer(), dSub_STAR_STAR.Buffer(), - wVector.data(), w.ColComm(), + ( int(n), d_STAR_STAR.Buffer(), dSub_STAR_STAR.Buffer(), + wVector.data(), w.ColComm(), ctrl.subset.lowerBound, ctrl.subset.upperBound ); else if( ctrl.subset.indexSubset ) rangeInfo = herm_tridiag_eig::Eig - ( int(n), d_STAR_STAR.Buffer(), dSub_STAR_STAR.Buffer(), - wVector.data(), w.ColComm(), + ( int(n), d_STAR_STAR.Buffer(), dSub_STAR_STAR.Buffer(), + wVector.data(), w.ColComm(), int(ctrl.subset.lowerIndex), int(ctrl.subset.upperIndex) ); else rangeInfo = herm_tridiag_eig::Eig - ( int(n), d_STAR_STAR.Buffer(), dSub_STAR_STAR.Buffer(), + ( int(n), d_STAR_STAR.Buffer(), dSub_STAR_STAR.Buffer(), wVector.data(), w.ColComm() ); w.Resize( rangeInfo.numGlobalEigenvalues, 1 ); for( Int iLoc=0; iLoc>> +template>> HermitianTridiagEigInfo Helper ( const AbstractDistMatrix& d, @@ -629,12 +641,14 @@ Helper } } -template>,typename=void> +template>, + typename=void> HermitianTridiagEigInfo Helper ( const AbstractDistMatrix& d, const AbstractDistMatrix& dSub, - AbstractDistMatrix& w, + AbstractDistMatrix& w, const HermitianTridiagEigCtrl& ctrl ) { EL_DEBUG_CSE @@ -648,12 +662,13 @@ Helper } } -template>> +template>> HermitianTridiagEigInfo MRRRHelper ( const AbstractDistMatrix& d, const AbstractDistMatrix>& dSub, - AbstractDistMatrix& wPre, + AbstractDistMatrix& wPre, const HermitianTridiagEigCtrl& ctrl ) { EL_DEBUG_CSE @@ -711,12 +726,13 @@ MRRRHelper return info; } -template>> +template>> HermitianTridiagEigInfo Helper ( const AbstractDistMatrix& d, const AbstractDistMatrix>& dSub, - AbstractDistMatrix& wPre, + AbstractDistMatrix& wPre, const HermitianTridiagEigCtrl& ctrl ) { EL_DEBUG_CSE @@ -734,12 +750,14 @@ Helper } } -template>,typename=void> +template>, + typename=void> HermitianTridiagEigInfo Helper ( const AbstractDistMatrix& d, const AbstractDistMatrix>& dSub, - AbstractDistMatrix& w, + AbstractDistMatrix& w, const HermitianTridiagEigCtrl& ctrl ) { EL_DEBUG_CSE @@ -760,7 +778,7 @@ HermitianTridiagEigInfo HermitianTridiagEig ( const AbstractDistMatrix>& d, const AbstractDistMatrix& dSub, - AbstractDistMatrix>& w, + AbstractDistMatrix>& w, const HermitianTridiagEigCtrl>& ctrl ) { EL_DEBUG_CSE @@ -783,7 +801,7 @@ QRHelper { EL_DEBUG_CSE HermitianTridiagEigInfo info; - w = d; + w = d; info.qrInfo = QRAlg( w, dSub, Q, ctrl ); herm_eig::SortAndFilter( w, Q, ctrl ); return info; @@ -806,13 +824,17 @@ DCHelper } else { - info.dcInfo = DivideAndConquer( d, dSub, w, Q, ctrl ); + auto ctrlMod( ctrl ); + ctrlMod.subset.indexSubset = false; + ctrlMod.subset.rangeSubset = false; + info.dcInfo = DivideAndConquer( d, dSub, w, Q, ctrlMod ); herm_eig::SortAndFilter( w, Q, ctrl ); } return info; } -template>> +template>> HermitianTridiagEigInfo LAPACKHelper ( Matrix& d, @@ -831,7 +853,7 @@ LAPACKHelper { Q.Resize( n, n ); const Int k = lapack::SymmetricTridiagEig - ( BlasInt(n), d.Buffer(), dSub.Buffer(), w.Buffer(), + ( BlasInt(n), d.Buffer(), dSub.Buffer(), w.Buffer(), Q.Buffer(), BlasInt(Q.LDim()), ctrl.subset.lowerBound, ctrl.subset.upperBound ); w.Resize( k, 1 ); @@ -842,7 +864,7 @@ LAPACKHelper const Int numEig = ctrl.subset.upperIndex-ctrl.subset.lowerIndex+1; Q.Resize( n, numEig ); lapack::SymmetricTridiagEig - ( BlasInt(n), d.Buffer(), dSub.Buffer(), w.Buffer(), + ( BlasInt(n), d.Buffer(), dSub.Buffer(), w.Buffer(), Q.Buffer(), BlasInt(Q.LDim()), BlasInt(ctrl.subset.lowerIndex), BlasInt(ctrl.subset.upperIndex) ); @@ -852,7 +874,7 @@ LAPACKHelper { Q.Resize( n, n ); lapack::SymmetricTridiagEig - ( BlasInt(n), d.Buffer(), dSub.Buffer(), w.Buffer(), + ( BlasInt(n), d.Buffer(), dSub.Buffer(), w.Buffer(), Q.Buffer(), BlasInt(Q.LDim()) ); } auto sortPairs = TaggedSort( w, ctrl.sort ); @@ -863,7 +885,8 @@ LAPACKHelper return info; } -template>> +template>> HermitianTridiagEigInfo Helper ( const Matrix& d, @@ -894,7 +917,9 @@ Helper } } -template>,typename=void> +template>, + typename=void> HermitianTridiagEigInfo Helper ( const Matrix& d, @@ -925,7 +950,7 @@ HermitianTridiagEigInfo Helper ( const Matrix& d, const Matrix>& dSub, - Matrix& w, + Matrix& w, Matrix>& Q, const HermitianTridiagEigCtrl& ctrl ) { @@ -953,7 +978,7 @@ HermitianTridiagEig ( const Matrix>& d, const Matrix& dSub, Matrix>& w, - Matrix& Q, + Matrix& Q, const HermitianTridiagEigCtrl>& ctrl ) { EL_DEBUG_CSE @@ -967,8 +992,8 @@ HermitianTridiagEigInfo QRHelper ( const AbstractDistMatrix& d, const AbstractDistMatrix& dSub, - AbstractDistMatrix& w, - AbstractDistMatrix& QPre, + AbstractDistMatrix& w, + AbstractDistMatrix& QPre, const HermitianTridiagEigCtrl& ctrl ) { EL_DEBUG_CSE @@ -1000,8 +1025,8 @@ HermitianTridiagEigInfo QRHelper ( const AbstractDistMatrix& d, const AbstractDistMatrix>& dSub, - AbstractDistMatrix& w, - AbstractDistMatrix>& QPre, + AbstractDistMatrix& w, + AbstractDistMatrix>& QPre, const HermitianTridiagEigCtrl& ctrl ) { EL_DEBUG_CSE @@ -1048,8 +1073,8 @@ HermitianTridiagEigInfo DCHelper ( const AbstractDistMatrix& d, const AbstractDistMatrix& dSub, - AbstractDistMatrix& wPre, - AbstractDistMatrix& QPre, + AbstractDistMatrix& wPre, + AbstractDistMatrix& QPre, const HermitianTridiagEigCtrl& ctrl ) { EL_DEBUG_CSE @@ -1066,9 +1091,12 @@ DCHelper auto& w = wProx.Get(); DistMatrixWriteProxy QProx( QPre ); auto& Q = QProx.Get(); + auto ctrlMod( ctrl ); + ctrlMod.subset.indexSubset = false; + ctrlMod.subset.rangeSubset = false; info.dcInfo = DivideAndConquer - ( d_STAR_STAR.Matrix(), dSub_STAR_STAR.Matrix(), w, Q, ctrl ); + ( d_STAR_STAR.Matrix(), dSub_STAR_STAR.Matrix(), w, Q, ctrlMod ); herm_eig::SortAndFilter( w, Q, ctrl ); } @@ -1107,9 +1135,12 @@ DCHelper DistMatrixWriteProxy wProx( wPre ); auto& w = wProx.Get(); + auto ctrlMod( ctrl ); + ctrlMod.subset.indexSubset = false; + ctrlMod.subset.rangeSubset = false; info.dcInfo = DivideAndConquer - ( d_STAR_STAR.Matrix(), dSubReal.Matrix(), w, QReal, ctrl ); + ( d_STAR_STAR.Matrix(), dSubReal.Matrix(), w, QReal, ctrlMod ); herm_eig::SortAndFilter( w, QReal, ctrl ); Copy( QReal, Q ); @@ -1119,13 +1150,14 @@ DCHelper return info; } -template>> +template>> HermitianTridiagEigInfo MRRRHelper ( const AbstractDistMatrix& d, const AbstractDistMatrix& dSub, - AbstractDistMatrix& wPre, - AbstractDistMatrix& QPre, + AbstractDistMatrix& wPre, + AbstractDistMatrix& QPre, const HermitianTridiagEigCtrl& ctrl ) { EL_DEBUG_CSE @@ -1174,17 +1206,17 @@ MRRRHelper vector wVector(n); if( ctrl.subset.rangeSubset ) rangeInfo = herm_tridiag_eig::Eig - ( int(n), d_STAR_STAR.Buffer(), dSub_STAR_STAR.Buffer(), + ( int(n), d_STAR_STAR.Buffer(), dSub_STAR_STAR.Buffer(), wVector.data(), Q.Buffer(), Q.LDim(), w.ColComm(), ctrl.subset.lowerBound, ctrl.subset.upperBound ); else if( ctrl.subset.indexSubset ) rangeInfo = herm_tridiag_eig::Eig - ( int(n), d_STAR_STAR.Buffer(), dSub_STAR_STAR.Buffer(), + ( int(n), d_STAR_STAR.Buffer(), dSub_STAR_STAR.Buffer(), wVector.data(), Q.Buffer(), Q.LDim(), w.ColComm(), int(ctrl.subset.lowerIndex), int(ctrl.subset.upperIndex) ); else rangeInfo = herm_tridiag_eig::Eig - ( int(n), d_STAR_STAR.Buffer(), dSub_STAR_STAR.Buffer(), + ( int(n), d_STAR_STAR.Buffer(), dSub_STAR_STAR.Buffer(), wVector.data(), Q.Buffer(), Q.LDim(), w.ColComm() ); w.Resize( rangeInfo.numGlobalEigenvalues, 1 ); Q.Resize( n, rangeInfo.numGlobalEigenvalues ); @@ -1199,13 +1231,14 @@ MRRRHelper return info; } -template>> +template>> HermitianTridiagEigInfo MRRRHelper ( const AbstractDistMatrix& d, const AbstractDistMatrix>& dSub, - AbstractDistMatrix& wPre, - AbstractDistMatrix>& QPre, + AbstractDistMatrix& wPre, + AbstractDistMatrix>& QPre, const HermitianTridiagEigCtrl& ctrl ) { EL_DEBUG_CSE @@ -1261,17 +1294,17 @@ MRRRHelper vector wVector(n); if( ctrl.subset.rangeSubset ) rangeInfo = herm_tridiag_eig::Eig - ( int(n), d_STAR_STAR.Buffer(), dSubReal.Buffer(), + ( int(n), d_STAR_STAR.Buffer(), dSubReal.Buffer(), wVector.data(), QReal.Buffer(), QReal.LDim(), w.ColComm(), ctrl.subset.lowerBound, ctrl.subset.upperBound ); else if( ctrl.subset.indexSubset ) rangeInfo = herm_tridiag_eig::Eig - ( int(n), d_STAR_STAR.Buffer(), dSubReal.Buffer(), + ( int(n), d_STAR_STAR.Buffer(), dSubReal.Buffer(), wVector.data(), QReal.Buffer(), QReal.LDim(), w.ColComm(), int(ctrl.subset.lowerIndex), int(ctrl.subset.upperIndex) ); else rangeInfo = herm_tridiag_eig::Eig - ( int(n), d_STAR_STAR.Buffer(), dSubReal.Buffer(), + ( int(n), d_STAR_STAR.Buffer(), dSubReal.Buffer(), wVector.data(), QReal.Buffer(), QReal.LDim(), w.ColComm() ); w.Resize( rangeInfo.numGlobalEigenvalues, 1 ); @@ -1287,18 +1320,19 @@ MRRRHelper ApplyTaggedSortToEachRow( sortPairs, QReal ); Copy( QReal, Q ); - DiagonalScale( LEFT, NORMAL, phase, Q ); + DiagonalScale( LEFT, NORMAL, phase, Q ); return info; } -template>> +template>> HermitianTridiagEigInfo Helper ( const AbstractDistMatrix& d, const AbstractDistMatrix& dSub, - AbstractDistMatrix& w, - AbstractDistMatrix& Q, + AbstractDistMatrix& w, + AbstractDistMatrix& Q, const HermitianTridiagEigCtrl& ctrl ) { EL_DEBUG_CSE @@ -1316,13 +1350,15 @@ Helper } } -template>,typename=void> +template>, + typename=void> HermitianTridiagEigInfo Helper ( const AbstractDistMatrix& d, const AbstractDistMatrix& dSub, - AbstractDistMatrix& w, - AbstractDistMatrix& Q, + AbstractDistMatrix& w, + AbstractDistMatrix& Q, const HermitianTridiagEigCtrl& ctrl ) { EL_DEBUG_CSE @@ -1336,13 +1372,14 @@ Helper } } -template>> +template>> HermitianTridiagEigInfo Helper ( const AbstractDistMatrix& d, const AbstractDistMatrix>& dSub, - AbstractDistMatrix& w, - AbstractDistMatrix>& Q, + AbstractDistMatrix& w, + AbstractDistMatrix>& Q, const HermitianTridiagEigCtrl& ctrl ) { EL_DEBUG_CSE @@ -1360,13 +1397,15 @@ Helper } } -template>,typename=void> +template>, + typename=void> HermitianTridiagEigInfo Helper ( const AbstractDistMatrix& d, const AbstractDistMatrix>& dSub, - AbstractDistMatrix& w, - AbstractDistMatrix>& QPre, + AbstractDistMatrix& w, + AbstractDistMatrix>& QPre, const HermitianTridiagEigCtrl& ctrl ) { EL_DEBUG_CSE @@ -1388,7 +1427,7 @@ HermitianTridiagEig ( const AbstractDistMatrix>& d, const AbstractDistMatrix& dSub, AbstractDistMatrix>& w, - AbstractDistMatrix& Q, + AbstractDistMatrix& Q, const HermitianTridiagEigCtrl>& ctrl ) { EL_DEBUG_CSE @@ -1397,7 +1436,8 @@ HermitianTridiagEig namespace herm_tridiag_eig { -template>> +template>> Int MRRREstimateHelper ( const AbstractDistMatrix& d, const AbstractDistMatrix& dSub, @@ -1421,7 +1461,9 @@ Int MRRREstimateHelper return estimate.numGlobalEigenvalues; } -template>,typename=void> +template>, + typename=void> Int MRRREstimateHelper ( const AbstractDistMatrix& d, const AbstractDistMatrix& dSub, @@ -1435,7 +1477,8 @@ Int MRRREstimateHelper } // Q is assumed to be sufficiently large and properly aligned -template>> +template>> HermitianTridiagEigInfo MRRRPostEstimateHelper ( const AbstractDistMatrix& d, @@ -1488,7 +1531,9 @@ MRRRPostEstimateHelper return info; } -template>,typename=void> +template>, + typename=void> HermitianTridiagEigInfo MRRRPostEstimateHelper ( const AbstractDistMatrix& d, diff --git a/src/lapack_like/spectral/HermitianTridiagEig/DivideAndConquer.hpp b/src/lapack_like/spectral/HermitianTridiagEig/DivideAndConquer.hpp index c729724852..a4797d12be 100644 --- a/src/lapack_like/spectral/HermitianTridiagEig/DivideAndConquer.hpp +++ b/src/lapack_like/spectral/HermitianTridiagEig/DivideAndConquer.hpp @@ -1,5 +1,5 @@ /* - Copyright (c) 2009-2016, Jack Poulson + Copyright (c) 2009-2017, Jack Poulson All rights reserved. This file is part of Elemental and is under the BSD 2-Clause License, @@ -1126,10 +1126,13 @@ DivideAndConquer EL_DEBUG_CSE const Int n = mainDiag.Height(); const auto& dcCtrl = ctrl.dcCtrl; + if( ctrl.subset.indexSubset || ctrl.subset.rangeSubset ) + LogicError + ("DivideAndConquer should not have been called directly for subset " + "computation"); DCInfo info; auto& secularInfo = info.secularInfo; - if( n <= Max(dcCtrl.cutoff,3) ) { auto ctrlMod( ctrl ); @@ -1230,11 +1233,15 @@ DivideAndConquer bool topLevel=true ) { EL_DEBUG_CSE - const Grid& grid = Q.Grid(); const Int n = mainDiag.Height(); const auto& dcCtrl = ctrl.dcCtrl; - DCInfo info; + const Grid& grid = Q.Grid(); + if( ctrl.subset.indexSubset || ctrl.subset.rangeSubset ) + LogicError + ("DivideAndConquer should not have been called directly for subset " + "computation"); + DCInfo info; if( n <= Max(dcCtrl.cutoff,3) ) { // Run the problem redundantly locally @@ -1327,8 +1334,7 @@ DivideAndConquer // TODO(poulson): A more intelligent split point. const Int split = (n/2) + 1; const Grid *leftGrid, *rightGrid; - const bool splitGrid = - SplitGrid( split, n-split, grid, leftGrid, rightGrid ); + SplitGrid( split, n-split, grid, leftGrid, rightGrid ); const Real& beta = superDiag(split-1); auto mainDiag0 = mainDiag( IR(0,split), ALL ); @@ -1425,12 +1431,6 @@ DivideAndConquer secularInfo.numSmallUpdateDeflations += info1.secularInfo.numSmallUpdateDeflations; } - if( splitGrid ) - { - // TODO(poulson): Introduce a smart pointer instead? - delete leftGrid; - delete rightGrid; - } if( topLevel ) {