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