glpluf.c

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

*  strategy, trying to choose such element v[p,q], which satisfies to
*  the stability condition (see above) and has smallest Markowitz cost
*  (nr[p]-1) * (nc[q]-1), where nr[p] and nc[q] are numbers of non-zero
*  elements, respectively, in the p-th row and in the q-th column of the
*  active submatrix.
* 
*  In order to reduce the search, i.e. not to walk through all elements
*  of the active submatrix, the routine exploits a technique proposed by
*  I.Duff. This technique is based on using the sets R[len] and C[len]
*  of active rows and columns.
* 
*  If the pivot element v[p,q] has been chosen, the routine stores its
*  indices to the locations *p and *q and returns zero. Otherwise, if
*  the active submatrix is empty and therefore the pivot element can't
*  be chosen, the routine returns non-zero. */

static int find_pivot(LUF *luf, int *_p, int *_q)
{     int n = luf->n;
      int *vr_ptr = luf->vr_ptr;
      int *vr_len = luf->vr_len;
      int *vc_ptr = luf->vc_ptr;
      int *vc_len = luf->vc_len;
      int *sv_ind = luf->sv_ind;
      double *sv_val = luf->sv_val;
      double *vr_max = luf->vr_max;
      int *rs_head = luf->rs_head;
      int *rs_next = luf->rs_next;
      int *cs_head = luf->cs_head;
      int *cs_prev = luf->cs_prev;
      int *cs_next = luf->cs_next;
      double piv_tol = luf->piv_tol;
      int piv_lim = luf->piv_lim;
      int suhl = luf->suhl;
      int p, q, len, i, i_beg, i_end, i_ptr, j, j_beg, j_end, j_ptr,
         ncand, next_j, min_p, min_q, min_len;
      double best, cost, big, temp;
      /* initially no pivot candidates have been found so far */
      p = q = 0, best = DBL_MAX, ncand = 0;
      /* if in the active submatrix there is a column that has the only
         non-zero (column singleton), choose it as pivot */
      j = cs_head[1];
      if (j != 0)
      {  xassert(vc_len[j] == 1);
         p = sv_ind[vc_ptr[j]], q = j;
         goto done;
      }
      /* if in the active submatrix there is a row that has the only
         non-zero (row singleton), choose it as pivot */
      i = rs_head[1];
      if (i != 0)
      {  xassert(vr_len[i] == 1);
         p = i, q = sv_ind[vr_ptr[i]];
         goto done;
      }
      /* there are no singletons in the active submatrix; walk through
         other non-empty rows and columns */
      for (len = 2; len <= n; len++)
      {  /* consider active columns that have len non-zeros */
         for (j = cs_head[len]; j != 0; j = next_j)
         {  /* the j-th column has len non-zeros */
            j_beg = vc_ptr[j];
            j_end = j_beg + vc_len[j] - 1;
            /* save pointer to the next column with the same length */
            next_j = cs_next[j];
            /* find an element in the j-th column, which is placed in a
               row with minimal number of non-zeros and satisfies to the
               stability condition (such element may not exist) */
            min_p = min_q = 0, min_len = INT_MAX;
            for (j_ptr = j_beg; j_ptr <= j_end; j_ptr++)
            {  /* get row index of v[i,j] */
               i = sv_ind[j_ptr];
               i_beg = vr_ptr[i];
               i_end = i_beg + vr_len[i] - 1;
               /* if the i-th row is not shorter than that one, where
                  minimal element is currently placed, skip v[i,j] */
               if (vr_len[i] >= min_len) continue;
               /* determine the largest of absolute values of elements
                  in the i-th row */
               big = vr_max[i];
               if (big < 0.0)
               {  /* the largest value is unknown yet; compute it */
                  for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++)
                  {  temp = sv_val[i_ptr];
                     if (temp < 0.0) temp = - temp;
                     if (big < temp) big = temp;
                  }
                  vr_max[i] = big;
               }
               /* find v[i,j] in the i-th row */
               for (i_ptr = vr_ptr[i]; sv_ind[i_ptr] != j; i_ptr++);
               xassert(i_ptr <= i_end);
               /* if v[i,j] doesn't satisfy to the stability condition,
                  skip it */
               temp = sv_val[i_ptr];
               if (temp < 0.0) temp = - temp;
               if (temp < piv_tol * big) continue;
               /* v[i,j] is better than the current minimal element */
               min_p = i, min_q = j, min_len = vr_len[i];
               /* if Markowitz cost of the current minimal element is
                  not greater than (len-1)**2, it can be chosen right
                  now; this heuristic reduces the search and works well
                  in many cases */
               if (min_len <= len)
               {  p = min_p, q = min_q;
                  goto done;
               }
            }
            /* the j-th column has been scanned */
            if (min_p != 0)
            {  /* the minimal element is a next pivot candidate */
               ncand++;
               /* compute its Markowitz cost */
               cost = (double)(min_len - 1) * (double)(len - 1);
               /* choose between the minimal element and the current
                  candidate */
               if (cost < best) p = min_p, q = min_q, best = cost;
               /* if piv_lim candidates have been considered, there are
                  doubts that a much better candidate exists; therefore
                  it's time to terminate the search */
               if (ncand == piv_lim) goto done;
            }
            else
            {  /* the j-th column has no elements, which satisfy to the
                  stability condition; Uwe Suhl suggests to exclude such
                  column from the further consideration until it becomes
                  a column singleton; in hard cases this significantly
                  reduces a time needed for pivot searching */
               if (suhl)
               {  /* remove the j-th column from the active set */
                  if (cs_prev[j] == 0)
                     cs_head[len] = cs_next[j];
                  else
                     cs_next[cs_prev[j]] = cs_next[j];
                  if (cs_next[j] == 0)
                     /* nop */;
                  else
                     cs_prev[cs_next[j]] = cs_prev[j];
                  /* the following assignment is used to avoid an error
                     when the routine eliminate (see below) will try to
                     remove the j-th column from the active set */
                  cs_prev[j] = cs_next[j] = j;
               }
            }
         }
         /* consider active rows that have len non-zeros */
         for (i = rs_head[len]; i != 0; i = rs_next[i])
         {  /* the i-th row has len non-zeros */
            i_beg = vr_ptr[i];
            i_end = i_beg + vr_len[i] - 1;
            /* determine the largest of absolute values of elements in
               the i-th row */
            big = vr_max[i];
            if (big < 0.0)
            {  /* the largest value is unknown yet; compute it */
               for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++)
               {  temp = sv_val[i_ptr];
                  if (temp < 0.0) temp = - temp;
                  if (big < temp) big = temp;
               }
               vr_max[i] = big;
            }
            /* find an element in the i-th row, which is placed in a
               column with minimal number of non-zeros and satisfies to
               the stability condition (such element always exists) */
            min_p = min_q = 0, min_len = INT_MAX;
            for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++)
            {  /* get column index of v[i,j] */
               j = sv_ind[i_ptr];
               /* if the j-th column is not shorter than that one, where
                  minimal element is currently placed, skip v[i,j] */
               if (vc_len[j] >= min_len) continue;
               /* if v[i,j] doesn't satisfy to the stability condition,
                  skip it */
               temp = sv_val[i_ptr];
               if (temp < 0.0) temp = - temp;
               if (temp < piv_tol * big) continue;
               /* v[i,j] is better than the current minimal element */
               min_p = i, min_q = j, min_len = vc_len[j];
               /* if Markowitz cost of the current minimal element is
                  not greater than (len-1)**2, it can be chosen right
                  now; this heuristic reduces the search and works well
                  in many cases */
               if (min_len <= len)
               {  p = min_p, q = min_q;
                  goto done;
               }
            }
            /* the i-th row has been scanned */
            if (min_p != 0)
            {  /* the minimal element is a next pivot candidate */
               ncand++;
               /* compute its Markowitz cost */
               cost = (double)(len - 1) * (double)(min_len - 1);
               /* choose between the minimal element and the current
                  candidate */
               if (cost < best) p = min_p, q = min_q, best = cost;
               /* if piv_lim candidates have been considered, there are
                  doubts that a much better candidate exists; therefore
                  it's time to terminate the search */
               if (ncand == piv_lim) goto done;
            }
            else
            {  /* this can't be because this can never be */
               xassert(min_p != min_p);
            }
         }
      }
done: /* bring the pivot to the factorizing routine */
      *_p = p, *_q = q;
      return (p == 0);
}

/***********************************************************************
*  eliminate - perform gaussian elimination.
* 
*  This routine performs elementary gaussian transformations in order
*  to eliminate subdiagonal elements in the k-th column of the matrix
*  U = P*V*Q using the pivot element u[k,k], where k is the number of
*  the current elimination step.
* 
*  The parameters p and q are, respectively, row and column indices of
*  the element v[p,q], which corresponds to the element u[k,k].
* 
*  Each time when the routine applies the elementary transformation to
*  a non-pivot row of the matrix V, it stores the corresponding element
*  to the matrix F in order to keep the main equality A = F*V.
* 
*  The routine assumes that on entry the matrices L = P*F*inv(P) and
*  U = P*V*Q are the following:
* 
*        1       k                  1       k         n
*     1  1 . . . . . . . . .     1  x x x x x x x x x x
*        x 1 . . . . . . . .        . x x x x x x x x x
*        x x 1 . . . . . . .        . . x x x x x x x x
*        x x x 1 . . . . . .        . . . x x x x x x x
*     k  x x x x 1 . . . . .     k  . . . . * * * * * *
*        x x x x _ 1 . . . .        . . . . # * * * * *
*        x x x x _ . 1 . . .        . . . . # * * * * *
*        x x x x _ . . 1 . .        . . . . # * * * * *
*        x x x x _ . . . 1 .        . . . . # * * * * *
*     n  x x x x _ . . . . 1     n  . . . . # * * * * *
* 
*             matrix L                   matrix U
* 
*  where rows and columns of the matrix U with numbers k, k+1, ..., n
*  form the active submatrix (eliminated elements are marked by '#' and
*  other elements of the active submatrix are marked by '*'). Note that
*  each eliminated non-zero element u[i,k] of the matrix U gives the
*  corresponding element l[i,k] of the matrix L (marked by '_').
* 
*  Actually all operations are performed on the matrix V. Should note
*  that the row-wise representation corresponds to the matrix V, but the
*  column-wise representation corresponds to the active submatrix of the
*  matrix V, i.e. elements of the matrix V, which doesn't belong to the
*  active submatrix, are missing from the column linked lists.
* 
*  Let u[k,k] = v[p,q] be the pivot. In order to eliminate subdiagonal
*  elements u[i',k] = v[i,q], i' = k+1, k+2, ..., n, the routine applies
*  the following elementary gaussian transformations:
* 
*     (i-th row of V) := (i-th row of V) - f[i,p] * (p-th row of V),
* 
*  where f[i,p] = v[i,q] / v[p,q] is a gaussian multiplier.
*
*  Additionally, in order to keep the main equality A = F*V, each time
*  when the routine applies the transformation to i-th row of the matrix
*  V, it also adds f[i,p] as a new element to the matrix F.
*
*  IMPORTANT: On entry the working arrays flag and work should contain
*  zeros. This status is provided by the routine on exit.
*
*  If no error occured, the routine returns zero. Otherwise, in case of
*  overflow of the sparse vector area, the routine returns non-zero. */

static int eliminate(LUF *luf, int p, int q)
{     int n = luf->n;
      int *fc_ptr = luf->fc_ptr;
      int *fc_len = luf->fc_len;
      int *vr_ptr = luf->vr_ptr;
      int *vr_len = luf->vr_len;
      int *vr_cap = luf->vr_cap;
      double *vr_piv = luf->vr_piv;
      int *vc_ptr = luf->vc_ptr;
      int *vc_len = luf->vc_len;
      int *vc_cap = luf->vc_cap;
      int *sv_ind = luf->sv_ind;
      double *sv_val = luf->sv_val;
      int *sv_prev = luf->sv_prev;
      int *sv_next = luf->sv_next;
      double *vr_max = luf->vr_max;
      int *rs_head = luf->rs_head;
      int *rs_prev = luf->rs_prev;
      int *rs_next = luf->rs_next;
      int *cs_head = luf->cs_head;
      int *cs_prev = luf->cs_prev;
      int *cs_next = luf->cs_next;
      int *flag = luf->flag;
      double *work = luf->work;
      double eps_tol = luf->eps_tol;
      /* at this stage the row-wise representation of the matrix F is
         not used, so fr_len can be used as a working array */
      int *ndx = luf->fr_len;
      int ret = 0;
      int len, fill, i, i_beg, i_end, i_ptr, j, j_beg, j_end, j_ptr, k,
         p_beg, p_end, p_ptr, q_beg, q_end, q_ptr;
      double fip, val, vpq, temp;
      xassert(1 <= p && p <= n);
      xassert(1 <= q && q <= n);
      /* remove the p-th (pivot) row from the active set; this row will
         never return there */
      if (rs_prev[p] == 0)
         rs_head[vr_len[p]] = rs_next[p];
      else
         rs_next[rs_prev[p]] = rs_next[p];
      if (rs_next[p] == 0)
         ;
      else
         rs_prev[rs_next[p]] = rs_prev[p];
      /* remove the q-th (pivot) column from the active set; this column
         will never return there */
      if (cs_prev[q] == 0)
         cs_head[vc_len[q]] = cs_next[q];
      else
         cs_next[cs_prev[q]] = cs_next[q];
      if (cs_next[q] == 0)
         ;
      else
         cs_prev[cs_next[q]] = cs_prev[q];
      /* find the pivot v[p,q] = u[k,k] in the p-th row */
      p_beg = vr_ptr[p];
      p_end = p_beg + vr_len[p] - 1;
      for (p_ptr = p_beg; sv_ind[p_ptr] != q; p_ptr++) /* nop */;
      xassert(p_ptr <= p_end);