|
3431
|
1 ## Copyright (C) 1996, 1998 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 |
|
7016
|
6 ## under the terms of the GNU General Public License as published by |
|
|
7 ## the Free Software Foundation; either version 3 of the License, or (at |
|
|
8 ## your option) any later version. |
|
3431
|
9 ## |
|
7016
|
10 ## Octave is distributed in the hope that it will be useful, but |
|
|
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
|
|
13 ## General Public License for more details. |
|
3431
|
14 ## |
|
|
15 ## You should have received a copy of the GNU General Public License |
|
7016
|
16 ## along with Octave; see the file COPYING. If not, see |
|
|
17 ## <http://www.gnu.org/licenses/>. |
|
3431
|
18 |
|
|
19 ## -*- texinfo -*- |
|
5016
|
20 ## @deftypefn {Function File} {@var{zz} =} zginit (@var{a}, @var{b}, @var{c}, @var{d}) |
|
|
21 ## Construct right hand side vector @var{zz} |
|
3431
|
22 ## for the zero-computation generalized eigenvalue problem |
|
5016
|
23 ## balancing procedure. Called by @command{zgepbal}. |
|
3431
|
24 ## @end deftypefn |
|
|
25 |
|
|
26 ## References: |
|
|
27 ## ZGEP: Hodel, "Computation of Zeros with Balancing," 1992, submitted to LAA |
|
|
28 ## Generalized CG: Golub and Van Loan, "Matrix Computations, 2nd ed" 1989 |
|
|
29 |
|
|
30 ## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu> |
|
|
31 ## Created: July 24, 1992 |
|
|
32 ## Conversion to Octave by R. Bruce Tenison, July 3, 1994 |
|
|
33 |
|
|
34 function zz = zginit (a, b, c, d) |
|
|
35 |
|
|
36 [nn,mm] = size(b); |
|
|
37 [pp,mm] = size(d); |
|
|
38 |
|
|
39 nmp = nn+mm+pp; |
|
|
40 |
|
|
41 ## set up log vector zz |
|
|
42 zz = zeros(nmp,1); |
|
|
43 |
|
|
44 ## zz part 1: |
|
|
45 for i=1:nn |
|
|
46 ## nonzero off diagonal entries of a |
|
|
47 if(nn > 1) |
|
|
48 nidx = complement(i,1:nn); |
|
|
49 a_row_i = a(i,nidx); a_col_i = a(nidx,i); |
|
|
50 arnz = a_row_i(find(a_row_i != 0)); acnz = a_col_i(find(a_col_i != 0)); |
|
|
51 else |
|
|
52 arnz = acnz = []; |
|
|
53 endif |
|
|
54 |
|
|
55 ## row of b |
|
|
56 bidx = find(b(i,:) != 0); |
|
|
57 b_row_i = b(i,bidx); |
|
|
58 |
|
|
59 ## column of c |
|
|
60 cidx = find(c(:,i) != 0); |
|
|
61 c_col_i = c(cidx,i); |
|
|
62 |
|
|
63 ## sum the entries |
|
|
64 zz(i) = sum(log(abs(acnz))) - sum(log(abs(arnz))) ... |
|
|
65 - sum(log(abs(b_row_i))) + sum(log(abs(c_col_i))); |
|
|
66 endfor |
|
|
67 |
|
|
68 ## zz part 2: |
|
|
69 bd = [b;d]; |
|
|
70 for i=1:mm |
|
|
71 i1 = i+nn; |
|
|
72 |
|
|
73 ## column of [b;d] |
|
|
74 bdidx = find(bd(:,i) != 0); |
|
|
75 bd_col_i = bd(bdidx,i); |
|
|
76 zz(i1) = sum(log(abs(bd_col_i))); |
|
|
77 endfor |
|
|
78 |
|
|
79 ## zz part 3: |
|
|
80 cd = [c, d]; |
|
|
81 for i=1:pp |
|
|
82 i1 = i+nn+mm; |
|
|
83 cdidx = find(cd(i,:) != 0); |
|
|
84 cd_row_i = cd(i,cdidx); |
|
|
85 zz(i1) = -sum(log(abs(cd_row_i))); |
|
|
86 endfor |
|
|
87 |
|
|
88 ## now set zz as log base 2 |
|
|
89 zz = zz*(1/log(2)); |
|
|
90 endfunction |
