Mercurial > octave-nkf
comparison scripts/control/util/zgscal.m @ 7132:b01db194c526
[project @ 2007-11-08 16:17:34 by jwe]
author | jwe |
---|---|
date | Thu, 08 Nov 2007 16:17:34 +0000 |
parents | f084ba47812b |
children |
comparison
equal
deleted
inserted
replaced
7131:a184bc985c37 | 7132:b01db194c526 |
---|---|
42 ## Givens rotations, diagonalized 2x2 block of F, gcg vector initialization | 42 ## Givens rotations, diagonalized 2x2 block of F, gcg vector initialization |
43 | 43 |
44 nmp = n+m+p; | 44 nmp = n+m+p; |
45 | 45 |
46 ## x_0 = x_{-1} = 0, r_0 = z | 46 ## x_0 = x_{-1} = 0, r_0 = z |
47 x = zeros(nmp,1); | 47 x = zeros (nmp, 1); |
48 xk1 = x; | 48 xk1 = x; |
49 xk2 = x; | 49 xk2 = x; |
50 rk1 = z; | 50 rk1 = z; |
51 k = 0; | 51 k = 0; |
52 | 52 |
53 ## construct balancing least squares problem | 53 ## construct balancing least squares problem |
54 F = eye(nmp); | 54 F = eye (nmp); |
55 for kk=1:nmp | 55 for kk = 1:nmp |
56 F(1:nmp,kk) = zgfmul(a,b,c,d,F(:,kk)); | 56 F(1:nmp,kk) = zgfmul (a, b, c, d, F(:,kk)); |
57 endfor | 57 endfor |
58 | 58 |
59 [U,H,k1] = krylov(F,z,nmp,1e-12,1); | 59 [U, H, k1] = krylov (F, z, nmp, 1e-12, 1); |
60 if(!issquare(H)) | 60 if (! issquare (H)) |
61 if(columns(H) != k1) | 61 if (columns (H) != k1) |
62 error("zgscal(tzero): k1=%d, columns(H)=%d",k1,columns(H)); | 62 error ("zgscal(tzero): k1=%d, columns(H)=%d", k1, columns (H)); |
63 elseif(rows(H) != k1+1) | 63 elseif (rows (H) != k1+1) |
64 error("zgscal: k1=%d, rows(H) = %d",k1,rows(H)); | 64 error ("zgscal: k1=%d, rows(H) = %d", k1, rows (H)); |
65 elseif ( norm(H(k1+1,:)) > 1e-12*norm(H,"inf") ) | 65 elseif (norm (H(k1+1,:)) > 1e-12*norm (H, "inf")) |
66 zgscal_last_row_of_H = H(k1+1,:) | 66 zgscal_last_row_of_H = H(k1+1,:) |
67 error("zgscal: last row of H nonzero (norm(H)=%e)",norm(H,"inf")) | 67 error ("zgscal: last row of H nonzero (norm(H)=%e)", norm (H, "inf")) |
68 endif | 68 endif |
69 H = H(1:k1,1:k1); | 69 H = H(1:k1,1:k1); |
70 U = U(:,1:k1); | 70 U = U(:,1:k1); |
71 endif | 71 endif |
72 | 72 |
73 ## tridiagonal H can still be rank deficient, so do permuted qr | 73 ## tridiagonal H can still be rank deficient, so do permuted qr |
74 ## factorization | 74 ## factorization |
75 [qq,rr,pp] = qr(H); # H = qq*rr*pp' | 75 [qq, rr, pp] = qr (H); # H = qq*rr*pp' |
76 nn = rank(rr); | 76 nn = rank (rr); |
77 qq = qq(:,1:nn); | 77 qq = qq(:,1:nn); |
78 rr = rr(1:nn,:); # rr may not be square, but "\" does least | 78 rr = rr(1:nn,:); # rr may not be square, but "\" does least |
79 xx = U*pp*(rr\qq'*(U'*z)); # squares solution, so this works | 79 xx = U*pp*(rr\qq'*(U'*z)); # squares solution, so this works |
80 ## xx1 = pinv(F)*z; | 80 ## xx1 = pinv(F)*z; |
81 ## zgscal_x_xx1_err = [xx,xx1,xx-xx1] | 81 ## zgscal_x_xx1_err = [xx,xx1,xx-xx1] |
84 ## the rest of this is left from the original zgscal; | 84 ## the rest of this is left from the original zgscal; |
85 ## I've had some numerical problems with the GCG algorithm, | 85 ## I've had some numerical problems with the GCG algorithm, |
86 ## so for now I'm solving it with the krylov routine. | 86 ## so for now I'm solving it with the krylov routine. |
87 | 87 |
88 ## initialize residual error norm | 88 ## initialize residual error norm |
89 rnorm = norm(rk1,1); | 89 rnorm = norm (rk1, 1); |
90 | 90 |
91 xnorm = 0; | 91 xnorm = 0; |
92 fnorm = 1e-12 * norm([a,b;c,d],1); | 92 fnorm = 1e-12 * norm ([a, b; c, d], 1); |
93 | 93 |
94 ## dummy defines for MATHTOOLS compiler | 94 gamk2 = 0; |
95 gamk2 = 0; omega1 = 0; ztmz2 = 0; | 95 omega1 = 0; |
96 ztmz2 = 0; | |
96 | 97 |
97 ## do until small changes to x | 98 ## do until small changes to x |
98 len_x = length(x); | 99 len_x = length(x); |
99 while ((k < 2*len_x) & (xnorm> 0.5) & (rnorm>fnorm))|(k == 0) | 100 while ((k < 2*len_x && xnorm > 0.5 && rnorm > fnorm) || k == 0) |
100 k = k+1; | 101 k++; |
101 | 102 |
102 ## solve F_d z_{k-1} = r_{k-1} | 103 ## solve F_d z_{k-1} = r_{k-1} |
103 zk1= zgfslv(n,m,p,rk1); | 104 zk1= zgfslv (n, m, p, rk1); |
104 | 105 |
105 ## Generalized CG iteration | 106 ## Generalized CG iteration |
106 ## gamk1 = (zk1'*F_d*zk1)/(zk1'*F*zk1); | 107 ## gamk1 = (zk1'*F_d*zk1)/(zk1'*F*zk1); |
107 ztMz1 = zk1'*rk1; | 108 ztMz1 = zk1'*rk1; |
108 gamk1 = ztMz1/(zk1'*zgfmul(a,b,c,d,zk1)); | 109 gamk1 = ztMz1/(zk1'*zgfmul (a, b, c, d, zk1)); |
109 | 110 |
110 if(rem(k,len_x) == 1) omega = 1; | 111 if (rem (k, len_x) == 1) |
111 else omega = 1/(1-gamk1*ztMz1/(gamk2*omega1*ztmz2)); | 112 omega = 1; |
113 else | |
114 omega = 1/(1-gamk1*ztMz1/(gamk2*omega1*ztmz2)); | |
112 endif | 115 endif |
113 | 116 |
114 ## store x in xk2 to save space | 117 ## store x in xk2 to save space |
115 xk2 = xk2 + omega*(gamk1*zk1 + xk1 - xk2); | 118 xk2 = xk2 + omega*(gamk1*zk1 + xk1 - xk2); |
116 | 119 |
117 ## compute new residual error: rk = z - F xk, check end conditions | 120 ## compute new residual error: rk = z - F xk, check end conditions |
118 rk1 = z - zgfmul(a,b,c,d,xk2); | 121 rk1 = z - zgfmul (a, b, c, d, xk2); |
119 rnorm = norm(rk1); | 122 rnorm = norm (rk1); |
120 xnorm = max(abs(xk1 - xk2)); | 123 xnorm = max (abs (xk1 - xk2)); |
121 | 124 |
122 ## printf("zgscal: k=%d, gamk1=%e, gamk2=%e, \nztMz1=%e ztmz2=%e\n", ... | 125 ## printf("zgscal: k=%d, gamk1=%e, gamk2=%e, \nztMz1=%e ztmz2=%e\n", ... |
123 ## k,gamk1, gamk2, ztMz1, ztmz2); | 126 ## k,gamk1, gamk2, ztMz1, ztmz2); |
124 ## xk2_1_zk1 = [xk2 xk1 zk1] | 127 ## xk2_1_zk1 = [xk2 xk1 zk1] |
125 ## ABCD = [a,b;c,d] | 128 ## ABCD = [a,b;c,d] |
127 | 130 |
128 ## get ready for next iteration | 131 ## get ready for next iteration |
129 gamk2 = gamk1; | 132 gamk2 = gamk1; |
130 omega1 = omega; | 133 omega1 = omega; |
131 ztmz2 = ztMz1; | 134 ztmz2 = ztMz1; |
132 [xk1,xk2] = swap(xk1,xk2); | 135 [xk1, xk2] = swap (xk1, xk2); |
133 endwhile | 136 endwhile |
134 x = xk2; | 137 x = xk2; |
135 | 138 |
136 ## check convergence | 139 ## check convergence |
137 if (xnorm> 0.5 & rnorm>fnorm) | 140 if (xnorm> 0.5 && rnorm > fnorm) |
138 warning("zgscal(tzero): GCG iteration failed; solving with pinv"); | 141 warning ("zgscal(tzero): GCG iteration failed; solving with pinv"); |
139 | 142 |
140 ## perform brute force least squares; construct F | 143 ## perform brute force least squares; construct F |
141 Am = eye(nmp); | 144 Am = eye (nmp); |
142 for ii=1:nmp | 145 for ii = 1:nmp |
143 Am(:,ii) = zgfmul(a,b,c,d,Am(:,ii)); | 146 Am(:,ii) = zgfmul (a, b, c, d, Am(:,ii)); |
144 endfor | 147 endfor |
145 | 148 |
146 ## now solve with qr factorization | 149 ## now solve with qr factorization |
147 x = pinv(Am)*z; | 150 x = pinv (Am) * z; |
148 endif | 151 endif |
152 | |
149 endfunction | 153 endfunction |