home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 2.ddi / MUTOOLS2.DI$ / GENLFT.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  1.3 KB  |  39 lines
  1. % function [out] = genlft(top,bot,topfeedout,botfeedout)
  2. %   forms Redheffer star-product of two CONSTANT matrices.
  3. %   main subroutine for STARP.
  4.  
  5. function [out] = genlft(top,bot,dim1,dim2)
  6.    if nargin ~= 4
  7.      disp('usage: [out] = genlft(top,bot,topfeedout,botfeedout)')
  8.      return
  9.    end
  10.    [nrtop,nctop] = size(top);
  11.    [nrbot,ncbot] = size(bot);
  12.    dimout1 = nrtop - dim1;
  13.    dimin1 = nctop - dim2;
  14.    dimout2 = nrbot - dim2;
  15.    dimin2 = ncbot - dim1;
  16.    top11 = top(1:dimout1,1:dimin1);
  17.    top12 = top(1:dimout1,dimin1+1:nctop);
  18.    top21 = top(dimout1+1:nrtop,1:dimin1);
  19.    top22 = top(dimout1+1:nrtop,dimin1+1:nctop);
  20.    bot11 = bot(1:dim2,1:dim1);
  21.    bot12 = bot(1:dim2,dim1+1:ncbot);
  22.    bot21 = bot(dim2+1:nrbot,1:dim1);
  23.    bot22 = bot(dim2+1:nrbot,dim1+1:ncbot);
  24.    tb = top22 * bot11;
  25.    [nrtb nctb] = size(tb);
  26.    imtb = eye(nrtb) - tb;
  27. %  should check the invertibility of imtb at this point, but
  28. %  we have had trouble with COND on large, sparse matrices
  29.    bt = bot11*top22;
  30.    [nrbt,ncbt] = size(bt);
  31.    imbt = eye(nrbt) - bt;
  32.    x = imbt\bot12;
  33.    y = imtb\top21;
  34.    upper = [ (top11 + top12 * bot11 * y) top12*x ];
  35.    lower = [ bot21*y (bot22 + bot21 * top22 * x) ];
  36.    out = [upper ; lower];
  37. %
  38. % Copyright MUSYN INC 1991,  All Rights Reserved
  39.