PMR^3: templated eigensolver for Elemental

CMakeLists.txt

@@ -410,6 +410,7 @@ endif()
 # Elemental's mod's of Parallel Multiple Relatively Robust Representations
 # ------------------------------------------------------------------------
 add_subdirectory(external/pmrrr)
+include_directories(external/pmrrr/include)
 if(EL_BUILT_SCALAPACK)
   add_dependencies(pmrrr project_scalapack)
 else()

examples/lapack_like/HermitianEig.cpp

@@ -7,12 +7,88 @@
    http://opensource.org/licenses/BSD-2-Clause
 */
 #include <El.hpp>
+
 using namespace std;
 using namespace El;
 
-// Typedef our real and complex types to 'Real' and 'C' for convenience
-typedef double Real;
-typedef Complex<Real> C;
+template<typename Real>
+void run_example(Int n, bool print)
+{
+    typedef Complex<Real> C;
+
+    // Create a 2d process grid from a communicator. In our case, it is
+    // mpi::COMM_WORLD. There is another constructor that allows you to 
+    // specify the grid dimensions, Grid g( comm, r ), which creates a
+    // grid of height r.
+    Grid g( mpi::COMM_WORLD );
+
+    // Create an n x n complex distributed matrix, 
+    // We distribute the matrix using grid 'g'.
+    //
+    // There are quite a few available constructors, including ones that 
+    // allow you to pass in your own local buffer and to specify the 
+    // distribution alignments (i.e., which process row and column owns the
+    // top-left element)
+    DistMatrix<C> H( n, n, g );
+
+    // Fill the matrix since we did not pass in a buffer. 
+    //
+    // We will fill entry (i,j) with the complex value (i+j,i-j) so that 
+    // the global matrix is Hermitian. However, only one triangle of the 
+    // matrix actually needs to be filled, the symmetry can be implicit.
+    //
+    const Int localHeight = H.LocalHeight();
+    const Int localWidth = H.LocalWidth();
+    for( Int jLoc=0; jLoc<localWidth; ++jLoc )
+    {
+        // Our process owns the rows colShift:colStride:n,
+        //           and the columns rowShift:rowStride:n
+        const Int j = H.GlobalCol(jLoc);
+        for( Int iLoc=0; iLoc<localHeight; ++iLoc )
+        {
+            const Int i = H.GlobalRow(iLoc);
+            H.SetLocal( iLoc, jLoc, C(i+j,i-j) );
+        }
+    }
+
+    // Make a backup of H before we overwrite it within the eigensolver
+    auto HCopy( H );
+
+    // Call the eigensolver. We first create an empty complex eigenvector 
+    // matrix, Q[MC,MR], and an eigenvalue column vector, w[VR,* ]
+    //
+    // Optional: set blocksizes and algorithmic choices here. See the 
+    //           'Tuning' section of the README for details.
+    DistMatrix<Real,VR,STAR> w( g );
+    DistMatrix<C> Q( g );
+    HermitianEig( LOWER, H, w, Q ); 
+
+    if( print )
+    {
+        Print( HCopy, "H" );
+        Print( Q, "Q" );
+        Print( w, "w" );
+    }
+
+    // Check the residual, || H Q - Omega Q ||_F
+    const Real frobH = HermitianFrobeniusNorm( LOWER, HCopy );
+    auto E( Q );
+    DiagonalScale( RIGHT, NORMAL, w, E );
+    Hemm( LEFT, LOWER, C(-1), HCopy, Q, C(1), E );
+    const Real frobResid = FrobeniusNorm( E );
+
+    // Check the orthogonality of Q
+    Identity( E, n, n );
+    Herk( LOWER, ADJOINT, Real(-1), Q, Real(1), E );
+    const Real frobOrthog = HermitianFrobeniusNorm( LOWER, E );
+
+    if( mpi::Rank() == 0 )
+        Output
+        ("|| H ||_F = ",frobH,"\n",
+         "|| H Q - Q Omega ||_F / || A ||_F = ",frobResid/frobH,"\n",
+         "|| Q' Q - I ||_F = ",frobOrthog,"\n");
+}
+
 
 int
 main( int argc, char* argv[] )
@@ -27,80 +103,14 @@ main( int argc, char* argv[] )
     {
         const Int n = Input("--size","size of matrix",100);
         const bool print = Input("--print","print matrices?",false);
+        const bool single_precision = Input("--single", "single precision?", false);
         ProcessInput();
         PrintInputReport();
-
-        // Create a 2d process grid from a communicator. In our case, it is
-        // mpi::COMM_WORLD. There is another constructor that allows you to 
-        // specify the grid dimensions, Grid g( comm, r ), which creates a
-        // grid of height r.
-        Grid g( mpi::COMM_WORLD );
-
-        // Create an n x n complex distributed matrix, 
-        // We distribute the matrix using grid 'g'.
-        //
-        // There are quite a few available constructors, including ones that 
-        // allow you to pass in your own local buffer and to specify the 
-        // distribution alignments (i.e., which process row and column owns the
-        // top-left element)
-        DistMatrix<C> H( n, n, g );
-
-        // Fill the matrix since we did not pass in a buffer. 
-        //
-        // We will fill entry (i,j) with the complex value (i+j,i-j) so that 
-        // the global matrix is Hermitian. However, only one triangle of the 
-        // matrix actually needs to be filled, the symmetry can be implicit.
-        //
-        const Int localHeight = H.LocalHeight();
-        const Int localWidth = H.LocalWidth();
-        for( Int jLoc=0; jLoc<localWidth; ++jLoc )
-        {
-            // Our process owns the rows colShift:colStride:n,
-            //           and the columns rowShift:rowStride:n
-            const Int j = H.GlobalCol(jLoc);
-            for( Int iLoc=0; iLoc<localHeight; ++iLoc )
-            {
-                const Int i = H.GlobalRow(iLoc);
-                H.SetLocal( iLoc, jLoc, C(i+j,i-j) );
-            }
+        if(single_precision) {
+            run_example<float>(n, print);
+        } else {
+            run_example<double>(n, print);
         }
-
-        // Make a backup of H before we overwrite it within the eigensolver
-        auto HCopy( H );
-
-        // Call the eigensolver. We first create an empty complex eigenvector 
-        // matrix, Q[MC,MR], and an eigenvalue column vector, w[VR,* ]
-        //
-        // Optional: set blocksizes and algorithmic choices here. See the 
-        //           'Tuning' section of the README for details.
-        DistMatrix<Real,VR,STAR> w( g );
-        DistMatrix<C> Q( g );
-        HermitianEig( LOWER, H, w, Q ); 
-
-        if( print )
-        {
-            Print( HCopy, "H" );
-            Print( Q, "Q" );
-            Print( w, "w" );
-        }
-
-        // Check the residual, || H Q - Omega Q ||_F
-        const Real frobH = HermitianFrobeniusNorm( LOWER, HCopy );
-        auto E( Q );
-        DiagonalScale( RIGHT, NORMAL, w, E );
-        Hemm( LEFT, LOWER, C(-1), HCopy, Q, C(1), E );
-        const Real frobResid = FrobeniusNorm( E );
-
-        // Check the orthogonality of Q
-        Identity( E, n, n );
-        Herk( LOWER, ADJOINT, Real(-1), Q, Real(1), E );
-        const Real frobOrthog = HermitianFrobeniusNorm( LOWER, E );
-
-        if( mpi::Rank() == 0 )
-            Output
-            ("|| H ||_F = ",frobH,"\n",
-             "|| H Q - Q Omega ||_F / || A ||_F = ",frobResid/frobH,"\n",
-             "|| Q' Q - I ||_F = ",frobOrthog,"\n");
     }
     catch( exception& e ) { ReportException(e); }
 

external/pmrrr/CMakeLists.txt

@@ -10,10 +10,13 @@
 #  http://opensource.org/licenses/BSD-2-Clause
 #
 
+# Create a set of defines for including in templated code
+set(pmrrr_defines "\n")
+
 option(HAVE_SPINLOCKS "Enable if pthread lib supports spinlocks" OFF)
 MARK_AS_ADVANCED(HAVE_SPINLOCKS)
 if(NOT HAVE_SPINLOCKS)
-  add_definitions(-DNOSPINLOCKS)
+  set(pmrrr_defines "${pmrrr_defines}#define NOSPINLOCKS\n")
 endif()
 
 # Include the header directory
@@ -24,9 +27,6 @@ include_directories(BEFORE "${CMAKE_CURRENT_SOURCE_DIR}/include")
 
 # Define the header files installation rules
 # ------------------------------------------
-install(DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/include/
-	DESTINATION ${CMAKE_INSTALL_INCLUDEDIR} 
-  FILES_MATCHING PATTERN "*.h")
 if(WIN32 OR APPLE OR NOT EL_HYBRID)
   # Disable for Windows because of lack of POSIX support;
   # disable for Mac OS X because of deprecation of unnamed semaphores
@@ -36,9 +36,21 @@ else()
   option(DISABLE_PTHREADS "Disable pthreads?" OFF)
 endif()
 if(DISABLE_PTHREADS)
-  add_definitions(-DDISABLE_PTHREADS)
+  set(pmrrr_defines "${pmrrr_defines}#define DISABLE_PTHREADS\n")
 endif()
 
+#Build directory
+configure_file("${CMAKE_CURRENT_SOURCE_DIR}/cmake/global.h.in"
+               "${PROJECT_BINARY_DIR}/include/pmrrr/definitions/global.h"
+               @ONLY)
+#Install directory
+configure_file("${CMAKE_CURRENT_SOURCE_DIR}/cmake/global.h.in"
+               "${CMAKE_INSTALL_INCLUDEDIR}/pmrrr/definitions/lobal.h"
+               @ONLY)
+install(DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/include/
+  DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}
+  FILES_MATCHING PATTERN "*.h" PATTERN "*.hpp")
+
 # Ensure that an MPI C compiler was found
 if(NOT MPI_C_FOUND)
   message(FATAL_ERROR "No MPI C compiler was found, so PMRRR cannot be built")
@@ -74,7 +86,7 @@ set(CMAKE_THREAD_LIBS_INIT ${CMAKE_THREAD_LIBS_INIT} PARENT_SCOPE)
 
 # Define the main library and its link libraries
 # ----------------------------------------------
-file(GLOB_RECURSE PMRRR_SRC RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} "*.c" "*.h")
+file(GLOB_RECURSE PMRRR_SRC RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} "*.cpp" "*.h")
 add_library(pmrrr ${LIBRARY_TYPE} ${PMRRR_SRC})
 if(DISABLE_PTHREADS)
   target_link_libraries(pmrrr ${MPI_C_LIBRARIES} ${MATH_LIBS})