#  MAPLE.  OCtober 7, 2005.
#
#  Implement the affine scaling method.
#
#  Input: A  - an m x n matrix with full row rank.
#         b  - an m x 1 matrix 
#         c  - and n x 1 matrix
#         x0 - a strictly positive solution to Ax=b
#         r  - a real number  0 < r < 1, the step size.
#
#  procedure "affineit" performs one iteration of the affine scaling method.
#  Output: x1 - an n x 1 matrix which is the next strictly positive 
#               feasible solution after x0 in the affine scaling method.
#  The code prints most intermediate calculations to stdout
#  Be careul when using sqrt(.), etc. since MAPLE cannot always tell
#  if entries are non-zero. It is simplest to use floating point values.

with(linalg):




affineit := proc( oA, b, c, x0, r)
local  A,cc,  n,m,x,X,i,tA,P,d,beta,j,ones;
# Assume: Given A, b, c with A (m by n) having full row rank and 
#   a strictly positive solution x0 (n by 1).
n := coldim(oA);
m := rowdim(oA);

x := x0;


   print(`---------------------------------------------------`);
   print(`Beginning of iteration`);
   print(`---------------------------------------------------`);

   X := matrix(n,n,0):
   for i to n do X[i,i] := x[i,1]; od:
   print(`   X = `,op(X) );
   A := multiply(oA,X);
   cc := multiply(X,c);
   print(`   A-hat = `,op(A) ,`   c-hat = `,op(cc) );
   tA := transpose(A);
   P := evalm( 1 - tA &* (A &* tA)^(-1) &* A );
   print(` A * A^T = `,evalm( A &* tA) );
   print(` (A * A^T)^(-1) = `,inverse( evalm( A &* tA) )   );
   print(` A^T * (A * A^T)^(-1) = `, evalm( tA &* (A &* tA)^(-1) ) );
   print(` A^T * (A * A^T)^(-1) * A = `, evalm( tA &* (A &* tA)^(-1) &* A) );
   print(` P = `, op(P) );
   d := multiply(P,cc);
   print( ` d = `, op(d) );
   beta := infinity;
   for j to n do if d[j,1] < 0 then beta := min( beta, -1/d[j,1] ); fi od:
   print(` beta = `,beta);
   ones := matrix(n,1,1);  # all ones vector
   x := evalm( ones + r*beta*d );
   print(` new scaled solution x-hat = X x = `,op(x) );
   x := multiply( X, x);
   print(` new solution  x = `, op(x) ,` objective value zeta =`, 
      multiply( transpose(c), x)  );

RETURN( x );
end;