Skip to content

Commit d19ad7c

Browse files
ViralBShahViralBShah
committed
Merge pull request #2308 from JuliaLang/anj/linalg2
WIP:New design of linalg (factorizations)
2 parents cfb683f + f70cd9f commit d19ad7c

34 files changed

+3929
-2930
lines changed

‎Makefile

+1-1
Original file line numberDiff line numberDiff line change
@@ -32,7 +32,7 @@ $(BUILD)/share/julia/helpdb.jl: doc/helpdb.jl | $(BUILD)/share/julia
3232
@cp $< $@
3333

3434
# use sys.ji if it exists, otherwise run two stages
35-
$(BUILD)/$(JL_PRIVATE_LIBDIR)/sys.ji: VERSION base/*.jl base/pkg/*.jl $(BUILD)/share/julia/helpdb.jl
35+
$(BUILD)/$(JL_PRIVATE_LIBDIR)/sys.ji: VERSION base/*.jl base/pkg/*.jl base/linalg/*.jl $(BUILD)/share/julia/helpdb.jl
3636
@#echo `git rev-parse --short HEAD`-$(OS)-$(ARCH) \(`date +"%Y-%m-%d %H:%M:%S"`\) > COMMIT
3737
$(QUIET_JULIA) cd base && \
3838
(test -f $(BUILD)/$(JL_PRIVATE_LIBDIR)/sys.ji || $(JULIA_EXECUTABLE) -bf sysimg.jl) && $(JULIA_EXECUTABLE) -f sysimg.jl || echo "Note: this error is usually fixed by running 'make clean'. If the error persists, 'make cleanall' may help."

‎base/deprecated.jl

+1
Original file line numberDiff line numberDiff line change
@@ -152,6 +152,7 @@ export PipeString
152152
@deprecate localize localpart
153153
@deprecate expr(hd, a...) Expr(hd, a...)
154154
@deprecate expr(hd, a::Array{Any,1}) Expr(hd, a...)
155+
155156
@deprecate logb exponent
156157
@deprecate ilogb exponent
157158
@deprecate ref_shape index_shape

‎base/exports.jl

+15-4
Original file line numberDiff line numberDiff line change
@@ -3,8 +3,10 @@ export
33
PCRE,
44
FFTW,
55
DSP,
6-
LAPACK,
6+
LinAlg,
77
BLAS,
8+
LAPACK,
9+
ARPACK,
810
LibRandom,
911
Random,
1012
Math,
@@ -108,13 +110,17 @@ export
108110
BunchKaufman,
109111
CholeskyDense,
110112
CholeskyPivotedDense,
113+
Eigen,
114+
GSVDDense,
115+
Hessenberg,
111116
LUDense,
112117
LUTridiagonal,
113118
LDLTTridiagonal,
114119
QRDense,
115120
QRPivotedDense,
116121
SVDDense,
117-
GSVDDense,
122+
Hermitian,
123+
Triangular,
118124
InsertionSort,
119125
QuickSort,
120126
MergeSort,
@@ -234,6 +240,7 @@ export
234240
A_rdiv_Bt,
235241
Ac_ldiv_B,
236242
Ac_ldiv_Bc,
243+
Ac_mul_b_RFP,
237244
Ac_mul_B,
238245
Ac_mul_Bc,
239246
Ac_rdiv_B,
@@ -458,7 +465,6 @@ export
458465
cumsum_kbn,
459466
cummin,
460467
cummax,
461-
diff,
462468
fill,
463469
fill!,
464470
find,
@@ -544,6 +550,7 @@ export
544550
chol,
545551
cholfact,
546552
cholfact!,
553+
cholp,
547554
cholpfact,
548555
cholpfact!,
549556
cond,
@@ -554,8 +561,12 @@ export
554561
diagm,
555562
diagmm,
556563
diagmm!,
564+
diff,
557565
dot,
558566
eig,
567+
eigenfact,
568+
eigenfact!,
569+
eigs,
559570
eigvals,
560571
expm,
561572
sqrtm,
@@ -598,7 +609,7 @@ export
598609
svd,
599610
svdfact!,
600611
svdfact,
601-
svdt,
612+
svds,
602613
svdvals!,
603614
svdvals,
604615
symmetrize!,

‎base/expr.jl

+1
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@ symbol(s::ASCIIString) = symbol(s.data)
55
symbol(s::UTF8String) = symbol(s.data)
66
symbol(a::Array{Uint8,1}) =
77
ccall(:jl_symbol_n, Any, (Ptr{Uint8}, Int32), a, length(a))::Symbol
8+
symbol(x::Char) = symbol(string(x))
89

910
gensym() = ccall(:jl_gensym, Any, ())::Symbol
1011

‎base/linalg/arnoldi.jl

+167
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,167 @@
1+
using ARPACK
2+
3+
## eigs
4+
5+
function eigs{T <: BlasFloat}(A::AbstractMatrix{T}, nev::Integer, evtype::ASCIIString, rvec::Bool)
6+
(m, n) = size(A)
7+
if m != n; error("Input must be square"); end
8+
sym = issym(A)
9+
cmplx = iscomplex(A)
10+
bmat = "I"
11+
12+
# Compute the Ritz values and Ritz vectors
13+
(select, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, rwork) =
14+
aupd_wrapper(T, n, sym, cmplx, bmat, nev, evtype, (x) -> A * x)
15+
16+
# Postprocessing to get eigenvalues and eigenvectors
17+
return eupd_wrapper(T, n, sym, cmplx, bmat, nev, evtype, rvec,
18+
select, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, rwork)
19+
20+
end
21+
22+
eigs(A::AbstractMatrix, nev::Integer, typ::ASCIIString) = eigs(A, nev, which, true)
23+
eigs(A::AbstractMatrix, nev::Integer, rvec::Bool) = eigs(A, nev, "LM", rvec)
24+
eigs(A::AbstractMatrix, rvec::Bool) = eigs(A, 6, "LM", rvec)
25+
eigs(A::AbstractMatrix, nev::Integer) = eigs(A, nev, "LM", true)
26+
eigs(A::AbstractMatrix) = eigs(A, 6, "LM", true)
27+
28+
## svds
29+
30+
# For a dense matrix A is ignored and At is actually A'*A
31+
sarupdate{T}(A::StridedMatrix{T}, At::StridedMatrix{T}, X::StridedVector{T}) = BLAS.symv('U', one(T), At, X)
32+
sarupdate{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, At::SparseMatrixCSC{Tv,Ti}, X::StridedVector{Tv}) = At*(A*X)
33+
34+
function svds{T <: Union(Float64,Float32)}(A::AbstractMatrix{T}, nev::Integer,
35+
which::ASCIIString, rvec::Bool)
36+
37+
(m, n) = size(A)
38+
if m < n error("m = $m, n = $n and only the m >= n case is implemented") end
39+
sym = true
40+
cmplx = false
41+
bmat = "I"
42+
At = isa(A, StridedMatrix) ? BLAS.syrk('U','T',1.0,A) : A'
43+
44+
# Compute the Ritz values and Ritz vectors
45+
(select, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, rwork) =
46+
aupd_wrapper(T, n, sym, cmplx, bmat, nev, which, (x) -> sarupdate(A, At, x))
47+
48+
# Postprocessing to get eigenvalues and eigenvectors
49+
(svals, svecs) = eupd_wrapper(T, n, sym, cmplx, bmat, nev, which, rvec,
50+
select, tol, resid, ncv, v, ldv, iparam, ipntr,
51+
workd, workl, lworkl, rwork)
52+
53+
svals = sqrt(svals)
54+
rvec ? (A*svecs*diagm(1./svals), svals, v.') : svals
55+
end
56+
57+
svds(A::AbstractMatrix, nev::Integer, which::ASCIIString) = svds(A, nev, which, true)
58+
svds(A::AbstractMatrix, nev::Integer, rvec::Bool) = svds(A, nev, "LA", rvec)
59+
svds(A::AbstractMatrix, rvec::Bool) = svds(A, 6, "LA", rvec)
60+
svds(A::AbstractMatrix, nev::Integer) = svds(A, nev, "LA", true)
61+
svds(A::AbstractMatrix) = svds(A, 6, "LA", true)
62+
63+
## aupd and eupd wrappers
64+
65+
function aupd_wrapper(T, n::Integer, sym::Bool, cmplx::Bool, bmat::ASCIIString,
66+
nev::Integer, evtype::ASCIIString, linop::Function)
67+
68+
ncv = min(max(nev*2, 20), n)
69+
70+
bmat = "I"
71+
lworkl = cmplx ? ncv * (3*ncv + 5) : ( lworkl = sym ? ncv * (ncv + 8) : ncv * (3*ncv + 6) )
72+
TR = cmplx ? T.types[1] : T
73+
74+
v = Array(T, n, ncv)
75+
workd = Array(T, 3*n)
76+
workl = Array(T, lworkl)
77+
rwork = cmplx ? Array(TR, ncv) : Array(TR, 0)
78+
resid = Array(T, n)
79+
select = Array(BlasInt, ncv)
80+
iparam = zeros(BlasInt, 11)
81+
ipntr = zeros(BlasInt, 14)
82+
83+
tol = zeros(TR, 1)
84+
ido = zeros(BlasInt, 1)
85+
info = zeros(BlasInt, 1)
86+
87+
iparam[1] = blas_int(1) # ishifts
88+
iparam[3] = blas_int(1000) # maxitr
89+
iparam[7] = blas_int(1) # mode 1
90+
91+
zernm1 = 0:(n-1)
92+
93+
while true
94+
if cmplx
95+
naupd(ido, bmat, n, evtype, nev, tol, resid, ncv, v, n,
96+
iparam, ipntr, workd, workl, lworkl, rwork, info)
97+
elseif sym
98+
saupd(ido, bmat, n, evtype, nev, tol, resid, ncv, v, n,
99+
iparam, ipntr, workd, workl, lworkl, info)
100+
else
101+
naupd(ido, bmat, n, evtype, nev, tol, resid, ncv, v, n,
102+
iparam, ipntr, workd, workl, lworkl, info)
103+
end
104+
if info[1] != 0; error("error code $(info[1]) from ARPACK aupd"); end
105+
if (ido[1] != -1 && ido[1] != 1); break; end
106+
workd[ipntr[2]+zernm1] = linop(getindex(workd, ipntr[1]+zernm1))
107+
end
108+
109+
return (select, tol, resid, ncv, v, n, iparam, ipntr, workd, workl, lworkl, rwork)
110+
end
111+
112+
function eupd_wrapper(T, n::Integer, sym::Bool, cmplx::Bool, bmat::ASCIIString,
113+
nev::Integer, evtype::ASCIIString, rvec::Bool,
114+
select, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, rwork)
115+
116+
howmny = "A"
117+
info = zeros(BlasInt, 1)
118+
119+
if cmplx
120+
121+
d = Array(T, nev+1)
122+
sigma = zeros(T, 1)
123+
workev = Array(T, 2ncv)
124+
neupd(rvec, howmny, select, d, v, ldv, workev, sigma,
125+
bmat, n, evtype, nev, tol, resid, ncv, v, ldv,
126+
iparam, ipntr, workd, workl, lworkl, rwork, info)
127+
if info[1] != 0; error("error code $(info[1]) from ARPACK eupd"); end
128+
return rvec ? (d, v[1:n, 1:nev]) : d
129+
130+
elseif sym
131+
132+
d = Array(T, nev)
133+
sigma = zeros(T, 1)
134+
seupd(rvec, howmny, select, d, v, ldv, sigma,
135+
bmat, n, evtype, nev, tol, resid, ncv, v, ldv,
136+
iparam, ipntr, workd, workl, lworkl, info)
137+
if info[1] != 0; error("error code $(info[1]) from ARPACK eupd"); end
138+
return rvec ? (d, v[1:n, 1:nev]) : d
139+
140+
else
141+
142+
dr = Array(T, nev+1)
143+
di = Array(T, nev+1)
144+
sigmar = zeros(T, 1)
145+
sigmai = zeros(T, 1)
146+
workev = Array(T, 3*ncv)
147+
neupd(rvec, howmny, select, dr, di, v, ldv, sigmar, sigmai,
148+
workev, bmat, n, evtype, nev, tol, resid, ncv, v, ldv,
149+
iparam, ipntr, workd, workl, lworkl, info)
150+
if info[1] != 0; error("error code $(info[1]) from ARPACK eupd"); end
151+
evec = complex(zeros(T, n, nev+1), zeros(T, n, nev+1))
152+
j = 1
153+
while j <= nev
154+
if di[j] == 0.0
155+
evec[:,j] = v[:,j]
156+
else
157+
evec[:,j] = v[:,j] + im*v[:,j+1]
158+
evec[:,j+1] = v[:,j] - im*v[:,j+1]
159+
j += 1
160+
end
161+
j += 1
162+
end
163+
d = complex(dr[1:nev],di[1:nev])
164+
return rvec ? (d, evec[1:n, 1:nev]) : d
165+
end
166+
167+
end

‎base/linalg/arpack.jl

+108
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,108 @@
1+
module ARPACK
2+
3+
const libarpack = "libarpack"
4+
5+
export naupd, neupd, saupd, seupd
6+
7+
import LinAlg.BlasInt
8+
import LinAlg.blas_int
9+
10+
for (T, saupd_name, seupd_name, naupd_name, neupd_name) in
11+
((:Float64, :dsaupd_, :dseupd_, :dnaupd_, :dneupd_),
12+
(:Float32, :ssaupd_, :sseupd_, :snaupd_, :sneupd_))
13+
@eval begin
14+
15+
function naupd(ido, bmat, n, evtype, nev, tol, resid, ncv, v, ldv,
16+
iparam, ipntr, workd::Array{$T}, workl::Array{$T}, lworkl, info)
17+
18+
ccall(($(string(naupd_name)), libarpack), Void,
19+
(Ptr{BlasInt}, Ptr{Uint8}, Ptr{BlasInt}, Ptr{Uint8}, Ptr{BlasInt},
20+
Ptr{$T}, Ptr{$T}, Ptr{BlasInt}, Ptr{$T}, Ptr{BlasInt},
21+
Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$T}, Ptr{$T}, Ptr{BlasInt}, Ptr{BlasInt}),
22+
ido, bmat, &n, evtype, &nev, tol, resid, &ncv, v, &ldv,
23+
iparam, ipntr, workd, workl, &lworkl, info)
24+
end
25+
26+
function neupd(rvec, howmny, select, dr, di, z, ldz, sigmar, sigmai,
27+
workev, bmat, n, evtype, nev, tol, resid, ncv, v, ldv,
28+
iparam, ipntr, workd::Array{$T}, workl::Array{$T}, lworkl, info)
29+
30+
ccall(($(string(neupd_name)), libarpack), Void,
31+
(Ptr{BlasInt}, Ptr{Uint8}, Ptr{BlasInt}, Ptr{$T}, Ptr{$T}, Ptr{$T},
32+
Ptr{BlasInt}, Ptr{$T}, Ptr{$T}, Ptr{$T}, Ptr{Uint8}, Ptr{BlasInt},
33+
Ptr{Uint8}, Ptr{BlasInt}, Ptr{$T}, Ptr{$T}, Ptr{BlasInt}, Ptr{$T},
34+
Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$T}, Ptr{$T},
35+
Ptr{BlasInt}, Ptr{BlasInt}),
36+
&rvec, howmny, select, dr, di, z, &ldz, sigmar, sigmai,
37+
workev, bmat, &n, evtype, &nev, tol, resid, &ncv, v, &ldv,
38+
iparam, ipntr, workd, workl, &lworkl, info)
39+
end
40+
41+
function saupd(ido, bmat, n, which, nev, tol, resid, ncv, v, ldv,
42+
iparam, ipntr, workd::Array{$T}, workl::Array{$T}, lworkl, info)
43+
44+
ccall(($(string(saupd_name)), libarpack), Void,
45+
(Ptr{BlasInt}, Ptr{Uint8}, Ptr{BlasInt}, Ptr{Uint8}, Ptr{BlasInt},
46+
Ptr{$T}, Ptr{$T}, Ptr{BlasInt}, Ptr{$T}, Ptr{BlasInt},
47+
Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$T}, Ptr{$T}, Ptr{BlasInt}, Ptr{BlasInt}),
48+
ido, bmat, &n, which, &nev, tol, resid, &ncv, v, &ldv,
49+
iparam, ipntr, workd, workl, &lworkl, info)
50+
51+
end
52+
53+
function seupd(rvec, howmny, select, d, z, ldz, sigma,
54+
bmat, n, evtype, nev, tol, resid, ncv, v, ldv,
55+
iparam, ipntr, workd::Array{$T}, workl::Array{$T}, lworkl, info)
56+
57+
ccall(($(string(seupd_name)), libarpack), Void,
58+
(Ptr{BlasInt}, Ptr{Uint8}, Ptr{BlasInt}, Ptr{$T}, Ptr{$T}, Ptr{BlasInt}, Ptr{$T},
59+
Ptr{Uint8}, Ptr{BlasInt}, Ptr{Uint8}, Ptr{BlasInt},
60+
Ptr{$T}, Ptr{$T}, Ptr{BlasInt}, Ptr{$T}, Ptr{BlasInt}, Ptr{BlasInt},
61+
Ptr{BlasInt}, Ptr{$T}, Ptr{$T}, Ptr{BlasInt}, Ptr{BlasInt}),
62+
&rvec, howmny, select, d, z, &ldz, sigma,
63+
bmat, &n, evtype, &nev, tol, resid, &ncv, v, &ldv,
64+
iparam, ipntr, workd, workl, &lworkl, info)
65+
end
66+
67+
end
68+
end
69+
70+
for (T, TR, naupd_name, neupd_name) in
71+
((:Complex128, :Float64, :znaupd_, :zneupd_),
72+
(:Complex64, :Float32, :cnaupd_, :cneupd_))
73+
@eval begin
74+
75+
function naupd(ido, bmat, n, evtype, nev, tol, resid, ncv, v, ldv,
76+
iparam, ipntr, workd::Array{$T}, workl::Array{$T}, lworkl,
77+
rwork::Array{$TR}, info)
78+
79+
ccall(($(string(naupd_name)), libarpack), Void,
80+
(Ptr{BlasInt}, Ptr{Uint8}, Ptr{BlasInt}, Ptr{Uint8}, Ptr{BlasInt},
81+
Ptr{$TR}, Ptr{$T}, Ptr{BlasInt}, Ptr{$T}, Ptr{BlasInt},
82+
Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$T}, Ptr{$T}, Ptr{BlasInt},
83+
Ptr{$TR}, Ptr{BlasInt}),
84+
ido, bmat, &n, evtype, &nev, tol, resid, &ncv, v, &ldv,
85+
iparam, ipntr, workd, workl, &lworkl, rwork, info)
86+
87+
end
88+
89+
function neupd(rvec, howmny, select, d, z, ldz, workev, sigma,
90+
bmat, n, evtype, nev, tol, resid, ncv, v, ldv,
91+
iparam, ipntr, workd::Array{$T}, workl::Array{$T}, lworkl,
92+
rwork::Array{$TR}, info)
93+
94+
ccall(($(string(neupd_name)), libarpack), Void,
95+
(Ptr{BlasInt}, Ptr{Uint8}, Ptr{BlasInt}, Ptr{$T}, Ptr{$T}, Ptr{BlasInt},
96+
Ptr{$T}, Ptr{$T}, Ptr{Uint8}, Ptr{BlasInt}, Ptr{Uint8}, Ptr{BlasInt},
97+
Ptr{$TR}, Ptr{$T}, Ptr{BlasInt}, Ptr{$T}, Ptr{BlasInt}, Ptr{BlasInt},
98+
Ptr{BlasInt}, Ptr{$T}, Ptr{$T}, Ptr{BlasInt}, Ptr{$TR}, Ptr{BlasInt}),
99+
&rvec, howmny, select, d, z, &ldz, workev, sigma,
100+
bmat, &n, evtype, &nev, tol, resid, &ncv, v, &ldv,
101+
iparam, ipntr, workd, workl, &lworkl, rwork, info)
102+
103+
end
104+
105+
end
106+
end
107+
108+
end # module ARPACK

0 commit comments

Comments
 (0)