CHOLMOD makeover

base/linalg/diagonal.jl

@@ -99,27 +99,28 @@ expm(D::Diagonal) = Diagonal(exp(D.diag))
 sqrtm(D::Diagonal) = Diagonal(sqrt(D.diag))
 
 #Linear solver
-function \{TD<:Number,TA<:Number}(D::Diagonal{TD}, A::AbstractArray{TA,1})
-    m, n = size(A,2)==1 ? (size(A,1),1) : size(A)
-    m==length(D.diag) || throw(DimensionMismatch())
-    (m == 0 || n == 0) && return A
-    C = Array(typeof(one(TD)/one(TA)),size(A))
+function A_ldiv_B!(D::Diagonal, B::StridedVecOrMat)
+    m, n = size(B, 1), size(B, 2)
+    m == length(D.diag) || throw(DimensionMismatch("diagonal matrix is $(length(D.diag)) by $(length(D.diag)) but right hand side has $m rows"))
+    (m == 0 || n == 0) && return B
     for j = 1:n
         for i = 1:m
             di = D.diag[i]
-            di==0 && throw(SingularException(i))
-            C[i,j] = A[i,j] / di
+            di == 0 && throw(SingularException(i))
+            B[i,j] /= di
         end
     end
-    return C
+    return B
 end
+\(D::Diagonal, B::StridedMatrix) = scale(1 ./ D.diag, B)
+\(D::Diagonal, b::StridedVector) = reshape(scale(1 ./ D.diag, reshape(b, length(b), 1)), length(b))
 \(Da::Diagonal, Db::Diagonal) = Diagonal(Db.diag ./ Da.diag)
 
 function inv{T}(D::Diagonal{T})
     Di = similar(D.diag)
     for i = 1:length(D.diag)
-        D.diag[i]==zero(T) && throw(SingularException(i))
-        Di[i]=inv(D.diag[i])
+        D.diag[i] == zero(T) && throw(SingularException(i))
+        Di[i] = inv(D.diag[i])
     end
     Diagonal(Di)
 end

base/linalg/generic.jl

@@ -240,7 +240,7 @@ end
 function \{TA,TB}(A::AbstractMatrix{TA}, B::AbstractVecOrMat{TB})
     TC = typeof(one(TA)/one(TB))
     size(A,1) == size(B,1) || throw(DimensionMismatch("LHS and RHS should have the same number of rows. LHS has $(size(A,1)) rows, but RHS has $(size(B,1)) rows."))
-    A_ldiv_B!(factorize(TA == TC ? copy(A) : convert(AbstractMatrix{TC}, A)), TB == TC ? copy(B) : convert(AbstractArray{TC}, B))
+    \(factorize(TA == TC ? A : convert(AbstractMatrix{TC}, A)), TB == TC ? copy(B) : convert(AbstractArray{TC}, B))
 end
 \(a::AbstractVector, b::AbstractArray) = reshape(a, length(a), 1) \ b
 /(A::AbstractVecOrMat, B::AbstractVecOrMat) = (B' \ A')'

base/linalg/symmetric.jl

@@ -13,6 +13,7 @@ typealias HermOrSym{T,S} Union(Hermitian{T,S}, Symmetric{T,S})
 typealias RealHermSymComplexHerm{T<:Real,S} Union(Hermitian{T,S}, Symmetric{T,S}, Hermitian{Complex{T},S})
 
 size(A::HermOrSym, args...) = size(A.data, args...)
+getindex(A::HermOrSym, i::Integer) = ((q, r) = divrem(i - 1, size(A, 1)); A[r + 1, q + 1])
 getindex(A::Symmetric, i::Integer, j::Integer) = (A.uplo == 'U') == (i < j) ? getindex(A.data, i, j) : getindex(A.data, j, i)
 getindex(A::Hermitian, i::Integer, j::Integer) = (A.uplo == 'U') == (i < j) ? getindex(A.data, i, j) : conj(getindex(A.data, j, i))
 full(A::Symmetric) = copytri!(copy(A.data), A.uplo)
@@ -45,7 +46,7 @@ A_mul_B!{T<:BlasComplex,S<:AbstractMatrix}(y::StridedMatrix{T}, A::Hermitian{T,S
 *(A::StridedMatrix, B::HermOrSym) = A*full(B)
 
 factorize(A::HermOrSym) = bkfact(A.data, symbol(A.uplo), issym(A))
-\(A::HermOrSym, B::StridedVecOrMat) = \(bkfact(A.data, symbol(A.uplo), issym(A)), B)
+\{T,S<:StridedMatrix}(A::HermOrSym{T,S}, B::StridedVecOrMat) = \(bkfact(A.data, symbol(A.uplo), issym(A)), B)
 inv{T<:BlasFloat,S<:StridedMatrix}(A::Hermitian{T,S}) = Hermitian{T,S}(inv(bkfact(A.data, symbol(A.uplo))), A.uplo)
 inv{T<:BlasFloat,S<:StridedMatrix}(A::Symmetric{T,S}) = Symmetric{T,S}(inv(bkfact(A.data, symbol(A.uplo), true)), A.uplo)
 

base/sparse.jl

@@ -1,21 +1,18 @@
-include("sparse/abstractsparse.jl")
-
 module SparseMatrix
 
+using Base: NonTupleType
+using Base.Sort: Forward
+using Base.LinAlg: AbstractTriangular
+
 importall Base
 importall Base.LinAlg
 
-import Base.NonTupleType, Base.float, Base.Order, Base.Sort.Forward
-import Base.transpose!, Base.ctranspose!
-import Base.LinAlg.AbstractTriangular
-
 export SparseMatrixCSC,
-       blkdiag, dense, diag, diagm, droptol!, dropzeros!, etree, full,
-       getindex, ishermitian, issparse, issym, istril, istriu, nnz,
-       setindex!, sparse, sparsevec, spdiagm, speye, spones,
-       sprand, sprandbool, sprandn, spzeros, symperm, trace, tril, tril!,
-       triu, triu!, nonzeros, rowvals, nzrange
+       blkdiag, dense, droptol!, dropzeros!, etree, issparse, nnz, nonzeros, nzrange,
+       rowvals, sparse, sparsevec, spdiagm, speye, spones, sprand, sprandbool, sprandn,
+       spzeros, symperm
 
+include("sparse/abstractsparse.jl")
 include("sparse/sparsematrix.jl")
 include("sparse/csparse.jl")