lemon-project-template-glpk
diff deps/glpk/src/glpini02.c @ 9:33de93886c88
Import GLPK 4.47
author | Alpar Juttner <alpar@cs.elte.hu> |
---|---|
date | Sun, 06 Nov 2011 20:59:10 +0100 |
parents | |
children |
line diff
1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 1.2 +++ b/deps/glpk/src/glpini02.c Sun Nov 06 20:59:10 2011 +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, 2011 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 */