TOC | Abstract | (1) Intro | (2) Simplex | (3) Impl | (4) Test | (5) Results | (6) Concl | Scripts | Biblio | Author | Chairman |
Forrest-Tomlin Method
function solution = ft (c, A, b, eps1, eps2, eps3, bfs) % % Solves: minimize cx subject to Ax <= b & x >= 0 % m number of rows in A % n number of columns in A % B_indices vector of columns in A comprising the solution basis % V_indices vector of columns in A not in solution basis [m n] = size(A); B_indices = find(bfs); V_indices = find(ones(1,n) - abs(sign(bfs))); ft_nnz = zeros(5000,2); refactors=0; iters=0; while (1 == 1) refactors=refactors + 1; % L lower triangular factor of the basis % U upper triangular factor of the basis % R matrix of R vectors % numR number of entries in R % Q column permutation matrix B_indices = B_indices(:,colmmd(A(:,B_indices))); [L U pt] = lu(A(:,B_indices)); P = [1:m] * pt'; Q=[1:m]; R = zeros(m, m+1); numR = 0; % Simplex method loops continuously until solution is found or discovered % to be impossible. ok = 1; while (nnz(R) < (.65 * m) & ok) iters=iters+1; ft_nnz(iters,1) = nnz(A(:,B_indices)); ft_nnz(iters,2) = nnz(L) + nnz(U) + nnz(R) - numR; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Step 1 % compute B^-1 % We won't explicitly calculate it here, but it will be solved when % necessary to take advantage of matrix structures: % B = L * R(1) * ... * R(t) * U * Q' % B^-1 = Q * U^-1 * Q * R(t)^-1 * ... * R(1)^-1 * L^-1 % Qinv reverse column permutation of Q Qinv(Q) = [1:m]; Pinv(P) = [1:m]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Step 2 % compute d = B^-1 * b = B \ b % dtemp accumulating value of d during R computation % z loop variable % d current solution vector dtemp = L \ b(P); for z = 1:numR rtemp = [-R(z, 2:R(z,1)) 1 -R(z,R(z,1)+2:m+1)]; dtemp(R(z,1)) = rtemp * dtemp; end; %while dtemp = U(Q,:) \ dtemp(Q,:); d = dtemp(Qinv,:); if (norm(A(:,B_indices)*d - b, inf) > eps3) ok = 0; end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Step 3/Step 4/Step 5 % compute c_tilde = c_V - c_B * B^-1 * V % ctemp accumulating value during R computation % c_tilde modified cost vector c_tilde = zeros(1,n); ctemp = L \ A(P,V_indices); for z = 1:numR rtemp = [-R(z, 2:R(z,1)) 1 -R(z,R(z,1)+2:m+1)]; ctemp(R(z,1),:) = rtemp * ctemp; end; %while ctemp = U(Q,:) \ ctemp(Q,:); c_tilde(:,V_indices) = c(:,V_indices) - c(:,B_indices) * ctemp(Qinv,:); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Step 6 % compute j s.t. c_tilde[j] <= c_tilde[k] for all k in V_indices % cj minimum cost value (negative) of non-basic columns % j column in A corresponding to minimum cost value [cj j] = min(c_tilde); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Step 7 % if cj >= 0 , then we're done -- return solution which is optimal if cj >= -eps1 solution = zeros(n,1); solution(B_indices,:) = d; return; end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Step 8 % otherwise, compute w = B^-1 * a[j] % wtemp accumulating value of w during R computation % w relative weight (vector) of column entering the basis wtemp = L \ A(P,j); for z = 1:numR rtemp = [-R(z, 2:R(z,1)) 1 -R(z,R(z,1)+2:m+1)]; wtemp(R(z,1)) = rtemp * wtemp; end; %while wtemp = U(Q,:) \ wtemp(Q,:); w = wtemp(Qinv,:); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Step 9 % compute i s.t. w[i]>0 and d[i]/w[i] is minimized % mn minimum of d[i]/w[i] when w[i] > 0 % i row corresponding to mn -- detemines outgoing column % k temporary storage variable mn = inf; i=0; zz = find (w > eps1)' ; [yy, ii] = min (d(zz) ./ w (zz)) ; i = zz(ii(1)) ; k = B_indices(i); B_indices(i) = j; V_indices(j == V_indices) = k; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%% Update U %%%%%%% ** Forrest-Tomlin ** %%%%%%%%%%%%%%%%%%%%%%%% %% U is a row-permuted upper-triangular matrix %% Q*U is the proper upper-triangular matrix % uc spiked column of U resulting from new basis column % ur row corresponding to uc uc = Qinv(i); ur = Q(uc); % r new entry in R that eliminates the nonzero values % in ur r = [zeros(1,uc) U(ur, uc+1:m)]; r = r / U; numR = numR+1; R(numR,:) = [ur r]; % We already know the effects of r, so there is no need to explicitly % calculate them. if (uc < m) U(ur,uc+1:m) = zeros(1,m-uc); end; % Update the incoming column according to the operations already done on U incoming = L \ A(P,j); for z = 1:numR rtemp = [-R(z, 2:R(z,1)) 1 -R(z,R(z,1)+2:m+1)]; incoming(R(z,1)) = rtemp * incoming; end; %while if (abs(incoming(ur) < eps2)) ok = 0; end; U = [U(:,[1:uc-1 uc+1:m]) incoming]; % And then update Q Q = Q(:,[1:uc-1 uc+1:m uc]); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Step 10 % REPEAT end; % while end; % while