glpspx01.c

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

      csa->bfd = lp->bfd, lp->bfd = NULL;
      /* matrix N (by rows) */
      alloc_N(csa);
      build_N(csa);
      /* working parameters */
      csa->phase = 0;
      csa->tm_beg = xtime();
      csa->it_beg = csa->it_cnt = lp->it_cnt;
      csa->it_dpy = -1;
      /* reference space and steepest edge coefficients */
      csa->refct = 0;
      memset(&refsp[1], 0, (m+n) * sizeof(char));
      for (j = 1; j <= n; j++) gamma[j] = 1.0;
      return;
}

/***********************************************************************
*  invert_B - compute factorization of the basis matrix
*
*  This routine computes factorization of the current basis matrix B.
*
*  If the operation is successful, the routine returns zero, otherwise
*  non-zero. */

static int inv_col(void *info, int i, int ind[], double val[])
{     /* this auxiliary routine returns row indices and numeric values
         of non-zero elements of i-th column of the basis matrix */
      struct csa *csa = info;
      int m = csa->m;
#ifdef GLP_DEBUG
      int n = csa->n;
#endif
      int *A_ptr = csa->A_ptr;
      int *A_ind = csa->A_ind;
      double *A_val = csa->A_val;
      int *head = csa->head;
      int k, len, ptr, t;
#ifdef GLP_DEBUG
      xassert(1 <= i && i <= m);
#endif
      k = head[i]; /* B[i] is k-th column of (I|-A) */
#ifdef GLP_DEBUG
      xassert(1 <= k && k <= m+n);
#endif
      if (k <= m)
      {  /* B[i] is k-th column of submatrix I */
         len = 1;
         ind[1] = k;
         val[1] = 1.0;
      }
      else
      {  /* B[i] is (k-m)-th column of submatrix (-A) */
         ptr = A_ptr[k-m];
         len = A_ptr[k-m+1] - ptr;
         memcpy(&ind[1], &A_ind[ptr], len * sizeof(int));
         memcpy(&val[1], &A_val[ptr], len * sizeof(double));
         for (t = 1; t <= len; t++) val[t] = - val[t];
      }
      return len;
}

static int invert_B(struct csa *csa)
{     int ret;
      ret = bfd_factorize(csa->bfd, csa->m, NULL, inv_col, csa);
      csa->valid = (ret == 0);
      return ret;
}

/***********************************************************************
*  update_B - update factorization of the basis matrix
*
*  This routine replaces i-th column of the basis matrix B by k-th
*  column of the augmented constraint matrix (I|-A) and then updates
*  the factorization of B.
*
*  If the factorization has been successfully updated, the routine
*  returns zero, otherwise non-zero. */

static int update_B(struct csa *csa, int i, int k)
{     int m = csa->m;
#ifdef GLP_DEBUG
      int n = csa->n;
#endif
      int ret;
#ifdef GLP_DEBUG
      xassert(1 <= i && i <= m);
      xassert(1 <= k && k <= m+n);
#endif
      if (k <= m)
      {  /* new i-th column of B is k-th column of I */
         int ind[1+1];
         double val[1+1];
         ind[1] = k;
         val[1] = 1.0;
         xassert(csa->valid);
         ret = bfd_update_it(csa->bfd, i, 0, 1, ind, val);
      }
      else
      {  /* new i-th column of B is (k-m)-th column of (-A) */
         int *A_ptr = csa->A_ptr;
         int *A_ind = csa->A_ind;
         double *A_val = csa->A_val;
         double *val = csa->work1;
         int beg, end, ptr, len;
         beg = A_ptr[k-m];
         end = A_ptr[k-m+1];
         len = 0;
         for (ptr = beg; ptr < end; ptr++)
            val[++len] = - A_val[ptr];
         xassert(csa->valid);
         ret = bfd_update_it(csa->bfd, i, 0, len, &A_ind[beg-1], val);
      }
      csa->valid = (ret == 0);
      return ret;
}

/***********************************************************************
*  error_ftran - compute residual vector r = h - B * x
*
*  This routine computes the residual vector r = h - B * x, where B is
*  the current basis matrix, h is the vector of right-hand sides, x is
*  the solution vector. */

static void error_ftran(struct csa *csa, double h[], double x[],
      double r[])
{     int m = csa->m;
#ifdef GLP_DEBUG
      int n = csa->n;
#endif
      int *A_ptr = csa->A_ptr;
      int *A_ind = csa->A_ind;
      double *A_val = csa->A_val;
      int *head = csa->head;
      int i, k, beg, end, ptr;
      double temp;
      /* compute the residual vector:
         r = h - B * x = h - B[1] * x[1] - ... - B[m] * x[m],
         where B[1], ..., B[m] are columns of matrix B */
      memcpy(&r[1], &h[1], m * sizeof(double));
      for (i = 1; i <= m; i++)
      {  temp = x[i];
         if (temp == 0.0) continue;
         k = head[i]; /* B[i] is k-th column of (I|-A) */
#ifdef GLP_DEBUG
         xassert(1 <= k && k <= m+n);
#endif
         if (k <= m)
         {  /* B[i] is k-th column of submatrix I */
            r[k] -= temp;
         }
         else
         {  /* B[i] is (k-m)-th column of submatrix (-A) */
            beg = A_ptr[k-m];
            end = A_ptr[k-m+1];
            for (ptr = beg; ptr < end; ptr++)
               r[A_ind[ptr]] += A_val[ptr] * temp;
         }
      }
      return;
}

/***********************************************************************
*  refine_ftran - refine solution of B * x = h
*
*  This routine performs one iteration to refine the solution of
*  the system B * x = h, where B is the current basis matrix, h is the
*  vector of right-hand sides, x is the solution vector. */

static void refine_ftran(struct csa *csa, double h[], double x[])
{     int m = csa->m;
      double *r = csa->work1;
      double *d = csa->work1;
      int i;
      /* compute the residual vector r = h - B * x */
      error_ftran(csa, h, x, r);
      /* compute the correction vector d = inv(B) * r */
      xassert(csa->valid);
      bfd_ftran(csa->bfd, d);
      /* refine the solution vector (new x) = (old x) + d */
      for (i = 1; i <= m; i++) x[i] += d[i];
      return;
}

/***********************************************************************
*  error_btran - compute residual vector r = h - B'* x
*
*  This routine computes the residual vector r = h - B'* x, where B'
*  is a matrix transposed to the current basis matrix, h is the vector
*  of right-hand sides, x is the solution vector. */

static void error_btran(struct csa *csa, double h[], double x[],
      double r[])
{     int m = csa->m;
#ifdef GLP_DEBUG
      int n = csa->n;
#endif
      int *A_ptr = csa->A_ptr;
      int *A_ind = csa->A_ind;
      double *A_val = csa->A_val;
      int *head = csa->head;
      int i, k, beg, end, ptr;
      double temp;
      /* compute the residual vector r = b - B'* x */
      for (i = 1; i <= m; i++)
      {  /* r[i] := b[i] - (i-th column of B)'* x */
         k = head[i]; /* B[i] is k-th column of (I|-A) */
#ifdef GLP_DEBUG
         xassert(1 <= k && k <= m+n);
#endif
         temp = h[i];
         if (k <= m)
         {  /* B[i] is k-th column of submatrix I */
            temp -= x[k];
         }
         else
         {  /* B[i] is (k-m)-th column of submatrix (-A) */
            beg = A_ptr[k-m];
            end = A_ptr[k-m+1];
            for (ptr = beg; ptr < end; ptr++)
               temp += A_val[ptr] * x[A_ind[ptr]];
         }
         r[i] = temp;
      }
      return;
}

/***********************************************************************
*  refine_btran - refine solution of B'* x = h
*
*  This routine performs one iteration to refine the solution of the
*  system B'* x = h, where B' is a matrix transposed to the current
*  basis matrix, h is the vector of right-hand sides, x is the solution
*  vector. */

static void refine_btran(struct csa *csa, double h[], double x[])
{     int m = csa->m;
      double *r = csa->work1;
      double *d = csa->work1;
      int i;
      /* compute the residual vector r = h - B'* x */
      error_btran(csa, h, x, r);
      /* compute the correction vector d = inv(B') * r */
      xassert(csa->valid);
      bfd_btran(csa->bfd, d);
      /* refine the solution vector (new x) = (old x) + d */
      for (i = 1; i <= m; i++) x[i] += d[i];
      return;
}

/***********************************************************************
*  alloc_N - allocate matrix N
*
*  This routine determines maximal row lengths of matrix N, sets its
*  row pointers, and then allocates arrays N_ind and N_val.
*
*  Note that some fixed structural variables may temporarily become
*  double-bounded, so corresponding columns of matrix A should not be
*  ignored on calculating maximal row lengths of matrix N. */

static void alloc_N(struct csa *csa)
{     int m = csa->m;
      int n = csa->n;
      int *A_ptr = csa->A_ptr;
      int *A_ind = csa->A_ind;
      int *N_ptr = csa->N_ptr;
      int *N_len = csa->N_len;
      int i, j, beg, end, ptr;
      /* determine number of non-zeros in each row of the augmented
         constraint matrix (I|-A) */
      for (i = 1; i <= m; i++)
         N_len[i] = 1;
      for (j = 1; j <= n; j++)
      {  beg = A_ptr[j];
         end = A_ptr[j+1];
         for (ptr = beg; ptr < end; ptr++)
            N_len[A_ind[ptr]]++;
      }
      /* determine maximal row lengths of matrix N and set its row
         pointers */
      N_ptr[1] = 1;
      for (i = 1; i <= m; i++)
      {  /* row of matrix N cannot have more than n non-zeros */
         if (N_len[i] > n) N_len[i] = n;
         N_ptr[i+1] = N_ptr[i] + N_len[i];
      }
      /* now maximal number of non-zeros in matrix N is known */
      csa->N_ind = xcalloc(N_ptr[m+1], sizeof(int));
      csa->N_val = xcalloc(N_ptr[m+1], sizeof(double));
      return;
}

/***********************************************************************
*  add_N_col - add column of matrix (I|-A) to matrix N
*
*  This routine adds j-th column to matrix N which is k-th column of
*  the augmented constraint matrix (I|-A). (It is assumed that old j-th
*  column was previously removed from matrix N.) */

static void add_N_col(struct csa *csa, int j, int k)
{     int m = csa->m;
#ifdef GLP_DEBUG
      int n = csa->n;
#endif
      int *N_ptr = csa->N_ptr;
      int *N_len = csa->N_len;
      int *N_ind = csa->N_ind;
      double *N_val = csa->N_val;
      int pos;
#ifdef GLP_DEBUG
      xassert(1 <= j && j <= n);
      xassert(1 <= k && k <= m+n);
#endif
      if (k <= m)
      {  /* N[j] is k-th column of submatrix I */
         pos = N_ptr[k] + (N_len[k]++);
#ifdef GLP_DEBUG
         xassert(pos < N_ptr[k+1]);
#endif
         N_ind[pos] = j;
         N_val[pos] = 1.0;
      }
      else
      {  /* N[j] is (k-m)-th column of submatrix (-A) */
         int *A_ptr = csa->A_ptr;
         int *A_ind = csa->A_ind;
         double *A_val = csa->A_val;
         int i, beg, end, ptr;
         beg = A_ptr[k-m];
         end = A_ptr[k-m+1];
         for (ptr = beg; ptr < end; ptr++)
         {  i = A_ind[ptr]; /* row number */
            pos = N_ptr[i] + (N_len[i]++);
#ifdef GLP_DEBUG
            xassert(pos < N_ptr[i+1]);
#endif
            N_ind[pos] = j;
            N_val[pos] = - A_val[ptr];
         }
      }
      return;
}

/***********************************************************************
*  del_N_col - remove column of matrix (I|-A) from matrix N
*
*  This routine removes j-th column from matrix N which is k-th column
*  of the augmented constraint matrix (I|-A). */

static void del_N_col(struct csa *csa, int j, int k)
{     int m = csa->m;
#ifdef GLP_DEBUG
      int n = csa->n;
#endif
      int *N_ptr = csa->N_ptr;
      int *N_len = csa->N_len;
      int *N_ind = csa->N_ind;
      double *N_val = csa->N_val;
      int pos, head, tail;
#ifdef GLP_DEBUG
      xassert(1 <= j && j <= n);
      xassert(1 <= k && k <= m+n);
#endif
      if (k <= m)
      {  /* N[j] is k-th column of submatrix I */
         /* find element in k-th row of N */
         head = N_ptr[k];
         for (pos = head; N_ind[pos] != j; pos++) /* nop */;
         /* and remove it from the row list */
         tail = head + (--N_len[k]);
#ifdef GLP_DEBUG
         xassert(pos <= tail);
#endif
         N_ind[pos] = N_ind[tail];
         N_val[pos] = N_val[tail];
      }
      else
      {  /* N[j] is (k-m)-th column of submatrix (-A) */
         int *A_ptr = csa->A_ptr;
         int *A_ind = csa->A_ind;
         int i, beg, end, ptr;
         beg = A_ptr[k-m];
         end = A_ptr[k-m+1];
         for (ptr = beg; ptr < end; ptr++)
         {  i = A_ind[ptr]; /* row number */