#
# Julia interface to BASICLU
#
module basiclu
using LinearAlgebra
using LinearAlgebra: checksquare
using Printf
using SparseArrays
using SparseArrays: getcolptr
using Libdl
if haskey(ENV, "JULIA_BASICLU_LIBRARY_PATH")
const libbasiclu = joinpath(ENV["JULIA_BASICLU_LIBRARY_PATH"], "libbasiclu.$dlext")
else
using basiclu_jll
end
# comment out if the MAT module is installed
# include("test.jl")
# status codes
const BASICLU_OK = 0
const BASICLU_WARNING_SINGULAR_MATRIX = 2
const BASICLU_ERROR_INVALID_STORE = -1
const BASICLU_ERROR_INVALID_CALL = -2
const BASICLU_ERROR_ARGUMENT_MISSING = -3
const BASICLU_ERROR_INVALID_ARGUMENT = -4
const BASICLU_ERROR_MAXIMUM_UPDATES = -5
const BASICLU_ERROR_SINGULAR_UPDATE = -6
const BASICLU_ERROR_INVALID_OBJECT = -8
const BASICLU_ERROR_OUT_OF_MEMORY = -9
# user parameters in xstore (zero based)
const BASICLU_DROP_TOLERANCE = 4
const BASICLU_ABS_PIVOT_TOLERANCE = 5
const BASICLU_REL_PIVOT_TOLERANCE = 6
const BASICLU_BIAS_NONZEROS = 7
const BASICLU_MAXN_SEARCH_PIVOT = 8
const BASICLU_PAD = 9
const BASICLU_STRETCH = 10
const BASICLU_COMPRESSION_THRESHOLD = 11
const BASICLU_SPARSE_THRESHOLD = 12
const BASICLU_REMOVE_COLUMNS = 13
const BASICLU_SEARCH_ROWS = 14
# user readable from xstore (zero based)
const BASICLU_DIM = 64
const BASICLU_STATUS = 65
const BASICLU_NUPDATE = 70
const BASICLU_NFORREST = 71
const BASICLU_NFACTORIZE = 72
const BASICLU_NUPDATE_TOTAL = 73
const BASICLU_NFORREST_TOTAL = 74
const BASICLU_NSYMPERM_TOTAL = 75
const BASICLU_LNZ = 76
const BASICLU_UNZ = 77
const BASICLU_RNZ = 78
const BASICLU_MIN_PIVOT = 79
const BASICLU_MAX_PIVOT = 80
const BASICLU_UPDATE_COST = 81
const BASICLU_TIME_FACTORIZE = 82
const BASICLU_TIME_SOLVE = 83
const BASICLU_TIME_UPDATE = 84
const BASICLU_TIME_FACTORIZE_TOTAL = 85
const BASICLU_TIME_SOLVE_TOTAL = 86
const BASICLU_TIME_UPDATE_TOTAL = 87
const BASICLU_LFLOPS = 88
const BASICLU_UFLOPS = 89
const BASICLU_RFLOPS = 90
const BASICLU_CONDEST_L = 91
const BASICLU_CONDEST_U = 92
const BASICLU_MAX_ETA = 93
const BASICLU_NORM_L = 94
const BASICLU_NORM_U = 95
const BASICLU_NORMEST_LINV = 96
const BASICLU_NORMEST_UINV = 97
const BASICLU_MATRIX_ONENORM = 98
const BASICLU_MATRIX_INFNORM = 99
const BASICLU_RESIDUAL_TEST = 111
const BASICLU_MATRIX_NZ = 100
const BASICLU_RANK = 101
const BASICLU_BUMP_SIZE = 102
const BASICLU_BUMP_NZ = 103
const BASICLU_NSEARCH_PIVOT = 104
const BASICLU_NEXPAND = 105
const BASICLU_NGARBAGE = 106
const BASICLU_FACTOR_FLOPS = 107
const BASICLU_TIME_SINGLETONS = 108
const BASICLU_TIME_SEARCH_PIVOT = 109
const BASICLU_TIME_ELIM_PIVOT = 110
const BASICLU_PIVOT_ERROR = 120
mutable struct Param
# Factorization parameters
absPivotTol::Float64
relPivotTol::Float64
dropTol::Float64
biasNz::Int64
maxSearch::Int64
searchRows::Int64
removeCols::Int64
memPad::Int64
memStretch::Float64
memCompress::Float64
# Solve parameters
sparseThres::Float64
function Param()
param = new()
for f in fieldnames(Param)
setfield!(param, f, convert(fieldtype(Param, f), 0))
end
param
end
end
mutable struct Info
dim::Int64
status::Int64
# Info from last factorization
nPivot::Int64
nMatElem::Int64
nBumpCol::Int64
nBumpElem::Int64
nSearchPivot::Int64
nExpand::Int64
nGarbage::Int64
nFactorFlop::Int64
timeFactorize::Float64
timeSingletons::Float64
timeSearchPivot::Float64
timeElimPivot::Float64
residualTest::Float64
onenormMatrix::Float64
infnormMatrix::Float64
normL::Float64
normU::Float64
normestLinv::Float64
normestUinv::Float64
condestL::Float64
condestU::Float64
# Info from solves and updates since last factorization
nUpdate::Int64
nForrest::Int64
nElemL::Int64
nElemU::Int64
nElemR::Int64
nFlopL::Int64
nFlopU::Int64
nFlopR::Int64
updateCost::Float64
minPivot::Float64
maxPivot::Float64
maxEta::Float64
pivotError::Float64
timeSolve::Float64
timeUpdate::Float64
# Accumulated info since object was created
nFactorize::Int64
nUpdateTotal::Int64
nForrestTotal::Int64
nSympermTotal::Int64
timeFactorizeTotal::Float64
timeSolveTotal::Float64
timeUpdateTotal::Float64
function Info()
info = new()
for f in fieldnames(Info)
setfield!(info, f, convert(fieldtype(Info, f), 0))
end
info
end
end
function paramnumber(name::Symbol)
if name == :absPivotTol
return BASICLU_ABS_PIVOT_TOLERANCE
elseif name == :relPivotTol
return BASICLU_REL_PIVOT_TOLERANCE
elseif name == :dropTol
return BASICLU_DROP_TOLERANCE
elseif name == :biasNz
return BASICLU_BIAS_NONZEROS
elseif name == :maxSearch
return BASICLU_MAXN_SEARCH_PIVOT
elseif name == :searchRows
return BASICLU_SEARCH_ROWS
elseif name == :removeCols
return BASICLU_REMOVE_COLUMNS
elseif name == :memPad
return BASICLU_PAD
elseif name == :memStretch
return BASICLU_STRETCH
elseif name == :memCompress
return BASICLU_COMPRESSION_THRESHOLD
elseif name == :sparseThres
return BASICLU_SPARSE_THRESHOLD
end
throw(ErrorException("unknown parameter name"))
end
function infonumber(name::Symbol)
if name == :dim
return BASICLU_DIM
elseif name == :status
return BASICLU_STATUS
elseif name == :nPivot
return BASICLU_RANK
elseif name == :nMatElem
return BASICLU_MATRIX_NZ
elseif name == :nBumpCol
return BASICLU_BUMP_SIZE
elseif name == :nBumpElem
return BASICLU_BUMP_NZ
elseif name == :nSearchPivot
return BASICLU_NSEARCH_PIVOT
elseif name == :nExpand
return BASICLU_NEXPAND
elseif name == :nGarbage
return BASICLU_NGARBAGE
elseif name == :nFactorFlop
return BASICLU_FACTOR_FLOPS
elseif name == :timeFactorize
return BASICLU_TIME_FACTORIZE
elseif name == :timeSingletons
return BASICLU_TIME_SINGLETONS
elseif name == :timeSearchPivot
return BASICLU_TIME_SEARCH_PIVOT
elseif name == :timeElimPivot
return BASICLU_TIME_ELIM_PIVOT
elseif name == :residualTest
return BASICLU_RESIDUAL_TEST
elseif name == :onenormMatrix
return BASICLU_MATRIX_ONENORM
elseif name == :infnormMatrix
return BASICLU_MATRIX_INFNORM
elseif name == :normL
return BASICLU_NORM_L
elseif name == :normU
return BASICLU_NORM_U
elseif name == :normestLinv
return BASICLU_NORMEST_LINV
elseif name == :normestUinv
return BASICLU_NORMEST_UINV
elseif name == :condestL
return BASICLU_CONDEST_L
elseif name == :condestU
return BASICLU_CONDEST_U
elseif name == :nUpdate
return BASICLU_NUPDATE
elseif name == :nForrest
return BASICLU_NFORREST
elseif name == :nElemL
return BASICLU_LNZ
elseif name == :nElemU
return BASICLU_UNZ
elseif name == :nElemR
return BASICLU_RNZ
elseif name == :nFlopL
return BASICLU_LFLOPS
elseif name == :nFlopU
return BASICLU_UFLOPS
elseif name == :nFlopR
return BASICLU_RFLOPS
elseif name == :updateCost
return BASICLU_UPDATE_COST
elseif name == :minPivot
return BASICLU_MIN_PIVOT
elseif name == :maxPivot
return BASICLU_MAX_PIVOT
elseif name == :maxEta
return BASICLU_MAX_ETA
elseif name == :pivotError
return BASICLU_PIVOT_ERROR
elseif name == :timeSolve
return BASICLU_TIME_SOLVE
elseif name == :timeUpdate
return BASICLU_TIME_UPDATE
elseif name == :nFactorize
return BASICLU_NFACTORIZE
elseif name == :nUpdateTotal
return BASICLU_NUPDATE_TOTAL
elseif name == :nForrestTotal
return BASICLU_NFORREST_TOTAL
elseif name == :nSympermTotal
return BASICLU_NSYMPERM_TOTAL
elseif name == :timeFactorizeTotal
return BASICLU_TIME_FACTORIZE_TOTAL
elseif name == :timeSolveTotal
return BASICLU_TIME_SOLVE_TOTAL
elseif name == :timeUpdateTotal
return BASICLU_TIME_UPDATE_TOTAL
end
throw(ErrorException("unknown info name"))
end
mutable struct basiclu_object
istore::Ptr{Int64}
xstore::Ptr{Float64}
Li::Ptr{Int64}
Ui::Ptr{Int64}
Wi::Ptr{Int64}
Lx::Ptr{Float64}
Ux::Ptr{Float64}
Wx::Ptr{Float64}
lhs::Ptr{Float64}
ilhs::Ptr{Int64}
nzlhs::Int64
realloc_factor::Float64
function basiclu_object(dim::Int64)
obj = new()
retcode = ccall((:basiclu_obj_initialize, libbasiclu), Int64,
(Ptr{basiclu_object}, Int64),
pointer_from_objref(obj), dim)
checkretcode("basiclu_obj_initialize", retcode)
finalizer(obj) do x
ccall((:basiclu_obj_free, libbasiclu), Cvoid,
(Ptr{basiclu_object},), pointer_from_objref(x))
end
end
end
function getdim(obj::basiclu_object)
ccall((:basiclu_obj_get_dim, libbasiclu), Int64,
(Ptr{basiclu_object},), pointer_from_objref(obj))
end
function getparam(obj::basiclu_object, name::Symbol)
val = 0.0
if obj.xstore != C_NULL
val = unsafe_load(obj.xstore, paramnumber(name) + 1)
end
convert(fieldtype(Param, name), val)
end
function getparam(obj::basiclu_object)
param = Param()
for f in fieldnames(Param)
setfield!(param, f, getparam(obj, f))
end
param
end
function setparam!(obj::basiclu_object, name::Symbol, val)
if obj.xstore != C_NULL
unsafe_store!(obj.xstore, val, paramnumber(name) + 1)
end
end
function setparam!(obj::basiclu_object, param::Param)
for f in fieldnames(Param)
setparam!(obj, f, getfield(param, f))
end
end
function getinfo(obj::basiclu_object, name::Symbol)
val = 0.0
if obj.xstore != C_NULL
val = unsafe_load(obj.xstore, infonumber(name) + 1)
end
convert(fieldtype(Info, name), val)
end
function getinfo(obj::basiclu_object)
info = Info()
for f in fieldnames(Info)
setfield!(info, f, getinfo(obj, f))
end
info
end
"""
factorize(obj::basiclu_object, B::SparseMatrixCSC{Float64, Int64}; check::Bool=true)
Factorize sparse matrix `B`, which must be square and have the dimension for
which `obj` was created.
The factorization method stops when all elements of the active submatrix are
below the absolute pivot tolerance. In this case the `L` and `U` matrix are
padded with columns of the identity matrix, yielding the factorization of a
matrix `B*` which equals `B` except that the dependent columns are replaced by
columns of the identity matrix. The number of actual pivot steps performed can
be obtained with `getinfo(obj, :nPivot)`.
When `check = true`, an error is thrown if the number of pivot steps is less
than the dimension of `B`.
"""
function factorize(obj::basiclu_object, B::SparseMatrixCSC{Float64, Int64}; check::Bool=true)
dim = getdim(obj)
m = checksquare(B)
if dim != m
throw(DimensionMismatch("matrix dimension does not match basiclu object"))
end
if dim != 0
Bp = getcolptr(B) .- 1
Bi = rowvals(B) .- 1
Bx = nonzeros(B) # don't need a copy
retcode = ccall((:basiclu_obj_factorize, libbasiclu), Int64,
(Ptr{basiclu_object}, Ptr{Int64}, Ptr{Int64}, Ptr{Int64}, Ptr{Float64}),
pointer_from_objref(obj), Bp, pointer(Bp, 2), Bi, Bx)
if retcode == BASICLU_WARNING_SINGULAR_MATRIX && !check
retcode = BASICLU_OK
end
checkretcode("basiclu_obj_factorize", retcode)
end
end
"""
getfactors(obj::basiclu_object) -> L::SparseMatrixCSC{Float64, Int64},
U::SparseMatrixCSC{Float64, Int64},
p::Vector{Int64},
q::Vector{Int64}
Extract LU factors after fresh factorization. `L` is unit lower triangular, `U`
is upper triangular and `p` and `q` are permutation vectors such that (ignoring
round-off errors) `B[p,q]=LU` when matrix `B` was factorizied.
"""
function getfactors(obj::basiclu_object)
if getinfo(obj, :nUpdate) != 0
throw(ErrorException("cannot extract LU factors after factorization has been updated"))
end
dim = getdim(obj)
if dim == 0
L = spzeros(0, 0)
U = spzeros(0, 0)
rowperm = zeros(Int64, 0)
colperm = zeros(Int64, 0)
return L, U, rowperm, colperm
end
info = getinfo(obj)
rowperm = Vector{Int64}(undef, dim)
colperm = Vector{Int64}(undef, dim)
Lcolptr = Vector{Int64}(undef, dim + 1)
Lrowidx = Vector{Int64}(undef, info.nElemL + dim)
Lvalue = Vector{Float64}(undef, info.nElemL + dim)
Ucolptr = Vector{Int64}(undef, dim + 1)
Urowidx = Vector{Int64}(undef, info.nElemU + dim)
Uvalue = Vector{Float64}(undef, info.nElemU + dim)
retcode = ccall((:basiclu_obj_get_factors, libbasiclu), Int64,
(Ptr{basiclu_object}, Ptr{Int64}, Ptr{Int64},
Ptr{Int64}, Ptr{Int64}, Ptr{Float64},
Ptr{Int64}, Ptr{Int64}, Ptr{Float64}),
pointer_from_objref(obj), rowperm, colperm,
Lcolptr, Lrowidx, Lvalue,
Ucolptr, Urowidx, Uvalue)
if retcode == BASICLU_ERROR_INVALID_CALL
throw(ErrorException("no fresh factorization available"))
end
checkretcode("basiclu_get_factors", retcode)
L = SparseMatrixCSC{Float64, Int64}(dim, dim, Lcolptr .+ 1, Lrowidx .+ 1, Lvalue)
U = SparseMatrixCSC{Float64, Int64}(dim, dim, Ucolptr .+ 1, Urowidx .+ 1, Uvalue)
return L, U, rowperm .+ 1, colperm .+ 1
end
"""
solve(obj::basiclu_object, rhs::Vector{Float64}, trans::AbstractChar) -> Vector{Float64}
Solve linear system with dense right-hand side. `rhs` is not modified. `trans`
must be 'T' for transposed solve or 'N' for forward solve.
"""
function solve(obj::basiclu_object, rhs::Vector{Float64}, trans::AbstractChar)
lhs = copy(rhs)
solve!(obj, lhs, trans)
end
"""
solve!(obj::basiclu_object, rhs::Vector{Float64}, trans::AbstractChar) -> Vector{Float64}
Solve linear system with dense right-hand side. Solution overwrites `rhs`.
`trans` must be 'T' for transposed solve or 'N' for forward solve.
"""
function solve!(obj::basiclu_object, rhs::Vector{Float64}, trans::AbstractChar)
checktrans(trans)
dim = getdim(obj)
if length(rhs) != dim
throw(DimensionMismatch("dimension of right-hand side does not match basiclu object"))
end
if dim != 0
retcode = ccall((:basiclu_obj_solve_dense, libbasiclu), Int64,
(Ptr{basiclu_object}, Ptr{Float64}, Ptr{Float64}, Cchar),
pointer_from_objref(obj), rhs, rhs, trans)
if retcode == BASICLU_ERROR_INVALID_CALL
throw(ErrorException("no factorization available"))
end
checkretcode("basiclu_obj_solve_dense", retcode)
end
rhs
end
"""
solve(obj::basiclu_object, rhs::SparseVector{Float64, Int64}, trans::AbstractChar) -> SparseVector{Float64, Int64}
Solve linear system with sparse right-hand side. `rhs` is not modified. `trans`
must be 'T' for transposed solve or 'N' for forward solve.
"""
function solve(obj::basiclu_object, rhs::SparseVector{Float64, Int64}, trans::AbstractChar)
lhs = copy(rhs)
solve!(obj, lhs, trans)
end
"""
solve!(obj::basiclu_object, rhs::SparseVector{Float64, Int64}, trans::AbstractChar) -> SparseVector{Float64, Int64}
Solve linear system with sparse right-hand side. Solution overwrites `rhs`.
`trans` must be 'T' for transposed solve or 'N' for forward solve.
"""
function solve!(obj::basiclu_object, rhs::SparseVector{Float64, Int64}, trans::AbstractChar)
checktrans(trans)
dim = getdim(obj)
if length(rhs) != dim
throw(DimensionMismatch("dimension of right-hand side does not match basiclu object"))
end
if dim == 0
return spzeros(dim)
end
retcode = ccall((:basiclu_obj_solve_sparse, libbasiclu), Int64,
(Ptr{basiclu_object}, Int64, Ptr{Int64}, Ptr{Float64}, Cchar),
pointer_from_objref(obj), nnz(rhs), rowvals(rhs) .- 1, nonzeros(rhs), trans)
if retcode == BASICLU_ERROR_INVALID_CALL
throw(ErrorException("no factorization available"))
end
checkretcode("basiclu_obj_solve_sparse", retcode)
gathersol!(obj, rhs)
end
"""
solve_for_update(obj::basiclu_object, rhs::SparseVector{Float64, Int64}; getsol::Bool=false) -> SparseVector{Float64, Int64}
Solve forward system in preparation to update the factorization. `rhs` holds the
column to be inserted into the factorized matrix in the next call to
[`update`](@ref). When `getsol = true`, then the solution from the forward
solve with right-hand side `rhs` is returned. Otherwise only the update is
prepared.
"""
function solve_for_update(obj::basiclu_object, rhs::SparseVector{Float64, Int64}; getsol::Bool=false)
dim = getdim(obj)
if length(rhs) != dim
throw(DimensionMismatch("dimension of right-hand side does not match basiclu object"))
end
if dim == 0
return getsol ? spzeros(dim) : nothing
end
retcode = ccall((:basiclu_obj_solve_for_update, libbasiclu), Int64,
(Ptr{basiclu_object}, Int64, Ptr{Int64}, Ptr{Float64}, Cchar, Int64),
pointer_from_objref(obj), nnz(rhs), rowvals(rhs) .- 1, nonzeros(rhs), 'N', getsol ? 1 : 0)
if retcode == BASICLU_ERROR_INVALID_CALL
throw(ErrorException("no factorization available"))
end
if retcode == BASICLU_ERROR_MAXIMUM_UPDATES
throw(ErrorException("maximum number of updates reached"))
end
checkretcode("basiclu_obj_solve_for_update", retcode)
return getsol ? gathersol(obj) : nothing
end
"""
solve_for_update(obj::basiclu_object, pos::Int64; getsol::Bool=false) -> SparseVector{Float64, Int64}
Solve transposed system in preparation to update the factorization. `pos` holds
the column index of the factorized matrix to be replaced in the next call to
[`update`](@ref). When `getsol = true`, then the solution from the transposed
solve with a unit vector as right-hand side is returned. Otherwise only the
update is prepared.
"""
function solve_for_update(obj::basiclu_object, pos::Int64; getsol::Bool=false)
dim = getdim(obj)
if pos < 1 || pos > dim
throw(DimensionMismatch("column index outside dimension of basiclu object"))
end
if dim == 0
return getsol ? spzeros(dim) : nothing
end
retcode = ccall((:basiclu_obj_solve_for_update, libbasiclu), Int64,
(Ptr{basiclu_object}, Int64, Ptr{Int64}, Ptr{Float64}, Cchar, Int64),
pointer_from_objref(obj), 0, Ref{Int64}(pos-1), C_NULL, 'T', getsol ? 1 : 0)
if retcode == BASICLU_ERROR_INVALID_CALL
throw(ErrorException("no factorization available"))
end
if retcode == BASICLU_ERROR_MAXIMUM_UPDATES
throw(ErrorException("maximum number of updates reached"))
end
checkretcode("basiclu_obj_solve_for_update", retcode)
return getsol ? gathersol(obj) : nothing
end
"""
update(obj::basiclu_object, pivot::Float64) -> Float64
Update the factorization after a column modification. The column position and
the new column must have been set in previous calls to
[`solve_for_update`](@ref).
`pivot` is the pivot element corresponding to the update operation; i.e. when
column `j` of `B` is to be replaced by vector `v`, then `pivot = (B\\v)[j]`. The
absolute difference between `pivot` and a recomputed version can be obtained
with `getinfo(obj, :pivotError)`; this is also the return value. A pivot error
larger than 1e-8, say, indicates numerical instability and suggests
refactorization.
An error is thrown when the recomputed pivot element is below the absolute pivot
tolerance. In this case no update is performed and the old factorization remains
valid.
"""
function update(obj::basiclu_object, pivot::Float64)
dim = getdim(obj)
if dim != 0
retcode = ccall((:basiclu_obj_update, libbasiclu), Int64,
(Ptr{basiclu_object}, Float64),
pointer_from_objref(obj), pivot)
if retcode == BASICLU_ERROR_INVALID_CALL
throw(ErrorException("not prepared for update"))
end
checkretcode("basiclu_update", retcode)
end
getinfo(obj, :pivotError)
end
"""
maxvolume(obj::basiclu_object, A::SparseMatrixCSC{Float64, Int64}, basis::Vector{Int64}, volumetol::Float64=2.0) -> Int64
Given an initial basis such that `A[:,basis]` is square and nonsingular, make
one pass over the nonbasic columns of `A` and pivot each column into the basis
when it increases the absolute value of the determinant of the basis matrix by
more than a factor `volumetol`. On return `basis` has been updated. Return the
number of basis updates performed.
"""
function maxvolume(obj::basiclu_object, A::SparseMatrixCSC{Float64, Int64}, basis::Vector{Int64},
volumetol::Float64=2.0)
m, n = size(A)
if m != getdim(obj)
throw(DimensionMismatch("number of rows of A does not match dimension of basiclu object"))
end
if length(basis) != getdim(obj)
throw(DimensionMismatch("basis has incorrect length"))
end
isbasic = zeros(Int64, n)
isbasic[basis] .= 1
if sum(isbasic) != m
throw(ErrorException("duplicate index in basis"))
end
cbasis = basis .- 1
Ap = getcolptr(A) .- 1
Ai = rowvals(A) .- 1
Ax = nonzeros(A) # don't need a copy
p_nupdate = Ref{Int64}(0)
retcode = ccall((:basiclu_obj_maxvolume, libbasiclu), Int64,
(Ptr{basiclu_object}, Int64, Ptr{Int64}, Ptr{Int64}, Ptr{Float64},
Ptr{Int64}, Ptr{Int64}, Float64, Ptr{Int64}),
pointer_from_objref(obj), n, Ap, Ai, Ax, cbasis, isbasic, volumetol, p_nupdate)
basis[:] = cbasis .+ 1
checkretcode("basiclu_obj_maxvolume", retcode)
return p_nupdate[]
end
"""
maxvolbasis(A::SparseMatrixCSC{Float64, Int64}; lindeptol::Float64=1e-8,
volumetol::Float64=2.0, maxpass::Int64=2, verbose::Bool=true) -> Vector{Int64}, basiclu_object
Find a set of column indices for the matrix `AI = [A I]` such that `AI[:,basis]`
is square and nonsingular and the number of slack columns in the basis is
minimum (this is the row rank deficiency of `A`). Return the vector of column
indices of `AI` which form the basis matrix and a `basiclu_object` which holds a
factorization of the basis matrix.
Method: Scale the slack columns of `AI` by `lindeptol` and try to find a maximum
volume basis for this matrix by making at most `maxpass` calls to
[`maxvolume`](@ref). If `verbose` is true, then print the number of basis
updates after each call.
"""
function maxvolbasis(A::SparseMatrixCSC{Float64, Int64}; lindeptol::Float64=1e-8,
volumetol::Float64=2.0, maxpass::Int64=2,
verbose::Bool=true)
m, n = size(A)
rowmax = maximum(abs.(A), dims=2)[:]
rowmax = max.(rowmax, 1.)
AI = [A spdiagm(0 => lindeptol.*rowmax)]
basis = collect(n+1:n+m)
obj = basiclu_object(m)
for pass = 1:maxpass
nupdate = maxvolume(obj, AI, basis, volumetol)
if verbose
@printf("pass %d: %d updates\n", pass, nupdate)
end
if nupdate == 0 break; end
end
return basis, obj
end
function checktrans(trans::AbstractChar)
if !(trans == 'N' || trans == 'T')
throw(ArgumentError("trans argument must be 'N' (no transpose) or 'T' (transpose), got '$trans'"))
end
trans
end
function checkretcode(funcname::String, retcode::Int64)
if retcode == BASICLU_ERROR_OUT_OF_MEMORY
throw(OutOfMemoryError())
elseif retcode == BASICLU_WARNING_SINGULAR_MATRIX
msg = @sprintf("%s() encountered singular matrix", funcname)
throw(ErrorException(msg))
elseif retcode == BASICLU_ERROR_SINGULAR_UPDATE
msg = @sprintf("%s() encountered singular update", funcname)
throw(ErrorException(msg))
elseif retcode != BASICLU_OK
msg = @sprintf("%s() failed with return code %d", funcname, retcode)
throw(ErrorException(msg))
end
end
function gathersol(obj::basiclu_object)
dim = getdim(obj)
nzind = Vector{Int64}(undef, obj.nzlhs)
nzval = Vector{Float64}(undef, obj.nzlhs)
lhs = SparseVector{Float64, Int64}(dim, nzind, nzval)
gathersol!(obj, lhs)
end
function gathersol!(obj::basiclu_object, lhs::SparseVector{Float64, Int64})
dim = getdim(obj)
@assert length(lhs) == dim
if nnz(lhs) != obj.nzlhs
resize!(lhs.nzind, obj.nzlhs)
resize!(lhs.nzval, obj.nzlhs)
end
for i = 1:obj.nzlhs
idx = unsafe_load(obj.ilhs, i) + 1
@assert idx >= 1 && idx <= dim
lhs.nzind[i] = idx
lhs.nzval[i] = unsafe_load(obj.lhs, idx)
end
p = sortperm(lhs.nzind)
lhs.nzind[:] = lhs.nzind[p]
lhs.nzval[:] = lhs.nzval[p]
lhs
end
end