glpluf.c

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

            sv_val[j_ptr] = sv_val[i_ptr];
            /* increase length of the j-th column */
            vc_len[j]++;
         }
      }
      /* now columns are placed in the sparse vector area behind rows
         in the order n+1, n+2, ..., n+n; so insert column nodes in the
         addressing list using this order */
      for (k = n+1; k <= n+n; k++)
      {  sv_prev[k] = k-1;
         sv_next[k] = k+1;
      }
      sv_prev[n+1] = luf->sv_tail;
      sv_next[luf->sv_tail] = n+1;
      sv_next[n+n] = 0;
      luf->sv_tail = n+n;
done: /* return to the factorizing routine */
      return ret;
}

/***********************************************************************
*  build_f_rows - build the matrix F in row-wise format
*
*  This routine builds the row-wise representation of the matrix F using
*  its column-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_f_rows(LUF *luf)
{     int n = luf->n;
      int *fr_ptr = luf->fr_ptr;
      int *fr_len = luf->fr_len;
      int *fc_ptr = luf->fc_ptr;
      int *fc_len = luf->fc_len;
      int *sv_ind = luf->sv_ind;
      double *sv_val = luf->sv_val;
      int ret = 0;
      int i, j, j_beg, j_end, j_ptr, ptr, nnz;
      /* clear rows of the matrix F */
      for (i = 1; i <= n; i++) fr_len[i] = 0;
      /* count non-zeros in rows of the matrix F; count total number of
         non-zeros in this matrix */
      nnz = 0;
      for (j = 1; j <= n; j++)
      {  /* walk through elements of the j-th column and count non-zeros
            in the corresponding rows */
         j_beg = fc_ptr[j];
         j_end = j_beg + fc_len[j] - 1;
         for (j_ptr = j_beg; j_ptr <= j_end; j_ptr++)
            fr_len[sv_ind[j_ptr]]++;
         /* increase total number of non-zeros */
         nnz += fc_len[j];
      }
      /* store total number of non-zeros */
      luf->nnz_f = nnz;
      /* check for free locations */
      if (luf->sv_end - luf->sv_beg < nnz)
      {  /* overflow of the sparse vector area */
         ret = 1;
         goto done;
      }
      /* allocate rows of the matrix F */
      for (i = 1; i <= n; i++)
      {  /* set pointer to the end of the i-th row; later this pointer
            will be set to the beginning of the i-th row */
         fr_ptr[i] = luf->sv_end;
         /* reserve locations for the i-th row */
         luf->sv_end -= fr_len[i];
      }
      /* build the matrix F in row-wise format using this matrix in
         column-wise format */
      for (j = 1; j <= n; j++)
      {  /* walk through elements of the j-th column */
         j_beg = fc_ptr[j];
         j_end = j_beg + fc_len[j] - 1;
         for (j_ptr = j_beg; j_ptr <= j_end; j_ptr++)
         {  /* get row index */
            i = sv_ind[j_ptr];
            /* store element in the i-th row */
            ptr = --fr_ptr[i];
            sv_ind[ptr] = j;
            sv_val[ptr] = sv_val[j_ptr];
         }
      }
done: /* return to the factorizing routine */
      return ret;
}

/***********************************************************************
*  NAME
*
*  luf_factorize - compute LU-factorization
*
*  SYNOPSIS
*
*  #include "glpluf.h"
*  int luf_factorize(LUF *luf, int n, int (*col)(void *info, int j,
*     int ind[], double val[]), void *info);
*
*  DESCRIPTION
*
*  The routine luf_factorize computes LU-factorization of a specified
*  square matrix A.
*
*  The parameter luf specifies LU-factorization program object created
*  by the routine luf_create_it.
*
*  The parameter n specifies the order of A, n > 0.
*
*  The formal routine col specifies the matrix A to be factorized. To
*  obtain j-th column of A the routine luf_factorize calls the routine
*  col with the parameter j (1 <= j <= n). In response the routine col
*  should store row indices and numerical values of non-zero elements
*  of j-th column of A to locations ind[1,...,len] and val[1,...,len],
*  respectively, where len is the number of non-zeros in j-th column
*  returned on exit. Neither zero nor duplicate elements are allowed.
*
*  The parameter info is a transit pointer passed to the routine col.
*
*  RETURNS
*
*  0  LU-factorization has been successfully computed.
*
*  LUF_ESING
*     The specified matrix is singular within the working precision.
*     (On some elimination step the active submatrix is exactly zero,
*     so no pivot can be chosen.)
*
*  LUF_ECOND
*     The specified matrix is ill-conditioned.
*     (On some elimination step too intensive growth of elements of the
*     active submatix has been detected.)
*
*  If matrix A is well scaled, the return code LUF_ECOND may also mean
*  that the threshold pivoting tolerance piv_tol should be increased.
*
*  In case of non-zero return code the factorization becomes invalid.
*  It should not be used in other operations until the cause of failure
*  has been eliminated and the factorization has been recomputed again
*  with the routine luf_factorize.
*
*  REPAIRING SINGULAR MATRIX
*
*  If the routine luf_factorize returns non-zero code, 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 this routine is used for factorizing the basis matrix
*  within the simplex method procedure.
*
*  On exit linearly dependent columns of the (partially transformed)
*  matrix U have numbers rank+1, rank+2, ..., n, where rank is estimated
*  rank of the matrix A stored by the routine to the member luf->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 numbers qq_col[rank+1], qq_col[rank+2], ..., qq_col[n], where
*  qq_col is the column-like representation of the permutation matrix Q.
*  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 pp_row[j], where pp_row is
*  the row-like representation of the permutation matrix P. Thus, each
*  j-th linearly dependent column of the matrix U should be replaced by
*  column of the unity matrix with the number pp_row[j].
*
*  The code that repairs the matrix A may look like follows:
*
*     for (j = rank+1; j <= n; j++)
*     {  replace the column qq_col[j] of the matrix A by the column
*        pp_row[j] of the unity matrix;
*     }
*
*  where rank, pp_row, and qq_col are members of the structure LUF. */

int luf_factorize(LUF *luf, int n, int (*col)(void *info, int j,
      int ind[], double val[]), void *info)
{     int *pp_row, *pp_col, *qq_row, *qq_col;
      double max_gro = luf->max_gro;
      int i, j, k, p, q, t, ret;
      if (n < 1)
         xfault("luf_factorize: n = %d; invalid parameter\n", n);
      if (n > N_MAX)
         xfault("luf_factorize: n = %d; matrix too big\n", n);
      /* invalidate the factorization */
      luf->valid = 0;
      /* reallocate arrays, if necessary */
      reallocate(luf, n);
      pp_row = luf->pp_row;
      pp_col = luf->pp_col;
      qq_row = luf->qq_row;
      qq_col = luf->qq_col;
      /* estimate initial size of the SVA, if not specified */
      if (luf->sv_size == 0 && luf->new_sva == 0)
         luf->new_sva = 5 * (n + 10);
more: /* reallocate the sparse vector area, if required */
      if (luf->new_sva > 0)
      {  if (luf->sv_ind != NULL) xfree(luf->sv_ind);
         if (luf->sv_val != NULL) xfree(luf->sv_val);
         luf->sv_size = luf->new_sva;
         luf->sv_ind = xcalloc(1+luf->sv_size, sizeof(int));
         luf->sv_val = xcalloc(1+luf->sv_size, sizeof(double));
         luf->new_sva = 0;
      }
      /* initialize LU-factorization data structures */
      if (initialize(luf, col, info))
      {  /* overflow of the sparse vector area */
         luf->new_sva = luf->sv_size + luf->sv_size;
         xassert(luf->new_sva > luf->sv_size);
         goto more;
      }
      /* main elimination loop */
      for (k = 1; k <= n; k++)
      {  /* choose a pivot element v[p,q] */
         if (find_pivot(luf, &p, &q))
         {  /* no pivot can be chosen, because the active submatrix is
               exactly zero */
            luf->rank = k - 1;
            ret = LUF_ESING;
            goto done;
         }
         /* 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 = pp_col[p], j = qq_row[q];
         xassert(k <= i && i <= n && k <= j && j <= n);
         /* permute k-th and i-th rows of the matrix U */
         t = pp_row[k];
         pp_row[i] = t, pp_col[t] = i;
         pp_row[k] = p, pp_col[p] = k;
         /* permute k-th and j-th columns of the matrix U */
         t = qq_col[k];
         qq_col[j] = t, qq_row[t] = j;
         qq_col[k] = q, qq_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] */
         if (eliminate(luf, p, q))
         {  /* overflow of the sparse vector area */
            luf->new_sva = luf->sv_size + luf->sv_size;
            xassert(luf->new_sva > luf->sv_size);
            goto more;
         }
         /* check relative growth of elements of the matrix V */
         if (luf->big_v > max_gro * luf->max_a)
         {  /* the growth is too intensive, therefore most probably the
               matrix A is ill-conditioned */
            luf->rank = k - 1;
            ret = LUF_ECOND;
            goto done;
         }
      }
      /* now the matrix U = P*V*Q is upper triangular, the matrix V has
         been built in row-wise format, and the matrix F has been built
         in column-wise format */
      /* defragment the sparse vector area in order to merge all free
         locations in one continuous extent */
      luf_defrag_sva(luf);
      /* build the matrix V in column-wise format */
      if (build_v_cols(luf))
      {  /* overflow of the sparse vector area */
         luf->new_sva = luf->sv_size + luf->sv_size;
         xassert(luf->new_sva > luf->sv_size);
         goto more;
      }
      /* build the matrix F in row-wise format */
      if (build_f_rows(luf))
      {  /* overflow of the sparse vector area */
         luf->new_sva = luf->sv_size + luf->sv_size;
         xassert(luf->new_sva > luf->sv_size);
         goto more;
      }
      /* the LU-factorization has been successfully computed */
      luf->valid = 1;
      luf->rank = n;
      ret = 0;
      /* if there are few free locations in the sparse vector area, try
         increasing its size in the future */
      t = 3 * (n + luf->nnz_v) + 2 * luf->nnz_f;
      if (luf->sv_size < t)
      {  luf->new_sva = luf->sv_size;
         while (luf->new_sva < t)
         {  k = luf->new_sva;
            luf->new_sva = k + k;
            xassert(luf->new_sva > k);
         }
      }
done: /* return to the calling program */
      return ret;
}

/***********************************************************************
*  NAME
*
*  luf_f_solve - solve system F*x = b or F'*x = b
*
*  SYNOPSIS
*
*  #include "glpluf.h"
*  void luf_f_solve(LUF *luf, int tr, double x[]);
*
*  DESCRIPTION
*
*  The routine luf_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 luf, 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 luf_f_solve(LUF *luf, int tr, double x[])
{     int n = luf->n;
      int *fr_ptr = luf->fr_ptr;
      int *fr_len = luf->fr_len;
      int *fc_ptr = luf->fc_ptr;
      int *fc_len = luf->fc_len;
      int *pp_row = luf->pp_row;
      int *sv_ind = luf->sv_ind;
      double *sv_val = luf->sv_val;
      int i, j, k, beg, end, ptr;
      double xk;
      if (!luf->valid)
         xfault("luf_f_solve: LU-factorization is not valid\n");
      if (!tr)
      {  /* solve the system F*x = b */
         for (j = 1; j <= n; j++)
         {  k = pp_ro