glplpf.c

著名的大规模线性规划求解器源码GLPK.C语言版本,可以修剪.内有详细帮助文档.

      lpf->v_size = v_size;
      lpf->v_ind = xcalloc(1+v_size, sizeof(int));
      lpf->v_val = xcalloc(1+v_size, sizeof(double));
      xassert(used >= 0);
      memcpy(&lpf->v_ind[1], &v_ind[1], used * sizeof(int));
      memcpy(&lpf->v_val[1], &v_val[1], used * sizeof(double));
      xfree(v_ind);
      xfree(v_val);
      return;
}

/***********************************************************************
*  NAME
*
*  lpf_update_it - update LP basis factorization
*
*  SYNOPSIS
*
*  #include "glplpf.h"
*  int lpf_update_it(LPF *lpf, int j, int bh, int len, const int ind[],
*     const double val[]);
*
*  DESCRIPTION
*
*  The routine lpf_update_it updates the factorization of the basis
*  matrix B after replacing its j-th column by a new vector.
*
*  The parameter j specifies the number of column of B, which has been
*  replaced, 1 <= j <= m, where m is the order of B.
*
*  The parameter bh specifies the basis header entry for the new column
*  of B, which is the number of the new column in some original matrix.
*  This parameter is optional and can be specified as 0.
*
*  Row indices and numerical values of non-zero elements of the new
*  column of B should be placed in locations ind[1], ..., ind[len] and
*  val[1], ..., val[len], resp., where len is the number of non-zeros
*  in the column. Neither zero nor duplicate elements are allowed.
*
*  RETURNS
*
*  0  The factorization has been successfully updated.
*
*  LPF_ESING
*     New basis B is singular within the working precision.
*
*  LPF_ELIMIT
*     Maximal number of additional rows and columns has been reached.
*
*  BACKGROUND
*
*  Let j-th column of the current basis matrix B have to be replaced by
*  a new column a. This replacement is equivalent to removing the old
*  j-th column by fixing it at zero and introducing the new column as
*  follows:
*
*                   ( B   F^| a )
*     ( B  F^)      (       |   )
*     (      ) ---> ( G^  H^| 0 )
*     ( G^ H^)      (-------+---)
*                   ( e'j 0 | 0 )
*
*  where ej is a unit vector with 1 in j-th position which used to fix
*  the old j-th column of B (at zero). Then using the main equality we
*  have:
*
*     ( B   F^| a )            ( B0  F | f )
*     (       |   )   ( P  0 ) (       |   ) ( Q  0 )
*     ( G^  H^| 0 ) = (      ) ( G   H | g ) (      ) =
*     (-------+---)   ( 0  1 ) (-------+---) ( 0  1 )
*     ( e'j 0 | 0 )            ( v'  w'| 0 )
*
*       [   ( B0  F )|   ( f ) ]            [   ( B0 F )  |   ( f ) ]
*       [ P (       )| P (   ) ] ( Q  0 )   [ P (      ) Q| P (   ) ]
*     = [   ( G   H )|   ( g ) ] (      ) = [   ( G  H )  |   ( g ) ]
*       [------------+-------- ] ( 0  1 )   [-------------+---------]
*       [   ( v'  w')|     0   ]            [   ( v' w') Q|    0    ]
*
*  where:
*
*     ( a )     ( f )      ( f )        ( a )
*     (   ) = P (   )  =>  (   ) = P' * (   )
*     ( 0 )     ( g )      ( g )        ( 0 )
*
*                                 ( ej )      ( v )    ( v )     ( ej )
*     ( e'j  0 ) = ( v' w' ) Q => (    ) = Q' (   ) => (   ) = Q (    )
*                                 ( 0  )      ( w )    ( w )     ( 0  )
*
*  On the other hand:
*
*              ( B0| F  f )
*     ( P  0 ) (---+------) ( Q  0 )         ( B0    new F )
*     (      ) ( G | H  g ) (      ) = new P (             ) new Q
*     ( 0  1 ) (   |      ) ( 0  1 )         ( new G new H )
*              ( v'| w' 0 )
*
*  where:
*                               ( G )           ( H  g )
*     new F = ( F  f ), new G = (   ),  new H = (      ),
*                               ( v')           ( w' 0 )
*
*             ( P  0 )           ( Q  0 )
*     new P = (      ) , new Q = (      ) .
*             ( 0  1 )           ( 0  1 )
*
*  The factorization structure for the new augmented matrix remains the
*  same, therefore:
*
*           ( B0    new F )         ( L0     0 ) ( U0 new R )
*     new P (             ) new Q = (          ) (          )
*           ( new G new H )         ( new S  I ) ( 0  new C )
*
*  where:
*
*     new F = L0 * new R  =>
*
*     new R = inv(L0) * new F = inv(L0) * (F  f) = ( R  inv(L0)*f )
*
*     new G = new S * U0  =>
*
*                               ( G )             (     S      )
*     new S = new G * inv(U0) = (   ) * inv(U0) = (            )
*                               ( v')             ( v'*inv(U0) )
*
*     new H = new S * new R + new C  =>
*
*     new C = new H - new S * new R =
*
*             ( H  g )   (     S      )
*           = (      ) - (            ) * ( R  inv(L0)*f ) =
*             ( w' 0 )   ( v'*inv(U0) )
*
*             ( H - S*R           g - S*inv(L0)*f      )   ( C  x )
*           = (                                        ) = (      )
*             ( w'- v'*inv(U0)*R  -v'*inv(U0)*inv(L0)*f)   ( y' z )
*
*  Note that new C is resulted by expanding old C with new column x,
*  row y', and diagonal element z, where:
*
*     x = g - S * inv(L0) * f = g - S * (new column of R)
*
*     y = w - R'* inv(U'0)* v = w - R'* (new row of S)
*
*     z = - (new row of S) * (new column of R)
*
*  Finally, to replace old B by new B we have to permute j-th and last
*  (just added) columns of the matrix
*
*     ( B   F^| a )
*     (       |   )
*     ( G^  H^| 0 )
*     (-------+---)
*     ( e'j 0 | 0 )
*
*  and to keep the main equality do the same for matrix Q. */

int lpf_update_it(LPF *lpf, int j, int bh, int len, const int ind[],
      const double val[])
{     int m0 = lpf->m0;
      int m = lpf->m;
#if _GLPLPF_DEBUG
      double *B = lpf->B;
#endif
      int n = lpf->n;
      int *R_ptr = lpf->R_ptr;
      int *R_len = lpf->R_len;
      int *S_ptr = lpf->S_ptr;
      int *S_len = lpf->S_len;
      int *P_row = lpf->P_row;
      int *P_col = lpf->P_col;
      int *Q_row = lpf->Q_row;
      int *Q_col = lpf->Q_col;
      int v_ptr = lpf->v_ptr;
      int *v_ind = lpf->v_ind;
      double *v_val = lpf->v_val;
      double *a = lpf->work2; /* new column */
      double *fg = lpf->work1, *f = fg, *g = fg + m0;
      double *vw = lpf->work2, *v = vw, *w = vw + m0;
      double *x = g, *y = w, z;
      int i, ii, k, ret;
      xassert(bh == bh);
      if (!lpf->valid)
         xfault("lpf_update_it: the factorization is not valid\n");
      if (!(1 <= j && j <= m))
         xfault("lpf_update_it: j = %d; column number out of range\n",
            j);
      xassert(0 <= m && m <= m0 + n);
      /* check if the basis factorization can be expanded */
      if (n == lpf->n_max)
      {  lpf->valid = 0;
         ret = LPF_ELIMIT;
         goto done;
      }
      /* convert new j-th column of B to dense format */
      for (i = 1; i <= m; i++)
         a[i] = 0.0;
      for (k = 1; k <= len; k++)
      {  i = ind[k];
         if (!(1 <= i && i <= m))
            xfault("lpf_update_it: ind[%d] = %d; row number out of rang"
               "e\n", k, i);
         if (a[i] != 0.0)
            xfault("lpf_update_it: ind[%d] = %d; duplicate row index no"
               "t allowed\n", k, i);
         if (val[k] == 0.0)
            xfault("lpf_update_it: val[%d] = %g; zero element not allow"
               "ed\n", k, val[k]);
         a[i] = val[k];
      }
#if _GLPLPF_DEBUG
      /* change column in the basis matrix for debugging */
      for (i = 1; i <= m; i++)
         B[(i - 1) * m + j] = a[i];
#endif
      /* (f g) := inv(P) * (a 0) */
      for (i = 1; i <= m0+n; i++)
         fg[i] = ((ii = P_col[i]) <= m ? a[ii] : 0.0);
      /* (v w) := Q * (ej 0) */
      for (i = 1; i <= m0+n; i++) vw[i] = 0.0;
      vw[Q_col[j]] = 1.0;
      /* f1 := inv(L0) * f (new column of R) */
      luf_f_solve(lpf->luf, 0, f);
      /* v1 := inv(U'0) * v (new row of S) */
      luf_v_solve(lpf->luf, 1, v);
      /* we need at most 2 * m0 available locations in the SVA to store
         new column of matrix R and new row of matrix S */
      if (lpf->v_size < v_ptr + m0 + m0)
      {  enlarge_sva(lpf, v_ptr + m0 + m0);
         v_ind = lpf->v_ind;
         v_val = lpf->v_val;
      }
      /* store new column of R */
      R_ptr[n+1] = v_ptr;
      for (i = 1; i <= m0; i++)
      {  if (f[i] != 0.0)
            v_ind[v_ptr] = i, v_val[v_ptr] = f[i], v_ptr++;
      }
      R_len[n+1] = v_ptr - lpf->v_ptr;
      lpf->v_ptr = v_ptr;
      /* store new row of S */
      S_ptr[n+1] = v_ptr;
      for (i = 1; i <= m0; i++)
      {  if (v[i] != 0.0)
            v_ind[v_ptr] = i, v_val[v_ptr] = v[i], v_ptr++;
      }
      S_len[n+1] = v_ptr - lpf->v_ptr;
      lpf->v_ptr = v_ptr;
      /* x := g - S * f1 (new column of C) */
      s_prod(lpf, x, -1.0, f);
      /* y := w - R' * v1 (new row of C) */
      rt_prod(lpf, y, -1.0, v);
      /* z := - v1 * f1 (new diagonal element of C) */
      z = 0.0;
      for (i = 1; i <= m0; i++) z -= v[i] * f[i];
      /* update factorization of new matrix C */
      switch (scf_update_exp(lpf->scf, x, y, z))
      {  case 0:
            break;
         case SCF_ESING:
            lpf->valid = 0;
            ret = LPF_ESING;
            goto done;
         case SCF_ELIMIT:
            xassert(lpf != lpf);
         default:
            xassert(lpf != lpf);
      }
      /* expand matrix P */
      P_row[m0+n+1] = P_col[m0+n+1] = m0+n+1;
      /* expand matrix Q */
      Q_row[m0+n+1] = Q_col[m0+n+1] = m0+n+1;
      /* permute j-th and last (just added) column of matrix Q */
      i = Q_col[j], ii = Q_col[m0+n+1];
      Q_row[i] = m0+n+1, Q_col[m0+n+1] = i;
      Q_row[ii] = j, Q_col[j] = ii;
      /* increase the number of additional rows and columns */
      lpf->n++;
      xassert(lpf->n <= lpf->n_max);
      /* the factorization has been successfully updated */
      ret = 0;
done: /* return to the calling program */
      return ret;
}

/***********************************************************************
*  NAME
*
*  lpf_delete_it - delete LP basis factorization
*
*  SYNOPSIS
*
*  #include "glplpf.h"
*  void lpf_delete_it(LPF *lpf)
*
*  DESCRIPTION
*
*  The routine lpf_delete_it deletes LP basis factorization specified
*  by the parameter lpf and frees all memory allocated to this program
*  object. */

void lpf_delete_it(LPF *lpf)
{     luf_delete_it(lpf->luf);
#if _GLPLPF_DEBUG
      if (lpf->B != NULL) xfree(lpf->B);
#else
      xassert(lpf->B == NULL);
#endif
      if (lpf->R_ptr != NULL) xfree(lpf->R_ptr);
      if (lpf->R_len != NULL) xfree(lpf->R_len);
      if (lpf->S_ptr != NULL) xfree(lpf->S_ptr);
      if (lpf->S_len != NULL) xfree(lpf->S_len);
      if (lpf->scf != NULL) scf_delete_it(lpf->scf);
      if (lpf->P_row != NULL) xfree(lpf->P_row);
      if (lpf->P_col != NULL) xfree(lpf->P_col);
      if (lpf->Q_row != NULL) xfree(lpf->Q_row);
      if (lpf->Q_col != NULL) xfree(lpf->Q_col);
      if (lpf->v_ind != NULL) xfree(lpf->v_ind);
      if (lpf->v_val != NULL) xfree(lpf->v_val);
      if (lpf->work1 != NULL) xfree(lpf->work1);
      if (lpf->work2 != NULL) xfree(lpf->work2);
      xfree(lpf);
      return;
}

/* eof */