alpar@1: /* glpapi12.c (basis factorization and simplex tableau routines) */ alpar@1: alpar@1: /*********************************************************************** alpar@1: * This code is part of GLPK (GNU Linear Programming Kit). alpar@1: * alpar@1: * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, alpar@1: * 2009, 2010 Andrew Makhorin, Department for Applied Informatics, alpar@1: * Moscow Aviation Institute, Moscow, Russia. All rights reserved. alpar@1: * E-mail: . alpar@1: * alpar@1: * GLPK is free software: you can redistribute it and/or modify it alpar@1: * under the terms of the GNU General Public License as published by alpar@1: * the Free Software Foundation, either version 3 of the License, or alpar@1: * (at your option) any later version. alpar@1: * alpar@1: * GLPK is distributed in the hope that it will be useful, but WITHOUT alpar@1: * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY alpar@1: * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public alpar@1: * License for more details. alpar@1: * alpar@1: * You should have received a copy of the GNU General Public License alpar@1: * along with GLPK. If not, see . alpar@1: ***********************************************************************/ alpar@1: alpar@1: #include "glpapi.h" alpar@1: alpar@1: /*********************************************************************** alpar@1: * NAME alpar@1: * alpar@1: * glp_bf_exists - check if the basis factorization exists alpar@1: * alpar@1: * SYNOPSIS alpar@1: * alpar@1: * int glp_bf_exists(glp_prob *lp); alpar@1: * alpar@1: * RETURNS alpar@1: * alpar@1: * If the basis factorization for the current basis associated with alpar@1: * the specified problem object exists and therefore is available for alpar@1: * computations, the routine glp_bf_exists returns non-zero. Otherwise alpar@1: * the routine returns zero. */ alpar@1: alpar@1: int glp_bf_exists(glp_prob *lp) alpar@1: { int ret; alpar@1: ret = (lp->m == 0 || lp->valid); alpar@1: return ret; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * NAME alpar@1: * alpar@1: * glp_factorize - compute the basis factorization alpar@1: * alpar@1: * SYNOPSIS alpar@1: * alpar@1: * int glp_factorize(glp_prob *lp); alpar@1: * alpar@1: * DESCRIPTION alpar@1: * alpar@1: * The routine glp_factorize computes the basis factorization for the alpar@1: * current basis associated with the specified problem object. alpar@1: * alpar@1: * RETURNS alpar@1: * alpar@1: * 0 The basis factorization has been successfully computed. alpar@1: * alpar@1: * GLP_EBADB alpar@1: * The basis matrix is invalid, i.e. the number of basic (auxiliary alpar@1: * and structural) variables differs from the number of rows in the alpar@1: * problem object. alpar@1: * alpar@1: * GLP_ESING alpar@1: * The basis matrix is singular within the working precision. alpar@1: * alpar@1: * GLP_ECOND alpar@1: * The basis matrix is ill-conditioned. */ alpar@1: alpar@1: static int b_col(void *info, int j, int ind[], double val[]) alpar@1: { glp_prob *lp = info; alpar@1: int m = lp->m; alpar@1: GLPAIJ *aij; alpar@1: int k, len; alpar@1: xassert(1 <= j && j <= m); alpar@1: /* determine the ordinal number of basic auxiliary or structural alpar@1: variable x[k] corresponding to basic variable xB[j] */ alpar@1: k = lp->head[j]; alpar@1: /* build j-th column of the basic matrix, which is k-th column of alpar@1: the scaled augmented matrix (I | -R*A*S) */ alpar@1: if (k <= m) alpar@1: { /* x[k] is auxiliary variable */ alpar@1: len = 1; alpar@1: ind[1] = k; alpar@1: val[1] = 1.0; alpar@1: } alpar@1: else alpar@1: { /* x[k] is structural variable */ alpar@1: len = 0; alpar@1: for (aij = lp->col[k-m]->ptr; aij != NULL; aij = aij->c_next) alpar@1: { len++; alpar@1: ind[len] = aij->row->i; alpar@1: val[len] = - aij->row->rii * aij->val * aij->col->sjj; alpar@1: } alpar@1: } alpar@1: return len; alpar@1: } alpar@1: alpar@1: static void copy_bfcp(glp_prob *lp); alpar@1: alpar@1: int glp_factorize(glp_prob *lp) alpar@1: { int m = lp->m; alpar@1: int n = lp->n; alpar@1: GLPROW **row = lp->row; alpar@1: GLPCOL **col = lp->col; alpar@1: int *head = lp->head; alpar@1: int j, k, stat, ret; alpar@1: /* invalidate the basis factorization */ alpar@1: lp->valid = 0; alpar@1: /* build the basis header */ alpar@1: j = 0; alpar@1: for (k = 1; k <= m+n; k++) alpar@1: { if (k <= m) alpar@1: { stat = row[k]->stat; alpar@1: row[k]->bind = 0; alpar@1: } alpar@1: else alpar@1: { stat = col[k-m]->stat; alpar@1: col[k-m]->bind = 0; alpar@1: } alpar@1: if (stat == GLP_BS) alpar@1: { j++; alpar@1: if (j > m) alpar@1: { /* too many basic variables */ alpar@1: ret = GLP_EBADB; alpar@1: goto fini; alpar@1: } alpar@1: head[j] = k; alpar@1: if (k <= m) alpar@1: row[k]->bind = j; alpar@1: else alpar@1: col[k-m]->bind = j; alpar@1: } alpar@1: } alpar@1: if (j < m) alpar@1: { /* too few basic variables */ alpar@1: ret = GLP_EBADB; alpar@1: goto fini; alpar@1: } alpar@1: /* try to factorize the basis matrix */ alpar@1: if (m > 0) alpar@1: { if (lp->bfd == NULL) alpar@1: { lp->bfd = bfd_create_it(); alpar@1: copy_bfcp(lp); alpar@1: } alpar@1: switch (bfd_factorize(lp->bfd, m, lp->head, b_col, lp)) alpar@1: { case 0: alpar@1: /* ok */ alpar@1: break; alpar@1: case BFD_ESING: alpar@1: /* singular matrix */ alpar@1: ret = GLP_ESING; alpar@1: goto fini; alpar@1: case BFD_ECOND: alpar@1: /* ill-conditioned matrix */ alpar@1: ret = GLP_ECOND; alpar@1: goto fini; alpar@1: default: alpar@1: xassert(lp != lp); alpar@1: } alpar@1: lp->valid = 1; alpar@1: } alpar@1: /* factorization successful */ alpar@1: ret = 0; alpar@1: fini: /* bring the return code to the calling program */ alpar@1: return ret; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * NAME alpar@1: * alpar@1: * glp_bf_updated - check if the basis factorization has been updated alpar@1: * alpar@1: * SYNOPSIS alpar@1: * alpar@1: * int glp_bf_updated(glp_prob *lp); alpar@1: * alpar@1: * RETURNS alpar@1: * alpar@1: * If the basis factorization has been just computed from scratch, the alpar@1: * routine glp_bf_updated returns zero. Otherwise, if the factorization alpar@1: * has been updated one or more times, the routine returns non-zero. */ alpar@1: alpar@1: int glp_bf_updated(glp_prob *lp) alpar@1: { int cnt; alpar@1: if (!(lp->m == 0 || lp->valid)) alpar@1: xerror("glp_bf_update: basis factorization does not exist\n"); alpar@1: #if 0 /* 15/XI-2009 */ alpar@1: cnt = (lp->m == 0 ? 0 : lp->bfd->upd_cnt); alpar@1: #else alpar@1: cnt = (lp->m == 0 ? 0 : bfd_get_count(lp->bfd)); alpar@1: #endif alpar@1: return cnt; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * NAME alpar@1: * alpar@1: * glp_get_bfcp - retrieve basis factorization control parameters alpar@1: * alpar@1: * SYNOPSIS alpar@1: * alpar@1: * void glp_get_bfcp(glp_prob *lp, glp_bfcp *parm); alpar@1: * alpar@1: * DESCRIPTION alpar@1: * alpar@1: * The routine glp_get_bfcp retrieves control parameters, which are alpar@1: * used on computing and updating the basis factorization associated alpar@1: * with the specified problem object. alpar@1: * alpar@1: * Current values of control parameters are stored by the routine in alpar@1: * a glp_bfcp structure, which the parameter parm points to. */ alpar@1: alpar@1: void glp_get_bfcp(glp_prob *lp, glp_bfcp *parm) alpar@1: { glp_bfcp *bfcp = lp->bfcp; alpar@1: if (bfcp == NULL) alpar@1: { parm->type = GLP_BF_FT; alpar@1: parm->lu_size = 0; alpar@1: parm->piv_tol = 0.10; alpar@1: parm->piv_lim = 4; alpar@1: parm->suhl = GLP_ON; alpar@1: parm->eps_tol = 1e-15; alpar@1: parm->max_gro = 1e+10; alpar@1: parm->nfs_max = 100; alpar@1: parm->upd_tol = 1e-6; alpar@1: parm->nrs_max = 100; alpar@1: parm->rs_size = 0; alpar@1: } alpar@1: else alpar@1: memcpy(parm, bfcp, sizeof(glp_bfcp)); alpar@1: return; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * NAME alpar@1: * alpar@1: * glp_set_bfcp - change basis factorization control parameters alpar@1: * alpar@1: * SYNOPSIS alpar@1: * alpar@1: * void glp_set_bfcp(glp_prob *lp, const glp_bfcp *parm); alpar@1: * alpar@1: * DESCRIPTION alpar@1: * alpar@1: * The routine glp_set_bfcp changes control parameters, which are used alpar@1: * by internal GLPK routines in computing and updating the basis alpar@1: * factorization associated with the specified problem object. alpar@1: * alpar@1: * New values of the control parameters should be passed in a structure alpar@1: * glp_bfcp, which the parameter parm points to. alpar@1: * alpar@1: * The parameter parm can be specified as NULL, in which case all alpar@1: * control parameters are reset to their default values. */ alpar@1: alpar@1: #if 0 /* 15/XI-2009 */ alpar@1: static void copy_bfcp(glp_prob *lp) alpar@1: { glp_bfcp _parm, *parm = &_parm; alpar@1: BFD *bfd = lp->bfd; alpar@1: glp_get_bfcp(lp, parm); alpar@1: xassert(bfd != NULL); alpar@1: bfd->type = parm->type; alpar@1: bfd->lu_size = parm->lu_size; alpar@1: bfd->piv_tol = parm->piv_tol; alpar@1: bfd->piv_lim = parm->piv_lim; alpar@1: bfd->suhl = parm->suhl; alpar@1: bfd->eps_tol = parm->eps_tol; alpar@1: bfd->max_gro = parm->max_gro; alpar@1: bfd->nfs_max = parm->nfs_max; alpar@1: bfd->upd_tol = parm->upd_tol; alpar@1: bfd->nrs_max = parm->nrs_max; alpar@1: bfd->rs_size = parm->rs_size; alpar@1: return; alpar@1: } alpar@1: #else alpar@1: static void copy_bfcp(glp_prob *lp) alpar@1: { glp_bfcp _parm, *parm = &_parm; alpar@1: glp_get_bfcp(lp, parm); alpar@1: bfd_set_parm(lp->bfd, parm); alpar@1: return; alpar@1: } alpar@1: #endif alpar@1: alpar@1: void glp_set_bfcp(glp_prob *lp, const glp_bfcp *parm) alpar@1: { glp_bfcp *bfcp = lp->bfcp; alpar@1: if (parm == NULL) alpar@1: { /* reset to default values */ alpar@1: if (bfcp != NULL) alpar@1: xfree(bfcp), lp->bfcp = NULL; alpar@1: } alpar@1: else alpar@1: { /* set to specified values */ alpar@1: if (bfcp == NULL) alpar@1: bfcp = lp->bfcp = xmalloc(sizeof(glp_bfcp)); alpar@1: memcpy(bfcp, parm, sizeof(glp_bfcp)); alpar@1: if (!(bfcp->type == GLP_BF_FT || bfcp->type == GLP_BF_BG || alpar@1: bfcp->type == GLP_BF_GR)) alpar@1: xerror("glp_set_bfcp: type = %d; invalid parameter\n", alpar@1: bfcp->type); alpar@1: if (bfcp->lu_size < 0) alpar@1: xerror("glp_set_bfcp: lu_size = %d; invalid parameter\n", alpar@1: bfcp->lu_size); alpar@1: if (!(0.0 < bfcp->piv_tol && bfcp->piv_tol < 1.0)) alpar@1: xerror("glp_set_bfcp: piv_tol = %g; invalid parameter\n", alpar@1: bfcp->piv_tol); alpar@1: if (bfcp->piv_lim < 1) alpar@1: xerror("glp_set_bfcp: piv_lim = %d; invalid parameter\n", alpar@1: bfcp->piv_lim); alpar@1: if (!(bfcp->suhl == GLP_ON || bfcp->suhl == GLP_OFF)) alpar@1: xerror("glp_set_bfcp: suhl = %d; invalid parameter\n", alpar@1: bfcp->suhl); alpar@1: if (!(0.0 <= bfcp->eps_tol && bfcp->eps_tol <= 1e-6)) alpar@1: xerror("glp_set_bfcp: eps_tol = %g; invalid parameter\n", alpar@1: bfcp->eps_tol); alpar@1: if (bfcp->max_gro < 1.0) alpar@1: xerror("glp_set_bfcp: max_gro = %g; invalid parameter\n", alpar@1: bfcp->max_gro); alpar@1: if (!(1 <= bfcp->nfs_max && bfcp->nfs_max <= 32767)) alpar@1: xerror("glp_set_bfcp: nfs_max = %d; invalid parameter\n", alpar@1: bfcp->nfs_max); alpar@1: if (!(0.0 < bfcp->upd_tol && bfcp->upd_tol < 1.0)) alpar@1: xerror("glp_set_bfcp: upd_tol = %g; invalid parameter\n", alpar@1: bfcp->upd_tol); alpar@1: if (!(1 <= bfcp->nrs_max && bfcp->nrs_max <= 32767)) alpar@1: xerror("glp_set_bfcp: nrs_max = %d; invalid parameter\n", alpar@1: bfcp->nrs_max); alpar@1: if (bfcp->rs_size < 0) alpar@1: xerror("glp_set_bfcp: rs_size = %d; invalid parameter\n", alpar@1: bfcp->nrs_max); alpar@1: if (bfcp->rs_size == 0) alpar@1: bfcp->rs_size = 20 * bfcp->nrs_max; alpar@1: } alpar@1: if (lp->bfd != NULL) copy_bfcp(lp); alpar@1: return; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * NAME alpar@1: * alpar@1: * glp_get_bhead - retrieve the basis header information alpar@1: * alpar@1: * SYNOPSIS alpar@1: * alpar@1: * int glp_get_bhead(glp_prob *lp, int k); alpar@1: * alpar@1: * DESCRIPTION alpar@1: * alpar@1: * The routine glp_get_bhead returns the basis header information for alpar@1: * the current basis associated with the specified problem object. alpar@1: * alpar@1: * RETURNS alpar@1: * alpar@1: * If xB[k], 1 <= k <= m, is i-th auxiliary variable (1 <= i <= m), the alpar@1: * routine returns i. Otherwise, if xB[k] is j-th structural variable alpar@1: * (1 <= j <= n), the routine returns m+j. Here m is the number of rows alpar@1: * and n is the number of columns in the problem object. */ alpar@1: alpar@1: int glp_get_bhead(glp_prob *lp, int k) alpar@1: { if (!(lp->m == 0 || lp->valid)) alpar@1: xerror("glp_get_bhead: basis factorization does not exist\n"); alpar@1: if (!(1 <= k && k <= lp->m)) alpar@1: xerror("glp_get_bhead: k = %d; index out of range\n", k); alpar@1: return lp->head[k]; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * NAME alpar@1: * alpar@1: * glp_get_row_bind - retrieve row index in the basis header alpar@1: * alpar@1: * SYNOPSIS alpar@1: * alpar@1: * int glp_get_row_bind(glp_prob *lp, int i); alpar@1: * alpar@1: * RETURNS alpar@1: * alpar@1: * The routine glp_get_row_bind returns the index k of basic variable alpar@1: * xB[k], 1 <= k <= m, which is i-th auxiliary variable, 1 <= i <= m, alpar@1: * in the current basis associated with the specified problem object, alpar@1: * where m is the number of rows. However, if i-th auxiliary variable alpar@1: * is non-basic, the routine returns zero. */ alpar@1: alpar@1: int glp_get_row_bind(glp_prob *lp, int i) alpar@1: { if (!(lp->m == 0 || lp->valid)) alpar@1: xerror("glp_get_row_bind: basis factorization does not exist\n" alpar@1: ); alpar@1: if (!(1 <= i && i <= lp->m)) alpar@1: xerror("glp_get_row_bind: i = %d; row number out of range\n", alpar@1: i); alpar@1: return lp->row[i]->bind; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * NAME alpar@1: * alpar@1: * glp_get_col_bind - retrieve column index in the basis header alpar@1: * alpar@1: * SYNOPSIS alpar@1: * alpar@1: * int glp_get_col_bind(glp_prob *lp, int j); alpar@1: * alpar@1: * RETURNS alpar@1: * alpar@1: * The routine glp_get_col_bind returns the index k of basic variable alpar@1: * xB[k], 1 <= k <= m, which is j-th structural variable, 1 <= j <= n, alpar@1: * in the current basis associated with the specified problem object, alpar@1: * where m is the number of rows, n is the number of columns. However, alpar@1: * if j-th structural variable is non-basic, the routine returns zero.*/ alpar@1: alpar@1: int glp_get_col_bind(glp_prob *lp, int j) alpar@1: { if (!(lp->m == 0 || lp->valid)) alpar@1: xerror("glp_get_col_bind: basis factorization does not exist\n" alpar@1: ); alpar@1: if (!(1 <= j && j <= lp->n)) alpar@1: xerror("glp_get_col_bind: j = %d; column number out of range\n" alpar@1: , j); alpar@1: return lp->col[j]->bind; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * NAME alpar@1: * alpar@1: * glp_ftran - perform forward transformation (solve system B*x = b) alpar@1: * alpar@1: * SYNOPSIS alpar@1: * alpar@1: * void glp_ftran(glp_prob *lp, double x[]); alpar@1: * alpar@1: * DESCRIPTION alpar@1: * alpar@1: * The routine glp_ftran performs forward transformation, i.e. solves alpar@1: * the system B*x = b, where B is the basis matrix corresponding to the alpar@1: * current basis for the specified problem object, x is the vector of alpar@1: * unknowns to be computed, b is the vector of right-hand sides. alpar@1: * alpar@1: * On entry elements of the vector b should be stored in dense format alpar@1: * in locations x[1], ..., x[m], where m is the number of rows. On exit alpar@1: * the routine stores elements of the vector x in the same locations. alpar@1: * alpar@1: * SCALING/UNSCALING alpar@1: * alpar@1: * Let A~ = (I | -A) is the augmented constraint matrix of the original alpar@1: * (unscaled) problem. In the scaled LP problem instead the matrix A the alpar@1: * scaled matrix A" = R*A*S is actually used, so alpar@1: * alpar@1: * A~" = (I | A") = (I | R*A*S) = (R*I*inv(R) | R*A*S) = alpar@1: * (1) alpar@1: * = R*(I | A)*S~ = R*A~*S~, alpar@1: * alpar@1: * is the scaled augmented constraint matrix, where R and S are diagonal alpar@1: * scaling matrices used to scale rows and columns of the matrix A, and alpar@1: * alpar@1: * S~ = diag(inv(R) | S) (2) alpar@1: * alpar@1: * is an augmented diagonal scaling matrix. alpar@1: * alpar@1: * By definition: alpar@1: * alpar@1: * A~ = (B | N), (3) alpar@1: * alpar@1: * where B is the basic matrix, which consists of basic columns of the alpar@1: * augmented constraint matrix A~, and N is a matrix, which consists of alpar@1: * non-basic columns of A~. From (1) it follows that: alpar@1: * alpar@1: * A~" = (B" | N") = (R*B*SB | R*N*SN), (4) alpar@1: * alpar@1: * where SB and SN are parts of the augmented scaling matrix S~, which alpar@1: * correspond to basic and non-basic variables, respectively. Therefore alpar@1: * alpar@1: * B" = R*B*SB, (5) alpar@1: * alpar@1: * which is the scaled basis matrix. */ alpar@1: alpar@1: void glp_ftran(glp_prob *lp, double x[]) alpar@1: { int m = lp->m; alpar@1: GLPROW **row = lp->row; alpar@1: GLPCOL **col = lp->col; alpar@1: int i, k; alpar@1: /* B*x = b ===> (R*B*SB)*(inv(SB)*x) = R*b ===> alpar@1: B"*x" = b", where b" = R*b, x = SB*x" */ alpar@1: if (!(m == 0 || lp->valid)) alpar@1: xerror("glp_ftran: basis factorization does not exist\n"); alpar@1: /* b" := R*b */ alpar@1: for (i = 1; i <= m; i++) alpar@1: x[i] *= row[i]->rii; alpar@1: /* x" := inv(B")*b" */ alpar@1: if (m > 0) bfd_ftran(lp->bfd, x); alpar@1: /* x := SB*x" */ alpar@1: for (i = 1; i <= m; i++) alpar@1: { k = lp->head[i]; alpar@1: if (k <= m) alpar@1: x[i] /= row[k]->rii; alpar@1: else alpar@1: x[i] *= col[k-m]->sjj; alpar@1: } alpar@1: return; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * NAME alpar@1: * alpar@1: * glp_btran - perform backward transformation (solve system B'*x = b) alpar@1: * alpar@1: * SYNOPSIS alpar@1: * alpar@1: * void glp_btran(glp_prob *lp, double x[]); alpar@1: * alpar@1: * DESCRIPTION alpar@1: * alpar@1: * The routine glp_btran performs backward transformation, i.e. solves alpar@1: * the system B'*x = b, where B' is a matrix transposed to the basis alpar@1: * matrix corresponding to the current basis for the specified problem alpar@1: * problem object, x is the vector of unknowns to be computed, b is the alpar@1: * vector of right-hand sides. alpar@1: * alpar@1: * On entry elements of the vector b should be stored in dense format alpar@1: * in locations x[1], ..., x[m], where m is the number of rows. On exit alpar@1: * the routine stores elements of the vector x in the same locations. alpar@1: * alpar@1: * SCALING/UNSCALING alpar@1: * alpar@1: * See comments to the routine glp_ftran. */ alpar@1: alpar@1: void glp_btran(glp_prob *lp, double x[]) alpar@1: { int m = lp->m; alpar@1: GLPROW **row = lp->row; alpar@1: GLPCOL **col = lp->col; alpar@1: int i, k; alpar@1: /* B'*x = b ===> (SB*B'*R)*(inv(R)*x) = SB*b ===> alpar@1: (B")'*x" = b", where b" = SB*b, x = R*x" */ alpar@1: if (!(m == 0 || lp->valid)) alpar@1: xerror("glp_btran: basis factorization does not exist\n"); alpar@1: /* b" := SB*b */ alpar@1: for (i = 1; i <= m; i++) alpar@1: { k = lp->head[i]; alpar@1: if (k <= m) alpar@1: x[i] /= row[k]->rii; alpar@1: else alpar@1: x[i] *= col[k-m]->sjj; alpar@1: } alpar@1: /* x" := inv[(B")']*b" */ alpar@1: if (m > 0) bfd_btran(lp->bfd, x); alpar@1: /* x := R*x" */ alpar@1: for (i = 1; i <= m; i++) alpar@1: x[i] *= row[i]->rii; alpar@1: return; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * NAME alpar@1: * alpar@1: * glp_warm_up - "warm up" LP basis alpar@1: * alpar@1: * SYNOPSIS alpar@1: * alpar@1: * int glp_warm_up(glp_prob *P); alpar@1: * alpar@1: * DESCRIPTION alpar@1: * alpar@1: * The routine glp_warm_up "warms up" the LP basis for the specified alpar@1: * problem object using current statuses assigned to rows and columns alpar@1: * (that is, to auxiliary and structural variables). alpar@1: * alpar@1: * This operation includes computing factorization of the basis matrix alpar@1: * (if it does not exist), computing primal and dual components of basic alpar@1: * solution, and determining the solution status. alpar@1: * alpar@1: * RETURNS alpar@1: * alpar@1: * 0 The operation has been successfully performed. alpar@1: * alpar@1: * GLP_EBADB alpar@1: * The basis matrix is invalid, i.e. the number of basic (auxiliary alpar@1: * and structural) variables differs from the number of rows in the alpar@1: * problem object. alpar@1: * alpar@1: * GLP_ESING alpar@1: * The basis matrix is singular within the working precision. alpar@1: * alpar@1: * GLP_ECOND alpar@1: * The basis matrix is ill-conditioned. */ alpar@1: alpar@1: int glp_warm_up(glp_prob *P) alpar@1: { GLPROW *row; alpar@1: GLPCOL *col; alpar@1: GLPAIJ *aij; alpar@1: int i, j, type, ret; alpar@1: double eps, temp, *work; alpar@1: /* invalidate basic solution */ alpar@1: P->pbs_stat = P->dbs_stat = GLP_UNDEF; alpar@1: P->obj_val = 0.0; alpar@1: P->some = 0; alpar@1: for (i = 1; i <= P->m; i++) alpar@1: { row = P->row[i]; alpar@1: row->prim = row->dual = 0.0; alpar@1: } alpar@1: for (j = 1; j <= P->n; j++) alpar@1: { col = P->col[j]; alpar@1: col->prim = col->dual = 0.0; alpar@1: } alpar@1: /* compute the basis factorization, if necessary */ alpar@1: if (!glp_bf_exists(P)) alpar@1: { ret = glp_factorize(P); alpar@1: if (ret != 0) goto done; alpar@1: } alpar@1: /* allocate working array */ alpar@1: work = xcalloc(1+P->m, sizeof(double)); alpar@1: /* determine and store values of non-basic variables, compute alpar@1: vector (- N * xN) */ alpar@1: for (i = 1; i <= P->m; i++) alpar@1: work[i] = 0.0; alpar@1: for (i = 1; i <= P->m; i++) alpar@1: { row = P->row[i]; alpar@1: if (row->stat == GLP_BS) alpar@1: continue; alpar@1: else if (row->stat == GLP_NL) alpar@1: row->prim = row->lb; alpar@1: else if (row->stat == GLP_NU) alpar@1: row->prim = row->ub; alpar@1: else if (row->stat == GLP_NF) alpar@1: row->prim = 0.0; alpar@1: else if (row->stat == GLP_NS) alpar@1: row->prim = row->lb; alpar@1: else alpar@1: xassert(row != row); alpar@1: /* N[j] is i-th column of matrix (I|-A) */ alpar@1: work[i] -= row->prim; alpar@1: } alpar@1: for (j = 1; j <= P->n; j++) alpar@1: { col = P->col[j]; alpar@1: if (col->stat == GLP_BS) alpar@1: continue; alpar@1: else if (col->stat == GLP_NL) alpar@1: col->prim = col->lb; alpar@1: else if (col->stat == GLP_NU) alpar@1: col->prim = col->ub; alpar@1: else if (col->stat == GLP_NF) alpar@1: col->prim = 0.0; alpar@1: else if (col->stat == GLP_NS) alpar@1: col->prim = col->lb; alpar@1: else alpar@1: xassert(col != col); alpar@1: /* N[j] is (m+j)-th column of matrix (I|-A) */ alpar@1: if (col->prim != 0.0) alpar@1: { for (aij = col->ptr; aij != NULL; aij = aij->c_next) alpar@1: work[aij->row->i] += aij->val * col->prim; alpar@1: } alpar@1: } alpar@1: /* compute vector of basic variables xB = - inv(B) * N * xN */ alpar@1: glp_ftran(P, work); alpar@1: /* store values of basic variables, check primal feasibility */ alpar@1: P->pbs_stat = GLP_FEAS; alpar@1: for (i = 1; i <= P->m; i++) alpar@1: { row = P->row[i]; alpar@1: if (row->stat != GLP_BS) alpar@1: continue; alpar@1: row->prim = work[row->bind]; alpar@1: type = row->type; alpar@1: if (type == GLP_LO || type == GLP_DB || type == GLP_FX) alpar@1: { eps = 1e-6 + 1e-9 * fabs(row->lb); alpar@1: if (row->prim < row->lb - eps) alpar@1: P->pbs_stat = GLP_INFEAS; alpar@1: } alpar@1: if (type == GLP_UP || type == GLP_DB || type == GLP_FX) alpar@1: { eps = 1e-6 + 1e-9 * fabs(row->ub); alpar@1: if (row->prim > row->ub + eps) alpar@1: P->pbs_stat = GLP_INFEAS; alpar@1: } alpar@1: } alpar@1: for (j = 1; j <= P->n; j++) alpar@1: { col = P->col[j]; alpar@1: if (col->stat != GLP_BS) alpar@1: continue; alpar@1: col->prim = work[col->bind]; alpar@1: type = col->type; alpar@1: if (type == GLP_LO || type == GLP_DB || type == GLP_FX) alpar@1: { eps = 1e-6 + 1e-9 * fabs(col->lb); alpar@1: if (col->prim < col->lb - eps) alpar@1: P->pbs_stat = GLP_INFEAS; alpar@1: } alpar@1: if (type == GLP_UP || type == GLP_DB || type == GLP_FX) alpar@1: { eps = 1e-6 + 1e-9 * fabs(col->ub); alpar@1: if (col->prim > col->ub + eps) alpar@1: P->pbs_stat = GLP_INFEAS; alpar@1: } alpar@1: } alpar@1: /* compute value of the objective function */ alpar@1: P->obj_val = P->c0; alpar@1: for (j = 1; j <= P->n; j++) alpar@1: { col = P->col[j]; alpar@1: P->obj_val += col->coef * col->prim; alpar@1: } alpar@1: /* build vector cB of objective coefficients at basic variables */ alpar@1: for (i = 1; i <= P->m; i++) alpar@1: work[i] = 0.0; alpar@1: for (j = 1; j <= P->n; j++) alpar@1: { col = P->col[j]; alpar@1: if (col->stat == GLP_BS) alpar@1: work[col->bind] = col->coef; alpar@1: } alpar@1: /* compute vector of simplex multipliers pi = inv(B') * cB */ alpar@1: glp_btran(P, work); alpar@1: /* compute and store reduced costs of non-basic variables d[j] = alpar@1: c[j] - N'[j] * pi, check dual feasibility */ alpar@1: P->dbs_stat = GLP_FEAS; alpar@1: for (i = 1; i <= P->m; i++) alpar@1: { row = P->row[i]; alpar@1: if (row->stat == GLP_BS) alpar@1: { row->dual = 0.0; alpar@1: continue; alpar@1: } alpar@1: /* N[j] is i-th column of matrix (I|-A) */ alpar@1: row->dual = - work[i]; alpar@1: type = row->type; alpar@1: temp = (P->dir == GLP_MIN ? + row->dual : - row->dual); alpar@1: if ((type == GLP_FR || type == GLP_LO) && temp < -1e-5 || alpar@1: (type == GLP_FR || type == GLP_UP) && temp > +1e-5) alpar@1: P->dbs_stat = GLP_INFEAS; alpar@1: } alpar@1: for (j = 1; j <= P->n; j++) alpar@1: { col = P->col[j]; alpar@1: if (col->stat == GLP_BS) alpar@1: { col->dual = 0.0; alpar@1: continue; alpar@1: } alpar@1: /* N[j] is (m+j)-th column of matrix (I|-A) */ alpar@1: col->dual = col->coef; alpar@1: for (aij = col->ptr; aij != NULL; aij = aij->c_next) alpar@1: col->dual += aij->val * work[aij->row->i]; alpar@1: type = col->type; alpar@1: temp = (P->dir == GLP_MIN ? + col->dual : - col->dual); alpar@1: if ((type == GLP_FR || type == GLP_LO) && temp < -1e-5 || alpar@1: (type == GLP_FR || type == GLP_UP) && temp > +1e-5) alpar@1: P->dbs_stat = GLP_INFEAS; alpar@1: } alpar@1: /* free working array */ alpar@1: xfree(work); alpar@1: ret = 0; alpar@1: done: return ret; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * NAME alpar@1: * alpar@1: * glp_eval_tab_row - compute row of the simplex tableau alpar@1: * alpar@1: * SYNOPSIS alpar@1: * alpar@1: * int glp_eval_tab_row(glp_prob *lp, int k, int ind[], double val[]); alpar@1: * alpar@1: * DESCRIPTION alpar@1: * alpar@1: * The routine glp_eval_tab_row computes a row of the current simplex alpar@1: * tableau for the basic variable, which is specified by the number k: alpar@1: * if 1 <= k <= m, x[k] is k-th auxiliary variable; if m+1 <= k <= m+n, alpar@1: * x[k] is (k-m)-th structural variable, where m is number of rows, and alpar@1: * n is number of columns. The current basis must be available. alpar@1: * alpar@1: * The routine stores column indices and numerical values of non-zero alpar@1: * elements of the computed row using sparse format to the locations alpar@1: * ind[1], ..., ind[len] and val[1], ..., val[len], respectively, where alpar@1: * 0 <= len <= n is number of non-zeros returned on exit. alpar@1: * alpar@1: * Element indices stored in the array ind have the same sense as the alpar@1: * index k, i.e. indices 1 to m denote auxiliary variables and indices alpar@1: * m+1 to m+n denote structural ones (all these variables are obviously alpar@1: * non-basic by definition). alpar@1: * alpar@1: * The computed row shows how the specified basic variable x[k] = xB[i] alpar@1: * depends on non-basic variables: alpar@1: * alpar@1: * xB[i] = alfa[i,1]*xN[1] + alfa[i,2]*xN[2] + ... + alfa[i,n]*xN[n], alpar@1: * alpar@1: * where alfa[i,j] are elements of the simplex table row, xN[j] are alpar@1: * non-basic (auxiliary and structural) variables. alpar@1: * alpar@1: * RETURNS alpar@1: * alpar@1: * The routine returns number of non-zero elements in the simplex table alpar@1: * row stored in the arrays ind and val. alpar@1: * alpar@1: * BACKGROUND alpar@1: * alpar@1: * The system of equality constraints of the LP problem is: alpar@1: * alpar@1: * xR = A * xS, (1) alpar@1: * alpar@1: * where xR is the vector of auxliary variables, xS is the vector of alpar@1: * structural variables, A is the matrix of constraint coefficients. alpar@1: * alpar@1: * The system (1) can be written in homogenous form as follows: alpar@1: * alpar@1: * A~ * x = 0, (2) alpar@1: * alpar@1: * where A~ = (I | -A) is the augmented constraint matrix (has m rows alpar@1: * and m+n columns), x = (xR | xS) is the vector of all (auxiliary and alpar@1: * structural) variables. alpar@1: * alpar@1: * By definition for the current basis we have: alpar@1: * alpar@1: * A~ = (B | N), (3) alpar@1: * alpar@1: * where B is the basis matrix. Thus, the system (2) can be written as: alpar@1: * alpar@1: * B * xB + N * xN = 0. (4) alpar@1: * alpar@1: * From (4) it follows that: alpar@1: * alpar@1: * xB = A^ * xN, (5) alpar@1: * alpar@1: * where the matrix alpar@1: * alpar@1: * A^ = - inv(B) * N (6) alpar@1: * alpar@1: * is called the simplex table. alpar@1: * alpar@1: * It is understood that i-th row of the simplex table is: alpar@1: * alpar@1: * e * A^ = - e * inv(B) * N, (7) alpar@1: * alpar@1: * where e is a unity vector with e[i] = 1. alpar@1: * alpar@1: * To compute i-th row of the simplex table the routine first computes alpar@1: * i-th row of the inverse: alpar@1: * alpar@1: * rho = inv(B') * e, (8) alpar@1: * alpar@1: * where B' is a matrix transposed to B, and then computes elements of alpar@1: * i-th row of the simplex table as scalar products: alpar@1: * alpar@1: * alfa[i,j] = - rho * N[j] for all j, (9) alpar@1: * alpar@1: * where N[j] is a column of the augmented constraint matrix A~, which alpar@1: * corresponds to some non-basic auxiliary or structural variable. */ alpar@1: alpar@1: int glp_eval_tab_row(glp_prob *lp, int k, int ind[], double val[]) alpar@1: { int m = lp->m; alpar@1: int n = lp->n; alpar@1: int i, t, len, lll, *iii; alpar@1: double alfa, *rho, *vvv; alpar@1: if (!(m == 0 || lp->valid)) alpar@1: xerror("glp_eval_tab_row: basis factorization does not exist\n" alpar@1: ); alpar@1: if (!(1 <= k && k <= m+n)) alpar@1: xerror("glp_eval_tab_row: k = %d; variable number out of range" alpar@1: , k); alpar@1: /* determine xB[i] which corresponds to x[k] */ alpar@1: if (k <= m) alpar@1: i = glp_get_row_bind(lp, k); alpar@1: else alpar@1: i = glp_get_col_bind(lp, k-m); alpar@1: if (i == 0) alpar@1: xerror("glp_eval_tab_row: k = %d; variable must be basic", k); alpar@1: xassert(1 <= i && i <= m); alpar@1: /* allocate working arrays */ alpar@1: rho = xcalloc(1+m, sizeof(double)); alpar@1: iii = xcalloc(1+m, sizeof(int)); alpar@1: vvv = xcalloc(1+m, sizeof(double)); alpar@1: /* compute i-th row of the inverse; see (8) */ alpar@1: for (t = 1; t <= m; t++) rho[t] = 0.0; alpar@1: rho[i] = 1.0; alpar@1: glp_btran(lp, rho); alpar@1: /* compute i-th row of the simplex table */ alpar@1: len = 0; alpar@1: for (k = 1; k <= m+n; k++) alpar@1: { if (k <= m) alpar@1: { /* x[k] is auxiliary variable, so N[k] is a unity column */ alpar@1: if (glp_get_row_stat(lp, k) == GLP_BS) continue; alpar@1: /* compute alfa[i,j]; see (9) */ alpar@1: alfa = - rho[k]; alpar@1: } alpar@1: else alpar@1: { /* x[k] is structural variable, so N[k] is a column of the alpar@1: original constraint matrix A with negative sign */ alpar@1: if (glp_get_col_stat(lp, k-m) == GLP_BS) continue; alpar@1: /* compute alfa[i,j]; see (9) */ alpar@1: lll = glp_get_mat_col(lp, k-m, iii, vvv); alpar@1: alfa = 0.0; alpar@1: for (t = 1; t <= lll; t++) alfa += rho[iii[t]] * vvv[t]; alpar@1: } alpar@1: /* store alfa[i,j] */ alpar@1: if (alfa != 0.0) len++, ind[len] = k, val[len] = alfa; alpar@1: } alpar@1: xassert(len <= n); alpar@1: /* free working arrays */ alpar@1: xfree(rho); alpar@1: xfree(iii); alpar@1: xfree(vvv); alpar@1: /* return to the calling program */ alpar@1: return len; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * NAME alpar@1: * alpar@1: * glp_eval_tab_col - compute column of the simplex tableau alpar@1: * alpar@1: * SYNOPSIS alpar@1: * alpar@1: * int glp_eval_tab_col(glp_prob *lp, int k, int ind[], double val[]); alpar@1: * alpar@1: * DESCRIPTION alpar@1: * alpar@1: * The routine glp_eval_tab_col computes a column of the current simplex alpar@1: * table for the non-basic variable, which is specified by the number k: alpar@1: * if 1 <= k <= m, x[k] is k-th auxiliary variable; if m+1 <= k <= m+n, alpar@1: * x[k] is (k-m)-th structural variable, where m is number of rows, and alpar@1: * n is number of columns. The current basis must be available. alpar@1: * alpar@1: * The routine stores row indices and numerical values of non-zero alpar@1: * elements of the computed column using sparse format to the locations alpar@1: * ind[1], ..., ind[len] and val[1], ..., val[len] respectively, where alpar@1: * 0 <= len <= m is number of non-zeros returned on exit. alpar@1: * alpar@1: * Element indices stored in the array ind have the same sense as the alpar@1: * index k, i.e. indices 1 to m denote auxiliary variables and indices alpar@1: * m+1 to m+n denote structural ones (all these variables are obviously alpar@1: * basic by the definition). alpar@1: * alpar@1: * The computed column shows how basic variables depend on the specified alpar@1: * non-basic variable x[k] = xN[j]: alpar@1: * alpar@1: * xB[1] = ... + alfa[1,j]*xN[j] + ... alpar@1: * xB[2] = ... + alfa[2,j]*xN[j] + ... alpar@1: * . . . . . . alpar@1: * xB[m] = ... + alfa[m,j]*xN[j] + ... alpar@1: * alpar@1: * where alfa[i,j] are elements of the simplex table column, xB[i] are alpar@1: * basic (auxiliary and structural) variables. alpar@1: * alpar@1: * RETURNS alpar@1: * alpar@1: * The routine returns number of non-zero elements in the simplex table alpar@1: * column stored in the arrays ind and val. alpar@1: * alpar@1: * BACKGROUND alpar@1: * alpar@1: * As it was explained in comments to the routine glp_eval_tab_row (see alpar@1: * above) the simplex table is the following matrix: alpar@1: * alpar@1: * A^ = - inv(B) * N. (1) alpar@1: * alpar@1: * Therefore j-th column of the simplex table is: alpar@1: * alpar@1: * A^ * e = - inv(B) * N * e = - inv(B) * N[j], (2) alpar@1: * alpar@1: * where e is a unity vector with e[j] = 1, B is the basis matrix, N[j] alpar@1: * is a column of the augmented constraint matrix A~, which corresponds alpar@1: * to the given non-basic auxiliary or structural variable. */ alpar@1: alpar@1: int glp_eval_tab_col(glp_prob *lp, int k, int ind[], double val[]) alpar@1: { int m = lp->m; alpar@1: int n = lp->n; alpar@1: int t, len, stat; alpar@1: double *col; alpar@1: if (!(m == 0 || lp->valid)) alpar@1: xerror("glp_eval_tab_col: basis factorization does not exist\n" alpar@1: ); alpar@1: if (!(1 <= k && k <= m+n)) alpar@1: xerror("glp_eval_tab_col: k = %d; variable number out of range" alpar@1: , k); alpar@1: if (k <= m) alpar@1: stat = glp_get_row_stat(lp, k); alpar@1: else alpar@1: stat = glp_get_col_stat(lp, k-m); alpar@1: if (stat == GLP_BS) alpar@1: xerror("glp_eval_tab_col: k = %d; variable must be non-basic", alpar@1: k); alpar@1: /* obtain column N[k] with negative sign */ alpar@1: col = xcalloc(1+m, sizeof(double)); alpar@1: for (t = 1; t <= m; t++) col[t] = 0.0; alpar@1: if (k <= m) alpar@1: { /* x[k] is auxiliary variable, so N[k] is a unity column */ alpar@1: col[k] = -1.0; alpar@1: } alpar@1: else alpar@1: { /* x[k] is structural variable, so N[k] is a column of the alpar@1: original constraint matrix A with negative sign */ alpar@1: len = glp_get_mat_col(lp, k-m, ind, val); alpar@1: for (t = 1; t <= len; t++) col[ind[t]] = val[t]; alpar@1: } alpar@1: /* compute column of the simplex table, which corresponds to the alpar@1: specified non-basic variable x[k] */ alpar@1: glp_ftran(lp, col); alpar@1: len = 0; alpar@1: for (t = 1; t <= m; t++) alpar@1: { if (col[t] != 0.0) alpar@1: { len++; alpar@1: ind[len] = glp_get_bhead(lp, t); alpar@1: val[len] = col[t]; alpar@1: } alpar@1: } alpar@1: xfree(col); alpar@1: /* return to the calling program */ alpar@1: return len; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * NAME alpar@1: * alpar@1: * glp_transform_row - transform explicitly specified row alpar@1: * alpar@1: * SYNOPSIS alpar@1: * alpar@1: * int glp_transform_row(glp_prob *P, int len, int ind[], double val[]); alpar@1: * alpar@1: * DESCRIPTION alpar@1: * alpar@1: * The routine glp_transform_row performs the same operation as the alpar@1: * routine glp_eval_tab_row with exception that the row to be alpar@1: * transformed is specified explicitly as a sparse vector. alpar@1: * alpar@1: * The explicitly specified row may be thought as a linear form: alpar@1: * alpar@1: * x = a[1]*x[m+1] + a[2]*x[m+2] + ... + a[n]*x[m+n], (1) alpar@1: * alpar@1: * where x is an auxiliary variable for this row, a[j] are coefficients alpar@1: * of the linear form, x[m+j] are structural variables. alpar@1: * alpar@1: * On entry column indices and numerical values of non-zero elements of alpar@1: * the row should be stored in locations ind[1], ..., ind[len] and alpar@1: * val[1], ..., val[len], where len is the number of non-zero elements. alpar@1: * alpar@1: * This routine uses the system of equality constraints and the current alpar@1: * basis in order to express the auxiliary variable x in (1) through the alpar@1: * current non-basic variables (as if the transformed row were added to alpar@1: * the problem object and its auxiliary variable were basic), i.e. the alpar@1: * resultant row has the form: alpar@1: * alpar@1: * x = alfa[1]*xN[1] + alfa[2]*xN[2] + ... + alfa[n]*xN[n], (2) alpar@1: * alpar@1: * where xN[j] are non-basic (auxiliary or structural) variables, n is alpar@1: * the number of columns in the LP problem object. alpar@1: * alpar@1: * On exit the routine stores indices and numerical values of non-zero alpar@1: * elements of the resultant row (2) in locations ind[1], ..., ind[len'] alpar@1: * and val[1], ..., val[len'], where 0 <= len' <= n is the number of alpar@1: * non-zero elements in the resultant row returned by the routine. Note alpar@1: * that indices (numbers) of non-basic variables stored in the array ind alpar@1: * correspond to original ordinal numbers of variables: indices 1 to m alpar@1: * mean auxiliary variables and indices m+1 to m+n mean structural ones. alpar@1: * alpar@1: * RETURNS alpar@1: * alpar@1: * The routine returns len', which is the number of non-zero elements in alpar@1: * the resultant row stored in the arrays ind and val. alpar@1: * alpar@1: * BACKGROUND alpar@1: * alpar@1: * The explicitly specified row (1) is transformed in the same way as it alpar@1: * were the objective function row. alpar@1: * alpar@1: * From (1) it follows that: alpar@1: * alpar@1: * x = aB * xB + aN * xN, (3) alpar@1: * alpar@1: * where xB is the vector of basic variables, xN is the vector of alpar@1: * non-basic variables. alpar@1: * alpar@1: * The simplex table, which corresponds to the current basis, is: alpar@1: * alpar@1: * xB = [-inv(B) * N] * xN. (4) alpar@1: * alpar@1: * Therefore substituting xB from (4) to (3) we have: alpar@1: * alpar@1: * x = aB * [-inv(B) * N] * xN + aN * xN = alpar@1: * (5) alpar@1: * = rho * (-N) * xN + aN * xN = alfa * xN, alpar@1: * alpar@1: * where: alpar@1: * alpar@1: * rho = inv(B') * aB, (6) alpar@1: * alpar@1: * and alpar@1: * alpar@1: * alfa = aN + rho * (-N) (7) alpar@1: * alpar@1: * is the resultant row computed by the routine. */ alpar@1: alpar@1: int glp_transform_row(glp_prob *P, int len, int ind[], double val[]) alpar@1: { int i, j, k, m, n, t, lll, *iii; alpar@1: double alfa, *a, *aB, *rho, *vvv; alpar@1: if (!glp_bf_exists(P)) alpar@1: xerror("glp_transform_row: basis factorization does not exist " alpar@1: "\n"); alpar@1: m = glp_get_num_rows(P); alpar@1: n = glp_get_num_cols(P); alpar@1: /* unpack the row to be transformed to the array a */ alpar@1: a = xcalloc(1+n, sizeof(double)); alpar@1: for (j = 1; j <= n; j++) a[j] = 0.0; alpar@1: if (!(0 <= len && len <= n)) alpar@1: xerror("glp_transform_row: len = %d; invalid row length\n", alpar@1: len); alpar@1: for (t = 1; t <= len; t++) alpar@1: { j = ind[t]; alpar@1: if (!(1 <= j && j <= n)) alpar@1: xerror("glp_transform_row: ind[%d] = %d; column index out o" alpar@1: "f range\n", t, j); alpar@1: if (val[t] == 0.0) alpar@1: xerror("glp_transform_row: val[%d] = 0; zero coefficient no" alpar@1: "t allowed\n", t); alpar@1: if (a[j] != 0.0) alpar@1: xerror("glp_transform_row: ind[%d] = %d; duplicate column i" alpar@1: "ndices not allowed\n", t, j); alpar@1: a[j] = val[t]; alpar@1: } alpar@1: /* construct the vector aB */ alpar@1: aB = xcalloc(1+m, sizeof(double)); alpar@1: for (i = 1; i <= m; i++) alpar@1: { k = glp_get_bhead(P, i); alpar@1: /* xB[i] is k-th original variable */ alpar@1: xassert(1 <= k && k <= m+n); alpar@1: aB[i] = (k <= m ? 0.0 : a[k-m]); alpar@1: } alpar@1: /* solve the system B'*rho = aB to compute the vector rho */ alpar@1: rho = aB, glp_btran(P, rho); alpar@1: /* compute coefficients at non-basic auxiliary variables */ alpar@1: len = 0; alpar@1: for (i = 1; i <= m; i++) alpar@1: { if (glp_get_row_stat(P, i) != GLP_BS) alpar@1: { alfa = - rho[i]; alpar@1: if (alfa != 0.0) alpar@1: { len++; alpar@1: ind[len] = i; alpar@1: val[len] = alfa; alpar@1: } alpar@1: } alpar@1: } alpar@1: /* compute coefficients at non-basic structural variables */ alpar@1: iii = xcalloc(1+m, sizeof(int)); alpar@1: vvv = xcalloc(1+m, sizeof(double)); alpar@1: for (j = 1; j <= n; j++) alpar@1: { if (glp_get_col_stat(P, j) != GLP_BS) alpar@1: { alfa = a[j]; alpar@1: lll = glp_get_mat_col(P, j, iii, vvv); alpar@1: for (t = 1; t <= lll; t++) alfa += vvv[t] * rho[iii[t]]; alpar@1: if (alfa != 0.0) alpar@1: { len++; alpar@1: ind[len] = m+j; alpar@1: val[len] = alfa; alpar@1: } alpar@1: } alpar@1: } alpar@1: xassert(len <= n); alpar@1: xfree(iii); alpar@1: xfree(vvv); alpar@1: xfree(aB); alpar@1: xfree(a); alpar@1: return len; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * NAME alpar@1: * alpar@1: * glp_transform_col - transform explicitly specified column alpar@1: * alpar@1: * SYNOPSIS alpar@1: * alpar@1: * int glp_transform_col(glp_prob *P, int len, int ind[], double val[]); alpar@1: * alpar@1: * DESCRIPTION alpar@1: * alpar@1: * The routine glp_transform_col performs the same operation as the alpar@1: * routine glp_eval_tab_col with exception that the column to be alpar@1: * transformed is specified explicitly as a sparse vector. alpar@1: * alpar@1: * The explicitly specified column may be thought as if it were added alpar@1: * to the original system of equality constraints: alpar@1: * alpar@1: * x[1] = a[1,1]*x[m+1] + ... + a[1,n]*x[m+n] + a[1]*x alpar@1: * x[2] = a[2,1]*x[m+1] + ... + a[2,n]*x[m+n] + a[2]*x (1) alpar@1: * . . . . . . . . . . . . . . . alpar@1: * x[m] = a[m,1]*x[m+1] + ... + a[m,n]*x[m+n] + a[m]*x alpar@1: * alpar@1: * where x[i] are auxiliary variables, x[m+j] are structural variables, alpar@1: * x is a structural variable for the explicitly specified column, a[i] alpar@1: * are constraint coefficients for x. alpar@1: * alpar@1: * On entry row indices and numerical values of non-zero elements of alpar@1: * the column should be stored in locations ind[1], ..., ind[len] and alpar@1: * val[1], ..., val[len], where len is the number of non-zero elements. alpar@1: * alpar@1: * This routine uses the system of equality constraints and the current alpar@1: * basis in order to express the current basic variables through the alpar@1: * structural variable x in (1) (as if the transformed column were added alpar@1: * to the problem object and the variable x were non-basic), i.e. the alpar@1: * resultant column has the form: alpar@1: * alpar@1: * xB[1] = ... + alfa[1]*x alpar@1: * xB[2] = ... + alfa[2]*x (2) alpar@1: * . . . . . . alpar@1: * xB[m] = ... + alfa[m]*x alpar@1: * alpar@1: * where xB are basic (auxiliary and structural) variables, m is the alpar@1: * number of rows in the problem object. alpar@1: * alpar@1: * On exit the routine stores indices and numerical values of non-zero alpar@1: * elements of the resultant column (2) in locations ind[1], ..., alpar@1: * ind[len'] and val[1], ..., val[len'], where 0 <= len' <= m is the alpar@1: * number of non-zero element in the resultant column returned by the alpar@1: * routine. Note that indices (numbers) of basic variables stored in alpar@1: * the array ind correspond to original ordinal numbers of variables: alpar@1: * indices 1 to m mean auxiliary variables and indices m+1 to m+n mean alpar@1: * structural ones. alpar@1: * alpar@1: * RETURNS alpar@1: * alpar@1: * The routine returns len', which is the number of non-zero elements alpar@1: * in the resultant column stored in the arrays ind and val. alpar@1: * alpar@1: * BACKGROUND alpar@1: * alpar@1: * The explicitly specified column (1) is transformed in the same way alpar@1: * as any other column of the constraint matrix using the formula: alpar@1: * alpar@1: * alfa = inv(B) * a, (3) alpar@1: * alpar@1: * where alfa is the resultant column computed by the routine. */ alpar@1: alpar@1: int glp_transform_col(glp_prob *P, int len, int ind[], double val[]) alpar@1: { int i, m, t; alpar@1: double *a, *alfa; alpar@1: if (!glp_bf_exists(P)) alpar@1: xerror("glp_transform_col: basis factorization does not exist " alpar@1: "\n"); alpar@1: m = glp_get_num_rows(P); alpar@1: /* unpack the column to be transformed to the array a */ alpar@1: a = xcalloc(1+m, sizeof(double)); alpar@1: for (i = 1; i <= m; i++) a[i] = 0.0; alpar@1: if (!(0 <= len && len <= m)) alpar@1: xerror("glp_transform_col: len = %d; invalid column length\n", alpar@1: len); alpar@1: for (t = 1; t <= len; t++) alpar@1: { i = ind[t]; alpar@1: if (!(1 <= i && i <= m)) alpar@1: xerror("glp_transform_col: ind[%d] = %d; row index out of r" alpar@1: "ange\n", t, i); alpar@1: if (val[t] == 0.0) alpar@1: xerror("glp_transform_col: val[%d] = 0; zero coefficient no" alpar@1: "t allowed\n", t); alpar@1: if (a[i] != 0.0) alpar@1: xerror("glp_transform_col: ind[%d] = %d; duplicate row indi" alpar@1: "ces not allowed\n", t, i); alpar@1: a[i] = val[t]; alpar@1: } alpar@1: /* solve the system B*a = alfa to compute the vector alfa */ alpar@1: alfa = a, glp_ftran(P, alfa); alpar@1: /* store resultant coefficients */ alpar@1: len = 0; alpar@1: for (i = 1; i <= m; i++) alpar@1: { if (alfa[i] != 0.0) alpar@1: { len++; alpar@1: ind[len] = glp_get_bhead(P, i); alpar@1: val[len] = alfa[i]; alpar@1: } alpar@1: } alpar@1: xfree(a); alpar@1: return len; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * NAME alpar@1: * alpar@1: * glp_prim_rtest - perform primal ratio test alpar@1: * alpar@1: * SYNOPSIS alpar@1: * alpar@1: * int glp_prim_rtest(glp_prob *P, int len, const int ind[], alpar@1: * const double val[], int dir, double eps); alpar@1: * alpar@1: * DESCRIPTION alpar@1: * alpar@1: * The routine glp_prim_rtest performs the primal ratio test using an alpar@1: * explicitly specified column of the simplex table. alpar@1: * alpar@1: * The current basic solution associated with the LP problem object alpar@1: * must be primal feasible. alpar@1: * alpar@1: * The explicitly specified column of the simplex table shows how the alpar@1: * basic variables xB depend on some non-basic variable x (which is not alpar@1: * necessarily presented in the problem object): alpar@1: * alpar@1: * xB[1] = ... + alfa[1] * x + ... alpar@1: * xB[2] = ... + alfa[2] * x + ... (*) alpar@1: * . . . . . . . . alpar@1: * xB[m] = ... + alfa[m] * x + ... alpar@1: * alpar@1: * The column (*) is specifed on entry to the routine using the sparse alpar@1: * format. Ordinal numbers of basic variables xB[i] should be placed in alpar@1: * locations ind[1], ..., ind[len], where ordinal number 1 to m denote alpar@1: * auxiliary variables, and ordinal numbers m+1 to m+n denote structural alpar@1: * variables. The corresponding non-zero coefficients alfa[i] should be alpar@1: * placed in locations val[1], ..., val[len]. The arrays ind and val are alpar@1: * not changed on exit. alpar@1: * alpar@1: * The parameter dir specifies direction in which the variable x changes alpar@1: * on entering the basis: +1 means increasing, -1 means decreasing. alpar@1: * alpar@1: * The parameter eps is an absolute tolerance (small positive number) alpar@1: * used by the routine to skip small alfa[j] of the row (*). alpar@1: * alpar@1: * The routine determines which basic variable (among specified in alpar@1: * ind[1], ..., ind[len]) should leave the basis in order to keep primal alpar@1: * feasibility. alpar@1: * alpar@1: * RETURNS alpar@1: * alpar@1: * The routine glp_prim_rtest returns the index piv in the arrays ind alpar@1: * and val corresponding to the pivot element chosen, 1 <= piv <= len. alpar@1: * If the adjacent basic solution is primal unbounded and therefore the alpar@1: * choice cannot be made, the routine returns zero. alpar@1: * alpar@1: * COMMENTS alpar@1: * alpar@1: * If the non-basic variable x is presented in the LP problem object, alpar@1: * the column (*) can be computed with the routine glp_eval_tab_col; alpar@1: * otherwise it can be computed with the routine glp_transform_col. */ alpar@1: alpar@1: int glp_prim_rtest(glp_prob *P, int len, const int ind[], alpar@1: const double val[], int dir, double eps) alpar@1: { int k, m, n, piv, t, type, stat; alpar@1: double alfa, big, beta, lb, ub, temp, teta; alpar@1: if (glp_get_prim_stat(P) != GLP_FEAS) alpar@1: xerror("glp_prim_rtest: basic solution is not primal feasible " alpar@1: "\n"); alpar@1: if (!(dir == +1 || dir == -1)) alpar@1: xerror("glp_prim_rtest: dir = %d; invalid parameter\n", dir); alpar@1: if (!(0.0 < eps && eps < 1.0)) alpar@1: xerror("glp_prim_rtest: eps = %g; invalid parameter\n", eps); alpar@1: m = glp_get_num_rows(P); alpar@1: n = glp_get_num_cols(P); alpar@1: /* initial settings */ alpar@1: piv = 0, teta = DBL_MAX, big = 0.0; alpar@1: /* walk through the entries of the specified column */ alpar@1: for (t = 1; t <= len; t++) alpar@1: { /* get the ordinal number of basic variable */ alpar@1: k = ind[t]; alpar@1: if (!(1 <= k && k <= m+n)) alpar@1: xerror("glp_prim_rtest: ind[%d] = %d; variable number out o" alpar@1: "f range\n", t, k); alpar@1: /* determine type, bounds, status and primal value of basic alpar@1: variable xB[i] = x[k] in the current basic solution */ alpar@1: if (k <= m) alpar@1: { type = glp_get_row_type(P, k); alpar@1: lb = glp_get_row_lb(P, k); alpar@1: ub = glp_get_row_ub(P, k); alpar@1: stat = glp_get_row_stat(P, k); alpar@1: beta = glp_get_row_prim(P, k); alpar@1: } alpar@1: else alpar@1: { type = glp_get_col_type(P, k-m); alpar@1: lb = glp_get_col_lb(P, k-m); alpar@1: ub = glp_get_col_ub(P, k-m); alpar@1: stat = glp_get_col_stat(P, k-m); alpar@1: beta = glp_get_col_prim(P, k-m); alpar@1: } alpar@1: if (stat != GLP_BS) alpar@1: xerror("glp_prim_rtest: ind[%d] = %d; non-basic variable no" alpar@1: "t allowed\n", t, k); alpar@1: /* determine influence coefficient at basic variable xB[i] alpar@1: in the explicitly specified column and turn to the case of alpar@1: increasing the variable x in order to simplify the program alpar@1: logic */ alpar@1: alfa = (dir > 0 ? + val[t] : - val[t]); alpar@1: /* analyze main cases */ alpar@1: if (type == GLP_FR) alpar@1: { /* xB[i] is free variable */ alpar@1: continue; alpar@1: } alpar@1: else if (type == GLP_LO) alpar@1: lo: { /* xB[i] has an lower bound */ alpar@1: if (alfa > - eps) continue; alpar@1: temp = (lb - beta) / alfa; alpar@1: } alpar@1: else if (type == GLP_UP) alpar@1: up: { /* xB[i] has an upper bound */ alpar@1: if (alfa < + eps) continue; alpar@1: temp = (ub - beta) / alfa; alpar@1: } alpar@1: else if (type == GLP_DB) alpar@1: { /* xB[i] has both lower and upper bounds */ alpar@1: if (alfa < 0.0) goto lo; else goto up; alpar@1: } alpar@1: else if (type == GLP_FX) alpar@1: { /* xB[i] is fixed variable */ alpar@1: if (- eps < alfa && alfa < + eps) continue; alpar@1: temp = 0.0; alpar@1: } alpar@1: else alpar@1: xassert(type != type); alpar@1: /* if the value of the variable xB[i] violates its lower or alpar@1: upper bound (slightly, because the current basis is assumed alpar@1: to be primal feasible), temp is negative; we can think this alpar@1: happens due to round-off errors and the value is exactly on alpar@1: the bound; this allows replacing temp by zero */ alpar@1: if (temp < 0.0) temp = 0.0; alpar@1: /* apply the minimal ratio test */ alpar@1: if (teta > temp || teta == temp && big < fabs(alfa)) alpar@1: piv = t, teta = temp, big = fabs(alfa); alpar@1: } alpar@1: /* return index of the pivot element chosen */ alpar@1: return piv; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * NAME alpar@1: * alpar@1: * glp_dual_rtest - perform dual ratio test alpar@1: * alpar@1: * SYNOPSIS alpar@1: * alpar@1: * int glp_dual_rtest(glp_prob *P, int len, const int ind[], alpar@1: * const double val[], int dir, double eps); alpar@1: * alpar@1: * DESCRIPTION alpar@1: * alpar@1: * The routine glp_dual_rtest performs the dual ratio test using an alpar@1: * explicitly specified row of the simplex table. alpar@1: * alpar@1: * The current basic solution associated with the LP problem object alpar@1: * must be dual feasible. alpar@1: * alpar@1: * The explicitly specified row of the simplex table is a linear form alpar@1: * that shows how some basic variable x (which is not necessarily alpar@1: * presented in the problem object) depends on non-basic variables xN: alpar@1: * alpar@1: * x = alfa[1] * xN[1] + alfa[2] * xN[2] + ... + alfa[n] * xN[n]. (*) alpar@1: * alpar@1: * The row (*) is specified on entry to the routine using the sparse alpar@1: * format. Ordinal numbers of non-basic variables xN[j] should be placed alpar@1: * in locations ind[1], ..., ind[len], where ordinal numbers 1 to m alpar@1: * denote auxiliary variables, and ordinal numbers m+1 to m+n denote alpar@1: * structural variables. The corresponding non-zero coefficients alfa[j] alpar@1: * should be placed in locations val[1], ..., val[len]. The arrays ind alpar@1: * and val are not changed on exit. alpar@1: * alpar@1: * The parameter dir specifies direction in which the variable x changes alpar@1: * on leaving the basis: +1 means that x goes to its lower bound, and -1 alpar@1: * means that x goes to its upper bound. alpar@1: * alpar@1: * The parameter eps is an absolute tolerance (small positive number) alpar@1: * used by the routine to skip small alfa[j] of the row (*). alpar@1: * alpar@1: * The routine determines which non-basic variable (among specified in alpar@1: * ind[1], ..., ind[len]) should enter the basis in order to keep dual alpar@1: * feasibility. alpar@1: * alpar@1: * RETURNS alpar@1: * alpar@1: * The routine glp_dual_rtest returns the index piv in the arrays ind alpar@1: * and val corresponding to the pivot element chosen, 1 <= piv <= len. alpar@1: * If the adjacent basic solution is dual unbounded and therefore the alpar@1: * choice cannot be made, the routine returns zero. alpar@1: * alpar@1: * COMMENTS alpar@1: * alpar@1: * If the basic variable x is presented in the LP problem object, the alpar@1: * row (*) can be computed with the routine glp_eval_tab_row; otherwise alpar@1: * it can be computed with the routine glp_transform_row. */ alpar@1: alpar@1: int glp_dual_rtest(glp_prob *P, int len, const int ind[], alpar@1: const double val[], int dir, double eps) alpar@1: { int k, m, n, piv, t, stat; alpar@1: double alfa, big, cost, obj, temp, teta; alpar@1: if (glp_get_dual_stat(P) != GLP_FEAS) alpar@1: xerror("glp_dual_rtest: basic solution is not dual feasible\n") alpar@1: ; alpar@1: if (!(dir == +1 || dir == -1)) alpar@1: xerror("glp_dual_rtest: dir = %d; invalid parameter\n", dir); alpar@1: if (!(0.0 < eps && eps < 1.0)) alpar@1: xerror("glp_dual_rtest: eps = %g; invalid parameter\n", eps); alpar@1: m = glp_get_num_rows(P); alpar@1: n = glp_get_num_cols(P); alpar@1: /* take into account optimization direction */ alpar@1: obj = (glp_get_obj_dir(P) == GLP_MIN ? +1.0 : -1.0); alpar@1: /* initial settings */ alpar@1: piv = 0, teta = DBL_MAX, big = 0.0; alpar@1: /* walk through the entries of the specified row */ alpar@1: for (t = 1; t <= len; t++) alpar@1: { /* get ordinal number of non-basic variable */ alpar@1: k = ind[t]; alpar@1: if (!(1 <= k && k <= m+n)) alpar@1: xerror("glp_dual_rtest: ind[%d] = %d; variable number out o" alpar@1: "f range\n", t, k); alpar@1: /* determine status and reduced cost of non-basic variable alpar@1: x[k] = xN[j] in the current basic solution */ alpar@1: if (k <= m) alpar@1: { stat = glp_get_row_stat(P, k); alpar@1: cost = glp_get_row_dual(P, k); alpar@1: } alpar@1: else alpar@1: { stat = glp_get_col_stat(P, k-m); alpar@1: cost = glp_get_col_dual(P, k-m); alpar@1: } alpar@1: if (stat == GLP_BS) alpar@1: xerror("glp_dual_rtest: ind[%d] = %d; basic variable not al" alpar@1: "lowed\n", t, k); alpar@1: /* determine influence coefficient at non-basic variable xN[j] alpar@1: in the explicitly specified row and turn to the case of alpar@1: increasing the variable x in order to simplify the program alpar@1: logic */ alpar@1: alfa = (dir > 0 ? + val[t] : - val[t]); alpar@1: /* analyze main cases */ alpar@1: if (stat == GLP_NL) alpar@1: { /* xN[j] is on its lower bound */ alpar@1: if (alfa < + eps) continue; alpar@1: temp = (obj * cost) / alfa; alpar@1: } alpar@1: else if (stat == GLP_NU) alpar@1: { /* xN[j] is on its upper bound */ alpar@1: if (alfa > - eps) continue; alpar@1: temp = (obj * cost) / alfa; alpar@1: } alpar@1: else if (stat == GLP_NF) alpar@1: { /* xN[j] is non-basic free variable */ alpar@1: if (- eps < alfa && alfa < + eps) continue; alpar@1: temp = 0.0; alpar@1: } alpar@1: else if (stat == GLP_NS) alpar@1: { /* xN[j] is non-basic fixed variable */ alpar@1: continue; alpar@1: } alpar@1: else alpar@1: xassert(stat != stat); alpar@1: /* if the reduced cost of the variable xN[j] violates its zero alpar@1: bound (slightly, because the current basis is assumed to be alpar@1: dual feasible), temp is negative; we can think this happens alpar@1: due to round-off errors and the reduced cost is exact zero; alpar@1: this allows replacing temp by zero */ alpar@1: if (temp < 0.0) temp = 0.0; alpar@1: /* apply the minimal ratio test */ alpar@1: if (teta > temp || teta == temp && big < fabs(alfa)) alpar@1: piv = t, teta = temp, big = fabs(alfa); alpar@1: } alpar@1: /* return index of the pivot element chosen */ alpar@1: return piv; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * NAME alpar@1: * alpar@1: * glp_analyze_row - simulate one iteration of dual simplex method alpar@1: * alpar@1: * SYNOPSIS alpar@1: * alpar@1: * int glp_analyze_row(glp_prob *P, int len, const int ind[], alpar@1: * const double val[], int type, double rhs, double eps, int *piv, alpar@1: * double *x, double *dx, double *y, double *dy, double *dz); alpar@1: * alpar@1: * DESCRIPTION alpar@1: * alpar@1: * Let the current basis be optimal or dual feasible, and there be alpar@1: * specified a row (constraint), which is violated by the current basic alpar@1: * solution. The routine glp_analyze_row simulates one iteration of the alpar@1: * dual simplex method to determine some information on the adjacent alpar@1: * basis (see below), where the specified row becomes active constraint alpar@1: * (i.e. its auxiliary variable becomes non-basic). alpar@1: * alpar@1: * The current basic solution associated with the problem object passed alpar@1: * to the routine must be dual feasible, and its primal components must alpar@1: * be defined. alpar@1: * alpar@1: * The row to be analyzed must be previously transformed either with alpar@1: * the routine glp_eval_tab_row (if the row is in the problem object) alpar@1: * or with the routine glp_transform_row (if the row is external, i.e. alpar@1: * not in the problem object). This is needed to express the row only alpar@1: * through (auxiliary and structural) variables, which are non-basic in alpar@1: * the current basis: alpar@1: * alpar@1: * y = alfa[1] * xN[1] + alfa[2] * xN[2] + ... + alfa[n] * xN[n], alpar@1: * alpar@1: * where y is an auxiliary variable of the row, alfa[j] is an influence alpar@1: * coefficient, xN[j] is a non-basic variable. alpar@1: * alpar@1: * The row is passed to the routine in sparse format. Ordinal numbers alpar@1: * of non-basic variables are stored in locations ind[1], ..., ind[len], alpar@1: * where numbers 1 to m denote auxiliary variables while numbers m+1 to alpar@1: * m+n denote structural variables. Corresponding non-zero coefficients alpar@1: * alfa[j] are stored in locations val[1], ..., val[len]. The arrays alpar@1: * ind and val are ot changed on exit. alpar@1: * alpar@1: * The parameters type and rhs specify the row type and its right-hand alpar@1: * side as follows: alpar@1: * alpar@1: * type = GLP_LO: y = sum alfa[j] * xN[j] >= rhs alpar@1: * alpar@1: * type = GLP_UP: y = sum alfa[j] * xN[j] <= rhs alpar@1: * alpar@1: * The parameter eps is an absolute tolerance (small positive number) alpar@1: * used by the routine to skip small coefficients alfa[j] on performing alpar@1: * the dual ratio test. alpar@1: * alpar@1: * If the operation was successful, the routine stores the following alpar@1: * information to corresponding location (if some parameter is NULL, alpar@1: * its value is not stored): alpar@1: * alpar@1: * piv index in the array ind and val, 1 <= piv <= len, determining alpar@1: * the non-basic variable, which would enter the adjacent basis; alpar@1: * alpar@1: * x value of the non-basic variable in the current basis; alpar@1: * alpar@1: * dx difference between values of the non-basic variable in the alpar@1: * adjacent and current bases, dx = x.new - x.old; alpar@1: * alpar@1: * y value of the row (i.e. of its auxiliary variable) in the alpar@1: * current basis; alpar@1: * alpar@1: * dy difference between values of the row in the adjacent and alpar@1: * current bases, dy = y.new - y.old; alpar@1: * alpar@1: * dz difference between values of the objective function in the alpar@1: * adjacent and current bases, dz = z.new - z.old. Note that in alpar@1: * case of minimization dz >= 0, and in case of maximization alpar@1: * dz <= 0, i.e. in the adjacent basis the objective function alpar@1: * always gets worse (degrades). */ alpar@1: alpar@1: int _glp_analyze_row(glp_prob *P, int len, const int ind[], alpar@1: const double val[], int type, double rhs, double eps, int *_piv, alpar@1: double *_x, double *_dx, double *_y, double *_dy, double *_dz) alpar@1: { int t, k, dir, piv, ret = 0; alpar@1: double x, dx, y, dy, dz; alpar@1: if (P->pbs_stat == GLP_UNDEF) alpar@1: xerror("glp_analyze_row: primal basic solution components are " alpar@1: "undefined\n"); alpar@1: if (P->dbs_stat != GLP_FEAS) alpar@1: xerror("glp_analyze_row: basic solution is not dual feasible\n" alpar@1: ); alpar@1: /* compute the row value y = sum alfa[j] * xN[j] in the current alpar@1: basis */ alpar@1: if (!(0 <= len && len <= P->n)) alpar@1: xerror("glp_analyze_row: len = %d; invalid row length\n", len); alpar@1: y = 0.0; alpar@1: for (t = 1; t <= len; t++) alpar@1: { /* determine value of x[k] = xN[j] in the current basis */ alpar@1: k = ind[t]; alpar@1: if (!(1 <= k && k <= P->m+P->n)) alpar@1: xerror("glp_analyze_row: ind[%d] = %d; row/column index out" alpar@1: " of range\n", t, k); alpar@1: if (k <= P->m) alpar@1: { /* x[k] is auxiliary variable */ alpar@1: if (P->row[k]->stat == GLP_BS) alpar@1: xerror("glp_analyze_row: ind[%d] = %d; basic auxiliary v" alpar@1: "ariable is not allowed\n", t, k); alpar@1: x = P->row[k]->prim; alpar@1: } alpar@1: else alpar@1: { /* x[k] is structural variable */ alpar@1: if (P->col[k-P->m]->stat == GLP_BS) alpar@1: xerror("glp_analyze_row: ind[%d] = %d; basic structural " alpar@1: "variable is not allowed\n", t, k); alpar@1: x = P->col[k-P->m]->prim; alpar@1: } alpar@1: y += val[t] * x; alpar@1: } alpar@1: /* check if the row is primal infeasible in the current basis, alpar@1: i.e. the constraint is violated at the current point */ alpar@1: if (type == GLP_LO) alpar@1: { if (y >= rhs) alpar@1: { /* the constraint is not violated */ alpar@1: ret = 1; alpar@1: goto done; alpar@1: } alpar@1: /* in the adjacent basis y goes to its lower bound */ alpar@1: dir = +1; alpar@1: } alpar@1: else if (type == GLP_UP) alpar@1: { if (y <= rhs) alpar@1: { /* the constraint is not violated */ alpar@1: ret = 1; alpar@1: goto done; alpar@1: } alpar@1: /* in the adjacent basis y goes to its upper bound */ alpar@1: dir = -1; alpar@1: } alpar@1: else alpar@1: xerror("glp_analyze_row: type = %d; invalid parameter\n", alpar@1: type); alpar@1: /* compute dy = y.new - y.old */ alpar@1: dy = rhs - y; alpar@1: /* perform dual ratio test to determine which non-basic variable alpar@1: should enter the adjacent basis to keep it dual feasible */ alpar@1: piv = glp_dual_rtest(P, len, ind, val, dir, eps); alpar@1: if (piv == 0) alpar@1: { /* no dual feasible adjacent basis exists */ alpar@1: ret = 2; alpar@1: goto done; alpar@1: } alpar@1: /* non-basic variable x[k] = xN[j] should enter the basis */ alpar@1: k = ind[piv]; alpar@1: xassert(1 <= k && k <= P->m+P->n); alpar@1: /* determine its value in the current basis */ alpar@1: if (k <= P->m) alpar@1: x = P->row[k]->prim; alpar@1: else alpar@1: x = P->col[k-P->m]->prim; alpar@1: /* compute dx = x.new - x.old = dy / alfa[j] */ alpar@1: xassert(val[piv] != 0.0); alpar@1: dx = dy / val[piv]; alpar@1: /* compute dz = z.new - z.old = d[j] * dx, where d[j] is reduced alpar@1: cost of xN[j] in the current basis */ alpar@1: if (k <= P->m) alpar@1: dz = P->row[k]->dual * dx; alpar@1: else alpar@1: dz = P->col[k-P->m]->dual * dx; alpar@1: /* store the analysis results */ alpar@1: if (_piv != NULL) *_piv = piv; alpar@1: if (_x != NULL) *_x = x; alpar@1: if (_dx != NULL) *_dx = dx; alpar@1: if (_y != NULL) *_y = y; alpar@1: if (_dy != NULL) *_dy = dy; alpar@1: if (_dz != NULL) *_dz = dz; alpar@1: done: return ret; alpar@1: } alpar@1: alpar@1: #if 0 alpar@1: int main(void) alpar@1: { /* example program for the routine glp_analyze_row */ alpar@1: glp_prob *P; alpar@1: glp_smcp parm; alpar@1: int i, k, len, piv, ret, ind[1+100]; alpar@1: double rhs, x, dx, y, dy, dz, val[1+100]; alpar@1: P = glp_create_prob(); alpar@1: /* read plan.mps (see glpk/examples) */ alpar@1: ret = glp_read_mps(P, GLP_MPS_DECK, NULL, "plan.mps"); alpar@1: glp_assert(ret == 0); alpar@1: /* and solve it to optimality */ alpar@1: ret = glp_simplex(P, NULL); alpar@1: glp_assert(ret == 0); alpar@1: glp_assert(glp_get_status(P) == GLP_OPT); alpar@1: /* the optimal objective value is 296.217 */ alpar@1: /* we would like to know what happens if we would add a new row alpar@1: (constraint) to plan.mps: alpar@1: .01 * bin1 + .01 * bin2 + .02 * bin4 + .02 * bin5 <= 12 */ alpar@1: /* first, we specify this new row */ alpar@1: glp_create_index(P); alpar@1: len = 0; alpar@1: ind[++len] = glp_find_col(P, "BIN1"), val[len] = .01; alpar@1: ind[++len] = glp_find_col(P, "BIN2"), val[len] = .01; alpar@1: ind[++len] = glp_find_col(P, "BIN4"), val[len] = .02; alpar@1: ind[++len] = glp_find_col(P, "BIN5"), val[len] = .02; alpar@1: rhs = 12; alpar@1: /* then we can compute value of the row (i.e. of its auxiliary alpar@1: variable) in the current basis to see if the constraint is alpar@1: violated */ alpar@1: y = 0.0; alpar@1: for (k = 1; k <= len; k++) alpar@1: y += val[k] * glp_get_col_prim(P, ind[k]); alpar@1: glp_printf("y = %g\n", y); alpar@1: /* this prints y = 15.1372, so the constraint is violated, since alpar@1: we require that y <= rhs = 12 */ alpar@1: /* now we transform the row to express it only through non-basic alpar@1: (auxiliary and artificial) variables */ alpar@1: len = glp_transform_row(P, len, ind, val); alpar@1: /* finally, we simulate one step of the dual simplex method to alpar@1: obtain necessary information for the adjacent basis */ alpar@1: ret = _glp_analyze_row(P, len, ind, val, GLP_UP, rhs, 1e-9, &piv, alpar@1: &x, &dx, &y, &dy, &dz); alpar@1: glp_assert(ret == 0); alpar@1: glp_printf("k = %d, x = %g; dx = %g; y = %g; dy = %g; dz = %g\n", alpar@1: ind[piv], x, dx, y, dy, dz); alpar@1: /* this prints dz = 5.64418 and means that in the adjacent basis alpar@1: the objective function would be 296.217 + 5.64418 = 301.861 */ alpar@1: /* now we actually include the row into the problem object; note alpar@1: that the arrays ind and val are clobbered, so we need to build alpar@1: them once again */ alpar@1: len = 0; alpar@1: ind[++len] = glp_find_col(P, "BIN1"), val[len] = .01; alpar@1: ind[++len] = glp_find_col(P, "BIN2"), val[len] = .01; alpar@1: ind[++len] = glp_find_col(P, "BIN4"), val[len] = .02; alpar@1: ind[++len] = glp_find_col(P, "BIN5"), val[len] = .02; alpar@1: rhs = 12; alpar@1: i = glp_add_rows(P, 1); alpar@1: glp_set_row_bnds(P, i, GLP_UP, 0, rhs); alpar@1: glp_set_mat_row(P, i, len, ind, val); alpar@1: /* and perform one dual simplex iteration */ alpar@1: glp_init_smcp(&parm); alpar@1: parm.meth = GLP_DUAL; alpar@1: parm.it_lim = 1; alpar@1: glp_simplex(P, &parm); alpar@1: /* the current objective value is 301.861 */ alpar@1: return 0; alpar@1: } alpar@1: #endif alpar@1: alpar@1: /*********************************************************************** alpar@1: * NAME alpar@1: * alpar@1: * glp_analyze_bound - analyze active bound of non-basic variable alpar@1: * alpar@1: * SYNOPSIS alpar@1: * alpar@1: * void glp_analyze_bound(glp_prob *P, int k, double *limit1, int *var1, alpar@1: * double *limit2, int *var2); alpar@1: * alpar@1: * DESCRIPTION alpar@1: * alpar@1: * The routine glp_analyze_bound analyzes the effect of varying the alpar@1: * active bound of specified non-basic variable. alpar@1: * alpar@1: * The non-basic variable is specified by the parameter k, where alpar@1: * 1 <= k <= m means auxiliary variable of corresponding row while alpar@1: * m+1 <= k <= m+n means structural variable (column). alpar@1: * alpar@1: * Note that the current basic solution must be optimal, and the basis alpar@1: * factorization must exist. alpar@1: * alpar@1: * Results of the analysis have the following meaning. alpar@1: * alpar@1: * value1 is the minimal value of the active bound, at which the basis alpar@1: * still remains primal feasible and thus optimal. -DBL_MAX means that alpar@1: * the active bound has no lower limit. alpar@1: * alpar@1: * var1 is the ordinal number of an auxiliary (1 to m) or structural alpar@1: * (m+1 to n) basic variable, which reaches its bound first and thereby alpar@1: * limits further decreasing the active bound being analyzed. alpar@1: * if value1 = -DBL_MAX, var1 is set to 0. alpar@1: * alpar@1: * value2 is the maximal value of the active bound, at which the basis alpar@1: * still remains primal feasible and thus optimal. +DBL_MAX means that alpar@1: * the active bound has no upper limit. alpar@1: * alpar@1: * var2 is the ordinal number of an auxiliary (1 to m) or structural alpar@1: * (m+1 to n) basic variable, which reaches its bound first and thereby alpar@1: * limits further increasing the active bound being analyzed. alpar@1: * if value2 = +DBL_MAX, var2 is set to 0. */ alpar@1: alpar@1: void glp_analyze_bound(glp_prob *P, int k, double *value1, int *var1, alpar@1: double *value2, int *var2) alpar@1: { GLPROW *row; alpar@1: GLPCOL *col; alpar@1: int m, n, stat, kase, p, len, piv, *ind; alpar@1: double x, new_x, ll, uu, xx, delta, *val; alpar@1: /* sanity checks */ alpar@1: if (P == NULL || P->magic != GLP_PROB_MAGIC) alpar@1: xerror("glp_analyze_bound: P = %p; invalid problem object\n", alpar@1: P); alpar@1: m = P->m, n = P->n; alpar@1: if (!(P->pbs_stat == GLP_FEAS && P->dbs_stat == GLP_FEAS)) alpar@1: xerror("glp_analyze_bound: optimal basic solution required\n"); alpar@1: if (!(m == 0 || P->valid)) alpar@1: xerror("glp_analyze_bound: basis factorization required\n"); alpar@1: if (!(1 <= k && k <= m+n)) alpar@1: xerror("glp_analyze_bound: k = %d; variable number out of rang" alpar@1: "e\n", k); alpar@1: /* retrieve information about the specified non-basic variable alpar@1: x[k] whose active bound is to be analyzed */ alpar@1: if (k <= m) alpar@1: { row = P->row[k]; alpar@1: stat = row->stat; alpar@1: x = row->prim; alpar@1: } alpar@1: else alpar@1: { col = P->col[k-m]; alpar@1: stat = col->stat; alpar@1: x = col->prim; alpar@1: } alpar@1: if (stat == GLP_BS) alpar@1: xerror("glp_analyze_bound: k = %d; basic variable not allowed " alpar@1: "\n", k); alpar@1: /* allocate working arrays */ alpar@1: ind = xcalloc(1+m, sizeof(int)); alpar@1: val = xcalloc(1+m, sizeof(double)); alpar@1: /* compute column of the simplex table corresponding to the alpar@1: non-basic variable x[k] */ alpar@1: len = glp_eval_tab_col(P, k, ind, val); alpar@1: xassert(0 <= len && len <= m); alpar@1: /* perform analysis */ alpar@1: for (kase = -1; kase <= +1; kase += 2) alpar@1: { /* kase < 0 means active bound of x[k] is decreasing; alpar@1: kase > 0 means active bound of x[k] is increasing */ alpar@1: /* use the primal ratio test to determine some basic variable alpar@1: x[p] which reaches its bound first */ alpar@1: piv = glp_prim_rtest(P, len, ind, val, kase, 1e-9); alpar@1: if (piv == 0) alpar@1: { /* nothing limits changing the active bound of x[k] */ alpar@1: p = 0; alpar@1: new_x = (kase < 0 ? -DBL_MAX : +DBL_MAX); alpar@1: goto store; alpar@1: } alpar@1: /* basic variable x[p] limits changing the active bound of alpar@1: x[k]; determine its value in the current basis */ alpar@1: xassert(1 <= piv && piv <= len); alpar@1: p = ind[piv]; alpar@1: if (p <= m) alpar@1: { row = P->row[p]; alpar@1: ll = glp_get_row_lb(P, row->i); alpar@1: uu = glp_get_row_ub(P, row->i); alpar@1: stat = row->stat; alpar@1: xx = row->prim; alpar@1: } alpar@1: else alpar@1: { col = P->col[p-m]; alpar@1: ll = glp_get_col_lb(P, col->j); alpar@1: uu = glp_get_col_ub(P, col->j); alpar@1: stat = col->stat; alpar@1: xx = col->prim; alpar@1: } alpar@1: xassert(stat == GLP_BS); alpar@1: /* determine delta x[p] = bound of x[p] - value of x[p] */ alpar@1: if (kase < 0 && val[piv] > 0.0 || alpar@1: kase > 0 && val[piv] < 0.0) alpar@1: { /* delta x[p] < 0, so x[p] goes toward its lower bound */ alpar@1: xassert(ll != -DBL_MAX); alpar@1: delta = ll - xx; alpar@1: } alpar@1: else alpar@1: { /* delta x[p] > 0, so x[p] goes toward its upper bound */ alpar@1: xassert(uu != +DBL_MAX); alpar@1: delta = uu - xx; alpar@1: } alpar@1: /* delta x[p] = alfa[p,k] * delta x[k], so new x[k] = x[k] + alpar@1: delta x[k] = x[k] + delta x[p] / alfa[p,k] is the value of alpar@1: x[k] in the adjacent basis */ alpar@1: xassert(val[piv] != 0.0); alpar@1: new_x = x + delta / val[piv]; alpar@1: store: /* store analysis results */ alpar@1: if (kase < 0) alpar@1: { if (value1 != NULL) *value1 = new_x; alpar@1: if (var1 != NULL) *var1 = p; alpar@1: } alpar@1: else alpar@1: { if (value2 != NULL) *value2 = new_x; alpar@1: if (var2 != NULL) *var2 = p; alpar@1: } alpar@1: } alpar@1: /* free working arrays */ alpar@1: xfree(ind); alpar@1: xfree(val); alpar@1: return; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * NAME alpar@1: * alpar@1: * glp_analyze_coef - analyze objective coefficient at basic variable alpar@1: * alpar@1: * SYNOPSIS alpar@1: * alpar@1: * void glp_analyze_coef(glp_prob *P, int k, double *coef1, int *var1, alpar@1: * double *value1, double *coef2, int *var2, double *value2); alpar@1: * alpar@1: * DESCRIPTION alpar@1: * alpar@1: * The routine glp_analyze_coef analyzes the effect of varying the alpar@1: * objective coefficient at specified basic variable. alpar@1: * alpar@1: * The basic variable is specified by the parameter k, where alpar@1: * 1 <= k <= m means auxiliary variable of corresponding row while alpar@1: * m+1 <= k <= m+n means structural variable (column). alpar@1: * alpar@1: * Note that the current basic solution must be optimal, and the basis alpar@1: * factorization must exist. alpar@1: * alpar@1: * Results of the analysis have the following meaning. alpar@1: * alpar@1: * coef1 is the minimal value of the objective coefficient, at which alpar@1: * the basis still remains dual feasible and thus optimal. -DBL_MAX alpar@1: * means that the objective coefficient has no lower limit. alpar@1: * alpar@1: * var1 is the ordinal number of an auxiliary (1 to m) or structural alpar@1: * (m+1 to n) non-basic variable, whose reduced cost reaches its zero alpar@1: * bound first and thereby limits further decreasing the objective alpar@1: * coefficient being analyzed. If coef1 = -DBL_MAX, var1 is set to 0. alpar@1: * alpar@1: * value1 is value of the basic variable being analyzed in an adjacent alpar@1: * basis, which is defined as follows. Let the objective coefficient alpar@1: * reaches its minimal value (coef1) and continues decreasing. Then the alpar@1: * reduced cost of the limiting non-basic variable (var1) becomes dual alpar@1: * infeasible and the current basis becomes non-optimal that forces the alpar@1: * limiting non-basic variable to enter the basis replacing there some alpar@1: * basic variable that leaves the basis to keep primal feasibility. alpar@1: * Should note that on determining the adjacent basis current bounds alpar@1: * of the basic variable being analyzed are ignored as if it were free alpar@1: * (unbounded) variable, so it cannot leave the basis. It may happen alpar@1: * that no dual feasible adjacent basis exists, in which case value1 is alpar@1: * set to -DBL_MAX or +DBL_MAX. alpar@1: * alpar@1: * coef2 is the maximal value of the objective coefficient, at which alpar@1: * the basis still remains dual feasible and thus optimal. +DBL_MAX alpar@1: * means that the objective coefficient has no upper limit. alpar@1: * alpar@1: * var2 is the ordinal number of an auxiliary (1 to m) or structural alpar@1: * (m+1 to n) non-basic variable, whose reduced cost reaches its zero alpar@1: * bound first and thereby limits further increasing the objective alpar@1: * coefficient being analyzed. If coef2 = +DBL_MAX, var2 is set to 0. alpar@1: * alpar@1: * value2 is value of the basic variable being analyzed in an adjacent alpar@1: * basis, which is defined exactly in the same way as value1 above with alpar@1: * exception that now the objective coefficient is increasing. */ alpar@1: alpar@1: void glp_analyze_coef(glp_prob *P, int k, double *coef1, int *var1, alpar@1: double *value1, double *coef2, int *var2, double *value2) alpar@1: { GLPROW *row; GLPCOL *col; alpar@1: int m, n, type, stat, kase, p, q, dir, clen, cpiv, rlen, rpiv, alpar@1: *cind, *rind; alpar@1: double lb, ub, coef, x, lim_coef, new_x, d, delta, ll, uu, xx, alpar@1: *rval, *cval; alpar@1: /* sanity checks */ alpar@1: if (P == NULL || P->magic != GLP_PROB_MAGIC) alpar@1: xerror("glp_analyze_coef: P = %p; invalid problem object\n", alpar@1: P); alpar@1: m = P->m, n = P->n; alpar@1: if (!(P->pbs_stat == GLP_FEAS && P->dbs_stat == GLP_FEAS)) alpar@1: xerror("glp_analyze_coef: optimal basic solution required\n"); alpar@1: if (!(m == 0 || P->valid)) alpar@1: xerror("glp_analyze_coef: basis factorization required\n"); alpar@1: if (!(1 <= k && k <= m+n)) alpar@1: xerror("glp_analyze_coef: k = %d; variable number out of range" alpar@1: "\n", k); alpar@1: /* retrieve information about the specified basic variable x[k] alpar@1: whose objective coefficient c[k] is to be analyzed */ alpar@1: if (k <= m) alpar@1: { row = P->row[k]; alpar@1: type = row->type; alpar@1: lb = row->lb; alpar@1: ub = row->ub; alpar@1: coef = 0.0; alpar@1: stat = row->stat; alpar@1: x = row->prim; alpar@1: } alpar@1: else alpar@1: { col = P->col[k-m]; alpar@1: type = col->type; alpar@1: lb = col->lb; alpar@1: ub = col->ub; alpar@1: coef = col->coef; alpar@1: stat = col->stat; alpar@1: x = col->prim; alpar@1: } alpar@1: if (stat != GLP_BS) alpar@1: xerror("glp_analyze_coef: k = %d; non-basic variable not allow" alpar@1: "ed\n", k); alpar@1: /* allocate working arrays */ alpar@1: cind = xcalloc(1+m, sizeof(int)); alpar@1: cval = xcalloc(1+m, sizeof(double)); alpar@1: rind = xcalloc(1+n, sizeof(int)); alpar@1: rval = xcalloc(1+n, sizeof(double)); alpar@1: /* compute row of the simplex table corresponding to the basic alpar@1: variable x[k] */ alpar@1: rlen = glp_eval_tab_row(P, k, rind, rval); alpar@1: xassert(0 <= rlen && rlen <= n); alpar@1: /* perform analysis */ alpar@1: for (kase = -1; kase <= +1; kase += 2) alpar@1: { /* kase < 0 means objective coefficient c[k] is decreasing; alpar@1: kase > 0 means objective coefficient c[k] is increasing */ alpar@1: /* note that decreasing c[k] is equivalent to increasing dual alpar@1: variable lambda[k] and vice versa; we need to correctly set alpar@1: the dir flag as required by the routine glp_dual_rtest */ alpar@1: if (P->dir == GLP_MIN) alpar@1: dir = - kase; alpar@1: else if (P->dir == GLP_MAX) alpar@1: dir = + kase; alpar@1: else alpar@1: xassert(P != P); alpar@1: /* use the dual ratio test to determine non-basic variable alpar@1: x[q] whose reduced cost d[q] reaches zero bound first */ alpar@1: rpiv = glp_dual_rtest(P, rlen, rind, rval, dir, 1e-9); alpar@1: if (rpiv == 0) alpar@1: { /* nothing limits changing c[k] */ alpar@1: lim_coef = (kase < 0 ? -DBL_MAX : +DBL_MAX); alpar@1: q = 0; alpar@1: /* x[k] keeps its current value */ alpar@1: new_x = x; alpar@1: goto store; alpar@1: } alpar@1: /* non-basic variable x[q] limits changing coefficient c[k]; alpar@1: determine its status and reduced cost d[k] in the current alpar@1: basis */ alpar@1: xassert(1 <= rpiv && rpiv <= rlen); alpar@1: q = rind[rpiv]; alpar@1: xassert(1 <= q && q <= m+n); alpar@1: if (q <= m) alpar@1: { row = P->row[q]; alpar@1: stat = row->stat; alpar@1: d = row->dual; alpar@1: } alpar@1: else alpar@1: { col = P->col[q-m]; alpar@1: stat = col->stat; alpar@1: d = col->dual; alpar@1: } alpar@1: /* note that delta d[q] = new d[q] - d[q] = - d[q], because alpar@1: new d[q] = 0; delta d[q] = alfa[k,q] * delta c[k], so alpar@1: delta c[k] = delta d[q] / alfa[k,q] = - d[q] / alfa[k,q] */ alpar@1: xassert(rval[rpiv] != 0.0); alpar@1: delta = - d / rval[rpiv]; alpar@1: /* compute new c[k] = c[k] + delta c[k], which is the limiting alpar@1: value of the objective coefficient c[k] */ alpar@1: lim_coef = coef + delta; alpar@1: /* let c[k] continue decreasing/increasing that makes d[q] alpar@1: dual infeasible and forces x[q] to enter the basis; alpar@1: to perform the primal ratio test we need to know in which alpar@1: direction x[q] changes on entering the basis; we determine alpar@1: that analyzing the sign of delta d[q] (see above), since alpar@1: d[q] may be close to zero having wrong sign */ alpar@1: /* let, for simplicity, the problem is minimization */ alpar@1: if (kase < 0 && rval[rpiv] > 0.0 || alpar@1: kase > 0 && rval[rpiv] < 0.0) alpar@1: { /* delta d[q] < 0, so d[q] being non-negative will become alpar@1: negative, so x[q] will increase */ alpar@1: dir = +1; alpar@1: } alpar@1: else alpar@1: { /* delta d[q] > 0, so d[q] being non-positive will become alpar@1: positive, so x[q] will decrease */ alpar@1: dir = -1; alpar@1: } alpar@1: /* if the problem is maximization, correct the direction */ alpar@1: if (P->dir == GLP_MAX) dir = - dir; alpar@1: /* check that we didn't make a silly mistake */ alpar@1: if (dir > 0) alpar@1: xassert(stat == GLP_NL || stat == GLP_NF); alpar@1: else alpar@1: xassert(stat == GLP_NU || stat == GLP_NF); alpar@1: /* compute column of the simplex table corresponding to the alpar@1: non-basic variable x[q] */ alpar@1: clen = glp_eval_tab_col(P, q, cind, cval); alpar@1: /* make x[k] temporarily free (unbounded) */ alpar@1: if (k <= m) alpar@1: { row = P->row[k]; alpar@1: row->type = GLP_FR; alpar@1: row->lb = row->ub = 0.0; alpar@1: } alpar@1: else alpar@1: { col = P->col[k-m]; alpar@1: col->type = GLP_FR; alpar@1: col->lb = col->ub = 0.0; alpar@1: } alpar@1: /* use the primal ratio test to determine some basic variable alpar@1: which leaves the basis */ alpar@1: cpiv = glp_prim_rtest(P, clen, cind, cval, dir, 1e-9); alpar@1: /* restore original bounds of the basic variable x[k] */ alpar@1: if (k <= m) alpar@1: { row = P->row[k]; alpar@1: row->type = type; alpar@1: row->lb = lb, row->ub = ub; alpar@1: } alpar@1: else alpar@1: { col = P->col[k-m]; alpar@1: col->type = type; alpar@1: col->lb = lb, col->ub = ub; alpar@1: } alpar@1: if (cpiv == 0) alpar@1: { /* non-basic variable x[q] can change unlimitedly */ alpar@1: if (dir < 0 && rval[rpiv] > 0.0 || alpar@1: dir > 0 && rval[rpiv] < 0.0) alpar@1: { /* delta x[k] = alfa[k,q] * delta x[q] < 0 */ alpar@1: new_x = -DBL_MAX; alpar@1: } alpar@1: else alpar@1: { /* delta x[k] = alfa[k,q] * delta x[q] > 0 */ alpar@1: new_x = +DBL_MAX; alpar@1: } alpar@1: goto store; alpar@1: } alpar@1: /* some basic variable x[p] limits changing non-basic variable alpar@1: x[q] in the adjacent basis */ alpar@1: xassert(1 <= cpiv && cpiv <= clen); alpar@1: p = cind[cpiv]; alpar@1: xassert(1 <= p && p <= m+n); alpar@1: xassert(p != k); alpar@1: if (p <= m) alpar@1: { row = P->row[p]; alpar@1: xassert(row->stat == GLP_BS); alpar@1: ll = glp_get_row_lb(P, row->i); alpar@1: uu = glp_get_row_ub(P, row->i); alpar@1: xx = row->prim; alpar@1: } alpar@1: else alpar@1: { col = P->col[p-m]; alpar@1: xassert(col->stat == GLP_BS); alpar@1: ll = glp_get_col_lb(P, col->j); alpar@1: uu = glp_get_col_ub(P, col->j); alpar@1: xx = col->prim; alpar@1: } alpar@1: /* determine delta x[p] = new x[p] - x[p] */ alpar@1: if (dir < 0 && cval[cpiv] > 0.0 || alpar@1: dir > 0 && cval[cpiv] < 0.0) alpar@1: { /* delta x[p] < 0, so x[p] goes toward its lower bound */ alpar@1: xassert(ll != -DBL_MAX); alpar@1: delta = ll - xx; alpar@1: } alpar@1: else alpar@1: { /* delta x[p] > 0, so x[p] goes toward its upper bound */ alpar@1: xassert(uu != +DBL_MAX); alpar@1: delta = uu - xx; alpar@1: } alpar@1: /* compute new x[k] = x[k] + alfa[k,q] * delta x[q], where alpar@1: delta x[q] = delta x[p] / alfa[p,q] */ alpar@1: xassert(cval[cpiv] != 0.0); alpar@1: new_x = x + (rval[rpiv] / cval[cpiv]) * delta; alpar@1: store: /* store analysis results */ alpar@1: if (kase < 0) alpar@1: { if (coef1 != NULL) *coef1 = lim_coef; alpar@1: if (var1 != NULL) *var1 = q; alpar@1: if (value1 != NULL) *value1 = new_x; alpar@1: } alpar@1: else alpar@1: { if (coef2 != NULL) *coef2 = lim_coef; alpar@1: if (var2 != NULL) *var2 = q; alpar@1: if (value2 != NULL) *value2 = new_x; alpar@1: } alpar@1: } alpar@1: /* free working arrays */ alpar@1: xfree(cind); alpar@1: xfree(cval); alpar@1: xfree(rind); alpar@1: xfree(rval); alpar@1: return; alpar@1: } alpar@1: alpar@1: /* eof */