Mercurial > hg > octave-lyh

comparison src/balance.cc @ 636:fae2bd91c027

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1994-08-23 18:39:50 by jwe]
author jwe
date Tue, 23 Aug 1994 18:39:50 +0000
parents 14b2a186a5c0
children 0a81458ef677
comparison
equal deleted inserted replaced
635:5338832d2cf6 636:fae2bd91c027
37 37
38 #include "tree-const.h" 38 #include "tree-const.h"
39 #include "user-prefs.h" 39 #include "user-prefs.h"
40 #include "gripes.h" 40 #include "gripes.h"
41 #include "error.h" 41 #include "error.h"
42 #include "utils.h"
42 #include "help.h" 43 #include "help.h"
43 #include "defun-dld.h" 44 #include "defun-dld.h"
44 45
45 DEFUN_DLD ("balance", Fbalance, Sbalance, 4, 4, 46 DEFUN_DLD ("balance", Fbalance, Sbalance, 4, 4,
46 "AA = balance (A [, OPT]) or [[DD,] AA] = balance (A [, OPT])\n\ 47 "AA = balance (A [, OPT]) or [[DD,] AA] = balance (A [, OPT])\n\
72 } 73 }
73 74
74 char *bal_job; 75 char *bal_job;
75 int my_nargin; // # args w/o optional string arg 76 int my_nargin; // # args w/o optional string arg
76 77
77 // determine if balancing option is listed 78 // Determine if balancing option is listed. Set my_nargin to the
78 // set my_nargin to the number of matrix inputs 79 // number of matrix inputs.
80
79 if (args(nargin-1).is_string ()) 81 if (args(nargin-1).is_string ())
80 { 82 {
81 bal_job = args(nargin-1).string_value (); 83 bal_job = args(nargin-1).string_value ();
82 my_nargin = nargin-2; 84 my_nargin = nargin-2;
83 } 85 }
85 { 87 {
86 bal_job = "B"; 88 bal_job = "B";
87 my_nargin = nargin-1; 89 my_nargin = nargin-1;
88 } 90 }
89 91
90 tree_constant arg = args(1).make_numeric (); 92 tree_constant arg_a = args(1);
91 int a_nr = arg.rows (); 93
92 int a_nc = arg.columns (); 94 int a_nr = arg_a.rows ();
95 int a_nc = arg_a.columns ();
93 96
94 // Check argument 1 dimensions. 97 // Check argument 1 dimensions.
95 98
96 if (a_nr == 0 || a_nc == 0) 99 if (empty_arg ("balance", a_nr, a_nc) < 0)
97 { 100 return retval;
98 int flag = user_pref.propagate_empty_matrices;
99 if (flag != 0)
100 {
101 if (flag < 0)
102 warning ("balance: argument is empty matrix");
103
104 Matrix m;
105 retval.resize (2);
106 retval(0) = m;
107 retval(1) = m;
108 }
109 else
110 error ("balance: empty matrix is invalid as argument");
111
112 return retval;
113 }
114 101
115 if (a_nr != a_nc) 102 if (a_nr != a_nc)
116 { 103 {
117 gripe_square_matrix_required ("balance"); 104 gripe_square_matrix_required ("balance");
118 return retval; 105 return retval;
120 107
121 // Extract argument 1 parameter for both AEP and GEP. 108 // Extract argument 1 parameter for both AEP and GEP.
122 109
123 Matrix aa; 110 Matrix aa;
124 ComplexMatrix caa; 111 ComplexMatrix caa;
125 if (arg.is_complex_type ()) 112 if (arg_a.is_complex_type ())
126 { 113 caa = arg_a.complex_matrix_value ();
127 if (arg.is_matrix_type ())
128 caa=arg.complex_matrix_value ();
129 else
130 {
131 caa.resize (1, 1);
132 caa.elem (0, 0) = arg.complex_value ();
133 }
134 }
135 else 114 else
136 { 115 aa = arg_a.matrix_value ();
137 if (arg.is_matrix_type ()) 116
138 aa = arg.matrix_value (); 117 if (error_state)
139 else 118 return retval;
140 {
141 double d = arg.double_value ();
142 aa.resize (1, 1);
143 aa.elem (0, 0) = d;
144 }
145 }
146 119
147 // Treat AEP/ GEP cases. 120 // Treat AEP/ GEP cases.
148 121
149 switch (my_nargin) 122 switch (my_nargin)
150 { 123 {
151 case 1: 124 case 1:
152 125
153 // Algebraic eigenvalue problem. 126 // Algebraic eigenvalue problem.
154 127
155 retval.resize (nargout ? nargout : 1); 128 if (arg_a.is_complex_type ())
156 if (arg.is_complex_type ())
157 { 129 {
158 ComplexAEPBALANCE result (caa, bal_job); 130 ComplexAEPBALANCE result (caa, bal_job);
159 131
160 if (nargout == 0 || nargout == 1) 132 if (nargout == 0 || nargout == 1)
161 retval(0) = result.balanced_matrix (); 133 retval(0) = result.balanced_matrix ();
162 else 134 else
163 { 135 {
136 retval(1) = result.balanced_matrix ();
164 retval(0) = result.balancing_matrix (); 137 retval(0) = result.balancing_matrix ();
165 retval(1) = result.balanced_matrix ();
166 } 138 }
167 } 139 }
168 else 140 else
169 { 141 {
170 AEPBALANCE result (aa, bal_job); 142 AEPBALANCE result (aa, bal_job);
171 143
172 if (nargout == 0 || nargout == 1) 144 if (nargout == 0 || nargout == 1)
173 retval(0) = result.balanced_matrix (); 145 retval(0) = result.balanced_matrix ();
174 else 146 else
175 { 147 {
148 retval(1) = result.balanced_matrix ();
176 retval(0) = result.balancing_matrix (); 149 retval(0) = result.balancing_matrix ();
177 retval(1) = result.balanced_matrix ();
178 } 150 }
179 } 151 }
180 break; 152 break;
153
181 case 2: 154 case 2:
182 155 {
183 // Generalized eigenvalue problem. 156 // Generalized eigenvalue problem.
184 157
185 {
186 retval.resize (nargout ? nargout : 1);
187
188 // 1st we have to check argument 2 dimensions and type... 158 // 1st we have to check argument 2 dimensions and type...
189 159
190 tree_constant brg = args(2).make_numeric (); 160 tree_constant arg_b = args(2);
191 int b_nr = brg.rows (); 161
192 int b_nc = brg.columns (); 162 int b_nr = arg_b.rows ();
163 int b_nc = arg_b.columns ();
193 164
194 // Check argument 2 dimensions -- must match arg 1. 165 // Check argument 2 dimensions -- must match arg 1.
195 166
196 if ((b_nr != b_nc) || (b_nr != a_nr)) 167 if (b_nr != b_nc || b_nr != a_nr)
197 { 168 {
198 gripe_nonconformant (); 169 gripe_nonconformant ();
199 return retval; 170 return retval;
200 } 171 }
201 172
202 // Now, extract the second matrix... 173 // Now, extract the second matrix...
203 // Extract argument 1 parameter for both AEP and GEP. 174 // Extract argument 1 parameter for both AEP and GEP.
204 175
205 Matrix bb; 176 Matrix bb;
206 ComplexMatrix cbb; 177 ComplexMatrix cbb;
207 if (brg.is_complex_type ()) 178 if (arg_b.is_complex_type ())
208 { 179 cbb = arg_b.complex_matrix_value ();
209 if (brg.is_matrix_type ())
210 cbb = brg.complex_matrix_value ();
211 else
212 {
213 cbb.resize (1, 1);
214 cbb.elem (0, 0) = brg.complex_value ();
215 }
216 }
217 else 180 else
218 { 181 bb = arg_b.matrix_value ();
219 if (brg.is_matrix_type ()) 182
220 bb = brg.matrix_value (); 183 if (error_state)
221 else 184 return retval;
222 { 185
223 double d = brg.double_value ();
224 bb.resize (1, 1);
225 bb.elem (0, 0) = d;
226 }
227 }
228
229 // Both matrices loaded, now let's check what kind of arithmetic: 186 // Both matrices loaded, now let's check what kind of arithmetic:
230 187
231 if (arg.is_complex_type () || brg.is_complex_type ()) 188 if (arg_a.is_complex_type () || arg_b.is_complex_type ())
232 { 189 {
233 if (arg.is_real_type ()) 190 if (arg_a.is_real_type ())
234 caa = aa; 191 caa = aa;
235 else if (brg.is_real_type ()) 192
193 if (arg_b.is_real_type ())
236 cbb = bb; 194 cbb = bb;
237 195
238 // Compute magnitudes of elements for balancing purposes. 196 // Compute magnitudes of elements for balancing purposes.
239 // Surely there's a function I can call someplace! 197 // Surely there's a function I can call someplace!
240 198
241 for (int i = 0; i < a_nr; i++) 199 for (int i = 0; i < a_nr; i++)
242 for (int j = 0; j < a_nr; j++) 200 for (int j = 0; j < a_nc; j++)
243 { 201 {
244 aa.elem (i, j) = abs (caa.elem (i, j)); 202 aa.elem (i, j) = abs (caa.elem (i, j));
245 bb.elem (i, j) = abs (cbb.elem (i, j)); 203 bb.elem (i, j) = abs (cbb.elem (i, j));
246 } 204 }
247 } 205 }
248 206
249 GEPBALANCE result (aa, bb, bal_job); 207 GEPBALANCE result (aa, bb, bal_job);
250 208
251 if (arg.is_complex_type () || brg.is_complex_type ()) 209 if (arg_a.is_complex_type () || arg_b.is_complex_type ())
252 { 210 {
253 caa = result.left_balancing_matrix () * caa 211 caa = result.left_balancing_matrix () * caa
254 * result.right_balancing_matrix (); 212 * result.right_balancing_matrix ();
255 213
256 cbb = result.left_balancing_matrix () * cbb 214 cbb = result.left_balancing_matrix () * cbb
261 case 0: 219 case 0:
262 case 1: 220 case 1:
263 warning ("balance: should use two output arguments"); 221 warning ("balance: should use two output arguments");
264 retval(0) = caa; 222 retval(0) = caa;
265 break; 223 break;
224
266 case 2: 225 case 2:
226 retval(1) = cbb;
267 retval(0) = caa; 227 retval(0) = caa;
268 retval(1) = cbb; 228 break;
269 break; 229
270 case 4: 230 case 4:
231 retval(3) = cbb;
232 retval(2) = caa;
233 retval(1) = result.right_balancing_matrix ();
271 retval(0) = result.left_balancing_matrix (); 234 retval(0) = result.left_balancing_matrix ();
272 retval(1) = result.right_balancing_matrix (); 235 break;
273 retval(2) = caa; 236
274 retval(3) = cbb;
275 break;
276 default: 237 default:
277 error ("balance: invalid number of output arguments"); 238 error ("balance: invalid number of output arguments");
278 break; 239 break;
279 } 240 }
280 } 241 }
285 case 0: 246 case 0:
286 case 1: 247 case 1:
287 warning ("balance: should use two output arguments"); 248 warning ("balance: should use two output arguments");
288 retval(0) = result.balanced_a_matrix (); 249 retval(0) = result.balanced_a_matrix ();
289 break; 250 break;
251
290 case 2: 252 case 2:
253 retval(1) = result.balanced_b_matrix ();
291 retval(0) = result.balanced_a_matrix (); 254 retval(0) = result.balanced_a_matrix ();
292 retval(1) = result.balanced_b_matrix (); 255 break;
293 break; 256
294 case 4: 257 case 4:
258 retval(3) = result.balanced_b_matrix ();
259 retval(2) = result.balanced_a_matrix ();
260 retval(1) = result.right_balancing_matrix ();
295 retval(0) = result.left_balancing_matrix (); 261 retval(0) = result.left_balancing_matrix ();
296 retval(1) = result.right_balancing_matrix (); 262 break;
297 retval(2) = result.balanced_a_matrix (); 263
298 retval(3) = result.balanced_b_matrix ();
299 break;
300 default: 264 default:
301 error ("balance: invalid number of output arguments"); 265 error ("balance: invalid number of output arguments");
302 break; 266 break;
303 } 267 }
304 } 268 }
305 } 269 }
306 break; 270 break;
271
307 default: 272 default:
308 error ("balance requires one (AEP) or two (GEP) numeric arguments"); 273 error ("balance requires one (AEP) or two (GEP) numeric arguments");
309 break; 274 break;
310 } 275 }
276
311 return retval; 277 return retval;
312 } 278 }
313 279
314 /* 280 /*
315 ;;; Local Variables: *** 281 ;;; Local Variables: ***