%function [cm,dl,dri,sens,err,tprog,iter] = twoblu(mat,blk,ld)

%mod 3/10
function [cm,dl,dri,sens,err,tprog,cprog,iter] = twoblu(mat,blk,ld)
 blkp = ptrs(blk);
 [nr,nc]=size(mat);
 err = 0;
 tprog = [];
 cprog = [];
 ml = mat;
 mm = mat;
 mr = mat;
 ur = mat(1:blk(1,2),blk(1,1)+1:nc);
 ll = mat(blk(1,2)+1:nr,1:blk(1,1));
 mt = mat;
 if nargin == 2
   tl = -8;
   tm = 0;
   tr = 8;
   ml(1:blk(1,2),blk(1,1)+1:nc) = exp(tl)*ur;
   ml(blk(1,2)+1:nr,1:blk(1,1)) = exp(-tl)*ll;
   mr(1:blk(1,2),blk(1,1)+1:nc) = exp(tr)*ur;
   mr(blk(1,2)+1:nr,1:blk(1,1)) = exp(-tr)*ll;
   cl = norm(ml);
   cr = norm(mr);
   cm = norm(mm);
   if cl < cm
     mu = cl;
     dl = [exp(tl)*ones(1,blk(1,2)) ones(1,blk(2,2))];
     dri = [exp(-tl)*ones(1,blk(1,1)) ones(1,blk(2,1))];
     dl = diag(dl);
     dri = diag(dri);
%     [lb,dummy,rowp] = cheapl(ml,blk,blkp);
     err = 1;
%     disp('INF IS NOT ACHIEVeD')
%    inf is not achieved
     return
   elseif cr < cm
     mu = cr;
     dl = [exp(tr)*ones(1,blk(1,2)) ones(1,blk(2,2))];
     dri = [exp(-tr)*ones(1,blk(1,1)) ones(1,blk(2,1))];
     dl = diag(dl);
     dri = diag(dri);
%     [lb,dummy,rowp] = cheapl(ml,blk,blkp);
     err = 1;
%     disp('INF IS NOT ACHIEVeD')
%    inf is not achieved
     return
   end
 else
   [tl,tm,tr,cl,cm,cr,nop] = initt(mat,blk,ld);
   if nop == -1
     mu = cr;
     dl = [exp(tr)*ones(1,blk(1,2)) ones(1,blk(2,2))];
     dri = [exp(-tr)*ones(1,blk(1,1)) ones(1,blk(2,1))];
     dl = diag(dl);
     dri = diag(dri);
%     [lb,dummy,rowp] = cheapl(ml,blk,blkp);
     err = 1;
%     disp('INF IS NOT ACHIEVD')
%    inf is not achieved
     return
   end
 end
 itermax = 13;
 rel = 1;
 iter = 0;
 while rel > 0.00000001 & iter<6*itermax
   val = [tr;tm;tl];
   c3 = [ 1;1;1];
   co = [(val .^2) val c3];
   ve = [cr;cm;cl];
   param = co \ ve;
   t = -param(2)/(2*param(1));
   if t < tm
     t = (tl+tm)/2;
   else
     t = (tr+tm)/2;
   end
   mt(1:blk(1,2),blk(1,1)+1:nc) = exp(t)*ur;
   mt(blk(1,2)+1:nr,1:blk(1,1)) = exp(-t)*ll;
   c = norm(mt);
   if t < tm
     if c < cm
       tr = tm;
       tm = t;
       cr = cm;
       cm = c;
     elseif c > cm
       tl = t;
       cl = c;
     else
       tl = t;
       tr = tm;
       cl = c;
       cr = cm;
       tm = (tl+tr)/2;
       mt(1:blk(1,2),blk(1,1)+1:nc) = exp(tm)*ur;
       mt(blk(1,2)+1:nr,1:blk(1,1)) = exp(-tm)*ll;
       cm = norm(mt);
     end
   elseif t > tm
     if c < cm
       tl = tm;
       cl = cm;
       tm = t;
       cm = c;
     elseif c > cm
       tr = t;
       cr = c;
     else
       tl = tm;
       tr = t;
       cl = cm;
       cr = c;
       tm = (tl+tr)/2;
       mt(1:blk(1,2),blk(1,1)+1:nc) = exp(tm)*ur;
       mt(blk(1,2)+1:nr,1:blk(1,1)) = exp(-tm)*ll;
       cm = norm(mt);
     end
   end
   tprog = [tprog ; iter tl tm tr cm];
   cprog = [cprog ; iter cl cm cr cm];
   rel = (max([cr cl]) - cm)/abs(cm);
   iter = iter+1;
 end
 dl = [exp(tm)*ones(1,blk(1,2)) ones(1,blk(2,2))];
 dri = [exp(-tm)*ones(1,blk(1,1)) ones(1,blk(2,1))];
 dl = diag(dl);
 dri = diag(dri);
% add sensitivity calculation
 scaledm = dl*mat*dri;
 ur = scaledm(1:blk(1,2),blk(1,1)+1:nc);
 ll = scaledm(blk(1,2)+1:nr,1:blk(1,1));
 val = norm(ur) + norm(ll);
 sens = [val val];

%
% Copyright MUSYN INC 1991,  All Rights Reserved
