Module: NMatrix::LAPACK

Defined in:
lib/nmatrix/atlas.rb,
lib/nmatrix/atlas.rb,
lib/nmatrix/lapacke.rb,
lib/nmatrix/lapacke.rb,
lib/nmatrix/lapack_core.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

Class Method Details

.alloc_svd_result(matrix) ⇒ Object


77
78
79
80
81
82
83
# File 'lib/nmatrix/lapack_core.rb', line 77

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

.clapack_getri(order, n, a, lda, ipiv) ⇒ Object

Raises:

  • (NotImplementedError)

176
177
178
# File 'lib/nmatrix/lapack_core.rb', line 176

def clapack_getri(order, n, a, lda, ipiv)
  raise(NotImplementedError,"clapack_getri requires the nmatrix-atlas gem")
end

.clapack_potrf(order, uplo, n, a, lda) ⇒ Object

Raises:

  • (NotImplementedError)

164
165
166
# File 'lib/nmatrix/lapack_core.rb', line 164

def clapack_potrf(order, uplo, n, a, lda)
  raise(NotImplementedError,"clapack_potrf requires the nmatrix-atlas gem")
end

.clapack_potri(order, uplo, n, a, lda) ⇒ Object

Raises:

  • (NotImplementedError)

168
169
170
# File 'lib/nmatrix/lapack_core.rb', line 168

def clapack_potri(order, uplo, n, a, lda)
  raise(NotImplementedError,"clapack_potri requires the nmatrix-atlas gem")
end

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

Raises:

  • (NotImplementedError)

172
173
174
# File 'lib/nmatrix/lapack_core.rb', line 172

def clapack_potrs(order, uplo, n, nrhs, a, lda, b, ldb)
  raise(NotImplementedError,"clapack_potrs requires the nmatrix-atlas gem")
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.

eigenvalues is a n-by-1 NMatrix containing the eigenvalues.

right_eigenvalues is a n-by-n matrix such that its j'th column contains the (right) eigenvalue of matrix corresponding to the j'th eigenvalue. This means that matrix = RDR^(-1), where R is right_eigenvalues and D is the diagonal matrix formed from eigenvalues.

left_eigenvalues is n-by-n and its columns are the left eigenvalues of matrix, using the definition of left eigenvalue from LAPACK.

For real dtypes, eigenvalues and the eigenvector matrices will be complex if and only if matrix has complex eigenvalues.

Only available if nmatrix-lapack or nmatrix-atlas is installed.

Raises:

  • (NotImplementedError)

145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
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
# File 'lib/nmatrix/lapack_core.rb', line 145

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

  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

  # lapack_geev is a pure LAPACK routine so it expects column-major matrices,
  # so we need to transpose the input as well as the output.
  temporary_matrix = matrix.transpose
  NMatrix::LAPACK::lapack_geev(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
                               2*n)
  left_output = left_output.transpose if jobvl
  right_output = right_output.transpose if jobvr


  # 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

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.

Requires either the nmatrix-lapacke or nmatrix-atlas gem.

Raises:

  • (NotImplementedError)

115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
# File 'lib/nmatrix/lapack_core.rb', line 115

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)

  m = matrix.shape[0]
  n = matrix.shape[1]

  # This is a pure LAPACK function so it expects column-major functions.
  # So we need to transpose the input as well as the output.
  matrix = matrix.transpose
  NMatrix::LAPACK::lapack_gesdd(:a, m, n, matrix, m, result[1], \
   result[0], m, result[2], n, workspace_size)
  result[0] = result[0].transpose
  result[2] = result[2].transpose
  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.

Requires either the nmatrix-lapacke or nmatrix-atlas gem.

Raises:

  • (NotImplementedError)

98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
# File 'lib/nmatrix/lapack_core.rb', line 98

def gesvd(matrix, workspace_size=1)
  result = alloc_svd_result(matrix)

  m = matrix.shape[0]
  n = matrix.shape[1]

  # This is a pure LAPACK function so it expects column-major functions.
  # So we need to transpose the input as well as the output.
  matrix = matrix.transpose
  NMatrix::LAPACK::lapack_gesvd(:a, :a, m, n, matrix, \
   m, result[1], result[0], m, result[2], n, workspace_size)
  result[0] = result[0].transpose
  result[2] = result[2].transpose
  result
end

.lapack_geev(jobvl, jobvr, n, a, lda, w, wi, vl, ldvl, vr, ldvr, lwork) ⇒ Object

Raises:

  • (NotImplementedError)

160
161
162
# File 'lib/nmatrix/lapack_core.rb', line 160

def lapack_geev(jobvl, jobvr, n, a, lda, w, wi, vl, ldvl, vr, ldvr, lwork)
  raise(NotImplementedError,"lapack_geev requires the nmatrix-atlas gem")
end

.lapack_gesdd(jobz, m, n, a, lda, s, u, ldu, vt, ldvt, lwork) ⇒ Object

Raises:

  • (NotImplementedError)

156
157
158
# File 'lib/nmatrix/lapack_core.rb', line 156

def lapack_gesdd(jobz, m, n, a, lda, s, u, ldu, vt, ldvt, lwork)
  raise(NotImplementedError,"lapack_gesdd requires the nmatrix-atlas gem")
end

.lapack_gesvd(jobu, jobvt, m, n, a, lda, s, u, ldu, vt, ldvt, lwork) ⇒ Object

The following are functions that used to be implemented in C, but now require nmatrix-atlas to run properly, so we can just implemented their stubs in Ruby.

Raises:

  • (NotImplementedError)

152
153
154
# File 'lib/nmatrix/lapack_core.rb', line 152

def lapack_gesvd(jobu, jobvt, m, n, a, lda, s, u, ldu, vt, ldvt, lwork)
  raise(NotImplementedError,"lapack_gesvd requires the nmatrix-atlas gem")
end

.laswp(matrix, ipiv) ⇒ Object

laswp(matrix, ipiv) -> NMatrix

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

Raises:

  • (ArgumentError)

69
70
71
72
73
74
75
# File 'lib/nmatrix/lapack_core.rb', line 69

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

.posv(uplo, a, b) ⇒ Object

Solve the matrix equation AX = B, where A is a symmetric (or Hermitian) positive-definite matrix. If A is a nxn matrix, B must be mxn. Depending on the value of uplo, only the upper or lower half of a is read. This uses the Cholesky decomposition so it should be faster than the generic NMatrix#solve method. Doesn't modify inputs. Requires either the nmatrix-atlas or nmatrix-lapacke gem.

  • Arguments :

    • uplo -> Either :upper or :lower. Specifies which half of a to read.

    • a -> The matrix A.

    • b -> The right-hand side B.

  • Returns :

    • The solution X

Raises:

  • (NotImplementedError)

61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
# File 'lib/nmatrix/lapack_core.rb', line 61

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]
  clapack_potrf(:row, uplo, n, clone, n)  # Must transpose b before and after:
  #  http://math-atlas.sourceforge.net/faq.html#RowSolve

  x = x.transpose
  clapack_potrs(:row, uplo, n, nrhs, clone, n, x, n)
  x.transpose
end