/* glpini02.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 "glpapi.h"

struct var

{     /* structural variable */

      int j;

      /* ordinal number */

      double q;

      /* penalty value */

};

static int fcmp(const void *ptr1, const void *ptr2)

{     /* this routine is passed to the qsort() function */

      struct var *col1 = (void *)ptr1, *col2 = (void *)ptr2;

      if (col1->q < col2->q) return -1;

      if (col1->q > col2->q) return +1;

      return 0;

}

static int get_column(glp_prob *lp, int j, int ind[], double val[])

{     /* Bixby's algorithm assumes that the constraint matrix is scaled

         such that the maximum absolute value in every non-zero row and

         column is 1 */

      int k, len;

      double big;

      len = glp_get_mat_col(lp, j, ind, val);

      big = 0.0;

      for (k = 1; k <= len; k++)

         if (big < fabs(val[k])) big = fabs(val[k]);

      if (big == 0.0) big = 1.0;

      for (k = 1; k <= len; k++) val[k] /= big;

      return len;

}

static void cpx_basis(glp_prob *lp)

{     /* main routine */

      struct var *C, *C2, *C3, *C4;

      int m, n, i, j, jk, k, l, ll, t, n2, n3, n4, type, len, *I, *r,

         *ind;

      double alpha, gamma, cmax, temp, *v, *val;

      xprintf("Constructing initial basis...\n");

      /* determine the number of rows and columns */

      m = glp_get_num_rows(lp);

      n = glp_get_num_cols(lp);

      /* allocate working arrays */

      C = xcalloc(1+n, sizeof(struct var));

      I = xcalloc(1+m, sizeof(int));

      r = xcalloc(1+m, sizeof(int));

      v = xcalloc(1+m, sizeof(double));

      ind = xcalloc(1+m, sizeof(int));

      val = xcalloc(1+m, sizeof(double));

      /* make all auxiliary variables non-basic */

      for (i = 1; i <= m; i++)

      {  if (glp_get_row_type(lp, i) != GLP_DB)

            glp_set_row_stat(lp, i, GLP_NS);

         else if (fabs(glp_get_row_lb(lp, i)) <=

                  fabs(glp_get_row_ub(lp, i)))

            glp_set_row_stat(lp, i, GLP_NL);

         else

            glp_set_row_stat(lp, i, GLP_NU);

      }

      /* make all structural variables non-basic */

      for (j = 1; j <= n; j++)

      {  if (glp_get_col_type(lp, j) != GLP_DB)

            glp_set_col_stat(lp, j, GLP_NS);

         else if (fabs(glp_get_col_lb(lp, j)) <=

                  fabs(glp_get_col_ub(lp, j)))

            glp_set_col_stat(lp, j, GLP_NL);

         else

            glp_set_col_stat(lp, j, GLP_NU);

      }

      /* C2 is a set of free structural variables */

      n2 = 0, C2 = C + 0;

      for (j = 1; j <= n; j++)

      {  type = glp_get_col_type(lp, j);

         if (type == GLP_FR)

         {  n2++;

            C2[n2].j = j;

            C2[n2].q = 0.0;

         }

      }

      /* C3 is a set of structural variables having excatly one (lower

         or upper) bound */

      n3 = 0, C3 = C2 + n2;

      for (j = 1; j <= n; j++)

      {  type = glp_get_col_type(lp, j);

         if (type == GLP_LO)

         {  n3++;

            C3[n3].j = j;

            C3[n3].q = + glp_get_col_lb(lp, j);

         }

         else if (type == GLP_UP)

         {  n3++;

            C3[n3].j = j;

            C3[n3].q = - glp_get_col_ub(lp, j);

         }

      }

      /* C4 is a set of structural variables having both (lower and

         upper) bounds */

      n4 = 0, C4 = C3 + n3;

      for (j = 1; j <= n; j++)

      {  type = glp_get_col_type(lp, j);

         if (type == GLP_DB)

         {  n4++;

            C4[n4].j = j;

            C4[n4].q = glp_get_col_lb(lp, j) - glp_get_col_ub(lp, j);

         }

      }

      /* compute gamma = max{|c[j]|: 1 <= j <= n} */

      gamma = 0.0;

      for (j = 1; j <= n; j++)

      {  temp = fabs(glp_get_obj_coef(lp, j));

         if (gamma < temp) gamma = temp;

      }

      /* compute cmax */

      cmax = (gamma == 0.0 ? 1.0 : 1000.0 * gamma);

      /* compute final penalty for all structural variables within sets

         C2, C3, and C4 */

      switch (glp_get_obj_dir(lp))

      {  case GLP_MIN: temp = +1.0; break;

         case GLP_MAX: temp = -1.0; break;

         default: xassert(lp != lp);

      }

      for (k = 1; k <= n2+n3+n4; k++)

      {  j = C[k].j;

         C[k].q += (temp * glp_get_obj_coef(lp, j)) / cmax;

      }

      /* sort structural variables within C2, C3, and C4 in ascending

         order of penalty value */

      qsort(C2+1, n2, sizeof(struct var), fcmp);

      for (k = 1; k < n2; k++) xassert(C2[k].q <= C2[k+1].q);

      qsort(C3+1, n3, sizeof(struct var), fcmp);

      for (k = 1; k < n3; k++) xassert(C3[k].q <= C3[k+1].q);

      qsort(C4+1, n4, sizeof(struct var), fcmp);

      for (k = 1; k < n4; k++) xassert(C4[k].q <= C4[k+1].q);

      /*** STEP 1 ***/

      for (i = 1; i <= m; i++)

      {  type = glp_get_row_type(lp, i);

         if (type != GLP_FX)

         {  /* row i is either free or inequality constraint */

            glp_set_row_stat(lp, i, GLP_BS);

            I[i] = 1;

            r[i] = 1;

         }

         else

         {  /* row i is equality constraint */

            I[i] = 0;

            r[i] = 0;

         }

         v[i] = +DBL_MAX;

      }

      /*** STEP 2 ***/

      for (k = 1; k <= n2+n3+n4; k++)

      {  jk = C[k].j;

         len = get_column(lp, jk, ind, val);

         /* let alpha = max{|A[l,jk]|: r[l] = 0} and let l' be such

            that alpha = |A[l',jk]| */

         alpha = 0.0, ll = 0;

         for (t = 1; t <= len; t++)

         {  l = ind[t];

            if (r[l] == 0 && alpha < fabs(val[t]))

               alpha = fabs(val[t]), ll = l;

         }

         if (alpha >= 0.99)

         {  /* B := B union {jk} */

            glp_set_col_stat(lp, jk, GLP_BS);

            I[ll] = 1;

            v[ll] = alpha;

            /* r[l] := r[l] + 1 for all l such that |A[l,jk]| != 0 */

            for (t = 1; t <= len; t++)

            {  l = ind[t];

               if (val[t] != 0.0) r[l]++;

            }

            /* continue to the next k */

            continue;

         }

         /* if |A[l,jk]| > 0.01 * v[l] for some l, continue to the

            next k */

         for (t = 1; t <= len; t++)

         {  l = ind[t];

            if (fabs(val[t]) > 0.01 * v[l]) break;

         }

         if (t <= len) continue;

         /* otherwise, let alpha = max{|A[l,jk]|: I[l] = 0} and let l'

            be such that alpha = |A[l',jk]| */

         alpha = 0.0, ll = 0;

         for (t = 1; t <= len; t++)

         {  l = ind[t];

            if (I[l] == 0 && alpha < fabs(val[t]))

               alpha = fabs(val[t]), ll = l;

         }

         /* if alpha = 0, continue to the next k */

         if (alpha == 0.0) continue;

         /* B := B union {jk} */

         glp_set_col_stat(lp, jk, GLP_BS);

         I[ll] = 1;

         v[ll] = alpha;

         /* r[l] := r[l] + 1 for all l such that |A[l,jk]| != 0 */

         for (t = 1; t <= len; t++)

         {  l = ind[t];

            if (val[t] != 0.0) r[l]++;

         }

      }

      /*** STEP 3 ***/

      /* add an artificial variable (auxiliary variable for equality

         constraint) to cover each remaining uncovered row */

      for (i = 1; i <= m; i++)

         if (I[i] == 0) glp_set_row_stat(lp, i, GLP_BS);

      /* free working arrays */

      xfree(C);

      xfree(I);

      xfree(r);

      xfree(v);

      xfree(ind);

      xfree(val);

      return;

}

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

*  NAME

*

*  glp_cpx_basis - construct Bixby's initial LP basis

*

*  SYNOPSIS

*

*  void glp_cpx_basis(glp_prob *lp);

*

*  DESCRIPTION

*

*  The routine glp_cpx_basis constructs an advanced initial basis for

*  the specified problem object.

*

*  The routine is based on Bixby's algorithm described in the paper:

*

*  Robert E. Bixby. Implementing the Simplex Method: The Initial Basis.

*  ORSA Journal on Computing, Vol. 4, No. 3, 1992, pp. 267-84. */

void glp_cpx_basis(glp_prob *lp)

{     if (lp->m == 0 || lp->n == 0)

         glp_std_basis(lp);

      else

         cpx_basis(lp);

      return;

}

/* eof */