lemon-project-template-glpk

annotate deps/glpk/src/glpnpp01.c @ 9:33de93886c88

find changesets by author, revision, files, or words in the commit message
Import GLPK 4.47
author Alpar Juttner <alpar@cs.elte.hu>
date Sun, 06 Nov 2011 20:59:10 +0100
parents
children
rev   line source
alpar@9 1 /* glpnpp01.c */
alpar@9 2
alpar@9 3 /***********************************************************************
alpar@9 4 * This code is part of GLPK (GNU Linear Programming Kit).
alpar@9 5 *
alpar@9 6 * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
alpar@9 7 * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
alpar@9 8 * Moscow Aviation Institute, Moscow, Russia. All rights reserved.
alpar@9 9 * E-mail: <mao@gnu.org>.
alpar@9 10 *
alpar@9 11 * GLPK is free software: you can redistribute it and/or modify it
alpar@9 12 * under the terms of the GNU General Public License as published by
alpar@9 13 * the Free Software Foundation, either version 3 of the License, or
alpar@9 14 * (at your option) any later version.
alpar@9 15 *
alpar@9 16 * GLPK is distributed in the hope that it will be useful, but WITHOUT
alpar@9 17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
alpar@9 18 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
alpar@9 19 * License for more details.
alpar@9 20 *
alpar@9 21 * You should have received a copy of the GNU General Public License
alpar@9 22 * along with GLPK. If not, see <http://www.gnu.org/licenses/>.
alpar@9 23 ***********************************************************************/
alpar@9 24
alpar@9 25 #include "glpnpp.h"
alpar@9 26
alpar@9 27 NPP *npp_create_wksp(void)
alpar@9 28 { /* create LP/MIP preprocessor workspace */
alpar@9 29 NPP *npp;
alpar@9 30 npp = xmalloc(sizeof(NPP));
alpar@9 31 npp->orig_dir = 0;
alpar@9 32 npp->orig_m = npp->orig_n = npp->orig_nnz = 0;
alpar@9 33 npp->pool = dmp_create_pool();
alpar@9 34 npp->name = npp->obj = NULL;
alpar@9 35 npp->c0 = 0.0;
alpar@9 36 npp->nrows = npp->ncols = 0;
alpar@9 37 npp->r_head = npp->r_tail = NULL;
alpar@9 38 npp->c_head = npp->c_tail = NULL;
alpar@9 39 npp->stack = dmp_create_pool();
alpar@9 40 npp->top = NULL;
alpar@9 41 #if 0 /* 16/XII-2009 */
alpar@9 42 memset(&npp->count, 0, sizeof(npp->count));
alpar@9 43 #endif
alpar@9 44 npp->m = npp->n = npp->nnz = 0;
alpar@9 45 npp->row_ref = npp->col_ref = NULL;
alpar@9 46 npp->sol = npp->scaling = 0;
alpar@9 47 npp->p_stat = npp->d_stat = npp->t_stat = npp->i_stat = 0;
alpar@9 48 npp->r_stat = NULL;
alpar@9 49 /*npp->r_prim =*/ npp->r_pi = NULL;
alpar@9 50 npp->c_stat = NULL;
alpar@9 51 npp->c_value = /*npp->c_dual =*/ NULL;
alpar@9 52 return npp;
alpar@9 53 }
alpar@9 54
alpar@9 55 void npp_insert_row(NPP *npp, NPPROW *row, int where)
alpar@9 56 { /* insert row to the row list */
alpar@9 57 if (where == 0)
alpar@9 58 { /* insert row to the beginning of the row list */
alpar@9 59 row->prev = NULL;
alpar@9 60 row->next = npp->r_head;
alpar@9 61 if (row->next == NULL)
alpar@9 62 npp->r_tail = row;
alpar@9 63 else
alpar@9 64 row->next->prev = row;
alpar@9 65 npp->r_head = row;
alpar@9 66 }
alpar@9 67 else
alpar@9 68 { /* insert row to the end of the row list */
alpar@9 69 row->prev = npp->r_tail;
alpar@9 70 row->next = NULL;
alpar@9 71 if (row->prev == NULL)
alpar@9 72 npp->r_head = row;
alpar@9 73 else
alpar@9 74 row->prev->next = row;
alpar@9 75 npp->r_tail = row;
alpar@9 76 }
alpar@9 77 return;
alpar@9 78 }
alpar@9 79
alpar@9 80 void npp_remove_row(NPP *npp, NPPROW *row)
alpar@9 81 { /* remove row from the row list */
alpar@9 82 if (row->prev == NULL)
alpar@9 83 npp->r_head = row->next;
alpar@9 84 else
alpar@9 85 row->prev->next = row->next;
alpar@9 86 if (row->next == NULL)
alpar@9 87 npp->r_tail = row->prev;
alpar@9 88 else
alpar@9 89 row->next->prev = row->prev;
alpar@9 90 return;
alpar@9 91 }
alpar@9 92
alpar@9 93 void npp_activate_row(NPP *npp, NPPROW *row)
alpar@9 94 { /* make row active */
alpar@9 95 if (!row->temp)
alpar@9 96 { row->temp = 1;
alpar@9 97 /* move the row to the beginning of the row list */
alpar@9 98 npp_remove_row(npp, row);
alpar@9 99 npp_insert_row(npp, row, 0);
alpar@9 100 }
alpar@9 101 return;
alpar@9 102 }
alpar@9 103
alpar@9 104 void npp_deactivate_row(NPP *npp, NPPROW *row)
alpar@9 105 { /* make row inactive */
alpar@9 106 if (row->temp)
alpar@9 107 { row->temp = 0;
alpar@9 108 /* move the row to the end of the row list */
alpar@9 109 npp_remove_row(npp, row);
alpar@9 110 npp_insert_row(npp, row, 1);
alpar@9 111 }
alpar@9 112 return;
alpar@9 113 }
alpar@9 114
alpar@9 115 void npp_insert_col(NPP *npp, NPPCOL *col, int where)
alpar@9 116 { /* insert column to the column list */
alpar@9 117 if (where == 0)
alpar@9 118 { /* insert column to the beginning of the column list */
alpar@9 119 col->prev = NULL;
alpar@9 120 col->next = npp->c_head;
alpar@9 121 if (col->next == NULL)
alpar@9 122 npp->c_tail = col;
alpar@9 123 else
alpar@9 124 col->next->prev = col;
alpar@9 125 npp->c_head = col;
alpar@9 126 }
alpar@9 127 else
alpar@9 128 { /* insert column to the end of the column list */
alpar@9 129 col->prev = npp->c_tail;
alpar@9 130 col->next = NULL;
alpar@9 131 if (col->prev == NULL)
alpar@9 132 npp->c_head = col;
alpar@9 133 else
alpar@9 134 col->prev->next = col;
alpar@9 135 npp->c_tail = col;
alpar@9 136 }
alpar@9 137 return;
alpar@9 138 }
alpar@9 139
alpar@9 140 void npp_remove_col(NPP *npp, NPPCOL *col)
alpar@9 141 { /* remove column from the column list */
alpar@9 142 if (col->prev == NULL)
alpar@9 143 npp->c_head = col->next;
alpar@9 144 else
alpar@9 145 col->prev->next = col->next;
alpar@9 146 if (col->next == NULL)
alpar@9 147 npp->c_tail = col->prev;
alpar@9 148 else
alpar@9 149 col->next->prev = col->prev;
alpar@9 150 return;
alpar@9 151 }
alpar@9 152
alpar@9 153 void npp_activate_col(NPP *npp, NPPCOL *col)
alpar@9 154 { /* make column active */
alpar@9 155 if (!col->temp)
alpar@9 156 { col->temp = 1;
alpar@9 157 /* move the column to the beginning of the column list */
alpar@9 158 npp_remove_col(npp, col);
alpar@9 159 npp_insert_col(npp, col, 0);
alpar@9 160 }
alpar@9 161 return;
alpar@9 162 }
alpar@9 163
alpar@9 164 void npp_deactivate_col(NPP *npp, NPPCOL *col)
alpar@9 165 { /* make column inactive */
alpar@9 166 if (col->temp)
alpar@9 167 { col->temp = 0;
alpar@9 168 /* move the column to the end of the column list */
alpar@9 169 npp_remove_col(npp, col);
alpar@9 170 npp_insert_col(npp, col, 1);
alpar@9 171 }
alpar@9 172 return;
alpar@9 173 }
alpar@9 174
alpar@9 175 NPPROW *npp_add_row(NPP *npp)
alpar@9 176 { /* add new row to the current problem */
alpar@9 177 NPPROW *row;
alpar@9 178 row = dmp_get_atom(npp->pool, sizeof(NPPROW));
alpar@9 179 row->i = ++(npp->nrows);
alpar@9 180 row->name = NULL;
alpar@9 181 row->lb = -DBL_MAX, row->ub = +DBL_MAX;
alpar@9 182 row->ptr = NULL;
alpar@9 183 row->temp = 0;
alpar@9 184 npp_insert_row(npp, row, 1);
alpar@9 185 return row;
alpar@9 186 }
alpar@9 187
alpar@9 188 NPPCOL *npp_add_col(NPP *npp)
alpar@9 189 { /* add new column to the current problem */
alpar@9 190 NPPCOL *col;
alpar@9 191 col = dmp_get_atom(npp->pool, sizeof(NPPCOL));
alpar@9 192 col->j = ++(npp->ncols);
alpar@9 193 col->name = NULL;
alpar@9 194 #if 0
alpar@9 195 col->kind = GLP_CV;
alpar@9 196 #else
alpar@9 197 col->is_int = 0;
alpar@9 198 #endif
alpar@9 199 col->lb = col->ub = col->coef = 0.0;
alpar@9 200 col->ptr = NULL;
alpar@9 201 col->temp = 0;
alpar@9 202 npp_insert_col(npp, col, 1);
alpar@9 203 return col;
alpar@9 204 }
alpar@9 205
alpar@9 206 NPPAIJ *npp_add_aij(NPP *npp, NPPROW *row, NPPCOL *col, double val)
alpar@9 207 { /* add new element to the constraint matrix */
alpar@9 208 NPPAIJ *aij;
alpar@9 209 aij = dmp_get_atom(npp->pool, sizeof(NPPAIJ));
alpar@9 210 aij->row = row;
alpar@9 211 aij->col = col;
alpar@9 212 aij->val = val;
alpar@9 213 aij->r_prev = NULL;
alpar@9 214 aij->r_next = row->ptr;
alpar@9 215 aij->c_prev = NULL;
alpar@9 216 aij->c_next = col->ptr;
alpar@9 217 if (aij->r_next != NULL)
alpar@9 218 aij->r_next->r_prev = aij;
alpar@9 219 if (aij->c_next != NULL)
alpar@9 220 aij->c_next->c_prev = aij;
alpar@9 221 row->ptr = col->ptr = aij;
alpar@9 222 return aij;
alpar@9 223 }
alpar@9 224
alpar@9 225 int npp_row_nnz(NPP *npp, NPPROW *row)
alpar@9 226 { /* count number of non-zero coefficients in row */
alpar@9 227 NPPAIJ *aij;
alpar@9 228 int nnz;
alpar@9 229 xassert(npp == npp);
alpar@9 230 nnz = 0;
alpar@9 231 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
alpar@9 232 nnz++;
alpar@9 233 return nnz;
alpar@9 234 }
alpar@9 235
alpar@9 236 int npp_col_nnz(NPP *npp, NPPCOL *col)
alpar@9 237 { /* count number of non-zero coefficients in column */
alpar@9 238 NPPAIJ *aij;
alpar@9 239 int nnz;
alpar@9 240 xassert(npp == npp);
alpar@9 241 nnz = 0;
alpar@9 242 for (aij = col->ptr; aij != NULL; aij = aij->c_next)
alpar@9 243 nnz++;
alpar@9 244 return nnz;
alpar@9 245 }
alpar@9 246
alpar@9 247 void *npp_push_tse(NPP *npp, int (*func)(NPP *npp, void *info),
alpar@9 248 int size)
alpar@9 249 { /* push new entry to the transformation stack */
alpar@9 250 NPPTSE *tse;
alpar@9 251 tse = dmp_get_atom(npp->stack, sizeof(NPPTSE));
alpar@9 252 tse->func = func;
alpar@9 253 tse->info = dmp_get_atom(npp->stack, size);
alpar@9 254 tse->link = npp->top;
alpar@9 255 npp->top = tse;
alpar@9 256 return tse->info;
alpar@9 257 }
alpar@9 258
alpar@9 259 #if 1 /* 23/XII-2009 */
alpar@9 260 void npp_erase_row(NPP *npp, NPPROW *row)
alpar@9 261 { /* erase row content to make it empty */
alpar@9 262 NPPAIJ *aij;
alpar@9 263 while (row->ptr != NULL)
alpar@9 264 { aij = row->ptr;
alpar@9 265 row->ptr = aij->r_next;
alpar@9 266 if (aij->c_prev == NULL)
alpar@9 267 aij->col->ptr = aij->c_next;
alpar@9 268 else
alpar@9 269 aij->c_prev->c_next = aij->c_next;
alpar@9 270 if (aij->c_next == NULL)
alpar@9 271 ;
alpar@9 272 else
alpar@9 273 aij->c_next->c_prev = aij->c_prev;
alpar@9 274 dmp_free_atom(npp->pool, aij, sizeof(NPPAIJ));
alpar@9 275 }
alpar@9 276 return;
alpar@9 277 }
alpar@9 278 #endif
alpar@9 279
alpar@9 280 void npp_del_row(NPP *npp, NPPROW *row)
alpar@9 281 { /* remove row from the current problem */
alpar@9 282 #if 0 /* 23/XII-2009 */
alpar@9 283 NPPAIJ *aij;
alpar@9 284 #endif
alpar@9 285 if (row->name != NULL)
alpar@9 286 dmp_free_atom(npp->pool, row->name, strlen(row->name)+1);
alpar@9 287 #if 0 /* 23/XII-2009 */
alpar@9 288 while (row->ptr != NULL)
alpar@9 289 { aij = row->ptr;
alpar@9 290 row->ptr = aij->r_next;
alpar@9 291 if (aij->c_prev == NULL)
alpar@9 292 aij->col->ptr = aij->c_next;
alpar@9 293 else
alpar@9 294 aij->c_prev->c_next = aij->c_next;
alpar@9 295 if (aij->c_next == NULL)
alpar@9 296 ;
alpar@9 297 else
alpar@9 298 aij->c_next->c_prev = aij->c_prev;
alpar@9 299 dmp_free_atom(npp->pool, aij, sizeof(NPPAIJ));
alpar@9 300 }
alpar@9 301 #else
alpar@9 302 npp_erase_row(npp, row);
alpar@9 303 #endif
alpar@9 304 npp_remove_row(npp, row);
alpar@9 305 dmp_free_atom(npp->pool, row, sizeof(NPPROW));
alpar@9 306 return;
alpar@9 307 }
alpar@9 308
alpar@9 309 void npp_del_col(NPP *npp, NPPCOL *col)
alpar@9 310 { /* remove column from the current problem */
alpar@9 311 NPPAIJ *aij;
alpar@9 312 if (col->name != NULL)
alpar@9 313 dmp_free_atom(npp->pool, col->name, strlen(col->name)+1);
alpar@9 314 while (col->ptr != NULL)
alpar@9 315 { aij = col->ptr;
alpar@9 316 col->ptr = aij->c_next;
alpar@9 317 if (aij->r_prev == NULL)
alpar@9 318 aij->row->ptr = aij->r_next;
alpar@9 319 else
alpar@9 320 aij->r_prev->r_next = aij->r_next;
alpar@9 321 if (aij->r_next == NULL)
alpar@9 322 ;
alpar@9 323 else
alpar@9 324 aij->r_next->r_prev = aij->r_prev;
alpar@9 325 dmp_free_atom(npp->pool, aij, sizeof(NPPAIJ));
alpar@9 326 }
alpar@9 327 npp_remove_col(npp, col);
alpar@9 328 dmp_free_atom(npp->pool, col, sizeof(NPPCOL));
alpar@9 329 return;
alpar@9 330 }
alpar@9 331
alpar@9 332 void npp_del_aij(NPP *npp, NPPAIJ *aij)
alpar@9 333 { /* remove element from the constraint matrix */
alpar@9 334 if (aij->r_prev == NULL)
alpar@9 335 aij->row->ptr = aij->r_next;
alpar@9 336 else
alpar@9 337 aij->r_prev->r_next = aij->r_next;
alpar@9 338 if (aij->r_next == NULL)
alpar@9 339 ;
alpar@9 340 else
alpar@9 341 aij->r_next->r_prev = aij->r_prev;
alpar@9 342 if (aij->c_prev == NULL)
alpar@9 343 aij->col->ptr = aij->c_next;
alpar@9 344 else
alpar@9 345 aij->c_prev->c_next = aij->c_next;
alpar@9 346 if (aij->c_next == NULL)
alpar@9 347 ;
alpar@9 348 else
alpar@9 349 aij->c_next->c_prev = aij->c_prev;
alpar@9 350 dmp_free_atom(npp->pool, aij, sizeof(NPPAIJ));
alpar@9 351 return;
alpar@9 352 }
alpar@9 353
alpar@9 354 void npp_load_prob(NPP *npp, glp_prob *orig, int names, int sol,
alpar@9 355 int scaling)
alpar@9 356 { /* load original problem into the preprocessor workspace */
alpar@9 357 int m = orig->m;
alpar@9 358 int n = orig->n;
alpar@9 359 NPPROW **link;
alpar@9 360 int i, j;
alpar@9 361 double dir;
alpar@9 362 xassert(names == GLP_OFF || names == GLP_ON);
alpar@9 363 xassert(sol == GLP_SOL || sol == GLP_IPT || sol == GLP_MIP);
alpar@9 364 xassert(scaling == GLP_OFF || scaling == GLP_ON);
alpar@9 365 if (sol == GLP_MIP) xassert(!scaling);
alpar@9 366 npp->orig_dir = orig->dir;
alpar@9 367 if (npp->orig_dir == GLP_MIN)
alpar@9 368 dir = +1.0;
alpar@9 369 else if (npp->orig_dir == GLP_MAX)
alpar@9 370 dir = -1.0;
alpar@9 371 else
alpar@9 372 xassert(npp != npp);
alpar@9 373 npp->orig_m = m;
alpar@9 374 npp->orig_n = n;
alpar@9 375 npp->orig_nnz = orig->nnz;
alpar@9 376 if (names && orig->name != NULL)
alpar@9 377 { npp->name = dmp_get_atom(npp->pool, strlen(orig->name)+1);
alpar@9 378 strcpy(npp->name, orig->name);
alpar@9 379 }
alpar@9 380 if (names && orig->obj != NULL)
alpar@9 381 { npp->obj = dmp_get_atom(npp->pool, strlen(orig->obj)+1);
alpar@9 382 strcpy(npp->obj, orig->obj);
alpar@9 383 }
alpar@9 384 npp->c0 = dir * orig->c0;
alpar@9 385 /* load rows */
alpar@9 386 link = xcalloc(1+m, sizeof(NPPROW *));
alpar@9 387 for (i = 1; i <= m; i++)
alpar@9 388 { GLPROW *rrr = orig->row[i];
alpar@9 389 NPPROW *row;
alpar@9 390 link[i] = row = npp_add_row(npp);
alpar@9 391 xassert(row->i == i);
alpar@9 392 if (names && rrr->name != NULL)
alpar@9 393 { row->name = dmp_get_atom(npp->pool, strlen(rrr->name)+1);
alpar@9 394 strcpy(row->name, rrr->name);
alpar@9 395 }
alpar@9 396 if (!scaling)
alpar@9 397 { if (rrr->type == GLP_FR)
alpar@9 398 row->lb = -DBL_MAX, row->ub = +DBL_MAX;
alpar@9 399 else if (rrr->type == GLP_LO)
alpar@9 400 row->lb = rrr->lb, row->ub = +DBL_MAX;
alpar@9 401 else if (rrr->type == GLP_UP)
alpar@9 402 row->lb = -DBL_MAX, row->ub = rrr->ub;
alpar@9 403 else if (rrr->type == GLP_DB)
alpar@9 404 row->lb = rrr->lb, row->ub = rrr->ub;
alpar@9 405 else if (rrr->type == GLP_FX)
alpar@9 406 row->lb = row->ub = rrr->lb;
alpar@9 407 else
alpar@9 408 xassert(rrr != rrr);
alpar@9 409 }
alpar@9 410 else
alpar@9 411 { double rii = rrr->rii;
alpar@9 412 if (rrr->type == GLP_FR)
alpar@9 413 row->lb = -DBL_MAX, row->ub = +DBL_MAX;
alpar@9 414 else if (rrr->type == GLP_LO)
alpar@9 415 row->lb = rrr->lb * rii, row->ub = +DBL_MAX;
alpar@9 416 else if (rrr->type == GLP_UP)
alpar@9 417 row->lb = -DBL_MAX, row->ub = rrr->ub * rii;
alpar@9 418 else if (rrr->type == GLP_DB)
alpar@9 419 row->lb = rrr->lb * rii, row->ub = rrr->ub * rii;
alpar@9 420 else if (rrr->type == GLP_FX)
alpar@9 421 row->lb = row->ub = rrr->lb * rii;
alpar@9 422 else
alpar@9 423 xassert(rrr != rrr);
alpar@9 424 }
alpar@9 425 }
alpar@9 426 /* load columns and constraint coefficients */
alpar@9 427 for (j = 1; j <= n; j++)
alpar@9 428 { GLPCOL *ccc = orig->col[j];
alpar@9 429 GLPAIJ *aaa;
alpar@9 430 NPPCOL *col;
alpar@9 431 col = npp_add_col(npp);
alpar@9 432 xassert(col->j == j);
alpar@9 433 if (names && ccc->name != NULL)
alpar@9 434 { col->name = dmp_get_atom(npp->pool, strlen(ccc->name)+1);
alpar@9 435 strcpy(col->name, ccc->name);
alpar@9 436 }
alpar@9 437 if (sol == GLP_MIP)
alpar@9 438 #if 0
alpar@9 439 col->kind = ccc->kind;
alpar@9 440 #else
alpar@9 441 col->is_int = (char)(ccc->kind == GLP_IV);
alpar@9 442 #endif
alpar@9 443 if (!scaling)
alpar@9 444 { if (ccc->type == GLP_FR)
alpar@9 445 col->lb = -DBL_MAX, col->ub = +DBL_MAX;
alpar@9 446 else if (ccc->type == GLP_LO)
alpar@9 447 col->lb = ccc->lb, col->ub = +DBL_MAX;
alpar@9 448 else if (ccc->type == GLP_UP)
alpar@9 449 col->lb = -DBL_MAX, col->ub = ccc->ub;
alpar@9 450 else if (ccc->type == GLP_DB)
alpar@9 451 col->lb = ccc->lb, col->ub = ccc->ub;
alpar@9 452 else if (ccc->type == GLP_FX)
alpar@9 453 col->lb = col->ub = ccc->lb;
alpar@9 454 else
alpar@9 455 xassert(ccc != ccc);
alpar@9 456 col->coef = dir * ccc->coef;
alpar@9 457 for (aaa = ccc->ptr; aaa != NULL; aaa = aaa->c_next)
alpar@9 458 npp_add_aij(npp, link[aaa->row->i], col, aaa->val);
alpar@9 459 }
alpar@9 460 else
alpar@9 461 { double sjj = ccc->sjj;
alpar@9 462 if (ccc->type == GLP_FR)
alpar@9 463 col->lb = -DBL_MAX, col->ub = +DBL_MAX;
alpar@9 464 else if (ccc->type == GLP_LO)
alpar@9 465 col->lb = ccc->lb / sjj, col->ub = +DBL_MAX;
alpar@9 466 else if (ccc->type == GLP_UP)
alpar@9 467 col->lb = -DBL_MAX, col->ub = ccc->ub / sjj;
alpar@9 468 else if (ccc->type == GLP_DB)
alpar@9 469 col->lb = ccc->lb / sjj, col->ub = ccc->ub / sjj;
alpar@9 470 else if (ccc->type == GLP_FX)
alpar@9 471 col->lb = col->ub = ccc->lb / sjj;
alpar@9 472 else
alpar@9 473 xassert(ccc != ccc);
alpar@9 474 col->coef = dir * ccc->coef * sjj;
alpar@9 475 for (aaa = ccc->ptr; aaa != NULL; aaa = aaa->c_next)
alpar@9 476 npp_add_aij(npp, link[aaa->row->i], col,
alpar@9 477 aaa->row->rii * aaa->val * sjj);
alpar@9 478 }
alpar@9 479 }
alpar@9 480 xfree(link);
alpar@9 481 /* keep solution indicator and scaling option */
alpar@9 482 npp->sol = sol;
alpar@9 483 npp->scaling = scaling;
alpar@9 484 return;
alpar@9 485 }
alpar@9 486
alpar@9 487 void npp_build_prob(NPP *npp, glp_prob *prob)
alpar@9 488 { /* build resultant (preprocessed) problem */
alpar@9 489 NPPROW *row;
alpar@9 490 NPPCOL *col;
alpar@9 491 NPPAIJ *aij;
alpar@9 492 int i, j, type, len, *ind;
alpar@9 493 double dir, *val;
alpar@9 494 glp_erase_prob(prob);
alpar@9 495 glp_set_prob_name(prob, npp->name);
alpar@9 496 glp_set_obj_name(prob, npp->obj);
alpar@9 497 glp_set_obj_dir(prob, npp->orig_dir);
alpar@9 498 if (npp->orig_dir == GLP_MIN)
alpar@9 499 dir = +1.0;
alpar@9 500 else if (npp->orig_dir == GLP_MAX)
alpar@9 501 dir = -1.0;
alpar@9 502 else
alpar@9 503 xassert(npp != npp);
alpar@9 504 glp_set_obj_coef(prob, 0, dir * npp->c0);
alpar@9 505 /* build rows */
alpar@9 506 for (row = npp->r_head; row != NULL; row = row->next)
alpar@9 507 { row->temp = i = glp_add_rows(prob, 1);
alpar@9 508 glp_set_row_name(prob, i, row->name);
alpar@9 509 if (row->lb == -DBL_MAX && row->ub == +DBL_MAX)
alpar@9 510 type = GLP_FR;
alpar@9 511 else if (row->ub == +DBL_MAX)
alpar@9 512 type = GLP_LO;
alpar@9 513 else if (row->lb == -DBL_MAX)
alpar@9 514 type = GLP_UP;
alpar@9 515 else if (row->lb != row->ub)
alpar@9 516 type = GLP_DB;
alpar@9 517 else
alpar@9 518 type = GLP_FX;
alpar@9 519 glp_set_row_bnds(prob, i, type, row->lb, row->ub);
alpar@9 520 }
alpar@9 521 /* build columns and the constraint matrix */
alpar@9 522 ind = xcalloc(1+prob->m, sizeof(int));
alpar@9 523 val = xcalloc(1+prob->m, sizeof(double));
alpar@9 524 for (col = npp->c_head; col != NULL; col = col->next)
alpar@9 525 { j = glp_add_cols(prob, 1);
alpar@9 526 glp_set_col_name(prob, j, col->name);
alpar@9 527 #if 0
alpar@9 528 glp_set_col_kind(prob, j, col->kind);
alpar@9 529 #else
alpar@9 530 glp_set_col_kind(prob, j, col->is_int ? GLP_IV : GLP_CV);
alpar@9 531 #endif
alpar@9 532 if (col->lb == -DBL_MAX && col->ub == +DBL_MAX)
alpar@9 533 type = GLP_FR;
alpar@9 534 else if (col->ub == +DBL_MAX)
alpar@9 535 type = GLP_LO;
alpar@9 536 else if (col->lb == -DBL_MAX)
alpar@9 537 type = GLP_UP;
alpar@9 538 else if (col->lb != col->ub)
alpar@9 539 type = GLP_DB;
alpar@9 540 else
alpar@9 541 type = GLP_FX;
alpar@9 542 glp_set_col_bnds(prob, j, type, col->lb, col->ub);
alpar@9 543 glp_set_obj_coef(prob, j, dir * col->coef);
alpar@9 544 len = 0;
alpar@9 545 for (aij = col->ptr; aij != NULL; aij = aij->c_next)
alpar@9 546 { len++;
alpar@9 547 ind[len] = aij->row->temp;
alpar@9 548 val[len] = aij->val;
alpar@9 549 }
alpar@9 550 glp_set_mat_col(prob, j, len, ind, val);
alpar@9 551 }
alpar@9 552 xfree(ind);
alpar@9 553 xfree(val);
alpar@9 554 /* resultant problem has been built */
alpar@9 555 npp->m = prob->m;
alpar@9 556 npp->n = prob->n;
alpar@9 557 npp->nnz = prob->nnz;
alpar@9 558 npp->row_ref = xcalloc(1+npp->m, sizeof(int));
alpar@9 559 npp->col_ref = xcalloc(1+npp->n, sizeof(int));
alpar@9 560 for (row = npp->r_head, i = 0; row != NULL; row = row->next)
alpar@9 561 npp->row_ref[++i] = row->i;
alpar@9 562 for (col = npp->c_head, j = 0; col != NULL; col = col->next)
alpar@9 563 npp->col_ref[++j] = col->j;
alpar@9 564 /* transformed problem segment is no longer needed */
alpar@9 565 dmp_delete_pool(npp->pool), npp->pool = NULL;
alpar@9 566 npp->name = npp->obj = NULL;
alpar@9 567 npp->c0 = 0.0;
alpar@9 568 npp->r_head = npp->r_tail = NULL;
alpar@9 569 npp->c_head = npp->c_tail = NULL;
alpar@9 570 return;
alpar@9 571 }
alpar@9 572
alpar@9 573 void npp_postprocess(NPP *npp, glp_prob *prob)
alpar@9 574 { /* postprocess solution from the resultant problem */
alpar@9 575 GLPROW *row;
alpar@9 576 GLPCOL *col;
alpar@9 577 NPPTSE *tse;
alpar@9 578 int i, j, k;
alpar@9 579 double dir;
alpar@9 580 xassert(npp->orig_dir == prob->dir);
alpar@9 581 if (npp->orig_dir == GLP_MIN)
alpar@9 582 dir = +1.0;
alpar@9 583 else if (npp->orig_dir == GLP_MAX)
alpar@9 584 dir = -1.0;
alpar@9 585 else
alpar@9 586 xassert(npp != npp);
alpar@9 587 xassert(npp->m == prob->m);
alpar@9 588 xassert(npp->n == prob->n);
alpar@9 589 xassert(npp->nnz == prob->nnz);
alpar@9 590 /* copy solution status */
alpar@9 591 if (npp->sol == GLP_SOL)
alpar@9 592 { npp->p_stat = prob->pbs_stat;
alpar@9 593 npp->d_stat = prob->dbs_stat;
alpar@9 594 }
alpar@9 595 else if (npp->sol == GLP_IPT)
alpar@9 596 npp->t_stat = prob->ipt_stat;
alpar@9 597 else if (npp->sol == GLP_MIP)
alpar@9 598 npp->i_stat = prob->mip_stat;
alpar@9 599 else
alpar@9 600 xassert(npp != npp);
alpar@9 601 /* allocate solution arrays */
alpar@9 602 if (npp->sol == GLP_SOL)
alpar@9 603 { if (npp->r_stat == NULL)
alpar@9 604 npp->r_stat = xcalloc(1+npp->nrows, sizeof(char));
alpar@9 605 for (i = 1; i <= npp->nrows; i++)
alpar@9 606 npp->r_stat[i] = 0;
alpar@9 607 if (npp->c_stat == NULL)
alpar@9 608 npp->c_stat = xcalloc(1+npp->ncols, sizeof(char));
alpar@9 609 for (j = 1; j <= npp->ncols; j++)
alpar@9 610 npp->c_stat[j] = 0;
alpar@9 611 }
alpar@9 612 #if 0
alpar@9 613 if (npp->r_prim == NULL)
alpar@9 614 npp->r_prim = xcalloc(1+npp->nrows, sizeof(double));
alpar@9 615 for (i = 1; i <= npp->nrows; i++)
alpar@9 616 npp->r_prim[i] = DBL_MAX;
alpar@9 617 #endif
alpar@9 618 if (npp->c_value == NULL)
alpar@9 619 npp->c_value = xcalloc(1+npp->ncols, sizeof(double));
alpar@9 620 for (j = 1; j <= npp->ncols; j++)
alpar@9 621 npp->c_value[j] = DBL_MAX;
alpar@9 622 if (npp->sol != GLP_MIP)
alpar@9 623 { if (npp->r_pi == NULL)
alpar@9 624 npp->r_pi = xcalloc(1+npp->nrows, sizeof(double));
alpar@9 625 for (i = 1; i <= npp->nrows; i++)
alpar@9 626 npp->r_pi[i] = DBL_MAX;
alpar@9 627 #if 0
alpar@9 628 if (npp->c_dual == NULL)
alpar@9 629 npp->c_dual = xcalloc(1+npp->ncols, sizeof(double));
alpar@9 630 for (j = 1; j <= npp->ncols; j++)
alpar@9 631 npp->c_dual[j] = DBL_MAX;
alpar@9 632 #endif
alpar@9 633 }
alpar@9 634 /* copy solution components from the resultant problem */
alpar@9 635 if (npp->sol == GLP_SOL)
alpar@9 636 { for (i = 1; i <= npp->m; i++)
alpar@9 637 { row = prob->row[i];
alpar@9 638 k = npp->row_ref[i];
alpar@9 639 npp->r_stat[k] = (char)row->stat;
alpar@9 640 /*npp->r_prim[k] = row->prim;*/
alpar@9 641 npp->r_pi[k] = dir * row->dual;
alpar@9 642 }
alpar@9 643 for (j = 1; j <= npp->n; j++)
alpar@9 644 { col = prob->col[j];
alpar@9 645 k = npp->col_ref[j];
alpar@9 646 npp->c_stat[k] = (char)col->stat;
alpar@9 647 npp->c_value[k] = col->prim;
alpar@9 648 /*npp->c_dual[k] = dir * col->dual;*/
alpar@9 649 }
alpar@9 650 }
alpar@9 651 else if (npp->sol == GLP_IPT)
alpar@9 652 { for (i = 1; i <= npp->m; i++)
alpar@9 653 { row = prob->row[i];
alpar@9 654 k = npp->row_ref[i];
alpar@9 655 /*npp->r_prim[k] = row->pval;*/
alpar@9 656 npp->r_pi[k] = dir * row->dval;
alpar@9 657 }
alpar@9 658 for (j = 1; j <= npp->n; j++)
alpar@9 659 { col = prob->col[j];
alpar@9 660 k = npp->col_ref[j];
alpar@9 661 npp->c_value[k] = col->pval;
alpar@9 662 /*npp->c_dual[k] = dir * col->dval;*/
alpar@9 663 }
alpar@9 664 }
alpar@9 665 else if (npp->sol == GLP_MIP)
alpar@9 666 {
alpar@9 667 #if 0
alpar@9 668 for (i = 1; i <= npp->m; i++)
alpar@9 669 { row = prob->row[i];
alpar@9 670 k = npp->row_ref[i];
alpar@9 671 /*npp->r_prim[k] = row->mipx;*/
alpar@9 672 }
alpar@9 673 #endif
alpar@9 674 for (j = 1; j <= npp->n; j++)
alpar@9 675 { col = prob->col[j];
alpar@9 676 k = npp->col_ref[j];
alpar@9 677 npp->c_value[k] = col->mipx;
alpar@9 678 }
alpar@9 679 }
alpar@9 680 else
alpar@9 681 xassert(npp != npp);
alpar@9 682 /* perform postprocessing to construct solution to the original
alpar@9 683 problem */
alpar@9 684 for (tse = npp->top; tse != NULL; tse = tse->link)
alpar@9 685 { xassert(tse->func != NULL);
alpar@9 686 xassert(tse->func(npp, tse->info) == 0);
alpar@9 687 }
alpar@9 688 return;
alpar@9 689 }
alpar@9 690
alpar@9 691 void npp_unload_sol(NPP *npp, glp_prob *orig)
alpar@9 692 { /* store solution to the original problem */
alpar@9 693 GLPROW *row;
alpar@9 694 GLPCOL *col;
alpar@9 695 int i, j;
alpar@9 696 double dir;
alpar@9 697 xassert(npp->orig_dir == orig->dir);
alpar@9 698 if (npp->orig_dir == GLP_MIN)
alpar@9 699 dir = +1.0;
alpar@9 700 else if (npp->orig_dir == GLP_MAX)
alpar@9 701 dir = -1.0;
alpar@9 702 else
alpar@9 703 xassert(npp != npp);
alpar@9 704 xassert(npp->orig_m == orig->m);
alpar@9 705 xassert(npp->orig_n == orig->n);
alpar@9 706 xassert(npp->orig_nnz == orig->nnz);
alpar@9 707 if (npp->sol == GLP_SOL)
alpar@9 708 { /* store basic solution */
alpar@9 709 orig->valid = 0;
alpar@9 710 orig->pbs_stat = npp->p_stat;
alpar@9 711 orig->dbs_stat = npp->d_stat;
alpar@9 712 orig->obj_val = orig->c0;
alpar@9 713 orig->some = 0;
alpar@9 714 for (i = 1; i <= orig->m; i++)
alpar@9 715 { row = orig->row[i];
alpar@9 716 row->stat = npp->r_stat[i];
alpar@9 717 if (!npp->scaling)
alpar@9 718 { /*row->prim = npp->r_prim[i];*/
alpar@9 719 row->dual = dir * npp->r_pi[i];
alpar@9 720 }
alpar@9 721 else
alpar@9 722 { /*row->prim = npp->r_prim[i] / row->rii;*/
alpar@9 723 row->dual = dir * npp->r_pi[i] * row->rii;
alpar@9 724 }
alpar@9 725 if (row->stat == GLP_BS)
alpar@9 726 row->dual = 0.0;
alpar@9 727 else if (row->stat == GLP_NL)
alpar@9 728 { xassert(row->type == GLP_LO || row->type == GLP_DB);
alpar@9 729 row->prim = row->lb;
alpar@9 730 }
alpar@9 731 else if (row->stat == GLP_NU)
alpar@9 732 { xassert(row->type == GLP_UP || row->type == GLP_DB);
alpar@9 733 row->prim = row->ub;
alpar@9 734 }
alpar@9 735 else if (row->stat == GLP_NF)
alpar@9 736 { xassert(row->type == GLP_FR);
alpar@9 737 row->prim = 0.0;
alpar@9 738 }
alpar@9 739 else if (row->stat == GLP_NS)
alpar@9 740 { xassert(row->type == GLP_FX);
alpar@9 741 row->prim = row->lb;
alpar@9 742 }
alpar@9 743 else
alpar@9 744 xassert(row != row);
alpar@9 745 }
alpar@9 746 for (j = 1; j <= orig->n; j++)
alpar@9 747 { col = orig->col[j];
alpar@9 748 col->stat = npp->c_stat[j];
alpar@9 749 if (!npp->scaling)
alpar@9 750 { col->prim = npp->c_value[j];
alpar@9 751 /*col->dual = dir * npp->c_dual[j];*/
alpar@9 752 }
alpar@9 753 else
alpar@9 754 { col->prim = npp->c_value[j] * col->sjj;
alpar@9 755 /*col->dual = dir * npp->c_dual[j] / col->sjj;*/
alpar@9 756 }
alpar@9 757 if (col->stat == GLP_BS)
alpar@9 758 col->dual = 0.0;
alpar@9 759 #if 1
alpar@9 760 else if (col->stat == GLP_NL)
alpar@9 761 { xassert(col->type == GLP_LO || col->type == GLP_DB);
alpar@9 762 col->prim = col->lb;
alpar@9 763 }
alpar@9 764 else if (col->stat == GLP_NU)
alpar@9 765 { xassert(col->type == GLP_UP || col->type == GLP_DB);
alpar@9 766 col->prim = col->ub;
alpar@9 767 }
alpar@9 768 else if (col->stat == GLP_NF)
alpar@9 769 { xassert(col->type == GLP_FR);
alpar@9 770 col->prim = 0.0;
alpar@9 771 }
alpar@9 772 else if (col->stat == GLP_NS)
alpar@9 773 { xassert(col->type == GLP_FX);
alpar@9 774 col->prim = col->lb;
alpar@9 775 }
alpar@9 776 else
alpar@9 777 xassert(col != col);
alpar@9 778 #endif
alpar@9 779 orig->obj_val += col->coef * col->prim;
alpar@9 780 }
alpar@9 781 #if 1
alpar@9 782 /* compute primal values of inactive rows */
alpar@9 783 for (i = 1; i <= orig->m; i++)
alpar@9 784 { row = orig->row[i];
alpar@9 785 if (row->stat == GLP_BS)
alpar@9 786 { GLPAIJ *aij;
alpar@9 787 double temp;
alpar@9 788 temp = 0.0;
alpar@9 789 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
alpar@9 790 temp += aij->val * aij->col->prim;
alpar@9 791 row->prim = temp;
alpar@9 792 }
alpar@9 793 }
alpar@9 794 /* compute reduced costs of active columns */
alpar@9 795 for (j = 1; j <= orig->n; j++)
alpar@9 796 { col = orig->col[j];
alpar@9 797 if (col->stat != GLP_BS)
alpar@9 798 { GLPAIJ *aij;
alpar@9 799 double temp;
alpar@9 800 temp = col->coef;
alpar@9 801 for (aij = col->ptr; aij != NULL; aij = aij->c_next)
alpar@9 802 temp -= aij->val * aij->row->dual;
alpar@9 803 col->dual = temp;
alpar@9 804 }
alpar@9 805 }
alpar@9 806 #endif
alpar@9 807 }
alpar@9 808 else if (npp->sol == GLP_IPT)
alpar@9 809 { /* store interior-point solution */
alpar@9 810 orig->ipt_stat = npp->t_stat;
alpar@9 811 orig->ipt_obj = orig->c0;
alpar@9 812 for (i = 1; i <= orig->m; i++)
alpar@9 813 { row = orig->row[i];
alpar@9 814 if (!npp->scaling)
alpar@9 815 { /*row->pval = npp->r_prim[i];*/
alpar@9 816 row->dval = dir * npp->r_pi[i];
alpar@9 817 }
alpar@9 818 else
alpar@9 819 { /*row->pval = npp->r_prim[i] / row->rii;*/
alpar@9 820 row->dval = dir * npp->r_pi[i] * row->rii;
alpar@9 821 }
alpar@9 822 }
alpar@9 823 for (j = 1; j <= orig->n; j++)
alpar@9 824 { col = orig->col[j];
alpar@9 825 if (!npp->scaling)
alpar@9 826 { col->pval = npp->c_value[j];
alpar@9 827 /*col->dval = dir * npp->c_dual[j];*/
alpar@9 828 }
alpar@9 829 else
alpar@9 830 { col->pval = npp->c_value[j] * col->sjj;
alpar@9 831 /*col->dval = dir * npp->c_dual[j] / col->sjj;*/
alpar@9 832 }
alpar@9 833 orig->ipt_obj += col->coef * col->pval;
alpar@9 834 }
alpar@9 835 #if 1
alpar@9 836 /* compute row primal values */
alpar@9 837 for (i = 1; i <= orig->m; i++)
alpar@9 838 { row = orig->row[i];
alpar@9 839 { GLPAIJ *aij;
alpar@9 840 double temp;
alpar@9 841 temp = 0.0;
alpar@9 842 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
alpar@9 843 temp += aij->val * aij->col->pval;
alpar@9 844 row->pval = temp;
alpar@9 845 }
alpar@9 846 }
alpar@9 847 /* compute column dual values */
alpar@9 848 for (j = 1; j <= orig->n; j++)
alpar@9 849 { col = orig->col[j];
alpar@9 850 { GLPAIJ *aij;
alpar@9 851 double temp;
alpar@9 852 temp = col->coef;
alpar@9 853 for (aij = col->ptr; aij != NULL; aij = aij->c_next)
alpar@9 854 temp -= aij->val * aij->row->dval;
alpar@9 855 col->dval = temp;
alpar@9 856 }
alpar@9 857 }
alpar@9 858 #endif
alpar@9 859 }
alpar@9 860 else if (npp->sol == GLP_MIP)
alpar@9 861 { /* store MIP solution */
alpar@9 862 xassert(!npp->scaling);
alpar@9 863 orig->mip_stat = npp->i_stat;
alpar@9 864 orig->mip_obj = orig->c0;
alpar@9 865 #if 0
alpar@9 866 for (i = 1; i <= orig->m; i++)
alpar@9 867 { row = orig->row[i];
alpar@9 868 /*row->mipx = npp->r_prim[i];*/
alpar@9 869 }
alpar@9 870 #endif
alpar@9 871 for (j = 1; j <= orig->n; j++)
alpar@9 872 { col = orig->col[j];
alpar@9 873 col->mipx = npp->c_value[j];
alpar@9 874 if (col->kind == GLP_IV)
alpar@9 875 xassert(col->mipx == floor(col->mipx));
alpar@9 876 orig->mip_obj += col->coef * col->mipx;
alpar@9 877 }
alpar@9 878 #if 1
alpar@9 879 /* compute row primal values */
alpar@9 880 for (i = 1; i <= orig->m; i++)
alpar@9 881 { row = orig->row[i];
alpar@9 882 { GLPAIJ *aij;
alpar@9 883 double temp;
alpar@9 884 temp = 0.0;
alpar@9 885 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
alpar@9 886 temp += aij->val * aij->col->mipx;
alpar@9 887 row->mipx = temp;
alpar@9 888 }
alpar@9 889 }
alpar@9 890 #endif
alpar@9 891 }
alpar@9 892 else
alpar@9 893 xassert(npp != npp);
alpar@9 894 return;
alpar@9 895 }
alpar@9 896
alpar@9 897 void npp_delete_wksp(NPP *npp)
alpar@9 898 { /* delete LP/MIP preprocessor workspace */
alpar@9 899 if (npp->pool != NULL)
alpar@9 900 dmp_delete_pool(npp->pool);
alpar@9 901 if (npp->stack != NULL)
alpar@9 902 dmp_delete_pool(npp->stack);
alpar@9 903 if (npp->row_ref != NULL)
alpar@9 904 xfree(npp->row_ref);
alpar@9 905 if (npp->col_ref != NULL)
alpar@9 906 xfree(npp->col_ref);
alpar@9 907 if (npp->r_stat != NULL)
alpar@9 908 xfree(npp->r_stat);
alpar@9 909 #if 0
alpar@9 910 if (npp->r_prim != NULL)
alpar@9 911 xfree(npp->r_prim);
alpar@9 912 #endif
alpar@9 913 if (npp->r_pi != NULL)
alpar@9 914 xfree(npp->r_pi);
alpar@9 915 if (npp->c_stat != NULL)
alpar@9 916 xfree(npp->c_stat);
alpar@9 917 if (npp->c_value != NULL)
alpar@9 918 xfree(npp->c_value);
alpar@9 919 #if 0
alpar@9 920 if (npp->c_dual != NULL)
alpar@9 921 xfree(npp->c_dual);
alpar@9 922 #endif
alpar@9 923 xfree(npp);
alpar@9 924 return;
alpar@9 925 }
alpar@9 926
alpar@9 927 /* eof */