Example 11: Hermitian Eigenvalue Problems

This example demonstrates computing eigenvalues and eigenvectors for Hermitian (or symmetric) matrices.

Key Concepts

  1. Eigenvalues Only: Computing just the eigenvalues \(\lambda\) using eig_vals (heev).

  2. Eigenvectors: Computing eigenvalues and eigenvectors \(Z\) such that \(A Z = Z \Lambda\).

  3. Simplified vs Traditional: Using slate::eig vs slate::heev.

C++ Example

Setup (Lines 28-32)

slate::HermitianMatrix<scalar_type> A( ... );
slate::Matrix<scalar_type> Z( ... );
std::vector<real_t> Lambda( n );

  • A: Input Hermitian matrix n by n.

  • Z: Matrix to store eigenvectors. Must be n by n.

  • Lambda: Vector to store eigenvalues. Must be size n.

Eigenvalues Only (Lines 43-49)

slate::eig_vals( A, Lambda );  // simplified API
// or
slate::eig( A, Lambda );       // simplified API
// or
slate::heev( A, Lambda );      // traditional API

Computes only the eigenvalues. A is overwritten by the tridiagonal factors (if heev logic is followed, though conceptually A is destroyed). Lambda contains the eigenvalues in ascending order.

  • Note: heev stands for Hermitian EigenValues.

Eigenvalues and Eigenvectors (Lines 63-68)

slate::eig( A, Lambda, Z );    // simplified API
// or
slate::heev( A, Lambda, Z );   // traditional API

Computes both eigenvalues and eigenvectors.

  • Z is overwritten with the eigenvectors.

  • The columns of Z correspond to the eigenvalues in Lambda.

  1// ex11_hermitian_eig.cc
  2// Solve Hermitian eigenvalues A = Z Lambda Z^H
  3
  4/// !!!   Lines between `//---------- begin label`          !!!
  5/// !!!             and `//---------- end label`            !!!
  6/// !!!   are included in the SLATE Users' Guide.           !!!
  7
  8#include <slate/slate.hh>
  9
 10#include "util.hh"
 11
 12int mpi_size = 0;
 13int mpi_rank = 0;
 14int grid_p = 0;
 15int grid_q = 0;
 16
 17//------------------------------------------------------------------------------
 18template <typename scalar_type>
 19void test_hermitian_eig()
 20{
 21    using real_t = blas::real_type<scalar_type>;
 22
 23    print_func( mpi_rank );
 24
 25    int64_t n=1000, nb=256;
 26
 27    //---------- begin eig1
 28    slate::HermitianMatrix<scalar_type>
 29        A( slate::Uplo::Lower, n, nb, grid_p, grid_q, MPI_COMM_WORLD );
 30    slate::Matrix<scalar_type>
 31        Z( n, n, nb, grid_p, grid_q, MPI_COMM_WORLD );
 32    std::vector<real_t> Lambda( n );
 33    // ...
 34    //---------- end eig1
 35
 36    A.insertLocalTiles();
 37    Z.insertLocalTiles();
 38    random_matrix( A );
 39
 40    //----------------------------------------
 41    //---------- begin eig2
 42    // A = Z Lambda Z^H, eigenvalues only
 43    slate::eig_vals( A, Lambda );  // simplified API, or
 44    //---------- end eig2
 45
 46    random_matrix( A );
 47
 48    //---------- begin eig3
 49    slate::eig( A, Lambda );       // simplified API
 50    //---------- end eig3
 51
 52    random_matrix( A );
 53
 54    //---------- begin eig4
 55    slate::heev( A, Lambda );      // traditional API
 56    //---------- end eig4
 57
 58    random_matrix( A );
 59
 60    //----------------------------------------
 61    //---------- begin eig5
 62    // A = Z Lambda Z^H, eigenvalues and eigenvectors
 63    slate::eig( A, Lambda, Z );    // simplified API
 64    //---------- end eig5
 65    random_matrix( A );
 66
 67    //---------- begin eig6
 68    slate::heev( A, Lambda, Z );   // traditional API
 69    //---------- end eig6
 70}
 71
 72//------------------------------------------------------------------------------
 73int main( int argc, char** argv )
 74{
 75    try {
 76        // Parse command line to set types for s, d, c, z precisions.
 77        bool types[ 4 ];
 78        parse_args( argc, argv, types );
 79
 80        int provided = 0;
 81        slate_mpi_call(
 82            MPI_Init_thread( &argc, &argv, MPI_THREAD_MULTIPLE, &provided ) );
 83        assert( provided == MPI_THREAD_MULTIPLE );
 84
 85        slate_mpi_call(
 86            MPI_Comm_size( MPI_COMM_WORLD, &mpi_size ) );
 87
 88        slate_mpi_call(
 89            MPI_Comm_rank( MPI_COMM_WORLD, &mpi_rank ) );
 90
 91        // Determine p-by-q grid for this MPI size.
 92        // Hermitian eig requires square MPI grid.
 93        grid_size_square( mpi_size, &grid_p, &grid_q );
 94        if (mpi_rank == 0) {
 95            printf( "mpi_size %d, grid_p %d, grid_q %d\n",
 96                    mpi_size, grid_p, grid_q );
 97        }
 98
 99        // so random_matrix is different on different ranks.
100        srand( 100 * mpi_rank );
101
102        if (types[ 0 ]) {
103            test_hermitian_eig< float >();
104        }
105
106        if (types[ 1 ]) {
107            test_hermitian_eig< double >();
108        }
109
110        if (types[ 2 ]) {
111            test_hermitian_eig< std::complex<float> >();
112        }
113
114        if (types[ 3 ]) {
115            test_hermitian_eig< std::complex<double> >();
116        }
117
118        slate_mpi_call(
119            MPI_Finalize() );
120    }
121    catch (std::exception const& ex) {
122        fprintf( stderr, "%s", ex.what() );
123        return 1;
124    }
125    return 0;
126}