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

Class Method Details

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

Raises:

  • (StorageTypeError) 


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

Raises:

  • (ShapeError) 


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