Mercurial > octave-nkf

comparison scripts/control/util/zgfslv.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 f084ba47812b
children
comparison
equal deleted inserted replaced
7131:a184bc985c37 7132:b01db194c526
30 if (nargin != 4) 30 if (nargin != 4)
31 print_usage (); 31 print_usage ();
32 endif 32 endif
33 33
34 nmp = n+m+p; 34 nmp = n+m+p;
35 gam1 = (2*n)+m+p; gam2 = n+p; gam3 = n+m; 35 gam1 = (2*n)+m+p;
36 gam2 = n+p;
37 gam3 = n+m;
36 38
37 G1 = givens(sqrt(m),-sqrt(p))'; 39 G1 = givens (sqrt (m), -sqrt (p))';
38 G2 = givens(m+p,sqrt(n*(m+p)))'; 40 G2 = givens (m+p, sqrt (n*(m+p)))';
39 41
40 x = b; 42 x = b;
41 43
42 ## 1) U1 e^n = sqrt(n)e_1^n 44 ## 1) U1 e^n = sqrt(n)e_1^n
43 ## 2) U2 e^m = sqrt(m)e_1^m 45 ## 2) U2 e^m = sqrt(m)e_1^m
44 ## 3) U3 e^p = sqrt(p)e_1^p 46 ## 3) U3 e^p = sqrt(p)e_1^p
45 xdx1 = 1:n; xdx2 = n+(1:m); xdx3 = n+m+(1:p); 47 xdx1 = 1:n;
46 x(xdx1,1) = zgshsr(x(xdx1,1)); 48 xdx2 = n+(1:m);
47 x(xdx2,1) = zgshsr(x(xdx2,1)); 49 xdx3 = n+m+(1:p);
48 x(xdx3,1) = zgshsr(x(xdx3,1)); 50
51 x(xdx1,1) = zgshsr (x(xdx1,1));
52 x(xdx2,1) = zgshsr (x(xdx2,1));
53 x(xdx3,1) = zgshsr (x(xdx3,1));
49 54
50 ## 4) Givens rotations to reduce stray non-zero elements 55 ## 4) Givens rotations to reduce stray non-zero elements
51 idx1 = [n+1,n+m+1]; idx2 = [1,n+1]; 56 idx1 = [n+1, n+m+1];
57 idx2 = [1, n+1];
58
52 x(idx1) = G1'*x(idx1); 59 x(idx1) = G1'*x(idx1);
53 x(idx2) = G2'*x(idx2); 60 x(idx2) = G2'*x(idx2);
54 61
55 ## 6) Scale x, then back-transform to get x 62 ## 6) Scale x, then back-transform to get x
56 en = ones(n,1); em = ones(m,1); ep = ones(p,1); 63 en = ones (n, 1);
57 lam = [gam1*en;gam2*em;gam3*ep]; 64 em = ones (m, 1);
65 ep = ones (p, 1);
66 lam = [gam1*en; gam2*em; gam3*ep];
58 lam(1) = n+m+p; 67 lam(1) = n+m+p;
59 lam(n+1) = 1; # dummy value to avoid divide by zero 68 lam(n+1) = 1; # dummy value to avoid divide by zero
60 lam(n+m+1)=n+m+p; 69 lam(n+m+1) = n+m+p;
61 70
62 x = x ./ lam; x(n+1) = 0; # minimum norm solution 71 x = x ./ lam;
72 x(n+1) = 0; # minimum norm solution
63 73
64 ## back transform now. 74 ## back transform now.
65 x(idx2) = G2*x(idx2); 75 x(idx2) = G2*x(idx2);
66 x(idx1) = G1*x(idx1); 76 x(idx1) = G1*x(idx1);
67 x(xdx3,1) = zgshsr(x(xdx3,1)); 77 x(xdx3,1) = zgshsr (x(xdx3,1));
68 x(xdx2,1) = zgshsr(x(xdx2,1)); 78 x(xdx2,1) = zgshsr (x(xdx2,1));
69 x(xdx1,1) = zgshsr(x(xdx1,1)); 79 x(xdx1,1) = zgshsr (x(xdx1,1));
70 80
71 endfunction 81 endfunction
72 82