forked from JuliaLang/julia
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathspqr.jl
174 lines (148 loc) · 6.78 KB
/
spqr.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
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
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
# This file is a part of Julia. License is MIT: http://julialang.org/license
module SPQR
import Base: \
# ordering options */
const ORDERING_FIXED = Int32(0)
const ORDERING_NATURAL = Int32(1)
const ORDERING_COLAMD = Int32(2)
const ORDERING_GIVEN = Int32(3) # only used for C/C++ interface
const ORDERING_CHOLMOD = Int32(4) # CHOLMOD best-effort (COLAMD, METIS,...)
const ORDERING_AMD = Int32(5) # AMD(A'*A)
const ORDERING_METIS = Int32(6) # metis(A'*A)
const ORDERING_DEFAULT = Int32(7) # SuiteSparseQR default ordering
const ORDERING_BEST = Int32(8) # try COLAMD, AMD, and METIS; pick best
const ORDERING_BESTAMD = Int32(9) # try COLAMD and AMD; pick best#
# Let [m n] = size of the matrix after pruning singletons. The default
# ordering strategy is to use COLAMD if m <= 2*n. Otherwise, AMD(A'A) is
# tried. If there is a high fill-in with AMD then try METIS(A'A) and take
# the best of AMD and METIS. METIS is not tried if it isn't installed.
# tol options
const DEFAULT_TOL = Int32(-2) # if tol <= -2, the default tol is used
const NO_TOL = Int32(-1) # if -2 < tol < 0, then no tol is used
# for qmult, method can be 0,1,2,3:
const QTX = Int32(0)
const QX = Int32(1)
const XQT = Int32(2)
const XQ = Int32(3)
# system can be 0,1,2,3: Given Q*R=A*E from SuiteSparseQR_factorize:
const RX_EQUALS_B = Int32(0) # solve R*X=B or X = R\B
const RETX_EQUALS_B = Int32(1) # solve R*E'*X=B or X = E*(R\B)
const RTX_EQUALS_B = Int32(2) # solve R'*X=B or X = R'\B
const RTX_EQUALS_ETB = Int32(3) # solve R'*X=E'*B or X = R'\(E'*B)
using ..SparseArrays: SparseMatrixCSC
using ..SparseArrays.CHOLMOD: C_Dense, C_Sparse, Dense, ITypes, Sparse, SuiteSparseStruct, VTypes, common
import Base: size
import Base.LinAlg: qrfact
import ..SparseArrays.CHOLMOD: convert, free!
immutable C_Factorization{Tv<:VTypes} <: SuiteSparseStruct
xtype::Cint
factors::Ptr{Tv}
end
type Factorization{Tv<:VTypes} <: Base.LinAlg.Factorization{Tv}
m::Int
n::Int
p::Ptr{C_Factorization{Tv}}
function Factorization(m::Integer, n::Integer, p::Ptr{C_Factorization{Tv}})
if p == C_NULL
throw(ArgumentError("factorization failed for unknown reasons. Please submit a bug report."))
end
new(m, n, p)
end
end
Factorization{Tv<:VTypes}(m::Integer, n::Integer, p::Ptr{C_Factorization{Tv}}) = Factorization{Tv}(m, n, p)
size(F::Factorization) = (F.m, F.n)
function size(F::Factorization, i::Integer)
if i < 1
throw(ArgumentError("dimension must be positive"))
end
if i == 1
return F.m
elseif i == 2
return F.n
end
return 1
end
function free!{Tv<:VTypes}(F::Factorization{Tv})
ccall((:SuiteSparseQR_C_free, :libspqr), Cint,
(Ptr{Ptr{C_Factorization{Tv}}}, Ptr{Void}),
&F.p, common()) == 1
end
function backslash{Tv<:VTypes}(ordering::Integer, tol::Real, A::Sparse{Tv}, B::Dense{Tv})
m, n = size(A)
if m != size(B, 1)
throw(DimensionMismatch("left hand side and right hand side must have same number of rows"))
end
d = Dense(ccall((:SuiteSparseQR_C_backslash, :libspqr), Ptr{C_Dense{Tv}},
(Cint, Cdouble, Ptr{C_Sparse{Tv}}, Ptr{C_Dense{Tv}}, Ptr{Void}),
ordering, tol, get(A.p), get(B.p), common()))
finalizer(d, free!)
d
end
function factorize{Tv<:VTypes}(ordering::Integer, tol::Real, A::Sparse{Tv})
s = unsafe_load(A.p)
if s.stype != 0
throw(ArgumentError("stype must be zero"))
end
f = Factorization(size(A)..., ccall((:SuiteSparseQR_C_factorize, :libspqr), Ptr{C_Factorization{Tv}},
(Cint, Cdouble, Ptr{Sparse{Tv}}, Ptr{Void}),
ordering, tol, get(A.p), common()))
finalizer(f, free!)
f
end
function solve{Tv<:VTypes}(system::Integer, QR::Factorization{Tv}, B::Dense{Tv})
m, n = size(QR)
mB = size(B, 1)
if (system == RX_EQUALS_B || system == RETX_EQUALS_B) && m != mB
throw(DimensionMismatch("number of rows in factorized matrix must equal number of rows in right hand side"))
elseif (system == RTX_EQUALS_ETB || system == RTX_EQUALS_B) && n != mB
throw(DimensionMismatch("number of columns in factorized matrix must equal number of rows in right hand side"))
end
d = Dense(ccall((:SuiteSparseQR_C_solve, :libspqr), Ptr{C_Dense{Tv}},
(Cint, Ptr{C_Factorization{Tv}}, Ptr{C_Dense{Tv}}, Ptr{Void}),
system, get(QR.p), get(B.p), common()))
finalizer(d, free!)
d
end
function qmult{Tv<:VTypes}(method::Integer, QR::Factorization{Tv}, X::Dense{Tv})
mQR, nQR = size(QR)
mX, nX = size(X)
if (method == QTX || method == QX) && mQR != mX
throw(DimensionMismatch("Q matrix size $mQR and dense matrix has $mX rows"))
elseif (method == XQT || method == XQ) && mQR != nX
throw(DimensionMismatch("Q matrix size $mQR and dense matrix has $nX columns"))
end
d = Dense(ccall((:SuiteSparseQR_C_qmult, :libspqr), Ptr{C_Dense{Tv}},
(Cint, Ptr{C_Factorization{Tv}}, Ptr{C_Dense{Tv}}, Ptr{Void}),
method, get(QR.p), get(X.p), common()))
finalizer(d, free!)
d
end
qrfact(A::SparseMatrixCSC, ::Type{Val{true}}) = factorize(ORDERING_DEFAULT, DEFAULT_TOL, Sparse(A, 0))
qrfact(A::SparseMatrixCSC) = qrfact(A, Val{true})
# With a real lhs and complex rhs with the same precision, we can reinterpret
# the complex rhs as a real rhs with twice the number of columns
#
# This definition is similar to the definition in factorization.jl except the we here have to use
# \ instead of A_ldiv_B! because of limitations in SPQR
## Two helper methods
_ret_size(F::Factorization, b::AbstractVector) = (size(F, 2),)
_ret_size(F::Factorization, B::AbstractMatrix) = (size(F, 2), size(B, 2))
function (\)(F::Factorization{Float64}, B::VecOrMat{Complex{Float64}})
# | z | z | reinterpret | x | x | x | x | transpose | x | y | reshape | x | y | x | y |
# | z | z | -> | y | y | y | y | -> | x | y | -> | x | y | x | y |
# | x | y |
# | x | y |
c2r = reshape(transpose(reinterpret(Float64, B, (2, length(B)))), size(B, 1), 2*size(B, 2))
x = F\c2r
# | z | z | reinterpret | x | x | x | x | transpose | x | y | reshape | x | y | x | y |
# | z | z | <- | y | y | y | y | <- | x | y | <- | x | y | x | y |
# | x | y |
# | x | y |
return reinterpret(Complex{Float64}, transpose(reshape(x, (length(x) >> 1), 2)), _ret_size(F, B))
end
function (\){T<:VTypes}(F::Factorization{T}, B::StridedVecOrMat{T})
QtB = qmult(QTX, F, Dense(B))
convert(typeof(B), solve(RETX_EQUALS_B, F, QtB))
end
(\)(F::Factorization, B::StridedVecOrMat) = F\convert(AbstractArray{eltype(F)}, B)
end # module