Adding multi-bulge Hessenberg QR iteration; uniting HessQRCtrl and HessenbergSchurCtrl; fixing a mistake in the previous, overzealous commit to fix an Intel compilation issue reported in Issue #177

Author: poulson

examples/lapack_like/ChunkedPseudospectra.cpp

@@ -241,9 +241,9 @@ main( int argc, char* argv[] )
         DistMatrix<C> QCpx(g);
         SchurCtrl<Real> ctrl;
 #ifdef EL_HAVE_SCALAPACK
-        ctrl.qrCtrl.blockHeight = nbDist;
-        ctrl.qrCtrl.blockWidth = nbDist;
-        ctrl.qrCtrl.distAED = false;
+        ctrl.hessSchurCtrl.blockHeight = nbDist;
+        ctrl.hessSchurCtrl.blockWidth = nbDist;
+        ctrl.hessSchurCtrl.distAED = false;
 #else
         ctrl.useSDC = true;
         ctrl.sdcCtrl.cutoff = cutoff;

examples/lapack_like/Pseudospectra.cpp

@@ -249,9 +249,9 @@ main( int argc, char* argv[] )
         psCtrl.basisSize = basisSize;
         psCtrl.progress = progress;
 #ifdef EL_HAVE_SCALAPACK
-        psCtrl.schurCtrl.qrCtrl.blockHeight = nbDist;
-        psCtrl.schurCtrl.qrCtrl.blockWidth = nbDist;
-        psCtrl.schurCtrl.qrCtrl.distAED = false;
+        psCtrl.schurCtrl.hessSchurCtrl.blockHeight = nbDist;
+        psCtrl.schurCtrl.hessSchurCtrl.blockWidth = nbDist;
+        psCtrl.schurCtrl.hessSchurCtrl.distAED = false;
 #else
         psCtrl.schurCtrl.sdcCtrl.cutoff = cutoff;
         psCtrl.schurCtrl.sdcCtrl.maxInnerIts = maxInnerIts;

examples/lapack_like/RealSchur.cpp

@@ -51,7 +51,7 @@ main( int argc, char* argv[] )
         DistMatrix<Complex<Real>,VR,STAR> w;
         SchurCtrl<Real> ctrl;
 #ifdef EL_HAVE_SCALAPACK
-        ctrl.qrCtrl.distAED = aed;
+        ctrl.hessSchurCtrl.distAED = aed;
         // TODO: distribution block size
 #else
         ctrl.useSDC = true;

examples/lapack_like/Schur.cpp

@@ -55,9 +55,9 @@ main( int argc, char* argv[] )
         DistMatrix<C,VR,STAR> w(g);
         SchurCtrl<Real> ctrl;
 #ifdef EL_HAVE_SCALAPACK
-        //ctrl.qrCtrl.blockHeight = nbDist;
-        //ctrl.qrCtrl.blockWidth = nbDist;
-        ctrl.qrCtrl.distAED = false;
+        //ctrl.hessSchurCtrl.blockHeight = nbDist;
+        //ctrl.hessSchurCtrl.blockWidth = nbDist;
+        ctrl.hessSchurCtrl.distAED = false;
 #else
         ctrl.useSDC = true;
         ctrl.sdcCtrl.cutoff = cutoff;

include/El/lapack_like/spectral.h

@@ -80,9 +80,9 @@ EL_EXPORT ElError ElSecularSVDCtrlDefault_d
 /* Hermitian tridiagonal eigensolvers
    ================================== */
 typedef enum {
-  EL_HERM_TRIDIAG_EIG_QR = 0,
-  EL_HERM_TRIDIAG_EIG_DC = 1,
-  EL_HERM_TRIDIAG_EIG_MRRR = 2
+  EL_HERM_TRIDIAG_EIG_QR=0,
+  EL_HERM_TRIDIAG_EIG_DC=1,
+  EL_HERM_TRIDIAG_EIG_MRRR=2
 } ElHermitianTridiagEigAlg;
 
 /* HermitianEigSubset */
@@ -587,12 +587,36 @@ EL_EXPORT ElError ElHermitianPolarDecompDist_z
 
 /* Schur decomposition
    =================== */
-/* HessQRCtrl */
+typedef enum {
+  EL_HESSENBERG_SCHUR_AED=0,
+  EL_HESSENBERG_SCHUR_MULTIBULGE=1,
+  EL_HESSENBERG_SCHUR_SIMPLE=2
+} ElHessenbergSchurAlg;
+
+/* HessenbergSchurCtrl */
 typedef struct {
+  ElInt winBeg;
+  ElInt winEnd;
+  bool fullTriangle;
+  bool wantSchurVecs;
+  bool demandConverged;
+
+  ElHessenbergSchurAlg alg;
+  bool recursiveAED;
+  bool accumulateReflections;
+
+  bool progress;
+
+  ElInt minMultiBulgeSize;
+  
+  ElInt (*numShifts)(ElInt,ElInt);
+  ElInt (*deflationSize)(ElInt,ElInt,ElInt);
+  ElInt (*sufficientDeflation)(ElInt);
+
   bool distAED;
   ElInt blockHeight, blockWidth;
-} ElHessQRCtrl;
-EL_EXPORT ElError ElHessQRCtrlDefault( ElHessQRCtrl* ctrl );
+} ElHessenbergSchurCtrl;
+EL_EXPORT ElError ElHessenbergSchurCtrlDefault( ElHessenbergSchurCtrl* ctrl );
 
 /* SDCCtrl */
 typedef struct {
@@ -620,15 +644,15 @@ EL_EXPORT ElError ElSDCCtrlDefault_d( ElSDCCtrl_d* ctrl );
 /* SchurCtrl */
 typedef struct {
   bool useSDC;
-  ElHessQRCtrl qrCtrl;
+  ElHessenbergSchurCtrl hessSchurCtrl;
   ElSDCCtrl_s sdcCtrl;
   bool time;
 } ElSchurCtrl_s;
 EL_EXPORT ElError ElSchurCtrlDefault_s( ElSchurCtrl_s* ctrl );
 
 typedef struct {
   bool useSDC;
-  ElHessQRCtrl qrCtrl;
+  ElHessenbergSchurCtrl hessSchurCtrl;
   ElSDCCtrl_d sdcCtrl;
   bool time;
 } ElSchurCtrl_d;

include/El/lapack_like/spectral.hpp

@@ -682,6 +682,12 @@ inline Int SufficientDeflation( Int deflationSize )
 } // namespace aed
 } // namespace hess_schur
 
+enum HessenbergSchurAlg {
+  HESSENBERG_SCHUR_AED=0,
+  HESSENBERG_SCHUR_MULTIBULGE=1,
+  HESSENBERG_SCHUR_SIMPLE=2
+};
+
 struct HessenbergSchurCtrl
 {
     Int winBeg=0;
@@ -690,15 +696,15 @@ struct HessenbergSchurCtrl
     bool wantSchurVecs=false;
     bool demandConverged=true;
 
-    bool useAED=true;
+    HessenbergSchurAlg alg=HESSENBERG_SCHUR_AED;
     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;
+    Int minMultiBulgeSize = 75;
 
     function<Int(Int,Int)> numShifts =
       function<Int(Int,Int)>(hess_schur::aed::NumShifts);
@@ -708,6 +714,11 @@ struct HessenbergSchurCtrl
 
     function<Int(Int)> sufficientDeflation =
       function<Int(Int)>(hess_schur::aed::SufficientDeflation);
+
+    // For the distributed Hessenberg QR algorithm
+    // TODO(poulson): Move this into a substructure?
+    bool distAED=false;
+    Int blockHeight=DefaultBlockHeight(), blockWidth=DefaultBlockWidth();
 };
 
 template<typename Real>
@@ -756,18 +767,11 @@ struct SDCCtrl
     SignCtrl<Real> signCtrl;
 };
 
-// TODO: Combine with HessenbergSchurCtrl
-struct HessQRCtrl 
-{
-    bool distAED=false;
-    Int blockHeight=DefaultBlockHeight(), blockWidth=DefaultBlockWidth();
-};
-
 template<typename Real>
 struct SchurCtrl 
 {
     bool useSDC=false;
-    HessQRCtrl qrCtrl;
+    HessenbergSchurCtrl hessSchurCtrl;
     SDCCtrl<Real> sdcCtrl;    
     bool time=false;
 };

include/El/lapack_like/spectral/CReflect.hpp

@@ -724,19 +724,75 @@ inline ElSVDCtrl_d CReflect( const SVDCtrl<double>& ctrl )
     return ctrlC;
 }
 
-/* HessQRCtrl */
-inline ElHessQRCtrl CReflect( const HessQRCtrl& ctrl )
-{
-    ElHessQRCtrl ctrlC;
+/* HessenbergSchurCtrl */
+inline ElHessenbergSchurAlg CReflect( const HessenbergSchurAlg& alg )
+{ return static_cast<ElHessenbergSchurAlg>(alg); }
+inline HessenbergSchurAlg CReflect( const ElHessenbergSchurAlg& alg )
+{ return static_cast<HessenbergSchurAlg>(alg); }
+
+inline ElHessenbergSchurCtrl CReflect( const HessenbergSchurCtrl& ctrl )
+{
+    ElHessenbergSchurCtrl ctrlC;
+    ctrlC.winBeg = ctrl.winBeg;
+    ctrlC.winEnd = ctrl.winEnd;
+    ctrlC.fullTriangle = ctrl.fullTriangle;
+    ctrlC.wantSchurVecs = ctrl.wantSchurVecs;
+    ctrlC.demandConverged = ctrl.demandConverged;
+
+    ctrlC.alg = CReflect(ctrl.alg);
+    ctrlC.recursiveAED = ctrl.recursiveAED;
+    ctrlC.accumulateReflections = ctrl.accumulateReflections;
+
+    ctrlC.progress = ctrl.progress;
+
+    ctrlC.minMultiBulgeSize = ctrl.minMultiBulgeSize;
+    auto numShiftsRes = ctrl.numShifts.target<ElInt(*)(ElInt,ElInt)>(); 
+    if( numShiftsRes )
+        ctrlC.numShifts = *numShiftsRes;
+    else
+        RuntimeError("Could not convert numShifts to C function pointer");
+
+    auto deflationSizeRes =
+      ctrl.deflationSize.target<ElInt(*)(ElInt,ElInt,ElInt)>();
+    if( deflationSizeRes )
+        ctrlC.deflationSize = *deflationSizeRes;
+    else
+        RuntimeError("Could not convert deflationSize to C function pointer");
+
+    auto sufficientDeflationRes =
+      ctrl.sufficientDeflation.target<ElInt(*)(ElInt)>();
+    if( sufficientDeflationRes )
+        ctrlC.sufficientDeflation = *sufficientDeflationRes;
+    else
+        RuntimeError
+        ("Could not convert sufficientDeflation to C function pointer");
+
     ctrlC.distAED = ctrl.distAED;
     ctrlC.blockHeight = ctrl.blockHeight;
     ctrlC.blockWidth = ctrl.blockWidth;
     return ctrlC;
 }
 
-inline HessQRCtrl CReflect( const ElHessQRCtrl& ctrlC )
+inline HessenbergSchurCtrl CReflect( const ElHessenbergSchurCtrl& ctrlC )
 {
-    HessQRCtrl ctrl;
+    HessenbergSchurCtrl ctrl;
+    ctrl.winBeg = ctrlC.winBeg;
+    ctrl.winEnd = ctrlC.winEnd;
+    ctrl.fullTriangle = ctrlC.fullTriangle;
+    ctrl.wantSchurVecs = ctrlC.wantSchurVecs;
+    ctrl.demandConverged = ctrlC.demandConverged;
+
+    ctrl.alg = CReflect(ctrlC.alg);
+    ctrl.recursiveAED = ctrlC.recursiveAED;
+    ctrl.accumulateReflections = ctrlC.accumulateReflections;
+
+    ctrl.progress = ctrlC.progress;
+
+    ctrl.minMultiBulgeSize = ctrlC.minMultiBulgeSize;
+    ctrl.numShifts = ctrlC.numShifts;
+    ctrl.deflationSize = ctrlC.deflationSize;
+    ctrl.sufficientDeflation = ctrlC.sufficientDeflation;
+
     ctrl.distAED = ctrlC.distAED;
     ctrl.blockHeight = ctrlC.blockHeight;
     ctrl.blockWidth = ctrlC.blockWidth;
@@ -799,7 +855,7 @@ inline ElSchurCtrl_s CReflect( const SchurCtrl<float>& ctrl )
 {
     ElSchurCtrl_s ctrlC;
     ctrlC.useSDC = ctrl.useSDC;
-    ctrlC.qrCtrl = CReflect( ctrl.qrCtrl );
+    ctrlC.hessSchurCtrl = CReflect( ctrl.hessSchurCtrl );
     ctrlC.sdcCtrl = CReflect( ctrl.sdcCtrl );
     ctrlC.time = ctrl.time;
     return ctrlC;
@@ -808,7 +864,7 @@ inline ElSchurCtrl_d CReflect( const SchurCtrl<double>& ctrl )
 {
     ElSchurCtrl_d ctrlC;
     ctrlC.useSDC = ctrl.useSDC;
-    ctrlC.qrCtrl = CReflect( ctrl.qrCtrl );
+    ctrlC.hessSchurCtrl = CReflect( ctrl.hessSchurCtrl );
     ctrlC.sdcCtrl = CReflect( ctrl.sdcCtrl );
     ctrlC.time = ctrl.time;
     return ctrlC;
@@ -818,7 +874,7 @@ inline SchurCtrl<float> CReflect( const ElSchurCtrl_s& ctrlC )
 {
     SchurCtrl<float> ctrl;
     ctrl.useSDC = ctrlC.useSDC;
-    ctrl.qrCtrl = CReflect( ctrlC.qrCtrl );
+    ctrl.hessSchurCtrl = CReflect( ctrlC.hessSchurCtrl );
     ctrl.sdcCtrl = CReflect( ctrlC.sdcCtrl );
     ctrl.time = ctrlC.time;
     return ctrl;
@@ -827,7 +883,7 @@ inline SchurCtrl<double> CReflect( const ElSchurCtrl_d& ctrlC )
 {
     SchurCtrl<double> ctrl;
     ctrl.useSDC = ctrlC.useSDC;
-    ctrl.qrCtrl = CReflect( ctrlC.qrCtrl );
+    ctrl.hessSchurCtrl = CReflect( ctrlC.hessSchurCtrl );
     ctrl.sdcCtrl = CReflect( ctrlC.sdcCtrl );
     ctrl.time = ctrlC.time;
     return ctrl;

python/lapack_like/spectral.py

@@ -10,6 +10,7 @@
 from ..blas_like import Copy, EntrywiseMap
 from ..io import ProcessEvents
 import ctypes
+from ctypes import CFUNCTYPE
 
 # Cubic secular
 # =============
@@ -550,12 +551,28 @@ class SignCtrl_d(ctypes.Structure):
   def __init__(self):
     lib.ElSignCtrlDefault_d(pointer(self))
 
-lib.ElHessQRCtrlDefault.argtypes = [c_void_p]
-class HessQRCtrl(ctypes.Structure):
-  _fields_ = [("distAED",bType),
+(HESSENBERG_SCHUR_AED,HESSENBERG_SCHUR_MULTIBULGE,HESSENBERG_SCHUR_SIMPLE)= \
+(0,1,2)
+
+lib.ElHessenbergSchurCtrlDefault.argtypes = [c_void_p]
+class HessenbergSchurCtrl(ctypes.Structure):
+  _fields_ = [("winBeg",iType),
+              ("winEnd",iType),
+              ("fullTriangle",bType),
+              ("wantSchurVecs",bType),
+              ("demandConverged",bType),
+              ("alg",c_uint),
+              ("recursiveAED",bType),
+              ("accumulateReflections",bType),
+              ("progress",bType),
+              ("minMultiBulgeSize",iType),
+              ("numShifts",CFUNCTYPE(iType,iType,iType)),
+              ("deflationSize",CFUNCTYPE(iType,iType,iType,iType)),
+              ("sufficientDeflation",CFUNCTYPE(iType,iType)),
+              ("distAED",bType),
               ("blockHeight",iType),("blockWidth",iType)]
   def __init__(self):
-    lib.ElHessQRCtrlDefault(pointer(self))
+    lib.ElHessenbergSchurCtrlDefault(pointer(self))
 
 lib.ElSDCCtrlDefault_s.argtypes = [c_void_p]
 class SDCCtrl_s(ctypes.Structure):
@@ -584,7 +601,7 @@ def __init__(self):
 lib.ElSchurCtrlDefault_s.argtypes = [c_void_p]
 class SchurCtrl_s(ctypes.Structure):
   _fields_ = [("useSDC",bType),
-              ("qrCtrl",HessQRCtrl),
+              ("hessSchurCtrl",HessenbergSchurCtrl),
               ("sdcCtrl",SDCCtrl_s),
               ("time",bType)]
   def __init__(self):
@@ -593,7 +610,7 @@ def __init__(self):
 lib.ElSchurCtrlDefault_d.argtypes = [c_void_p]
 class SchurCtrl_d(ctypes.Structure):
   _fields_ = [("useSDC",bType),
-              ("qrCtrl",HessQRCtrl),
+              ("hessSchurCtrl",HessenbergSchurCtrl),
               ("sdcCtrl",SDCCtrl_d),
               ("time",bType)]
   def __init__(self):