Skip to content

Commit

Permalink
Pulling new Hessenberg QR algorith implementations into the main libr…
Browse files Browse the repository at this point in the history
…ary as HessenbergSchur
  • Loading branch information
Jack Poulson committed Jun 25, 2016
1 parent cd61fb5 commit a95d621
Show file tree
Hide file tree
Showing 16 changed files with 464 additions and 326 deletions.
9 changes: 5 additions & 4 deletions README.md
Expand Up @@ -58,7 +58,8 @@ datatypes:
- Tikhonov (and ridge) regression
- Equality-constrained Least Squares
- General (Gauss-Markov) Linear Models
* High-performance pseudospectral computation and visualization
* (Arbitrary-precision) High-performance pseudospectral computation and visualization
* (Arbitrary-precision) Aggressive Early Deflation Schur decompositions (currently sequential only)
* (Arbitrary-precision) blocked column-pivoted QR via Johnson-Lindenstrauss
* (Arbitrary-precision) quadratic-time low-rank Cholesky and LU modifications
* (Arbitrary-precision) Bunch-Kaufman and Bunch-Parlett for accurate symmetric
Expand All @@ -75,10 +76,10 @@ datatypes:
* Sign-based Lyapunov/Ricatti/Sylvester solvers

**Lattice reduction**:
* An extension of [Householder-based LLL](http://perso.ens-lyon.fr/damien.stehle/HLLL.html) to real and complex linearly-dependent bases
* Generalizations of [BKZ 2.0](http://link.springer.com/chapter/10.1007%2F978-3-642-25385-0_1) to complex bases
* An extension of [Householder-based LLL](http://perso.ens-lyon.fr/damien.stehle/HLLL.html) to real and complex linearly-dependent bases (currently sequential only)
* Generalizations of [BKZ 2.0](http://link.springer.com/chapter/10.1007%2F978-3-642-25385-0_1) to complex bases (currently sequential only)
incorporating ["y-sparse" enumeration](https://eprint.iacr.org/2014/980)
* Integer images/kernels and relation-finding
* Integer images/kernels and relation-finding (currently sequential only)

### The current development roadmap

Expand Down
13 changes: 13 additions & 0 deletions include/El/lapack_like/condense.hpp
Expand Up @@ -170,6 +170,19 @@ void ApplyQ
const ElementalMatrix<F>& phase,
ElementalMatrix<F>& B );

template<typename F>
void FormQ
( UpperOrLower uplo,
const Matrix<F>& A,
const Matrix<F>& phase,
Matrix<F>& Q );
template<typename F>
void FormQ
( UpperOrLower uplo,
const ElementalMatrix<F>& A,
const ElementalMatrix<F>& phase,
ElementalMatrix<F>& Q );

} // namespace hessenberg

} // namespace El
Expand Down
120 changes: 120 additions & 0 deletions include/El/lapack_like/spectral.hpp
Expand Up @@ -302,6 +302,126 @@ PolarInfo HermitianPolar
ElementalMatrix<F>& P,
const PolarCtrl& ctrl=PolarCtrl() );

// Hessenberg Schur decomposition
// ==============================
struct HessenbergSchurInfo
{
Int numUnconverged=0;
Int numIterations=0;
};

namespace hess_schur {
namespace aed {

// Cf. LAPACK's IPARMQ for these choices. The primary difference here is that
// we do not use a fixed value (of 256) for windows of size at least 6000.
inline Int NumShifts( Int n, Int winSize )
{
Int numShifts;
if( winSize < 30 )
numShifts = 2;
else if( winSize < 60 )
numShifts = 4;
else if( winSize < 150 )
numShifts = 10;
else if( winSize < 590 )
numShifts = Max( 10, winSize/Int(Log2(double(winSize))) );
else if( winSize < 3000 )
numShifts = 64;
else if( winSize < 6000 )
numShifts = 128;
else
numShifts = Max( 256, winSize/Int(2*Log2(double(winSize))) );

numShifts = Min( numShifts, winSize );
numShifts = Max( 2, numShifts-Mod(numShifts,2) );

return numShifts;
}

// Cf. LAPACK's IPARMQ for these deflation window sizes
inline Int DeflationSize( Int n, Int winSize, Int numShifts )
{
Int deflationSize;
if( winSize <= 500 )
deflationSize = numShifts;
else
deflationSize = (3*numShifts) / 2;

deflationSize = Min( deflationSize, winSize );
deflationSize = Min( deflationSize, (n-1)/3 );
deflationSize = Max( 2, deflationSize-Mod(deflationSize,2) );

return deflationSize;
}

// Cf. LAPACK's IPARMQ for the choice of skipping a QR sweep if at least
// 14% of the eigenvalues in a window deflated
inline Int SufficientDeflation( Int deflationSize )
{
const Int nibble = 14;
return (nibble*deflationSize) / 100;
}

} // namespace aed
} // namespace hess_schur

struct HessenbergSchurCtrl
{
Int winBeg=0;
Int winEnd=END;
bool fullTriangle=true;
bool wantSchurVecs=false;
bool demandConverged=true;

bool useAED=true;
bool recursiveAED=true;
bool accumulateReflections=true;

bool progress=false;

// Cf. LAPACK's IPARMQ for this choice;
// note that LAPACK's hard minimum of 12 does not apply to us
Int minAEDSize = 75;

function<Int(Int,Int)> numShifts =
function<Int(Int,Int)>(hess_schur::aed::NumShifts);

function<Int(Int,Int,Int)> deflationSize =
function<Int(Int,Int,Int)>(hess_schur::aed::DeflationSize);

function<Int(Int)> sufficientDeflation =
function<Int(Int)>(hess_schur::aed::SufficientDeflation);
};

template<typename Real>
HessenbergSchurInfo
HessenbergSchur
( Matrix<Real>& H,
Matrix<Complex<Real>>& w,
const HessenbergSchurCtrl& ctrl=HessenbergSchurCtrl() );
template<typename Real>
HessenbergSchurInfo
HessenbergSchur
( Matrix<Real>& H,
Matrix<Complex<Real>>& w,
Matrix<Real>& Z,
const HessenbergSchurCtrl& ctrl=HessenbergSchurCtrl() );

template<typename Real>
HessenbergSchurInfo
HessenbergSchur
( Matrix<Complex<Real>>& H,
Matrix<Complex<Real>>& w,
const HessenbergSchurCtrl& ctrl=HessenbergSchurCtrl() );
template<typename Real>
HessenbergSchurInfo
HessenbergSchur
( Matrix<Complex<Real>>& H,
Matrix<Complex<Real>>& w,
Matrix<Complex<Real>>& Z,
const HessenbergSchurCtrl& ctrl=HessenbergSchurCtrl() );

// Schur decomposition
// ===================
// Forward declaration
Expand Down
13 changes: 12 additions & 1 deletion src/lapack_like/condense/Hessenberg.cpp
Expand Up @@ -11,6 +11,7 @@
#include "./Hessenberg/L.hpp"
#include "./Hessenberg/U.hpp"
#include "./Hessenberg/ApplyQ.hpp"
#include "./Hessenberg/FormQ.hpp"

namespace El {

Expand Down Expand Up @@ -81,7 +82,17 @@ void ExplicitCondensed( UpperOrLower uplo, ElementalMatrix<F>& A )
( LeftOrRight side, UpperOrLower uplo, Orientation orientation, \
const ElementalMatrix<F>& A, \
const ElementalMatrix<F>& phase, \
ElementalMatrix<F>& B );
ElementalMatrix<F>& B ); \
template void hessenberg::FormQ \
( UpperOrLower uplo, \
const Matrix<F>& A, \
const Matrix<F>& phase, \
Matrix<F>& Q ); \
template void hessenberg::FormQ \
( UpperOrLower uplo, \
const ElementalMatrix<F>& A, \
const ElementalMatrix<F>& phase, \
ElementalMatrix<F>& Q );

#define EL_NO_INT_PROTO
#define EL_ENABLE_DOUBLEDOUBLE
Expand Down
46 changes: 46 additions & 0 deletions src/lapack_like/condense/Hessenberg/FormQ.hpp
@@ -0,0 +1,46 @@
/*
Copyright (c) 2009-2016, 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
http://opensource.org/licenses/BSD-2-Clause
*/
#ifndef EL_HESSENBERG_FORMQ_HPP
#define EL_HESSENBERG_FORMQ_HPP

namespace El {
namespace hessenberg {

template<typename F>
void FormQ
( UpperOrLower uplo,
const Matrix<F>& A,
const Matrix<F>& phase,
Matrix<F>& Q )
{
DEBUG_CSE
// TODO: Make this smarter
const Int n = A.Height();
Identity( Q, n, n );
ApplyQ( LEFT, uplo, NORMAL, A, phase, Q );
}

template<typename F>
void FormQ
( UpperOrLower uplo,
const ElementalMatrix<F>& A,
const ElementalMatrix<F>& phase,
ElementalMatrix<F>& Q )
{
DEBUG_CSE
// TODO: Make this smarter
const Int n = A.Height();
Identity( Q, n, n );
ApplyQ( LEFT, uplo, NORMAL, A, phase, Q );
}

} // namespace hessenberg
} // namespace El

#endif // ifndef EL_HESSENBERG_FORMQ_HPP
131 changes: 131 additions & 0 deletions src/lapack_like/spectral/HessenbergSchur.cpp
@@ -0,0 +1,131 @@
/*
Copyright (c) 2009-2016, 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
http://opensource.org/licenses/BSD-2-Clause
*/
#include <El.hpp>

#include "./HessenbergSchur/SingleShift.hpp"
#include "./HessenbergSchur/DoubleShift.hpp"
#include "./HessenbergSchur/AED.hpp"

namespace El {

template<typename Real>
HessenbergSchurInfo
HessenbergSchur
( Matrix<Real>& H,
Matrix<Complex<Real>>& w,
const HessenbergSchurCtrl& ctrl )
{
const Int n = H.Height();
auto ctrlMod( ctrl );
ctrlMod.winBeg = ( ctrl.winBeg==END ? n : ctrl.winBeg );
ctrlMod.winEnd = ( ctrl.winEnd==END ? n : ctrl.winEnd );
ctrlMod.wantSchurVecs = false;

Matrix<Real> Z;
if( ctrl.useAED )
{
return hess_schur::AED( H, w, Z, ctrlMod );
}
else
{
return hess_schur::DoubleShift( H, w, Z, ctrlMod );
}
}

template<typename Real>
HessenbergSchurInfo
HessenbergSchur
( Matrix<Real>& H,
Matrix<Complex<Real>>& w,
Matrix<Real>& Z,
const HessenbergSchurCtrl& ctrl )
{
const Int n = H.Height();
auto ctrlMod( ctrl );
ctrlMod.winBeg = ( ctrl.winBeg==END ? n : ctrl.winBeg );
ctrlMod.winEnd = ( ctrl.winEnd==END ? n : ctrl.winEnd );
ctrlMod.wantSchurVecs = true;

if( ctrl.useAED )
{
return hess_schur::AED( H, w, Z, ctrlMod );
}
else
{
return hess_schur::DoubleShift( H, w, Z, ctrlMod );
}
}

template<typename Real>
HessenbergSchurInfo
HessenbergSchur
( Matrix<Complex<Real>>& H,
Matrix<Complex<Real>>& w,
const HessenbergSchurCtrl& ctrl )
{
const Int n = H.Height();
auto ctrlMod( ctrl );
ctrlMod.winBeg = ( ctrl.winBeg==END ? n : ctrl.winBeg );
ctrlMod.winEnd = ( ctrl.winEnd==END ? n : ctrl.winEnd );
ctrlMod.wantSchurVecs = false;

Matrix<Complex<Real>> Z;
if( ctrl.useAED )
{
return hess_schur::AED( H, w, Z, ctrlMod );
}
else
{
return hess_schur::SingleShift( H, w, Z, ctrlMod );
}
}

template<typename Real>
HessenbergSchurInfo
HessenbergSchur
( Matrix<Complex<Real>>& H,
Matrix<Complex<Real>>& w,
Matrix<Complex<Real>>& Z,
const HessenbergSchurCtrl& ctrl )
{
const Int n = H.Height();
auto ctrlMod( ctrl );
ctrlMod.winBeg = ( ctrl.winBeg==END ? n : ctrl.winBeg );
ctrlMod.winEnd = ( ctrl.winEnd==END ? n : ctrl.winEnd );
ctrlMod.wantSchurVecs = true;

if( ctrl.useAED )
{
return hess_schur::AED( H, w, Z, ctrlMod );
}
else
{
return hess_schur::SingleShift( H, w, Z, ctrlMod );
}
}

#define PROTO(F) \
template HessenbergSchurInfo HessenbergSchur \
( Matrix<F>& H, \
Matrix<Complex<Base<F>>>& w, \
const HessenbergSchurCtrl& ctrl ); \
template HessenbergSchurInfo HessenbergSchur \
( Matrix<F>& H, \
Matrix<Complex<Base<F>>>& w, \
Matrix<F>& Z, \
const HessenbergSchurCtrl& ctrl );

#define EL_NO_INT_PROTO
#define EL_ENABLE_DOUBLEDOUBLE
#define EL_ENABLE_QUADDOUBLE
#define EL_ENABLE_QUAD
#define EL_ENABLE_BIGFLOAT
#include <El/macros/Instantiate.h>

} // namespace El

0 comments on commit a95d621

Please sign in to comment.