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