glpluf.c

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

      /* store value of the pivot */
      vpq = (vr_piv[p] = sv_val[p_ptr]);
      /* remove the pivot from the p-th row */
      sv_ind[p_ptr] = sv_ind[p_end];
      sv_val[p_ptr] = sv_val[p_end];
      vr_len[p]--;
      p_end--;
      /* find the pivot v[p,q] = u[k,k] in the q-th column */
      q_beg = vc_ptr[q];
      q_end = q_beg + vc_len[q] - 1;
      for (q_ptr = q_beg; sv_ind[q_ptr] != p; q_ptr++) /* nop */;
      xassert(q_ptr <= q_end);
      /* remove the pivot from the q-th column */
      sv_ind[q_ptr] = sv_ind[q_end];
      vc_len[q]--;
      q_end--;
      /* walk through the p-th (pivot) row, which doesn't contain the
         pivot v[p,q] already, and do the following... */
      for (p_ptr = p_beg; p_ptr <= p_end; p_ptr++)
      {  /* get column index of v[p,j] */
         j = sv_ind[p_ptr];
         /* store v[p,j] to the working array */
         flag[j] = 1;
         work[j] = sv_val[p_ptr];
         /* remove the j-th column from the active set; this column will
            return there later with new length */
         if (cs_prev[j] == 0)
            cs_head[vc_len[j]] = cs_next[j];
         else
            cs_next[cs_prev[j]] = cs_next[j];
         if (cs_next[j] == 0)
            ;
         else
            cs_prev[cs_next[j]] = cs_prev[j];
         /* find v[p,j] in the j-th column */
         j_beg = vc_ptr[j];
         j_end = j_beg + vc_len[j] - 1;
         for (j_ptr = j_beg; sv_ind[j_ptr] != p; j_ptr++) /* nop */;
         xassert(j_ptr <= j_end);
         /* since v[p,j] leaves the active submatrix, remove it from the
            j-th column; however, v[p,j] is kept in the p-th row */
         sv_ind[j_ptr] = sv_ind[j_end];
         vc_len[j]--;
      }
      /* walk through the q-th (pivot) column, which doesn't contain the
         pivot v[p,q] already, and perform gaussian elimination */
      while (q_beg <= q_end)
      {  /* element v[i,q] should be eliminated */
         /* get row index of v[i,q] */
         i = sv_ind[q_beg];
         /* remove the i-th row from the active set; later this row will
            return there with new length */
         if (rs_prev[i] == 0)
            rs_head[vr_len[i]] = rs_next[i];
         else
            rs_next[rs_prev[i]] = rs_next[i];
         if (rs_next[i] == 0)
            ;
         else
            rs_prev[rs_next[i]] = rs_prev[i];
         /* find v[i,q] in the i-th row */
         i_beg = vr_ptr[i];
         i_end = i_beg + vr_len[i] - 1;
         for (i_ptr = i_beg; sv_ind[i_ptr] != q; i_ptr++) /* nop */;
         xassert(i_ptr <= i_end);
         /* compute gaussian multiplier f[i,p] = v[i,q] / v[p,q] */
         fip = sv_val[i_ptr] / vpq;
         /* since v[i,q] should be eliminated, remove it from the i-th
            row */
         sv_ind[i_ptr] = sv_ind[i_end];
         sv_val[i_ptr] = sv_val[i_end];
         vr_len[i]--;
         i_end--;
         /* and from the q-th column */
         sv_ind[q_beg] = sv_ind[q_end];
         vc_len[q]--;
         q_end--;
         /* perform gaussian transformation:
            (i-th row) := (i-th row) - f[i,p] * (p-th row)
            note that now the p-th row, which is in the working array,
            doesn't contain the pivot v[p,q], and the i-th row doesn't
            contain the eliminated element v[i,q] */
         /* walk through the i-th row and transform existing non-zero
            elements */
         fill = vr_len[p];
         for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++)
         {  /* get column index of v[i,j] */
            j = sv_ind[i_ptr];
            /* v[i,j] := v[i,j] - f[i,p] * v[p,j] */
            if (flag[j])
            {  /* v[p,j] != 0 */
               temp = (sv_val[i_ptr] -= fip * work[j]);
               if (temp < 0.0) temp = - temp;
               flag[j] = 0;
               fill--; /* since both v[i,j] and v[p,j] exist */
               if (temp == 0.0 || temp < eps_tol)
               {  /* new v[i,j] is closer to zero; replace it by exact
                     zero, i.e. remove it from the active submatrix */
                  /* remove v[i,j] from the i-th row */
                  sv_ind[i_ptr] = sv_ind[i_end];
                  sv_val[i_ptr] = sv_val[i_end];
                  vr_len[i]--;
                  i_ptr--;
                  i_end--;
                  /* find v[i,j] in the j-th column */
                  j_beg = vc_ptr[j];
                  j_end = j_beg + vc_len[j] - 1;
                  for (j_ptr = j_beg; sv_ind[j_ptr] != i; j_ptr++);
                  xassert(j_ptr <= j_end);
                  /* remove v[i,j] from the j-th column */
                  sv_ind[j_ptr] = sv_ind[j_end];
                  vc_len[j]--;
               }
               else
               {  /* v_big := max(v_big, |v[i,j]|) */
                  if (luf->big_v < temp) luf->big_v = temp;
               }
            }
         }
         /* now flag is the pattern of the set v[p,*] \ v[i,*], and fill
            is number of non-zeros in this set; therefore up to fill new
            non-zeros may appear in the i-th row */
         if (vr_len[i] + fill > vr_cap[i])
         {  /* enlarge the i-th row */
            if (luf_enlarge_row(luf, i, vr_len[i] + fill))
            {  /* overflow of the sparse vector area */
               ret = 1;
               goto done;
            }
            /* defragmentation may change row and column pointers of the
               matrix V */
            p_beg = vr_ptr[p];
            p_end = p_beg + vr_len[p] - 1;
            q_beg = vc_ptr[q];
            q_end = q_beg + vc_len[q] - 1;
         }
         /* walk through the p-th (pivot) row and create new elements
            of the i-th row that appear due to fill-in; column indices
            of these new elements are accumulated in the array ndx */
         len = 0;
         for (p_ptr = p_beg; p_ptr <= p_end; p_ptr++)
         {  /* get column index of v[p,j], which may cause fill-in */
            j = sv_ind[p_ptr];
            if (flag[j])
            {  /* compute new non-zero v[i,j] = 0 - f[i,p] * v[p,j] */
               temp = (val = - fip * work[j]);
               if (temp < 0.0) temp = - temp;
               if (temp == 0.0 || temp < eps_tol)
                  /* if v[i,j] is closer to zero; just ignore it */;
               else
               {  /* add v[i,j] to the i-th row */
                  i_ptr = vr_ptr[i] + vr_len[i];
                  sv_ind[i_ptr] = j;
                  sv_val[i_ptr] = val;
                  vr_len[i]++;
                  /* remember column index of v[i,j] */
                  ndx[++len] = j;
                  /* big_v := max(big_v, |v[i,j]|) */
                  if (luf->big_v < temp) luf->big_v = temp;
               }
            }
            else
            {  /* there is no fill-in, because v[i,j] already exists in
                  the i-th row; restore the flag of the element v[p,j],
                  which was reset before */
               flag[j] = 1;
            }
         }
         /* add new non-zeros v[i,j] to the corresponding columns */
         for (k = 1; k <= len; k++)
         {  /* get column index of new non-zero v[i,j] */
            j = ndx[k];
            /* one free location is needed in the j-th column */
            if (vc_len[j] + 1 > vc_cap[j])
            {  /* enlarge the j-th column */
               if (luf_enlarge_col(luf, j, vc_len[j] + 10))
               {  /* overflow of the sparse vector area */
                  ret = 1;
                  goto done;
               }
               /* defragmentation may change row and column pointers of
                  the matrix V */
               p_beg = vr_ptr[p];
               p_end = p_beg + vr_len[p] - 1;
               q_beg = vc_ptr[q];
               q_end = q_beg + vc_len[q] - 1;
            }
            /* add new non-zero v[i,j] to the j-th column */
            j_ptr = vc_ptr[j] + vc_len[j];
            sv_ind[j_ptr] = i;
            vc_len[j]++;
         }
         /* now the i-th row has been completely transformed, therefore
            it can return to the active set with new length */
         rs_prev[i] = 0;
         rs_next[i] = rs_head[vr_len[i]];
         if (rs_next[i] != 0) rs_prev[rs_next[i]] = i;
         rs_head[vr_len[i]] = i;
         /* the largest of absolute values of elements in the i-th row
            is currently unknown */
         vr_max[i] = -1.0;
         /* at least one free location is needed to store the gaussian
            multiplier */
         if (luf->sv_end - luf->sv_beg < 1)
         {  /* there are no free locations at all; defragment SVA */
            luf_defrag_sva(luf);
            if (luf->sv_end - luf->sv_beg < 1)
            {  /* overflow of the sparse vector area */
               ret = 1;
               goto done;
            }
            /* defragmentation may change row and column pointers of the
               matrix V */
            p_beg = vr_ptr[p];
            p_end = p_beg + vr_len[p] - 1;
            q_beg = vc_ptr[q];
            q_end = q_beg + vc_len[q] - 1;
         }
         /* add the element f[i,p], which is the gaussian multiplier,
            to the matrix F */
         luf->sv_end--;
         sv_ind[luf->sv_end] = i;
         sv_val[luf->sv_end] = fip;
         fc_len[p]++;
         /* end of elimination loop */
      }
      /* at this point the q-th (pivot) column should be empty */
      xassert(vc_len[q] == 0);
      /* reset capacity of the q-th column */
      vc_cap[q] = 0;
      /* remove node of the q-th column from the addressing list */
      k = n + q;
      if (sv_prev[k] == 0)
         luf->sv_head = sv_next[k];
      else
         sv_next[sv_prev[k]] = sv_next[k];
      if (sv_next[k] == 0)
         luf->sv_tail = sv_prev[k];
      else
         sv_prev[sv_next[k]] = sv_prev[k];
      /* the p-th column of the matrix F has been completely built; set
         its pointer */
      fc_ptr[p] = luf->sv_end;
      /* walk through the p-th (pivot) row and do the following... */
      for (p_ptr = p_beg; p_ptr <= p_end; p_ptr++)
      {  /* get column index of v[p,j] */
         j = sv_ind[p_ptr];
         /* erase v[p,j] from the working array */
         flag[j] = 0;
         work[j] = 0.0;
         /* the j-th column has been completely transformed, therefore
            it can return to the active set with new length; however
            the special case c_prev[j] = c_next[j] = j means that the
            routine find_pivot excluded the j-th column from the active
            set due to Uwe Suhl's rule, and therefore in this case the
            column can return to the active set only if it is a column
            singleton */
         if (!(vc_len[j] != 1 && cs_prev[j] == j && cs_next[j] == j))
         {  cs_prev[j] = 0;
            cs_next[j] = cs_head[vc_len[j]];
            if (cs_next[j] != 0) cs_prev[cs_next[j]] = j;
            cs_head[vc_len[j]] = j;
         }
      }
done: /* return to the factorizing routine */
      return ret;
}

/***********************************************************************
*  build_v_cols - build the matrix V in column-wise format
*
*  This routine builds the column-wise representation of the matrix V
*  using its row-wise representation.
*
*  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 build_v_cols(LUF *luf)
{     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 *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;
      int ret = 0;
      int i, i_beg, i_end, i_ptr, j, j_ptr, k, nnz;
      /* it is assumed that on entry all columns of the matrix V are
         empty, i.e. vc_len[j] = vc_cap[j] = 0 for all j = 1, ..., n,
         and have been removed from the addressing list */
      /* count non-zeros in columns of the matrix V; count total number
         of non-zeros in this matrix */
      nnz = 0;
      for (i = 1; i <= n; i++)
      {  /* walk through elements of the i-th row and count non-zeros
            in the corresponding columns */
         i_beg = vr_ptr[i];
         i_end = i_beg + vr_len[i] - 1;
         for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++)
            vc_cap[sv_ind[i_ptr]]++;
         /* count total number of non-zeros */
         nnz += vr_len[i];
      }
      /* store total number of non-zeros */
      luf->nnz_v = nnz;
      /* check for free locations */
      if (luf->sv_end - luf->sv_beg < nnz)
      {  /* overflow of the sparse vector area */
         ret = 1;
         goto done;
      }
      /* allocate columns of the matrix V */
      for (j = 1; j <= n; j++)
      {  /* set pointer to the j-th column */
         vc_ptr[j] = luf->sv_beg;
         /* reserve locations for the j-th column */
         luf->sv_beg += vc_cap[j];
      }
      /* build the matrix V in column-wise format using this matrix in
         row-wise format */
      for (i = 1; i <= n; i++)
      {  /* walk through elements of the i-th row */
         i_beg = vr_ptr[i];
         i_end = i_beg + vr_len[i] - 1;
         for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++)
         {  /* get column index */
            j = sv_ind[i_ptr];
            /* store element in the j-th column */
            j_ptr = vc_ptr[j] + vc_len[j];
            sv_ind[j_ptr] = i;