Mercurial > octave-nkf

comparison scripts/control/util/__zgpbal__.m @ 7132:b01db194c526

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-11-08 16:17:34 by jwe]
author jwe
date Thu, 08 Nov 2007 16:17:34 +0000
parents a1dbe9d80eee
children
comparison
equal deleted inserted replaced
7131:a184bc985c37 7132:b01db194c526
45 ## Created: July 24, 1992 45 ## Created: July 24, 1992
46 ## Conversion to Octave by R. Bruce Tenison July 3, 1994 46 ## Conversion to Octave by R. Bruce Tenison July 3, 1994
47 47
48 function retsys = __zgpbal__ (Asys) 48 function retsys = __zgpbal__ (Asys)
49 49
50 if( (nargin != 1) | (!isstruct(Asys))) 50 if (nargin != 1 || ! isstruct (Asys))
51 print_usage (); 51 print_usage ();
52 endif 52 endif
53 53
54 Asys = sysupdate(Asys,"ss"); 54 Asys = sysupdate (Asys, "ss");
55 [a,b,c,d] = sys2ss(Asys); 55 [a, b, c, d] = sys2ss (Asys);
56 56
57 [nn,mm,pp] = abcddim(a,b,c,d); 57 [nn, mm, pp] = abcddim (a, b, c, d);
58 58
59 np1 = nn+1; 59 np1 = nn+1;
60 nmp = nn+mm+pp; 60 nmp = nn+mm+pp;
61 61
62 ## set up log vector zz, incidence matrix ff 62 ## set up log vector zz, incidence matrix ff
63 zz = zginit(a,b,c,d); 63 zz = zginit (a, b, c, d);
64 64
65 ## disp("__zgpbal__: zginit returns") 65 ## disp("__zgpbal__: zginit returns")
66 ## zz 66 ## zz
67 ## disp("/__zgpbal__") 67 ## disp("/__zgpbal__")
68 68
69 if (norm(zz)) 69 if (norm (zz))
70 ## generalized conjugate gradient approach 70 ## generalized conjugate gradient approach
71 xx = zgscal(a,b,c,d,zz,nn,mm,pp); 71 xx = zgscal (a, b, c, d, zz, nn, mm, pp);
72 72
73 for i=1:nmp 73 for i = 1:nmp
74 xx(i) = floor(xx(i)+0.5); 74 xx(i) = floor (xx(i)+0.5);
75 xx(i) = 2.0^xx(i); 75 xx(i) = 2.0^xx(i);
76 endfor 76 endfor
77 77
78 ## now scale a 78 ## now scale a
79 ## block 1: a = sigma a inv(sigma) 79 ## block 1: a = sigma a inv(sigma)
80 for i=1:nn 80 for i = 1:nn
81 a(i,1:nn) = a(i,1:nn)*xx(i); 81 a(i,1:nn) = a(i,1:nn)*xx(i);
82 a(1:nn,i) = a(1:nn,i)/xx(i); 82 a(1:nn,i) = a(1:nn,i)/xx(i);
83 endfor 83 endfor
84 ## block 2: b= sigma a phi 84 ## block 2: b= sigma a phi
85 for j=1:mm 85 for j = 1:mm
86 j1 = j+nn; 86 j1 = j+nn;
87 b(1:nn,j) = b(1:nn,j)*xx(j1); 87 b(1:nn,j) = b(1:nn,j)*xx(j1);
88 endfor 88 endfor
89 for i=1:nn 89 for i = 1:nn
90 b(i,1:mm) = b(i,1:mm)*xx(i); 90 b(i,1:mm) = b(i,1:mm)*xx(i);
91 endfor 91 endfor
92 for i=1:pp 92 for i = 1:pp
93 i1 = i+nn+mm; 93 i1 = i+nn+mm;
94 ## block 3: c = psi C inv(sigma) 94 ## block 3: c = psi C inv(sigma)
95 c(i,1:nn) = c(i,1:nn)*xx(i1); 95 c(i,1:nn) = c(i,1:nn)*xx(i1);
96 endfor 96 endfor
97 for j=1:nn 97 for j = 1:nn
98 c(1:pp,j) = c(1:pp,j)/xx(j); 98 c(1:pp,j) = c(1:pp,j)/xx(j);
99 endfor 99 endfor
100 ## block 4: d = psi D phi 100 ## block 4: d = psi D phi
101 for j=1:mm 101 for j = 1:mm
102 j1 = j+nn; 102 j1 = j+nn;
103 d(1:pp,j) = d(1:pp,j)*xx(j1); 103 d(1:pp,j) = d(1:pp,j)*xx(j1);
104 endfor 104 endfor
105 for i=1:pp 105 for i = 1:pp
106 i1 = i + nn + mm; 106 i1 = i + nn + mm;
107 d(i,1:mm) = d(i,1:mm)*xx(i1); 107 d(i,1:mm) = d(i,1:mm)*xx(i1);
108 endfor 108 endfor
109 endif 109 endif
110 110
111 retsys = ss(a,b,c,d); 111 retsys = ss (a, b, c, d);
112
112 endfunction 113 endfunction
113