home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 2.ddi / MUTOOLS2.DI$ / FMLINP.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.2 KB  |  64 lines
  1. %function [basic,sol,cost,lambda,tnpiv,flopcnt] = fmlinp(a,b,c,startbasic)
  2. %
  3. %  *****  UNTESTED SUBROUTINE  *****
  4. %
  5. % Simplex routine for linear programs arising in FITMAGLP.
  6. %  A modification of LINP that exploits some structure in FITMAGLP
  7. %  and tailors some constants.  Also tries to avoid degeneracy in 
  8. %  the special structure of FITMAGLP.
  9.  
  10. function [basic,sol,cost,lambda,tnpiv,flopcnt] = fmlinp(a,b,c,startbasic)
  11.  if nargin < 4
  12.    disp('usage: [basic,sol,cost,lambda] = lpmagfit(a,b,c,startbasic)')
  13.    return
  14.  end
  15.  flopin = flops;tnpiv=0;
  16.  [con,var] = size(a);
  17.  [nrb,ncb] = size(b);
  18.  [nrc,ncc] = size(c);
  19.  if con ~= nrb | var ~= ncc
  20.    disp('error dimensioning')
  21.    return
  22.  end
  23. %check conditioning of startbasic and change if necessary
  24. startbasic=sort(startbasic); stdiff=diff(startbasic);
  25. if  any([length(startbasic)~=con,stdiff==0,max(startbasic)>var]),
  26.    [q2,r2,e2]=qr(a); for i=1:con, startbasic(i)=find(e2(:,i));end;
  27.  else [q1,r1,e1]=qr(a(:,startbasic));
  28.       [q2,r2,e2]=qr(a);
  29.      if r1(con,con)<(1e-4)*r2(con,con), 
  30.        for i=1:con, startbasic(i)=find(e2(:,i));end;end
  31.   end
  32.  sol1=a(:,startbasic)\b;
  33.  negind=find(sol1<0);
  34.  p=length(negind);
  35.  if p>0,posind=startbasic(comple(negind,con));
  36.      else posind=startbasic(con);  end;%if p>0
  37.  feasbasic=startbasic;
  38.  cost1=-sum(sol1(negind));
  39.  if p>0,   %find an initial basic feasible solution
  40.    newstartbasic=[posind var+1:var+p];
  41.    auga = [a -a(:,startbasic(negind))];
  42.    augc = [2*cost1 zeros(1,var-1) ones(1,p)];
  43.    [feasbasic,augsol,cost2,auglambda,piv] = fmlp(auga,b,augc,newstartbasic);
  44.    tnpiv=tnpiv+piv;
  45.    if max(feasbasic)>var, 
  46.      augc(1,1)=0;
  47.      [feasbasic,augsol,cost3,auglambda,piv] = fmlp(auga,b,augc,feasbasic);
  48.      tnpiv=tnpiv+piv;
  49.      [feasbasic,err] = fmfixbas(a,feasbasic,augsol,var,con);
  50.      if err==1,  %problem with feasible solution so form degenerate one.
  51.           feasbasic=sort(feasbasic);augsol=[1 ; zeros(con-1,1)];
  52.           [feasbasic,err] = fmfixbas(a,feasbasic,augsol,var,con);
  53.       end %if err==1
  54.     end %if max(feasbasic)>var
  55.   end %if p>0
  56.  
  57. [basic,sol,cost,lambda,tnpiv] = fmlp(a,b,c,feasbasic);
  58. tnpiv=tnpiv+piv;
  59.  
  60.  flopout = flops;
  61.  flopcnt = flopout - flopin;
  62. %
  63. % Copyright MUSYN INC 1991,  All Rights Reserved
  64.