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