1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
/*
* basiclu_obj_maxvolume.c
*
* Copyright (C) 2016-2017 ERGO-Code
*
*/
#include "lu_internal.h"
/*
* factorize() - factorize A[:,basis]
*/
static lu_int factorize(struct basiclu_object *obj,
const lu_int *Ap,
const lu_int *Ai,
const double *Ax,
const lu_int *basis)
{
double *xstore = obj->xstore;
const lu_int m = xstore[BASICLU_DIM];
lu_int *begin = NULL;
lu_int *end = NULL;
lu_int i, status = BASICLU_OK;
begin = malloc(m*sizeof(lu_int));
end = malloc(m*sizeof(lu_int));
if (!begin || !end) {
status = BASICLU_ERROR_out_of_memory;
goto cleanup;
}
for (i = 0; i < m; i++)
{
begin[i] = Ap[basis[i]];
end[i] = Ap[basis[i]+1];
}
status = basiclu_obj_factorize(obj, begin, end, Ai, Ax);
cleanup:
if (begin) free(begin);
if (end) free(end);
return status;
}
/*
* refactorize_if_needed() - refactorize the basis if required or favourable
*
* The basis matrix is refactorized if
* - the maximum number of updates is reached, or
* - the previous update had a large pivot error, or
* - it is favourable for performance
*
* factorize() is called for the actual factorization.
*
* Note: refactorize_if_needed() will not do an initial factorization.
*/
static lu_int refactorize_if_needed(struct basiclu_object *obj,
const lu_int *Ap,
const lu_int *Ai,
const double *Ax,
const lu_int *basis)
{
lu_int status = BASICLU_OK;
const double piverr_tol = 1e-8;
double *xstore = obj->xstore;
if (xstore[BASICLU_NFORREST] == xstore[BASICLU_DIM] ||
xstore[BASICLU_PIVOT_ERROR] > piverr_tol ||
xstore[BASICLU_UPDATE_COST] > 1.0)
status = factorize(obj, Ap, Ai, Ax, basis);
return status;
}
/*
* basiclu_obj_maxvolume() - one pass over columns of A doing basis updates
*
* For each column a_j not in B, compute lhs = B^{-1}*a_j and find the maximum
* entry lhs[imax]. If it is bigger than @volumetol in absolute value, then
* replace position imax of the basis by index j. On return *p_nupdate is the
* number of basis updates done.
*/
lu_int basiclu_obj_maxvolume
(
struct basiclu_object *obj,
lu_int ncol,
const lu_int Ap[],
const lu_int Ai[],
const double Ax[],
lu_int basis[],
lu_int isbasic[],
double volumetol,
lu_int *p_nupdate
)
{
lu_int i, j, k;
lu_int nzrhs, imax, begin, nupdate = 0;
double xtbl, xmax;
lu_int status = BASICLU_OK;
if (volumetol < 1.0)
{
status = BASICLU_ERROR_invalid_argument;
goto cleanup;
}
/* Compute initial factorization. */
status = factorize(obj, Ap, Ai, Ax, basis);
if (status != BASICLU_OK)
goto cleanup;
for (j = 0; j < ncol; j++)
{
if (isbasic[j])
continue;
/* compute B^{-1}*a_j */
nzrhs = Ap[j+1] - Ap[j];
begin = Ap[j];
status = basiclu_obj_solve_for_update(obj, nzrhs, Ai+begin, Ax+begin,
'N', 1);
if (status != BASICLU_OK)
goto cleanup;
/* Find the maximum entry. */
xmax = 0.0;
xtbl = 0.0;
imax = 0;
for (k = 0; k < obj->nzlhs; k++) {
i = obj->ilhs[k];
if (fabs(obj->lhs[i]) > xmax) {
xtbl = obj->lhs[i];
xmax = fabs(xtbl);
imax = i;
}
}
if (xmax <= volumetol)
continue;
/* Update basis. */
isbasic[basis[imax]] = 0;
isbasic[j] = 1;
basis[imax] = j;
nupdate++;
/* Prepare to update factorization. */
status = basiclu_obj_solve_for_update(obj, 0, &imax, NULL, 'T', 0);
if (status != BASICLU_OK)
goto cleanup;
status = basiclu_obj_update(obj, xtbl);
if (status != BASICLU_OK)
goto cleanup;
status = refactorize_if_needed(obj, Ap, Ai, Ax, basis);
if (status != BASICLU_OK)
goto cleanup;
}
cleanup:
if (p_nupdate) *p_nupdate = nupdate;
return status;
}