Module: NMatrix::LAPACK

Defined in:
lib/nmatrix/lapack.rb

Class Method Summary collapse

Class Method Details

.alloc_evd_result(matrix) ⇒ Object 



179
180
181
182
183
184
185
186
 
# File 'lib/nmatrix/lapack.rb', line 179

def alloc_evd_result(matrix)
  [
    NMatrix.new(matrix.shape[0], dtype: matrix.dtype),
    NMatrix.new(matrix.shape[0], dtype: matrix.dtype),
    NMatrix.new([matrix.shape[0],1], dtype: matrix.dtype),
    NMatrix.new([matrix.shape[0],1], dtype: matrix.dtype),
  ]
end

.alloc_svd_result(matrix) ⇒ Object 



135
136
137
138
139
140
141
 
# File 'lib/nmatrix/lapack.rb', line 135

def alloc_svd_result(matrix)
  [
    NMatrix.new(matrix.shape[0], dtype: matrix.dtype),
    NMatrix.new([matrix.shape[0],1], dtype: matrix.dtype),
    NMatrix.new(matrix.shape[1], dtype: matrix.dtype)
  ]
end

.clapack_gesv(order, n, nrhs, a, lda, b, ldb, ipiv = nil) ⇒ Object

call-seq:

clapack_gesv(order, n, nrhs, a, lda, ipiv, b, ldb) -> NMatrix

Computes the solution to a system of linear equations

A * X = B,

where A is an N-by-N matrix and X and B are N-by-NRHS matrices.

The LU factorization used to factor A is dependent on the order parameter, as detailed in the leading comments of clapack_getrf.

The factored form of A is then used solve the system of equations A * X = B.

A is overwritten with the appropriate LU factorization, and B, which contains B on input, is overwritten with the solution X on output.

From ATLAS 3.8.0.

Note: Because this function is implemented in Ruby, the ATLAS lib version is never called! For float32, float64, complex64, and complex128, the ATLAS lib versions of getrf and getrs will be called.

  • Arguments :

    • order ->

    • n ->

    • nrhs ->

    • a ->

    • lda ->

    • b ->

    • ldb ->

    • ipiv -> A pivot array (if nil, one will be generated with clapack_getrf)

  • Returns : -

  • Raises :

    • ++ ->



75
76
77
78
 
# File 'lib/nmatrix/lapack.rb', line 75

def clapack_gesv(order, n, nrhs, a, lda, b, ldb, ipiv=nil)
  ipiv ||= clapack_getrf(order, n, n, a, lda)
  clapack_getrs(order, :no_transpose, n, nrhs, a, lda, ipiv, b, ldb)
end

.clapack_posv(order, uplo, n, nrhs, a, lda, b, ldb) ⇒ Object

call-seq:

clapack_posv(order, uplo, n ,nrhs, a, lda, b, ldb) -> ...

TODO Complete this description.

Computes the solution to a real system of linear equations

A * X = B,

where A is an N-by-N symmetric positive definite matrix and X and B are N-by-NRHS matrices.

The Cholesky decomposition is used to factor A as

A = U**T* U,  if UPLO = 'U', or
A = L * L**T,  if UPLO = 'L',

where U is an upper triangular matrix and L is a lower triangular matrix. The factored form of A is then used to solve the system of equations A * X = B.

From ATLAS 3.8.0.

Note: Because this function is implemented in Ruby, the ATLAS lib version is never called! For float32, float64, complex64, and complex128, the ATLAS lib versions of potrf and potrs will be called.

  • Arguments :

    • order ->

    • uplo ->

    • n ->

    • nrhs ->

    • a ->

    • lda ->

    • b ->

    • ldb ->

  • Returns : -

  • Raises :

    • ++ ->



118
119
120
121
 
# File 'lib/nmatrix/lapack.rb', line 118

def clapack_posv(order, uplo, n, nrhs, a, lda, b, ldb)
  clapack_potrf(order, uplo, n, a, lda)
  clapack_potrs(order, uplo, n, nrhs, a, lda, b, ldb)
end

.geev(matrix, which = :both) ⇒ Object

call-seq:

geev(matrix) -> [eigenvalues, left_eigenvectors, right_eigenvectors]
geev(matrix, :left) -> [eigenvalues, left_eigenvectors]
geev(matrix, :right) -> [eigenvalues, right_eigenvectors]

Perform eigenvalue decomposition on a matrix using LAPACK’s xGEEV function.



197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
 
# File 'lib/nmatrix/lapack.rb', line 197

def geev(matrix, which=:both)
  jobvl = (which == :both || which == :left) ? :left : false
  jobvr = (which == :both || which == :right) ? :right : false

  # Copy the matrix so it doesn't get overwritten.
  temporary_matrix = matrix.clone

  # Outputs
  real_eigenvalues = NMatrix.new([matrix.shape[0], 1], dtype: matrix.dtype)
  imag_eigenvalues = NMatrix.new([matrix.shape[0], 1], dtype: matrix.dtype)

  left_output      = jobvl == :left ? matrix.clone_structure : NMatrix.new(1, dtype: matrix.dtype)
  right_output     = jobvr == :right ? matrix.clone_structure : NMatrix.new(1, dtype: matrix.dtype)

  NMatrix::LAPACK::lapack_geev(jobvl, # compute left eigenvectors of A?
                               jobvr, # compute right eigenvectors of A? (left eigenvectors of A**T)
                               matrix.shape[0], # order of the matrix
                               temporary_matrix,# input matrix (used as work)
                               matrix.shape[0], # leading dimension of matrix
                               real_eigenvalues,# real part of computed eigenvalues
                               imag_eigenvalues,# imag part of computed eigenvalues
                               left_output,     # left eigenvectors, if applicable
                               left_output.shape[0], # leading dimension of left_output
                               right_output,    # right eigenvectors, if applicable
                               right_output.shape[0], # leading dimension of right_output
                               2*matrix.shape[0])

  # Put the real and imaginary parts together
  eigenvalues = real_eigenvalues.to_a.flatten.map.with_index do |real,i|
    imag_eigenvalues[i] != 0 ? Complex(real, imag_eigenvalues[i]) : real
  end

  if which == :both
    return [eigenvalues, left_output.transpose, right_output.transpose]
  elsif which == :left
    return [eigenvalues, left_output.transpose]
  else
    return [eigenvalues, right_output]
  end
end

.gesdd(matrix, workspace_size = nil) ⇒ Object

call-seq:

gesdd(matrix) -> [u, sigma, v_transpose]
gesdd(matrix) -> [u, sigma, v_conjugate_transpose] # complex

Compute the singular value decomposition of a matrix using LAPACK’s GESDD function. This uses a divide-and-conquer strategy. See also #gesvd.

Optionally accepts a workspace_size parameter, which will be honored only if it is larger than what LAPACK requires.



171
172
173
174
175
176
177
 
# File 'lib/nmatrix/lapack.rb', line 171

def gesdd(matrix, workspace_size=nil)
  min_workspace_size = matrix.shape.min * (6 + 4 * matrix.shape.min) + matrix.shape.max
  workspace_size = min_workspace_size if workspace_size.nil? || workspace_size < min_workspace_size
  result = alloc_svd_result(matrix)
  NMatrix::LAPACK::lapack_gesdd(:a, matrix.shape[0], matrix.shape[1], matrix, matrix.shape[0], result[1], result[0], matrix.shape[0], result[2], matrix.shape[1], workspace_size)
  result
end

.gesvd(matrix, workspace_size = 1) ⇒ Object

call-seq:

gesvd(matrix) -> [u, sigma, v_transpose]
gesvd(matrix) -> [u, sigma, v_conjugate_transpose] # complex

Compute the singular value decomposition of a matrix using LAPACK’s GESVD function.

Optionally accepts a workspace_size parameter, which will be honored only if it is larger than what LAPACK requires.



154
155
156
157
158
 
# File 'lib/nmatrix/lapack.rb', line 154

def gesvd(matrix, workspace_size=1)
  result = alloc_svd_result(matrix)
  NMatrix::LAPACK::lapack_gesvd(:a, :a, matrix.shape[0], matrix.shape[1], matrix, matrix.shape[0], result[1], result[0], matrix.shape[0], result[2], matrix.shape[1], workspace_size)
  result
end

.laswp(matrix, ipiv) ⇒ Object

laswp(matrix, ipiv) -> NMatrix

Permute the columns of a matrix (in-place) according to the Array ipiv.

Raises:

  • (ArgumentError) 


127
128
129
130
131
132
133
 
# File 'lib/nmatrix/lapack.rb', line 127

def laswp(matrix, ipiv)
  raise(ArgumentError, "expected NMatrix for argument 0") unless matrix.is_a?(NMatrix)
  raise(StorageTypeError, "LAPACK functions only work on :dense NMatrix instances") unless matrix.stype == :dense
  raise(ArgumentError, "expected Array ipiv to have no more entries than NMatrix a has columns") if ipiv.size > matrix.shape[1]

  clapack_laswp(matrix.shape[0], matrix, matrix.shape[1], 0, ipiv.size-1, ipiv, 1)
end