lemon-project-template-glpk
comparison deps/glpk/src/glplpf.c @ 11:4fc6ad2fb8a6
Test GLPK in src/main.cc
author | Alpar Juttner <alpar@cs.elte.hu> |
---|---|
date | Sun, 06 Nov 2011 21:43:29 +0100 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:2dc206db8947 |
---|---|
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, 2011 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 */ |