glplpf.c

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

*  Since matrix R is available by rows, the product is computed as a
*  linear combination:
*
*     y := y + alpha * (S'[1] * x[1] + ... + S'[n] * x[n]),
*
*  where S'[i] is i-th row of S. */

static void st_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++)
      {  if (x[i] == 0.0) continue;
         /* y := y + alpha * S'[i] * x[i] */
         t = a * x[i];
         beg = S_ptr[i];
         end = beg + S_len[i];
         for (ptr = beg; ptr < end; ptr++)
            y[v_ind[ptr]] += t * v_val[ptr];
      }
      return;
}

#if _GLPLPF_DEBUG
/***********************************************************************
*  The routine check_error computes the maximal relative error between
*  left- and right-hand sides for the system B * x = b (if tr is zero)
*  or B' * x = b (if tr is non-zero), where B' is a matrix transposed
*  to B. (This routine is intended for debugging only.) */

static void check_error(LPF *lpf, int tr, const double x[],
      const double b[])
{     int m = lpf->m;
      double *B = lpf->B;
      int i, j;
      double  d, dmax = 0.0, s, t, tmax;
      for (i = 1; i <= m; i++)
      {  s = 0.0;
         tmax = 1.0;
         for (j = 1; j <= m; j++)
         {  if (!tr)
               t = B[m * (i - 1) + j] * x[j];
            else
               t = B[m * (j - 1) + i] * x[j];
            if (tmax < fabs(t)) tmax = fabs(t);
            s += t;
         }
         d = fabs(s - b[i]) / tmax;
         if (dmax < d) dmax = d;
      }
      if (dmax > 1e-8)
         xprintf("%s: dmax = %g; relative error too large\n",
            !tr ? "lpf_ftran" : "lpf_btran", dmax);
      return;
}
#endif

/***********************************************************************
*  NAME
*
*  lpf_ftran - perform forward transformation (solve system B*x = b)
*
*  SYNOPSIS
*
*  #include "glplpf.h"
*  void lpf_ftran(LPF *lpf, double x[]);
*
*  DESCRIPTION
*
*  The routine lpf_ftran performs forward transformation, i.e. solves
*  the system B*x = b, where B is the basis matrix, x is the vector of
*  unknowns to be computed, b is the vector of right-hand sides.
*
*  On entry elements of the vector b should be stored in dense format
*  in locations x[1], ..., x[m], where m is the number of rows. On exit
*  the routine stores elements of the vector x in the same locations.
*
*  BACKGROUND
*
*  Solution of the system B * x = b can be obtained by solving the
*  following augmented system:
*
*     ( B  F^) ( x )   ( b )
*     (      ) (   ) = (   )
*     ( G^ H^) ( y )   ( 0 )
*
*  which, using the main equality, can be written as follows:
*
*       ( L0 0 ) ( U0 R )   ( x )   ( b )
*     P (      ) (      ) Q (   ) = (   )
*       ( S  I ) ( 0  C )   ( y )   ( 0 )
*
*  therefore,
*
*     ( x )      ( U0 R )-1 ( L0 0 )-1    ( b )
*     (   ) = Q' (      )   (      )   P' (   )
*     ( y )      ( 0  C )   ( S  I )      ( 0 )
*
*  Thus, computing the solution includes the following steps:
*
*  1. Compute
*
*     ( f )      ( b )
*     (   ) = P' (   )
*     ( g )      ( 0 )
*
*  2. Solve the system
*
*     ( f1 )   ( L0 0 )-1 ( f )      ( L0 0 ) ( f1 )   ( f )
*     (    ) = (      )   (   )  =>  (      ) (    ) = (   )
*     ( g1 )   ( S  I )   ( g )      ( S  I ) ( g1 )   ( g )
*
*     from which it follows that:
*
*     { L0 * f1      = f      f1 = inv(L0) * f
*     {                   =>
*     {  S * f1 + g1 = g      g1 = g - S * f1
*
*  3. Solve the system
*
*     ( f2 )   ( U0 R )-1 ( f1 )      ( U0 R ) ( f2 )   ( f1 )
*     (    ) = (      )   (    )  =>  (      ) (    ) = (    )
*     ( g2 )   ( 0  C )   ( g1 )      ( 0  C ) ( g2 )   ( g1 )
*
*     from which it follows that:
*
*     { U0 * f2 + R * g2 = f1      f2 = inv(U0) * (f1 - R * g2)
*     {                        =>
*     {           C * g2 = g1      g2 = inv(C) * g1
*
*  4. Compute
*
*     ( x )      ( f2 )
*     (   ) = Q' (    )
*     ( y )      ( g2 )                                               */

void lpf_ftran(LPF *lpf, double x[])
{     int m0 = lpf->m0;
      int m = lpf->m;
      int n  = lpf->n;
      int *P_col = lpf->P_col;
      int *Q_col = lpf->Q_col;
      double *fg = lpf->work1;
      double *f = fg;
      double *g = fg + m0;
      int i, ii;
#if _GLPLPF_DEBUG
      double *b;
#endif
      if (!lpf->valid)
         xfault("lpf_ftran: the factorization is not valid\n");
      xassert(0 <= m && m <= m0 + n);
#if _GLPLPF_DEBUG
      /* save the right-hand side vector */
      b = xcalloc(1+m, sizeof(double));
      for (i = 1; i <= m; i++) b[i] = x[i];
#endif
      /* (f g) := inv(P) * (b 0) */
      for (i = 1; i <= m0 + n; i++)
         fg[i] = ((ii = P_col[i]) <= m ? x[ii] : 0.0);
      /* f1 := inv(L0) * f */
      luf_f_solve(lpf->luf, 0, f);
      /* g1 := g - S * f1 */
      s_prod(lpf, g, -1.0, f);
      /* g2 := inv(C) * g1 */
      scf_solve_it(lpf->scf, 0, g);
      /* f2 := inv(U0) * (f1 - R * g2) */
      r_prod(lpf, f, -1.0, g);
      luf_v_solve(lpf->luf, 0, f);
      /* (x y) := inv(Q) * (f2 g2) */
      for (i = 1; i <= m; i++)
         x[i] = fg[Q_col[i]];
#if _GLPLPF_DEBUG
      /* check relative error in solution */
      check_error(lpf, 0, x, b);
      xfree(b);
#endif
      return;
}

/***********************************************************************
*  NAME
*
*  lpf_btran - perform backward transformation (solve system B'*x = b)
*
*  SYNOPSIS
*
*  #include "glplpf.h"
*  void lpf_btran(LPF *lpf, double x[]);
*
*  DESCRIPTION
*
*  The routine lpf_btran performs backward transformation, i.e. solves
*  the system B'*x = b, where B' is a matrix transposed to the basis
*  matrix B, x is the vector of unknowns to be computed, b is the vector
*  of right-hand sides.
*
*  On entry elements of the vector b should be stored in dense format
*  in locations x[1], ..., x[m], where m is the number of rows. On exit
*  the routine stores elements of the vector x in the same locations.
*
*  BACKGROUND
*
*  Solution of the system B' * x = b, where B' is a matrix transposed
*  to B, can be obtained by solving the following augmented system:
*
*     ( B  F^)T ( x )   ( b )
*     (      )  (   ) = (   )
*     ( G^ H^)  ( y )   ( 0 )
*
*  which, using the main equality, can be written as follows:
*
*      T ( U0 R )T ( L0 0 )T  T ( x )   ( b )
*     Q  (      )  (      )  P  (   ) = (   )
*        ( 0  C )  ( S  I )     ( y )   ( 0 )
*
*  or, equivalently, as follows:
*
*        ( U'0 0 ) ( L'0 S')    ( x )   ( b )
*     Q' (       ) (       ) P' (   ) = (   )
*        ( R'  C') ( 0   I )    ( y )   ( 0 )
*
*  therefore,
*
*     ( x )     ( L'0 S')-1 ( U'0 0 )-1   ( b )
*     (   ) = P (       )   (       )   Q (   )
*     ( y )     ( 0   I )   ( R'  C')     ( 0 )
*
*  Thus, computing the solution includes the following steps:
*
*  1. Compute
*
*     ( f )     ( b )
*     (   ) = Q (   )
*     ( g )     ( 0 )
*
*  2. Solve the system
*
*     ( f1 )   ( U'0 0 )-1 ( f )      ( U'0 0 ) ( f1 )   ( f )
*     (    ) = (       )   (   )  =>  (       ) (    ) = (   )
*     ( g1 )   ( R'  C')   ( g )      ( R'  C') ( g1 )   ( g )
*
*     from which it follows that:
*
*     { U'0 * f1           = f      f1 = inv(U'0) * f
*     {                         =>
*     { R'  * f1 + C' * g1 = g      g1 = inv(C') * (g - R' * f1)
*
*  3. Solve the system
*
*     ( f2 )   ( L'0 S')-1 ( f1 )      ( L'0 S') ( f2 )   ( f1 )
*     (    ) = (       )   (    )  =>  (       ) (    ) = (    )
*     ( g2 )   ( 0   I )   ( g1 )      ( 0   I ) ( g2 )   ( g1 )
*
*     from which it follows that:
*
*     { L'0 * f2 + S' * g2 = f1
*     {                          =>  f2 = inv(L'0) * ( f1 - S' * g2)
*     {                 g2 = g1
*
*  4. Compute
*
*     ( x )     ( f2 )
*     (   ) = P (    )
*     ( y )     ( g2 )                                                */

void lpf_btran(LPF *lpf, double x[])
{     int m0 = lpf->m0;
      int m = lpf->m;
      int n = lpf->n;
      int *P_row = lpf->P_row;
      int *Q_row = lpf->Q_row;
      double *fg = lpf->work1;
      double *f = fg;
      double *g = fg + m0;
      int i, ii;
#if _GLPLPF_DEBUG
      double *b;
#endif
      if (!lpf->valid)
         xfault("lpf_btran: the factorization is not valid\n");
      xassert(0 <= m && m <= m0 + n);
#if _GLPLPF_DEBUG
      /* save the right-hand side vector */
      b = xcalloc(1+m, sizeof(double));
      for (i = 1; i <= m; i++) b[i] = x[i];
#endif
      /* (f g) := Q * (b 0) */
      for (i = 1; i <= m0 + n; i++)
         fg[i] = ((ii = Q_row[i]) <= m ? x[ii] : 0.0);
      /* f1 := inv(U'0) * f */
      luf_v_solve(lpf->luf, 1, f);
      /* g1 := inv(C') * (g - R' * f1) */
      rt_prod(lpf, g, -1.0, f);
      scf_solve_it(lpf->scf, 1, g);
      /* g2 := g1 */
      g = g;
      /* f2 := inv(L'0) * (f1 - S' * g2) */
      st_prod(lpf, f, -1.0, g);
      luf_f_solve(lpf->luf, 1, f);
      /* (x y) := P * (f2 g2) */
      for (i = 1; i <= m; i++)
         x[i] = fg[P_row[i]];
#if _GLPLPF_DEBUG
      /* check relative error in solution */
      check_error(lpf, 1, x, b);
      xfree(b);
#endif
      return;
}

/***********************************************************************
*  The routine enlarge_sva enlarges the Sparse Vector Area to new_size
*  locations by reallocating the arrays v_ind and v_val. */

static void enlarge_sva(LPF *lpf, int new_size)
{     int v_size = lpf->v_size;
      int used = lpf->v_ptr - 1;
      int *v_ind = lpf->v_ind;
      double *v_val = lpf->v_val;
      xassert(v_size < new_size);
      while (v_size < new_size) v_size += v_size;