glplux.c

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

            /* compute its Markowitz cost */
            cost = (double)(len - 1) * (double)(min_len - 1);
            /* choose between the current candidate and this element */
            if (cost < best) piv = some, best = cost;
            /* if piv_lim candidates have been considered, there is a
               doubt that a much better candidate exists; therefore it
               is the time to terminate the search */
            if (ncand == piv_lim) goto done;
         }
      }
done: /* bring the pivot v[p,q] to the factorizing routine */
      return piv;
}

/*----------------------------------------------------------------------
// 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 parameter piv specifies the pivot element v[p,q] = 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. */

static void eliminate(LUX *lux, LUXWKA *wka, LUXELM *piv, int flag[],
      mpq_t work[])
{     DMP *pool = lux->pool;
      LUXELM **F_row = lux->F_row;
      LUXELM **F_col = lux->F_col;
      mpq_t *V_piv = lux->V_piv;
      LUXELM **V_row = lux->V_row;
      LUXELM **V_col = lux->V_col;
      int *R_len = wka->R_len;
      int *R_head = wka->R_head;
      int *R_prev = wka->R_prev;
      int *R_next = wka->R_next;
      int *C_len = wka->C_len;
      int *C_head = wka->C_head;
      int *C_prev = wka->C_prev;
      int *C_next = wka->C_next;
      LUXELM *fip, *vij, *vpj, *viq, *next;
      mpq_t temp;
      int i, j, p, q;
      mpq_init(temp);
      /* determine row and column indices of the pivot v[p,q] */
      xassert(piv != NULL);
      p = piv->i, q = piv->j;
      /* remove p-th (pivot) row from the active set; it will never
         return there */
      if (R_prev[p] == 0)
         R_head[R_len[p]] = R_next[p];
      else
         R_next[R_prev[p]] = R_next[p];
      if (R_next[p] == 0)
         ;
      else
         R_prev[R_next[p]] = R_prev[p];
      /* remove q-th (pivot) column from the active set; it will never
         return there */
      if (C_prev[q] == 0)
         C_head[C_len[q]] = C_next[q];
      else
         C_next[C_prev[q]] = C_next[q];
      if (C_next[q] == 0)
         ;
      else
         C_prev[C_next[q]] = C_prev[q];
      /* store the pivot value in a separate array */
      mpq_set(V_piv[p], piv->val);
      /* remove the pivot from p-th row */
      if (piv->r_prev == NULL)
         V_row[p] = piv->r_next;
      else
         piv->r_prev->r_next = piv->r_next;
      if (piv->r_next == NULL)
         ;
      else
         piv->r_next->r_prev = piv->r_prev;
      R_len[p]--;
      /* remove the pivot from q-th column */
      if (piv->c_prev == NULL)
         V_col[q] = piv->c_next;
      else
         piv->c_prev->c_next = piv->c_next;
      if (piv->c_next == NULL)
         ;
      else
         piv->c_next->c_prev = piv->c_prev;
      C_len[q]--;
      /* free the space occupied by the pivot */
      mpq_clear(piv->val);
      dmp_free_atom(pool, piv, sizeof(LUXELM));
      /* walk through p-th (pivot) row, which already does not contain
         the pivot v[p,q], and do the following... */
      for (vpj = V_row[p]; vpj != NULL; vpj = vpj->r_next)
      {  /* get column index of v[p,j] */
         j = vpj->j;
         /* store v[p,j] in the working array */
         flag[j] = 1;
         mpq_set(work[j], vpj->val);
         /* remove j-th column from the active set; it will return there
            later with a new length */
         if (C_prev[j] == 0)
            C_head[C_len[j]] = C_next[j];
         else
            C_next[C_prev[j]] = C_next[j];
         if (C_next[j] == 0)
            ;
         else
            C_prev[C_next[j]] = C_prev[j];
         /* v[p,j] leaves the active submatrix, so remove it from j-th
            column; however, v[p,j] is kept in p-th row */
         if (vpj->c_prev == NULL)
            V_col[j] = vpj->c_next;
         else
            vpj->c_prev->c_next = vpj->c_next;
         if (vpj->c_next == NULL)
            ;
         else
            vpj->c_next->c_prev = vpj->c_prev;
         C_len[j]--;
      }
      /* now walk through q-th (pivot) column, which already does not
         contain the pivot v[p,q], and perform gaussian elimination */
      while (V_col[q] != NULL)
      {  /* element v[i,q] has to be eliminated */
         viq = V_col[q];
         /* get row index of v[i,q] */
         i = viq->i;
         /* remove i-th row from the active set; later it will return
            there with a new length */
         if (R_prev[i] == 0)
            R_head[R_len[i]] = R_next[i];
         else
            R_next[R_prev[i]] = R_next[i];
         if (R_next[i] == 0)
            ;
         else
            R_prev[R_next[i]] = R_prev[i];
         /* compute gaussian multiplier f[i,p] = v[i,q] / v[p,q] and
            store it in the matrix F */
         fip = dmp_get_atom(pool, sizeof(LUXELM));
         fip->i = i, fip->j = p;
         mpq_init(fip->val);
         mpq_div(fip->val, viq->val, V_piv[p]);
         fip->r_prev = NULL;
         fip->r_next = F_row[i];
         fip->c_prev = NULL;
         fip->c_next = F_col[p];
         if (fip->r_next != NULL) fip->r_next->r_prev = fip;
         if (fip->c_next != NULL) fip->c_next->c_prev = fip;
         F_row[i] = F_col[p] = fip;
         /* v[i,q] has to be eliminated, so remove it from i-th row */
         if (viq->r_prev == NULL)
            V_row[i] = viq->r_next;
         else
            viq->r_prev->r_next = viq->r_next;
         if (viq->r_next == NULL)
            ;
         else
            viq->r_next->r_prev = viq->r_prev;
         R_len[i]--;
         /* and also from q-th column */
         V_col[q] = viq->c_next;
         C_len[q]--;
         /* free the space occupied by v[i,q] */
         mpq_clear(viq->val);
         dmp_free_atom(pool, viq, sizeof(LUXELM));
         /* perform gaussian transformation:
            (i-th row) := (i-th row) - f[i,p] * (p-th row)
            note that now p-th row, which is in the working array,
            does not contain the pivot v[p,q], and i-th row does not
            contain the element v[i,q] to be eliminated */
         /* walk through i-th row and transform existing non-zero
            elements */
         for (vij = V_row[i]; vij != NULL; vij = next)
         {  next = vij->r_next;
            /* get column index of v[i,j] */
            j = vij->j;
            /* v[i,j] := v[i,j] - f[i,p] * v[p,j] */
            if (flag[j])
            {  /* v[p,j] != 0 */
               flag[j] = 0;
               mpq_mul(temp, fip->val, work[j]);
               mpq_sub(vij->val, vij->val, temp);
               if (mpq_sgn(vij->val) == 0)
               {  /* new v[i,j] is zero, so remove it from the active
                     submatrix */
                  /* remove v[i,j] from i-th row */
                  if (vij->r_prev == NULL)
                     V_row[i] = vij->r_next;
                  else
                     vij->r_prev->r_next = vij->r_next;
                  if (vij->r_next == NULL)
                     ;
                  else
                     vij->r_next->r_prev = vij->r_prev;
                  R_len[i]--;
                  /* remove v[i,j] from j-th column */
                  if (vij->c_prev == NULL)
                     V_col[j] = vij->c_next;
                  else
                     vij->c_prev->c_next = vij->c_next;
                  if (vij->c_next == NULL)
                     ;
                  else
                     vij->c_next->c_prev = vij->c_prev;
                  C_len[j]--;
                  /* free the space occupied by v[i,j] */
                  mpq_clear(vij->val);
                  dmp_free_atom(pool, vij, sizeof(LUXELM));
               }
            }
         }
         /* now flag is the pattern of the set v[p,*] \ v[i,*] */
         /* walk through p-th (pivot) row and create new elements in
            i-th row, which appear due to fill-in */
         for (vpj = V_row[p]; vpj != NULL; vpj = vpj->r_next)
         {  j = vpj->j;
            if (flag[j])
            {  /* create new non-zero v[i,j] = 0 - f[i,p] * v[p,j] and
                  add it to i-th row and j-th column */
               vij = dmp_get_atom(pool, sizeof(LUXELM));
               vij->i = i, vij->j = j;
               mpq_init(vij->val);
               mpq_mul(vij->val, fip->val, work[j]);
               mpq_neg(vij->val, vij->val);
               vij->r_prev = NULL;
               vij->r_next = V_row[i];
               vij->c_prev = NULL;
               vij->c_next = V_col[j];
               if (vij->r_next != NULL) vij->r_next->r_prev = vij;
               if (vij->c_next != NULL) vij->c_next->c_prev = vij;
               V_row[i] = V_col[j] = vij;
               R_len[i]++, C_len[j]++;
            }
            else
            {  /* there is no fill-in, because v[i,j] already exists in
                  i-th row; restore the flag, which was reset before */
               flag[j] = 1;
            }
         }
         /* now i-th row has been completely transformed and can return
            to the active set with a new length */
         R_prev[i] = 0;
         R_next[i] = R_head[R_len[i]];
         if (R_next[i] != 0) R_prev[R_next[i]] = i;
         R_head[R_len[i]] = i;
      }
      /* at this point q-th (pivot) column must be empty */
      xassert(C_len[q] == 0);
      /* walk through p-th (pivot) row again and do the following... */
      for (vpj = V_row[p]; vpj != NULL; vpj = vpj->r_next)
      {  /* get column index of v[p,j] */
         j = vpj->j;
         /* erase v[p,j] from the working array */
         flag[j] = 0;
         mpq_set_si(work[j], 0, 1);
         /* now j-th column has been completely transformed, so it can
            return to the active list with a new length */
         C_prev[j] = 0;
         C_next[j] = C_head[C_len[j]];
         if (C_next[j] != 0) C_prev[C_next[j]] = j;
         C_head[C_len[j]] = j;
      }
      mpq_clear(temp);
      /* return to the factorizing routine */
      return;
}

/*----------------------------------------------------------------------
// lux_decomp - compute LU-factorization.
//
// SYNOPSIS
//
// #include "glplux.h"
// int lux_decomp(LUX *lux, int (*col)(void *info, int j, int ind[],
//    mpq_t val[]), void *info);
//
// DESCRIPTION
//
// The routine lux_decomp computes LU-factorization of a given square
// matrix A.
//
// The parameter lux specifies LU-factorization data structure built by
// means of the routine lux_create.
//
// The formal routine col specifies the original matrix A. In order to
// obtain j-th column of the matrix A the routine lux_decomp calls the
// routine col with the parameter j (1 <= j <= n, where n is the order
// of A). In response the routine col should store row indices and
// numerical values of non-zero elements of j-th column of A to the
// locations ind[1], ..., ind[len] and val[1], ..., val[len], resp.,
// where len is the number of non-zeros in j-th column, which should be
// returned on exit. Neiter zero nor duplicate elements are allowed.