alpar@1: /* glpspx02.c (dual simplex method) */ 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 "glpspx.h" alpar@1: alpar@1: #define GLP_DEBUG 1 alpar@1: alpar@1: #if 0 alpar@1: #define GLP_LONG_STEP 1 alpar@1: #endif alpar@1: alpar@1: struct csa alpar@1: { /* common storage area */ alpar@1: /*--------------------------------------------------------------*/ alpar@1: /* LP data */ alpar@1: int m; alpar@1: /* number of rows (auxiliary variables), m > 0 */ alpar@1: int n; alpar@1: /* number of columns (structural variables), n > 0 */ alpar@1: char *type; /* char type[1+m+n]; */ alpar@1: /* type[0] is not used; alpar@1: type[k], 1 <= k <= m+n, is the type of variable x[k]: alpar@1: GLP_FR - free variable alpar@1: GLP_LO - variable with lower bound alpar@1: GLP_UP - variable with upper bound alpar@1: GLP_DB - double-bounded variable alpar@1: GLP_FX - fixed variable */ alpar@1: double *lb; /* double lb[1+m+n]; */ alpar@1: /* lb[0] is not used; alpar@1: lb[k], 1 <= k <= m+n, is an lower bound of variable x[k]; alpar@1: if x[k] has no lower bound, lb[k] is zero */ alpar@1: double *ub; /* double ub[1+m+n]; */ alpar@1: /* ub[0] is not used; alpar@1: ub[k], 1 <= k <= m+n, is an upper bound of variable x[k]; alpar@1: if x[k] has no upper bound, ub[k] is zero; alpar@1: if x[k] is of fixed type, ub[k] is the same as lb[k] */ alpar@1: double *coef; /* double coef[1+m+n]; */ alpar@1: /* coef[0] is not used; alpar@1: coef[k], 1 <= k <= m+n, is an objective coefficient at alpar@1: variable x[k] */ alpar@1: /*--------------------------------------------------------------*/ alpar@1: /* original bounds of variables */ alpar@1: char *orig_type; /* char orig_type[1+m+n]; */ alpar@1: double *orig_lb; /* double orig_lb[1+m+n]; */ alpar@1: double *orig_ub; /* double orig_ub[1+m+n]; */ alpar@1: /*--------------------------------------------------------------*/ alpar@1: /* original objective function */ alpar@1: double *obj; /* double obj[1+n]; */ alpar@1: /* obj[0] is a constant term of the original objective function; alpar@1: obj[j], 1 <= j <= n, is an original objective coefficient at alpar@1: structural variable x[m+j] */ alpar@1: double zeta; alpar@1: /* factor used to scale original objective coefficients; its alpar@1: sign defines original optimization direction: zeta > 0 means alpar@1: minimization, zeta < 0 means maximization */ alpar@1: /*--------------------------------------------------------------*/ alpar@1: /* constraint matrix A; it has m rows and n columns and is stored alpar@1: by columns */ alpar@1: int *A_ptr; /* int A_ptr[1+n+1]; */ alpar@1: /* A_ptr[0] is not used; alpar@1: A_ptr[j], 1 <= j <= n, is starting position of j-th column in alpar@1: arrays A_ind and A_val; note that A_ptr[1] is always 1; alpar@1: A_ptr[n+1] indicates the position after the last element in alpar@1: arrays A_ind and A_val */ alpar@1: int *A_ind; /* int A_ind[A_ptr[n+1]]; */ alpar@1: /* row indices */ alpar@1: double *A_val; /* double A_val[A_ptr[n+1]]; */ alpar@1: /* non-zero element values */ alpar@1: #if 1 /* 06/IV-2009 */ alpar@1: /* constraint matrix A stored by rows */ alpar@1: int *AT_ptr; /* int AT_ptr[1+m+1]; alpar@1: /* AT_ptr[0] is not used; alpar@1: AT_ptr[i], 1 <= i <= m, is starting position of i-th row in alpar@1: arrays AT_ind and AT_val; note that AT_ptr[1] is always 1; alpar@1: AT_ptr[m+1] indicates the position after the last element in alpar@1: arrays AT_ind and AT_val */ alpar@1: int *AT_ind; /* int AT_ind[AT_ptr[m+1]]; */ alpar@1: /* column indices */ alpar@1: double *AT_val; /* double AT_val[AT_ptr[m+1]]; */ alpar@1: /* non-zero element values */ alpar@1: #endif alpar@1: /*--------------------------------------------------------------*/ alpar@1: /* basis header */ alpar@1: int *head; /* int head[1+m+n]; */ alpar@1: /* head[0] is not used; alpar@1: head[i], 1 <= i <= m, is the ordinal number of basic variable alpar@1: xB[i]; head[i] = k means that xB[i] = x[k] and i-th column of alpar@1: matrix B is k-th column of matrix (I|-A); alpar@1: head[m+j], 1 <= j <= n, is the ordinal number of non-basic alpar@1: variable xN[j]; head[m+j] = k means that xN[j] = x[k] and j-th alpar@1: column of matrix N is k-th column of matrix (I|-A) */ alpar@1: #if 1 /* 06/IV-2009 */ alpar@1: int *bind; /* int bind[1+m+n]; */ alpar@1: /* bind[0] is not used; alpar@1: bind[k], 1 <= k <= m+n, is the position of k-th column of the alpar@1: matrix (I|-A) in the matrix (B|N); that is, bind[k] = k' means alpar@1: that head[k'] = k */ alpar@1: #endif alpar@1: char *stat; /* char stat[1+n]; */ alpar@1: /* stat[0] is not used; alpar@1: stat[j], 1 <= j <= n, is the status of non-basic variable alpar@1: xN[j], which defines its active bound: alpar@1: GLP_NL - lower bound is active alpar@1: GLP_NU - upper bound is active alpar@1: GLP_NF - free variable alpar@1: GLP_NS - fixed variable */ alpar@1: /*--------------------------------------------------------------*/ alpar@1: /* matrix B is the basis matrix; it is composed from columns of alpar@1: the augmented constraint matrix (I|-A) corresponding to basic alpar@1: variables and stored in a factorized (invertable) form */ alpar@1: int valid; alpar@1: /* factorization is valid only if this flag is set */ alpar@1: BFD *bfd; /* BFD bfd[1:m,1:m]; */ alpar@1: /* factorized (invertable) form of the basis matrix */ alpar@1: #if 0 /* 06/IV-2009 */ alpar@1: /*--------------------------------------------------------------*/ alpar@1: /* matrix N is a matrix composed from columns of the augmented alpar@1: constraint matrix (I|-A) corresponding to non-basic variables alpar@1: except fixed ones; it is stored by rows and changes every time alpar@1: the basis changes */ alpar@1: int *N_ptr; /* int N_ptr[1+m+1]; */ alpar@1: /* N_ptr[0] is not used; alpar@1: N_ptr[i], 1 <= i <= m, is starting position of i-th row in alpar@1: arrays N_ind and N_val; note that N_ptr[1] is always 1; alpar@1: N_ptr[m+1] indicates the position after the last element in alpar@1: arrays N_ind and N_val */ alpar@1: int *N_len; /* int N_len[1+m]; */ alpar@1: /* N_len[0] is not used; alpar@1: N_len[i], 1 <= i <= m, is length of i-th row (0 to n) */ alpar@1: int *N_ind; /* int N_ind[N_ptr[m+1]]; */ alpar@1: /* column indices */ alpar@1: double *N_val; /* double N_val[N_ptr[m+1]]; */ alpar@1: /* non-zero element values */ alpar@1: #endif alpar@1: /*--------------------------------------------------------------*/ alpar@1: /* working parameters */ alpar@1: int phase; alpar@1: /* search phase: alpar@1: 0 - not determined yet alpar@1: 1 - search for dual feasible solution alpar@1: 2 - search for optimal solution */ alpar@1: glp_long tm_beg; alpar@1: /* time value at the beginning of the search */ alpar@1: int it_beg; alpar@1: /* simplex iteration count at the beginning of the search */ alpar@1: int it_cnt; alpar@1: /* simplex iteration count; it increases by one every time the alpar@1: basis changes */ alpar@1: int it_dpy; alpar@1: /* simplex iteration count at the most recent display output */ alpar@1: /*--------------------------------------------------------------*/ alpar@1: /* basic solution components */ alpar@1: double *bbar; /* double bbar[1+m]; */ alpar@1: /* bbar[0] is not used on phase I; on phase II it is the current alpar@1: value of the original objective function; alpar@1: bbar[i], 1 <= i <= m, is primal value of basic variable xB[i] alpar@1: (if xB[i] is free, its primal value is not updated) */ alpar@1: double *cbar; /* double cbar[1+n]; */ alpar@1: /* cbar[0] is not used; alpar@1: cbar[j], 1 <= j <= n, is reduced cost of non-basic variable alpar@1: xN[j] (if xN[j] is fixed, its reduced cost is not updated) */ alpar@1: /*--------------------------------------------------------------*/ alpar@1: /* the following pricing technique options may be used: alpar@1: GLP_PT_STD - standard ("textbook") pricing; alpar@1: GLP_PT_PSE - projected steepest edge; alpar@1: GLP_PT_DVX - Devex pricing (not implemented yet); alpar@1: in case of GLP_PT_STD the reference space is not used, and all alpar@1: steepest edge coefficients are set to 1 */ alpar@1: int refct; alpar@1: /* this count is set to an initial value when the reference space alpar@1: is defined and decreases by one every time the basis changes; alpar@1: once this count reaches zero, the reference space is redefined alpar@1: again */ alpar@1: char *refsp; /* char refsp[1+m+n]; */ alpar@1: /* refsp[0] is not used; alpar@1: refsp[k], 1 <= k <= m+n, is the flag which means that variable alpar@1: x[k] belongs to the current reference space */ alpar@1: double *gamma; /* double gamma[1+m]; */ alpar@1: /* gamma[0] is not used; alpar@1: gamma[i], 1 <= i <= n, is the steepest edge coefficient for alpar@1: basic variable xB[i]; if xB[i] is free, gamma[i] is not used alpar@1: and just set to 1 */ alpar@1: /*--------------------------------------------------------------*/ alpar@1: /* basic variable xB[p] chosen to leave the basis */ alpar@1: int p; alpar@1: /* index of the basic variable xB[p] chosen, 1 <= p <= m; alpar@1: if the set of eligible basic variables is empty (i.e. if the alpar@1: current basic solution is primal feasible within a tolerance) alpar@1: and thus no variable has been chosen, p is set to 0 */ alpar@1: double delta; alpar@1: /* change of xB[p] in the adjacent basis; alpar@1: delta > 0 means that xB[p] violates its lower bound and will alpar@1: increase to achieve it in the adjacent basis; alpar@1: delta < 0 means that xB[p] violates its upper bound and will alpar@1: decrease to achieve it in the adjacent basis */ alpar@1: /*--------------------------------------------------------------*/ alpar@1: /* pivot row of the simplex table corresponding to basic variable alpar@1: xB[p] chosen is the following vector: alpar@1: T' * e[p] = - N' * inv(B') * e[p] = - N' * rho, alpar@1: where B' is a matrix transposed to the current basis matrix, alpar@1: N' is a matrix, whose rows are columns of the matrix (I|-A) alpar@1: corresponding to non-basic non-fixed variables */ alpar@1: int trow_nnz; alpar@1: /* number of non-zero components, 0 <= nnz <= n */ alpar@1: int *trow_ind; /* int trow_ind[1+n]; */ alpar@1: /* trow_ind[0] is not used; alpar@1: trow_ind[t], 1 <= t <= nnz, is an index of non-zero component, alpar@1: i.e. trow_ind[t] = j means that trow_vec[j] != 0 */ alpar@1: double *trow_vec; /* int trow_vec[1+n]; */ alpar@1: /* trow_vec[0] is not used; alpar@1: trow_vec[j], 1 <= j <= n, is a numeric value of j-th component alpar@1: of the row */ alpar@1: double trow_max; alpar@1: /* infinity (maximum) norm of the row (max |trow_vec[j]|) */ alpar@1: int trow_num; alpar@1: /* number of significant non-zero components, which means that: alpar@1: |trow_vec[j]| >= eps for j in trow_ind[1,...,num], alpar@1: |tcol_vec[j]| < eps for j in trow_ind[num+1,...,nnz], alpar@1: where eps is a pivot tolerance */ alpar@1: /*--------------------------------------------------------------*/ alpar@1: #ifdef GLP_LONG_STEP /* 07/IV-2009 */ alpar@1: int nbps; alpar@1: /* number of breakpoints, 0 <= nbps <= n */ alpar@1: struct bkpt alpar@1: { int j; alpar@1: /* index of non-basic variable xN[j], 1 <= j <= n */ alpar@1: double t; alpar@1: /* value of dual ray parameter at breakpoint, t >= 0 */ alpar@1: double dz; alpar@1: /* dz = zeta(t = t[k]) - zeta(t = 0) */ alpar@1: } *bkpt; /* struct bkpt bkpt[1+n]; */ alpar@1: /* bkpt[0] is not used; alpar@1: bkpt[k], 1 <= k <= nbps, is k-th breakpoint of the dual alpar@1: objective */ alpar@1: #endif alpar@1: /*--------------------------------------------------------------*/ alpar@1: /* non-basic variable xN[q] chosen to enter the basis */ alpar@1: int q; alpar@1: /* index of the non-basic variable xN[q] chosen, 1 <= q <= n; alpar@1: if no variable has been chosen, q is set to 0 */ alpar@1: double new_dq; alpar@1: /* reduced cost of xN[q] in the adjacent basis (it is the change alpar@1: of lambdaB[p]) */ alpar@1: /*--------------------------------------------------------------*/ alpar@1: /* pivot column of the simplex table corresponding to non-basic alpar@1: variable xN[q] chosen is the following vector: alpar@1: T * e[q] = - inv(B) * N * e[q] = - inv(B) * N[q], alpar@1: where B is the current basis matrix, N[q] is a column of the alpar@1: matrix (I|-A) corresponding to xN[q] */ alpar@1: int tcol_nnz; alpar@1: /* number of non-zero components, 0 <= nnz <= m */ alpar@1: int *tcol_ind; /* int tcol_ind[1+m]; */ alpar@1: /* tcol_ind[0] is not used; alpar@1: tcol_ind[t], 1 <= t <= nnz, is an index of non-zero component, alpar@1: i.e. tcol_ind[t] = i means that tcol_vec[i] != 0 */ alpar@1: double *tcol_vec; /* double tcol_vec[1+m]; */ alpar@1: /* tcol_vec[0] is not used; alpar@1: tcol_vec[i], 1 <= i <= m, is a numeric value of i-th component alpar@1: of the column */ alpar@1: /*--------------------------------------------------------------*/ alpar@1: /* working arrays */ alpar@1: double *work1; /* double work1[1+m]; */ alpar@1: double *work2; /* double work2[1+m]; */ alpar@1: double *work3; /* double work3[1+m]; */ alpar@1: double *work4; /* double work4[1+m]; */ alpar@1: }; alpar@1: alpar@1: static const double kappa = 0.10; alpar@1: alpar@1: /*********************************************************************** alpar@1: * alloc_csa - allocate common storage area alpar@1: * alpar@1: * This routine allocates all arrays in the common storage area (CSA) alpar@1: * and returns a pointer to the CSA. */ alpar@1: alpar@1: static struct csa *alloc_csa(glp_prob *lp) alpar@1: { struct csa *csa; alpar@1: int m = lp->m; alpar@1: int n = lp->n; alpar@1: int nnz = lp->nnz; alpar@1: csa = xmalloc(sizeof(struct csa)); alpar@1: xassert(m > 0 && n > 0); alpar@1: csa->m = m; alpar@1: csa->n = n; alpar@1: csa->type = xcalloc(1+m+n, sizeof(char)); alpar@1: csa->lb = xcalloc(1+m+n, sizeof(double)); alpar@1: csa->ub = xcalloc(1+m+n, sizeof(double)); alpar@1: csa->coef = xcalloc(1+m+n, sizeof(double)); alpar@1: csa->orig_type = xcalloc(1+m+n, sizeof(char)); alpar@1: csa->orig_lb = xcalloc(1+m+n, sizeof(double)); alpar@1: csa->orig_ub = xcalloc(1+m+n, sizeof(double)); alpar@1: csa->obj = xcalloc(1+n, sizeof(double)); alpar@1: csa->A_ptr = xcalloc(1+n+1, sizeof(int)); alpar@1: csa->A_ind = xcalloc(1+nnz, sizeof(int)); alpar@1: csa->A_val = xcalloc(1+nnz, sizeof(double)); alpar@1: #if 1 /* 06/IV-2009 */ alpar@1: csa->AT_ptr = xcalloc(1+m+1, sizeof(int)); alpar@1: csa->AT_ind = xcalloc(1+nnz, sizeof(int)); alpar@1: csa->AT_val = xcalloc(1+nnz, sizeof(double)); alpar@1: #endif alpar@1: csa->head = xcalloc(1+m+n, sizeof(int)); alpar@1: #if 1 /* 06/IV-2009 */ alpar@1: csa->bind = xcalloc(1+m+n, sizeof(int)); alpar@1: #endif alpar@1: csa->stat = xcalloc(1+n, sizeof(char)); alpar@1: #if 0 /* 06/IV-2009 */ alpar@1: csa->N_ptr = xcalloc(1+m+1, sizeof(int)); alpar@1: csa->N_len = xcalloc(1+m, sizeof(int)); alpar@1: csa->N_ind = NULL; /* will be allocated later */ alpar@1: csa->N_val = NULL; /* will be allocated later */ alpar@1: #endif alpar@1: csa->bbar = xcalloc(1+m, sizeof(double)); alpar@1: csa->cbar = xcalloc(1+n, sizeof(double)); alpar@1: csa->refsp = xcalloc(1+m+n, sizeof(char)); alpar@1: csa->gamma = xcalloc(1+m, sizeof(double)); alpar@1: csa->trow_ind = xcalloc(1+n, sizeof(int)); alpar@1: csa->trow_vec = xcalloc(1+n, sizeof(double)); alpar@1: #ifdef GLP_LONG_STEP /* 07/IV-2009 */ alpar@1: csa->bkpt = xcalloc(1+n, sizeof(struct bkpt)); alpar@1: #endif alpar@1: csa->tcol_ind = xcalloc(1+m, sizeof(int)); alpar@1: csa->tcol_vec = xcalloc(1+m, sizeof(double)); alpar@1: csa->work1 = xcalloc(1+m, sizeof(double)); alpar@1: csa->work2 = xcalloc(1+m, sizeof(double)); alpar@1: csa->work3 = xcalloc(1+m, sizeof(double)); alpar@1: csa->work4 = xcalloc(1+m, sizeof(double)); alpar@1: return csa; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * init_csa - initialize common storage area alpar@1: * alpar@1: * This routine initializes all data structures in the common storage alpar@1: * area (CSA). */ alpar@1: alpar@1: static void init_csa(struct csa *csa, glp_prob *lp) alpar@1: { int m = csa->m; alpar@1: int n = csa->n; alpar@1: char *type = csa->type; alpar@1: double *lb = csa->lb; alpar@1: double *ub = csa->ub; alpar@1: double *coef = csa->coef; alpar@1: char *orig_type = csa->orig_type; alpar@1: double *orig_lb = csa->orig_lb; alpar@1: double *orig_ub = csa->orig_ub; alpar@1: double *obj = csa->obj; alpar@1: int *A_ptr = csa->A_ptr; alpar@1: int *A_ind = csa->A_ind; alpar@1: double *A_val = csa->A_val; alpar@1: #if 1 /* 06/IV-2009 */ alpar@1: int *AT_ptr = csa->AT_ptr; alpar@1: int *AT_ind = csa->AT_ind; alpar@1: double *AT_val = csa->AT_val; alpar@1: #endif alpar@1: int *head = csa->head; alpar@1: #if 1 /* 06/IV-2009 */ alpar@1: int *bind = csa->bind; alpar@1: #endif alpar@1: char *stat = csa->stat; alpar@1: char *refsp = csa->refsp; alpar@1: double *gamma = csa->gamma; alpar@1: int i, j, k, loc; alpar@1: double cmax; alpar@1: /* auxiliary variables */ alpar@1: for (i = 1; i <= m; i++) alpar@1: { GLPROW *row = lp->row[i]; alpar@1: type[i] = (char)row->type; alpar@1: lb[i] = row->lb * row->rii; alpar@1: ub[i] = row->ub * row->rii; alpar@1: coef[i] = 0.0; alpar@1: } alpar@1: /* structural variables */ alpar@1: for (j = 1; j <= n; j++) alpar@1: { GLPCOL *col = lp->col[j]; alpar@1: type[m+j] = (char)col->type; alpar@1: lb[m+j] = col->lb / col->sjj; alpar@1: ub[m+j] = col->ub / col->sjj; alpar@1: coef[m+j] = col->coef * col->sjj; alpar@1: } alpar@1: /* original bounds of variables */ alpar@1: memcpy(&orig_type[1], &type[1], (m+n) * sizeof(char)); alpar@1: memcpy(&orig_lb[1], &lb[1], (m+n) * sizeof(double)); alpar@1: memcpy(&orig_ub[1], &ub[1], (m+n) * sizeof(double)); alpar@1: /* original objective function */ alpar@1: obj[0] = lp->c0; alpar@1: memcpy(&obj[1], &coef[m+1], n * sizeof(double)); alpar@1: /* factor used to scale original objective coefficients */ alpar@1: cmax = 0.0; alpar@1: for (j = 1; j <= n; j++) alpar@1: if (cmax < fabs(obj[j])) cmax = fabs(obj[j]); alpar@1: if (cmax == 0.0) cmax = 1.0; alpar@1: switch (lp->dir) alpar@1: { case GLP_MIN: alpar@1: csa->zeta = + 1.0 / cmax; alpar@1: break; alpar@1: case GLP_MAX: alpar@1: csa->zeta = - 1.0 / cmax; alpar@1: break; alpar@1: default: alpar@1: xassert(lp != lp); alpar@1: } alpar@1: #if 1 alpar@1: if (fabs(csa->zeta) < 1.0) csa->zeta *= 1000.0; alpar@1: #endif alpar@1: /* scale working objective coefficients */ alpar@1: for (j = 1; j <= n; j++) coef[m+j] *= csa->zeta; alpar@1: /* matrix A (by columns) */ alpar@1: loc = 1; alpar@1: for (j = 1; j <= n; j++) alpar@1: { GLPAIJ *aij; alpar@1: A_ptr[j] = loc; alpar@1: for (aij = lp->col[j]->ptr; aij != NULL; aij = aij->c_next) alpar@1: { A_ind[loc] = aij->row->i; alpar@1: A_val[loc] = aij->row->rii * aij->val * aij->col->sjj; alpar@1: loc++; alpar@1: } alpar@1: } alpar@1: A_ptr[n+1] = loc; alpar@1: xassert(loc-1 == lp->nnz); alpar@1: #if 1 /* 06/IV-2009 */ alpar@1: /* matrix A (by rows) */ alpar@1: loc = 1; alpar@1: for (i = 1; i <= m; i++) alpar@1: { GLPAIJ *aij; alpar@1: AT_ptr[i] = loc; alpar@1: for (aij = lp->row[i]->ptr; aij != NULL; aij = aij->r_next) alpar@1: { AT_ind[loc] = aij->col->j; alpar@1: AT_val[loc] = aij->row->rii * aij->val * aij->col->sjj; alpar@1: loc++; alpar@1: } alpar@1: } alpar@1: AT_ptr[m+1] = loc; alpar@1: xassert(loc-1 == lp->nnz); alpar@1: #endif alpar@1: /* basis header */ alpar@1: xassert(lp->valid); alpar@1: memcpy(&head[1], &lp->head[1], m * sizeof(int)); alpar@1: k = 0; alpar@1: for (i = 1; i <= m; i++) alpar@1: { GLPROW *row = lp->row[i]; alpar@1: if (row->stat != GLP_BS) alpar@1: { k++; alpar@1: xassert(k <= n); alpar@1: head[m+k] = i; alpar@1: stat[k] = (char)row->stat; alpar@1: } alpar@1: } alpar@1: for (j = 1; j <= n; j++) alpar@1: { GLPCOL *col = lp->col[j]; alpar@1: if (col->stat != GLP_BS) alpar@1: { k++; alpar@1: xassert(k <= n); alpar@1: head[m+k] = m + j; alpar@1: stat[k] = (char)col->stat; alpar@1: } alpar@1: } alpar@1: xassert(k == n); alpar@1: #if 1 /* 06/IV-2009 */ alpar@1: for (k = 1; k <= m+n; k++) alpar@1: bind[head[k]] = k; alpar@1: #endif alpar@1: /* factorization of matrix B */ alpar@1: csa->valid = 1, lp->valid = 0; alpar@1: csa->bfd = lp->bfd, lp->bfd = NULL; alpar@1: #if 0 /* 06/IV-2009 */ alpar@1: /* matrix N (by rows) */ alpar@1: alloc_N(csa); alpar@1: build_N(csa); alpar@1: #endif alpar@1: /* working parameters */ alpar@1: csa->phase = 0; alpar@1: csa->tm_beg = xtime(); alpar@1: csa->it_beg = csa->it_cnt = lp->it_cnt; alpar@1: csa->it_dpy = -1; alpar@1: /* reference space and steepest edge coefficients */ alpar@1: csa->refct = 0; alpar@1: memset(&refsp[1], 0, (m+n) * sizeof(char)); alpar@1: for (i = 1; i <= m; i++) gamma[i] = 1.0; alpar@1: return; alpar@1: } alpar@1: alpar@1: #if 1 /* copied from primal */ alpar@1: /*********************************************************************** alpar@1: * invert_B - compute factorization of the basis matrix alpar@1: * alpar@1: * This routine computes factorization of the current basis matrix B. alpar@1: * alpar@1: * If the operation is successful, the routine returns zero, otherwise alpar@1: * non-zero. */ alpar@1: alpar@1: static int inv_col(void *info, int i, int ind[], double val[]) alpar@1: { /* this auxiliary routine returns row indices and numeric values alpar@1: of non-zero elements of i-th column of the basis matrix */ alpar@1: struct csa *csa = info; alpar@1: int m = csa->m; alpar@1: #ifdef GLP_DEBUG alpar@1: int n = csa->n; alpar@1: #endif alpar@1: int *A_ptr = csa->A_ptr; alpar@1: int *A_ind = csa->A_ind; alpar@1: double *A_val = csa->A_val; alpar@1: int *head = csa->head; alpar@1: int k, len, ptr, t; alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= i && i <= m); alpar@1: #endif alpar@1: k = head[i]; /* B[i] is k-th column of (I|-A) */ alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= k && k <= m+n); alpar@1: #endif alpar@1: if (k <= m) alpar@1: { /* B[i] is k-th column of submatrix I */ alpar@1: len = 1; alpar@1: ind[1] = k; alpar@1: val[1] = 1.0; alpar@1: } alpar@1: else alpar@1: { /* B[i] is (k-m)-th column of submatrix (-A) */ alpar@1: ptr = A_ptr[k-m]; alpar@1: len = A_ptr[k-m+1] - ptr; alpar@1: memcpy(&ind[1], &A_ind[ptr], len * sizeof(int)); alpar@1: memcpy(&val[1], &A_val[ptr], len * sizeof(double)); alpar@1: for (t = 1; t <= len; t++) val[t] = - val[t]; alpar@1: } alpar@1: return len; alpar@1: } alpar@1: alpar@1: static int invert_B(struct csa *csa) alpar@1: { int ret; alpar@1: ret = bfd_factorize(csa->bfd, csa->m, NULL, inv_col, csa); alpar@1: csa->valid = (ret == 0); alpar@1: return ret; alpar@1: } alpar@1: #endif alpar@1: alpar@1: #if 1 /* copied from primal */ alpar@1: /*********************************************************************** alpar@1: * update_B - update factorization of the basis matrix alpar@1: * alpar@1: * This routine replaces i-th column of the basis matrix B by k-th alpar@1: * column of the augmented constraint matrix (I|-A) and then updates alpar@1: * the factorization of B. alpar@1: * alpar@1: * If the factorization has been successfully updated, the routine alpar@1: * returns zero, otherwise non-zero. */ alpar@1: alpar@1: static int update_B(struct csa *csa, int i, int k) alpar@1: { int m = csa->m; alpar@1: #ifdef GLP_DEBUG alpar@1: int n = csa->n; alpar@1: #endif alpar@1: int ret; alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= i && i <= m); alpar@1: xassert(1 <= k && k <= m+n); alpar@1: #endif alpar@1: if (k <= m) alpar@1: { /* new i-th column of B is k-th column of I */ alpar@1: int ind[1+1]; alpar@1: double val[1+1]; alpar@1: ind[1] = k; alpar@1: val[1] = 1.0; alpar@1: xassert(csa->valid); alpar@1: ret = bfd_update_it(csa->bfd, i, 0, 1, ind, val); alpar@1: } alpar@1: else alpar@1: { /* new i-th column of B is (k-m)-th column of (-A) */ alpar@1: int *A_ptr = csa->A_ptr; alpar@1: int *A_ind = csa->A_ind; alpar@1: double *A_val = csa->A_val; alpar@1: double *val = csa->work1; alpar@1: int beg, end, ptr, len; alpar@1: beg = A_ptr[k-m]; alpar@1: end = A_ptr[k-m+1]; alpar@1: len = 0; alpar@1: for (ptr = beg; ptr < end; ptr++) alpar@1: val[++len] = - A_val[ptr]; alpar@1: xassert(csa->valid); alpar@1: ret = bfd_update_it(csa->bfd, i, 0, len, &A_ind[beg-1], val); alpar@1: } alpar@1: csa->valid = (ret == 0); alpar@1: return ret; alpar@1: } alpar@1: #endif alpar@1: alpar@1: #if 1 /* copied from primal */ alpar@1: /*********************************************************************** alpar@1: * error_ftran - compute residual vector r = h - B * x alpar@1: * alpar@1: * This routine computes the residual vector r = h - B * x, where B is alpar@1: * the current basis matrix, h is the vector of right-hand sides, x is alpar@1: * the solution vector. */ alpar@1: alpar@1: static void error_ftran(struct csa *csa, double h[], double x[], alpar@1: double r[]) alpar@1: { int m = csa->m; alpar@1: #ifdef GLP_DEBUG alpar@1: int n = csa->n; alpar@1: #endif alpar@1: int *A_ptr = csa->A_ptr; alpar@1: int *A_ind = csa->A_ind; alpar@1: double *A_val = csa->A_val; alpar@1: int *head = csa->head; alpar@1: int i, k, beg, end, ptr; alpar@1: double temp; alpar@1: /* compute the residual vector: alpar@1: r = h - B * x = h - B[1] * x[1] - ... - B[m] * x[m], alpar@1: where B[1], ..., B[m] are columns of matrix B */ alpar@1: memcpy(&r[1], &h[1], m * sizeof(double)); alpar@1: for (i = 1; i <= m; i++) alpar@1: { temp = x[i]; alpar@1: if (temp == 0.0) continue; alpar@1: k = head[i]; /* B[i] is k-th column of (I|-A) */ alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= k && k <= m+n); alpar@1: #endif alpar@1: if (k <= m) alpar@1: { /* B[i] is k-th column of submatrix I */ alpar@1: r[k] -= temp; alpar@1: } alpar@1: else alpar@1: { /* B[i] is (k-m)-th column of submatrix (-A) */ alpar@1: beg = A_ptr[k-m]; alpar@1: end = A_ptr[k-m+1]; alpar@1: for (ptr = beg; ptr < end; ptr++) alpar@1: r[A_ind[ptr]] += A_val[ptr] * temp; alpar@1: } alpar@1: } alpar@1: return; alpar@1: } alpar@1: #endif alpar@1: alpar@1: #if 1 /* copied from primal */ alpar@1: /*********************************************************************** alpar@1: * refine_ftran - refine solution of B * x = h alpar@1: * alpar@1: * This routine performs one iteration to refine the solution of alpar@1: * the system B * x = h, where B is the current basis matrix, h is the alpar@1: * vector of right-hand sides, x is the solution vector. */ alpar@1: alpar@1: static void refine_ftran(struct csa *csa, double h[], double x[]) alpar@1: { int m = csa->m; alpar@1: double *r = csa->work1; alpar@1: double *d = csa->work1; alpar@1: int i; alpar@1: /* compute the residual vector r = h - B * x */ alpar@1: error_ftran(csa, h, x, r); alpar@1: /* compute the correction vector d = inv(B) * r */ alpar@1: xassert(csa->valid); alpar@1: bfd_ftran(csa->bfd, d); alpar@1: /* refine the solution vector (new x) = (old x) + d */ alpar@1: for (i = 1; i <= m; i++) x[i] += d[i]; alpar@1: return; alpar@1: } alpar@1: #endif alpar@1: alpar@1: #if 1 /* copied from primal */ alpar@1: /*********************************************************************** alpar@1: * error_btran - compute residual vector r = h - B'* x alpar@1: * alpar@1: * This routine computes the residual vector r = h - B'* x, where B' alpar@1: * is a matrix transposed to the current basis matrix, h is the vector alpar@1: * of right-hand sides, x is the solution vector. */ alpar@1: alpar@1: static void error_btran(struct csa *csa, double h[], double x[], alpar@1: double r[]) alpar@1: { int m = csa->m; alpar@1: #ifdef GLP_DEBUG alpar@1: int n = csa->n; alpar@1: #endif alpar@1: int *A_ptr = csa->A_ptr; alpar@1: int *A_ind = csa->A_ind; alpar@1: double *A_val = csa->A_val; alpar@1: int *head = csa->head; alpar@1: int i, k, beg, end, ptr; alpar@1: double temp; alpar@1: /* compute the residual vector r = b - B'* x */ alpar@1: for (i = 1; i <= m; i++) alpar@1: { /* r[i] := b[i] - (i-th column of B)'* x */ alpar@1: k = head[i]; /* B[i] is k-th column of (I|-A) */ alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= k && k <= m+n); alpar@1: #endif alpar@1: temp = h[i]; alpar@1: if (k <= m) alpar@1: { /* B[i] is k-th column of submatrix I */ alpar@1: temp -= x[k]; alpar@1: } alpar@1: else alpar@1: { /* B[i] is (k-m)-th column of submatrix (-A) */ alpar@1: beg = A_ptr[k-m]; alpar@1: end = A_ptr[k-m+1]; alpar@1: for (ptr = beg; ptr < end; ptr++) alpar@1: temp += A_val[ptr] * x[A_ind[ptr]]; alpar@1: } alpar@1: r[i] = temp; alpar@1: } alpar@1: return; alpar@1: } alpar@1: #endif alpar@1: alpar@1: #if 1 /* copied from primal */ alpar@1: /*********************************************************************** alpar@1: * refine_btran - refine solution of B'* x = h alpar@1: * alpar@1: * This routine performs one iteration to refine the solution of the alpar@1: * system B'* x = h, where B' is a matrix transposed to the current alpar@1: * basis matrix, h is the vector of right-hand sides, x is the solution alpar@1: * vector. */ alpar@1: alpar@1: static void refine_btran(struct csa *csa, double h[], double x[]) alpar@1: { int m = csa->m; alpar@1: double *r = csa->work1; alpar@1: double *d = csa->work1; alpar@1: int i; alpar@1: /* compute the residual vector r = h - B'* x */ alpar@1: error_btran(csa, h, x, r); alpar@1: /* compute the correction vector d = inv(B') * r */ alpar@1: xassert(csa->valid); alpar@1: bfd_btran(csa->bfd, d); alpar@1: /* refine the solution vector (new x) = (old x) + d */ alpar@1: for (i = 1; i <= m; i++) x[i] += d[i]; alpar@1: return; alpar@1: } alpar@1: #endif alpar@1: alpar@1: #if 1 /* copied from primal */ alpar@1: /*********************************************************************** alpar@1: * get_xN - determine current value of non-basic variable xN[j] alpar@1: * alpar@1: * This routine returns the current value of non-basic variable xN[j], alpar@1: * which is a value of its active bound. */ alpar@1: alpar@1: static double get_xN(struct csa *csa, int j) alpar@1: { int m = csa->m; alpar@1: #ifdef GLP_DEBUG alpar@1: int n = csa->n; alpar@1: #endif alpar@1: double *lb = csa->lb; alpar@1: double *ub = csa->ub; alpar@1: int *head = csa->head; alpar@1: char *stat = csa->stat; alpar@1: int k; alpar@1: double xN; alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= j && j <= n); alpar@1: #endif alpar@1: k = head[m+j]; /* x[k] = xN[j] */ alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= k && k <= m+n); alpar@1: #endif alpar@1: switch (stat[j]) alpar@1: { case GLP_NL: alpar@1: /* x[k] is on its lower bound */ alpar@1: xN = lb[k]; break; alpar@1: case GLP_NU: alpar@1: /* x[k] is on its upper bound */ alpar@1: xN = ub[k]; break; alpar@1: case GLP_NF: alpar@1: /* x[k] is free non-basic variable */ alpar@1: xN = 0.0; break; alpar@1: case GLP_NS: alpar@1: /* x[k] is fixed non-basic variable */ alpar@1: xN = lb[k]; break; alpar@1: default: alpar@1: xassert(stat != stat); alpar@1: } alpar@1: return xN; alpar@1: } alpar@1: #endif alpar@1: alpar@1: #if 1 /* copied from primal */ alpar@1: /*********************************************************************** alpar@1: * eval_beta - compute primal values of basic variables alpar@1: * alpar@1: * This routine computes current primal values of all basic variables: alpar@1: * alpar@1: * beta = - inv(B) * N * xN, alpar@1: * alpar@1: * where B is the current basis matrix, N is a matrix built of columns alpar@1: * of matrix (I|-A) corresponding to non-basic variables, and xN is the alpar@1: * vector of current values of non-basic variables. */ alpar@1: alpar@1: static void eval_beta(struct csa *csa, double beta[]) alpar@1: { int m = csa->m; alpar@1: int n = csa->n; alpar@1: int *A_ptr = csa->A_ptr; alpar@1: int *A_ind = csa->A_ind; alpar@1: double *A_val = csa->A_val; alpar@1: int *head = csa->head; alpar@1: double *h = csa->work2; alpar@1: int i, j, k, beg, end, ptr; alpar@1: double xN; alpar@1: /* compute the right-hand side vector: alpar@1: h := - N * xN = - N[1] * xN[1] - ... - N[n] * xN[n], alpar@1: where N[1], ..., N[n] are columns of matrix N */ alpar@1: for (i = 1; i <= m; i++) alpar@1: h[i] = 0.0; alpar@1: for (j = 1; j <= n; j++) alpar@1: { k = head[m+j]; /* x[k] = xN[j] */ alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= k && k <= m+n); alpar@1: #endif alpar@1: /* determine current value of xN[j] */ alpar@1: xN = get_xN(csa, j); alpar@1: if (xN == 0.0) continue; alpar@1: if (k <= m) alpar@1: { /* N[j] is k-th column of submatrix I */ alpar@1: h[k] -= xN; alpar@1: } alpar@1: else alpar@1: { /* N[j] is (k-m)-th column of submatrix (-A) */ alpar@1: beg = A_ptr[k-m]; alpar@1: end = A_ptr[k-m+1]; alpar@1: for (ptr = beg; ptr < end; ptr++) alpar@1: h[A_ind[ptr]] += xN * A_val[ptr]; alpar@1: } alpar@1: } alpar@1: /* solve system B * beta = h */ alpar@1: memcpy(&beta[1], &h[1], m * sizeof(double)); alpar@1: xassert(csa->valid); alpar@1: bfd_ftran(csa->bfd, beta); alpar@1: /* and refine the solution */ alpar@1: refine_ftran(csa, h, beta); alpar@1: return; alpar@1: } alpar@1: #endif alpar@1: alpar@1: #if 1 /* copied from primal */ alpar@1: /*********************************************************************** alpar@1: * eval_pi - compute vector of simplex multipliers alpar@1: * alpar@1: * This routine computes the vector of current simplex multipliers: alpar@1: * alpar@1: * pi = inv(B') * cB, alpar@1: * alpar@1: * where B' is a matrix transposed to the current basis matrix, cB is alpar@1: * a subvector of objective coefficients at basic variables. */ alpar@1: alpar@1: static void eval_pi(struct csa *csa, double pi[]) alpar@1: { int m = csa->m; alpar@1: double *c = csa->coef; alpar@1: int *head = csa->head; alpar@1: double *cB = csa->work2; alpar@1: int i; alpar@1: /* construct the right-hand side vector cB */ alpar@1: for (i = 1; i <= m; i++) alpar@1: cB[i] = c[head[i]]; alpar@1: /* solve system B'* pi = cB */ alpar@1: memcpy(&pi[1], &cB[1], m * sizeof(double)); alpar@1: xassert(csa->valid); alpar@1: bfd_btran(csa->bfd, pi); alpar@1: /* and refine the solution */ alpar@1: refine_btran(csa, cB, pi); alpar@1: return; alpar@1: } alpar@1: #endif alpar@1: alpar@1: #if 1 /* copied from primal */ alpar@1: /*********************************************************************** alpar@1: * eval_cost - compute reduced cost of non-basic variable xN[j] alpar@1: * alpar@1: * This routine computes the current reduced cost of non-basic variable alpar@1: * xN[j]: alpar@1: * alpar@1: * d[j] = cN[j] - N'[j] * pi, alpar@1: * alpar@1: * where cN[j] is the objective coefficient at variable xN[j], N[j] is alpar@1: * a column of the augmented constraint matrix (I|-A) corresponding to alpar@1: * xN[j], pi is the vector of simplex multipliers. */ alpar@1: alpar@1: static double eval_cost(struct csa *csa, double pi[], int j) alpar@1: { int m = csa->m; alpar@1: #ifdef GLP_DEBUG alpar@1: int n = csa->n; alpar@1: #endif alpar@1: double *coef = csa->coef; alpar@1: int *head = csa->head; alpar@1: int k; alpar@1: double dj; alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= j && j <= n); alpar@1: #endif alpar@1: k = head[m+j]; /* x[k] = xN[j] */ alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= k && k <= m+n); alpar@1: #endif alpar@1: dj = coef[k]; alpar@1: if (k <= m) alpar@1: { /* N[j] is k-th column of submatrix I */ alpar@1: dj -= pi[k]; alpar@1: } alpar@1: else alpar@1: { /* N[j] is (k-m)-th column of submatrix (-A) */ alpar@1: int *A_ptr = csa->A_ptr; alpar@1: int *A_ind = csa->A_ind; alpar@1: double *A_val = csa->A_val; alpar@1: int beg, end, ptr; alpar@1: beg = A_ptr[k-m]; alpar@1: end = A_ptr[k-m+1]; alpar@1: for (ptr = beg; ptr < end; ptr++) alpar@1: dj += A_val[ptr] * pi[A_ind[ptr]]; alpar@1: } alpar@1: return dj; alpar@1: } alpar@1: #endif alpar@1: alpar@1: #if 1 /* copied from primal */ alpar@1: /*********************************************************************** alpar@1: * eval_bbar - compute and store primal values of basic variables alpar@1: * alpar@1: * This routine computes primal values of all basic variables and then alpar@1: * stores them in the solution array. */ alpar@1: alpar@1: static void eval_bbar(struct csa *csa) alpar@1: { eval_beta(csa, csa->bbar); alpar@1: return; alpar@1: } alpar@1: #endif alpar@1: alpar@1: #if 1 /* copied from primal */ alpar@1: /*********************************************************************** alpar@1: * eval_cbar - compute and store reduced costs of non-basic variables alpar@1: * alpar@1: * This routine computes reduced costs of all non-basic variables and alpar@1: * then stores them in the solution array. */ alpar@1: alpar@1: static void eval_cbar(struct csa *csa) alpar@1: { alpar@1: #ifdef GLP_DEBUG alpar@1: int m = csa->m; alpar@1: #endif alpar@1: int n = csa->n; alpar@1: #ifdef GLP_DEBUG alpar@1: int *head = csa->head; alpar@1: #endif alpar@1: double *cbar = csa->cbar; alpar@1: double *pi = csa->work3; alpar@1: int j; alpar@1: #ifdef GLP_DEBUG alpar@1: int k; alpar@1: #endif alpar@1: /* compute simplex multipliers */ alpar@1: eval_pi(csa, pi); alpar@1: /* compute and store reduced costs */ alpar@1: for (j = 1; j <= n; j++) alpar@1: { alpar@1: #ifdef GLP_DEBUG alpar@1: k = head[m+j]; /* x[k] = xN[j] */ alpar@1: xassert(1 <= k && k <= m+n); alpar@1: #endif alpar@1: cbar[j] = eval_cost(csa, pi, j); alpar@1: } alpar@1: return; alpar@1: } alpar@1: #endif alpar@1: alpar@1: /*********************************************************************** alpar@1: * reset_refsp - reset the reference space alpar@1: * alpar@1: * This routine resets (redefines) the reference space used in the alpar@1: * projected steepest edge pricing algorithm. */ alpar@1: alpar@1: static void reset_refsp(struct csa *csa) alpar@1: { int m = csa->m; alpar@1: int n = csa->n; alpar@1: int *head = csa->head; alpar@1: char *refsp = csa->refsp; alpar@1: double *gamma = csa->gamma; alpar@1: int i, k; alpar@1: xassert(csa->refct == 0); alpar@1: csa->refct = 1000; alpar@1: memset(&refsp[1], 0, (m+n) * sizeof(char)); alpar@1: for (i = 1; i <= m; i++) alpar@1: { k = head[i]; /* x[k] = xB[i] */ alpar@1: refsp[k] = 1; alpar@1: gamma[i] = 1.0; alpar@1: } alpar@1: return; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * eval_gamma - compute steepest edge coefficients alpar@1: * alpar@1: * This routine computes the vector of steepest edge coefficients for alpar@1: * all basic variables (except free ones) using its direct definition: alpar@1: * alpar@1: * gamma[i] = eta[i] + sum alfa[i,j]^2, i = 1,...,m, alpar@1: * j in C alpar@1: * alpar@1: * where eta[i] = 1 means that xB[i] is in the current reference space, alpar@1: * and 0 otherwise; C is a set of non-basic non-fixed variables xN[j], alpar@1: * which are in the current reference space; alfa[i,j] are elements of alpar@1: * the current simplex table. alpar@1: * alpar@1: * NOTE: The routine is intended only for debugginig purposes. */ alpar@1: alpar@1: static void eval_gamma(struct csa *csa, double gamma[]) alpar@1: { int m = csa->m; alpar@1: int n = csa->n; alpar@1: char *type = csa->type; alpar@1: int *head = csa->head; alpar@1: char *refsp = csa->refsp; alpar@1: double *alfa = csa->work3; alpar@1: double *h = csa->work3; alpar@1: int i, j, k; alpar@1: /* gamma[i] := eta[i] (or 1, if xB[i] is free) */ alpar@1: for (i = 1; i <= m; i++) alpar@1: { k = head[i]; /* x[k] = xB[i] */ alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= k && k <= m+n); alpar@1: #endif alpar@1: if (type[k] == GLP_FR) alpar@1: gamma[i] = 1.0; alpar@1: else alpar@1: gamma[i] = (refsp[k] ? 1.0 : 0.0); alpar@1: } alpar@1: /* compute columns of the current simplex table */ alpar@1: for (j = 1; j <= n; j++) alpar@1: { k = head[m+j]; /* x[k] = xN[j] */ alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= k && k <= m+n); alpar@1: #endif alpar@1: /* skip column, if xN[j] is not in C */ alpar@1: if (!refsp[k]) continue; alpar@1: #ifdef GLP_DEBUG alpar@1: /* set C must not contain fixed variables */ alpar@1: xassert(type[k] != GLP_FX); alpar@1: #endif alpar@1: /* construct the right-hand side vector h = - N[j] */ alpar@1: for (i = 1; i <= m; i++) alpar@1: h[i] = 0.0; alpar@1: if (k <= m) alpar@1: { /* N[j] is k-th column of submatrix I */ alpar@1: h[k] = -1.0; alpar@1: } alpar@1: else alpar@1: { /* N[j] is (k-m)-th column of submatrix (-A) */ alpar@1: int *A_ptr = csa->A_ptr; alpar@1: int *A_ind = csa->A_ind; alpar@1: double *A_val = csa->A_val; alpar@1: int beg, end, ptr; alpar@1: beg = A_ptr[k-m]; alpar@1: end = A_ptr[k-m+1]; alpar@1: for (ptr = beg; ptr < end; ptr++) alpar@1: h[A_ind[ptr]] = A_val[ptr]; alpar@1: } alpar@1: /* solve system B * alfa = h */ alpar@1: xassert(csa->valid); alpar@1: bfd_ftran(csa->bfd, alfa); alpar@1: /* gamma[i] := gamma[i] + alfa[i,j]^2 */ alpar@1: for (i = 1; i <= m; i++) alpar@1: { k = head[i]; /* x[k] = xB[i] */ alpar@1: if (type[k] != GLP_FR) alpar@1: gamma[i] += alfa[i] * alfa[i]; alpar@1: } alpar@1: } alpar@1: return; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * chuzr - choose basic variable (row of the simplex table) alpar@1: * alpar@1: * This routine chooses basic variable xB[p] having largest weighted alpar@1: * bound violation: alpar@1: * alpar@1: * |r[p]| / sqrt(gamma[p]) = max |r[i]| / sqrt(gamma[i]), alpar@1: * i in I alpar@1: * alpar@1: * / lB[i] - beta[i], if beta[i] < lB[i] alpar@1: * | alpar@1: * r[i] = < 0, if lB[i] <= beta[i] <= uB[i] alpar@1: * | alpar@1: * \ uB[i] - beta[i], if beta[i] > uB[i] alpar@1: * alpar@1: * where beta[i] is primal value of xB[i] in the current basis, lB[i] alpar@1: * and uB[i] are lower and upper bounds of xB[i], I is a subset of alpar@1: * eligible basic variables, which significantly violates their bounds, alpar@1: * gamma[i] is the steepest edge coefficient. alpar@1: * alpar@1: * If |r[i]| is less than a specified tolerance, xB[i] is not included alpar@1: * in I and therefore ignored. alpar@1: * alpar@1: * If I is empty and no variable has been chosen, p is set to 0. */ alpar@1: alpar@1: static void chuzr(struct csa *csa, double tol_bnd) alpar@1: { int m = csa->m; alpar@1: #ifdef GLP_DEBUG alpar@1: int n = csa->n; alpar@1: #endif alpar@1: char *type = csa->type; alpar@1: double *lb = csa->lb; alpar@1: double *ub = csa->ub; alpar@1: int *head = csa->head; alpar@1: double *bbar = csa->bbar; alpar@1: double *gamma = csa->gamma; alpar@1: int i, k, p; alpar@1: double delta, best, eps, ri, temp; alpar@1: /* nothing is chosen so far */ alpar@1: p = 0, delta = 0.0, best = 0.0; alpar@1: /* look through the list of basic variables */ alpar@1: for (i = 1; i <= m; i++) alpar@1: { k = head[i]; /* x[k] = xB[i] */ alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= k && k <= m+n); alpar@1: #endif alpar@1: /* determine bound violation ri[i] */ alpar@1: ri = 0.0; alpar@1: if (type[k] == GLP_LO || type[k] == GLP_DB || alpar@1: type[k] == GLP_FX) alpar@1: { /* xB[i] has lower bound */ alpar@1: eps = tol_bnd * (1.0 + kappa * fabs(lb[k])); alpar@1: if (bbar[i] < lb[k] - eps) alpar@1: { /* and significantly violates it */ alpar@1: ri = lb[k] - bbar[i]; alpar@1: } alpar@1: } alpar@1: if (type[k] == GLP_UP || type[k] == GLP_DB || alpar@1: type[k] == GLP_FX) alpar@1: { /* xB[i] has upper bound */ alpar@1: eps = tol_bnd * (1.0 + kappa * fabs(ub[k])); alpar@1: if (bbar[i] > ub[k] + eps) alpar@1: { /* and significantly violates it */ alpar@1: ri = ub[k] - bbar[i]; alpar@1: } alpar@1: } alpar@1: /* if xB[i] is not eligible, skip it */ alpar@1: if (ri == 0.0) continue; alpar@1: /* xB[i] is eligible basic variable; choose one with largest alpar@1: weighted bound violation */ alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(gamma[i] >= 0.0); alpar@1: #endif alpar@1: temp = gamma[i]; alpar@1: if (temp < DBL_EPSILON) temp = DBL_EPSILON; alpar@1: temp = (ri * ri) / temp; alpar@1: if (best < temp) alpar@1: p = i, delta = ri, best = temp; alpar@1: } alpar@1: /* store the index of basic variable xB[p] chosen and its change alpar@1: in the adjacent basis */ alpar@1: csa->p = p; alpar@1: csa->delta = delta; alpar@1: return; alpar@1: } alpar@1: alpar@1: #if 1 /* copied from primal */ alpar@1: /*********************************************************************** alpar@1: * eval_rho - compute pivot row of the inverse alpar@1: * alpar@1: * This routine computes the pivot (p-th) row of the inverse inv(B), alpar@1: * which corresponds to basic variable xB[p] chosen: alpar@1: * alpar@1: * rho = inv(B') * e[p], alpar@1: * alpar@1: * where B' is a matrix transposed to the current basis matrix, e[p] alpar@1: * is unity vector. */ alpar@1: alpar@1: static void eval_rho(struct csa *csa, double rho[]) alpar@1: { int m = csa->m; alpar@1: int p = csa->p; alpar@1: double *e = rho; alpar@1: int i; alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= p && p <= m); alpar@1: #endif alpar@1: /* construct the right-hand side vector e[p] */ alpar@1: for (i = 1; i <= m; i++) alpar@1: e[i] = 0.0; alpar@1: e[p] = 1.0; alpar@1: /* solve system B'* rho = e[p] */ alpar@1: xassert(csa->valid); alpar@1: bfd_btran(csa->bfd, rho); alpar@1: return; alpar@1: } alpar@1: #endif alpar@1: alpar@1: #if 1 /* copied from primal */ alpar@1: /*********************************************************************** alpar@1: * refine_rho - refine pivot row of the inverse alpar@1: * alpar@1: * This routine refines the pivot row of the inverse inv(B) assuming alpar@1: * that it was previously computed by the routine eval_rho. */ alpar@1: alpar@1: static void refine_rho(struct csa *csa, double rho[]) alpar@1: { int m = csa->m; alpar@1: int p = csa->p; alpar@1: double *e = csa->work3; alpar@1: int i; alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= p && p <= m); alpar@1: #endif alpar@1: /* construct the right-hand side vector e[p] */ alpar@1: for (i = 1; i <= m; i++) alpar@1: e[i] = 0.0; alpar@1: e[p] = 1.0; alpar@1: /* refine solution of B'* rho = e[p] */ alpar@1: refine_btran(csa, e, rho); alpar@1: return; alpar@1: } alpar@1: #endif alpar@1: alpar@1: #if 1 /* 06/IV-2009 */ alpar@1: /*********************************************************************** alpar@1: * eval_trow - compute pivot row of the simplex table alpar@1: * alpar@1: * This routine computes the pivot row of the simplex table, which alpar@1: * corresponds to basic variable xB[p] chosen. alpar@1: * alpar@1: * The pivot row is the following vector: alpar@1: * alpar@1: * trow = T'* e[p] = - N'* inv(B') * e[p] = - N' * rho, alpar@1: * alpar@1: * where rho is the pivot row of the inverse inv(B) previously computed alpar@1: * by the routine eval_rho. alpar@1: * alpar@1: * Note that elements of the pivot row corresponding to fixed non-basic alpar@1: * variables are not computed. alpar@1: * alpar@1: * NOTES alpar@1: * alpar@1: * Computing pivot row of the simplex table is one of the most time alpar@1: * consuming operations, and for some instances it may take more than alpar@1: * 50% of the total solution time. alpar@1: * alpar@1: * In the current implementation there are two routines to compute the alpar@1: * pivot row. The routine eval_trow1 computes elements of the pivot row alpar@1: * as inner products of columns of the matrix N and the vector rho; it alpar@1: * is used when the vector rho is relatively dense. The routine alpar@1: * eval_trow2 computes the pivot row as a linear combination of rows of alpar@1: * the matrix N; it is used when the vector rho is relatively sparse. */ alpar@1: alpar@1: static void eval_trow1(struct csa *csa, double rho[]) alpar@1: { int m = csa->m; alpar@1: int n = csa->n; alpar@1: int *A_ptr = csa->A_ptr; alpar@1: int *A_ind = csa->A_ind; alpar@1: double *A_val = csa->A_val; alpar@1: int *head = csa->head; alpar@1: char *stat = csa->stat; alpar@1: int *trow_ind = csa->trow_ind; alpar@1: double *trow_vec = csa->trow_vec; alpar@1: int j, k, beg, end, ptr, nnz; alpar@1: double temp; alpar@1: /* compute the pivot row as inner products of columns of the alpar@1: matrix N and vector rho: trow[j] = - rho * N[j] */ alpar@1: nnz = 0; alpar@1: for (j = 1; j <= n; j++) alpar@1: { if (stat[j] == GLP_NS) alpar@1: { /* xN[j] is fixed */ alpar@1: trow_vec[j] = 0.0; alpar@1: continue; alpar@1: } alpar@1: k = head[m+j]; /* x[k] = xN[j] */ alpar@1: if (k <= m) alpar@1: { /* N[j] is k-th column of submatrix I */ alpar@1: temp = - rho[k]; alpar@1: } alpar@1: else alpar@1: { /* N[j] is (k-m)-th column of submatrix (-A) */ alpar@1: beg = A_ptr[k-m], end = A_ptr[k-m+1]; alpar@1: temp = 0.0; alpar@1: for (ptr = beg; ptr < end; ptr++) alpar@1: temp += rho[A_ind[ptr]] * A_val[ptr]; alpar@1: } alpar@1: if (temp != 0.0) alpar@1: trow_ind[++nnz] = j; alpar@1: trow_vec[j] = temp; alpar@1: } alpar@1: csa->trow_nnz = nnz; alpar@1: return; alpar@1: } alpar@1: alpar@1: static void eval_trow2(struct csa *csa, double rho[]) alpar@1: { int m = csa->m; alpar@1: int n = csa->n; alpar@1: int *AT_ptr = csa->AT_ptr; alpar@1: int *AT_ind = csa->AT_ind; alpar@1: double *AT_val = csa->AT_val; alpar@1: int *bind = csa->bind; alpar@1: char *stat = csa->stat; alpar@1: int *trow_ind = csa->trow_ind; alpar@1: double *trow_vec = csa->trow_vec; alpar@1: int i, j, beg, end, ptr, nnz; alpar@1: double temp; alpar@1: /* clear the pivot row */ alpar@1: for (j = 1; j <= n; j++) alpar@1: trow_vec[j] = 0.0; alpar@1: /* compute the pivot row as a linear combination of rows of the alpar@1: matrix N: trow = - rho[1] * N'[1] - ... - rho[m] * N'[m] */ alpar@1: for (i = 1; i <= m; i++) alpar@1: { temp = rho[i]; alpar@1: if (temp == 0.0) continue; alpar@1: /* trow := trow - rho[i] * N'[i] */ alpar@1: j = bind[i] - m; /* x[i] = xN[j] */ alpar@1: if (j >= 1 && stat[j] != GLP_NS) alpar@1: trow_vec[j] -= temp; alpar@1: beg = AT_ptr[i], end = AT_ptr[i+1]; alpar@1: for (ptr = beg; ptr < end; ptr++) alpar@1: { j = bind[m + AT_ind[ptr]] - m; /* x[k] = xN[j] */ alpar@1: if (j >= 1 && stat[j] != GLP_NS) alpar@1: trow_vec[j] += temp * AT_val[ptr]; alpar@1: } alpar@1: } alpar@1: /* construct sparse pattern of the pivot row */ alpar@1: nnz = 0; alpar@1: for (j = 1; j <= n; j++) alpar@1: { if (trow_vec[j] != 0.0) alpar@1: trow_ind[++nnz] = j; alpar@1: } alpar@1: csa->trow_nnz = nnz; alpar@1: return; alpar@1: } alpar@1: alpar@1: static void eval_trow(struct csa *csa, double rho[]) alpar@1: { int m = csa->m; alpar@1: int i, nnz; alpar@1: double dens; alpar@1: /* determine the density of the vector rho */ alpar@1: nnz = 0; alpar@1: for (i = 1; i <= m; i++) alpar@1: if (rho[i] != 0.0) nnz++; alpar@1: dens = (double)nnz / (double)m; alpar@1: if (dens >= 0.20) alpar@1: { /* rho is relatively dense */ alpar@1: eval_trow1(csa, rho); alpar@1: } alpar@1: else alpar@1: { /* rho is relatively sparse */ alpar@1: eval_trow2(csa, rho); alpar@1: } alpar@1: return; alpar@1: } alpar@1: #endif alpar@1: alpar@1: /*********************************************************************** alpar@1: * sort_trow - sort pivot row of the simplex table alpar@1: * alpar@1: * This routine reorders the list of non-zero elements of the pivot alpar@1: * row to put significant elements, whose magnitude is not less than alpar@1: * a specified tolerance, in front of the list, and stores the number alpar@1: * of significant elements in trow_num. */ alpar@1: alpar@1: static void sort_trow(struct csa *csa, double tol_piv) alpar@1: { alpar@1: #ifdef GLP_DEBUG alpar@1: int n = csa->n; alpar@1: char *stat = csa->stat; alpar@1: #endif alpar@1: int nnz = csa->trow_nnz; alpar@1: int *trow_ind = csa->trow_ind; alpar@1: double *trow_vec = csa->trow_vec; alpar@1: int j, num, pos; alpar@1: double big, eps, temp; alpar@1: /* compute infinity (maximum) norm of the row */ alpar@1: big = 0.0; alpar@1: for (pos = 1; pos <= nnz; pos++) alpar@1: { alpar@1: #ifdef GLP_DEBUG alpar@1: j = trow_ind[pos]; alpar@1: xassert(1 <= j && j <= n); alpar@1: xassert(stat[j] != GLP_NS); alpar@1: #endif alpar@1: temp = fabs(trow_vec[trow_ind[pos]]); alpar@1: if (big < temp) big = temp; alpar@1: } alpar@1: csa->trow_max = big; alpar@1: /* determine absolute pivot tolerance */ alpar@1: eps = tol_piv * (1.0 + 0.01 * big); alpar@1: /* move significant row components to the front of the list */ alpar@1: for (num = 0; num < nnz; ) alpar@1: { j = trow_ind[nnz]; alpar@1: if (fabs(trow_vec[j]) < eps) alpar@1: nnz--; alpar@1: else alpar@1: { num++; alpar@1: trow_ind[nnz] = trow_ind[num]; alpar@1: trow_ind[num] = j; alpar@1: } alpar@1: } alpar@1: csa->trow_num = num; alpar@1: return; alpar@1: } alpar@1: alpar@1: #ifdef GLP_LONG_STEP /* 07/IV-2009 */ alpar@1: static int ls_func(const void *p1_, const void *p2_) alpar@1: { const struct bkpt *p1 = p1_, *p2 = p2_; alpar@1: if (p1->t < p2->t) return -1; alpar@1: if (p1->t > p2->t) return +1; alpar@1: return 0; alpar@1: } alpar@1: alpar@1: static int ls_func1(const void *p1_, const void *p2_) alpar@1: { const struct bkpt *p1 = p1_, *p2 = p2_; alpar@1: if (p1->dz < p2->dz) return -1; alpar@1: if (p1->dz > p2->dz) return +1; alpar@1: return 0; alpar@1: } alpar@1: alpar@1: static void long_step(struct csa *csa) alpar@1: { int m = csa->m; alpar@1: #ifdef GLP_DEBUG alpar@1: int n = csa->n; alpar@1: #endif alpar@1: char *type = csa->type; alpar@1: double *lb = csa->lb; alpar@1: double *ub = csa->ub; alpar@1: int *head = csa->head; alpar@1: char *stat = csa->stat; alpar@1: double *cbar = csa->cbar; alpar@1: double delta = csa->delta; alpar@1: int *trow_ind = csa->trow_ind; alpar@1: double *trow_vec = csa->trow_vec; alpar@1: int trow_num = csa->trow_num; alpar@1: struct bkpt *bkpt = csa->bkpt; alpar@1: int j, k, kk, nbps, pos; alpar@1: double alfa, s, slope, dzmax; alpar@1: /* delta > 0 means that xB[p] violates its lower bound, so to alpar@1: increase the dual objective lambdaB[p] must increase; alpar@1: delta < 0 means that xB[p] violates its upper bound, so to alpar@1: increase the dual objective lambdaB[p] must decrease */ alpar@1: /* s := sign(delta) */ alpar@1: s = (delta > 0.0 ? +1.0 : -1.0); alpar@1: /* determine breakpoints of the dual objective */ alpar@1: nbps = 0; alpar@1: for (pos = 1; pos <= trow_num; pos++) alpar@1: { j = trow_ind[pos]; alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= j && j <= n); alpar@1: xassert(stat[j] != GLP_NS); alpar@1: #endif alpar@1: /* if there is free non-basic variable, switch to the standard alpar@1: ratio test */ alpar@1: if (stat[j] == GLP_NF) alpar@1: { nbps = 0; alpar@1: goto done; alpar@1: } alpar@1: /* lambdaN[j] = ... - alfa * t - ..., where t = s * lambdaB[i] alpar@1: is the dual ray parameter, t >= 0 */ alpar@1: alfa = s * trow_vec[j]; alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(alfa != 0.0); alpar@1: xassert(stat[j] == GLP_NL || stat[j] == GLP_NU); alpar@1: #endif alpar@1: if (alfa > 0.0 && stat[j] == GLP_NL || alpar@1: alfa < 0.0 && stat[j] == GLP_NU) alpar@1: { /* either lambdaN[j] >= 0 (if stat = GLP_NL) and decreases alpar@1: or lambdaN[j] <= 0 (if stat = GLP_NU) and increases; in alpar@1: both cases we have a breakpoint */ alpar@1: nbps++; alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(nbps <= n); alpar@1: #endif alpar@1: bkpt[nbps].j = j; alpar@1: bkpt[nbps].t = cbar[j] / alfa; alpar@1: /* alpar@1: if (stat[j] == GLP_NL && cbar[j] < 0.0 || alpar@1: stat[j] == GLP_NU && cbar[j] > 0.0) alpar@1: xprintf("%d %g\n", stat[j], cbar[j]); alpar@1: */ alpar@1: /* if t is negative, replace it by exact zero (see comments alpar@1: in the routine chuzc) */ alpar@1: if (bkpt[nbps].t < 0.0) bkpt[nbps].t = 0.0; alpar@1: } alpar@1: } alpar@1: /* if there are less than two breakpoints, switch to the standard alpar@1: ratio test */ alpar@1: if (nbps < 2) alpar@1: { nbps = 0; alpar@1: goto done; alpar@1: } alpar@1: /* sort breakpoints by ascending the dual ray parameter, t */ alpar@1: qsort(&bkpt[1], nbps, sizeof(struct bkpt), ls_func); alpar@1: /* determine last breakpoint, at which the dual objective still alpar@1: greater than at t = 0 */ alpar@1: dzmax = 0.0; alpar@1: slope = fabs(delta); /* initial slope */ alpar@1: for (kk = 1; kk <= nbps; kk++) alpar@1: { if (kk == 1) alpar@1: bkpt[kk].dz = alpar@1: 0.0 + slope * (bkpt[kk].t - 0.0); alpar@1: else alpar@1: bkpt[kk].dz = alpar@1: bkpt[kk-1].dz + slope * (bkpt[kk].t - bkpt[kk-1].t); alpar@1: if (dzmax < bkpt[kk].dz) alpar@1: dzmax = bkpt[kk].dz; alpar@1: else if (bkpt[kk].dz < 0.05 * (1.0 + dzmax)) alpar@1: { nbps = kk - 1; alpar@1: break; alpar@1: } alpar@1: j = bkpt[kk].j; alpar@1: k = head[m+j]; /* x[k] = xN[j] */ alpar@1: if (type[k] == GLP_DB) alpar@1: slope -= fabs(trow_vec[j]) * (ub[k] - lb[k]); alpar@1: else alpar@1: { nbps = kk; alpar@1: break; alpar@1: } alpar@1: } alpar@1: /* if there are less than two breakpoints, switch to the standard alpar@1: ratio test */ alpar@1: if (nbps < 2) alpar@1: { nbps = 0; alpar@1: goto done; alpar@1: } alpar@1: /* sort breakpoints by ascending the dual change, dz */ alpar@1: qsort(&bkpt[1], nbps, sizeof(struct bkpt), ls_func1); alpar@1: /* alpar@1: for (kk = 1; kk <= nbps; kk++) alpar@1: xprintf("%d; t = %g; dz = %g\n", kk, bkpt[kk].t, bkpt[kk].dz); alpar@1: */ alpar@1: done: csa->nbps = nbps; alpar@1: return; alpar@1: } alpar@1: #endif alpar@1: alpar@1: /*********************************************************************** alpar@1: * chuzc - choose non-basic variable (column of the simplex table) alpar@1: * alpar@1: * This routine chooses non-basic variable xN[q], which being entered alpar@1: * in the basis keeps dual feasibility of the basic solution. alpar@1: * alpar@1: * The parameter rtol is a relative tolerance used to relax zero bounds alpar@1: * of reduced costs of non-basic variables. If rtol = 0, the routine alpar@1: * implements the standard ratio test. Otherwise, if rtol > 0, the alpar@1: * routine implements Harris' two-pass ratio test. In the latter case alpar@1: * rtol should be about three times less than a tolerance used to check alpar@1: * dual feasibility. */ alpar@1: alpar@1: static void chuzc(struct csa *csa, double rtol) alpar@1: { alpar@1: #ifdef GLP_DEBUG alpar@1: int m = csa->m; alpar@1: int n = csa->n; alpar@1: #endif alpar@1: char *stat = csa->stat; alpar@1: double *cbar = csa->cbar; alpar@1: #ifdef GLP_DEBUG alpar@1: int p = csa->p; alpar@1: #endif alpar@1: double delta = csa->delta; alpar@1: int *trow_ind = csa->trow_ind; alpar@1: double *trow_vec = csa->trow_vec; alpar@1: int trow_num = csa->trow_num; alpar@1: int j, pos, q; alpar@1: double alfa, big, s, t, teta, tmax; alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= p && p <= m); alpar@1: #endif alpar@1: /* delta > 0 means that xB[p] violates its lower bound and goes alpar@1: to it in the adjacent basis, so lambdaB[p] is increasing from alpar@1: its lower zero bound; alpar@1: delta < 0 means that xB[p] violates its upper bound and goes alpar@1: to it in the adjacent basis, so lambdaB[p] is decreasing from alpar@1: its upper zero bound */ alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(delta != 0.0); alpar@1: #endif alpar@1: /* s := sign(delta) */ alpar@1: s = (delta > 0.0 ? +1.0 : -1.0); alpar@1: /*** FIRST PASS ***/ alpar@1: /* nothing is chosen so far */ alpar@1: q = 0, teta = DBL_MAX, big = 0.0; alpar@1: /* walk through significant elements of the pivot row */ alpar@1: for (pos = 1; pos <= trow_num; pos++) alpar@1: { j = trow_ind[pos]; alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= j && j <= n); alpar@1: #endif alpar@1: alfa = s * trow_vec[j]; alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(alfa != 0.0); alpar@1: #endif alpar@1: /* lambdaN[j] = ... - alfa * lambdaB[p] - ..., and due to s we alpar@1: need to consider only increasing lambdaB[p] */ alpar@1: if (alfa > 0.0) alpar@1: { /* lambdaN[j] is decreasing */ alpar@1: if (stat[j] == GLP_NL || stat[j] == GLP_NF) alpar@1: { /* lambdaN[j] has zero lower bound */ alpar@1: t = (cbar[j] + rtol) / alfa; alpar@1: } alpar@1: else alpar@1: { /* lambdaN[j] has no lower bound */ alpar@1: continue; alpar@1: } alpar@1: } alpar@1: else alpar@1: { /* lambdaN[j] is increasing */ alpar@1: if (stat[j] == GLP_NU || stat[j] == GLP_NF) alpar@1: { /* lambdaN[j] has zero upper bound */ alpar@1: t = (cbar[j] - rtol) / alfa; alpar@1: } alpar@1: else alpar@1: { /* lambdaN[j] has no upper bound */ alpar@1: continue; alpar@1: } alpar@1: } alpar@1: /* t is a change of lambdaB[p], on which lambdaN[j] reaches alpar@1: its zero bound (possibly relaxed); since the basic solution alpar@1: is assumed to be dual feasible, t has to be non-negative by alpar@1: definition; however, it may happen that lambdaN[j] slightly alpar@1: (i.e. within a tolerance) violates its zero bound, that alpar@1: leads to negative t; in the latter case, if xN[j] is chosen, alpar@1: negative t means that lambdaB[p] changes in wrong direction alpar@1: that may cause wrong results on updating reduced costs; alpar@1: thus, if t is negative, we should replace it by exact zero alpar@1: assuming that lambdaN[j] is exactly on its zero bound, and alpar@1: violation appears due to round-off errors */ alpar@1: if (t < 0.0) t = 0.0; alpar@1: /* apply minimal ratio test */ alpar@1: if (teta > t || teta == t && big < fabs(alfa)) alpar@1: q = j, teta = t, big = fabs(alfa); alpar@1: } alpar@1: /* the second pass is skipped in the following cases: */ alpar@1: /* if the standard ratio test is used */ alpar@1: if (rtol == 0.0) goto done; alpar@1: /* if no non-basic variable has been chosen on the first pass */ alpar@1: if (q == 0) goto done; alpar@1: /* if lambdaN[q] prevents lambdaB[p] from any change */ alpar@1: if (teta == 0.0) goto done; alpar@1: /*** SECOND PASS ***/ alpar@1: /* here tmax is a maximal change of lambdaB[p], on which the alpar@1: solution remains dual feasible within a tolerance */ alpar@1: #if 0 alpar@1: tmax = (1.0 + 10.0 * DBL_EPSILON) * teta; alpar@1: #else alpar@1: tmax = teta; alpar@1: #endif alpar@1: /* nothing is chosen so far */ alpar@1: q = 0, teta = DBL_MAX, big = 0.0; alpar@1: /* walk through significant elements of the pivot row */ alpar@1: for (pos = 1; pos <= trow_num; pos++) alpar@1: { j = trow_ind[pos]; alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= j && j <= n); alpar@1: #endif alpar@1: alfa = s * trow_vec[j]; alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(alfa != 0.0); alpar@1: #endif alpar@1: /* lambdaN[j] = ... - alfa * lambdaB[p] - ..., and due to s we alpar@1: need to consider only increasing lambdaB[p] */ alpar@1: if (alfa > 0.0) alpar@1: { /* lambdaN[j] is decreasing */ alpar@1: if (stat[j] == GLP_NL || stat[j] == GLP_NF) alpar@1: { /* lambdaN[j] has zero lower bound */ alpar@1: t = cbar[j] / alfa; alpar@1: } alpar@1: else alpar@1: { /* lambdaN[j] has no lower bound */ alpar@1: continue; alpar@1: } alpar@1: } alpar@1: else alpar@1: { /* lambdaN[j] is increasing */ alpar@1: if (stat[j] == GLP_NU || stat[j] == GLP_NF) alpar@1: { /* lambdaN[j] has zero upper bound */ alpar@1: t = cbar[j] / alfa; alpar@1: } alpar@1: else alpar@1: { /* lambdaN[j] has no upper bound */ alpar@1: continue; alpar@1: } alpar@1: } alpar@1: /* (see comments for the first pass) */ alpar@1: if (t < 0.0) t = 0.0; alpar@1: /* t is a change of lambdaB[p], on which lambdaN[j] reaches alpar@1: its zero (lower or upper) bound; if t <= tmax, all reduced alpar@1: costs can violate their zero bounds only within relaxation alpar@1: tolerance rtol, so we can choose non-basic variable having alpar@1: largest influence coefficient to avoid possible numerical alpar@1: instability */ alpar@1: if (t <= tmax && big < fabs(alfa)) alpar@1: q = j, teta = t, big = fabs(alfa); alpar@1: } alpar@1: /* something must be chosen on the second pass */ alpar@1: xassert(q != 0); alpar@1: done: /* store the index of non-basic variable xN[q] chosen */ alpar@1: csa->q = q; alpar@1: /* store reduced cost of xN[q] in the adjacent basis */ alpar@1: csa->new_dq = s * teta; alpar@1: return; alpar@1: } alpar@1: alpar@1: #if 1 /* copied from primal */ alpar@1: /*********************************************************************** alpar@1: * eval_tcol - compute pivot column of the simplex table alpar@1: * alpar@1: * This routine computes the pivot column of the simplex table, which alpar@1: * corresponds to non-basic variable xN[q] chosen. alpar@1: * alpar@1: * The pivot column is the following vector: alpar@1: * alpar@1: * tcol = T * e[q] = - inv(B) * N * e[q] = - inv(B) * N[q], alpar@1: * alpar@1: * where B is the current basis matrix, N[q] is a column of the matrix alpar@1: * (I|-A) corresponding to variable xN[q]. */ alpar@1: alpar@1: static void eval_tcol(struct csa *csa) alpar@1: { int m = csa->m; alpar@1: #ifdef GLP_DEBUG alpar@1: int n = csa->n; alpar@1: #endif alpar@1: int *head = csa->head; alpar@1: int q = csa->q; alpar@1: int *tcol_ind = csa->tcol_ind; alpar@1: double *tcol_vec = csa->tcol_vec; alpar@1: double *h = csa->tcol_vec; alpar@1: int i, k, nnz; alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= q && q <= n); alpar@1: #endif alpar@1: k = head[m+q]; /* x[k] = xN[q] */ alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= k && k <= m+n); alpar@1: #endif alpar@1: /* construct the right-hand side vector h = - N[q] */ alpar@1: for (i = 1; i <= m; i++) alpar@1: h[i] = 0.0; alpar@1: if (k <= m) alpar@1: { /* N[q] is k-th column of submatrix I */ alpar@1: h[k] = -1.0; alpar@1: } alpar@1: else alpar@1: { /* N[q] is (k-m)-th column of submatrix (-A) */ alpar@1: int *A_ptr = csa->A_ptr; alpar@1: int *A_ind = csa->A_ind; alpar@1: double *A_val = csa->A_val; alpar@1: int beg, end, ptr; alpar@1: beg = A_ptr[k-m]; alpar@1: end = A_ptr[k-m+1]; alpar@1: for (ptr = beg; ptr < end; ptr++) alpar@1: h[A_ind[ptr]] = A_val[ptr]; alpar@1: } alpar@1: /* solve system B * tcol = h */ alpar@1: xassert(csa->valid); alpar@1: bfd_ftran(csa->bfd, tcol_vec); alpar@1: /* construct sparse pattern of the pivot column */ alpar@1: nnz = 0; alpar@1: for (i = 1; i <= m; i++) alpar@1: { if (tcol_vec[i] != 0.0) alpar@1: tcol_ind[++nnz] = i; alpar@1: } alpar@1: csa->tcol_nnz = nnz; alpar@1: return; alpar@1: } alpar@1: #endif alpar@1: alpar@1: #if 1 /* copied from primal */ alpar@1: /*********************************************************************** alpar@1: * refine_tcol - refine pivot column of the simplex table alpar@1: * alpar@1: * This routine refines the pivot column of the simplex table assuming alpar@1: * that it was previously computed by the routine eval_tcol. */ alpar@1: alpar@1: static void refine_tcol(struct csa *csa) alpar@1: { int m = csa->m; alpar@1: #ifdef GLP_DEBUG alpar@1: int n = csa->n; alpar@1: #endif alpar@1: int *head = csa->head; alpar@1: int q = csa->q; alpar@1: int *tcol_ind = csa->tcol_ind; alpar@1: double *tcol_vec = csa->tcol_vec; alpar@1: double *h = csa->work3; alpar@1: int i, k, nnz; alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= q && q <= n); alpar@1: #endif alpar@1: k = head[m+q]; /* x[k] = xN[q] */ alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= k && k <= m+n); alpar@1: #endif alpar@1: /* construct the right-hand side vector h = - N[q] */ alpar@1: for (i = 1; i <= m; i++) alpar@1: h[i] = 0.0; alpar@1: if (k <= m) alpar@1: { /* N[q] is k-th column of submatrix I */ alpar@1: h[k] = -1.0; alpar@1: } alpar@1: else alpar@1: { /* N[q] is (k-m)-th column of submatrix (-A) */ alpar@1: int *A_ptr = csa->A_ptr; alpar@1: int *A_ind = csa->A_ind; alpar@1: double *A_val = csa->A_val; alpar@1: int beg, end, ptr; alpar@1: beg = A_ptr[k-m]; alpar@1: end = A_ptr[k-m+1]; alpar@1: for (ptr = beg; ptr < end; ptr++) alpar@1: h[A_ind[ptr]] = A_val[ptr]; alpar@1: } alpar@1: /* refine solution of B * tcol = h */ alpar@1: refine_ftran(csa, h, tcol_vec); alpar@1: /* construct sparse pattern of the pivot column */ alpar@1: nnz = 0; alpar@1: for (i = 1; i <= m; i++) alpar@1: { if (tcol_vec[i] != 0.0) alpar@1: tcol_ind[++nnz] = i; alpar@1: } alpar@1: csa->tcol_nnz = nnz; alpar@1: return; alpar@1: } alpar@1: #endif alpar@1: alpar@1: /*********************************************************************** alpar@1: * update_cbar - update reduced costs of non-basic variables alpar@1: * alpar@1: * This routine updates reduced costs of all (except fixed) non-basic alpar@1: * variables for the adjacent basis. */ alpar@1: alpar@1: static void update_cbar(struct csa *csa) alpar@1: { alpar@1: #ifdef GLP_DEBUG alpar@1: int n = csa->n; alpar@1: #endif alpar@1: double *cbar = csa->cbar; alpar@1: int trow_nnz = csa->trow_nnz; alpar@1: int *trow_ind = csa->trow_ind; alpar@1: double *trow_vec = csa->trow_vec; alpar@1: int q = csa->q; alpar@1: double new_dq = csa->new_dq; alpar@1: int j, pos; alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= q && q <= n); alpar@1: #endif alpar@1: /* set new reduced cost of xN[q] */ alpar@1: cbar[q] = new_dq; alpar@1: /* update reduced costs of other non-basic variables */ alpar@1: if (new_dq == 0.0) goto done; alpar@1: for (pos = 1; pos <= trow_nnz; pos++) alpar@1: { j = trow_ind[pos]; alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= j && j <= n); alpar@1: #endif alpar@1: if (j != q) alpar@1: cbar[j] -= trow_vec[j] * new_dq; alpar@1: } alpar@1: done: return; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * update_bbar - update values of basic variables alpar@1: * alpar@1: * This routine updates values of all basic variables for the adjacent alpar@1: * basis. */ alpar@1: alpar@1: static void update_bbar(struct csa *csa) alpar@1: { alpar@1: #ifdef GLP_DEBUG alpar@1: int m = csa->m; alpar@1: int n = csa->n; alpar@1: #endif alpar@1: double *bbar = csa->bbar; alpar@1: int p = csa->p; alpar@1: double delta = csa->delta; alpar@1: int q = csa->q; alpar@1: int tcol_nnz = csa->tcol_nnz; alpar@1: int *tcol_ind = csa->tcol_ind; alpar@1: double *tcol_vec = csa->tcol_vec; alpar@1: int i, pos; alpar@1: double teta; alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= p && p <= m); alpar@1: xassert(1 <= q && q <= n); alpar@1: #endif alpar@1: /* determine the change of xN[q] in the adjacent basis */ alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(tcol_vec[p] != 0.0); alpar@1: #endif alpar@1: teta = delta / tcol_vec[p]; alpar@1: /* set new primal value of xN[q] */ alpar@1: bbar[p] = get_xN(csa, q) + teta; alpar@1: /* update primal values of other basic variables */ alpar@1: if (teta == 0.0) goto done; alpar@1: for (pos = 1; pos <= tcol_nnz; pos++) alpar@1: { i = tcol_ind[pos]; alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= i && i <= m); alpar@1: #endif alpar@1: if (i != p) alpar@1: bbar[i] += tcol_vec[i] * teta; alpar@1: } alpar@1: done: return; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * update_gamma - update steepest edge coefficients alpar@1: * alpar@1: * This routine updates steepest-edge coefficients for the adjacent alpar@1: * basis. */ alpar@1: alpar@1: static void update_gamma(struct csa *csa) alpar@1: { int m = csa->m; alpar@1: #ifdef GLP_DEBUG alpar@1: int n = csa->n; alpar@1: #endif alpar@1: char *type = csa->type; alpar@1: int *head = csa->head; alpar@1: char *refsp = csa->refsp; alpar@1: double *gamma = csa->gamma; alpar@1: int p = csa->p; alpar@1: int trow_nnz = csa->trow_nnz; alpar@1: int *trow_ind = csa->trow_ind; alpar@1: double *trow_vec = csa->trow_vec; alpar@1: int q = csa->q; alpar@1: int tcol_nnz = csa->tcol_nnz; alpar@1: int *tcol_ind = csa->tcol_ind; alpar@1: double *tcol_vec = csa->tcol_vec; alpar@1: double *u = csa->work3; alpar@1: int i, j, k,pos; alpar@1: double gamma_p, eta_p, pivot, t, t1, t2; alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= p && p <= m); alpar@1: xassert(1 <= q && q <= n); alpar@1: #endif alpar@1: /* the basis changes, so decrease the count */ alpar@1: xassert(csa->refct > 0); alpar@1: csa->refct--; alpar@1: /* recompute gamma[p] for the current basis more accurately and alpar@1: compute auxiliary vector u */ alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(type[head[p]] != GLP_FR); alpar@1: #endif alpar@1: gamma_p = eta_p = (refsp[head[p]] ? 1.0 : 0.0); alpar@1: for (i = 1; i <= m; i++) u[i] = 0.0; alpar@1: for (pos = 1; pos <= trow_nnz; pos++) alpar@1: { j = trow_ind[pos]; alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= j && j <= n); alpar@1: #endif alpar@1: k = head[m+j]; /* x[k] = xN[j] */ alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= k && k <= m+n); alpar@1: xassert(type[k] != GLP_FX); alpar@1: #endif alpar@1: if (!refsp[k]) continue; alpar@1: t = trow_vec[j]; alpar@1: gamma_p += t * t; alpar@1: /* u := u + N[j] * delta[j] * trow[j] */ alpar@1: if (k <= m) alpar@1: { /* N[k] = k-j stolbec submatrix I */ alpar@1: u[k] += t; alpar@1: } alpar@1: else alpar@1: { /* N[k] = k-m-k stolbec (-A) */ alpar@1: int *A_ptr = csa->A_ptr; alpar@1: int *A_ind = csa->A_ind; alpar@1: double *A_val = csa->A_val; alpar@1: int beg, end, ptr; alpar@1: beg = A_ptr[k-m]; alpar@1: end = A_ptr[k-m+1]; alpar@1: for (ptr = beg; ptr < end; ptr++) alpar@1: u[A_ind[ptr]] -= t * A_val[ptr]; alpar@1: } alpar@1: } alpar@1: xassert(csa->valid); alpar@1: bfd_ftran(csa->bfd, u); alpar@1: /* update gamma[i] for other basic variables (except xB[p] and alpar@1: free variables) */ alpar@1: pivot = tcol_vec[p]; alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(pivot != 0.0); alpar@1: #endif alpar@1: for (pos = 1; pos <= tcol_nnz; pos++) alpar@1: { i = tcol_ind[pos]; alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= i && i <= m); alpar@1: #endif alpar@1: k = head[i]; alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= k && k <= m+n); alpar@1: #endif alpar@1: /* skip xB[p] */ alpar@1: if (i == p) continue; alpar@1: /* skip free basic variable */ alpar@1: if (type[head[i]] == GLP_FR) alpar@1: { alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(gamma[i] == 1.0); alpar@1: #endif alpar@1: continue; alpar@1: } alpar@1: /* compute gamma[i] for the adjacent basis */ alpar@1: t = tcol_vec[i] / pivot; alpar@1: t1 = gamma[i] + t * t * gamma_p + 2.0 * t * u[i]; alpar@1: t2 = (refsp[k] ? 1.0 : 0.0) + eta_p * t * t; alpar@1: gamma[i] = (t1 >= t2 ? t1 : t2); alpar@1: /* (though gamma[i] can be exact zero, because the reference alpar@1: space does not include non-basic fixed variables) */ alpar@1: if (gamma[i] < DBL_EPSILON) gamma[i] = DBL_EPSILON; alpar@1: } alpar@1: /* compute gamma[p] for the adjacent basis */ alpar@1: if (type[head[m+q]] == GLP_FR) alpar@1: gamma[p] = 1.0; alpar@1: else alpar@1: { gamma[p] = gamma_p / (pivot * pivot); alpar@1: if (gamma[p] < DBL_EPSILON) gamma[p] = DBL_EPSILON; alpar@1: } alpar@1: /* if xB[p], which becomes xN[q] in the adjacent basis, is fixed alpar@1: and belongs to the reference space, remove it from there, and alpar@1: change all gamma's appropriately */ alpar@1: k = head[p]; alpar@1: if (type[k] == GLP_FX && refsp[k]) alpar@1: { refsp[k] = 0; alpar@1: for (pos = 1; pos <= tcol_nnz; pos++) alpar@1: { i = tcol_ind[pos]; alpar@1: if (i == p) alpar@1: { if (type[head[m+q]] == GLP_FR) continue; alpar@1: t = 1.0 / tcol_vec[p]; alpar@1: } alpar@1: else alpar@1: { if (type[head[i]] == GLP_FR) continue; alpar@1: t = tcol_vec[i] / tcol_vec[p]; alpar@1: } alpar@1: gamma[i] -= t * t; alpar@1: if (gamma[i] < DBL_EPSILON) gamma[i] = DBL_EPSILON; alpar@1: } alpar@1: } alpar@1: return; alpar@1: } alpar@1: alpar@1: #if 1 /* copied from primal */ alpar@1: /*********************************************************************** alpar@1: * err_in_bbar - compute maximal relative error in primal solution alpar@1: * alpar@1: * This routine returns maximal relative error: alpar@1: * alpar@1: * max |beta[i] - bbar[i]| / (1 + |beta[i]|), alpar@1: * alpar@1: * where beta and bbar are, respectively, directly computed and the alpar@1: * current (updated) values of basic variables. alpar@1: * alpar@1: * NOTE: The routine is intended only for debugginig purposes. */ alpar@1: alpar@1: static double err_in_bbar(struct csa *csa) alpar@1: { int m = csa->m; alpar@1: double *bbar = csa->bbar; alpar@1: int i; alpar@1: double e, emax, *beta; alpar@1: beta = xcalloc(1+m, sizeof(double)); alpar@1: eval_beta(csa, beta); alpar@1: emax = 0.0; alpar@1: for (i = 1; i <= m; i++) alpar@1: { e = fabs(beta[i] - bbar[i]) / (1.0 + fabs(beta[i])); alpar@1: if (emax < e) emax = e; alpar@1: } alpar@1: xfree(beta); alpar@1: return emax; alpar@1: } alpar@1: #endif alpar@1: alpar@1: #if 1 /* copied from primal */ alpar@1: /*********************************************************************** alpar@1: * err_in_cbar - compute maximal relative error in dual solution alpar@1: * alpar@1: * This routine returns maximal relative error: alpar@1: * alpar@1: * max |cost[j] - cbar[j]| / (1 + |cost[j]|), alpar@1: * alpar@1: * where cost and cbar are, respectively, directly computed and the alpar@1: * current (updated) reduced costs of non-basic non-fixed variables. alpar@1: * alpar@1: * NOTE: The routine is intended only for debugginig purposes. */ alpar@1: alpar@1: static double err_in_cbar(struct csa *csa) alpar@1: { int m = csa->m; alpar@1: int n = csa->n; alpar@1: char *stat = csa->stat; alpar@1: double *cbar = csa->cbar; alpar@1: int j; alpar@1: double e, emax, cost, *pi; alpar@1: pi = xcalloc(1+m, sizeof(double)); alpar@1: eval_pi(csa, pi); alpar@1: emax = 0.0; alpar@1: for (j = 1; j <= n; j++) alpar@1: { if (stat[j] == GLP_NS) continue; alpar@1: cost = eval_cost(csa, pi, j); alpar@1: e = fabs(cost - cbar[j]) / (1.0 + fabs(cost)); alpar@1: if (emax < e) emax = e; alpar@1: } alpar@1: xfree(pi); alpar@1: return emax; alpar@1: } alpar@1: #endif alpar@1: alpar@1: /*********************************************************************** alpar@1: * err_in_gamma - compute maximal relative error in steepest edge cff. alpar@1: * alpar@1: * This routine returns maximal relative error: alpar@1: * alpar@1: * max |gamma'[j] - gamma[j]| / (1 + |gamma'[j]), alpar@1: * alpar@1: * where gamma'[j] and gamma[j] are, respectively, directly computed alpar@1: * and the current (updated) steepest edge coefficients for non-basic alpar@1: * non-fixed variable x[j]. alpar@1: * alpar@1: * NOTE: The routine is intended only for debugginig purposes. */ alpar@1: alpar@1: static double err_in_gamma(struct csa *csa) alpar@1: { int m = csa->m; alpar@1: char *type = csa->type; alpar@1: int *head = csa->head; alpar@1: double *gamma = csa->gamma; alpar@1: double *exact = csa->work4; alpar@1: int i; alpar@1: double e, emax, temp; alpar@1: eval_gamma(csa, exact); alpar@1: emax = 0.0; alpar@1: for (i = 1; i <= m; i++) alpar@1: { if (type[head[i]] == GLP_FR) alpar@1: { xassert(gamma[i] == 1.0); alpar@1: xassert(exact[i] == 1.0); alpar@1: continue; alpar@1: } alpar@1: temp = exact[i]; alpar@1: e = fabs(temp - gamma[i]) / (1.0 + fabs(temp)); alpar@1: if (emax < e) emax = e; alpar@1: } alpar@1: return emax; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * change_basis - change basis header alpar@1: * alpar@1: * This routine changes the basis header to make it corresponding to alpar@1: * the adjacent basis. */ alpar@1: alpar@1: static void change_basis(struct csa *csa) alpar@1: { int m = csa->m; alpar@1: #ifdef GLP_DEBUG alpar@1: int n = csa->n; alpar@1: #endif alpar@1: char *type = csa->type; alpar@1: int *head = csa->head; alpar@1: #if 1 /* 06/IV-2009 */ alpar@1: int *bind = csa->bind; alpar@1: #endif alpar@1: char *stat = csa->stat; alpar@1: int p = csa->p; alpar@1: double delta = csa->delta; alpar@1: int q = csa->q; alpar@1: int k; alpar@1: /* xB[p] leaves the basis, xN[q] enters the basis */ alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= p && p <= m); alpar@1: xassert(1 <= q && q <= n); alpar@1: #endif alpar@1: /* xB[p] <-> xN[q] */ alpar@1: k = head[p], head[p] = head[m+q], head[m+q] = k; alpar@1: #if 1 /* 06/IV-2009 */ alpar@1: bind[head[p]] = p, bind[head[m+q]] = m + q; alpar@1: #endif alpar@1: if (type[k] == GLP_FX) alpar@1: stat[q] = GLP_NS; alpar@1: else if (delta > 0.0) alpar@1: { alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(type[k] == GLP_LO || type[k] == GLP_DB); alpar@1: #endif alpar@1: stat[q] = GLP_NL; alpar@1: } alpar@1: else /* delta < 0.0 */ alpar@1: { alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(type[k] == GLP_UP || type[k] == GLP_DB); alpar@1: #endif alpar@1: stat[q] = GLP_NU; alpar@1: } alpar@1: return; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * check_feas - check dual feasibility of basic solution alpar@1: * alpar@1: * If the current basic solution is dual feasible within a tolerance, alpar@1: * this routine returns zero, otherwise it returns non-zero. */ alpar@1: alpar@1: static int check_feas(struct csa *csa, double tol_dj) alpar@1: { int m = csa->m; alpar@1: int n = csa->n; alpar@1: char *orig_type = csa->orig_type; alpar@1: int *head = csa->head; alpar@1: double *cbar = csa->cbar; alpar@1: int j, k; alpar@1: for (j = 1; j <= n; j++) alpar@1: { k = head[m+j]; /* x[k] = xN[j] */ alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= k && k <= m+n); alpar@1: #endif alpar@1: if (cbar[j] < - tol_dj) alpar@1: if (orig_type[k] == GLP_LO || orig_type[k] == GLP_FR) alpar@1: return 1; alpar@1: if (cbar[j] > + tol_dj) alpar@1: if (orig_type[k] == GLP_UP || orig_type[k] == GLP_FR) alpar@1: return 1; alpar@1: } alpar@1: return 0; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * set_aux_bnds - assign auxiliary bounds to variables alpar@1: * alpar@1: * This routine assigns auxiliary bounds to variables to construct an alpar@1: * LP problem solved on phase I. */ alpar@1: alpar@1: static void set_aux_bnds(struct csa *csa) alpar@1: { int m = csa->m; alpar@1: int n = csa->n; alpar@1: char *type = csa->type; alpar@1: double *lb = csa->lb; alpar@1: double *ub = csa->ub; alpar@1: char *orig_type = csa->orig_type; alpar@1: int *head = csa->head; alpar@1: char *stat = csa->stat; alpar@1: double *cbar = csa->cbar; alpar@1: int j, k; alpar@1: for (k = 1; k <= m+n; k++) alpar@1: { switch (orig_type[k]) alpar@1: { case GLP_FR: alpar@1: #if 0 alpar@1: type[k] = GLP_DB, lb[k] = -1.0, ub[k] = +1.0; alpar@1: #else alpar@1: /* to force free variables to enter the basis */ alpar@1: type[k] = GLP_DB, lb[k] = -1e3, ub[k] = +1e3; alpar@1: #endif alpar@1: break; alpar@1: case GLP_LO: alpar@1: type[k] = GLP_DB, lb[k] = 0.0, ub[k] = +1.0; alpar@1: break; alpar@1: case GLP_UP: alpar@1: type[k] = GLP_DB, lb[k] = -1.0, ub[k] = 0.0; alpar@1: break; alpar@1: case GLP_DB: alpar@1: case GLP_FX: alpar@1: type[k] = GLP_FX, lb[k] = ub[k] = 0.0; alpar@1: break; alpar@1: default: alpar@1: xassert(orig_type != orig_type); alpar@1: } alpar@1: } alpar@1: for (j = 1; j <= n; j++) alpar@1: { k = head[m+j]; /* x[k] = xN[j] */ alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= k && k <= m+n); alpar@1: #endif alpar@1: if (type[k] == GLP_FX) alpar@1: stat[j] = GLP_NS; alpar@1: else if (cbar[j] >= 0.0) alpar@1: stat[j] = GLP_NL; alpar@1: else alpar@1: stat[j] = GLP_NU; alpar@1: } alpar@1: return; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * set_orig_bnds - restore original bounds of variables alpar@1: * alpar@1: * This routine restores original types and bounds of variables and alpar@1: * determines statuses of non-basic variables assuming that the current alpar@1: * basis is dual feasible. */ alpar@1: alpar@1: static void set_orig_bnds(struct csa *csa) alpar@1: { int m = csa->m; alpar@1: int n = csa->n; alpar@1: char *type = csa->type; alpar@1: double *lb = csa->lb; alpar@1: double *ub = csa->ub; alpar@1: char *orig_type = csa->orig_type; alpar@1: double *orig_lb = csa->orig_lb; alpar@1: double *orig_ub = csa->orig_ub; alpar@1: int *head = csa->head; alpar@1: char *stat = csa->stat; alpar@1: double *cbar = csa->cbar; alpar@1: int j, k; alpar@1: memcpy(&type[1], &orig_type[1], (m+n) * sizeof(char)); alpar@1: memcpy(&lb[1], &orig_lb[1], (m+n) * sizeof(double)); alpar@1: memcpy(&ub[1], &orig_ub[1], (m+n) * sizeof(double)); alpar@1: for (j = 1; j <= n; j++) alpar@1: { k = head[m+j]; /* x[k] = xN[j] */ alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= k && k <= m+n); alpar@1: #endif alpar@1: switch (type[k]) alpar@1: { case GLP_FR: alpar@1: stat[j] = GLP_NF; alpar@1: break; alpar@1: case GLP_LO: alpar@1: stat[j] = GLP_NL; alpar@1: break; alpar@1: case GLP_UP: alpar@1: stat[j] = GLP_NU; alpar@1: break; alpar@1: case GLP_DB: alpar@1: if (cbar[j] >= +DBL_EPSILON) alpar@1: stat[j] = GLP_NL; alpar@1: else if (cbar[j] <= -DBL_EPSILON) alpar@1: stat[j] = GLP_NU; alpar@1: else if (fabs(lb[k]) <= fabs(ub[k])) alpar@1: stat[j] = GLP_NL; alpar@1: else alpar@1: stat[j] = GLP_NU; alpar@1: break; alpar@1: case GLP_FX: alpar@1: stat[j] = GLP_NS; alpar@1: break; alpar@1: default: alpar@1: xassert(type != type); alpar@1: } alpar@1: } alpar@1: return; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * check_stab - check numerical stability of basic solution alpar@1: * alpar@1: * If the current basic solution is dual feasible within a tolerance, alpar@1: * this routine returns zero, otherwise it returns non-zero. */ alpar@1: alpar@1: static int check_stab(struct csa *csa, double tol_dj) alpar@1: { int n = csa->n; alpar@1: char *stat = csa->stat; alpar@1: double *cbar = csa->cbar; alpar@1: int j; alpar@1: for (j = 1; j <= n; j++) alpar@1: { if (cbar[j] < - tol_dj) alpar@1: if (stat[j] == GLP_NL || stat[j] == GLP_NF) return 1; alpar@1: if (cbar[j] > + tol_dj) alpar@1: if (stat[j] == GLP_NU || stat[j] == GLP_NF) return 1; alpar@1: } alpar@1: return 0; alpar@1: } alpar@1: alpar@1: #if 1 /* copied from primal */ alpar@1: /*********************************************************************** alpar@1: * eval_obj - compute original objective function alpar@1: * alpar@1: * This routine computes the current value of the original objective alpar@1: * function. */ alpar@1: alpar@1: static double eval_obj(struct csa *csa) alpar@1: { int m = csa->m; alpar@1: int n = csa->n; alpar@1: double *obj = csa->obj; alpar@1: int *head = csa->head; alpar@1: double *bbar = csa->bbar; alpar@1: int i, j, k; alpar@1: double sum; alpar@1: sum = obj[0]; alpar@1: /* walk through the list of basic variables */ alpar@1: for (i = 1; i <= m; i++) alpar@1: { k = head[i]; /* x[k] = xB[i] */ alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= k && k <= m+n); alpar@1: #endif alpar@1: if (k > m) alpar@1: sum += obj[k-m] * bbar[i]; alpar@1: } alpar@1: /* walk through the list of non-basic variables */ alpar@1: for (j = 1; j <= n; j++) alpar@1: { k = head[m+j]; /* x[k] = xN[j] */ alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= k && k <= m+n); alpar@1: #endif alpar@1: if (k > m) alpar@1: sum += obj[k-m] * get_xN(csa, j); alpar@1: } alpar@1: return sum; alpar@1: } alpar@1: #endif alpar@1: alpar@1: /*********************************************************************** alpar@1: * display - display the search progress alpar@1: * alpar@1: * This routine displays some information about the search progress. */ alpar@1: alpar@1: static void display(struct csa *csa, const glp_smcp *parm, int spec) alpar@1: { int m = csa->m; alpar@1: int n = csa->n; alpar@1: double *coef = csa->coef; alpar@1: char *orig_type = csa->orig_type; alpar@1: int *head = csa->head; alpar@1: char *stat = csa->stat; alpar@1: int phase = csa->phase; alpar@1: double *bbar = csa->bbar; alpar@1: double *cbar = csa->cbar; alpar@1: int i, j, cnt; alpar@1: double sum; alpar@1: if (parm->msg_lev < GLP_MSG_ON) goto skip; alpar@1: if (parm->out_dly > 0 && alpar@1: 1000.0 * xdifftime(xtime(), csa->tm_beg) < parm->out_dly) alpar@1: goto skip; alpar@1: if (csa->it_cnt == csa->it_dpy) goto skip; alpar@1: if (!spec && csa->it_cnt % parm->out_frq != 0) goto skip; alpar@1: /* compute the sum of dual infeasibilities */ alpar@1: sum = 0.0; alpar@1: if (phase == 1) alpar@1: { for (i = 1; i <= m; i++) alpar@1: sum -= coef[head[i]] * bbar[i]; alpar@1: for (j = 1; j <= n; j++) alpar@1: sum -= coef[head[m+j]] * get_xN(csa, j); alpar@1: } alpar@1: else alpar@1: { for (j = 1; j <= n; j++) alpar@1: { if (cbar[j] < 0.0) alpar@1: if (stat[j] == GLP_NL || stat[j] == GLP_NF) alpar@1: sum -= cbar[j]; alpar@1: if (cbar[j] > 0.0) alpar@1: if (stat[j] == GLP_NU || stat[j] == GLP_NF) alpar@1: sum += cbar[j]; alpar@1: } alpar@1: } alpar@1: /* determine the number of basic fixed variables */ alpar@1: cnt = 0; alpar@1: for (i = 1; i <= m; i++) alpar@1: if (orig_type[head[i]] == GLP_FX) cnt++; alpar@1: if (csa->phase == 1) alpar@1: xprintf(" %6d: %24s infeas = %10.3e (%d)\n", alpar@1: csa->it_cnt, "", sum, cnt); alpar@1: else alpar@1: xprintf("|%6d: obj = %17.9e infeas = %10.3e (%d)\n", alpar@1: csa->it_cnt, eval_obj(csa), sum, cnt); alpar@1: csa->it_dpy = csa->it_cnt; alpar@1: skip: return; alpar@1: } alpar@1: alpar@1: #if 1 /* copied from primal */ alpar@1: /*********************************************************************** alpar@1: * store_sol - store basic solution back to the problem object alpar@1: * alpar@1: * This routine stores basic solution components back to the problem alpar@1: * object. */ alpar@1: alpar@1: static void store_sol(struct csa *csa, glp_prob *lp, int p_stat, alpar@1: int d_stat, int ray) alpar@1: { int m = csa->m; alpar@1: int n = csa->n; alpar@1: double zeta = csa->zeta; alpar@1: int *head = csa->head; alpar@1: char *stat = csa->stat; alpar@1: double *bbar = csa->bbar; alpar@1: double *cbar = csa->cbar; alpar@1: int i, j, k; alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(lp->m == m); alpar@1: xassert(lp->n == n); alpar@1: #endif alpar@1: /* basis factorization */ alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(!lp->valid && lp->bfd == NULL); alpar@1: xassert(csa->valid && csa->bfd != NULL); alpar@1: #endif alpar@1: lp->valid = 1, csa->valid = 0; alpar@1: lp->bfd = csa->bfd, csa->bfd = NULL; alpar@1: memcpy(&lp->head[1], &head[1], m * sizeof(int)); alpar@1: /* basic solution status */ alpar@1: lp->pbs_stat = p_stat; alpar@1: lp->dbs_stat = d_stat; alpar@1: /* objective function value */ alpar@1: lp->obj_val = eval_obj(csa); alpar@1: /* simplex iteration count */ alpar@1: lp->it_cnt = csa->it_cnt; alpar@1: /* unbounded ray */ alpar@1: lp->some = ray; alpar@1: /* basic variables */ alpar@1: for (i = 1; i <= m; i++) alpar@1: { k = head[i]; /* x[k] = xB[i] */ alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= k && k <= m+n); alpar@1: #endif alpar@1: if (k <= m) alpar@1: { GLPROW *row = lp->row[k]; alpar@1: row->stat = GLP_BS; alpar@1: row->bind = i; alpar@1: row->prim = bbar[i] / row->rii; alpar@1: row->dual = 0.0; alpar@1: } alpar@1: else alpar@1: { GLPCOL *col = lp->col[k-m]; alpar@1: col->stat = GLP_BS; alpar@1: col->bind = i; alpar@1: col->prim = bbar[i] * col->sjj; alpar@1: col->dual = 0.0; alpar@1: } alpar@1: } alpar@1: /* non-basic variables */ alpar@1: for (j = 1; j <= n; j++) alpar@1: { k = head[m+j]; /* x[k] = xN[j] */ alpar@1: #ifdef GLP_DEBUG alpar@1: xassert(1 <= k && k <= m+n); alpar@1: #endif alpar@1: if (k <= m) alpar@1: { GLPROW *row = lp->row[k]; alpar@1: row->stat = stat[j]; alpar@1: row->bind = 0; alpar@1: #if 0 alpar@1: row->prim = get_xN(csa, j) / row->rii; alpar@1: #else alpar@1: switch (stat[j]) alpar@1: { case GLP_NL: alpar@1: row->prim = row->lb; break; alpar@1: case GLP_NU: alpar@1: row->prim = row->ub; break; alpar@1: case GLP_NF: alpar@1: row->prim = 0.0; break; alpar@1: case GLP_NS: alpar@1: row->prim = row->lb; break; alpar@1: default: alpar@1: xassert(stat != stat); alpar@1: } alpar@1: #endif alpar@1: row->dual = (cbar[j] * row->rii) / zeta; alpar@1: } alpar@1: else alpar@1: { GLPCOL *col = lp->col[k-m]; alpar@1: col->stat = stat[j]; alpar@1: col->bind = 0; alpar@1: #if 0 alpar@1: col->prim = get_xN(csa, j) * col->sjj; alpar@1: #else alpar@1: switch (stat[j]) alpar@1: { case GLP_NL: alpar@1: col->prim = col->lb; break; alpar@1: case GLP_NU: alpar@1: col->prim = col->ub; break; alpar@1: case GLP_NF: alpar@1: col->prim = 0.0; break; alpar@1: case GLP_NS: alpar@1: col->prim = col->lb; break; alpar@1: default: alpar@1: xassert(stat != stat); alpar@1: } alpar@1: #endif alpar@1: col->dual = (cbar[j] / col->sjj) / zeta; alpar@1: } alpar@1: } alpar@1: return; alpar@1: } alpar@1: #endif alpar@1: alpar@1: /*********************************************************************** alpar@1: * free_csa - deallocate common storage area alpar@1: * alpar@1: * This routine frees all the memory allocated to arrays in the common alpar@1: * storage area (CSA). */ alpar@1: alpar@1: static void free_csa(struct csa *csa) alpar@1: { xfree(csa->type); alpar@1: xfree(csa->lb); alpar@1: xfree(csa->ub); alpar@1: xfree(csa->coef); alpar@1: xfree(csa->orig_type); alpar@1: xfree(csa->orig_lb); alpar@1: xfree(csa->orig_ub); alpar@1: xfree(csa->obj); alpar@1: xfree(csa->A_ptr); alpar@1: xfree(csa->A_ind); alpar@1: xfree(csa->A_val); alpar@1: #if 1 /* 06/IV-2009 */ alpar@1: xfree(csa->AT_ptr); alpar@1: xfree(csa->AT_ind); alpar@1: xfree(csa->AT_val); alpar@1: #endif alpar@1: xfree(csa->head); alpar@1: #if 1 /* 06/IV-2009 */ alpar@1: xfree(csa->bind); alpar@1: #endif alpar@1: xfree(csa->stat); alpar@1: #if 0 /* 06/IV-2009 */ alpar@1: xfree(csa->N_ptr); alpar@1: xfree(csa->N_len); alpar@1: xfree(csa->N_ind); alpar@1: xfree(csa->N_val); alpar@1: #endif alpar@1: xfree(csa->bbar); alpar@1: xfree(csa->cbar); alpar@1: xfree(csa->refsp); alpar@1: xfree(csa->gamma); alpar@1: xfree(csa->trow_ind); alpar@1: xfree(csa->trow_vec); alpar@1: #ifdef GLP_LONG_STEP /* 07/IV-2009 */ alpar@1: xfree(csa->bkpt); alpar@1: #endif alpar@1: xfree(csa->tcol_ind); alpar@1: xfree(csa->tcol_vec); alpar@1: xfree(csa->work1); alpar@1: xfree(csa->work2); alpar@1: xfree(csa->work3); alpar@1: xfree(csa->work4); alpar@1: xfree(csa); alpar@1: return; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * spx_dual - core LP solver based on the dual simplex method alpar@1: * alpar@1: * SYNOPSIS alpar@1: * alpar@1: * #include "glpspx.h" alpar@1: * int spx_dual(glp_prob *lp, const glp_smcp *parm); alpar@1: * alpar@1: * DESCRIPTION alpar@1: * alpar@1: * The routine spx_dual is a core LP solver based on the two-phase dual alpar@1: * simplex method. alpar@1: * alpar@1: * RETURNS alpar@1: * alpar@1: * 0 LP instance has been successfully solved. alpar@1: * alpar@1: * GLP_EOBJLL alpar@1: * Objective lower limit has been reached (maximization). alpar@1: * alpar@1: * GLP_EOBJUL alpar@1: * Objective upper limit has been reached (minimization). alpar@1: * alpar@1: * GLP_EITLIM alpar@1: * Iteration limit has been exhausted. alpar@1: * alpar@1: * GLP_ETMLIM alpar@1: * Time limit has been exhausted. alpar@1: * alpar@1: * GLP_EFAIL alpar@1: * The solver failed to solve LP instance. */ alpar@1: alpar@1: int spx_dual(glp_prob *lp, const glp_smcp *parm) alpar@1: { struct csa *csa; alpar@1: int binv_st = 2; alpar@1: /* status of basis matrix factorization: alpar@1: 0 - invalid; 1 - just computed; 2 - updated */ alpar@1: int bbar_st = 0; alpar@1: /* status of primal values of basic variables: alpar@1: 0 - invalid; 1 - just computed; 2 - updated */ alpar@1: int cbar_st = 0; alpar@1: /* status of reduced costs of non-basic variables: alpar@1: 0 - invalid; 1 - just computed; 2 - updated */ alpar@1: int rigorous = 0; alpar@1: /* rigorous mode flag; this flag is used to enable iterative alpar@1: refinement on computing pivot rows and columns of the simplex alpar@1: table */ alpar@1: int check = 0; alpar@1: int p_stat, d_stat, ret; alpar@1: /* allocate and initialize the common storage area */ alpar@1: csa = alloc_csa(lp); alpar@1: init_csa(csa, lp); alpar@1: if (parm->msg_lev >= GLP_MSG_DBG) alpar@1: xprintf("Objective scale factor = %g\n", csa->zeta); alpar@1: loop: /* main loop starts here */ alpar@1: /* compute factorization of the basis matrix */ alpar@1: if (binv_st == 0) alpar@1: { ret = invert_B(csa); alpar@1: if (ret != 0) alpar@1: { if (parm->msg_lev >= GLP_MSG_ERR) alpar@1: { xprintf("Error: unable to factorize the basis matrix (%d" alpar@1: ")\n", ret); alpar@1: xprintf("Sorry, basis recovery procedure not implemented" alpar@1: " yet\n"); alpar@1: } alpar@1: xassert(!lp->valid && lp->bfd == NULL); alpar@1: lp->bfd = csa->bfd, csa->bfd = NULL; alpar@1: lp->pbs_stat = lp->dbs_stat = GLP_UNDEF; alpar@1: lp->obj_val = 0.0; alpar@1: lp->it_cnt = csa->it_cnt; alpar@1: lp->some = 0; alpar@1: ret = GLP_EFAIL; alpar@1: goto done; alpar@1: } alpar@1: csa->valid = 1; alpar@1: binv_st = 1; /* just computed */ alpar@1: /* invalidate basic solution components */ alpar@1: bbar_st = cbar_st = 0; alpar@1: } alpar@1: /* compute reduced costs of non-basic variables */ alpar@1: if (cbar_st == 0) alpar@1: { eval_cbar(csa); alpar@1: cbar_st = 1; /* just computed */ alpar@1: /* determine the search phase, if not determined yet */ alpar@1: if (csa->phase == 0) alpar@1: { if (check_feas(csa, 0.90 * parm->tol_dj) != 0) alpar@1: { /* current basic solution is dual infeasible */ alpar@1: /* start searching for dual feasible solution */ alpar@1: csa->phase = 1; alpar@1: set_aux_bnds(csa); alpar@1: } alpar@1: else alpar@1: { /* current basic solution is dual feasible */ alpar@1: /* start searching for optimal solution */ alpar@1: csa->phase = 2; alpar@1: set_orig_bnds(csa); alpar@1: } alpar@1: xassert(check_stab(csa, parm->tol_dj) == 0); alpar@1: /* some non-basic double-bounded variables might become alpar@1: fixed (on phase I) or vice versa (on phase II) */ alpar@1: #if 0 /* 06/IV-2009 */ alpar@1: build_N(csa); alpar@1: #endif alpar@1: csa->refct = 0; alpar@1: /* bounds of non-basic variables have been changed, so alpar@1: invalidate primal values */ alpar@1: bbar_st = 0; alpar@1: } alpar@1: /* make sure that the current basic solution remains dual alpar@1: feasible */ alpar@1: if (check_stab(csa, parm->tol_dj) != 0) alpar@1: { if (parm->msg_lev >= GLP_MSG_ERR) alpar@1: xprintf("Warning: numerical instability (dual simplex, p" alpar@1: "hase %s)\n", csa->phase == 1 ? "I" : "II"); alpar@1: #if 1 alpar@1: if (parm->meth == GLP_DUALP) alpar@1: { store_sol(csa, lp, GLP_UNDEF, GLP_UNDEF, 0); alpar@1: ret = GLP_EFAIL; alpar@1: goto done; alpar@1: } alpar@1: #endif alpar@1: /* restart the search */ alpar@1: csa->phase = 0; alpar@1: binv_st = 0; alpar@1: rigorous = 5; alpar@1: goto loop; alpar@1: } alpar@1: } alpar@1: xassert(csa->phase == 1 || csa->phase == 2); alpar@1: /* on phase I we do not need to wait until the current basic alpar@1: solution becomes primal feasible; it is sufficient to make alpar@1: sure that all reduced costs have correct signs */ alpar@1: if (csa->phase == 1 && check_feas(csa, parm->tol_dj) == 0) alpar@1: { /* the current basis is dual feasible; switch to phase II */ alpar@1: display(csa, parm, 1); alpar@1: csa->phase = 2; alpar@1: if (cbar_st != 1) alpar@1: { eval_cbar(csa); alpar@1: cbar_st = 1; alpar@1: } alpar@1: set_orig_bnds(csa); alpar@1: #if 0 /* 06/IV-2009 */ alpar@1: build_N(csa); alpar@1: #endif alpar@1: csa->refct = 0; alpar@1: bbar_st = 0; alpar@1: } alpar@1: /* compute primal values of basic variables */ alpar@1: if (bbar_st == 0) alpar@1: { eval_bbar(csa); alpar@1: if (csa->phase == 2) alpar@1: csa->bbar[0] = eval_obj(csa); alpar@1: bbar_st = 1; /* just computed */ alpar@1: } alpar@1: /* redefine the reference space, if required */ alpar@1: switch (parm->pricing) alpar@1: { case GLP_PT_STD: alpar@1: break; alpar@1: case GLP_PT_PSE: alpar@1: if (csa->refct == 0) reset_refsp(csa); alpar@1: break; alpar@1: default: alpar@1: xassert(parm != parm); alpar@1: } alpar@1: /* at this point the basis factorization and all basic solution alpar@1: components are valid */ alpar@1: xassert(binv_st && bbar_st && cbar_st); alpar@1: /* check accuracy of current basic solution components (only for alpar@1: debugging) */ alpar@1: if (check) alpar@1: { double e_bbar = err_in_bbar(csa); alpar@1: double e_cbar = err_in_cbar(csa); alpar@1: double e_gamma = alpar@1: (parm->pricing == GLP_PT_PSE ? err_in_gamma(csa) : 0.0); alpar@1: xprintf("e_bbar = %10.3e; e_cbar = %10.3e; e_gamma = %10.3e\n", alpar@1: e_bbar, e_cbar, e_gamma); alpar@1: xassert(e_bbar <= 1e-5 && e_cbar <= 1e-5 && e_gamma <= 1e-3); alpar@1: } alpar@1: /* if the objective has to be maximized, check if it has reached alpar@1: its lower limit */ alpar@1: if (csa->phase == 2 && csa->zeta < 0.0 && alpar@1: parm->obj_ll > -DBL_MAX && csa->bbar[0] <= parm->obj_ll) alpar@1: { if (bbar_st != 1 || cbar_st != 1) alpar@1: { if (bbar_st != 1) bbar_st = 0; alpar@1: if (cbar_st != 1) cbar_st = 0; alpar@1: goto loop; alpar@1: } alpar@1: display(csa, parm, 1); alpar@1: if (parm->msg_lev >= GLP_MSG_ALL) alpar@1: xprintf("OBJECTIVE LOWER LIMIT REACHED; SEARCH TERMINATED\n" alpar@1: ); alpar@1: store_sol(csa, lp, GLP_INFEAS, GLP_FEAS, 0); alpar@1: ret = GLP_EOBJLL; alpar@1: goto done; alpar@1: } alpar@1: /* if the objective has to be minimized, check if it has reached alpar@1: its upper limit */ alpar@1: if (csa->phase == 2 && csa->zeta > 0.0 && alpar@1: parm->obj_ul < +DBL_MAX && csa->bbar[0] >= parm->obj_ul) alpar@1: { if (bbar_st != 1 || cbar_st != 1) alpar@1: { if (bbar_st != 1) bbar_st = 0; alpar@1: if (cbar_st != 1) cbar_st = 0; alpar@1: goto loop; alpar@1: } alpar@1: display(csa, parm, 1); alpar@1: if (parm->msg_lev >= GLP_MSG_ALL) alpar@1: xprintf("OBJECTIVE UPPER LIMIT REACHED; SEARCH TERMINATED\n" alpar@1: ); alpar@1: store_sol(csa, lp, GLP_INFEAS, GLP_FEAS, 0); alpar@1: ret = GLP_EOBJUL; alpar@1: goto done; alpar@1: } alpar@1: /* check if the iteration limit has been exhausted */ alpar@1: if (parm->it_lim < INT_MAX && alpar@1: csa->it_cnt - csa->it_beg >= parm->it_lim) alpar@1: { if (csa->phase == 2 && bbar_st != 1 || cbar_st != 1) alpar@1: { if (csa->phase == 2 && bbar_st != 1) bbar_st = 0; alpar@1: if (cbar_st != 1) cbar_st = 0; alpar@1: goto loop; alpar@1: } alpar@1: display(csa, parm, 1); alpar@1: if (parm->msg_lev >= GLP_MSG_ALL) alpar@1: xprintf("ITERATION LIMIT EXCEEDED; SEARCH TERMINATED\n"); alpar@1: switch (csa->phase) alpar@1: { case 1: alpar@1: d_stat = GLP_INFEAS; alpar@1: set_orig_bnds(csa); alpar@1: eval_bbar(csa); alpar@1: break; alpar@1: case 2: alpar@1: d_stat = GLP_FEAS; alpar@1: break; alpar@1: default: alpar@1: xassert(csa != csa); alpar@1: } alpar@1: store_sol(csa, lp, GLP_INFEAS, d_stat, 0); alpar@1: ret = GLP_EITLIM; alpar@1: goto done; alpar@1: } alpar@1: /* check if the time limit has been exhausted */ alpar@1: if (parm->tm_lim < INT_MAX && alpar@1: 1000.0 * xdifftime(xtime(), csa->tm_beg) >= parm->tm_lim) alpar@1: { if (csa->phase == 2 && bbar_st != 1 || cbar_st != 1) alpar@1: { if (csa->phase == 2 && bbar_st != 1) bbar_st = 0; alpar@1: if (cbar_st != 1) cbar_st = 0; alpar@1: goto loop; alpar@1: } alpar@1: display(csa, parm, 1); alpar@1: if (parm->msg_lev >= GLP_MSG_ALL) alpar@1: xprintf("TIME LIMIT EXCEEDED; SEARCH TERMINATED\n"); alpar@1: switch (csa->phase) alpar@1: { case 1: alpar@1: d_stat = GLP_INFEAS; alpar@1: set_orig_bnds(csa); alpar@1: eval_bbar(csa); alpar@1: break; alpar@1: case 2: alpar@1: d_stat = GLP_FEAS; alpar@1: break; alpar@1: default: alpar@1: xassert(csa != csa); alpar@1: } alpar@1: store_sol(csa, lp, GLP_INFEAS, d_stat, 0); alpar@1: ret = GLP_ETMLIM; alpar@1: goto done; alpar@1: } alpar@1: /* display the search progress */ alpar@1: display(csa, parm, 0); alpar@1: /* choose basic variable xB[p] */ alpar@1: chuzr(csa, parm->tol_bnd); alpar@1: if (csa->p == 0) alpar@1: { if (bbar_st != 1 || cbar_st != 1) alpar@1: { if (bbar_st != 1) bbar_st = 0; alpar@1: if (cbar_st != 1) cbar_st = 0; alpar@1: goto loop; alpar@1: } alpar@1: display(csa, parm, 1); alpar@1: switch (csa->phase) alpar@1: { case 1: alpar@1: if (parm->msg_lev >= GLP_MSG_ALL) alpar@1: xprintf("PROBLEM HAS NO DUAL FEASIBLE SOLUTION\n"); alpar@1: set_orig_bnds(csa); alpar@1: eval_bbar(csa); alpar@1: p_stat = GLP_INFEAS, d_stat = GLP_NOFEAS; alpar@1: break; alpar@1: case 2: alpar@1: if (parm->msg_lev >= GLP_MSG_ALL) alpar@1: xprintf("OPTIMAL SOLUTION FOUND\n"); alpar@1: p_stat = d_stat = GLP_FEAS; alpar@1: break; alpar@1: default: alpar@1: xassert(csa != csa); alpar@1: } alpar@1: store_sol(csa, lp, p_stat, d_stat, 0); alpar@1: ret = 0; alpar@1: goto done; alpar@1: } alpar@1: /* compute pivot row of the simplex table */ alpar@1: { double *rho = csa->work4; alpar@1: eval_rho(csa, rho); alpar@1: if (rigorous) refine_rho(csa, rho); alpar@1: eval_trow(csa, rho); alpar@1: sort_trow(csa, parm->tol_bnd); alpar@1: } alpar@1: /* unlike primal simplex there is no need to check accuracy of alpar@1: the primal value of xB[p] (which might be computed using the alpar@1: pivot row), since bbar is a result of FTRAN */ alpar@1: #ifdef GLP_LONG_STEP /* 07/IV-2009 */ alpar@1: long_step(csa); alpar@1: if (csa->nbps > 0) alpar@1: { csa->q = csa->bkpt[csa->nbps].j; alpar@1: if (csa->delta > 0.0) alpar@1: csa->new_dq = + csa->bkpt[csa->nbps].t; alpar@1: else alpar@1: csa->new_dq = - csa->bkpt[csa->nbps].t; alpar@1: } alpar@1: else alpar@1: #endif alpar@1: /* choose non-basic variable xN[q] */ alpar@1: switch (parm->r_test) alpar@1: { case GLP_RT_STD: alpar@1: chuzc(csa, 0.0); alpar@1: break; alpar@1: case GLP_RT_HAR: alpar@1: chuzc(csa, 0.30 * parm->tol_dj); alpar@1: break; alpar@1: default: alpar@1: xassert(parm != parm); alpar@1: } alpar@1: if (csa->q == 0) alpar@1: { if (bbar_st != 1 || cbar_st != 1 || !rigorous) alpar@1: { if (bbar_st != 1) bbar_st = 0; alpar@1: if (cbar_st != 1) cbar_st = 0; alpar@1: rigorous = 1; alpar@1: goto loop; alpar@1: } alpar@1: display(csa, parm, 1); alpar@1: switch (csa->phase) alpar@1: { case 1: alpar@1: if (parm->msg_lev >= GLP_MSG_ERR) alpar@1: xprintf("Error: unable to choose basic variable on ph" alpar@1: "ase I\n"); alpar@1: xassert(!lp->valid && lp->bfd == NULL); alpar@1: lp->bfd = csa->bfd, csa->bfd = NULL; alpar@1: lp->pbs_stat = lp->dbs_stat = GLP_UNDEF; alpar@1: lp->obj_val = 0.0; alpar@1: lp->it_cnt = csa->it_cnt; alpar@1: lp->some = 0; alpar@1: ret = GLP_EFAIL; alpar@1: break; alpar@1: case 2: alpar@1: if (parm->msg_lev >= GLP_MSG_ALL) alpar@1: xprintf("PROBLEM HAS NO FEASIBLE SOLUTION\n"); alpar@1: store_sol(csa, lp, GLP_NOFEAS, GLP_FEAS, alpar@1: csa->head[csa->p]); alpar@1: ret = 0; alpar@1: break; alpar@1: default: alpar@1: xassert(csa != csa); alpar@1: } alpar@1: goto done; alpar@1: } alpar@1: /* check if the pivot element is acceptable */ alpar@1: { double piv = csa->trow_vec[csa->q]; alpar@1: double eps = 1e-5 * (1.0 + 0.01 * csa->trow_max); alpar@1: if (fabs(piv) < eps) alpar@1: { if (parm->msg_lev >= GLP_MSG_DBG) alpar@1: xprintf("piv = %.12g; eps = %g\n", piv, eps); alpar@1: if (!rigorous) alpar@1: { rigorous = 5; alpar@1: goto loop; alpar@1: } alpar@1: } alpar@1: } alpar@1: /* now xN[q] and xB[p] have been chosen anyhow */ alpar@1: /* compute pivot column of the simplex table */ alpar@1: eval_tcol(csa); alpar@1: if (rigorous) refine_tcol(csa); alpar@1: /* accuracy check based on the pivot element */ alpar@1: { double piv1 = csa->tcol_vec[csa->p]; /* more accurate */ alpar@1: double piv2 = csa->trow_vec[csa->q]; /* less accurate */ alpar@1: xassert(piv1 != 0.0); alpar@1: if (fabs(piv1 - piv2) > 1e-8 * (1.0 + fabs(piv1)) || alpar@1: !(piv1 > 0.0 && piv2 > 0.0 || piv1 < 0.0 && piv2 < 0.0)) alpar@1: { if (parm->msg_lev >= GLP_MSG_DBG) alpar@1: xprintf("piv1 = %.12g; piv2 = %.12g\n", piv1, piv2); alpar@1: if (binv_st != 1 || !rigorous) alpar@1: { if (binv_st != 1) binv_st = 0; alpar@1: rigorous = 5; alpar@1: goto loop; alpar@1: } alpar@1: /* (not a good idea; should be revised later) */ alpar@1: if (csa->tcol_vec[csa->p] == 0.0) alpar@1: { csa->tcol_nnz++; alpar@1: xassert(csa->tcol_nnz <= csa->m); alpar@1: csa->tcol_ind[csa->tcol_nnz] = csa->p; alpar@1: } alpar@1: csa->tcol_vec[csa->p] = piv2; alpar@1: } alpar@1: } alpar@1: /* update primal values of basic variables */ alpar@1: #ifdef GLP_LONG_STEP /* 07/IV-2009 */ alpar@1: if (csa->nbps > 0) alpar@1: { int kk, j, k; alpar@1: for (kk = 1; kk < csa->nbps; kk++) alpar@1: { if (csa->bkpt[kk].t >= csa->bkpt[csa->nbps].t) continue; alpar@1: j = csa->bkpt[kk].j; alpar@1: k = csa->head[csa->m + j]; alpar@1: xassert(csa->type[k] == GLP_DB); alpar@1: if (csa->stat[j] == GLP_NL) alpar@1: csa->stat[j] = GLP_NU; alpar@1: else alpar@1: csa->stat[j] = GLP_NL; alpar@1: } alpar@1: } alpar@1: bbar_st = 0; alpar@1: #else alpar@1: update_bbar(csa); alpar@1: if (csa->phase == 2) alpar@1: csa->bbar[0] += (csa->cbar[csa->q] / csa->zeta) * alpar@1: (csa->delta / csa->tcol_vec[csa->p]); alpar@1: bbar_st = 2; /* updated */ alpar@1: #endif alpar@1: /* update reduced costs of non-basic variables */ alpar@1: update_cbar(csa); alpar@1: cbar_st = 2; /* updated */ alpar@1: /* update steepest edge coefficients */ alpar@1: switch (parm->pricing) alpar@1: { case GLP_PT_STD: alpar@1: break; alpar@1: case GLP_PT_PSE: alpar@1: if (csa->refct > 0) update_gamma(csa); alpar@1: break; alpar@1: default: alpar@1: xassert(parm != parm); alpar@1: } alpar@1: /* update factorization of the basis matrix */ alpar@1: ret = update_B(csa, csa->p, csa->head[csa->m+csa->q]); alpar@1: if (ret == 0) alpar@1: binv_st = 2; /* updated */ alpar@1: else alpar@1: { csa->valid = 0; alpar@1: binv_st = 0; /* invalid */ alpar@1: } alpar@1: #if 0 /* 06/IV-2009 */ alpar@1: /* update matrix N */ alpar@1: del_N_col(csa, csa->q, csa->head[csa->m+csa->q]); alpar@1: if (csa->type[csa->head[csa->p]] != GLP_FX) alpar@1: add_N_col(csa, csa->q, csa->head[csa->p]); alpar@1: #endif alpar@1: /* change the basis header */ alpar@1: change_basis(csa); alpar@1: /* iteration complete */ alpar@1: csa->it_cnt++; alpar@1: if (rigorous > 0) rigorous--; alpar@1: goto loop; alpar@1: done: /* deallocate the common storage area */ alpar@1: free_csa(csa); alpar@1: /* return to the calling program */ alpar@1: return ret; alpar@1: } alpar@1: alpar@1: /* eof */