Yichao,
That's great, thanks very much! It's almost working now (revised code
below). I still have two issues, which are perhaps best answered by
someone familiar with SPQR. First, although the comments in the SPQR
package say that E is "column permutation", in fact I am getting back a
vector containing many 1's (and a few other integers in the range 2...n),
so it is apparently not a permutation. Anyone know what E is? I couldn't
find anything about E in the PDF manuals for SPQR.
Also, I don't know how to free the memory associated with E (a C array of
integers).
-- Steve
function qr_fac_get_r(A, tol, ordering =
Base.SparseMatrix.SPQR.ORDERING_DEFAULT)
Aw = Sparse(A,0)
s = unsafe_load(Aw.p)
if s.stype != 0
throw(ArgumentError("stype must be zero"))
end
ptr_R = Ref{Ptr{C_Sparse{Float64}}}()
ptr_E = Ref{Ptr{Clong}}()
rk = ccall((:SuiteSparseQR_C, :libspqr), Clong,
(Cint, #ordering
Cdouble, #tol
Clong, #econ
Cint, #getCTX
Ptr{C_Sparse{Float64}}, # A
Ptr{Void}, #Bsparse
Ptr{Void}, #Bdense
Ptr{Void}, #Zsparse
Ptr{Void}, #Zdense
Ptr{Void}, #R
Ptr{Void}, #E
Ptr{Void}, #H
Ptr{Void}, #HPInv
Ptr{Void}, #HTau
Ptr{Void}), #cc
ordering, #ordering
tol, #tol
1000000000, #econ
0, #getCTX
get(Aw.p), # A
C_NULL, #Bsparse
C_NULL, #Bdense
C_NULL, #Zsparse
C_NULL, #Zdense
ptr_R, #R
ptr_E, #E
C_NULL, #H
C_NULL, #HPInv
C_NULL, #HTau
common()) #cc
println("rk = ", rk)
R = ptr_R[]
E = ptr_E[]
if E != C_NULL
colprm = pointer_to_array(E, size(A,2), false) + 1
else
colprm = collect(1:size(A,2))
end
R1 = unsafe_load(R)
if R1.stype != 0
throw(ArgumentError("matrix has stype != 0. Convert to matrix with
stype == 0 before converting to SparseMatrixCSC"))
end
R_cvt = SparseMatrixCSC(R1.nrow,
R1.ncol,
increment(pointer_to_array(R1.p, (R1.ncol +
1,), false)),
increment(pointer_to_array(R1.i, (R1.nzmax,),
false)),
copy(pointer_to_array(R1.x, (R1.nzmax,),
false)))
free_sparse!(R)
R_cvt, colprm
end
On Thursday, August 4, 2016 at 11:00:55 PM UTC-4, Yichao Yu wrote:
>
>
> On Fri, Aug 5, 2016 at 10:47 AM, <[email protected] <javascript:>>
> wrote:
>
>> In an earlier post, I asked how to use spqr (C language subroutines that
>> are part of the Julia base) to obtain the R factor of a sparse QR
>> factorization. Based on the lack of response, I assume that this has not
>> been tried before. So now I'm requesting help on how to use ccall, which I
>> have not used previously. The C-routine that I think I need is below (in
>> the Julia comments). This C function, as far as I understand, returns its
>> outputs via pointers-to-pointers. In particular, the caller is supposed to
>> pass in a location in which a pointer to the answer is to be stored, or
>> else 0 if the caller discards that output.
>>
>> My wrapper function below did not work. The C-function appears to have
>> run, because rk has the correct value. But the values printed for nrow,
>> ncol, nzmax are garbage, i.e., the code does not correctly retrieve R.
>> Furthermore, the call to finalizer gave an error message ("Objects of type
>> Ptr{...} cannot be finalized"). I cribbed the finalizer invocations from
>> spqr.jl but evidently not correctly.
>>
>
> It's really unclear from your code below what you are trying to do. (maybe
> someone with more experience with suitsparse would). You'll be better off
> trying ccall on a few simple functions.
>
> From your description above, I'm guessing you have a c function that looks
> like
>
> ```
> void f(void **p)
> {
> *p = new object;
> }
> ```
>
> The way to call this is
>
> ```
> pp = Ref{Ptr{Void}}()
> ccall(:f, Void, (Ptr{Ptr{Void}},), pp)
> p = pp[]
> ```
>
> You can't add finalizer to a bare pointer. There are mutiple ways to
> manage memory returned by C functions.
>
> 1. If the pointer is pointing to a buffer of known size and can be free
> with `free`, use `unsafe_wrap(Array)` to get an array
> 2. If the pointer is merely used to return a small amount of information,
> you can eagerly convert that info to a julia object
> 3. If the pointer is an opaque object or if you need exactly the same
> pointer later for interacting with the C library, put the pointer in a
> mutable type. Define `unsafe_convert` and `cconvert` on that type to
> convert the julia object to pointer argument and add finalizer on the julia
> object that frees the pointer.
>
>
>
>> -- Steve Vavasis
>>
>>
>>
>> module test_qr_fac
>>
>>
>> import Base.SparseMatrix.CHOLMOD.common
>> import Base.SparseMatrix.CHOLMOD.C_Sparse
>> import Base.SparseMatrix.CHOLMOD.Sparse
>> import Base.SparseMatrix.CHOLMOD.free!
>>
>>
>> #SuiteSparse_long SuiteSparseQR_C /* returns rank(A) estimate, (-1) if
>> failure */
>> #(
>> # /* inputs: */
>> # int ordering, /* all, except 3:given treated as
>> 0:fixed */
>> # double tol, /* columns with 2-norm <= tol treated as
>> 0 */
>> # SuiteSparse_long econ, /* e = max(min(m,econ),rank(A)) */
>> # int getCTX, /* 0: Z=C (e-by-k), 1: Z=C', 2: Z=X
>> (e-by-k) */
>> # cholmod_sparse *A, /* m-by-n sparse matrix to factorize */
>> # cholmod_sparse *Bsparse, /* sparse m-by-k B */
>> # cholmod_dense *Bdense, /* dense m-by-k B */
>> # /* outputs: */
>> # cholmod_sparse **Zsparse, /* sparse Z */
>> # cholmod_dense **Zdense, /* dense Z */
>> # cholmod_sparse **R, /* e-by-n sparse matrix */
>> # SuiteSparse_long **E, /* size n column perm, NULL if identity
>> */
>> # cholmod_sparse **H, /* m-by-nh Householder vectors */
>> # SuiteSparse_long **HPinv, /* size m row permutation */
>> # cholmod_dense **HTau, /* 1-by-nh Householder coefficients */
>> # cholmod_common *cc /* workspace and parameters */
>> #) ;
>>
>> function qr_fac_get_r(A, tol, ordering =
>> Base.SparseMatrix.SPQR.ORDERING_DEFAULT)
>> Aw = Sparse(A,0)
>> s = unsafe_load(Aw.p)
>> if s.stype != 0
>> throw(ArgumentError("stype must be zero"))
>> end
>> R = Ptr{C_Sparse{Float64}}(C_NULL)
>> ptr_R = Ptr{Void}(R)
>> E = Ptr{Clong}(C_NULL)
>> ptr_E = Ptr{Void}(E)
>>
>> rk = ccall((:SuiteSparseQR_C, :libspqr), Clong,
>> (Cint, #ordering
>> Cdouble, #tol
>> Clong, #econ
>> Cint, #getCTX
>> Ptr{C_Sparse{Float64}}, # A
>> Ptr{Void}, #Bsparse
>> Ptr{Void}, #Bdense
>> Ptr{Void}, #Zsparse
>> Ptr{Void}, #Zdense
>> Ptr{Void}, #R
>> Ptr{Void}, #E
>> Ptr{Void}, #H
>> Ptr{Void}, #HPInv
>> Ptr{Void}, #HTau
>> Ptr{Void}), #cc
>> ordering, #ordering
>> tol, #tol
>> 1000000000, #econ
>> 0, #getCTX
>> get(Aw.p), # A
>> C_NULL, #Bsparse
>> C_NULL, #Bdense
>> C_NULL, #Zsparse
>> C_NULL, #Zdense
>> ptr_R, #R
>> ptr_E, #E
>> C_NULL, #H
>> C_NULL, #HPInv
>> C_NULL, #HTau
>> common()) #cc
>> println("rk = ", rk)
>> R1 = unsafe_load(R)
>> println("nrow = ", R1.nrow, " ncol = ", R1.ncol, " nzmax = ",
>> R1.nzmax)
>> finalizer(R, free!)
>> finalizer(E, free!)
>> nothing
>> end
>>
>>
>> end
>>
>>
>
