Skip to content

Commit ba6ab0c

Browse files
authored
Merge pull request #30 from Uditgulati/linalg
Implement NumRuby::Linalg
2 parents 5428e46 + 4e4daf2 commit ba6ab0c

2 files changed

Lines changed: 109 additions & 15 deletions

File tree

‎.gitignore‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -3,4 +3,6 @@
33
Gemfile.lock
44
/.yardoc
55
.rake_tasks~
6-
*.so
6+
*.so
7+
.vscode
8+

‎lib/nmatrix/lapack.rb‎

Lines changed: 106 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,20 @@
11
module NumRuby::Linalg
2-
def self.inv(obj)
3-
if obj.is_a?(NMatrix)
4-
return obj.invert
2+
def self.inv(matrix)
3+
if not matrix.is_a?(NMatrix)
4+
raise("Invalid matrix. Not of type NMatrix.")
5+
end
6+
if matrix.dim != 2
7+
raise("Invalid shape of matrix. Should be 2.")
8+
end
9+
if matrix.shape[0] != matrix.shape[1]
10+
raise("Invalid shape. Expected square matrix.")
511
end
12+
m, n = matrix.shape
13+
14+
lu, ipiv = NumRuby::Lapack.getrf(matrix)
15+
inv_a = NumRuby::Lapack.getri(lu, ipiv)
16+
17+
return inv_a
618
end
719

820
def self.dot(lha, rha)
@@ -13,12 +25,22 @@ def self.norm
1325

1426
end
1527

16-
def self.solve
17-
28+
def self.solve(a, b, sym_pos: False, lower: False, assume_a: "gen", transposed: False)
29+
# TODO: implement this and remove NMatrix.solve
1830
end
1931

20-
def self.det
32+
def self.det(matrix)
33+
if not matrix.is_a?(NMatrix)
34+
raise("Invalid matrix. Not of type NMatrix.")
35+
end
36+
if matrix.dim != 2
37+
raise("Invalid shape of matrix. Should be 2.")
38+
end
39+
if matrix.shape[0] != matrix.shape[1]
40+
raise("Invalid shape. Expected square matrix.")
41+
end
2142

43+
return matrix.det
2244
end
2345

2446
def self.least_square
@@ -48,23 +70,47 @@ def self.eigvalsh
4870
# Matrix Decomposition
4971

5072

51-
def self.lu(matrix)
73+
def self.lu(matrix, permute_l: False)
74+
if not matrix.is_a?(NMatrix)
75+
raise("Invalid matrix. Not of type NMatrix.")
76+
end
77+
if matrix.dim != 2
78+
raise("Invalid shape of matrix. Should be 2.")
79+
end
80+
81+
lu, ipiv = NumRuby::Linalg.getrf(matrix)
5282

83+
# TODO: calulate p, l, u
5384
end
5485

5586
def self.lu_factor(matrix)
87+
if not matrix.is_a?(NMatrix)
88+
raise("Invalid matrix. Not of type NMatrix.")
89+
end
90+
if matrix.dim != 2
91+
raise("Invalid shape of matrix. Should be 2.")
92+
end
93+
if matrix.shape[0] != matrix.shape[1]
94+
raise("Invalid shape. Expected square matrix.")
95+
end
5696

97+
lu, ipiv = NumRuby::Linalg.getrf(matrix)
98+
99+
return [lu, ipiv]
57100
end
58101

59-
def self.lu_solve(matrix, rhs_val)
102+
def self.lu_solve(lu, ipiv, b, trans: 0)
103+
if lu.shape[0] != b.shape[0]
104+
raise("Incompatibel dimensions.")
105+
end
60106

107+
x = NumRuby::Lapack.getrs(lu, ipiv, b, trans)
108+
return x
61109
end
62110

63-
# Computes the QR decomposition of matrix.
111+
# Computes the SVD decomposition of matrix.
64112
# Args:
65113
# - input matrix, type: NMatrix
66-
# - mode, type: String
67-
# - pivoting, type: Boolean
68114
def self.svd(matrix)
69115

70116
end
@@ -89,11 +135,24 @@ def self.cholesky_solve(matrix)
89135

90136
end
91137

92-
# Computes the QR decomposition of matrix.
138+
# Computes QR decomposition of a matrix.
139+
#
140+
# Calculates the decomposition A = Q*R where Q is unitary/orthogonal and R is upper triangular.
141+
#
93142
# Args:
94-
# - input matrix, type: NMatrix
143+
# - matrix, type: NMatrix
144+
# Matrix to be decomposed
95145
# - mode, type: String
146+
# Determines what information is to be returned: either both Q and R
147+
# ('full', default), only R ('r') or both Q and R but computed in
148+
# economy-size ('economic', see Notes). The final option 'raw'
149+
# (added in Scipy 0.11) makes the function return two matrices
150+
# (Q, TAU) in the internal format used by LAPACK.
96151
# - pivoting, type: Boolean
152+
# Whether or not factorization should include pivoting for rank-revealing
153+
# qr decomposition. If pivoting, compute the decomposition
154+
# A*P = Q*R as above, but where P is chosen such that the diagonal
155+
# of R is non-increasing.
97156
def self.qr(matrix, mode: "full", pivoting: false)
98157
if not ['full', 'r', 'economic', 'raw'].include?(mode.downcase)
99158
raise("Invalid mode. Should be one of ['full', 'r', 'economic', 'raw']")
@@ -106,9 +165,42 @@ def self.qr(matrix, mode: "full", pivoting: false)
106165
end
107166
m, n = matrix.shape
108167

109-
if pivoting == false
168+
if pivoting == true
169+
qr, tau, jpvt = NumRuby::Lapack.geqp3(matrix)
170+
jpvt -= 1
171+
else
110172
qr, tau = NumRuby::Lapack.geqrf(matrix)
173+
end
174+
175+
# calculate R here for both pivot true & false
176+
177+
if ['economic', 'raw'].include?(mode.downcase) or m < n
178+
r = NumRuby.triu(matrix)
179+
else
180+
r = NumRuby.triu(matrix[0...n, 0...n])
181+
end
182+
183+
if pivoting == true
184+
rj = r, jpvt
185+
else
186+
rj = r
187+
end
188+
189+
if mode == 'r'
190+
return rj
191+
elsif mode == 'raw'
111192
return [qr, tau]
112193
end
194+
195+
if m < n
196+
q = NumRuby::Lapack.orgqr(qr[0...m, 0...m], tau)
197+
elsif mode == 'economic'
198+
q = NumRuby::Lapack.orgqr(qr, tau)
199+
else
200+
# TODO: Implement slice view and set slice
201+
q = NumRuby::Lapack.orgqr(qr, tau)
202+
end
203+
204+
return q, rj
113205
end
114206
end

0 commit comments

Comments
 (0)