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.
-- 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
