alpar@1: /* glpipm.c */ alpar@1: alpar@1: /*********************************************************************** alpar@1: * This code is part of GLPK (GNU Linear Programming Kit). alpar@1: * alpar@1: * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, alpar@1: * 2009, 2010 Andrew Makhorin, Department for Applied Informatics, alpar@1: * Moscow Aviation Institute, Moscow, Russia. All rights reserved. alpar@1: * E-mail: . alpar@1: * alpar@1: * GLPK is free software: you can redistribute it and/or modify it alpar@1: * under the terms of the GNU General Public License as published by alpar@1: * the Free Software Foundation, either version 3 of the License, or alpar@1: * (at your option) any later version. alpar@1: * alpar@1: * GLPK is distributed in the hope that it will be useful, but WITHOUT alpar@1: * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY alpar@1: * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public alpar@1: * License for more details. alpar@1: * alpar@1: * You should have received a copy of the GNU General Public License alpar@1: * along with GLPK. If not, see . alpar@1: ***********************************************************************/ alpar@1: alpar@1: #include "glpipm.h" alpar@1: #include "glpmat.h" alpar@1: alpar@1: #define ITER_MAX 100 alpar@1: /* maximal number of iterations */ alpar@1: alpar@1: struct csa alpar@1: { /* common storage area */ alpar@1: /*--------------------------------------------------------------*/ alpar@1: /* LP data */ alpar@1: int m; alpar@1: /* number of rows (equality constraints) */ alpar@1: int n; alpar@1: /* number of columns (structural variables) */ alpar@1: int *A_ptr; /* int A_ptr[1+m+1]; */ alpar@1: int *A_ind; /* int A_ind[A_ptr[m+1]]; */ alpar@1: double *A_val; /* double A_val[A_ptr[m+1]]; */ alpar@1: /* mxn-matrix A in storage-by-rows format */ alpar@1: double *b; /* double b[1+m]; */ alpar@1: /* m-vector b of right-hand sides */ alpar@1: double *c; /* double c[1+n]; */ alpar@1: /* n-vector c of objective coefficients; c[0] is constant term of alpar@1: the objective function */ alpar@1: /*--------------------------------------------------------------*/ alpar@1: /* LP solution */ alpar@1: double *x; /* double x[1+n]; */ alpar@1: double *y; /* double y[1+m]; */ alpar@1: double *z; /* double z[1+n]; */ alpar@1: /* current point in primal-dual space; the best point on exit */ alpar@1: /*--------------------------------------------------------------*/ alpar@1: /* control parameters */ alpar@1: const glp_iptcp *parm; alpar@1: /*--------------------------------------------------------------*/ alpar@1: /* working arrays and variables */ alpar@1: double *D; /* double D[1+n]; */ alpar@1: /* diagonal nxn-matrix D = X*inv(Z), where X = diag(x[j]) and alpar@1: Z = diag(z[j]) */ alpar@1: int *P; /* int P[1+m+m]; */ alpar@1: /* permutation mxm-matrix P used to minimize fill-in in Cholesky alpar@1: factorization */ alpar@1: int *S_ptr; /* int S_ptr[1+m+1]; */ alpar@1: int *S_ind; /* int S_ind[S_ptr[m+1]]; */ alpar@1: double *S_val; /* double S_val[S_ptr[m+1]]; */ alpar@1: double *S_diag; /* double S_diag[1+m]; */ alpar@1: /* symmetric mxm-matrix S = P*A*D*A'*P' whose upper triangular alpar@1: part without diagonal elements is stored in S_ptr, S_ind, and alpar@1: S_val in storage-by-rows format, diagonal elements are stored alpar@1: in S_diag */ alpar@1: int *U_ptr; /* int U_ptr[1+m+1]; */ alpar@1: int *U_ind; /* int U_ind[U_ptr[m+1]]; */ alpar@1: double *U_val; /* double U_val[U_ptr[m+1]]; */ alpar@1: double *U_diag; /* double U_diag[1+m]; */ alpar@1: /* upper triangular mxm-matrix U defining Cholesky factorization alpar@1: S = U'*U; its non-diagonal elements are stored in U_ptr, U_ind, alpar@1: U_val in storage-by-rows format, diagonal elements are stored alpar@1: in U_diag */ alpar@1: int iter; alpar@1: /* iteration number (0, 1, 2, ...); iter = 0 corresponds to the alpar@1: initial point */ alpar@1: double obj; alpar@1: /* current value of the objective function */ alpar@1: double rpi; alpar@1: /* relative primal infeasibility rpi = ||A*x-b||/(1+||b||) */ alpar@1: double rdi; alpar@1: /* relative dual infeasibility rdi = ||A'*y+z-c||/(1+||c||) */ alpar@1: double gap; alpar@1: /* primal-dual gap = |c'*x-b'*y|/(1+|c'*x|) which is a relative alpar@1: difference between primal and dual objective functions */ alpar@1: double phi; alpar@1: /* merit function phi = ||A*x-b||/max(1,||b||) + alpar@1: + ||A'*y+z-c||/max(1,||c||) + alpar@1: + |c'*x-b'*y|/max(1,||b||,||c||) */ alpar@1: double mu; alpar@1: /* duality measure mu = x'*z/n (used as barrier parameter) */ alpar@1: double rmu; alpar@1: /* rmu = max(||A*x-b||,||A'*y+z-c||)/mu */ alpar@1: double rmu0; alpar@1: /* the initial value of rmu on iteration 0 */ alpar@1: double *phi_min; /* double phi_min[1+ITER_MAX]; */ alpar@1: /* phi_min[k] = min(phi[k]), where phi[k] is the value of phi on alpar@1: k-th iteration, 0 <= k <= iter */ alpar@1: int best_iter; alpar@1: /* iteration number, on which the value of phi reached its best alpar@1: (minimal) value */ alpar@1: double *best_x; /* double best_x[1+n]; */ alpar@1: double *best_y; /* double best_y[1+m]; */ alpar@1: double *best_z; /* double best_z[1+n]; */ alpar@1: /* best point (in the sense of the merit function phi) which has alpar@1: been reached on iteration iter_best */ alpar@1: double best_obj; alpar@1: /* objective value at the best point */ alpar@1: double *dx_aff; /* double dx_aff[1+n]; */ alpar@1: double *dy_aff; /* double dy_aff[1+m]; */ alpar@1: double *dz_aff; /* double dz_aff[1+n]; */ alpar@1: /* affine scaling direction */ alpar@1: double alfa_aff_p, alfa_aff_d; alpar@1: /* maximal primal and dual stepsizes in affine scaling direction, alpar@1: on which x and z are still non-negative */ alpar@1: double mu_aff; alpar@1: /* duality measure mu_aff = x_aff'*z_aff/n in the boundary point alpar@1: x_aff' = x+alfa_aff_p*dx_aff, z_aff' = z+alfa_aff_d*dz_aff */ alpar@1: double sigma; alpar@1: /* Mehrotra's heuristic parameter (0 <= sigma <= 1) */ alpar@1: double *dx_cc; /* double dx_cc[1+n]; */ alpar@1: double *dy_cc; /* double dy_cc[1+m]; */ alpar@1: double *dz_cc; /* double dz_cc[1+n]; */ alpar@1: /* centering corrector direction */ alpar@1: double *dx; /* double dx[1+n]; */ alpar@1: double *dy; /* double dy[1+m]; */ alpar@1: double *dz; /* double dz[1+n]; */ alpar@1: /* final combined direction dx = dx_aff+dx_cc, dy = dy_aff+dy_cc, alpar@1: dz = dz_aff+dz_cc */ alpar@1: double alfa_max_p; alpar@1: double alfa_max_d; alpar@1: /* maximal primal and dual stepsizes in combined direction, on alpar@1: which x and z are still non-negative */ alpar@1: }; alpar@1: alpar@1: /*********************************************************************** alpar@1: * initialize - allocate and initialize common storage area alpar@1: * alpar@1: * This routine allocates and initializes the common storage area (CSA) alpar@1: * used by interior-point method routines. */ alpar@1: alpar@1: static void initialize(struct csa *csa) alpar@1: { int m = csa->m; alpar@1: int n = csa->n; alpar@1: int i; alpar@1: if (csa->parm->msg_lev >= GLP_MSG_ALL) alpar@1: xprintf("Matrix A has %d non-zeros\n", csa->A_ptr[m+1]-1); alpar@1: csa->D = xcalloc(1+n, sizeof(double)); alpar@1: /* P := I */ alpar@1: csa->P = xcalloc(1+m+m, sizeof(int)); alpar@1: for (i = 1; i <= m; i++) csa->P[i] = csa->P[m+i] = i; alpar@1: /* S := A*A', symbolically */ alpar@1: csa->S_ptr = xcalloc(1+m+1, sizeof(int)); alpar@1: csa->S_ind = adat_symbolic(m, n, csa->P, csa->A_ptr, csa->A_ind, alpar@1: csa->S_ptr); alpar@1: if (csa->parm->msg_lev >= GLP_MSG_ALL) alpar@1: xprintf("Matrix S = A*A' has %d non-zeros (upper triangle)\n", alpar@1: csa->S_ptr[m+1]-1 + m); alpar@1: /* determine P using specified ordering algorithm */ alpar@1: if (csa->parm->ord_alg == GLP_ORD_NONE) alpar@1: { if (csa->parm->msg_lev >= GLP_MSG_ALL) alpar@1: xprintf("Original ordering is being used\n"); alpar@1: for (i = 1; i <= m; i++) alpar@1: csa->P[i] = csa->P[m+i] = i; alpar@1: } alpar@1: else if (csa->parm->ord_alg == GLP_ORD_QMD) alpar@1: { if (csa->parm->msg_lev >= GLP_MSG_ALL) alpar@1: xprintf("Minimum degree ordering (QMD)...\n"); alpar@1: min_degree(m, csa->S_ptr, csa->S_ind, csa->P); alpar@1: } alpar@1: else if (csa->parm->ord_alg == GLP_ORD_AMD) alpar@1: { if (csa->parm->msg_lev >= GLP_MSG_ALL) alpar@1: xprintf("Approximate minimum degree ordering (AMD)...\n"); alpar@1: amd_order1(m, csa->S_ptr, csa->S_ind, csa->P); alpar@1: } alpar@1: else if (csa->parm->ord_alg == GLP_ORD_SYMAMD) alpar@1: { if (csa->parm->msg_lev >= GLP_MSG_ALL) alpar@1: xprintf("Approximate minimum degree ordering (SYMAMD)...\n") alpar@1: ; alpar@1: symamd_ord(m, csa->S_ptr, csa->S_ind, csa->P); alpar@1: } alpar@1: else alpar@1: xassert(csa != csa); alpar@1: /* S := P*A*A'*P', symbolically */ alpar@1: xfree(csa->S_ind); alpar@1: csa->S_ind = adat_symbolic(m, n, csa->P, csa->A_ptr, csa->A_ind, alpar@1: csa->S_ptr); alpar@1: csa->S_val = xcalloc(csa->S_ptr[m+1], sizeof(double)); alpar@1: csa->S_diag = xcalloc(1+m, sizeof(double)); alpar@1: /* compute Cholesky factorization S = U'*U, symbolically */ alpar@1: if (csa->parm->msg_lev >= GLP_MSG_ALL) alpar@1: xprintf("Computing Cholesky factorization S = L*L'...\n"); alpar@1: csa->U_ptr = xcalloc(1+m+1, sizeof(int)); alpar@1: csa->U_ind = chol_symbolic(m, csa->S_ptr, csa->S_ind, csa->U_ptr); alpar@1: if (csa->parm->msg_lev >= GLP_MSG_ALL) alpar@1: xprintf("Matrix L has %d non-zeros\n", csa->U_ptr[m+1]-1 + m); alpar@1: csa->U_val = xcalloc(csa->U_ptr[m+1], sizeof(double)); alpar@1: csa->U_diag = xcalloc(1+m, sizeof(double)); alpar@1: csa->iter = 0; alpar@1: csa->obj = 0.0; alpar@1: csa->rpi = 0.0; alpar@1: csa->rdi = 0.0; alpar@1: csa->gap = 0.0; alpar@1: csa->phi = 0.0; alpar@1: csa->mu = 0.0; alpar@1: csa->rmu = 0.0; alpar@1: csa->rmu0 = 0.0; alpar@1: csa->phi_min = xcalloc(1+ITER_MAX, sizeof(double)); alpar@1: csa->best_iter = 0; alpar@1: csa->best_x = xcalloc(1+n, sizeof(double)); alpar@1: csa->best_y = xcalloc(1+m, sizeof(double)); alpar@1: csa->best_z = xcalloc(1+n, sizeof(double)); alpar@1: csa->best_obj = 0.0; alpar@1: csa->dx_aff = xcalloc(1+n, sizeof(double)); alpar@1: csa->dy_aff = xcalloc(1+m, sizeof(double)); alpar@1: csa->dz_aff = xcalloc(1+n, sizeof(double)); alpar@1: csa->alfa_aff_p = 0.0; alpar@1: csa->alfa_aff_d = 0.0; alpar@1: csa->mu_aff = 0.0; alpar@1: csa->sigma = 0.0; alpar@1: csa->dx_cc = xcalloc(1+n, sizeof(double)); alpar@1: csa->dy_cc = xcalloc(1+m, sizeof(double)); alpar@1: csa->dz_cc = xcalloc(1+n, sizeof(double)); alpar@1: csa->dx = csa->dx_aff; alpar@1: csa->dy = csa->dy_aff; alpar@1: csa->dz = csa->dz_aff; alpar@1: csa->alfa_max_p = 0.0; alpar@1: csa->alfa_max_d = 0.0; alpar@1: return; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * A_by_vec - compute y = A*x alpar@1: * alpar@1: * This routine computes matrix-vector product y = A*x, where A is the alpar@1: * constraint matrix. */ alpar@1: alpar@1: static void A_by_vec(struct csa *csa, double x[], double y[]) alpar@1: { /* compute y = A*x */ alpar@1: int m = csa->m; alpar@1: int *A_ptr = csa->A_ptr; alpar@1: int *A_ind = csa->A_ind; alpar@1: double *A_val = csa->A_val; alpar@1: int i, t, beg, end; alpar@1: double temp; alpar@1: for (i = 1; i <= m; i++) alpar@1: { temp = 0.0; alpar@1: beg = A_ptr[i], end = A_ptr[i+1]; alpar@1: for (t = beg; t < end; t++) temp += A_val[t] * x[A_ind[t]]; alpar@1: y[i] = temp; alpar@1: } alpar@1: return; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * AT_by_vec - compute y = A'*x alpar@1: * alpar@1: * This routine computes matrix-vector product y = A'*x, where A' is a alpar@1: * matrix transposed to the constraint matrix A. */ alpar@1: alpar@1: static void AT_by_vec(struct csa *csa, double x[], double y[]) alpar@1: { /* compute y = A'*x, where A' is transposed to A */ alpar@1: int m = csa->m; alpar@1: int n = csa->n; alpar@1: int *A_ptr = csa->A_ptr; alpar@1: int *A_ind = csa->A_ind; alpar@1: double *A_val = csa->A_val; alpar@1: int i, j, t, beg, end; alpar@1: double temp; alpar@1: for (j = 1; j <= n; j++) y[j] = 0.0; alpar@1: for (i = 1; i <= m; i++) alpar@1: { temp = x[i]; alpar@1: if (temp == 0.0) continue; alpar@1: beg = A_ptr[i], end = A_ptr[i+1]; alpar@1: for (t = beg; t < end; t++) y[A_ind[t]] += A_val[t] * temp; alpar@1: } alpar@1: return; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * decomp_NE - numeric factorization of matrix S = P*A*D*A'*P' alpar@1: * alpar@1: * This routine implements numeric phase of Cholesky factorization of alpar@1: * the matrix S = P*A*D*A'*P', which is a permuted matrix of the normal alpar@1: * equation system. Matrix D is assumed to be already computed. */ alpar@1: alpar@1: static void decomp_NE(struct csa *csa) alpar@1: { adat_numeric(csa->m, csa->n, csa->P, csa->A_ptr, csa->A_ind, alpar@1: csa->A_val, csa->D, csa->S_ptr, csa->S_ind, csa->S_val, alpar@1: csa->S_diag); alpar@1: chol_numeric(csa->m, csa->S_ptr, csa->S_ind, csa->S_val, alpar@1: csa->S_diag, csa->U_ptr, csa->U_ind, csa->U_val, csa->U_diag); alpar@1: return; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * solve_NE - solve normal equation system alpar@1: * alpar@1: * This routine solves the normal equation system: alpar@1: * alpar@1: * A*D*A'*y = h. alpar@1: * alpar@1: * It is assumed that the matrix A*D*A' has been previously factorized alpar@1: * by the routine decomp_NE. alpar@1: * alpar@1: * On entry the array y contains the vector of right-hand sides h. On alpar@1: * exit this array contains the computed vector of unknowns y. alpar@1: * alpar@1: * Once the vector y has been computed the routine checks for numeric alpar@1: * stability. If the residual vector: alpar@1: * alpar@1: * r = A*D*A'*y - h alpar@1: * alpar@1: * is relatively small, the routine returns zero, otherwise non-zero is alpar@1: * returned. */ alpar@1: alpar@1: static int solve_NE(struct csa *csa, double y[]) alpar@1: { int m = csa->m; alpar@1: int n = csa->n; alpar@1: int *P = csa->P; alpar@1: int i, j, ret = 0; alpar@1: double *h, *r, *w; alpar@1: /* save vector of right-hand sides h */ alpar@1: h = xcalloc(1+m, sizeof(double)); alpar@1: for (i = 1; i <= m; i++) h[i] = y[i]; alpar@1: /* solve normal equation system (A*D*A')*y = h */ alpar@1: /* since S = P*A*D*A'*P' = U'*U, then A*D*A' = P'*U'*U*P, so we alpar@1: have inv(A*D*A') = P'*inv(U)*inv(U')*P */ alpar@1: /* w := P*h */ alpar@1: w = xcalloc(1+m, sizeof(double)); alpar@1: for (i = 1; i <= m; i++) w[i] = y[P[i]]; alpar@1: /* w := inv(U')*w */ alpar@1: ut_solve(m, csa->U_ptr, csa->U_ind, csa->U_val, csa->U_diag, w); alpar@1: /* w := inv(U)*w */ alpar@1: u_solve(m, csa->U_ptr, csa->U_ind, csa->U_val, csa->U_diag, w); alpar@1: /* y := P'*w */ alpar@1: for (i = 1; i <= m; i++) y[i] = w[P[m+i]]; alpar@1: xfree(w); alpar@1: /* compute residual vector r = A*D*A'*y - h */ alpar@1: r = xcalloc(1+m, sizeof(double)); alpar@1: /* w := A'*y */ alpar@1: w = xcalloc(1+n, sizeof(double)); alpar@1: AT_by_vec(csa, y, w); alpar@1: /* w := D*w */ alpar@1: for (j = 1; j <= n; j++) w[j] *= csa->D[j]; alpar@1: /* r := A*w */ alpar@1: A_by_vec(csa, w, r); alpar@1: xfree(w); alpar@1: /* r := r - h */ alpar@1: for (i = 1; i <= m; i++) r[i] -= h[i]; alpar@1: /* check for numeric stability */ alpar@1: for (i = 1; i <= m; i++) alpar@1: { if (fabs(r[i]) / (1.0 + fabs(h[i])) > 1e-4) alpar@1: { ret = 1; alpar@1: break; alpar@1: } alpar@1: } alpar@1: xfree(h); alpar@1: xfree(r); alpar@1: return ret; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * solve_NS - solve Newtonian system alpar@1: * alpar@1: * This routine solves the Newtonian system: alpar@1: * alpar@1: * A*dx = p alpar@1: * alpar@1: * A'*dy + dz = q alpar@1: * alpar@1: * Z*dx + X*dz = r alpar@1: * alpar@1: * where X = diag(x[j]), Z = diag(z[j]), by reducing it to the normal alpar@1: * equation system: alpar@1: * alpar@1: * (A*inv(Z)*X*A')*dy = A*inv(Z)*(X*q-r)+p alpar@1: * alpar@1: * (it is assumed that the matrix A*inv(Z)*X*A' has been factorized by alpar@1: * the routine decomp_NE). alpar@1: * alpar@1: * Once vector dy has been computed the routine computes vectors dx and alpar@1: * dz as follows: alpar@1: * alpar@1: * dx = inv(Z)*(X*(A'*dy-q)+r) alpar@1: * alpar@1: * dz = inv(X)*(r-Z*dx) alpar@1: * alpar@1: * The routine solve_NS returns the same code which was reported by the alpar@1: * routine solve_NE (see above). */ alpar@1: alpar@1: static int solve_NS(struct csa *csa, double p[], double q[], double r[], alpar@1: double dx[], double dy[], double dz[]) alpar@1: { int m = csa->m; alpar@1: int n = csa->n; alpar@1: double *x = csa->x; alpar@1: double *z = csa->z; alpar@1: int i, j, ret; alpar@1: double *w = dx; alpar@1: /* compute the vector of right-hand sides A*inv(Z)*(X*q-r)+p for alpar@1: the normal equation system */ alpar@1: for (j = 1; j <= n; j++) alpar@1: w[j] = (x[j] * q[j] - r[j]) / z[j]; alpar@1: A_by_vec(csa, w, dy); alpar@1: for (i = 1; i <= m; i++) dy[i] += p[i]; alpar@1: /* solve the normal equation system to compute vector dy */ alpar@1: ret = solve_NE(csa, dy); alpar@1: /* compute vectors dx and dz */ alpar@1: AT_by_vec(csa, dy, dx); alpar@1: for (j = 1; j <= n; j++) alpar@1: { dx[j] = (x[j] * (dx[j] - q[j]) + r[j]) / z[j]; alpar@1: dz[j] = (r[j] - z[j] * dx[j]) / x[j]; alpar@1: } alpar@1: return ret; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * initial_point - choose initial point using Mehrotra's heuristic alpar@1: * alpar@1: * This routine chooses a starting point using a heuristic proposed in alpar@1: * the paper: alpar@1: * alpar@1: * S. Mehrotra. On the implementation of a primal-dual interior point alpar@1: * method. SIAM J. on Optim., 2(4), pp. 575-601, 1992. alpar@1: * alpar@1: * The starting point x in the primal space is chosen as a solution of alpar@1: * the following least squares problem: alpar@1: * alpar@1: * minimize ||x|| alpar@1: * alpar@1: * subject to A*x = b alpar@1: * alpar@1: * which can be computed explicitly as follows: alpar@1: * alpar@1: * x = A'*inv(A*A')*b alpar@1: * alpar@1: * Similarly, the starting point (y, z) in the dual space is chosen as alpar@1: * a solution of the following least squares problem: alpar@1: * alpar@1: * minimize ||z|| alpar@1: * alpar@1: * subject to A'*y + z = c alpar@1: * alpar@1: * which can be computed explicitly as follows: alpar@1: * alpar@1: * y = inv(A*A')*A*c alpar@1: * alpar@1: * z = c - A'*y alpar@1: * alpar@1: * However, some components of the vectors x and z may be non-positive alpar@1: * or close to zero, so the routine uses a Mehrotra's heuristic to find alpar@1: * a more appropriate starting point. */ alpar@1: alpar@1: static void initial_point(struct csa *csa) alpar@1: { int m = csa->m; alpar@1: int n = csa->n; alpar@1: double *b = csa->b; alpar@1: double *c = csa->c; alpar@1: double *x = csa->x; alpar@1: double *y = csa->y; alpar@1: double *z = csa->z; alpar@1: double *D = csa->D; alpar@1: int i, j; alpar@1: double dp, dd, ex, ez, xz; alpar@1: /* factorize A*A' */ alpar@1: for (j = 1; j <= n; j++) D[j] = 1.0; alpar@1: decomp_NE(csa); alpar@1: /* x~ = A'*inv(A*A')*b */ alpar@1: for (i = 1; i <= m; i++) y[i] = b[i]; alpar@1: solve_NE(csa, y); alpar@1: AT_by_vec(csa, y, x); alpar@1: /* y~ = inv(A*A')*A*c */ alpar@1: A_by_vec(csa, c, y); alpar@1: solve_NE(csa, y); alpar@1: /* z~ = c - A'*y~ */ alpar@1: AT_by_vec(csa, y,z); alpar@1: for (j = 1; j <= n; j++) z[j] = c[j] - z[j]; alpar@1: /* use Mehrotra's heuristic in order to choose more appropriate alpar@1: starting point with positive components of vectors x and z */ alpar@1: dp = dd = 0.0; alpar@1: for (j = 1; j <= n; j++) alpar@1: { if (dp < -1.5 * x[j]) dp = -1.5 * x[j]; alpar@1: if (dd < -1.5 * z[j]) dd = -1.5 * z[j]; alpar@1: } alpar@1: /* note that b = 0 involves x = 0, and c = 0 involves y = 0 and alpar@1: z = 0, so we need to be careful */ alpar@1: if (dp == 0.0) dp = 1.5; alpar@1: if (dd == 0.0) dd = 1.5; alpar@1: ex = ez = xz = 0.0; alpar@1: for (j = 1; j <= n; j++) alpar@1: { ex += (x[j] + dp); alpar@1: ez += (z[j] + dd); alpar@1: xz += (x[j] + dp) * (z[j] + dd); alpar@1: } alpar@1: dp += 0.5 * (xz / ez); alpar@1: dd += 0.5 * (xz / ex); alpar@1: for (j = 1; j <= n; j++) alpar@1: { x[j] += dp; alpar@1: z[j] += dd; alpar@1: xassert(x[j] > 0.0 && z[j] > 0.0); alpar@1: } alpar@1: return; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * basic_info - perform basic computations at the current point alpar@1: * alpar@1: * This routine computes the following quantities at the current point: alpar@1: * alpar@1: * 1) value of the objective function: alpar@1: * alpar@1: * F = c'*x + c[0] alpar@1: * alpar@1: * 2) relative primal infeasibility: alpar@1: * alpar@1: * rpi = ||A*x-b|| / (1+||b||) alpar@1: * alpar@1: * 3) relative dual infeasibility: alpar@1: * alpar@1: * rdi = ||A'*y+z-c|| / (1+||c||) alpar@1: * alpar@1: * 4) primal-dual gap (relative difference between the primal and the alpar@1: * dual objective function values): alpar@1: * alpar@1: * gap = |c'*x-b'*y| / (1+|c'*x|) alpar@1: * alpar@1: * 5) merit function: alpar@1: * alpar@1: * phi = ||A*x-b|| / max(1,||b||) + ||A'*y+z-c|| / max(1,||c||) + alpar@1: * alpar@1: * + |c'*x-b'*y| / max(1,||b||,||c||) alpar@1: * alpar@1: * 6) duality measure: alpar@1: * alpar@1: * mu = x'*z / n alpar@1: * alpar@1: * 7) the ratio of infeasibility to mu: alpar@1: * alpar@1: * rmu = max(||A*x-b||,||A'*y+z-c||) / mu alpar@1: * alpar@1: * where ||*|| denotes euclidian norm, *' denotes transposition. */ alpar@1: alpar@1: static void basic_info(struct csa *csa) alpar@1: { int m = csa->m; alpar@1: int n = csa->n; alpar@1: double *b = csa->b; alpar@1: double *c = csa->c; alpar@1: double *x = csa->x; alpar@1: double *y = csa->y; alpar@1: double *z = csa->z; alpar@1: int i, j; alpar@1: double norm1, bnorm, norm2, cnorm, cx, by, *work, temp; alpar@1: /* compute value of the objective function */ alpar@1: temp = c[0]; alpar@1: for (j = 1; j <= n; j++) temp += c[j] * x[j]; alpar@1: csa->obj = temp; alpar@1: /* norm1 = ||A*x-b|| */ alpar@1: work = xcalloc(1+m, sizeof(double)); alpar@1: A_by_vec(csa, x, work); alpar@1: norm1 = 0.0; alpar@1: for (i = 1; i <= m; i++) alpar@1: norm1 += (work[i] - b[i]) * (work[i] - b[i]); alpar@1: norm1 = sqrt(norm1); alpar@1: xfree(work); alpar@1: /* bnorm = ||b|| */ alpar@1: bnorm = 0.0; alpar@1: for (i = 1; i <= m; i++) bnorm += b[i] * b[i]; alpar@1: bnorm = sqrt(bnorm); alpar@1: /* compute relative primal infeasibility */ alpar@1: csa->rpi = norm1 / (1.0 + bnorm); alpar@1: /* norm2 = ||A'*y+z-c|| */ alpar@1: work = xcalloc(1+n, sizeof(double)); alpar@1: AT_by_vec(csa, y, work); alpar@1: norm2 = 0.0; alpar@1: for (j = 1; j <= n; j++) alpar@1: norm2 += (work[j] + z[j] - c[j]) * (work[j] + z[j] - c[j]); alpar@1: norm2 = sqrt(norm2); alpar@1: xfree(work); alpar@1: /* cnorm = ||c|| */ alpar@1: cnorm = 0.0; alpar@1: for (j = 1; j <= n; j++) cnorm += c[j] * c[j]; alpar@1: cnorm = sqrt(cnorm); alpar@1: /* compute relative dual infeasibility */ alpar@1: csa->rdi = norm2 / (1.0 + cnorm); alpar@1: /* by = b'*y */ alpar@1: by = 0.0; alpar@1: for (i = 1; i <= m; i++) by += b[i] * y[i]; alpar@1: /* cx = c'*x */ alpar@1: cx = 0.0; alpar@1: for (j = 1; j <= n; j++) cx += c[j] * x[j]; alpar@1: /* compute primal-dual gap */ alpar@1: csa->gap = fabs(cx - by) / (1.0 + fabs(cx)); alpar@1: /* compute merit function */ alpar@1: csa->phi = 0.0; alpar@1: csa->phi += norm1 / (bnorm > 1.0 ? bnorm : 1.0); alpar@1: csa->phi += norm2 / (cnorm > 1.0 ? cnorm : 1.0); alpar@1: temp = 1.0; alpar@1: if (temp < bnorm) temp = bnorm; alpar@1: if (temp < cnorm) temp = cnorm; alpar@1: csa->phi += fabs(cx - by) / temp; alpar@1: /* compute duality measure */ alpar@1: temp = 0.0; alpar@1: for (j = 1; j <= n; j++) temp += x[j] * z[j]; alpar@1: csa->mu = temp / (double)n; alpar@1: /* compute the ratio of infeasibility to mu */ alpar@1: csa->rmu = (norm1 > norm2 ? norm1 : norm2) / csa->mu; alpar@1: return; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * make_step - compute next point using Mehrotra's technique alpar@1: * alpar@1: * This routine computes the next point using the predictor-corrector alpar@1: * technique proposed in the paper: alpar@1: * alpar@1: * S. Mehrotra. On the implementation of a primal-dual interior point alpar@1: * method. SIAM J. on Optim., 2(4), pp. 575-601, 1992. alpar@1: * alpar@1: * At first, the routine computes so called affine scaling (predictor) alpar@1: * direction (dx_aff,dy_aff,dz_aff) which is a solution of the system: alpar@1: * alpar@1: * A*dx_aff = b - A*x alpar@1: * alpar@1: * A'*dy_aff + dz_aff = c - A'*y - z alpar@1: * alpar@1: * Z*dx_aff + X*dz_aff = - X*Z*e alpar@1: * alpar@1: * where (x,y,z) is the current point, X = diag(x[j]), Z = diag(z[j]), alpar@1: * e = (1,...,1)'. alpar@1: * alpar@1: * Then, the routine computes the centering parameter sigma, using the alpar@1: * following Mehrotra's heuristic: alpar@1: * alpar@1: * alfa_aff_p = inf{0 <= alfa <= 1 | x+alfa*dx_aff >= 0} alpar@1: * alpar@1: * alfa_aff_d = inf{0 <= alfa <= 1 | z+alfa*dz_aff >= 0} alpar@1: * alpar@1: * mu_aff = (x+alfa_aff_p*dx_aff)'*(z+alfa_aff_d*dz_aff)/n alpar@1: * alpar@1: * sigma = (mu_aff/mu)^3 alpar@1: * alpar@1: * where alfa_aff_p is the maximal stepsize along the affine scaling alpar@1: * direction in the primal space, alfa_aff_d is the maximal stepsize alpar@1: * along the same direction in the dual space. alpar@1: * alpar@1: * After determining sigma the routine computes so called centering alpar@1: * (corrector) direction (dx_cc,dy_cc,dz_cc) which is the solution of alpar@1: * the system: alpar@1: * alpar@1: * A*dx_cc = 0 alpar@1: * alpar@1: * A'*dy_cc + dz_cc = 0 alpar@1: * alpar@1: * Z*dx_cc + X*dz_cc = sigma*mu*e - X*Z*e alpar@1: * alpar@1: * Finally, the routine computes the combined direction alpar@1: * alpar@1: * (dx,dy,dz) = (dx_aff,dy_aff,dz_aff) + (dx_cc,dy_cc,dz_cc) alpar@1: * alpar@1: * and determines maximal primal and dual stepsizes along the combined alpar@1: * direction: alpar@1: * alpar@1: * alfa_max_p = inf{0 <= alfa <= 1 | x+alfa*dx >= 0} alpar@1: * alpar@1: * alfa_max_d = inf{0 <= alfa <= 1 | z+alfa*dz >= 0} alpar@1: * alpar@1: * In order to prevent the next point to be too close to the boundary alpar@1: * of the positive ortant, the routine decreases maximal stepsizes: alpar@1: * alpar@1: * alfa_p = gamma_p * alfa_max_p alpar@1: * alpar@1: * alfa_d = gamma_d * alfa_max_d alpar@1: * alpar@1: * where gamma_p and gamma_d are scaling factors, and computes the next alpar@1: * point: alpar@1: * alpar@1: * x_new = x + alfa_p * dx alpar@1: * alpar@1: * y_new = y + alfa_d * dy alpar@1: * alpar@1: * z_new = z + alfa_d * dz alpar@1: * alpar@1: * which becomes the current point on the next iteration. */ alpar@1: alpar@1: static int make_step(struct csa *csa) alpar@1: { int m = csa->m; alpar@1: int n = csa->n; alpar@1: double *b = csa->b; alpar@1: double *c = csa->c; alpar@1: double *x = csa->x; alpar@1: double *y = csa->y; alpar@1: double *z = csa->z; alpar@1: double *dx_aff = csa->dx_aff; alpar@1: double *dy_aff = csa->dy_aff; alpar@1: double *dz_aff = csa->dz_aff; alpar@1: double *dx_cc = csa->dx_cc; alpar@1: double *dy_cc = csa->dy_cc; alpar@1: double *dz_cc = csa->dz_cc; alpar@1: double *dx = csa->dx; alpar@1: double *dy = csa->dy; alpar@1: double *dz = csa->dz; alpar@1: int i, j, ret = 0; alpar@1: double temp, gamma_p, gamma_d, *p, *q, *r; alpar@1: /* allocate working arrays */ alpar@1: p = xcalloc(1+m, sizeof(double)); alpar@1: q = xcalloc(1+n, sizeof(double)); alpar@1: r = xcalloc(1+n, sizeof(double)); alpar@1: /* p = b - A*x */ alpar@1: A_by_vec(csa, x, p); alpar@1: for (i = 1; i <= m; i++) p[i] = b[i] - p[i]; alpar@1: /* q = c - A'*y - z */ alpar@1: AT_by_vec(csa, y,q); alpar@1: for (j = 1; j <= n; j++) q[j] = c[j] - q[j] - z[j]; alpar@1: /* r = - X * Z * e */ alpar@1: for (j = 1; j <= n; j++) r[j] = - x[j] * z[j]; alpar@1: /* solve the first Newtonian system */ alpar@1: if (solve_NS(csa, p, q, r, dx_aff, dy_aff, dz_aff)) alpar@1: { ret = 1; alpar@1: goto done; alpar@1: } alpar@1: /* alfa_aff_p = inf{0 <= alfa <= 1 | x + alfa*dx_aff >= 0} */ alpar@1: /* alfa_aff_d = inf{0 <= alfa <= 1 | z + alfa*dz_aff >= 0} */ alpar@1: csa->alfa_aff_p = csa->alfa_aff_d = 1.0; alpar@1: for (j = 1; j <= n; j++) alpar@1: { if (dx_aff[j] < 0.0) alpar@1: { temp = - x[j] / dx_aff[j]; alpar@1: if (csa->alfa_aff_p > temp) csa->alfa_aff_p = temp; alpar@1: } alpar@1: if (dz_aff[j] < 0.0) alpar@1: { temp = - z[j] / dz_aff[j]; alpar@1: if (csa->alfa_aff_d > temp) csa->alfa_aff_d = temp; alpar@1: } alpar@1: } alpar@1: /* mu_aff = (x+alfa_aff_p*dx_aff)' * (z+alfa_aff_d*dz_aff) / n */ alpar@1: temp = 0.0; alpar@1: for (j = 1; j <= n; j++) alpar@1: temp += (x[j] + csa->alfa_aff_p * dx_aff[j]) * alpar@1: (z[j] + csa->alfa_aff_d * dz_aff[j]); alpar@1: csa->mu_aff = temp / (double)n; alpar@1: /* sigma = (mu_aff/mu)^3 */ alpar@1: temp = csa->mu_aff / csa->mu; alpar@1: csa->sigma = temp * temp * temp; alpar@1: /* p = 0 */ alpar@1: for (i = 1; i <= m; i++) p[i] = 0.0; alpar@1: /* q = 0 */ alpar@1: for (j = 1; j <= n; j++) q[j] = 0.0; alpar@1: /* r = sigma * mu * e - X * Z * e */ alpar@1: for (j = 1; j <= n; j++) alpar@1: r[j] = csa->sigma * csa->mu - dx_aff[j] * dz_aff[j]; alpar@1: /* solve the second Newtonian system with the same coefficients alpar@1: but with altered right-hand sides */ alpar@1: if (solve_NS(csa, p, q, r, dx_cc, dy_cc, dz_cc)) alpar@1: { ret = 1; alpar@1: goto done; alpar@1: } alpar@1: /* (dx,dy,dz) = (dx_aff,dy_aff,dz_aff) + (dx_cc,dy_cc,dz_cc) */ alpar@1: for (j = 1; j <= n; j++) dx[j] = dx_aff[j] + dx_cc[j]; alpar@1: for (i = 1; i <= m; i++) dy[i] = dy_aff[i] + dy_cc[i]; alpar@1: for (j = 1; j <= n; j++) dz[j] = dz_aff[j] + dz_cc[j]; alpar@1: /* alfa_max_p = inf{0 <= alfa <= 1 | x + alfa*dx >= 0} */ alpar@1: /* alfa_max_d = inf{0 <= alfa <= 1 | z + alfa*dz >= 0} */ alpar@1: csa->alfa_max_p = csa->alfa_max_d = 1.0; alpar@1: for (j = 1; j <= n; j++) alpar@1: { if (dx[j] < 0.0) alpar@1: { temp = - x[j] / dx[j]; alpar@1: if (csa->alfa_max_p > temp) csa->alfa_max_p = temp; alpar@1: } alpar@1: if (dz[j] < 0.0) alpar@1: { temp = - z[j] / dz[j]; alpar@1: if (csa->alfa_max_d > temp) csa->alfa_max_d = temp; alpar@1: } alpar@1: } alpar@1: /* determine scale factors (not implemented yet) */ alpar@1: gamma_p = 0.90; alpar@1: gamma_d = 0.90; alpar@1: /* compute the next point */ alpar@1: for (j = 1; j <= n; j++) alpar@1: { x[j] += gamma_p * csa->alfa_max_p * dx[j]; alpar@1: xassert(x[j] > 0.0); alpar@1: } alpar@1: for (i = 1; i <= m; i++) alpar@1: y[i] += gamma_d * csa->alfa_max_d * dy[i]; alpar@1: for (j = 1; j <= n; j++) alpar@1: { z[j] += gamma_d * csa->alfa_max_d * dz[j]; alpar@1: xassert(z[j] > 0.0); alpar@1: } alpar@1: done: /* free working arrays */ alpar@1: xfree(p); alpar@1: xfree(q); alpar@1: xfree(r); alpar@1: return ret; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * terminate - deallocate common storage area alpar@1: * alpar@1: * This routine frees all memory allocated to the common storage area alpar@1: * used by interior-point method routines. */ alpar@1: alpar@1: static void terminate(struct csa *csa) alpar@1: { xfree(csa->D); alpar@1: xfree(csa->P); alpar@1: xfree(csa->S_ptr); alpar@1: xfree(csa->S_ind); alpar@1: xfree(csa->S_val); alpar@1: xfree(csa->S_diag); alpar@1: xfree(csa->U_ptr); alpar@1: xfree(csa->U_ind); alpar@1: xfree(csa->U_val); alpar@1: xfree(csa->U_diag); alpar@1: xfree(csa->phi_min); alpar@1: xfree(csa->best_x); alpar@1: xfree(csa->best_y); alpar@1: xfree(csa->best_z); alpar@1: xfree(csa->dx_aff); alpar@1: xfree(csa->dy_aff); alpar@1: xfree(csa->dz_aff); alpar@1: xfree(csa->dx_cc); alpar@1: xfree(csa->dy_cc); alpar@1: xfree(csa->dz_cc); alpar@1: return; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * ipm_main - main interior-point method routine alpar@1: * alpar@1: * This is a main routine of the primal-dual interior-point method. alpar@1: * alpar@1: * The routine ipm_main returns one of the following codes: alpar@1: * alpar@1: * 0 - optimal solution found; alpar@1: * 1 - problem has no feasible (primal or dual) solution; alpar@1: * 2 - no convergence; alpar@1: * 3 - iteration limit exceeded; alpar@1: * 4 - numeric instability on solving Newtonian system. alpar@1: * alpar@1: * In case of non-zero return code the routine returns the best point, alpar@1: * which has been reached during optimization. */ alpar@1: alpar@1: static int ipm_main(struct csa *csa) alpar@1: { int m = csa->m; alpar@1: int n = csa->n; alpar@1: int i, j, status; alpar@1: double temp; alpar@1: /* choose initial point using Mehrotra's heuristic */ alpar@1: if (csa->parm->msg_lev >= GLP_MSG_ALL) alpar@1: xprintf("Guessing initial point...\n"); alpar@1: initial_point(csa); alpar@1: /* main loop starts here */ alpar@1: if (csa->parm->msg_lev >= GLP_MSG_ALL) alpar@1: xprintf("Optimization begins...\n"); alpar@1: for (;;) alpar@1: { /* perform basic computations at the current point */ alpar@1: basic_info(csa); alpar@1: /* save initial value of rmu */ alpar@1: if (csa->iter == 0) csa->rmu0 = csa->rmu; alpar@1: /* accumulate values of min(phi[k]) and save the best point */ alpar@1: xassert(csa->iter <= ITER_MAX); alpar@1: if (csa->iter == 0 || csa->phi_min[csa->iter-1] > csa->phi) alpar@1: { csa->phi_min[csa->iter] = csa->phi; alpar@1: csa->best_iter = csa->iter; alpar@1: for (j = 1; j <= n; j++) csa->best_x[j] = csa->x[j]; alpar@1: for (i = 1; i <= m; i++) csa->best_y[i] = csa->y[i]; alpar@1: for (j = 1; j <= n; j++) csa->best_z[j] = csa->z[j]; alpar@1: csa->best_obj = csa->obj; alpar@1: } alpar@1: else alpar@1: csa->phi_min[csa->iter] = csa->phi_min[csa->iter-1]; alpar@1: /* display information at the current point */ alpar@1: if (csa->parm->msg_lev >= GLP_MSG_ON) alpar@1: xprintf("%3d: obj = %17.9e; rpi = %8.1e; rdi = %8.1e; gap =" alpar@1: " %8.1e\n", csa->iter, csa->obj, csa->rpi, csa->rdi, alpar@1: csa->gap); alpar@1: /* check if the current point is optimal */ alpar@1: if (csa->rpi < 1e-8 && csa->rdi < 1e-8 && csa->gap < 1e-8) alpar@1: { if (csa->parm->msg_lev >= GLP_MSG_ALL) alpar@1: xprintf("OPTIMAL SOLUTION FOUND\n"); alpar@1: status = 0; alpar@1: break; alpar@1: } alpar@1: /* check if the problem has no feasible solution */ alpar@1: temp = 1e5 * csa->phi_min[csa->iter]; alpar@1: if (temp < 1e-8) temp = 1e-8; alpar@1: if (csa->phi >= temp) alpar@1: { if (csa->parm->msg_lev >= GLP_MSG_ALL) alpar@1: xprintf("PROBLEM HAS NO FEASIBLE PRIMAL/DUAL SOLUTION\n") alpar@1: ; alpar@1: status = 1; alpar@1: break; alpar@1: } alpar@1: /* check for very slow convergence or divergence */ alpar@1: if (((csa->rpi >= 1e-8 || csa->rdi >= 1e-8) && csa->rmu / alpar@1: csa->rmu0 >= 1e6) || alpar@1: (csa->iter >= 30 && csa->phi_min[csa->iter] >= 0.5 * alpar@1: csa->phi_min[csa->iter - 30])) alpar@1: { if (csa->parm->msg_lev >= GLP_MSG_ALL) alpar@1: xprintf("NO CONVERGENCE; SEARCH TERMINATED\n"); alpar@1: status = 2; alpar@1: break; alpar@1: } alpar@1: /* check for maximal number of iterations */ alpar@1: if (csa->iter == ITER_MAX) alpar@1: { if (csa->parm->msg_lev >= GLP_MSG_ALL) alpar@1: xprintf("ITERATION LIMIT EXCEEDED; SEARCH TERMINATED\n"); alpar@1: status = 3; alpar@1: break; alpar@1: } alpar@1: /* start the next iteration */ alpar@1: csa->iter++; alpar@1: /* factorize normal equation system */ alpar@1: for (j = 1; j <= n; j++) csa->D[j] = csa->x[j] / csa->z[j]; alpar@1: decomp_NE(csa); alpar@1: /* compute the next point using Mehrotra's predictor-corrector alpar@1: technique */ alpar@1: if (make_step(csa)) alpar@1: { if (csa->parm->msg_lev >= GLP_MSG_ALL) alpar@1: xprintf("NUMERIC INSTABILITY; SEARCH TERMINATED\n"); alpar@1: status = 4; alpar@1: break; alpar@1: } alpar@1: } alpar@1: /* restore the best point */ alpar@1: if (status != 0) alpar@1: { for (j = 1; j <= n; j++) csa->x[j] = csa->best_x[j]; alpar@1: for (i = 1; i <= m; i++) csa->y[i] = csa->best_y[i]; alpar@1: for (j = 1; j <= n; j++) csa->z[j] = csa->best_z[j]; alpar@1: if (csa->parm->msg_lev >= GLP_MSG_ALL) alpar@1: xprintf("Best point %17.9e was reached on iteration %d\n", alpar@1: csa->best_obj, csa->best_iter); alpar@1: } alpar@1: /* return to the calling program */ alpar@1: return status; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * NAME alpar@1: * alpar@1: * ipm_solve - core LP solver based on the interior-point method alpar@1: * alpar@1: * SYNOPSIS alpar@1: * alpar@1: * #include "glpipm.h" alpar@1: * int ipm_solve(glp_prob *P, const glp_iptcp *parm); alpar@1: * alpar@1: * DESCRIPTION alpar@1: * alpar@1: * The routine ipm_solve is a core LP solver based on the primal-dual alpar@1: * interior-point method. alpar@1: * alpar@1: * The routine assumes the following standard formulation of LP problem alpar@1: * to be solved: alpar@1: * alpar@1: * minimize alpar@1: * alpar@1: * F = c[0] + c[1]*x[1] + c[2]*x[2] + ... + c[n]*x[n] alpar@1: * alpar@1: * subject to linear constraints alpar@1: * alpar@1: * a[1,1]*x[1] + a[1,2]*x[2] + ... + a[1,n]*x[n] = b[1] alpar@1: * alpar@1: * a[2,1]*x[1] + a[2,2]*x[2] + ... + a[2,n]*x[n] = b[2] alpar@1: * alpar@1: * . . . . . . alpar@1: * alpar@1: * a[m,1]*x[1] + a[m,2]*x[2] + ... + a[m,n]*x[n] = b[m] alpar@1: * alpar@1: * and non-negative variables alpar@1: * alpar@1: * x[1] >= 0, x[2] >= 0, ..., x[n] >= 0 alpar@1: * alpar@1: * where: alpar@1: * F is the objective function; alpar@1: * x[1], ..., x[n] are (structural) variables; alpar@1: * c[0] is a constant term of the objective function; alpar@1: * c[1], ..., c[n] are objective coefficients; alpar@1: * a[1,1], ..., a[m,n] are constraint coefficients; alpar@1: * b[1], ..., b[n] are right-hand sides. alpar@1: * alpar@1: * The solution is three vectors x, y, and z, which are stored by the alpar@1: * routine in the arrays x, y, and z, respectively. These vectors alpar@1: * correspond to the best primal-dual point found during optimization. alpar@1: * They are approximate solution of the following system (which is the alpar@1: * Karush-Kuhn-Tucker optimality conditions): alpar@1: * alpar@1: * A*x = b (primal feasibility condition) alpar@1: * alpar@1: * A'*y + z = c (dual feasibility condition) alpar@1: * alpar@1: * x'*z = 0 (primal-dual complementarity condition) alpar@1: * alpar@1: * x >= 0, z >= 0 (non-negativity condition) alpar@1: * alpar@1: * where: alpar@1: * x[1], ..., x[n] are primal (structural) variables; alpar@1: * y[1], ..., y[m] are dual variables (Lagrange multipliers) for alpar@1: * equality constraints; alpar@1: * z[1], ..., z[n] are dual variables (Lagrange multipliers) for alpar@1: * non-negativity constraints. alpar@1: * alpar@1: * RETURNS alpar@1: * alpar@1: * 0 LP has been successfully solved. alpar@1: * alpar@1: * GLP_ENOCVG alpar@1: * No convergence. alpar@1: * alpar@1: * GLP_EITLIM alpar@1: * Iteration limit exceeded. alpar@1: * alpar@1: * GLP_EINSTAB alpar@1: * Numeric instability on solving Newtonian system. alpar@1: * alpar@1: * In case of non-zero return code the routine returns the best point, alpar@1: * which has been reached during optimization. */ alpar@1: alpar@1: int ipm_solve(glp_prob *P, const glp_iptcp *parm) alpar@1: { struct csa _dsa, *csa = &_dsa; alpar@1: int m = P->m; alpar@1: int n = P->n; alpar@1: int nnz = P->nnz; alpar@1: GLPROW *row; alpar@1: GLPCOL *col; alpar@1: GLPAIJ *aij; alpar@1: int i, j, loc, ret, *A_ind, *A_ptr; alpar@1: double dir, *A_val, *b, *c, *x, *y, *z; alpar@1: xassert(m > 0); alpar@1: xassert(n > 0); alpar@1: /* allocate working arrays */ alpar@1: A_ptr = xcalloc(1+m+1, sizeof(int)); alpar@1: A_ind = xcalloc(1+nnz, sizeof(int)); alpar@1: A_val = xcalloc(1+nnz, sizeof(double)); alpar@1: b = xcalloc(1+m, sizeof(double)); alpar@1: c = xcalloc(1+n, sizeof(double)); alpar@1: x = xcalloc(1+n, sizeof(double)); alpar@1: y = xcalloc(1+m, sizeof(double)); alpar@1: z = xcalloc(1+n, sizeof(double)); alpar@1: /* prepare rows and constraint coefficients */ alpar@1: loc = 1; alpar@1: for (i = 1; i <= m; i++) alpar@1: { row = P->row[i]; alpar@1: xassert(row->type == GLP_FX); alpar@1: b[i] = row->lb * row->rii; alpar@1: A_ptr[i] = loc; alpar@1: for (aij = row->ptr; aij != NULL; aij = aij->r_next) alpar@1: { A_ind[loc] = aij->col->j; alpar@1: A_val[loc] = row->rii * aij->val * aij->col->sjj; alpar@1: loc++; alpar@1: } alpar@1: } alpar@1: A_ptr[m+1] = loc; alpar@1: xassert(loc-1 == nnz); alpar@1: /* prepare columns and objective coefficients */ alpar@1: if (P->dir == GLP_MIN) alpar@1: dir = +1.0; alpar@1: else if (P->dir == GLP_MAX) alpar@1: dir = -1.0; alpar@1: else alpar@1: xassert(P != P); alpar@1: c[0] = dir * P->c0; alpar@1: for (j = 1; j <= n; j++) alpar@1: { col = P->col[j]; alpar@1: xassert(col->type == GLP_LO && col->lb == 0.0); alpar@1: c[j] = dir * col->coef * col->sjj; alpar@1: } alpar@1: /* allocate and initialize the common storage area */ alpar@1: csa->m = m; alpar@1: csa->n = n; alpar@1: csa->A_ptr = A_ptr; alpar@1: csa->A_ind = A_ind; alpar@1: csa->A_val = A_val; alpar@1: csa->b = b; alpar@1: csa->c = c; alpar@1: csa->x = x; alpar@1: csa->y = y; alpar@1: csa->z = z; alpar@1: csa->parm = parm; alpar@1: initialize(csa); alpar@1: /* solve LP with the interior-point method */ alpar@1: ret = ipm_main(csa); alpar@1: /* deallocate the common storage area */ alpar@1: terminate(csa); alpar@1: /* determine solution status */ alpar@1: if (ret == 0) alpar@1: { /* optimal solution found */ alpar@1: P->ipt_stat = GLP_OPT; alpar@1: ret = 0; alpar@1: } alpar@1: else if (ret == 1) alpar@1: { /* problem has no feasible (primal or dual) solution */ alpar@1: P->ipt_stat = GLP_NOFEAS; alpar@1: ret = 0; alpar@1: } alpar@1: else if (ret == 2) alpar@1: { /* no convergence */ alpar@1: P->ipt_stat = GLP_INFEAS; alpar@1: ret = GLP_ENOCVG; alpar@1: } alpar@1: else if (ret == 3) alpar@1: { /* iteration limit exceeded */ alpar@1: P->ipt_stat = GLP_INFEAS; alpar@1: ret = GLP_EITLIM; alpar@1: } alpar@1: else if (ret == 4) alpar@1: { /* numeric instability on solving Newtonian system */ alpar@1: P->ipt_stat = GLP_INFEAS; alpar@1: ret = GLP_EINSTAB; alpar@1: } alpar@1: else alpar@1: xassert(ret != ret); alpar@1: /* store row solution components */ alpar@1: for (i = 1; i <= m; i++) alpar@1: { row = P->row[i]; alpar@1: row->pval = row->lb; alpar@1: row->dval = dir * y[i] * row->rii; alpar@1: } alpar@1: /* store column solution components */ alpar@1: P->ipt_obj = P->c0; alpar@1: for (j = 1; j <= n; j++) alpar@1: { col = P->col[j]; alpar@1: col->pval = x[j] * col->sjj; alpar@1: col->dval = dir * z[j] / col->sjj; alpar@1: P->ipt_obj += col->coef * col->pval; alpar@1: } alpar@1: /* free working arrays */ alpar@1: xfree(A_ptr); alpar@1: xfree(A_ind); alpar@1: xfree(A_val); alpar@1: xfree(b); alpar@1: xfree(c); alpar@1: xfree(x); alpar@1: xfree(y); alpar@1: xfree(z); alpar@1: return ret; alpar@1: } alpar@1: alpar@1: /* eof */