src/glpapi12.c
author Alpar Juttner <alpar@cs.elte.hu>
Mon, 06 Dec 2010 13:09:21 +0100
changeset 1 c445c931472f
permissions -rw-r--r--
Import glpk-4.45

- Generated files and doc/notes are removed
alpar@1
     1
/* glpapi12.c (basis factorization and simplex tableau routines) */
alpar@1
     2
alpar@1
     3
/***********************************************************************
alpar@1
     4
*  This code is part of GLPK (GNU Linear Programming Kit).
alpar@1
     5
*
alpar@1
     6
*  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
alpar@1
     7
*  2009, 2010 Andrew Makhorin, Department for Applied Informatics,
alpar@1
     8
*  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
alpar@1
     9
*  E-mail: <mao@gnu.org>.
alpar@1
    10
*
alpar@1
    11
*  GLPK is free software: you can redistribute it and/or modify it
alpar@1
    12
*  under the terms of the GNU General Public License as published by
alpar@1
    13
*  the Free Software Foundation, either version 3 of the License, or
alpar@1
    14
*  (at your option) any later version.
alpar@1
    15
*
alpar@1
    16
*  GLPK is distributed in the hope that it will be useful, but WITHOUT
alpar@1
    17
*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
alpar@1
    18
*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
alpar@1
    19
*  License for more details.
alpar@1
    20
*
alpar@1
    21
*  You should have received a copy of the GNU General Public License
alpar@1
    22
*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
alpar@1
    23
***********************************************************************/
alpar@1
    24
alpar@1
    25
#include "glpapi.h"
alpar@1
    26
alpar@1
    27
/***********************************************************************
alpar@1
    28
*  NAME
alpar@1
    29
*
alpar@1
    30
*  glp_bf_exists - check if the basis factorization exists
alpar@1
    31
*
alpar@1
    32
*  SYNOPSIS
alpar@1
    33
*
alpar@1
    34
*  int glp_bf_exists(glp_prob *lp);
alpar@1
    35
*
alpar@1
    36
*  RETURNS
alpar@1
    37
*
alpar@1
    38
*  If the basis factorization for the current basis associated with
alpar@1
    39
*  the specified problem object exists and therefore is available for
alpar@1
    40
*  computations, the routine glp_bf_exists returns non-zero. Otherwise
alpar@1
    41
*  the routine returns zero. */
alpar@1
    42
alpar@1
    43
int glp_bf_exists(glp_prob *lp)
alpar@1
    44
{     int ret;
alpar@1
    45
      ret = (lp->m == 0 || lp->valid);
alpar@1
    46
      return ret;
alpar@1
    47
}
alpar@1
    48
alpar@1
    49
/***********************************************************************
alpar@1
    50
*  NAME
alpar@1
    51
*
alpar@1
    52
*  glp_factorize - compute the basis factorization
alpar@1
    53
*
alpar@1
    54
*  SYNOPSIS
alpar@1
    55
*
alpar@1
    56
*  int glp_factorize(glp_prob *lp);
alpar@1
    57
*
alpar@1
    58
*  DESCRIPTION
alpar@1
    59
*
alpar@1
    60
*  The routine glp_factorize computes the basis factorization for the
alpar@1
    61
*  current basis associated with the specified problem object.
alpar@1
    62
*
alpar@1
    63
*  RETURNS
alpar@1
    64
*
alpar@1
    65
*  0  The basis factorization has been successfully computed.
alpar@1
    66
*
alpar@1
    67
*  GLP_EBADB
alpar@1
    68
*     The basis matrix is invalid, i.e. the number of basic (auxiliary
alpar@1
    69
*     and structural) variables differs from the number of rows in the
alpar@1
    70
*     problem object.
alpar@1
    71
*
alpar@1
    72
*  GLP_ESING
alpar@1
    73
*     The basis matrix is singular within the working precision.
alpar@1
    74
*
alpar@1
    75
*  GLP_ECOND
alpar@1
    76
*     The basis matrix is ill-conditioned. */
alpar@1
    77
alpar@1
    78
static int b_col(void *info, int j, int ind[], double val[])
alpar@1
    79
{     glp_prob *lp = info;
alpar@1
    80
      int m = lp->m;
alpar@1
    81
      GLPAIJ *aij;
alpar@1
    82
      int k, len;
alpar@1
    83
      xassert(1 <= j && j <= m);
alpar@1
    84
      /* determine the ordinal number of basic auxiliary or structural
alpar@1
    85
         variable x[k] corresponding to basic variable xB[j] */
alpar@1
    86
      k = lp->head[j];
alpar@1
    87
      /* build j-th column of the basic matrix, which is k-th column of
alpar@1
    88
         the scaled augmented matrix (I | -R*A*S) */
alpar@1
    89
      if (k <= m)
alpar@1
    90
      {  /* x[k] is auxiliary variable */
alpar@1
    91
         len = 1;
alpar@1
    92
         ind[1] = k;
alpar@1
    93
         val[1] = 1.0;
alpar@1
    94
      }
alpar@1
    95
      else
alpar@1
    96
      {  /* x[k] is structural variable */
alpar@1
    97
         len = 0;
alpar@1
    98
         for (aij = lp->col[k-m]->ptr; aij != NULL; aij = aij->c_next)
alpar@1
    99
         {  len++;
alpar@1
   100
            ind[len] = aij->row->i;
alpar@1
   101
            val[len] = - aij->row->rii * aij->val * aij->col->sjj;
alpar@1
   102
         }
alpar@1
   103
      }
alpar@1
   104
      return len;
alpar@1
   105
}
alpar@1
   106
alpar@1
   107
static void copy_bfcp(glp_prob *lp);
alpar@1
   108
alpar@1
   109
int glp_factorize(glp_prob *lp)
alpar@1
   110
{     int m = lp->m;
alpar@1
   111
      int n = lp->n;
alpar@1
   112
      GLPROW **row = lp->row;
alpar@1
   113
      GLPCOL **col = lp->col;
alpar@1
   114
      int *head = lp->head;
alpar@1
   115
      int j, k, stat, ret;
alpar@1
   116
      /* invalidate the basis factorization */
alpar@1
   117
      lp->valid = 0;
alpar@1
   118
      /* build the basis header */
alpar@1
   119
      j = 0;
alpar@1
   120
      for (k = 1; k <= m+n; k++)
alpar@1
   121
      {  if (k <= m)
alpar@1
   122
         {  stat = row[k]->stat;
alpar@1
   123
            row[k]->bind = 0;
alpar@1
   124
         }
alpar@1
   125
         else
alpar@1
   126
         {  stat = col[k-m]->stat;
alpar@1
   127
            col[k-m]->bind = 0;
alpar@1
   128
         }
alpar@1
   129
         if (stat == GLP_BS)
alpar@1
   130
         {  j++;
alpar@1
   131
            if (j > m)
alpar@1
   132
            {  /* too many basic variables */
alpar@1
   133
               ret = GLP_EBADB;
alpar@1
   134
               goto fini;
alpar@1
   135
            }
alpar@1
   136
            head[j] = k;
alpar@1
   137
            if (k <= m)
alpar@1
   138
               row[k]->bind = j;
alpar@1
   139
            else
alpar@1
   140
               col[k-m]->bind = j;
alpar@1
   141
         }
alpar@1
   142
      }
alpar@1
   143
      if (j < m)
alpar@1
   144
      {  /* too few basic variables */
alpar@1
   145
         ret = GLP_EBADB;
alpar@1
   146
         goto fini;
alpar@1
   147
      }
alpar@1
   148
      /* try to factorize the basis matrix */
alpar@1
   149
      if (m > 0)
alpar@1
   150
      {  if (lp->bfd == NULL)
alpar@1
   151
         {  lp->bfd = bfd_create_it();
alpar@1
   152
            copy_bfcp(lp);
alpar@1
   153
         }
alpar@1
   154
         switch (bfd_factorize(lp->bfd, m, lp->head, b_col, lp))
alpar@1
   155
         {  case 0:
alpar@1
   156
               /* ok */
alpar@1
   157
               break;
alpar@1
   158
            case BFD_ESING:
alpar@1
   159
               /* singular matrix */
alpar@1
   160
               ret = GLP_ESING;
alpar@1
   161
               goto fini;
alpar@1
   162
            case BFD_ECOND:
alpar@1
   163
               /* ill-conditioned matrix */
alpar@1
   164
               ret = GLP_ECOND;
alpar@1
   165
               goto fini;
alpar@1
   166
            default:
alpar@1
   167
               xassert(lp != lp);
alpar@1
   168
         }
alpar@1
   169
         lp->valid = 1;
alpar@1
   170
      }
alpar@1
   171
      /* factorization successful */
alpar@1
   172
      ret = 0;
alpar@1
   173
fini: /* bring the return code to the calling program */
alpar@1
   174
      return ret;
alpar@1
   175
}
alpar@1
   176
alpar@1
   177
/***********************************************************************
alpar@1
   178
*  NAME
alpar@1
   179
*
alpar@1
   180
*  glp_bf_updated - check if the basis factorization has been updated
alpar@1
   181
*
alpar@1
   182
*  SYNOPSIS
alpar@1
   183
*
alpar@1
   184
*  int glp_bf_updated(glp_prob *lp);
alpar@1
   185
*
alpar@1
   186
*  RETURNS
alpar@1
   187
*
alpar@1
   188
*  If the basis factorization has been just computed from scratch, the
alpar@1
   189
*  routine glp_bf_updated returns zero. Otherwise, if the factorization
alpar@1
   190
*  has been updated one or more times, the routine returns non-zero. */
alpar@1
   191
alpar@1
   192
int glp_bf_updated(glp_prob *lp)
alpar@1
   193
{     int cnt;
alpar@1
   194
      if (!(lp->m == 0 || lp->valid))
alpar@1
   195
         xerror("glp_bf_update: basis factorization does not exist\n");
alpar@1
   196
#if 0 /* 15/XI-2009 */
alpar@1
   197
      cnt = (lp->m == 0 ? 0 : lp->bfd->upd_cnt);
alpar@1
   198
#else
alpar@1
   199
      cnt = (lp->m == 0 ? 0 : bfd_get_count(lp->bfd));
alpar@1
   200
#endif
alpar@1
   201
      return cnt;
alpar@1
   202
}
alpar@1
   203
alpar@1
   204
/***********************************************************************
alpar@1
   205
*  NAME
alpar@1
   206
*
alpar@1
   207
*  glp_get_bfcp - retrieve basis factorization control parameters
alpar@1
   208
*
alpar@1
   209
*  SYNOPSIS
alpar@1
   210
*
alpar@1
   211
*  void glp_get_bfcp(glp_prob *lp, glp_bfcp *parm);
alpar@1
   212
*
alpar@1
   213
*  DESCRIPTION
alpar@1
   214
*
alpar@1
   215
*  The routine glp_get_bfcp retrieves control parameters, which are
alpar@1
   216
*  used on computing and updating the basis factorization associated
alpar@1
   217
*  with the specified problem object.
alpar@1
   218
*
alpar@1
   219
*  Current values of control parameters are stored by the routine in
alpar@1
   220
*  a glp_bfcp structure, which the parameter parm points to. */
alpar@1
   221
alpar@1
   222
void glp_get_bfcp(glp_prob *lp, glp_bfcp *parm)
alpar@1
   223
{     glp_bfcp *bfcp = lp->bfcp;
alpar@1
   224
      if (bfcp == NULL)
alpar@1
   225
      {  parm->type = GLP_BF_FT;
alpar@1
   226
         parm->lu_size = 0;
alpar@1
   227
         parm->piv_tol = 0.10;
alpar@1
   228
         parm->piv_lim = 4;
alpar@1
   229
         parm->suhl = GLP_ON;
alpar@1
   230
         parm->eps_tol = 1e-15;
alpar@1
   231
         parm->max_gro = 1e+10;
alpar@1
   232
         parm->nfs_max = 100;
alpar@1
   233
         parm->upd_tol = 1e-6;
alpar@1
   234
         parm->nrs_max = 100;
alpar@1
   235
         parm->rs_size = 0;
alpar@1
   236
      }
alpar@1
   237
      else
alpar@1
   238
         memcpy(parm, bfcp, sizeof(glp_bfcp));
alpar@1
   239
      return;
alpar@1
   240
}
alpar@1
   241
alpar@1
   242
/***********************************************************************
alpar@1
   243
*  NAME
alpar@1
   244
*
alpar@1
   245
*  glp_set_bfcp - change basis factorization control parameters
alpar@1
   246
*
alpar@1
   247
*  SYNOPSIS
alpar@1
   248
*
alpar@1
   249
*  void glp_set_bfcp(glp_prob *lp, const glp_bfcp *parm);
alpar@1
   250
*
alpar@1
   251
*  DESCRIPTION
alpar@1
   252
*
alpar@1
   253
*  The routine glp_set_bfcp changes control parameters, which are used
alpar@1
   254
*  by internal GLPK routines in computing and updating the basis
alpar@1
   255
*  factorization associated with the specified problem object.
alpar@1
   256
*
alpar@1
   257
*  New values of the control parameters should be passed in a structure
alpar@1
   258
*  glp_bfcp, which the parameter parm points to.
alpar@1
   259
*
alpar@1
   260
*  The parameter parm can be specified as NULL, in which case all
alpar@1
   261
*  control parameters are reset to their default values. */
alpar@1
   262
alpar@1
   263
#if 0 /* 15/XI-2009 */
alpar@1
   264
static void copy_bfcp(glp_prob *lp)
alpar@1
   265
{     glp_bfcp _parm, *parm = &_parm;
alpar@1
   266
      BFD *bfd = lp->bfd;
alpar@1
   267
      glp_get_bfcp(lp, parm);
alpar@1
   268
      xassert(bfd != NULL);
alpar@1
   269
      bfd->type = parm->type;
alpar@1
   270
      bfd->lu_size = parm->lu_size;
alpar@1
   271
      bfd->piv_tol = parm->piv_tol;
alpar@1
   272
      bfd->piv_lim = parm->piv_lim;
alpar@1
   273
      bfd->suhl = parm->suhl;
alpar@1
   274
      bfd->eps_tol = parm->eps_tol;
alpar@1
   275
      bfd->max_gro = parm->max_gro;
alpar@1
   276
      bfd->nfs_max = parm->nfs_max;
alpar@1
   277
      bfd->upd_tol = parm->upd_tol;
alpar@1
   278
      bfd->nrs_max = parm->nrs_max;
alpar@1
   279
      bfd->rs_size = parm->rs_size;
alpar@1
   280
      return;
alpar@1
   281
}
alpar@1
   282
#else
alpar@1
   283
static void copy_bfcp(glp_prob *lp)
alpar@1
   284
{     glp_bfcp _parm, *parm = &_parm;
alpar@1
   285
      glp_get_bfcp(lp, parm);
alpar@1
   286
      bfd_set_parm(lp->bfd, parm);
alpar@1
   287
      return;
alpar@1
   288
}
alpar@1
   289
#endif
alpar@1
   290
alpar@1
   291
void glp_set_bfcp(glp_prob *lp, const glp_bfcp *parm)
alpar@1
   292
{     glp_bfcp *bfcp = lp->bfcp;
alpar@1
   293
      if (parm == NULL)
alpar@1
   294
      {  /* reset to default values */
alpar@1
   295
         if (bfcp != NULL)
alpar@1
   296
            xfree(bfcp), lp->bfcp = NULL;
alpar@1
   297
      }
alpar@1
   298
      else
alpar@1
   299
      {  /* set to specified values */
alpar@1
   300
         if (bfcp == NULL)
alpar@1
   301
            bfcp = lp->bfcp = xmalloc(sizeof(glp_bfcp));
alpar@1
   302
         memcpy(bfcp, parm, sizeof(glp_bfcp));
alpar@1
   303
         if (!(bfcp->type == GLP_BF_FT || bfcp->type == GLP_BF_BG ||
alpar@1
   304
               bfcp->type == GLP_BF_GR))
alpar@1
   305
            xerror("glp_set_bfcp: type = %d; invalid parameter\n",
alpar@1
   306
               bfcp->type);
alpar@1
   307
         if (bfcp->lu_size < 0)
alpar@1
   308
            xerror("glp_set_bfcp: lu_size = %d; invalid parameter\n",
alpar@1
   309
               bfcp->lu_size);
alpar@1
   310
         if (!(0.0 < bfcp->piv_tol && bfcp->piv_tol < 1.0))
alpar@1
   311
            xerror("glp_set_bfcp: piv_tol = %g; invalid parameter\n",
alpar@1
   312
               bfcp->piv_tol);
alpar@1
   313
         if (bfcp->piv_lim < 1)
alpar@1
   314
            xerror("glp_set_bfcp: piv_lim = %d; invalid parameter\n",
alpar@1
   315
               bfcp->piv_lim);
alpar@1
   316
         if (!(bfcp->suhl == GLP_ON || bfcp->suhl == GLP_OFF))
alpar@1
   317
            xerror("glp_set_bfcp: suhl = %d; invalid parameter\n",
alpar@1
   318
               bfcp->suhl);
alpar@1
   319
         if (!(0.0 <= bfcp->eps_tol && bfcp->eps_tol <= 1e-6))
alpar@1
   320
            xerror("glp_set_bfcp: eps_tol = %g; invalid parameter\n",
alpar@1
   321
               bfcp->eps_tol);
alpar@1
   322
         if (bfcp->max_gro < 1.0)
alpar@1
   323
            xerror("glp_set_bfcp: max_gro = %g; invalid parameter\n",
alpar@1
   324
               bfcp->max_gro);
alpar@1
   325
         if (!(1 <= bfcp->nfs_max && bfcp->nfs_max <= 32767))
alpar@1
   326
            xerror("glp_set_bfcp: nfs_max = %d; invalid parameter\n",
alpar@1
   327
               bfcp->nfs_max);
alpar@1
   328
         if (!(0.0 < bfcp->upd_tol && bfcp->upd_tol < 1.0))
alpar@1
   329
            xerror("glp_set_bfcp: upd_tol = %g; invalid parameter\n",
alpar@1
   330
               bfcp->upd_tol);
alpar@1
   331
         if (!(1 <= bfcp->nrs_max && bfcp->nrs_max <= 32767))
alpar@1
   332
            xerror("glp_set_bfcp: nrs_max = %d; invalid parameter\n",
alpar@1
   333
               bfcp->nrs_max);
alpar@1
   334
         if (bfcp->rs_size < 0)
alpar@1
   335
            xerror("glp_set_bfcp: rs_size = %d; invalid parameter\n",
alpar@1
   336
               bfcp->nrs_max);
alpar@1
   337
         if (bfcp->rs_size == 0)
alpar@1
   338
            bfcp->rs_size = 20 * bfcp->nrs_max;
alpar@1
   339
      }
alpar@1
   340
      if (lp->bfd != NULL) copy_bfcp(lp);
alpar@1
   341
      return;
alpar@1
   342
}
alpar@1
   343
alpar@1
   344
/***********************************************************************
alpar@1
   345
*  NAME
alpar@1
   346
*
alpar@1
   347
*  glp_get_bhead - retrieve the basis header information
alpar@1
   348
*
alpar@1
   349
*  SYNOPSIS
alpar@1
   350
*
alpar@1
   351
*  int glp_get_bhead(glp_prob *lp, int k);
alpar@1
   352
*
alpar@1
   353
*  DESCRIPTION
alpar@1
   354
*
alpar@1
   355
*  The routine glp_get_bhead returns the basis header information for
alpar@1
   356
*  the current basis associated with the specified problem object.
alpar@1
   357
*
alpar@1
   358
*  RETURNS
alpar@1
   359
*
alpar@1
   360
*  If xB[k], 1 <= k <= m, is i-th auxiliary variable (1 <= i <= m), the
alpar@1
   361
*  routine returns i. Otherwise, if xB[k] is j-th structural variable
alpar@1
   362
*  (1 <= j <= n), the routine returns m+j. Here m is the number of rows
alpar@1
   363
*  and n is the number of columns in the problem object. */
alpar@1
   364
alpar@1
   365
int glp_get_bhead(glp_prob *lp, int k)
alpar@1
   366
{     if (!(lp->m == 0 || lp->valid))
alpar@1
   367
         xerror("glp_get_bhead: basis factorization does not exist\n");
alpar@1
   368
      if (!(1 <= k && k <= lp->m))
alpar@1
   369
         xerror("glp_get_bhead: k = %d; index out of range\n", k);
alpar@1
   370
      return lp->head[k];
alpar@1
   371
}
alpar@1
   372
alpar@1
   373
/***********************************************************************
alpar@1
   374
*  NAME
alpar@1
   375
*
alpar@1
   376
*  glp_get_row_bind - retrieve row index in the basis header
alpar@1
   377
*
alpar@1
   378
*  SYNOPSIS
alpar@1
   379
*
alpar@1
   380
*  int glp_get_row_bind(glp_prob *lp, int i);
alpar@1
   381
*
alpar@1
   382
*  RETURNS
alpar@1
   383
*
alpar@1
   384
*  The routine glp_get_row_bind returns the index k of basic variable
alpar@1
   385
*  xB[k], 1 <= k <= m, which is i-th auxiliary variable, 1 <= i <= m,
alpar@1
   386
*  in the current basis associated with the specified problem object,
alpar@1
   387
*  where m is the number of rows. However, if i-th auxiliary variable
alpar@1
   388
*  is non-basic, the routine returns zero. */
alpar@1
   389
alpar@1
   390
int glp_get_row_bind(glp_prob *lp, int i)
alpar@1
   391
{     if (!(lp->m == 0 || lp->valid))
alpar@1
   392
         xerror("glp_get_row_bind: basis factorization does not exist\n"
alpar@1
   393
            );
alpar@1
   394
      if (!(1 <= i && i <= lp->m))
alpar@1
   395
         xerror("glp_get_row_bind: i = %d; row number out of range\n",
alpar@1
   396
            i);
alpar@1
   397
      return lp->row[i]->bind;
alpar@1
   398
}
alpar@1
   399
alpar@1
   400
/***********************************************************************
alpar@1
   401
*  NAME
alpar@1
   402
*
alpar@1
   403
*  glp_get_col_bind - retrieve column index in the basis header
alpar@1
   404
*
alpar@1
   405
*  SYNOPSIS
alpar@1
   406
*
alpar@1
   407
*  int glp_get_col_bind(glp_prob *lp, int j);
alpar@1
   408
*
alpar@1
   409
*  RETURNS
alpar@1
   410
*
alpar@1
   411
*  The routine glp_get_col_bind returns the index k of basic variable
alpar@1
   412
*  xB[k], 1 <= k <= m, which is j-th structural variable, 1 <= j <= n,
alpar@1
   413
*  in the current basis associated with the specified problem object,
alpar@1
   414
*  where m is the number of rows, n is the number of columns. However,
alpar@1
   415
*  if j-th structural variable is non-basic, the routine returns zero.*/
alpar@1
   416
alpar@1
   417
int glp_get_col_bind(glp_prob *lp, int j)
alpar@1
   418
{     if (!(lp->m == 0 || lp->valid))
alpar@1
   419
         xerror("glp_get_col_bind: basis factorization does not exist\n"
alpar@1
   420
            );
alpar@1
   421
      if (!(1 <= j && j <= lp->n))
alpar@1
   422
         xerror("glp_get_col_bind: j = %d; column number out of range\n"
alpar@1
   423
            , j);
alpar@1
   424
      return lp->col[j]->bind;
alpar@1
   425
}
alpar@1
   426
alpar@1
   427
/***********************************************************************
alpar@1
   428
*  NAME
alpar@1
   429
*
alpar@1
   430
*  glp_ftran - perform forward transformation (solve system B*x = b)
alpar@1
   431
*
alpar@1
   432
*  SYNOPSIS
alpar@1
   433
*
alpar@1
   434
*  void glp_ftran(glp_prob *lp, double x[]);
alpar@1
   435
*
alpar@1
   436
*  DESCRIPTION
alpar@1
   437
*
alpar@1
   438
*  The routine glp_ftran performs forward transformation, i.e. solves
alpar@1
   439
*  the system B*x = b, where B is the basis matrix corresponding to the
alpar@1
   440
*  current basis for the specified problem object, x is the vector of
alpar@1
   441
*  unknowns to be computed, b is the vector of right-hand sides.
alpar@1
   442
*
alpar@1
   443
*  On entry elements of the vector b should be stored in dense format
alpar@1
   444
*  in locations x[1], ..., x[m], where m is the number of rows. On exit
alpar@1
   445
*  the routine stores elements of the vector x in the same locations.
alpar@1
   446
*
alpar@1
   447
*  SCALING/UNSCALING
alpar@1
   448
*
alpar@1
   449
*  Let A~ = (I | -A) is the augmented constraint matrix of the original
alpar@1
   450
*  (unscaled) problem. In the scaled LP problem instead the matrix A the
alpar@1
   451
*  scaled matrix A" = R*A*S is actually used, so
alpar@1
   452
*
alpar@1
   453
*     A~" = (I | A") = (I | R*A*S) = (R*I*inv(R) | R*A*S) =
alpar@1
   454
*                                                                    (1)
alpar@1
   455
*         = R*(I | A)*S~ = R*A~*S~,
alpar@1
   456
*
alpar@1
   457
*  is the scaled augmented constraint matrix, where R and S are diagonal
alpar@1
   458
*  scaling matrices used to scale rows and columns of the matrix A, and
alpar@1
   459
*
alpar@1
   460
*     S~ = diag(inv(R) | S)                                          (2)
alpar@1
   461
*
alpar@1
   462
*  is an augmented diagonal scaling matrix.
alpar@1
   463
*
alpar@1
   464
*  By definition:
alpar@1
   465
*
alpar@1
   466
*     A~ = (B | N),                                                  (3)
alpar@1
   467
*
alpar@1
   468
*  where B is the basic matrix, which consists of basic columns of the
alpar@1
   469
*  augmented constraint matrix A~, and N is a matrix, which consists of
alpar@1
   470
*  non-basic columns of A~. From (1) it follows that:
alpar@1
   471
*
alpar@1
   472
*     A~" = (B" | N") = (R*B*SB | R*N*SN),                           (4)
alpar@1
   473
*
alpar@1
   474
*  where SB and SN are parts of the augmented scaling matrix S~, which
alpar@1
   475
*  correspond to basic and non-basic variables, respectively. Therefore
alpar@1
   476
*
alpar@1
   477
*     B" = R*B*SB,                                                   (5)
alpar@1
   478
*
alpar@1
   479
*  which is the scaled basis matrix. */
alpar@1
   480
alpar@1
   481
void glp_ftran(glp_prob *lp, double x[])
alpar@1
   482
{     int m = lp->m;
alpar@1
   483
      GLPROW **row = lp->row;
alpar@1
   484
      GLPCOL **col = lp->col;
alpar@1
   485
      int i, k;
alpar@1
   486
      /* B*x = b ===> (R*B*SB)*(inv(SB)*x) = R*b ===>
alpar@1
   487
         B"*x" = b", where b" = R*b, x = SB*x" */
alpar@1
   488
      if (!(m == 0 || lp->valid))
alpar@1
   489
         xerror("glp_ftran: basis factorization does not exist\n");
alpar@1
   490
      /* b" := R*b */
alpar@1
   491
      for (i = 1; i <= m; i++)
alpar@1
   492
         x[i] *= row[i]->rii;
alpar@1
   493
      /* x" := inv(B")*b" */
alpar@1
   494
      if (m > 0) bfd_ftran(lp->bfd, x);
alpar@1
   495
      /* x := SB*x" */
alpar@1
   496
      for (i = 1; i <= m; i++)
alpar@1
   497
      {  k = lp->head[i];
alpar@1
   498
         if (k <= m)
alpar@1
   499
            x[i] /= row[k]->rii;
alpar@1
   500
         else
alpar@1
   501
            x[i] *= col[k-m]->sjj;
alpar@1
   502
      }
alpar@1
   503
      return;
alpar@1
   504
}
alpar@1
   505
alpar@1
   506
/***********************************************************************
alpar@1
   507
*  NAME
alpar@1
   508
*
alpar@1
   509
*  glp_btran - perform backward transformation (solve system B'*x = b)
alpar@1
   510
*
alpar@1
   511
*  SYNOPSIS
alpar@1
   512
*
alpar@1
   513
*  void glp_btran(glp_prob *lp, double x[]);
alpar@1
   514
*
alpar@1
   515
*  DESCRIPTION
alpar@1
   516
*
alpar@1
   517
*  The routine glp_btran performs backward transformation, i.e. solves
alpar@1
   518
*  the system B'*x = b, where B' is a matrix transposed to the basis
alpar@1
   519
*  matrix corresponding to the current basis for the specified problem
alpar@1
   520
*  problem object, x is the vector of unknowns to be computed, b is the
alpar@1
   521
*  vector of right-hand sides.
alpar@1
   522
*
alpar@1
   523
*  On entry elements of the vector b should be stored in dense format
alpar@1
   524
*  in locations x[1], ..., x[m], where m is the number of rows. On exit
alpar@1
   525
*  the routine stores elements of the vector x in the same locations.
alpar@1
   526
*
alpar@1
   527
*  SCALING/UNSCALING
alpar@1
   528
*
alpar@1
   529
*  See comments to the routine glp_ftran. */
alpar@1
   530
alpar@1
   531
void glp_btran(glp_prob *lp, double x[])
alpar@1
   532
{     int m = lp->m;
alpar@1
   533
      GLPROW **row = lp->row;
alpar@1
   534
      GLPCOL **col = lp->col;
alpar@1
   535
      int i, k;
alpar@1
   536
      /* B'*x = b ===> (SB*B'*R)*(inv(R)*x) = SB*b ===>
alpar@1
   537
         (B")'*x" = b", where b" = SB*b, x = R*x" */
alpar@1
   538
      if (!(m == 0 || lp->valid))
alpar@1
   539
         xerror("glp_btran: basis factorization does not exist\n");
alpar@1
   540
      /* b" := SB*b */
alpar@1
   541
      for (i = 1; i <= m; i++)
alpar@1
   542
      {  k = lp->head[i];
alpar@1
   543
         if (k <= m)
alpar@1
   544
            x[i] /= row[k]->rii;
alpar@1
   545
         else
alpar@1
   546
            x[i] *= col[k-m]->sjj;
alpar@1
   547
      }
alpar@1
   548
      /* x" := inv[(B")']*b" */
alpar@1
   549
      if (m > 0) bfd_btran(lp->bfd, x);
alpar@1
   550
      /* x := R*x" */
alpar@1
   551
      for (i = 1; i <= m; i++)
alpar@1
   552
         x[i] *= row[i]->rii;
alpar@1
   553
      return;
alpar@1
   554
}
alpar@1
   555
alpar@1
   556
/***********************************************************************
alpar@1
   557
*  NAME
alpar@1
   558
*
alpar@1
   559
*  glp_warm_up - "warm up" LP basis
alpar@1
   560
*
alpar@1
   561
*  SYNOPSIS
alpar@1
   562
*
alpar@1
   563
*  int glp_warm_up(glp_prob *P);
alpar@1
   564
*
alpar@1
   565
*  DESCRIPTION
alpar@1
   566
*
alpar@1
   567
*  The routine glp_warm_up "warms up" the LP basis for the specified
alpar@1
   568
*  problem object using current statuses assigned to rows and columns
alpar@1
   569
*  (that is, to auxiliary and structural variables).
alpar@1
   570
*
alpar@1
   571
*  This operation includes computing factorization of the basis matrix
alpar@1
   572
*  (if it does not exist), computing primal and dual components of basic
alpar@1
   573
*  solution, and determining the solution status.
alpar@1
   574
*
alpar@1
   575
*  RETURNS
alpar@1
   576
*
alpar@1
   577
*  0  The operation has been successfully performed.
alpar@1
   578
*
alpar@1
   579
*  GLP_EBADB
alpar@1
   580
*     The basis matrix is invalid, i.e. the number of basic (auxiliary
alpar@1
   581
*     and structural) variables differs from the number of rows in the
alpar@1
   582
*     problem object.
alpar@1
   583
*
alpar@1
   584
*  GLP_ESING
alpar@1
   585
*     The basis matrix is singular within the working precision.
alpar@1
   586
*
alpar@1
   587
*  GLP_ECOND
alpar@1
   588
*     The basis matrix is ill-conditioned. */
alpar@1
   589
alpar@1
   590
int glp_warm_up(glp_prob *P)
alpar@1
   591
{     GLPROW *row;
alpar@1
   592
      GLPCOL *col;
alpar@1
   593
      GLPAIJ *aij;
alpar@1
   594
      int i, j, type, ret;
alpar@1
   595
      double eps, temp, *work;
alpar@1
   596
      /* invalidate basic solution */
alpar@1
   597
      P->pbs_stat = P->dbs_stat = GLP_UNDEF;
alpar@1
   598
      P->obj_val = 0.0;
alpar@1
   599
      P->some = 0;
alpar@1
   600
      for (i = 1; i <= P->m; i++)
alpar@1
   601
      {  row = P->row[i];
alpar@1
   602
         row->prim = row->dual = 0.0;
alpar@1
   603
      }
alpar@1
   604
      for (j = 1; j <= P->n; j++)
alpar@1
   605
      {  col = P->col[j];
alpar@1
   606
         col->prim = col->dual = 0.0;
alpar@1
   607
      }
alpar@1
   608
      /* compute the basis factorization, if necessary */
alpar@1
   609
      if (!glp_bf_exists(P))
alpar@1
   610
      {  ret = glp_factorize(P);
alpar@1
   611
         if (ret != 0) goto done;
alpar@1
   612
      }
alpar@1
   613
      /* allocate working array */
alpar@1
   614
      work = xcalloc(1+P->m, sizeof(double));
alpar@1
   615
      /* determine and store values of non-basic variables, compute
alpar@1
   616
         vector (- N * xN) */
alpar@1
   617
      for (i = 1; i <= P->m; i++)
alpar@1
   618
         work[i] = 0.0;
alpar@1
   619
      for (i = 1; i <= P->m; i++)
alpar@1
   620
      {  row = P->row[i];
alpar@1
   621
         if (row->stat == GLP_BS)
alpar@1
   622
            continue;
alpar@1
   623
         else if (row->stat == GLP_NL)
alpar@1
   624
            row->prim = row->lb;
alpar@1
   625
         else if (row->stat == GLP_NU)
alpar@1
   626
            row->prim = row->ub;
alpar@1
   627
         else if (row->stat == GLP_NF)
alpar@1
   628
            row->prim = 0.0;
alpar@1
   629
         else if (row->stat == GLP_NS)
alpar@1
   630
            row->prim = row->lb;
alpar@1
   631
         else
alpar@1
   632
            xassert(row != row);
alpar@1
   633
         /* N[j] is i-th column of matrix (I|-A) */
alpar@1
   634
         work[i] -= row->prim;
alpar@1
   635
      }
alpar@1
   636
      for (j = 1; j <= P->n; j++)
alpar@1
   637
      {  col = P->col[j];
alpar@1
   638
         if (col->stat == GLP_BS)
alpar@1
   639
            continue;
alpar@1
   640
         else if (col->stat == GLP_NL)
alpar@1
   641
            col->prim = col->lb;
alpar@1
   642
         else if (col->stat == GLP_NU)
alpar@1
   643
            col->prim = col->ub;
alpar@1
   644
         else if (col->stat == GLP_NF)
alpar@1
   645
            col->prim = 0.0;
alpar@1
   646
         else if (col->stat == GLP_NS)
alpar@1
   647
            col->prim = col->lb;
alpar@1
   648
         else
alpar@1
   649
            xassert(col != col);
alpar@1
   650
         /* N[j] is (m+j)-th column of matrix (I|-A) */
alpar@1
   651
         if (col->prim != 0.0)
alpar@1
   652
         {  for (aij = col->ptr; aij != NULL; aij = aij->c_next)
alpar@1
   653
               work[aij->row->i] += aij->val * col->prim;
alpar@1
   654
         }
alpar@1
   655
      }
alpar@1
   656
      /* compute vector of basic variables xB = - inv(B) * N * xN */
alpar@1
   657
      glp_ftran(P, work);
alpar@1
   658
      /* store values of basic variables, check primal feasibility */
alpar@1
   659
      P->pbs_stat = GLP_FEAS;
alpar@1
   660
      for (i = 1; i <= P->m; i++)
alpar@1
   661
      {  row = P->row[i];
alpar@1
   662
         if (row->stat != GLP_BS)
alpar@1
   663
            continue;
alpar@1
   664
         row->prim = work[row->bind];
alpar@1
   665
         type = row->type;
alpar@1
   666
         if (type == GLP_LO || type == GLP_DB || type == GLP_FX)
alpar@1
   667
         {  eps = 1e-6 + 1e-9 * fabs(row->lb);
alpar@1
   668
            if (row->prim < row->lb - eps)
alpar@1
   669
               P->pbs_stat = GLP_INFEAS;
alpar@1
   670
         }
alpar@1
   671
         if (type == GLP_UP || type == GLP_DB || type == GLP_FX)
alpar@1
   672
         {  eps = 1e-6 + 1e-9 * fabs(row->ub);
alpar@1
   673
            if (row->prim > row->ub + eps)
alpar@1
   674
               P->pbs_stat = GLP_INFEAS;
alpar@1
   675
         }
alpar@1
   676
      }
alpar@1
   677
      for (j = 1; j <= P->n; j++)
alpar@1
   678
      {  col = P->col[j];
alpar@1
   679
         if (col->stat != GLP_BS)
alpar@1
   680
            continue;
alpar@1
   681
         col->prim = work[col->bind];
alpar@1
   682
         type = col->type;
alpar@1
   683
         if (type == GLP_LO || type == GLP_DB || type == GLP_FX)
alpar@1
   684
         {  eps = 1e-6 + 1e-9 * fabs(col->lb);
alpar@1
   685
            if (col->prim < col->lb - eps)
alpar@1
   686
               P->pbs_stat = GLP_INFEAS;
alpar@1
   687
         }
alpar@1
   688
         if (type == GLP_UP || type == GLP_DB || type == GLP_FX)
alpar@1
   689
         {  eps = 1e-6 + 1e-9 * fabs(col->ub);
alpar@1
   690
            if (col->prim > col->ub + eps)
alpar@1
   691
               P->pbs_stat = GLP_INFEAS;
alpar@1
   692
         }
alpar@1
   693
      }
alpar@1
   694
      /* compute value of the objective function */
alpar@1
   695
      P->obj_val = P->c0;
alpar@1
   696
      for (j = 1; j <= P->n; j++)
alpar@1
   697
      {  col = P->col[j];
alpar@1
   698
         P->obj_val += col->coef * col->prim;
alpar@1
   699
      }
alpar@1
   700
      /* build vector cB of objective coefficients at basic variables */
alpar@1
   701
      for (i = 1; i <= P->m; i++)
alpar@1
   702
         work[i] = 0.0;
alpar@1
   703
      for (j = 1; j <= P->n; j++)
alpar@1
   704
      {  col = P->col[j];
alpar@1
   705
         if (col->stat == GLP_BS)
alpar@1
   706
            work[col->bind] = col->coef;
alpar@1
   707
      }
alpar@1
   708
      /* compute vector of simplex multipliers pi = inv(B') * cB */
alpar@1
   709
      glp_btran(P, work);
alpar@1
   710
      /* compute and store reduced costs of non-basic variables d[j] =
alpar@1
   711
         c[j] - N'[j] * pi, check dual feasibility */
alpar@1
   712
      P->dbs_stat = GLP_FEAS;
alpar@1
   713
      for (i = 1; i <= P->m; i++)
alpar@1
   714
      {  row = P->row[i];
alpar@1
   715
         if (row->stat == GLP_BS)
alpar@1
   716
         {  row->dual = 0.0;
alpar@1
   717
            continue;
alpar@1
   718
         }
alpar@1
   719
         /* N[j] is i-th column of matrix (I|-A) */
alpar@1
   720
         row->dual = - work[i];
alpar@1
   721
         type = row->type;
alpar@1
   722
         temp = (P->dir == GLP_MIN ? + row->dual : - row->dual);
alpar@1
   723
         if ((type == GLP_FR || type == GLP_LO) && temp < -1e-5 ||
alpar@1
   724
             (type == GLP_FR || type == GLP_UP) && temp > +1e-5)
alpar@1
   725
            P->dbs_stat = GLP_INFEAS;
alpar@1
   726
      }
alpar@1
   727
      for (j = 1; j <= P->n; j++)
alpar@1
   728
      {  col = P->col[j];
alpar@1
   729
         if (col->stat == GLP_BS)
alpar@1
   730
         {  col->dual = 0.0;
alpar@1
   731
            continue;
alpar@1
   732
         }
alpar@1
   733
         /* N[j] is (m+j)-th column of matrix (I|-A) */
alpar@1
   734
         col->dual = col->coef;
alpar@1
   735
         for (aij = col->ptr; aij != NULL; aij = aij->c_next)
alpar@1
   736
            col->dual += aij->val * work[aij->row->i];
alpar@1
   737
         type = col->type;
alpar@1
   738
         temp = (P->dir == GLP_MIN ? + col->dual : - col->dual);
alpar@1
   739
         if ((type == GLP_FR || type == GLP_LO) && temp < -1e-5 ||
alpar@1
   740
             (type == GLP_FR || type == GLP_UP) && temp > +1e-5)
alpar@1
   741
            P->dbs_stat = GLP_INFEAS;
alpar@1
   742
      }
alpar@1
   743
      /* free working array */
alpar@1
   744
      xfree(work);
alpar@1
   745
      ret = 0;
alpar@1
   746
done: return ret;
alpar@1
   747
}
alpar@1
   748
alpar@1
   749
/***********************************************************************
alpar@1
   750
*  NAME
alpar@1
   751
*
alpar@1
   752
*  glp_eval_tab_row - compute row of the simplex tableau
alpar@1
   753
*
alpar@1
   754
*  SYNOPSIS
alpar@1
   755
*
alpar@1
   756
*  int glp_eval_tab_row(glp_prob *lp, int k, int ind[], double val[]);
alpar@1
   757
*
alpar@1
   758
*  DESCRIPTION
alpar@1
   759
*
alpar@1
   760
*  The routine glp_eval_tab_row computes a row of the current simplex
alpar@1
   761
*  tableau for the basic variable, which is specified by the number k:
alpar@1
   762
*  if 1 <= k <= m, x[k] is k-th auxiliary variable; if m+1 <= k <= m+n,
alpar@1
   763
*  x[k] is (k-m)-th structural variable, where m is number of rows, and
alpar@1
   764
*  n is number of columns. The current basis must be available.
alpar@1
   765
*
alpar@1
   766
*  The routine stores column indices and numerical values of non-zero
alpar@1
   767
*  elements of the computed row using sparse format to the locations
alpar@1
   768
*  ind[1], ..., ind[len] and val[1], ..., val[len], respectively, where
alpar@1
   769
*  0 <= len <= n is number of non-zeros returned on exit.
alpar@1
   770
*
alpar@1
   771
*  Element indices stored in the array ind have the same sense as the
alpar@1
   772
*  index k, i.e. indices 1 to m denote auxiliary variables and indices
alpar@1
   773
*  m+1 to m+n denote structural ones (all these variables are obviously
alpar@1
   774
*  non-basic by definition).
alpar@1
   775
*
alpar@1
   776
*  The computed row shows how the specified basic variable x[k] = xB[i]
alpar@1
   777
*  depends on non-basic variables:
alpar@1
   778
*
alpar@1
   779
*     xB[i] = alfa[i,1]*xN[1] + alfa[i,2]*xN[2] + ... + alfa[i,n]*xN[n],
alpar@1
   780
*
alpar@1
   781
*  where alfa[i,j] are elements of the simplex table row, xN[j] are
alpar@1
   782
*  non-basic (auxiliary and structural) variables.
alpar@1
   783
*
alpar@1
   784
*  RETURNS
alpar@1
   785
*
alpar@1
   786
*  The routine returns number of non-zero elements in the simplex table
alpar@1
   787
*  row stored in the arrays ind and val.
alpar@1
   788
*
alpar@1
   789
*  BACKGROUND
alpar@1
   790
*
alpar@1
   791
*  The system of equality constraints of the LP problem is:
alpar@1
   792
*
alpar@1
   793
*     xR = A * xS,                                                   (1)
alpar@1
   794
*
alpar@1
   795
*  where xR is the vector of auxliary variables, xS is the vector of
alpar@1
   796
*  structural variables, A is the matrix of constraint coefficients.
alpar@1
   797
*
alpar@1
   798
*  The system (1) can be written in homogenous form as follows:
alpar@1
   799
*
alpar@1
   800
*     A~ * x = 0,                                                    (2)
alpar@1
   801
*
alpar@1
   802
*  where A~ = (I | -A) is the augmented constraint matrix (has m rows
alpar@1
   803
*  and m+n columns), x = (xR | xS) is the vector of all (auxiliary and
alpar@1
   804
*  structural) variables.
alpar@1
   805
*
alpar@1
   806
*  By definition for the current basis we have:
alpar@1
   807
*
alpar@1
   808
*     A~ = (B | N),                                                  (3)
alpar@1
   809
*
alpar@1
   810
*  where B is the basis matrix. Thus, the system (2) can be written as:
alpar@1
   811
*
alpar@1
   812
*     B * xB + N * xN = 0.                                           (4)
alpar@1
   813
*
alpar@1
   814
*  From (4) it follows that:
alpar@1
   815
*
alpar@1
   816
*     xB = A^ * xN,                                                  (5)
alpar@1
   817
*
alpar@1
   818
*  where the matrix
alpar@1
   819
*
alpar@1
   820
*     A^ = - inv(B) * N                                              (6)
alpar@1
   821
*
alpar@1
   822
*  is called the simplex table.
alpar@1
   823
*
alpar@1
   824
*  It is understood that i-th row of the simplex table is:
alpar@1
   825
*
alpar@1
   826
*     e * A^ = - e * inv(B) * N,                                     (7)
alpar@1
   827
*
alpar@1
   828
*  where e is a unity vector with e[i] = 1.
alpar@1
   829
*
alpar@1
   830
*  To compute i-th row of the simplex table the routine first computes
alpar@1
   831
*  i-th row of the inverse:
alpar@1
   832
*
alpar@1
   833
*     rho = inv(B') * e,                                             (8)
alpar@1
   834
*
alpar@1
   835
*  where B' is a matrix transposed to B, and then computes elements of
alpar@1
   836
*  i-th row of the simplex table as scalar products:
alpar@1
   837
*
alpar@1
   838
*     alfa[i,j] = - rho * N[j]   for all j,                          (9)
alpar@1
   839
*
alpar@1
   840
*  where N[j] is a column of the augmented constraint matrix A~, which
alpar@1
   841
*  corresponds to some non-basic auxiliary or structural variable. */
alpar@1
   842
alpar@1
   843
int glp_eval_tab_row(glp_prob *lp, int k, int ind[], double val[])
alpar@1
   844
{     int m = lp->m;
alpar@1
   845
      int n = lp->n;
alpar@1
   846
      int i, t, len, lll, *iii;
alpar@1
   847
      double alfa, *rho, *vvv;
alpar@1
   848
      if (!(m == 0 || lp->valid))
alpar@1
   849
         xerror("glp_eval_tab_row: basis factorization does not exist\n"
alpar@1
   850
            );
alpar@1
   851
      if (!(1 <= k && k <= m+n))
alpar@1
   852
         xerror("glp_eval_tab_row: k = %d; variable number out of range"
alpar@1
   853
            , k);
alpar@1
   854
      /* determine xB[i] which corresponds to x[k] */
alpar@1
   855
      if (k <= m)
alpar@1
   856
         i = glp_get_row_bind(lp, k);
alpar@1
   857
      else
alpar@1
   858
         i = glp_get_col_bind(lp, k-m);
alpar@1
   859
      if (i == 0)
alpar@1
   860
         xerror("glp_eval_tab_row: k = %d; variable must be basic", k);
alpar@1
   861
      xassert(1 <= i && i <= m);
alpar@1
   862
      /* allocate working arrays */
alpar@1
   863
      rho = xcalloc(1+m, sizeof(double));
alpar@1
   864
      iii = xcalloc(1+m, sizeof(int));
alpar@1
   865
      vvv = xcalloc(1+m, sizeof(double));
alpar@1
   866
      /* compute i-th row of the inverse; see (8) */
alpar@1
   867
      for (t = 1; t <= m; t++) rho[t] = 0.0;
alpar@1
   868
      rho[i] = 1.0;
alpar@1
   869
      glp_btran(lp, rho);
alpar@1
   870
      /* compute i-th row of the simplex table */
alpar@1
   871
      len = 0;
alpar@1
   872
      for (k = 1; k <= m+n; k++)
alpar@1
   873
      {  if (k <= m)
alpar@1
   874
         {  /* x[k] is auxiliary variable, so N[k] is a unity column */
alpar@1
   875
            if (glp_get_row_stat(lp, k) == GLP_BS) continue;
alpar@1
   876
            /* compute alfa[i,j]; see (9) */
alpar@1
   877
            alfa = - rho[k];
alpar@1
   878
         }
alpar@1
   879
         else
alpar@1
   880
         {  /* x[k] is structural variable, so N[k] is a column of the
alpar@1
   881
               original constraint matrix A with negative sign */
alpar@1
   882
            if (glp_get_col_stat(lp, k-m) == GLP_BS) continue;
alpar@1
   883
            /* compute alfa[i,j]; see (9) */
alpar@1
   884
            lll = glp_get_mat_col(lp, k-m, iii, vvv);
alpar@1
   885
            alfa = 0.0;
alpar@1
   886
            for (t = 1; t <= lll; t++) alfa += rho[iii[t]] * vvv[t];
alpar@1
   887
         }
alpar@1
   888
         /* store alfa[i,j] */
alpar@1
   889
         if (alfa != 0.0) len++, ind[len] = k, val[len] = alfa;
alpar@1
   890
      }
alpar@1
   891
      xassert(len <= n);
alpar@1
   892
      /* free working arrays */
alpar@1
   893
      xfree(rho);
alpar@1
   894
      xfree(iii);
alpar@1
   895
      xfree(vvv);
alpar@1
   896
      /* return to the calling program */
alpar@1
   897
      return len;
alpar@1
   898
}
alpar@1
   899
alpar@1
   900
/***********************************************************************
alpar@1
   901
*  NAME
alpar@1
   902
*
alpar@1
   903
*  glp_eval_tab_col - compute column of the simplex tableau
alpar@1
   904
*
alpar@1
   905
*  SYNOPSIS
alpar@1
   906
*
alpar@1
   907
*  int glp_eval_tab_col(glp_prob *lp, int k, int ind[], double val[]);
alpar@1
   908
*
alpar@1
   909
*  DESCRIPTION
alpar@1
   910
*
alpar@1
   911
*  The routine glp_eval_tab_col computes a column of the current simplex
alpar@1
   912
*  table for the non-basic variable, which is specified by the number k:
alpar@1
   913
*  if 1 <= k <= m, x[k] is k-th auxiliary variable; if m+1 <= k <= m+n,
alpar@1
   914
*  x[k] is (k-m)-th structural variable, where m is number of rows, and
alpar@1
   915
*  n is number of columns. The current basis must be available.
alpar@1
   916
*
alpar@1
   917
*  The routine stores row indices and numerical values of non-zero
alpar@1
   918
*  elements of the computed column using sparse format to the locations
alpar@1
   919
*  ind[1], ..., ind[len] and val[1], ..., val[len] respectively, where
alpar@1
   920
*  0 <= len <= m is number of non-zeros returned on exit.
alpar@1
   921
*
alpar@1
   922
*  Element indices stored in the array ind have the same sense as the
alpar@1
   923
*  index k, i.e. indices 1 to m denote auxiliary variables and indices
alpar@1
   924
*  m+1 to m+n denote structural ones (all these variables are obviously
alpar@1
   925
*  basic by the definition).
alpar@1
   926
*
alpar@1
   927
*  The computed column shows how basic variables depend on the specified
alpar@1
   928
*  non-basic variable x[k] = xN[j]:
alpar@1
   929
*
alpar@1
   930
*     xB[1] = ... + alfa[1,j]*xN[j] + ...
alpar@1
   931
*     xB[2] = ... + alfa[2,j]*xN[j] + ...
alpar@1
   932
*              . . . . . .
alpar@1
   933
*     xB[m] = ... + alfa[m,j]*xN[j] + ...
alpar@1
   934
*
alpar@1
   935
*  where alfa[i,j] are elements of the simplex table column, xB[i] are
alpar@1
   936
*  basic (auxiliary and structural) variables.
alpar@1
   937
*
alpar@1
   938
*  RETURNS
alpar@1
   939
*
alpar@1
   940
*  The routine returns number of non-zero elements in the simplex table
alpar@1
   941
*  column stored in the arrays ind and val.
alpar@1
   942
*
alpar@1
   943
*  BACKGROUND
alpar@1
   944
*
alpar@1
   945
*  As it was explained in comments to the routine glp_eval_tab_row (see
alpar@1
   946
*  above) the simplex table is the following matrix:
alpar@1
   947
*
alpar@1
   948
*     A^ = - inv(B) * N.                                             (1)
alpar@1
   949
*
alpar@1
   950
*  Therefore j-th column of the simplex table is:
alpar@1
   951
*
alpar@1
   952
*     A^ * e = - inv(B) * N * e = - inv(B) * N[j],                   (2)
alpar@1
   953
*
alpar@1
   954
*  where e is a unity vector with e[j] = 1, B is the basis matrix, N[j]
alpar@1
   955
*  is a column of the augmented constraint matrix A~, which corresponds
alpar@1
   956
*  to the given non-basic auxiliary or structural variable. */
alpar@1
   957
alpar@1
   958
int glp_eval_tab_col(glp_prob *lp, int k, int ind[], double val[])
alpar@1
   959
{     int m = lp->m;
alpar@1
   960
      int n = lp->n;
alpar@1
   961
      int t, len, stat;
alpar@1
   962
      double *col;
alpar@1
   963
      if (!(m == 0 || lp->valid))
alpar@1
   964
         xerror("glp_eval_tab_col: basis factorization does not exist\n"
alpar@1
   965
            );
alpar@1
   966
      if (!(1 <= k && k <= m+n))
alpar@1
   967
         xerror("glp_eval_tab_col: k = %d; variable number out of range"
alpar@1
   968
            , k);
alpar@1
   969
      if (k <= m)
alpar@1
   970
         stat = glp_get_row_stat(lp, k);
alpar@1
   971
      else
alpar@1
   972
         stat = glp_get_col_stat(lp, k-m);
alpar@1
   973
      if (stat == GLP_BS)
alpar@1
   974
         xerror("glp_eval_tab_col: k = %d; variable must be non-basic",
alpar@1
   975
            k);
alpar@1
   976
      /* obtain column N[k] with negative sign */
alpar@1
   977
      col = xcalloc(1+m, sizeof(double));
alpar@1
   978
      for (t = 1; t <= m; t++) col[t] = 0.0;
alpar@1
   979
      if (k <= m)
alpar@1
   980
      {  /* x[k] is auxiliary variable, so N[k] is a unity column */
alpar@1
   981
         col[k] = -1.0;
alpar@1
   982
      }
alpar@1
   983
      else
alpar@1
   984
      {  /* x[k] is structural variable, so N[k] is a column of the
alpar@1
   985
            original constraint matrix A with negative sign */
alpar@1
   986
         len = glp_get_mat_col(lp, k-m, ind, val);
alpar@1
   987
         for (t = 1; t <= len; t++) col[ind[t]] = val[t];
alpar@1
   988
      }
alpar@1
   989
      /* compute column of the simplex table, which corresponds to the
alpar@1
   990
         specified non-basic variable x[k] */
alpar@1
   991
      glp_ftran(lp, col);
alpar@1
   992
      len = 0;
alpar@1
   993
      for (t = 1; t <= m; t++)
alpar@1
   994
      {  if (col[t] != 0.0)
alpar@1
   995
         {  len++;
alpar@1
   996
            ind[len] = glp_get_bhead(lp, t);
alpar@1
   997
            val[len] = col[t];
alpar@1
   998
         }
alpar@1
   999
      }
alpar@1
  1000
      xfree(col);
alpar@1
  1001
      /* return to the calling program */
alpar@1
  1002
      return len;
alpar@1
  1003
}
alpar@1
  1004
alpar@1
  1005
/***********************************************************************
alpar@1
  1006
*  NAME
alpar@1
  1007
*
alpar@1
  1008
*  glp_transform_row - transform explicitly specified row
alpar@1
  1009
*
alpar@1
  1010
*  SYNOPSIS
alpar@1
  1011
*
alpar@1
  1012
*  int glp_transform_row(glp_prob *P, int len, int ind[], double val[]);
alpar@1
  1013
*
alpar@1
  1014
*  DESCRIPTION
alpar@1
  1015
*
alpar@1
  1016
*  The routine glp_transform_row performs the same operation as the
alpar@1
  1017
*  routine glp_eval_tab_row with exception that the row to be
alpar@1
  1018
*  transformed is specified explicitly as a sparse vector.
alpar@1
  1019
*
alpar@1
  1020
*  The explicitly specified row may be thought as a linear form:
alpar@1
  1021
*
alpar@1
  1022
*     x = a[1]*x[m+1] + a[2]*x[m+2] + ... + a[n]*x[m+n],             (1)
alpar@1
  1023
*
alpar@1
  1024
*  where x is an auxiliary variable for this row, a[j] are coefficients
alpar@1
  1025
*  of the linear form, x[m+j] are structural variables.
alpar@1
  1026
*
alpar@1
  1027
*  On entry column indices and numerical values of non-zero elements of
alpar@1
  1028
*  the row should be stored in locations ind[1], ..., ind[len] and
alpar@1
  1029
*  val[1], ..., val[len], where len is the number of non-zero elements.
alpar@1
  1030
*
alpar@1
  1031
*  This routine uses the system of equality constraints and the current
alpar@1
  1032
*  basis in order to express the auxiliary variable x in (1) through the
alpar@1
  1033
*  current non-basic variables (as if the transformed row were added to
alpar@1
  1034
*  the problem object and its auxiliary variable were basic), i.e. the
alpar@1
  1035
*  resultant row has the form:
alpar@1
  1036
*
alpar@1
  1037
*     x = alfa[1]*xN[1] + alfa[2]*xN[2] + ... + alfa[n]*xN[n],       (2)
alpar@1
  1038
*
alpar@1
  1039
*  where xN[j] are non-basic (auxiliary or structural) variables, n is
alpar@1
  1040
*  the number of columns in the LP problem object.
alpar@1
  1041
*
alpar@1
  1042
*  On exit the routine stores indices and numerical values of non-zero
alpar@1
  1043
*  elements of the resultant row (2) in locations ind[1], ..., ind[len']
alpar@1
  1044
*  and val[1], ..., val[len'], where 0 <= len' <= n is the number of
alpar@1
  1045
*  non-zero elements in the resultant row returned by the routine. Note
alpar@1
  1046
*  that indices (numbers) of non-basic variables stored in the array ind
alpar@1
  1047
*  correspond to original ordinal numbers of variables: indices 1 to m
alpar@1
  1048
*  mean auxiliary variables and indices m+1 to m+n mean structural ones.
alpar@1
  1049
*
alpar@1
  1050
*  RETURNS
alpar@1
  1051
*
alpar@1
  1052
*  The routine returns len', which is the number of non-zero elements in
alpar@1
  1053
*  the resultant row stored in the arrays ind and val.
alpar@1
  1054
*
alpar@1
  1055
*  BACKGROUND
alpar@1
  1056
*
alpar@1
  1057
*  The explicitly specified row (1) is transformed in the same way as it
alpar@1
  1058
*  were the objective function row.
alpar@1
  1059
*
alpar@1
  1060
*  From (1) it follows that:
alpar@1
  1061
*
alpar@1
  1062
*     x = aB * xB + aN * xN,                                         (3)
alpar@1
  1063
*
alpar@1
  1064
*  where xB is the vector of basic variables, xN is the vector of
alpar@1
  1065
*  non-basic variables.
alpar@1
  1066
*
alpar@1
  1067
*  The simplex table, which corresponds to the current basis, is:
alpar@1
  1068
*
alpar@1
  1069
*     xB = [-inv(B) * N] * xN.                                       (4)
alpar@1
  1070
*
alpar@1
  1071
*  Therefore substituting xB from (4) to (3) we have:
alpar@1
  1072
*
alpar@1
  1073
*     x = aB * [-inv(B) * N] * xN + aN * xN =
alpar@1
  1074
*                                                                    (5)
alpar@1
  1075
*       = rho * (-N) * xN + aN * xN = alfa * xN,
alpar@1
  1076
*
alpar@1
  1077
*  where:
alpar@1
  1078
*
alpar@1
  1079
*     rho = inv(B') * aB,                                            (6)
alpar@1
  1080
*
alpar@1
  1081
*  and
alpar@1
  1082
*
alpar@1
  1083
*     alfa = aN + rho * (-N)                                         (7)
alpar@1
  1084
*
alpar@1
  1085
*  is the resultant row computed by the routine. */
alpar@1
  1086
alpar@1
  1087
int glp_transform_row(glp_prob *P, int len, int ind[], double val[])
alpar@1
  1088
{     int i, j, k, m, n, t, lll, *iii;
alpar@1
  1089
      double alfa, *a, *aB, *rho, *vvv;
alpar@1
  1090
      if (!glp_bf_exists(P))
alpar@1
  1091
         xerror("glp_transform_row: basis factorization does not exist "
alpar@1
  1092
            "\n");
alpar@1
  1093
      m = glp_get_num_rows(P);
alpar@1
  1094
      n = glp_get_num_cols(P);
alpar@1
  1095
      /* unpack the row to be transformed to the array a */
alpar@1
  1096
      a = xcalloc(1+n, sizeof(double));
alpar@1
  1097
      for (j = 1; j <= n; j++) a[j] = 0.0;
alpar@1
  1098
      if (!(0 <= len && len <= n))
alpar@1
  1099
         xerror("glp_transform_row: len = %d; invalid row length\n",
alpar@1
  1100
            len);
alpar@1
  1101
      for (t = 1; t <= len; t++)
alpar@1
  1102
      {  j = ind[t];
alpar@1
  1103
         if (!(1 <= j && j <= n))
alpar@1
  1104
            xerror("glp_transform_row: ind[%d] = %d; column index out o"
alpar@1
  1105
               "f range\n", t, j);
alpar@1
  1106
         if (val[t] == 0.0)
alpar@1
  1107
            xerror("glp_transform_row: val[%d] = 0; zero coefficient no"
alpar@1
  1108
               "t allowed\n", t);
alpar@1
  1109
         if (a[j] != 0.0)
alpar@1
  1110
            xerror("glp_transform_row: ind[%d] = %d; duplicate column i"
alpar@1
  1111
               "ndices not allowed\n", t, j);
alpar@1
  1112
         a[j] = val[t];
alpar@1
  1113
      }
alpar@1
  1114
      /* construct the vector aB */
alpar@1
  1115
      aB = xcalloc(1+m, sizeof(double));
alpar@1
  1116
      for (i = 1; i <= m; i++)
alpar@1
  1117
      {  k = glp_get_bhead(P, i);
alpar@1
  1118
         /* xB[i] is k-th original variable */
alpar@1
  1119
         xassert(1 <= k && k <= m+n);
alpar@1
  1120
         aB[i] = (k <= m ? 0.0 : a[k-m]);
alpar@1
  1121
      }
alpar@1
  1122
      /* solve the system B'*rho = aB to compute the vector rho */
alpar@1
  1123
      rho = aB, glp_btran(P, rho);
alpar@1
  1124
      /* compute coefficients at non-basic auxiliary variables */
alpar@1
  1125
      len = 0;
alpar@1
  1126
      for (i = 1; i <= m; i++)
alpar@1
  1127
      {  if (glp_get_row_stat(P, i) != GLP_BS)
alpar@1
  1128
         {  alfa = - rho[i];
alpar@1
  1129
            if (alfa != 0.0)
alpar@1
  1130
            {  len++;
alpar@1
  1131
               ind[len] = i;
alpar@1
  1132
               val[len] = alfa;
alpar@1
  1133
            }
alpar@1
  1134
         }
alpar@1
  1135
      }
alpar@1
  1136
      /* compute coefficients at non-basic structural variables */
alpar@1
  1137
      iii = xcalloc(1+m, sizeof(int));
alpar@1
  1138
      vvv = xcalloc(1+m, sizeof(double));
alpar@1
  1139
      for (j = 1; j <= n; j++)
alpar@1
  1140
      {  if (glp_get_col_stat(P, j) != GLP_BS)
alpar@1
  1141
         {  alfa = a[j];
alpar@1
  1142
            lll = glp_get_mat_col(P, j, iii, vvv);
alpar@1
  1143
            for (t = 1; t <= lll; t++) alfa += vvv[t] * rho[iii[t]];
alpar@1
  1144
            if (alfa != 0.0)
alpar@1
  1145
            {  len++;
alpar@1
  1146
               ind[len] = m+j;
alpar@1
  1147
               val[len] = alfa;
alpar@1
  1148
            }
alpar@1
  1149
         }
alpar@1
  1150
      }
alpar@1
  1151
      xassert(len <= n);
alpar@1
  1152
      xfree(iii);
alpar@1
  1153
      xfree(vvv);
alpar@1
  1154
      xfree(aB);
alpar@1
  1155
      xfree(a);
alpar@1
  1156
      return len;
alpar@1
  1157
}
alpar@1
  1158
alpar@1
  1159
/***********************************************************************
alpar@1
  1160
*  NAME
alpar@1
  1161
*
alpar@1
  1162
*  glp_transform_col - transform explicitly specified column
alpar@1
  1163
*
alpar@1
  1164
*  SYNOPSIS
alpar@1
  1165
*
alpar@1
  1166
*  int glp_transform_col(glp_prob *P, int len, int ind[], double val[]);
alpar@1
  1167
*
alpar@1
  1168
*  DESCRIPTION
alpar@1
  1169
*
alpar@1
  1170
*  The routine glp_transform_col performs the same operation as the
alpar@1
  1171
*  routine glp_eval_tab_col with exception that the column to be
alpar@1
  1172
*  transformed is specified explicitly as a sparse vector.
alpar@1
  1173
*
alpar@1
  1174
*  The explicitly specified column may be thought as if it were added
alpar@1
  1175
*  to the original system of equality constraints:
alpar@1
  1176
*
alpar@1
  1177
*     x[1] = a[1,1]*x[m+1] + ... + a[1,n]*x[m+n] + a[1]*x
alpar@1
  1178
*     x[2] = a[2,1]*x[m+1] + ... + a[2,n]*x[m+n] + a[2]*x            (1)
alpar@1
  1179
*        .  .  .  .  .  .  .  .  .  .  .  .  .  .  .
alpar@1
  1180
*     x[m] = a[m,1]*x[m+1] + ... + a[m,n]*x[m+n] + a[m]*x
alpar@1
  1181
*
alpar@1
  1182
*  where x[i] are auxiliary variables, x[m+j] are structural variables,
alpar@1
  1183
*  x is a structural variable for the explicitly specified column, a[i]
alpar@1
  1184
*  are constraint coefficients for x.
alpar@1
  1185
*
alpar@1
  1186
*  On entry row indices and numerical values of non-zero elements of
alpar@1
  1187
*  the column should be stored in locations ind[1], ..., ind[len] and
alpar@1
  1188
*  val[1], ..., val[len], where len is the number of non-zero elements.
alpar@1
  1189
*
alpar@1
  1190
*  This routine uses the system of equality constraints and the current
alpar@1
  1191
*  basis in order to express the current basic variables through the
alpar@1
  1192
*  structural variable x in (1) (as if the transformed column were added
alpar@1
  1193
*  to the problem object and the variable x were non-basic), i.e. the
alpar@1
  1194
*  resultant column has the form:
alpar@1
  1195
*
alpar@1
  1196
*     xB[1] = ... + alfa[1]*x
alpar@1
  1197
*     xB[2] = ... + alfa[2]*x                                        (2)
alpar@1
  1198
*        .  .  .  .  .  .
alpar@1
  1199
*     xB[m] = ... + alfa[m]*x
alpar@1
  1200
*
alpar@1
  1201
*  where xB are basic (auxiliary and structural) variables, m is the
alpar@1
  1202
*  number of rows in the problem object.
alpar@1
  1203
*
alpar@1
  1204
*  On exit the routine stores indices and numerical values of non-zero
alpar@1
  1205
*  elements of the resultant column (2) in locations ind[1], ...,
alpar@1
  1206
*  ind[len'] and val[1], ..., val[len'], where 0 <= len' <= m is the
alpar@1
  1207
*  number of non-zero element in the resultant column returned by the
alpar@1
  1208
*  routine. Note that indices (numbers) of basic variables stored in
alpar@1
  1209
*  the array ind correspond to original ordinal numbers of variables:
alpar@1
  1210
*  indices 1 to m mean auxiliary variables and indices m+1 to m+n mean
alpar@1
  1211
*  structural ones.
alpar@1
  1212
*
alpar@1
  1213
*  RETURNS
alpar@1
  1214
*
alpar@1
  1215
*  The routine returns len', which is the number of non-zero elements
alpar@1
  1216
*  in the resultant column stored in the arrays ind and val.
alpar@1
  1217
*
alpar@1
  1218
*  BACKGROUND
alpar@1
  1219
*
alpar@1
  1220
*  The explicitly specified column (1) is transformed in the same way
alpar@1
  1221
*  as any other column of the constraint matrix using the formula:
alpar@1
  1222
*
alpar@1
  1223
*     alfa = inv(B) * a,                                             (3)
alpar@1
  1224
*
alpar@1
  1225
*  where alfa is the resultant column computed by the routine. */
alpar@1
  1226
alpar@1
  1227
int glp_transform_col(glp_prob *P, int len, int ind[], double val[])
alpar@1
  1228
{     int i, m, t;
alpar@1
  1229
      double *a, *alfa;
alpar@1
  1230
      if (!glp_bf_exists(P))
alpar@1
  1231
         xerror("glp_transform_col: basis factorization does not exist "
alpar@1
  1232
            "\n");
alpar@1
  1233
      m = glp_get_num_rows(P);
alpar@1
  1234
      /* unpack the column to be transformed to the array a */
alpar@1
  1235
      a = xcalloc(1+m, sizeof(double));
alpar@1
  1236
      for (i = 1; i <= m; i++) a[i] = 0.0;
alpar@1
  1237
      if (!(0 <= len && len <= m))
alpar@1
  1238
         xerror("glp_transform_col: len = %d; invalid column length\n",
alpar@1
  1239
            len);
alpar@1
  1240
      for (t = 1; t <= len; t++)
alpar@1
  1241
      {  i = ind[t];
alpar@1
  1242
         if (!(1 <= i && i <= m))
alpar@1
  1243
            xerror("glp_transform_col: ind[%d] = %d; row index out of r"
alpar@1
  1244
               "ange\n", t, i);
alpar@1
  1245
         if (val[t] == 0.0)
alpar@1
  1246
            xerror("glp_transform_col: val[%d] = 0; zero coefficient no"
alpar@1
  1247
               "t allowed\n", t);
alpar@1
  1248
         if (a[i] != 0.0)
alpar@1
  1249
            xerror("glp_transform_col: ind[%d] = %d; duplicate row indi"
alpar@1
  1250
               "ces not allowed\n", t, i);
alpar@1
  1251
         a[i] = val[t];
alpar@1
  1252
      }
alpar@1
  1253
      /* solve the system B*a = alfa to compute the vector alfa */
alpar@1
  1254
      alfa = a, glp_ftran(P, alfa);
alpar@1
  1255
      /* store resultant coefficients */
alpar@1
  1256
      len = 0;
alpar@1
  1257
      for (i = 1; i <= m; i++)
alpar@1
  1258
      {  if (alfa[i] != 0.0)
alpar@1
  1259
         {  len++;
alpar@1
  1260
            ind[len] = glp_get_bhead(P, i);
alpar@1
  1261
            val[len] = alfa[i];
alpar@1
  1262
         }
alpar@1
  1263
      }
alpar@1
  1264
      xfree(a);
alpar@1
  1265
      return len;
alpar@1
  1266
}
alpar@1
  1267
alpar@1
  1268
/***********************************************************************
alpar@1
  1269
*  NAME
alpar@1
  1270
*
alpar@1
  1271
*  glp_prim_rtest - perform primal ratio test
alpar@1
  1272
*
alpar@1
  1273
*  SYNOPSIS
alpar@1
  1274
*
alpar@1
  1275
*  int glp_prim_rtest(glp_prob *P, int len, const int ind[],
alpar@1
  1276
*     const double val[], int dir, double eps);
alpar@1
  1277
*
alpar@1
  1278
*  DESCRIPTION
alpar@1
  1279
*
alpar@1
  1280
*  The routine glp_prim_rtest performs the primal ratio test using an
alpar@1
  1281
*  explicitly specified column of the simplex table.
alpar@1
  1282
*
alpar@1
  1283
*  The current basic solution associated with the LP problem object
alpar@1
  1284
*  must be primal feasible.
alpar@1
  1285
*
alpar@1
  1286
*  The explicitly specified column of the simplex table shows how the
alpar@1
  1287
*  basic variables xB depend on some non-basic variable x (which is not
alpar@1
  1288
*  necessarily presented in the problem object):
alpar@1
  1289
*
alpar@1
  1290
*     xB[1] = ... + alfa[1] * x + ...
alpar@1
  1291
*     xB[2] = ... + alfa[2] * x + ...                                (*)
alpar@1
  1292
*         .  .  .  .  .  .  .  .
alpar@1
  1293
*     xB[m] = ... + alfa[m] * x + ...
alpar@1
  1294
*
alpar@1
  1295
*  The column (*) is specifed on entry to the routine using the sparse
alpar@1
  1296
*  format. Ordinal numbers of basic variables xB[i] should be placed in
alpar@1
  1297
*  locations ind[1], ..., ind[len], where ordinal number 1 to m denote
alpar@1
  1298
*  auxiliary variables, and ordinal numbers m+1 to m+n denote structural
alpar@1
  1299
*  variables. The corresponding non-zero coefficients alfa[i] should be
alpar@1
  1300
*  placed in locations val[1], ..., val[len]. The arrays ind and val are
alpar@1
  1301
*  not changed on exit.
alpar@1
  1302
*
alpar@1
  1303
*  The parameter dir specifies direction in which the variable x changes
alpar@1
  1304
*  on entering the basis: +1 means increasing, -1 means decreasing.
alpar@1
  1305
*
alpar@1
  1306
*  The parameter eps is an absolute tolerance (small positive number)
alpar@1
  1307
*  used by the routine to skip small alfa[j] of the row (*).
alpar@1
  1308
*
alpar@1
  1309
*  The routine determines which basic variable (among specified in
alpar@1
  1310
*  ind[1], ..., ind[len]) should leave the basis in order to keep primal
alpar@1
  1311
*  feasibility.
alpar@1
  1312
*
alpar@1
  1313
*  RETURNS
alpar@1
  1314
*
alpar@1
  1315
*  The routine glp_prim_rtest returns the index piv in the arrays ind
alpar@1
  1316
*  and val corresponding to the pivot element chosen, 1 <= piv <= len.
alpar@1
  1317
*  If the adjacent basic solution is primal unbounded and therefore the
alpar@1
  1318
*  choice cannot be made, the routine returns zero.
alpar@1
  1319
*
alpar@1
  1320
*  COMMENTS
alpar@1
  1321
*
alpar@1
  1322
*  If the non-basic variable x is presented in the LP problem object,
alpar@1
  1323
*  the column (*) can be computed with the routine glp_eval_tab_col;
alpar@1
  1324
*  otherwise it can be computed with the routine glp_transform_col. */
alpar@1
  1325
alpar@1
  1326
int glp_prim_rtest(glp_prob *P, int len, const int ind[],
alpar@1
  1327
      const double val[], int dir, double eps)
alpar@1
  1328
{     int k, m, n, piv, t, type, stat;
alpar@1
  1329
      double alfa, big, beta, lb, ub, temp, teta;
alpar@1
  1330
      if (glp_get_prim_stat(P) != GLP_FEAS)
alpar@1
  1331
         xerror("glp_prim_rtest: basic solution is not primal feasible "
alpar@1
  1332
            "\n");
alpar@1
  1333
      if (!(dir == +1 || dir == -1))
alpar@1
  1334
         xerror("glp_prim_rtest: dir = %d; invalid parameter\n", dir);
alpar@1
  1335
      if (!(0.0 < eps && eps < 1.0))
alpar@1
  1336
         xerror("glp_prim_rtest: eps = %g; invalid parameter\n", eps);
alpar@1
  1337
      m = glp_get_num_rows(P);
alpar@1
  1338
      n = glp_get_num_cols(P);
alpar@1
  1339
      /* initial settings */
alpar@1
  1340
      piv = 0, teta = DBL_MAX, big = 0.0;
alpar@1
  1341
      /* walk through the entries of the specified column */
alpar@1
  1342
      for (t = 1; t <= len; t++)
alpar@1
  1343
      {  /* get the ordinal number of basic variable */
alpar@1
  1344
         k = ind[t];
alpar@1
  1345
         if (!(1 <= k && k <= m+n))
alpar@1
  1346
            xerror("glp_prim_rtest: ind[%d] = %d; variable number out o"
alpar@1
  1347
               "f range\n", t, k);
alpar@1
  1348
         /* determine type, bounds, status and primal value of basic
alpar@1
  1349
            variable xB[i] = x[k] in the current basic solution */
alpar@1
  1350
         if (k <= m)
alpar@1
  1351
         {  type = glp_get_row_type(P, k);
alpar@1
  1352
            lb = glp_get_row_lb(P, k);
alpar@1
  1353
            ub = glp_get_row_ub(P, k);
alpar@1
  1354
            stat = glp_get_row_stat(P, k);
alpar@1
  1355
            beta = glp_get_row_prim(P, k);
alpar@1
  1356
         }
alpar@1
  1357
         else
alpar@1
  1358
         {  type = glp_get_col_type(P, k-m);
alpar@1
  1359
            lb = glp_get_col_lb(P, k-m);
alpar@1
  1360
            ub = glp_get_col_ub(P, k-m);
alpar@1
  1361
            stat = glp_get_col_stat(P, k-m);
alpar@1
  1362
            beta = glp_get_col_prim(P, k-m);
alpar@1
  1363
         }
alpar@1
  1364
         if (stat != GLP_BS)
alpar@1
  1365
            xerror("glp_prim_rtest: ind[%d] = %d; non-basic variable no"
alpar@1
  1366
               "t allowed\n", t, k);
alpar@1
  1367
         /* determine influence coefficient at basic variable xB[i]
alpar@1
  1368
            in the explicitly specified column and turn to the case of
alpar@1
  1369
            increasing the variable x in order to simplify the program
alpar@1
  1370
            logic */
alpar@1
  1371
         alfa = (dir > 0 ? + val[t] : - val[t]);
alpar@1
  1372
         /* analyze main cases */
alpar@1
  1373
         if (type == GLP_FR)
alpar@1
  1374
         {  /* xB[i] is free variable */
alpar@1
  1375
            continue;
alpar@1
  1376
         }
alpar@1
  1377
         else if (type == GLP_LO)
alpar@1
  1378
lo:      {  /* xB[i] has an lower bound */
alpar@1
  1379
            if (alfa > - eps) continue;
alpar@1
  1380
            temp = (lb - beta) / alfa;
alpar@1
  1381
         }
alpar@1
  1382
         else if (type == GLP_UP)
alpar@1
  1383
up:      {  /* xB[i] has an upper bound */
alpar@1
  1384
            if (alfa < + eps) continue;
alpar@1
  1385
            temp = (ub - beta) / alfa;
alpar@1
  1386
         }
alpar@1
  1387
         else if (type == GLP_DB)
alpar@1
  1388
         {  /* xB[i] has both lower and upper bounds */
alpar@1
  1389
            if (alfa < 0.0) goto lo; else goto up;
alpar@1
  1390
         }
alpar@1
  1391
         else if (type == GLP_FX)
alpar@1
  1392
         {  /* xB[i] is fixed variable */
alpar@1
  1393
            if (- eps < alfa && alfa < + eps) continue;
alpar@1
  1394
            temp = 0.0;
alpar@1
  1395
         }
alpar@1
  1396
         else
alpar@1
  1397
            xassert(type != type);
alpar@1
  1398
         /* if the value of the variable xB[i] violates its lower or
alpar@1
  1399
            upper bound (slightly, because the current basis is assumed
alpar@1
  1400
            to be primal feasible), temp is negative; we can think this
alpar@1
  1401
            happens due to round-off errors and the value is exactly on
alpar@1
  1402
            the bound; this allows replacing temp by zero */
alpar@1
  1403
         if (temp < 0.0) temp = 0.0;
alpar@1
  1404
         /* apply the minimal ratio test */
alpar@1
  1405
         if (teta > temp || teta == temp && big < fabs(alfa))
alpar@1
  1406
            piv = t, teta = temp, big = fabs(alfa);
alpar@1
  1407
      }
alpar@1
  1408
      /* return index of the pivot element chosen */
alpar@1
  1409
      return piv;
alpar@1
  1410
}
alpar@1
  1411
alpar@1
  1412
/***********************************************************************
alpar@1
  1413
*  NAME
alpar@1
  1414
*
alpar@1
  1415
*  glp_dual_rtest - perform dual ratio test
alpar@1
  1416
*
alpar@1
  1417
*  SYNOPSIS
alpar@1
  1418
*
alpar@1
  1419
*  int glp_dual_rtest(glp_prob *P, int len, const int ind[],
alpar@1
  1420
*     const double val[], int dir, double eps);
alpar@1
  1421
*
alpar@1
  1422
*  DESCRIPTION
alpar@1
  1423
*
alpar@1
  1424
*  The routine glp_dual_rtest performs the dual ratio test using an
alpar@1
  1425
*  explicitly specified row of the simplex table.
alpar@1
  1426
*
alpar@1
  1427
*  The current basic solution associated with the LP problem object
alpar@1
  1428
*  must be dual feasible.
alpar@1
  1429
*
alpar@1
  1430
*  The explicitly specified row of the simplex table is a linear form
alpar@1
  1431
*  that shows how some basic variable x (which is not necessarily
alpar@1
  1432
*  presented in the problem object) depends on non-basic variables xN:
alpar@1
  1433
*
alpar@1
  1434
*     x = alfa[1] * xN[1] + alfa[2] * xN[2] + ... + alfa[n] * xN[n]. (*)
alpar@1
  1435
*
alpar@1
  1436
*  The row (*) is specified on entry to the routine using the sparse
alpar@1
  1437
*  format. Ordinal numbers of non-basic variables xN[j] should be placed
alpar@1
  1438
*  in locations ind[1], ..., ind[len], where ordinal numbers 1 to m
alpar@1
  1439
*  denote auxiliary variables, and ordinal numbers m+1 to m+n denote
alpar@1
  1440
*  structural variables. The corresponding non-zero coefficients alfa[j]
alpar@1
  1441
*  should be placed in locations val[1], ..., val[len]. The arrays ind
alpar@1
  1442
*  and val are not changed on exit.
alpar@1
  1443
*
alpar@1
  1444
*  The parameter dir specifies direction in which the variable x changes
alpar@1
  1445
*  on leaving the basis: +1 means that x goes to its lower bound, and -1
alpar@1
  1446
*  means that x goes to its upper bound.
alpar@1
  1447
*
alpar@1
  1448
*  The parameter eps is an absolute tolerance (small positive number)
alpar@1
  1449
*  used by the routine to skip small alfa[j] of the row (*).
alpar@1
  1450
*
alpar@1
  1451
*  The routine determines which non-basic variable (among specified in
alpar@1
  1452
*  ind[1], ..., ind[len]) should enter the basis in order to keep dual
alpar@1
  1453
*  feasibility.
alpar@1
  1454
*
alpar@1
  1455
*  RETURNS
alpar@1
  1456
*
alpar@1
  1457
*  The routine glp_dual_rtest returns the index piv in the arrays ind
alpar@1
  1458
*  and val corresponding to the pivot element chosen, 1 <= piv <= len.
alpar@1
  1459
*  If the adjacent basic solution is dual unbounded and therefore the
alpar@1
  1460
*  choice cannot be made, the routine returns zero.
alpar@1
  1461
*
alpar@1
  1462
*  COMMENTS
alpar@1
  1463
*
alpar@1
  1464
*  If the basic variable x is presented in the LP problem object, the
alpar@1
  1465
*  row (*) can be computed with the routine glp_eval_tab_row; otherwise
alpar@1
  1466
*  it can be computed with the routine glp_transform_row. */
alpar@1
  1467
alpar@1
  1468
int glp_dual_rtest(glp_prob *P, int len, const int ind[],
alpar@1
  1469
      const double val[], int dir, double eps)
alpar@1
  1470
{     int k, m, n, piv, t, stat;
alpar@1
  1471
      double alfa, big, cost, obj, temp, teta;
alpar@1
  1472
      if (glp_get_dual_stat(P) != GLP_FEAS)
alpar@1
  1473
         xerror("glp_dual_rtest: basic solution is not dual feasible\n")
alpar@1
  1474
            ;
alpar@1
  1475
      if (!(dir == +1 || dir == -1))
alpar@1
  1476
         xerror("glp_dual_rtest: dir = %d; invalid parameter\n", dir);
alpar@1
  1477
      if (!(0.0 < eps && eps < 1.0))
alpar@1
  1478
         xerror("glp_dual_rtest: eps = %g; invalid parameter\n", eps);
alpar@1
  1479
      m = glp_get_num_rows(P);
alpar@1
  1480
      n = glp_get_num_cols(P);
alpar@1
  1481
      /* take into account optimization direction */
alpar@1
  1482
      obj = (glp_get_obj_dir(P) == GLP_MIN ? +1.0 : -1.0);
alpar@1
  1483
      /* initial settings */
alpar@1
  1484
      piv = 0, teta = DBL_MAX, big = 0.0;
alpar@1
  1485
      /* walk through the entries of the specified row */
alpar@1
  1486
      for (t = 1; t <= len; t++)
alpar@1
  1487
      {  /* get ordinal number of non-basic variable */
alpar@1
  1488
         k = ind[t];
alpar@1
  1489
         if (!(1 <= k && k <= m+n))
alpar@1
  1490
            xerror("glp_dual_rtest: ind[%d] = %d; variable number out o"
alpar@1
  1491
               "f range\n", t, k);
alpar@1
  1492
         /* determine status and reduced cost of non-basic variable
alpar@1
  1493
            x[k] = xN[j] in the current basic solution */
alpar@1
  1494
         if (k <= m)
alpar@1
  1495
         {  stat = glp_get_row_stat(P, k);
alpar@1
  1496
            cost = glp_get_row_dual(P, k);
alpar@1
  1497
         }
alpar@1
  1498
         else
alpar@1
  1499
         {  stat = glp_get_col_stat(P, k-m);
alpar@1
  1500
            cost = glp_get_col_dual(P, k-m);
alpar@1
  1501
         }
alpar@1
  1502
         if (stat == GLP_BS)
alpar@1
  1503
            xerror("glp_dual_rtest: ind[%d] = %d; basic variable not al"
alpar@1
  1504
               "lowed\n", t, k);
alpar@1
  1505
         /* determine influence coefficient at non-basic variable xN[j]
alpar@1
  1506
            in the explicitly specified row and turn to the case of
alpar@1
  1507
            increasing the variable x in order to simplify the program
alpar@1
  1508
            logic */
alpar@1
  1509
         alfa = (dir > 0 ? + val[t] : - val[t]);
alpar@1
  1510
         /* analyze main cases */
alpar@1
  1511
         if (stat == GLP_NL)
alpar@1
  1512
         {  /* xN[j] is on its lower bound */
alpar@1
  1513
            if (alfa < + eps) continue;
alpar@1
  1514
            temp = (obj * cost) / alfa;
alpar@1
  1515
         }
alpar@1
  1516
         else if (stat == GLP_NU)
alpar@1
  1517
         {  /* xN[j] is on its upper bound */
alpar@1
  1518
            if (alfa > - eps) continue;
alpar@1
  1519
            temp = (obj * cost) / alfa;
alpar@1
  1520
         }
alpar@1
  1521
         else if (stat == GLP_NF)
alpar@1
  1522
         {  /* xN[j] is non-basic free variable */
alpar@1
  1523
            if (- eps < alfa && alfa < + eps) continue;
alpar@1
  1524
            temp = 0.0;
alpar@1
  1525
         }
alpar@1
  1526
         else if (stat == GLP_NS)
alpar@1
  1527
         {  /* xN[j] is non-basic fixed variable */
alpar@1
  1528
            continue;
alpar@1
  1529
         }
alpar@1
  1530
         else
alpar@1
  1531
            xassert(stat != stat);
alpar@1
  1532
         /* if the reduced cost of the variable xN[j] violates its zero
alpar@1
  1533
            bound (slightly, because the current basis is assumed to be
alpar@1
  1534
            dual feasible), temp is negative; we can think this happens
alpar@1
  1535
            due to round-off errors and the reduced cost is exact zero;
alpar@1
  1536
            this allows replacing temp by zero */
alpar@1
  1537
         if (temp < 0.0) temp = 0.0;
alpar@1
  1538
         /* apply the minimal ratio test */
alpar@1
  1539
         if (teta > temp || teta == temp && big < fabs(alfa))
alpar@1
  1540
            piv = t, teta = temp, big = fabs(alfa);
alpar@1
  1541
      }
alpar@1
  1542
      /* return index of the pivot element chosen */
alpar@1
  1543
      return piv;
alpar@1
  1544
}
alpar@1
  1545
alpar@1
  1546
/***********************************************************************
alpar@1
  1547
*  NAME
alpar@1
  1548
*
alpar@1
  1549
*  glp_analyze_row - simulate one iteration of dual simplex method
alpar@1
  1550
*
alpar@1
  1551
*  SYNOPSIS
alpar@1
  1552
*
alpar@1
  1553
*  int glp_analyze_row(glp_prob *P, int len, const int ind[],
alpar@1
  1554
*     const double val[], int type, double rhs, double eps, int *piv,
alpar@1
  1555
*     double *x, double *dx, double *y, double *dy, double *dz);
alpar@1
  1556
*
alpar@1
  1557
*  DESCRIPTION
alpar@1
  1558
*
alpar@1
  1559
*  Let the current basis be optimal or dual feasible, and there be
alpar@1
  1560
*  specified a row (constraint), which is violated by the current basic
alpar@1
  1561
*  solution. The routine glp_analyze_row simulates one iteration of the
alpar@1
  1562
*  dual simplex method to determine some information on the adjacent
alpar@1
  1563
*  basis (see below), where the specified row becomes active constraint
alpar@1
  1564
*  (i.e. its auxiliary variable becomes non-basic).
alpar@1
  1565
*
alpar@1
  1566
*  The current basic solution associated with the problem object passed
alpar@1
  1567
*  to the routine must be dual feasible, and its primal components must
alpar@1
  1568
*  be defined.
alpar@1
  1569
*
alpar@1
  1570
*  The row to be analyzed must be previously transformed either with
alpar@1
  1571
*  the routine glp_eval_tab_row (if the row is in the problem object)
alpar@1
  1572
*  or with the routine glp_transform_row (if the row is external, i.e.
alpar@1
  1573
*  not in the problem object). This is needed to express the row only
alpar@1
  1574
*  through (auxiliary and structural) variables, which are non-basic in
alpar@1
  1575
*  the current basis:
alpar@1
  1576
*
alpar@1
  1577
*     y = alfa[1] * xN[1] + alfa[2] * xN[2] + ... + alfa[n] * xN[n],
alpar@1
  1578
*
alpar@1
  1579
*  where y is an auxiliary variable of the row, alfa[j] is an influence
alpar@1
  1580
*  coefficient, xN[j] is a non-basic variable.
alpar@1
  1581
*
alpar@1
  1582
*  The row is passed to the routine in sparse format. Ordinal numbers
alpar@1
  1583
*  of non-basic variables are stored in locations ind[1], ..., ind[len],
alpar@1
  1584
*  where numbers 1 to m denote auxiliary variables while numbers m+1 to
alpar@1
  1585
*  m+n denote structural variables. Corresponding non-zero coefficients
alpar@1
  1586
*  alfa[j] are stored in locations val[1], ..., val[len]. The arrays
alpar@1
  1587
*  ind and val are ot changed on exit.
alpar@1
  1588
*
alpar@1
  1589
*  The parameters type and rhs specify the row type and its right-hand
alpar@1
  1590
*  side as follows:
alpar@1
  1591
*
alpar@1
  1592
*     type = GLP_LO: y = sum alfa[j] * xN[j] >= rhs
alpar@1
  1593
*
alpar@1
  1594
*     type = GLP_UP: y = sum alfa[j] * xN[j] <= rhs
alpar@1
  1595
*
alpar@1
  1596
*  The parameter eps is an absolute tolerance (small positive number)
alpar@1
  1597
*  used by the routine to skip small coefficients alfa[j] on performing
alpar@1
  1598
*  the dual ratio test.
alpar@1
  1599
*
alpar@1
  1600
*  If the operation was successful, the routine stores the following
alpar@1
  1601
*  information to corresponding location (if some parameter is NULL,
alpar@1
  1602
*  its value is not stored):
alpar@1
  1603
*
alpar@1
  1604
*  piv   index in the array ind and val, 1 <= piv <= len, determining
alpar@1
  1605
*        the non-basic variable, which would enter the adjacent basis;
alpar@1
  1606
*
alpar@1
  1607
*  x     value of the non-basic variable in the current basis;
alpar@1
  1608
*
alpar@1
  1609
*  dx    difference between values of the non-basic variable in the
alpar@1
  1610
*        adjacent and current bases, dx = x.new - x.old;
alpar@1
  1611
*
alpar@1
  1612
*  y     value of the row (i.e. of its auxiliary variable) in the
alpar@1
  1613
*        current basis;
alpar@1
  1614
*
alpar@1
  1615
*  dy    difference between values of the row in the adjacent and
alpar@1
  1616
*        current bases, dy = y.new - y.old;
alpar@1
  1617
*
alpar@1
  1618
*  dz    difference between values of the objective function in the
alpar@1
  1619
*        adjacent and current bases, dz = z.new - z.old. Note that in
alpar@1
  1620
*        case of minimization dz >= 0, and in case of maximization
alpar@1
  1621
*        dz <= 0, i.e. in the adjacent basis the objective function
alpar@1
  1622
*        always gets worse (degrades). */
alpar@1
  1623
alpar@1
  1624
int _glp_analyze_row(glp_prob *P, int len, const int ind[],
alpar@1
  1625
      const double val[], int type, double rhs, double eps, int *_piv,
alpar@1
  1626
      double *_x, double *_dx, double *_y, double *_dy, double *_dz)
alpar@1
  1627
{     int t, k, dir, piv, ret = 0;
alpar@1
  1628
      double x, dx, y, dy, dz;
alpar@1
  1629
      if (P->pbs_stat == GLP_UNDEF)
alpar@1
  1630
         xerror("glp_analyze_row: primal basic solution components are "
alpar@1
  1631
            "undefined\n");
alpar@1
  1632
      if (P->dbs_stat != GLP_FEAS)
alpar@1
  1633
         xerror("glp_analyze_row: basic solution is not dual feasible\n"
alpar@1
  1634
            );
alpar@1
  1635
      /* compute the row value y = sum alfa[j] * xN[j] in the current
alpar@1
  1636
         basis */
alpar@1
  1637
      if (!(0 <= len && len <= P->n))
alpar@1
  1638
         xerror("glp_analyze_row: len = %d; invalid row length\n", len);
alpar@1
  1639
      y = 0.0;
alpar@1
  1640
      for (t = 1; t <= len; t++)
alpar@1
  1641
      {  /* determine value of x[k] = xN[j] in the current basis */
alpar@1
  1642
         k = ind[t];
alpar@1
  1643
         if (!(1 <= k && k <= P->m+P->n))
alpar@1
  1644
            xerror("glp_analyze_row: ind[%d] = %d; row/column index out"
alpar@1
  1645
               " of range\n", t, k);
alpar@1
  1646
         if (k <= P->m)
alpar@1
  1647
         {  /* x[k] is auxiliary variable */
alpar@1
  1648
            if (P->row[k]->stat == GLP_BS)
alpar@1
  1649
               xerror("glp_analyze_row: ind[%d] = %d; basic auxiliary v"
alpar@1
  1650
                  "ariable is not allowed\n", t, k);
alpar@1
  1651
            x = P->row[k]->prim;
alpar@1
  1652
         }
alpar@1
  1653
         else
alpar@1
  1654
         {  /* x[k] is structural variable */
alpar@1
  1655
            if (P->col[k-P->m]->stat == GLP_BS)
alpar@1
  1656
               xerror("glp_analyze_row: ind[%d] = %d; basic structural "
alpar@1
  1657
                  "variable is not allowed\n", t, k);
alpar@1
  1658
            x = P->col[k-P->m]->prim;
alpar@1
  1659
         }
alpar@1
  1660
         y += val[t] * x;
alpar@1
  1661
      }
alpar@1
  1662
      /* check if the row is primal infeasible in the current basis,
alpar@1
  1663
         i.e. the constraint is violated at the current point */
alpar@1
  1664
      if (type == GLP_LO)
alpar@1
  1665
      {  if (y >= rhs)
alpar@1
  1666
         {  /* the constraint is not violated */
alpar@1
  1667
            ret = 1;
alpar@1
  1668
            goto done;
alpar@1
  1669
         }
alpar@1
  1670
         /* in the adjacent basis y goes to its lower bound */
alpar@1
  1671
         dir = +1;
alpar@1
  1672
      }
alpar@1
  1673
      else if (type == GLP_UP)
alpar@1
  1674
      {  if (y <= rhs)
alpar@1
  1675
         {  /* the constraint is not violated */
alpar@1
  1676
            ret = 1;
alpar@1
  1677
            goto done;
alpar@1
  1678
         }
alpar@1
  1679
         /* in the adjacent basis y goes to its upper bound */
alpar@1
  1680
         dir = -1;
alpar@1
  1681
      }
alpar@1
  1682
      else
alpar@1
  1683
         xerror("glp_analyze_row: type = %d; invalid parameter\n",
alpar@1
  1684
            type);
alpar@1
  1685
      /* compute dy = y.new - y.old */
alpar@1
  1686
      dy = rhs - y;
alpar@1
  1687
      /* perform dual ratio test to determine which non-basic variable
alpar@1
  1688
         should enter the adjacent basis to keep it dual feasible */
alpar@1
  1689
      piv = glp_dual_rtest(P, len, ind, val, dir, eps);
alpar@1
  1690
      if (piv == 0)
alpar@1
  1691
      {  /* no dual feasible adjacent basis exists */
alpar@1
  1692
         ret = 2;
alpar@1
  1693
         goto done;
alpar@1
  1694
      }
alpar@1
  1695
      /* non-basic variable x[k] = xN[j] should enter the basis */
alpar@1
  1696
      k = ind[piv];
alpar@1
  1697
      xassert(1 <= k && k <= P->m+P->n);
alpar@1
  1698
      /* determine its value in the current basis */
alpar@1
  1699
      if (k <= P->m)
alpar@1
  1700
         x = P->row[k]->prim;
alpar@1
  1701
      else
alpar@1
  1702
         x = P->col[k-P->m]->prim;
alpar@1
  1703
      /* compute dx = x.new - x.old = dy / alfa[j] */
alpar@1
  1704
      xassert(val[piv] != 0.0);
alpar@1
  1705
      dx = dy / val[piv];
alpar@1
  1706
      /* compute dz = z.new - z.old = d[j] * dx, where d[j] is reduced
alpar@1
  1707
         cost of xN[j] in the current basis */
alpar@1
  1708
      if (k <= P->m)
alpar@1
  1709
         dz = P->row[k]->dual * dx;
alpar@1
  1710
      else
alpar@1
  1711
         dz = P->col[k-P->m]->dual * dx;
alpar@1
  1712
      /* store the analysis results */
alpar@1
  1713
      if (_piv != NULL) *_piv = piv;
alpar@1
  1714
      if (_x   != NULL) *_x   = x;
alpar@1
  1715
      if (_dx  != NULL) *_dx  = dx;
alpar@1
  1716
      if (_y   != NULL) *_y   = y;
alpar@1
  1717
      if (_dy  != NULL) *_dy  = dy;
alpar@1
  1718
      if (_dz  != NULL) *_dz  = dz;
alpar@1
  1719
done: return ret;
alpar@1
  1720
}
alpar@1
  1721
alpar@1
  1722
#if 0
alpar@1
  1723
int main(void)
alpar@1
  1724
{     /* example program for the routine glp_analyze_row */
alpar@1
  1725
      glp_prob *P;
alpar@1
  1726
      glp_smcp parm;
alpar@1
  1727
      int i, k, len, piv, ret, ind[1+100];
alpar@1
  1728
      double rhs, x, dx, y, dy, dz, val[1+100];
alpar@1
  1729
      P = glp_create_prob();
alpar@1
  1730
      /* read plan.mps (see glpk/examples) */
alpar@1
  1731
      ret = glp_read_mps(P, GLP_MPS_DECK, NULL, "plan.mps");
alpar@1
  1732
      glp_assert(ret == 0);
alpar@1
  1733
      /* and solve it to optimality */
alpar@1
  1734
      ret = glp_simplex(P, NULL);
alpar@1
  1735
      glp_assert(ret == 0);
alpar@1
  1736
      glp_assert(glp_get_status(P) == GLP_OPT);
alpar@1
  1737
      /* the optimal objective value is 296.217 */
alpar@1
  1738
      /* we would like to know what happens if we would add a new row
alpar@1
  1739
         (constraint) to plan.mps:
alpar@1
  1740
         .01 * bin1 + .01 * bin2 + .02 * bin4 + .02 * bin5 <= 12 */
alpar@1
  1741
      /* first, we specify this new row */
alpar@1
  1742
      glp_create_index(P);
alpar@1
  1743
      len = 0;
alpar@1
  1744
      ind[++len] = glp_find_col(P, "BIN1"), val[len] = .01;
alpar@1
  1745
      ind[++len] = glp_find_col(P, "BIN2"), val[len] = .01;
alpar@1
  1746
      ind[++len] = glp_find_col(P, "BIN4"), val[len] = .02;
alpar@1
  1747
      ind[++len] = glp_find_col(P, "BIN5"), val[len] = .02;
alpar@1
  1748
      rhs = 12;
alpar@1
  1749
      /* then we can compute value of the row (i.e. of its auxiliary
alpar@1
  1750
         variable) in the current basis to see if the constraint is
alpar@1
  1751
         violated */
alpar@1
  1752
      y = 0.0;
alpar@1
  1753
      for (k = 1; k <= len; k++)
alpar@1
  1754
         y += val[k] * glp_get_col_prim(P, ind[k]);
alpar@1
  1755
      glp_printf("y = %g\n", y);
alpar@1
  1756
      /* this prints y = 15.1372, so the constraint is violated, since
alpar@1
  1757
         we require that y <= rhs = 12 */
alpar@1
  1758
      /* now we transform the row to express it only through non-basic
alpar@1
  1759
         (auxiliary and artificial) variables */
alpar@1
  1760
      len = glp_transform_row(P, len, ind, val);
alpar@1
  1761
      /* finally, we simulate one step of the dual simplex method to
alpar@1
  1762
         obtain necessary information for the adjacent basis */
alpar@1
  1763
      ret = _glp_analyze_row(P, len, ind, val, GLP_UP, rhs, 1e-9, &piv,
alpar@1
  1764
         &x, &dx, &y, &dy, &dz);
alpar@1
  1765
      glp_assert(ret == 0);
alpar@1
  1766
      glp_printf("k = %d, x = %g; dx = %g; y = %g; dy = %g; dz = %g\n",
alpar@1
  1767
         ind[piv], x, dx, y, dy, dz);
alpar@1
  1768
      /* this prints dz = 5.64418 and means that in the adjacent basis
alpar@1
  1769
         the objective function would be 296.217 + 5.64418 = 301.861 */
alpar@1
  1770
      /* now we actually include the row into the problem object; note
alpar@1
  1771
         that the arrays ind and val are clobbered, so we need to build
alpar@1
  1772
         them once again */
alpar@1
  1773
      len = 0;
alpar@1
  1774
      ind[++len] = glp_find_col(P, "BIN1"), val[len] = .01;
alpar@1
  1775
      ind[++len] = glp_find_col(P, "BIN2"), val[len] = .01;
alpar@1
  1776
      ind[++len] = glp_find_col(P, "BIN4"), val[len] = .02;
alpar@1
  1777
      ind[++len] = glp_find_col(P, "BIN5"), val[len] = .02;
alpar@1
  1778
      rhs = 12;
alpar@1
  1779
      i = glp_add_rows(P, 1);
alpar@1
  1780
      glp_set_row_bnds(P, i, GLP_UP, 0, rhs);
alpar@1
  1781
      glp_set_mat_row(P, i, len, ind, val);
alpar@1
  1782
      /* and perform one dual simplex iteration */
alpar@1
  1783
      glp_init_smcp(&parm);
alpar@1
  1784
      parm.meth = GLP_DUAL;
alpar@1
  1785
      parm.it_lim = 1;
alpar@1
  1786
      glp_simplex(P, &parm);
alpar@1
  1787
      /* the current objective value is 301.861 */
alpar@1
  1788
      return 0;
alpar@1
  1789
}
alpar@1
  1790
#endif
alpar@1
  1791
alpar@1
  1792
/***********************************************************************
alpar@1
  1793
*  NAME
alpar@1
  1794
*
alpar@1
  1795
*  glp_analyze_bound - analyze active bound of non-basic variable
alpar@1
  1796
*
alpar@1
  1797
*  SYNOPSIS
alpar@1
  1798
*
alpar@1
  1799
*  void glp_analyze_bound(glp_prob *P, int k, double *limit1, int *var1,
alpar@1
  1800
*     double *limit2, int *var2);
alpar@1
  1801
*
alpar@1
  1802
*  DESCRIPTION
alpar@1
  1803
*
alpar@1
  1804
*  The routine glp_analyze_bound analyzes the effect of varying the
alpar@1
  1805
*  active bound of specified non-basic variable.
alpar@1
  1806
*
alpar@1
  1807
*  The non-basic variable is specified by the parameter k, where
alpar@1
  1808
*  1 <= k <= m means auxiliary variable of corresponding row while
alpar@1
  1809
*  m+1 <= k <= m+n means structural variable (column).
alpar@1
  1810
*
alpar@1
  1811
*  Note that the current basic solution must be optimal, and the basis
alpar@1
  1812
*  factorization must exist.
alpar@1
  1813
*
alpar@1
  1814
*  Results of the analysis have the following meaning.
alpar@1
  1815
*
alpar@1
  1816
*  value1 is the minimal value of the active bound, at which the basis
alpar@1
  1817
*  still remains primal feasible and thus optimal. -DBL_MAX means that
alpar@1
  1818
*  the active bound has no lower limit.
alpar@1
  1819
*
alpar@1
  1820
*  var1 is the ordinal number of an auxiliary (1 to m) or structural
alpar@1
  1821
*  (m+1 to n) basic variable, which reaches its bound first and thereby
alpar@1
  1822
*  limits further decreasing the active bound being analyzed.
alpar@1
  1823
*  if value1 = -DBL_MAX, var1 is set to 0.
alpar@1
  1824
*
alpar@1
  1825
*  value2 is the maximal value of the active bound, at which the basis
alpar@1
  1826
*  still remains primal feasible and thus optimal. +DBL_MAX means that
alpar@1
  1827
*  the active bound has no upper limit.
alpar@1
  1828
*
alpar@1
  1829
*  var2 is the ordinal number of an auxiliary (1 to m) or structural
alpar@1
  1830
*  (m+1 to n) basic variable, which reaches its bound first and thereby
alpar@1
  1831
*  limits further increasing the active bound being analyzed.
alpar@1
  1832
*  if value2 = +DBL_MAX, var2 is set to 0. */
alpar@1
  1833
alpar@1
  1834
void glp_analyze_bound(glp_prob *P, int k, double *value1, int *var1,
alpar@1
  1835
      double *value2, int *var2)
alpar@1
  1836
{     GLPROW *row;
alpar@1
  1837
      GLPCOL *col;
alpar@1
  1838
      int m, n, stat, kase, p, len, piv, *ind;
alpar@1
  1839
      double x, new_x, ll, uu, xx, delta, *val;
alpar@1
  1840
      /* sanity checks */
alpar@1
  1841
      if (P == NULL || P->magic != GLP_PROB_MAGIC)
alpar@1
  1842
         xerror("glp_analyze_bound: P = %p; invalid problem object\n",
alpar@1
  1843
            P);
alpar@1
  1844
      m = P->m, n = P->n;
alpar@1
  1845
      if (!(P->pbs_stat == GLP_FEAS && P->dbs_stat == GLP_FEAS))
alpar@1
  1846
         xerror("glp_analyze_bound: optimal basic solution required\n");
alpar@1
  1847
      if (!(m == 0 || P->valid))
alpar@1
  1848
         xerror("glp_analyze_bound: basis factorization required\n");
alpar@1
  1849
      if (!(1 <= k && k <= m+n))
alpar@1
  1850
         xerror("glp_analyze_bound: k = %d; variable number out of rang"
alpar@1
  1851
            "e\n", k);
alpar@1
  1852
      /* retrieve information about the specified non-basic variable
alpar@1
  1853
         x[k] whose active bound is to be analyzed */
alpar@1
  1854
      if (k <= m)
alpar@1
  1855
      {  row = P->row[k];
alpar@1
  1856
         stat = row->stat;
alpar@1
  1857
         x = row->prim;
alpar@1
  1858
      }
alpar@1
  1859
      else
alpar@1
  1860
      {  col = P->col[k-m];
alpar@1
  1861
         stat = col->stat;
alpar@1
  1862
         x = col->prim;
alpar@1
  1863
      }
alpar@1
  1864
      if (stat == GLP_BS)
alpar@1
  1865
         xerror("glp_analyze_bound: k = %d; basic variable not allowed "
alpar@1
  1866
            "\n", k);
alpar@1
  1867
      /* allocate working arrays */
alpar@1
  1868
      ind = xcalloc(1+m, sizeof(int));
alpar@1
  1869
      val = xcalloc(1+m, sizeof(double));
alpar@1
  1870
      /* compute column of the simplex table corresponding to the
alpar@1
  1871
         non-basic variable x[k] */
alpar@1
  1872
      len = glp_eval_tab_col(P, k, ind, val);
alpar@1
  1873
      xassert(0 <= len && len <= m);
alpar@1
  1874
      /* perform analysis */
alpar@1
  1875
      for (kase = -1; kase <= +1; kase += 2)
alpar@1
  1876
      {  /* kase < 0 means active bound of x[k] is decreasing;
alpar@1
  1877
            kase > 0 means active bound of x[k] is increasing */
alpar@1
  1878
         /* use the primal ratio test to determine some basic variable
alpar@1
  1879
            x[p] which reaches its bound first */
alpar@1
  1880
         piv = glp_prim_rtest(P, len, ind, val, kase, 1e-9);
alpar@1
  1881
         if (piv == 0)
alpar@1
  1882
         {  /* nothing limits changing the active bound of x[k] */
alpar@1
  1883
            p = 0;
alpar@1
  1884
            new_x = (kase < 0 ? -DBL_MAX : +DBL_MAX);
alpar@1
  1885
            goto store;
alpar@1
  1886
         }
alpar@1
  1887
         /* basic variable x[p] limits changing the active bound of
alpar@1
  1888
            x[k]; determine its value in the current basis */
alpar@1
  1889
         xassert(1 <= piv && piv <= len);
alpar@1
  1890
         p = ind[piv];
alpar@1
  1891
         if (p <= m)
alpar@1
  1892
         {  row = P->row[p];
alpar@1
  1893
            ll = glp_get_row_lb(P, row->i);
alpar@1
  1894
            uu = glp_get_row_ub(P, row->i);
alpar@1
  1895
            stat = row->stat;
alpar@1
  1896
            xx = row->prim;
alpar@1
  1897
         }
alpar@1
  1898
         else
alpar@1
  1899
         {  col = P->col[p-m];
alpar@1
  1900
            ll = glp_get_col_lb(P, col->j);
alpar@1
  1901
            uu = glp_get_col_ub(P, col->j);
alpar@1
  1902
            stat = col->stat;
alpar@1
  1903
            xx = col->prim;
alpar@1
  1904
         }
alpar@1
  1905
         xassert(stat == GLP_BS);
alpar@1
  1906
         /* determine delta x[p] = bound of x[p] - value of x[p] */
alpar@1
  1907
         if (kase < 0 && val[piv] > 0.0 ||
alpar@1
  1908
             kase > 0 && val[piv] < 0.0)
alpar@1
  1909
         {  /* delta x[p] < 0, so x[p] goes toward its lower bound */
alpar@1
  1910
            xassert(ll != -DBL_MAX);
alpar@1
  1911
            delta = ll - xx;
alpar@1
  1912
         }
alpar@1
  1913
         else
alpar@1
  1914
         {  /* delta x[p] > 0, so x[p] goes toward its upper bound */
alpar@1
  1915
            xassert(uu != +DBL_MAX);
alpar@1
  1916
            delta = uu - xx;
alpar@1
  1917
         }
alpar@1
  1918
         /* delta x[p] = alfa[p,k] * delta x[k], so new x[k] = x[k] +
alpar@1
  1919
            delta x[k] = x[k] + delta x[p] / alfa[p,k] is the value of
alpar@1
  1920
            x[k] in the adjacent basis */
alpar@1
  1921
         xassert(val[piv] != 0.0);
alpar@1
  1922
         new_x = x + delta / val[piv];
alpar@1
  1923
store:   /* store analysis results */
alpar@1
  1924
         if (kase < 0)
alpar@1
  1925
         {  if (value1 != NULL) *value1 = new_x;
alpar@1
  1926
            if (var1 != NULL) *var1 = p;
alpar@1
  1927
         }
alpar@1
  1928
         else
alpar@1
  1929
         {  if (value2 != NULL) *value2 = new_x;
alpar@1
  1930
            if (var2 != NULL) *var2 = p;
alpar@1
  1931
         }
alpar@1
  1932
      }
alpar@1
  1933
      /* free working arrays */
alpar@1
  1934
      xfree(ind);
alpar@1
  1935
      xfree(val);
alpar@1
  1936
      return;
alpar@1
  1937
}
alpar@1
  1938
alpar@1
  1939
/***********************************************************************
alpar@1
  1940
*  NAME
alpar@1
  1941
*
alpar@1
  1942
*  glp_analyze_coef - analyze objective coefficient at basic variable
alpar@1
  1943
*
alpar@1
  1944
*  SYNOPSIS
alpar@1
  1945
*
alpar@1
  1946
*  void glp_analyze_coef(glp_prob *P, int k, double *coef1, int *var1,
alpar@1
  1947
*     double *value1, double *coef2, int *var2, double *value2);
alpar@1
  1948
*
alpar@1
  1949
*  DESCRIPTION
alpar@1
  1950
*
alpar@1
  1951
*  The routine glp_analyze_coef analyzes the effect of varying the
alpar@1
  1952
*  objective coefficient at specified basic variable.
alpar@1
  1953
*
alpar@1
  1954
*  The basic variable is specified by the parameter k, where
alpar@1
  1955
*  1 <= k <= m means auxiliary variable of corresponding row while
alpar@1
  1956
*  m+1 <= k <= m+n means structural variable (column).
alpar@1
  1957
*
alpar@1
  1958
*  Note that the current basic solution must be optimal, and the basis
alpar@1
  1959
*  factorization must exist.
alpar@1
  1960
*
alpar@1
  1961
*  Results of the analysis have the following meaning.
alpar@1
  1962
*
alpar@1
  1963
*  coef1 is the minimal value of the objective coefficient, at which
alpar@1
  1964
*  the basis still remains dual feasible and thus optimal. -DBL_MAX
alpar@1
  1965
*  means that the objective coefficient has no lower limit.
alpar@1
  1966
*
alpar@1
  1967
*  var1 is the ordinal number of an auxiliary (1 to m) or structural
alpar@1
  1968
*  (m+1 to n) non-basic variable, whose reduced cost reaches its zero
alpar@1
  1969
*  bound first and thereby limits further decreasing the objective
alpar@1
  1970
*  coefficient being analyzed. If coef1 = -DBL_MAX, var1 is set to 0.
alpar@1
  1971
*
alpar@1
  1972
*  value1 is value of the basic variable being analyzed in an adjacent
alpar@1
  1973
*  basis, which is defined as follows. Let the objective coefficient
alpar@1
  1974
*  reaches its minimal value (coef1) and continues decreasing. Then the
alpar@1
  1975
*  reduced cost of the limiting non-basic variable (var1) becomes dual
alpar@1
  1976
*  infeasible and the current basis becomes non-optimal that forces the
alpar@1
  1977
*  limiting non-basic variable to enter the basis replacing there some
alpar@1
  1978
*  basic variable that leaves the basis to keep primal feasibility.
alpar@1
  1979
*  Should note that on determining the adjacent basis current bounds
alpar@1
  1980
*  of the basic variable being analyzed are ignored as if it were free
alpar@1
  1981
*  (unbounded) variable, so it cannot leave the basis. It may happen
alpar@1
  1982
*  that no dual feasible adjacent basis exists, in which case value1 is
alpar@1
  1983
*  set to -DBL_MAX or +DBL_MAX.
alpar@1
  1984
*
alpar@1
  1985
*  coef2 is the maximal value of the objective coefficient, at which
alpar@1
  1986
*  the basis still remains dual feasible and thus optimal. +DBL_MAX
alpar@1
  1987
*  means that the objective coefficient has no upper limit.
alpar@1
  1988
*
alpar@1
  1989
*  var2 is the ordinal number of an auxiliary (1 to m) or structural
alpar@1
  1990
*  (m+1 to n) non-basic variable, whose reduced cost reaches its zero
alpar@1
  1991
*  bound first and thereby limits further increasing the objective
alpar@1
  1992
*  coefficient being analyzed. If coef2 = +DBL_MAX, var2 is set to 0.
alpar@1
  1993
*
alpar@1
  1994
*  value2 is value of the basic variable being analyzed in an adjacent
alpar@1
  1995
*  basis, which is defined exactly in the same way as value1 above with
alpar@1
  1996
*  exception that now the objective coefficient is increasing. */
alpar@1
  1997
alpar@1
  1998
void glp_analyze_coef(glp_prob *P, int k, double *coef1, int *var1,
alpar@1
  1999
      double *value1, double *coef2, int *var2, double *value2)
alpar@1
  2000
{     GLPROW *row; GLPCOL *col;
alpar@1
  2001
      int m, n, type, stat, kase, p, q, dir, clen, cpiv, rlen, rpiv,
alpar@1
  2002
         *cind, *rind;
alpar@1
  2003
      double lb, ub, coef, x, lim_coef, new_x, d, delta, ll, uu, xx,
alpar@1
  2004
         *rval, *cval;
alpar@1
  2005
      /* sanity checks */
alpar@1
  2006
      if (P == NULL || P->magic != GLP_PROB_MAGIC)
alpar@1
  2007
         xerror("glp_analyze_coef: P = %p; invalid problem object\n",
alpar@1
  2008
            P);
alpar@1
  2009
      m = P->m, n = P->n;
alpar@1
  2010
      if (!(P->pbs_stat == GLP_FEAS && P->dbs_stat == GLP_FEAS))
alpar@1
  2011
         xerror("glp_analyze_coef: optimal basic solution required\n");
alpar@1
  2012
      if (!(m == 0 || P->valid))
alpar@1
  2013
         xerror("glp_analyze_coef: basis factorization required\n");
alpar@1
  2014
      if (!(1 <= k && k <= m+n))
alpar@1
  2015
         xerror("glp_analyze_coef: k = %d; variable number out of range"
alpar@1
  2016
            "\n", k);
alpar@1
  2017
      /* retrieve information about the specified basic variable x[k]
alpar@1
  2018
         whose objective coefficient c[k] is to be analyzed */
alpar@1
  2019
      if (k <= m)
alpar@1
  2020
      {  row = P->row[k];
alpar@1
  2021
         type = row->type;
alpar@1
  2022
         lb = row->lb;
alpar@1
  2023
         ub = row->ub;
alpar@1
  2024
         coef = 0.0;
alpar@1
  2025
         stat = row->stat;
alpar@1
  2026
         x = row->prim;
alpar@1
  2027
      }
alpar@1
  2028
      else
alpar@1
  2029
      {  col = P->col[k-m];
alpar@1
  2030
         type = col->type;
alpar@1
  2031
         lb = col->lb;
alpar@1
  2032
         ub = col->ub;
alpar@1
  2033
         coef = col->coef;
alpar@1
  2034
         stat = col->stat;
alpar@1
  2035
         x = col->prim;
alpar@1
  2036
      }
alpar@1
  2037
      if (stat != GLP_BS)
alpar@1
  2038
         xerror("glp_analyze_coef: k = %d; non-basic variable not allow"
alpar@1
  2039
            "ed\n", k);
alpar@1
  2040
      /* allocate working arrays */
alpar@1
  2041
      cind = xcalloc(1+m, sizeof(int));
alpar@1
  2042
      cval = xcalloc(1+m, sizeof(double));
alpar@1
  2043
      rind = xcalloc(1+n, sizeof(int));
alpar@1
  2044
      rval = xcalloc(1+n, sizeof(double));
alpar@1
  2045
      /* compute row of the simplex table corresponding to the basic
alpar@1
  2046
         variable x[k] */
alpar@1
  2047
      rlen = glp_eval_tab_row(P, k, rind, rval);
alpar@1
  2048
      xassert(0 <= rlen && rlen <= n);
alpar@1
  2049
      /* perform analysis */
alpar@1
  2050
      for (kase = -1; kase <= +1; kase += 2)
alpar@1
  2051
      {  /* kase < 0 means objective coefficient c[k] is decreasing;
alpar@1
  2052
            kase > 0 means objective coefficient c[k] is increasing */
alpar@1
  2053
         /* note that decreasing c[k] is equivalent to increasing dual
alpar@1
  2054
            variable lambda[k] and vice versa; we need to correctly set
alpar@1
  2055
            the dir flag as required by the routine glp_dual_rtest */
alpar@1
  2056
         if (P->dir == GLP_MIN)
alpar@1
  2057
            dir = - kase;
alpar@1
  2058
         else if (P->dir == GLP_MAX)
alpar@1
  2059
            dir = + kase;
alpar@1
  2060
         else
alpar@1
  2061
            xassert(P != P);
alpar@1
  2062
         /* use the dual ratio test to determine non-basic variable
alpar@1
  2063
            x[q] whose reduced cost d[q] reaches zero bound first */
alpar@1
  2064
         rpiv = glp_dual_rtest(P, rlen, rind, rval, dir, 1e-9);
alpar@1
  2065
         if (rpiv == 0)
alpar@1
  2066
         {  /* nothing limits changing c[k] */
alpar@1
  2067
            lim_coef = (kase < 0 ? -DBL_MAX : +DBL_MAX);
alpar@1
  2068
            q = 0;
alpar@1
  2069
            /* x[k] keeps its current value */
alpar@1
  2070
            new_x = x;
alpar@1
  2071
            goto store;
alpar@1
  2072
         }
alpar@1
  2073
         /* non-basic variable x[q] limits changing coefficient c[k];
alpar@1
  2074
            determine its status and reduced cost d[k] in the current
alpar@1
  2075
            basis */
alpar@1
  2076
         xassert(1 <= rpiv && rpiv <= rlen);
alpar@1
  2077
         q = rind[rpiv];
alpar@1
  2078
         xassert(1 <= q && q <= m+n);
alpar@1
  2079
         if (q <= m)
alpar@1
  2080
         {  row = P->row[q];
alpar@1
  2081
            stat = row->stat;
alpar@1
  2082
            d = row->dual;
alpar@1
  2083
         }
alpar@1
  2084
         else
alpar@1
  2085
         {  col = P->col[q-m];
alpar@1
  2086
            stat = col->stat;
alpar@1
  2087
            d = col->dual;
alpar@1
  2088
         }
alpar@1
  2089
         /* note that delta d[q] = new d[q] - d[q] = - d[q], because
alpar@1
  2090
            new d[q] = 0; delta d[q] = alfa[k,q] * delta c[k], so
alpar@1
  2091
            delta c[k] = delta d[q] / alfa[k,q] = - d[q] / alfa[k,q] */
alpar@1
  2092
         xassert(rval[rpiv] != 0.0);
alpar@1
  2093
         delta = - d / rval[rpiv];
alpar@1
  2094
         /* compute new c[k] = c[k] + delta c[k], which is the limiting
alpar@1
  2095
            value of the objective coefficient c[k] */
alpar@1
  2096
         lim_coef = coef + delta;
alpar@1
  2097
         /* let c[k] continue decreasing/increasing that makes d[q]
alpar@1
  2098
            dual infeasible and forces x[q] to enter the basis;
alpar@1
  2099
            to perform the primal ratio test we need to know in which
alpar@1
  2100
            direction x[q] changes on entering the basis; we determine
alpar@1
  2101
            that analyzing the sign of delta d[q] (see above), since
alpar@1
  2102
            d[q] may be close to zero having wrong sign */
alpar@1
  2103
         /* let, for simplicity, the problem is minimization */
alpar@1
  2104
         if (kase < 0 && rval[rpiv] > 0.0 ||
alpar@1
  2105
             kase > 0 && rval[rpiv] < 0.0)
alpar@1
  2106
         {  /* delta d[q] < 0, so d[q] being non-negative will become
alpar@1
  2107
               negative, so x[q] will increase */
alpar@1
  2108
            dir = +1;
alpar@1
  2109
         }
alpar@1
  2110
         else
alpar@1
  2111
         {  /* delta d[q] > 0, so d[q] being non-positive will become
alpar@1
  2112
               positive, so x[q] will decrease */
alpar@1
  2113
            dir = -1;
alpar@1
  2114
         }
alpar@1
  2115
         /* if the problem is maximization, correct the direction */
alpar@1
  2116
         if (P->dir == GLP_MAX) dir = - dir;
alpar@1
  2117
         /* check that we didn't make a silly mistake */
alpar@1
  2118
         if (dir > 0)
alpar@1
  2119
            xassert(stat == GLP_NL || stat == GLP_NF);
alpar@1
  2120
         else
alpar@1
  2121
            xassert(stat == GLP_NU || stat == GLP_NF);
alpar@1
  2122
         /* compute column of the simplex table corresponding to the
alpar@1
  2123
            non-basic variable x[q] */
alpar@1
  2124
         clen = glp_eval_tab_col(P, q, cind, cval);
alpar@1
  2125
         /* make x[k] temporarily free (unbounded) */
alpar@1
  2126
         if (k <= m)
alpar@1
  2127
         {  row = P->row[k];
alpar@1
  2128
            row->type = GLP_FR;
alpar@1
  2129
            row->lb = row->ub = 0.0;
alpar@1
  2130
         }
alpar@1
  2131
         else
alpar@1
  2132
         {  col = P->col[k-m];
alpar@1
  2133
            col->type = GLP_FR;
alpar@1
  2134
            col->lb = col->ub = 0.0;
alpar@1
  2135
         }
alpar@1
  2136
         /* use the primal ratio test to determine some basic variable
alpar@1
  2137
            which leaves the basis */
alpar@1
  2138
         cpiv = glp_prim_rtest(P, clen, cind, cval, dir, 1e-9);
alpar@1
  2139
         /* restore original bounds of the basic variable x[k] */
alpar@1
  2140
         if (k <= m)
alpar@1
  2141
         {  row = P->row[k];
alpar@1
  2142
            row->type = type;
alpar@1
  2143
            row->lb = lb, row->ub = ub;
alpar@1
  2144
         }
alpar@1
  2145
         else
alpar@1
  2146
         {  col = P->col[k-m];
alpar@1
  2147
            col->type = type;
alpar@1
  2148
            col->lb = lb, col->ub = ub;
alpar@1
  2149
         }
alpar@1
  2150
         if (cpiv == 0)
alpar@1
  2151
         {  /* non-basic variable x[q] can change unlimitedly */
alpar@1
  2152
            if (dir < 0 && rval[rpiv] > 0.0 ||
alpar@1
  2153
                dir > 0 && rval[rpiv] < 0.0)
alpar@1
  2154
            {  /* delta x[k] = alfa[k,q] * delta x[q] < 0 */
alpar@1
  2155
               new_x = -DBL_MAX;
alpar@1
  2156
            }
alpar@1
  2157
            else
alpar@1
  2158
            {  /* delta x[k] = alfa[k,q] * delta x[q] > 0 */
alpar@1
  2159
               new_x = +DBL_MAX;
alpar@1
  2160
            }
alpar@1
  2161
            goto store;
alpar@1
  2162
         }
alpar@1
  2163
         /* some basic variable x[p] limits changing non-basic variable
alpar@1
  2164
            x[q] in the adjacent basis */
alpar@1
  2165
         xassert(1 <= cpiv && cpiv <= clen);
alpar@1
  2166
         p = cind[cpiv];
alpar@1
  2167
         xassert(1 <= p && p <= m+n);
alpar@1
  2168
         xassert(p != k);
alpar@1
  2169
         if (p <= m)
alpar@1
  2170
         {  row = P->row[p];
alpar@1
  2171
            xassert(row->stat == GLP_BS);
alpar@1
  2172
            ll = glp_get_row_lb(P, row->i);
alpar@1
  2173
            uu = glp_get_row_ub(P, row->i);
alpar@1
  2174
            xx = row->prim;
alpar@1
  2175
         }
alpar@1
  2176
         else
alpar@1
  2177
         {  col = P->col[p-m];
alpar@1
  2178
            xassert(col->stat == GLP_BS);
alpar@1
  2179
            ll = glp_get_col_lb(P, col->j);
alpar@1
  2180
            uu = glp_get_col_ub(P, col->j);
alpar@1
  2181
            xx = col->prim;
alpar@1
  2182
         }
alpar@1
  2183
         /* determine delta x[p] = new x[p] - x[p] */
alpar@1
  2184
         if (dir < 0 && cval[cpiv] > 0.0 ||
alpar@1
  2185
             dir > 0 && cval[cpiv] < 0.0)
alpar@1
  2186
         {  /* delta x[p] < 0, so x[p] goes toward its lower bound */
alpar@1
  2187
            xassert(ll != -DBL_MAX);
alpar@1
  2188
            delta = ll - xx;
alpar@1
  2189
         }
alpar@1
  2190
         else
alpar@1
  2191
         {  /* delta x[p] > 0, so x[p] goes toward its upper bound */
alpar@1
  2192
            xassert(uu != +DBL_MAX);
alpar@1
  2193
            delta = uu - xx;
alpar@1
  2194
         }
alpar@1
  2195
         /* compute new x[k] = x[k] + alfa[k,q] * delta x[q], where
alpar@1
  2196
            delta x[q] = delta x[p] / alfa[p,q] */
alpar@1
  2197
         xassert(cval[cpiv] != 0.0);
alpar@1
  2198
         new_x = x + (rval[rpiv] / cval[cpiv]) * delta;
alpar@1
  2199
store:   /* store analysis results */
alpar@1
  2200
         if (kase < 0)
alpar@1
  2201
         {  if (coef1 != NULL) *coef1 = lim_coef;
alpar@1
  2202
            if (var1 != NULL) *var1 = q;
alpar@1
  2203
            if (value1 != NULL) *value1 = new_x;
alpar@1
  2204
         }
alpar@1
  2205
         else
alpar@1
  2206
         {  if (coef2 != NULL) *coef2 = lim_coef;
alpar@1
  2207
            if (var2 != NULL) *var2 = q;
alpar@1
  2208
            if (value2 != NULL) *value2 = new_x;
alpar@1
  2209
         }
alpar@1
  2210
      }
alpar@1
  2211
      /* free working arrays */
alpar@1
  2212
      xfree(cind);
alpar@1
  2213
      xfree(cval);
alpar@1
  2214
      xfree(rind);
alpar@1
  2215
      xfree(rval);
alpar@1
  2216
      return;
alpar@1
  2217
}
alpar@1
  2218
alpar@1
  2219
/* eof */