home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 2.ddi / MUTOOLS2.DI$ / TWOBL.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  6.1 KB  |  245 lines
  1. %function [cm,dl,dri,sens,rowp,lb,err,tprog,iter] = twobl(mat,blk,ld)
  2.  
  3. %mod 3/10
  4. function [cm,dl,dri,sens,rowp,lb,err,tprog,cprog,iter] = twobl(mat,blk,ld)
  5.  blkp = ptrs(blk);
  6.  [nr,nc]=size(mat);
  7.  err = 0;
  8.  tprog = [];
  9.  cprog = [];
  10.  ml = mat;
  11.  mm = mat;
  12.  mr = mat;
  13.  ur = mat(1:blk(1,2),blk(1,1)+1:nc);
  14.  ll = mat(blk(1,2)+1:nr,1:blk(1,1));
  15.  mt = mat;
  16.  if nargin == 2
  17.    tl = -8;
  18.    tm = 0;
  19.    tr = 8;
  20.    ml(1:blk(1,2),blk(1,1)+1:nc) = exp(tl)*ur;
  21.    ml(blk(1,2)+1:nr,1:blk(1,1)) = exp(-tl)*ll;
  22.    mr(1:blk(1,2),blk(1,1)+1:nc) = exp(tr)*ur;
  23.    mr(blk(1,2)+1:nr,1:blk(1,1)) = exp(-tr)*ll;
  24.    cl = norm(ml);
  25.    cr = norm(mr);
  26.    cm = norm(mm);
  27.    if cl < cm
  28.      mu = cl;
  29.      dl = [exp(tl)*ones(1,blk(1,2)) ones(1,blk(2,2))];
  30.      dri = [exp(-tl)*ones(1,blk(1,1)) ones(1,blk(2,1))];
  31.      dl = diag(dl);
  32.      dri = diag(dri);
  33.      [lb,dummy,rowp] = cheapl(ml,blk,blkp);
  34.      scaledm = dl*mat*dri;
  35.      ur = scaledm(1:blk(1,2),blk(1,1)+1:nc);
  36.      ll = scaledm(blk(1,2)+1:nr,1:blk(1,1));
  37.      val = norm(ur) + norm(ll);
  38.      sens = [val val];
  39.      err = 1;
  40. %     disp('INF IS NOT ACHIEVeD')
  41. %    inf is not achieved
  42.      return
  43.    elseif cr < cm
  44.      mu = cr;
  45.      dl = [exp(tr)*ones(1,blk(1,2)) ones(1,blk(2,2))];
  46.      dri = [exp(-tr)*ones(1,blk(1,1)) ones(1,blk(2,1))];
  47.      dl = diag(dl);
  48.      dri = diag(dri);
  49.      [lb,dummy,rowp] = cheapl(ml,blk,blkp);
  50.      err = 1;
  51. %     disp('INF IS NOT ACHIEVeD')
  52. %    inf is not achieved
  53.      scaledm = dl*mat*dri;
  54.      ur = scaledm(1:blk(1,2),blk(1,1)+1:nc);
  55.      ll = scaledm(blk(1,2)+1:nr,1:blk(1,1));
  56.      val = norm(ur) + norm(ll);
  57.      sens = [val val];
  58.      return
  59.    end
  60.  else
  61.    [tl,tm,tr,cl,cm,cr,nop] = initt(mat,blk,ld);
  62.    if nop == -1
  63.      mu = cr;
  64.      dl = [exp(tr)*ones(1,blk(1,2)) ones(1,blk(2,2))];
  65.      dri = [exp(-tr)*ones(1,blk(1,1)) ones(1,blk(2,1))];
  66.      dl = diag(dl);
  67.      dri = diag(dri);
  68.      [lb,dummy,rowp] = cheapl(ml,blk,blkp);
  69.      err = 1;
  70. %     disp('INF IS NOT ACHIEVD')
  71. %    inf is not achieved
  72.      scaledm = dl*mat*dri;
  73.      ur = scaledm(1:blk(1,2),blk(1,1)+1:nc);
  74.      ll = scaledm(blk(1,2)+1:nr,1:blk(1,1));
  75.      val = norm(ur) + norm(ll);
  76.      sens = [val val];
  77.      return
  78.    end
  79.  end
  80.  itermax = 13;
  81.  rel = 1;
  82.  iter = 0;
  83.  while rel > 0.00000001 & iter<6*itermax
  84.    val = [tr;tm;tl];
  85.    c3 = [ 1;1;1];
  86.    co = [(val .^2) val c3];
  87.    ve = [cr;cm;cl];
  88.    param = co \ ve;
  89.    t = -param(2)/(2*param(1));
  90.    if t < tm
  91.      t = (tl+tm)/2;
  92.    else
  93.      t = (tr+tm)/2;
  94.    end
  95.    mt(1:blk(1,2),blk(1,1)+1:nc) = exp(t)*ur;
  96.    mt(blk(1,2)+1:nr,1:blk(1,1)) = exp(-t)*ll;
  97.    c = norm(mt);
  98.    if t < tm
  99.      if c < cm
  100.        tr = tm;
  101.        tm = t;
  102.        cr = cm;
  103.        cm = c;
  104.      elseif c > cm
  105.        tl = t;
  106.        cl = c;
  107.      else
  108.        tl = t;
  109.        tr = tm;
  110.        cl = c;
  111.        cr = cm;
  112.        tm = (tl+tr)/2;
  113.        mt(1:blk(1,2),blk(1,1)+1:nc) = exp(tm)*ur;
  114.        mt(blk(1,2)+1:nr,1:blk(1,1)) = exp(-tm)*ll;
  115.        cm = norm(mt);
  116.      end
  117.    elseif t > tm
  118.      if c < cm
  119.        tl = tm;
  120.        cl = cm;
  121.        tm = t;
  122.        cm = c;
  123.      elseif c > cm
  124.        tr = t;
  125.        cr = c;
  126.      else
  127.        tl = tm;
  128.        tr = t;
  129.        cl = cm;
  130.        cr = c;
  131.        tm = (tl+tr)/2;
  132.        mt(1:blk(1,2),blk(1,1)+1:nc) = exp(tm)*ur;
  133.        mt(blk(1,2)+1:nr,1:blk(1,1)) = exp(-tm)*ll;
  134.        cm = norm(mt);
  135.      end
  136.    end
  137.    tprog = [tprog ; iter tl tm tr cm];
  138.    cprog = [cprog ; iter cl cm cr cm];
  139.    rel = (max([cr cl]) - cm)/abs(cm);
  140.    iter = iter+1;
  141.  end
  142.  dl = [exp(tm)*ones(1,blk(1,2)) ones(1,blk(2,2))];
  143.  dri = [exp(-tm)*ones(1,blk(1,1)) ones(1,blk(2,1))];
  144.  dl = diag(dl);
  145.  dri = diag(dri);
  146.  [u,s,v] = svd(mt);
  147.  s = diag(s);
  148. % changed coal to .999  from .9995
  149.  coal = length(find(s>(0.999*s(1))));
  150.  e1 = u(1:blk(1,2),1:coal);
  151.  f1 = v(1:blk(1,1),1:coal);
  152.  e2 = u(blk(1,2)+1:nr,1:coal);
  153.  f2 = v(blk(1,1)+1:nc,1:coal);
  154.  if coal == 1
  155.    del1 = abs(e1'*e1 - f1'*f1);
  156.    del2 = abs(e2'*e2 - f2'*f2);
  157.    if max(del1,del2) > 1e-4
  158. %    added 13
  159.      pee1 = f1*e1';
  160.      pee2 = f2*e2';
  161.      if norm(pee1) > 1e-4 & norm(pee2) > 1e-4
  162.        peee = daug((1/(norm(f1)*norm(e1)))*pee1,(1/(norm(f2)*norm(e2)))*pee2);
  163.        lb = max(abs(eig(mt*peee)));
  164.        rowp = [];
  165.        for i=1:blk(1,1)
  166.          rowp = [rowp peee(i,1:blk(1,2))];
  167.        end
  168.        for i=1:blk(2,1)
  169.          rowp = [rowp peee(i+blk(1,1),blk(1,2)+1:nr)];
  170.        end
  171.        rowp = (1/lb) * rowp;
  172. %      till above
  173.      else
  174.        err = -max(del1,del2);
  175.        [lb,dummy,rowp] = cheapl(mt,blk,blkp);
  176.      end
  177. %    keep sens calc outside
  178.      scaledm = dl*mat*dri;
  179.      ur = scaledm(1:blk(1,2),blk(1,1)+1:nc);
  180.      ll = scaledm(blk(1,2)+1:nr,1:blk(1,1));
  181.      val = norm(ur) + norm(ll);
  182.      sens = [val val];
  183.      return
  184.    else
  185.      eta = 1;
  186.    end
  187.  else
  188.    [wvec,wval]= eig(e1'*e1 - f1'*f1);
  189.    wval = diag(real(wval));
  190.    [wval,index] = sort(wval);
  191.    wvec = wvec(:,index);
  192.    if wval(1) > 0 | wval(coal) < 0
  193.      err = -1;
  194.      [lb,dummy,rowp] = cheapl(mt,blk,blkp);
  195.      scaledm = dl*mat*dri;
  196.      ur = scaledm(1:blk(1,2),blk(1,1)+1:nc);
  197.      ll = scaledm(blk(1,2)+1:nr,1:blk(1,1));
  198.      val = norm(ur) + norm(ll);
  199.      sens = [val val];
  200.      return
  201.    else
  202.      eta = zeros(coal,1);
  203.      eta(1) = -wval(coal);
  204.      eta(coal) = wval(1);
  205.      if norm(eta) == 0
  206.        eta(1) = 1; 
  207.      else
  208.        eta = eta/norm(eta);
  209.      end
  210.      eta = wvec*eta;
  211.    end
  212.  end
  213.  e1e = e1*eta;
  214.  e2e = e2*eta;
  215.  f1e = f1*eta;
  216.  f2e = f2*eta;
  217.  rowp = [];
  218.  if norm(e1e) == 0
  219.    rowp = zeros(1,blk(1,1)*blk(1,2));
  220.  else
  221.    f1e = f1e/norm(f1e);
  222.    e1e = e1e/norm(e1e);
  223.    for i=1:blk(1,1)
  224.      rowp = [rowp f1e(i)*e1e'];
  225.    end
  226.  end
  227.  if norm(e2e) == 0
  228.    rowp = [rowp zeros(1,blk(2,1)*blk(2,2))];
  229.  else
  230.    f2e = f2e/norm(f2e);
  231.    e2e = e2e/norm(e2e);
  232.    for i=1:blk(2,1)
  233.      rowp = [rowp f2e(i)*e2e'];
  234.    end
  235.  end
  236.  rowp = rowp/cm;
  237. % add sensitivity calculation
  238.  scaledm = dl*mat*dri;
  239.  ur = scaledm(1:blk(1,2),blk(1,1)+1:nc);
  240.  ll = scaledm(blk(1,2)+1:nr,1:blk(1,1));
  241.  val = norm(ur) + norm(ll);
  242.  sens = [val val];
  243. %
  244. % Copyright MUSYN INC 1991,  All Rights Reserved
  245.