#-- # = 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: http://narray.rubyforge.org # # == Copyright Information # # SciRuby is Copyright (c) 2010 - 2016, Ruby Science Foundation # NMatrix is Copyright (c) 2012 - 2016, 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: # # * https://github.com/SciRuby/sciruby/wiki/Contributor-Agreement # # == atlas.rb # # ruby file for the nmatrix-atlas gem. Loads the C extension and defines # nice ruby interfaces for ATLAS functions. #++ require 'nmatrix/nmatrix.rb' #need to have nmatrix required first or else bad things will happen require_relative 'lapack_ext_common' NMatrix.register_lapack_extension("nmatrix-atlas") require "nmatrix_atlas.so" class NMatrix #Add functions from the ATLAS C extension to the main LAPACK and BLAS modules. #This will overwrite the original functions where applicable. module LAPACK class << self NMatrix::ATLAS::LAPACK.singleton_methods.each do |m| define_method m, NMatrix::ATLAS::LAPACK.method(m).to_proc end end end module BLAS class << self NMatrix::ATLAS::BLAS.singleton_methods.each do |m| define_method m, NMatrix::ATLAS::BLAS.method(m).to_proc end end end module LAPACK class << self 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 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 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 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 end end 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? # Even though we are using the ATLAS plugin, we still might be missing # CLAPACK (and thus clapack_getri) if we are on OS X. if NMatrix.has_clapack? # Get the pivot array; factor the matrix # We can't used getrf! here since it doesn't have the clapack behavior, # so it doesn't play nicely with clapack_getri n = self.shape[0] pivot = NMatrix::LAPACK::clapack_getrf(:row, n, n, self, n) # Now calculate the inverse using the pivot array NMatrix::LAPACK::clapack_getri(:row, n, self, n, pivot) self else __inverse__(self,true) end end def potrf!(which) raise(StorageTypeError, "ATLAS 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::clapack_potrf(:row, which, self.shape[0], self, self.shape[1]) end end