From df8f86e23e465682011a127698c63dffed446237 Mon Sep 17 00:00:00 2001
From: "Viral B. Shah" <viral@mayin.org>
Date: Sun, 3 Feb 2013 20:23:14 -0500
Subject: [PATCH 01/11] Reinstate all the commits from #2069. Relevant comments
 are in #2062.

---
 base/blas.jl    | 254 +++++++++++++++++++++++-------------------------
 test/Makefile   |   2 +-
 test/blas.jl    |  78 +++++++++++++++
 test/default.jl |   1 +
 test/linalg.jl  | 103 +++++++++-----------
 5 files changed, 252 insertions(+), 186 deletions(-)
 create mode 100644 test/blas.jl

diff --git a/base/blas.jl b/base/blas.jl
index c14970d613c84..9b88042097d31 100644
--- a/base/blas.jl
+++ b/base/blas.jl
@@ -166,18 +166,16 @@ function axpy!{T,Ta<:Number,Ti<:Integer}(alpha::Ta, x::Array{T}, rx::Union(Range
     return axpy!(length(rx), convert(T, alpha), pointer(x)+(first(rx)-1)*sizeof(T), step(rx), pointer(y)+(first(ry)-1)*sizeof(T), step(ry))
 end
 
-
-# SUBROUTINE DSYRK(UPLO,TRANS,N,K,ALPHA,A,LDA,BETA,C,LDC)
-# *     .. Scalar Arguments ..
-#       REAL ALPHA,BETA
-#       INTEGER K,LDA,LDC,N
-#       CHARACTER TRANS,UPLO
-# *     ..
-# *     .. Array Arguments ..
-#       REAL A(LDA,*),C(LDC,*)
 for (fname, elty) in ((:dsyrk_,:Float64), (:ssyrk_,:Float32),
                       (:zsyrk_,:Complex128), (:csyrk_,:Complex64))
    @eval begin
+       # SUBROUTINE DSYRK(UPLO,TRANS,N,K,ALPHA,A,LDA,BETA,C,LDC)
+       # *     .. Scalar Arguments ..
+       #       REAL ALPHA,BETA
+       #       INTEGER K,LDA,LDC,N
+       #       CHARACTER TRANS,UPLO
+       # *     .. Array Arguments ..
+       #       REAL A(LDA,*),C(LDC,*)
        function syrk!(uplo::BlasChar, trans::BlasChar, alpha::($elty), A::StridedVecOrMat{$elty},
                       beta::($elty), C::StridedMatrix{$elty})
            m, n = size(C)
@@ -193,14 +191,9 @@ for (fname, elty) in ((:dsyrk_,:Float64), (:ssyrk_,:Float32),
        end
        function syrk(uplo::BlasChar, trans::BlasChar, alpha::($elty), A::StridedVecOrMat{$elty})
            n = size(A, trans == 'N' ? 1 : 2)
-           k = size(A, trans == 'N' ? 2 : 1)
-           C = Array($elty, (n, n)) 
-           ccall(($(string(fname)),libblas), Void,
-                 (Ptr{Uint8}, Ptr{Uint8}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty},
-                  Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}),
-                 &uplo, &trans, &n, &k, &alpha, A, &stride(A,2), &0., C, &stride(C,2))
-           C
+           syrk!(uplo, trans, alpha, A, zero($elty), Array($elty, (n, n)))
        end
+       syrk(uplo::BlasChar, trans::BlasChar, A::StridedVecOrMat{$elty}) = syrk(uplo, trans, one($elty), A)
    end
 end
 
@@ -214,7 +207,7 @@ end
 #       COMPLEX A(LDA,*),C(LDC,*)
 for (fname, elty) in ((:zherk_,:Complex128), (:cherk_,:Complex64))
    @eval begin
-       function herk!(uplo::BlasChar, trans, alpha::($elty), A::StridedVecOrMat{$elty},
+       function herk!(uplo::BlasChar, trans::BlasChar, alpha::($elty), A::StridedVecOrMat{$elty},
                       beta::($elty), C::StridedMatrix{$elty})
            m, n = size(C)
            if m != n error("syrk!: matrix C must be square") end
@@ -227,16 +220,11 @@ for (fname, elty) in ((:zherk_,:Complex128), (:cherk_,:Complex64))
                  &uplo, &trans, &n, &k, &alpha, A, &stride(A,2), &beta, C, &stride(C,2))
            C
        end
-       function herk(uplo::BlasChar, trans, alpha::($elty), A::StridedVecOrMat{$elty})
+       function herk(uplo::BlasChar, trans::BlasChar, alpha::($elty), A::StridedVecOrMat{$elty})
            n = size(A, trans == 'N' ? 1 : 2)
-           k = size(A, trans == 'N' ? 2 : 1)
-           C = Array($elty, (n, n)) 
-           ccall(($(string(fname)),libblas), Void,
-                 (Ptr{Uint8}, Ptr{Uint8}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty},
-                  Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}),
-                 &uplo, &trans, &n, &k, &alpha, A, &stride(A,2), &0., C, &stride(C,2))
-           C
+           herk!(uplo, trans, alpha, A, zero($elty), Array($elty, (n,n)))
        end
+       herk(uplo::BlasChar, trans::BlasChar, A::StridedVecOrMat{$elty}) = herk(uplo, trans, one($elty), A)
    end
 end
 
@@ -266,16 +254,12 @@ for (fname, elty) in ((:dgbmv_,:Float64), (:sgbmv_,:Float32),
        function gbmv(trans::BlasChar, m::Integer, kl::Integer, ku::Integer,
                      alpha::($elty), A::StridedMatrix{$elty}, x::StridedVector{$elty})
            n = stride(A,2)
-           y = Array($elty, n)
-           ccall(($(string(fname)),libblas), Void,
-                 (Ptr{Uint8}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt},
-                  Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
-                  Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}),
-                 &trans, &m, &n, &kl, &ku, &alpha, A, &stride(A,2),
-                 x, &stride(x,1), &0., y, &1)
-           y
+           gbmv!(trans, m, kl, ku, alpha, A, x, zero($elty), Array($elty, n))
+       end
+       function gbmv(trans::BlasChar, m::Integer, kl::Integer, ku::Integer,
+                     A::StridedMatrix{$elty}, x::StridedVector{$elty})
+           gbmv(trans, m, kl, ku, one($elty), A, x)
        end
-
    end
 end
 
@@ -303,34 +287,42 @@ for (fname, elty) in ((:dsbmv_,:Float64), (:ssbmv_,:Float32),
        function sbmv(uplo::BlasChar, k::Integer, alpha::($elty), A::StridedMatrix{$elty},
                      x::StridedVector{$elty})
            n = size(A,2)
-           y = Array($elty, n)
-           ccall(($(string(fname)),libblas), Void,
-                 (Ptr{Uint8}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt},
-                 Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}),
-                 &uplo, &size(A,2), &k, &alpha, A, &stride(A,2), x, &stride(x,1), &0., y, &1)
-           y
+           sbmv!(uplo, k, alpha, A, x, zero($elty), Array($elty, n))
+       end
+       function sbmv(uplo::BlasChar, k::Integer, A::StridedMatrix{$elty}, x::StridedVector{$elty})
+           sbmv(uplo, k, one($elty), A, x)
        end
    end
 end
 
-# (GE) general matrix-matrix multiplication
-# SUBROUTINE DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
-# *     .. Scalar Arguments ..
-#       DOUBLE PRECISION ALPHA,BETA
-#       INTEGER K,LDA,LDB,LDC,M,N
-#       CHARACTER TRANSA,TRANSB
-# *     .. Array Arguments ..
-#       DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
-for (fname, elty) in ((:dgemm_,:Float64), (:sgemm_,:Float32),
-                      (:zgemm_,:Complex128), (:cgemm_,:Complex64))
+# (GE) general matrix-matrix and matrix-vector multiplication
+for (gemm, gemv, elty) in
+    ((:dgemm_,:dgemv_,:Float64),
+     (:sgemm_,:sgemv_,:Float32),
+     (:zgemm_,:zgemv_,:Complex128),
+     (:cgemm_,:cgemv_,:Complex64))
    @eval begin
-       function gemm!(transA::BlasChar, transB::BlasChar, alpha::($elty), A::StridedMatrix{$elty},
-                      B::StridedMatrix{$elty}, beta::($elty), C::StridedMatrix{$elty})
+       # SUBROUTINE DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
+       # *     .. Scalar Arguments ..
+       #       DOUBLE PRECISION ALPHA,BETA
+       #       INTEGER K,LDA,LDB,LDC,M,N
+       #       CHARACTER TRANSA,TRANSB
+       # *     .. Array Arguments ..
+       #       DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
+       function gemm!(transA::BlasChar, transB::BlasChar,
+                      alpha::($elty), A::StridedMatrix{$elty},
+                      B::StridedMatrix{$elty},
+                      beta::($elty), C::StridedMatrix{$elty})
+#           if any([stride(A,1), stride(B,1), stride(C,1)] .!= 1)
+#               error("gemm!: BLAS module requires contiguous matrix columns")
+#           end  # should this be checked on every call?
            m = size(A, transA == 'N' ? 1 : 2)
            k = size(A, transA == 'N' ? 2 : 1)
            n = size(B, transB == 'N' ? 2 : 1)
-           if m != size(C,1) || n != size(C,2) error("gemm!: mismatched dimensions") end
-           ccall(($(string(fname)),libblas), Void,
+           if m != size(C,1) || n != size(C,2)
+               error("gemm!: mismatched dimensions")
+           end
+           ccall(($(string(gemm)),libblas), Void,
                  (Ptr{Uint8}, Ptr{Uint8}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt},
                   Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
                   Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}),
@@ -338,92 +330,66 @@ for (fname, elty) in ((:dgemm_,:Float64), (:sgemm_,:Float32),
                  B, &stride(B,2), &beta, C, &stride(C,2))
            C
        end
-       function gemm(transA::BlasChar, transB::BlasChar, alpha::($elty), A::StridedMatrix{$elty}, B::StridedMatrix{$elty})
-           m = size(A, transA == 'N' ? 1 : 2)
-           k = size(A, transA == 'N' ? 2 : 1)
-           if k != size(B, transB == 'N' ? 1 : 2) error("gemm!: mismatched dimensions") end
-           n = size(B, transB == 'N' ? 2 : 1)
-           C = Array($elty, (m, n))
-           ccall(($(string(fname)),libblas), Void,
-                 (Ptr{Uint8}, Ptr{Uint8}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt},
-                  Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
-                  Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}),
-                 &transA, &transB, &m, &n, &k, &alpha, A, &stride(A,2),
-                 B, &stride(B,2), &0., C, &stride(C,2))
-           C
+       function gemm(transA::BlasChar, transB::BlasChar,
+                     alpha::($elty), A::StridedMatrix{$elty},
+                     B::StridedMatrix{$elty})
+           gemm!(transA, transB, alpha, A, B, zero($elty),
+                 Array($elty, (size(A, transA == 'N' ? 1 : 2),
+                               size(B, transB == 'N' ? 2 : 1))))
        end
-   end
-end
-
-#SUBROUTINE DGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
-#*     .. Scalar Arguments ..
-#      DOUBLE PRECISION ALPHA,BETA
-#      INTEGER INCX,INCY,LDA,M,N
-#      CHARACTER TRANS
-#*     .. Array Arguments ..
-#      DOUBLE PRECISION A(LDA,*),X(*),Y(*)
-
-for (fname, elty) in ((:dgemv_,:Float64), (:sgemv_,:Float32),
-                      (:zgemv_,:Complex128), (:cgemv_,:Complex64))
-   @eval begin
-       function gemv!(trans::BlasChar, alpha::($elty), A::StridedMatrix{$elty},
-                      X::StridedVector{$elty}, beta::($elty), Y::StridedVector{$elty})
-           ccall(($(string(fname)),libblas), Void,
-                 (Ptr{Uint8}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt},
-                  Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}),
+       function gemm(transA::BlasChar, transB::BlasChar,
+                     A::StridedMatrix{$elty}, B::StridedMatrix{$elty})
+           gemm(transA, transB, one($elty), A, B)
+       end
+       #SUBROUTINE DGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
+       #*     .. Scalar Arguments ..
+       #      DOUBLE PRECISION ALPHA,BETA
+       #      INTEGER INCX,INCY,LDA,M,N
+       #      CHARACTER TRANS
+       #*     .. Array Arguments ..
+       #      DOUBLE PRECISION A(LDA,*),X(*),Y(*)
+       function gemv!(trans::BlasChar,
+                      alpha::($elty), A::StridedMatrix{$elty},
+                      X::StridedVector{$elty},
+                      beta::($elty), Y::StridedVector{$elty})
+           ccall(($(string(gemv)),libblas), Void,
+                 (Ptr{Uint8}, Ptr{BlasInt}, Ptr{BlasInt},
+                  Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt},
+                  Ptr{$elty}, Ptr{BlasInt},
+                  Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}),
                  &trans, &size(A,1), &size(A,2), &alpha, A, &stride(A,2),
                  X, &stride(X,1), &beta, Y, &stride(Y,1))
            Y
        end
-       function gemv(trans::BlasChar, alpha::($elty), A::StridedMatrix{$elty}, X::StridedVector{$elty})
-           Y = Array($elty, size(A,1))
-           gemv!(trans, alpha, A, X, zero($elty), Y)
-           Y
+       function gemv(trans::BlasChar,
+                     alpha::($elty), A::StridedMatrix{$elty},
+                     X::StridedVector{$elty})
+           gemv!(trans, alpha, A, X, zero($elty),
+                 Array($elty, size(A, (trans == 'N' ? 1 : 2))))
+       end
+       function gemv(trans::BlasChar, A::StridedMatrix{$elty}, X::StridedVector{$elty})
+           gemv!(trans, one($elty), A, X, zero($elty),
+                 Array($elty, size(A, (trans == 'N' ? 1 : 2))))
        end
    end
 end
 
 # (SY) symmetric matrix-matrix and matrix-vector multiplication
-
-#     SUBROUTINE DSYMM(SIDE,UPLO,M,N,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
-#     .. Scalar Arguments ..
-#     DOUBLE PRECISION ALPHA,BETA
-#     INTEGER LDA,LDB,LDC,M,N
-#     CHARACTER SIDE,UPLO
-#     .. Array Arguments ..
-#     DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
-
-#      SUBROUTINE DSYMV(UPLO,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
-#     .. Scalar Arguments ..
-#      DOUBLE PRECISION ALPHA,BETA
-#      INTEGER INCX,INCY,LDA,N
-#      CHARACTER UPLO
-#     .. Array Arguments ..
-#      DOUBLE PRECISION A(LDA,*),X(*),Y(*)
-
-for (vfname, mfname, elty) in
-    ((:dsymv_,:dsymm_,:Float64),
-     (:ssymv_,:ssymm_,:Float32),
-     (:zsymv_,:zsymm_,:Complex128),
-     (:csymv_,:csymm_,:Complex64))
+for (mfname, vfname, elty) in
+    ((:dsymm_,:dsymv_,:Float64),
+     (:ssymm_,:ssymv_,:Float32),
+     (:zsymm_,:zsymv_,:Complex128),
+     (:csymm_,:csymv_,:Complex64))
    @eval begin
-       function symv!(uplo::BlasChar, alpha::($elty), A::StridedMatrix{$elty}, X::StridedVector{$elty},
-                      beta::($elty), Y::StridedVector{$elty})
-           m, n = size(A)
-           if m != n error("symm!: matrix A is $m by $n but must be square") end
-           if m != length(X) || m != length(Y) error("symm!: dimension mismatch") end
-           ccall(($(string(vfname)),libblas), Void,
-                 (Ptr{Uint8}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt},
-                 Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}),
-                 &uplo, &n, &alpha, A, &stride(A,2), X, &stride(X,1), &beta, Y, &stride(Y,1))
-           Y
-       end
-       function symv(uplo::BlasChar, alpha::($elty), A::StridedMatrix{$elty}, X::StridedVector{$elty})
-           symv!(uplo, alpha, A, X, zero($elty), similar(X))
-       end
-       function symm!(side::BlasChar, uplo::BlasChar, alpha::($elty), A::StridedMatrix{$elty}, B::StridedMatrix{$elty},
-                      beta::($elty), C::StridedMatrix{$elty})
-           side = uppercase(convert(Char, side))
+       #     SUBROUTINE DSYMM(SIDE,UPLO,M,N,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
+       #     .. Scalar Arguments ..
+       #     DOUBLE PRECISION ALPHA,BETA
+       #     INTEGER LDA,LDB,LDC,M,N
+       #     CHARACTER SIDE,UPLO
+       #     .. Array Arguments ..
+       #     DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
+       function symm!(side::BlasChar, uplo::BlasChar, alpha::($elty), A::StridedMatrix{$elty},
+                      B::StridedMatrix{$elty}, beta::($elty), C::StridedMatrix{$elty})
            m, n = size(C)
            k, j = size(A)
            if k != j error("symm!: matrix A is $k by $j but must be square") end
@@ -435,9 +401,37 @@ for (vfname, mfname, elty) in
                  &beta, C, &stride(C,2))
            C
        end
-       function symm(side::BlasChar, uplo::BlasChar, alpha::($elty), A::StridedMatrix{$elty}, B::StridedMatrix{$elty})
+       function symm(side::BlasChar, uplo::BlasChar, alpha::($elty), A::StridedMatrix{$elty},
+                     B::StridedMatrix{$elty})
            symm!(side, uplo, alpha, A, B, zero($elty), similar(B))
        end
+       function symm(side::BlasChar, uplo::BlasChar, A::StridedMatrix{$elty}, B::StridedMatrix{$elty})
+           symm(side, uplo, one($elty), A, B)
+       end
+       #      SUBROUTINE DSYMV(UPLO,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
+       #     .. Scalar Arguments ..
+       #      DOUBLE PRECISION ALPHA,BETA
+       #      INTEGER INCX,INCY,LDA,N
+       #      CHARACTER UPLO
+       #     .. Array Arguments ..
+       #      DOUBLE PRECISION A(LDA,*),X(*),Y(*)
+       function symv!(uplo::BlasChar, alpha::($elty), A::StridedMatrix{$elty}, x::StridedVector{$elty},
+                      beta::($elty), y::StridedVector{$elty})
+           m, n = size(A)
+           if m != n error("symm!: matrix A is $m by $n but must be square") end
+           if m != length(x) || m != length(y) error("symm!: dimension mismatch") end
+           ccall(($(string(vfname)),libblas), Void,
+                 (Ptr{Uint8}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt},
+                 Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}),
+                 &uplo, &n, &alpha, A, &stride(A,2), x, &stride(x,1), &beta, y, &stride(y,1))
+           Y
+       end
+       function symv(uplo::BlasChar, alpha::($elty), A::StridedMatrix{$elty}, x::StridedVector{$elty})
+           symv!(uplo, alpha, A, x, zero($elty), similar(x))
+       end
+       function symv(uplo::BlasChar, A::StridedMatrix{$elty}, x::StridedVector{$elty})
+           symv(uplo, one($elty), A, x)
+       end
    end
 end
 
diff --git a/test/Makefile b/test/Makefile
index 4e655a0287404..81d13ff4e00e2 100644
--- a/test/Makefile
+++ b/test/Makefile
@@ -8,7 +8,7 @@ core numbers strings unicode corelib hashing remote \
 arrayops linalg fft dct sparse bitarray suitesparse arpack \
 random math functional bigint bigfloat sorting \
 statistics glpk linprog poly file Rmath remote zlib image \
-iostring gzip integers spawn ccall parallel
+iostring gzip integers spawn ccall blas parallel
 
 $(TESTS) ::
 	$(QUIET_JULIA) $(JULIA_EXECUTABLE) ./runtests.jl $@
diff --git a/test/blas.jl b/test/blas.jl
new file mode 100644
index 0000000000000..7c58a753ca312
--- /dev/null
+++ b/test/blas.jl
@@ -0,0 +1,78 @@
+## BLAS tests - testing the interface code to BLAS routines
+for elty in (Float32, Float64, Complex64, Complex128)
+
+    o4 = ones(elty, 4)
+    z4 = zeros(elty, 4)
+
+    I4 = eye(elty, 4)
+    L4 = tril(ones(elty, (4,4)))    
+    U4 = triu(ones(elty, (4,4)))
+    Z4 = zeros(elty, (4,4))
+
+    elm1 = convert(elty, -1)    
+    el2 = convert(elty, 2)
+    v14 = convert(Vector{elty}, [1:4])
+    v41 = convert(Vector{elty}, [4:-1:1])    
+                                        # gemv
+    @assert all(BLAS.gemv('N', I4, o4) .== o4)
+    @assert all(BLAS.gemv('T', I4, o4) .== o4)
+    @assert all(BLAS.gemv('N', el2, I4, o4) .== el2 * o4)
+    @assert all(BLAS.gemv('T', el2, I4, o4) .== el2 * o4)
+    o4cp = copy(o4)
+    @assert all(BLAS.gemv!('N', one(elty), I4, o4, elm1, o4cp) .== z4)
+    @assert all(o4cp .== z4)
+    o4cp[:] = o4
+    @assert all(BLAS.gemv!('T', one(elty), I4, o4, elm1, o4cp) .== z4)
+    @assert all(o4cp .== z4)
+    @assert all(BLAS.gemv('N', U4, o4) .== v41)
+    @assert all(BLAS.gemv('N', U4, o4) .== v41)
+                                        # gemm
+    @assert all(BLAS.gemm('N', 'N', I4, I4) .== I4)
+    @assert all(BLAS.gemm('N', 'T', I4, I4) .== I4)
+    @assert all(BLAS.gemm('T', 'N', I4, I4) .== I4)
+    @assert all(BLAS.gemm('T', 'T', I4, I4) .== I4)
+    @assert all(BLAS.gemm('N', 'N', el2, I4, I4) .== el2 * I4)    
+    @assert all(BLAS.gemm('N', 'T', el2, I4, I4) .== el2 * I4)    
+    @assert all(BLAS.gemm('T', 'N', el2, I4, I4) .== el2 * I4)    
+    @assert all(BLAS.gemm('T', 'T', el2, I4, I4) .== el2 * I4)
+    I4cp = copy(I4)
+    @assert all(BLAS.gemm!('N', 'N', one(elty), I4, I4, elm1, I4cp) .== Z4)
+    @assert all(I4cp .== Z4)
+    I4cp[:] = I4
+    @assert all(BLAS.gemm!('N', 'T', one(elty), I4, I4, elm1, I4cp) .== Z4)
+    @assert all(I4cp .== Z4)
+    I4cp[:] = I4
+    @assert all(BLAS.gemm!('T', 'N', one(elty), I4, I4, elm1, I4cp) .== Z4)
+    @assert all(I4cp .== Z4)
+    I4cp[:] = I4
+    @assert all(BLAS.gemm!('T', 'T', one(elty), I4, I4, elm1, I4cp) .== Z4)
+    @assert all(I4cp .== Z4)
+    @assert all(BLAS.gemm('N', 'N', I4, U4) .== U4)
+    @assert all(BLAS.gemm('N', 'T', I4, U4) .== L4)
+                                        # gemm compared to (sy)(he)rk
+    if iscomplex(elm1)
+        @assert all(triu(BLAS.herk('U', 'N', U4)) .== triu(BLAS.gemm('N', 'T', U4, U4)))
+        @assert all(tril(BLAS.herk('L', 'N', U4)) .== tril(BLAS.gemm('N', 'T', U4, U4)))
+        @assert all(triu(BLAS.herk('U', 'N', L4)) .== triu(BLAS.gemm('N', 'T', L4, L4)))
+        @assert all(tril(BLAS.herk('L', 'N', L4)) .== tril(BLAS.gemm('N', 'T', L4, L4)))
+        @assert all(triu(BLAS.herk('U', 'T', U4)) .== triu(BLAS.gemm('T', 'N', U4, U4)))
+        @assert all(tril(BLAS.herk('L', 'T', U4)) .== tril(BLAS.gemm('T', 'N', U4, U4)))
+        @assert all(triu(BLAS.herk('U', 'T', L4)) .== triu(BLAS.gemm('T', 'N', L4, L4)))
+        @assert all(tril(BLAS.herk('L', 'T', L4)) .== tril(BLAS.gemm('T', 'N', L4, L4)))
+        ans = similar(L4)
+        @assert all(tril(BLAS.herk('L','T', L4)) .== tril(BLAS.herk!('L', 'T', one(elty), L4, zero(elty), ans)))
+        @assert all(symmetrize!(ans, 'L') .== BLAS.gemm('T', 'N', L4, L4))
+    else
+        @assert all(triu(BLAS.syrk('U', 'N', U4)) .== triu(BLAS.gemm('N', 'T', U4, U4)))
+        @assert all(tril(BLAS.syrk('L', 'N', U4)) .== tril(BLAS.gemm('N', 'T', U4, U4)))
+        @assert all(triu(BLAS.syrk('U', 'N', L4)) .== triu(BLAS.gemm('N', 'T', L4, L4)))
+        @assert all(tril(BLAS.syrk('L', 'N', L4)) .== tril(BLAS.gemm('N', 'T', L4, L4)))
+        @assert all(triu(BLAS.syrk('U', 'T', U4)) .== triu(BLAS.gemm('T', 'N', U4, U4)))
+        @assert all(tril(BLAS.syrk('L', 'T', U4)) .== tril(BLAS.gemm('T', 'N', U4, U4)))
+        @assert all(triu(BLAS.syrk('U', 'T', L4)) .== triu(BLAS.gemm('T', 'N', L4, L4)))
+        @assert all(tril(BLAS.syrk('L', 'T', L4)) .== tril(BLAS.gemm('T', 'N', L4, L4)))
+        ans = similar(L4)
+        @assert all(tril(BLAS.syrk('L','T', L4)) .== tril(BLAS.syrk!('L', 'T', one(elty), L4, zero(elty), ans)))
+        @assert all(symmetrize!(ans, 'L') .== BLAS.gemm('T', 'N', L4, L4))
+    end
+end
diff --git a/test/default.jl b/test/default.jl
index 06eb626c75ace..3843bffbc8061 100644
--- a/test/default.jl
+++ b/test/default.jl
@@ -11,6 +11,7 @@ runtests("iostring")
 # array/matrix tests
 runtests("arrayops")
 runtests("linalg")
+runtests("blas")
 runtests("fft")
 runtests("dct")
 runtests("sparse")
diff --git a/test/linalg.jl b/test/linalg.jl
index 25420fde36626..04209a3da571e 100644
--- a/test/linalg.jl
+++ b/test/linalg.jl
@@ -4,35 +4,35 @@ a     = rand(n,n)
 b     = rand(n)
 for elty in (Float32, Float64, Complex64, Complex128)
         a     = convert(Matrix{elty}, a)
-        asym  = a' + a                          # symmetric indefinite
-        apd   = a'*a                            # symmetric positive-definite
+        asym  = a' + a                  # symmetric indefinite
+        apd   = a'*a                    # symmetric positive-definite
         b     = convert(Vector{elty}, b)
         
-        capd  = chold(apd)                      # upper Cholesky factor
+        capd  = chold(apd)              # upper Cholesky factor
         r     = factors(capd)
         @assert_approx_eq r'*r apd
         @assert_approx_eq b apd * (capd\b)
         @assert_approx_eq apd * inv(capd) eye(elty, n)
-        @assert_approx_eq a*(capd\(a'*b)) b          # least squares soln for square a
+        @assert_approx_eq a*(capd\(a'*b)) b # least squares soln for square a
         @assert_approx_eq det(capd) det(apd)
 
-        l     = factors(chold(apd, 'L'))      # lower Cholesky factor
+        l     = factors(chold(apd, 'L')) # lower Cholesky factor
         @assert_approx_eq l*l' apd
 
-        cpapd = cholpd(apd)                     # pivoted Choleksy decomposition
+        cpapd = cholpd(apd)           # pivoted Choleksy decomposition
         @test rank(cpapd) == n
-        @test all(diff(diag(real(cpapd.LR))).<=0.)  # diagonal show be non-increasing
+        @test all(diff(diag(real(cpapd.LR))).<=0.) # diagonal should be non-increasing
         @assert_approx_eq b apd * (cpapd\b)
         @assert_approx_eq apd * inv(cpapd) eye(elty, n)
 
-        bc1   = BunchKaufman(asym)              # Bunch-Kaufman factor of indefinite matrix
+        bc1   = BunchKaufman(asym) # Bunch-Kaufman factor of indefinite matrix
         @assert_approx_eq inv(bc1) * asym eye(elty, n)
         @assert_approx_eq asym * (bc1\b) b
-        bc2   = BunchKaufman(apd)               # Bunch-Kaufman factors of a pos-def matrix
+        bc2   = BunchKaufman(apd) # Bunch-Kaufman factors of a pos-def matrix
         @assert_approx_eq inv(bc2) * apd eye(elty, n)
         @assert_approx_eq apd * (bc2\b) b
 
-        lua   = lud(a)                          # LU decomposition
+        lua   = lud(a)                  # LU decomposition
         l,u,p = lu(a)
         L,U,P = factors(lua)
         @test l == L && u == U && p == P
@@ -41,7 +41,7 @@ for elty in (Float32, Float64, Complex64, Complex128)
         @assert_approx_eq a * inv(lua) eye(elty, n)
         @assert_approx_eq a*(lua\b) b
 
-        qra   = qrd(a)                          # QR decomposition
+        qra   = qrd(a)                  # QR decomposition
         q,r   = factors(qra)
         @assert_approx_eq q'*q eye(elty, n)
         @assert_approx_eq q*q' eye(elty, n)
@@ -52,7 +52,7 @@ for elty in (Float32, Float64, Complex64, Complex128)
         @assert_approx_eq qra'*b Q'*b
         @assert_approx_eq a*(qra\b) b
 
-        qrpa  = qrpd(a)                         # pivoted QR decomposition
+        qrpa  = qrpd(a)                 # pivoted QR decomposition
         q,r,p = factors(qrpa)
         @assert_approx_eq q'*q eye(elty, n)
         @assert_approx_eq q*q' eye(elty, n)
@@ -62,20 +62,20 @@ for elty in (Float32, Float64, Complex64, Complex128)
         @assert_approx_eq q*r[:,invperm(p)] a
         @assert_approx_eq a*(qrpa\b) b
 
-        d,v   = eig(asym)                       # symmetric eigen-decomposition
+        d,v   = eig(asym)              # symmetric eigen-decomposition
         @assert_approx_eq asym*v[:,1] d[1]*v[:,1]
         @assert_approx_eq v*diagmm(d,v') asym
 
-        d,v   = eig(a)                          # non-symmetric eigen decomposition
+        d,v   = eig(a)                 # non-symmetric eigen decomposition
         for i in 1:size(a,2) @assert_approx_eq a*v[:,i] d[i]*v[:,i] end
     
-        u, q, v = schur(a)                         # Schur
+        u, q, v = schur(a)             # Schur
         @assert_approx_eq q*u*q' a
         @assert_approx_eq sort(real(v)) sort(real(d))
         @assert_approx_eq sort(imag(v)) sort(imag(d))
         @test istriu(u) || isreal(a)
 
-        u,s,vt = svdt(a)                        # singular value decomposition
+        u,s,vt = svdt(a)               # singular value decomposition
         @assert_approx_eq u*diagmm(s,vt) a
     
         x = a \ b
@@ -87,63 +87,60 @@ for elty in (Float32, Float64, Complex64, Complex128)
         x = tril(a)\b
         @assert_approx_eq tril(a)*x b
 
-        # Test null
+                                        # Test null
         a15null = null(a[:,1:5]')
         @assert_approx_eq_eps norm(a[:,1:5]'a15null) zero(elty) n*eps(one(elty))
         @assert_approx_eq_eps norm(a15null'a[:,1:5]) zero(elty) n*eps(one(elty))
         @test size(null(b), 2) == 0
 
-        # Test pinv
+                                        # Test pinv
         pinva15 = pinv(a[:,1:5])
         @assert_approx_eq a[:,1:5]*pinva15*a[:,1:5] a[:,1:5]
         @assert_approx_eq pinva15*a[:,1:5]*pinva15 pinva15
     
-        # Complex complex rhs real lhs
+                                        # Complex vector rhs
         x = a\complex(b)
         @assert_approx_eq a*x complex(b)
 
-        # Test cond
+                                        # Test cond
         @assert_approx_eq_eps cond(a, 1) 4.837320054554436e+02 0.01
         @assert_approx_eq_eps cond(a, 2) 1.960057871514615e+02 0.01
         @assert_approx_eq_eps cond(a, Inf) 3.757017682707787e+02 0.01
         @assert_approx_eq_eps cond(a[:,1:5]) 10.233059337453463 0.01
 
-        # Matrix square root
+                                        # Matrix square root
         asq = sqrtm(a)
         @assert_approx_eq asq*asq a
 end
+
+## Least squares solutions
 a = [ones(20) 1:20 1:20]
 b = reshape(eye(8, 5), 20, 2)
 for elty in (Float32, Float64, Complex64, Complex128)
         a = convert(Matrix{elty}, a)
         b = convert(Matrix{elty}, b)
 
-        # Matrix and vector multiplication
-        @assert_approx_eq b'b convert(Matrix{elty}, [3 0; 0 2])
-        @assert_approx_eq b'b[:,1] convert(Vector{elty}, [3, 0])
-        @assert_approx_eq dot(b[:,1], b[:,1]) convert(elty, 3.0)
-
-        # Least squares        
-        x = a[:,1:2]\b[:,1]                             # Vector rhs
+        x = a[:,1:2]\b[:,1]             # Vector rhs
         @assert_approx_eq ((a[:,1:2]*x-b[:,1])'*(a[:,1:2]*x-b[:,1]))[1] convert(elty, 2.546616541353384)
     
-        x = a[:,1:2]\b                                  # Matrix rhs
+        x = a[:,1:2]\b                  # Matrix rhs
         @assert_approx_eq det((a[:,1:2]*x-b)'*(a[:,1:2]*x-b)) convert(elty, 4.437969924812031)
     
-        x = a\b                            # Rank deficient
+        x = a\b                         # Rank deficient
         @assert_approx_eq det((a*x-b)'*(a*x-b)) convert(elty, 4.437969924812031)
 
-        x = convert(Matrix{elty}, [1 0 0; 0 1 -1]) \ convert(Vector{elty}, [1,1]) # Underdetermined minimum norm
+                                        # Underdetermined minimum norm
+        x = convert(Matrix{elty}, [1 0 0; 0 1 -1]) \ convert(Vector{elty}, [1,1]) 
         @assert_approx_eq x convert(Vector{elty}, [1, 0.5, -0.5])
 
-        # symmetric, positive definite
+                                        # symmetric, positive definite
         @assert_approx_eq inv(convert(Matrix{elty}, [6. 2; 2 1])) convert(Matrix{elty}, [0.5 -1; -1 3])
-        # symmetric, negative definite
+                                        # symmetric, negative definite
         @assert_approx_eq inv(convert(Matrix{elty}, [1. 2; 2 1])) convert(Matrix{elty}, [-1. 2; 2 -1]/3)               
 end
 
 ## Test Julia fallbacks to BLAS routines
-# matrices with zero dimensions
+                                        # matrices with zero dimensions
 @test ones(0,5)*ones(5,3) == zeros(0,3)
 @test ones(3,5)*ones(5,0) == zeros(3,0)
 @test ones(3,0)*ones(0,4) == zeros(3,4)
@@ -151,7 +148,7 @@ end
 @test ones(0,0)*ones(0,4) == zeros(0,4)
 @test ones(3,0)*ones(0,0) == zeros(3,0)
 @test ones(0,0)*ones(0,0) == zeros(0,0)
-# 2x2
+                                        # 2x2
 A = [1 2; 3 4]
 B = [5 6; 7 8]
 @test A*B == [19 22; 43 50]
@@ -164,7 +161,7 @@ Bi = B+(2.5*im).*A[[2,1],[2,1]]
 @test Ac_mul_B(Ai, Bi) == [68.5-12im 57.5-28im; 88-3im 76.5-25im]
 @test A_mul_Bc(Ai, Bi) == [64.5+5.5im 43+31.5im; 104-18.5im 80.5+31.5im]
 @test Ac_mul_Bc(Ai, Bi) == [-28.25-66im 9.75-58im; -26-89im 21-73im]
-# 3x3
+                                        # 3x3
 A = [1 2 3; 4 5 6; 7 8 9]-5
 B = [1 0 5; 6 -10 3; 2 -4 -1]
 @test A*B == [-26 38 -27; 1 -4 -6; 28 -46 15]
@@ -177,7 +174,7 @@ Bi = B+(2.5*im).*A[[2,1,3],[2,3,1]]
 @test Ac_mul_B(Ai, Bi) == [-21+2im -1.75+49im -51.25+19.5im; 25.5+56.5im -7-35.5im 22+35.5im; -3+12im -32.25+43im -34.75-2.5im]
 @test A_mul_Bc(Ai, Bi) == [-20.25+15.5im -28.75-54.5im 22.25+68.5im; -12.25+13im -15.5+75im -23+27im; 18.25+im 1.5+94.5im -27-54.5im]
 @test Ac_mul_Bc(Ai, Bi) == [1+2im 20.75+9im -44.75+42im; 19.5+17.5im -54-36.5im 51-14.5im; 13+7.5im 11.25+31.5im -43.25-14.5im]
-# Generic integer matrix multiplication
+                                        # Generic integer matrix multiplication
 A = [1 2 3; 4 5 6] - 3
 B = [2 -2; 3 -5; -4 7]
 @test A*B == [-7 9; -4 9]
@@ -189,13 +186,13 @@ A = rand(1:20, 5, 5) - 10
 B = rand(1:20, 5, 5) - 10
 @test At_mul_B(A, B) == A'*B
 @test A_mul_Bt(A, B) == A*B'
-# Preallocated
+                                        # Preallocated
 C = Array(Int, size(A, 1), size(B, 2))
 @test A_mul_B(C, A, B) == A*B
 @test At_mul_B(C, A, B) == A'*B
 @test A_mul_Bt(C, A, B) == A*B'
 @test At_mul_Bt(C, A, B) == A'*B'
-# matrix algebra with subarrays of floats (stride != 1)
+                                        # matrix algebra with subarrays of floats (stride != 1)
 A = reshape(float64(1:20),5,4)
 Aref = A[1:2:end,1:2:end]
 Asub = sub(A, 1:2:5, 1:2:4)
@@ -208,7 +205,7 @@ Aref = Ai[1:2:end,1:2:end]
 Asub = sub(Ai, 1:2:5, 1:2:4)
 @test Ac_mul_B(Asub, Asub) == Ac_mul_B(Aref, Aref)
 @test A_mul_Bc(Asub, Asub) == A_mul_Bc(Aref, Aref)
-# syrk & herk
+                                        # syrk & herk
 A = reshape(1:1503, 501, 3)-750.0
 res = float64([135228751 9979252 -115270247; 9979252 10481254 10983256; -115270247 10983256 137236759])
 @test At_mul_B(A, A) == res
@@ -223,7 +220,7 @@ Asub = sub(Ai, 1:2:2*cutoff, 1:3)
 Aref = Ai[1:2:2*cutoff, 1:3]
 @test Ac_mul_B(Asub, Asub) == Ac_mul_B(Aref, Aref)
 
-# Matrix exponential and Hessenberg
+                                        # Matrix exponential
 for elty in (Float32, Float64, Complex64, Complex128)
         A1  = convert(Matrix{elty}, [4 2 0; 1 4 1; 1 1 4])
         eA1 = convert(Matrix{elty}, [147.866622446369 127.781085523181  127.781085523182;
@@ -247,27 +244,27 @@ for elty in (Float32, Float64, Complex64, Complex128)
         0.135335281175235 0.406005843524598 0.541341126763207]')
         @assert_approx_eq expm(A3) eA3
 
-        # Hessenberg
+                                        # Hessenberg
         @assert_approx_eq hess(A1) convert(Matrix{elty}, 
                         [4.000000000000000  -1.414213562373094  -1.414213562373095
                         -1.414213562373095   4.999999999999996  -0.000000000000000
                                          0  -0.000000000000002   3.000000000000000])
 end
 
-# matmul for types w/o sizeof (issue #1282)
+                                        # matmul for types w/o sizeof (issue #1282)
 A = Array(ComplexPair{Int},10,10)
 A[:] = complex(1,1)
 A2 = A^2
 @test A2[1,1] == 20im
 
-# basic tridiagonal operations
+                                        # basic tridiagonal operations
 n = 5
 d = 1 + rand(n)
 dl = -rand(n-1)
 du = -rand(n-1)
 v = randn(n)
 B = randn(n,2)
-# Woodbury
+                                        # Woodbury
 U = randn(n,2)
 V = randn(2,n)
 C = randn(2,2)
@@ -285,8 +282,7 @@ for elty in (Float32, Float64, Complex64, Complex128)
             F[i+1,i] = dl[i]
         end
         @test full(T) == F
-
-        # tridiagonal linear algebra
+                                        # tridiagonal linear algebra
         v = convert(Vector{elty}, v)
         @assert_approx_eq T*v F*v
         invFv = F\v
@@ -298,24 +294,21 @@ for elty in (Float32, Float64, Complex64, Complex128)
         x = Tlu\v
         @assert_approx_eq x invFv
         @assert_approx_eq det(T) det(F)
-
-        # symmetric tridiagonal
+                                        # symmetric tridiagonal
         Ts = SymTridiagonal(d, dl)
         Fs = full(Ts)
         invFsv = Fs\v
         Tldlt = ldltd(Ts)
         x = Tldlt\v
         @assert_approx_eq x invFsv
-
-        # eigenvalues/eigenvectors of symmetric tridiagonal
+                                        # eigenvalues/eigenvectors of symmetric tridiagonal
         if elty === Float32 || elty === Float64
             DT, VT = eig(Ts)
             D, Vecs = eig(Fs)
             @assert_approx_eq DT D
             @assert_approx_eq abs(VT'Vecs) eye(elty, n)
         end
-        
-        # Woodbury
+                                        # Woodbury
         U = convert(Matrix{elty}, U)
         V = convert(Matrix{elty}, V)
         C = convert(Matrix{elty}, C)
@@ -355,10 +348,10 @@ for elty in (Float32, Float64, Complex64, Complex128)
         @assert_approx_eq_eps det(ones(elty, 3,3)) zero(elty) 3*eps(one(elty))
 end
 
-# LAPACK tests
+                                        # LAPACK tests
 Ainit = randn(5,5)
 for elty in (Float32, Float64, Complex64, Complex128)
-        # syevr!
+                                        # syevr!
         A = convert(Array{elty, 2}, Ainit)
         Asym = A'A
         Z = Array(elty, 5, 5)

From bc7dec79035656a92203e8e7651fec3b755dd092 Mon Sep 17 00:00:00 2001
From: "Viral B. Shah" <viral@mayin.org>
Date: Wed, 6 Feb 2013 18:24:13 -0500
Subject: [PATCH 02/11] Get rid of chold. Add documentation for chol.

---
 base/lapack.jl       |  4 ++++
 base/linalg_dense.jl | 33 +++++++++++++++++----------------
 doc/stdlib/base.rst  |  8 ++++++--
 3 files changed, 27 insertions(+), 18 deletions(-)

diff --git a/base/lapack.jl b/base/lapack.jl
index 9824911db7452..e7d8828cd2022 100644
--- a/base/lapack.jl
+++ b/base/lapack.jl
@@ -12,6 +12,10 @@ type LapackException <: Exception
     info::BlasInt
 end
 
+type PosDefException <: Exception
+    info::BlasInt
+end
+
 type SingularException <: Exception
     info::BlasInt
 end
diff --git a/base/linalg_dense.jl b/base/linalg_dense.jl
index bc6c56658c22d..562549fdae0d6 100644
--- a/base/linalg_dense.jl
+++ b/base/linalg_dense.jl
@@ -389,13 +389,15 @@ end
 
 \{T<:BlasFloat}(B::BunchKaufman{T}, R::StridedVecOrMat{T}) =
     LAPACK.sytrs!(B.UL, B.LD, B.ipiv, copy(R))
-    
+
+# Dense Cholesky factorization
+
 type CholeskyDense{T<:BlasFloat} <: Factorization{T}
     LR::Matrix{T}
     UL::BlasChar
-    function CholeskyDense(A::Matrix{T}, UL::BlasChar)
+    function CholeskyDense{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar)
         A, info = LAPACK.potrf!(UL, A)
-        if info != 0 error("Matrix A not positive-definite") end
+        if info != 0; throw(LAPACK.PosDefException(info)); end
         if UL == 'U'
             new(triu!(A), UL)
         elseif UL == 'L'
@@ -407,7 +409,7 @@ type CholeskyDense{T<:BlasFloat} <: Factorization{T}
 end
 
 size(C::CholeskyDense) = size(C.LR)
-size(C::CholeskyDense,d::Integer) = size(C.LR,d)
+size(C::CholeskyDense, d::Integer) = size(C.LR, d)
 
 factors(C::CholeskyDense) = C.LR
 
@@ -417,28 +419,27 @@ factors(C::CholeskyDense) = C.LR
 function det{T}(C::CholeskyDense{T})
     ff = C.LR
     dd = one(T)
-    for i in 1:size(ff,1) dd *= abs2(ff[i,i]) end
+    for i in 1:size(ff,1); dd *= abs2(ff[i,i]); end
     dd
 end
-    
+
 function inv(C::CholeskyDense)
     Ci, info = LAPACK.potri!(C.UL, copy(C.LR))
-    if info != 0 error("Matrix singular") end 
+    if info != 0; throw(LAPACK.SingularException(info)); end 
     symmetrize!(Ci, C.UL)
 end
 
 ## Should these functions check that the matrix is Hermitian?
-chold!{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar) = CholeskyDense{T}(A, UL)
-chold!{T<:BlasFloat}(A::Matrix{T}) = chold!(A, 'U')
-chold{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar) = chold!(copy(A), UL)
-chold{T<:Number}(A::Matrix{T}, UL::BlasChar) = chold(float64(A), UL)
-chold{T<:Number}(A::Matrix{T}) = chold(A, 'U')
-
-## Matlab (and R) compatible
-chol(A::Matrix, UL::BlasChar) = factors(chold(A, UL))
-chol(A::Matrix) = chol(A, 'U')
+chol!{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar) = CholeskyDense{T}(A, UL)
+chol!{T<:BlasFloat}(A::Matrix{T}) = chol!(A, 'U')
+chol{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar) = chol!(copy(A), UL)
+chol{T<:Number}(A::Matrix{T}, UL::BlasChar) = chol(float64(A), UL)
+chol{T<:Number}(A::Matrix{T}) = chol(A, 'U')
+
 chol(x::Number) = imag(x) == 0 && real(x) > 0 ? sqrt(x) : error("Argument not positive-definite")
 
+# Pivoted Cholesky factorization
+
 type CholeskyDensePivoted{T<:BlasFloat} <: Factorization{T}
     LR::Matrix{T}
     UL::BlasChar
diff --git a/doc/stdlib/base.rst b/doc/stdlib/base.rst
index 87dee32e81e1c..993f5d5aa0a4c 100644
--- a/doc/stdlib/base.rst
+++ b/doc/stdlib/base.rst
@@ -1720,9 +1720,13 @@ Linear algebra functions in Julia are largely implemented by calling functions f
 
    Compute LU factorization. LU is an "LU factorization" type that can be used as an ordinary matrix.
 
-.. function:: chol(A)
+.. function:: chol(A, [LU])
 
-   Compute Cholesky factorization
+   Compute the Cholesky factorization of ``A`` and return a ``CholeskyDense`` object. ``LU`` may be 'L' for using the lower part or 'U' for the upper part. The default is to use 'U'. The ``factors(chol(A))`` returns the triangular matrix containing the factorization. The following functions are available for ``CholeskyDense`` objects: ``size``, ``factors``, ``\``, ``inv``, ``det``.
+
+.. function:: chol!(A, [LU])
+
+   ``chol!`` is the same as ``chol``, but overwrites the input matrix A with the factorization.
 
 .. function:: qr(A)
 

From fe13a03ed63539281644093d119515ec932d4f30 Mon Sep 17 00:00:00 2001
From: "Viral B. Shah" <viral@mayin.org>
Date: Wed, 6 Feb 2013 18:36:10 -0500
Subject: [PATCH 03/11] Rename cholpd to cholpivot. Add docs for cholpivot.

---
 base/exports.jl      |  7 +++----
 base/lapack.jl       |  4 ++++
 base/linalg_dense.jl | 32 ++++++++++++++++++--------------
 doc/stdlib/base.rst  |  8 ++++++++
 test/linalg.jl       |  6 +++---
 5 files changed, 36 insertions(+), 21 deletions(-)

diff --git a/base/exports.jl b/base/exports.jl
index 4e5efa1380e5c..a38d9fa1fea28 100644
--- a/base/exports.jl
+++ b/base/exports.jl
@@ -523,10 +523,9 @@ export
 
 # linear algebra
     chol,
-    chold,
-    chold!,
-    cholpd,
-    cholpd!,
+    chol!,
+    cholpivot,
+    cholpivot!,
     cond,
     cross,
     ctranspose,
diff --git a/base/lapack.jl b/base/lapack.jl
index e7d8828cd2022..5e6627d4561fa 100644
--- a/base/lapack.jl
+++ b/base/lapack.jl
@@ -16,6 +16,10 @@ type PosDefException <: Exception
     info::BlasInt
 end
 
+type RankDeficientException <: Exception
+    info::BlasInt
+end
+
 type SingularException <: Exception
     info::BlasInt
 end
diff --git a/base/linalg_dense.jl b/base/linalg_dense.jl
index 562549fdae0d6..068cc0ef3e3fc 100644
--- a/base/linalg_dense.jl
+++ b/base/linalg_dense.jl
@@ -448,6 +448,7 @@ type CholeskyDensePivoted{T<:BlasFloat} <: Factorization{T}
     tol::Real
     function CholeskyDensePivoted(A::Matrix{T}, UL::BlasChar, tol::Real)
         A, piv, rank, info = LAPACK.pstrf!(UL, A, tol)
+        if info != 0; throw(LAPACK.RankDeficientException(info)); end
         if UL == 'U'
             new(triu!(A), UL, piv, rank, tol)
         elseif UL == 'L'
@@ -464,11 +465,12 @@ size(C::CholeskyDensePivoted,d::Integer) = size(C.LR,d)
 factors(C::CholeskyDensePivoted) = C.LR, C.piv
 
 function \{T<:BlasFloat}(C::CholeskyDensePivoted{T}, B::StridedVector{T})
-    if C.rank < size(C.LR, 1) error("Matrix is not positive-definite") end
+    if C.rank < size(C.LR, 1); throw(LAPACK.RankDeficientException(info)); end
     LAPACK.potrs!(C.UL, C.LR, copy(B)[C.piv])[invperm(C.piv)]
 end
+
 function \{T<:BlasFloat}(C::CholeskyDensePivoted{T}, B::StridedMatrix{T})
-    if C.rank < size(C.LR, 1) error("Matrix is not positive-definite") end
+    if C.rank < size(C.LR, 1); throw(LAPACK.RankDeficientException(info)); end
     LAPACK.potrs!(C.UL, C.LR, copy(B)[C.piv,:])[invperm(C.piv),:]
 end
 
@@ -491,18 +493,20 @@ function inv(C::CholeskyDensePivoted)
 end
 
 ## Should these functions check that the matrix is Hermitian?
-cholpd!{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar, tol::Real) = CholeskyDensePivoted{T}(A, UL, tol)
-cholpd!{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar) = cholpd!(A, UL, -1.)
-cholpd!{T<:BlasFloat}(A::Matrix{T}, tol::Real) = cholpd!(A, 'U', tol)
-cholpd!{T<:BlasFloat}(A::Matrix{T}) = cholpd!(A, 'U', -1.)
-cholpd{T<:Number}(A::Matrix{T}, UL::BlasChar, tol::Real) = cholpd(float64(A), UL, tol)
-cholpd{T<:Number}(A::Matrix{T}, UL::BlasChar) = cholpd(float64(A), UL, -1.)
-cholpd{T<:Number}(A::Matrix{T}, tol::Real) = cholpd(float64(A), true, tol)
-cholpd{T<:Number}(A::Matrix{T}) = cholpd(float64(A), 'U', -1.)
-cholpd{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar, tol::Real) = cholpd!(copy(A), UL, tol)
-cholpd{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar) = cholpd!(copy(A), UL, -1.)
-cholpd{T<:BlasFloat}(A::Matrix{T}, tol::Real) = cholpd!(copy(A), 'U', tol)
-cholpd{T<:BlasFloat}(A::Matrix{T}) = cholpd!(copy(A), 'U', -1.)
+cholpivot!{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar, tol::Real) = CholeskyDensePivoted{T}(A, UL, tol)
+cholpivot!{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar) = cholpivot!(A, UL, -1.)
+cholpivot!{T<:BlasFloat}(A::Matrix{T}, tol::Real) = cholpivot!(A, 'U', tol)
+cholpivot!{T<:BlasFloat}(A::Matrix{T}) = cholpivot!(A, 'U', -1.)
+cholpivot{T<:Number}(A::Matrix{T}, UL::BlasChar, tol::Real) = cholpivot(float64(A), UL, tol)
+cholpivot{T<:Number}(A::Matrix{T}, UL::BlasChar) = cholpivot(float64(A), UL, -1.)
+cholpivot{T<:Number}(A::Matrix{T}, tol::Real) = cholpivot(float64(A), true, tol)
+cholpivot{T<:Number}(A::Matrix{T}) = cholpivot(float64(A), 'U', -1.)
+cholpivot{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar, tol::Real) = cholpivot!(copy(A), UL, tol)
+cholpivot{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar) = cholpivot!(copy(A), UL, -1.)
+cholpivot{T<:BlasFloat}(A::Matrix{T}, tol::Real) = cholpivot!(copy(A), 'U', tol)
+cholpivot{T<:BlasFloat}(A::Matrix{T}) = cholpivot!(copy(A), 'U', -1.)
+
+# LU Factorization
 
 type LUDense{T} <: Factorization{T}
     lu::Matrix{T}
diff --git a/doc/stdlib/base.rst b/doc/stdlib/base.rst
index 993f5d5aa0a4c..d4b4ce8d6a967 100644
--- a/doc/stdlib/base.rst
+++ b/doc/stdlib/base.rst
@@ -1728,6 +1728,14 @@ Linear algebra functions in Julia are largely implemented by calling functions f
 
    ``chol!`` is the same as ``chol``, but overwrites the input matrix A with the factorization.
 
+.. function:: cholpivot(A, [LU])
+
+   Compute the pivoted Cholesky factorization of ``A`` and return a ``CholeskyDensePivoted`` object. ``LU`` may be 'L' for using the lower part or 'U' for the upper part. The default is to use 'U'. The ``factors(cholpivot(A))`` returns the triangular matrix containing the factorization. The following functions are available for ``CholeskyDensePivoted`` objects: ``size``, ``factors``, ``\``, ``inv``, ``det``.
+
+.. function:: cholpivot!(A, [LU])
+
+   ``cholpivot!`` is the same as ``cholpivot``, but overwrites the input matrix A with the factorization.
+
 .. function:: qr(A)
 
    Compute QR factorization
diff --git a/test/linalg.jl b/test/linalg.jl
index 04209a3da571e..5918c1d1d38e5 100644
--- a/test/linalg.jl
+++ b/test/linalg.jl
@@ -8,7 +8,7 @@ for elty in (Float32, Float64, Complex64, Complex128)
         apd   = a'*a                    # symmetric positive-definite
         b     = convert(Vector{elty}, b)
         
-        capd  = chold(apd)              # upper Cholesky factor
+        capd  = chol(apd)              # upper Cholesky factor
         r     = factors(capd)
         @assert_approx_eq r'*r apd
         @assert_approx_eq b apd * (capd\b)
@@ -16,10 +16,10 @@ for elty in (Float32, Float64, Complex64, Complex128)
         @assert_approx_eq a*(capd\(a'*b)) b # least squares soln for square a
         @assert_approx_eq det(capd) det(apd)
 
-        l     = factors(chold(apd, 'L')) # lower Cholesky factor
+        l     = factors(chol(apd, 'L')) # lower Cholesky factor
         @assert_approx_eq l*l' apd
 
-        cpapd = cholpd(apd)           # pivoted Choleksy decomposition
+        cpapd = cholpivot(apd)           # pivoted Choleksy decomposition
         @test rank(cpapd) == n
         @test all(diff(diag(real(cpapd.LR))).<=0.) # diagonal should be non-increasing
         @assert_approx_eq b apd * (cpapd\b)

From ccf2024a00d2325644dfb9c179a1d94b1b2d2533 Mon Sep 17 00:00:00 2001
From: Douglas Bates <dmbates@gmail.com>
Date: Wed, 6 Feb 2013 17:37:00 -0600
Subject: [PATCH 04/11] Improve UMFpack interface

- use lud and lud! methods for the decomposition
- allocate umf_conf and umf_info and use them throughout
- preserve the pointers to the symbolic and numeric factorizations in
  the UmfpackLU object.  Also keep around the decremented
  index/pointer vectors.
---
 extras/suitesparse.jl | 732 ++++++++++++++++++++----------------------
 test/suitesparse.jl   |   1 +
 2 files changed, 355 insertions(+), 378 deletions(-)

diff --git a/extras/suitesparse.jl b/extras/suitesparse.jl
index 7c89c6f0b2f36..81dd1b2cb62e6 100644
--- a/extras/suitesparse.jl
+++ b/extras/suitesparse.jl
@@ -1,7 +1,8 @@
 module SuiteSparse
 
 import Base.SparseMatrixCSC, Base.size, Base.nnz, Base.eltype, Base.show
-import Base.triu, Base.norm, Base.solve, Base.(\), Base.ctranspose, Base.transpose
+import Base.triu, Base.norm, Base.solve, Base.(\), Base.lu
+import Base.Ac_ldiv_B, Base.At_ldiv_B, Base.lud, Base.lud!
 
 import Base.BlasInt
 import Base.blas_int
@@ -13,93 +14,341 @@ export                                  # types
     CholmodFactor,
     CholmodDense,
     CholmodSparseOut,
-    UmfpackPtr,
-    UmfpackLU,
+    UmfpackLU,                          # call these lud and lud! instead?
     UmfpackLU!,
-    UmfpackLUTrans,
                                         # methods
-    chm_aat, # drop prefix?
-    eltype,  #? maybe not
-    indtype, #? maybe not
-    nnz,
-    show,
-    size,
-    solve,
-    \,
-    At_ldiv_B,
-    Ac_ldiv_B
+    decrement,
+    decrement!,
+    increment,
+    increment!,
+    indtype,
+    show_umf_ctrl,
+    show_umf_info
 
 include("suitesparse_h.jl")
 
-const libsuitesparse_wrapper = "libsuitesparse_wrapper"
-const libcholmod = "libcholmod"
-const libumfpack = "libumfpack"
-const libspqr = "libspqr"
-
-const _chm_aat       = (:cholmod_aat, libcholmod)
-const _chm_amd       = (:cholmod_amd, libcholmod)
-const _chm_analyze   = (:cholmod_analyze, libcholmod)
-const _chm_colamd    = (:cholmod_colamd, libcholmod)
-const _chm_copy      = (:cholmod_copy, libcholmod)
-const _chm_factorize = (:cholmod_factorize, libcholmod)
-const _chm_free_dn   = (:cholmod_free_dense, libcholmod)
-const _chm_free_fa   = (:cholmod_free_factor, libcholmod)
-const _chm_free_sp   = (:cholmod_free_sparse, libcholmod)
-const _chm_print_dn  = (:cholmod_print_dense, libcholmod)
-const _chm_print_fa  = (:cholmod_print_factor, libcholmod)
-const _chm_print_sp  = (:cholmod_print_sparse, libcholmod)
-const _chm_solve     = (:cholmod_solve, libcholmod)
-const _chm_sort      = (:cholmod_sort, libcholmod)
-const _chm_submatrix = (:cholmod_submatrix, libcholmod)
-
-const _spqr_C_QR                = (:SuiteSparseQR_C_QR, libspqr)
-const _spqr_C_backslash         = (:SuiteSparseQR_C_backslash, libspqr)
-const _spqr_C_backslash_default = (:SuiteSparseQR_C_backslash_default, libspqr)
-const _spqr_C_backslash_sparse  = (:SuiteSparseQR_C_backslash_sparse, libspqr)
-const _spqr_C_factorize         = (:SuiteSparseQR_C_factorize, libspqr)
-const _spqr_C_symbolic          = (:SuiteSparseQR_C_symbolic, libspqr)
-const _spqr_C_numeric           = (:SuiteSparseQR_C_numeric, libspqr)
-const _spqr_C_free              = (:SuiteSparseQR_C_free, libspqr)
-const _spqr_C_solve             = (:SuiteSparseQR_C_solve, libspqr)
-const _spqr_C_qmult             = (:SuiteSparseQR_C_qmult, libspqr)
-
 type MatrixIllConditionedException <: Exception end
 
-function convert_to_0_based_indexing!(S::SparseMatrixCSC)
-    for i=1:(S.colptr[end]-1); S.rowval[i] -= 1; end
-    for i=1:length(S.colptr); S.colptr[i] -= 1; end
-    return S
+function decrement!{T<:Integer}(A::AbstractArray{T})
+    for i in 1:length(A) A[i] -= one(T) end
+    A
+end
+function decrement{T<:Integer}(A::AbstractArray{T})
+    B = similar(A)
+    for i in 1:length(B) B[i] = A[i] - one(T) end
+    B
+end
+function increment!{T<:Integer}(A::AbstractArray{T})
+    for i in 1:length(A) A[i] += one(T) end
+    A
+end
+function increment{T<:Integer}(A::AbstractArray{T})
+    increment!(copy(A))
 end
 
-function convert_to_1_based_indexing!(S::SparseMatrixCSC)
-    for i=1:length(S.colptr); S.colptr[i] += 1; end
-    for i=1:(S.colptr[end]-1); S.rowval[i] += 1; end
-    return S
+typealias CHMITypes Union(Int32,Int64)       # also ITypes for UMFPACK
+typealias CHMVTypes Union(Complex64, Complex128, Float32, Float64)
+typealias UMFVTypes Union(Float64,Complex128)
+
+## UMFPACK
+
+# the control and info arrays
+const umf_ctrl = Array(Float64, UMFPACK_CONTROL)
+ccall((:umfpack_dl_defaults, :libumfpack), Void, (Ptr{Float64},), umf_ctrl)
+const umf_info = Array(Float64, UMFPACK_INFO)
+
+function show_umf_ctrl(level::Real)
+    old_prt::Float64 = umf_ctrl[1]
+    umf_ctrl[1] = float64(level)
+    ccall((:umfpack_dl_report_control, :libumfpack), Void, (Ptr{Float64},), umf_ctrl)
+    umf_ctrl[1] = old_prt
 end
+show_umf_ctrl() = show_umf_ctrl(2.)
 
-convert_to_0_based_indexing(S) = convert_to_0_based_indexing!(copy(S))
-convert_to_1_based_indexing(S) = convert_to_1_based_indexing!(copy(S))
+function show_umf_info(level::Real)
+    old_prt::Float64 = umf_ctrl[1]
+    umf_ctrl[1] = float64(level)
+    ccall((:umfpack_dl_report_info, :libumfpack), Void,
+          (Ptr{Float64}, Ptr{Float64}), umf_ctrl, umf_info)
+    umf_ctrl[1] = old_prt
+end
+show_umf_info() = show_umf_info(2.)
 
-## CHOLMOD
+# Wrapper for memory allocated by umfpack. Carry along the value and index types.
+## type UmfpackPtr{Tv<:UMFVTypes,Ti<:CHMITypes}
+##     val::Vector{Ptr{Void}} 
+## end
 
-typealias CHMVTypes Union(Complex64, Complex128, Float32, Float64)
-typealias CHMITypes Union(Int32, Int64)
+type UmfpackLU{Tv<:UMFVTypes,Ti<:CHMITypes} <: Factorization{Tv}
+    symbolic::Ptr{Void}
+    numeric::Ptr{Void}
+    m::Int
+    n::Int
+    colptr::Vector{Ti}                  # 0-based column pointers
+    rowval::Vector{Ti}                  # 0-based row indices
+    nzval::Vector{Tv}
+end
+
+function lud{Tv<:UMFVTypes,Ti<:CHMITypes}(S::SparseMatrixCSC{Tv,Ti})
+    zerobased = S.colptr[1] == 0
+    lu = UmfpackLU(C_NULL, C_NULL, S.m, S.n,
+                   zerobased ? copy(S.colptr) : decrement(S.colptr),
+                   zerobased ? copy(S.rowval) : decrement(S.rowval),
+                   copy(S.nzval))
+    umfpack_numeric!(lu)
+end
+
+function lud!{Tv<:UMFVTypes,Ti<:CHMITypes}(S::SparseMatrixCSC{Tv,Ti})
+    zerobased = S.colptr[1] == 0
+    UmfpackLU(C_NULL, C_NULL, S.m, S.n,
+              zerobased ? S.colptr : decrement!(S.colptr),
+              zerobased ? S.rowval : decrement!(S.rowval),
+              S.nzval)
+end
 
-function chm_itype{Tv,Ti}(S::SparseMatrixCSC{Tv,Ti})
-    if !(Ti<:CHMITypes) error("chm_itype: indtype(S) must be in CHMITypes") end
-    Ti == Int32 ? CHOLMOD_INT : CHOLMOD_LONG
+function show(io::IO, f::UmfpackLU)
+    @printf(io, "UMFPACK LU Factorization of a %d-by-%d sparse matrix\n",
+           f.m, f.n)
+    if f.numeric != C_NULL println(f.numeric) end
 end
 
-function chm_xtype{T}(S::SparseMatrixCSC{T})
-    if !(T<:CHMVTypes) error("chm_xtype: eltype(S) must be in CHMVTypes") end
-    T <: Complex ? CHOLMOD_COMPLEX : CHOLMOD_REAL
+# Solve with Factorization
+
+(\){T<:UMFVTypes}(fact::UmfpackLU{T}, b::Vector{T}) = umfpack_solve(fact, b)
+(\){Ts<:UMFVTypes,Tb<:Number}(fact::UmfpackLU{Ts}, b::Vector{Tb}) = fact\convert(Vector{Ts},b)
+
+# Solve directly with matrix
+
+(\)(S::SparseMatrixCSC, b::Vector) = lud(S) \ b
+At_ldiv_B{T<:UMFVTypes}(S::SparseMatrixCSC{T}, b::Vector{T}) = umfpack_solve(lud(S), b, UMFPACK_Aat)
+function At_ldiv_B{Ts<:UMFVTypes,Tb<:Number}(S::SparseMatrixCSC{Ts}, b::Vector{Tb})
+    ## should be more careful here in case Ts<:Real and Tb<:Complex
+    At_ldiv_B(S, convert(Vector{Ts}, b))
+end
+Ac_ldiv_B{T<:UMFVTypes}(S::SparseMatrixCSC{T}, b::Vector{T}) = umfpack_solve(lud(S), b, UMFPACK_At)
+function Ac_ldiv_B{Ts<:UMFVTypes,Tb<:Number}(S::SparseMatrixCSC{Ts}, b::Vector{Tb})
+    ## should be more careful here in case Ts<:Real and Tb<:Complex
+    Ac_ldiv_B(S, convert(Vector{Ts}, b))
 end
 
-function chm_dtype{T}(S::SparseMatrixCSC{T})
-    if !(T<:CHMVTypes) error("chm_dtype: eltype(S) must be in CHMVTypes") end
-    T <: Union(Float32, Complex64) ? CHOLMOD_SINGLE : CHOLMOD_DOUBLE
+## Wrappers around UMFPACK routines
+
+for (f_sym_r, f_num_r, f_sym_c, f_num_c, itype) in
+    (("umfpack_di_symbolic","umfpack_di_numeric","umfpack_zi_symbolic","umfpack_zi_numeric",:Int32),
+     ("umfpack_dl_symbolic","umfpack_dl_numeric","umfpack_zl_symbolic","umfpack_zl_numeric",:Int64))
+    @eval begin
+        function umfpack_symbolic!{Tv<:Float64,Ti<:$itype}(U::UmfpackLU{Tv,Ti})
+            if U.symbolic != C_NULL return U end
+            tmp = Array(Ptr{Void},1)
+            status = ccall(($f_sym_r, :libumfpack), Ti,
+                           (Ti, Ti, Ptr{Ti}, Ptr{Ti}, Ptr{Tv}, Ptr{Void},
+                            Ptr{Float64}, Ptr{Float64}),
+                           U.m, U.n, U.colptr, U.rowval, U.nzval, tmp,
+                           umf_ctrl, umf_info)
+            if status != UMFPACK_OK; error("Error code $status from symbolic factorization"); end
+            U.symbolic = tmp[1]
+            finalizer(U.symbolic,umfpack_free_symbolic)
+            U
+        end
+
+        function umfpack_symbolic!{Tv<:Complex128,Ti<:$itype}(U::UmfpackLU{Tv,Ti})
+            if U.symbolic != C_NULL return U end
+            tmp = Array(Ptr{Void},1)
+            status = ccall(($f_sym_r, :libumfpack), Ti,
+                           (Ti, Ti, Ptr{Ti}, Ptr{Ti}, Ptr{Float64}, Ptr{Float64}, Ptr{Void},
+                            Ptr{Float64}, Ptr{Float64}),
+                           U.m, U.n, U.colptr, U.rowval, real(U.nzval), imag(U.nzval), tmp,
+                           umf_ctrl, umf_info)
+            if status != UMFPACK_OK; error("Error code $status from symbolic factorization"); end
+            U.symbolic = tmp[1]
+            finalizer(U.symbolic,umfpack_free_symbolic)
+            U
+        end
+        
+        function umfpack_numeric!{Tv<:Float64,Ti<:$itype}(U::UmfpackLU{Tv,Ti})
+            if U.numeric != C_NULL return U end
+            if U.symbolic == C_NULL umfpack_symbolic!(U) end
+            tmp = Array(Ptr{Void}, 1)
+            status = ccall(($f_num_r, :libumfpack), Ti,
+                           (Ptr{Ti}, Ptr{Ti}, Ptr{Float64}, Ptr{Void}, Ptr{Void}, 
+                            Ptr{Float64}, Ptr{Float64}),
+                           U.colptr, U.rowval, U.nzval, U.symbolic, tmp,
+                           umf_ctrl, umf_info)
+            if status > 0; throw(MatrixIllConditionedException); end
+            if status != UMFPACK_OK; error("Error code $status from numeric factorization"); end
+            U.numeric = tmp[1]
+            finalizer(U.numeric,umfpack_free_numeric)
+            U
+        end
+
+        function umfpack_numeric!{Tv<:Complex128,Ti<:$itype}(U::UmfpackLU{Tv,Ti})
+            if U.numeric != C_NULL return U end
+            if U.symbolic == C_NULL umfpack_symbolic!(U) end
+            tmp = Array(Ptr{Void}, 1)
+            status = ccall(($f_num_r, :libumfpack), Ti,
+                           (Ptr{Ti}, Ptr{Ti}, Ptr{Float64}, Ptr{Float64}, Ptr{Void}, Ptr{Void}, 
+                            Ptr{Float64}, Ptr{Float64}),
+                           U.colptr, U.rowval, real(U.nzval), imag(U.nzval), U.symbolic, tmp,
+                           umf_ctrl, umf_info)
+            if status > 0; throw(MatrixIllConditionedException); end
+            if status != UMFPACK_OK; error("Error code $status from numeric factorization"); end
+            U.numeric = tmp[1]
+            finalizer(U.numeric,umfpack_free_numeric)
+            U
+        end
+    end
 end
 
+for (f_sol_r, f_sol_c, inttype) in
+    (("umfpack_di_solve","umfpack_zi_solve",:Int32),
+     ("umfpack_dl_solve","umfpack_zl_solve",:Int64))
+    @eval begin
+        function umfpack_solve{Tv<:Float64,Ti<:$inttype}(lu::UmfpackLU{Tv,Ti}, b::Vector{Tv}, typ::Integer)
+            umfpack_numeric!(lu)
+            x = similar(b)
+            status = ccall(($f_sol_r, :libumfpack), Ti,
+                           (Ti, Ptr{Ti}, Ptr{Ti}, Ptr{Float64}, Ptr{Float64},
+                            Ptr{Float64}, Ptr{Void}, Ptr{Float64}, Ptr{Float64}),
+                           typ, lu.colptr, lu.rowval, lu.nzval, x, b, lu.numeric, umf_ctrl, umf_info)
+            if status != UMFPACK_OK; error("Error code $status in umfpack_solve"); end
+            return x
+        end
+
+        function umfpack_solve{Tv<:Complex128,Ti<:$inttype}(lu::UmfpackLU{Tv,Ti}, b::Vector{Tv}, typ::Integer)
+            umfpack_numeric!(lu)
+            xr = similar(b, Float64)
+            xi = similar(b, Float64)
+            status = ccall(($f_sol_c, :libumfpack),
+                           Ti,
+                           (Ti, Ptr{Ti}, Ptr{Ti}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, 
+                            Ptr{Float64}, Ptr{Float64}, Ptr{Void}, Ptr{Float64}, Ptr{Float64}),
+                           typ, lu.colptr, lu.rowval, real(lu.nzval), imag(lu.nzval), 
+                           xr, xi, real(b), imag(b), lu.num, umf_ctrl, umf_info)
+            if status != UMFPACK_OK; error("Error code $status from umfpack_solve"); end
+            return complex(xr,xi)
+        end
+    end
+end
+
+umfpack_solve(lu::UmfpackLU, b::Vector) = umfpack_solve(lu, b, UMFPACK_A)
+ 
+## The C functions called by these Julia functions do not depend on
+## the numeric and index types, even though the umfpack names indicate
+## they do.  The umfpack_free_* functions can be called on C_NULL without harm.
+function umfpack_free_symbolic(symb::Ptr{Void})
+    tmp = [symb]
+    ccall((:umfpack_dl_free_symbolic, :libumfpack), Void, (Ptr{Void},), tmp)
+end
+
+function umfpack_free_symbolic(lu::UmfpackLU)
+    if lu.symbolic == C_NULL return lu end
+    umfpack_free_numeric(lu)
+    umfpack_free_symbolic(lu.symbolic)
+    lu.symbolic = C_NULL
+    lu
+end
+
+function umfpack_free_numeric(num::Ptr{Void})
+    tmp = [num]
+    ccall((:umfpack_dl_free_numeric, :libumfpack), Void, (Ptr{Void},), tmp)
+end
+
+function umfpack_free_symbolic(lu::UmfpackLU)
+    if lu.numeric == C_NULL return lu end
+    umfpack_free_numeric(lu.numeric)
+    lu.numeric = C_NULL
+    lu
+end
+ 
+function umfpack_report_symbolic(symb::Ptr{Void}, level::Real)
+    old_prl::Float64 = umf_ctrl[UMFPACK_PRL]
+    umf_ctrl[UMFPACK_PRL] = float64(level)
+    status = ccall((:umfpack_dl_report_symbolic, :libumfpack), Int,
+                   (Ptr{Void}, Ptr{Float64}), symb, umf_ctrl)
+    umf_ctrl[UMFPACK_PRL] = old_prl
+    if status != 0
+        error("Error code $status from umfpack_report_symbolic")
+    end
+end
+
+umfpack_report_symbolic(symb::Ptr{Void}) = umfpack_report_symbolic(symb, 4.)
+
+function umfpack_report_symbolic(lu::UmfpackLU, level::Real)
+    umfpack_report_symbolic(umfpack_symbolic!(lu).symbolic, level)
+end
+
+umfpack_report_symbolic(lu::UmfpackLU) = umfpack_report_symbolic(lu.symbolic,4.)
+function umfpack_report_numeric(num::Ptr{Void}, level::Real)
+    old_prl::Float64 = umf_ctrl[UMFPACK_PRL]
+    umf_ctrl[UMFPACK_PRL] = float64(level)
+    status = ccall((:umfpack_dl_report_numeric, :libumfpack), Int,
+                   (Ptr{Void}, Ptr{Float64}), num, umf_ctrl)
+    umf_ctrl[UMFPACK_PRL] = old_prl
+    if status != 0
+        error("Error code $status from umfpack_report_numeric")
+    end
+end
+
+umfpack_report_numeric(num::Ptr{Void}) = umfpack_report_numeric(num, 4.)
+function umfpack_report_numeric(lu::UmfpackLU, level::Real)
+    umfpack_report_numeric(umfpack_numeric!(lu).symbolic, level)
+end
+
+umfpack_report_numeric(lu::UmfpackLU) = umfpack_report_numeric(lu.symbolic,4.)
+
+## CHOLMOD
+
+const chm_com_sz = ccall((:jl_cholmod_common_size, :libsuitesparse_wrapper), Int, ())
+
+### Need one cholmod_common struct for Int32 and one for Int64
+const chm_com    = Array(Uint8, chm_com_sz)
+ccall((:cholmod_start, :libcholmod), Int32, (Ptr{Uint8},), chm_com)
+
+const chm_l_com  = Array(Uint8, chm_com_sz)
+ccall((:cholmod_l_start, :libcholmod), Int32, (Ptr{Uint8},), chm_l_com)
+
+### These offsets should be reconfigured to be less error-prone in matches
+const chm_com_offsets = Array(Int, 20)
+ccall((:jl_cholmod_common_offsets, :libsuitesparse_wrapper),
+      Void, (Ptr{Uint},), chm_com_offsets)
+                                               
+### A way of examining some of the fields in chm_com or chm_l_com
+type ChmCommon
+    dbound::Float64
+    maxrank::Int
+    supernodal_switch::Float64
+    supernodal::Int32
+    final_asis::Int32
+    final_super::Int32
+    final_ll::Int32
+    final_pack::Int32
+    final_monotonic::Int32
+    final_resymbol::Int32
+    prefer_zomplex::Int32               # should always be false
+    prefer_upper::Int32
+    print::Int32                        # print level. Default: 3
+    precise::Int32                      # print 16 digits, otherwise 5
+    nmethods::Int32                     # number of ordering methods
+    selected::Int32
+    postorder::Int32
+    itype::Int32
+    dtype::Int32
+end
+
+### there must be an easier way but at least this works.
+function ChmCommon(aa::Array{Uint8,1})
+    eval(Expr(:call,
+              unshift!({reinterpret(ChmCommon.types[i],
+                                    aa[(1:sizeof(ChmCommon.types[i])) +
+                                       chm_com_offsets[i]])[1]
+                        for i in 1:length(ChmCommon.types)}, :ChmCommon),
+              Any))
+end
+
+chm_itype{Tv,Ti<:CHMITypes}(S::SparseMatrixCSC{Tv,Ti}) = Ti == Int32 ? CHOLMOD_INT : CHOLMOD_LONG
+chm_xtype{Tv<:CHMVTypes}(S::SparseMatrixCSC{Tv}) = Tv <: Complex ? CHOLMOD_COMPLEX : CHOLMOD_REAL
+chm_dtype{Tv<:CHMVTypes}(S::SparseMatrixCSC{Tv}) = T <: Union(Float32, Complex64) ? CHOLMOD_SINGLE : CHOLMOD_DOUBLE
+
 # Wrapper for memory allocated by CHOLMOD. Carry along the value and index types.
 ## FIXME: CholmodPtr and UmfpackPtr should be amalgamated
 type CholmodPtr{Tv<:CHMVTypes,Ti<:CHMITypes}
@@ -110,7 +359,7 @@ eltype{Tv,Ti}(P::CholmodPtr{Tv,Ti}) = Tv
 indtype{Tv,Ti}(P::CholmodPtr{Tv,Ti}) = Ti
 
 function cholmod_common_finalizer(x::Vector{Ptr{Void}})
-    st = ccall((:cholmod_finish, libcholmod), BlasInt, (Ptr{Void},), x[1])
+    st = ccall((:cholmod_finish, :libcholmod), BlasInt, (Ptr{Void},), x[1])
     if st != CHOLMOD_TRUE error("Error calling cholmod_finish") end
     c_free(x[1])
 end
@@ -119,9 +368,9 @@ type CholmodCommon
     pt::Vector{Ptr{Void}}
     function CholmodCommon()
         pt = Array(Ptr{Void}, 1)
-        ccall((:jl_cholmod_common, libsuitesparse_wrapper), Void,
+        ccall((:jl_cholmod_common, :libsuitesparse_wrapper), Void,
               (Ptr{Void},), pt)
-        st = ccall((:cholmod_start, libcholmod), BlasInt, (Ptr{Void}, ), pt[1])
+        st = ccall((:cholmod_start, :libcholmod), BlasInt, (Ptr{Void}, ), pt[1])
         if st != CHOLMOD_TRUE error("Error calling cholmod_start") end
         finalizer(pt, cholmod_common_finalizer)
         new(pt)
@@ -129,7 +378,7 @@ type CholmodCommon
 end
 
 function show(io::IO, cm::CholmodCommon)
-    st = ccall((:cholmod_print_common, libcholmod), BlasInt,
+    st = ccall((:cholmod_print_common, :libcholmod), BlasInt,
                (Ptr{Uint8},Ptr{Void}), "", cm.pt[1])
     if st != CHOLMOD_TRUE error("Error calling cholmod_print_common") end
 end
@@ -144,7 +393,7 @@ type CholmodSparse{Tv<:CHMVTypes,Ti<:CHMITypes}
         pt = CholmodPtr{Tv,Ti}(Array(Ptr{Void}, 1))
         cp = convert_to_0_based_indexing(S)
         
-        ccall((:jl_cholmod_sparse, libsuitesparse_wrapper), Void,
+        ccall((:jl_cholmod_sparse, :libsuitesparse_wrapper), Void,
               (Ptr{Void}, Uint, Uint, Uint, Ptr{Void}, Ptr{Void}, Ptr{Void},
                Ptr{Void}, Ptr{Void}, BlasInt, BlasInt, BlasInt, BlasInt, BlasInt, Int),
               pt.val, S.m, S.n, nnz(S), cp.colptr, cp.rowval, C_NULL,
@@ -179,7 +428,7 @@ SparseMatrixCSC(cs::CholmodSparse) = convert_to_1_based_indexing(cs.cp)
 ## For testing only.  The infinity and 1 norms of a sparse matrix are simply
 ## the same norm applied to its nzval field.
 function norm(cs::CholmodSparse, p::Number)
-    ccall((:cholmod_norm_sparse, libcholmod), Float64,
+    ccall((:cholmod_norm_sparse, :libcholmod), Float64,
           (Ptr{Void}, BlasInt, Ptr{Void}), cs.pt.val[1], p == Inf ? 0 : 1, cs.cm.pt[1])
 end
 
@@ -188,10 +437,10 @@ norm(cs::CholmodSparse) = norm(cs, Inf)
 ## Approximate minimal degree ordering
 function chm_amd(cs::CholmodSparse)
     aa = Array(BlasInt, cs.cp.m)
-    st = cs.stype == 0 ? ccall(_chm_colamd, BlasInt,
+    st = cs.stype == 0 ? ccall(:cholmod_colamd, BlasInt,
                                (Ptr{Void}, Ptr{Void}, Uint, BlasInt, Ptr{BlasInt}, Ptr{Void}),
                                cs.pt.val[1], C_NULL, 0, 1, aa, cs.cm.pt[1]) :
-                         ccall(_chm_amd, BlasInt, (Ptr{Void}, Ptr{Void}, Uint, Ptr{BlasInt}, Ptr{Void}),
+                         ccall(:cholmod_amd, BlasInt, (Ptr{Void}, Ptr{Void}, Uint, Ptr{BlasInt}, Ptr{Void}),
                                cs.pt.val[1], C_NULL, 0, aa, cs.cm.pt[1])
     if st != CHOLMOD_TRUE error("Error in cholmod_amd") end
     aa
@@ -208,7 +457,7 @@ type CholmodFactor{Tv<:CHMVTypes,Ti<:CHMITypes} <: Factorization{Tv}
 end
 
 function cholmod_factor_finalizer(x::CholmodFactor)
-    if ccall(_chm_free_fa, BlasInt, (Ptr{Void}, Ptr{Void}), x.pt.val, x.cs.cm[1]) != CHOLMOD_TRUE
+    if ccall((:cholmod_free_factor, :libcholmod), BlasInt, (Ptr{Void}, Ptr{Void}), x.pt.val, x.cs.cm[1]) != CHOLMOD_TRUE
         error("CHOLMOD error in cholmod_free_factor")
     end
 end
@@ -223,16 +472,16 @@ indtype{Tv,Ti}(F::CholmodFactor{Tv,Ti}) = Ti
 
 function CholmodFactor{Tv,Ti}(cs::CholmodSparse{Tv,Ti})
     pt = CholmodPtr{Tv,Ti}(Array(Ptr{Void}, 1))
-    pt.val[1] = ccall(_chm_analyze, Ptr{Void},
+    pt.val[1] = ccall(:cholmod_analyze, Ptr{Void},
                       (Ptr{Void}, Ptr{Void}), cs.pt.val[1], cs.cm.pt[1])
-    st = ccall(_chm_factorize, BlasInt,
+    st = ccall(:cholmod_factorize, BlasInt,
                (Ptr{Void}, Ptr{Void}, Ptr{Void}), cs.pt.val[1], pt.val[1], cs.cm.pt[1])
     if st != CHOLMOD_TRUE error("CHOLMOD failure in factorize") end
     CholmodFactor{Tv,Ti}(pt, cs)
 end
 
 function show(io::IO, cf::CholmodFactor)
-    st = ccall(_chm_print_fa, BlasInt, (Ptr{Void}, Ptr{Uint8}, Ptr{Void}), cf.pt.val[1], "", cf.cs.cm.pt[1])
+    st = ccall(:cholmod_print_fa, BlasInt, (Ptr{Void}, Ptr{Uint8}, Ptr{Void}), cf.pt.val[1], "", cf.cs.cm.pt[1])
     if st != CHOLMOD_TRUE error("Cholmod error in print_factor") end
 end
 
@@ -253,7 +502,7 @@ function CholmodDense{T<:CHMVTypes}(b::VecOrMat{T}, cm::CholmodCommon)
 
     pt = Array(Ptr{Void}, 1)
 
-    ccall((:jl_cholmod_dense, libsuitesparse_wrapper), Void,
+    ccall((:jl_cholmod_dense, :libsuitesparse_wrapper), Void,
           (Ptr{Void}, Uint, Uint, Uint, Uint, Ptr{Void}, Ptr{Void}, BlasInt, Int),
           pt, m, n, length(b), m, b, C_NULL, xtype, dtype)
     finalizer(pt, x->c_free(pt[1]))
@@ -265,7 +514,7 @@ CholmodDense{T<:Integer}(B::VecOrMat{T}, cm::CholmodCommon) = CholmodDense(float
 size(cd::CholmodDense) = (cd.m, cd.n)
 
 function show(io::IO, cd::CholmodDense)
-    st = ccall(_chm_print_dn, BlasInt, (Ptr{Void},Ptr{Uint8},Ptr{Void}), cd.pt[1], "", cd.cm.pt[1])
+    st = ccall(:cholmod_print_dn, BlasInt, (Ptr{Void},Ptr{Uint8},Ptr{Void}), cd.pt[1], "", cd.cm.pt[1])
     if st != CHOLMOD_TRUE error("Cholmod error in print_dense") end
 end
 
@@ -282,7 +531,7 @@ type CholmodDenseOut{Tv<:CHMVTypes,Ti<:CHMITypes}
 end
 
 function cholmod_denseout_finalizer(cd::CholmodDenseOut)
-    st = ccall(_chm_free_dn, BlasInt, (Ptr{Void}, Ptr{Void}), cd.pt.val, cd.cm.pt[1])
+    st = ccall((:cholmod_free_dense, :libcholmod), BlasInt, (Ptr{Void}, Ptr{Void}), cd.pt.val, cd.cm.pt[1])
     if st != CHOLMOD_TRUE error("Error in cholmod_free_dense") end
 end
 
@@ -292,7 +541,7 @@ size(cd::CholmodDenseOut) = (cd.m, cd.n)
 
 function convert{T}(::Type{Array{T}}, cdo::CholmodDenseOut{T})
     mm = Array(T, size(cdo))
-    ccall((:jl_cholmod_dense_copy_out, libsuitesparse_wrapper), Void,
+    ccall((:jl_cholmod_dense_copy_out, :libsuitesparse_wrapper), Void,
           (Ptr{Void}, Ptr{T}), cdo.pt.val[1], mm)
     mm
 end
@@ -300,7 +549,7 @@ end
 function solve{Tv,Ti}(cf::CholmodFactor{Tv,Ti}, B::CholmodDense{Tv}, solv::Integer)
     m, n = size(B)
     cdo = CholmodPtr{Tv,Ti}(Array(Ptr{Void},1))
-    cdo.val[1] = ccall(_chm_solve, Ptr{Void},
+    cdo.val[1] = ccall(:cholmod_solve, Ptr{Void},
                        (BlasInt, Ptr{Void}, Ptr{Void}, Ptr{Void}),
                        solv, cf.pt.val[1], B.pt[1], cf.cs.cm.pt[1])
     return cdo, m, n, cf.cs.cm
@@ -324,13 +573,13 @@ type CholmodSparseOut{Tv<:CHMVTypes,Ti<:CHMITypes}
 end
 
 function cholmod_sparseout_finalizer(cso::CholmodSparseOut)
-    st = ccall(_chm_free_sp, BlasInt,
+    st = ccall((:cholmod_free_sparse, :libcholmod), BlasInt,
                (Ptr{Void}, Ptr{Void}), cso.pt.val, cso.cm.pt[1])
     if st != CHOLMOD_TRUE error("Error in cholmod_free_sparse") end
 end
 
 function nnz(cso::CholmodSparseOut)
-    ccall((:cholmod_nnz, libcholmod), BlasInt,
+    ccall((:cholmod_nnz, :libcholmod), BlasInt,
           (Ptr{Void}, Ptr{Void}), cso.pt.val[1], cso.cm.pt[1])
 end
 size(cso::CholmodSparseOut) = (cso.m, cso.n)
@@ -340,7 +589,7 @@ indtype{Tv,Ti}(cso::CholmodSparseOut{Tv,Ti}) = Ti
 function solve{Tv,Ti}(cf::CholmodFactor{Tv,Ti}, B::CholmodSparse{Tv,Ti}, solv::Integer)
     m, n = size(B)
     cso = CholmodPtr{Tv,Ti}(Array(Ptr{Void},1))
-    cso.val[1] = ccall((:cholmod_spsolve, libcholmod), Ptr{Void},
+    cso.val[1] = ccall((:cholmod_spsolve, :libcholmod), Ptr{Void},
                        (BlasInt, Ptr{Void}, Ptr{Void}, Ptr{Void}),
                        solv, cf.pt.val[1], B.pt[1], B.cm.pt[1])
     CholmodSparseOut{Tv,Ti}(cso, m, n, cf.cs.cm)
@@ -349,7 +598,7 @@ end
 function CholmodSparseOut{Tv,Ti}(cf::CholmodFactor{Tv,Ti})
     n = size(cf.cs)[1]
     cso = CholmodPtr{Tv,Ti}(Array(Ptr{Void},1))
-    cso.val[1] = ccall((:cholmod_factor_to_sparse, libcholmod), Ptr{Void},
+    cso.val[1] = ccall((:cholmod_factor_to_sparse, :libcholmod), Ptr{Void},
                        (Ptr{Void}, Ptr{Void}), cf.pt.val[1], cf.cs.cm.pt[1])
     CholmodSparseOut{Tv,Ti}(cso, n, n, cf.cs.cm)
 end
@@ -357,7 +606,7 @@ end
 function SparseMatrixCSC{Tv,Ti}(cso::CholmodSparseOut{Tv,Ti})
     nz = nnz(cso)
     sp = SparseMatrixCSC{Tv,Ti}(cso.m, cso.n, Array(Ti, cso.n + 1), Array(Ti, nz), Array(Tv, nz))
-    st = ccall((:jl_cholmod_sparse_copy_out, libsuitesparse_wrapper), BlasInt,
+    st = ccall((:jl_cholmod_sparse_copy_out, :libsuitesparse_wrapper), BlasInt,
                 (Ptr{Void}, Ptr{Ti}, Ptr{Ti}, Ptr{Tv}),
                 cso.pt.val[1], sp.colptr, sp.rowval, sp.nzval)
     if st == 1 error("CholmodSparseOut object is not packed") end
@@ -367,22 +616,23 @@ function SparseMatrixCSC{Tv,Ti}(cso::CholmodSparseOut{Tv,Ti})
 end
 
 function show(io::IO, cso::CholmodSparseOut)
-    sp = ccall(_chm_print_sp, BlasInt, (Ptr{Void}, Ptr{Uint8},Ptr{Void}), cso.pt.val[1], "", cso.cm.pt[1])
+    sp = ccall(:cholmod_print_sp, BlasInt, (Ptr{Void}, Ptr{Uint8},Ptr{Void}), cso.pt.val[1], "", cso.cm.pt[1])
     if sp != CHOLMOD_TRUE error("Cholmod error in print_sparse") end
 end
 
 function chm_aat{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, symm::Bool)
     cs = CholmodSparse(A, 0)
     aa = CholmodPtr{Tv,Ti}(Array(Ptr{Void},1))
-    aa.val[1] = ccall(_chm_aat, Ptr{Void}, (Ptr{Void},Ptr{BlasInt},BlasInt,BlasInt,Ptr{Void}),
+    aa.val[1] = ccall(:cholmod_aat, Ptr{Void}, (Ptr{Void},Ptr{BlasInt},BlasInt,BlasInt,Ptr{Void}),
                       cs.pt.val[1], C_NULL, 0, 1, cs.cm.pt[1])
-    if ccall(_chm_sort, BlasInt, (Ptr{Void}, Ptr{Void}), aa.val[1], cs.cm.pt[1]) != CHOLMOD_TRUE
+    if ccall(:cholmod_sort, BlasInt, (Ptr{Void}, Ptr{Void}), aa.val[1], cs.cm.pt[1]) != CHOLMOD_TRUE
         error("Cholmod error in sort")
     end
     if symm
-        pt = ccall(_chm_copy, Ptr{Void}, (Ptr{Void}, BlasInt, BlasInt, Ptr{Void}),
+        pt = ccall(:cholmod_copy, Ptr{Void}, (Ptr{Void}, BlasInt, BlasInt, Ptr{Void}),
                    aa.val[1], 1, 1, cs.cm.pt[1])
-        if ccall(_chm_free_sp, BlasInt, (Ptr{Void}, Ptr{Void}), aa.val, cs.cm.pt[1]) != CHOLMOD_TRUE
+        if ccall((:cholmod_free_sparse, :libcholmod), BlasInt,
+                 (Ptr{Void}, Ptr{Void}), aa.val, cs.cm.pt[1]) != CHOLMOD_TRUE
             error("Cholmod error in free_sparse")
         end
         aa.val[1] = pt
@@ -408,7 +658,7 @@ function cholmod_sparse{Tv,Ti}(S::SparseMatrixCSC{Tv,Ti}, stype::Int)
     if     Tv == Float64 || Tv == Complex128; dtype = CHOLMOD_DOUBLE; 
     elseif Tv == Float32 || Tv == Complex64 ; dtype = CHOLMOD_SINGLE; end
 
-    ccall((:jl_cholmod_sparse, libsuitesparse_wrapper),
+    ccall((:jl_cholmod_sparse, :libsuitesparse_wrapper),
           Ptr{Void},
           (Ptr{Void}, BlasInt, BlasInt, BlasInt, Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void},
            BlasInt, BlasInt, BlasInt, BlasInt, BlasInt, Int),
@@ -432,17 +682,16 @@ function cholmod_dense{T}(B::VecOrMat{T})
     if     T == Float64 || T == Complex128; dtype = CHOLMOD_DOUBLE; 
     elseif T == Float32 || T == Complex64 ; dtype = CHOLMOD_SINGLE; end
 
-    ccall((:jl_cholmod_dense, libsuitesparse_wrapper),
-          Ptr{Void},
+    ccall((:jl_cholmod_dense, :libsuitesparse_wrapper), Ptr{Void},
           (Ptr{Void}, BlasInt, BlasInt, BlasInt, BlasInt, Ptr{T}, Ptr{Void}, BlasInt, Int),
-          cd, m, n, length(B), m, B, C_NULL, xtype, dtype
+          cd, size(B,1), size(B,2), length(B), stride(B,2), B, C_NULL, xtype, dtype
           )
 
     return cd
 end
 
 function cholmod_dense_copy_out{T}(x::Ptr{Void}, sol::VecOrMat{T})
-    ccall((:jl_cholmod_dense_copy_out, libsuitesparse_wrapper),
+    ccall((:jl_cholmod_dense_copy_out, :libsuitesparse_wrapper),
           Void,
           (Ptr{Void}, Ptr{T}),
           x, sol
@@ -469,291 +718,18 @@ function cholmod_transpose_unsym{Tv,Ti}(S::SparseMatrixCSC{Tv,Ti}, cm::Array{Ptr
     return S_t
 end
 
-function cholmod_analyze{Tv<:Union(Float64,Complex128), Ti<:CHMITypes}(cs::Array{Ptr{Void},1}, cm::Array{Ptr{Void},1})
-    ccall(_chm_analyze, Ptr{Void}, (Ptr{Void}, Ptr{Void}), cs[1], cm[1])
+function cholmod_analyze{Tv<:CHMVTypes, Ti<:CHMITypes}(cs::Array{Ptr{Void},1}, cm::Array{Ptr{Void},1})
+    ccall(:cholmod_analyze, Ptr{Void}, (Ptr{Void}, Ptr{Void}), cs[1], cm[1])
 end
 
-function cholmod_factorize{Tv<:Union(Float64,Complex128), Ti<:CHMITypes}(cs::Array{Ptr{Void},1}, cs_factor::Ptr{Void}, cm::Array{Ptr{Void},1})
-    st = ccall(_chm_factorize, BlasInt, (Ptr{Void}, Ptr{Void}, Ptr{Void}), cs[1], cs_factor, cm[1])
+function cholmod_factorize{Tv<:CHMVTypes, Ti<:CHMITypes}(cs::Array{Ptr{Void},1}, cs_factor::Ptr{Void}, cm::Array{Ptr{Void},1})
+    st = ccall(:cholmod_factorize, BlasInt, (Ptr{Void}, Ptr{Void}, Ptr{Void}), cs[1], cs_factor, cm[1])
     if st != CHOLMOD_TRUE error("CHOLMOD could not factorize the matrix") end
 end
 
 function cholmod_solve(cs_factor::Ptr{Void}, cd_rhs::Array{Ptr{Void},1}, cm::Array{Ptr{Void},1})
-    ccall(_chm_solve, Ptr{Void}, (BlasInt, Ptr{Void}, Ptr{Void}, Ptr{Void}),
+    ccall(:cholmod_solve, Ptr{Void}, (BlasInt, Ptr{Void}, Ptr{Void}, Ptr{Void}),
           CHOLMOD_A, cs_factor, cd_rhs[1], cm[1])
 end
 
-## UMFPACK
-
-# Wrapper for memory allocated by umfpack. Carry along the value and index types.
-type UmfpackPtr{Tv<:Union(Float64,Complex128),Ti<:CHMITypes}
-    val::Vector{Ptr{Void}} 
-end
-
-type UmfpackLU{Tv<:Union(Float64,Complex128),Ti<:CHMITypes} <: Factorization{Tv}
-    numeric::UmfpackPtr{Tv,Ti}
-    mat::SparseMatrixCSC{Tv,Ti}
-end
-
-function show(io::IO, f::UmfpackLU)
-    @printf(io, "UMFPACK LU Factorization of a %d-by-%d sparse matrix\n",
-           size(f.mat,1), size(f.mat,2))
-    println(f.numeric)
-    umfpack_report(f)
-end
-
-type UmfpackLUTrans{Tv<:Union(Float64,Complex128),Ti<:CHMITypes} <: Factorization{Tv}
-    numeric::UmfpackPtr{Tv,Ti}
-    mat::SparseMatrixCSC{Tv,Ti}
-end
-
-function show(io::IO, f::UmfpackLUTrans)
-    @printf(io, "UMFPACK LU Factorization of a transposed %d-by-%d sparse matrix\n",
-           size(f.mat,1), size(f.mat,2))
-    println(f.numeric)
-    umfpack_report(f)
-end
-
-function UmfpackLU{Tv<:Union(Float64,Complex128),Ti<:CHMITypes}(S::SparseMatrixCSC{Tv,Ti})
-    Scopy = copy(S) 
-    Scopy = convert_to_0_based_indexing!(Scopy)
-    numeric = []
-
-    try
-        symbolic = umfpack_symbolic(Scopy)
-        numeric = umfpack_numeric(Scopy, symbolic)
-    catch e
-        if is(e,MatrixIllConditionedException)
-            error("Input matrix is ill conditioned or singular");
-        else
-            error("Error calling UMFPACK")
-        end
-    end
-    
-    return UmfpackLU(numeric,Scopy) 
-end
-
-function UmfpackLU!{Tv<:Union(Float64,Complex128),Ti<:CHMITypes}(S::SparseMatrixCSC{Tv,Ti})
-    Sshallow = SparseMatrixCSC(S.m,S.n,S.colptr,S.rowval,S.nzval)
-    Sshallow = convert_to_0_based_indexing!(Sshallow)
-    numeric = []
-
-    try
-        symbolic = umfpack_symbolic(Sshallow)
-        numeric = umfpack_numeric(Sshallow, symbolic)
-    catch e
-        Sshallow = convert_to_1_based_indexing!(Sshallow)
-        if is(e,MatrixIllConditionedException)
-            error("Input matrix is ill conditioned or singular");
-        else
-            error("Error calling UMFPACK")
-        end
-    end
-
-    S.rowval = []
-    S.nzval = []
-    S.colptr = ones(S.n+1)
-    
-    return UmfpackLU(numeric,Sshallow)
-end
-
-function UmfpackLUTrans(S::SparseMatrixCSC) 
-    x = UmfpackLU(S)
-    return UmfpackLUTrans(x.numeric, x.mat)
-end
-
-# Solve with Factorization
-
-(\){T}(fact::UmfpackLU{T}, b::Vector) = fact \ convert(Array{T,1}, b)
-(\){T}(fact::UmfpackLU{T}, b::Vector{T}) = umfpack_solve(fact.mat,b,fact.numeric)
-
-(\){T}(fact::UmfpackLUTrans{T}, b::Vector) = fact \ convert(Array{T,1}, b)
-(\){T}(fact::UmfpackLUTrans{T}, b::Vector{T}) = umfpack_transpose_solve(fact.mat,b,fact.numeric)
-
-ctranspose(fact::UmfpackLU) = UmfpackLUTrans(fact.numeric, fact.mat)
-
-# Solve directly with matrix
-
-(\)(S::SparseMatrixCSC, b::Vector) = UmfpackLU(S) \ b
-At_ldiv_B(S::SparseMatrixCSC, b::Vector) = UmfpackLUTrans(S) \ b
-Ac_ldiv_B(S::SparseMatrixCSC, b::Vector) = UmfpackLUTrans(S) \ b
-
-## Wrappers around UMFPACK routines
-
-for (f_sym_r, f_sym_c, inttype) in
-    (("umfpack_di_symbolic","umfpack_zi_symbolic",:Int32),
-     ("umfpack_dl_symbolic","umfpack_zl_symbolic",:Int64))
-    @eval begin
-
-        function umfpack_symbolic{Tv<:Float64,Ti<:$inttype}(S::SparseMatrixCSC{Tv,Ti})
-            # Pointer to store the symbolic factorization returned by UMFPACK
-            Symbolic = UmfpackPtr{Tv,Ti}(Array(Ptr{Void},1))
-            status = ccall(($f_sym_r, libumfpack),
-                           Ti,
-                           (Ti, Ti, 
-                            Ptr{Ti}, Ptr{Ti}, Ptr{Float64}, Ptr{Void}, Ptr{Float64}, Ptr{Float64}),
-                           S.m, S.n, 
-                           S.colptr, S.rowval, S.nzval, Symbolic.val, C_NULL, C_NULL)
-            if status != UMFPACK_OK; error("Error in symbolic factorization"); end
-            finalizer(Symbolic,umfpack_free_symbolic)
-            return Symbolic
-        end
-
-        function umfpack_symbolic{Tv<:Complex128,Ti<:$inttype}(S::SparseMatrixCSC{Tv,Ti})
-            # Pointer to store the symbolic factorization returned by UMFPACK
-            Symbolic = UmfpackPtr{Tv,Ti}(Array(Ptr{Void},1))
-            status = ccall(($f_sym_c, libumfpack),
-                           Ti,
-                           (Ti, Ti, 
-                            Ptr{Ti}, Ptr{Ti}, Ptr{Float64}, Ptr{Float64}, Ptr{Void}, 
-                            Ptr{Float64}, Ptr{Float64}),
-                           S.m, S.n, 
-                           S.colptr, S.rowval, real(S.nzval), imag(S.nzval), Symbolic.val, 
-                           C_NULL, C_NULL)
-            if status != UMFPACK_OK; error("Error in symbolic factorization"); end
-            finalizer(Symbolic,umfpack_free_symbolic) # Check: do we need to free if there was an error?
-            return Symbolic
-        end
-
-    end
-end
-
-for (f_num_r, f_num_c, inttype) in
-    (("umfpack_di_numeric","umfpack_zi_numeric",:Int32),
-     ("umfpack_dl_numeric","umfpack_zl_numeric",:Int64))
-    @eval begin
-
-        function umfpack_numeric{Tv<:Float64,Ti<:$inttype}(S::SparseMatrixCSC{Tv,Ti}, Symbolic)
-            # Pointer to store the numeric factorization returned by UMFPACK
-            Numeric = UmfpackPtr{Tv,Ti}(Array(Ptr{Void},1))
-            status = ccall(($f_num_r, libumfpack),
-                           Ti,
-                           (Ptr{Ti}, Ptr{Ti}, Ptr{Float64}, Ptr{Void}, Ptr{Void}, 
-                            Ptr{Float64}, Ptr{Float64}),
-                           S.colptr, S.rowval, S.nzval, Symbolic.val[1], Numeric.val, 
-                           C_NULL, C_NULL)
-            if status > 0; throw(MatrixIllConditionedException); end
-            if status != UMFPACK_OK; error("Error in numeric factorization"); end
-            finalizer(Numeric,umfpack_free_numeric)
-            return Numeric
-        end
-
-        function umfpack_numeric{Tv<:Complex128,Ti<:$inttype}(S::SparseMatrixCSC{Tv,Ti}, Symbolic)
-            # Pointer to store the numeric factorization returned by UMFPACK
-            Numeric = UmfpackPtr{Tv,Ti}(Array(Ptr{Void},1))
-            status = ccall(($f_num_c, libumfpack),
-                           Ti,
-                           (Ptr{Ti}, Ptr{Ti}, Ptr{Float64}, Ptr{Float64}, Ptr{Void}, Ptr{Void}, 
-                            Ptr{Float64}, Ptr{Float64}),
-                           S.colptr, S.rowval, real(S.nzval), imag(S.nzval), Symbolic.val[1], Numeric.val, 
-                           C_NULL, C_NULL)
-            if status > 0; throw(MatrixIllConditionedException); end
-            if status != UMFPACK_OK; error("Error in numeric factorization"); end
-            finalizer(Numeric,umfpack_free_numeric)
-            return Numeric
-        end
-
-    end
-end
-
-for (f_sol_r, f_sol_c, inttype) in
-    (("umfpack_di_solve","umfpack_zi_solve",:Int32),
-     ("umfpack_dl_solve","umfpack_zl_solve",:Int64))
-    @eval begin
-
-        function umfpack_solve{Tv<:Float64,Ti<:$inttype}(S::SparseMatrixCSC{Tv,Ti}, 
-                                                             b::Vector{Tv}, Numeric::UmfpackPtr{Tv,Ti})
-            x = similar(b)
-            status = ccall(($f_sol_r, libumfpack),
-                           Ti,
-                           (Ti, Ptr{Ti}, Ptr{Ti}, Ptr{Float64}, 
-                            Ptr{Float64}, Ptr{Float64}, Ptr{Void}, Ptr{Float64}, Ptr{Float64}),
-                           UMFPACK_A, S.colptr, S.rowval, S.nzval, 
-                           x, b, Numeric.val[1], C_NULL, C_NULL)
-            if status != UMFPACK_OK; error("Error in solve"); end
-            return x
-        end
-
-        function umfpack_solve{Tv<:Complex128,Ti<:$inttype}(S::SparseMatrixCSC{Tv,Ti}, 
-                                                                b::Vector{Tv}, Numeric::UmfpackPtr{Tv,Ti})
-            xr = similar(b, Float64)
-            xi = similar(b, Float64)
-            status = ccall(($f_sol_c, libumfpack),
-                           Ti,
-                           (Ti, Ptr{Ti}, Ptr{Ti}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, 
-                            Ptr{Float64}, Ptr{Float64}, Ptr{Void}, Ptr{Float64}, Ptr{Float64}),
-                           UMFPACK_A, S.colptr, S.rowval, real(S.nzval), imag(S.nzval), 
-                           xr, xi, real(b), imag(b), Numeric.val[1], C_NULL, C_NULL)
-            if status != UMFPACK_OK; error("Error in solve"); end
-            return complex(xr,xi)
-        end
-
-        function umfpack_transpose_solve{Tv<:Float64,Ti<:$inttype}(S::SparseMatrixCSC{Tv,Ti}, 
-                                                             b::Vector{Tv}, Numeric::UmfpackPtr{Tv,Ti})
-            x = similar(b)
-            status = ccall(($f_sol_r, libumfpack),
-                           Ti,
-                           (Ti, Ptr{Ti}, Ptr{Ti}, Ptr{Float64}, 
-                            Ptr{Float64}, Ptr{Float64}, Ptr{Void}, Ptr{Float64}, Ptr{Float64}),
-                           UMFPACK_At, S.colptr, S.rowval, S.nzval, 
-                           x, b, Numeric.val[1], C_NULL, C_NULL)
-            if status != UMFPACK_OK; error("Error in solve"); end
-            return x
-        end
-
-        function umfpack_transpose_solve{Tv<:Complex128,Ti<:$inttype}(S::SparseMatrixCSC{Tv,Ti}, 
-                                                                b::Vector{Tv}, Numeric::UmfpackPtr{Tv,Ti})
-            xr = similar(b, Float64)
-            xi = similar(b, Float64)
-            status = ccall(($f_sol_c, libumfpack),
-                           Ti,
-                           (Ti, Ptr{Ti}, Ptr{Ti}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, 
-                            Ptr{Float64}, Ptr{Float64}, Ptr{Void}, Ptr{Float64}, Ptr{Float64}),
-                           UMFPACK_At, S.colptr, S.rowval, real(S.nzval), imag(S.nzval), 
-                           xr, xi, real(b), imag(b), Numeric.val[1], C_NULL, C_NULL)
-            if status != UMFPACK_OK; error("Error in solve"); end
-            return complex(xr,xi)
-        end
-
-    end
-end
-
-for (f_report, elty, inttype) in
-    (("umfpack_di_report_numeric", :Float64,    :Int),
-     ("umfpack_zi_report_numeric", :Complex128, :Int),
-     ("umfpack_dl_report_numeric", :Float64,    :Int64),
-     ("umfpack_zl_report_numeric", :Complex128, :Int64))
-     @eval begin
-
-         function umfpack_report{Tv<:$elty,Ti<:$inttype}(slu::UmfpackLU{Tv,Ti})
-
-             control = zeros(Float64, UMFPACK_CONTROL)
-             control[UMFPACK_PRL] = 4
-         
-             ccall(($f_report, libumfpack),
-                   Ti,
-                   (Ptr{Void}, Ptr{Float64}),
-                   slu.numeric.val[1], control)
-         end
-
-    end
-end
-
-
-for (f_symfree, f_numfree, elty, inttype) in
-    (("umfpack_di_free_symbolic","umfpack_di_free_numeric",:Float64,:Int32),
-     ("umfpack_zi_free_symbolic","umfpack_zi_free_numeric",:Complex128,:Int32),
-     ("umfpack_dl_free_symbolic","umfpack_dl_free_numeric",:Float64,:Int64),
-     ("umfpack_zl_free_symbolic","umfpack_zl_free_numeric",:Complex128,:Int64))
-    @eval begin
-
-        umfpack_free_symbolic{Tv<:$elty,Ti<:$inttype}(Symbolic::UmfpackPtr{Tv,Ti}) =
-        ccall(($f_symfree, libumfpack), Void, (Ptr{Void},), Symbolic.val)
-        
-        umfpack_free_numeric{Tv<:$elty,Ti<:$inttype}(Numeric::UmfpackPtr{Tv,Ti}) =
-        ccall(($f_numfree, libumfpack), Void, (Ptr{Void},), Numeric.val)
-        
-    end
-end
-
 end #module
diff --git a/test/suitesparse.jl b/test/suitesparse.jl
index c2756429c93f1..d9ce338f15a07 100644
--- a/test/suitesparse.jl
+++ b/test/suitesparse.jl
@@ -5,3 +5,4 @@ using SuiteSparse
 se33 = speye(3)
 do33 = ones(3)
 @test isequal(se33 \ do33, do33)
+

From c5b0fe03805b0a0c9360f1123b35894c36fb3bd5 Mon Sep 17 00:00:00 2001
From: "Viral B. Shah" <viral@mayin.org>
Date: Wed, 6 Feb 2013 19:26:12 -0500
Subject: [PATCH 05/11] Rename lud to lu. Update docs.

---
 base/exports.jl      |  3 +--
 base/linalg_dense.jl | 22 +++++++++-------------
 doc/stdlib/base.rst  | 16 ++++++++++------
 test/linalg.jl       |  8 +++-----
 4 files changed, 23 insertions(+), 26 deletions(-)

diff --git a/base/exports.jl b/base/exports.jl
index a38d9fa1fea28..e194b1ef28fad 100644
--- a/base/exports.jl
+++ b/base/exports.jl
@@ -553,8 +553,7 @@ export
     ldltd,
     linreg,
     lu,
-    lud,
-    lud!,
+    lu!,
     norm,
     normfro,
     null,
diff --git a/base/linalg_dense.jl b/base/linalg_dense.jl
index 068cc0ef3e3fc..e4f0c123307f3 100644
--- a/base/linalg_dense.jl
+++ b/base/linalg_dense.jl
@@ -536,17 +536,13 @@ function factors{T}(lu::LUDense{T})
     L, U, P
 end
 
-function lud!{T<:BlasFloat}(A::Matrix{T})
+function lu!{T<:BlasFloat}(A::Matrix{T})
     lu, ipiv, info = LAPACK.getrf!(A)
     LUDense{T}(lu, ipiv, info)
 end
 
-lud{T<:BlasFloat}(A::Matrix{T}) = lud!(copy(A))
-lud{T<:Number}(A::Matrix{T}) = lud(float64(A))
-
-## Matlab-compatible
-lu{T<:Number}(A::Matrix{T}) = factors(lud(A))
-lu(x::Number) = (one(x), x, [1])
+lu{T<:BlasFloat}(A::Matrix{T}) = lu!(copy(A))
+lu{T<:Number}(A::Matrix{T}) = lu(float64(A))
 
 function det{T}(lu::LUDense{T})
     m, n = size(lu)
@@ -680,7 +676,7 @@ function det(A::Matrix)
     m, n = size(A)
     if m != n; error("det only defined for square matrices"); end
     if istriu(A) | istril(A); return prod(diag(A)); end
-    return det(lud(A))
+    return det(lu(A))
 end
 det(x::Number) = x
 
@@ -894,7 +890,7 @@ function cond(A::StridedMatrix, p)
     elseif p == 1 || p == Inf
         m, n = size(A)
         if m != n; error("Use 2-norm for non-square matrices"); end
-        cnd = 1 / LAPACK.gecon!(p == 1 ? '1' : 'I', lud(A).lu, norm(A, p))
+        cnd = 1 / LAPACK.gecon!(p == 1 ? '1' : 'I', lu(A).lu, norm(A, p))
     else
         error("Norm type must be 1, 2 or Inf")
     end
@@ -1190,16 +1186,16 @@ end
 
 #show(io, lu::LUTridiagonal) = print(io, "LU decomposition of ", summary(lu.lu))
 
-lud!{T}(A::Tridiagonal{T}) = LUTridiagonal{T}(LAPACK.gttrf!(A.dl,A.d,A.du)...)
-lud{T}(A::Tridiagonal{T}) = LUTridiagonal{T}(LAPACK.gttrf!(copy(A.dl),copy(A.d),copy(A.du))...)
-lu(A::Tridiagonal) = factors(lud(A))
+lu!{T}(A::Tridiagonal{T}) = LUTridiagonal{T}(LAPACK.gttrf!(A.dl,A.d,A.du)...)
+lu{T}(A::Tridiagonal{T}) = LUTridiagonal{T}(LAPACK.gttrf!(copy(A.dl),copy(A.d),copy(A.du))...)
+lu(A::Tridiagonal) = factors(lu(A))
 
 function det{T}(lu::LUTridiagonal{T})
     n = length(lu.d)
     prod(lu.d) * (bool(sum(lu.ipiv .!= 1:n) % 2) ? -one(T) : one(T))
 end
 
-det(A::Tridiagonal) = det(lud(A))
+det(A::Tridiagonal) = det(lu(A))
 
 (\){T<:BlasFloat}(lu::LUTridiagonal{T}, B::StridedVecOrMat{T}) =
     LAPACK.gttrs!('N', lu.dl, lu.d, lu.du, lu.du2, lu.ipiv, copy(B))
diff --git a/doc/stdlib/base.rst b/doc/stdlib/base.rst
index d4b4ce8d6a967..27bf5ecd59d9d 100644
--- a/doc/stdlib/base.rst
+++ b/doc/stdlib/base.rst
@@ -1716,13 +1716,17 @@ Linear algebra functions in Julia are largely implemented by calling functions f
 
    Compute the norm of a ``Vector`` or a ``Matrix``
 
-.. function:: lu(A) -> LU
+.. function:: lu(A)
 
-   Compute LU factorization. LU is an "LU factorization" type that can be used as an ordinary matrix.
+   Compute the LU factorization of ``A`` and return a ``LUDense`` object. ``factors(lu(A))`` returns the triangular matrices containing the factorization. The following functions are available for ``LUDense`` objects: ``size``, ``factors``, ``\``, ``inv``, ``det``.
+
+.. function:: lu!(A)
+
+   ``lu!`` is the same as ``lu``, but overwrites the input matrix A with the factorization.
 
 .. function:: chol(A, [LU])
 
-   Compute the Cholesky factorization of ``A`` and return a ``CholeskyDense`` object. ``LU`` may be 'L' for using the lower part or 'U' for the upper part. The default is to use 'U'. The ``factors(chol(A))`` returns the triangular matrix containing the factorization. The following functions are available for ``CholeskyDense`` objects: ``size``, ``factors``, ``\``, ``inv``, ``det``.
+   Compute the Cholesky factorization of ``A`` and return a ``CholeskyDense`` object. ``LU`` may be 'L' for using the lower part or 'U' for the upper part. The default is to use 'U'. ``factors(chol(A))`` returns the triangular matrix containing the factorization. The following functions are available for ``CholeskyDense`` objects: ``size``, ``factors``, ``\``, ``inv``, ``det``. A ``LAPACK.PosDefException`` error is thrown in case the matrix is not positive definite.
 
 .. function:: chol!(A, [LU])
 
@@ -1730,7 +1734,7 @@ Linear algebra functions in Julia are largely implemented by calling functions f
 
 .. function:: cholpivot(A, [LU])
 
-   Compute the pivoted Cholesky factorization of ``A`` and return a ``CholeskyDensePivoted`` object. ``LU`` may be 'L' for using the lower part or 'U' for the upper part. The default is to use 'U'. The ``factors(cholpivot(A))`` returns the triangular matrix containing the factorization. The following functions are available for ``CholeskyDensePivoted`` objects: ``size``, ``factors``, ``\``, ``inv``, ``det``.
+   Compute the pivoted Cholesky factorization of ``A`` and return a ``CholeskyDensePivoted`` object. ``LU`` may be 'L' for using the lower part or 'U' for the upper part. The default is to use 'U'. ``factors(cholpivot(A))`` returns the triangular matrix containing the factorization. The following functions are available for ``CholeskyDensePivoted`` objects: ``size``, ``factors``, ``\``, ``inv``, ``det``. A ``LAPACK.RankDeficientException`` error is thrown in case the matrix is rank deficient.
 
 .. function:: cholpivot!(A, [LU])
 
@@ -1744,9 +1748,9 @@ Linear algebra functions in Julia are largely implemented by calling functions f
 
    Compute QR factorization with pivoting
 
-.. function:: factors(D)
+.. function:: factors(F)
 
-   Return the factors of a decomposition D. For an LU decomposition, factors(LU) -> L, U, p
+   Return the factors from a ``Factorization`` object.
 
 .. function:: eig(A) -> D, V
 
diff --git a/test/linalg.jl b/test/linalg.jl
index 5918c1d1d38e5..a0917e5bbaae4 100644
--- a/test/linalg.jl
+++ b/test/linalg.jl
@@ -32,10 +32,8 @@ for elty in (Float32, Float64, Complex64, Complex128)
         @assert_approx_eq inv(bc2) * apd eye(elty, n)
         @assert_approx_eq apd * (bc2\b) b
 
-        lua   = lud(a)                  # LU decomposition
-        l,u,p = lu(a)
-        L,U,P = factors(lua)
-        @test l == L && u == U && p == P
+        lua   = lu(a)                  # LU decomposition
+        l,u,p = factors(lua)
         @assert_approx_eq l*u a[p,:]
         @assert_approx_eq l[invperm(p),:]*u a
         @assert_approx_eq a * inv(lua) eye(elty, n)
@@ -290,7 +288,7 @@ for elty in (Float32, Float64, Complex64, Complex128)
         @assert_approx_eq solve(T,v) invFv
         B = convert(Matrix{elty}, B)
         @assert_approx_eq solve(T, B) F\B
-        Tlu = lud(T)
+        Tlu = lu(T)
         x = Tlu\v
         @assert_approx_eq x invFv
         @assert_approx_eq det(T) det(F)

From 2fe3c04f0442fcad1b636955e77b28a2b062490b Mon Sep 17 00:00:00 2001
From: "Viral B. Shah" <viral@mayin.org>
Date: Wed, 6 Feb 2013 19:59:02 -0500
Subject: [PATCH 06/11] Rename QRPDense to QRDensePivoted. Use qr and qrpivot
 instead of qrd and qrpd. Update docs for qr.

---
 base/exports.jl      | 10 ++++------
 base/linalg_dense.jl | 42 +++++++++++++++++++++---------------------
 doc/stdlib/base.rst  | 14 +++++++++++---
 test/linalg.jl       | 12 ++++--------
 4 files changed, 40 insertions(+), 38 deletions(-)

diff --git a/base/exports.jl b/base/exports.jl
index e194b1ef28fad..4e6b88e95830f 100644
--- a/base/exports.jl
+++ b/base/exports.jl
@@ -109,7 +109,7 @@ export
     LUTridiagonal,
     LDLTTridiagonal,
     QRDense,
-    QRPDense,
+    QRDensePivoted,
     InsertionSort,
     QuickSort,
     MergeSort,
@@ -559,11 +559,9 @@ export
     null,
     pinv,
     qr,
-    qrd!,
-    qrd,
-    qrp,
-    qrpd!,
-    qrpd,
+    qr!,
+    qrpivot,
+    qrpivot!,
     randsym,
     rank,
     rref,
diff --git a/base/linalg_dense.jl b/base/linalg_dense.jl
index e4f0c123307f3..b303f7026ddee 100644
--- a/base/linalg_dense.jl
+++ b/base/linalg_dense.jl
@@ -573,17 +573,16 @@ end
 size(A::QRDense) = size(A.hh)
 size(A::QRDense,n) = size(A.hh,n)
 
-qrd!{T<:BlasFloat}(A::StridedMatrix{T}) = QRDense{T}(LAPACK.geqrf!(A)...)
-qrd{T<:BlasFloat}(A::StridedMatrix{T}) = qrd!(copy(A))
-qrd{T<:Real}(A::StridedMatrix{T}) = qrd(float64(A))
+qr!{T<:BlasFloat}(A::StridedMatrix{T}) = QRDense{T}(LAPACK.geqrf!(A)...)
+qr{T<:BlasFloat}(A::StridedMatrix{T}) = qr!(copy(A))
+qr{T<:Real}(A::StridedMatrix{T}) = qr(float64(A))
 
-function factors{T<:BlasFloat}(qrd::QRDense{T})
-    aa  = copy(qrd.hh)
+function factors{T<:BlasFloat}(qr::QRDense{T})
+    aa  = copy(qr.hh)
     R   = triu(aa[1:min(size(aa)),:])   # must be *before* call to orgqr!
-    LAPACK.orgqr!(aa, qrd.tau, size(aa,2)), R
+    LAPACK.orgqr!(aa, qr.tau, size(aa,2)), R
 end
 
-qr{T<:Number}(x::StridedMatrix{T}) = factors(qrd(x))
 qr(x::Number) = (one(x), x)
 
 ## Multiplication by Q from the QR decomposition
@@ -602,41 +601,42 @@ function (\){T<:BlasFloat}(A::QRDense{T}, B::StridedVecOrMat{T})
     return ans
 end
 
-type QRPDense{T} <: Factorization{T}
+## Pivoted QR decomposition
+
+type QRDensePivoted{T} <: Factorization{T}
     hh::Matrix{T}
     tau::Vector{T}
     jpvt::Vector{BlasInt}
-    function QRPDense(hh::Matrix{T}, tau::Vector{T}, jpvt::Vector{BlasInt})
+    function QRDensePivoted(hh::Matrix{T}, tau::Vector{T}, jpvt::Vector{BlasInt})
         m, n = size(hh)
         if length(tau) != min(m,n) || length(jpvt) != n
-            error("QRPDense: mismatched dimensions")
+            error("QRDensePivoted: mismatched dimensions")
         end
         new(hh,tau,jpvt)
     end
 end
-size(x::QRPDense)   = size(x.hh)
-size(x::QRPDense,d) = size(x.hh,d)
+size(x::QRDensePivoted)   = size(x.hh)
+size(x::QRDensePivoted,d) = size(x.hh,d)
 ## Multiplication by Q from the QR decomposition
-(*){T<:BlasFloat}(A::QRPDense{T}, B::StridedVecOrMat{T}) =
+(*){T<:BlasFloat}(A::QRDensePivoted{T}, B::StridedVecOrMat{T}) =
     LAPACK.ormqr!('L', 'N', A.hh, size(A,2), A.tau, copy(B))
 ## Multiplication by Q' from the QR decomposition
-Ac_mul_B{T<:BlasFloat}(A::QRPDense{T}, B::StridedVecOrMat{T}) =
+Ac_mul_B{T<:BlasFloat}(A::QRDensePivoted{T}, B::StridedVecOrMat{T}) =
     LAPACK.ormqr!('L', iscomplex(A.tau)?'C':'T', A.hh, size(A,2), A.tau, copy(B))
 
-qrpd!{T<:BlasFloat}(A::StridedMatrix{T}) = QRPDense{T}(LAPACK.geqp3!(A)...)
-qrpd{T<:BlasFloat}(A::StridedMatrix{T}) = qrpd!(copy(A))
-qrpd{T<:Real}(x::StridedMatrix{T}) = qrpd(float64(x))
+qrpivot!{T<:BlasFloat}(A::StridedMatrix{T}) = QRDensePivoted{T}(LAPACK.geqp3!(A)...)
+qrpivot{T<:BlasFloat}(A::StridedMatrix{T}) = qrpivot!(copy(A))
+qrpivot{T<:Real}(x::StridedMatrix{T}) = qrpivot(float64(x))
 
-function factors{T<:BlasFloat}(x::QRPDense{T})
+function factors{T<:BlasFloat}(x::QRDensePivoted{T})
     aa = copy(x.hh)
     R  = triu(aa[1:min(size(aa)),:])
     LAPACK.orgqr!(aa, x.tau, size(aa,2)), R, x.jpvt
 end
 
-qrp{T<:BlasFloat}(x::StridedMatrix{T}) = factors(qrpd(x))
-qrp{T<:Real}(x::StridedMatrix{T}) = qrp(float64(x))
+qrpivot{T<:Real}(x::StridedMatrix{T}) = qrpivot(float64(x))
 
-function (\){T<:BlasFloat}(A::QRPDense{T}, B::StridedVecOrMat{T})
+function (\){T<:BlasFloat}(A::QRDensePivoted{T}, B::StridedVecOrMat{T})
     n = length(A.tau)
     x, info = LAPACK.trtrs!('U','N','N',A.hh[1:n,:],(A'*B)[1:n,:])
     if info > 0; throw(LAPACK.SingularException(info)); end
diff --git a/doc/stdlib/base.rst b/doc/stdlib/base.rst
index 27bf5ecd59d9d..445062096ba63 100644
--- a/doc/stdlib/base.rst
+++ b/doc/stdlib/base.rst
@@ -1742,11 +1742,19 @@ Linear algebra functions in Julia are largely implemented by calling functions f
 
 .. function:: qr(A)
 
-   Compute QR factorization
+   Compute the QR factorization of ``A`` and return a ``QRDense`` object. ``factors(qr(A))`` returns ``Q`` and ``R``. The following functions are available for ``QRDense`` objects: ``size``, ``factors``, ``*``, ``Ac_mul_B``, ``\``. 
 
-.. function:: qrp(A)
+.. function:: qr!(A)
 
-   Compute QR factorization with pivoting
+   ``qr!`` is the same as ``qr``, but overwrites the input matrix A with the factorization.
+
+.. function:: qrpivot(A)
+
+   Compute the QR factorization of ``A`` and return a ``QRDensePivoted`` object. ``factors(qr(A))`` returns ``Q`` and ``R``. The following functions are available for ``QRDensePivoted`` objects: ``size``, ``factors``, ``*``, ``Ac_mul_B``, ``\``. 
+
+.. function:: qrpivot!(A)
+
+   ``qrpivot!`` is the same as ``qrpivot``, but overwrites the input matrix A with the factorization.
 
 .. function:: factors(F)
 
diff --git a/test/linalg.jl b/test/linalg.jl
index a0917e5bbaae4..caccfbb733158 100644
--- a/test/linalg.jl
+++ b/test/linalg.jl
@@ -39,23 +39,19 @@ for elty in (Float32, Float64, Complex64, Complex128)
         @assert_approx_eq a * inv(lua) eye(elty, n)
         @assert_approx_eq a*(lua\b) b
 
-        qra   = qrd(a)                  # QR decomposition
+        qra   = qr(a)                  # QR decomposition
         q,r   = factors(qra)
         @assert_approx_eq q'*q eye(elty, n)
         @assert_approx_eq q*q' eye(elty, n)
-        Q,R   = qr(a)
-        @test q == Q && r == R
         @assert_approx_eq q*r a
-        @assert_approx_eq qra*b Q*b
-        @assert_approx_eq qra'*b Q'*b
+        @assert_approx_eq qra*b q*b
+        @assert_approx_eq qra'*b q'*b
         @assert_approx_eq a*(qra\b) b
 
-        qrpa  = qrpd(a)                 # pivoted QR decomposition
+        qrpa  = qrpivot(a)                 # pivoted QR decomposition
         q,r,p = factors(qrpa)
         @assert_approx_eq q'*q eye(elty, n)
         @assert_approx_eq q*q' eye(elty, n)
-        Q,R,P = qrp(a)
-        @test q == Q && r == R && p == P
         @assert_approx_eq q*r a[:,p]
         @assert_approx_eq q*r[:,invperm(p)] a
         @assert_approx_eq a*(qrpa\b) b

From b259333030fdc659eda47715eb0bbdcadd8bea42 Mon Sep 17 00:00:00 2001
From: "Viral B. Shah" <viral@mayin.org>
Date: Wed, 6 Feb 2013 20:06:06 -0500
Subject: [PATCH 07/11] Deprecate old dense matrix factorization names.

---
 base/deprecated.jl | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/base/deprecated.jl b/base/deprecated.jl
index 8307f702ff9bf..51fc0b2921a86 100644
--- a/base/deprecated.jl
+++ b/base/deprecated.jl
@@ -59,6 +59,11 @@ end
 @deprecate  chars         collect
 @deprecate  elements      collect
 @deprecate  strcat        string
+@deprecate  chold         chol
+@deprecate  cholpd        cholpivot
+@deprecate  lud           lu
+@deprecate  qrd           qr
+@deprecate  qrpd          qrpivot
 
 export randi
 function randi(n,x...)

From bf427cdcfd9e031d978e19d39e8efe7e14f815dc Mon Sep 17 00:00:00 2001
From: "Viral B. Shah" <viral@mayin.org>
Date: Wed, 6 Feb 2013 20:14:52 -0500
Subject: [PATCH 08/11] Replace lud with lu in suitesparse.jl

---
 extras/suitesparse.jl | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/extras/suitesparse.jl b/extras/suitesparse.jl
index 81dd1b2cb62e6..9106889c81060 100644
--- a/extras/suitesparse.jl
+++ b/extras/suitesparse.jl
@@ -1,8 +1,8 @@
 module SuiteSparse
 
 import Base.SparseMatrixCSC, Base.size, Base.nnz, Base.eltype, Base.show
-import Base.triu, Base.norm, Base.solve, Base.(\), Base.lu
-import Base.Ac_ldiv_B, Base.At_ldiv_B, Base.lud, Base.lud!
+import Base.triu, Base.norm, Base.solve, Base.(\), Base.lu, Base.lu!
+import Base.Ac_ldiv_B, Base.At_ldiv_B
 
 import Base.BlasInt
 import Base.blas_int
@@ -14,7 +14,7 @@ export                                  # types
     CholmodFactor,
     CholmodDense,
     CholmodSparseOut,
-    UmfpackLU,                          # call these lud and lud! instead?
+    UmfpackLU,                          # call these lu and lu! instead?
     UmfpackLU!,
                                         # methods
     decrement,
@@ -89,7 +89,7 @@ type UmfpackLU{Tv<:UMFVTypes,Ti<:CHMITypes} <: Factorization{Tv}
     nzval::Vector{Tv}
 end
 
-function lud{Tv<:UMFVTypes,Ti<:CHMITypes}(S::SparseMatrixCSC{Tv,Ti})
+function lu{Tv<:UMFVTypes,Ti<:CHMITypes}(S::SparseMatrixCSC{Tv,Ti})
     zerobased = S.colptr[1] == 0
     lu = UmfpackLU(C_NULL, C_NULL, S.m, S.n,
                    zerobased ? copy(S.colptr) : decrement(S.colptr),
@@ -98,7 +98,7 @@ function lud{Tv<:UMFVTypes,Ti<:CHMITypes}(S::SparseMatrixCSC{Tv,Ti})
     umfpack_numeric!(lu)
 end
 
-function lud!{Tv<:UMFVTypes,Ti<:CHMITypes}(S::SparseMatrixCSC{Tv,Ti})
+function lu!{Tv<:UMFVTypes,Ti<:CHMITypes}(S::SparseMatrixCSC{Tv,Ti})
     zerobased = S.colptr[1] == 0
     UmfpackLU(C_NULL, C_NULL, S.m, S.n,
               zerobased ? S.colptr : decrement!(S.colptr),
@@ -119,13 +119,13 @@ end
 
 # Solve directly with matrix
 
-(\)(S::SparseMatrixCSC, b::Vector) = lud(S) \ b
-At_ldiv_B{T<:UMFVTypes}(S::SparseMatrixCSC{T}, b::Vector{T}) = umfpack_solve(lud(S), b, UMFPACK_Aat)
+(\)(S::SparseMatrixCSC, b::Vector) = lu(S) \ b
+At_ldiv_B{T<:UMFVTypes}(S::SparseMatrixCSC{T}, b::Vector{T}) = umfpack_solve(lu(S), b, UMFPACK_Aat)
 function At_ldiv_B{Ts<:UMFVTypes,Tb<:Number}(S::SparseMatrixCSC{Ts}, b::Vector{Tb})
     ## should be more careful here in case Ts<:Real and Tb<:Complex
     At_ldiv_B(S, convert(Vector{Ts}, b))
 end
-Ac_ldiv_B{T<:UMFVTypes}(S::SparseMatrixCSC{T}, b::Vector{T}) = umfpack_solve(lud(S), b, UMFPACK_At)
+Ac_ldiv_B{T<:UMFVTypes}(S::SparseMatrixCSC{T}, b::Vector{T}) = umfpack_solve(lu(S), b, UMFPACK_At)
 function Ac_ldiv_B{Ts<:UMFVTypes,Tb<:Number}(S::SparseMatrixCSC{Ts}, b::Vector{Tb})
     ## should be more careful here in case Ts<:Real and Tb<:Complex
     Ac_ldiv_B(S, convert(Vector{Ts}, b))

From 5ddf16b76f4fb5d7d120d4736db1bdf4865d3f5a Mon Sep 17 00:00:00 2001
From: Douglas Bates <dmbates@gmail.com>
Date: Thu, 7 Feb 2013 16:10:54 -0600
Subject: [PATCH 09/11] Got the basic CHOLMOD functions associated with a dense
 matrix working by defining a Julia type

Well, working on a 64-bit system and these are only the calls to
cholmod_check_dense and cholmod_print_dense, as in

julia> aa = randn(10,2)
10x2 Float64 Array:
 -1.1131     -0.276926
 -1.18138    -2.17424
  0.615017   -0.341257
  0.0714101  -0.0933698
 -0.493773   -0.520607
  0.989568   -0.0106482
  1.33313    -0.204813
 -0.199893   -0.983869
  0.934582   -0.0652311
  0.792507   -0.369899

julia> cda = CholmodDense(aa)
CholmodDense{Float64}(10,2,20,10,Ptr{Float64} @0x0000000003186070,Ptr{Void} @0x0000000000000000,1,0,10x2 Float64 Array:
 -1.1131     -0.276926
 -1.18138    -2.17424
  0.615017   -0.341257
  0.0714101  -0.0933698
 -0.493773   -0.520607
  0.989568   -0.0106482
  1.33313    -0.204813
 -0.199893   -0.983869
  0.934582   -0.0652311
  0.792507   -0.369899 ,false)

julia> SuiteSparse.chm_chk_dn(cda)
1

julia> SuiteSparse.chm_prt_dn(cda)

CHOLMOD dense:   :  10-by-2,
  leading dimension 10, nzmax 20, real, double
  col 0:
         0: -1.1131
         1: -1.1814
         2: 0.61502
         3: 0.07141
         4: -0.49377
         5: 0.98957
         6: 1.3331
         7: -0.19989
    ...
  col 1:
         0: -0.27693
         1: -2.1742
         2: -0.34126
         3: -0.09337
    ...
         6: -0.20481
         7: -0.98387
         8: -0.065231
         9: -0.3699
  OK

4-element Uint8 Array:
 0x03
 0x00
 0x00
 0x00

Not exactly earth-shattering but at least it is a proof-of-concept for
using a Julia type to pass a pointer to a struct to one of the CHOLMOD
functions.
---
 extras/suitesparse.jl | 143 ++++++++++++++++++++++++++++++------------
 1 file changed, 103 insertions(+), 40 deletions(-)

diff --git a/extras/suitesparse.jl b/extras/suitesparse.jl
index 9106889c81060..35ebf67077ef6 100644
--- a/extras/suitesparse.jl
+++ b/extras/suitesparse.jl
@@ -1,19 +1,19 @@
 module SuiteSparse
 
 import Base.SparseMatrixCSC, Base.size, Base.nnz, Base.eltype, Base.show
-import Base.triu, Base.norm, Base.solve, Base.(\), Base.lu, Base.lu!
+import Base.triu, Base.norm, Base.solve, Base.lu, Base.lu!, Base.(\)
 import Base.Ac_ldiv_B, Base.At_ldiv_B
 
 import Base.BlasInt
 import Base.blas_int
 
 export                                  # types
-    CholmodPtr,
-    CholmodCommon,
-    CholmodSparse,
-    CholmodFactor,
+#    CholmodPtr,
+#    CholmodCommon,
+#    CholmodSparse,
+#    CholmodFactor,
     CholmodDense,
-    CholmodSparseOut,
+#    CholmodSparseOut,
     UmfpackLU,                          # call these lu and lu! instead?
     UmfpackLU!,
                                         # methods
@@ -338,17 +338,80 @@ end
 ### there must be an easier way but at least this works.
 function ChmCommon(aa::Array{Uint8,1})
     eval(Expr(:call,
-              unshift!({reinterpret(ChmCommon.types[i],
+              unshift!({
+                        reinterpret(ChmCommon.types[i],
                                     aa[(1:sizeof(ChmCommon.types[i])) +
                                        chm_com_offsets[i]])[1]
-                        for i in 1:length(ChmCommon.types)}, :ChmCommon),
-              Any))
+                        for i in 1:length(ChmCommon.types)
+                            },
+                            :ChmCommon),
+                            Any)
+                            )
 end
 
-chm_itype{Tv,Ti<:CHMITypes}(S::SparseMatrixCSC{Tv,Ti}) = Ti == Int32 ? CHOLMOD_INT : CHOLMOD_LONG
-chm_xtype{Tv<:CHMVTypes}(S::SparseMatrixCSC{Tv}) = Tv <: Complex ? CHOLMOD_COMPLEX : CHOLMOD_REAL
-chm_dtype{Tv<:CHMVTypes}(S::SparseMatrixCSC{Tv}) = T <: Union(Float32, Complex64) ? CHOLMOD_SINGLE : CHOLMOD_DOUBLE
+function chm_itype{Tv,Ti<:CHMITypes}(S::SparseMatrixCSC{Tv,Ti})
+    int32(Ti == Int32 ? CHOLMOD_INT : CHOLMOD_LONG)
+end
+function chm_xtype{Tv<:CHMVTypes}(S::SparseMatrixCSC{Tv})
+    int32(Tv <: Complex ? CHOLMOD_COMPLEX : CHOLMOD_REAL)
+end
+function chm_dtype{Tv<:CHMVTypes}(S::SparseMatrixCSC{Tv})
+    int32(T <: Union(Float32, Complex64) ? CHOLMOD_SINGLE : CHOLMOD_DOUBLE)
+end
+
+function set_chm_prt_lev(cm::Array{Uint8}, lev::Integer)
+    cm[(1:4) + chm_com_offsets[13]] = reinterpret(Uint8, [int32(lev)])
+end
+
+type CholmodDense{T<:CHMVTypes}
+    m::Int
+    n::Int
+    nnz::Int
+    lda::Int
+    xpt::Ptr{T}
+    zpt::Ptr{Void}
+    xtype::Int32
+    dtype::Int32
+    aa::Array{T}
+    chm::Bool                       # was storage allocated by Cholmod?
+end
+
+function CholmodDense{T<:CHMVTypes}(aa::VecOrMat{T})
+    m = size(aa,1); n = size(aa,2)
+    typs = CholmodDense{T}.types[1:8]   # doesn't actually depend on T
+    sz = [sizeof(t) for t in typs]
+    inds = zeros(Int,9)
+    inds[2:9] = cumsum(sz)
+    CholmodDense{T}(m, n, m*n, stride(aa,2), convert(Ptr{T}, aa), C_NULL,
+                    int32(T<:Complex ? CHOLMOD_COMPLEX : CHOLMOD_REAL),
+                    int32(T<:Union(Float32,Complex64) ? CHOLMOD_SINGLE : CHOLMOD_DOUBLE),
+                    aa, false)
+end
+
+function chm_chk_dn{T<:CHMVTypes}(cd::CholmodDense{T})
+    ccall((:cholmod_check_dense, :libcholmod), Int,
+          (Ptr{CholmodDense{T}}, Ptr{Uint8}), &cd, chm_com)
+end
 
+function chm_chk_dn{T<:CHMVTypes}(cd::CholmodDense{T})
+    ccall((:cholmod_check_dense, :libcholmod), Int,
+          (Ptr{CholmodDense{T}}, Ptr{Uint8}), &cd, chm_com)
+end
+
+function chm_prt_dn{T<:CHMVTypes}(cd::CholmodDense{T}, lev::Integer, nm::ASCIIString)
+    inds = (1:4) + chm_com_offsets[13]
+    orig = chm_com[inds]
+    chm_com[inds] = reinterpret(Uint8, [int32(lev)])
+    status = ccall((:cholmod_print_dense, :libcholmod), Int32,
+                   (Ptr{CholmodDense{T}}, Ptr{Uint8}, Ptr{Uint8}),
+                   &cd, nm, chm_com)
+    chm_com[inds] = orig
+end
+
+chm_prt_dn{T<:CHMVTypes}(cd::CholmodDense{T}, lev::Integer) = chm_prt_dn(cd, lev, "")    
+chm_prt_dn{T<:CHMVTypes}(cd::CholmodDense{T}) = chm_prt_dn(cd, int32(4), "")
+    
+if false    
 # Wrapper for memory allocated by CHOLMOD. Carry along the value and index types.
 ## FIXME: CholmodPtr and UmfpackPtr should be amalgamated
 type CholmodPtr{Tv<:CHMVTypes,Ti<:CHMITypes}
@@ -358,38 +421,38 @@ end
 eltype{Tv,Ti}(P::CholmodPtr{Tv,Ti}) = Tv
 indtype{Tv,Ti}(P::CholmodPtr{Tv,Ti}) = Ti
 
-function cholmod_common_finalizer(x::Vector{Ptr{Void}})
-    st = ccall((:cholmod_finish, :libcholmod), BlasInt, (Ptr{Void},), x[1])
-    if st != CHOLMOD_TRUE error("Error calling cholmod_finish") end
-    c_free(x[1])
-end
+## function cholmod_common_finalizer(x::Vector{Ptr{Void}})
+##     st = ccall((:cholmod_finish, :libcholmod), BlasInt, (Ptr{Void},), x[1])
+##     if st != CHOLMOD_TRUE error("Error calling cholmod_finish") end
+##     c_free(x[1])
+## end
 
-type CholmodCommon
-    pt::Vector{Ptr{Void}}
-    function CholmodCommon()
-        pt = Array(Ptr{Void}, 1)
-        ccall((:jl_cholmod_common, :libsuitesparse_wrapper), Void,
-              (Ptr{Void},), pt)
-        st = ccall((:cholmod_start, :libcholmod), BlasInt, (Ptr{Void}, ), pt[1])
-        if st != CHOLMOD_TRUE error("Error calling cholmod_start") end
-        finalizer(pt, cholmod_common_finalizer)
-        new(pt)
-    end
-end
+## type CholmodCommon
+##     pt::Vector{Ptr{Void}}
+##     function CholmodCommon()
+##         pt = Array(Ptr{Void}, 1)
+##         ccall((:jl_cholmod_common, :libsuitesparse_wrapper), Void,
+##               (Ptr{Void},), pt)
+##         st = ccall((:cholmod_start, :libcholmod), BlasInt, (Ptr{Void}, ), pt[1])
+##         if st != CHOLMOD_TRUE error("Error calling cholmod_start") end
+##         finalizer(pt, cholmod_common_finalizer)
+##         new(pt)
+##     end
+## end
 
-function show(io::IO, cm::CholmodCommon)
-    st = ccall((:cholmod_print_common, :libcholmod), BlasInt,
-               (Ptr{Uint8},Ptr{Void}), "", cm.pt[1])
-    if st != CHOLMOD_TRUE error("Error calling cholmod_print_common") end
-end
+## function show(io::IO, cm::CholmodCommon)
+##     st = ccall((:cholmod_print_common, :libcholmod), BlasInt,
+##                (Ptr{Uint8},Ptr{Void}), "", cm.pt[1])
+##     if st != CHOLMOD_TRUE error("Error calling cholmod_print_common") end
+## end
 
 type CholmodSparse{Tv<:CHMVTypes,Ti<:CHMITypes}
     pt::CholmodPtr{Tv,Ti}
     ## cp contains a copy of the original matrix but with 0-based indices
     cp::SparseMatrixCSC{Tv,Ti}
     stype::Int
-    cm::CholmodCommon
-    function CholmodSparse(S::SparseMatrixCSC{Tv,Ti}, stype::BlasInt, cm::CholmodCommon)
+#    cm::CholmodCommon
+    function CholmodSparse(S::SparseMatrixCSC{Tv,Ti}, stype::BlasInt)#, cm::CholmodCommon)
         pt = CholmodPtr{Tv,Ti}(Array(Ptr{Void}, 1))
         cp = convert_to_0_based_indexing(S)
         
@@ -406,12 +469,12 @@ end
 
 CholmodSparse{Tv,Ti}(S::SparseMatrixCSC{Tv,Ti}, stype::Int) = CholmodSparse{Tv,Ti}(S, stype, CholmodCommon())
 
-function CholmodSparse{Tv,Ti}(S::SparseMatrixCSC{Tv,Ti}, cm::CholmodCommon)
+function CholmodSparse{Tv,Ti}(S::SparseMatrixCSC{Tv,Ti})#, cm::CholmodCommon)
     stype = S.m == S.n && ishermitian(S)
     CholmodSparse{Tv,Ti}(stype ? triu(S) : S, blas_int(stype), cm)
 end
 
-CholmodSparse(S::SparseMatrixCSC) = CholmodSparse(S, CholmodCommon())
+CholmodSparse(S::SparseMatrixCSC) = CholmodSparse(S)#, CholmodCommon())
 
 function show(io::IO, cs::CholmodSparse)
     ccall(_chm_print_sp,
@@ -490,7 +553,7 @@ type CholmodDense{T<:CHMVTypes}
     m::Int
     n::Int
     aa::VecOrMat{T}                     # original array
-    cm::CholmodCommon
+#    cm::CholmodCommon
 end
 
 function CholmodDense{T<:CHMVTypes}(b::VecOrMat{T}, cm::CholmodCommon)
@@ -731,5 +794,5 @@ function cholmod_solve(cs_factor::Ptr{Void}, cd_rhs::Array{Ptr{Void},1}, cm::Arr
     ccall(:cholmod_solve, Ptr{Void}, (BlasInt, Ptr{Void}, Ptr{Void}, Ptr{Void}),
           CHOLMOD_A, cs_factor, cd_rhs[1], cm[1])
 end
-
+end # code omission
 end #module

From 2ce6abcb4a7aefd6658f4b65cbb03057a3e02fbe Mon Sep 17 00:00:00 2001
From: Douglas Bates <dmbates@gmail.com>
Date: Thu, 7 Feb 2013 17:07:07 -0600
Subject: [PATCH 10/11] Defined a CholmodSparse type

but my external constructor from a SparseMatrix{Tv,Ti} is not working
as hoped.  I'm pushing this commit in case someone can try
```julia
require("suitesparse")
using SuiteSparse
A = sprand(8,4,0.25);
ca = CholmodSparse(A)
```
and tell me why I get a no method error.
---
 extras/suitesparse.jl | 51 ++++++++++++++++++++++++++++++++++---------
 1 file changed, 41 insertions(+), 10 deletions(-)

diff --git a/extras/suitesparse.jl b/extras/suitesparse.jl
index 35ebf67077ef6..6a6e397a35f1d 100644
--- a/extras/suitesparse.jl
+++ b/extras/suitesparse.jl
@@ -10,7 +10,7 @@ import Base.blas_int
 export                                  # types
 #    CholmodPtr,
 #    CholmodCommon,
-#    CholmodSparse,
+    CholmodSparse,
 #    CholmodFactor,
     CholmodDense,
 #    CholmodSparseOut,
@@ -352,11 +352,11 @@ end
 function chm_itype{Tv,Ti<:CHMITypes}(S::SparseMatrixCSC{Tv,Ti})
     int32(Ti == Int32 ? CHOLMOD_INT : CHOLMOD_LONG)
 end
-function chm_xtype{Tv<:CHMVTypes}(S::SparseMatrixCSC{Tv})
-    int32(Tv <: Complex ? CHOLMOD_COMPLEX : CHOLMOD_REAL)
+function chm_xtype{T<:CHMVTypes}(S::SparseMatrixCSC{T})
+    int32(T<:Complex ? CHOLMOD_COMPLEX : CHOLMOD_REAL)
 end
-function chm_dtype{Tv<:CHMVTypes}(S::SparseMatrixCSC{Tv})
-    int32(T <: Union(Float32, Complex64) ? CHOLMOD_SINGLE : CHOLMOD_DOUBLE)
+function chm_dtype{T<:CHMVTypes}(S::SparseMatrixCSC{T})
+    int32(T<:Union(Float32, Complex64) ? CHOLMOD_SINGLE : CHOLMOD_DOUBLE)
 end
 
 function set_chm_prt_lev(cm::Array{Uint8}, lev::Integer)
@@ -366,7 +366,7 @@ end
 type CholmodDense{T<:CHMVTypes}
     m::Int
     n::Int
-    nnz::Int
+    nzmax::Int
     lda::Int
     xpt::Ptr{T}
     zpt::Ptr{Void}
@@ -378,10 +378,6 @@ end
 
 function CholmodDense{T<:CHMVTypes}(aa::VecOrMat{T})
     m = size(aa,1); n = size(aa,2)
-    typs = CholmodDense{T}.types[1:8]   # doesn't actually depend on T
-    sz = [sizeof(t) for t in typs]
-    inds = zeros(Int,9)
-    inds[2:9] = cumsum(sz)
     CholmodDense{T}(m, n, m*n, stride(aa,2), convert(Ptr{T}, aa), C_NULL,
                     int32(T<:Complex ? CHOLMOD_COMPLEX : CHOLMOD_REAL),
                     int32(T<:Union(Float32,Complex64) ? CHOLMOD_SINGLE : CHOLMOD_DOUBLE),
@@ -410,6 +406,41 @@ end
 
 chm_prt_dn{T<:CHMVTypes}(cd::CholmodDense{T}, lev::Integer) = chm_prt_dn(cd, lev, "")    
 chm_prt_dn{T<:CHMVTypes}(cd::CholmodDense{T}) = chm_prt_dn(cd, int32(4), "")
+
+type CholmodSparse{Tv<:CHMVTypes,Ti<:CHMITypes}
+    m::Int
+    n::Int
+    nzmax::Int
+    ppt::Ptr{Ti}
+    ipt::Ptr{Ti}
+    nzpt::Ptr{Void}
+    xpt::Ptr{Tv}
+    zpt::Ptr{Void}
+    stype::Int32
+    itype::Int32
+    xtype::Int32
+    dtype::Int32
+    sorted::Int32
+    packed::Int32
+    colptr0::Vector{Ti} # 0-based column pointers
+    rowval0::Vector{Ti} # 0-based row indices
+    nzval::Vector{Tv}   # a copy of the non-zero values
+    chm_free::Bool      # was storage allocated by Cholmod?
+    structpt::Ptr{Void} # pointer to the C struct (when created from Cholmod results)
+end
+
+function CholmodSparse{Tv<:CHMVTypes,Ti<:CHMITypes}(A::SparseMatrixCSC{Tv,Ti})
+    colptr0 = decrement(A.colptr)
+    rowval0 = decrement(A.rowval)
+    nzval = copy(A.nzval)
+    CholmodSparse{Tv,Ti}(size(A,1),size(A,2),int(colptr0[end]),
+                         convert(Ptr{Ti}, colptr0), convert(Ptr{Ti}, rowval0),
+                         convert(Ptr{Tv}, nzval), C_NULL, int32(0),
+                         chm_itype(A), chm_xtype(A), chm_dtype(A),
+                         CHOLMOD_TRUE, CHOLMOD_TRUE,
+                         colptr0, rowval0, nzval, false, C_NULL)
+end
+
     
 if false    
 # Wrapper for memory allocated by CHOLMOD. Carry along the value and index types.

From 614530e2d6023100974e9827d1e5275009ed004b Mon Sep 17 00:00:00 2001
From: Douglas Bates <dmbates@gmail.com>
Date: Fri, 8 Feb 2013 17:43:20 -0600
Subject: [PATCH 11/11] Incorporate suggestion from Jameson on keeping pointer
 contents alive

Defined types c_CholmodDense which mirrors the cholmod_dense C struct and a wrapper
type CholmodDense that retains the Julia matrix when constructed from one and retains
the pointer itself when constructed from *cholmod_dense returned from C.

Still need to add a finalizer for the CholmodDense type and create
similar types for CholmodSparse, CholmodFactor and CholmodTriplet.
---
 extras/suitesparse.jl | 529 +++++++-----------------------------------
 1 file changed, 90 insertions(+), 439 deletions(-)

diff --git a/extras/suitesparse.jl b/extras/suitesparse.jl
index 6a6e397a35f1d..ad06f2cb6c6ec 100644
--- a/extras/suitesparse.jl
+++ b/extras/suitesparse.jl
@@ -1,33 +1,33 @@
 module SuiteSparse
 
 import Base.SparseMatrixCSC, Base.size, Base.nnz, Base.eltype, Base.show
-import Base.triu, Base.norm, Base.solve, Base.lu, Base.lu!, Base.(\)
-import Base.Ac_ldiv_B, Base.At_ldiv_B
+import Base.triu, Base.norm, Base.solve, Base.lu, Base.lu!
+#, Base.(\)
+import Base.Ac_ldiv_B, Base.At_ldiv_B, Base.convert
 
 import Base.BlasInt
 import Base.blas_int
 
 export                                  # types
-#    CholmodPtr,
-#    CholmodCommon,
-    CholmodSparse,
+CholmodSparse,
 #    CholmodFactor,
-    CholmodDense,
+CholmodDense,
 #    CholmodSparseOut,
-    UmfpackLU,                          # call these lu and lu! instead?
-    UmfpackLU!,
-                                        # methods
-    decrement,
-    decrement!,
-    increment,
-    increment!,
-    indtype,
-    show_umf_ctrl,
-    show_umf_info
+UmfpackLU,                          # call these lu and lu! instead?
+UmfpackLU!,
+# methods
+decrement,
+decrement!,
+increment,
+increment!,
+indtype,
+show_umf_ctrl,
+show_umf_info
 
 include("suitesparse_h.jl")
 
 type MatrixIllConditionedException <: Exception end
+type CholmodException <: Exception end
 
 function decrement!{T<:Integer}(A::AbstractArray{T})
     for i in 1:length(A) A[i] -= one(T) end
@@ -108,18 +108,18 @@ end
 
 function show(io::IO, f::UmfpackLU)
     @printf(io, "UMFPACK LU Factorization of a %d-by-%d sparse matrix\n",
-           f.m, f.n)
+            f.m, f.n)
     if f.numeric != C_NULL println(f.numeric) end
 end
 
-# Solve with Factorization
+### Solve with Factorization
 
-(\){T<:UMFVTypes}(fact::UmfpackLU{T}, b::Vector{T}) = umfpack_solve(fact, b)
-(\){Ts<:UMFVTypes,Tb<:Number}(fact::UmfpackLU{Ts}, b::Vector{Tb}) = fact\convert(Vector{Ts},b)
+#(\){T<:UMFVTypes}(fact::UmfpackLU{T}, b::Vector{T}) = umfpack_solve(fact, b)
+#(\){Ts<:UMFVTypes,Tb<:Number}(fact::UmfpackLU{Ts}, b::Vector{Tb}) = fact\convert(Vector{Ts},b)
+  
+### Solve directly with matrix
 
-# Solve directly with matrix
-
-(\)(S::SparseMatrixCSC, b::Vector) = lu(S) \ b
+#(\)(S::SparseMatrixCSC, b::Vector) = lu(S) \ b
 At_ldiv_B{T<:UMFVTypes}(S::SparseMatrixCSC{T}, b::Vector{T}) = umfpack_solve(lu(S), b, UMFPACK_Aat)
 function At_ldiv_B{Ts<:UMFVTypes,Tb<:Number}(S::SparseMatrixCSC{Ts}, b::Vector{Tb})
     ## should be more careful here in case Ts<:Real and Tb<:Complex
@@ -150,7 +150,7 @@ for (f_sym_r, f_num_r, f_sym_c, f_num_c, itype) in
             finalizer(U.symbolic,umfpack_free_symbolic)
             U
         end
-
+        
         function umfpack_symbolic!{Tv<:Complex128,Ti<:$itype}(U::UmfpackLU{Tv,Ti})
             if U.symbolic != C_NULL return U end
             tmp = Array(Ptr{Void},1)
@@ -180,7 +180,7 @@ for (f_sym_r, f_num_r, f_sym_c, f_num_c, itype) in
             finalizer(U.numeric,umfpack_free_numeric)
             U
         end
-
+        
         function umfpack_numeric!{Tv<:Complex128,Ti<:$itype}(U::UmfpackLU{Tv,Ti})
             if U.numeric != C_NULL return U end
             if U.symbolic == C_NULL umfpack_symbolic!(U) end
@@ -213,7 +213,7 @@ for (f_sol_r, f_sol_c, inttype) in
             if status != UMFPACK_OK; error("Error code $status in umfpack_solve"); end
             return x
         end
-
+        
         function umfpack_solve{Tv<:Complex128,Ti<:$inttype}(lu::UmfpackLU{Tv,Ti}, b::Vector{Tv}, typ::Integer)
             umfpack_numeric!(lu)
             xr = similar(b, Float64)
@@ -259,7 +259,7 @@ function umfpack_free_symbolic(lu::UmfpackLU)
     lu.numeric = C_NULL
     lu
 end
- 
+
 function umfpack_report_symbolic(symb::Ptr{Void}, level::Real)
     old_prl::Float64 = umf_ctrl[UMFPACK_PRL]
     umf_ctrl[UMFPACK_PRL] = float64(level)
@@ -306,12 +306,6 @@ ccall((:cholmod_start, :libcholmod), Int32, (Ptr{Uint8},), chm_com)
 
 const chm_l_com  = Array(Uint8, chm_com_sz)
 ccall((:cholmod_l_start, :libcholmod), Int32, (Ptr{Uint8},), chm_l_com)
-
-### These offsets should be reconfigured to be less error-prone in matches
-const chm_com_offsets = Array(Int, 20)
-ccall((:jl_cholmod_common_offsets, :libsuitesparse_wrapper),
-      Void, (Ptr{Uint},), chm_com_offsets)
-                                               
 ### A way of examining some of the fields in chm_com or chm_l_com
 type ChmCommon
     dbound::Float64
@@ -335,18 +329,17 @@ type ChmCommon
     dtype::Int32
 end
 
+### These offsets should be reconfigured to be less error-prone in matches
+const chm_com_offsets = Array(Int, length(ChmCommon.types))
+ccall((:jl_cholmod_common_offsets, :libsuitesparse_wrapper),
+      Void, (Ptr{Uint8},), chm_com_offsets)
+                                               
 ### there must be an easier way but at least this works.
 function ChmCommon(aa::Array{Uint8,1})
-    eval(Expr(:call,
-              unshift!({
-                        reinterpret(ChmCommon.types[i],
-                                    aa[(1:sizeof(ChmCommon.types[i])) +
-                                       chm_com_offsets[i]])[1]
-                        for i in 1:length(ChmCommon.types)
-                            },
-                            :ChmCommon),
-                            Any)
-                            )
+    typs = ChmCommon.types
+    sz = map(sizeof, typs)
+    args = map(i->reinterpret(typs[i], aa[chm_com_offsets[i] + (1:sz[i])])[1], 1:length(sz))
+    eval(Expr(:call, unshift!(args, :ChmCommon), Any))
 end
 
 function chm_itype{Tv,Ti<:CHMITypes}(S::SparseMatrixCSC{Tv,Ti})
@@ -363,7 +356,10 @@ function set_chm_prt_lev(cm::Array{Uint8}, lev::Integer)
     cm[(1:4) + chm_com_offsets[13]] = reinterpret(Uint8, [int32(lev)])
 end
 
-type CholmodDense{T<:CHMVTypes}
+## cholmod_dense pointers passed to or returned from C functions are of Julia type
+## Ptr{c_CholmodDense}.  The CholmodDense type contains a c_CholmodDense object and other
+## fields then ensure the memory pointed to is freed when it should be and not before.
+type c_CholmodDense{T<:CHMVTypes}
     m::Int
     n::Int
     nzmax::Int
@@ -372,28 +368,68 @@ type CholmodDense{T<:CHMVTypes}
     zpt::Ptr{Void}
     xtype::Int32
     dtype::Int32
-    aa::Array{T}
-    chm::Bool                       # was storage allocated by Cholmod?
 end
 
+type CholmodDense{T<:CHMVTypes}
+    c::c_CholmodDense
+    m::Matrix{T} # Array(T,(0,0)) when created from Ptr{c_CholmodDense{T}}
+    ptr::Ptr{c_CholmodDense{T}} # null pointer when created from Julia array
+end
+
+convert(::Type{c_CholmodDense}, d::CholmodDense) = d.c
+
 function CholmodDense{T<:CHMVTypes}(aa::VecOrMat{T})
     m = size(aa,1); n = size(aa,2)
-    CholmodDense{T}(m, n, m*n, stride(aa,2), convert(Ptr{T}, aa), C_NULL,
-                    int32(T<:Complex ? CHOLMOD_COMPLEX : CHOLMOD_REAL),
-                    int32(T<:Union(Float32,Complex64) ? CHOLMOD_SINGLE : CHOLMOD_DOUBLE),
-                    aa, false)
+    CholmodDense{T}(c_CholmodDense{T}(m, n, m*n, stride(aa,2), convert(Ptr{T}, aa), C_NULL,
+                                      int32(T<:Complex ? CHOLMOD_COMPLEX : CHOLMOD_REAL),
+                                      int32(T<:Union(Float32,Complex64) ? CHOLMOD_SINGLE : CHOLMOD_DOUBLE)),
+                    aa, convert(Ptr{c_CholmodDense{T}}, C_NULL))
 end
 
-function chm_chk_dn{T<:CHMVTypes}(cd::CholmodDense{T})
-    ccall((:cholmod_check_dense, :libcholmod), Int,
-          (Ptr{CholmodDense{T}}, Ptr{Uint8}), &cd, chm_com)
+function CholmodDense{T<:CHMVTypes}(c::Ptr{c_CholmodDense{T}})
+    cp = unsafe_ref(c)
+    if cp.lda != cp.m || cp.nzmax != cp.m * cp.n
+        error("overallocated cholmod_sparse returned object of size $(cp.m) by $(cp.n) with leading dim $(cp.lda) and nzmax $(cp.nzmax)")
+    end
+    CholmodDense{T}(cp, Array(T, (0,0)), c)
 end
 
+function convert{T<:CHMVTypes}(::Type{Array}, chmpt::Ptr{c_CholmodDense{T}}, own::Bool)
+    cp = unsafe_ref(chmpt)
+    if cp.lda != cp.m || cp.nzmax != cp.m * cp.n
+        error("overallocated cholmod_sparse returned object of size $(cp.m) by $(cp.n) with leading dim $(cp.lda) and nzmax $(cp.nzmax)")
+    end
+    arr = pointer_to_array(cp.xpt, (cp.m, cp.n), own)
+    if !own return copy(arr) end
+    c_free(chmpt)
+end
+convert{T<:CHMVTypes}(::Type{Array}, chmpt::Ptr{c_CholmodDense{T}}) = convert(Array, chmpt, false)
+    
 function chm_chk_dn{T<:CHMVTypes}(cd::CholmodDense{T})
     ccall((:cholmod_check_dense, :libcholmod), Int,
           (Ptr{CholmodDense{T}}, Ptr{Uint8}), &cd, chm_com)
 end
 
+function chm_ones(m::Integer, n::Integer, t::Float64)
+    ccall((:cholmod_ones, :libcholmod), Ptr{c_CholmodDense{Float64}},
+                  (Int, Int, Int32, Ptr{Uint8}),
+                  m, n, CHOLMOD_REAL, chm_com)
+end
+function chm_ones(m::Integer, n::Integer, t::Complex128)
+    ccall((:cholmod_ones, :libcholmod), Ptr{c_CholmodDense{Complex128}},
+          (Int, Int, Int32, Ptr{Uint8}),
+          m, n, CHOLMOD_COMPLEX, chm_com)
+end
+chm_ones(m::Integer, n::Integer) = chm_ones(m, n, 1.)
+
+function chm_free_dn{T<:UMFVTypes}(A::Ptr{c_CholmodDense{T}})
+    aa = [A]
+    status = ccall((:cholmod_free_dense, :libcholmod), Int32,
+                   (Ptr{Ptr{c_CholmodDense{T}}}, Ptr{Uint8}), aa, chm_com)
+    if status != CHOLMOD_TRUE throw(CholmodException) end
+end
+
+## This should be converted to a show method    
 function chm_prt_dn{T<:CHMVTypes}(cd::CholmodDense{T}, lev::Integer, nm::ASCIIString)
     inds = (1:4) + chm_com_offsets[13]
     orig = chm_com[inds]
@@ -435,395 +471,10 @@ function CholmodSparse{Tv<:CHMVTypes,Ti<:CHMITypes}(A::SparseMatrixCSC{Tv,Ti})
     nzval = copy(A.nzval)
     CholmodSparse{Tv,Ti}(size(A,1),size(A,2),int(colptr0[end]),
                          convert(Ptr{Ti}, colptr0), convert(Ptr{Ti}, rowval0),
-                         convert(Ptr{Tv}, nzval), C_NULL, int32(0),
+                         C_NULL, convert(Ptr{Tv}, nzval), C_NULL, int32(0),
                          chm_itype(A), chm_xtype(A), chm_dtype(A),
                          CHOLMOD_TRUE, CHOLMOD_TRUE,
                          colptr0, rowval0, nzval, false, C_NULL)
 end
 
-    
-if false    
-# Wrapper for memory allocated by CHOLMOD. Carry along the value and index types.
-## FIXME: CholmodPtr and UmfpackPtr should be amalgamated
-type CholmodPtr{Tv<:CHMVTypes,Ti<:CHMITypes}
-    val::Vector{Ptr{Void}} 
-end
-
-eltype{Tv,Ti}(P::CholmodPtr{Tv,Ti}) = Tv
-indtype{Tv,Ti}(P::CholmodPtr{Tv,Ti}) = Ti
-
-## function cholmod_common_finalizer(x::Vector{Ptr{Void}})
-##     st = ccall((:cholmod_finish, :libcholmod), BlasInt, (Ptr{Void},), x[1])
-##     if st != CHOLMOD_TRUE error("Error calling cholmod_finish") end
-##     c_free(x[1])
-## end
-
-## type CholmodCommon
-##     pt::Vector{Ptr{Void}}
-##     function CholmodCommon()
-##         pt = Array(Ptr{Void}, 1)
-##         ccall((:jl_cholmod_common, :libsuitesparse_wrapper), Void,
-##               (Ptr{Void},), pt)
-##         st = ccall((:cholmod_start, :libcholmod), BlasInt, (Ptr{Void}, ), pt[1])
-##         if st != CHOLMOD_TRUE error("Error calling cholmod_start") end
-##         finalizer(pt, cholmod_common_finalizer)
-##         new(pt)
-##     end
-## end
-
-## function show(io::IO, cm::CholmodCommon)
-##     st = ccall((:cholmod_print_common, :libcholmod), BlasInt,
-##                (Ptr{Uint8},Ptr{Void}), "", cm.pt[1])
-##     if st != CHOLMOD_TRUE error("Error calling cholmod_print_common") end
-## end
-
-type CholmodSparse{Tv<:CHMVTypes,Ti<:CHMITypes}
-    pt::CholmodPtr{Tv,Ti}
-    ## cp contains a copy of the original matrix but with 0-based indices
-    cp::SparseMatrixCSC{Tv,Ti}
-    stype::Int
-#    cm::CholmodCommon
-    function CholmodSparse(S::SparseMatrixCSC{Tv,Ti}, stype::BlasInt)#, cm::CholmodCommon)
-        pt = CholmodPtr{Tv,Ti}(Array(Ptr{Void}, 1))
-        cp = convert_to_0_based_indexing(S)
-        
-        ccall((:jl_cholmod_sparse, :libsuitesparse_wrapper), Void,
-              (Ptr{Void}, Uint, Uint, Uint, Ptr{Void}, Ptr{Void}, Ptr{Void},
-               Ptr{Void}, Ptr{Void}, BlasInt, BlasInt, BlasInt, BlasInt, BlasInt, Int),
-              pt.val, S.m, S.n, nnz(S), cp.colptr, cp.rowval, C_NULL,
-              cp.nzval, C_NULL, stype, chm_itype(S), chm_xtype(S), chm_dtype(S),
-              CHOLMOD_TRUE, CHOLMOD_TRUE)
-        finalizer(pt, x->c_free(x.val[1]))
-        new(pt, cp, blas_int(stype), cm)
-    end
-end
-
-CholmodSparse{Tv,Ti}(S::SparseMatrixCSC{Tv,Ti}, stype::Int) = CholmodSparse{Tv,Ti}(S, stype, CholmodCommon())
-
-function CholmodSparse{Tv,Ti}(S::SparseMatrixCSC{Tv,Ti})#, cm::CholmodCommon)
-    stype = S.m == S.n && ishermitian(S)
-    CholmodSparse{Tv,Ti}(stype ? triu(S) : S, blas_int(stype), cm)
-end
-
-CholmodSparse(S::SparseMatrixCSC) = CholmodSparse(S)#, CholmodCommon())
-
-function show(io::IO, cs::CholmodSparse)
-    ccall(_chm_print_sp,
-          BlasInt, (Ptr{Void}, Ptr{Uint8},Ptr{Void}), cs.pt.val[1], "", cs.cm.pt[1])
-end
-
-size(cs::CholmodSparse) = size(cs.cp)
-nnz(cs::CholmodSparse) = cs.cp.colptr[end]
-eltype{T}(cs::CholmodSparse{T}) = T
-indtype{Tv,Ti}(cs::CholmodSparse{Tv,Ti}) = Ti
-
-SparseMatrixCSC(cs::CholmodSparse) = convert_to_1_based_indexing(cs.cp)
-
-## For testing only.  The infinity and 1 norms of a sparse matrix are simply
-## the same norm applied to its nzval field.
-function norm(cs::CholmodSparse, p::Number)
-    ccall((:cholmod_norm_sparse, :libcholmod), Float64,
-          (Ptr{Void}, BlasInt, Ptr{Void}), cs.pt.val[1], p == Inf ? 0 : 1, cs.cm.pt[1])
-end
-
-norm(cs::CholmodSparse) = norm(cs, Inf)
-
-## Approximate minimal degree ordering
-function chm_amd(cs::CholmodSparse)
-    aa = Array(BlasInt, cs.cp.m)
-    st = cs.stype == 0 ? ccall(:cholmod_colamd, BlasInt,
-                               (Ptr{Void}, Ptr{Void}, Uint, BlasInt, Ptr{BlasInt}, Ptr{Void}),
-                               cs.pt.val[1], C_NULL, 0, 1, aa, cs.cm.pt[1]) :
-                         ccall(:cholmod_amd, BlasInt, (Ptr{Void}, Ptr{Void}, Uint, Ptr{BlasInt}, Ptr{Void}),
-                               cs.pt.val[1], C_NULL, 0, aa, cs.cm.pt[1])
-    if st != CHOLMOD_TRUE error("Error in cholmod_amd") end
-    aa
-end
-
-type CholmodFactor{Tv<:CHMVTypes,Ti<:CHMITypes} <: Factorization{Tv}
-    pt::CholmodPtr{Tv,Ti}
-    cs::CholmodSparse{Tv,Ti}
-    function CholmodFactor(pt::CholmodPtr{Tv,Ti}, cs::CholmodSparse{Tv,Ti})
-        ff = new(pt, cs)
-        finalizer(ff, cholmod_factor_finalizer)
-        ff
-    end
-end
-
-function cholmod_factor_finalizer(x::CholmodFactor)
-    if ccall((:cholmod_free_factor, :libcholmod), BlasInt, (Ptr{Void}, Ptr{Void}), x.pt.val, x.cs.cm[1]) != CHOLMOD_TRUE
-        error("CHOLMOD error in cholmod_free_factor")
-    end
-end
-
-function size(F::CholmodFactor)
-    n = size(F.cs,1)
-    (n, n)
-end
-
-eltype{T}(F::CholmodFactor{T}) = T
-indtype{Tv,Ti}(F::CholmodFactor{Tv,Ti}) = Ti
-
-function CholmodFactor{Tv,Ti}(cs::CholmodSparse{Tv,Ti})
-    pt = CholmodPtr{Tv,Ti}(Array(Ptr{Void}, 1))
-    pt.val[1] = ccall(:cholmod_analyze, Ptr{Void},
-                      (Ptr{Void}, Ptr{Void}), cs.pt.val[1], cs.cm.pt[1])
-    st = ccall(:cholmod_factorize, BlasInt,
-               (Ptr{Void}, Ptr{Void}, Ptr{Void}), cs.pt.val[1], pt.val[1], cs.cm.pt[1])
-    if st != CHOLMOD_TRUE error("CHOLMOD failure in factorize") end
-    CholmodFactor{Tv,Ti}(pt, cs)
-end
-
-function show(io::IO, cf::CholmodFactor)
-    st = ccall(:cholmod_print_fa, BlasInt, (Ptr{Void}, Ptr{Uint8}, Ptr{Void}), cf.pt.val[1], "", cf.cs.cm.pt[1])
-    if st != CHOLMOD_TRUE error("Cholmod error in print_factor") end
-end
-
-type CholmodDense{T<:CHMVTypes}
-    pt::Vector{Ptr{Void}}
-    m::Int
-    n::Int
-    aa::VecOrMat{T}                     # original array
-#    cm::CholmodCommon
-end
-
-function CholmodDense{T<:CHMVTypes}(b::VecOrMat{T}, cm::CholmodCommon)
-    m = size(b, 1)
-    n = isa(b, Matrix) ? size(b, 2) : 1
-
-    xtype = T <: Complex ? CHOLMOD_COMPLEX : CHOLMOD_REAL
-    dtype = T <: Float32 || T == Complex64 ? CHOLMOD_SINGLE : CHOLMOD_DOUBLE
-
-    pt = Array(Ptr{Void}, 1)
-
-    ccall((:jl_cholmod_dense, :libsuitesparse_wrapper), Void,
-          (Ptr{Void}, Uint, Uint, Uint, Uint, Ptr{Void}, Ptr{Void}, BlasInt, Int),
-          pt, m, n, length(b), m, b, C_NULL, xtype, dtype)
-    finalizer(pt, x->c_free(pt[1]))
-    CholmodDense{T}(pt, m, n, copy(b), cm)
-end
-
-CholmodDense{T<:Integer}(B::VecOrMat{T}, cm::CholmodCommon) = CholmodDense(float64(B), cm)
-    
-size(cd::CholmodDense) = (cd.m, cd.n)
-
-function show(io::IO, cd::CholmodDense)
-    st = ccall(:cholmod_print_dn, BlasInt, (Ptr{Void},Ptr{Uint8},Ptr{Void}), cd.pt[1], "", cd.cm.pt[1])
-    if st != CHOLMOD_TRUE error("Cholmod error in print_dense") end
-end
-
-type CholmodDenseOut{Tv<:CHMVTypes,Ti<:CHMITypes}
-    pt::CholmodPtr{Tv,Ti}
-    m::Int
-    n::Int
-    cm::CholmodCommon
-    function CholmodDenseOut(pt::CholmodPtr{Tv,Ti}, m::BlasInt, n::BlasInt, cm::CholmodCommon)
-        dd = new(pt, m, n, cm)
-        finalizer(dd, cholmod_denseout_finalizer)
-        dd
-    end
-end
-
-function cholmod_denseout_finalizer(cd::CholmodDenseOut)
-    st = ccall((:cholmod_free_dense, :libcholmod), BlasInt, (Ptr{Void}, Ptr{Void}), cd.pt.val, cd.cm.pt[1])
-    if st != CHOLMOD_TRUE error("Error in cholmod_free_dense") end
-end
-
-eltype{T}(cdo::CholmodDenseOut{T}) = T
-indtype{Tv,Ti}(cdo::CholmodDenseOut{Tv,Ti}) = Ti
-size(cd::CholmodDenseOut) = (cd.m, cd.n)
-
-function convert{T}(::Type{Array{T}}, cdo::CholmodDenseOut{T})
-    mm = Array(T, size(cdo))
-    ccall((:jl_cholmod_dense_copy_out, :libsuitesparse_wrapper), Void,
-          (Ptr{Void}, Ptr{T}), cdo.pt.val[1], mm)
-    mm
-end
-
-function solve{Tv,Ti}(cf::CholmodFactor{Tv,Ti}, B::CholmodDense{Tv}, solv::Integer)
-    m, n = size(B)
-    cdo = CholmodPtr{Tv,Ti}(Array(Ptr{Void},1))
-    cdo.val[1] = ccall(:cholmod_solve, Ptr{Void},
-                       (BlasInt, Ptr{Void}, Ptr{Void}, Ptr{Void}),
-                       solv, cf.pt.val[1], B.pt[1], cf.cs.cm.pt[1])
-    return cdo, m, n, cf.cs.cm
-    CholmodDenseOut(cdo, m, n, cf.cs.cm)
-end
-
-solve(cf::CholmodFactor, B::CholmodDense) = solve(cf, B, CHOLMOD_A)
-
-(\){Tf,Tb}(cf::CholmodFactor{Tf}, b::VecOrMat{Tb}) = solve(cf, CholmodDense{Tf}(convert(Array{Tf},b), cf.cs.cm), CHOLMOD_A)
-
-type CholmodSparseOut{Tv<:CHMVTypes,Ti<:CHMITypes}
-    pt::CholmodPtr{Tv,Ti}
-    m::Int
-    n::Int
-    cm::CholmodCommon
-    function CholmodSparseOut(pt::CholmodPtr{Tv,Ti}, m::BlasInt, n::BlasInt, cm::CholmodCommon)
-        cso = new(pt, m, n, cm)
-        finalizer(cso, cholmod_sparseout_finalizer)
-        cso
-    end
-end
-
-function cholmod_sparseout_finalizer(cso::CholmodSparseOut)
-    st = ccall((:cholmod_free_sparse, :libcholmod), BlasInt,
-               (Ptr{Void}, Ptr{Void}), cso.pt.val, cso.cm.pt[1])
-    if st != CHOLMOD_TRUE error("Error in cholmod_free_sparse") end
-end
-
-function nnz(cso::CholmodSparseOut)
-    ccall((:cholmod_nnz, :libcholmod), BlasInt,
-          (Ptr{Void}, Ptr{Void}), cso.pt.val[1], cso.cm.pt[1])
-end
-size(cso::CholmodSparseOut) = (cso.m, cso.n)
-eltype{T}(cso::CholmodSparseOut{T}) = T
-indtype{Tv,Ti}(cso::CholmodSparseOut{Tv,Ti}) = Ti
-
-function solve{Tv,Ti}(cf::CholmodFactor{Tv,Ti}, B::CholmodSparse{Tv,Ti}, solv::Integer)
-    m, n = size(B)
-    cso = CholmodPtr{Tv,Ti}(Array(Ptr{Void},1))
-    cso.val[1] = ccall((:cholmod_spsolve, :libcholmod), Ptr{Void},
-                       (BlasInt, Ptr{Void}, Ptr{Void}, Ptr{Void}),
-                       solv, cf.pt.val[1], B.pt[1], B.cm.pt[1])
-    CholmodSparseOut{Tv,Ti}(cso, m, n, cf.cs.cm)
-end
-
-function CholmodSparseOut{Tv,Ti}(cf::CholmodFactor{Tv,Ti})
-    n = size(cf.cs)[1]
-    cso = CholmodPtr{Tv,Ti}(Array(Ptr{Void},1))
-    cso.val[1] = ccall((:cholmod_factor_to_sparse, :libcholmod), Ptr{Void},
-                       (Ptr{Void}, Ptr{Void}), cf.pt.val[1], cf.cs.cm.pt[1])
-    CholmodSparseOut{Tv,Ti}(cso, n, n, cf.cs.cm)
-end
-
-function SparseMatrixCSC{Tv,Ti}(cso::CholmodSparseOut{Tv,Ti})
-    nz = nnz(cso)
-    sp = SparseMatrixCSC{Tv,Ti}(cso.m, cso.n, Array(Ti, cso.n + 1), Array(Ti, nz), Array(Tv, nz))
-    st = ccall((:jl_cholmod_sparse_copy_out, :libsuitesparse_wrapper), BlasInt,
-                (Ptr{Void}, Ptr{Ti}, Ptr{Ti}, Ptr{Tv}),
-                cso.pt.val[1], sp.colptr, sp.rowval, sp.nzval)
-    if st == 1 error("CholmodSparseOut object is not packed") end
-    if st == 2 error("CholmodSparseOut object is not sorted") end # Should not occur
-    if st == 3 error("CholmodSparseOut object has INTLONG itype") end
-    convert_to_1_based_indexing!(sp)
-end
-
-function show(io::IO, cso::CholmodSparseOut)
-    sp = ccall(:cholmod_print_sp, BlasInt, (Ptr{Void}, Ptr{Uint8},Ptr{Void}), cso.pt.val[1], "", cso.cm.pt[1])
-    if sp != CHOLMOD_TRUE error("Cholmod error in print_sparse") end
-end
-
-function chm_aat{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, symm::Bool)
-    cs = CholmodSparse(A, 0)
-    aa = CholmodPtr{Tv,Ti}(Array(Ptr{Void},1))
-    aa.val[1] = ccall(:cholmod_aat, Ptr{Void}, (Ptr{Void},Ptr{BlasInt},BlasInt,BlasInt,Ptr{Void}),
-                      cs.pt.val[1], C_NULL, 0, 1, cs.cm.pt[1])
-    if ccall(:cholmod_sort, BlasInt, (Ptr{Void}, Ptr{Void}), aa.val[1], cs.cm.pt[1]) != CHOLMOD_TRUE
-        error("Cholmod error in sort")
-    end
-    if symm
-        pt = ccall(:cholmod_copy, Ptr{Void}, (Ptr{Void}, BlasInt, BlasInt, Ptr{Void}),
-                   aa.val[1], 1, 1, cs.cm.pt[1])
-        if ccall((:cholmod_free_sparse, :libcholmod), BlasInt,
-                 (Ptr{Void}, Ptr{Void}), aa.val, cs.cm.pt[1]) != CHOLMOD_TRUE
-            error("Cholmod error in free_sparse")
-        end
-        aa.val[1] = pt
-    end
-    m  = size(A, 1)
-    CholmodSparseOut{Tv,Ti}(aa, m, m, cs.cm)
-end
-
-chm_aat{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}) = chm_aat(A, false)
-
-## call wrapper function to create cholmod_sparse objects
-cholmod_sparse(S) = cholmod_sparse(S, 0)
-
-function cholmod_sparse{Tv,Ti}(S::SparseMatrixCSC{Tv,Ti}, stype::Int)
-    cs = Array(Ptr{Void}, 1)
-
-    if     Ti == Int; itype = CHOLMOD_INT;
-    elseif Ti == Int64; itype = CHOLMOD_LONG; end
-
-    if     Tv == Float64    || Tv == Float32;    xtype = CHOLMOD_REAL;
-    elseif Tv == Complex128 || Tv == Complex64 ; xtype = CHOLMOD_COMPLEX; end
-
-    if     Tv == Float64 || Tv == Complex128; dtype = CHOLMOD_DOUBLE; 
-    elseif Tv == Float32 || Tv == Complex64 ; dtype = CHOLMOD_SINGLE; end
-
-    ccall((:jl_cholmod_sparse, :libsuitesparse_wrapper),
-          Ptr{Void},
-          (Ptr{Void}, BlasInt, BlasInt, BlasInt, Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void},
-           BlasInt, BlasInt, BlasInt, BlasInt, BlasInt, Int),
-          cs, blas_int(S.m), blas_int(S.n), blas_int(length(S.nzval)), S.colptr, S.rowval, C_NULL, S.nzval, C_NULL,
-          int32(stype), itype, xtype, dtype, CHOLMOD_TRUE, CHOLMOD_TRUE
-          )
-
-    return cs
-end
-
-## Call wrapper function to create cholmod_dense objects
-function cholmod_dense{T}(B::VecOrMat{T})
-    m = size(B, 1)
-    n = isa(B, Matrix) ? size(B, 2) : 1
-
-    cd = Array(Ptr{Void}, 1)
-
-    if     T == Float64    || T == Float32;    xtype = CHOLMOD_REAL;
-    elseif T == Complex128 || T == Complex64 ; xtype = CHOLMOD_COMPLEX; end
-
-    if     T == Float64 || T == Complex128; dtype = CHOLMOD_DOUBLE; 
-    elseif T == Float32 || T == Complex64 ; dtype = CHOLMOD_SINGLE; end
-
-    ccall((:jl_cholmod_dense, :libsuitesparse_wrapper), Ptr{Void},
-          (Ptr{Void}, BlasInt, BlasInt, BlasInt, BlasInt, Ptr{T}, Ptr{Void}, BlasInt, Int),
-          cd, size(B,1), size(B,2), length(B), stride(B,2), B, C_NULL, xtype, dtype
-          )
-
-    return cd
-end
-
-function cholmod_dense_copy_out{T}(x::Ptr{Void}, sol::VecOrMat{T})
-    ccall((:jl_cholmod_dense_copy_out, :libsuitesparse_wrapper),
-          Void,
-          (Ptr{Void}, Ptr{T}),
-          x, sol
-          )
-    return sol
-end
-
-function cholmod_transpose_unsym{Tv,Ti}(S::SparseMatrixCSC{Tv,Ti}, cm::Array{Ptr{Void}, 1})
-    S_t = SparseMatrixCSC(Tv, S.n, S.m, nnz(S)+1)
-
-    # Allocate space for a cholmod_sparse object
-    cs = cholmod_sparse(S)
-    cs_t = cholmod_sparse(S_t)
-    
-    status = ccall((:cholmod_transpose_unsym),
-                   Int32,
-                   (Ptr{Void}, BlasInt, Ptr{BlasInt}, Ptr{BlasInt}, BlasInt, Ptr{Void}, Ptr{Void}),
-                   cs[1], int32(1), C_NULL, C_NULL, int32(-1), cs_t[1], cm[1]);
-
-    # Deallocate space for cholmod_sparse objects
-    c_free(cs[1])
-    c_free(cs_t[1])
-
-    return S_t
-end
-
-function cholmod_analyze{Tv<:CHMVTypes, Ti<:CHMITypes}(cs::Array{Ptr{Void},1}, cm::Array{Ptr{Void},1})
-    ccall(:cholmod_analyze, Ptr{Void}, (Ptr{Void}, Ptr{Void}), cs[1], cm[1])
-end
-
-function cholmod_factorize{Tv<:CHMVTypes, Ti<:CHMITypes}(cs::Array{Ptr{Void},1}, cs_factor::Ptr{Void}, cm::Array{Ptr{Void},1})
-    st = ccall(:cholmod_factorize, BlasInt, (Ptr{Void}, Ptr{Void}, Ptr{Void}), cs[1], cs_factor, cm[1])
-    if st != CHOLMOD_TRUE error("CHOLMOD could not factorize the matrix") end
-end
-
-function cholmod_solve(cs_factor::Ptr{Void}, cd_rhs::Array{Ptr{Void},1}, cm::Array{Ptr{Void},1})
-    ccall(:cholmod_solve, Ptr{Void}, (BlasInt, Ptr{Void}, Ptr{Void}, Ptr{Void}),
-          CHOLMOD_A, cs_factor, cd_rhs[1], cm[1])
-end
-end # code omission
 end #module