! LINGO MODEL for interior point/barrier primal-dual method. (IntPtAlg.lg4)
   Original version by Kipp Martin. See article by
    Marsten, et. al.(1990) in INTERFACES, vol. 20, no. 4.;
 ! Keywords: Interior Point Method, Barrier Method;
SETS:
   ROW  : B, Y, DY;
   COL  : C, DX, X, DZ, Z;
   MATRIX( ROW, COL) : A;
ENDSETS

 ! Solve the optimization problem 
     Min CX
      s.t.
      AX = B,
       X >= 0,
 by solving approximately the system:
   1)  AX = B,     (Primal conditions)
   2)  A'Y - C = Z,(Dual conditions)
   3)  X >= 0,     (Non-negativity)
   4)  Z >= 0
   5)  ZX = 0.     (Complementarity)
 We start with a solution to the linear constraints, (1-4), and
 then use iterative approximation to satisfy (5);
  DATA:
   COL = A__ C__ D__ S_1 S_2 S_3;
   B = 60, 50, 120;       ! RHS;
   C = -20, -30, -47,  0,  0,  0; ! Objective coefficients, cost minimize form;
   A =   1,   0,   1,  1,  0,  0, ! Coefficient matrix for the;
         0,   1,   1,  0,  1,  0, ! Astro/Cosmo;
         1,   2,   3,  0,  0,  1; ! Problem in Schrage;

   X =   1    1    1  58  48 114; ! Initial feasible primal soln;
   Z =  40   60  103  30  30  30; ! Initial feasible reduced costs;
   Y =  30 30 30;      ! Initial dual prices;
   TOLZR  = .0000001;  ! Zero tolerance;
   TOLGAP = .00000005; ! Termination tolerance;
  ENDDATA

 SUBMODEL DIRECTION:
  !  Given a current point:
       X : feasible (nonnegative) primal point,
       Y : feasible dual prices,
       Z : feasible (nonnegative) reduced costs,
    Choose a direction: DX, DY and DZ that
     tries to better satisfy the complementary slackness constraint;

  ! The Primal equations: A*( X + DX) = B;
  @FOR( ROW( I):
   [P] @SUM( COL( J) : A( I, J)*( X( J) + DX( J))) = B( I);
       );

  ! The Dual equations: A*( Y + DY) -(Z + DZ) = -C;
  @FOR( COL( J):
   [D] @SUM( ROW( I): A( I, J) * ( Y( I) + DY( I))) -
    (Z( J) + DZ( J)) = -C( J)
      );

  ! The Complementary slackness equation:
    ( Z + DZ) * ( X + DX) =  0, or
    Z*X + Z*DX + DZ*X + DZ*DX = 0, 
    or the linear approximation
    Z*X + Z*DX + DZ*X ~ 0, ;
  ! and the linear approximation to it;
  @FOR( COL( J):
   [CM] ( Z( J) + DZ( J)) * X( J) + Z( J) * DX( J) = MU
       );

  ! Allow for movements in the negative direction;
  @FOR( ROW: @FREE( DY));
  @FOR( COL: @FREE( DX));
  @FOR( COL: @FREE( DZ));
 ENDSUBMODEL


 CALC:
   N = @SIZE( COL);
   ! Turn off all messages;
   @SET('TERSEO',3); ! Set to <= 2 for debugging;
   maxiter = 30; ! Seldom need more than this # iterations;
   more = 1;     ! Set to 0 when done;
   iter = 0;
   @WHILE( iter #LE# maxiter #AND# more:
    iter = iter + 1; 
    gap = @SUM( COL: X * Z); ! Duality gap;
    @WRITE(iter,') Current gap= ', gap,@NEWLINE(1));
! Marsten's rule for error target, MU;
    MU = gap / ( N * N);
    @SOLVE(  DIRECTION); ! Get the move direction;

   ! Calculate step size in the move direction so that
     X and Z remain > 0;
   !  Do not let any variable get too close to zero;
    ALPHAP = .995 * ( @MIN(COL( J) | DX( J) #LT# -TOLZR:
      -X(J) / DX(J)));
    ALPHAD = .995 * ( @MIN(COL( J) | DZ( J) #LT# -TOLZR:
      -Z(J) / DZ(J)));

   ! Move to the new point;
    @FOR( COL( J): X( J) = X( J) + ALPHAP * DX( J));
    @FOR( COL( J): Z( J) = Z( J) + ALPHAD * DZ( J));
    @FOR( ROW( I): Y( I) = Y( I) + ALPHAD * DY( I));

   ! Calculate primal and dual objectives;
    POBJ = -@SUM( COL: X * C);
    DOBJ = @SUM( ROW: Y * B);
    gap = DOBJ - POBJ ;

  ! Write out details of new point;
    @WRITE(' New Primal, Dual, Gap= ', @FORMAT(POBJ,'12.4f'),' ',@FORMAT(DOBJ,'12.4f'),' ',gap,@NEWLINE(1));
    @FOR( COL(j): @WRITE('        ',COL(j)));
    @WRITE(@NEWLINE(1),'X= ');
    @FOR( COL(j):
      @WRITE( @FORMAT(X(j),'10.3f'),' ');
        );
    @WRITE(@NEWLINE(1));
    @WRITE('Z= ');
    @FOR( COL(j):
      @WRITE( @FORMAT(Z(j),'10.3f'),' ');
        );
    @WRITE(@NEWLINE(1));
    @WRITE('Y= ');
    @FOR( ROW(i):
      @WRITE( @FORMAT(Y(i),'10.3f'),' ');
        );
    @WRITE(@NEWLINE(1));

    ! Exit if gap is small enough;
    @IFC( gap #LE# TOLGAP *( 1 + @ABS( @SUM( COL: C * X))): more = 0); 
    );
 ENDCALC