glpipm.c

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

/* glpipm.c */

/***********************************************************************
*  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 "glpipm.h"
#include "glplib.h"
#include "glpmat.h"

/*----------------------------------------------------------------------
-- ipm_main - solve LP with primal-dual interior-point method.
--
-- *Synopsis*
--
-- #include "glpipm.h"
-- int ipm_main(int m, int n, int A_ptr[], int A_ind[], double A_val[],
--    double b[], double c[], double x[], double y[], double z[]);
--
-- *Description*
--
-- The routine ipm_main is a *tentative* implementation of primal-dual
-- interior point method for solving linear programming problems.
--
-- The routine ipm_main assumes the following *standard* formulation of
-- LP problem to be solved:
--
--    minimize
--
--       F = c[0] + c[1]*x[1] + c[2]*x[2] + ... + c[n]*x[n]
--
--    subject to linear constraints
--
--       a[1,1]*x[1] + a[1,2]*x[2] + ... + a[1,n]*x[n] = b[1]
--       a[2,1]*x[1] + a[2,2]*x[2] + ... + a[2,n]*x[n] = b[2]
--             . . . . . .
--       a[m,1]*x[1] + a[m,2]*x[2] + ... + a[m,n]*x[n] = b[m]
--
--    and non-negative variables
--
--       x[1] >= 0, x[2] >= 0, ..., x[n] >= 0
--
-- where:
-- F                    is objective function;
-- x[1], ..., x[n]      are (structural) variables;
-- c[0]                 is constant term of the objective function;
-- c[1], ..., c[n]      are objective coefficients;
-- a[1,1], ..., a[m,n]  are constraint coefficients;
-- b[1], ..., b[n]      are right-hand sides.
--
-- The parameter m is the number of rows (constraints).
--
-- The parameter n is the number of columns (variables).
--
-- The arrays A_ptr, A_ind, and A_val specify the mxn constraint matrix
-- A in storage-by-rows format. These arrays are not changed on exit.
--
-- The array b specifies the vector of right-hand sides b, which should
-- be stored in locations b[1], ..., b[m]. This array is not changed on
-- exit.
--
-- The array c specifies the vector of objective coefficients c, which
-- should be stored in locations c[1], ..., c[n], and the constant term
-- of the objective function, which should be stored in location c[0].
-- This array is not changed on exit.
--
-- The solution is three vectors x, y, and z, which are stored by the
-- routine in the arrays x, y, and z, respectively. These vectors
-- correspond to the best primal-dual point found during optimization.
-- They are approximate solution of the following system (which is the
-- Karush-Kuhn-Tucker optimality conditions):
--
--    A*x      = b      (primal feasibility condition)
--    A'*y + z = c      (dual feasibility condition)
--    x'*z     = 0      (primal-dual complementarity condition)
--    x >= 0, z >= 0    (non-negativity condition)
--
-- where:
-- x[1], ..., x[n]      are primal (structural) variables;
-- y[1], ..., y[m]      are dual variables (Lagrange multipliers) for
--                      equality constraints;
-- z[1], ..., z[n]      are dual variables (Lagrange multipliers) for
--                      non-negativity constraints.
--
-- *Returns*
--
-- The routine ipm_main returns one of the following codes:
--
-- 0 - optimal solution found;
-- 1 - problem has no feasible (primal or dual) solution;
-- 2 - no convergence;
-- 3 - iteration limit exceeded;
-- 4 - numeric instability on solving Newtonian system.
--
-- In case of non-zero return code the routine returns the best point,
-- which has been reached during optimization. */

#define ITER_MAX 100
/* maximal number of iterations */

struct dsa
{     /* working area used by interior-point routines */
      int m;
      /* number of rows (equality constraints) */
      int n;
      /* number of columns (structural variables) */
      int *A_ptr; /* int A_ptr[1+m+1]; */
      int *A_ind; /* int A_ind[A_ptr[m+1]]; */
      double *A_val; /* double A_val[A_ptr[m+1]]; */
      /* mxn matrix A in storage-by-rows format */
      double *b; /* double b[1+m]; */
      /* m-vector b of right hand-sides */
      double *c; /* double c[1+n]; */
      /* n-vector c of objective coefficients; c[0] is constant term of
         the objective function */
      double *x; /* double x[1+n]; */
      double *y; /* double y[1+m]; */
      double *z; /* double z[1+n]; */
      /* current point in primal-dual space; the best point on exit */
      double *D; /* double D[1+n]; */
      /* nxn diagonal matrix D = X*inv(Z), where X = diag(x[j]) and
         Z = diag(z[j]) */
      int *P; /* int P[1+m+m]; */
      /* mxm permutation matrix P used to minimize fill-in in Cholesky
         factorization */
      int *S_ptr; /* int S_ptr[1+m+1]; */
      int *S_ind; /* int S_ind[S_ptr[m+1]]; */
      double *S_val; /* double S_val[S_ptr[m+1]]; */
      double *S_diag; /* double S_diag[1+m]; */
      /* mxm symmetric matrix S = P*A*D*A'*P' whose upper triangular
         part without diagonal elements is stored in S_ptr, S_ind, and
         S_val in storage-by-rows format, diagonal elements are stored
         in S_diag */
      int *U_ptr; /* int U_ptr[1+m+1]; */
      int *U_ind; /* int U_ind[U_ptr[m+1]]; */
      double *U_val; /* double U_val[U_ptr[m+1]]; */
      double *U_diag; /* double U_diag[1+m]; */
      /* mxm upper triangular matrix U defining Cholesky factorization
         S = U'*U; its non-diagonal elements are stored in U_ptr, U_ind,
         U_val in storage-by-rows format, diagonal elements are stored
         in U_diag */
      int iter;
      /* iteration number (0, 1, 2, ...); iter = 0 corresponds to the
         initial point */
      double obj;
      /* current value of the objective function */
      double rpi;
      /* relative primal infeasibility rpi = ||A*x-b||/(1+||b||) */
      double rdi;
      /* relative dual infeasibility rdi = ||A'*y+z-c||/(1+||c||) */
      double gap;
      /* primal-dual gap = |c'*x-b'*y|/(1+|c'*x|) which is a relative
         difference between primal and dual objective functions */
      double phi;
      /* merit function phi = ||A*x-b||/max(1,||b||) +
                            + ||A'*y+z-c||/max(1,||c||) +
                            + |c'*x-b'*y|/max(1,||b||,||c||) */
      double mu;
      /* duality measure mu = x'*z/n (used as barrier parameter) */
      double rmu;
      /* rmu = max(||A*x-b||,||A'*y+z-c||)/mu */
      double rmu0;
      /* the initial value of rmu on iteration 0 */
      double *phi_min; /* double phi_min[1+ITER_MAX]; */
      /* phi_min[k] = min(phi[k]), where phi[k] is the value of phi on
         k-th iteration, 0 <= k <= iter */
      int best_iter;
      /* iteration number, on which the value of phi reached its best
         (minimal) value */
      double *best_x; /* double best_x[1+n]; */
      double *best_y; /* double best_y[1+m]; */
      double *best_z; /* double best_z[1+n]; */
      /* best point (in the sense of the merit function phi) which has
         been reached on iteration iter_best */
      double best_obj;
      /* objective value at the best point */
      double *dx_aff; /* double dx_aff[1+n]; */
      double *dy_aff; /* double dy_aff[1+m]; */
      double *dz_aff; /* double dz_aff[1+n]; */
      /* affine scaling direction */
      double alfa_aff_p, alfa_aff_d;
      /* maximal primal and dual stepsizes in affine scaling direction,
         on which x and z are still non-negative */
      double mu_aff;
      /* duality measure mu_aff = x_aff'*z_aff/n in the boundary point
         x_aff' = x+alfa_aff_p*dx_aff, z_aff' = z+alfa_aff_d*dz_aff */
      double sigma;
      /* Mehrotra's heuristic parameter (0 <= sigma <= 1) */
      double *dx_cc; /* double dx_cc[1+n]; */
      double *dy_cc; /* double dy_cc[1+m]; */
      double *dz_cc; /* double dz_cc[1+n]; */
      /* centering corrector direction */
      double *dx; /* double dx[1+n]; */
      double *dy; /* double dy[1+m]; */
      double *dz; /* double dz[1+n]; */
      /* final combined direction dx = dx_aff+dx_cc, dy = dy_aff+dy_cc,
         dz = dz_aff+dz_cc */
      double alfa_max_p;
      double alfa_max_d;
      /* maximal primal and dual stepsizes in combined direction, on
         which x and z are still non-negative */
};

/*----------------------------------------------------------------------
-- initialize - allocate and initialize working area.
--
-- This routine allocates and initializes the working area used by all
-- interior-point method routines. */

static void initialize(struct dsa *dsa, int m, int n, int A_ptr[],
      int A_ind[], double A_val[], double b[], double c[], double x[],
      double y[], double z[])
{     int i;
      dsa->m = m;
      dsa->n = n;
      dsa->A_ptr = A_ptr;
      dsa->A_ind = A_ind;
      dsa->A_val = A_val;
      xprintf("lpx_interior: A has %d non-zeros\n", A_ptr[m+1]-1);
      dsa->b = b;
      dsa->c = c;
      dsa->x = x;
      dsa->y = y;
      dsa->z = z;
      dsa->D = xcalloc(1+n, sizeof(double));
      /* P := I */
      dsa->P = xcalloc(1+m+m, sizeof(int));
      for (i = 1; i <= m; i++) dsa->P[i] = dsa->P[m+i] = i;
      /* S := A*A', symbolically */
      dsa->S_ptr = xcalloc(1+m+1, sizeof(int));
      dsa->S_ind = adat_symbolic(m, n, dsa->P, dsa->A_ptr, dsa->A_ind,
         dsa->S_ptr);
      xprintf("lpx_interior: S has %d non-zeros (upper triangle)\n",
         dsa->S_ptr[m+1]-1 + m);
      /* determine P using minimal degree algorithm */
      xprintf("lpx_interior: minimal degree ordering...\n");
      min_degree(m, dsa->S_ptr, dsa->S_ind, dsa->P);
      /* S = P*A*A'*P', symbolically */
      xfree(dsa->S_ind);
      dsa->S_ind = adat_symbolic(m, n, dsa->P, dsa->A_ptr, dsa->A_ind,
         dsa->S_ptr);
      dsa->S_val = xcalloc(dsa->S_ptr[m+1], sizeof(double));
      dsa->S_diag = xcalloc(1+m, sizeof(double));
      /* compute Cholesky factorization S = U'*U, symbolically */
      xprintf("lpx_interior: computing Cholesky factorization...\n");
      dsa->U_ptr = xcalloc(1+m+1, sizeof(int));
      dsa->U_ind = chol_symbolic(m, dsa->S_ptr, dsa->S_ind, dsa->U_ptr);
      xprintf("lpx_interior: U has %d non-zeros\n",
         dsa->U_ptr[m+1]-1 + m);
      dsa->U_val = xcalloc(dsa->U_ptr[m+1], sizeof(double));
      dsa->U_diag = xcalloc(1+m, sizeof(double));
      dsa->iter = 0;
      dsa->obj = 0.0;
      dsa->rpi = 0.0;
      dsa->rdi = 0.0;
      dsa->gap = 0.0;
      dsa->phi = 0.0;
      dsa->mu = 0.0;
      dsa->rmu = 0.0;
      dsa->rmu0 = 0.0;
      dsa->phi_min = xcalloc(1+ITER_MAX, sizeof(double));
      dsa->best_iter = 0;
      dsa->best_x = xcalloc(1+n, sizeof(double));
      dsa->best_y = xcalloc(1+m, sizeof(double));
      dsa->best_z = xcalloc(1+n, sizeof(double));
      dsa->best_obj = 0.0;
      dsa->dx_aff = xcalloc(1+n, sizeof(double));
      dsa->dy_aff = xcalloc(1+m, sizeof(double));
      dsa->dz_aff = xcalloc(1+n, sizeof(double));
      dsa->alfa_aff_p = 0.0;
      dsa->alfa_aff_d = 0.0;
      dsa->mu_aff = 0.0;
      dsa->sigma = 0.0;
      dsa->dx_cc = xcalloc(1+n, sizeof(double));
      dsa->dy_cc = xcalloc(1+m, sizeof(double));
      dsa->dz_cc = xcalloc(1+n, sizeof(double));
      dsa->dx = dsa->dx_aff;
      dsa->dy = dsa->dy_aff;
      dsa->dz = dsa->dz_aff;
      dsa->alfa_max_p = 0.0;
      dsa->alfa_max_d = 0.0;
      return;
}

/*----------------------------------------------------------------------
-- A_by_vec - compute y = A*x.
--
-- This routine computes the matrix-vector product y = A*x, where A is
-- the constraint matrix. */

static void A_by_vec(struct dsa *dsa, double x[], double y[])
{     /* compute y = A*x */
      int m = dsa->m;
      int *A_ptr = dsa->A_ptr;
      int *A_ind = dsa->A_ind;
      double *A_val = dsa->A_val;
      int i, t, beg, end;
      double temp;
      for (i = 1; i <= m; i++)
      {  temp = 0.0;
         beg = A_ptr[i], end = A_ptr[i+1];
         for (t = beg; t < end; t++) temp += A_val[t] * x[A_ind[t]];
         y[i] = temp;
      }
      return;
}

/*----------------------------------------------------------------------
-- AT_by_vec - compute y = A'*x.
--
-- This routine computes the matrix-vector product y = A'*x, where A' is