1 /*
2 * lu_markowitz.c
3 *
4 * Copyright (C) 2016-2019 ERGO-Code
5 *
6 * Search for pivot element with small Markowitz cost. An eligible pivot
7 * must be nonzero and satisfy
8 *
9 * (1) abs(piv) >= abstol,
10 * (2) abs(piv) >= reltol * max[pivot column].
11 *
12 * From all eligible pivots search for one that minimizes
13 *
14 * mc := (nnz[pivot row] - 1) * (nnz[pivot column] - 1).
15 *
16 * The search is terminated when maxsearch rows or columns with eligible pivots
17 * have been searched (if not before). The row and column of the cheapest one
18 * found is stored in pivot_row and pivot_col.
19 *
20 * When the active submatrix contains columns with column count = 0, then such a
21 * column is chosen immediately and pivot_row = -1 is returned. Otherwise, when
22 * the Markowitz search does not find a pivot that is nonzero and >= abstol,
23 * then pivot_col = pivot_row = -1 is returned. (The latter cannot happen in the
24 * current version of the code because lu_pivot() erases columns of the active
25 * submatrix whose maximum absolute value drops below abstol.)
26 *
27 * The Markowitz search is implemented as described in [1].
28 *
29 * [1] U. Suhl, L. Suhl, "Computing Sparse LU Factorizations for Large-Scale
30 * Linear Programming Bases", ORSA Journal on Computing (1990)
31 *
32 */
33
34 #include "lu_internal.h"
35 #include "lu_list.h"
36
lu_markowitz(struct lu * this)37 lu_int lu_markowitz(struct lu *this)
38 {
39 const lu_int m = this->m;
40 const lu_int *Wbegin = this->Wbegin;
41 const lu_int *Wend = this->Wend;
42 const lu_int *Windex = this->Windex;
43 const double *Wvalue = this->Wvalue;
44 const lu_int *colcount_flink = this->colcount_flink;
45 lu_int *rowcount_flink = this->rowcount_flink;
46 lu_int *rowcount_blink = this->rowcount_blink;
47 const double *colmax = this->col_pivot;
48 const double abstol = this->abstol;
49 const double reltol = this->reltol;
50 const lu_int maxsearch = this->maxsearch;
51 const lu_int search_rows = this->search_rows;
52 const lu_int nz_start = search_rows ?
53 MIN(this->min_colnz, this->min_rownz) : this->min_colnz;
54
55 lu_int i, j, pos, where, inext, nz, pivot_row, pivot_col;
56 lu_int nsearch, cheap, found, min_colnz, min_rownz;
57 double cmx, x, tol, tic[2];
58
59 /* integers for Markowitz cost must be 64 bit to prevent overflow */
60 const int_least64_t M = m;
61 int_least64_t nz1, nz2, mc, MC;
62
63 pivot_row = -1; /* row of best pivot so far */
64 pivot_col = -1; /* col of best pivot so far */
65 MC = M*M; /* Markowitz cost of best pivot so far */
66 nsearch = 0; /* count rows/columns searched */
67 min_colnz = -1; /* minimum col count in active submatrix */
68 min_rownz = -1; /* minimum row count in active submatrix */
69 assert(nz_start >= 1);
70
71 /* If the active submatrix contains empty columns, choose one and return
72 with pivot_row = -1. */
73 if (colcount_flink[m] != m)
74 {
75 pivot_col = colcount_flink[m];
76 assert(pivot_col >= 0 && pivot_col < m);
77 assert(Wend[pivot_col] == Wbegin[pivot_col]);
78 goto done;
79 }
80
81 for (nz = nz_start; nz <= m; nz++)
82 {
83 /* Search columns with nz nonzeros. */
84 for (j = colcount_flink[m+nz]; j < m; j = colcount_flink[j])
85 {
86 if (min_colnz == -1)
87 min_colnz = nz;
88 assert(Wend[j] - Wbegin[j] == nz);
89 cmx = colmax[j];
90 assert(cmx >= 0);
91 if (!cmx || cmx < abstol)
92 continue;
93 tol = fmax(abstol, reltol*cmx);
94 for (pos = Wbegin[j]; pos < Wend[j]; pos++)
95 {
96 x = fabs(Wvalue[pos]);
97 if (!x || x < tol)
98 continue;
99 i = Windex[pos];
100 assert(i >= 0 && i < m);
101 nz1 = nz;
102 nz2 = Wend[m+i] - Wbegin[m+i];
103 assert(nz2 >= 1);
104 mc = (nz1-1) * (nz2-1);
105 if (mc < MC)
106 {
107 MC = mc;
108 pivot_row = i;
109 pivot_col = j;
110 if (search_rows && MC <= (nz1-1)*(nz1-1))
111 goto done;
112 }
113 }
114 /* We have seen at least one eligible pivot in column j. */
115 assert(MC < M*M);
116 if (++nsearch >= maxsearch)
117 goto done;
118 }
119 assert(j == m+nz);
120
121 if (!search_rows)
122 continue;
123
124 /* Search rows with nz nonzeros. */
125 for (i = rowcount_flink[m+nz]; i < m; i = inext)
126 {
127 if (min_rownz == -1)
128 min_rownz = nz;
129 /* rowcount_flink[i] might be changed below, so keep a copy */
130 inext = rowcount_flink[i];
131 assert(Wend[m+i] - Wbegin[m+i] == nz);
132 cheap = 0; /* row has entries with Markowitz cost < MC? */
133 found = 0; /* eligible pivot found? */
134 for (pos = Wbegin[m+i]; pos < Wend[m+i]; pos++)
135 {
136 j = Windex[pos];
137 assert(j >= 0 && j < m);
138 nz1 = nz;
139 nz2 = Wend[j] - Wbegin[j];
140 assert(nz2 >= 1);
141 mc = (nz1-1) * (nz2-1);
142 if (mc >= MC)
143 continue;
144 cheap = 1;
145 cmx = colmax[j];
146 assert(cmx >= 0);
147 if (!cmx || cmx < abstol)
148 continue;
149 /* find position of pivot in column file */
150 for (where = Wbegin[j]; Windex[where] != i; where++)
151 assert(where < Wend[j] - 1);
152 x = fabs(Wvalue[where]);
153 if (x >= abstol && x >= reltol*cmx)
154 {
155 found = 1;
156 MC = mc;
157 pivot_row = i;
158 pivot_col = j;
159 if (MC <= nz1*(nz1-1))
160 goto done;
161 }
162 }
163 /* If row i has cheap entries but none of them is numerically
164 * acceptable, then don't search the row again until updated. */
165 if (cheap && !found)
166 {
167 lu_list_move(i, m+1, rowcount_flink, rowcount_blink, m, NULL);
168 }
169 else
170 {
171 assert(MC < M*M);
172 if (++nsearch >= maxsearch)
173 goto done;
174 }
175 }
176 assert(i == m+nz);
177 }
178
179 done:
180 this->pivot_row = pivot_row;
181 this->pivot_col = pivot_col;
182 this->nsearch_pivot += nsearch;
183 if (min_colnz >= 0)
184 this->min_colnz = min_colnz;
185 if (min_rownz >= 0)
186 this->min_rownz = min_rownz;
187 return BASICLU_OK;
188 }
189