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
lpf_create_it(void)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
lpf_factorize(LPF * lpf,int m,const int bh[],int (* col)(void * info,int j,int ind[],double val[]),void * info)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
r_prod(LPF * lpf,double y[],double a,const double x[])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
rt_prod(LPF * lpf,double y[],double a,const double x[])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
s_prod(LPF * lpf,double y[],double a,const double x[])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
st_prod(LPF * lpf,double y[],double a,const double x[])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
check_error(LPF * lpf,int tr,const double x[],const double b[])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
lpf_ftran(LPF * lpf,double x[])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
lpf_btran(LPF * lpf,double x[])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
enlarge_sva(LPF * lpf,int new_size)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
lpf_update_it(LPF * lpf,int j,int bh,int len,const int ind[],const double val[])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
lpf_delete_it(LPF * lpf)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 */
979