1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glpnpp01.c Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,927 @@
1.4 +/* glpnpp01.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 "glpnpp.h"
1.29 +
1.30 +NPP *npp_create_wksp(void)
1.31 +{ /* create LP/MIP preprocessor workspace */
1.32 + NPP *npp;
1.33 + npp = xmalloc(sizeof(NPP));
1.34 + npp->orig_dir = 0;
1.35 + npp->orig_m = npp->orig_n = npp->orig_nnz = 0;
1.36 + npp->pool = dmp_create_pool();
1.37 + npp->name = npp->obj = NULL;
1.38 + npp->c0 = 0.0;
1.39 + npp->nrows = npp->ncols = 0;
1.40 + npp->r_head = npp->r_tail = NULL;
1.41 + npp->c_head = npp->c_tail = NULL;
1.42 + npp->stack = dmp_create_pool();
1.43 + npp->top = NULL;
1.44 +#if 0 /* 16/XII-2009 */
1.45 + memset(&npp->count, 0, sizeof(npp->count));
1.46 +#endif
1.47 + npp->m = npp->n = npp->nnz = 0;
1.48 + npp->row_ref = npp->col_ref = NULL;
1.49 + npp->sol = npp->scaling = 0;
1.50 + npp->p_stat = npp->d_stat = npp->t_stat = npp->i_stat = 0;
1.51 + npp->r_stat = NULL;
1.52 + /*npp->r_prim =*/ npp->r_pi = NULL;
1.53 + npp->c_stat = NULL;
1.54 + npp->c_value = /*npp->c_dual =*/ NULL;
1.55 + return npp;
1.56 +}
1.57 +
1.58 +void npp_insert_row(NPP *npp, NPPROW *row, int where)
1.59 +{ /* insert row to the row list */
1.60 + if (where == 0)
1.61 + { /* insert row to the beginning of the row list */
1.62 + row->prev = NULL;
1.63 + row->next = npp->r_head;
1.64 + if (row->next == NULL)
1.65 + npp->r_tail = row;
1.66 + else
1.67 + row->next->prev = row;
1.68 + npp->r_head = row;
1.69 + }
1.70 + else
1.71 + { /* insert row to the end of the row list */
1.72 + row->prev = npp->r_tail;
1.73 + row->next = NULL;
1.74 + if (row->prev == NULL)
1.75 + npp->r_head = row;
1.76 + else
1.77 + row->prev->next = row;
1.78 + npp->r_tail = row;
1.79 + }
1.80 + return;
1.81 +}
1.82 +
1.83 +void npp_remove_row(NPP *npp, NPPROW *row)
1.84 +{ /* remove row from the row list */
1.85 + if (row->prev == NULL)
1.86 + npp->r_head = row->next;
1.87 + else
1.88 + row->prev->next = row->next;
1.89 + if (row->next == NULL)
1.90 + npp->r_tail = row->prev;
1.91 + else
1.92 + row->next->prev = row->prev;
1.93 + return;
1.94 +}
1.95 +
1.96 +void npp_activate_row(NPP *npp, NPPROW *row)
1.97 +{ /* make row active */
1.98 + if (!row->temp)
1.99 + { row->temp = 1;
1.100 + /* move the row to the beginning of the row list */
1.101 + npp_remove_row(npp, row);
1.102 + npp_insert_row(npp, row, 0);
1.103 + }
1.104 + return;
1.105 +}
1.106 +
1.107 +void npp_deactivate_row(NPP *npp, NPPROW *row)
1.108 +{ /* make row inactive */
1.109 + if (row->temp)
1.110 + { row->temp = 0;
1.111 + /* move the row to the end of the row list */
1.112 + npp_remove_row(npp, row);
1.113 + npp_insert_row(npp, row, 1);
1.114 + }
1.115 + return;
1.116 +}
1.117 +
1.118 +void npp_insert_col(NPP *npp, NPPCOL *col, int where)
1.119 +{ /* insert column to the column list */
1.120 + if (where == 0)
1.121 + { /* insert column to the beginning of the column list */
1.122 + col->prev = NULL;
1.123 + col->next = npp->c_head;
1.124 + if (col->next == NULL)
1.125 + npp->c_tail = col;
1.126 + else
1.127 + col->next->prev = col;
1.128 + npp->c_head = col;
1.129 + }
1.130 + else
1.131 + { /* insert column to the end of the column list */
1.132 + col->prev = npp->c_tail;
1.133 + col->next = NULL;
1.134 + if (col->prev == NULL)
1.135 + npp->c_head = col;
1.136 + else
1.137 + col->prev->next = col;
1.138 + npp->c_tail = col;
1.139 + }
1.140 + return;
1.141 +}
1.142 +
1.143 +void npp_remove_col(NPP *npp, NPPCOL *col)
1.144 +{ /* remove column from the column list */
1.145 + if (col->prev == NULL)
1.146 + npp->c_head = col->next;
1.147 + else
1.148 + col->prev->next = col->next;
1.149 + if (col->next == NULL)
1.150 + npp->c_tail = col->prev;
1.151 + else
1.152 + col->next->prev = col->prev;
1.153 + return;
1.154 +}
1.155 +
1.156 +void npp_activate_col(NPP *npp, NPPCOL *col)
1.157 +{ /* make column active */
1.158 + if (!col->temp)
1.159 + { col->temp = 1;
1.160 + /* move the column to the beginning of the column list */
1.161 + npp_remove_col(npp, col);
1.162 + npp_insert_col(npp, col, 0);
1.163 + }
1.164 + return;
1.165 +}
1.166 +
1.167 +void npp_deactivate_col(NPP *npp, NPPCOL *col)
1.168 +{ /* make column inactive */
1.169 + if (col->temp)
1.170 + { col->temp = 0;
1.171 + /* move the column to the end of the column list */
1.172 + npp_remove_col(npp, col);
1.173 + npp_insert_col(npp, col, 1);
1.174 + }
1.175 + return;
1.176 +}
1.177 +
1.178 +NPPROW *npp_add_row(NPP *npp)
1.179 +{ /* add new row to the current problem */
1.180 + NPPROW *row;
1.181 + row = dmp_get_atom(npp->pool, sizeof(NPPROW));
1.182 + row->i = ++(npp->nrows);
1.183 + row->name = NULL;
1.184 + row->lb = -DBL_MAX, row->ub = +DBL_MAX;
1.185 + row->ptr = NULL;
1.186 + row->temp = 0;
1.187 + npp_insert_row(npp, row, 1);
1.188 + return row;
1.189 +}
1.190 +
1.191 +NPPCOL *npp_add_col(NPP *npp)
1.192 +{ /* add new column to the current problem */
1.193 + NPPCOL *col;
1.194 + col = dmp_get_atom(npp->pool, sizeof(NPPCOL));
1.195 + col->j = ++(npp->ncols);
1.196 + col->name = NULL;
1.197 +#if 0
1.198 + col->kind = GLP_CV;
1.199 +#else
1.200 + col->is_int = 0;
1.201 +#endif
1.202 + col->lb = col->ub = col->coef = 0.0;
1.203 + col->ptr = NULL;
1.204 + col->temp = 0;
1.205 + npp_insert_col(npp, col, 1);
1.206 + return col;
1.207 +}
1.208 +
1.209 +NPPAIJ *npp_add_aij(NPP *npp, NPPROW *row, NPPCOL *col, double val)
1.210 +{ /* add new element to the constraint matrix */
1.211 + NPPAIJ *aij;
1.212 + aij = dmp_get_atom(npp->pool, sizeof(NPPAIJ));
1.213 + aij->row = row;
1.214 + aij->col = col;
1.215 + aij->val = val;
1.216 + aij->r_prev = NULL;
1.217 + aij->r_next = row->ptr;
1.218 + aij->c_prev = NULL;
1.219 + aij->c_next = col->ptr;
1.220 + if (aij->r_next != NULL)
1.221 + aij->r_next->r_prev = aij;
1.222 + if (aij->c_next != NULL)
1.223 + aij->c_next->c_prev = aij;
1.224 + row->ptr = col->ptr = aij;
1.225 + return aij;
1.226 +}
1.227 +
1.228 +int npp_row_nnz(NPP *npp, NPPROW *row)
1.229 +{ /* count number of non-zero coefficients in row */
1.230 + NPPAIJ *aij;
1.231 + int nnz;
1.232 + xassert(npp == npp);
1.233 + nnz = 0;
1.234 + for (aij = row->ptr; aij != NULL; aij = aij->r_next)
1.235 + nnz++;
1.236 + return nnz;
1.237 +}
1.238 +
1.239 +int npp_col_nnz(NPP *npp, NPPCOL *col)
1.240 +{ /* count number of non-zero coefficients in column */
1.241 + NPPAIJ *aij;
1.242 + int nnz;
1.243 + xassert(npp == npp);
1.244 + nnz = 0;
1.245 + for (aij = col->ptr; aij != NULL; aij = aij->c_next)
1.246 + nnz++;
1.247 + return nnz;
1.248 +}
1.249 +
1.250 +void *npp_push_tse(NPP *npp, int (*func)(NPP *npp, void *info),
1.251 + int size)
1.252 +{ /* push new entry to the transformation stack */
1.253 + NPPTSE *tse;
1.254 + tse = dmp_get_atom(npp->stack, sizeof(NPPTSE));
1.255 + tse->func = func;
1.256 + tse->info = dmp_get_atom(npp->stack, size);
1.257 + tse->link = npp->top;
1.258 + npp->top = tse;
1.259 + return tse->info;
1.260 +}
1.261 +
1.262 +#if 1 /* 23/XII-2009 */
1.263 +void npp_erase_row(NPP *npp, NPPROW *row)
1.264 +{ /* erase row content to make it empty */
1.265 + NPPAIJ *aij;
1.266 + while (row->ptr != NULL)
1.267 + { aij = row->ptr;
1.268 + row->ptr = aij->r_next;
1.269 + if (aij->c_prev == NULL)
1.270 + aij->col->ptr = aij->c_next;
1.271 + else
1.272 + aij->c_prev->c_next = aij->c_next;
1.273 + if (aij->c_next == NULL)
1.274 + ;
1.275 + else
1.276 + aij->c_next->c_prev = aij->c_prev;
1.277 + dmp_free_atom(npp->pool, aij, sizeof(NPPAIJ));
1.278 + }
1.279 + return;
1.280 +}
1.281 +#endif
1.282 +
1.283 +void npp_del_row(NPP *npp, NPPROW *row)
1.284 +{ /* remove row from the current problem */
1.285 +#if 0 /* 23/XII-2009 */
1.286 + NPPAIJ *aij;
1.287 +#endif
1.288 + if (row->name != NULL)
1.289 + dmp_free_atom(npp->pool, row->name, strlen(row->name)+1);
1.290 +#if 0 /* 23/XII-2009 */
1.291 + while (row->ptr != NULL)
1.292 + { aij = row->ptr;
1.293 + row->ptr = aij->r_next;
1.294 + if (aij->c_prev == NULL)
1.295 + aij->col->ptr = aij->c_next;
1.296 + else
1.297 + aij->c_prev->c_next = aij->c_next;
1.298 + if (aij->c_next == NULL)
1.299 + ;
1.300 + else
1.301 + aij->c_next->c_prev = aij->c_prev;
1.302 + dmp_free_atom(npp->pool, aij, sizeof(NPPAIJ));
1.303 + }
1.304 +#else
1.305 + npp_erase_row(npp, row);
1.306 +#endif
1.307 + npp_remove_row(npp, row);
1.308 + dmp_free_atom(npp->pool, row, sizeof(NPPROW));
1.309 + return;
1.310 +}
1.311 +
1.312 +void npp_del_col(NPP *npp, NPPCOL *col)
1.313 +{ /* remove column from the current problem */
1.314 + NPPAIJ *aij;
1.315 + if (col->name != NULL)
1.316 + dmp_free_atom(npp->pool, col->name, strlen(col->name)+1);
1.317 + while (col->ptr != NULL)
1.318 + { aij = col->ptr;
1.319 + col->ptr = aij->c_next;
1.320 + if (aij->r_prev == NULL)
1.321 + aij->row->ptr = aij->r_next;
1.322 + else
1.323 + aij->r_prev->r_next = aij->r_next;
1.324 + if (aij->r_next == NULL)
1.325 + ;
1.326 + else
1.327 + aij->r_next->r_prev = aij->r_prev;
1.328 + dmp_free_atom(npp->pool, aij, sizeof(NPPAIJ));
1.329 + }
1.330 + npp_remove_col(npp, col);
1.331 + dmp_free_atom(npp->pool, col, sizeof(NPPCOL));
1.332 + return;
1.333 +}
1.334 +
1.335 +void npp_del_aij(NPP *npp, NPPAIJ *aij)
1.336 +{ /* remove element from the constraint matrix */
1.337 + if (aij->r_prev == NULL)
1.338 + aij->row->ptr = aij->r_next;
1.339 + else
1.340 + aij->r_prev->r_next = aij->r_next;
1.341 + if (aij->r_next == NULL)
1.342 + ;
1.343 + else
1.344 + aij->r_next->r_prev = aij->r_prev;
1.345 + if (aij->c_prev == NULL)
1.346 + aij->col->ptr = aij->c_next;
1.347 + else
1.348 + aij->c_prev->c_next = aij->c_next;
1.349 + if (aij->c_next == NULL)
1.350 + ;
1.351 + else
1.352 + aij->c_next->c_prev = aij->c_prev;
1.353 + dmp_free_atom(npp->pool, aij, sizeof(NPPAIJ));
1.354 + return;
1.355 +}
1.356 +
1.357 +void npp_load_prob(NPP *npp, glp_prob *orig, int names, int sol,
1.358 + int scaling)
1.359 +{ /* load original problem into the preprocessor workspace */
1.360 + int m = orig->m;
1.361 + int n = orig->n;
1.362 + NPPROW **link;
1.363 + int i, j;
1.364 + double dir;
1.365 + xassert(names == GLP_OFF || names == GLP_ON);
1.366 + xassert(sol == GLP_SOL || sol == GLP_IPT || sol == GLP_MIP);
1.367 + xassert(scaling == GLP_OFF || scaling == GLP_ON);
1.368 + if (sol == GLP_MIP) xassert(!scaling);
1.369 + npp->orig_dir = orig->dir;
1.370 + if (npp->orig_dir == GLP_MIN)
1.371 + dir = +1.0;
1.372 + else if (npp->orig_dir == GLP_MAX)
1.373 + dir = -1.0;
1.374 + else
1.375 + xassert(npp != npp);
1.376 + npp->orig_m = m;
1.377 + npp->orig_n = n;
1.378 + npp->orig_nnz = orig->nnz;
1.379 + if (names && orig->name != NULL)
1.380 + { npp->name = dmp_get_atom(npp->pool, strlen(orig->name)+1);
1.381 + strcpy(npp->name, orig->name);
1.382 + }
1.383 + if (names && orig->obj != NULL)
1.384 + { npp->obj = dmp_get_atom(npp->pool, strlen(orig->obj)+1);
1.385 + strcpy(npp->obj, orig->obj);
1.386 + }
1.387 + npp->c0 = dir * orig->c0;
1.388 + /* load rows */
1.389 + link = xcalloc(1+m, sizeof(NPPROW *));
1.390 + for (i = 1; i <= m; i++)
1.391 + { GLPROW *rrr = orig->row[i];
1.392 + NPPROW *row;
1.393 + link[i] = row = npp_add_row(npp);
1.394 + xassert(row->i == i);
1.395 + if (names && rrr->name != NULL)
1.396 + { row->name = dmp_get_atom(npp->pool, strlen(rrr->name)+1);
1.397 + strcpy(row->name, rrr->name);
1.398 + }
1.399 + if (!scaling)
1.400 + { if (rrr->type == GLP_FR)
1.401 + row->lb = -DBL_MAX, row->ub = +DBL_MAX;
1.402 + else if (rrr->type == GLP_LO)
1.403 + row->lb = rrr->lb, row->ub = +DBL_MAX;
1.404 + else if (rrr->type == GLP_UP)
1.405 + row->lb = -DBL_MAX, row->ub = rrr->ub;
1.406 + else if (rrr->type == GLP_DB)
1.407 + row->lb = rrr->lb, row->ub = rrr->ub;
1.408 + else if (rrr->type == GLP_FX)
1.409 + row->lb = row->ub = rrr->lb;
1.410 + else
1.411 + xassert(rrr != rrr);
1.412 + }
1.413 + else
1.414 + { double rii = rrr->rii;
1.415 + if (rrr->type == GLP_FR)
1.416 + row->lb = -DBL_MAX, row->ub = +DBL_MAX;
1.417 + else if (rrr->type == GLP_LO)
1.418 + row->lb = rrr->lb * rii, row->ub = +DBL_MAX;
1.419 + else if (rrr->type == GLP_UP)
1.420 + row->lb = -DBL_MAX, row->ub = rrr->ub * rii;
1.421 + else if (rrr->type == GLP_DB)
1.422 + row->lb = rrr->lb * rii, row->ub = rrr->ub * rii;
1.423 + else if (rrr->type == GLP_FX)
1.424 + row->lb = row->ub = rrr->lb * rii;
1.425 + else
1.426 + xassert(rrr != rrr);
1.427 + }
1.428 + }
1.429 + /* load columns and constraint coefficients */
1.430 + for (j = 1; j <= n; j++)
1.431 + { GLPCOL *ccc = orig->col[j];
1.432 + GLPAIJ *aaa;
1.433 + NPPCOL *col;
1.434 + col = npp_add_col(npp);
1.435 + xassert(col->j == j);
1.436 + if (names && ccc->name != NULL)
1.437 + { col->name = dmp_get_atom(npp->pool, strlen(ccc->name)+1);
1.438 + strcpy(col->name, ccc->name);
1.439 + }
1.440 + if (sol == GLP_MIP)
1.441 +#if 0
1.442 + col->kind = ccc->kind;
1.443 +#else
1.444 + col->is_int = (char)(ccc->kind == GLP_IV);
1.445 +#endif
1.446 + if (!scaling)
1.447 + { if (ccc->type == GLP_FR)
1.448 + col->lb = -DBL_MAX, col->ub = +DBL_MAX;
1.449 + else if (ccc->type == GLP_LO)
1.450 + col->lb = ccc->lb, col->ub = +DBL_MAX;
1.451 + else if (ccc->type == GLP_UP)
1.452 + col->lb = -DBL_MAX, col->ub = ccc->ub;
1.453 + else if (ccc->type == GLP_DB)
1.454 + col->lb = ccc->lb, col->ub = ccc->ub;
1.455 + else if (ccc->type == GLP_FX)
1.456 + col->lb = col->ub = ccc->lb;
1.457 + else
1.458 + xassert(ccc != ccc);
1.459 + col->coef = dir * ccc->coef;
1.460 + for (aaa = ccc->ptr; aaa != NULL; aaa = aaa->c_next)
1.461 + npp_add_aij(npp, link[aaa->row->i], col, aaa->val);
1.462 + }
1.463 + else
1.464 + { double sjj = ccc->sjj;
1.465 + if (ccc->type == GLP_FR)
1.466 + col->lb = -DBL_MAX, col->ub = +DBL_MAX;
1.467 + else if (ccc->type == GLP_LO)
1.468 + col->lb = ccc->lb / sjj, col->ub = +DBL_MAX;
1.469 + else if (ccc->type == GLP_UP)
1.470 + col->lb = -DBL_MAX, col->ub = ccc->ub / sjj;
1.471 + else if (ccc->type == GLP_DB)
1.472 + col->lb = ccc->lb / sjj, col->ub = ccc->ub / sjj;
1.473 + else if (ccc->type == GLP_FX)
1.474 + col->lb = col->ub = ccc->lb / sjj;
1.475 + else
1.476 + xassert(ccc != ccc);
1.477 + col->coef = dir * ccc->coef * sjj;
1.478 + for (aaa = ccc->ptr; aaa != NULL; aaa = aaa->c_next)
1.479 + npp_add_aij(npp, link[aaa->row->i], col,
1.480 + aaa->row->rii * aaa->val * sjj);
1.481 + }
1.482 + }
1.483 + xfree(link);
1.484 + /* keep solution indicator and scaling option */
1.485 + npp->sol = sol;
1.486 + npp->scaling = scaling;
1.487 + return;
1.488 +}
1.489 +
1.490 +void npp_build_prob(NPP *npp, glp_prob *prob)
1.491 +{ /* build resultant (preprocessed) problem */
1.492 + NPPROW *row;
1.493 + NPPCOL *col;
1.494 + NPPAIJ *aij;
1.495 + int i, j, type, len, *ind;
1.496 + double dir, *val;
1.497 + glp_erase_prob(prob);
1.498 + glp_set_prob_name(prob, npp->name);
1.499 + glp_set_obj_name(prob, npp->obj);
1.500 + glp_set_obj_dir(prob, npp->orig_dir);
1.501 + if (npp->orig_dir == GLP_MIN)
1.502 + dir = +1.0;
1.503 + else if (npp->orig_dir == GLP_MAX)
1.504 + dir = -1.0;
1.505 + else
1.506 + xassert(npp != npp);
1.507 + glp_set_obj_coef(prob, 0, dir * npp->c0);
1.508 + /* build rows */
1.509 + for (row = npp->r_head; row != NULL; row = row->next)
1.510 + { row->temp = i = glp_add_rows(prob, 1);
1.511 + glp_set_row_name(prob, i, row->name);
1.512 + if (row->lb == -DBL_MAX && row->ub == +DBL_MAX)
1.513 + type = GLP_FR;
1.514 + else if (row->ub == +DBL_MAX)
1.515 + type = GLP_LO;
1.516 + else if (row->lb == -DBL_MAX)
1.517 + type = GLP_UP;
1.518 + else if (row->lb != row->ub)
1.519 + type = GLP_DB;
1.520 + else
1.521 + type = GLP_FX;
1.522 + glp_set_row_bnds(prob, i, type, row->lb, row->ub);
1.523 + }
1.524 + /* build columns and the constraint matrix */
1.525 + ind = xcalloc(1+prob->m, sizeof(int));
1.526 + val = xcalloc(1+prob->m, sizeof(double));
1.527 + for (col = npp->c_head; col != NULL; col = col->next)
1.528 + { j = glp_add_cols(prob, 1);
1.529 + glp_set_col_name(prob, j, col->name);
1.530 +#if 0
1.531 + glp_set_col_kind(prob, j, col->kind);
1.532 +#else
1.533 + glp_set_col_kind(prob, j, col->is_int ? GLP_IV : GLP_CV);
1.534 +#endif
1.535 + if (col->lb == -DBL_MAX && col->ub == +DBL_MAX)
1.536 + type = GLP_FR;
1.537 + else if (col->ub == +DBL_MAX)
1.538 + type = GLP_LO;
1.539 + else if (col->lb == -DBL_MAX)
1.540 + type = GLP_UP;
1.541 + else if (col->lb != col->ub)
1.542 + type = GLP_DB;
1.543 + else
1.544 + type = GLP_FX;
1.545 + glp_set_col_bnds(prob, j, type, col->lb, col->ub);
1.546 + glp_set_obj_coef(prob, j, dir * col->coef);
1.547 + len = 0;
1.548 + for (aij = col->ptr; aij != NULL; aij = aij->c_next)
1.549 + { len++;
1.550 + ind[len] = aij->row->temp;
1.551 + val[len] = aij->val;
1.552 + }
1.553 + glp_set_mat_col(prob, j, len, ind, val);
1.554 + }
1.555 + xfree(ind);
1.556 + xfree(val);
1.557 + /* resultant problem has been built */
1.558 + npp->m = prob->m;
1.559 + npp->n = prob->n;
1.560 + npp->nnz = prob->nnz;
1.561 + npp->row_ref = xcalloc(1+npp->m, sizeof(int));
1.562 + npp->col_ref = xcalloc(1+npp->n, sizeof(int));
1.563 + for (row = npp->r_head, i = 0; row != NULL; row = row->next)
1.564 + npp->row_ref[++i] = row->i;
1.565 + for (col = npp->c_head, j = 0; col != NULL; col = col->next)
1.566 + npp->col_ref[++j] = col->j;
1.567 + /* transformed problem segment is no longer needed */
1.568 + dmp_delete_pool(npp->pool), npp->pool = NULL;
1.569 + npp->name = npp->obj = NULL;
1.570 + npp->c0 = 0.0;
1.571 + npp->r_head = npp->r_tail = NULL;
1.572 + npp->c_head = npp->c_tail = NULL;
1.573 + return;
1.574 +}
1.575 +
1.576 +void npp_postprocess(NPP *npp, glp_prob *prob)
1.577 +{ /* postprocess solution from the resultant problem */
1.578 + GLPROW *row;
1.579 + GLPCOL *col;
1.580 + NPPTSE *tse;
1.581 + int i, j, k;
1.582 + double dir;
1.583 + xassert(npp->orig_dir == prob->dir);
1.584 + if (npp->orig_dir == GLP_MIN)
1.585 + dir = +1.0;
1.586 + else if (npp->orig_dir == GLP_MAX)
1.587 + dir = -1.0;
1.588 + else
1.589 + xassert(npp != npp);
1.590 + xassert(npp->m == prob->m);
1.591 + xassert(npp->n == prob->n);
1.592 + xassert(npp->nnz == prob->nnz);
1.593 + /* copy solution status */
1.594 + if (npp->sol == GLP_SOL)
1.595 + { npp->p_stat = prob->pbs_stat;
1.596 + npp->d_stat = prob->dbs_stat;
1.597 + }
1.598 + else if (npp->sol == GLP_IPT)
1.599 + npp->t_stat = prob->ipt_stat;
1.600 + else if (npp->sol == GLP_MIP)
1.601 + npp->i_stat = prob->mip_stat;
1.602 + else
1.603 + xassert(npp != npp);
1.604 + /* allocate solution arrays */
1.605 + if (npp->sol == GLP_SOL)
1.606 + { if (npp->r_stat == NULL)
1.607 + npp->r_stat = xcalloc(1+npp->nrows, sizeof(char));
1.608 + for (i = 1; i <= npp->nrows; i++)
1.609 + npp->r_stat[i] = 0;
1.610 + if (npp->c_stat == NULL)
1.611 + npp->c_stat = xcalloc(1+npp->ncols, sizeof(char));
1.612 + for (j = 1; j <= npp->ncols; j++)
1.613 + npp->c_stat[j] = 0;
1.614 + }
1.615 +#if 0
1.616 + if (npp->r_prim == NULL)
1.617 + npp->r_prim = xcalloc(1+npp->nrows, sizeof(double));
1.618 + for (i = 1; i <= npp->nrows; i++)
1.619 + npp->r_prim[i] = DBL_MAX;
1.620 +#endif
1.621 + if (npp->c_value == NULL)
1.622 + npp->c_value = xcalloc(1+npp->ncols, sizeof(double));
1.623 + for (j = 1; j <= npp->ncols; j++)
1.624 + npp->c_value[j] = DBL_MAX;
1.625 + if (npp->sol != GLP_MIP)
1.626 + { if (npp->r_pi == NULL)
1.627 + npp->r_pi = xcalloc(1+npp->nrows, sizeof(double));
1.628 + for (i = 1; i <= npp->nrows; i++)
1.629 + npp->r_pi[i] = DBL_MAX;
1.630 +#if 0
1.631 + if (npp->c_dual == NULL)
1.632 + npp->c_dual = xcalloc(1+npp->ncols, sizeof(double));
1.633 + for (j = 1; j <= npp->ncols; j++)
1.634 + npp->c_dual[j] = DBL_MAX;
1.635 +#endif
1.636 + }
1.637 + /* copy solution components from the resultant problem */
1.638 + if (npp->sol == GLP_SOL)
1.639 + { for (i = 1; i <= npp->m; i++)
1.640 + { row = prob->row[i];
1.641 + k = npp->row_ref[i];
1.642 + npp->r_stat[k] = (char)row->stat;
1.643 + /*npp->r_prim[k] = row->prim;*/
1.644 + npp->r_pi[k] = dir * row->dual;
1.645 + }
1.646 + for (j = 1; j <= npp->n; j++)
1.647 + { col = prob->col[j];
1.648 + k = npp->col_ref[j];
1.649 + npp->c_stat[k] = (char)col->stat;
1.650 + npp->c_value[k] = col->prim;
1.651 + /*npp->c_dual[k] = dir * col->dual;*/
1.652 + }
1.653 + }
1.654 + else if (npp->sol == GLP_IPT)
1.655 + { for (i = 1; i <= npp->m; i++)
1.656 + { row = prob->row[i];
1.657 + k = npp->row_ref[i];
1.658 + /*npp->r_prim[k] = row->pval;*/
1.659 + npp->r_pi[k] = dir * row->dval;
1.660 + }
1.661 + for (j = 1; j <= npp->n; j++)
1.662 + { col = prob->col[j];
1.663 + k = npp->col_ref[j];
1.664 + npp->c_value[k] = col->pval;
1.665 + /*npp->c_dual[k] = dir * col->dval;*/
1.666 + }
1.667 + }
1.668 + else if (npp->sol == GLP_MIP)
1.669 + {
1.670 +#if 0
1.671 + for (i = 1; i <= npp->m; i++)
1.672 + { row = prob->row[i];
1.673 + k = npp->row_ref[i];
1.674 + /*npp->r_prim[k] = row->mipx;*/
1.675 + }
1.676 +#endif
1.677 + for (j = 1; j <= npp->n; j++)
1.678 + { col = prob->col[j];
1.679 + k = npp->col_ref[j];
1.680 + npp->c_value[k] = col->mipx;
1.681 + }
1.682 + }
1.683 + else
1.684 + xassert(npp != npp);
1.685 + /* perform postprocessing to construct solution to the original
1.686 + problem */
1.687 + for (tse = npp->top; tse != NULL; tse = tse->link)
1.688 + { xassert(tse->func != NULL);
1.689 + xassert(tse->func(npp, tse->info) == 0);
1.690 + }
1.691 + return;
1.692 +}
1.693 +
1.694 +void npp_unload_sol(NPP *npp, glp_prob *orig)
1.695 +{ /* store solution to the original problem */
1.696 + GLPROW *row;
1.697 + GLPCOL *col;
1.698 + int i, j;
1.699 + double dir;
1.700 + xassert(npp->orig_dir == orig->dir);
1.701 + if (npp->orig_dir == GLP_MIN)
1.702 + dir = +1.0;
1.703 + else if (npp->orig_dir == GLP_MAX)
1.704 + dir = -1.0;
1.705 + else
1.706 + xassert(npp != npp);
1.707 + xassert(npp->orig_m == orig->m);
1.708 + xassert(npp->orig_n == orig->n);
1.709 + xassert(npp->orig_nnz == orig->nnz);
1.710 + if (npp->sol == GLP_SOL)
1.711 + { /* store basic solution */
1.712 + orig->valid = 0;
1.713 + orig->pbs_stat = npp->p_stat;
1.714 + orig->dbs_stat = npp->d_stat;
1.715 + orig->obj_val = orig->c0;
1.716 + orig->some = 0;
1.717 + for (i = 1; i <= orig->m; i++)
1.718 + { row = orig->row[i];
1.719 + row->stat = npp->r_stat[i];
1.720 + if (!npp->scaling)
1.721 + { /*row->prim = npp->r_prim[i];*/
1.722 + row->dual = dir * npp->r_pi[i];
1.723 + }
1.724 + else
1.725 + { /*row->prim = npp->r_prim[i] / row->rii;*/
1.726 + row->dual = dir * npp->r_pi[i] * row->rii;
1.727 + }
1.728 + if (row->stat == GLP_BS)
1.729 + row->dual = 0.0;
1.730 + else if (row->stat == GLP_NL)
1.731 + { xassert(row->type == GLP_LO || row->type == GLP_DB);
1.732 + row->prim = row->lb;
1.733 + }
1.734 + else if (row->stat == GLP_NU)
1.735 + { xassert(row->type == GLP_UP || row->type == GLP_DB);
1.736 + row->prim = row->ub;
1.737 + }
1.738 + else if (row->stat == GLP_NF)
1.739 + { xassert(row->type == GLP_FR);
1.740 + row->prim = 0.0;
1.741 + }
1.742 + else if (row->stat == GLP_NS)
1.743 + { xassert(row->type == GLP_FX);
1.744 + row->prim = row->lb;
1.745 + }
1.746 + else
1.747 + xassert(row != row);
1.748 + }
1.749 + for (j = 1; j <= orig->n; j++)
1.750 + { col = orig->col[j];
1.751 + col->stat = npp->c_stat[j];
1.752 + if (!npp->scaling)
1.753 + { col->prim = npp->c_value[j];
1.754 + /*col->dual = dir * npp->c_dual[j];*/
1.755 + }
1.756 + else
1.757 + { col->prim = npp->c_value[j] * col->sjj;
1.758 + /*col->dual = dir * npp->c_dual[j] / col->sjj;*/
1.759 + }
1.760 + if (col->stat == GLP_BS)
1.761 + col->dual = 0.0;
1.762 +#if 1
1.763 + else if (col->stat == GLP_NL)
1.764 + { xassert(col->type == GLP_LO || col->type == GLP_DB);
1.765 + col->prim = col->lb;
1.766 + }
1.767 + else if (col->stat == GLP_NU)
1.768 + { xassert(col->type == GLP_UP || col->type == GLP_DB);
1.769 + col->prim = col->ub;
1.770 + }
1.771 + else if (col->stat == GLP_NF)
1.772 + { xassert(col->type == GLP_FR);
1.773 + col->prim = 0.0;
1.774 + }
1.775 + else if (col->stat == GLP_NS)
1.776 + { xassert(col->type == GLP_FX);
1.777 + col->prim = col->lb;
1.778 + }
1.779 + else
1.780 + xassert(col != col);
1.781 +#endif
1.782 + orig->obj_val += col->coef * col->prim;
1.783 + }
1.784 +#if 1
1.785 + /* compute primal values of inactive rows */
1.786 + for (i = 1; i <= orig->m; i++)
1.787 + { row = orig->row[i];
1.788 + if (row->stat == GLP_BS)
1.789 + { GLPAIJ *aij;
1.790 + double temp;
1.791 + temp = 0.0;
1.792 + for (aij = row->ptr; aij != NULL; aij = aij->r_next)
1.793 + temp += aij->val * aij->col->prim;
1.794 + row->prim = temp;
1.795 + }
1.796 + }
1.797 + /* compute reduced costs of active columns */
1.798 + for (j = 1; j <= orig->n; j++)
1.799 + { col = orig->col[j];
1.800 + if (col->stat != GLP_BS)
1.801 + { GLPAIJ *aij;
1.802 + double temp;
1.803 + temp = col->coef;
1.804 + for (aij = col->ptr; aij != NULL; aij = aij->c_next)
1.805 + temp -= aij->val * aij->row->dual;
1.806 + col->dual = temp;
1.807 + }
1.808 + }
1.809 +#endif
1.810 + }
1.811 + else if (npp->sol == GLP_IPT)
1.812 + { /* store interior-point solution */
1.813 + orig->ipt_stat = npp->t_stat;
1.814 + orig->ipt_obj = orig->c0;
1.815 + for (i = 1; i <= orig->m; i++)
1.816 + { row = orig->row[i];
1.817 + if (!npp->scaling)
1.818 + { /*row->pval = npp->r_prim[i];*/
1.819 + row->dval = dir * npp->r_pi[i];
1.820 + }
1.821 + else
1.822 + { /*row->pval = npp->r_prim[i] / row->rii;*/
1.823 + row->dval = dir * npp->r_pi[i] * row->rii;
1.824 + }
1.825 + }
1.826 + for (j = 1; j <= orig->n; j++)
1.827 + { col = orig->col[j];
1.828 + if (!npp->scaling)
1.829 + { col->pval = npp->c_value[j];
1.830 + /*col->dval = dir * npp->c_dual[j];*/
1.831 + }
1.832 + else
1.833 + { col->pval = npp->c_value[j] * col->sjj;
1.834 + /*col->dval = dir * npp->c_dual[j] / col->sjj;*/
1.835 + }
1.836 + orig->ipt_obj += col->coef * col->pval;
1.837 + }
1.838 +#if 1
1.839 + /* compute row primal values */
1.840 + for (i = 1; i <= orig->m; i++)
1.841 + { row = orig->row[i];
1.842 + { GLPAIJ *aij;
1.843 + double temp;
1.844 + temp = 0.0;
1.845 + for (aij = row->ptr; aij != NULL; aij = aij->r_next)
1.846 + temp += aij->val * aij->col->pval;
1.847 + row->pval = temp;
1.848 + }
1.849 + }
1.850 + /* compute column dual values */
1.851 + for (j = 1; j <= orig->n; j++)
1.852 + { col = orig->col[j];
1.853 + { GLPAIJ *aij;
1.854 + double temp;
1.855 + temp = col->coef;
1.856 + for (aij = col->ptr; aij != NULL; aij = aij->c_next)
1.857 + temp -= aij->val * aij->row->dval;
1.858 + col->dval = temp;
1.859 + }
1.860 + }
1.861 +#endif
1.862 + }
1.863 + else if (npp->sol == GLP_MIP)
1.864 + { /* store MIP solution */
1.865 + xassert(!npp->scaling);
1.866 + orig->mip_stat = npp->i_stat;
1.867 + orig->mip_obj = orig->c0;
1.868 +#if 0
1.869 + for (i = 1; i <= orig->m; i++)
1.870 + { row = orig->row[i];
1.871 + /*row->mipx = npp->r_prim[i];*/
1.872 + }
1.873 +#endif
1.874 + for (j = 1; j <= orig->n; j++)
1.875 + { col = orig->col[j];
1.876 + col->mipx = npp->c_value[j];
1.877 + if (col->kind == GLP_IV)
1.878 + xassert(col->mipx == floor(col->mipx));
1.879 + orig->mip_obj += col->coef * col->mipx;
1.880 + }
1.881 +#if 1
1.882 + /* compute row primal values */
1.883 + for (i = 1; i <= orig->m; i++)
1.884 + { row = orig->row[i];
1.885 + { GLPAIJ *aij;
1.886 + double temp;
1.887 + temp = 0.0;
1.888 + for (aij = row->ptr; aij != NULL; aij = aij->r_next)
1.889 + temp += aij->val * aij->col->mipx;
1.890 + row->mipx = temp;
1.891 + }
1.892 + }
1.893 +#endif
1.894 + }
1.895 + else
1.896 + xassert(npp != npp);
1.897 + return;
1.898 +}
1.899 +
1.900 +void npp_delete_wksp(NPP *npp)
1.901 +{ /* delete LP/MIP preprocessor workspace */
1.902 + if (npp->pool != NULL)
1.903 + dmp_delete_pool(npp->pool);
1.904 + if (npp->stack != NULL)
1.905 + dmp_delete_pool(npp->stack);
1.906 + if (npp->row_ref != NULL)
1.907 + xfree(npp->row_ref);
1.908 + if (npp->col_ref != NULL)
1.909 + xfree(npp->col_ref);
1.910 + if (npp->r_stat != NULL)
1.911 + xfree(npp->r_stat);
1.912 +#if 0
1.913 + if (npp->r_prim != NULL)
1.914 + xfree(npp->r_prim);
1.915 +#endif
1.916 + if (npp->r_pi != NULL)
1.917 + xfree(npp->r_pi);
1.918 + if (npp->c_stat != NULL)
1.919 + xfree(npp->c_stat);
1.920 + if (npp->c_value != NULL)
1.921 + xfree(npp->c_value);
1.922 +#if 0
1.923 + if (npp->c_dual != NULL)
1.924 + xfree(npp->c_dual);
1.925 +#endif
1.926 + xfree(npp);
1.927 + return;
1.928 +}
1.929 +
1.930 +/* eof */