1 /*
2 * lu_singletons.c
3 *
4 * Copyright (C) 2016-2018 ERGO-Code
5 *
6 * Build initial triangular factors
7 *
8 */
9
10 #include "lu_internal.h"
11
12 static lu_int singleton_cols
13 (
14 const lu_int m,
15 const lu_int *Bbegin, /* B columnwise */
16 const lu_int *Bend,
17 const lu_int *Bi,
18 const double *Bx,
19 const lu_int *Btp, /* B rowwise */
20 const lu_int *Bti,
21 const double *Btx,
22 lu_int *Up,
23 lu_int *Ui,
24 double *Ux,
25 lu_int *Lp,
26 lu_int *Li,
27 double *Lx,
28 double *col_pivot,
29 lu_int *pinv,
30 lu_int *qinv,
31 lu_int *iset, /* size m workspace */
32 lu_int *queue, /* size m workspace */
33 lu_int rank,
34 double abstol
35 );
36
37 static lu_int singleton_rows
38 (
39 const lu_int m,
40 const lu_int *Bbegin, /* B columnwise */
41 const lu_int *Bend,
42 const lu_int *Bi,
43 const double *Bx,
44 const lu_int *Btp, /* B rowwise */
45 const lu_int *Bti,
46 const double *Btx,
47 lu_int *Up,
48 lu_int *Ui,
49 double *Ux,
50 lu_int *Lp,
51 lu_int *Li,
52 double *Lx,
53 double *col_pivot,
54 lu_int *pinv,
55 lu_int *qinv,
56 lu_int *iset, /* size m workspace */
57 lu_int *queue, /* size m workspace */
58 lu_int rank,
59 double abstol
60 );
61
62
63 /*
64 * lu_singletons()
65 *
66 * Initialize the data structures which store the LU factors during
67 * factorization and eliminate pivots with Markowitz cost zero.
68 *
69 * During factorization the inverse pivot sequence is recorded in pinv, qinv:
70 *
71 * pinv[i] >= 0 if row i was pivot row in stage pinv[i]
72 * pinv[i] == -1 if row i has not been pivot row yet
73 * qinv[j] >= 0 if col j was pivot col in stage qinv[j]
74 * qinv[j] == -1 if col j has not been pivot col yet
75 *
76 * The lower triangular factor is composed columnwise in Lindex, Lvalue.
77 * The upper triangular factor is composed rowwise in Uindex, Uvalue.
78 * After rank steps of factorization:
79 *
80 * Lbegin_p[rank] is the next unused position in Lindex, Lvalue.
81 *
82 * Lindex[Lbegin_p[k]..], Lvalue[Lbegin_p[k]..] for 0 <= k < rank
83 * stores the column of L computed in stage k without the unit diagonal.
84 * The column is terminated by a negative index.
85 *
86 * Ubegin[rank] is the next unused position in Uindex, Uvalue.
87 *
88 * Uindex[Ubegin[k]..Ubegin[k+1]-1], Uvalue[Ubegin[k]..Ubegin[k+1]-1]
89 * stores the row of U computed in stage k without the pivot element.
90 *
91 * lu_singletons() does rank >= 0 steps of factorization until no singletons are
92 * left. We can either eliminate singleton columns before singleton rows or vice
93 * versa. When nzbias >= 0, then eliminate singleton columns first to keep L
94 * sparse. Otherwise eliminate singleton rows first. The resulting permutations
95 * P, Q (stored in inverse form) make PBQ' of the form
96 *
97 * \uuuuuuuuuuuuuuuuuuuuuuu
98 * \u u
99 * \u u
100 * \u u
101 * \u u
102 * PBQ' = \uuuuuuu__________u singleton columns before
103 * \ | | singleton rows
104 * l\ | |
105 * ll\ | |
106 * l l\ | BUMP |
107 * l l\ | |
108 * lllll\|__________|
109 *
110 * \
111 * l\
112 * ll\
113 * l l\
114 * l l\
115 * l l\ __________
116 * PBQ' = l l\uuuuu| | singleton rows before
117 * l l \u u| | singleton columns
118 * l l \u u| |
119 * l l \uu| BUMP |
120 * l l \u| |
121 * llllll \|__________|
122 *
123 * Off-diagonals from singleton columns (u) are stored in U, off-diagonals from
124 * singleton rows (l) are stored in L and divided by the diagonal. Diagonals (\)
125 * are stored in col_pivot.
126 *
127 * Do not pivot on elements which are zero or less than abstol in magnitude.
128 * When such pivots occur, the row/column remains in the active submatrix and
129 * the bump factorization will detect the singularity.
130 *
131 * Return:
132 *
133 * BASICLU_REALLOCATE less than nnz(B) memory in L, U or W
134 * BASICLU_ERROR_invalid_argument matrix B is invalid (negative number of
135 * entries in column, index out of range,
136 * duplicates)
137 * BASICLU_OK
138 */
139
lu_singletons(struct lu * this,const lu_int * Bbegin,const lu_int * Bend,const lu_int * Bi,const double * Bx)140 lu_int lu_singletons(
141 struct lu *this, const lu_int *Bbegin, const lu_int *Bend, const lu_int *Bi,
142 const double *Bx)
143 {
144 const lu_int m = this->m;
145 const lu_int Lmem = this->Lmem;
146 const lu_int Umem = this->Umem;
147 const lu_int Wmem = this->Wmem;
148 const double abstol = this->abstol;
149 const lu_int nzbias = this->nzbias;
150 lu_int *pinv = this->pinv;
151 lu_int *qinv = this->qinv;
152 lu_int *Lbegin_p = this->Lbegin_p;
153 lu_int *Ubegin = this->Ubegin;
154 double *col_pivot = this->col_pivot;
155 lu_int *Lindex = this->Lindex;
156 double *Lvalue = this->Lvalue;
157 lu_int *Uindex = this->Uindex;
158 double *Uvalue = this->Uvalue;
159 lu_int *iwork1 = this->iwork1;
160 lu_int *iwork2 = iwork1 + m;
161
162 lu_int *Btp = this->Wbegin; /* build B rowwise in W */
163 lu_int *Bti = this->Windex;
164 double *Btx = this->Wvalue;
165
166 lu_int i, j, pos, put, rank, Bnz, ok;
167 double tic[2];
168
169 /* -------------------------------- */
170 /* Check matrix and build transpose */
171 /* -------------------------------- */
172
173 /* Check pointers and count nnz(B). */
174 Bnz = 0;
175 ok = 1;
176 for (j = 0; j < m && ok; j++)
177 {
178 if (Bend[j] < Bbegin[j])
179 ok = 0;
180 else
181 Bnz += Bend[j] - Bbegin[j];
182 }
183 if (!ok)
184 return BASICLU_ERROR_invalid_argument;
185
186 /* Check if sufficient memory in L, U, W. */
187 ok = 1;
188 if (Lmem < Bnz) { this->addmemL = Bnz-Lmem; ok = 0; }
189 if (Umem < Bnz) { this->addmemU = Bnz-Umem; ok = 0; }
190 if (Wmem < Bnz) { this->addmemW = Bnz-Wmem; ok = 0; }
191 if (!ok)
192 return BASICLU_REALLOCATE;
193
194 /* Count nz per row, check indices. */
195 memset(iwork1, 0, m*sizeof(lu_int)); /* row counts */
196 ok = 1;
197 for (j = 0; j < m && ok; j++)
198 {
199 for (pos = Bbegin[j]; pos < Bend[j] && ok; pos++)
200 {
201 i = Bi[pos];
202 if (i < 0 || i >= m)
203 ok = 0;
204 else
205 iwork1[i]++;
206 }
207 }
208 if (!ok)
209 return BASICLU_ERROR_invalid_argument;
210
211 /* Pack matrix rowwise, check for duplicates. */
212 put = 0;
213 for (i = 0; i < m; i++) /* set row pointers */
214 {
215 Btp[i] = put;
216 put += iwork1[i];
217 iwork1[i] = Btp[i];
218 }
219 Btp[m] = put;
220 assert(put == Bnz);
221 ok = 1;
222 for (j = 0; j < m; j++) /* fill rows */
223 {
224 for (pos = Bbegin[j]; pos < Bend[j]; pos++)
225 {
226 i = Bi[pos];
227 put = iwork1[i]++;
228 Bti[put] = j;
229 Btx[put] = Bx [pos];
230 if (put > Btp[i] && Bti[put-1] == j)
231 ok = 0;
232 }
233 }
234 if (!ok)
235 return BASICLU_ERROR_invalid_argument;
236
237 /* ---------------- */
238 /* Pivot singletons */
239 /* ---------------- */
240
241 /* No pivot rows or pivot columns so far. */
242 for (i = 0; i < m; i++)
243 pinv[i] = -1;
244 for (j = 0; j < m; j++)
245 qinv[j] = -1;
246
247 if (nzbias >= 0) /* put more in U */
248 {
249 Lbegin_p[0] = Ubegin[0] = rank = 0;
250
251 rank = singleton_cols(m, Bbegin, Bend, Bi, Bx, Btp, Bti, Btx,
252 Ubegin, Uindex, Uvalue, Lbegin_p, Lindex, Lvalue,
253 col_pivot, pinv, qinv, iwork1, iwork2, rank,
254 abstol);
255
256 rank = singleton_rows(m, Bbegin, Bend, Bi, Bx, Btp, Bti, Btx,
257 Ubegin, Uindex, Uvalue, Lbegin_p, Lindex, Lvalue,
258 col_pivot, pinv, qinv, iwork1, iwork2, rank,
259 abstol);
260 }
261 else /* put more in L */
262 {
263 Lbegin_p[0] = Ubegin[0] = rank = 0;
264
265 rank = singleton_rows(m, Bbegin, Bend, Bi, Bx, Btp, Bti, Btx,
266 Ubegin, Uindex, Uvalue, Lbegin_p, Lindex, Lvalue,
267 col_pivot, pinv, qinv, iwork1, iwork2, rank,
268 abstol);
269
270 rank = singleton_cols(m, Bbegin, Bend, Bi, Bx, Btp, Bti, Btx,
271 Ubegin, Uindex, Uvalue, Lbegin_p, Lindex, Lvalue,
272 col_pivot, pinv, qinv, iwork1, iwork2, rank,
273 abstol);
274 }
275
276 /* pinv, qinv were used as nonzero counters. Reset to -1 if not pivoted. */
277 for (i = 0; i < m; i++)
278 if (pinv[i] < 0)
279 pinv[i] = -1;
280 for (j = 0; j < m; j++)
281 if (qinv[j] < 0)
282 qinv[j] = -1;
283
284 this->matrix_nz = Bnz;
285 this->rank = rank;
286 return BASICLU_OK;
287 }
288
289
290 /*
291 * singleton_cols()
292 *
293 * The method successively removes singleton cols from an active submatrix.
294 * The active submatrix is composed of columns j for which qinv[j] < 0 and
295 * rows i for which pinv[i] < 0. When removing a singleton column and its
296 * associated row generates new singleton columns, these are appended to a
297 * queue. The method stops when the active submatrix has no more singleton
298 * columns.
299 *
300 * For each active column j iset[j] is the XOR of row indices in the column
301 * in the active submatrix. For a singleton column, this is its single row
302 * index. The technique is due to J. Gilbert and described in [1], ex 3.7.
303 *
304 * For each eliminated column its associated row is stored in U without the
305 * pivot element. The pivot elements are stored in col_pivot. For each
306 * eliminated pivot an empty column is appended to L.
307 *
308 * Pivot elements which are zero or less than abstol, and empty columns in
309 * the active submatrix are not eliminated. In these cases the matrix is
310 * numerically or structurally singular and the bump factorization handles
311 * it. (We want singularities at the end of the pivot sequence.)
312 *
313 * [1] T. Davis, "Direct methods for sparse linear systems"
314 */
315
singleton_cols(const lu_int m,const lu_int * Bbegin,const lu_int * Bend,const lu_int * Bi,const double * Bx,const lu_int * Btp,const lu_int * Bti,const double * Btx,lu_int * Up,lu_int * Ui,double * Ux,lu_int * Lp,lu_int * Li,double * Lx,double * col_pivot,lu_int * pinv,lu_int * qinv,lu_int * iset,lu_int * queue,lu_int rank,double abstol)316 static lu_int singleton_cols
317 (
318 const lu_int m,
319 const lu_int *Bbegin, /* B columnwise */
320 const lu_int *Bend,
321 const lu_int *Bi,
322 const double *Bx,
323 const lu_int *Btp, /* B rowwise */
324 const lu_int *Bti,
325 const double *Btx,
326 lu_int *Up,
327 lu_int *Ui,
328 double *Ux,
329 lu_int *Lp,
330 lu_int *Li,
331 double *Lx,
332 double *col_pivot,
333 lu_int *pinv,
334 lu_int *qinv,
335 lu_int *iset, /* size m workspace */
336 lu_int *queue, /* size m workspace */
337 lu_int rank,
338 double abstol
339 )
340 {
341 lu_int i, j, j2, nz, pos, put, end, front, tail, rk = rank;
342 double piv;
343
344 /* Build index sets and initialize queue. */
345 tail = 0;
346 for (j = 0; j < m; j++)
347 {
348 if (qinv[j] < 0)
349 {
350 nz = Bend[j] - Bbegin[j];
351 i = 0;
352 for (pos = Bbegin[j]; pos < Bend[j]; pos++)
353 i ^= Bi[pos]; /* put row into set j */
354 iset[j] = i;
355 qinv[j] = -nz-1; /* use as nonzero counter */
356 if (nz == 1)
357 queue[tail++] = j;
358 }
359 }
360
361 /* Eliminate singleton columns. */
362 put = Up [rank];
363 for (front = 0; front < tail; front++)
364 {
365 j = queue[front];
366 assert(qinv[j] == -2 || qinv[j] == -1);
367 if (qinv[j] == -1)
368 continue; /* empty column in active submatrix */
369 i = iset[j];
370 assert(i >= 0 && i < m);
371 assert(pinv[i] < 0);
372 end = Btp[i+1];
373 for (pos = Btp[i]; Bti[pos] != j; pos++) /* find pivot */
374 assert(pos < end-1);
375 piv = Btx[pos];
376 if (!piv || fabs(piv) < abstol)
377 continue; /* skip singularity */
378
379 /* Eliminate pivot. */
380 qinv[j] = rank;
381 pinv[i] = rank;
382 for (pos = Btp[i]; pos < end; pos++)
383 {
384 j2 = Bti[pos];
385 if (qinv[j2] < 0)
386 /* test is mandatory because the initial active submatrix may
387 not be the entire matrix (rows eliminated before) */
388 {
389 Ui[put] = j2;
390 Ux[put++] = Btx[pos];
391 iset[j2] ^= i; /* remove i from set j2 */
392 if (++qinv[j2] == -2)
393 queue[tail++] = j2; /* new singleton */
394 }
395 }
396 Up[rank+1] = put;
397 col_pivot[j] = piv;
398 rank++;
399 }
400
401 /* Put empty columns into L. */
402 pos = Lp[rk];
403 for ( ; rk < rank; rk++)
404 {
405 Li[pos++] = -1;
406 Lp[rk+1] = pos;
407 }
408 return rank;
409 }
410
411
412 /*
413 * singleton_rows()
414 *
415 * Analogeous singleton_cols except that for each singleton row the
416 * associated column is stored in L and divided by the pivot element. The
417 * pivot element is stored in col_pivot.
418 */
419
singleton_rows(const lu_int m,const lu_int * Bbegin,const lu_int * Bend,const lu_int * Bi,const double * Bx,const lu_int * Btp,const lu_int * Bti,const double * Btx,lu_int * Up,lu_int * Ui,double * Ux,lu_int * Lp,lu_int * Li,double * Lx,double * col_pivot,lu_int * pinv,lu_int * qinv,lu_int * iset,lu_int * queue,lu_int rank,double abstol)420 static lu_int singleton_rows
421 (
422 const lu_int m,
423 const lu_int *Bbegin, /* B columnwise */
424 const lu_int *Bend,
425 const lu_int *Bi,
426 const double *Bx,
427 const lu_int *Btp, /* B rowwise */
428 const lu_int *Bti,
429 const double *Btx,
430 lu_int *Up,
431 lu_int *Ui,
432 double *Ux,
433 lu_int *Lp,
434 lu_int *Li,
435 double *Lx,
436 double *col_pivot,
437 lu_int *pinv,
438 lu_int *qinv,
439 lu_int *iset, /* size m workspace */
440 lu_int *queue, /* size m workspace */
441 lu_int rank,
442 double abstol
443 )
444 {
445 lu_int i, j, i2, nz, pos, put, end, front, tail, rk = rank;
446 double piv;
447
448 /* Build index sets and initialize queue. */
449 tail = 0;
450 for (i = 0; i < m; i++)
451 {
452 if (pinv[i] < 0)
453 {
454 nz = Btp[i+1] - Btp[i];
455 j = 0;
456 for (pos = Btp[i]; pos < Btp[i+1]; pos++)
457 j ^= Bti[pos]; /* put column into set i */
458 iset[i] = j;
459 pinv[i] = -nz-1; /* use as nonzero counter */
460 if (nz == 1)
461 queue[tail++] = i;
462 }
463 }
464
465 /* Eliminate singleton rows. */
466 put = Lp[rank];
467 for (front = 0; front < tail; front++)
468 {
469 i = queue[front];
470 assert(pinv[i] == -2 || pinv[i] == -1);
471 if (pinv[i] == -1)
472 continue; /* empty column in active submatrix */
473 j = iset [i];
474 assert(j >= 0 && j < m);
475 assert(qinv[j] < 0);
476 end = Bend[j];
477 for (pos = Bbegin[j]; Bi[pos] != i; pos++) /* find pivot */
478 assert(pos < end-1);
479 piv = Bx[pos];
480 if (!piv || fabs(piv) < abstol)
481 continue; /* skip singularity */
482
483 /* Eliminate pivot. */
484 qinv[j] = rank;
485 pinv[i] = rank;
486 for (pos = Bbegin[j]; pos < end; pos++)
487 {
488 i2 = Bi[pos];
489 if (pinv[i2] < 0)
490 /* test is mandatory because the initial active submatrix may
491 not be the entire matrix (columns eliminated before) */
492 {
493 Li[put] = i2;
494 Lx[put++] = Bx[pos] / piv;
495 iset[i2] ^= j; /* remove j from set i2 */
496 if (++pinv[i2] == -2)
497 queue[tail++] = i2; /* new singleton */
498 }
499 }
500 Li[put++] = -1; /* terminate column */
501 Lp[rank+1] = put;
502 col_pivot[j] = piv;
503 rank++;
504 }
505
506 /* Put empty rows into U. */
507 pos = Up[rk];
508 for ( ; rk < rank; rk++)
509 Up[rk+1] = pos;
510
511 return rank;
512 }
513