1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glpini02.c Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,269 @@
1.4 +/* glpini02.c */
1.5 +
1.6 +/***********************************************************************
1.7 +* This code is part of GLPK (GNU Linear Programming Kit).
1.8 +*
1.9 +* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
1.10 +* 2009, 2010 Andrew Makhorin, Department for Applied Informatics,
1.11 +* Moscow Aviation Institute, Moscow, Russia. All rights reserved.
1.12 +* E-mail: <mao@gnu.org>.
1.13 +*
1.14 +* GLPK is free software: you can redistribute it and/or modify it
1.15 +* under the terms of the GNU General Public License as published by
1.16 +* the Free Software Foundation, either version 3 of the License, or
1.17 +* (at your option) any later version.
1.18 +*
1.19 +* GLPK is distributed in the hope that it will be useful, but WITHOUT
1.20 +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.21 +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
1.22 +* License for more details.
1.23 +*
1.24 +* You should have received a copy of the GNU General Public License
1.25 +* along with GLPK. If not, see <http://www.gnu.org/licenses/>.
1.26 +***********************************************************************/
1.27 +
1.28 +#include "glpapi.h"
1.29 +
1.30 +struct var
1.31 +{ /* structural variable */
1.32 + int j;
1.33 + /* ordinal number */
1.34 + double q;
1.35 + /* penalty value */
1.36 +};
1.37 +
1.38 +static int fcmp(const void *ptr1, const void *ptr2)
1.39 +{ /* this routine is passed to the qsort() function */
1.40 + struct var *col1 = (void *)ptr1, *col2 = (void *)ptr2;
1.41 + if (col1->q < col2->q) return -1;
1.42 + if (col1->q > col2->q) return +1;
1.43 + return 0;
1.44 +}
1.45 +
1.46 +static int get_column(glp_prob *lp, int j, int ind[], double val[])
1.47 +{ /* Bixby's algorithm assumes that the constraint matrix is scaled
1.48 + such that the maximum absolute value in every non-zero row and
1.49 + column is 1 */
1.50 + int k, len;
1.51 + double big;
1.52 + len = glp_get_mat_col(lp, j, ind, val);
1.53 + big = 0.0;
1.54 + for (k = 1; k <= len; k++)
1.55 + if (big < fabs(val[k])) big = fabs(val[k]);
1.56 + if (big == 0.0) big = 1.0;
1.57 + for (k = 1; k <= len; k++) val[k] /= big;
1.58 + return len;
1.59 +}
1.60 +
1.61 +static void cpx_basis(glp_prob *lp)
1.62 +{ /* main routine */
1.63 + struct var *C, *C2, *C3, *C4;
1.64 + int m, n, i, j, jk, k, l, ll, t, n2, n3, n4, type, len, *I, *r,
1.65 + *ind;
1.66 + double alpha, gamma, cmax, temp, *v, *val;
1.67 + xprintf("Constructing initial basis...\n");
1.68 + /* determine the number of rows and columns */
1.69 + m = glp_get_num_rows(lp);
1.70 + n = glp_get_num_cols(lp);
1.71 + /* allocate working arrays */
1.72 + C = xcalloc(1+n, sizeof(struct var));
1.73 + I = xcalloc(1+m, sizeof(int));
1.74 + r = xcalloc(1+m, sizeof(int));
1.75 + v = xcalloc(1+m, sizeof(double));
1.76 + ind = xcalloc(1+m, sizeof(int));
1.77 + val = xcalloc(1+m, sizeof(double));
1.78 + /* make all auxiliary variables non-basic */
1.79 + for (i = 1; i <= m; i++)
1.80 + { if (glp_get_row_type(lp, i) != GLP_DB)
1.81 + glp_set_row_stat(lp, i, GLP_NS);
1.82 + else if (fabs(glp_get_row_lb(lp, i)) <=
1.83 + fabs(glp_get_row_ub(lp, i)))
1.84 + glp_set_row_stat(lp, i, GLP_NL);
1.85 + else
1.86 + glp_set_row_stat(lp, i, GLP_NU);
1.87 + }
1.88 + /* make all structural variables non-basic */
1.89 + for (j = 1; j <= n; j++)
1.90 + { if (glp_get_col_type(lp, j) != GLP_DB)
1.91 + glp_set_col_stat(lp, j, GLP_NS);
1.92 + else if (fabs(glp_get_col_lb(lp, j)) <=
1.93 + fabs(glp_get_col_ub(lp, j)))
1.94 + glp_set_col_stat(lp, j, GLP_NL);
1.95 + else
1.96 + glp_set_col_stat(lp, j, GLP_NU);
1.97 + }
1.98 + /* C2 is a set of free structural variables */
1.99 + n2 = 0, C2 = C + 0;
1.100 + for (j = 1; j <= n; j++)
1.101 + { type = glp_get_col_type(lp, j);
1.102 + if (type == GLP_FR)
1.103 + { n2++;
1.104 + C2[n2].j = j;
1.105 + C2[n2].q = 0.0;
1.106 + }
1.107 + }
1.108 + /* C3 is a set of structural variables having excatly one (lower
1.109 + or upper) bound */
1.110 + n3 = 0, C3 = C2 + n2;
1.111 + for (j = 1; j <= n; j++)
1.112 + { type = glp_get_col_type(lp, j);
1.113 + if (type == GLP_LO)
1.114 + { n3++;
1.115 + C3[n3].j = j;
1.116 + C3[n3].q = + glp_get_col_lb(lp, j);
1.117 + }
1.118 + else if (type == GLP_UP)
1.119 + { n3++;
1.120 + C3[n3].j = j;
1.121 + C3[n3].q = - glp_get_col_ub(lp, j);
1.122 + }
1.123 + }
1.124 + /* C4 is a set of structural variables having both (lower and
1.125 + upper) bounds */
1.126 + n4 = 0, C4 = C3 + n3;
1.127 + for (j = 1; j <= n; j++)
1.128 + { type = glp_get_col_type(lp, j);
1.129 + if (type == GLP_DB)
1.130 + { n4++;
1.131 + C4[n4].j = j;
1.132 + C4[n4].q = glp_get_col_lb(lp, j) - glp_get_col_ub(lp, j);
1.133 + }
1.134 + }
1.135 + /* compute gamma = max{|c[j]|: 1 <= j <= n} */
1.136 + gamma = 0.0;
1.137 + for (j = 1; j <= n; j++)
1.138 + { temp = fabs(glp_get_obj_coef(lp, j));
1.139 + if (gamma < temp) gamma = temp;
1.140 + }
1.141 + /* compute cmax */
1.142 + cmax = (gamma == 0.0 ? 1.0 : 1000.0 * gamma);
1.143 + /* compute final penalty for all structural variables within sets
1.144 + C2, C3, and C4 */
1.145 + switch (glp_get_obj_dir(lp))
1.146 + { case GLP_MIN: temp = +1.0; break;
1.147 + case GLP_MAX: temp = -1.0; break;
1.148 + default: xassert(lp != lp);
1.149 + }
1.150 + for (k = 1; k <= n2+n3+n4; k++)
1.151 + { j = C[k].j;
1.152 + C[k].q += (temp * glp_get_obj_coef(lp, j)) / cmax;
1.153 + }
1.154 + /* sort structural variables within C2, C3, and C4 in ascending
1.155 + order of penalty value */
1.156 + qsort(C2+1, n2, sizeof(struct var), fcmp);
1.157 + for (k = 1; k < n2; k++) xassert(C2[k].q <= C2[k+1].q);
1.158 + qsort(C3+1, n3, sizeof(struct var), fcmp);
1.159 + for (k = 1; k < n3; k++) xassert(C3[k].q <= C3[k+1].q);
1.160 + qsort(C4+1, n4, sizeof(struct var), fcmp);
1.161 + for (k = 1; k < n4; k++) xassert(C4[k].q <= C4[k+1].q);
1.162 + /*** STEP 1 ***/
1.163 + for (i = 1; i <= m; i++)
1.164 + { type = glp_get_row_type(lp, i);
1.165 + if (type != GLP_FX)
1.166 + { /* row i is either free or inequality constraint */
1.167 + glp_set_row_stat(lp, i, GLP_BS);
1.168 + I[i] = 1;
1.169 + r[i] = 1;
1.170 + }
1.171 + else
1.172 + { /* row i is equality constraint */
1.173 + I[i] = 0;
1.174 + r[i] = 0;
1.175 + }
1.176 + v[i] = +DBL_MAX;
1.177 + }
1.178 + /*** STEP 2 ***/
1.179 + for (k = 1; k <= n2+n3+n4; k++)
1.180 + { jk = C[k].j;
1.181 + len = get_column(lp, jk, ind, val);
1.182 + /* let alpha = max{|A[l,jk]|: r[l] = 0} and let l' be such
1.183 + that alpha = |A[l',jk]| */
1.184 + alpha = 0.0, ll = 0;
1.185 + for (t = 1; t <= len; t++)
1.186 + { l = ind[t];
1.187 + if (r[l] == 0 && alpha < fabs(val[t]))
1.188 + alpha = fabs(val[t]), ll = l;
1.189 + }
1.190 + if (alpha >= 0.99)
1.191 + { /* B := B union {jk} */
1.192 + glp_set_col_stat(lp, jk, GLP_BS);
1.193 + I[ll] = 1;
1.194 + v[ll] = alpha;
1.195 + /* r[l] := r[l] + 1 for all l such that |A[l,jk]| != 0 */
1.196 + for (t = 1; t <= len; t++)
1.197 + { l = ind[t];
1.198 + if (val[t] != 0.0) r[l]++;
1.199 + }
1.200 + /* continue to the next k */
1.201 + continue;
1.202 + }
1.203 + /* if |A[l,jk]| > 0.01 * v[l] for some l, continue to the
1.204 + next k */
1.205 + for (t = 1; t <= len; t++)
1.206 + { l = ind[t];
1.207 + if (fabs(val[t]) > 0.01 * v[l]) break;
1.208 + }
1.209 + if (t <= len) continue;
1.210 + /* otherwise, let alpha = max{|A[l,jk]|: I[l] = 0} and let l'
1.211 + be such that alpha = |A[l',jk]| */
1.212 + alpha = 0.0, ll = 0;
1.213 + for (t = 1; t <= len; t++)
1.214 + { l = ind[t];
1.215 + if (I[l] == 0 && alpha < fabs(val[t]))
1.216 + alpha = fabs(val[t]), ll = l;
1.217 + }
1.218 + /* if alpha = 0, continue to the next k */
1.219 + if (alpha == 0.0) continue;
1.220 + /* B := B union {jk} */
1.221 + glp_set_col_stat(lp, jk, GLP_BS);
1.222 + I[ll] = 1;
1.223 + v[ll] = alpha;
1.224 + /* r[l] := r[l] + 1 for all l such that |A[l,jk]| != 0 */
1.225 + for (t = 1; t <= len; t++)
1.226 + { l = ind[t];
1.227 + if (val[t] != 0.0) r[l]++;
1.228 + }
1.229 + }
1.230 + /*** STEP 3 ***/
1.231 + /* add an artificial variable (auxiliary variable for equality
1.232 + constraint) to cover each remaining uncovered row */
1.233 + for (i = 1; i <= m; i++)
1.234 + if (I[i] == 0) glp_set_row_stat(lp, i, GLP_BS);
1.235 + /* free working arrays */
1.236 + xfree(C);
1.237 + xfree(I);
1.238 + xfree(r);
1.239 + xfree(v);
1.240 + xfree(ind);
1.241 + xfree(val);
1.242 + return;
1.243 +}
1.244 +
1.245 +/***********************************************************************
1.246 +* NAME
1.247 +*
1.248 +* glp_cpx_basis - construct Bixby's initial LP basis
1.249 +*
1.250 +* SYNOPSIS
1.251 +*
1.252 +* void glp_cpx_basis(glp_prob *lp);
1.253 +*
1.254 +* DESCRIPTION
1.255 +*
1.256 +* The routine glp_cpx_basis constructs an advanced initial basis for
1.257 +* the specified problem object.
1.258 +*
1.259 +* The routine is based on Bixby's algorithm described in the paper:
1.260 +*
1.261 +* Robert E. Bixby. Implementing the Simplex Method: The Initial Basis.
1.262 +* ORSA Journal on Computing, Vol. 4, No. 3, 1992, pp. 267-84. */
1.263 +
1.264 +void glp_cpx_basis(glp_prob *lp)
1.265 +{ if (lp->m == 0 || lp->n == 0)
1.266 + glp_std_basis(lp);
1.267 + else
1.268 + cpx_basis(lp);
1.269 + return;
1.270 +}
1.271 +
1.272 +/* eof */