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

%mod 3/10
function [cm,dl,dri,sens,rowp,lb,err,tprog,cprog,iter] = twobl(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);
     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];
     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
     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];
     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
     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];
     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);
 [u,s,v] = svd(mt);
 s = diag(s);
% changed coal to .999  from .9995
 coal = length(find(s>(0.999*s(1))));
 e1 = u(1:blk(1,2),1:coal);
 f1 = v(1:blk(1,1),1:coal);
 e2 = u(blk(1,2)+1:nr,1:coal);
 f2 = v(blk(1,1)+1:nc,1:coal);
 if coal == 1
   del1 = abs(e1'*e1 - f1'*f1);
   del2 = abs(e2'*e2 - f2'*f2);
   if max(del1,del2) > 1e-4
%    added 13
     pee1 = f1*e1';
     pee2 = f2*e2';
     if norm(pee1) > 1e-4 & norm(pee2) > 1e-4
       peee = daug((1/(norm(f1)*norm(e1)))*pee1,(1/(norm(f2)*norm(e2)))*pee2);
       lb = max(abs(eig(mt*peee)));
       rowp = [];
       for i=1:blk(1,1)
         rowp = [rowp peee(i,1:blk(1,2))];
       end
       for i=1:blk(2,1)
         rowp = [rowp peee(i+blk(1,1),blk(1,2)+1:nr)];
       end
       rowp = (1/lb) * rowp;
%      till above
     else
       err = -max(del1,del2);
       [lb,dummy,rowp] = cheapl(mt,blk,blkp);
     end
%    keep sens calc outside
     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];
     return
   else
     eta = 1;
   end
 else
   [wvec,wval]= eig(e1'*e1 - f1'*f1);
   wval = diag(real(wval));
   [wval,index] = sort(wval);
   wvec = wvec(:,index);
   if wval(1) > 0 | wval(coal) < 0
     err = -1;
     [lb,dummy,rowp] = cheapl(mt,blk,blkp);
     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];
     return
   else
     eta = zeros(coal,1);
     eta(1) = -wval(coal);
     eta(coal) = wval(1);
     if norm(eta) == 0
       eta(1) = 1; 
     else
       eta = eta/norm(eta);
     end
     eta = wvec*eta;
   end
 end
 e1e = e1*eta;
 e2e = e2*eta;
 f1e = f1*eta;
 f2e = f2*eta;
 rowp = [];
 if norm(e1e) == 0
   rowp = zeros(1,blk(1,1)*blk(1,2));
 else
   f1e = f1e/norm(f1e);
   e1e = e1e/norm(e1e);
   for i=1:blk(1,1)
     rowp = [rowp f1e(i)*e1e'];
   end
 end
 if norm(e2e) == 0
   rowp = [rowp zeros(1,blk(2,1)*blk(2,2))];
 else
   f2e = f2e/norm(f2e);
   e2e = e2e/norm(e2e);
   for i=1:blk(2,1)
     rowp = [rowp f2e(i)*e2e'];
   end
 end
 rowp = rowp/cm;
% 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
