On Fri, Aug 5, 2016 at 10:47 AM, <vava...@uwaterloo.ca> 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
>
>
