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