home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21eb.zip / octave / SCRIPTS.ZIP / scripts / control / zginit.m < prev    next >
Text File  |  1999-03-05  |  2KB  |  88 lines
  1. # Copyright (C) 1996,1998 A. Scottedward Hodel 
  2. #
  3. # This file is part of Octave. 
  4. #
  5. # Octave is free software; you can redistribute it and/or modify it 
  6. # under the terms of the GNU General Public License as published by the 
  7. # Free Software Foundation; either version 2, or (at your option) any 
  8. # later version. 
  9. # 
  10. # Octave is distributed in the hope that it will be useful, but WITHOUT 
  11. # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 
  12. # FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License 
  13. # for more details.
  14. # 
  15. # You should have received a copy of the GNU General Public License 
  16. # along with Octave; see the file COPYING.  If not, write to the Free 
  17. # Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. 
  18.  
  19. function zz = zginit(a,b,c,d)
  20.   # zz = zginit(a,b,c,d)
  21.   # construct right hand side vector zz
  22.   # for the zero-computation generalized eigenvalue problem
  23.   # balancing procedure
  24.   # called by zgepbal
  25.   # References:
  26.   # ZGEP: Hodel, "Computation of Zeros with Balancing," Linear Algebra and
  27.   #              its Applications, 1993
  28.   # Generalized CG: Golub and Van Loan, "Matrix Computations, 2nd ed" 1989
  29.   
  30.   # A. S. Hodel July 24 1992
  31.   # Conversion to Octave by R. Bruce Tenison, July 3, 1994
  32.  
  33.   [nn,mm] = size(b);
  34.   [pp,mm] = size(d);
  35.  
  36.   nmp = nn+mm+pp;
  37.  
  38.   # set up log vector zz
  39.   zz = zeros(nmp,1);
  40.  
  41.   # zz part 1:
  42.   for i=1:nn
  43.     # nonzero off diagonal entries of a
  44.     if(nn > 1)
  45.       nidx = complement(i,1:nn);
  46.       a_row_i = a(i,nidx);                 a_col_i = a(nidx,i);
  47.       arnz = a_row_i(find(a_row_i != 0));  acnz = a_col_i(find(a_col_i != 0));
  48.     else
  49.       arnz = acnz = [];
  50.     endif
  51.  
  52.     # row of b
  53.     bidx = find(b(i,:) != 0);
  54.     b_row_i = b(i,bidx);
  55.  
  56.     # column of c
  57.     cidx = find(c(:,i) != 0);
  58.     c_col_i = c(cidx,i);
  59.    
  60.     # sum the entries
  61.     zz(i) = sum(log(abs(acnz))) - sum(log(abs(arnz))) ...
  62.             - sum(log(abs(b_row_i))) + sum(log(abs(c_col_i)));
  63.   endfor
  64.  
  65.   # zz part 2:
  66.   bd = [b;d];
  67.   for i=1:mm
  68.     i1 = i+nn;
  69.  
  70.     # column of [b;d]
  71.     bdidx = find(bd(:,i) != 0);
  72.     bd_col_i = bd(bdidx,i);
  73.     zz(i1) = sum(log(abs(bd_col_i)));
  74.   endfor
  75.  
  76.   # zz part 3:
  77.   cd = [c, d];
  78.   for i=1:pp
  79.     i1 = i+nn+mm;
  80.     cdidx = find(cd(i,:) != 0);
  81.     cd_row_i = cd(i,cdidx);
  82.     zz(i1) = -sum(log(abs(cd_row_i)));
  83.   endfor
  84.  
  85.   # now set zz as log base 2
  86.   zz = zz*(1/log(2));
  87. endfunction
  88.