FeastKit.jl

A Julia implementation of the FeastKit eigenvalue solver for finding eigenvalues and eigenvectors of large-scale eigenvalue problems within a specified region.

Overview

FeastKit.jl is a pure Julia translation of the original FEAST library. FEAST is a numerical algorithm for solving both standard and generalized eigenvalue problems by computing eigenvalues located inside a given region in the complex plane.

Key Features

• Interval-based eigenvalue computation: Find eigenvalues in specified intervals or regions
• Multiple matrix formats: Support for dense, sparse, and banded matrices
• Generalized eigenvalue problems: Solve both Ax = λx and Ax = λBx problems
• Complex arithmetic: Handle both real symmetric/Hermitian and general complex matrices
• Precision flexibility: Works seamlessly with Float64/ComplexF64 and Float32/ComplexF32
• Parallel computation: Multi-threaded, distributed, and optional MPI-based contour integration
• Reverse Communication Interface (RCI): Advanced interface for custom linear solvers
• Polynomial eigenvalue problems: Support for polynomial eigenvalue problems

Installation

using Pkg
Pkg.add("FeastKit")

Quick Start

Basic Usage

using FeastKit
using LinearAlgebra

# Create a test matrix
n = 100
A = diagm(-1 => -ones(n-1), 0 => 2*ones(n), 1 => -ones(n-1))

# Find eigenvalues in the interval [0.5, 1.5]
result = feast(A, (0.5, 1.5), M0=10)

println("Found $(result.M) eigenvalues")
println("Eigenvalues: ", result.lambda)

Generalized Eigenvalue Problems

# For generalized problem Ax = λBx
B = diagm(0 => ones(n))
result = feast(A, B, (0.5, 1.5), M0=10)

Sparse Matrices

using SparseArrays

# Create sparse matrix
A_sparse = spdiagm(-1 => -ones(n-1), 0 => 2*ones(n), 1 => -ones(n-1))
result = feast(A_sparse, (0.5, 1.5), M0=10)

Complex Eigenvalue Problems

# For general complex matrices, use circular search region
A_complex = randn(ComplexF64, 50, 50)
B_complex = Matrix{ComplexF64}(I, 50, 50)

# Search in circle centered at origin with radius 2
center = 0.0 + 0.0im
radius = 2.0
result = feast_general(A_complex, B_complex, center, radius, M0=15)

Advanced Usage

Custom FeastKit Parameters

# Initialize Feast parameters
fpm = feastinit()

# Customize parameters
feast_set_defaults!(fpm.fpm, 
                   print_level=1,
                   integration_points=16,  # More integration points
                   tolerance_exp=14,       # Higher precision
                   max_refinement=30)      # More refinement loops

# Use custom parameters
result = feast(A, (0.5, 1.5), M0=10, fpm=fpm.fpm)

Banded Matrices

# For banded matrices stored in LAPACK format
n = 100
k = 2  # Number of super-diagonals
A_banded = zeros(k+1, n)
# ... fill A_banded with appropriate values ...

result = feast_banded(A_banded, k, (0.5, 1.5), M0=10)

MPI Acceleration (optional)

MPI wrappers mirror the pdfeast_* routines in the original FEAST library. They are only loaded when MPI.jl is available and FEASTKIT_ENABLE_MPI=true is set in the environment:

export FEASTKIT_ENABLE_MPI=true
mpiexec -n 4 julia --project -e 'using FeastKit; println(FeastKit.mpi_available())'

Once enabled you can call the MPI interfaces directly:

if FeastKit.mpi_available()
    using MPI
    MPI.Init()
    A = Matrix{Float64}(I, 50, 50)
    B = Matrix{Float64}(I, 50, 50)
    fpm = zeros(Int, 64)
    feastinit!(fpm)
    FeastKit.mpi_feast_sygv!(A, B, 0.0, 2.0, 16, fpm)
    MPI.Finalize()
end

Matrix-Free Operations

# Define matrix-vector operations
function A_mul!(y, x)
    # Implement y = A*x
    y .= A * x
end

function B_mul!(y, x)
    # Implement y = B*x  
    y .= B * x
end

# Use matrix-free interface
result = feast_matvec(A_mul!, B_mul!, n, (0.5, 1.5), M0=10)

Parallel Computing

FeastKit.jl supports parallel computation where each contour integration point is solved independently, leading to significant speedups for large problems.

Multi-threaded Execution

# Enable parallel computation using threads
result = feast(A, (0.5, 1.5), M0=10, parallel=:threads)

# Or use the explicit parallel interface
result = feast_parallel(A, B, (0.5, 1.5), M0=10, use_threads=true)

Distributed Computing

using Distributed

# Add worker processes
addprocs(4)

# Use distributed computing for contour integration
result = feast(A, (0.5, 1.5), M0=10, parallel=:distributed)

MPI Support for HPC Clusters

FeastKit.jl provides full MPI support for high-performance computing clusters:

using MPI
using FeastKit

# Initialize MPI (if not already done)
MPI.Init()

# Basic MPI FeastKit
result = feast(A, B, (0.5, 1.5), M0=10, parallel=:mpi)

# Explicit MPI interface with communicator
comm = MPI.COMM_WORLD
result = mpi_feast(A, B, (0.5, 1.5), M0=10, comm=comm)

# Hybrid MPI + threading (best for modern HPC)
result = feast_hybrid(A, B, (0.5, 1.5), M0=10, 
                     comm=comm, use_threads_per_rank=true)

Running on HPC Systems

Create a job script for SLURM/PBS:

#!/bin/bash
#SBATCH --nodes=4
#SBATCH --ntasks-per-node=8
#SBATCH --cpus-per-task=4
#SBATCH --time=01:00:00

module load julia mpi

# Run with MPI + threading
mpirun -np 32 julia --threads=4 feast_mpi_example.jl

Julia script (feast_mpi_example.jl):

using MPI
using FeastKit
using LinearAlgebra

MPI.Init()

# Create distributed problem
comm = MPI.COMM_WORLD
rank = MPI.Comm_rank(comm)
size = MPI.Comm_size(comm)

if rank == 0
    println("Running FeastKit on $size MPI processes with $(Threads.nthreads()) threads each")
end

# Large eigenvalue problem
n = 10000
A = spdiagm(-1 => -ones(n-1), 0 => 2*ones(n), 1 => -ones(n-1))
B = sparse(I, n, n)

# Hybrid MPI + threading execution
result = feast_hybrid(A, B, (0.5, 1.5), M0=20, 
                     comm=comm, use_threads_per_rank=true)

if rank == 0
    println("Found $(result.M) eigenvalues")
    println("Computation time measured on each rank")
end

MPI.Finalize()

Parallel RCI Interface

For advanced users, a parallel RCI interface is available:

using FeastKit

# Create parallel state
state = ParallelFeastState{Float64}(ne=8, M0=10, use_parallel=true, use_threads=true)

# Initialize workspace
N = size(A, 1)
work = randn(N, M0)
workc = zeros(ComplexF64, N, M0)
Aq = zeros(M0, M0)
Sq = zeros(M0, M0)
lambda = zeros(M0)
q = zeros(N, M0)
res = zeros(M0)

# RCI loop
while true
    pfeast_srci!(state, N, work, workc, Aq, Sq, fpm, 
                Emin, Emax, M0, lambda, q, res)
    
    if state.ijob == Feast_RCI_PARALLEL_SOLVE.value
        # Solve all contour points in parallel
        pfeast_compute_all_contour_points!(state, A, B, work, M0)
        
    elseif state.ijob == Feast_RCI_MULT_A.value
        # Compute A*q for residual calculation
        M = state.mode
        work[:, 1:M] .= A * q[:, 1:M]
        
    elseif state.ijob == Feast_RCI_DONE.value
        break
    end
end

Automatic Backend Selection

# FeastKit automatically selects the best available backend
result = feast(A, B, (0.5, 1.5), M0=10, parallel=:auto)

# Manual backend selection
result = feast(A, B, (0.5, 1.5), M0=10, parallel=:mpi)      # Force MPI
result = feast(A, B, (0.5, 1.5), M0=10, parallel=:threads)  # Force threading
result = feast(A, B, (0.5, 1.5), M0=10, parallel=:serial)   # Force serial

Algorithm Overview

FeastKit uses contour integration in the complex plane to compute eigenvalues. The key steps are:

1. Contour Definition: Define an integration contour enclosing the desired eigenvalues
2. Moment Computation: Compute spectral projector moments using numerical integration
3. Subspace Extraction: Extract eigenspace using computed moments
4. Refinement: Iteratively refine the solution until convergence

The algorithm is particularly effective for:

• Large sparse matrices
• Finding eigenvalues in specific intervals
• Parallel computation (via the original PFEAST routines)

Result Structure

FeastKit returns a FeastResult object containing:

struct FeastResult{T<:Real, VT}
    lambda::Vector{T}      # Computed eigenvalues
    q::Matrix{VT}          # Computed eigenvectors  
    M::Int                 # Number of eigenvalues found
    res::Vector{T}         # Residual norms
    info::Int              # Exit status (0 = success)
    epsout::T              # Final residual
    loop::Int              # Number of refinement loops
end

Error Codes

• info = 0: Successful convergence
• info = 1: Invalid matrix size
• info = 2: Invalid search subspace size
• info = 3: Invalid search interval
• info = 5: No convergence achieved
• info = 6: Memory allocation error
• info = 8: LAPACK error

Performance Tips

1. Choose appropriate M0: Set M0 slightly larger than expected number of eigenvalues
2. Integration points: More points (fpm[2]) improve accuracy but increase cost
3. Sparse matrices: Use sparse format for large problems with few non-zeros
4. Initial guess: Provide good initial guess when available (fpm[5] = 1)
5. Parallel execution: Use parallel=:auto to automatically select the best backend
6. HPC clusters: Use parallel=:mpi or feast_hybrid() for optimal cluster performance
7. Hybrid parallelism: Combine MPI processes with threading for maximum performance
8. Load balancing: FeastKit automatically distributes contour points for optimal load balancing

Limitations

This is a Julia translation of the original FEAST library. Some features from the original may not be fully implemented:

• Custom contour integration is partially implemented
• Some advanced PFEAST (parallel) features are not included
• Matrix-free interface is simplified compared to the original

References

1. E. Polizzi, “Density-matrix-based algorithm for solving eigenvalue problems”, Physical Review B 79, 115112 (2009)
2. Feast official website: http://www.feast-solver.org

Contributing

Contributions are welcome! Please feel free to submit issues, feature requests, or pull requests.

License

This project is licensed under the same terms as the original Feast library.

This site is open source. Improve this page.