home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21fb.zip / octave / SCRIPTS.ZIP / scripts.fat / control / zgpbal.m < prev    next >
Text File  |  1999-12-24  |  3KB  |  110 lines
  1. ## Copyright (C) 1996 Auburn University.  All Rights Reserved
  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, 59 Temple Place, Suite 330, Boston, MA 02111 USA. 
  18.  
  19. ## -*- texinfo -*-
  20. ## @deftypefn {Function File } {retsys =} zgpbal(Asys)
  21. ##
  22. ## used internally in @code{tzero}; minimal argument checking performed
  23. ##
  24. ## implementation of zero computation generalized eigenvalue problem 
  25. ## balancing method (Hodel and Tiller, Allerton Conference, 1991)
  26. ## Based on Ward's balancing algorithm (SIAM J. Sci Stat. Comput., 1981)
  27. ##
  28. ## zgpbal computes a state/input/output weighting that attempts to 
  29. ## reduced the range of the magnitudes of the nonzero elements of [a,b,c,d]
  30. ## The weighting uses scalar multiplication by powers of 2, so no roundoff
  31. ## will occur.  
  32. ##
  33. ## zgpbal should be followed by zgpred
  34. ##
  35. ## @end deftypefn
  36.  
  37. ## References:
  38. ## ZGEP: Hodel, "Computation of Zeros with Balancing," 1992, submitted to  LAA
  39. ## Generalized CG: Golub and Van Loan, "Matrix Computations, 2nd ed" 1989
  40.  
  41. function retsys = zgpbal (Asys)  
  42.  
  43.   ## A. S. Hodel July 24 1992
  44.   ## Conversion to Octave by R. Bruce Tenison July 3, 1994
  45.  
  46.   if( (nargin != 1) | (!is_struct(Asys)))
  47.     usage("retsys = zgpbal(Asys)");
  48.   endif
  49.  
  50.   Asys = sysupdat(Asys,"ss");
  51.   [a,b,c,d] = sys2ss(Asys);
  52.  
  53.   [nn,mm,pp] = abcddim(a,b,c,d);
  54.   
  55.   np1 = nn+1;
  56.   nmp = nn+mm+pp;
  57.  
  58.   ## set up log vector zz, incidence matrix ff
  59.   zz = zginit(a,b,c,d);
  60.  
  61.   ## disp("zgpbal: zginit returns")
  62.   ## zz
  63.   ## disp("/zgpbal")
  64.  
  65.   if (norm(zz))
  66.     ## generalized conjugate gradient approach
  67.     xx = zgscal(a,b,c,d,zz,nn,mm,pp);
  68.     
  69.     for i=1:nmp
  70.       xx(i) = floor(xx(i)+0.5);
  71.       xx(i) = 2.0^xx(i);
  72.     endfor
  73.     
  74.     ## now scale a
  75.     ## block 1: a = sigma a inv(sigma)
  76.     for i=1:nn
  77.       a(i,1:nn) = a(i,1:nn)*xx(i);
  78.       a(1:nn,i) = a(1:nn,i)/xx(i);
  79.     endfor
  80.     ## block 2: b= sigma a phi
  81.     for j=1:mm
  82.       j1 = j+nn;
  83.       b(1:nn,j) = b(1:nn,j)*xx(j1);
  84.     endfor
  85.     for i=1:nn
  86.       b(i,1:mm) = b(i,1:mm)*xx(i);
  87.     endfor
  88.     for i=1:pp
  89.       i1 = i+nn+mm;
  90.       ## block 3: c = psi C inv(sigma)
  91.       c(i,1:nn) = c(i,1:nn)*xx(i1);
  92.     endfor
  93.     for j=1:nn
  94.       c(1:pp,j) = c(1:pp,j)/xx(j);
  95.     endfor
  96.     ## block 4: d = psi D phi
  97.     for j=1:mm
  98.       j1 = j+nn;
  99.       d(1:pp,j) = d(1:pp,j)*xx(j1);
  100.     endfor
  101.     for i=1:pp
  102.       i1 = i + nn + mm;
  103.       d(i,1:mm) = d(i,1:mm)*xx(i1);
  104.     endfor
  105.   endif
  106.   
  107.   retsys = ss2sys(a,b,c,d);
  108. endfunction
  109.  
  110.