|
1 /* glplpf.c (LP basis factorization, Schur complement version) */ |
|
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 "glplpf.h" |
|
26 #include "glpenv.h" |
|
27 #define xfault xerror |
|
28 |
|
29 #define _GLPLPF_DEBUG 0 |
|
30 |
|
31 /* CAUTION: DO NOT CHANGE THE LIMIT BELOW */ |
|
32 |
|
33 #define M_MAX 100000000 /* = 100*10^6 */ |
|
34 /* maximal order of the basis matrix */ |
|
35 |
|
36 /*********************************************************************** |
|
37 * NAME |
|
38 * |
|
39 * lpf_create_it - create LP basis factorization |
|
40 * |
|
41 * SYNOPSIS |
|
42 * |
|
43 * #include "glplpf.h" |
|
44 * LPF *lpf_create_it(void); |
|
45 * |
|
46 * DESCRIPTION |
|
47 * |
|
48 * The routine lpf_create_it creates a program object, which represents |
|
49 * a factorization of LP basis. |
|
50 * |
|
51 * RETURNS |
|
52 * |
|
53 * The routine lpf_create_it returns a pointer to the object created. */ |
|
54 |
|
55 LPF *lpf_create_it(void) |
|
56 { LPF *lpf; |
|
57 #if _GLPLPF_DEBUG |
|
58 xprintf("lpf_create_it: warning: debug mode enabled\n"); |
|
59 #endif |
|
60 lpf = xmalloc(sizeof(LPF)); |
|
61 lpf->valid = 0; |
|
62 lpf->m0_max = lpf->m0 = 0; |
|
63 lpf->luf = luf_create_it(); |
|
64 lpf->m = 0; |
|
65 lpf->B = NULL; |
|
66 lpf->n_max = 50; |
|
67 lpf->n = 0; |
|
68 lpf->R_ptr = lpf->R_len = NULL; |
|
69 lpf->S_ptr = lpf->S_len = NULL; |
|
70 lpf->scf = NULL; |
|
71 lpf->P_row = lpf->P_col = NULL; |
|
72 lpf->Q_row = lpf->Q_col = NULL; |
|
73 lpf->v_size = 1000; |
|
74 lpf->v_ptr = 0; |
|
75 lpf->v_ind = NULL; |
|
76 lpf->v_val = NULL; |
|
77 lpf->work1 = lpf->work2 = NULL; |
|
78 return lpf; |
|
79 } |
|
80 |
|
81 /*********************************************************************** |
|
82 * NAME |
|
83 * |
|
84 * lpf_factorize - compute LP basis factorization |
|
85 * |
|
86 * SYNOPSIS |
|
87 * |
|
88 * #include "glplpf.h" |
|
89 * int lpf_factorize(LPF *lpf, int m, const int bh[], int (*col) |
|
90 * (void *info, int j, int ind[], double val[]), void *info); |
|
91 * |
|
92 * DESCRIPTION |
|
93 * |
|
94 * The routine lpf_factorize computes the factorization of the basis |
|
95 * matrix B specified by the routine col. |
|
96 * |
|
97 * The parameter lpf specified the basis factorization data structure |
|
98 * created with the routine lpf_create_it. |
|
99 * |
|
100 * The parameter m specifies the order of B, m > 0. |
|
101 * |
|
102 * The array bh specifies the basis header: bh[j], 1 <= j <= m, is the |
|
103 * number of j-th column of B in some original matrix. The array bh is |
|
104 * optional and can be specified as NULL. |
|
105 * |
|
106 * The formal routine col specifies the matrix B to be factorized. To |
|
107 * obtain j-th column of A the routine lpf_factorize calls the routine |
|
108 * col with the parameter j (1 <= j <= n). In response the routine col |
|
109 * should store row indices and numerical values of non-zero elements |
|
110 * of j-th column of B to locations ind[1,...,len] and val[1,...,len], |
|
111 * respectively, where len is the number of non-zeros in j-th column |
|
112 * returned on exit. Neither zero nor duplicate elements are allowed. |
|
113 * |
|
114 * The parameter info is a transit pointer passed to the routine col. |
|
115 * |
|
116 * RETURNS |
|
117 * |
|
118 * 0 The factorization has been successfully computed. |
|
119 * |
|
120 * LPF_ESING |
|
121 * The specified matrix is singular within the working precision. |
|
122 * |
|
123 * LPF_ECOND |
|
124 * The specified matrix is ill-conditioned. |
|
125 * |
|
126 * For more details see comments to the routine luf_factorize. */ |
|
127 |
|
128 int lpf_factorize(LPF *lpf, int m, const int bh[], int (*col) |
|
129 (void *info, int j, int ind[], double val[]), void *info) |
|
130 { int k, ret; |
|
131 #if _GLPLPF_DEBUG |
|
132 int i, j, len, *ind; |
|
133 double *B, *val; |
|
134 #endif |
|
135 xassert(bh == bh); |
|
136 if (m < 1) |
|
137 xfault("lpf_factorize: m = %d; invalid parameter\n", m); |
|
138 if (m > M_MAX) |
|
139 xfault("lpf_factorize: m = %d; matrix too big\n", m); |
|
140 lpf->m0 = lpf->m = m; |
|
141 /* invalidate the factorization */ |
|
142 lpf->valid = 0; |
|
143 /* allocate/reallocate arrays, if necessary */ |
|
144 if (lpf->R_ptr == NULL) |
|
145 lpf->R_ptr = xcalloc(1+lpf->n_max, sizeof(int)); |
|
146 if (lpf->R_len == NULL) |
|
147 lpf->R_len = xcalloc(1+lpf->n_max, sizeof(int)); |
|
148 if (lpf->S_ptr == NULL) |
|
149 lpf->S_ptr = xcalloc(1+lpf->n_max, sizeof(int)); |
|
150 if (lpf->S_len == NULL) |
|
151 lpf->S_len = xcalloc(1+lpf->n_max, sizeof(int)); |
|
152 if (lpf->scf == NULL) |
|
153 lpf->scf = scf_create_it(lpf->n_max); |
|
154 if (lpf->v_ind == NULL) |
|
155 lpf->v_ind = xcalloc(1+lpf->v_size, sizeof(int)); |
|
156 if (lpf->v_val == NULL) |
|
157 lpf->v_val = xcalloc(1+lpf->v_size, sizeof(double)); |
|
158 if (lpf->m0_max < m) |
|
159 { if (lpf->P_row != NULL) xfree(lpf->P_row); |
|
160 if (lpf->P_col != NULL) xfree(lpf->P_col); |
|
161 if (lpf->Q_row != NULL) xfree(lpf->Q_row); |
|
162 if (lpf->Q_col != NULL) xfree(lpf->Q_col); |
|
163 if (lpf->work1 != NULL) xfree(lpf->work1); |
|
164 if (lpf->work2 != NULL) xfree(lpf->work2); |
|
165 lpf->m0_max = m + 100; |
|
166 lpf->P_row = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(int)); |
|
167 lpf->P_col = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(int)); |
|
168 lpf->Q_row = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(int)); |
|
169 lpf->Q_col = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(int)); |
|
170 lpf->work1 = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(double)); |
|
171 lpf->work2 = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(double)); |
|
172 } |
|
173 /* try to factorize the basis matrix */ |
|
174 switch (luf_factorize(lpf->luf, m, col, info)) |
|
175 { case 0: |
|
176 break; |
|
177 case LUF_ESING: |
|
178 ret = LPF_ESING; |
|
179 goto done; |
|
180 case LUF_ECOND: |
|
181 ret = LPF_ECOND; |
|
182 goto done; |
|
183 default: |
|
184 xassert(lpf != lpf); |
|
185 } |
|
186 /* the basis matrix has been successfully factorized */ |
|
187 lpf->valid = 1; |
|
188 #if _GLPLPF_DEBUG |
|
189 /* store the basis matrix for debugging */ |
|
190 if (lpf->B != NULL) xfree(lpf->B); |
|
191 xassert(m <= 32767); |
|
192 lpf->B = B = xcalloc(1+m*m, sizeof(double)); |
|
193 ind = xcalloc(1+m, sizeof(int)); |
|
194 val = xcalloc(1+m, sizeof(double)); |
|
195 for (k = 1; k <= m * m; k++) |
|
196 B[k] = 0.0; |
|
197 for (j = 1; j <= m; j++) |
|
198 { len = col(info, j, ind, val); |
|
199 xassert(0 <= len && len <= m); |
|
200 for (k = 1; k <= len; k++) |
|
201 { i = ind[k]; |
|
202 xassert(1 <= i && i <= m); |
|
203 xassert(B[(i - 1) * m + j] == 0.0); |
|
204 xassert(val[k] != 0.0); |
|
205 B[(i - 1) * m + j] = val[k]; |
|
206 } |
|
207 } |
|
208 xfree(ind); |
|
209 xfree(val); |
|
210 #endif |
|
211 /* B = B0, so there are no additional rows/columns */ |
|
212 lpf->n = 0; |
|
213 /* reset the Schur complement factorization */ |
|
214 scf_reset_it(lpf->scf); |
|
215 /* P := Q := I */ |
|
216 for (k = 1; k <= m; k++) |
|
217 { lpf->P_row[k] = lpf->P_col[k] = k; |
|
218 lpf->Q_row[k] = lpf->Q_col[k] = k; |
|
219 } |
|
220 /* make all SVA locations free */ |
|
221 lpf->v_ptr = 1; |
|
222 ret = 0; |
|
223 done: /* return to the calling program */ |
|
224 return ret; |
|
225 } |
|
226 |
|
227 /*********************************************************************** |
|
228 * The routine r_prod computes the product y := y + alpha * R * x, |
|
229 * where x is a n-vector, alpha is a scalar, y is a m0-vector. |
|
230 * |
|
231 * Since matrix R is available by columns, the product is computed as |
|
232 * a linear combination: |
|
233 * |
|
234 * y := y + alpha * (R[1] * x[1] + ... + R[n] * x[n]), |
|
235 * |
|
236 * where R[j] is j-th column of R. */ |
|
237 |
|
238 static void r_prod(LPF *lpf, double y[], double a, const double x[]) |
|
239 { int n = lpf->n; |
|
240 int *R_ptr = lpf->R_ptr; |
|
241 int *R_len = lpf->R_len; |
|
242 int *v_ind = lpf->v_ind; |
|
243 double *v_val = lpf->v_val; |
|
244 int j, beg, end, ptr; |
|
245 double t; |
|
246 for (j = 1; j <= n; j++) |
|
247 { if (x[j] == 0.0) continue; |
|
248 /* y := y + alpha * R[j] * x[j] */ |
|
249 t = a * x[j]; |
|
250 beg = R_ptr[j]; |
|
251 end = beg + R_len[j]; |
|
252 for (ptr = beg; ptr < end; ptr++) |
|
253 y[v_ind[ptr]] += t * v_val[ptr]; |
|
254 } |
|
255 return; |
|
256 } |
|
257 |
|
258 /*********************************************************************** |
|
259 * The routine rt_prod computes the product y := y + alpha * R' * x, |
|
260 * where R' is a matrix transposed to R, x is a m0-vector, alpha is a |
|
261 * scalar, y is a n-vector. |
|
262 * |
|
263 * Since matrix R is available by columns, the product components are |
|
264 * computed as inner products: |
|
265 * |
|
266 * y[j] := y[j] + alpha * (j-th column of R) * x |
|
267 * |
|
268 * for j = 1, 2, ..., n. */ |
|
269 |
|
270 static void rt_prod(LPF *lpf, double y[], double a, const double x[]) |
|
271 { int n = lpf->n; |
|
272 int *R_ptr = lpf->R_ptr; |
|
273 int *R_len = lpf->R_len; |
|
274 int *v_ind = lpf->v_ind; |
|
275 double *v_val = lpf->v_val; |
|
276 int j, beg, end, ptr; |
|
277 double t; |
|
278 for (j = 1; j <= n; j++) |
|
279 { /* t := (j-th column of R) * x */ |
|
280 t = 0.0; |
|
281 beg = R_ptr[j]; |
|
282 end = beg + R_len[j]; |
|
283 for (ptr = beg; ptr < end; ptr++) |
|
284 t += v_val[ptr] * x[v_ind[ptr]]; |
|
285 /* y[j] := y[j] + alpha * t */ |
|
286 y[j] += a * t; |
|
287 } |
|
288 return; |
|
289 } |
|
290 |
|
291 /*********************************************************************** |
|
292 * The routine s_prod computes the product y := y + alpha * S * x, |
|
293 * where x is a m0-vector, alpha is a scalar, y is a n-vector. |
|
294 * |
|
295 * Since matrix S is available by rows, the product components are |
|
296 * computed as inner products: |
|
297 * |
|
298 * y[i] = y[i] + alpha * (i-th row of S) * x |
|
299 * |
|
300 * for i = 1, 2, ..., n. */ |
|
301 |
|
302 static void s_prod(LPF *lpf, double y[], double a, const double x[]) |
|
303 { int n = lpf->n; |
|
304 int *S_ptr = lpf->S_ptr; |
|
305 int *S_len = lpf->S_len; |
|
306 int *v_ind = lpf->v_ind; |
|
307 double *v_val = lpf->v_val; |
|
308 int i, beg, end, ptr; |
|
309 double t; |
|
310 for (i = 1; i <= n; i++) |
|
311 { /* t := (i-th row of S) * x */ |
|
312 t = 0.0; |
|
313 beg = S_ptr[i]; |
|
314 end = beg + S_len[i]; |
|
315 for (ptr = beg; ptr < end; ptr++) |
|
316 t += v_val[ptr] * x[v_ind[ptr]]; |
|
317 /* y[i] := y[i] + alpha * t */ |
|
318 y[i] += a * t; |
|
319 } |
|
320 return; |
|
321 } |
|
322 |
|
323 /*********************************************************************** |
|
324 * The routine st_prod computes the product y := y + alpha * S' * x, |
|
325 * where S' is a matrix transposed to S, x is a n-vector, alpha is a |
|
326 * scalar, y is m0-vector. |
|
327 * |
|
328 * Since matrix R is available by rows, the product is computed as a |
|
329 * linear combination: |
|
330 * |
|
331 * y := y + alpha * (S'[1] * x[1] + ... + S'[n] * x[n]), |
|
332 * |
|
333 * where S'[i] is i-th row of S. */ |
|
334 |
|
335 static void st_prod(LPF *lpf, double y[], double a, const double x[]) |
|
336 { int n = lpf->n; |
|
337 int *S_ptr = lpf->S_ptr; |
|
338 int *S_len = lpf->S_len; |
|
339 int *v_ind = lpf->v_ind; |
|
340 double *v_val = lpf->v_val; |
|
341 int i, beg, end, ptr; |
|
342 double t; |
|
343 for (i = 1; i <= n; i++) |
|
344 { if (x[i] == 0.0) continue; |
|
345 /* y := y + alpha * S'[i] * x[i] */ |
|
346 t = a * x[i]; |
|
347 beg = S_ptr[i]; |
|
348 end = beg + S_len[i]; |
|
349 for (ptr = beg; ptr < end; ptr++) |
|
350 y[v_ind[ptr]] += t * v_val[ptr]; |
|
351 } |
|
352 return; |
|
353 } |
|
354 |
|
355 #if _GLPLPF_DEBUG |
|
356 /*********************************************************************** |
|
357 * The routine check_error computes the maximal relative error between |
|
358 * left- and right-hand sides for the system B * x = b (if tr is zero) |
|
359 * or B' * x = b (if tr is non-zero), where B' is a matrix transposed |
|
360 * to B. (This routine is intended for debugging only.) */ |
|
361 |
|
362 static void check_error(LPF *lpf, int tr, const double x[], |
|
363 const double b[]) |
|
364 { int m = lpf->m; |
|
365 double *B = lpf->B; |
|
366 int i, j; |
|
367 double d, dmax = 0.0, s, t, tmax; |
|
368 for (i = 1; i <= m; i++) |
|
369 { s = 0.0; |
|
370 tmax = 1.0; |
|
371 for (j = 1; j <= m; j++) |
|
372 { if (!tr) |
|
373 t = B[m * (i - 1) + j] * x[j]; |
|
374 else |
|
375 t = B[m * (j - 1) + i] * x[j]; |
|
376 if (tmax < fabs(t)) tmax = fabs(t); |
|
377 s += t; |
|
378 } |
|
379 d = fabs(s - b[i]) / tmax; |
|
380 if (dmax < d) dmax = d; |
|
381 } |
|
382 if (dmax > 1e-8) |
|
383 xprintf("%s: dmax = %g; relative error too large\n", |
|
384 !tr ? "lpf_ftran" : "lpf_btran", dmax); |
|
385 return; |
|
386 } |
|
387 #endif |
|
388 |
|
389 /*********************************************************************** |
|
390 * NAME |
|
391 * |
|
392 * lpf_ftran - perform forward transformation (solve system B*x = b) |
|
393 * |
|
394 * SYNOPSIS |
|
395 * |
|
396 * #include "glplpf.h" |
|
397 * void lpf_ftran(LPF *lpf, double x[]); |
|
398 * |
|
399 * DESCRIPTION |
|
400 * |
|
401 * The routine lpf_ftran performs forward transformation, i.e. solves |
|
402 * the system B*x = b, where B is the basis matrix, x is the vector of |
|
403 * unknowns to be computed, b is the vector of right-hand sides. |
|
404 * |
|
405 * On entry elements of the vector b should be stored in dense format |
|
406 * in locations x[1], ..., x[m], where m is the number of rows. On exit |
|
407 * the routine stores elements of the vector x in the same locations. |
|
408 * |
|
409 * BACKGROUND |
|
410 * |
|
411 * Solution of the system B * x = b can be obtained by solving the |
|
412 * following augmented system: |
|
413 * |
|
414 * ( B F^) ( x ) ( b ) |
|
415 * ( ) ( ) = ( ) |
|
416 * ( G^ H^) ( y ) ( 0 ) |
|
417 * |
|
418 * which, using the main equality, can be written as follows: |
|
419 * |
|
420 * ( L0 0 ) ( U0 R ) ( x ) ( b ) |
|
421 * P ( ) ( ) Q ( ) = ( ) |
|
422 * ( S I ) ( 0 C ) ( y ) ( 0 ) |
|
423 * |
|
424 * therefore, |
|
425 * |
|
426 * ( x ) ( U0 R )-1 ( L0 0 )-1 ( b ) |
|
427 * ( ) = Q' ( ) ( ) P' ( ) |
|
428 * ( y ) ( 0 C ) ( S I ) ( 0 ) |
|
429 * |
|
430 * Thus, computing the solution includes the following steps: |
|
431 * |
|
432 * 1. Compute |
|
433 * |
|
434 * ( f ) ( b ) |
|
435 * ( ) = P' ( ) |
|
436 * ( g ) ( 0 ) |
|
437 * |
|
438 * 2. Solve the system |
|
439 * |
|
440 * ( f1 ) ( L0 0 )-1 ( f ) ( L0 0 ) ( f1 ) ( f ) |
|
441 * ( ) = ( ) ( ) => ( ) ( ) = ( ) |
|
442 * ( g1 ) ( S I ) ( g ) ( S I ) ( g1 ) ( g ) |
|
443 * |
|
444 * from which it follows that: |
|
445 * |
|
446 * { L0 * f1 = f f1 = inv(L0) * f |
|
447 * { => |
|
448 * { S * f1 + g1 = g g1 = g - S * f1 |
|
449 * |
|
450 * 3. Solve the system |
|
451 * |
|
452 * ( f2 ) ( U0 R )-1 ( f1 ) ( U0 R ) ( f2 ) ( f1 ) |
|
453 * ( ) = ( ) ( ) => ( ) ( ) = ( ) |
|
454 * ( g2 ) ( 0 C ) ( g1 ) ( 0 C ) ( g2 ) ( g1 ) |
|
455 * |
|
456 * from which it follows that: |
|
457 * |
|
458 * { U0 * f2 + R * g2 = f1 f2 = inv(U0) * (f1 - R * g2) |
|
459 * { => |
|
460 * { C * g2 = g1 g2 = inv(C) * g1 |
|
461 * |
|
462 * 4. Compute |
|
463 * |
|
464 * ( x ) ( f2 ) |
|
465 * ( ) = Q' ( ) |
|
466 * ( y ) ( g2 ) */ |
|
467 |
|
468 void lpf_ftran(LPF *lpf, double x[]) |
|
469 { int m0 = lpf->m0; |
|
470 int m = lpf->m; |
|
471 int n = lpf->n; |
|
472 int *P_col = lpf->P_col; |
|
473 int *Q_col = lpf->Q_col; |
|
474 double *fg = lpf->work1; |
|
475 double *f = fg; |
|
476 double *g = fg + m0; |
|
477 int i, ii; |
|
478 #if _GLPLPF_DEBUG |
|
479 double *b; |
|
480 #endif |
|
481 if (!lpf->valid) |
|
482 xfault("lpf_ftran: the factorization is not valid\n"); |
|
483 xassert(0 <= m && m <= m0 + n); |
|
484 #if _GLPLPF_DEBUG |
|
485 /* save the right-hand side vector */ |
|
486 b = xcalloc(1+m, sizeof(double)); |
|
487 for (i = 1; i <= m; i++) b[i] = x[i]; |
|
488 #endif |
|
489 /* (f g) := inv(P) * (b 0) */ |
|
490 for (i = 1; i <= m0 + n; i++) |
|
491 fg[i] = ((ii = P_col[i]) <= m ? x[ii] : 0.0); |
|
492 /* f1 := inv(L0) * f */ |
|
493 luf_f_solve(lpf->luf, 0, f); |
|
494 /* g1 := g - S * f1 */ |
|
495 s_prod(lpf, g, -1.0, f); |
|
496 /* g2 := inv(C) * g1 */ |
|
497 scf_solve_it(lpf->scf, 0, g); |
|
498 /* f2 := inv(U0) * (f1 - R * g2) */ |
|
499 r_prod(lpf, f, -1.0, g); |
|
500 luf_v_solve(lpf->luf, 0, f); |
|
501 /* (x y) := inv(Q) * (f2 g2) */ |
|
502 for (i = 1; i <= m; i++) |
|
503 x[i] = fg[Q_col[i]]; |
|
504 #if _GLPLPF_DEBUG |
|
505 /* check relative error in solution */ |
|
506 check_error(lpf, 0, x, b); |
|
507 xfree(b); |
|
508 #endif |
|
509 return; |
|
510 } |
|
511 |
|
512 /*********************************************************************** |
|
513 * NAME |
|
514 * |
|
515 * lpf_btran - perform backward transformation (solve system B'*x = b) |
|
516 * |
|
517 * SYNOPSIS |
|
518 * |
|
519 * #include "glplpf.h" |
|
520 * void lpf_btran(LPF *lpf, double x[]); |
|
521 * |
|
522 * DESCRIPTION |
|
523 * |
|
524 * The routine lpf_btran performs backward transformation, i.e. solves |
|
525 * the system B'*x = b, where B' is a matrix transposed to the basis |
|
526 * matrix B, x is the vector of unknowns to be computed, b is the vector |
|
527 * of right-hand sides. |
|
528 * |
|
529 * On entry elements of the vector b should be stored in dense format |
|
530 * in locations x[1], ..., x[m], where m is the number of rows. On exit |
|
531 * the routine stores elements of the vector x in the same locations. |
|
532 * |
|
533 * BACKGROUND |
|
534 * |
|
535 * Solution of the system B' * x = b, where B' is a matrix transposed |
|
536 * to B, can be obtained by solving the following augmented system: |
|
537 * |
|
538 * ( B F^)T ( x ) ( b ) |
|
539 * ( ) ( ) = ( ) |
|
540 * ( G^ H^) ( y ) ( 0 ) |
|
541 * |
|
542 * which, using the main equality, can be written as follows: |
|
543 * |
|
544 * T ( U0 R )T ( L0 0 )T T ( x ) ( b ) |
|
545 * Q ( ) ( ) P ( ) = ( ) |
|
546 * ( 0 C ) ( S I ) ( y ) ( 0 ) |
|
547 * |
|
548 * or, equivalently, as follows: |
|
549 * |
|
550 * ( U'0 0 ) ( L'0 S') ( x ) ( b ) |
|
551 * Q' ( ) ( ) P' ( ) = ( ) |
|
552 * ( R' C') ( 0 I ) ( y ) ( 0 ) |
|
553 * |
|
554 * therefore, |
|
555 * |
|
556 * ( x ) ( L'0 S')-1 ( U'0 0 )-1 ( b ) |
|
557 * ( ) = P ( ) ( ) Q ( ) |
|
558 * ( y ) ( 0 I ) ( R' C') ( 0 ) |
|
559 * |
|
560 * Thus, computing the solution includes the following steps: |
|
561 * |
|
562 * 1. Compute |
|
563 * |
|
564 * ( f ) ( b ) |
|
565 * ( ) = Q ( ) |
|
566 * ( g ) ( 0 ) |
|
567 * |
|
568 * 2. Solve the system |
|
569 * |
|
570 * ( f1 ) ( U'0 0 )-1 ( f ) ( U'0 0 ) ( f1 ) ( f ) |
|
571 * ( ) = ( ) ( ) => ( ) ( ) = ( ) |
|
572 * ( g1 ) ( R' C') ( g ) ( R' C') ( g1 ) ( g ) |
|
573 * |
|
574 * from which it follows that: |
|
575 * |
|
576 * { U'0 * f1 = f f1 = inv(U'0) * f |
|
577 * { => |
|
578 * { R' * f1 + C' * g1 = g g1 = inv(C') * (g - R' * f1) |
|
579 * |
|
580 * 3. Solve the system |
|
581 * |
|
582 * ( f2 ) ( L'0 S')-1 ( f1 ) ( L'0 S') ( f2 ) ( f1 ) |
|
583 * ( ) = ( ) ( ) => ( ) ( ) = ( ) |
|
584 * ( g2 ) ( 0 I ) ( g1 ) ( 0 I ) ( g2 ) ( g1 ) |
|
585 * |
|
586 * from which it follows that: |
|
587 * |
|
588 * { L'0 * f2 + S' * g2 = f1 |
|
589 * { => f2 = inv(L'0) * ( f1 - S' * g2) |
|
590 * { g2 = g1 |
|
591 * |
|
592 * 4. Compute |
|
593 * |
|
594 * ( x ) ( f2 ) |
|
595 * ( ) = P ( ) |
|
596 * ( y ) ( g2 ) */ |
|
597 |
|
598 void lpf_btran(LPF *lpf, double x[]) |
|
599 { int m0 = lpf->m0; |
|
600 int m = lpf->m; |
|
601 int n = lpf->n; |
|
602 int *P_row = lpf->P_row; |
|
603 int *Q_row = lpf->Q_row; |
|
604 double *fg = lpf->work1; |
|
605 double *f = fg; |
|
606 double *g = fg + m0; |
|
607 int i, ii; |
|
608 #if _GLPLPF_DEBUG |
|
609 double *b; |
|
610 #endif |
|
611 if (!lpf->valid) |
|
612 xfault("lpf_btran: the factorization is not valid\n"); |
|
613 xassert(0 <= m && m <= m0 + n); |
|
614 #if _GLPLPF_DEBUG |
|
615 /* save the right-hand side vector */ |
|
616 b = xcalloc(1+m, sizeof(double)); |
|
617 for (i = 1; i <= m; i++) b[i] = x[i]; |
|
618 #endif |
|
619 /* (f g) := Q * (b 0) */ |
|
620 for (i = 1; i <= m0 + n; i++) |
|
621 fg[i] = ((ii = Q_row[i]) <= m ? x[ii] : 0.0); |
|
622 /* f1 := inv(U'0) * f */ |
|
623 luf_v_solve(lpf->luf, 1, f); |
|
624 /* g1 := inv(C') * (g - R' * f1) */ |
|
625 rt_prod(lpf, g, -1.0, f); |
|
626 scf_solve_it(lpf->scf, 1, g); |
|
627 /* g2 := g1 */ |
|
628 g = g; |
|
629 /* f2 := inv(L'0) * (f1 - S' * g2) */ |
|
630 st_prod(lpf, f, -1.0, g); |
|
631 luf_f_solve(lpf->luf, 1, f); |
|
632 /* (x y) := P * (f2 g2) */ |
|
633 for (i = 1; i <= m; i++) |
|
634 x[i] = fg[P_row[i]]; |
|
635 #if _GLPLPF_DEBUG |
|
636 /* check relative error in solution */ |
|
637 check_error(lpf, 1, x, b); |
|
638 xfree(b); |
|
639 #endif |
|
640 return; |
|
641 } |
|
642 |
|
643 /*********************************************************************** |
|
644 * The routine enlarge_sva enlarges the Sparse Vector Area to new_size |
|
645 * locations by reallocating the arrays v_ind and v_val. */ |
|
646 |
|
647 static void enlarge_sva(LPF *lpf, int new_size) |
|
648 { int v_size = lpf->v_size; |
|
649 int used = lpf->v_ptr - 1; |
|
650 int *v_ind = lpf->v_ind; |
|
651 double *v_val = lpf->v_val; |
|
652 xassert(v_size < new_size); |
|
653 while (v_size < new_size) v_size += v_size; |
|
654 lpf->v_size = v_size; |
|
655 lpf->v_ind = xcalloc(1+v_size, sizeof(int)); |
|
656 lpf->v_val = xcalloc(1+v_size, sizeof(double)); |
|
657 xassert(used >= 0); |
|
658 memcpy(&lpf->v_ind[1], &v_ind[1], used * sizeof(int)); |
|
659 memcpy(&lpf->v_val[1], &v_val[1], used * sizeof(double)); |
|
660 xfree(v_ind); |
|
661 xfree(v_val); |
|
662 return; |
|
663 } |
|
664 |
|
665 /*********************************************************************** |
|
666 * NAME |
|
667 * |
|
668 * lpf_update_it - update LP basis factorization |
|
669 * |
|
670 * SYNOPSIS |
|
671 * |
|
672 * #include "glplpf.h" |
|
673 * int lpf_update_it(LPF *lpf, int j, int bh, int len, const int ind[], |
|
674 * const double val[]); |
|
675 * |
|
676 * DESCRIPTION |
|
677 * |
|
678 * The routine lpf_update_it updates the factorization of the basis |
|
679 * matrix B after replacing its j-th column by a new vector. |
|
680 * |
|
681 * The parameter j specifies the number of column of B, which has been |
|
682 * replaced, 1 <= j <= m, where m is the order of B. |
|
683 * |
|
684 * The parameter bh specifies the basis header entry for the new column |
|
685 * of B, which is the number of the new column in some original matrix. |
|
686 * This parameter is optional and can be specified as 0. |
|
687 * |
|
688 * Row indices and numerical values of non-zero elements of the new |
|
689 * column of B should be placed in locations ind[1], ..., ind[len] and |
|
690 * val[1], ..., val[len], resp., where len is the number of non-zeros |
|
691 * in the column. Neither zero nor duplicate elements are allowed. |
|
692 * |
|
693 * RETURNS |
|
694 * |
|
695 * 0 The factorization has been successfully updated. |
|
696 * |
|
697 * LPF_ESING |
|
698 * New basis B is singular within the working precision. |
|
699 * |
|
700 * LPF_ELIMIT |
|
701 * Maximal number of additional rows and columns has been reached. |
|
702 * |
|
703 * BACKGROUND |
|
704 * |
|
705 * Let j-th column of the current basis matrix B have to be replaced by |
|
706 * a new column a. This replacement is equivalent to removing the old |
|
707 * j-th column by fixing it at zero and introducing the new column as |
|
708 * follows: |
|
709 * |
|
710 * ( B F^| a ) |
|
711 * ( B F^) ( | ) |
|
712 * ( ) ---> ( G^ H^| 0 ) |
|
713 * ( G^ H^) (-------+---) |
|
714 * ( e'j 0 | 0 ) |
|
715 * |
|
716 * where ej is a unit vector with 1 in j-th position which used to fix |
|
717 * the old j-th column of B (at zero). Then using the main equality we |
|
718 * have: |
|
719 * |
|
720 * ( B F^| a ) ( B0 F | f ) |
|
721 * ( | ) ( P 0 ) ( | ) ( Q 0 ) |
|
722 * ( G^ H^| 0 ) = ( ) ( G H | g ) ( ) = |
|
723 * (-------+---) ( 0 1 ) (-------+---) ( 0 1 ) |
|
724 * ( e'j 0 | 0 ) ( v' w'| 0 ) |
|
725 * |
|
726 * [ ( B0 F )| ( f ) ] [ ( B0 F ) | ( f ) ] |
|
727 * [ P ( )| P ( ) ] ( Q 0 ) [ P ( ) Q| P ( ) ] |
|
728 * = [ ( G H )| ( g ) ] ( ) = [ ( G H ) | ( g ) ] |
|
729 * [------------+-------- ] ( 0 1 ) [-------------+---------] |
|
730 * [ ( v' w')| 0 ] [ ( v' w') Q| 0 ] |
|
731 * |
|
732 * where: |
|
733 * |
|
734 * ( a ) ( f ) ( f ) ( a ) |
|
735 * ( ) = P ( ) => ( ) = P' * ( ) |
|
736 * ( 0 ) ( g ) ( g ) ( 0 ) |
|
737 * |
|
738 * ( ej ) ( v ) ( v ) ( ej ) |
|
739 * ( e'j 0 ) = ( v' w' ) Q => ( ) = Q' ( ) => ( ) = Q ( ) |
|
740 * ( 0 ) ( w ) ( w ) ( 0 ) |
|
741 * |
|
742 * On the other hand: |
|
743 * |
|
744 * ( B0| F f ) |
|
745 * ( P 0 ) (---+------) ( Q 0 ) ( B0 new F ) |
|
746 * ( ) ( G | H g ) ( ) = new P ( ) new Q |
|
747 * ( 0 1 ) ( | ) ( 0 1 ) ( new G new H ) |
|
748 * ( v'| w' 0 ) |
|
749 * |
|
750 * where: |
|
751 * ( G ) ( H g ) |
|
752 * new F = ( F f ), new G = ( ), new H = ( ), |
|
753 * ( v') ( w' 0 ) |
|
754 * |
|
755 * ( P 0 ) ( Q 0 ) |
|
756 * new P = ( ) , new Q = ( ) . |
|
757 * ( 0 1 ) ( 0 1 ) |
|
758 * |
|
759 * The factorization structure for the new augmented matrix remains the |
|
760 * same, therefore: |
|
761 * |
|
762 * ( B0 new F ) ( L0 0 ) ( U0 new R ) |
|
763 * new P ( ) new Q = ( ) ( ) |
|
764 * ( new G new H ) ( new S I ) ( 0 new C ) |
|
765 * |
|
766 * where: |
|
767 * |
|
768 * new F = L0 * new R => |
|
769 * |
|
770 * new R = inv(L0) * new F = inv(L0) * (F f) = ( R inv(L0)*f ) |
|
771 * |
|
772 * new G = new S * U0 => |
|
773 * |
|
774 * ( G ) ( S ) |
|
775 * new S = new G * inv(U0) = ( ) * inv(U0) = ( ) |
|
776 * ( v') ( v'*inv(U0) ) |
|
777 * |
|
778 * new H = new S * new R + new C => |
|
779 * |
|
780 * new C = new H - new S * new R = |
|
781 * |
|
782 * ( H g ) ( S ) |
|
783 * = ( ) - ( ) * ( R inv(L0)*f ) = |
|
784 * ( w' 0 ) ( v'*inv(U0) ) |
|
785 * |
|
786 * ( H - S*R g - S*inv(L0)*f ) ( C x ) |
|
787 * = ( ) = ( ) |
|
788 * ( w'- v'*inv(U0)*R -v'*inv(U0)*inv(L0)*f) ( y' z ) |
|
789 * |
|
790 * Note that new C is resulted by expanding old C with new column x, |
|
791 * row y', and diagonal element z, where: |
|
792 * |
|
793 * x = g - S * inv(L0) * f = g - S * (new column of R) |
|
794 * |
|
795 * y = w - R'* inv(U'0)* v = w - R'* (new row of S) |
|
796 * |
|
797 * z = - (new row of S) * (new column of R) |
|
798 * |
|
799 * Finally, to replace old B by new B we have to permute j-th and last |
|
800 * (just added) columns of the matrix |
|
801 * |
|
802 * ( B F^| a ) |
|
803 * ( | ) |
|
804 * ( G^ H^| 0 ) |
|
805 * (-------+---) |
|
806 * ( e'j 0 | 0 ) |
|
807 * |
|
808 * and to keep the main equality do the same for matrix Q. */ |
|
809 |
|
810 int lpf_update_it(LPF *lpf, int j, int bh, int len, const int ind[], |
|
811 const double val[]) |
|
812 { int m0 = lpf->m0; |
|
813 int m = lpf->m; |
|
814 #if _GLPLPF_DEBUG |
|
815 double *B = lpf->B; |
|
816 #endif |
|
817 int n = lpf->n; |
|
818 int *R_ptr = lpf->R_ptr; |
|
819 int *R_len = lpf->R_len; |
|
820 int *S_ptr = lpf->S_ptr; |
|
821 int *S_len = lpf->S_len; |
|
822 int *P_row = lpf->P_row; |
|
823 int *P_col = lpf->P_col; |
|
824 int *Q_row = lpf->Q_row; |
|
825 int *Q_col = lpf->Q_col; |
|
826 int v_ptr = lpf->v_ptr; |
|
827 int *v_ind = lpf->v_ind; |
|
828 double *v_val = lpf->v_val; |
|
829 double *a = lpf->work2; /* new column */ |
|
830 double *fg = lpf->work1, *f = fg, *g = fg + m0; |
|
831 double *vw = lpf->work2, *v = vw, *w = vw + m0; |
|
832 double *x = g, *y = w, z; |
|
833 int i, ii, k, ret; |
|
834 xassert(bh == bh); |
|
835 if (!lpf->valid) |
|
836 xfault("lpf_update_it: the factorization is not valid\n"); |
|
837 if (!(1 <= j && j <= m)) |
|
838 xfault("lpf_update_it: j = %d; column number out of range\n", |
|
839 j); |
|
840 xassert(0 <= m && m <= m0 + n); |
|
841 /* check if the basis factorization can be expanded */ |
|
842 if (n == lpf->n_max) |
|
843 { lpf->valid = 0; |
|
844 ret = LPF_ELIMIT; |
|
845 goto done; |
|
846 } |
|
847 /* convert new j-th column of B to dense format */ |
|
848 for (i = 1; i <= m; i++) |
|
849 a[i] = 0.0; |
|
850 for (k = 1; k <= len; k++) |
|
851 { i = ind[k]; |
|
852 if (!(1 <= i && i <= m)) |
|
853 xfault("lpf_update_it: ind[%d] = %d; row number out of rang" |
|
854 "e\n", k, i); |
|
855 if (a[i] != 0.0) |
|
856 xfault("lpf_update_it: ind[%d] = %d; duplicate row index no" |
|
857 "t allowed\n", k, i); |
|
858 if (val[k] == 0.0) |
|
859 xfault("lpf_update_it: val[%d] = %g; zero element not allow" |
|
860 "ed\n", k, val[k]); |
|
861 a[i] = val[k]; |
|
862 } |
|
863 #if _GLPLPF_DEBUG |
|
864 /* change column in the basis matrix for debugging */ |
|
865 for (i = 1; i <= m; i++) |
|
866 B[(i - 1) * m + j] = a[i]; |
|
867 #endif |
|
868 /* (f g) := inv(P) * (a 0) */ |
|
869 for (i = 1; i <= m0+n; i++) |
|
870 fg[i] = ((ii = P_col[i]) <= m ? a[ii] : 0.0); |
|
871 /* (v w) := Q * (ej 0) */ |
|
872 for (i = 1; i <= m0+n; i++) vw[i] = 0.0; |
|
873 vw[Q_col[j]] = 1.0; |
|
874 /* f1 := inv(L0) * f (new column of R) */ |
|
875 luf_f_solve(lpf->luf, 0, f); |
|
876 /* v1 := inv(U'0) * v (new row of S) */ |
|
877 luf_v_solve(lpf->luf, 1, v); |
|
878 /* we need at most 2 * m0 available locations in the SVA to store |
|
879 new column of matrix R and new row of matrix S */ |
|
880 if (lpf->v_size < v_ptr + m0 + m0) |
|
881 { enlarge_sva(lpf, v_ptr + m0 + m0); |
|
882 v_ind = lpf->v_ind; |
|
883 v_val = lpf->v_val; |
|
884 } |
|
885 /* store new column of R */ |
|
886 R_ptr[n+1] = v_ptr; |
|
887 for (i = 1; i <= m0; i++) |
|
888 { if (f[i] != 0.0) |
|
889 v_ind[v_ptr] = i, v_val[v_ptr] = f[i], v_ptr++; |
|
890 } |
|
891 R_len[n+1] = v_ptr - lpf->v_ptr; |
|
892 lpf->v_ptr = v_ptr; |
|
893 /* store new row of S */ |
|
894 S_ptr[n+1] = v_ptr; |
|
895 for (i = 1; i <= m0; i++) |
|
896 { if (v[i] != 0.0) |
|
897 v_ind[v_ptr] = i, v_val[v_ptr] = v[i], v_ptr++; |
|
898 } |
|
899 S_len[n+1] = v_ptr - lpf->v_ptr; |
|
900 lpf->v_ptr = v_ptr; |
|
901 /* x := g - S * f1 (new column of C) */ |
|
902 s_prod(lpf, x, -1.0, f); |
|
903 /* y := w - R' * v1 (new row of C) */ |
|
904 rt_prod(lpf, y, -1.0, v); |
|
905 /* z := - v1 * f1 (new diagonal element of C) */ |
|
906 z = 0.0; |
|
907 for (i = 1; i <= m0; i++) z -= v[i] * f[i]; |
|
908 /* update factorization of new matrix C */ |
|
909 switch (scf_update_exp(lpf->scf, x, y, z)) |
|
910 { case 0: |
|
911 break; |
|
912 case SCF_ESING: |
|
913 lpf->valid = 0; |
|
914 ret = LPF_ESING; |
|
915 goto done; |
|
916 case SCF_ELIMIT: |
|
917 xassert(lpf != lpf); |
|
918 default: |
|
919 xassert(lpf != lpf); |
|
920 } |
|
921 /* expand matrix P */ |
|
922 P_row[m0+n+1] = P_col[m0+n+1] = m0+n+1; |
|
923 /* expand matrix Q */ |
|
924 Q_row[m0+n+1] = Q_col[m0+n+1] = m0+n+1; |
|
925 /* permute j-th and last (just added) column of matrix Q */ |
|
926 i = Q_col[j], ii = Q_col[m0+n+1]; |
|
927 Q_row[i] = m0+n+1, Q_col[m0+n+1] = i; |
|
928 Q_row[ii] = j, Q_col[j] = ii; |
|
929 /* increase the number of additional rows and columns */ |
|
930 lpf->n++; |
|
931 xassert(lpf->n <= lpf->n_max); |
|
932 /* the factorization has been successfully updated */ |
|
933 ret = 0; |
|
934 done: /* return to the calling program */ |
|
935 return ret; |
|
936 } |
|
937 |
|
938 /*********************************************************************** |
|
939 * NAME |
|
940 * |
|
941 * lpf_delete_it - delete LP basis factorization |
|
942 * |
|
943 * SYNOPSIS |
|
944 * |
|
945 * #include "glplpf.h" |
|
946 * void lpf_delete_it(LPF *lpf) |
|
947 * |
|
948 * DESCRIPTION |
|
949 * |
|
950 * The routine lpf_delete_it deletes LP basis factorization specified |
|
951 * by the parameter lpf and frees all memory allocated to this program |
|
952 * object. */ |
|
953 |
|
954 void lpf_delete_it(LPF *lpf) |
|
955 { luf_delete_it(lpf->luf); |
|
956 #if _GLPLPF_DEBUG |
|
957 if (lpf->B != NULL) xfree(lpf->B); |
|
958 #else |
|
959 xassert(lpf->B == NULL); |
|
960 #endif |
|
961 if (lpf->R_ptr != NULL) xfree(lpf->R_ptr); |
|
962 if (lpf->R_len != NULL) xfree(lpf->R_len); |
|
963 if (lpf->S_ptr != NULL) xfree(lpf->S_ptr); |
|
964 if (lpf->S_len != NULL) xfree(lpf->S_len); |
|
965 if (lpf->scf != NULL) scf_delete_it(lpf->scf); |
|
966 if (lpf->P_row != NULL) xfree(lpf->P_row); |
|
967 if (lpf->P_col != NULL) xfree(lpf->P_col); |
|
968 if (lpf->Q_row != NULL) xfree(lpf->Q_row); |
|
969 if (lpf->Q_col != NULL) xfree(lpf->Q_col); |
|
970 if (lpf->v_ind != NULL) xfree(lpf->v_ind); |
|
971 if (lpf->v_val != NULL) xfree(lpf->v_val); |
|
972 if (lpf->work1 != NULL) xfree(lpf->work1); |
|
973 if (lpf->work2 != NULL) xfree(lpf->work2); |
|
974 xfree(lpf); |
|
975 return; |
|
976 } |
|
977 |
|
978 /* eof */ |