@@ -18,12 +18,35 @@ immutable CholeskyPivoted{T,S<:AbstractMatrix} <: Factorization{T}
18
18
end
19
19
CholeskyPivoted {T} (UL:: AbstractMatrix{T} , uplo:: Char , piv:: Vector{BlasInt} , rank:: BlasInt , tol:: Real , info:: BlasInt ) = CholeskyPivoted {T,typeof(UL)} (UL, uplo, piv, rank, tol, info)
20
20
21
- function chol! {T<:BlasFloat} (A:: StridedMatrix{T} , uplo:: Symbol = :U )
21
+ function chol! {T<:BlasFloat} (A:: StridedMatrix{T} )
22
+ C, info = LAPACK. potrf! (' U' , A)
23
+ return @assertposdef UpperTriangular (C) info
24
+ end
25
+ function chol! {T<:BlasFloat} (A:: StridedMatrix{T} , uplo:: Symbol )
22
26
C, info = LAPACK. potrf! (char_uplo (uplo), A)
23
- return @assertposdef Triangular {eltype (C),typeof(C),uplo,false} (C) info
27
+ return @assertposdef uplo == :U ? UpperTriangular (C) : LowerTriangular (C) info
24
28
end
25
29
26
- function chol! {T} (A:: AbstractMatrix{T} , uplo:: Symbol = :U )
30
+ function chol! {T} (A:: AbstractMatrix{T} )
31
+ n = chksquare (A)
32
+ @inbounds begin
33
+ for k = 1 : n
34
+ for i = 1 : k - 1
35
+ A[k,k] -= A[i,k]' A[i,k]
36
+ end
37
+ A[k,k] = chol! (A[k,k], uplo)
38
+ AkkInv = inv (A[k,k])
39
+ for j = k + 1 : n
40
+ for i = 1 : k - 1
41
+ A[k,j] -= A[i,k]' A[i,j]
42
+ end
43
+ A[k,j] = A[k,k]' \ A[k,j]
44
+ end
45
+ end
46
+ end
47
+ return UpperTriangular (A)
48
+ end
49
+ function chol! {T} (A:: AbstractMatrix{T} , uplo:: Symbol )
27
50
n = chksquare (A)
28
51
@inbounds begin
29
52
if uplo == :L
@@ -58,7 +81,7 @@ function chol!{T}(A::AbstractMatrix{T}, uplo::Symbol=:U)
58
81
throw (ArgumentError (" uplo must be either :U or :L but was $(uplo) " ))
59
82
end
60
83
end
61
- return Triangular (A, uplo, false )
84
+ return uplo == :U ? UpperTriangular (A) : LowerTriangular (A )
62
85
end
63
86
64
87
function cholfact! {T<:BlasFloat} (A:: StridedMatrix{T} , uplo:: Symbol = :U ; pivot= false , tol= 0.0 )
@@ -118,14 +141,14 @@ size(C::Union(Cholesky, CholeskyPivoted)) = size(C.UL)
118
141
size (C:: Union(Cholesky, CholeskyPivoted) , d:: Integer ) = size (C. UL,d)
119
142
120
143
function getindex {T,S,UpLo} (C:: Cholesky{T,S,UpLo} , d:: Symbol )
121
- d == :U && return Triangular (UpLo == d ? C. UL : C. UL' , :U )
122
- d == :L && return Triangular (UpLo == d ? C. UL : C. UL' , :L )
123
- d == :UL && return Triangular (C. UL, UpLo )
144
+ d == :U && return UpperTriangular (UpLo == d ? C. UL : C. UL' )
145
+ d == :L && return LowerTriangular (UpLo == d ? C. UL : C. UL' )
146
+ d == :UL && return UpLo == :U ? UpperTriangular (C. UL) : LowerTriangular (C . UL )
124
147
throw (KeyError (d))
125
148
end
126
149
function getindex {T<:BlasFloat} (C:: CholeskyPivoted{T} , d:: Symbol )
127
- d == :U && return Triangular (symbol (C. uplo) == d ? C. UL : C. UL' , :U )
128
- d == :L && return Triangular (symbol (C. uplo) == d ? C. UL : C. UL' , :L )
150
+ d == :U && return UpperTriangular (symbol (C. uplo) == d ? C. UL : C. UL' )
151
+ d == :L && return LowerTriangular (symbol (C. uplo) == d ? C. UL : C. UL' )
129
152
d == :p && return C. piv
130
153
if d == :P
131
154
n = size (C, 1 )
@@ -142,8 +165,8 @@ show{T,S<:AbstractMatrix,UpLo}(io::IO, C::Cholesky{T,S,UpLo}) = (println("$(type
142
165
143
166
A_ldiv_B! {T<:BlasFloat,S<:AbstractMatrix} (C:: Cholesky{T,S,:U} , B:: StridedVecOrMat{T} ) = LAPACK. potrs! (' U' , C. UL, B)
144
167
A_ldiv_B! {T<:BlasFloat,S<:AbstractMatrix} (C:: Cholesky{T,S,:L} , B:: StridedVecOrMat{T} ) = LAPACK. potrs! (' L' , C. UL, B)
145
- A_ldiv_B! {T,S<:AbstractMatrix} (C:: Cholesky{T,S,:L} , B:: StridedVecOrMat ) = Ac_ldiv_B! (Triangular (C. UL, :L , false ), A_ldiv_B! (Triangular (C. UL, :L , false ), B))
146
- A_ldiv_B! {T,S<:AbstractMatrix} (C:: Cholesky{T,S,:U} , B:: StridedVecOrMat ) = A_ldiv_B! (Triangular (C. UL, :U , false ), Ac_ldiv_B! (Triangular (C. UL, :U , false ), B))
168
+ A_ldiv_B! {T,S<:AbstractMatrix} (C:: Cholesky{T,S,:L} , B:: StridedVecOrMat ) = Ac_ldiv_B! (LowerTriangular (C. UL), A_ldiv_B! (LowerTriangular (C. UL), B))
169
+ A_ldiv_B! {T,S<:AbstractMatrix} (C:: Cholesky{T,S,:U} , B:: StridedVecOrMat ) = A_ldiv_B! (UpperTriangular (C. UL), Ac_ldiv_B! (UpperTriangular (C. UL), B))
147
170
148
171
function A_ldiv_B! {T<:BlasFloat} (C:: CholeskyPivoted{T} , B:: StridedVector{T} )
149
172
chkfullrank (C)
@@ -161,8 +184,8 @@ function A_ldiv_B!{T<:BlasFloat}(C::CholeskyPivoted{T}, B::StridedMatrix{T})
161
184
end
162
185
B
163
186
end
164
- A_ldiv_B! (C:: CholeskyPivoted , B:: StridedVector ) = C. uplo== ' L' ? Ac_ldiv_B! (Triangular (C. UL, symbol (C . uplo ), false ), A_ldiv_B! (Triangular (C. UL, symbol (C . uplo ), false ), B[C. piv]))[invperm (C. piv)] : A_ldiv_B! (Triangular (C. UL, symbol (C . uplo ), false ), Ac_ldiv_B! (Triangular (C. UL, symbol (C . uplo), false ), B[C. piv]))[invperm (C. piv)]
165
- A_ldiv_B! (C:: CholeskyPivoted , B:: StridedMatrix ) = C. uplo== ' L' ? Ac_ldiv_B! (Triangular (C. UL, symbol (C . uplo ), false ), A_ldiv_B! (Triangular (C. UL, symbol (C . uplo ), false ), B[C. piv,:]))[invperm (C. piv),:] : A_ldiv_B! (Triangular (C. UL, symbol (C . uplo ), false ), Ac_ldiv_B! (Triangular (C. UL, symbol (C . uplo), false ), B[C. piv,:]))[invperm (C. piv),:]
187
+ A_ldiv_B! (C:: CholeskyPivoted , B:: StridedVector ) = C. uplo== ' L' ? Ac_ldiv_B! (LowerTriangular (C. UL), A_ldiv_B! (LowerTriangular (C. UL), B[C. piv]))[invperm (C. piv)] : A_ldiv_B! (UpperTriangular (C. UL), Ac_ldiv_B! (UpperTriangular (C. UL), B[C. piv]))[invperm (C. piv)]
188
+ A_ldiv_B! (C:: CholeskyPivoted , B:: StridedMatrix ) = C. uplo== ' L' ? Ac_ldiv_B! (LowerTriangular (C. UL), A_ldiv_B! (LowerTriangular (C. UL), B[C. piv,:]))[invperm (C. piv),:] : A_ldiv_B! (UpperTriangular (C. UL), Ac_ldiv_B! (UpperTriangular (C. UL), B[C. piv,:]))[invperm (C. piv),:]
166
189
167
190
function det {T,S,UpLo} (C:: Cholesky{T,S,UpLo} )
168
191
dd = one (T)
0 commit comments