27
27
# nice ruby interfaces for ATLAS functions.
28
28
#++
29
29
30
- require ' nmatrix/nmatrix.rb'
31
- # need to have nmatrix required first or else bad things will happen
32
- require_relative ' lapack_ext_common'
30
+ require " nmatrix/nmatrix.rb"
31
+ # need to have nmatrix required first or else bad things will happen
32
+ require_relative " lapack_ext_common"
33
33
34
34
NMatrix . register_lapack_extension ( "nmatrix-atlas" )
35
35
36
36
require "nmatrix_atlas.so"
37
37
38
38
class NMatrix
39
-
40
- #Add functions from the ATLAS C extension to the main LAPACK and BLAS modules.
41
- #This will overwrite the original functions where applicable.
39
+ # Add functions from the ATLAS C extension to the main LAPACK and BLAS modules.
40
+ # This will overwrite the original functions where applicable.
42
41
module LAPACK
43
42
class << self
44
43
NMatrix ::ATLAS ::LAPACK . singleton_methods . each do |m |
@@ -68,8 +67,8 @@ def posv(uplo, a, b)
68
67
unless a . stype == :dense && b . stype == :dense
69
68
70
69
raise ( DataTypeError , "only works for non-integer, non-object dtypes" ) \
71
- if a . integer_dtype? || a . object_dtype? || \
72
- b . integer_dtype? || b . object_dtype?
70
+ if a . integer_dtype? || a . object_dtype? || \
71
+ b . integer_dtype? || b . object_dtype?
73
72
74
73
x = b . clone
75
74
clone = a . clone
@@ -83,78 +82,81 @@ def posv(uplo, a, b)
83
82
x . transpose
84
83
end
85
84
86
- def geev ( matrix , which = :both )
85
+ def geev ( matrix , which = :both )
87
86
raise ( StorageTypeError , "LAPACK functions only work on dense matrices" ) \
88
87
unless matrix . dense?
89
88
90
89
raise ( ShapeError , "eigenvalues can only be computed for square matrices" ) \
91
90
unless matrix . dim == 2 && matrix . shape [ 0 ] == matrix . shape [ 1 ]
92
91
93
- jobvl = ( which == :both || which == :left ) ? :t : false
94
- jobvr = ( which == :both || which == :right ) ? :t : false
92
+ jobvl = which == :both || which == :left ? :t : false
93
+ jobvr = which == :both || which == :right ? :t : false
95
94
96
95
n = matrix . shape [ 0 ]
97
96
98
97
# Outputs
99
98
eigenvalues = NMatrix . new ( [ n , 1 ] , dtype : matrix . dtype )
100
- # For real dtypes this holds only the real part of the eigenvalues.
99
+ # For real dtypes this holds only the real part of the eigenvalues.
101
100
imag_eigenvalues = matrix . complex_dtype? ? nil : NMatrix . new ( [ n , 1 ] , \
102
- dtype : matrix . dtype ) # For complex dtypes, this is unused.
101
+ dtype : matrix . dtype ) # For complex dtypes, this is unused.
103
102
left_output = jobvl ? matrix . clone_structure : nil
104
103
right_output = jobvr ? matrix . clone_structure : nil
105
104
106
105
# lapack_geev is a pure LAPACK routine so it expects column-major matrices,
107
106
# so we need to transpose the input as well as the output.
108
107
temporary_matrix = matrix . transpose
109
- NMatrix ::LAPACK :: lapack_geev ( jobvl , # compute left eigenvectors of A?
110
- jobvr , # compute right eigenvectors of A? (left eigenvectors of A**T)
111
- n , # order of the matrix
112
- temporary_matrix , # input matrix (used as work)
113
- n , # leading dimension of matrix
114
- eigenvalues , # real part of computed eigenvalues
115
- imag_eigenvalues , # imag part of computed eigenvalues
116
- left_output , # left eigenvectors, if applicable
117
- n , # leading dimension of left_output
118
- right_output , # right eigenvectors, if applicable
119
- n , # leading dimension of right_output
120
- 2 * n )
108
+ NMatrix ::LAPACK . lapack_geev ( jobvl , # compute left eigenvectors of A?
109
+ jobvr , # compute right eigenvectors of A? (left eigenvectors of A**T)
110
+ n , # order of the matrix
111
+ temporary_matrix , # input matrix (used as work)
112
+ n , # leading dimension of matrix
113
+ eigenvalues , # real part of computed eigenvalues
114
+ imag_eigenvalues , # imag part of computed eigenvalues
115
+ left_output , # left eigenvectors, if applicable
116
+ n , # leading dimension of left_output
117
+ right_output , # right eigenvectors, if applicable
118
+ n , # leading dimension of right_output
119
+ 2 * n )
121
120
left_output = left_output . transpose if jobvl
122
121
right_output = right_output . transpose if jobvr
123
122
124
-
125
123
# For real dtypes, transform left_output and right_output into correct forms.
126
124
# If the j'th and the (j+1)'th eigenvalues form a complex conjugate
127
125
# pair, then the j'th and (j+1)'th columns of the matrix are
128
126
# the real and imag parts of the eigenvector corresponding
129
127
# to the j'th eigenvalue.
130
- if ! matrix . complex_dtype?
128
+ unless matrix . complex_dtype?
131
129
complex_indices = [ ]
132
130
n . times do |i |
133
131
complex_indices << i if imag_eigenvalues [ i ] != 0.0
134
132
end
135
133
136
- if ! complex_indices . empty?
134
+ unless complex_indices . empty?
137
135
# For real dtypes, put the real and imaginary parts together
138
- eigenvalues = eigenvalues + imag_eigenvalues * \
139
- Complex ( 0.0 , 1.0 )
140
- left_output = left_output . cast ( dtype : \
141
- NMatrix . upcast ( :complex64 , matrix . dtype ) ) if left_output
142
- right_output = right_output . cast ( dtype : NMatrix . upcast ( :complex64 , \
143
- matrix . dtype ) ) if right_output
136
+ eigenvalues += imag_eigenvalues * \
137
+ Complex ( 0.0 , 1.0 )
138
+ if left_output
139
+ left_output = left_output . cast ( dtype : \
140
+ NMatrix . upcast ( :complex64 , matrix . dtype ) )
141
+ end
142
+ if right_output
143
+ right_output = right_output . cast ( dtype : NMatrix . upcast ( :complex64 , \
144
+ matrix . dtype ) )
145
+ end
144
146
end
145
147
146
148
complex_indices . each_slice ( 2 ) do |i , _ |
147
149
if right_output
148
- right_output [ 0 ...n , i ] = right_output [ 0 ...n , i ] + \
149
- right_output [ 0 ...n , i + 1 ] * Complex ( 0.0 , 1.0 )
150
- right_output [ 0 ...n , i + 1 ] = \
151
- right_output [ 0 ...n , i ] . complex_conjugate
150
+ right_output [ 0 ...n , i ] = right_output [ 0 ...n , i ] + \
151
+ right_output [ 0 ...n , i + 1 ] * Complex ( 0.0 , 1.0 )
152
+ right_output [ 0 ...n , i + 1 ] = \
153
+ right_output [ 0 ...n , i ] . complex_conjugate
152
154
end
153
155
154
156
if left_output
155
- left_output [ 0 ...n , i ] = left_output [ 0 ...n , i ] + \
156
- left_output [ 0 ...n , i + 1 ] * Complex ( 0.0 , 1.0 )
157
- left_output [ 0 ...n , i + 1 ] = left_output [ 0 ...n , i ] . complex_conjugate
157
+ left_output [ 0 ...n , i ] = left_output [ 0 ...n , i ] + \
158
+ left_output [ 0 ...n , i + 1 ] * Complex ( 0.0 , 1.0 )
159
+ left_output [ 0 ...n , i + 1 ] = left_output [ 0 ...n , i ] . complex_conjugate
158
160
end
159
161
end
160
162
end
@@ -168,7 +170,7 @@ def geev(matrix, which=:both)
168
170
end
169
171
end
170
172
171
- def gesvd ( matrix , workspace_size = 1 )
173
+ def gesvd ( matrix , workspace_size = 1 )
172
174
result = alloc_svd_result ( matrix )
173
175
174
176
m = matrix . shape [ 0 ]
@@ -177,16 +179,16 @@ def gesvd(matrix, workspace_size=1)
177
179
# This is a pure LAPACK function so it expects column-major functions.
178
180
# So we need to transpose the input as well as the output.
179
181
matrix = matrix . transpose
180
- NMatrix ::LAPACK :: lapack_gesvd ( :a , :a , m , n , matrix , \
181
- m , result [ 1 ] , result [ 0 ] , m , result [ 2 ] , n , workspace_size )
182
+ NMatrix ::LAPACK . lapack_gesvd ( :a , :a , m , n , matrix , \
183
+ m , result [ 1 ] , result [ 0 ] , m , result [ 2 ] , n , workspace_size )
182
184
result [ 0 ] = result [ 0 ] . transpose
183
185
result [ 2 ] = result [ 2 ] . transpose
184
186
result
185
187
end
186
188
187
- def gesdd ( matrix , workspace_size = nil )
189
+ def gesdd ( matrix , workspace_size = nil )
188
190
min_workspace_size = matrix . shape . min * \
189
- ( 6 + 4 * matrix . shape . min ) + matrix . shape . max
191
+ ( 6 + 4 * matrix . shape . min ) + matrix . shape . max
190
192
workspace_size = min_workspace_size if \
191
193
workspace_size . nil? || workspace_size < min_workspace_size
192
194
@@ -198,8 +200,8 @@ def gesdd(matrix, workspace_size=nil)
198
200
# This is a pure LAPACK function so it expects column-major functions.
199
201
# So we need to transpose the input as well as the output.
200
202
matrix = matrix . transpose
201
- NMatrix ::LAPACK :: lapack_gesdd ( :a , m , n , matrix , m , result [ 1 ] , \
202
- result [ 0 ] , m , result [ 2 ] , n , workspace_size )
203
+ NMatrix ::LAPACK . lapack_gesdd ( :a , m , n , matrix , m , result [ 1 ] , \
204
+ result [ 0 ] , m , result [ 2 ] , n , workspace_size )
203
205
result [ 0 ] = result [ 0 ] . transpose
204
206
result [ 2 ] = result [ 2 ] . transpose
205
207
result
@@ -209,36 +211,36 @@ def gesdd(matrix, workspace_size=nil)
209
211
210
212
def invert!
211
213
raise ( StorageTypeError , "invert only works on dense matrices currently" ) \
212
- unless self . dense?
214
+ unless dense?
213
215
214
216
raise ( ShapeError , "Cannot invert non-square matrix" ) \
215
217
unless shape [ 0 ] == shape [ 1 ]
216
218
217
219
raise ( DataTypeError , "Cannot invert an integer matrix in-place" ) \
218
- if self . integer_dtype?
220
+ if integer_dtype?
219
221
220
222
# Even though we are using the ATLAS plugin, we still might be missing
221
223
# CLAPACK (and thus clapack_getri) if we are on OS X.
222
224
if NMatrix . has_clapack?
223
225
# Get the pivot array; factor the matrix
224
226
# We can't used getrf! here since it doesn't have the clapack behavior,
225
227
# so it doesn't play nicely with clapack_getri
226
- n = self . shape [ 0 ]
227
- pivot = NMatrix ::LAPACK :: clapack_getrf ( :row , n , n , self , n )
228
+ n = shape [ 0 ]
229
+ pivot = NMatrix ::LAPACK . clapack_getrf ( :row , n , n , self , n )
228
230
# Now calculate the inverse using the pivot array
229
- NMatrix ::LAPACK :: clapack_getri ( :row , n , self , n , pivot )
231
+ NMatrix ::LAPACK . clapack_getri ( :row , n , self , n , pivot )
230
232
self
231
233
else
232
- __inverse__ ( self , true )
234
+ __inverse__ ( self , true )
233
235
end
234
236
end
235
237
236
238
def potrf! ( which )
237
239
raise ( StorageTypeError , "ATLAS functions only work on dense matrices" ) \
238
- unless self . dense?
240
+ unless dense?
239
241
raise ( ShapeError , "Cholesky decomposition only valid for square matrices" ) \
240
- unless self . dim == 2 && self . shape [ 0 ] == self . shape [ 1 ]
242
+ unless dim == 2 && shape [ 0 ] == shape [ 1 ]
241
243
242
- NMatrix ::LAPACK :: clapack_potrf ( :row , which , self . shape [ 0 ] , self , self . shape [ 1 ] )
244
+ NMatrix ::LAPACK . clapack_potrf ( :row , which , shape [ 0 ] , self , shape [ 1 ] )
243
245
end
244
246
end
0 commit comments