home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 2.ddi / MUTOOLS2.DI$ / TWOBLU.M < prev   
Encoding:
Text File  |  1993-03-11  |  3.3 KB  |  140 lines
  1. %function [cm,dl,dri,sens,err,tprog,iter] = twoblu(mat,blk,ld)
  2.  
  3. %mod 3/10
  4. function [cm,dl,dri,sens,err,tprog,cprog,iter] = twoblu(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.      err = 1;
  35. %     disp('INF IS NOT ACHIEVeD')
  36. %    inf is not achieved
  37.      return
  38.    elseif cr < cm
  39.      mu = cr;
  40.      dl = [exp(tr)*ones(1,blk(1,2)) ones(1,blk(2,2))];
  41.      dri = [exp(-tr)*ones(1,blk(1,1)) ones(1,blk(2,1))];
  42.      dl = diag(dl);
  43.      dri = diag(dri);
  44. %     [lb,dummy,rowp] = cheapl(ml,blk,blkp);
  45.      err = 1;
  46. %     disp('INF IS NOT ACHIEVeD')
  47. %    inf is not achieved
  48.      return
  49.    end
  50.  else
  51.    [tl,tm,tr,cl,cm,cr,nop] = initt(mat,blk,ld);
  52.    if nop == -1
  53.      mu = cr;
  54.      dl = [exp(tr)*ones(1,blk(1,2)) ones(1,blk(2,2))];
  55.      dri = [exp(-tr)*ones(1,blk(1,1)) ones(1,blk(2,1))];
  56.      dl = diag(dl);
  57.      dri = diag(dri);
  58. %     [lb,dummy,rowp] = cheapl(ml,blk,blkp);
  59.      err = 1;
  60. %     disp('INF IS NOT ACHIEVD')
  61. %    inf is not achieved
  62.      return
  63.    end
  64.  end
  65.  itermax = 13;
  66.  rel = 1;
  67.  iter = 0;
  68.  while rel > 0.00000001 & iter<6*itermax
  69.    val = [tr;tm;tl];
  70.    c3 = [ 1;1;1];
  71.    co = [(val .^2) val c3];
  72.    ve = [cr;cm;cl];
  73.    param = co \ ve;
  74.    t = -param(2)/(2*param(1));
  75.    if t < tm
  76.      t = (tl+tm)/2;
  77.    else
  78.      t = (tr+tm)/2;
  79.    end
  80.    mt(1:blk(1,2),blk(1,1)+1:nc) = exp(t)*ur;
  81.    mt(blk(1,2)+1:nr,1:blk(1,1)) = exp(-t)*ll;
  82.    c = norm(mt);
  83.    if t < tm
  84.      if c < cm
  85.        tr = tm;
  86.        tm = t;
  87.        cr = cm;
  88.        cm = c;
  89.      elseif c > cm
  90.        tl = t;
  91.        cl = c;
  92.      else
  93.        tl = t;
  94.        tr = tm;
  95.        cl = c;
  96.        cr = cm;
  97.        tm = (tl+tr)/2;
  98.        mt(1:blk(1,2),blk(1,1)+1:nc) = exp(tm)*ur;
  99.        mt(blk(1,2)+1:nr,1:blk(1,1)) = exp(-tm)*ll;
  100.        cm = norm(mt);
  101.      end
  102.    elseif t > tm
  103.      if c < cm
  104.        tl = tm;
  105.        cl = cm;
  106.        tm = t;
  107.        cm = c;
  108.      elseif c > cm
  109.        tr = t;
  110.        cr = c;
  111.      else
  112.        tl = tm;
  113.        tr = t;
  114.        cl = cm;
  115.        cr = c;
  116.        tm = (tl+tr)/2;
  117.        mt(1:blk(1,2),blk(1,1)+1:nc) = exp(tm)*ur;
  118.        mt(blk(1,2)+1:nr,1:blk(1,1)) = exp(-tm)*ll;
  119.        cm = norm(mt);
  120.      end
  121.    end
  122.    tprog = [tprog ; iter tl tm tr cm];
  123.    cprog = [cprog ; iter cl cm cr cm];
  124.    rel = (max([cr cl]) - cm)/abs(cm);
  125.    iter = iter+1;
  126.  end
  127.  dl = [exp(tm)*ones(1,blk(1,2)) ones(1,blk(2,2))];
  128.  dri = [exp(-tm)*ones(1,blk(1,1)) ones(1,blk(2,1))];
  129.  dl = diag(dl);
  130.  dri = diag(dri);
  131. % add sensitivity calculation
  132.  scaledm = dl*mat*dri;
  133.  ur = scaledm(1:blk(1,2),blk(1,1)+1:nc);
  134.  ll = scaledm(blk(1,2)+1:nr,1:blk(1,1));
  135.  val = norm(ur) + norm(ll);
  136.  sens = [val val];
  137.  
  138. %
  139. % Copyright MUSYN INC 1991,  All Rights Reserved
  140.