1 // Copyright (c) 2018 ERGO-Code. See license.txt for license.
2
3 #include "basiclu_kernel.h"
4 #include <cassert>
5 #include <new> // bad_alloc
6 #include <stdexcept>
7 #include "basiclu.h"
8
9 namespace ipx {
10
11 // Wraps the BASICLU object into a struct with contructor/destructor for RAII.
12 struct BasicLuHelper {
13 static_assert(sizeof(Int) == sizeof(lu_int),
14 "IPX integer type does not match BASICLU integer type");
15
BasicLuHelperipx::BasicLuHelper16 BasicLuHelper(Int dim) {
17 Int err = basiclu_obj_initialize(&obj, dim);
18 if (err == BASICLU_ERROR_out_of_memory)
19 throw std::bad_alloc();
20 if (err != BASICLU_OK)
21 throw std::logic_error("basiclu_obj_initialize failed");
22 }
23
~BasicLuHelperipx::BasicLuHelper24 ~BasicLuHelper() {
25 basiclu_obj_free(&obj);
26 }
27
28 struct basiclu_object obj;
29 };
30
_Factorize(Int dim,const Int * Bbegin,const Int * Bend,const Int * Bi,const double * Bx,double pivottol,bool strict_abs_pivottol,SparseMatrix * L,SparseMatrix * U,std::vector<Int> * rowperm,std::vector<Int> * colperm,std::vector<Int> * dependent_cols)31 void BasicLuKernel::_Factorize(Int dim, const Int* Bbegin, const Int* Bend,
32 const Int* Bi, const double* Bx, double pivottol,
33 bool strict_abs_pivottol,
34 SparseMatrix* L, SparseMatrix* U,
35 std::vector<Int>* rowperm,
36 std::vector<Int>* colperm,
37 std::vector<Int>* dependent_cols) {
38 BasicLuHelper lu(dim);
39 lu.obj.xstore[BASICLU_REL_PIVOT_TOLERANCE] = pivottol;
40 if (strict_abs_pivottol) {
41 lu.obj.xstore[BASICLU_ABS_PIVOT_TOLERANCE] = kLuDependencyTol;
42 lu.obj.xstore[BASICLU_REMOVE_COLUMNS] = 1;
43 }
44 Int err = basiclu_obj_factorize(&lu.obj, Bbegin, Bend, Bi, Bx);
45 if (err == BASICLU_ERROR_out_of_memory)
46 throw std::bad_alloc();
47 if (err != BASICLU_OK && err != BASICLU_WARNING_singular_matrix)
48 throw std::logic_error("basiclu_obj_factorize failed");
49
50 // Extract factors. Dependent columns are at the end of the BASICLU pivot
51 // sequence.
52 Int rank = lu.obj.xstore[BASICLU_RANK];
53 dependent_cols->clear();
54 for (Int k = rank; k < dim; k++)
55 dependent_cols->push_back(k);
56 L->resize(dim, dim, dim + lu.obj.xstore[BASICLU_LNZ]);
57 U->resize(dim, dim, dim + lu.obj.xstore[BASICLU_UNZ]);
58 rowperm->resize(dim);
59 colperm->resize(dim);
60 err = basiclu_obj_get_factors(&lu.obj, rowperm->data(), colperm->data(),
61 L->colptr(), L->rowidx(), L->values(),
62 U->colptr(), U->rowidx(), U->values());
63 if (err != BASICLU_OK)
64 throw std::logic_error("basiclu_obj_get_factors failed");
65
66 // Remove unit diagonal from L.
67 Int num_dropped = RemoveDiagonal(*L, nullptr);
68 assert(num_dropped == dim);
69 assert(L->entries() == lu.obj.xstore[BASICLU_LNZ]);
70 }
71
72 } // namespace ipx
73