glplux.c

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

//
// The parameter info is a transit pointer passed to the formal routine
// col; it can be used for various purposes.
//
// RETURNS
//
// The routine lux_decomp returns the singularity flag. Zero flag means
// that the original matrix A is non-singular while non-zero flag means
// that A is (exactly!) singular.
//
// Note that LU-factorization is valid in both cases, however, in case
// of singularity some rows of the matrix V (including pivot elements)
// will be empty.
//
// REPAIRING SINGULAR MATRIX
//
// If the routine lux_decomp returns non-zero flag, it provides all
// necessary information that can be used for "repairing" the matrix A,
// where "repairing" means replacing linearly dependent columns of the
// matrix A by appropriate columns of the unity matrix. This feature is
// needed when the routine lux_decomp is used for reinverting the basis
// matrix within the simplex method procedure.
//
// On exit linearly dependent columns of the matrix U have the numbers
// rank+1, rank+2, ..., n, where rank is the exact rank of the matrix A
// stored by the routine to the member lux->rank. The correspondence
// between columns of A and U is the same as between columns of V and U.
// Thus, linearly dependent columns of the matrix A have the numbers
// Q_col[rank+1], Q_col[rank+2], ..., Q_col[n], where Q_col is an array
// representing the permutation matrix Q in column-like format. It is
// understood that each j-th linearly dependent column of the matrix U
// should be replaced by the unity vector, where all elements are zero
// except the unity diagonal element u[j,j]. On the other hand j-th row
// of the matrix U corresponds to the row of the matrix V (and therefore
// of the matrix A) with the number P_row[j], where P_row is an array
// representing the permutation matrix P in row-like format. Thus, each
// j-th linearly dependent column of the matrix U should be replaced by
// a column of the unity matrix with the number P_row[j].
//
// The code that repairs the matrix A may look like follows:
//
//    for (j = rank+1; j <= n; j++)
//    {  replace column Q_col[j] of the matrix A by column P_row[j] of
//       the unity matrix;
//    }
//
// where rank, P_row, and Q_col are members of the structure LUX. */

int lux_decomp(LUX *lux, int (*col)(void *info, int j, int ind[],
      mpq_t val[]), void *info)
{     int n = lux->n;
      LUXELM **V_row = lux->V_row;
      LUXELM **V_col = lux->V_col;
      int *P_row = lux->P_row;
      int *P_col = lux->P_col;
      int *Q_row = lux->Q_row;
      int *Q_col = lux->Q_col;
      LUXELM *piv, *vij;
      LUXWKA *wka;
      int i, j, k, p, q, t, *flag;
      mpq_t *work;
      /* allocate working area */
      wka = xmalloc(sizeof(LUXWKA));
      wka->R_len = xcalloc(1+n, sizeof(int));
      wka->R_head = xcalloc(1+n, sizeof(int));
      wka->R_prev = xcalloc(1+n, sizeof(int));
      wka->R_next = xcalloc(1+n, sizeof(int));
      wka->C_len = xcalloc(1+n, sizeof(int));
      wka->C_head = xcalloc(1+n, sizeof(int));
      wka->C_prev = xcalloc(1+n, sizeof(int));
      wka->C_next = xcalloc(1+n, sizeof(int));
      /* initialize LU-factorization data structures */
      initialize(lux, col, info, wka);
      /* allocate working arrays */
      flag = xcalloc(1+n, sizeof(int));
      work = xcalloc(1+n, sizeof(mpq_t));
      for (k = 1; k <= n; k++)
      {  flag[k] = 0;
         mpq_init(work[k]);
      }
      /* main elimination loop */
      for (k = 1; k <= n; k++)
      {  /* choose a pivot element v[p,q] */
         piv = find_pivot(lux, wka);
         if (piv == NULL)
         {  /* no pivot can be chosen, because the active submatrix is
               empty */
            break;
         }
         /* determine row and column indices of the pivot element */
         p = piv->i, q = piv->j;
         /* let v[p,q] correspond to u[i',j']; permute k-th and i'-th
            rows and k-th and j'-th columns of the matrix U = P*V*Q to
            move the element u[i',j'] to the position u[k,k] */
         i = P_col[p], j = Q_row[q];
         xassert(k <= i && i <= n && k <= j && j <= n);
         /* permute k-th and i-th rows of the matrix U */
         t = P_row[k];
         P_row[i] = t, P_col[t] = i;
         P_row[k] = p, P_col[p] = k;
         /* permute k-th and j-th columns of the matrix U */
         t = Q_col[k];
         Q_col[j] = t, Q_row[t] = j;
         Q_col[k] = q, Q_row[q] = k;
         /* eliminate subdiagonal elements of k-th column of the matrix
            U = P*V*Q using the pivot element u[k,k] = v[p,q] */
         eliminate(lux, wka, piv, flag, work);
      }
      /* determine the rank of A (and V) */
      lux->rank = k - 1;
      /* free working arrays */
      xfree(flag);
      for (k = 1; k <= n; k++) mpq_clear(work[k]);
      xfree(work);
      /* build column lists of the matrix V using its row lists */
      for (j = 1; j <= n; j++)
         xassert(V_col[j] == NULL);
      for (i = 1; i <= n; i++)
      {  for (vij = V_row[i]; vij != NULL; vij = vij->r_next)
         {  j = vij->j;
            vij->c_prev = NULL;
            vij->c_next = V_col[j];
            if (vij->c_next != NULL) vij->c_next->c_prev = vij;
            V_col[j] = vij;
         }
      }
      /* free working area */
      xfree(wka->R_len);
      xfree(wka->R_head);
      xfree(wka->R_prev);
      xfree(wka->R_next);
      xfree(wka->C_len);
      xfree(wka->C_head);
      xfree(wka->C_prev);
      xfree(wka->C_next);
      xfree(wka);
      /* return to the calling program */
      return (lux->rank < n);
}

/*----------------------------------------------------------------------
// lux_f_solve - solve system F*x = b or F'*x = b.
//
// SYNOPSIS
//
// #include "glplux.h"
// void lux_f_solve(LUX *lux, int tr, mpq_t x[]);
//
// DESCRIPTION
//
// The routine lux_f_solve solves either the system F*x = b (if the
// flag tr is zero) or the system F'*x = b (if the flag tr is non-zero),
// where the matrix F is a component of LU-factorization specified by
// the parameter lux, F' is a matrix transposed to F.
//
// On entry the array x should contain elements of the right-hand side
// vector b in locations x[1], ..., x[n], where n is the order of the
// matrix F. On exit this array will contain elements of the solution
// vector x in the same locations. */

void lux_f_solve(LUX *lux, int tr, mpq_t x[])
{     int n = lux->n;
      LUXELM **F_row = lux->F_row;
      LUXELM **F_col = lux->F_col;
      int *P_row = lux->P_row;
      LUXELM *fik, *fkj;
      int i, j, k;
      mpq_t temp;
      mpq_init(temp);
      if (!tr)
      {  /* solve the system F*x = b */
         for (j = 1; j <= n; j++)
         {  k = P_row[j];
            if (mpq_sgn(x[k]) != 0)
            {  for (fik = F_col[k]; fik != NULL; fik = fik->c_next)
               {  mpq_mul(temp, fik->val, x[k]);
                  mpq_sub(x[fik->i], x[fik->i], temp);
               }
            }
         }
      }
      else
      {  /* solve the system F'*x = b */
         for (i = n; i >= 1; i--)
         {  k = P_row[i];
            if (mpq_sgn(x[k]) != 0)
            {  for (fkj = F_row[k]; fkj != NULL; fkj = fkj->r_next)
               {  mpq_mul(temp, fkj->val, x[k]);
                  mpq_sub(x[fkj->j], x[fkj->j], temp);
               }
            }
         }
      }
      mpq_clear(temp);
      return;
}

/*----------------------------------------------------------------------
// lux_v_solve - solve system V*x = b or V'*x = b.
//
// SYNOPSIS
//
// #include "glplux.h"
// void lux_v_solve(LUX *lux, int tr, double x[]);
//
// DESCRIPTION
//
// The routine lux_v_solve solves either the system V*x = b (if the
// flag tr is zero) or the system V'*x = b (if the flag tr is non-zero),
// where the matrix V is a component of LU-factorization specified by
// the parameter lux, V' is a matrix transposed to V.
//
// On entry the array x should contain elements of the right-hand side
// vector b in locations x[1], ..., x[n], where n is the order of the
// matrix V. On exit this array will contain elements of the solution
// vector x in the same locations. */

void lux_v_solve(LUX *lux, int tr, mpq_t x[])
{     int n = lux->n;
      mpq_t *V_piv = lux->V_piv;
      LUXELM **V_row = lux->V_row;
      LUXELM **V_col = lux->V_col;
      int *P_row = lux->P_row;
      int *Q_col = lux->Q_col;
      LUXELM *vij;
      int i, j, k;
      mpq_t *b, temp;
      b = xcalloc(1+n, sizeof(mpq_t));
      for (k = 1; k <= n; k++)
         mpq_init(b[k]), mpq_set(b[k], x[k]), mpq_set_si(x[k], 0, 1);
      mpq_init(temp);
      if (!tr)
      {  /* solve the system V*x = b */
         for (k = n; k >= 1; k--)
         {  i = P_row[k], j = Q_col[k];
            if (mpq_sgn(b[i]) != 0)
            {  mpq_set(x[j], b[i]);
               mpq_div(x[j], x[j], V_piv[i]);
               for (vij = V_col[j]; vij != NULL; vij = vij->c_next)
               {  mpq_mul(temp, vij->val, x[j]);
                  mpq_sub(b[vij->i], b[vij->i], temp);
               }
            }
         }
      }
      else
      {  /* solve the system V'*x = b */
         for (k = 1; k <= n; k++)
         {  i = P_row[k], j = Q_col[k];
            if (mpq_sgn(b[j]) != 0)
            {  mpq_set(x[i], b[j]);
               mpq_div(x[i], x[i], V_piv[i]);
               for (vij = V_row[i]; vij != NULL; vij = vij->r_next)
               {  mpq_mul(temp, vij->val, x[i]);
                  mpq_sub(b[vij->j], b[vij->j], temp);
               }
            }
         }
      }
      for (k = 1; k <= n; k++) mpq_clear(b[k]);
      mpq_clear(temp);
      xfree(b);
      return;
}

/*----------------------------------------------------------------------
// lux_solve - solve system A*x = b or A'*x = b.
//
// SYNOPSIS
//
// #include "glplux.h"
// void lux_solve(LUX *lux, int tr, mpq_t x[]);
//
// DESCRIPTION
//
// The routine lux_solve solves either the system A*x = b (if the flag
// tr is zero) or the system A'*x = b (if the flag tr is non-zero),
// where the parameter lux specifies LU-factorization of the matrix A,
// A' is a matrix transposed to A.
//
// On entry the array x should contain elements of the right-hand side
// vector b in locations x[1], ..., x[n], where n is the order of the
// matrix A. On exit this array will contain elements of the solution
// vector x in the same locations. */

void lux_solve(LUX *lux, int tr, mpq_t x[])
{     if (lux->rank < lux->n)
         xfault("lux_solve: LU-factorization has incomplete rank\n");
      if (!tr)
      {  /* A = F*V, therefore inv(A) = inv(V)*inv(F) */
         lux_f_solve(lux, 0, x);
         lux_v_solve(lux, 0, x);
      }
      else
      {  /* A' = V'*F', therefore inv(A') = inv(F')*inv(V') */
         lux_v_solve(lux, 1, x);
         lux_f_solve(lux, 1, x);
      }
      return;
}

/*----------------------------------------------------------------------
// lux_delete - delete LU-factorization.
//
// SYNOPSIS
//
// #include "glplux.h"
// void lux_delete(LUX *lux);
//
// DESCRIPTION
//
// The routine lux_delete deletes LU-factorization data structure,
// which the parameter lux points to, freeing all the memory allocated
// to this object. */

void lux_delete(LUX *lux)
{     int n = lux->n;
      LUXELM *fij, *vij;
      int i;
      for (i = 1; i <= n; i++)
      {  for (fij = lux->F_row[i]; fij != NULL; fij = fij->r_next)
            mpq_clear(fij->val);
         mpq_clear(lux->V_piv[i]);
         for (vij = lux->V_row[i]; vij != NULL; vij = vij->r_next)
            mpq_clear(vij->val);
      }
      dmp_delete_pool(lux->pool);
      xfree(lux->F_row);
      xfree(lux->F_col);
      xfree(lux->V_piv);
      xfree(lux->V_row);
      xfree(lux->V_col);
      xfree(lux->P_row);
      xfree(lux->P_col);
      xfree(lux->Q_row);
      xfree(lux->Q_col);
      xfree(lux);
      return;
}

/* eof */