disp('WARNING: please change the value of pt if you want pauses between iterations')
disp(' ')
pt=0;  % no pause in computations
disp('This matlab file solves assignment 5, fall 1994, co350')
disp('Note that the matlab program minimizes - adjust solution accordingly.')
disp('The calculations for the dual simplex method are not shown.')
probno=input('enter desired problem number, 1 or 2 or 3, ')
%%%%%%%%%%%%%%%%%%%%%%%prob 1%%%%%%%%%%%%%%%%%%%
if probno==1,
  co=-[-2 -1 -2 3];
  ao=[
  8 -4 -1 3
  2 3 1 -1 
  0 -2 -1 4
  ];
  [m,n]=size(ao);
  bo=[10 7 12]';
  m1=3;m2=0;
  a=ao;b=bo;c=co;
  simplex
  disp('Original part of prob 1 solved.')
  disp('Now solve the desired parts.')
  pause(1)
  disp('The numbers 1,2,3,4 correspond to parts a,b,c,d of problem 1.')
  disp('Enter the parts of prob 1 to solve as a vector,') 
  parti=input('e.g. 1 or  [1,3] or [1,2,3,4]  :  ')
 for i=1:length(parti),
  part=parti(i);
    if part==1,
      disp('for Part a of Prob 1:')
      disp('for the min problem c_1 can vary from -inf to c_1 + delta')
      disp('delta is c_1bar from the final tableau')
      disp('c_1 and delta are:')
      [-co(1) c(1)]
      'hit enter to continue'
      pause
    elseif part==2,
      disp('for Part b of Prob 1:')
      aoo=[ao eye(m)];   %  original matrix after adding slacks
      ab=aoo(:,bas);     %  basic columns at optimum
      an=aoo(:,nbas);    %  nonbasic columns at optimum
      abi=inv(ab);       %  inverse basis matrix (can also be obtained 
                            %  from final tableau)
      cv=abi(2,:)*an;    %  change in cost coefficient
      disp('following are the two vectors to check ratios')
      disp('cnbar')
      [c(nbas)]
      disp('second row of (inverse of basis matrix) times nonbasic cols')
      [cv]
      disp('Since cv > 0 and')
      disp('since we are minimizing, to maintain optimality:')
      disp('cnbar - delta cv  >= 0 so take ratios: delta <= min(cnbar./cv)')
      min(c(nbas)./cv)
      'hit enter to continue'
      pause
    elseif part==3,
      disp('for Part c of Prob 1:')
      deltab=[-10 0 0]';     % change in the rhs vector
      aoo=[ao eye(m)];   %  original matrix after adding slacks
      ab=aoo(:,bas);     %  basic columns at optimum
      an=aoo(:,nbas);    %  nonbasic columns at optimum
      abi=inv(ab);       %  inverse basis matrix (can also be obtained 
                            %  from final tableau)
      deltax=abi*deltab;  % change in optimum
      xnew=b+deltax;      %new optimum
      disp('new optimum vector of basic solutions is:');
      xnew
      disp('prob 1  c solved')
      'hit enter to continue'
      pause
    elseif part==4,
      disp('for Part d of Prob 1:')
      deltab=[ 0 -1 0]';     % change in the rhs vector
      aoo=[ao eye(m)];   %  original matrix after adding slacks
      ab=aoo(:,bas);     %  basic columns at optimum
      an=aoo(:,nbas);    %  nonbasic columns at optimum
      abi=inv(ab);       %  inverse basis matrix (can also be obtained 
                            %  from final tableau)
      deltax=abi*deltab;  % change in optimum
      xnew=b+deltax;      %new optimum
      disp('new optimum vector of basic solutions is:')
      xnew
      disp('prob 1 d  solved')
      'hit enter to continue'
    end
  end  % for part=....
%%%%%%%%%%%%%%%%%%%%%%%prob 2%%%%%%%%%%%%%%%%%%%
elseif probno==2,
  co=[-3 -2 3 6 -10 5];  % the original cost from the max problem
  ao=[
  1 2 0 1 0 -6
  0 1 1 3 -2 -1
  1 2 1 3 -1 -5
  ];
  bo=[11 6 13]';
  a=ao;b=bo;c=co;
  m1=0;m2=3;
  disp('The numbers 1,2,3,4,5 correspond to parts a,b,c,d,e of problem 2.')
  disp('Enter the parts of prob 2 to solve as a vector,') 
  parti=input('e.g. 1 or  [1,3] or [1,2,3,4,5]  :  ')
 for i=1:length(parti),
  part=parti(i);
  if part==1,
    disp('for Part a of Prob 2:')
    bas=[1 2 3];
    nbas=[4 5 6];
    ab=ao(:,bas);
    cb=co(bas);
    cn=co(nbas);
    y=cb/ab;
    disp('reduced costs found by pricing out are <= 0:')
    cnbar=cn-y*ao(:,nbas)
    disp('vector of basic variables is >=0:')
    ab\bo
    disp('above shows that the basis is optimal')
    disp('prob 2 a solved')
    disp('hit enter to continue')
    pause
  elseif part==2,
    disp('for Part b of Prob 2:')
    a=[ao [3 -1 1]' ];
    c=[co 4];
    aoo=a;coo=c;
    bas=[1 2 3];
    nbas=[4 5 6 7];    % there are 7 columns in a now
    ab=a(:,bas);
    cb=c(bas);
    cn=c(nbas);
    y=cb/ab;
    disp('pricing out shows that the new column enters the basis')
    cnbar=cn-y*a(:,nbas)
    disp('right hand side is')
    xb=ab\bo
    disp('column for ratio check is')
    ajbar=ab\a(:,7)
    disp('1 leaves and 7 enters the basis')
    bas=[7 2 3];
    nbas=[4 5 6 1];
    ab=a(:,bas);
    cb=c(bas);
    cn=c(nbas);
    y=cb/ab;
    disp('Pricing out shows that we still do not have an optimum.')
    cnbar=cn-y*a(:,nbas)
    disp('Let column 6 enter the basis.')
    disp('the column for the ratios is nonpositive ')
    ab\a(:,6)
    disp('Problem is UNBOUNDED')
    disp('prob 2 b  solved')
    disp('hit enter to continue')
    pause
  elseif part==3,
    disp('for Part c of Prob 2:')
    bas=[1 2 3];
    ab=ao(:,bas);
    xbas=ab\bo;             % basic variable 1,2,3
    disp('After adding a slack variable, we see that the current')
    disp('optimum violates the new constraint.')
    -7 - [1 -1 3]*xbas
    disp('Therefore continue with dual simplex method till feasibility')
    disp('is recovered.')
    disp('We now provide the optimum by using phase 1 and 2 on the primal.')
    a=ao;b=bo;c=-[co 0]; % use -co since we are minimizing
    a=[a [0 0 0]'
     -[1 -1 3 0 0 0 1] ];
    m2=4;
    b=[b
      7];
    simplex
    m2=3;
    disp('prob 2 c  solved')
    disp('hit enter to continue')
    pause
  elseif part==4,
    disp('for Part d of Prob 2:')
    bas=[1 2 3];
    nbas=[4 5 6];
    ab=ao(:,bas);
    c=co;
    c(1) = -5;
    cb=c(bas);
    cn=c(nbas);
    y=cb/ab;
    disp('New reduced costs found by pricing out are not <= 0:')
    cnbar=cn-y*ao(:,nbas)
    a=ao;b=bo;c=-c;   % -co because we are minimizing
    disp('hit return to continue with simplex method')
    pause
    % To avoid using phase 1, multiply through by inverse of basis matrix.
    %  Pricing out is down in reg.m.
    a=ab\a;b=ab\b;
    pt=26
    reg   % skip phase 1 - use previous optimal basis
    disp('prob 2 d  solved')
    disp('hit enter to continue')
    pause
  elseif part==5,
    disp('for Part e of Prob 2:')
    bas=[1 2 3];
    ab=ao(:,bas);
    xbas=ab\bo;
    deltab= [4 0 0]';
    disp('The new optimum is infeasible.')
    xbas+ab\deltab
    disp('Continue with the dual simplex method to recover optimality.')
    disp('We show that the new problem is infeasible using primal simplex.')
    a=ao;b=bo;c=co;b(1)=15;
    simplex
    disp('prob 2 e  solved')
    disp('hit enter to continue')
    pause
  end
 end
elseif probno==3,
  disp('We present the solution by applying the primal')
  disp('simplex method to the dual problem,')
  disp('obtained after changing the constraints to >=')
  co=-[3 2 -4];  % the - is used because our algor. minimizes
  ao= [
  -1 2 1 -1 -1  
  2 1 1 -1 -1
  -1 -1 2 -2 3 
  ]';      % negative and transposed A
  bo=[3 4 2 1 5]';
  a=ao;b=bo;c=co;
  m1=5;m2=0;
  simplex
end
