1 // Copyright (c) 2018 ERGO-Code. See license.txt for license.
2
3 #include "basiclu_wrapper.h"
4 #include <cassert>
5 #include <cmath>
6 #include <stdexcept>
7 #include "basiclu.h"
8
9 namespace ipx {
10
BasicLu(const Control & control,Int dim)11 BasicLu::BasicLu(const Control& control, Int dim) : control_(control) {
12 static_assert(sizeof(Int) == sizeof(lu_int),
13 "IPX integer type does not match BASICLU integer type");
14 istore_.resize(BASICLU_SIZE_ISTORE_1 + BASICLU_SIZE_ISTORE_M * dim);
15 xstore_.resize(BASICLU_SIZE_ISTORE_1 + BASICLU_SIZE_ISTORE_M * dim);
16
17 Int status = basiclu_initialize(dim, istore_.data(), xstore_.data());
18 if (status != BASICLU_OK)
19 throw std::logic_error("basiclu_initialize failed");
20
21 // Set initial size of BASICLU work arrays to 1 element, so that data()
22 // does not return NULL (not sure if it would if size is 0).
23 Li_.resize(1);
24 Lx_.resize(1);
25 Ui_.resize(1);
26 Ux_.resize(1);
27 Wi_.resize(1);
28 Wx_.resize(1);
29 xstore_[BASICLU_MEMORYL] = 1;
30 xstore_[BASICLU_MEMORYU] = 1;
31 xstore_[BASICLU_MEMORYW] = 1;
32 }
33
_Factorize(const Int * Bbegin,const Int * Bend,const Int * Bi,const double * Bx,bool strict_abs_pivottol)34 Int BasicLu::_Factorize(const Int* Bbegin, const Int* Bend, const Int* Bi,
35 const double* Bx, bool strict_abs_pivottol) {
36 Int status;
37 if (strict_abs_pivottol) {
38 xstore_[BASICLU_REMOVE_COLUMNS] = 1;
39 xstore_[BASICLU_ABS_PIVOT_TOLERANCE] = kLuDependencyTol;
40 } else {
41 xstore_[BASICLU_REMOVE_COLUMNS] = 0;
42 xstore_[BASICLU_ABS_PIVOT_TOLERANCE] = 1e-14; // BASICLU default
43 }
44 for (Int ncall = 0; ; ncall++) {
45 status = basiclu_factorize(istore_.data(), xstore_.data(),
46 Li_.data(), Lx_.data(),
47 Ui_.data(), Ux_.data(),
48 Wi_.data(), Wx_.data(),
49 Bbegin, Bend, Bi, Bx, ncall);
50 if (status != BASICLU_REALLOCATE)
51 break;
52 Reallocate();
53 }
54 if (status != BASICLU_OK && status != BASICLU_WARNING_singular_matrix)
55 throw std::logic_error("basiclu_factorize failed");
56
57 Int matrix_nz = xstore_[BASICLU_MATRIX_NZ];
58 Int lnz = xstore_[BASICLU_LNZ];
59 Int unz = xstore_[BASICLU_UNZ];
60 Int dim = xstore_[BASICLU_DIM];
61 fill_factor_ = 1.0 * (lnz+unz+dim) / matrix_nz;
62
63 double normLinv = xstore_[BASICLU_NORMEST_LINV];
64 double normUinv = xstore_[BASICLU_NORMEST_UINV];
65 double stability = xstore_[BASICLU_RESIDUAL_TEST];
66 control_.Debug(3)
67 << " normLinv = " << sci2(normLinv) << ','
68 << " normUinv = " << sci2(normUinv) << ','
69 << " stability = " << sci2(stability) << '\n';
70
71 Int ret = 0;
72 if (stability > kLuStabilityThreshold)
73 ret |= 1;
74 if (status == BASICLU_WARNING_singular_matrix)
75 ret |= 2;
76 return ret;
77 }
78
_GetFactors(SparseMatrix * L,SparseMatrix * U,Int * rowperm,Int * colperm,std::vector<Int> * dependent_cols)79 void BasicLu::_GetFactors(SparseMatrix* L, SparseMatrix* U, Int* rowperm,
80 Int* colperm, std::vector<Int>* dependent_cols) {
81 Int *Lbegin = nullptr, *Ubegin = nullptr;
82 Int *Lindex = nullptr, *Uindex = nullptr;
83 double *Lvalue = nullptr, *Uvalue = nullptr;
84 Int dim = xstore_[BASICLU_DIM];
85
86 if (L) {
87 Int lnz = xstore_[BASICLU_LNZ];
88 L->resize(dim, dim, dim+lnz);
89 Lbegin = L->colptr();
90 Lindex = L->rowidx();
91 Lvalue = L->values();
92 }
93 if (U) {
94 Int unz = xstore_[BASICLU_UNZ];
95 U->resize(dim, dim, dim+unz);
96 Ubegin = U->colptr();
97 Uindex = U->rowidx();
98 Uvalue = U->values();
99 }
100 Int status = basiclu_get_factors(istore_.data(), xstore_.data(),
101 Li_.data(), Lx_.data(),
102 Ui_.data(), Ux_.data(),
103 Wi_.data(), Wx_.data(),
104 rowperm, colperm,
105 Lbegin, Lindex, Lvalue,
106 Ubegin, Uindex, Uvalue);
107 if (status != BASICLU_OK)
108 throw std::logic_error("basiclu_get_factors failed");
109
110 if (L) {
111 // Remove unit diagonal from L.
112 Int num_dropped = RemoveDiagonal(*L, nullptr);
113 assert(num_dropped == dim);
114 (void)(num_dropped);
115 }
116 if (dependent_cols) {
117 // Dependent columns are at the end of the BASICLU pivot sequence.
118 Int rank = xstore_[BASICLU_RANK];
119 dependent_cols->clear();
120 for (Int k = rank; k < dim; k++)
121 dependent_cols->push_back(k);
122 }
123 }
124
_SolveDense(const Vector & rhs,Vector & lhs,char trans)125 void BasicLu::_SolveDense(const Vector& rhs, Vector& lhs, char trans) {
126 Int status = basiclu_solve_dense(istore_.data(), xstore_.data(),
127 Li_.data(), Lx_.data(),
128 Ui_.data(), Ux_.data(),
129 Wi_.data(), Wx_.data(),
130 &rhs[0], &lhs[0], trans);
131 if (status != BASICLU_OK)
132 throw std::logic_error("basiclu_solve_dense failed");
133 }
134
_FtranForUpdate(Int nzrhs,const Int * bi,const double * bx)135 void BasicLu::_FtranForUpdate(Int nzrhs, const Int* bi, const double* bx) {
136 Int status;
137 for (Int ncall = 0; ; ncall++) {
138 status = basiclu_solve_for_update(istore_.data(), xstore_.data(),
139 Li_.data(), Lx_.data(),
140 Ui_.data(), Ux_.data(),
141 Wi_.data(), Wx_.data(),
142 nzrhs, bi, bx,
143 nullptr, nullptr, nullptr, 'N');
144 if (status != BASICLU_REALLOCATE)
145 break;
146 Reallocate();
147 }
148 if (status != BASICLU_OK)
149 throw std::logic_error(
150 "basiclu_solve_for_update (ftran without lhs) failed");
151 }
152
_FtranForUpdate(Int nzrhs,const Int * bi,const double * bx,IndexedVector & lhs)153 void BasicLu::_FtranForUpdate(Int nzrhs, const Int* bi, const double* bx,
154 IndexedVector& lhs) {
155 Int status;
156 Int nzlhs = 0;
157 lhs.set_to_zero();
158 for (Int ncall = 0; ; ncall++) {
159 status = basiclu_solve_for_update(istore_.data(), xstore_.data(),
160 Li_.data(), Lx_.data(),
161 Ui_.data(), Ux_.data(),
162 Wi_.data(), Wx_.data(),
163 nzrhs, bi, bx,
164 &nzlhs, lhs.pattern(), lhs.elements(),
165 'N');
166 if (status != BASICLU_REALLOCATE)
167 break;
168 Reallocate();
169 }
170 if (status != BASICLU_OK)
171 throw std::logic_error(
172 "basiclu_solve_for_update (ftran with lhs) failed");
173 lhs.set_nnz(nzlhs);
174 }
175
_BtranForUpdate(Int j)176 void BasicLu::_BtranForUpdate(Int j) {
177 Int status;
178 for (Int ncall = 0; ; ncall++) {
179 status = basiclu_solve_for_update(istore_.data(), xstore_.data(),
180 Li_.data(), Lx_.data(),
181 Ui_.data(), Ux_.data(),
182 Wi_.data(), Wx_.data(),
183 0, &j, nullptr,
184 nullptr, nullptr, nullptr, 'T');
185 if (status != BASICLU_REALLOCATE)
186 break;
187 Reallocate();
188 }
189 if (status != BASICLU_OK)
190 throw std::logic_error(
191 "basiclu_solve_for_update (btran without lhs) failed");
192 }
193
_BtranForUpdate(Int j,IndexedVector & lhs)194 void BasicLu::_BtranForUpdate(Int j, IndexedVector& lhs) {
195 Int status;
196 Int nzlhs = 0;
197 lhs.set_to_zero();
198 for (Int ncall = 0; ; ncall++) {
199 status = basiclu_solve_for_update(istore_.data(), xstore_.data(),
200 Li_.data(), Lx_.data(),
201 Ui_.data(), Ux_.data(),
202 Wi_.data(), Wx_.data(),
203 0, &j, nullptr,
204 &nzlhs, lhs.pattern(), lhs.elements(),
205 'T');
206 if (status != BASICLU_REALLOCATE)
207 break;
208 Reallocate();
209 }
210 if (status != BASICLU_OK)
211 throw std::logic_error(
212 "basiclu_solve_for_update (btran with lhs) failed");
213 lhs.set_nnz(nzlhs);
214 }
215
_Update(double pivot)216 Int BasicLu::_Update(double pivot) {
217 double max_eta_old = xstore_[BASICLU_MAX_ETA];
218 Int status;
219 for (Int ncall = 0; ; ncall++) {
220 status = basiclu_update(istore_.data(), xstore_.data(),
221 Li_.data(), Lx_.data(),
222 Ui_.data(), Ux_.data(),
223 Wi_.data(), Wx_.data(), pivot);
224 if (status != BASICLU_REALLOCATE)
225 break;
226 Reallocate();
227 }
228 if (status != BASICLU_OK && status != BASICLU_ERROR_singular_update)
229 throw std::logic_error("basiclu_update failed");
230 if (status == BASICLU_ERROR_singular_update)
231 return -1;
232
233 // Print a debugging message if a new eta entry is large.
234 double max_eta = xstore_[BASICLU_MAX_ETA];
235 if (max_eta > 1e10 && max_eta > max_eta_old)
236 control_.Debug(3) << " max eta = " << sci2(max_eta) << '\n';
237
238 // stability check
239 double pivot_error = xstore_[BASICLU_PIVOT_ERROR];
240 if (pivot_error > kFtDiagErrorTol) {
241 control_.Debug(3)
242 << " relative error in new diagonal entry of U = "
243 << sci2(pivot_error) << '\n';
244 return 1;
245 }
246 return 0;
247 }
248
_NeedFreshFactorization()249 bool BasicLu::_NeedFreshFactorization() {
250 Int dim = xstore_[BASICLU_DIM];
251 Int nforrest = xstore_[BASICLU_NFORREST];
252 double update_cost = xstore_[BASICLU_UPDATE_COST];
253
254 return nforrest == dim || update_cost > 1.0;
255 }
256
_fill_factor() const257 double BasicLu::_fill_factor() const {
258 return fill_factor_;
259 }
260
_pivottol() const261 double BasicLu::_pivottol() const {
262 return xstore_[BASICLU_REL_PIVOT_TOLERANCE];
263 }
264
_pivottol(double new_pivottol)265 void BasicLu::_pivottol(double new_pivottol) {
266 xstore_[BASICLU_REL_PIVOT_TOLERANCE] = new_pivottol;
267 }
268
Reallocate()269 void BasicLu::Reallocate() {
270 assert(Li_.size() == xstore_[BASICLU_MEMORYL]);
271 assert(Lx_.size() == xstore_[BASICLU_MEMORYL]);
272 assert(Ui_.size() == xstore_[BASICLU_MEMORYU]);
273 assert(Ux_.size() == xstore_[BASICLU_MEMORYU]);
274 assert(Wi_.size() == xstore_[BASICLU_MEMORYW]);
275 assert(Wx_.size() == xstore_[BASICLU_MEMORYW]);
276
277 if (xstore_[BASICLU_ADD_MEMORYL] > 0) {
278 Int new_size = xstore_[BASICLU_MEMORYL] + xstore_[BASICLU_ADD_MEMORYL];
279 new_size *= kReallocFactor;
280 Li_.resize(new_size);
281 Lx_.resize(new_size);
282 xstore_[BASICLU_MEMORYL] = new_size;
283 }
284 if (xstore_[BASICLU_ADD_MEMORYU] > 0) {
285 Int new_size = xstore_[BASICLU_MEMORYU] + xstore_[BASICLU_ADD_MEMORYU];
286 new_size *= kReallocFactor;
287 Ui_.resize(new_size);
288 Ux_.resize(new_size);
289 xstore_[BASICLU_MEMORYU] = new_size;
290 }
291 if (xstore_[BASICLU_ADD_MEMORYW] > 0) {
292 Int new_size = xstore_[BASICLU_MEMORYW] + xstore_[BASICLU_ADD_MEMORYW];
293 new_size *= kReallocFactor;
294 Wi_.resize(new_size);
295 Wx_.resize(new_size);
296 xstore_[BASICLU_MEMORYW] = new_size;
297 }
298 }
299
300 } // namespace ipx
301