lemon-project-template-glpk

annotate deps/glpk/src/glpbfd.c @ 9:33de93886c88

find changesets by author, revision, files, or words in the commit message
Import GLPK 4.47
author Alpar Juttner <alpar@cs.elte.hu>
date Sun, 06 Nov 2011 20:59:10 +0100
parents
children
rev   line source
alpar@9 1 /* glpbfd.c (LP basis factorization driver) */
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 typedef struct BFD BFD;
alpar@9 26
alpar@9 27 #define GLPBFD_PRIVATE
alpar@9 28 #include "glpapi.h"
alpar@9 29 #include "glpfhv.h"
alpar@9 30 #include "glplpf.h"
alpar@9 31
alpar@9 32 /* CAUTION: DO NOT CHANGE THE LIMIT BELOW */
alpar@9 33
alpar@9 34 #define M_MAX 100000000 /* = 100*10^6 */
alpar@9 35 /* maximal order of the basis matrix */
alpar@9 36
alpar@9 37 struct BFD
alpar@9 38 { /* LP basis factorization */
alpar@9 39 int valid;
alpar@9 40 /* factorization is valid only if this flag is set */
alpar@9 41 int type;
alpar@9 42 /* factorization type:
alpar@9 43 GLP_BF_FT - LUF + Forrest-Tomlin
alpar@9 44 GLP_BF_BG - LUF + Schur compl. + Bartels-Golub
alpar@9 45 GLP_BF_GR - LUF + Schur compl. + Givens rotation */
alpar@9 46 FHV *fhv;
alpar@9 47 /* LP basis factorization (GLP_BF_FT) */
alpar@9 48 LPF *lpf;
alpar@9 49 /* LP basis factorization (GLP_BF_BG, GLP_BF_GR) */
alpar@9 50 int lu_size; /* luf.sv_size */
alpar@9 51 double piv_tol; /* luf.piv_tol */
alpar@9 52 int piv_lim; /* luf.piv_lim */
alpar@9 53 int suhl; /* luf.suhl */
alpar@9 54 double eps_tol; /* luf.eps_tol */
alpar@9 55 double max_gro; /* luf.max_gro */
alpar@9 56 int nfs_max; /* fhv.hh_max */
alpar@9 57 double upd_tol; /* fhv.upd_tol */
alpar@9 58 int nrs_max; /* lpf.n_max */
alpar@9 59 int rs_size; /* lpf.v_size */
alpar@9 60 /* internal control parameters */
alpar@9 61 int upd_lim;
alpar@9 62 /* the factorization update limit */
alpar@9 63 int upd_cnt;
alpar@9 64 /* the factorization update count */
alpar@9 65 };
alpar@9 66
alpar@9 67 /***********************************************************************
alpar@9 68 * NAME
alpar@9 69 *
alpar@9 70 * bfd_create_it - create LP basis factorization
alpar@9 71 *
alpar@9 72 * SYNOPSIS
alpar@9 73 *
alpar@9 74 * #include "glpbfd.h"
alpar@9 75 * BFD *bfd_create_it(void);
alpar@9 76 *
alpar@9 77 * DESCRIPTION
alpar@9 78 *
alpar@9 79 * The routine bfd_create_it creates a program object, which represents
alpar@9 80 * a factorization of LP basis.
alpar@9 81 *
alpar@9 82 * RETURNS
alpar@9 83 *
alpar@9 84 * The routine bfd_create_it returns a pointer to the object created. */
alpar@9 85
alpar@9 86 BFD *bfd_create_it(void)
alpar@9 87 { BFD *bfd;
alpar@9 88 bfd = xmalloc(sizeof(BFD));
alpar@9 89 bfd->valid = 0;
alpar@9 90 bfd->type = GLP_BF_FT;
alpar@9 91 bfd->fhv = NULL;
alpar@9 92 bfd->lpf = NULL;
alpar@9 93 bfd->lu_size = 0;
alpar@9 94 bfd->piv_tol = 0.10;
alpar@9 95 bfd->piv_lim = 4;
alpar@9 96 bfd->suhl = 1;
alpar@9 97 bfd->eps_tol = 1e-15;
alpar@9 98 bfd->max_gro = 1e+10;
alpar@9 99 bfd->nfs_max = 100;
alpar@9 100 bfd->upd_tol = 1e-6;
alpar@9 101 bfd->nrs_max = 100;
alpar@9 102 bfd->rs_size = 1000;
alpar@9 103 bfd->upd_lim = -1;
alpar@9 104 bfd->upd_cnt = 0;
alpar@9 105 return bfd;
alpar@9 106 }
alpar@9 107
alpar@9 108 /**********************************************************************/
alpar@9 109
alpar@9 110 void bfd_set_parm(BFD *bfd, const void *_parm)
alpar@9 111 { /* change LP basis factorization control parameters */
alpar@9 112 const glp_bfcp *parm = _parm;
alpar@9 113 xassert(bfd != NULL);
alpar@9 114 bfd->type = parm->type;
alpar@9 115 bfd->lu_size = parm->lu_size;
alpar@9 116 bfd->piv_tol = parm->piv_tol;
alpar@9 117 bfd->piv_lim = parm->piv_lim;
alpar@9 118 bfd->suhl = parm->suhl;
alpar@9 119 bfd->eps_tol = parm->eps_tol;
alpar@9 120 bfd->max_gro = parm->max_gro;
alpar@9 121 bfd->nfs_max = parm->nfs_max;
alpar@9 122 bfd->upd_tol = parm->upd_tol;
alpar@9 123 bfd->nrs_max = parm->nrs_max;
alpar@9 124 bfd->rs_size = parm->rs_size;
alpar@9 125 return;
alpar@9 126 }
alpar@9 127
alpar@9 128 /***********************************************************************
alpar@9 129 * NAME
alpar@9 130 *
alpar@9 131 * bfd_factorize - compute LP basis factorization
alpar@9 132 *
alpar@9 133 * SYNOPSIS
alpar@9 134 *
alpar@9 135 * #include "glpbfd.h"
alpar@9 136 * int bfd_factorize(BFD *bfd, int m, int bh[], int (*col)(void *info,
alpar@9 137 * int j, int ind[], double val[]), void *info);
alpar@9 138 *
alpar@9 139 * DESCRIPTION
alpar@9 140 *
alpar@9 141 * The routine bfd_factorize computes the factorization of the basis
alpar@9 142 * matrix B specified by the routine col.
alpar@9 143 *
alpar@9 144 * The parameter bfd specified the basis factorization data structure
alpar@9 145 * created with the routine bfd_create_it.
alpar@9 146 *
alpar@9 147 * The parameter m specifies the order of B, m > 0.
alpar@9 148 *
alpar@9 149 * The array bh specifies the basis header: bh[j], 1 <= j <= m, is the
alpar@9 150 * number of j-th column of B in some original matrix. The array bh is
alpar@9 151 * optional and can be specified as NULL.
alpar@9 152 *
alpar@9 153 * The formal routine col specifies the matrix B to be factorized. To
alpar@9 154 * obtain j-th column of A the routine bfd_factorize calls the routine
alpar@9 155 * col with the parameter j (1 <= j <= n). In response the routine col
alpar@9 156 * should store row indices and numerical values of non-zero elements
alpar@9 157 * of j-th column of B to locations ind[1,...,len] and val[1,...,len],
alpar@9 158 * respectively, where len is the number of non-zeros in j-th column
alpar@9 159 * returned on exit. Neither zero nor duplicate elements are allowed.
alpar@9 160 *
alpar@9 161 * The parameter info is a transit pointer passed to the routine col.
alpar@9 162 *
alpar@9 163 * RETURNS
alpar@9 164 *
alpar@9 165 * 0 The factorization has been successfully computed.
alpar@9 166 *
alpar@9 167 * BFD_ESING
alpar@9 168 * The specified matrix is singular within the working precision.
alpar@9 169 *
alpar@9 170 * BFD_ECOND
alpar@9 171 * The specified matrix is ill-conditioned.
alpar@9 172 *
alpar@9 173 * For more details see comments to the routine luf_factorize. */
alpar@9 174
alpar@9 175 int bfd_factorize(BFD *bfd, int m, const int bh[], int (*col)
alpar@9 176 (void *info, int j, int ind[], double val[]), void *info)
alpar@9 177 { LUF *luf;
alpar@9 178 int nov, ret;
alpar@9 179 xassert(bfd != NULL);
alpar@9 180 xassert(1 <= m && m <= M_MAX);
alpar@9 181 /* invalidate the factorization */
alpar@9 182 bfd->valid = 0;
alpar@9 183 /* create the factorization, if necessary */
alpar@9 184 nov = 0;
alpar@9 185 switch (bfd->type)
alpar@9 186 { case GLP_BF_FT:
alpar@9 187 if (bfd->lpf != NULL)
alpar@9 188 lpf_delete_it(bfd->lpf), bfd->lpf = NULL;
alpar@9 189 if (bfd->fhv == NULL)
alpar@9 190 bfd->fhv = fhv_create_it(), nov = 1;
alpar@9 191 break;
alpar@9 192 case GLP_BF_BG:
alpar@9 193 case GLP_BF_GR:
alpar@9 194 if (bfd->fhv != NULL)
alpar@9 195 fhv_delete_it(bfd->fhv), bfd->fhv = NULL;
alpar@9 196 if (bfd->lpf == NULL)
alpar@9 197 bfd->lpf = lpf_create_it(), nov = 1;
alpar@9 198 break;
alpar@9 199 default:
alpar@9 200 xassert(bfd != bfd);
alpar@9 201 }
alpar@9 202 /* set control parameters specific to LUF */
alpar@9 203 if (bfd->fhv != NULL)
alpar@9 204 luf = bfd->fhv->luf;
alpar@9 205 else if (bfd->lpf != NULL)
alpar@9 206 luf = bfd->lpf->luf;
alpar@9 207 else
alpar@9 208 xassert(bfd != bfd);
alpar@9 209 if (nov) luf->new_sva = bfd->lu_size;
alpar@9 210 luf->piv_tol = bfd->piv_tol;
alpar@9 211 luf->piv_lim = bfd->piv_lim;
alpar@9 212 luf->suhl = bfd->suhl;
alpar@9 213 luf->eps_tol = bfd->eps_tol;
alpar@9 214 luf->max_gro = bfd->max_gro;
alpar@9 215 /* set control parameters specific to FHV */
alpar@9 216 if (bfd->fhv != NULL)
alpar@9 217 { if (nov) bfd->fhv->hh_max = bfd->nfs_max;
alpar@9 218 bfd->fhv->upd_tol = bfd->upd_tol;
alpar@9 219 }
alpar@9 220 /* set control parameters specific to LPF */
alpar@9 221 if (bfd->lpf != NULL)
alpar@9 222 { if (nov) bfd->lpf->n_max = bfd->nrs_max;
alpar@9 223 if (nov) bfd->lpf->v_size = bfd->rs_size;
alpar@9 224 }
alpar@9 225 /* try to factorize the basis matrix */
alpar@9 226 if (bfd->fhv != NULL)
alpar@9 227 { switch (fhv_factorize(bfd->fhv, m, col, info))
alpar@9 228 { case 0:
alpar@9 229 break;
alpar@9 230 case FHV_ESING:
alpar@9 231 ret = BFD_ESING;
alpar@9 232 goto done;
alpar@9 233 case FHV_ECOND:
alpar@9 234 ret = BFD_ECOND;
alpar@9 235 goto done;
alpar@9 236 default:
alpar@9 237 xassert(bfd != bfd);
alpar@9 238 }
alpar@9 239 }
alpar@9 240 else if (bfd->lpf != NULL)
alpar@9 241 { switch (lpf_factorize(bfd->lpf, m, bh, col, info))
alpar@9 242 { case 0:
alpar@9 243 /* set the Schur complement update type */
alpar@9 244 switch (bfd->type)
alpar@9 245 { case GLP_BF_BG:
alpar@9 246 /* Bartels-Golub update */
alpar@9 247 bfd->lpf->scf->t_opt = SCF_TBG;
alpar@9 248 break;
alpar@9 249 case GLP_BF_GR:
alpar@9 250 /* Givens rotation update */
alpar@9 251 bfd->lpf->scf->t_opt = SCF_TGR;
alpar@9 252 break;
alpar@9 253 default:
alpar@9 254 xassert(bfd != bfd);
alpar@9 255 }
alpar@9 256 break;
alpar@9 257 case LPF_ESING:
alpar@9 258 ret = BFD_ESING;
alpar@9 259 goto done;
alpar@9 260 case LPF_ECOND:
alpar@9 261 ret = BFD_ECOND;
alpar@9 262 goto done;
alpar@9 263 default:
alpar@9 264 xassert(bfd != bfd);
alpar@9 265 }
alpar@9 266 }
alpar@9 267 else
alpar@9 268 xassert(bfd != bfd);
alpar@9 269 /* the basis matrix has been successfully factorized */
alpar@9 270 bfd->valid = 1;
alpar@9 271 bfd->upd_cnt = 0;
alpar@9 272 ret = 0;
alpar@9 273 done: /* return to the calling program */
alpar@9 274 return ret;
alpar@9 275 }
alpar@9 276
alpar@9 277 /***********************************************************************
alpar@9 278 * NAME
alpar@9 279 *
alpar@9 280 * bfd_ftran - perform forward transformation (solve system B*x = b)
alpar@9 281 *
alpar@9 282 * SYNOPSIS
alpar@9 283 *
alpar@9 284 * #include "glpbfd.h"
alpar@9 285 * void bfd_ftran(BFD *bfd, double x[]);
alpar@9 286 *
alpar@9 287 * DESCRIPTION
alpar@9 288 *
alpar@9 289 * The routine bfd_ftran performs forward transformation, i.e. solves
alpar@9 290 * the system B*x = b, where B is the basis matrix, x is the vector of
alpar@9 291 * unknowns to be computed, b is the vector of right-hand sides.
alpar@9 292 *
alpar@9 293 * On entry elements of the vector b should be stored in dense format
alpar@9 294 * in locations x[1], ..., x[m], where m is the number of rows. On exit
alpar@9 295 * the routine stores elements of the vector x in the same locations. */
alpar@9 296
alpar@9 297 void bfd_ftran(BFD *bfd, double x[])
alpar@9 298 { xassert(bfd != NULL);
alpar@9 299 xassert(bfd->valid);
alpar@9 300 if (bfd->fhv != NULL)
alpar@9 301 fhv_ftran(bfd->fhv, x);
alpar@9 302 else if (bfd->lpf != NULL)
alpar@9 303 lpf_ftran(bfd->lpf, x);
alpar@9 304 else
alpar@9 305 xassert(bfd != bfd);
alpar@9 306 return;
alpar@9 307 }
alpar@9 308
alpar@9 309 /***********************************************************************
alpar@9 310 * NAME
alpar@9 311 *
alpar@9 312 * bfd_btran - perform backward transformation (solve system B'*x = b)
alpar@9 313 *
alpar@9 314 * SYNOPSIS
alpar@9 315 *
alpar@9 316 * #include "glpbfd.h"
alpar@9 317 * void bfd_btran(BFD *bfd, double x[]);
alpar@9 318 *
alpar@9 319 * DESCRIPTION
alpar@9 320 *
alpar@9 321 * The routine bfd_btran performs backward transformation, i.e. solves
alpar@9 322 * the system B'*x = b, where B' is a matrix transposed to the basis
alpar@9 323 * matrix B, x is the vector of unknowns to be computed, b is the vector
alpar@9 324 * of right-hand sides.
alpar@9 325 *
alpar@9 326 * On entry elements of the vector b should be stored in dense format
alpar@9 327 * in locations x[1], ..., x[m], where m is the number of rows. On exit
alpar@9 328 * the routine stores elements of the vector x in the same locations. */
alpar@9 329
alpar@9 330 void bfd_btran(BFD *bfd, double x[])
alpar@9 331 { xassert(bfd != NULL);
alpar@9 332 xassert(bfd->valid);
alpar@9 333 if (bfd->fhv != NULL)
alpar@9 334 fhv_btran(bfd->fhv, x);
alpar@9 335 else if (bfd->lpf != NULL)
alpar@9 336 lpf_btran(bfd->lpf, x);
alpar@9 337 else
alpar@9 338 xassert(bfd != bfd);
alpar@9 339 return;
alpar@9 340 }
alpar@9 341
alpar@9 342 /***********************************************************************
alpar@9 343 * NAME
alpar@9 344 *
alpar@9 345 * bfd_update_it - update LP basis factorization
alpar@9 346 *
alpar@9 347 * SYNOPSIS
alpar@9 348 *
alpar@9 349 * #include "glpbfd.h"
alpar@9 350 * int bfd_update_it(BFD *bfd, int j, int bh, int len, const int ind[],
alpar@9 351 * const double val[]);
alpar@9 352 *
alpar@9 353 * DESCRIPTION
alpar@9 354 *
alpar@9 355 * The routine bfd_update_it updates the factorization of the basis
alpar@9 356 * matrix B after replacing its j-th column by a new vector.
alpar@9 357 *
alpar@9 358 * The parameter j specifies the number of column of B, which has been
alpar@9 359 * replaced, 1 <= j <= m, where m is the order of B.
alpar@9 360 *
alpar@9 361 * The parameter bh specifies the basis header entry for the new column
alpar@9 362 * of B, which is the number of the new column in some original matrix.
alpar@9 363 * This parameter is optional and can be specified as 0.
alpar@9 364 *
alpar@9 365 * Row indices and numerical values of non-zero elements of the new
alpar@9 366 * column of B should be placed in locations ind[1], ..., ind[len] and
alpar@9 367 * val[1], ..., val[len], resp., where len is the number of non-zeros
alpar@9 368 * in the column. Neither zero nor duplicate elements are allowed.
alpar@9 369 *
alpar@9 370 * RETURNS
alpar@9 371 *
alpar@9 372 * 0 The factorization has been successfully updated.
alpar@9 373 *
alpar@9 374 * BFD_ESING
alpar@9 375 * New basis matrix is singular within the working precision.
alpar@9 376 *
alpar@9 377 * BFD_ECHECK
alpar@9 378 * The factorization is inaccurate.
alpar@9 379 *
alpar@9 380 * BFD_ELIMIT
alpar@9 381 * Factorization update limit has been reached.
alpar@9 382 *
alpar@9 383 * BFD_EROOM
alpar@9 384 * Overflow of the sparse vector area.
alpar@9 385 *
alpar@9 386 * In case of non-zero return code the factorization becomes invalid.
alpar@9 387 * It should not be used until it has been recomputed with the routine
alpar@9 388 * bfd_factorize. */
alpar@9 389
alpar@9 390 int bfd_update_it(BFD *bfd, int j, int bh, int len, const int ind[],
alpar@9 391 const double val[])
alpar@9 392 { int ret;
alpar@9 393 xassert(bfd != NULL);
alpar@9 394 xassert(bfd->valid);
alpar@9 395 /* try to update the factorization */
alpar@9 396 if (bfd->fhv != NULL)
alpar@9 397 { switch (fhv_update_it(bfd->fhv, j, len, ind, val))
alpar@9 398 { case 0:
alpar@9 399 break;
alpar@9 400 case FHV_ESING:
alpar@9 401 bfd->valid = 0;
alpar@9 402 ret = BFD_ESING;
alpar@9 403 goto done;
alpar@9 404 case FHV_ECHECK:
alpar@9 405 bfd->valid = 0;
alpar@9 406 ret = BFD_ECHECK;
alpar@9 407 goto done;
alpar@9 408 case FHV_ELIMIT:
alpar@9 409 bfd->valid = 0;
alpar@9 410 ret = BFD_ELIMIT;
alpar@9 411 goto done;
alpar@9 412 case FHV_EROOM:
alpar@9 413 bfd->valid = 0;
alpar@9 414 ret = BFD_EROOM;
alpar@9 415 goto done;
alpar@9 416 default:
alpar@9 417 xassert(bfd != bfd);
alpar@9 418 }
alpar@9 419 }
alpar@9 420 else if (bfd->lpf != NULL)
alpar@9 421 { switch (lpf_update_it(bfd->lpf, j, bh, len, ind, val))
alpar@9 422 { case 0:
alpar@9 423 break;
alpar@9 424 case LPF_ESING:
alpar@9 425 bfd->valid = 0;
alpar@9 426 ret = BFD_ESING;
alpar@9 427 goto done;
alpar@9 428 case LPF_ELIMIT:
alpar@9 429 bfd->valid = 0;
alpar@9 430 ret = BFD_ELIMIT;
alpar@9 431 goto done;
alpar@9 432 default:
alpar@9 433 xassert(bfd != bfd);
alpar@9 434 }
alpar@9 435 }
alpar@9 436 else
alpar@9 437 xassert(bfd != bfd);
alpar@9 438 /* the factorization has been successfully updated */
alpar@9 439 /* increase the update count */
alpar@9 440 bfd->upd_cnt++;
alpar@9 441 ret = 0;
alpar@9 442 done: /* return to the calling program */
alpar@9 443 return ret;
alpar@9 444 }
alpar@9 445
alpar@9 446 /**********************************************************************/
alpar@9 447
alpar@9 448 int bfd_get_count(BFD *bfd)
alpar@9 449 { /* determine factorization update count */
alpar@9 450 xassert(bfd != NULL);
alpar@9 451 xassert(bfd->valid);
alpar@9 452 return bfd->upd_cnt;
alpar@9 453 }
alpar@9 454
alpar@9 455 /***********************************************************************
alpar@9 456 * NAME
alpar@9 457 *
alpar@9 458 * bfd_delete_it - delete LP basis factorization
alpar@9 459 *
alpar@9 460 * SYNOPSIS
alpar@9 461 *
alpar@9 462 * #include "glpbfd.h"
alpar@9 463 * void bfd_delete_it(BFD *bfd);
alpar@9 464 *
alpar@9 465 * DESCRIPTION
alpar@9 466 *
alpar@9 467 * The routine bfd_delete_it deletes LP basis factorization specified
alpar@9 468 * by the parameter fhv and frees all memory allocated to this program
alpar@9 469 * object. */
alpar@9 470
alpar@9 471 void bfd_delete_it(BFD *bfd)
alpar@9 472 { xassert(bfd != NULL);
alpar@9 473 if (bfd->fhv != NULL)
alpar@9 474 fhv_delete_it(bfd->fhv);
alpar@9 475 if (bfd->lpf != NULL)
alpar@9 476 lpf_delete_it(bfd->lpf);
alpar@9 477 xfree(bfd);
alpar@9 478 return;
alpar@9 479 }
alpar@9 480
alpar@9 481 /* eof */