Class: NMatrix

Inherits:
Object
  • Object
show all
Defined in:
lib/nmatrix/lapack_ext_common.rb,
lib/nmatrix/lapacke.rb

Overview

–

NMatrix

A linear algebra library for scientific computation in Ruby. NMatrix is part of SciRuby.

NMatrix was originally inspired by and derived from NArray, by Masahiro Tanaka: narray.rubyforge.org

Copyright Information

SciRuby is Copyright © 2010 - 2014, Ruby Science Foundation NMatrix is Copyright © 2012 - 2014, John Woods and the Ruby Science Foundation

Please see LICENSE.txt for additional copyright notices.

Contributing

By contributing source code to SciRuby, you agree to be bound by our Contributor Agreement:

lapack_ext_common.rb

Contains functions shared by nmatrix-atlas and nmatrix-lapacke gems. ++

Defined Under Namespace

Modules: BLAS, LAPACK

Class Method Summary collapse

Instance Method Summary collapse

Class Method Details

.register_lapack_extension(name) ⇒ Object 



30
31
32
33
34
35
36
 
# File 'lib/nmatrix/lapack_ext_common.rb', line 30

def NMatrix.register_lapack_extension(name)
  if (defined? @@lapack_extension)
    raise "Attempting to load #{name} when #{@@lapack_extension} is already loaded. You can only load one LAPACK extension."
  end

  @@lapack_extension = name
end

Instance Method Details

#dot(right_v) ⇒ Object 



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
 
# File 'lib/nmatrix/lapack_ext_common.rb', line 40

def dot(right_v)
  if (right_v.is_a?(NMatrix) && self.stype == :dense && right_v.stype == :dense &&
      self.dim == 2 && right_v.dim == 2 && self.shape[1] == right_v.shape[0])

    result_dtype = NMatrix.upcast(self.dtype,right_v.dtype)
    left = self.dtype == result_dtype ? self : self.cast(dtype: result_dtype)
    right = right_v.dtype == result_dtype ? right_v : right_v.cast(dtype: result_dtype)

    left = left.clone if left.is_ref?
    right = right.clone if right.is_ref?

    result_m = left.shape[0]
    result_n = right.shape[1]
    left_n = left.shape[1]
    vector = result_n == 1
    result = NMatrix.new([result_m,result_n], dtype: result_dtype)

    if vector
      NMatrix::BLAS.cblas_gemv(false, result_m, left_n, 1, left, left_n, right, 1, 0, result, 1)
    else
      NMatrix::BLAS.cblas_gemm(:row, false, false, result_m, result_n, left_n, 1, left, left_n, right, result_n, 0, result, result_n)
    end
    return result
  else
    #internal_dot will handle non-dense matrices (and also dot-products for NMatrix's with dim=1),
    #and also all error-handling if the input is not valid
    self.internal_dot(right_v)
  end
end

#geqrf! ⇒ Object

call-seq:

geqrf! -> shape.min x 1 NMatrix

QR factorization of a general M-by-N matrix A.

The QR factorization is A = QR, where Q is orthogonal and R is Upper Triangular A is overwritten with the elements of R and Q with Q being represented by the elements below A’s diagonal and an array of scalar factors in the output NMatrix.

The matrix Q is represented as a product of elementary reflectors

Q = H(1) H(2) . . . H(k), where k = min(m,n).

Each H(i) has the form

H(i) = I - tau * v * v'

www.netlib.org/lapack/explore-html/d3/d69/dgeqrf_8f.html

Only works for dense matrices.

  • Returns :

    • Vector TAU. Q and R are stored in A. Q is represented by TAU and A

  • Raises :

    • StorageTypeError -> LAPACK functions only work on dense matrices.

Raises:

  • (StorageTypeError) 


262
263
264
265
266
267
268
269
 
# File 'lib/nmatrix/lapacke.rb', line 262

def geqrf!
  raise(StorageTypeError, "LAPACK functions only work on dense matrices") unless self.dense?
  
  tau = NMatrix.new([self.shape.min,1], dtype: self.dtype)
  NMatrix::LAPACK::lapacke_geqrf(:row, self.shape[0], self.shape[1], self, self.shape[1], tau)
  
  tau
end

#getrf! ⇒ Object

Raises:

  • (StorageTypeError) 


169
170
171
172
173
174
175
 
# File 'lib/nmatrix/lapacke.rb', line 169

def getrf!
  raise(StorageTypeError, "LAPACK functions only work on dense matrices") unless self.dense?

  ipiv = NMatrix::LAPACK::lapacke_getrf(:row, self.shape[0], self.shape[1], self, self.shape[1])

  return ipiv
end

#internal_dot ⇒ Object 



38
 
# File 'lib/nmatrix/lapack_ext_common.rb', line 38

alias_method :internal_dot, :dot

#invert! ⇒ Object

Raises:

  • (StorageTypeError) 


177
178
179
180
181
182
183
184
185
186
187
188
189
 
# File 'lib/nmatrix/lapacke.rb', line 177

def invert!
  raise(StorageTypeError, "invert only works on dense matrices currently") unless self.dense?
  raise(ShapeError, "Cannot invert non-square matrix") unless shape[0] == shape[1]
  raise(DataTypeError, "Cannot invert an integer matrix in-place") if self.integer_dtype?

  # Get the pivot array; factor the matrix
  n = self.shape[0]
  pivot = NMatrix::LAPACK::lapacke_getrf(:row, n, n, self, n)
  # Now calculate the inverse using the pivot array
  NMatrix::LAPACK::lapacke_getri(:row, n, self, n, pivot)

  self
end

#ormqr(tau, side = :left, transpose = false, c = nil) ⇒ Object

call-seq:

ormqr(tau) -> NMatrix
ormqr(tau, side, transpose, c) -> NMatrix

Returns the product Q * c or c * Q after a call to geqrf! used in QR factorization. c is overwritten with the elements of the result NMatrix if supplied. Q is the orthogonal matrix represented by tau and the calling NMatrix

Only works on float types, use unmqr for complex types.

Arguments

  • tau - vector containing scalar factors of elementary reflectors

  • side - direction of multiplication [:left, :right]

  • transpose - apply Q with or without transpose [false, :transpose]

  • c - NMatrix multplication argument that is overwritten, no argument assumes c = identity

  • Returns :

    • Q * c or c * Q Where Q may be transposed before multiplication.

  • Raises :

    • StorageTypeError -> LAPACK functions only work on dense matrices.

    • TypeError -> Works only on floating point matrices, use unmqr for complex types

    • TypeError -> c must have the same dtype as the calling NMatrix

Raises:

  • (StorageTypeError) 


299
300
301
302
303
304
305
306
307
308
309
310
 
# File 'lib/nmatrix/lapacke.rb', line 299

def ormqr(tau, side=:left, transpose=false, c=nil)
  raise(StorageTypeError, "LAPACK functions only work on dense matrices") unless self.dense?
  raise(TypeError, "Works only on floating point matrices, use unmqr for complex types") if self.complex_dtype?
  raise(TypeError, "c must have the same dtype as the calling NMatrix") if c and c.dtype != self.dtype


  #Default behaviour produces Q * I  = Q if c is not supplied.
  result = c ? c.clone : NMatrix.identity(self.shape[0], dtype: self.dtype)
  NMatrix::LAPACK::lapacke_ormqr(:row, side, transpose, result.shape[0], result.shape[1], tau.shape[0], self, self.shape[1], tau, result, result.shape[1])
  
  result
end

#potrf!(which) ⇒ Object

Raises:

  • (StorageTypeError) 


191
192
193
194
195
196
 
# File 'lib/nmatrix/lapacke.rb', line 191

def potrf!(which)
  raise(StorageTypeError, "LAPACK functions only work on dense matrices") unless self.dense?
  raise(ShapeError, "Cholesky decomposition only valid for square matrices") unless self.dim == 2 && self.shape[0] == self.shape[1]

  NMatrix::LAPACK::lapacke_potrf(:row, which, self.shape[0], self, self.shape[1])
end

#solve(b, opts = {}) ⇒ Object

Raises:

  • (ShapeError) 


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
 
# File 'lib/nmatrix/lapacke.rb', line 198

def solve(b, opts = {})
  raise(ShapeError, "Must be called on square matrix") unless self.dim == 2 && self.shape[0] == self.shape[1]
  raise(ShapeError, "number of rows of b must equal number of cols of self") if 
    self.shape[1] != b.shape[0]
  raise(ArgumentError, "only works with dense matrices") if self.stype != :dense
  raise(ArgumentError, "only works for non-integer, non-object dtypes") if 
    integer_dtype? or object_dtype? or b.integer_dtype? or b.object_dtype?

  opts = { form: :general }.merge(opts)
  x    = b.clone
  n    = self.shape[0]
  nrhs = b.shape[1]

  case opts[:form] 
  when :general
    clone = self.clone
    ipiv = NMatrix::LAPACK.lapacke_getrf(:row, n, n, clone, n)
    NMatrix::LAPACK.lapacke_getrs(:row, :no_transpose, n, nrhs, clone, n, ipiv, x, nrhs)
    x
  when :upper_tri, :upper_triangular
    raise(ArgumentError, "upper triangular solver does not work with complex dtypes") if
      complex_dtype? or b.complex_dtype?
    NMatrix::BLAS::cblas_trsm(:row, :left, :upper, false, :nounit, n, nrhs, 1.0, self, n, x, nrhs)
    x
  when :lower_tri, :lower_triangular
    raise(ArgumentError, "lower triangular solver does not work with complex dtypes") if
      complex_dtype? or b.complex_dtype?
    NMatrix::BLAS::cblas_trsm(:row, :left, :lower, false, :nounit, n, nrhs, 1.0, self, n, x, nrhs)
    x
  when :pos_def, :positive_definite
    u, l = self.factorize_cholesky
    z = l.solve(b, form: :lower_tri)
    u.solve(z, form: :upper_tri)
  else
    raise(ArgumentError, "#{opts[:form]} is not a valid form option")
  end
end

#unmqr(tau, side = :left, transpose = false, c = nil) ⇒ Object

call-seq:

unmqr(tau) -> NMatrix
unmqr(tau, side, transpose, c) -> NMatrix

Returns the product Q * c or c * Q after a call to geqrf! used in QR factorization. c is overwritten with the elements of the result NMatrix if it is supplied. Q is the orthogonal matrix represented by tau and the calling NMatrix

Only works on complex types, use ormqr for float types.

Arguments

  • tau - vector containing scalar factors of elementary reflectors

  • side - direction of multiplication [:left, :right]

  • transpose - apply Q as Q or its complex conjugate [false, :complex_conjugate]

  • c - NMatrix multplication argument that is overwritten, no argument assumes c = identity

  • Returns :

    • Q * c or c * Q Where Q may be transformed to its complex conjugate before multiplication.

  • Raises :

    • StorageTypeError -> LAPACK functions only work on dense matrices.

    • TypeError -> Works only on floating point matrices, use unmqr for complex types

    • TypeError -> c must have the same dtype as the calling NMatrix

Raises:

  • (StorageTypeError) 


340
341
342
343
344
345
346
347
348
349
350
 
# File 'lib/nmatrix/lapacke.rb', line 340

def unmqr(tau, side=:left, transpose=false, c=nil)
  raise(StorageTypeError, "ATLAS functions only work on dense matrices") unless self.dense?
  raise(TypeError, "Works only on complex matrices, use ormqr for normal floating point matrices") unless self.complex_dtype?
  raise(TypeError, "c must have the same dtype as the calling NMatrix") if c and c.dtype != self.dtype

  #Default behaviour produces Q * I  = Q if c is not supplied.
  result = c ? c.clone : NMatrix.identity(self.shape[0], dtype: self.dtype)
  NMatrix::LAPACK::lapacke_unmqr(:row, side, transpose, result.shape[0], result.shape[1], tau.shape[0], self, self.shape[1], tau, result, result.shape[1])
  
  result
end