1 | /* glplux.c */ |
---|
2 | |
---|
3 | /*********************************************************************** |
---|
4 | * This code is part of GLPK (GNU Linear Programming Kit). |
---|
5 | * |
---|
6 | * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, |
---|
7 | * 2009, 2010 Andrew Makhorin, Department for Applied Informatics, |
---|
8 | * Moscow Aviation Institute, Moscow, Russia. All rights reserved. |
---|
9 | * E-mail: <mao@gnu.org>. |
---|
10 | * |
---|
11 | * GLPK is free software: you can redistribute it and/or modify it |
---|
12 | * under the terms of the GNU General Public License as published by |
---|
13 | * the Free Software Foundation, either version 3 of the License, or |
---|
14 | * (at your option) any later version. |
---|
15 | * |
---|
16 | * GLPK is distributed in the hope that it will be useful, but WITHOUT |
---|
17 | * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
---|
18 | * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public |
---|
19 | * License for more details. |
---|
20 | * |
---|
21 | * You should have received a copy of the GNU General Public License |
---|
22 | * along with GLPK. If not, see <http://www.gnu.org/licenses/>. |
---|
23 | ***********************************************************************/ |
---|
24 | |
---|
25 | #include "glplux.h" |
---|
26 | #define xfault xerror |
---|
27 | #define dmp_create_poolx(size) dmp_create_pool() |
---|
28 | |
---|
29 | /*---------------------------------------------------------------------- |
---|
30 | // lux_create - create LU-factorization. |
---|
31 | // |
---|
32 | // SYNOPSIS |
---|
33 | // |
---|
34 | // #include "glplux.h" |
---|
35 | // LUX *lux_create(int n); |
---|
36 | // |
---|
37 | // DESCRIPTION |
---|
38 | // |
---|
39 | // The routine lux_create creates LU-factorization data structure for |
---|
40 | // a matrix of the order n. Initially the factorization corresponds to |
---|
41 | // the unity matrix (F = V = P = Q = I, so A = I). |
---|
42 | // |
---|
43 | // RETURNS |
---|
44 | // |
---|
45 | // The routine returns a pointer to the created LU-factorization data |
---|
46 | // structure, which represents the unity matrix of the order n. */ |
---|
47 | |
---|
48 | LUX *lux_create(int n) |
---|
49 | { LUX *lux; |
---|
50 | int k; |
---|
51 | if (n < 1) |
---|
52 | xfault("lux_create: n = %d; invalid parameter\n", n); |
---|
53 | lux = xmalloc(sizeof(LUX)); |
---|
54 | lux->n = n; |
---|
55 | lux->pool = dmp_create_poolx(sizeof(LUXELM)); |
---|
56 | lux->F_row = xcalloc(1+n, sizeof(LUXELM *)); |
---|
57 | lux->F_col = xcalloc(1+n, sizeof(LUXELM *)); |
---|
58 | lux->V_piv = xcalloc(1+n, sizeof(mpq_t)); |
---|
59 | lux->V_row = xcalloc(1+n, sizeof(LUXELM *)); |
---|
60 | lux->V_col = xcalloc(1+n, sizeof(LUXELM *)); |
---|
61 | lux->P_row = xcalloc(1+n, sizeof(int)); |
---|
62 | lux->P_col = xcalloc(1+n, sizeof(int)); |
---|
63 | lux->Q_row = xcalloc(1+n, sizeof(int)); |
---|
64 | lux->Q_col = xcalloc(1+n, sizeof(int)); |
---|
65 | for (k = 1; k <= n; k++) |
---|
66 | { lux->F_row[k] = lux->F_col[k] = NULL; |
---|
67 | mpq_init(lux->V_piv[k]); |
---|
68 | mpq_set_si(lux->V_piv[k], 1, 1); |
---|
69 | lux->V_row[k] = lux->V_col[k] = NULL; |
---|
70 | lux->P_row[k] = lux->P_col[k] = k; |
---|
71 | lux->Q_row[k] = lux->Q_col[k] = k; |
---|
72 | } |
---|
73 | lux->rank = n; |
---|
74 | return lux; |
---|
75 | } |
---|
76 | |
---|
77 | /*---------------------------------------------------------------------- |
---|
78 | // initialize - initialize LU-factorization data structures. |
---|
79 | // |
---|
80 | // This routine initializes data structures for subsequent computing |
---|
81 | // the LU-factorization of a given matrix A, which is specified by the |
---|
82 | // formal routine col. On exit V = A and F = P = Q = I, where I is the |
---|
83 | // unity matrix. */ |
---|
84 | |
---|
85 | static void initialize(LUX *lux, int (*col)(void *info, int j, |
---|
86 | int ind[], mpq_t val[]), void *info, LUXWKA *wka) |
---|
87 | { int n = lux->n; |
---|
88 | DMP *pool = lux->pool; |
---|
89 | LUXELM **F_row = lux->F_row; |
---|
90 | LUXELM **F_col = lux->F_col; |
---|
91 | mpq_t *V_piv = lux->V_piv; |
---|
92 | LUXELM **V_row = lux->V_row; |
---|
93 | LUXELM **V_col = lux->V_col; |
---|
94 | int *P_row = lux->P_row; |
---|
95 | int *P_col = lux->P_col; |
---|
96 | int *Q_row = lux->Q_row; |
---|
97 | int *Q_col = lux->Q_col; |
---|
98 | int *R_len = wka->R_len; |
---|
99 | int *R_head = wka->R_head; |
---|
100 | int *R_prev = wka->R_prev; |
---|
101 | int *R_next = wka->R_next; |
---|
102 | int *C_len = wka->C_len; |
---|
103 | int *C_head = wka->C_head; |
---|
104 | int *C_prev = wka->C_prev; |
---|
105 | int *C_next = wka->C_next; |
---|
106 | LUXELM *fij, *vij; |
---|
107 | int i, j, k, len, *ind; |
---|
108 | mpq_t *val; |
---|
109 | /* F := I */ |
---|
110 | for (i = 1; i <= n; i++) |
---|
111 | { while (F_row[i] != NULL) |
---|
112 | { fij = F_row[i], F_row[i] = fij->r_next; |
---|
113 | mpq_clear(fij->val); |
---|
114 | dmp_free_atom(pool, fij, sizeof(LUXELM)); |
---|
115 | } |
---|
116 | } |
---|
117 | for (j = 1; j <= n; j++) F_col[j] = NULL; |
---|
118 | /* V := 0 */ |
---|
119 | for (k = 1; k <= n; k++) mpq_set_si(V_piv[k], 0, 1); |
---|
120 | for (i = 1; i <= n; i++) |
---|
121 | { while (V_row[i] != NULL) |
---|
122 | { vij = V_row[i], V_row[i] = vij->r_next; |
---|
123 | mpq_clear(vij->val); |
---|
124 | dmp_free_atom(pool, vij, sizeof(LUXELM)); |
---|
125 | } |
---|
126 | } |
---|
127 | for (j = 1; j <= n; j++) V_col[j] = NULL; |
---|
128 | /* V := A */ |
---|
129 | ind = xcalloc(1+n, sizeof(int)); |
---|
130 | val = xcalloc(1+n, sizeof(mpq_t)); |
---|
131 | for (k = 1; k <= n; k++) mpq_init(val[k]); |
---|
132 | for (j = 1; j <= n; j++) |
---|
133 | { /* obtain j-th column of matrix A */ |
---|
134 | len = col(info, j, ind, val); |
---|
135 | if (!(0 <= len && len <= n)) |
---|
136 | xfault("lux_decomp: j = %d: len = %d; invalid column length" |
---|
137 | "\n", j, len); |
---|
138 | /* copy elements of j-th column to matrix V */ |
---|
139 | for (k = 1; k <= len; k++) |
---|
140 | { /* get row index of a[i,j] */ |
---|
141 | i = ind[k]; |
---|
142 | if (!(1 <= i && i <= n)) |
---|
143 | xfault("lux_decomp: j = %d: i = %d; row index out of ran" |
---|
144 | "ge\n", j, i); |
---|
145 | /* check for duplicate indices */ |
---|
146 | if (V_row[i] != NULL && V_row[i]->j == j) |
---|
147 | xfault("lux_decomp: j = %d: i = %d; duplicate row indice" |
---|
148 | "s not allowed\n", j, i); |
---|
149 | /* check for zero value */ |
---|
150 | if (mpq_sgn(val[k]) == 0) |
---|
151 | xfault("lux_decomp: j = %d: i = %d; zero elements not al" |
---|
152 | "lowed\n", j, i); |
---|
153 | /* add new element v[i,j] = a[i,j] to V */ |
---|
154 | vij = dmp_get_atom(pool, sizeof(LUXELM)); |
---|
155 | vij->i = i, vij->j = j; |
---|
156 | mpq_init(vij->val); |
---|
157 | mpq_set(vij->val, val[k]); |
---|
158 | vij->r_prev = NULL; |
---|
159 | vij->r_next = V_row[i]; |
---|
160 | vij->c_prev = NULL; |
---|
161 | vij->c_next = V_col[j]; |
---|
162 | if (vij->r_next != NULL) vij->r_next->r_prev = vij; |
---|
163 | if (vij->c_next != NULL) vij->c_next->c_prev = vij; |
---|
164 | V_row[i] = V_col[j] = vij; |
---|
165 | } |
---|
166 | } |
---|
167 | xfree(ind); |
---|
168 | for (k = 1; k <= n; k++) mpq_clear(val[k]); |
---|
169 | xfree(val); |
---|
170 | /* P := Q := I */ |
---|
171 | for (k = 1; k <= n; k++) |
---|
172 | P_row[k] = P_col[k] = Q_row[k] = Q_col[k] = k; |
---|
173 | /* the rank of A and V is not determined yet */ |
---|
174 | lux->rank = -1; |
---|
175 | /* initially the entire matrix V is active */ |
---|
176 | /* determine its row lengths */ |
---|
177 | for (i = 1; i <= n; i++) |
---|
178 | { len = 0; |
---|
179 | for (vij = V_row[i]; vij != NULL; vij = vij->r_next) len++; |
---|
180 | R_len[i] = len; |
---|
181 | } |
---|
182 | /* build linked lists of active rows */ |
---|
183 | for (len = 0; len <= n; len++) R_head[len] = 0; |
---|
184 | for (i = 1; i <= n; i++) |
---|
185 | { len = R_len[i]; |
---|
186 | R_prev[i] = 0; |
---|
187 | R_next[i] = R_head[len]; |
---|
188 | if (R_next[i] != 0) R_prev[R_next[i]] = i; |
---|
189 | R_head[len] = i; |
---|
190 | } |
---|
191 | /* determine its column lengths */ |
---|
192 | for (j = 1; j <= n; j++) |
---|
193 | { len = 0; |
---|
194 | for (vij = V_col[j]; vij != NULL; vij = vij->c_next) len++; |
---|
195 | C_len[j] = len; |
---|
196 | } |
---|
197 | /* build linked lists of active columns */ |
---|
198 | for (len = 0; len <= n; len++) C_head[len] = 0; |
---|
199 | for (j = 1; j <= n; j++) |
---|
200 | { len = C_len[j]; |
---|
201 | C_prev[j] = 0; |
---|
202 | C_next[j] = C_head[len]; |
---|
203 | if (C_next[j] != 0) C_prev[C_next[j]] = j; |
---|
204 | C_head[len] = j; |
---|
205 | } |
---|
206 | return; |
---|
207 | } |
---|
208 | |
---|
209 | /*---------------------------------------------------------------------- |
---|
210 | // find_pivot - choose a pivot element. |
---|
211 | // |
---|
212 | // This routine chooses a pivot element v[p,q] in the active submatrix |
---|
213 | // of matrix U = P*V*Q. |
---|
214 | // |
---|
215 | // It is assumed that on entry the matrix U has the following partially |
---|
216 | // triangularized form: |
---|
217 | // |
---|
218 | // 1 k n |
---|
219 | // 1 x x x x x x x x x x |
---|
220 | // . x x x x x x x x x |
---|
221 | // . . x x x x x x x x |
---|
222 | // . . . x x x x x x x |
---|
223 | // k . . . . * * * * * * |
---|
224 | // . . . . * * * * * * |
---|
225 | // . . . . * * * * * * |
---|
226 | // . . . . * * * * * * |
---|
227 | // . . . . * * * * * * |
---|
228 | // n . . . . * * * * * * |
---|
229 | // |
---|
230 | // where rows and columns k, k+1, ..., n belong to the active submatrix |
---|
231 | // (elements of the active submatrix are marked by '*'). |
---|
232 | // |
---|
233 | // Since the matrix U = P*V*Q is not stored, the routine works with the |
---|
234 | // matrix V. It is assumed that the row-wise representation corresponds |
---|
235 | // to the matrix V, but the column-wise representation corresponds to |
---|
236 | // the active submatrix of the matrix V, i.e. elements of the matrix V, |
---|
237 | // which does not belong to the active submatrix, are missing from the |
---|
238 | // column linked lists. It is also assumed that each active row of the |
---|
239 | // matrix V is in the set R[len], where len is number of non-zeros in |
---|
240 | // the row, and each active column of the matrix V is in the set C[len], |
---|
241 | // where len is number of non-zeros in the column (in the latter case |
---|
242 | // only elements of the active submatrix are counted; such elements are |
---|
243 | // marked by '*' on the figure above). |
---|
244 | // |
---|
245 | // Due to exact arithmetic any non-zero element of the active submatrix |
---|
246 | // can be chosen as a pivot. However, to keep sparsity of the matrix V |
---|
247 | // the routine uses Markowitz strategy, trying to choose such element |
---|
248 | // v[p,q], which has smallest Markowitz cost (nr[p]-1) * (nc[q]-1), |
---|
249 | // where nr[p] and nc[q] are the number of non-zero elements, resp., in |
---|
250 | // p-th row and in q-th column of the active submatrix. |
---|
251 | // |
---|
252 | // In order to reduce the search, i.e. not to walk through all elements |
---|
253 | // of the active submatrix, the routine exploits a technique proposed by |
---|
254 | // I.Duff. This technique is based on using the sets R[len] and C[len] |
---|
255 | // of active rows and columns. |
---|
256 | // |
---|
257 | // On exit the routine returns a pointer to a pivot v[p,q] chosen, or |
---|
258 | // NULL, if the active submatrix is empty. */ |
---|
259 | |
---|
260 | static LUXELM *find_pivot(LUX *lux, LUXWKA *wka) |
---|
261 | { int n = lux->n; |
---|
262 | LUXELM **V_row = lux->V_row; |
---|
263 | LUXELM **V_col = lux->V_col; |
---|
264 | int *R_len = wka->R_len; |
---|
265 | int *R_head = wka->R_head; |
---|
266 | int *R_next = wka->R_next; |
---|
267 | int *C_len = wka->C_len; |
---|
268 | int *C_head = wka->C_head; |
---|
269 | int *C_next = wka->C_next; |
---|
270 | LUXELM *piv, *some, *vij; |
---|
271 | int i, j, len, min_len, ncand, piv_lim = 5; |
---|
272 | double best, cost; |
---|
273 | /* nothing is chosen so far */ |
---|
274 | piv = NULL, best = DBL_MAX, ncand = 0; |
---|
275 | /* if in the active submatrix there is a column that has the only |
---|
276 | non-zero (column singleton), choose it as a pivot */ |
---|
277 | j = C_head[1]; |
---|
278 | if (j != 0) |
---|
279 | { xassert(C_len[j] == 1); |
---|
280 | piv = V_col[j]; |
---|
281 | xassert(piv != NULL && piv->c_next == NULL); |
---|
282 | goto done; |
---|
283 | } |
---|
284 | /* if in the active submatrix there is a row that has the only |
---|
285 | non-zero (row singleton), choose it as a pivot */ |
---|
286 | i = R_head[1]; |
---|
287 | if (i != 0) |
---|
288 | { xassert(R_len[i] == 1); |
---|
289 | piv = V_row[i]; |
---|
290 | xassert(piv != NULL && piv->r_next == NULL); |
---|
291 | goto done; |
---|
292 | } |
---|
293 | /* there are no singletons in the active submatrix; walk through |
---|
294 | other non-empty rows and columns */ |
---|
295 | for (len = 2; len <= n; len++) |
---|
296 | { /* consider active columns having len non-zeros */ |
---|
297 | for (j = C_head[len]; j != 0; j = C_next[j]) |
---|
298 | { /* j-th column has len non-zeros */ |
---|
299 | /* find an element in the row of minimal length */ |
---|
300 | some = NULL, min_len = INT_MAX; |
---|
301 | for (vij = V_col[j]; vij != NULL; vij = vij->c_next) |
---|
302 | { if (min_len > R_len[vij->i]) |
---|
303 | some = vij, min_len = R_len[vij->i]; |
---|
304 | /* if Markowitz cost of this element is not greater than |
---|
305 | (len-1)**2, it can be chosen right now; this heuristic |
---|
306 | reduces the search and works well in many cases */ |
---|
307 | if (min_len <= len) |
---|
308 | { piv = some; |
---|
309 | goto done; |
---|
310 | } |
---|
311 | } |
---|
312 | /* j-th column has been scanned */ |
---|
313 | /* the minimal element found is a next pivot candidate */ |
---|
314 | xassert(some != NULL); |
---|
315 | ncand++; |
---|
316 | /* compute its Markowitz cost */ |
---|
317 | cost = (double)(min_len - 1) * (double)(len - 1); |
---|
318 | /* choose between the current candidate and this element */ |
---|
319 | if (cost < best) piv = some, best = cost; |
---|
320 | /* if piv_lim candidates have been considered, there is a |
---|
321 | doubt that a much better candidate exists; therefore it |
---|
322 | is the time to terminate the search */ |
---|
323 | if (ncand == piv_lim) goto done; |
---|
324 | } |
---|
325 | /* now consider active rows having len non-zeros */ |
---|
326 | for (i = R_head[len]; i != 0; i = R_next[i]) |
---|
327 | { /* i-th row has len non-zeros */ |
---|
328 | /* find an element in the column of minimal length */ |
---|
329 | some = NULL, min_len = INT_MAX; |
---|
330 | for (vij = V_row[i]; vij != NULL; vij = vij->r_next) |
---|
331 | { if (min_len > C_len[vij->j]) |
---|
332 | some = vij, min_len = C_len[vij->j]; |
---|
333 | /* if Markowitz cost of this element is not greater than |
---|
334 | (len-1)**2, it can be chosen right now; this heuristic |
---|
335 | reduces the search and works well in many cases */ |
---|
336 | if (min_len <= len) |
---|
337 | { piv = some; |
---|
338 | goto done; |
---|
339 | } |
---|
340 | } |
---|
341 | /* i-th row has been scanned */ |
---|
342 | /* the minimal element found is a next pivot candidate */ |
---|
343 | xassert(some != NULL); |
---|
344 | ncand++; |
---|
345 | /* compute its Markowitz cost */ |
---|
346 | cost = (double)(len - 1) * (double)(min_len - 1); |
---|
347 | /* choose between the current candidate and this element */ |
---|
348 | if (cost < best) piv = some, best = cost; |
---|
349 | /* if piv_lim candidates have been considered, there is a |
---|
350 | doubt that a much better candidate exists; therefore it |
---|
351 | is the time to terminate the search */ |
---|
352 | if (ncand == piv_lim) goto done; |
---|
353 | } |
---|
354 | } |
---|
355 | done: /* bring the pivot v[p,q] to the factorizing routine */ |
---|
356 | return piv; |
---|
357 | } |
---|
358 | |
---|
359 | /*---------------------------------------------------------------------- |
---|
360 | // eliminate - perform gaussian elimination. |
---|
361 | // |
---|
362 | // This routine performs elementary gaussian transformations in order |
---|
363 | // to eliminate subdiagonal elements in the k-th column of the matrix |
---|
364 | // U = P*V*Q using the pivot element u[k,k], where k is the number of |
---|
365 | // the current elimination step. |
---|
366 | // |
---|
367 | // The parameter piv specifies the pivot element v[p,q] = u[k,k]. |
---|
368 | // |
---|
369 | // Each time when the routine applies the elementary transformation to |
---|
370 | // a non-pivot row of the matrix V, it stores the corresponding element |
---|
371 | // to the matrix F in order to keep the main equality A = F*V. |
---|
372 | // |
---|
373 | // The routine assumes that on entry the matrices L = P*F*inv(P) and |
---|
374 | // U = P*V*Q are the following: |
---|
375 | // |
---|
376 | // 1 k 1 k n |
---|
377 | // 1 1 . . . . . . . . . 1 x x x x x x x x x x |
---|
378 | // x 1 . . . . . . . . . x x x x x x x x x |
---|
379 | // x x 1 . . . . . . . . . x x x x x x x x |
---|
380 | // x x x 1 . . . . . . . . . x x x x x x x |
---|
381 | // k x x x x 1 . . . . . k . . . . * * * * * * |
---|
382 | // x x x x _ 1 . . . . . . . . # * * * * * |
---|
383 | // x x x x _ . 1 . . . . . . . # * * * * * |
---|
384 | // x x x x _ . . 1 . . . . . . # * * * * * |
---|
385 | // x x x x _ . . . 1 . . . . . # * * * * * |
---|
386 | // n x x x x _ . . . . 1 n . . . . # * * * * * |
---|
387 | // |
---|
388 | // matrix L matrix U |
---|
389 | // |
---|
390 | // where rows and columns of the matrix U with numbers k, k+1, ..., n |
---|
391 | // form the active submatrix (eliminated elements are marked by '#' and |
---|
392 | // other elements of the active submatrix are marked by '*'). Note that |
---|
393 | // each eliminated non-zero element u[i,k] of the matrix U gives the |
---|
394 | // corresponding element l[i,k] of the matrix L (marked by '_'). |
---|
395 | // |
---|
396 | // Actually all operations are performed on the matrix V. Should note |
---|
397 | // that the row-wise representation corresponds to the matrix V, but the |
---|
398 | // column-wise representation corresponds to the active submatrix of the |
---|
399 | // matrix V, i.e. elements of the matrix V, which doesn't belong to the |
---|
400 | // active submatrix, are missing from the column linked lists. |
---|
401 | // |
---|
402 | // Let u[k,k] = v[p,q] be the pivot. In order to eliminate subdiagonal |
---|
403 | // elements u[i',k] = v[i,q], i' = k+1, k+2, ..., n, the routine applies |
---|
404 | // the following elementary gaussian transformations: |
---|
405 | // |
---|
406 | // (i-th row of V) := (i-th row of V) - f[i,p] * (p-th row of V), |
---|
407 | // |
---|
408 | // where f[i,p] = v[i,q] / v[p,q] is a gaussian multiplier. |
---|
409 | // |
---|
410 | // Additionally, in order to keep the main equality A = F*V, each time |
---|
411 | // when the routine applies the transformation to i-th row of the matrix |
---|
412 | // V, it also adds f[i,p] as a new element to the matrix F. |
---|
413 | // |
---|
414 | // IMPORTANT: On entry the working arrays flag and work should contain |
---|
415 | // zeros. This status is provided by the routine on exit. */ |
---|
416 | |
---|
417 | static void eliminate(LUX *lux, LUXWKA *wka, LUXELM *piv, int flag[], |
---|
418 | mpq_t work[]) |
---|
419 | { DMP *pool = lux->pool; |
---|
420 | LUXELM **F_row = lux->F_row; |
---|
421 | LUXELM **F_col = lux->F_col; |
---|
422 | mpq_t *V_piv = lux->V_piv; |
---|
423 | LUXELM **V_row = lux->V_row; |
---|
424 | LUXELM **V_col = lux->V_col; |
---|
425 | int *R_len = wka->R_len; |
---|
426 | int *R_head = wka->R_head; |
---|
427 | int *R_prev = wka->R_prev; |
---|
428 | int *R_next = wka->R_next; |
---|
429 | int *C_len = wka->C_len; |
---|
430 | int *C_head = wka->C_head; |
---|
431 | int *C_prev = wka->C_prev; |
---|
432 | int *C_next = wka->C_next; |
---|
433 | LUXELM *fip, *vij, *vpj, *viq, *next; |
---|
434 | mpq_t temp; |
---|
435 | int i, j, p, q; |
---|
436 | mpq_init(temp); |
---|
437 | /* determine row and column indices of the pivot v[p,q] */ |
---|
438 | xassert(piv != NULL); |
---|
439 | p = piv->i, q = piv->j; |
---|
440 | /* remove p-th (pivot) row from the active set; it will never |
---|
441 | return there */ |
---|
442 | if (R_prev[p] == 0) |
---|
443 | R_head[R_len[p]] = R_next[p]; |
---|
444 | else |
---|
445 | R_next[R_prev[p]] = R_next[p]; |
---|
446 | if (R_next[p] == 0) |
---|
447 | ; |
---|
448 | else |
---|
449 | R_prev[R_next[p]] = R_prev[p]; |
---|
450 | /* remove q-th (pivot) column from the active set; it will never |
---|
451 | return there */ |
---|
452 | if (C_prev[q] == 0) |
---|
453 | C_head[C_len[q]] = C_next[q]; |
---|
454 | else |
---|
455 | C_next[C_prev[q]] = C_next[q]; |
---|
456 | if (C_next[q] == 0) |
---|
457 | ; |
---|
458 | else |
---|
459 | C_prev[C_next[q]] = C_prev[q]; |
---|
460 | /* store the pivot value in a separate array */ |
---|
461 | mpq_set(V_piv[p], piv->val); |
---|
462 | /* remove the pivot from p-th row */ |
---|
463 | if (piv->r_prev == NULL) |
---|
464 | V_row[p] = piv->r_next; |
---|
465 | else |
---|
466 | piv->r_prev->r_next = piv->r_next; |
---|
467 | if (piv->r_next == NULL) |
---|
468 | ; |
---|
469 | else |
---|
470 | piv->r_next->r_prev = piv->r_prev; |
---|
471 | R_len[p]--; |
---|
472 | /* remove the pivot from q-th column */ |
---|
473 | if (piv->c_prev == NULL) |
---|
474 | V_col[q] = piv->c_next; |
---|
475 | else |
---|
476 | piv->c_prev->c_next = piv->c_next; |
---|
477 | if (piv->c_next == NULL) |
---|
478 | ; |
---|
479 | else |
---|
480 | piv->c_next->c_prev = piv->c_prev; |
---|
481 | C_len[q]--; |
---|
482 | /* free the space occupied by the pivot */ |
---|
483 | mpq_clear(piv->val); |
---|
484 | dmp_free_atom(pool, piv, sizeof(LUXELM)); |
---|
485 | /* walk through p-th (pivot) row, which already does not contain |
---|
486 | the pivot v[p,q], and do the following... */ |
---|
487 | for (vpj = V_row[p]; vpj != NULL; vpj = vpj->r_next) |
---|
488 | { /* get column index of v[p,j] */ |
---|
489 | j = vpj->j; |
---|
490 | /* store v[p,j] in the working array */ |
---|
491 | flag[j] = 1; |
---|
492 | mpq_set(work[j], vpj->val); |
---|
493 | /* remove j-th column from the active set; it will return there |
---|
494 | later with a new length */ |
---|
495 | if (C_prev[j] == 0) |
---|
496 | C_head[C_len[j]] = C_next[j]; |
---|
497 | else |
---|
498 | C_next[C_prev[j]] = C_next[j]; |
---|
499 | if (C_next[j] == 0) |
---|
500 | ; |
---|
501 | else |
---|
502 | C_prev[C_next[j]] = C_prev[j]; |
---|
503 | /* v[p,j] leaves the active submatrix, so remove it from j-th |
---|
504 | column; however, v[p,j] is kept in p-th row */ |
---|
505 | if (vpj->c_prev == NULL) |
---|
506 | V_col[j] = vpj->c_next; |
---|
507 | else |
---|
508 | vpj->c_prev->c_next = vpj->c_next; |
---|
509 | if (vpj->c_next == NULL) |
---|
510 | ; |
---|
511 | else |
---|
512 | vpj->c_next->c_prev = vpj->c_prev; |
---|
513 | C_len[j]--; |
---|
514 | } |
---|
515 | /* now walk through q-th (pivot) column, which already does not |
---|
516 | contain the pivot v[p,q], and perform gaussian elimination */ |
---|
517 | while (V_col[q] != NULL) |
---|
518 | { /* element v[i,q] has to be eliminated */ |
---|
519 | viq = V_col[q]; |
---|
520 | /* get row index of v[i,q] */ |
---|
521 | i = viq->i; |
---|
522 | /* remove i-th row from the active set; later it will return |
---|
523 | there with a new length */ |
---|
524 | if (R_prev[i] == 0) |
---|
525 | R_head[R_len[i]] = R_next[i]; |
---|
526 | else |
---|
527 | R_next[R_prev[i]] = R_next[i]; |
---|
528 | if (R_next[i] == 0) |
---|
529 | ; |
---|
530 | else |
---|
531 | R_prev[R_next[i]] = R_prev[i]; |
---|
532 | /* compute gaussian multiplier f[i,p] = v[i,q] / v[p,q] and |
---|
533 | store it in the matrix F */ |
---|
534 | fip = dmp_get_atom(pool, sizeof(LUXELM)); |
---|
535 | fip->i = i, fip->j = p; |
---|
536 | mpq_init(fip->val); |
---|
537 | mpq_div(fip->val, viq->val, V_piv[p]); |
---|
538 | fip->r_prev = NULL; |
---|
539 | fip->r_next = F_row[i]; |
---|
540 | fip->c_prev = NULL; |
---|
541 | fip->c_next = F_col[p]; |
---|
542 | if (fip->r_next != NULL) fip->r_next->r_prev = fip; |
---|
543 | if (fip->c_next != NULL) fip->c_next->c_prev = fip; |
---|
544 | F_row[i] = F_col[p] = fip; |
---|
545 | /* v[i,q] has to be eliminated, so remove it from i-th row */ |
---|
546 | if (viq->r_prev == NULL) |
---|
547 | V_row[i] = viq->r_next; |
---|
548 | else |
---|
549 | viq->r_prev->r_next = viq->r_next; |
---|
550 | if (viq->r_next == NULL) |
---|
551 | ; |
---|
552 | else |
---|
553 | viq->r_next->r_prev = viq->r_prev; |
---|
554 | R_len[i]--; |
---|
555 | /* and also from q-th column */ |
---|
556 | V_col[q] = viq->c_next; |
---|
557 | C_len[q]--; |
---|
558 | /* free the space occupied by v[i,q] */ |
---|
559 | mpq_clear(viq->val); |
---|
560 | dmp_free_atom(pool, viq, sizeof(LUXELM)); |
---|
561 | /* perform gaussian transformation: |
---|
562 | (i-th row) := (i-th row) - f[i,p] * (p-th row) |
---|
563 | note that now p-th row, which is in the working array, |
---|
564 | does not contain the pivot v[p,q], and i-th row does not |
---|
565 | contain the element v[i,q] to be eliminated */ |
---|
566 | /* walk through i-th row and transform existing non-zero |
---|
567 | elements */ |
---|
568 | for (vij = V_row[i]; vij != NULL; vij = next) |
---|
569 | { next = vij->r_next; |
---|
570 | /* get column index of v[i,j] */ |
---|
571 | j = vij->j; |
---|
572 | /* v[i,j] := v[i,j] - f[i,p] * v[p,j] */ |
---|
573 | if (flag[j]) |
---|
574 | { /* v[p,j] != 0 */ |
---|
575 | flag[j] = 0; |
---|
576 | mpq_mul(temp, fip->val, work[j]); |
---|
577 | mpq_sub(vij->val, vij->val, temp); |
---|
578 | if (mpq_sgn(vij->val) == 0) |
---|
579 | { /* new v[i,j] is zero, so remove it from the active |
---|
580 | submatrix */ |
---|
581 | /* remove v[i,j] from i-th row */ |
---|
582 | if (vij->r_prev == NULL) |
---|
583 | V_row[i] = vij->r_next; |
---|
584 | else |
---|
585 | vij->r_prev->r_next = vij->r_next; |
---|
586 | if (vij->r_next == NULL) |
---|
587 | ; |
---|
588 | else |
---|
589 | vij->r_next->r_prev = vij->r_prev; |
---|
590 | R_len[i]--; |
---|
591 | /* remove v[i,j] from j-th column */ |
---|
592 | if (vij->c_prev == NULL) |
---|
593 | V_col[j] = vij->c_next; |
---|
594 | else |
---|
595 | vij->c_prev->c_next = vij->c_next; |
---|
596 | if (vij->c_next == NULL) |
---|
597 | ; |
---|
598 | else |
---|
599 | vij->c_next->c_prev = vij->c_prev; |
---|
600 | C_len[j]--; |
---|
601 | /* free the space occupied by v[i,j] */ |
---|
602 | mpq_clear(vij->val); |
---|
603 | dmp_free_atom(pool, vij, sizeof(LUXELM)); |
---|
604 | } |
---|
605 | } |
---|
606 | } |
---|
607 | /* now flag is the pattern of the set v[p,*] \ v[i,*] */ |
---|
608 | /* walk through p-th (pivot) row and create new elements in |
---|
609 | i-th row, which appear due to fill-in */ |
---|
610 | for (vpj = V_row[p]; vpj != NULL; vpj = vpj->r_next) |
---|
611 | { j = vpj->j; |
---|
612 | if (flag[j]) |
---|
613 | { /* create new non-zero v[i,j] = 0 - f[i,p] * v[p,j] and |
---|
614 | add it to i-th row and j-th column */ |
---|
615 | vij = dmp_get_atom(pool, sizeof(LUXELM)); |
---|
616 | vij->i = i, vij->j = j; |
---|
617 | mpq_init(vij->val); |
---|
618 | mpq_mul(vij->val, fip->val, work[j]); |
---|
619 | mpq_neg(vij->val, vij->val); |
---|
620 | vij->r_prev = NULL; |
---|
621 | vij->r_next = V_row[i]; |
---|
622 | vij->c_prev = NULL; |
---|
623 | vij->c_next = V_col[j]; |
---|
624 | if (vij->r_next != NULL) vij->r_next->r_prev = vij; |
---|
625 | if (vij->c_next != NULL) vij->c_next->c_prev = vij; |
---|
626 | V_row[i] = V_col[j] = vij; |
---|
627 | R_len[i]++, C_len[j]++; |
---|
628 | } |
---|
629 | else |
---|
630 | { /* there is no fill-in, because v[i,j] already exists in |
---|
631 | i-th row; restore the flag, which was reset before */ |
---|
632 | flag[j] = 1; |
---|
633 | } |
---|
634 | } |
---|
635 | /* now i-th row has been completely transformed and can return |
---|
636 | to the active set with a new length */ |
---|
637 | R_prev[i] = 0; |
---|
638 | R_next[i] = R_head[R_len[i]]; |
---|
639 | if (R_next[i] != 0) R_prev[R_next[i]] = i; |
---|
640 | R_head[R_len[i]] = i; |
---|
641 | } |
---|
642 | /* at this point q-th (pivot) column must be empty */ |
---|
643 | xassert(C_len[q] == 0); |
---|
644 | /* walk through p-th (pivot) row again and do the following... */ |
---|
645 | for (vpj = V_row[p]; vpj != NULL; vpj = vpj->r_next) |
---|
646 | { /* get column index of v[p,j] */ |
---|
647 | j = vpj->j; |
---|
648 | /* erase v[p,j] from the working array */ |
---|
649 | flag[j] = 0; |
---|
650 | mpq_set_si(work[j], 0, 1); |
---|
651 | /* now j-th column has been completely transformed, so it can |
---|
652 | return to the active list with a new length */ |
---|
653 | C_prev[j] = 0; |
---|
654 | C_next[j] = C_head[C_len[j]]; |
---|
655 | if (C_next[j] != 0) C_prev[C_next[j]] = j; |
---|
656 | C_head[C_len[j]] = j; |
---|
657 | } |
---|
658 | mpq_clear(temp); |
---|
659 | /* return to the factorizing routine */ |
---|
660 | return; |
---|
661 | } |
---|
662 | |
---|
663 | /*---------------------------------------------------------------------- |
---|
664 | // lux_decomp - compute LU-factorization. |
---|
665 | // |
---|
666 | // SYNOPSIS |
---|
667 | // |
---|
668 | // #include "glplux.h" |
---|
669 | // int lux_decomp(LUX *lux, int (*col)(void *info, int j, int ind[], |
---|
670 | // mpq_t val[]), void *info); |
---|
671 | // |
---|
672 | // DESCRIPTION |
---|
673 | // |
---|
674 | // The routine lux_decomp computes LU-factorization of a given square |
---|
675 | // matrix A. |
---|
676 | // |
---|
677 | // The parameter lux specifies LU-factorization data structure built by |
---|
678 | // means of the routine lux_create. |
---|
679 | // |
---|
680 | // The formal routine col specifies the original matrix A. In order to |
---|
681 | // obtain j-th column of the matrix A the routine lux_decomp calls the |
---|
682 | // routine col with the parameter j (1 <= j <= n, where n is the order |
---|
683 | // of A). In response the routine col should store row indices and |
---|
684 | // numerical values of non-zero elements of j-th column of A to the |
---|
685 | // locations ind[1], ..., ind[len] and val[1], ..., val[len], resp., |
---|
686 | // where len is the number of non-zeros in j-th column, which should be |
---|
687 | // returned on exit. Neiter zero nor duplicate elements are allowed. |
---|
688 | // |
---|
689 | // The parameter info is a transit pointer passed to the formal routine |
---|
690 | // col; it can be used for various purposes. |
---|
691 | // |
---|
692 | // RETURNS |
---|
693 | // |
---|
694 | // The routine lux_decomp returns the singularity flag. Zero flag means |
---|
695 | // that the original matrix A is non-singular while non-zero flag means |
---|
696 | // that A is (exactly!) singular. |
---|
697 | // |
---|
698 | // Note that LU-factorization is valid in both cases, however, in case |
---|
699 | // of singularity some rows of the matrix V (including pivot elements) |
---|
700 | // will be empty. |
---|
701 | // |
---|
702 | // REPAIRING SINGULAR MATRIX |
---|
703 | // |
---|
704 | // If the routine lux_decomp returns non-zero flag, it provides all |
---|
705 | // necessary information that can be used for "repairing" the matrix A, |
---|
706 | // where "repairing" means replacing linearly dependent columns of the |
---|
707 | // matrix A by appropriate columns of the unity matrix. This feature is |
---|
708 | // needed when the routine lux_decomp is used for reinverting the basis |
---|
709 | // matrix within the simplex method procedure. |
---|
710 | // |
---|
711 | // On exit linearly dependent columns of the matrix U have the numbers |
---|
712 | // rank+1, rank+2, ..., n, where rank is the exact rank of the matrix A |
---|
713 | // stored by the routine to the member lux->rank. The correspondence |
---|
714 | // between columns of A and U is the same as between columns of V and U. |
---|
715 | // Thus, linearly dependent columns of the matrix A have the numbers |
---|
716 | // Q_col[rank+1], Q_col[rank+2], ..., Q_col[n], where Q_col is an array |
---|
717 | // representing the permutation matrix Q in column-like format. It is |
---|
718 | // understood that each j-th linearly dependent column of the matrix U |
---|
719 | // should be replaced by the unity vector, where all elements are zero |
---|
720 | // except the unity diagonal element u[j,j]. On the other hand j-th row |
---|
721 | // of the matrix U corresponds to the row of the matrix V (and therefore |
---|
722 | // of the matrix A) with the number P_row[j], where P_row is an array |
---|
723 | // representing the permutation matrix P in row-like format. Thus, each |
---|
724 | // j-th linearly dependent column of the matrix U should be replaced by |
---|
725 | // a column of the unity matrix with the number P_row[j]. |
---|
726 | // |
---|
727 | // The code that repairs the matrix A may look like follows: |
---|
728 | // |
---|
729 | // for (j = rank+1; j <= n; j++) |
---|
730 | // { replace column Q_col[j] of the matrix A by column P_row[j] of |
---|
731 | // the unity matrix; |
---|
732 | // } |
---|
733 | // |
---|
734 | // where rank, P_row, and Q_col are members of the structure LUX. */ |
---|
735 | |
---|
736 | int lux_decomp(LUX *lux, int (*col)(void *info, int j, int ind[], |
---|
737 | mpq_t val[]), void *info) |
---|
738 | { int n = lux->n; |
---|
739 | LUXELM **V_row = lux->V_row; |
---|
740 | LUXELM **V_col = lux->V_col; |
---|
741 | int *P_row = lux->P_row; |
---|
742 | int *P_col = lux->P_col; |
---|
743 | int *Q_row = lux->Q_row; |
---|
744 | int *Q_col = lux->Q_col; |
---|
745 | LUXELM *piv, *vij; |
---|
746 | LUXWKA *wka; |
---|
747 | int i, j, k, p, q, t, *flag; |
---|
748 | mpq_t *work; |
---|
749 | /* allocate working area */ |
---|
750 | wka = xmalloc(sizeof(LUXWKA)); |
---|
751 | wka->R_len = xcalloc(1+n, sizeof(int)); |
---|
752 | wka->R_head = xcalloc(1+n, sizeof(int)); |
---|
753 | wka->R_prev = xcalloc(1+n, sizeof(int)); |
---|
754 | wka->R_next = xcalloc(1+n, sizeof(int)); |
---|
755 | wka->C_len = xcalloc(1+n, sizeof(int)); |
---|
756 | wka->C_head = xcalloc(1+n, sizeof(int)); |
---|
757 | wka->C_prev = xcalloc(1+n, sizeof(int)); |
---|
758 | wka->C_next = xcalloc(1+n, sizeof(int)); |
---|
759 | /* initialize LU-factorization data structures */ |
---|
760 | initialize(lux, col, info, wka); |
---|
761 | /* allocate working arrays */ |
---|
762 | flag = xcalloc(1+n, sizeof(int)); |
---|
763 | work = xcalloc(1+n, sizeof(mpq_t)); |
---|
764 | for (k = 1; k <= n; k++) |
---|
765 | { flag[k] = 0; |
---|
766 | mpq_init(work[k]); |
---|
767 | } |
---|
768 | /* main elimination loop */ |
---|
769 | for (k = 1; k <= n; k++) |
---|
770 | { /* choose a pivot element v[p,q] */ |
---|
771 | piv = find_pivot(lux, wka); |
---|
772 | if (piv == NULL) |
---|
773 | { /* no pivot can be chosen, because the active submatrix is |
---|
774 | empty */ |
---|
775 | break; |
---|
776 | } |
---|
777 | /* determine row and column indices of the pivot element */ |
---|
778 | p = piv->i, q = piv->j; |
---|
779 | /* let v[p,q] correspond to u[i',j']; permute k-th and i'-th |
---|
780 | rows and k-th and j'-th columns of the matrix U = P*V*Q to |
---|
781 | move the element u[i',j'] to the position u[k,k] */ |
---|
782 | i = P_col[p], j = Q_row[q]; |
---|
783 | xassert(k <= i && i <= n && k <= j && j <= n); |
---|
784 | /* permute k-th and i-th rows of the matrix U */ |
---|
785 | t = P_row[k]; |
---|
786 | P_row[i] = t, P_col[t] = i; |
---|
787 | P_row[k] = p, P_col[p] = k; |
---|
788 | /* permute k-th and j-th columns of the matrix U */ |
---|
789 | t = Q_col[k]; |
---|
790 | Q_col[j] = t, Q_row[t] = j; |
---|
791 | Q_col[k] = q, Q_row[q] = k; |
---|
792 | /* eliminate subdiagonal elements of k-th column of the matrix |
---|
793 | U = P*V*Q using the pivot element u[k,k] = v[p,q] */ |
---|
794 | eliminate(lux, wka, piv, flag, work); |
---|
795 | } |
---|
796 | /* determine the rank of A (and V) */ |
---|
797 | lux->rank = k - 1; |
---|
798 | /* free working arrays */ |
---|
799 | xfree(flag); |
---|
800 | for (k = 1; k <= n; k++) mpq_clear(work[k]); |
---|
801 | xfree(work); |
---|
802 | /* build column lists of the matrix V using its row lists */ |
---|
803 | for (j = 1; j <= n; j++) |
---|
804 | xassert(V_col[j] == NULL); |
---|
805 | for (i = 1; i <= n; i++) |
---|
806 | { for (vij = V_row[i]; vij != NULL; vij = vij->r_next) |
---|
807 | { j = vij->j; |
---|
808 | vij->c_prev = NULL; |
---|
809 | vij->c_next = V_col[j]; |
---|
810 | if (vij->c_next != NULL) vij->c_next->c_prev = vij; |
---|
811 | V_col[j] = vij; |
---|
812 | } |
---|
813 | } |
---|
814 | /* free working area */ |
---|
815 | xfree(wka->R_len); |
---|
816 | xfree(wka->R_head); |
---|
817 | xfree(wka->R_prev); |
---|
818 | xfree(wka->R_next); |
---|
819 | xfree(wka->C_len); |
---|
820 | xfree(wka->C_head); |
---|
821 | xfree(wka->C_prev); |
---|
822 | xfree(wka->C_next); |
---|
823 | xfree(wka); |
---|
824 | /* return to the calling program */ |
---|
825 | return (lux->rank < n); |
---|
826 | } |
---|
827 | |
---|
828 | /*---------------------------------------------------------------------- |
---|
829 | // lux_f_solve - solve system F*x = b or F'*x = b. |
---|
830 | // |
---|
831 | // SYNOPSIS |
---|
832 | // |
---|
833 | // #include "glplux.h" |
---|
834 | // void lux_f_solve(LUX *lux, int tr, mpq_t x[]); |
---|
835 | // |
---|
836 | // DESCRIPTION |
---|
837 | // |
---|
838 | // The routine lux_f_solve solves either the system F*x = b (if the |
---|
839 | // flag tr is zero) or the system F'*x = b (if the flag tr is non-zero), |
---|
840 | // where the matrix F is a component of LU-factorization specified by |
---|
841 | // the parameter lux, F' is a matrix transposed to F. |
---|
842 | // |
---|
843 | // On entry the array x should contain elements of the right-hand side |
---|
844 | // vector b in locations x[1], ..., x[n], where n is the order of the |
---|
845 | // matrix F. On exit this array will contain elements of the solution |
---|
846 | // vector x in the same locations. */ |
---|
847 | |
---|
848 | void lux_f_solve(LUX *lux, int tr, mpq_t x[]) |
---|
849 | { int n = lux->n; |
---|
850 | LUXELM **F_row = lux->F_row; |
---|
851 | LUXELM **F_col = lux->F_col; |
---|
852 | int *P_row = lux->P_row; |
---|
853 | LUXELM *fik, *fkj; |
---|
854 | int i, j, k; |
---|
855 | mpq_t temp; |
---|
856 | mpq_init(temp); |
---|
857 | if (!tr) |
---|
858 | { /* solve the system F*x = b */ |
---|
859 | for (j = 1; j <= n; j++) |
---|
860 | { k = P_row[j]; |
---|
861 | if (mpq_sgn(x[k]) != 0) |
---|
862 | { for (fik = F_col[k]; fik != NULL; fik = fik->c_next) |
---|
863 | { mpq_mul(temp, fik->val, x[k]); |
---|
864 | mpq_sub(x[fik->i], x[fik->i], temp); |
---|
865 | } |
---|
866 | } |
---|
867 | } |
---|
868 | } |
---|
869 | else |
---|
870 | { /* solve the system F'*x = b */ |
---|
871 | for (i = n; i >= 1; i--) |
---|
872 | { k = P_row[i]; |
---|
873 | if (mpq_sgn(x[k]) != 0) |
---|
874 | { for (fkj = F_row[k]; fkj != NULL; fkj = fkj->r_next) |
---|
875 | { mpq_mul(temp, fkj->val, x[k]); |
---|
876 | mpq_sub(x[fkj->j], x[fkj->j], temp); |
---|
877 | } |
---|
878 | } |
---|
879 | } |
---|
880 | } |
---|
881 | mpq_clear(temp); |
---|
882 | return; |
---|
883 | } |
---|
884 | |
---|
885 | /*---------------------------------------------------------------------- |
---|
886 | // lux_v_solve - solve system V*x = b or V'*x = b. |
---|
887 | // |
---|
888 | // SYNOPSIS |
---|
889 | // |
---|
890 | // #include "glplux.h" |
---|
891 | // void lux_v_solve(LUX *lux, int tr, double x[]); |
---|
892 | // |
---|
893 | // DESCRIPTION |
---|
894 | // |
---|
895 | // The routine lux_v_solve solves either the system V*x = b (if the |
---|
896 | // flag tr is zero) or the system V'*x = b (if the flag tr is non-zero), |
---|
897 | // where the matrix V is a component of LU-factorization specified by |
---|
898 | // the parameter lux, V' is a matrix transposed to V. |
---|
899 | // |
---|
900 | // On entry the array x should contain elements of the right-hand side |
---|
901 | // vector b in locations x[1], ..., x[n], where n is the order of the |
---|
902 | // matrix V. On exit this array will contain elements of the solution |
---|
903 | // vector x in the same locations. */ |
---|
904 | |
---|
905 | void lux_v_solve(LUX *lux, int tr, mpq_t x[]) |
---|
906 | { int n = lux->n; |
---|
907 | mpq_t *V_piv = lux->V_piv; |
---|
908 | LUXELM **V_row = lux->V_row; |
---|
909 | LUXELM **V_col = lux->V_col; |
---|
910 | int *P_row = lux->P_row; |
---|
911 | int *Q_col = lux->Q_col; |
---|
912 | LUXELM *vij; |
---|
913 | int i, j, k; |
---|
914 | mpq_t *b, temp; |
---|
915 | b = xcalloc(1+n, sizeof(mpq_t)); |
---|
916 | for (k = 1; k <= n; k++) |
---|
917 | mpq_init(b[k]), mpq_set(b[k], x[k]), mpq_set_si(x[k], 0, 1); |
---|
918 | mpq_init(temp); |
---|
919 | if (!tr) |
---|
920 | { /* solve the system V*x = b */ |
---|
921 | for (k = n; k >= 1; k--) |
---|
922 | { i = P_row[k], j = Q_col[k]; |
---|
923 | if (mpq_sgn(b[i]) != 0) |
---|
924 | { mpq_set(x[j], b[i]); |
---|
925 | mpq_div(x[j], x[j], V_piv[i]); |
---|
926 | for (vij = V_col[j]; vij != NULL; vij = vij->c_next) |
---|
927 | { mpq_mul(temp, vij->val, x[j]); |
---|
928 | mpq_sub(b[vij->i], b[vij->i], temp); |
---|
929 | } |
---|
930 | } |
---|
931 | } |
---|
932 | } |
---|
933 | else |
---|
934 | { /* solve the system V'*x = b */ |
---|
935 | for (k = 1; k <= n; k++) |
---|
936 | { i = P_row[k], j = Q_col[k]; |
---|
937 | if (mpq_sgn(b[j]) != 0) |
---|
938 | { mpq_set(x[i], b[j]); |
---|
939 | mpq_div(x[i], x[i], V_piv[i]); |
---|
940 | for (vij = V_row[i]; vij != NULL; vij = vij->r_next) |
---|
941 | { mpq_mul(temp, vij->val, x[i]); |
---|
942 | mpq_sub(b[vij->j], b[vij->j], temp); |
---|
943 | } |
---|
944 | } |
---|
945 | } |
---|
946 | } |
---|
947 | for (k = 1; k <= n; k++) mpq_clear(b[k]); |
---|
948 | mpq_clear(temp); |
---|
949 | xfree(b); |
---|
950 | return; |
---|
951 | } |
---|
952 | |
---|
953 | /*---------------------------------------------------------------------- |
---|
954 | // lux_solve - solve system A*x = b or A'*x = b. |
---|
955 | // |
---|
956 | // SYNOPSIS |
---|
957 | // |
---|
958 | // #include "glplux.h" |
---|
959 | // void lux_solve(LUX *lux, int tr, mpq_t x[]); |
---|
960 | // |
---|
961 | // DESCRIPTION |
---|
962 | // |
---|
963 | // The routine lux_solve solves either the system A*x = b (if the flag |
---|
964 | // tr is zero) or the system A'*x = b (if the flag tr is non-zero), |
---|
965 | // where the parameter lux specifies LU-factorization of the matrix A, |
---|
966 | // A' is a matrix transposed to A. |
---|
967 | // |
---|
968 | // On entry the array x should contain elements of the right-hand side |
---|
969 | // vector b in locations x[1], ..., x[n], where n is the order of the |
---|
970 | // matrix A. On exit this array will contain elements of the solution |
---|
971 | // vector x in the same locations. */ |
---|
972 | |
---|
973 | void lux_solve(LUX *lux, int tr, mpq_t x[]) |
---|
974 | { if (lux->rank < lux->n) |
---|
975 | xfault("lux_solve: LU-factorization has incomplete rank\n"); |
---|
976 | if (!tr) |
---|
977 | { /* A = F*V, therefore inv(A) = inv(V)*inv(F) */ |
---|
978 | lux_f_solve(lux, 0, x); |
---|
979 | lux_v_solve(lux, 0, x); |
---|
980 | } |
---|
981 | else |
---|
982 | { /* A' = V'*F', therefore inv(A') = inv(F')*inv(V') */ |
---|
983 | lux_v_solve(lux, 1, x); |
---|
984 | lux_f_solve(lux, 1, x); |
---|
985 | } |
---|
986 | return; |
---|
987 | } |
---|
988 | |
---|
989 | /*---------------------------------------------------------------------- |
---|
990 | // lux_delete - delete LU-factorization. |
---|
991 | // |
---|
992 | // SYNOPSIS |
---|
993 | // |
---|
994 | // #include "glplux.h" |
---|
995 | // void lux_delete(LUX *lux); |
---|
996 | // |
---|
997 | // DESCRIPTION |
---|
998 | // |
---|
999 | // The routine lux_delete deletes LU-factorization data structure, |
---|
1000 | // which the parameter lux points to, freeing all the memory allocated |
---|
1001 | // to this object. */ |
---|
1002 | |
---|
1003 | void lux_delete(LUX *lux) |
---|
1004 | { int n = lux->n; |
---|
1005 | LUXELM *fij, *vij; |
---|
1006 | int i; |
---|
1007 | for (i = 1; i <= n; i++) |
---|
1008 | { for (fij = lux->F_row[i]; fij != NULL; fij = fij->r_next) |
---|
1009 | mpq_clear(fij->val); |
---|
1010 | mpq_clear(lux->V_piv[i]); |
---|
1011 | for (vij = lux->V_row[i]; vij != NULL; vij = vij->r_next) |
---|
1012 | mpq_clear(vij->val); |
---|
1013 | } |
---|
1014 | dmp_delete_pool(lux->pool); |
---|
1015 | xfree(lux->F_row); |
---|
1016 | xfree(lux->F_col); |
---|
1017 | xfree(lux->V_piv); |
---|
1018 | xfree(lux->V_row); |
---|
1019 | xfree(lux->V_col); |
---|
1020 | xfree(lux->P_row); |
---|
1021 | xfree(lux->P_col); |
---|
1022 | xfree(lux->Q_row); |
---|
1023 | xfree(lux->Q_col); |
---|
1024 | xfree(lux); |
---|
1025 | return; |
---|
1026 | } |
---|
1027 | |
---|
1028 | /* eof */ |
---|