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