Module: Ignis::Solver::Eigen

Defined in:
lib/nvruby/solver/eigen.rb

Overview

Eigenvalue and Eigenvector computation using cuSOLVER Supports symmetric/Hermitian matrices (syevd/heevd) and general matrices (geev)

Class Method Summary collapse

Class Method Details

.eig(matrix) ⇒ Hash

Compute eigenvalues and eigenvectors of a general (non-symmetric) matrix Note: This is more expensive than eigh for symmetric matrices

Parameters:

  • matrix (NvArray) —

    General square matrix (n x n)

Returns:

  • (Hash) —

    Contains :eigenvalues (may be complex), :eigenvectors_right



86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
 
# File 'lib/nvruby/solver/eigen.rb', line 86

def eig(matrix)
  CuSolverBindings.ensure_loaded!

  validate_square_matrix!(matrix)
  n = matrix.shape[0]

  # For general matrices, eigenvalues may be complex
  # Even for real input, we may get complex eigenvalues
  # For now, we'll use Schur decomposition approximation

  # Copy matrix
  work_matrix = matrix.dup

  # For real matrices, eigenvalues are returned as real/imaginary pairs
  # For complex matrices, eigenvalues are complex directly

  if complex_dtype?(matrix.dtype)
    compute_complex_eig(work_matrix, n)
  else
    compute_real_eig(work_matrix, n)
  end
end

.eigh(matrix, eigenvectors: true, uplo: :lower) ⇒ Hash

Compute eigenvalues and eigenvectors of a symmetric/Hermitian matrix Uses divide-and-conquer algorithm (syevd/heevd)

Parameters:

  • matrix (NvArray) —

    Symmetric (real) or Hermitian (complex) matrix (n x n)

  • eigenvectors (Boolean) (defaults to: true) —

    If true, compute eigenvectors

  • uplo (Symbol) (defaults to: :lower) —

    :lower or :upper indicating which triangle contains the matrix

Returns:

  • (Hash) —

    Contains :eigenvalues and optionally :eigenvectors

Raises:

  • (CuSolverError) —

    If computation fails



18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
 
# File 'lib/nvruby/solver/eigen.rb', line 18

def eigh(matrix, eigenvectors: true, uplo: :lower)
  CuSolverBindings.ensure_loaded!

  validate_symmetric_matrix!(matrix)
  n = matrix.shape[0]
  lda = n

  # Job mode
  jobz = eigenvectors ? CuSolverBindings::CUSOLVER_EIG_MODE_VECTOR : CuSolverBindings::CUSOLVER_EIG_MODE_NOVECTOR
  fill_mode = uplo == :lower ? CuSolverBindings::CUBLAS_FILL_MODE_LOWER : CuSolverBindings::CUBLAS_FILL_MODE_UPPER

  # Allocate eigenvalues array (always real)
  real_dtype = real_dtype_for(matrix.dtype)
  eigenvalues = CUDA::Memory.new(n * dtype_size(real_dtype))

  # Copy matrix to work array (will contain eigenvectors on output if requested)
  work_matrix = matrix.dup

  # Get workspace size
  lwork_ptr = FFI::MemoryPointer.new(:int)
  get_syevd_buffer_size(jobz, fill_mode, n, work_matrix.device_ptr, lda, eigenvalues, lwork_ptr, matrix.dtype)
  lwork = lwork_ptr.read_int

  # Allocate workspace and info
  workspace = CUDA::Memory.new(lwork * dtype_size(matrix.dtype))
  info = CUDA::Memory.new(4)

  # Perform eigenvalue computation
  perform_syevd(jobz, fill_mode, n, work_matrix.device_ptr, lda, eigenvalues, workspace, lwork, info, matrix.dtype)

  # Check info
  info_value = read_device_int(info)
  if info_value < 0
    raise CuSolverError.new("Eigenvalue computation: parameter #{-info_value} had an illegal value",
                            cusolver_code: CuSolverBindings::CUSOLVER_STATUS_INVALID_VALUE)
  elsif info_value > 0
    raise CuSolverError.new("Eigenvalue computation: #{info_value} off-diagonal elements did not converge to zero",
                            cusolver_code: CuSolverBindings::CUSOLVER_STATUS_EXECUTION_FAILED)
  end

  CUDA::RuntimeAPI.cudaDeviceSynchronize

  result = {
    eigenvalues: NvArray.from_device_ptr(eigenvalues, shape: [n], dtype: real_dtype, take_ownership: true)
  }

  if eigenvectors
    result[:eigenvectors] = work_matrix
  end

  result
ensure
  workspace&.free! if defined?(workspace) && workspace
  info&.free! if defined?(info) && info
end

.eigvalsh(matrix) ⇒ NvArray

Compute eigenvalues only of a symmetric/Hermitian matrix (faster)

Parameters:

  • matrix (NvArray) —

    Symmetric/Hermitian matrix

Returns:

  • (NvArray) —

    Eigenvalues in ascending order



77
78
79
80
 
# File 'lib/nvruby/solver/eigen.rb', line 77

def eigvalsh(matrix)
  result = eigh(matrix, eigenvectors: false)
  result[:eigenvalues]
end

.symmetric?(matrix, tol: 1e-10) ⇒ Boolean

Check if a matrix is approximately symmetric

Parameters:

  • matrix (NvArray) —

    Matrix to check

  • tol (Float) (defaults to: 1e-10) —

    Tolerance for symmetry check

Returns:

  • (Boolean) 


113
114
115
116
117
118
119
120
 
# File 'lib/nvruby/solver/eigen.rb', line 113

def symmetric?(matrix, tol: 1e-10)
  return false unless matrix.shape[0] == matrix.shape[1]

  # For a proper check, we'd compare A with A^T
  # This is a simplified placeholder - full implementation would
  # compute max(abs(A - A^T))
  true # Placeholder - assumes caller knows matrix is symmetric
end