praat/external/glpk/glpnpp05.c

/* glpnpp05.c */

/***********************************************************************
*  This code is part of GLPK (GNU Linear Programming Kit).
*
*  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
*  2009, 2010 Andrew Makhorin, Department for Applied Informatics,
*  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
*  E-mail: <mao@gnu.org>.
*
*  GLPK is free software: you can redistribute it and/or modify it
*  under the terms of the GNU General Public License as published by
*  the Free Software Foundation, either version 3 of the License, or
*  (at your option) any later version.
*
*  GLPK is distributed in the hope that it will be useful, but WITHOUT
*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
*  License for more details.
*
*  You should have received a copy of the GNU General Public License
*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
***********************************************************************/

#include "glpnpp.h"

/***********************************************************************
*  NAME
*
*  npp_clean_prob - perform initial LP/MIP processing
*
*  SYNOPSIS
*
*  #include "glpnpp.h"
*  void npp_clean_prob(NPP *npp);
*
*  DESCRIPTION
*
*  The routine npp_clean_prob performs initial LP/MIP processing that
*  currently includes:
*
*  1) removing free rows;
*
*  2) replacing double-sided constraint rows with almost identical
*     bounds, by equality constraint rows;
*
*  3) removing fixed columns;
*
*  4) replacing double-bounded columns with almost identical bounds by
*     fixed columns and removing those columns;
*
*  5) initial processing constraint coefficients (not implemented);
*
*  6) initial processing objective coefficients (not implemented). */

void npp_clean_prob(NPP *npp)
{     /* perform initial LP/MIP processing */
      NPPROW *row, *next_row;
      NPPCOL *col, *next_col;
      int ret;
      xassert(npp == npp);
      /* process rows which originally are free */
      for (row = npp->r_head; row != NULL; row = next_row)
      {  next_row = row->next;
         if (row->lb == -DBL_MAX && row->ub == +DBL_MAX)
         {  /* process free row */
#ifdef GLP_DEBUG
            xprintf("1");
#endif
            npp_free_row(npp, row);
            /* row was deleted */
         }
      }
      /* process rows which originally are double-sided inequalities */
      for (row = npp->r_head; row != NULL; row = next_row)
      {  next_row = row->next;
         if (row->lb != -DBL_MAX && row->ub != +DBL_MAX &&
             row->lb < row->ub)
         {  ret = npp_make_equality(npp, row);
            if (ret == 0)
               ;
            else if (ret == 1)
            {  /* row was replaced by equality constraint */
#ifdef GLP_DEBUG
               xprintf("2");
#endif
            }
            else
               xassert(ret != ret);
         }
      }
      /* process columns which are originally fixed */
      for (col = npp->c_head; col != NULL; col = next_col)
      {  next_col = col->next;
         if (col->lb == col->ub)
         {  /* process fixed column */
#ifdef GLP_DEBUG
            xprintf("3");
#endif
            npp_fixed_col(npp, col);
            /* column was deleted */
         }
      }
      /* process columns which are originally double-bounded */
      for (col = npp->c_head; col != NULL; col = next_col)
      {  next_col = col->next;
         if (col->lb != -DBL_MAX && col->ub != +DBL_MAX &&
             col->lb < col->ub)
         {  ret = npp_make_fixed(npp, col);
            if (ret == 0)
               ;
            else if (ret == 1)
            {  /* column was replaced by fixed column; process it */
#ifdef GLP_DEBUG
               xprintf("4");
#endif
               npp_fixed_col(npp, col);
               /* column was deleted */
            }
         }
      }
      return;
}

/***********************************************************************
*  NAME
*
*  npp_process_row - perform basic row processing
*
*  SYNOPSIS
*
*  #include "glpnpp.h"
*  int npp_process_row(NPP *npp, NPPROW *row, int hard);
*
*  DESCRIPTION
*
*  The routine npp_process_row performs basic row processing that
*  currently includes:
*
*  1) removing empty row;
*
*  2) removing equality constraint row singleton and corresponding
*     column;
*
*  3) removing inequality constraint row singleton and corresponding
*     column if it was fixed;
*
*  4) performing general row analysis;
*
*  5) removing redundant row bounds;
*
*  6) removing forcing row and corresponding columns;
*
*  7) removing row which becomes free due to redundant bounds;
*
*  8) computing implied bounds for all columns in the row and using
*     them to strengthen current column bounds (MIP only, optional,
*     performed if the flag hard is on).
*
*  Additionally the routine may activate affected rows and/or columns
*  for further processing.
*
*  RETURNS
*
*  0           success;
*
*  GLP_ENOPFS  primal/integer infeasibility detected;
*
*  GLP_ENODFS  dual infeasibility detected. */

int npp_process_row(NPP *npp, NPPROW *row, int hard)
{     /* perform basic row processing */
      NPPCOL *col;
      NPPAIJ *aij, *next_aij, *aaa;
      int ret;
      /* row must not be free */
      xassert(!(row->lb == -DBL_MAX && row->ub == +DBL_MAX));
      /* start processing row */
      if (row->ptr == NULL)
      {  /* empty row */
         ret = npp_empty_row(npp, row);
         if (ret == 0)
         {  /* row was deleted */
#ifdef GLP_DEBUG
            xprintf("A");
#endif
            return 0;
         }
         else if (ret == 1)
         {  /* primal infeasibility */
            return GLP_ENOPFS;
         }
         else
            xassert(ret != ret);
      }
      if (row->ptr->r_next == NULL)
      {  /* row singleton */
         col = row->ptr->col;
         if (row->lb == row->ub)
         {  /* equality constraint */
            ret = npp_eq_singlet(npp, row);
            if (ret == 0)
            {  /* column was fixed, row was deleted */
#ifdef GLP_DEBUG
               xprintf("B");
#endif
               /* activate rows affected by column */
               for (aij = col->ptr; aij != NULL; aij = aij->c_next)
                  npp_activate_row(npp, aij->row);
               /* process fixed column */
               npp_fixed_col(npp, col);
               /* column was deleted */
               return 0;
            }
            else if (ret == 1 || ret == 2)
            {  /* primal/integer infeasibility */
               return GLP_ENOPFS;
            }
            else
               xassert(ret != ret);
         }
         else
         {  /* inequality constraint */
            ret = npp_ineq_singlet(npp, row);
            if (0 <= ret && ret <= 3)
            {  /* row was deleted */
#ifdef GLP_DEBUG
               xprintf("C");
#endif
               /* activate column, since its length was changed due to
                  row deletion */
               npp_activate_col(npp, col);
               if (ret >= 2)
               {  /* column bounds changed significantly or column was
                     fixed */
                  /* activate rows affected by column */
                  for (aij = col->ptr; aij != NULL; aij = aij->c_next)
                     npp_activate_row(npp, aij->row);
               }
               if (ret == 3)
               {  /* column was fixed; process it */
#ifdef GLP_DEBUG
                  xprintf("D");
#endif
                  npp_fixed_col(npp, col);
                  /* column was deleted */
               }
               return 0;
            }
            else if (ret == 4)
            {  /* primal infeasibility */
               return GLP_ENOPFS;
            }
            else
               xassert(ret != ret);
         }
      }
#if 0
      /* sometimes this causes too large round-off errors; probably
         pivot coefficient should be chosen more carefully */
      if (row->ptr->r_next->r_next == NULL)
      {  /* row doubleton */
         if (row->lb == row->ub)
         {  /* equality constraint */
            if (!(row->ptr->col->is_int ||
                  row->ptr->r_next->col->is_int))
            {  /* both columns are continuous */
               NPPCOL *q;
               q = npp_eq_doublet(npp, row);
               if (q != NULL)
               {  /* column q was eliminated */
#ifdef GLP_DEBUG
                  xprintf("E");
#endif
                  /* now column q is singleton of type "implied slack
                     variable"; we process it here to make sure that on
                     recovering basic solution the row is always active
                     equality constraint (as required by the routine
                     rcv_eq_doublet) */
                  xassert(npp_process_col(npp, q) == 0);
                  /* column q was deleted; note that row p also may be
                     deleted */
                  return 0;
               }
            }
         }
      }
#endif
      /* general row analysis */
      ret = npp_analyze_row(npp, row);
      xassert(0x00 <= ret && ret <= 0xFF);
      if (ret == 0x33)
      {  /* row bounds are inconsistent with column bounds */
         return GLP_ENOPFS;
      }
      if ((ret & 0x0F) == 0x00)
      {  /* row lower bound does not exist or redundant */
         if (row->lb != -DBL_MAX)
         {  /* remove redundant row lower bound */
#ifdef GLP_DEBUG
            xprintf("F");
#endif
            npp_inactive_bound(npp, row, 0);
         }
      }
      else if ((ret & 0x0F) == 0x01)
      {  /* row lower bound can be active */
         /* see below */
      }
      else if ((ret & 0x0F) == 0x02)
      {  /* row lower bound is a forcing bound */
#ifdef GLP_DEBUG
         xprintf("G");
#endif
         /* process forcing row */
         if (npp_forcing_row(npp, row, 0) == 0)
fixup:   {  /* columns were fixed, row was made free */
            for (aij = row->ptr; aij != NULL; aij = next_aij)
            {  /* process column fixed by forcing row */
#ifdef GLP_DEBUG
               xprintf("H");
#endif
               col = aij->col;
               next_aij = aij->r_next;
               /* activate rows affected by column */
               for (aaa = col->ptr; aaa != NULL; aaa = aaa->c_next)
                  npp_activate_row(npp, aaa->row);
               /* process fixed column */
               npp_fixed_col(npp, col);
               /* column was deleted */
            }
            /* process free row (which now is empty due to deletion of
               all its columns) */
            npp_free_row(npp, row);
            /* row was deleted */
            return 0;
         }
      }
      else
         xassert(ret != ret);
      if ((ret & 0xF0) == 0x00)
      {  /* row upper bound does not exist or redundant */
         if (row->ub != +DBL_MAX)
         {  /* remove redundant row upper bound */
#ifdef GLP_DEBUG
            xprintf("I");
#endif
            npp_inactive_bound(npp, row, 1);
         }
      }
      else if ((ret & 0xF0) == 0x10)
      {  /* row upper bound can be active */
         /* see below */
      }
      else if ((ret & 0xF0) == 0x20)
      {  /* row upper bound is a forcing bound */
#ifdef GLP_DEBUG
         xprintf("J");
#endif
         /* process forcing row */
         if (npp_forcing_row(npp, row, 1) == 0) goto fixup;
      }
      else
         xassert(ret != ret);
      if (row->lb == -DBL_MAX && row->ub == +DBL_MAX)
      {  /* row became free due to redundant bounds removal */
#ifdef GLP_DEBUG
         xprintf("K");
#endif
         /* activate its columns, since their length will change due
            to row deletion */
         for (aij = row->ptr; aij != NULL; aij = aij->r_next)
            npp_activate_col(npp, aij->col);
         /* process free row */
         npp_free_row(npp, row);
         /* row was deleted */
         return 0;
      }
#if 1 /* 23/XII-2009 */
      /* row lower and/or upper bounds can be active */
      if (npp->sol == GLP_MIP && hard)
      {  /* improve current column bounds (optional) */
         if (npp_improve_bounds(npp, row, 1) < 0)
            return GLP_ENOPFS;
      }
#endif
      return 0;
}

/***********************************************************************
*  NAME
*
*  npp_improve_bounds - improve current column bounds
*
*  SYNOPSIS
*
*  #include "glpnpp.h"
*  int npp_improve_bounds(NPP *npp, NPPROW *row, int flag);
*
*  DESCRIPTION
*
*  The routine npp_improve_bounds analyzes specified row (inequality
*  or equality constraint) to determine implied column bounds and then
*  uses these bounds to improve (strengthen) current column bounds.
*
*  If the flag is on and current column bounds changed significantly
*  or the column was fixed, the routine activate rows affected by the
*  column for further processing. (This feature is intended to be used
*  in the main loop of the routine npp_process_row.)
*
*  NOTE: This operation can be used for MIP problem only.
*
*  RETURNS
*
*  The routine npp_improve_bounds returns the number of significantly
*  changed bounds plus the number of column having been fixed due to
*  bound improvements. However, if the routine detects primal/integer
*  infeasibility, it returns a negative value. */

int npp_improve_bounds(NPP *npp, NPPROW *row, int flag)
{     /* improve current column bounds */
      NPPCOL *col;
      NPPAIJ *aij, *next_aij, *aaa;
      int kase, ret, count = 0;
      double lb, ub;
      xassert(npp->sol == GLP_MIP);
      /* row must not be free */
      xassert(!(row->lb == -DBL_MAX && row->ub == +DBL_MAX));
      /* determine implied column bounds */
      npp_implied_bounds(npp, row);
      /* and use these bounds to strengthen current column bounds */
      for (aij = row->ptr; aij != NULL; aij = next_aij)
      {  col = aij->col;
         next_aij = aij->r_next;
         for (kase = 0; kase <= 1; kase++)
         {  /* save current column bounds */
            lb = col->lb, ub = col->ub;
            if (kase == 0)
            {  /* process implied column lower bound */
               if (col->ll.ll == -DBL_MAX) continue;
               ret = npp_implied_lower(npp, col, col->ll.ll);
            }
            else
            {  /* process implied column upper bound */
               if (col->uu.uu == +DBL_MAX) continue;
               ret = npp_implied_upper(npp, col, col->uu.uu);
            }
            if (ret == 0 || ret == 1)
            {  /* current column bounds did not change or changed, but
                  not significantly; restore current column bounds */
               col->lb = lb, col->ub = ub;
            }
            else if (ret == 2 || ret == 3)
            {  /* current column bounds changed significantly or column
                  was fixed */
#ifdef GLP_DEBUG
               xprintf("L");
#endif
               count++;
               /* activate other rows affected by column, if required */
               if (flag)
               {  for (aaa = col->ptr; aaa != NULL; aaa = aaa->c_next)
                  {  if (aaa->row != row)
                        npp_activate_row(npp, aaa->row);
                  }
               }
               if (ret == 3)
               {  /* process fixed column */
#ifdef GLP_DEBUG
                  xprintf("M");
#endif
                  npp_fixed_col(npp, col);
                  /* column was deleted */
                  break; /* for kase */
               }
            }
            else if (ret == 4)
            {  /* primal/integer infeasibility */
               return -1;
            }
            else
               xassert(ret != ret);
         }
      }
      return count;
}

/***********************************************************************
*  NAME
*
*  npp_process_col - perform basic column processing
*
*  SYNOPSIS
*
*  #include "glpnpp.h"
*  int npp_process_col(NPP *npp, NPPCOL *col);
*
*  DESCRIPTION
*
*  The routine npp_process_col performs basic column processing that
*  currently includes:
*
*  1) fixing and removing empty column;
*
*  2) removing column singleton, which is implied slack variable, and
*     corresponding row if it becomes free;
*
*  3) removing bounds of column, which is implied free variable, and
*     replacing corresponding row by equality constraint.
*
*  Additionally the routine may activate affected rows and/or columns
*  for further processing.
*
*  RETURNS
*
*  0           success;
*
*  GLP_ENOPFS  primal/integer infeasibility detected;
*
*  GLP_ENODFS  dual infeasibility detected. */

int npp_process_col(NPP *npp, NPPCOL *col)
{     /* perform basic column processing */
      NPPROW *row;
      NPPAIJ *aij;
      int ret;
      /* column must not be fixed */
      xassert(col->lb < col->ub);
      /* start processing column */
      if (col->ptr == NULL)
      {  /* empty column */
         ret = npp_empty_col(npp, col);
         if (ret == 0)
         {  /* column was fixed and deleted */
#ifdef GLP_DEBUG
            xprintf("N");
#endif
            return 0;
         }
         else if (ret == 1)
         {  /* dual infeasibility */
            return GLP_ENODFS;
         }
         else
            xassert(ret != ret);
      }
      if (col->ptr->c_next == NULL)
      {  /* column singleton */
         row = col->ptr->row;
         if (row->lb == row->ub)
         {  /* equality constraint */
            if (!col->is_int)
slack:      {  /* implied slack variable */
#ifdef GLP_DEBUG
               xprintf("O");
#endif
               npp_implied_slack(npp, col);
               /* column was deleted */
               if (row->lb == -DBL_MAX && row->ub == +DBL_MAX)
               {  /* row became free due to implied slack variable */
#ifdef GLP_DEBUG
                  xprintf("P");
#endif
                  /* activate columns affected by row */
                  for (aij = row->ptr; aij != NULL; aij = aij->r_next)
                     npp_activate_col(npp, aij->col);
                  /* process free row */
                  npp_free_row(npp, row);
                  /* row was deleted */
               }
               else
               {  /* row became inequality constraint; activate it
                     since its length changed due to column deletion */
                  npp_activate_row(npp, row);
               }
               return 0;
            }
         }
         else
         {  /* inequality constraint */
            if (!col->is_int)
            {  ret = npp_implied_free(npp, col);
               if (ret == 0)
               {  /* implied free variable */
#ifdef GLP_DEBUG
                  xprintf("Q");
#endif
                  /* column bounds were removed, row was replaced by
                     equality constraint */
                  goto slack;
               }
               else if (ret == 1)
               {  /* column is not implied free variable, because its
                     lower and/or upper bounds can be active */
               }
               else if (ret == 2)
               {  /* dual infeasibility */
                  return GLP_ENODFS;
               }
            }
         }
      }
      /* column still exists */
      return 0;
}

/***********************************************************************
*  NAME
*
*  npp_process_prob - perform basic LP/MIP processing
*
*  SYNOPSIS
*
*  #include "glpnpp.h"
*  int npp_process_prob(NPP *npp, int hard);
*
*  DESCRIPTION
*
*  The routine npp_process_prob performs basic LP/MIP processing that
*  currently includes:
*
*  1) initial LP/MIP processing (see the routine npp_clean_prob),
*
*  2) basic row processing (see the routine npp_process_row), and
*
*  3) basic column processing (see the routine npp_process_col).
*
*  If the flag hard is on, the routine attempts to improve current
*  column bounds multiple times within the main processing loop, in
*  which case this feature may take a time. Otherwise, if the flag hard
*  is off, improving column bounds is performed only once at the end of
*  the main loop. (Note that this feature is used for MIP only.)
*
*  The routine uses two sets: the set of active rows and the set of
*  active columns. Rows/columns are marked by a flag (the field temp in
*  NPPROW/NPPCOL). If the flag is non-zero, the row/column is active,
*  in which case it is placed in the beginning of the row/column list;
*  otherwise, if the flag is zero, the row/column is inactive, in which
*  case it is placed in the end of the row/column list. If a row/column
*  being currently processed may affect other rows/columns, the latters
*  are activated for further processing.
*
*  RETURNS
*
*  0           success;
*
*  GLP_ENOPFS  primal/integer infeasibility detected;
*
*  GLP_ENODFS  dual infeasibility detected. */

int npp_process_prob(NPP *npp, int hard)
{     /* perform basic LP/MIP processing */
      NPPROW *row;
      NPPCOL *col;
      int processing, ret;
      /* perform initial LP/MIP processing */
      npp_clean_prob(npp);
      /* activate all remaining rows and columns */
      for (row = npp->r_head; row != NULL; row = row->next)
         row->temp = 1;
      for (col = npp->c_head; col != NULL; col = col->next)
         col->temp = 1;
      /* main processing loop */
      processing = 1;
      while (processing)
      {  processing = 0;
         /* process all active rows */
         for (;;)
         {  row = npp->r_head;
            if (row == NULL || !row->temp) break;
            npp_deactivate_row(npp, row);
            ret = npp_process_row(npp, row, hard);
            if (ret != 0) goto done;
            processing = 1;
         }
         /* process all active columns */
         for (;;)
         {  col = npp->c_head;
            if (col == NULL || !col->temp) break;
            npp_deactivate_col(npp, col);
            ret = npp_process_col(npp, col);
            if (ret != 0) goto done;
            processing = 1;
         }
      }
#if 1 /* 23/XII-2009 */
      if (npp->sol == GLP_MIP && !hard)
      {  /* improve current column bounds (optional) */
         for (row = npp->r_head; row != NULL; row = row->next)
         {  if (npp_improve_bounds(npp, row, 0) < 0)
            {  ret = GLP_ENOPFS;
               goto done;
            }
         }
      }
#endif
      /* all seems ok */
      ret = 0;
done: xassert(ret == 0 || ret == GLP_ENOPFS || ret == GLP_ENODFS);
#ifdef GLP_DEBUG
      xprintf("\n");
#endif
      return ret;
}

/**********************************************************************/

int npp_simplex(NPP *npp, const glp_smcp *parm)
{     /* process LP prior to applying primal/dual simplex method */
      int ret;
      xassert(npp->sol == GLP_SOL);
      xassert(parm == parm);
      ret = npp_process_prob(npp, 0);
      return ret;
}

/**********************************************************************/

int npp_integer(NPP *npp, const glp_iocp *parm)
{     /* process MIP prior to applying branch-and-bound method */
      NPPROW *row, *prev_row;
      NPPCOL *col;
      NPPAIJ *aij;
      int count, ret;
      xassert(npp->sol == GLP_MIP);
      xassert(parm == parm);
      /*==============================================================*/
      /* perform basic MIP processing */
      ret = npp_process_prob(npp, 1);
      if (ret != 0) goto done;
      /*==============================================================*/
      /* binarize problem, if required */
      if (parm->binarize)
         npp_binarize_prob(npp);
      /*==============================================================*/
      /* identify hidden packing inequalities */
      count = 0;
      /* new rows will be added to the end of the row list, so we go
         from the end to beginning of the row list */
      for (row = npp->r_tail; row != NULL; row = prev_row)
      {  prev_row = row->prev;
         /* skip free row */
         if (row->lb == -DBL_MAX && row->ub == +DBL_MAX) continue;
         /* skip equality constraint */
         if (row->lb == row->ub) continue;
         /* skip row having less than two variables */
         if (row->ptr == NULL || row->ptr->r_next == NULL) continue;
         /* skip row having non-binary variables */
         for (aij = row->ptr; aij != NULL; aij = aij->r_next)
         {  col = aij->col;
            if (!(col->is_int && col->lb == 0.0 && col->ub == 1.0))
               break;
         }
         if (aij != NULL) continue;
         count += npp_hidden_packing(npp, row);
      }
      if (count > 0)
         xprintf("%d hidden packing inequaliti(es) were detected\n",
            count);
      /*==============================================================*/
      /* identify hidden covering inequalities */
      count = 0;
      /* new rows will be added to the end of the row list, so we go
         from the end to beginning of the row list */
      for (row = npp->r_tail; row != NULL; row = prev_row)
      {  prev_row = row->prev;
         /* skip free row */
         if (row->lb == -DBL_MAX && row->ub == +DBL_MAX) continue;
         /* skip equality constraint */
         if (row->lb == row->ub) continue;
         /* skip row having less than three variables */
         if (row->ptr == NULL || row->ptr->r_next == NULL ||
             row->ptr->r_next->r_next == NULL) continue;
         /* skip row having non-binary variables */
         for (aij = row->ptr; aij != NULL; aij = aij->r_next)
         {  col = aij->col;
            if (!(col->is_int && col->lb == 0.0 && col->ub == 1.0))
               break;
         }
         if (aij != NULL) continue;
         count += npp_hidden_covering(npp, row);
      }
      if (count > 0)
         xprintf("%d hidden covering inequaliti(es) were detected\n",
            count);
      /*==============================================================*/
      /* reduce inequality constraint coefficients */
      count = 0;
      /* new rows will be added to the end of the row list, so we go
         from the end to beginning of the row list */
      for (row = npp->r_tail; row != NULL; row = prev_row)
      {  prev_row = row->prev;
         /* skip equality constraint */
         if (row->lb == row->ub) continue;
         count += npp_reduce_ineq_coef(npp, row);
      }
      if (count > 0)
         xprintf("%d constraint coefficient(s) were reduced\n", count);
      /*==============================================================*/
#ifdef GLP_DEBUG
      routine(npp);
#endif
      /*==============================================================*/
      /* all seems ok */
      ret = 0;
done: return ret;
}

/* eof */