glplpf.c

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

/* glplpf.c (LP basis factorization, Schur complement version) */

/***********************************************************************
*  This code is part of GLPK (GNU Linear Programming Kit).
*
*  Copyright (C) 2000,01,02,03,04,05,06,07,08,2009 Andrew Makhorin,
*  Department for Applied Informatics, Moscow Aviation Institute,
*  Moscow, Russia. All rights reserved. E-mail: <[e-mail]>.
*
*  GLPK is free software: you can redistribute it and/or modify it
*  under the terms of the GNU General Public License as published by
*  the Free Software Foundation, either version 3 of the License, or
*  (at your option) any later version.
*
*  GLPK is distributed in the hope that it will be useful, but WITHOUT
*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
*  License for more details.
*
*  You should have received a copy of the GNU General Public License
*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
***********************************************************************/

#include "glplpf.h"
#include "glplib.h"
#define xfault xerror

#define _GLPLPF_DEBUG 0

/* CAUTION: DO NOT CHANGE THE LIMIT BELOW */

#define M_MAX 100000000 /* = 100*10^6 */
/* maximal order of the basis matrix */

/***********************************************************************
*  NAME
*
*  lpf_create_it - create LP basis factorization
*
*  SYNOPSIS
*
*  #include "glplpf.h"
*  LPF *lpf_create_it(void);
*
*  DESCRIPTION
*
*  The routine lpf_create_it creates a program object, which represents
*  a factorization of LP basis.
*
*  RETURNS
*
*  The routine lpf_create_it returns a pointer to the object created. */

LPF *lpf_create_it(void)
{     LPF *lpf;
#if _GLPLPF_DEBUG
      xprintf("lpf_create_it: warning: debug mode enabled\n");
#endif
      lpf = xmalloc(sizeof(LPF));
      lpf->valid = 0;
      lpf->m0_max = lpf->m0 = 0;
      lpf->luf = luf_create_it();
      lpf->m = 0;
      lpf->B = NULL;
      lpf->n_max = 50;
      lpf->n = 0;
      lpf->R_ptr = lpf->R_len = NULL;
      lpf->S_ptr = lpf->S_len = NULL;
      lpf->scf = NULL;
      lpf->P_row = lpf->P_col = NULL;
      lpf->Q_row = lpf->Q_col = NULL;
      lpf->v_size = 1000;
      lpf->v_ptr = 0;
      lpf->v_ind = NULL;
      lpf->v_val = NULL;
      lpf->work1 = lpf->work2 = NULL;
      return lpf;
}

/***********************************************************************
*  NAME
*
*  lpf_factorize - compute LP basis factorization
*
*  SYNOPSIS
*
*  #include "glplpf.h"
*  int lpf_factorize(LPF *lpf, int m, const int bh[], int (*col)
*     (void *info, int j, int ind[], double val[]), void *info);
*
*  DESCRIPTION
*
*  The routine lpf_factorize computes the factorization of the basis
*  matrix B specified by the routine col.
*
*  The parameter lpf specified the basis factorization data structure
*  created with the routine lpf_create_it.
*
*  The parameter m specifies the order of B, m > 0.
*
*  The array bh specifies the basis header: bh[j], 1 <= j <= m, is the
*  number of j-th column of B in some original matrix. The array bh is
*  optional and can be specified as NULL.
*
*  The formal routine col specifies the matrix B to be factorized. To
*  obtain j-th column of A the routine lpf_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 B 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  The factorization has been successfully computed.
*
*  LPF_ESING
*     The specified matrix is singular within the working precision.
*
*  LPF_ECOND
*     The specified matrix is ill-conditioned.
*
*  For more details see comments to the routine luf_factorize. */

int lpf_factorize(LPF *lpf, int m, const int bh[], int (*col)
      (void *info, int j, int ind[], double val[]), void *info)
{     int k, ret;
#if _GLPLPF_DEBUG
      int i, j, len, *ind;
      double *B, *val;
#endif
      xassert(bh == bh);
      if (m < 1)
         xfault("lpf_factorize: m = %d; invalid parameter\n", m);
      if (m > M_MAX)
         xfault("lpf_factorize: m = %d; matrix too big\n", m);
      lpf->m0 = lpf->m = m;
      /* invalidate the factorization */
      lpf->valid = 0;
      /* allocate/reallocate arrays, if necessary */
      if (lpf->R_ptr == NULL)
         lpf->R_ptr = xcalloc(1+lpf->n_max, sizeof(int));
      if (lpf->R_len == NULL)
         lpf->R_len = xcalloc(1+lpf->n_max, sizeof(int));
      if (lpf->S_ptr == NULL)
         lpf->S_ptr = xcalloc(1+lpf->n_max, sizeof(int));
      if (lpf->S_len == NULL)
         lpf->S_len = xcalloc(1+lpf->n_max, sizeof(int));
      if (lpf->scf == NULL)
         lpf->scf = scf_create_it(lpf->n_max);
      if (lpf->v_ind == NULL)
         lpf->v_ind = xcalloc(1+lpf->v_size, sizeof(int));
      if (lpf->v_val == NULL)
         lpf->v_val = xcalloc(1+lpf->v_size, sizeof(double));
      if (lpf->m0_max < m)
      {  if (lpf->P_row != NULL) xfree(lpf->P_row);
         if (lpf->P_col != NULL) xfree(lpf->P_col);
         if (lpf->Q_row != NULL) xfree(lpf->Q_row);
         if (lpf->Q_col != NULL) xfree(lpf->Q_col);
         if (lpf->work1 != NULL) xfree(lpf->work1);
         if (lpf->work2 != NULL) xfree(lpf->work2);
         lpf->m0_max = m + 100;
         lpf->P_row = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(int));
         lpf->P_col = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(int));
         lpf->Q_row = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(int));
         lpf->Q_col = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(int));
         lpf->work1 = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(double));
         lpf->work2 = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(double));
      }
      /* try to factorize the basis matrix */
      switch (luf_factorize(lpf->luf, m, col, info))
      {  case 0:
            break;
         case LUF_ESING:
            ret = LPF_ESING;
            goto done;
         case LUF_ECOND:
            ret = LPF_ECOND;
            goto done;
         default:
            xassert(lpf != lpf);
      }
      /* the basis matrix has been successfully factorized */
      lpf->valid = 1;
#if _GLPLPF_DEBUG
      /* store the basis matrix for debugging */
      if (lpf->B != NULL) xfree(lpf->B);
      xassert(m <= 32767);
      lpf->B = B = xcalloc(1+m*m, sizeof(double));
      ind = xcalloc(1+m, sizeof(int));
      val = xcalloc(1+m, sizeof(double));
      for (k = 1; k <= m * m; k++)
         B[k] = 0.0;
      for (j = 1; j <= m; j++)
      {  len = col(info, j, ind, val);
         xassert(0 <= len && len <= m);
         for (k = 1; k <= len; k++)
         {  i = ind[k];
            xassert(1 <= i && i <= m);
            xassert(B[(i - 1) * m + j] == 0.0);
            xassert(val[k] != 0.0);
            B[(i - 1) * m + j] = val[k];
         }
      }
      xfree(ind);
      xfree(val);
#endif
      /* B = B0, so there are no additional rows/columns */
      lpf->n = 0;
      /* reset the Schur complement factorization */
      scf_reset_it(lpf->scf);
      /* P := Q := I */
      for (k = 1; k <= m; k++)
      {  lpf->P_row[k] = lpf->P_col[k] = k;
         lpf->Q_row[k] = lpf->Q_col[k] = k;
      }
      /* make all SVA locations free */
      lpf->v_ptr = 1;
      ret = 0;
done: /* return to the calling program */
      return ret;
}

/***********************************************************************
*  The routine r_prod computes the product y := y + alpha * R * x,
*  where x is a n-vector, alpha is a scalar, y is a m0-vector.
*
*  Since matrix R is available by columns, the product is computed as
*  a linear combination:
*
*     y := y + alpha * (R[1] * x[1] + ... + R[n] * x[n]),
*
*  where R[j] is j-th column of R. */

static void r_prod(LPF *lpf, double y[], double a, const double x[])
{     int n = lpf->n;
      int *R_ptr = lpf->R_ptr;
      int *R_len = lpf->R_len;
      int *v_ind = lpf->v_ind;
      double *v_val = lpf->v_val;
      int j, beg, end, ptr;
      double t;
      for (j = 1; j <= n; j++)
      {  if (x[j] == 0.0) continue;
         /* y := y + alpha * R[j] * x[j] */
         t = a * x[j];
         beg = R_ptr[j];
         end = beg + R_len[j];
         for (ptr = beg; ptr < end; ptr++)
            y[v_ind[ptr]] += t * v_val[ptr];
      }
      return;
}

/***********************************************************************
*  The routine rt_prod computes the product y := y + alpha * R' * x,
*  where R' is a matrix transposed to R, x is a m0-vector, alpha is a
*  scalar, y is a n-vector.
*
*  Since matrix R is available by columns, the product components are
*  computed as inner products:
*
*     y[j] := y[j] + alpha * (j-th column of R) * x
*
*  for j = 1, 2, ..., n. */

static void rt_prod(LPF *lpf, double y[], double a, const double x[])
{     int n = lpf->n;
      int *R_ptr = lpf->R_ptr;
      int *R_len = lpf->R_len;
      int *v_ind = lpf->v_ind;
      double *v_val = lpf->v_val;
      int j, beg, end, ptr;
      double t;
      for (j = 1; j <= n; j++)
      {  /* t := (j-th column of R) * x */
         t = 0.0;
         beg = R_ptr[j];
         end = beg + R_len[j];
         for (ptr = beg; ptr < end; ptr++)
            t += v_val[ptr] * x[v_ind[ptr]];
         /* y[j] := y[j] + alpha * t */
         y[j] += a * t;
      }
      return;
}

/***********************************************************************
*  The routine s_prod computes the product y := y + alpha * S * x,
*  where x is a m0-vector, alpha is a scalar, y is a n-vector.
*
*  Since matrix S is available by rows, the product components are
*  computed as inner products:
*
*     y[i] = y[i] + alpha * (i-th row of S) * x
*
*  for i = 1, 2, ..., n. */

static void s_prod(LPF *lpf, double y[], double a, const double x[])
{     int n = lpf->n;
      int *S_ptr = lpf->S_ptr;
      int *S_len = lpf->S_len;
      int *v_ind = lpf->v_ind;
      double *v_val = lpf->v_val;
      int i, beg, end, ptr;
      double t;
      for (i = 1; i <= n; i++)
      {  /* t := (i-th row of S) * x */
         t = 0.0;
         beg = S_ptr[i];
         end = beg + S_len[i];
         for (ptr = beg; ptr < end; ptr++)
            t += v_val[ptr] * x[v_ind[ptr]];
         /* y[i] := y[i] + alpha * t */
         y[i] += a * t;
      }
      return;
}

/***********************************************************************
*  The routine st_prod computes the product y := y + alpha * S' * x,
*  where S' is a matrix transposed to S, x is a n-vector, alpha is a
*  scalar, y is m0-vector.
*