alpar@9: /* glpspx01.c (primal 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: 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] (note that on phase I auxiliary variables also alpar@9: may have non-zero objective coefficients) */ 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: /*--------------------------------------------------------------*/ 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: 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: /*--------------------------------------------------------------*/ 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: /*--------------------------------------------------------------*/ alpar@9: /* working parameters */ alpar@9: int phase; alpar@9: /* search phase: alpar@9: 0 - not determined yet alpar@9: 1 - search for primal 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 (including the case when a non-basic variable alpar@9: jumps to its opposite bound) */ 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; 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+n]; */ alpar@9: /* gamma[0] is not used; alpar@9: gamma[j], 1 <= j <= n, is the steepest edge coefficient for alpar@9: non-basic variable xN[j]; if xN[j] is fixed, gamma[j] is not alpar@9: used and just set to 1 */ 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 the set of eligible non-basic variables is empty and thus alpar@9: no variable has been chosen, q is set to 0 */ 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: double tcol_max; alpar@9: /* infinity (maximum) norm of the column (max |tcol_vec[i]|) */ alpar@9: int tcol_num; alpar@9: /* number of significant non-zero components, which means that: alpar@9: |tcol_vec[i]| >= eps for i in tcol_ind[1,...,num], alpar@9: |tcol_vec[i]| < eps for i in tcol_ind[num+1,...,nnz], alpar@9: where eps is a pivot tolerance */ 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: p = 0 means that no basic variable reaches its bound; alpar@9: p < 0 means that non-basic variable xN[q] reaches its opposite alpar@9: bound before any basic variable */ alpar@9: int p_stat; alpar@9: /* new status (GLP_NL, GLP_NU, or GLP_NS) to be assigned to xB[p] alpar@9: once it has left the basis */ alpar@9: double teta; alpar@9: /* change of non-basic variable xN[q] (see above), on which xB[p] alpar@9: (or, if p < 0, xN[q] itself) reaches its bound */ 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: /*--------------------------------------------------------------*/ 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->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: csa->head = xcalloc(1+m+n, sizeof(int)); alpar@9: csa->stat = xcalloc(1+n, sizeof(char)); 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: 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+n, sizeof(double)); alpar@9: csa->tcol_ind = xcalloc(1+m, sizeof(int)); alpar@9: csa->tcol_vec = 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: 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 alloc_N(struct csa *csa); alpar@9: static void build_N(struct csa *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: 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: int *head = csa->head; 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 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: /* 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 == lp->nnz+1); 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: /* factorization of matrix B */ alpar@9: csa->valid = 1, lp->valid = 0; alpar@9: csa->bfd = lp->bfd, lp->bfd = NULL; alpar@9: /* matrix N (by rows) */ alpar@9: alloc_N(csa); alpar@9: build_N(csa); 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 (j = 1; j <= n; j++) gamma[j] = 1.0; alpar@9: return; alpar@9: } alpar@9: 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: 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: 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: 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: 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: 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: alpar@9: /*********************************************************************** alpar@9: * alloc_N - allocate matrix N alpar@9: * alpar@9: * This routine determines maximal row lengths of matrix N, sets its alpar@9: * row pointers, and then allocates arrays N_ind and N_val. alpar@9: * alpar@9: * Note that some fixed structural variables may temporarily become alpar@9: * double-bounded, so corresponding columns of matrix A should not be alpar@9: * ignored on calculating maximal row lengths of matrix N. */ alpar@9: alpar@9: static void alloc_N(struct csa *csa) 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: int *N_ptr = csa->N_ptr; alpar@9: int *N_len = csa->N_len; alpar@9: int i, j, beg, end, ptr; alpar@9: /* determine number of non-zeros in each row of the augmented alpar@9: constraint matrix (I|-A) */ alpar@9: for (i = 1; i <= m; i++) alpar@9: N_len[i] = 1; alpar@9: for (j = 1; j <= n; j++) alpar@9: { beg = A_ptr[j]; alpar@9: end = A_ptr[j+1]; alpar@9: for (ptr = beg; ptr < end; ptr++) alpar@9: N_len[A_ind[ptr]]++; alpar@9: } alpar@9: /* determine maximal row lengths of matrix N and set its row alpar@9: pointers */ alpar@9: N_ptr[1] = 1; alpar@9: for (i = 1; i <= m; i++) alpar@9: { /* row of matrix N cannot have more than n non-zeros */ alpar@9: if (N_len[i] > n) N_len[i] = n; alpar@9: N_ptr[i+1] = N_ptr[i] + N_len[i]; alpar@9: } alpar@9: /* now maximal number of non-zeros in matrix N is known */ alpar@9: csa->N_ind = xcalloc(N_ptr[m+1], sizeof(int)); alpar@9: csa->N_val = xcalloc(N_ptr[m+1], sizeof(double)); alpar@9: return; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * add_N_col - add column of matrix (I|-A) to matrix N alpar@9: * alpar@9: * This routine adds j-th column to matrix N which is k-th column of alpar@9: * the augmented constraint matrix (I|-A). (It is assumed that old j-th alpar@9: * column was previously removed from matrix N.) */ alpar@9: alpar@9: static void add_N_col(struct csa *csa, int j, 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 *N_ptr = csa->N_ptr; alpar@9: int *N_len = csa->N_len; alpar@9: int *N_ind = csa->N_ind; alpar@9: double *N_val = csa->N_val; alpar@9: int pos; alpar@9: #ifdef GLP_DEBUG alpar@9: xassert(1 <= j && j <= n); alpar@9: xassert(1 <= k && k <= m+n); alpar@9: #endif alpar@9: if (k <= m) alpar@9: { /* N[j] is k-th column of submatrix I */ alpar@9: pos = N_ptr[k] + (N_len[k]++); alpar@9: #ifdef GLP_DEBUG alpar@9: xassert(pos < N_ptr[k+1]); alpar@9: #endif alpar@9: N_ind[pos] = j; alpar@9: N_val[pos] = 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 i, 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: { i = A_ind[ptr]; /* row number */ alpar@9: pos = N_ptr[i] + (N_len[i]++); alpar@9: #ifdef GLP_DEBUG alpar@9: xassert(pos < N_ptr[i+1]); alpar@9: #endif alpar@9: N_ind[pos] = j; alpar@9: N_val[pos] = - A_val[ptr]; alpar@9: } alpar@9: } alpar@9: return; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * del_N_col - remove column of matrix (I|-A) from matrix N alpar@9: * alpar@9: * This routine removes j-th column from matrix N which is k-th column alpar@9: * of the augmented constraint matrix (I|-A). */ alpar@9: alpar@9: static void del_N_col(struct csa *csa, int j, 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 *N_ptr = csa->N_ptr; alpar@9: int *N_len = csa->N_len; alpar@9: int *N_ind = csa->N_ind; alpar@9: double *N_val = csa->N_val; alpar@9: int pos, head, tail; alpar@9: #ifdef GLP_DEBUG alpar@9: xassert(1 <= j && j <= n); alpar@9: xassert(1 <= k && k <= m+n); alpar@9: #endif alpar@9: if (k <= m) alpar@9: { /* N[j] is k-th column of submatrix I */ alpar@9: /* find element in k-th row of N */ alpar@9: head = N_ptr[k]; alpar@9: for (pos = head; N_ind[pos] != j; pos++) /* nop */; alpar@9: /* and remove it from the row list */ alpar@9: tail = head + (--N_len[k]); alpar@9: #ifdef GLP_DEBUG alpar@9: xassert(pos <= tail); alpar@9: #endif alpar@9: N_ind[pos] = N_ind[tail]; alpar@9: N_val[pos] = N_val[tail]; 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: int i, 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: { i = A_ind[ptr]; /* row number */ alpar@9: /* find element in i-th row of N */ alpar@9: head = N_ptr[i]; alpar@9: for (pos = head; N_ind[pos] != j; pos++) /* nop */; alpar@9: /* and remove it from the row list */ alpar@9: tail = head + (--N_len[i]); alpar@9: #ifdef GLP_DEBUG alpar@9: xassert(pos <= tail); alpar@9: #endif alpar@9: N_ind[pos] = N_ind[tail]; alpar@9: N_val[pos] = N_val[tail]; alpar@9: } alpar@9: } alpar@9: return; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * build_N - build matrix N for current basis alpar@9: * alpar@9: * This routine builds matrix N for the current basis from columns alpar@9: * of the augmented constraint matrix (I|-A) corresponding to non-basic alpar@9: * non-fixed variables. */ alpar@9: alpar@9: static void build_N(struct csa *csa) alpar@9: { int m = csa->m; alpar@9: int n = csa->n; alpar@9: int *head = csa->head; alpar@9: char *stat = csa->stat; alpar@9: int *N_len = csa->N_len; alpar@9: int j, k; alpar@9: /* N := empty matrix */ alpar@9: memset(&N_len[1], 0, m * sizeof(int)); alpar@9: /* go through non-basic columns of matrix (I|-A) */ alpar@9: for (j = 1; j <= n; j++) alpar@9: { if (stat[j] != GLP_NS) alpar@9: { /* xN[j] is non-fixed; add j-th column to matrix N which is alpar@9: k-th column of matrix (I|-A) */ 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: add_N_col(csa, j, k); alpar@9: } alpar@9: } alpar@9: return; alpar@9: } alpar@9: 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: 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: 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: 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: 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: 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: 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 j, 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 (j = 1; j <= n; j++) alpar@9: { k = head[m+j]; /* x[k] = xN[j] */ alpar@9: refsp[k] = 1; alpar@9: gamma[j] = 1.0; alpar@9: } alpar@9: return; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * eval_gamma - compute steepest edge coefficient alpar@9: * alpar@9: * This routine computes the steepest edge coefficient for non-basic alpar@9: * variable xN[j] using its direct definition: alpar@9: * alpar@9: * gamma[j] = delta[j] + sum alfa[i,j]^2, alpar@9: * i in R alpar@9: * alpar@9: * where delta[j] = 1, if xN[j] is in the current reference space, alpar@9: * and 0 otherwise; R is a set of basic variables xB[i], which are in alpar@9: * the current reference space; alfa[i,j] are elements of the current alpar@9: * simplex table. alpar@9: * alpar@9: * NOTE: The routine is intended only for debugginig purposes. */ alpar@9: alpar@9: static double eval_gamma(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: 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, k; alpar@9: double gamma; 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: /* 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: /* compute gamma */ alpar@9: gamma = (refsp[k] ? 1.0 : 0.0); alpar@9: for (i = 1; i <= m; i++) alpar@9: { k = head[i]; alpar@9: #ifdef GLP_DEBUG alpar@9: xassert(1 <= k && k <= m+n); alpar@9: #endif alpar@9: if (refsp[k]) gamma += alfa[i] * alfa[i]; alpar@9: } alpar@9: return gamma; alpar@9: } 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 has largest alpar@9: * weighted reduced cost: alpar@9: * alpar@9: * |d[q]| / sqrt(gamma[q]) = max |d[j]| / sqrt(gamma[j]), alpar@9: * j in J alpar@9: * alpar@9: * where J is a subset of eligible non-basic variables xN[j], d[j] is alpar@9: * reduced cost of xN[j], gamma[j] is the steepest edge coefficient. alpar@9: * alpar@9: * The working objective function is always minimized, so the sign of alpar@9: * d[q] determines direction, in which xN[q] has to change: alpar@9: * alpar@9: * if d[q] < 0, xN[q] has to increase; alpar@9: * alpar@9: * if d[q] > 0, xN[q] has to decrease. alpar@9: * alpar@9: * If |d[j]| <= tol_dj, where tol_dj is a specified tolerance, xN[j] alpar@9: * is not included in J and therefore ignored. (It is assumed that the alpar@9: * working objective row is appropriately scaled, i.e. max|c[k]| = 1.) alpar@9: * alpar@9: * If J is empty and no variable has been chosen, q is set to 0. */ alpar@9: alpar@9: static void chuzc(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: double *gamma = csa->gamma; alpar@9: int j, q; alpar@9: double dj, best, temp; alpar@9: /* nothing is chosen so far */ alpar@9: q = 0, best = 0.0; alpar@9: /* look through the list of non-basic variables */ alpar@9: for (j = 1; j <= n; j++) alpar@9: { dj = cbar[j]; alpar@9: switch (stat[j]) alpar@9: { case GLP_NL: alpar@9: /* xN[j] can increase */ alpar@9: if (dj >= - tol_dj) continue; alpar@9: break; alpar@9: case GLP_NU: alpar@9: /* xN[j] can decrease */ alpar@9: if (dj <= + tol_dj) continue; alpar@9: break; alpar@9: case GLP_NF: alpar@9: /* xN[j] can change in any direction */ alpar@9: if (- tol_dj <= dj && dj <= + tol_dj) continue; alpar@9: break; alpar@9: case GLP_NS: alpar@9: /* xN[j] cannot change at all */ alpar@9: continue; alpar@9: default: alpar@9: xassert(stat != stat); alpar@9: } alpar@9: /* xN[j] is eligible non-basic variable; choose one which has alpar@9: largest weighted reduced cost */ alpar@9: #ifdef GLP_DEBUG alpar@9: xassert(gamma[j] > 0.0); alpar@9: #endif alpar@9: temp = (dj * dj) / gamma[j]; alpar@9: if (best < temp) alpar@9: q = j, best = temp; alpar@9: } alpar@9: /* store the index of non-basic variable xN[q] chosen */ alpar@9: csa->q = q; alpar@9: return; alpar@9: } alpar@9: 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: 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: alpar@9: /*********************************************************************** alpar@9: * sort_tcol - sort pivot column of the simplex table alpar@9: * alpar@9: * This routine reorders the list of non-zero elements of the pivot alpar@9: * column 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 tcol_num. */ alpar@9: alpar@9: static void sort_tcol(struct csa *csa, double tol_piv) alpar@9: { alpar@9: #ifdef GLP_DEBUG alpar@9: int m = csa->m; alpar@9: #endif alpar@9: int nnz = csa->tcol_nnz; alpar@9: int *tcol_ind = csa->tcol_ind; alpar@9: double *tcol_vec = csa->tcol_vec; alpar@9: int i, num, pos; alpar@9: double big, eps, temp; alpar@9: /* compute infinity (maximum) norm of the column */ alpar@9: big = 0.0; alpar@9: for (pos = 1; pos <= nnz; pos++) alpar@9: { alpar@9: #ifdef GLP_DEBUG alpar@9: i = tcol_ind[pos]; alpar@9: xassert(1 <= i && i <= m); alpar@9: #endif alpar@9: temp = fabs(tcol_vec[tcol_ind[pos]]); alpar@9: if (big < temp) big = temp; alpar@9: } alpar@9: csa->tcol_max = big; alpar@9: /* determine absolute pivot tolerance */ alpar@9: eps = tol_piv * (1.0 + 0.01 * big); alpar@9: /* move significant column components to front of the list */ alpar@9: for (num = 0; num < nnz; ) alpar@9: { i = tcol_ind[nnz]; alpar@9: if (fabs(tcol_vec[i]) < eps) alpar@9: nnz--; alpar@9: else alpar@9: { num++; alpar@9: tcol_ind[nnz] = tcol_ind[num]; alpar@9: tcol_ind[num] = i; alpar@9: } alpar@9: } alpar@9: csa->tcol_num = num; 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], which reaches its bound alpar@9: * first on changing non-basic variable xN[q] in valid direction. alpar@9: * alpar@9: * The parameter rtol is a relative tolerance used to relax bounds of alpar@9: * basic variables. If rtol = 0, the routine implements the standard alpar@9: * ratio test. Otherwise, if rtol > 0, the routine implements Harris' alpar@9: * two-pass ratio test. In the latter case rtol should be about three alpar@9: * times less than a tolerance used to check primal feasibility. */ alpar@9: alpar@9: static void chuzr(struct csa *csa, double rtol) 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: double *coef = csa->coef; alpar@9: int *head = csa->head; alpar@9: int phase = csa->phase; alpar@9: double *bbar = csa->bbar; alpar@9: double *cbar = csa->cbar; alpar@9: int q = csa->q; alpar@9: int *tcol_ind = csa->tcol_ind; alpar@9: double *tcol_vec = csa->tcol_vec; alpar@9: int tcol_num = csa->tcol_num; alpar@9: int i, i_stat, k, p, p_stat, pos; alpar@9: double alfa, big, delta, s, t, teta, tmax; alpar@9: #ifdef GLP_DEBUG alpar@9: xassert(1 <= q && q <= n); alpar@9: #endif alpar@9: /* s := - sign(d[q]), where d[q] is reduced cost of xN[q] */ alpar@9: #ifdef GLP_DEBUG alpar@9: xassert(cbar[q] != 0.0); alpar@9: #endif alpar@9: s = (cbar[q] > 0.0 ? -1.0 : +1.0); alpar@9: /*** FIRST PASS ***/ 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: if (type[k] == GLP_DB) alpar@9: { /* xN[q] has both lower and upper bounds */ alpar@9: p = -1, p_stat = 0, teta = ub[k] - lb[k], big = 1.0; alpar@9: } alpar@9: else alpar@9: { /* xN[q] has no opposite bound */ alpar@9: p = 0, p_stat = 0, teta = DBL_MAX, big = 0.0; alpar@9: } alpar@9: /* walk through significant elements of the pivot column */ alpar@9: for (pos = 1; pos <= tcol_num; 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]; /* x[k] = xB[i] */ alpar@9: #ifdef GLP_DEBUG alpar@9: xassert(1 <= k && k <= m+n); alpar@9: #endif alpar@9: alfa = s * tcol_vec[i]; alpar@9: #ifdef GLP_DEBUG alpar@9: xassert(alfa != 0.0); alpar@9: #endif alpar@9: /* xB[i] = ... + alfa * xN[q] + ..., and due to s we need to alpar@9: consider the only case when xN[q] is increasing */ alpar@9: if (alfa > 0.0) alpar@9: { /* xB[i] is increasing */ alpar@9: if (phase == 1 && coef[k] < 0.0) alpar@9: { /* xB[i] violates its lower bound, which plays the role alpar@9: of an upper bound on phase I */ alpar@9: delta = rtol * (1.0 + kappa * fabs(lb[k])); alpar@9: t = ((lb[k] + delta) - bbar[i]) / alfa; alpar@9: i_stat = GLP_NL; alpar@9: } alpar@9: else if (phase == 1 && coef[k] > 0.0) alpar@9: { /* xB[i] violates its upper bound, which plays the role alpar@9: of an lower bound on phase I */ alpar@9: continue; alpar@9: } alpar@9: else if (type[k] == GLP_UP || type[k] == GLP_DB || alpar@9: type[k] == GLP_FX) alpar@9: { /* xB[i] is within its bounds and has an upper bound */ alpar@9: delta = rtol * (1.0 + kappa * fabs(ub[k])); alpar@9: t = ((ub[k] + delta) - bbar[i]) / alfa; alpar@9: i_stat = GLP_NU; alpar@9: } alpar@9: else alpar@9: { /* xB[i] is within its bounds and has no upper bound */ alpar@9: continue; alpar@9: } alpar@9: } alpar@9: else alpar@9: { /* xB[i] is decreasing */ alpar@9: if (phase == 1 && coef[k] > 0.0) alpar@9: { /* xB[i] violates its upper bound, which plays the role alpar@9: of an lower bound on phase I */ alpar@9: delta = rtol * (1.0 + kappa * fabs(ub[k])); alpar@9: t = ((ub[k] - delta) - bbar[i]) / alfa; alpar@9: i_stat = GLP_NU; alpar@9: } alpar@9: else if (phase == 1 && coef[k] < 0.0) alpar@9: { /* xB[i] violates its lower bound, which plays the role alpar@9: of an upper bound on phase I */ alpar@9: continue; alpar@9: } alpar@9: else if (type[k] == GLP_LO || type[k] == GLP_DB || alpar@9: type[k] == GLP_FX) alpar@9: { /* xB[i] is within its bounds and has an lower bound */ alpar@9: delta = rtol * (1.0 + kappa * fabs(lb[k])); alpar@9: t = ((lb[k] - delta) - bbar[i]) / alfa; alpar@9: i_stat = GLP_NL; alpar@9: } alpar@9: else alpar@9: { /* xB[i] is within its bounds and has no lower bound */ alpar@9: continue; alpar@9: } alpar@9: } alpar@9: /* t is a change of xN[q], on which xB[i] reaches its bound alpar@9: (possibly relaxed); since the basic solution is assumed to alpar@9: be primal feasible (or pseudo feasible on phase I), t has alpar@9: to be non-negative by definition; however, it may happen alpar@9: that xB[i] slightly (i.e. within a tolerance) violates its alpar@9: bound, that leads to negative t; in the latter case, if alpar@9: xB[i] is chosen, negative t means that xN[q] changes in alpar@9: wrong direction; if pivot alfa[i,q] is close to zero, even alpar@9: small bound violation of xB[i] may lead to a large change alpar@9: of xN[q] in wrong direction; let, for example, xB[i] >= 0 alpar@9: and in the current basis its value be -5e-9; let also xN[q] alpar@9: be on its zero bound and should increase; from the ratio alpar@9: test rule it follows that the pivot alfa[i,q] < 0; however, alpar@9: if alfa[i,q] is, say, -1e-9, the change of xN[q] in wrong alpar@9: direction is 5e-9 / (-1e-9) = -5, and using it for updating alpar@9: values of other basic variables will give absolutely wrong alpar@9: results; therefore, if t is negative, we should replace it alpar@9: by exact zero assuming that xB[i] is exactly on its bound, alpar@9: and the 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: p = i, p_stat = i_stat, 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 xN[q] reaches its opposite bound or if no basic variable alpar@9: has been chosen on the first pass */ alpar@9: if (p <= 0) goto done; alpar@9: /* if xB[p] is a blocking variable, i.e. if it prevents xN[q] alpar@9: from any change */ alpar@9: if (teta == 0.0) goto done; alpar@9: /*** SECOND PASS ***/ alpar@9: /* here tmax is a maximal change of xN[q], on which the solution alpar@9: remains primal feasible (or pseudo feasible on phase I) within alpar@9: 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: p = 0, p_stat = 0, teta = DBL_MAX, big = 0.0; alpar@9: /* walk through significant elements of the pivot column */ alpar@9: for (pos = 1; pos <= tcol_num; 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]; /* x[k] = xB[i] */ alpar@9: #ifdef GLP_DEBUG alpar@9: xassert(1 <= k && k <= m+n); alpar@9: #endif alpar@9: alfa = s * tcol_vec[i]; alpar@9: #ifdef GLP_DEBUG alpar@9: xassert(alfa != 0.0); alpar@9: #endif alpar@9: /* xB[i] = ... + alfa * xN[q] + ..., and due to s we need to alpar@9: consider the only case when xN[q] is increasing */ alpar@9: if (alfa > 0.0) alpar@9: { /* xB[i] is increasing */ alpar@9: if (phase == 1 && coef[k] < 0.0) alpar@9: { /* xB[i] violates its lower bound, which plays the role alpar@9: of an upper bound on phase I */ alpar@9: t = (lb[k] - bbar[i]) / alfa; alpar@9: i_stat = GLP_NL; alpar@9: } alpar@9: else if (phase == 1 && coef[k] > 0.0) alpar@9: { /* xB[i] violates its upper bound, which plays the role alpar@9: of an lower bound on phase I */ alpar@9: continue; alpar@9: } alpar@9: else if (type[k] == GLP_UP || type[k] == GLP_DB || alpar@9: type[k] == GLP_FX) alpar@9: { /* xB[i] is within its bounds and has an upper bound */ alpar@9: t = (ub[k] - bbar[i]) / alfa; alpar@9: i_stat = GLP_NU; alpar@9: } alpar@9: else alpar@9: { /* xB[i] is within its bounds and has no upper bound */ alpar@9: continue; alpar@9: } alpar@9: } alpar@9: else alpar@9: { /* xB[i] is decreasing */ alpar@9: if (phase == 1 && coef[k] > 0.0) alpar@9: { /* xB[i] violates its upper bound, which plays the role alpar@9: of an lower bound on phase I */ alpar@9: t = (ub[k] - bbar[i]) / alfa; alpar@9: i_stat = GLP_NU; alpar@9: } alpar@9: else if (phase == 1 && coef[k] < 0.0) alpar@9: { /* xB[i] violates its lower bound, which plays the role alpar@9: of an upper bound on phase I */ alpar@9: continue; alpar@9: } alpar@9: else if (type[k] == GLP_LO || type[k] == GLP_DB || alpar@9: type[k] == GLP_FX) alpar@9: { /* xB[i] is within its bounds and has an lower bound */ alpar@9: t = (lb[k] - bbar[i]) / alfa; alpar@9: i_stat = GLP_NL; alpar@9: } alpar@9: else alpar@9: { /* xB[i] is within its bounds and has no lower 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 xN[q], on which xB[i] reaches its bound; alpar@9: if t <= tmax, all basic variables can violate their bounds alpar@9: only within relaxation tolerance delta; we can use this alpar@9: freedom and choose basic variable having largest influence alpar@9: coefficient to avoid possible numeric instability */ alpar@9: if (t <= tmax && big < fabs(alfa)) alpar@9: p = i, p_stat = i_stat, teta = t, big = fabs(alfa); alpar@9: } alpar@9: /* something must be chosen on the second pass */ alpar@9: xassert(p != 0); alpar@9: done: /* store the index and status of basic variable xB[p] chosen */ alpar@9: csa->p = p; alpar@9: if (p > 0 && type[head[p]] == GLP_FX) alpar@9: csa->p_stat = GLP_NS; alpar@9: else alpar@9: csa->p_stat = p_stat; alpar@9: /* store corresponding change of non-basic variable xN[q] */ alpar@9: #ifdef GLP_DEBUG alpar@9: xassert(teta >= 0.0); alpar@9: #endif alpar@9: csa->teta = s * teta; alpar@9: return; alpar@9: } alpar@9: 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: 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: 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: static void eval_trow(struct csa *csa, double rho[]) alpar@9: { int m = csa->m; alpar@9: int n = csa->n; alpar@9: #ifdef GLP_DEBUG alpar@9: char *stat = csa->stat; alpar@9: #endif alpar@9: int *N_ptr = csa->N_ptr; alpar@9: int *N_len = csa->N_len; alpar@9: int *N_ind = csa->N_ind; alpar@9: double *N_val = csa->N_val; 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: beg = N_ptr[i]; alpar@9: end = beg + N_len[i]; alpar@9: for (ptr = beg; ptr < end; ptr++) alpar@9: { alpar@9: #ifdef GLP_DEBUG alpar@9: j = N_ind[ptr]; alpar@9: xassert(1 <= j && j <= n); alpar@9: xassert(stat[j] != GLP_NS); alpar@9: #endif alpar@9: trow_vec[N_ind[ptr]] -= temp * N_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: /*********************************************************************** 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 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 p = csa->p; alpar@9: double teta = csa->teta; alpar@9: int i, pos; alpar@9: #ifdef GLP_DEBUG alpar@9: xassert(1 <= q && q <= n); alpar@9: xassert(p < 0 || 1 <= p && p <= m); alpar@9: #endif alpar@9: /* if xN[q] leaves the basis, compute its value in the adjacent alpar@9: basis, where it will replace xB[p] */ alpar@9: if (p > 0) alpar@9: bbar[p] = get_xN(csa, q) + teta; alpar@9: /* update values of other basic variables (except xB[p], because alpar@9: it will be replaced by xN[q]) */ 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: /* skip xB[p] */ alpar@9: if (i == p) continue; alpar@9: /* (change of xB[i]) = alfa[i,q] * (change of xN[q]) */ alpar@9: bbar[i] += tcol_vec[i] * teta; alpar@9: } alpar@9: done: return; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * reeval_cost - recompute reduced cost of non-basic variable xN[q] alpar@9: * alpar@9: * This routine recomputes reduced cost of non-basic variable xN[q] for alpar@9: * the current basis more accurately using its direct definition: alpar@9: * alpar@9: * d[q] = cN[q] - N'[q] * pi = alpar@9: * alpar@9: * = cN[q] - N'[q] * (inv(B') * cB) = alpar@9: * alpar@9: * = cN[q] - (cB' * inv(B) * N[q]) = alpar@9: * alpar@9: * = cN[q] + cB' * (pivot column). alpar@9: * alpar@9: * It is assumed that the pivot column of the simplex table is already alpar@9: * computed. */ alpar@9: alpar@9: static double reeval_cost(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: double *coef = csa->coef; alpar@9: int *head = csa->head; 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 dq; alpar@9: #ifdef GLP_DEBUG alpar@9: xassert(1 <= q && q <= n); alpar@9: #endif alpar@9: dq = coef[head[m+q]]; 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: dq += coef[head[i]] * tcol_vec[i]; alpar@9: } alpar@9: return dq; alpar@9: } 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 q = csa->q; 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 j, pos; alpar@9: double new_dq; alpar@9: #ifdef GLP_DEBUG alpar@9: xassert(1 <= q && q <= n); alpar@9: #endif alpar@9: /* compute reduced cost of xB[p] in the adjacent basis, where it alpar@9: will replace xN[q] */ alpar@9: #ifdef GLP_DEBUG alpar@9: xassert(trow_vec[q] != 0.0); alpar@9: #endif alpar@9: new_dq = (cbar[q] /= trow_vec[q]); alpar@9: /* update reduced costs of other non-basic variables (except alpar@9: xN[q], because it will be replaced by xB[p]) */ alpar@9: for (pos = 1; pos <= trow_nnz; pos++) alpar@9: { j = trow_ind[pos]; alpar@9: /* skip xN[q] */ alpar@9: if (j == q) continue; alpar@9: cbar[j] -= trow_vec[j] * new_dq; alpar@9: } alpar@9: 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 *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 *refsp = csa->refsp; alpar@9: double *gamma = csa->gamma; 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 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: double *u = csa->work3; alpar@9: int i, j, k, pos, beg, end, ptr; alpar@9: double gamma_q, delta_q, pivot, s, 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[q] for the current basis more accurately and alpar@9: compute auxiliary vector u */ alpar@9: gamma_q = delta_q = (refsp[head[m+q]] ? 1.0 : 0.0); alpar@9: for (i = 1; i <= m; i++) u[i] = 0.0; alpar@9: for (pos = 1; pos <= tcol_nnz; pos++) alpar@9: { i = tcol_ind[pos]; alpar@9: if (refsp[head[i]]) alpar@9: { u[i] = t = tcol_vec[i]; alpar@9: gamma_q += t * t; alpar@9: } alpar@9: else alpar@9: u[i] = 0.0; alpar@9: } alpar@9: xassert(csa->valid); alpar@9: bfd_btran(csa->bfd, u); alpar@9: /* update gamma[k] for other non-basic variables (except fixed alpar@9: variables and xN[q], because it will be replaced by xB[p]) */ alpar@9: pivot = trow_vec[q]; alpar@9: #ifdef GLP_DEBUG alpar@9: xassert(pivot != 0.0); alpar@9: #endif alpar@9: for (pos = 1; pos <= trow_nnz; pos++) alpar@9: { j = trow_ind[pos]; alpar@9: /* skip xN[q] */ alpar@9: if (j == q) continue; alpar@9: /* compute t */ alpar@9: t = trow_vec[j] / pivot; alpar@9: /* compute inner product s = N'[j] * u */ alpar@9: k = head[m+j]; /* x[k] = xN[j] */ alpar@9: if (k <= m) alpar@9: s = u[k]; alpar@9: else alpar@9: { s = 0.0; 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: s -= A_val[ptr] * u[A_ind[ptr]]; alpar@9: } alpar@9: /* compute gamma[k] for the adjacent basis */ alpar@9: t1 = gamma[j] + t * t * gamma_q + 2.0 * t * s; alpar@9: t2 = (refsp[k] ? 1.0 : 0.0) + delta_q * t * t; alpar@9: gamma[j] = (t1 >= t2 ? t1 : t2); alpar@9: if (gamma[j] < DBL_EPSILON) gamma[j] = DBL_EPSILON; alpar@9: } alpar@9: /* compute gamma[q] for the adjacent basis */ alpar@9: if (type[head[p]] == GLP_FX) alpar@9: gamma[q] = 1.0; alpar@9: else alpar@9: { gamma[q] = gamma_q / (pivot * pivot); alpar@9: if (gamma[q] < DBL_EPSILON) gamma[q] = DBL_EPSILON; alpar@9: } alpar@9: return; alpar@9: } alpar@9: 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: 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: 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 n = csa->n; alpar@9: char *stat = csa->stat; alpar@9: double *gamma = csa->gamma; alpar@9: int j; alpar@9: double e, emax, temp; alpar@9: emax = 0.0; alpar@9: for (j = 1; j <= n; j++) alpar@9: { if (stat[j] == GLP_NS) alpar@9: { xassert(gamma[j] == 1.0); alpar@9: continue; alpar@9: } alpar@9: temp = eval_gamma(csa, j); alpar@9: e = fabs(temp - gamma[j]) / (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: char *type = csa->type; alpar@9: #endif alpar@9: int *head = csa->head; alpar@9: char *stat = csa->stat; alpar@9: int q = csa->q; alpar@9: int p = csa->p; alpar@9: int p_stat = csa->p_stat; alpar@9: int k; alpar@9: #ifdef GLP_DEBUG alpar@9: xassert(1 <= q && q <= n); alpar@9: #endif alpar@9: if (p < 0) alpar@9: { /* xN[q] goes to its opposite bound */ alpar@9: #ifdef GLP_DEBUG alpar@9: k = head[m+q]; /* x[k] = xN[q] */ alpar@9: xassert(1 <= k && k <= m+n); alpar@9: xassert(type[k] == GLP_DB); alpar@9: #endif alpar@9: switch (stat[q]) alpar@9: { case GLP_NL: alpar@9: /* xN[q] increases */ alpar@9: stat[q] = GLP_NU; alpar@9: break; alpar@9: case GLP_NU: alpar@9: /* xN[q] decreases */ alpar@9: stat[q] = GLP_NL; alpar@9: break; alpar@9: default: alpar@9: xassert(stat != stat); alpar@9: } alpar@9: } alpar@9: else 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: k = head[p]; /* x[k] = xB[p] */ alpar@9: switch (p_stat) alpar@9: { case GLP_NL: alpar@9: /* xB[p] goes to its lower bound */ alpar@9: xassert(type[k] == GLP_LO || type[k] == GLP_DB); alpar@9: break; alpar@9: case GLP_NU: alpar@9: /* xB[p] goes to its upper bound */ alpar@9: xassert(type[k] == GLP_UP || type[k] == GLP_DB); alpar@9: break; alpar@9: case GLP_NS: alpar@9: /* xB[p] goes to its fixed value */ alpar@9: xassert(type[k] == GLP_NS); alpar@9: break; alpar@9: default: alpar@9: xassert(p_stat != p_stat); alpar@9: } 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: stat[q] = (char)p_stat; alpar@9: } alpar@9: return; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * set_aux_obj - construct auxiliary objective function alpar@9: * alpar@9: * The auxiliary objective function is a separable piecewise linear alpar@9: * convex function, which is the sum of primal infeasibilities: alpar@9: * alpar@9: * z = t[1] + ... + t[m+n] -> minimize, alpar@9: * alpar@9: * where: alpar@9: * alpar@9: * / lb[k] - x[k], if x[k] < lb[k] alpar@9: * | alpar@9: * t[k] = < 0, if lb[k] <= x[k] <= ub[k] alpar@9: * | alpar@9: * \ x[k] - ub[k], if x[k] > ub[k] alpar@9: * alpar@9: * This routine computes objective coefficients for the current basis alpar@9: * and returns the number of non-zero terms t[k]. */ alpar@9: alpar@9: static int set_aux_obj(struct csa *csa, double tol_bnd) 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: int *head = csa->head; alpar@9: double *bbar = csa->bbar; alpar@9: int i, k, cnt = 0; alpar@9: double eps; alpar@9: /* use a bit more restrictive tolerance */ alpar@9: tol_bnd *= 0.90; alpar@9: /* clear all objective coefficients */ alpar@9: for (k = 1; k <= m+n; k++) alpar@9: coef[k] = 0.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: if (type[k] == GLP_LO || type[k] == GLP_DB || alpar@9: type[k] == GLP_FX) alpar@9: { /* x[k] 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 violates it */ alpar@9: coef[k] = -1.0; alpar@9: cnt++; alpar@9: } alpar@9: } alpar@9: if (type[k] == GLP_UP || type[k] == GLP_DB || alpar@9: type[k] == GLP_FX) alpar@9: { /* x[k] 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 violates it */ alpar@9: coef[k] = +1.0; alpar@9: cnt++; alpar@9: } alpar@9: } alpar@9: } alpar@9: return cnt; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * set_orig_obj - restore original objective function alpar@9: * alpar@9: * This routine assigns scaled original objective coefficients to the alpar@9: * working objective function. */ alpar@9: alpar@9: static void set_orig_obj(struct csa *csa) alpar@9: { int m = csa->m; alpar@9: int n = csa->n; alpar@9: double *coef = csa->coef; alpar@9: double *obj = csa->obj; alpar@9: double zeta = csa->zeta; alpar@9: int i, j; alpar@9: for (i = 1; i <= m; i++) alpar@9: coef[i] = 0.0; alpar@9: for (j = 1; j <= n; j++) alpar@9: coef[m+j] = zeta * obj[j]; 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 primal feasible (or pseudo feasible alpar@9: * on phase I) within a tolerance, this routine returns zero, otherwise alpar@9: * it returns non-zero. */ alpar@9: alpar@9: static int check_stab(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: double *coef = csa->coef; alpar@9: int *head = csa->head; alpar@9: int phase = csa->phase; alpar@9: double *bbar = csa->bbar; alpar@9: int i, k; alpar@9: double eps; 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 (phase == 1 && coef[k] < 0.0) alpar@9: { /* x[k] must not be greater than its lower bound */ alpar@9: #ifdef GLP_DEBUG alpar@9: xassert(type[k] == GLP_LO || type[k] == GLP_DB || alpar@9: type[k] == GLP_FX); alpar@9: #endif alpar@9: eps = tol_bnd * (1.0 + kappa * fabs(lb[k])); alpar@9: if (bbar[i] > lb[k] + eps) return 1; alpar@9: } alpar@9: else if (phase == 1 && coef[k] > 0.0) alpar@9: { /* x[k] must not be less than its upper bound */ alpar@9: #ifdef GLP_DEBUG alpar@9: xassert(type[k] == GLP_UP || type[k] == GLP_DB || alpar@9: type[k] == GLP_FX); alpar@9: #endif alpar@9: eps = tol_bnd * (1.0 + kappa * fabs(ub[k])); alpar@9: if (bbar[i] < ub[k] - eps) return 1; alpar@9: } alpar@9: else alpar@9: { /* either phase = 1 and coef[k] = 0, or phase = 2 */ alpar@9: if (type[k] == GLP_LO || type[k] == GLP_DB || alpar@9: type[k] == GLP_FX) alpar@9: { /* x[k] must not be less than its lower bound */ alpar@9: eps = tol_bnd * (1.0 + kappa * fabs(lb[k])); alpar@9: if (bbar[i] < lb[k] - eps) return 1; alpar@9: } alpar@9: if (type[k] == GLP_UP || type[k] == GLP_DB || alpar@9: type[k] == GLP_FX) alpar@9: { /* x[k] must not be greater then its upper bound */ alpar@9: eps = tol_bnd * (1.0 + kappa * fabs(ub[k])); alpar@9: if (bbar[i] > ub[k] + eps) return 1; alpar@9: } alpar@9: } alpar@9: } alpar@9: /* basic solution is primal feasible within a tolerance */ alpar@9: return 0; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * check_feas - check primal feasibility of basic solution alpar@9: * alpar@9: * If the current basic solution is primal 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_bnd) alpar@9: { int m = csa->m; alpar@9: #ifdef GLP_DEBUG alpar@9: int n = csa->n; alpar@9: char *type = csa->type; alpar@9: #endif alpar@9: double *lb = csa->lb; alpar@9: double *ub = csa->ub; alpar@9: double *coef = csa->coef; alpar@9: int *head = csa->head; alpar@9: double *bbar = csa->bbar; alpar@9: int i, k; alpar@9: double eps; alpar@9: xassert(csa->phase == 1); 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 (coef[k] < 0.0) alpar@9: { /* check if x[k] still violates its lower bound */ alpar@9: #ifdef GLP_DEBUG alpar@9: xassert(type[k] == GLP_LO || type[k] == GLP_DB || alpar@9: type[k] == GLP_FX); alpar@9: #endif alpar@9: eps = tol_bnd * (1.0 + kappa * fabs(lb[k])); alpar@9: if (bbar[i] < lb[k] - eps) return 1; alpar@9: } alpar@9: else if (coef[k] > 0.0) alpar@9: { /* check if x[k] still violates its upper bound */ alpar@9: #ifdef GLP_DEBUG alpar@9: xassert(type[k] == GLP_UP || type[k] == GLP_DB || alpar@9: type[k] == GLP_FX); alpar@9: #endif alpar@9: eps = tol_bnd * (1.0 + kappa * fabs(ub[k])); alpar@9: if (bbar[i] > ub[k] + eps) return 1; alpar@9: } alpar@9: } alpar@9: /* basic solution is primal feasible within a tolerance */ alpar@9: return 0; alpar@9: } alpar@9: 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: alpar@9: /*********************************************************************** alpar@9: * display - display the search progress alpar@9: * alpar@9: * This routine displays some information about the search progress alpar@9: * that includes: alpar@9: * alpar@9: * the search phase; alpar@9: * alpar@9: * the number of simplex iterations performed by the solver; alpar@9: * alpar@9: * the original objective value; alpar@9: * alpar@9: * the sum of (scaled) primal infeasibilities; alpar@9: * alpar@9: * the number of basic fixed variables. */ alpar@9: alpar@9: static void display(struct csa *csa, const glp_smcp *parm, int spec) 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 phase = csa->phase; alpar@9: int *head = csa->head; alpar@9: double *bbar = csa->bbar; alpar@9: int i, k, 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 primal infeasibilities and determine the alpar@9: number of basic fixed variables */ alpar@9: sum = 0.0, cnt = 0; 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_LO || type[k] == GLP_DB || alpar@9: type[k] == GLP_FX) alpar@9: { /* x[k] has lower bound */ alpar@9: if (bbar[i] < lb[k]) alpar@9: sum += (lb[k] - bbar[i]); alpar@9: } alpar@9: if (type[k] == GLP_UP || type[k] == GLP_DB || alpar@9: type[k] == GLP_FX) alpar@9: { /* x[k] has upper bound */ alpar@9: if (bbar[i] > ub[k]) alpar@9: sum += (bbar[i] - ub[k]); alpar@9: } alpar@9: if (type[k] == GLP_FX) cnt++; alpar@9: } alpar@9: xprintf("%c%6d: obj = %17.9e infeas = %10.3e (%d)\n", alpar@9: phase == 1 ? ' ' : '*', 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: /*********************************************************************** 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: 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->obj); alpar@9: xfree(csa->A_ptr); alpar@9: xfree(csa->A_ind); alpar@9: xfree(csa->A_val); alpar@9: xfree(csa->head); alpar@9: xfree(csa->stat); 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: xfree(csa->bbar); alpar@9: xfree(csa->cbar); alpar@9: xfree(csa->refsp); alpar@9: xfree(csa->gamma); alpar@9: xfree(csa->tcol_ind); alpar@9: xfree(csa->tcol_vec); alpar@9: xfree(csa->trow_ind); alpar@9: xfree(csa->trow_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_primal - core LP solver based on the primal simplex method alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpspx.h" alpar@9: * int spx_primal(glp_prob *lp, const glp_smcp *parm); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine spx_primal is a core LP solver based on the two-phase alpar@9: * primal simplex method. alpar@9: * alpar@9: * RETURNS alpar@9: * alpar@9: * 0 LP instance has been successfully solved. 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_primal(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 primal values of basic variables */ alpar@9: if (bbar_st == 0) alpar@9: { eval_bbar(csa); alpar@9: bbar_st = 1; /* just computed */ alpar@9: /* determine the search phase, if not determined yet */ alpar@9: if (csa->phase == 0) alpar@9: { if (set_aux_obj(csa, parm->tol_bnd) > 0) alpar@9: { /* current basic solution is primal infeasible */ alpar@9: /* start to minimize the sum of infeasibilities */ alpar@9: csa->phase = 1; alpar@9: } alpar@9: else alpar@9: { /* current basic solution is primal feasible */ alpar@9: /* start to minimize the original objective function */ alpar@9: set_orig_obj(csa); alpar@9: csa->phase = 2; alpar@9: } alpar@9: xassert(check_stab(csa, parm->tol_bnd) == 0); alpar@9: /* working objective coefficients have been changed, so alpar@9: invalidate reduced costs */ alpar@9: cbar_st = 0; alpar@9: display(csa, parm, 1); alpar@9: } alpar@9: /* make sure that the current basic solution remains primal alpar@9: feasible (or pseudo feasible on phase I) */ alpar@9: if (check_stab(csa, parm->tol_bnd)) alpar@9: { /* there are excessive bound violations due to round-off alpar@9: errors */ alpar@9: if (parm->msg_lev >= GLP_MSG_ERR) alpar@9: xprintf("Warning: numerical instability (primal simplex," alpar@9: " phase %s)\n", csa->phase == 1 ? "I" : "II"); 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 dual feasible; it is sufficient to make sure alpar@9: that no basic variable violates its bounds */ alpar@9: if (csa->phase == 1 && !check_feas(csa, parm->tol_bnd)) alpar@9: { /* the current basis is primal feasible; switch to phase II */ alpar@9: csa->phase = 2; alpar@9: set_orig_obj(csa); alpar@9: cbar_st = 0; alpar@9: display(csa, parm, 1); 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: } 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: /* 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 (bbar_st != 1 || csa->phase == 2 && cbar_st != 1) alpar@9: { if (bbar_st != 1) bbar_st = 0; alpar@9: if (csa->phase == 2 && 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: p_stat = GLP_INFEAS; alpar@9: set_orig_obj(csa); alpar@9: eval_cbar(csa); alpar@9: break; alpar@9: case 2: alpar@9: p_stat = GLP_FEAS; alpar@9: break; alpar@9: default: alpar@9: xassert(csa != csa); alpar@9: } alpar@9: chuzc(csa, parm->tol_dj); alpar@9: d_stat = (csa->q == 0 ? GLP_FEAS : GLP_INFEAS); alpar@9: store_sol(csa, lp, p_stat, 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 (bbar_st != 1 || csa->phase == 2 && cbar_st != 1) alpar@9: { if (bbar_st != 1) bbar_st = 0; alpar@9: if (csa->phase == 2 && 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: p_stat = GLP_INFEAS; alpar@9: set_orig_obj(csa); alpar@9: eval_cbar(csa); alpar@9: break; alpar@9: case 2: alpar@9: p_stat = GLP_FEAS; alpar@9: break; alpar@9: default: alpar@9: xassert(csa != csa); alpar@9: } alpar@9: chuzc(csa, parm->tol_dj); alpar@9: d_stat = (csa->q == 0 ? GLP_FEAS : GLP_INFEAS); alpar@9: store_sol(csa, lp, p_stat, 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 non-basic variable xN[q] */ alpar@9: chuzc(csa, parm->tol_dj); alpar@9: if (csa->q == 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 FEASIBLE SOLUTION\n"); alpar@9: p_stat = GLP_NOFEAS; alpar@9: set_orig_obj(csa); alpar@9: eval_cbar(csa); alpar@9: chuzc(csa, parm->tol_dj); alpar@9: d_stat = (csa->q == 0 ? GLP_FEAS : GLP_INFEAS); 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 column of the simplex table */ alpar@9: eval_tcol(csa); alpar@9: if (rigorous) refine_tcol(csa); alpar@9: sort_tcol(csa, parm->tol_piv); alpar@9: /* check accuracy of the reduced cost of xN[q] */ alpar@9: { double d1 = csa->cbar[csa->q]; /* less accurate */ alpar@9: double d2 = reeval_cost(csa); /* more accurate */ alpar@9: xassert(d1 != 0.0); alpar@9: if (fabs(d1 - d2) > 1e-5 * (1.0 + fabs(d2)) || alpar@9: !(d1 < 0.0 && d2 < 0.0 || d1 > 0.0 && d2 > 0.0)) alpar@9: { if (parm->msg_lev >= GLP_MSG_DBG) alpar@9: xprintf("d1 = %.12g; d2 = %.12g\n", d1, d2); alpar@9: if (cbar_st != 1 || !rigorous) alpar@9: { if (cbar_st != 1) cbar_st = 0; alpar@9: rigorous = 5; alpar@9: goto loop; alpar@9: } alpar@9: } alpar@9: /* replace cbar[q] by more accurate value keeping its sign */ alpar@9: if (d1 > 0.0) alpar@9: csa->cbar[csa->q] = (d2 > 0.0 ? d2 : +DBL_EPSILON); alpar@9: else alpar@9: csa->cbar[csa->q] = (d2 < 0.0 ? d2 : -DBL_EPSILON); alpar@9: } alpar@9: /* choose basic variable xB[p] */ alpar@9: switch (parm->r_test) alpar@9: { case GLP_RT_STD: alpar@9: chuzr(csa, 0.0); alpar@9: break; alpar@9: case GLP_RT_HAR: alpar@9: chuzr(csa, 0.30 * parm->tol_bnd); alpar@9: break; alpar@9: default: alpar@9: xassert(parm != parm); alpar@9: } alpar@9: if (csa->p == 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 UNBOUNDED SOLUTION\n"); alpar@9: store_sol(csa, lp, GLP_FEAS, GLP_NOFEAS, alpar@9: csa->head[csa->m+csa->q]); 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: if (csa->p > 0) alpar@9: { double piv = csa->tcol_vec[csa->p]; alpar@9: double eps = 1e-5 * (1.0 + 0.01 * csa->tcol_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 row of the simplex table */ alpar@9: if (csa->p > 0) 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: } alpar@9: /* accuracy check based on the pivot element */ alpar@9: if (csa->p > 0) 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: /* use more accurate version in the pivot row */ alpar@9: if (csa->trow_vec[csa->q] == 0.0) alpar@9: { csa->trow_nnz++; alpar@9: xassert(csa->trow_nnz <= csa->n); alpar@9: csa->trow_ind[csa->trow_nnz] = csa->q; alpar@9: } alpar@9: csa->trow_vec[csa->q] = piv1; alpar@9: } alpar@9: } alpar@9: /* update primal values of basic variables */ alpar@9: update_bbar(csa); alpar@9: bbar_st = 2; /* updated */ alpar@9: /* update reduced costs of non-basic variables */ alpar@9: if (csa->p > 0) alpar@9: { update_cbar(csa); alpar@9: cbar_st = 2; /* updated */ alpar@9: /* on phase I objective coefficient of xB[p] in the adjacent alpar@9: basis becomes zero */ alpar@9: if (csa->phase == 1) alpar@9: { int k = csa->head[csa->p]; /* x[k] = xB[p] -> xN[q] */ alpar@9: csa->cbar[csa->q] -= csa->coef[k]; alpar@9: csa->coef[k] = 0.0; alpar@9: } alpar@9: } alpar@9: /* update steepest edge coefficients */ alpar@9: if (csa->p > 0) 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: } alpar@9: /* update factorization of the basis matrix */ alpar@9: if (csa->p > 0) 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: } alpar@9: /* update matrix N */ alpar@9: if (csa->p > 0) 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: } 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 */