1 /*
2 * lu_internal.c
3 *
4 * Copyright (C) 2016-2018 ERGO-Code
5 *
6 * Functions to load/save/reset struct lu objects
7 *
8 */
9
10 #include "lu_internal.h"
11
12 /* private entries in xstore */
13 #define BASICLU_TASK 256
14 #define BASICLU_FTCOLUMN_IN 257
15 #define BASICLU_FTCOLUMN_OUT 258
16 #define BASICLU_PIVOT_ROW 259
17 #define BASICLU_PIVOT_COL 260
18 #define BASICLU_RANKDEF 261
19 #define BASICLU_MIN_COLNZ 262
20 #define BASICLU_MIN_ROWNZ 263
21 #define BASICLU_MARKER 266
22 #define BASICLU_UPDATE_COST_NUMER 267
23 #define BASICLU_UPDATE_COST_DENOM 268
24 #define BASICLU_PIVOTLEN 269
25
26
27 /* ==========================================================================
28 * lu_load
29 *
30 * Initialize @this from @istore, @xstore if these are a valid BASICLU
31 * instance. The remaining arguments are copied only and can be NULL.
32 *
33 * Return BASICLU_OK or BASICLU_ERROR_invalid_store
34 * ========================================================================== */
35
lu_load(struct lu * this,lu_int * istore,double * xstore,lu_int * Li,double * Lx,lu_int * Ui,double * Ux,lu_int * Wi,double * Wx)36 lu_int lu_load(
37 struct lu *this, lu_int *istore, double *xstore, lu_int *Li, double *Lx,
38 lu_int *Ui, double *Ux, lu_int *Wi, double *Wx)
39 {
40 lu_int m, *iptr;
41 double *xptr;
42
43 if (!istore || istore[0] != BASICLU_HASH ||
44 !xstore || xstore[0] != BASICLU_HASH )
45 {
46 return BASICLU_ERROR_invalid_store;
47 }
48
49 /* user parameters */
50 this->Lmem = xstore[BASICLU_MEMORYL];
51 this->Umem = xstore[BASICLU_MEMORYU];
52 this->Wmem = xstore[BASICLU_MEMORYW];
53 this->droptol = xstore[BASICLU_DROP_TOLERANCE];
54 this->abstol = xstore[BASICLU_ABS_PIVOT_TOLERANCE];
55 this->reltol = xstore[BASICLU_REL_PIVOT_TOLERANCE];
56 this->reltol = fmin(this->reltol, 1.0);
57 this->nzbias = xstore[BASICLU_BIAS_NONZEROS];
58 this->maxsearch = xstore[BASICLU_MAXN_SEARCH_PIVOT];
59 this->pad = xstore[BASICLU_PAD];
60 this->stretch = xstore[BASICLU_STRETCH];
61 this->compress_thres = xstore[BASICLU_COMPRESSION_THRESHOLD];
62 this->sparse_thres = xstore[BASICLU_SPARSE_THRESHOLD];
63 this->search_rows = xstore[BASICLU_SEARCH_ROWS] != 0;
64
65 /* user readable */
66 this->m = m = xstore[BASICLU_DIM];
67 this->addmemL = 0;
68 this->addmemU = 0;
69 this->addmemW = 0;
70
71 this->nupdate = xstore[BASICLU_NUPDATE];
72 this->nforrest = xstore[BASICLU_NFORREST];
73 this->nfactorize = xstore[BASICLU_NFACTORIZE];
74 this->nupdate_total = xstore[BASICLU_NUPDATE_TOTAL];
75 this->nforrest_total = xstore[BASICLU_NFORREST_TOTAL];
76 this->nsymperm_total = xstore[BASICLU_NSYMPERM_TOTAL];
77 this->Lnz = xstore[BASICLU_LNZ];
78 this->Unz = xstore[BASICLU_UNZ];
79 this->Rnz = xstore[BASICLU_RNZ];
80 this->min_pivot = xstore[BASICLU_MIN_PIVOT];
81 this->max_pivot = xstore[BASICLU_MAX_PIVOT];
82 this->max_eta = xstore[BASICLU_MAX_ETA];
83 this->update_cost_numer = xstore[BASICLU_UPDATE_COST_NUMER];
84 this->update_cost_denom = xstore[BASICLU_UPDATE_COST_DENOM];
85 this->time_factorize = xstore[BASICLU_TIME_FACTORIZE];
86 this->time_solve = xstore[BASICLU_TIME_SOLVE];
87 this->time_update = xstore[BASICLU_TIME_UPDATE];
88 this->time_factorize_total = xstore[BASICLU_TIME_FACTORIZE_TOTAL];
89 this->time_solve_total = xstore[BASICLU_TIME_SOLVE_TOTAL];
90 this->time_update_total = xstore[BASICLU_TIME_UPDATE_TOTAL];
91 this->Lflops = xstore[BASICLU_LFLOPS];
92 this->Uflops = xstore[BASICLU_UFLOPS];
93 this->Rflops = xstore[BASICLU_RFLOPS];
94 this->condestL = xstore[BASICLU_CONDEST_L];
95 this->condestU = xstore[BASICLU_CONDEST_U];
96 this->normL = xstore[BASICLU_NORM_L];
97 this->normU = xstore[BASICLU_NORM_U];
98 this->normestLinv = xstore[BASICLU_NORMEST_LINV];
99 this->normestUinv = xstore[BASICLU_NORMEST_UINV];
100 this->onenorm = xstore[BASICLU_MATRIX_ONENORM];
101 this->infnorm = xstore[BASICLU_MATRIX_INFNORM];
102 this->residual_test = xstore[BASICLU_RESIDUAL_TEST];
103
104 this->matrix_nz = xstore[BASICLU_MATRIX_NZ];
105 this->rank = xstore[BASICLU_RANK];
106 this->bump_size = xstore[BASICLU_BUMP_SIZE];
107 this->bump_nz = xstore[BASICLU_BUMP_NZ];
108 this->nsearch_pivot = xstore[BASICLU_NSEARCH_PIVOT];
109 this->nexpand = xstore[BASICLU_NEXPAND];
110 this->ngarbage = xstore[BASICLU_NGARBAGE];
111 this->factor_flops = xstore[BASICLU_FACTOR_FLOPS];
112 this->time_singletons = xstore[BASICLU_TIME_SINGLETONS];
113 this->time_search_pivot = xstore[BASICLU_TIME_SEARCH_PIVOT];
114 this->time_elim_pivot = xstore[BASICLU_TIME_ELIM_PIVOT];
115
116 this->pivot_error = xstore[BASICLU_PIVOT_ERROR];
117
118 /* private */
119 this->task = xstore[BASICLU_TASK];
120 this->pivot_row = xstore[BASICLU_PIVOT_ROW];
121 this->pivot_col = xstore[BASICLU_PIVOT_COL];
122 this->ftran_for_update = xstore[BASICLU_FTCOLUMN_IN];
123 this->btran_for_update = xstore[BASICLU_FTCOLUMN_OUT];
124 this->marker = xstore[BASICLU_MARKER];
125 this->pivotlen = xstore[BASICLU_PIVOTLEN];
126 this->rankdef = xstore[BASICLU_RANKDEF];
127 this->min_colnz = xstore[BASICLU_MIN_COLNZ];
128 this->min_rownz = xstore[BASICLU_MIN_ROWNZ];
129
130 /* aliases to user arrays */
131 this->Lindex = Li; this->Lvalue = Lx;
132 this->Uindex = Ui; this->Uvalue = Ux;
133 this->Windex = Wi; this->Wvalue = Wx;
134
135 /* partition istore for factorize */
136 iptr = istore + 1;
137 this->colcount_flink = iptr; iptr += 2*m+2;
138 this->colcount_blink = iptr; iptr += 2*m+2;
139 this->rowcount_flink = iptr; iptr += 2*m+2;
140 this->rowcount_blink = iptr; iptr += 2*m+2;
141 this->Wbegin = iptr; iptr += 2*m+1;
142 this->Wend = iptr; iptr += 2*m+1;
143 this->Wflink = iptr; iptr += 2*m+1;
144 this->Wblink = iptr; iptr += 2*m+1;
145 this->pinv = iptr; iptr += m;
146 this->qinv = iptr; iptr += m;
147 this->Lbegin_p = iptr; iptr += m+1;
148 this->Ubegin = iptr; iptr += m+1;
149 this->iwork0 = iptr; iptr += m;
150
151 /* share istore memory for solve/update */
152 this->pivotcol = this->colcount_flink;
153 this->pivotrow = this->colcount_blink;
154 this->Rbegin = this->rowcount_flink;
155 this->eta_row = this->rowcount_flink + m+1;
156 this->iwork1 = this->rowcount_blink;
157 this->Lbegin = this->Wbegin + m+1;
158 this->Ltbegin = this->Wend + m+1;
159 this->Ltbegin_p = this->Wflink + m+1;
160 this->p = this->Wblink + m+1;
161 this->pmap = this->pinv;
162 this->qmap = this->qinv;
163 this->marked = this->iwork0;
164
165 /* partition xstore for factorize and update */
166 xptr = xstore + 512;
167 this->work0 = xptr; xptr += m;
168 this->work1 = xptr; xptr += m;
169 this->col_pivot = xptr; xptr += m;
170 this->row_pivot = xptr; xptr += m;
171
172 /*
173 * Reset @marked if increasing @marker by four causes overflow.
174 */
175 if (this->marker > LU_INT_MAX-4)
176 {
177 memset(this->marked, 0, m*sizeof(lu_int));
178 this->marker = 0;
179 }
180
181 /*
182 * One past the final position in @Wend must hold the file size.
183 * The file has 2*m lines while factorizing and m lines otherwise.
184 */
185 if (this->nupdate >= 0)
186 this->Wend[m] = this->Wmem;
187 else
188 this->Wend[2*m] = this->Wmem;
189
190 return BASICLU_OK;
191 }
192
193
194 /* ==========================================================================
195 * lu_save
196 *
197 * Copy scalar entries (except for user parameters) from @this to @istore,
198 * @xstore. Store status code.
199 *
200 * Return @status
201 * ========================================================================== */
202
lu_save(const struct lu * this,lu_int * istore,double * xstore,lu_int status)203 lu_int lu_save(
204 const struct lu *this, lu_int *istore, double *xstore, lu_int status)
205 {
206 /* user readable */
207 xstore[BASICLU_STATUS] = status;
208 xstore[BASICLU_ADD_MEMORYL] = this->addmemL;
209 xstore[BASICLU_ADD_MEMORYU] = this->addmemU;
210 xstore[BASICLU_ADD_MEMORYW] = this->addmemW;
211
212 xstore[BASICLU_NUPDATE] = this->nupdate;
213 xstore[BASICLU_NFORREST] = this->nforrest;
214 xstore[BASICLU_NFACTORIZE] = this->nfactorize;
215 xstore[BASICLU_NUPDATE_TOTAL] = this->nupdate_total;
216 xstore[BASICLU_NFORREST_TOTAL] = this->nforrest_total;
217 xstore[BASICLU_NSYMPERM_TOTAL] = this->nsymperm_total;
218 xstore[BASICLU_LNZ] = this->Lnz;
219 xstore[BASICLU_UNZ] = this->Unz;
220 xstore[BASICLU_RNZ] = this->Rnz;
221 xstore[BASICLU_MIN_PIVOT] = this->min_pivot;
222 xstore[BASICLU_MAX_PIVOT] = this->max_pivot;
223 xstore[BASICLU_MAX_ETA] = this->max_eta;
224 xstore[BASICLU_UPDATE_COST_NUMER] = this->update_cost_numer;
225 xstore[BASICLU_UPDATE_COST_DENOM] = this->update_cost_denom;
226 xstore[BASICLU_UPDATE_COST] =
227 this->update_cost_numer / this->update_cost_denom;
228 xstore[BASICLU_TIME_FACTORIZE] = this->time_factorize;
229 xstore[BASICLU_TIME_SOLVE] = this->time_solve;
230 xstore[BASICLU_TIME_UPDATE] = this->time_update;
231 xstore[BASICLU_TIME_FACTORIZE_TOTAL] = this->time_factorize_total;
232 xstore[BASICLU_TIME_SOLVE_TOTAL] = this->time_solve_total;
233 xstore[BASICLU_TIME_UPDATE_TOTAL] = this->time_update_total;
234 xstore[BASICLU_LFLOPS] = this->Lflops;
235 xstore[BASICLU_UFLOPS] = this->Uflops;
236 xstore[BASICLU_RFLOPS] = this->Rflops;
237 xstore[BASICLU_CONDEST_L] = this->condestL;
238 xstore[BASICLU_CONDEST_U] = this->condestU;
239 xstore[BASICLU_NORM_L] = this->normL;
240 xstore[BASICLU_NORM_U] = this->normU;
241 xstore[BASICLU_NORMEST_LINV] = this->normestLinv;
242 xstore[BASICLU_NORMEST_UINV] = this->normestUinv;
243 xstore[BASICLU_MATRIX_ONENORM] = this->onenorm;
244 xstore[BASICLU_MATRIX_INFNORM] = this->infnorm;
245 xstore[BASICLU_RESIDUAL_TEST] = this->residual_test;
246
247 xstore[BASICLU_MATRIX_NZ] = this->matrix_nz;
248 xstore[BASICLU_RANK] = this->rank;
249 xstore[BASICLU_BUMP_SIZE] = this->bump_size;
250 xstore[BASICLU_BUMP_NZ] = this->bump_nz;
251 xstore[BASICLU_NSEARCH_PIVOT] = this->nsearch_pivot;
252 xstore[BASICLU_NEXPAND] = this->nexpand;
253 xstore[BASICLU_NGARBAGE] = this->ngarbage;
254 xstore[BASICLU_FACTOR_FLOPS] = this->factor_flops;
255 xstore[BASICLU_TIME_SINGLETONS] = this->time_singletons;
256 xstore[BASICLU_TIME_SEARCH_PIVOT] = this->time_search_pivot;
257 xstore[BASICLU_TIME_ELIM_PIVOT] = this->time_elim_pivot;
258
259 xstore[BASICLU_PIVOT_ERROR] = this->pivot_error;
260
261 /* private */
262 xstore[BASICLU_TASK] = this->task;
263 xstore[BASICLU_PIVOT_ROW] = this->pivot_row;
264 xstore[BASICLU_PIVOT_COL] = this->pivot_col;
265 xstore[BASICLU_FTCOLUMN_IN] = this->ftran_for_update;
266 xstore[BASICLU_FTCOLUMN_OUT] = this->btran_for_update;
267 xstore[BASICLU_MARKER] = this->marker;
268 xstore[BASICLU_PIVOTLEN] = this->pivotlen;
269 xstore[BASICLU_RANKDEF] = this->rankdef;
270 xstore[BASICLU_MIN_COLNZ] = this->min_colnz;
271 xstore[BASICLU_MIN_ROWNZ] = this->min_rownz;
272
273 return status;
274 }
275
276
277 /* ==========================================================================
278 * lu_reset
279 *
280 * Reset @this for a new factorization. Invalidate current factorization.
281 * ========================================================================== */
282
lu_reset(struct lu * this)283 void lu_reset(struct lu *this)
284 {
285 /* user readable */
286 this->nupdate = -1; /* invalidate factorization */
287 this->nforrest = 0;
288 this->Lnz = 0;
289 this->Unz = 0;
290 this->Rnz = 0;
291 this->min_pivot = 0;
292 this->max_pivot = 0;
293 this->max_eta = 0;
294 this->update_cost_numer = 0;
295 this->update_cost_denom = 1;
296 this->time_factorize = 0;
297 this->time_solve = 0;
298 this->time_update = 0;
299 this->Lflops = 0;
300 this->Uflops = 0;
301 this->Rflops = 0;
302 this->condestL = 0;
303 this->condestU = 0;
304 this->normL = 0;
305 this->normU = 0;
306 this->normestLinv = 0;
307 this->normestUinv = 0;
308 this->onenorm = 0;
309 this->infnorm = 0;
310 this->residual_test = 0;
311
312 this->matrix_nz = 0;
313 this->rank = 0;
314 this->bump_size = 0;
315 this->bump_nz = 0;
316 this->nsearch_pivot = 0;
317 this->nexpand = 0;
318 this->ngarbage = 0;
319 this->factor_flops = 0;
320 this->time_singletons = 0;
321 this->time_search_pivot = 0;
322 this->time_elim_pivot = 0;
323
324 this->pivot_error = 0;
325
326 /* private */
327 this->task = NO_TASK;
328 this->pivot_row = -1;
329 this->pivot_col = -1;
330 this->ftran_for_update = -1;
331 this->btran_for_update = -1;
332 this->marker = 0;
333 this->pivotlen = 0;
334 this->rankdef = 0;
335 this->min_colnz = 1;
336 this->min_rownz = 1;
337
338 /*
339 * One past the final position in @Wend must hold the file size.
340 * The file has 2*m lines during factorization.
341 */
342 this->Wend[2*this->m] = this->Wmem;
343
344 /*
345 * The integer workspace iwork0 must be zeroed for a new factorization.
346 * The double workspace work0 actually needs only be zeroed once in the
347 * initialization of xstore. However, it is easier and more consistent
348 * to do that here as well.
349 */
350 memset(this->iwork0, 0, this->m * sizeof(lu_int));
351 memset(this->work0, 0, this->m * sizeof(lu_int));
352 }
353