Mercurial > octave-nkf

comparison scripts/control/util/zgfmul.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
35 35
36 if (nargin != 5) 36 if (nargin != 5)
37 print_usage (); 37 print_usage ();
38 endif 38 endif
39 39
40 [n,m] = size(b); 40 [n, m] = size (b);
41 [p,m1] = size(c); 41 [p, m1] = size (c);
42 nm = n+m; 42 nm = n+m;
43 y = zeros(nm+p,1); 43 y = zeros (nm+p, 1);
44 44
45 ## construct F column by column 45 ## construct F column by column
46 for jj=1:n 46 for jj = 1:n
47 Fj = zeros(nm+p,1); 47 Fj = zeros (nm+p, 1);
48 48
49 ## rows 1:n: F1 49 ## rows 1:n: F1
50 aridx = complement(jj,find(a(jj,:) != 0)); 50 aridx = complement (jj, find (a(jj,:) != 0));
51 acidx = complement(jj,find(a(:,jj) != 0)); 51 acidx = complement (jj, find (a(:,jj) != 0));
52 bidx = find(b(jj,:) != 0); 52 bidx = find (b(jj,:) != 0);
53 cidx = find(c(:,jj) != 0); 53 cidx = find (c(:,jj) != 0);
54 54
55 Fj(aridx) = Fj(aridx) - 1; # off diagonal entries of F1 55 Fj(aridx) = Fj(aridx) - 1; # off diagonal entries of F1
56 Fj(acidx) = Fj(acidx) - 1; 56 Fj(acidx) = Fj(acidx) - 1;
57 ## diagonal entry of F1 57 ## diagonal entry of F1
58 Fj(jj) = length(aridx)+length(acidx) + length(bidx) + length(cidx); 58 Fj(jj) = length (aridx) + length (acidx) + length (bidx) + length (cidx);
59 59
60 if(!isempty(bidx)) Fj(n+bidx) = 1; endif # B' incidence 60 ## B' incidence
61 if(!isempty(cidx)) Fj(n+m+cidx) = -1; endif # -C incidence 61 if (! isempty (bidx))
62 Fj(n+bidx) = 1;
63 endif
64
65 ## -C incidence
66 if (! isempty (cidx))
67 Fj(n+m+cidx) = -1;
68 endif
62 y = y + x(jj)*Fj; # multiply by corresponding entry of x 69 y = y + x(jj)*Fj; # multiply by corresponding entry of x
63 endfor 70 endfor
64 71
65 for jj=1:m 72 for jj = 1:m
66 Fj = zeros(nm+p,1); 73 Fj = zeros (nm+p, 1);
67 bidx = find(b(:,jj) != 0); 74 bidx = find (b(:,jj) != 0);
68 if(!isempty(bidx)) Fj(bidx) = 1; endif # B incidence 75 ## B incidence
69 didx = find(d(:,jj) != 0); 76 if (! isempty (bidx))
70 if(!isempty(didx)) Fj(n+m+didx) = 1; endif # D incidence 77 Fj(bidx) = 1;
78 endif
79 didx = find (d(:,jj) != 0);
80 ## D incidence
81 if (! isempty (didx))
82 Fj(n+m+didx) = 1;
83 endif
71 Fj(n+jj) = length(bidx) + length(didx); # F2 is diagonal 84 Fj(n+jj) = length(bidx) + length(didx); # F2 is diagonal
72 y = y + x(n+jj)*Fj; # multiply by corresponding entry of x 85 y = y + x(n+jj)*Fj; # multiply by corresponding entry of x
73 endfor 86 endfor
74 87
75 for jj=1:p 88 for jj = 1:p
76 Fj = zeros(nm+p,1); 89 Fj = zeros (nm+p, 1);
77 cidx = find(c(jj,:) != 0); 90 cidx = find (c(jj,:) != 0);
78 if(!isempty(cidx)) Fj(cidx) = -1; endif # -C' incidence 91 ## -C' incidence
92 if (! isempty (cidx))
93 Fj(cidx) = -1;
94 endif
79 didx = find(d(jj,:) != 0); 95 didx = find(d(jj,:) != 0);
80 if(!isempty(didx)) Fj(n+didx) = 1; endif # D' incidence 96 ## D' incidence
81 Fj(n+m+jj) = length(cidx) + length(didx); # F2 is diagonal 97 if (! isempty (didx))
98 Fj(n+didx) = 1;
99 endif
100 Fj(n+m+jj) = length (cidx) + length (didx); # F2 is diagonal
82 y = y + x(n+m+jj)*Fj; # multiply by corresponding entry of x 101 y = y + x(n+m+jj)*Fj; # multiply by corresponding entry of x
83 endfor 102 endfor
84 103
85 endfunction 104 endfunction