glpluf.c

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

            expense of old locations of the j-th column */
         kk = sv_prev[k];
         if (kk <= n) vr_cap[kk] += cur; else vc_cap[kk-n] += cur;
         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];
      /* insert the j-th column node to the end of the linked list */
      sv_prev[k] = luf->sv_tail;
      sv_next[k] = 0;
      if (sv_prev[k] == 0)
         luf->sv_head = k;
      else
         sv_next[sv_prev[k]] = k;
      luf->sv_tail = k;
done: return ret;
}

/***********************************************************************
*  reallocate - reallocate LU-factorization arrays
*
*  This routine reallocates arrays, whose size depends of n, the order
*  of the matrix A to be factorized. */

static void reallocate(LUF *luf, int n)
{     int n_max = luf->n_max;
      luf->n = n;
      if (n <= n_max) goto done;
      if (luf->fr_ptr != NULL) xfree(luf->fr_ptr);
      if (luf->fr_len != NULL) xfree(luf->fr_len);
      if (luf->fc_ptr != NULL) xfree(luf->fc_ptr);
      if (luf->fc_len != NULL) xfree(luf->fc_len);
      if (luf->vr_ptr != NULL) xfree(luf->vr_ptr);
      if (luf->vr_len != NULL) xfree(luf->vr_len);
      if (luf->vr_cap != NULL) xfree(luf->vr_cap);
      if (luf->vr_piv != NULL) xfree(luf->vr_piv);
      if (luf->vc_ptr != NULL) xfree(luf->vc_ptr);
      if (luf->vc_len != NULL) xfree(luf->vc_len);
      if (luf->vc_cap != NULL) xfree(luf->vc_cap);
      if (luf->pp_row != NULL) xfree(luf->pp_row);
      if (luf->pp_col != NULL) xfree(luf->pp_col);
      if (luf->qq_row != NULL) xfree(luf->qq_row);
      if (luf->qq_col != NULL) xfree(luf->qq_col);
      if (luf->sv_prev != NULL) xfree(luf->sv_prev);
      if (luf->sv_next != NULL) xfree(luf->sv_next);
      if (luf->vr_max != NULL) xfree(luf->vr_max);
      if (luf->rs_head != NULL) xfree(luf->rs_head);
      if (luf->rs_prev != NULL) xfree(luf->rs_prev);
      if (luf->rs_next != NULL) xfree(luf->rs_next);
      if (luf->cs_head != NULL) xfree(luf->cs_head);
      if (luf->cs_prev != NULL) xfree(luf->cs_prev);
      if (luf->cs_next != NULL) xfree(luf->cs_next);
      if (luf->flag != NULL) xfree(luf->flag);
      if (luf->work != NULL) xfree(luf->work);
      luf->n_max = n_max = n + 100;
      luf->fr_ptr = xcalloc(1+n_max, sizeof(int));
      luf->fr_len = xcalloc(1+n_max, sizeof(int));
      luf->fc_ptr = xcalloc(1+n_max, sizeof(int));
      luf->fc_len = xcalloc(1+n_max, sizeof(int));
      luf->vr_ptr = xcalloc(1+n_max, sizeof(int));
      luf->vr_len = xcalloc(1+n_max, sizeof(int));
      luf->vr_cap = xcalloc(1+n_max, sizeof(int));
      luf->vr_piv = xcalloc(1+n_max, sizeof(double));
      luf->vc_ptr = xcalloc(1+n_max, sizeof(int));
      luf->vc_len = xcalloc(1+n_max, sizeof(int));
      luf->vc_cap = xcalloc(1+n_max, sizeof(int));
      luf->pp_row = xcalloc(1+n_max, sizeof(int));
      luf->pp_col = xcalloc(1+n_max, sizeof(int));
      luf->qq_row = xcalloc(1+n_max, sizeof(int));
      luf->qq_col = xcalloc(1+n_max, sizeof(int));
      luf->sv_prev = xcalloc(1+n_max+n_max, sizeof(int));
      luf->sv_next = xcalloc(1+n_max+n_max, sizeof(int));
      luf->vr_max = xcalloc(1+n_max, sizeof(double));
      luf->rs_head = xcalloc(1+n_max, sizeof(int));
      luf->rs_prev = xcalloc(1+n_max, sizeof(int));
      luf->rs_next = xcalloc(1+n_max, sizeof(int));
      luf->cs_head = xcalloc(1+n_max, sizeof(int));
      luf->cs_prev = xcalloc(1+n_max, sizeof(int));
      luf->cs_next = xcalloc(1+n_max, sizeof(int));
      luf->flag = xcalloc(1+n_max, sizeof(int));
      luf->work = xcalloc(1+n_max, sizeof(double));
done: return;
}

/***********************************************************************
*  initialize - initialize LU-factorization data structures
*
*  This routine initializes data structures for subsequent computing
*  the LU-factorization of a given matrix A, which is specified by the
*  formal routine col. On exit V = A and F = P = Q = I, where I is the
*  unity matrix. (Row-wise representation of the matrix F is not used
*  at the factorization stage and therefore is not initialized.)
*
*  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 initialize(LUF *luf, int (*col)(void *info, int j, int rn[],
      double aj[]), void *info)
{     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;
      int *vc_ptr = luf->vc_ptr;
      int *vc_len = luf->vc_len;
      int *vc_cap = luf->vc_cap;
      int *pp_row = luf->pp_row;
      int *pp_col = luf->pp_col;
      int *qq_row = luf->qq_row;
      int *qq_col = luf->qq_col;
      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;
      int ret = 0;
      int i, i_ptr, j, j_beg, j_end, k, len, nnz, sv_beg, sv_end, ptr;
      double big, val;
      /* free all locations of the sparse vector area */
      sv_beg = 1;
      sv_end = luf->sv_size + 1;
      /* (row-wise representation of the matrix F is not initialized,
         because it is not used at the factorization stage) */
      /* build the matrix F in column-wise format (initially F = I) */
      for (j = 1; j <= n; j++)
      {  fc_ptr[j] = sv_end;
         fc_len[j] = 0;
      }
      /* clear rows of the matrix V; clear the flag array */
      for (i = 1; i <= n; i++)
         vr_len[i] = vr_cap[i] = 0, flag[i] = 0;
      /* build the matrix V in column-wise format (initially V = A);
         count non-zeros in rows of this matrix; count total number of
         non-zeros; compute largest of absolute values of elements */
      nnz = 0;
      big = 0.0;
      for (j = 1; j <= n; j++)
      {  int *rn = pp_row;
         double *aj = work;
         /* obtain j-th column of the matrix A */
         len = col(info, j, rn, aj);
         if (!(0 <= len && len <= n))
            xfault("luf_factorize: j = %d; len = %d; invalid column len"
               "gth\n", j, len);
         /* check for free locations */
         if (sv_end - sv_beg < len)
         {  /* overflow of the sparse vector area */
            ret = 1;
            goto done;
         }
         /* set pointer to the j-th column */
         vc_ptr[j] = sv_beg;
         /* set length of the j-th column */
         vc_len[j] = vc_cap[j] = len;
         /* count total number of non-zeros */
         nnz += len;
         /* walk through elements of the j-th column */
         for (ptr = 1; ptr <= len; ptr++)
         {  /* get row index and numerical value of a[i,j] */
            i = rn[ptr];
            val = aj[ptr];
            if (!(1 <= i && i <= n))
               xfault("luf_factorize: i = %d; j = %d; invalid row index"
                  "\n", i, j);
            if (flag[i])
               xfault("luf_factorize: i = %d; j = %d; duplicate element"
                  " not allowed\n", i, j);
            if (val == 0.0)
               xfault("luf_factorize: i = %d; j = %d; zero element not "
                  "allowed\n", i, j);
            /* add new element v[i,j] = a[i,j] to j-th column */
            sv_ind[sv_beg] = i;
            sv_val[sv_beg] = val;
            sv_beg++;
            /* big := max(big, |a[i,j]|) */
            if (val < 0.0) val = - val;
            if (big < val) big = val;
            /* mark non-zero in the i-th position of the j-th column */
            flag[i] = 1;
            /* increase length of the i-th row */
            vr_cap[i]++;
         }
         /* reset all non-zero marks */
         for (ptr = 1; ptr <= len; ptr++) flag[rn[ptr]] = 0;
      }
      /* allocate rows of the matrix V */
      for (i = 1; i <= n; i++)
      {  /* get length of the i-th row */
         len = vr_cap[i];
         /* check for free locations */
         if (sv_end - sv_beg < len)
         {  /* overflow of the sparse vector area */
            ret = 1;
            goto done;
         }
         /* set pointer to the i-th row */
         vr_ptr[i] = sv_beg;
         /* reserve locations for the i-th row */
         sv_beg += len;
      }
      /* build the matrix V in row-wise format using representation of
         this matrix in column-wise format */
      for (j = 1; j <= n; j++)
      {  /* walk through elements of the j-th column */
         j_beg = vc_ptr[j];
         j_end = j_beg + vc_len[j] - 1;
         for (k = j_beg; k <= j_end; k++)
         {  /* get row index and numerical value of v[i,j] */
            i = sv_ind[k];
            val = sv_val[k];
            /* store element in the i-th row */
            i_ptr = vr_ptr[i] + vr_len[i];
            sv_ind[i_ptr] = j;
            sv_val[i_ptr] = val;
            /* increase count of the i-th row */
            vr_len[i]++;
         }
      }
      /* initialize the matrices P and Q (initially P = Q = I) */
      for (k = 1; k <= n; k++)
         pp_row[k] = pp_col[k] = qq_row[k] = qq_col[k] = k;
      /* set sva partitioning pointers */
      luf->sv_beg = sv_beg;
      luf->sv_end = sv_end;
      /* the initial physical order of rows and columns of the matrix V
         is n+1, ..., n+n, 1, ..., n (firstly columns, then rows) */
      luf->sv_head = n+1;
      luf->sv_tail = n;
      for (i = 1; i <= n; i++)
      {  sv_prev[i] = i-1;
         sv_next[i] = i+1;
      }
      sv_prev[1] = n+n;
      sv_next[n] = 0;
      for (j = 1; j <= n; j++)
      {  sv_prev[n+j] = n+j-1;
         sv_next[n+j] = n+j+1;
      }
      sv_prev[n+1] = 0;
      sv_next[n+n] = 1;
      /* clear working arrays */
      for (k = 1; k <= n; k++)
      {  flag[k] = 0;
         work[k] = 0.0;
      }
      /* initialize some statistics */
      luf->nnz_a = nnz;
      luf->nnz_f = 0;
      luf->nnz_v = nnz;
      luf->max_a = big;
      luf->big_v = big;
      luf->rank = -1;
      /* initially the active submatrix is the entire matrix V */
      /* largest of absolute values of elements in each active row is
         unknown yet */
      for (i = 1; i <= n; i++) vr_max[i] = -1.0;
      /* build linked lists of active rows */
      for (len = 0; len <= n; len++) rs_head[len] = 0;
      for (i = 1; i <= n; i++)
      {  len = vr_len[i];
         rs_prev[i] = 0;
         rs_next[i] = rs_head[len];
         if (rs_next[i] != 0) rs_prev[rs_next[i]] = i;
         rs_head[len] = i;
      }
      /* build linked lists of active columns */
      for (len = 0; len <= n; len++) cs_head[len] = 0;
      for (j = 1; j <= n; j++)
      {  len = vc_len[j];
         cs_prev[j] = 0;
         cs_next[j] = cs_head[len];
         if (cs_next[j] != 0) cs_prev[cs_next[j]] = j;
         cs_head[len] = j;
      }
done: /* return to the factorizing routine */
      return ret;
}

/***********************************************************************
*  find_pivot - choose a pivot element
*
*  This routine chooses a pivot element in the active submatrix of the
*  matrix U = P*V*Q.
*
*  It is assumed that on entry the matrix U has the following partially
*  triangularized form:
* 
*        1       k         n
*     1  x x x x x x x x x x
*        . x x x x x x x x x
*        . . x x x x x x x x
*        . . . x x x x x x x
*     k  . . . . * * * * * *
*        . . . . * * * * * *
*        . . . . * * * * * *
*        . . . . * * * * * *
*        . . . . * * * * * *
*     n  . . . . * * * * * *
* 
*  where rows and columns k, k+1, ..., n belong to the active submatrix
*  (elements of the active submatrix are marked by '*').
*
*  Since the matrix U = P*V*Q is not stored, the routine works with the
*  matrix V. It is assumed 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. It is also assumed that each active row of the
*  matrix V is in the set R[len], where len is number of non-zeros in
*  the row, and each active column of the matrix V is in the set C[len],
*  where len is number of non-zeros in the column (in the latter case
*  only elements of the active submatrix are counted; such elements are
*  marked by '*' on the figure above).
* 
*  For the reason of numerical stability the routine applies so called
*  threshold pivoting proposed by J.Reid. It is assumed that an element
*  v[i,j] can be selected as a pivot candidate if it is not very small
*  (in absolute value) among other elements in the same row, i.e. if it
*  satisfies to the stability condition |v[i,j]| >= tol * max|v[i,*]|,
*  where 0 < tol < 1 is a given tolerance.
* 
*  In order to keep sparsity of the matrix V the routine uses Markowitz