/*******************  Quadratic Programming  ************************
Author: Rob Dittmar     dittmar@stls.frb.org
        Research Department
        Federal Reserve Bank of St. Louis
Date:   20 Sep 1996
This code is provided free, without performance guarantees,
for public, non-commercial use.


To solve a quadratic programming problem, you usually write the first order
conditions in the form of a what is referred to as a linear complementarity
problem and use a variant of the simplex algorithm to solve it.

If the quadratic programming problem to be solved is of the form

        Min c'x + x'Qx
        s. t. A x >= b
                x >= 0

then the first order conditions take the form

      v = 2 Q x - A' mu + c
      s = A x - b
      x, mu, v, s >= 0
      x'v + s'mu = 0

where v are the multipliers on the non-negativity constraints, the mus are
multipliers on the inequality constraints, and the s are slack variables for
the inequality contraints.

This problem is in the linear-complementarity problem form which is to find
vectors w and z to solve

     w = M z + q
     w, z, >= 0
     w'z = 0.

I have attached a short procedure that will solve the linear-complementarity
problem.  Write the FOC's for your quadratic problem in the above form and
solve the associated linear complementarity problem.
****************************************************************************/


/* Smallest allowable pivot element */
_pivtol = 1e-8;

/* Procedure to perform partial pivoting on a matrix.  The
   procedure takes a matrix argument, tableau, and two scalar
   arguments, r and c.  The procedure performs a row reduction
   using the [r,c]-element of the matrix as a pivot.  The returned
   matrix will have a 1 as its [r,c]-element, and zeros everywhere
   else in column c. */
PROC (1) = rowred(tableau,r,c);

LOCAL reduced,pivot,i;

reduced = tableau;
pivot = tableau[r,.]/tableau[r,c];

i = 1;
DO WHILE i <= ROWS(tableau);
  reduced[i,.] = tableau[i,.] - tableau[i,c]*pivot;
  i = i + 1;
ENDO;

reduced[r,.] = pivot;

RETP(reduced);

ENDP;

/* Procedure to solve the linear complementarity problem:

       w = M z + q
       w and z >= 0
       w'z = 0

   The procedure takes the matrix M and vector q as arguments.  The
   procedure has three returns.  The first and second returns are
   the final values of the vectors w and z found by complementary
   pivoting.  The third return is a 2 by 1 vector.  Its first
   component is a 1 if the algorithm was successful, and a 2 if a
   ray termination resulted.  The second component is the value of
   the artificial varible upon termination of the algorithm. */
PROC (3) = LCPSolve(M,q);

LOCAL dimen,tableau,basis,locat,cand,rayTerm,eMs,quots,
  oldVar,vars,w,z,retcode;

/* Dimension of LCP */
dimen = ROWS(M);

/* Create initial tableau */
tableau = EYE(dimen)~(-M)~(-ONES(dimen,1))~q;

/* Let artificial variable enter the basis */
basis = SEQA(1,1,dimen);
locat = MININDC(tableau[.,2*dimen+2]);
basis[locat] = 2*dimen + 1;
cand = locat + dimen;
tableau = rowred(tableau,locat,2*dimen+1);

/* Perform complementary pivoting */
rayTerm = 0;
DO WHILE MAXC(basis) == 2*dimen+1 AND NOT(rayTerm);
  eMs = tableau[.,cand].*MISS(tableau[.,cand] .> 0, 0);
  quots = tableau[.,2*dimen+2]./eMs;
  locat = MININDC(quots);
  IF abs(eMs[locat]) > _pivtol AND NOT ISMISS(PACKR(quots));
    @ Reduce tableau @
    tableau = rowred(tableau,locat,cand);
    oldVar = basis[locat];
    @ New variable enters the basis @
    basis[locat] = cand;
    @ Select next candidate for entering the basis @
    IF oldVar > dimen;
      cand = oldVar - dimen;
    ELSE;
      cand = oldVar + dimen;
    ENDIF;
  ELSE;
    rayTerm = 1;
  ENDIF;
ENDO;

@ Return solution to LCP @
vars = ZEROS(2*dimen+1,1);
vars[basis] = tableau[.,2*dimen+2];

w = vars[1:dimen];
z = vars[(dimen+1):2*dimen];

retcode = vars[(2*dimen+1)];
IF rayTerm;
  retcode = 2|retcode;
ELSE;
  retcode = 1|retcode;
ENDIF;

RETP(w,z,retcode);

ENDP;

