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