lemon-project-template-glpk

annotate deps/glpk/src/glplpf.c @ 11:4fc6ad2fb8a6

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