Module: NMatrix::LAPACK
- Defined in:
- lib/nmatrix/lapacke.rb,
lib/nmatrix/lapacke.rb
Overview
Add functions from the LAPACKE C extension to the main LAPACK and BLAS modules. This will overwrite the original functions where applicable.
Class Method Summary collapse
- .geev(matrix, which = :both) ⇒ Object
- .gesdd(matrix, workspace_size = nil) ⇒ Object
- .gesvd(matrix, workspace_size = 1) ⇒ Object
- .posv(uplo, a, b) ⇒ Object
Class Method Details
.geev(matrix, which = :both) ⇒ Object
74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 |
# File 'lib/nmatrix/lapacke.rb', line 74 def geev(matrix, which=:both) raise(StorageTypeError, "LAPACK functions only work on dense matrices") unless matrix.dense? raise(ShapeError, "eigenvalues can only be computed for square matrices") unless matrix.dim == 2 && matrix.shape[0] == matrix.shape[1] jobvl = (which == :both || which == :left) ? :t : false jobvr = (which == :both || which == :right) ? :t : false # Copy the matrix so it doesn't get overwritten. temporary_matrix = matrix.clone n = matrix.shape[0] # Outputs eigenvalues = NMatrix.new([n, 1], dtype: matrix.dtype) # For real dtypes this holds only the real part of the eigenvalues. imag_eigenvalues = matrix.complex_dtype? ? nil : NMatrix.new([n, 1], dtype: matrix.dtype) # For complex dtypes, this is unused. left_output = jobvl ? matrix.clone_structure : nil right_output = jobvr ? matrix.clone_structure : nil NMatrix::LAPACK::lapacke_geev(:row, jobvl, # compute left eigenvectors of A? jobvr, # compute right eigenvectors of A? (left eigenvectors of A**T) n, # order of the matrix temporary_matrix,# input matrix (used as work) n, # leading dimension of matrix eigenvalues,# real part of computed eigenvalues imag_eigenvalues,# imag part of computed eigenvalues left_output, # left eigenvectors, if applicable n, # leading dimension of left_output right_output, # right eigenvectors, if applicable n) # leading dimension of right_output # For real dtypes, transform left_output and right_output into correct forms. # If the j'th and the (j+1)'th eigenvalues form a complex conjugate # pair, then the j'th and (j+1)'th columns of the matrix are # the real and imag parts of the eigenvector corresponding # to the j'th eigenvalue. if !matrix.complex_dtype? complex_indices = [] n.times do |i| complex_indices << i if imag_eigenvalues[i] != 0.0 end if !complex_indices.empty? # For real dtypes, put the real and imaginary parts together eigenvalues = eigenvalues + imag_eigenvalues*Complex(0.0,1.0) left_output = left_output.cast(dtype: NMatrix.upcast(:complex64, matrix.dtype)) if left_output right_output = right_output.cast(dtype: NMatrix.upcast(:complex64, matrix.dtype)) if right_output end complex_indices.each_slice(2) do |i, _| if right_output right_output[0...n,i] = right_output[0...n,i] + right_output[0...n,i+1]*Complex(0.0,1.0) right_output[0...n,i+1] = right_output[0...n,i].complex_conjugate end if left_output left_output[0...n,i] = left_output[0...n,i] + left_output[0...n,i+1]*Complex(0.0,1.0) left_output[0...n,i+1] = left_output[0...n,i].complex_conjugate end end end if which == :both return [eigenvalues, left_output, right_output] elsif which == :left return [eigenvalues, left_output] else return [eigenvalues, right_output] end end |
.gesdd(matrix, workspace_size = nil) ⇒ Object
157 158 159 160 161 162 163 164 165 |
# File 'lib/nmatrix/lapacke.rb', line 157 def gesdd(matrix, workspace_size=nil) result = alloc_svd_result(matrix) m = matrix.shape[0] n = matrix.shape[1] NMatrix::LAPACK::lapacke_gesdd(:row, :a, m, n, matrix, n, result[1], result[0], m, result[2], n) result end |
.gesvd(matrix, workspace_size = 1) ⇒ Object
145 146 147 148 149 150 151 152 153 154 155 |
# File 'lib/nmatrix/lapacke.rb', line 145 def gesvd(matrix, workspace_size=1) result = alloc_svd_result(matrix) m = matrix.shape[0] n = matrix.shape[1] superb = NMatrix.new([[m,n].min], dtype: matrix.abs_dtype) NMatrix::LAPACK::lapacke_gesvd(:row, :a, :a, m, n, matrix, n, result[1], result[0], m, result[2], n, superb) result end |
.posv(uplo, a, b) ⇒ Object
58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 |
# File 'lib/nmatrix/lapacke.rb', line 58 def posv(uplo, a, b) raise(ShapeError, "a must be square") unless a.dim == 2 && a.shape[0] == a.shape[1] raise(ShapeError, "number of rows of b must equal number of cols of a") unless a.shape[1] == b.shape[0] raise(StorageTypeError, "only works with dense matrices") unless a.stype == :dense && b.stype == :dense raise(DataTypeError, "only works for non-integer, non-object dtypes") if a.integer_dtype? || a.object_dtype? || b.integer_dtype? || b.object_dtype? x = b.clone clone = a.clone n = a.shape[0] nrhs = b.shape[1] lapacke_potrf(:row, uplo, n, clone, n) lapacke_potrs(:row, uplo, n, nrhs, clone, n, x, b.shape[1]) x end |