|
1
|
1 // The constants for the tree class. -*- C++ -*- |
|
|
2 /* |
|
|
3 |
|
296
|
4 Copyright (C) 1992, 1993, 1994 John W. Eaton |
|
1
|
5 |
|
|
6 This file is part of Octave. |
|
|
7 |
|
|
8 Octave is free software; you can redistribute it and/or modify it |
|
|
9 under the terms of the GNU General Public License as published by the |
|
|
10 Free Software Foundation; either version 2, or (at your option) any |
|
|
11 later version. |
|
|
12 |
|
|
13 Octave is distributed in the hope that it will be useful, but WITHOUT |
|
|
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
|
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
|
16 for more details. |
|
|
17 |
|
|
18 You should have received a copy of the GNU General Public License |
|
|
19 along with Octave; see the file COPYING. If not, write to the Free |
|
|
20 Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. |
|
|
21 |
|
|
22 */ |
|
|
23 |
|
240
|
24 #ifdef HAVE_CONFIG_H |
|
|
25 #include "config.h" |
|
1
|
26 #endif |
|
|
27 |
|
455
|
28 #if defined (__GNUG__) |
|
|
29 #pragma implementation |
|
|
30 #endif |
|
|
31 |
|
493
|
32 #include <strstream.h> |
|
1
|
33 #include <iostream.h> |
|
493
|
34 #include <fstream.h> |
|
|
35 #include <string.h> |
|
|
36 #include <ctype.h> |
|
|
37 |
|
|
38 #include "EIG.h" |
|
1
|
39 |
|
493
|
40 #include "unwind-prot.h" |
|
|
41 #include "tree-const.h" |
|
|
42 #include "user-prefs.h" |
|
1
|
43 #include "variables.h" |
|
493
|
44 #include "octave.h" |
|
1
|
45 #include "error.h" |
|
|
46 #include "gripes.h" |
|
455
|
47 #include "input.h" |
|
|
48 #include "octave-hist.h" |
|
493
|
49 #include "pager.h" |
|
|
50 #include "utils.h" |
|
455
|
51 #include "parse.h" |
|
|
52 #include "lex.h" |
|
|
53 |
|
493
|
54 #ifndef MAX |
|
|
55 #define MAX(a,b) ((a) > (b) ? (a) : (b)) |
|
1
|
56 #endif |
|
|
57 |
|
493
|
58 #ifndef MIN |
|
|
59 #define MIN(a,b) ((a) < (b) ? (a) : (b)) |
|
|
60 #endif |
|
1
|
61 |
|
|
62 tree_constant::~tree_constant (void) |
|
|
63 { |
|
|
64 #if defined (MDEBUG) |
|
|
65 cerr << "~tree_constant: rep: " << rep |
|
|
66 << " rep->count: " << rep->count << "\n"; |
|
|
67 #endif |
|
|
68 |
|
|
69 if (--rep->count <= 0) |
|
|
70 { |
|
|
71 delete rep; |
|
|
72 rep = (tree_constant_rep *) NULL; |
|
|
73 } |
|
|
74 } |
|
|
75 |
|
|
76 #if defined (MDEBUG) |
|
|
77 void * |
|
|
78 tree_constant::operator new (size_t size) |
|
|
79 { |
|
|
80 tree_constant *p = ::new tree_constant; |
|
|
81 cerr << "tree_constant::new(): " << p << "\n"; |
|
|
82 return p; |
|
|
83 } |
|
|
84 |
|
|
85 void |
|
|
86 tree_constant::operator delete (void *p, size_t size) |
|
|
87 { |
|
|
88 cerr << "tree_constant::delete(): " << p << "\n"; |
|
|
89 ::delete p; |
|
|
90 } |
|
|
91 #endif |
|
|
92 |
|
|
93 /* |
|
96
|
94 * Construct return vector of empty matrices. Return empty matrices |
|
|
95 * and/or gripe when appropriate. |
|
|
96 */ |
|
500
|
97 Octave_object |
|
164
|
98 vector_of_empties (int nargout, const char *fcn_name) |
|
96
|
99 { |
|
500
|
100 Octave_object retval; |
|
96
|
101 |
|
|
102 // Got an empty argument, check if should gripe/return empty values. |
|
|
103 |
|
|
104 int flag = user_pref.propagate_empty_matrices; |
|
|
105 if (flag != 0) |
|
|
106 { |
|
|
107 if (flag < 0) |
|
|
108 gripe_empty_arg (fcn_name, 0); |
|
|
109 |
|
|
110 Matrix m; |
|
500
|
111 retval.resize (nargout ? nargout : 1); |
|
96
|
112 for (int i = 0; i < nargout; i++) |
|
500
|
113 retval(i) = tree_constant (m); |
|
96
|
114 } |
|
|
115 else |
|
|
116 gripe_empty_arg (fcn_name, 1); |
|
|
117 |
|
|
118 return retval; |
|
|
119 } |
|
|
120 |
|
493
|
121 Matrix |
|
|
122 max (const Matrix& a, const Matrix& b) |
|
455
|
123 { |
|
493
|
124 int nr = a.rows (); |
|
|
125 int nc = a.columns (); |
|
|
126 if (nr != b.rows () || nc != b.columns ()) |
|
|
127 { |
|
|
128 error ("two-arg max expecting args of same size"); |
|
|
129 return Matrix (); |
|
|
130 } |
|
455
|
131 |
|
493
|
132 Matrix result (nr, nc); |
|
|
133 |
|
|
134 for (int j = 0; j < nc; j++) |
|
|
135 for (int i = 0; i < nr; i++) |
|
|
136 { |
|
|
137 double a_elem = a.elem (i, j); |
|
|
138 double b_elem = b.elem (i, j); |
|
|
139 result.elem (i, j) = MAX (a_elem, b_elem); |
|
|
140 } |
|
|
141 |
|
|
142 return result; |
|
455
|
143 } |
|
|
144 |
|
493
|
145 ComplexMatrix |
|
|
146 max (const ComplexMatrix& a, const ComplexMatrix& b) |
|
455
|
147 { |
|
493
|
148 int nr = a.rows (); |
|
|
149 int nc = a.columns (); |
|
|
150 if (nr != b.rows () || nc != b.columns ()) |
|
455
|
151 { |
|
493
|
152 error ("two-arg max expecting args of same size"); |
|
|
153 return ComplexMatrix (); |
|
455
|
154 } |
|
|
155 |
|
493
|
156 ComplexMatrix result (nr, nc); |
|
|
157 |
|
|
158 for (int j = 0; j < nc; j++) |
|
|
159 for (int i = 0; i < nr; i++) |
|
|
160 { |
|
|
161 double abs_a_elem = abs (a.elem (i, j)); |
|
|
162 double abs_b_elem = abs (b.elem (i, j)); |
|
|
163 if (abs_a_elem > abs_b_elem) |
|
|
164 result.elem (i, j) = a.elem (i, j); |
|
|
165 else |
|
|
166 result.elem (i, j) = b.elem (i, j); |
|
|
167 } |
|
|
168 |
|
|
169 return result; |
|
455
|
170 } |
|
|
171 |
|
493
|
172 Matrix |
|
|
173 min (const Matrix& a, const Matrix& b) |
|
455
|
174 { |
|
493
|
175 int nr = a.rows (); |
|
|
176 int nc = a.columns (); |
|
|
177 if (nr != b.rows () || nc != b.columns ()) |
|
455
|
178 { |
|
493
|
179 error ("two-arg min expecting args of same size"); |
|
|
180 return Matrix (); |
|
|
181 } |
|
|
182 |
|
|
183 Matrix result (nr, nc); |
|
455
|
184 |
|
493
|
185 for (int j = 0; j < nc; j++) |
|
|
186 for (int i = 0; i < nr; i++) |
|
|
187 { |
|
|
188 double a_elem = a.elem (i, j); |
|
|
189 double b_elem = b.elem (i, j); |
|
|
190 result.elem (i, j) = MIN (a_elem, b_elem); |
|
|
191 } |
|
455
|
192 |
|
493
|
193 return result; |
|
455
|
194 } |
|
|
195 |
|
493
|
196 ComplexMatrix |
|
|
197 min (const ComplexMatrix& a, const ComplexMatrix& b) |
|
455
|
198 { |
|
493
|
199 int nr = a.rows (); |
|
|
200 int nc = a.columns (); |
|
|
201 if (nr != b.rows () || nc != b.columns ()) |
|
455
|
202 { |
|
493
|
203 error ("two-arg min expecting args of same size"); |
|
|
204 return ComplexMatrix (); |
|
|
205 } |
|
455
|
206 |
|
493
|
207 ComplexMatrix result (nr, nc); |
|
455
|
208 |
|
493
|
209 for (int j = 0; j < nc; j++) |
|
|
210 for (int i = 0; i < nr; i++) |
|
|
211 { |
|
|
212 double abs_a_elem = abs (a.elem (i, j)); |
|
|
213 double abs_b_elem = abs (b.elem (i, j)); |
|
|
214 if (abs_a_elem < abs_b_elem) |
|
|
215 result.elem (i, j) = a.elem (i, j); |
|
|
216 else |
|
|
217 result.elem (i, j) = b.elem (i, j); |
|
|
218 } |
|
455
|
219 |
|
493
|
220 return result; |
|
455
|
221 } |
|
|
222 |
|
493
|
223 static void |
|
|
224 get_dimensions (const tree_constant& a, const char *warn_for, |
|
|
225 int& nr, int& nc) |
|
455
|
226 { |
|
493
|
227 tree_constant tmpa = a.make_numeric (); |
|
455
|
228 |
|
493
|
229 if (tmpa.is_scalar_type ()) |
|
455
|
230 { |
|
493
|
231 double tmp = tmpa.double_value (); |
|
|
232 nr = nc = NINT (tmp); |
|
455
|
233 } |
|
|
234 else |
|
|
235 { |
|
493
|
236 nr = tmpa.rows (); |
|
|
237 nc = tmpa.columns (); |
|
|
238 |
|
|
239 if ((nr == 1 && nc == 2) || (nr == 2 && nc == 1)) |
|
|
240 { |
|
|
241 ColumnVector v = tmpa.to_vector (); |
|
|
242 |
|
|
243 nr = NINT (v.elem (0)); |
|
|
244 nc = NINT (v.elem (1)); |
|
|
245 } |
|
|
246 else |
|
|
247 warning ("%s (A): use %s (size (A)) instead", warn_for, warn_for); |
|
|
248 } |
|
|
249 |
|
|
250 check_dimensions (nr, nc, warn_for); // May set error_state. |
|
|
251 } |
|
|
252 |
|
|
253 static void |
|
|
254 get_dimensions (const tree_constant& a, const tree_constant& b, |
|
|
255 const char *warn_for, int& nr, int& nc) |
|
|
256 { |
|
|
257 tree_constant tmpa = a.make_numeric (); |
|
|
258 tree_constant tmpb = b.make_numeric (); |
|
|
259 |
|
|
260 if (tmpa.is_scalar_type () && tmpb.is_scalar_type ()) |
|
|
261 { |
|
|
262 nr = NINT (tmpa.double_value ()); |
|
|
263 nc = NINT (tmpb.double_value ()); |
|
|
264 |
|
|
265 check_dimensions (nr, nc, warn_for); // May set error_state. |
|
455
|
266 } |
|
493
|
267 else |
|
|
268 error ("%s: expecting two scalar arguments", warn_for); |
|
|
269 } |
|
|
270 |
|
|
271 tree_constant |
|
|
272 fill_matrix (const tree_constant& a, double val, const char *warn_for) |
|
|
273 { |
|
|
274 int nr, nc; |
|
|
275 get_dimensions (a, warn_for, nr, nc); |
|
|
276 |
|
|
277 if (error_state) |
|
|
278 return tree_constant (); |
|
|
279 |
|
|
280 Matrix m (nr, nc, val); |
|
|
281 |
|
|
282 return tree_constant (m); |
|
|
283 } |
|
|
284 |
|
|
285 tree_constant |
|
|
286 fill_matrix (const tree_constant& a, const tree_constant& b, |
|
|
287 double val, const char *warn_for) |
|
|
288 { |
|
|
289 int nr, nc; |
|
|
290 get_dimensions (a, b, warn_for, nr, nc); // May set error_state. |
|
|
291 |
|
|
292 if (error_state) |
|
|
293 return tree_constant (); |
|
|
294 |
|
|
295 Matrix m (nr, nc, val); |
|
|
296 |
|
|
297 return tree_constant (m); |
|
455
|
298 } |
|
|
299 |
|
493
|
300 tree_constant |
|
|
301 identity_matrix (const tree_constant& a) |
|
455
|
302 { |
|
493
|
303 int nr, nc; |
|
|
304 get_dimensions (a, "eye", nr, nc); // May set error_state. |
|
|
305 |
|
|
306 if (error_state) |
|
|
307 return tree_constant (); |
|
|
308 |
|
|
309 Matrix m (nr, nc, 0.0); |
|
455
|
310 |
|
493
|
311 if (nr > 0 && nc > 0) |
|
|
312 { |
|
|
313 int n = MIN (nr, nc); |
|
|
314 for (int i = 0; i < n; i++) |
|
|
315 m.elem (i, i) = 1.0; |
|
|
316 } |
|
455
|
317 |
|
493
|
318 return tree_constant (m); |
|
|
319 } |
|
455
|
320 |
|
493
|
321 tree_constant |
|
|
322 identity_matrix (const tree_constant& a, const tree_constant& b) |
|
|
323 { |
|
|
324 int nr, nc; |
|
|
325 get_dimensions (a, b, "eye", nr, nc); // May set error_state. |
|
|
326 |
|
|
327 if (error_state) |
|
|
328 return tree_constant (); |
|
|
329 |
|
|
330 Matrix m (nr, nc, 0.0); |
|
|
331 |
|
|
332 if (nr > 0 && nc > 0) |
|
455
|
333 { |
|
493
|
334 int n = MIN (nr, nc); |
|
|
335 for (int i = 0; i < n; i++) |
|
|
336 m.elem (i, i) = 1.0; |
|
|
337 } |
|
|
338 |
|
|
339 return tree_constant (m); |
|
|
340 } |
|
455
|
341 |
|
511
|
342 static Octave_object |
|
|
343 find_to_fortran_idx (const ColumnVector i_idx, const ColumnVector j_idx, |
|
|
344 const tree_constant& val, int nr, int nc, int nargout) |
|
493
|
345 { |
|
511
|
346 Octave_object retval (nargout); |
|
493
|
347 |
|
511
|
348 switch (nargout) |
|
493
|
349 { |
|
511
|
350 case 1: |
|
|
351 { |
|
|
352 int count = i_idx.length (); |
|
|
353 ColumnVector tmp (count); |
|
|
354 for (int i = 0; i < count; i++) |
|
|
355 tmp (i) = nr * (j_idx (i) - 1.0) + i_idx (i); |
|
512
|
356 retval(0) = tree_constant (tmp, 1); |
|
|
357 // If you want this to work more like Matlab, use the following line |
|
|
358 // instead of the previous one. |
|
|
359 // retval(0) = tree_constant (tmp, (nr != 1)); |
|
511
|
360 } |
|
|
361 break; |
|
|
362 case 3: |
|
513
|
363 retval(2) = val; |
|
511
|
364 case 2: |
|
513
|
365 retval(0) = tree_constant (i_idx, 1); |
|
512
|
366 // If you want this to work more like Matlab, use the following line |
|
|
367 // instead of the previous one. |
|
513
|
368 // retval(0) = tree_constant (i_idx, (nr != 1)); |
|
511
|
369 retval(1) = tree_constant (j_idx, 1); |
|
|
370 break; |
|
|
371 default: |
|
|
372 panic_impossible (); |
|
|
373 break; |
|
493
|
374 } |
|
511
|
375 |
|
|
376 return retval; |
|
455
|
377 } |
|
|
378 |
|
511
|
379 static Octave_object |
|
|
380 find_nonzero_elem_idx (const Matrix& m, int nargout) |
|
455
|
381 { |
|
493
|
382 int count = 0; |
|
|
383 int m_nr = m.rows (); |
|
|
384 int m_nc = m.columns (); |
|
455
|
385 |
|
511
|
386 int i, j; |
|
|
387 for (j = 0; j < m_nc; j++) |
|
|
388 for (i = 0; i < m_nr; i++) |
|
|
389 if (m.elem (i, j) != 0.0) |
|
|
390 count++; |
|
455
|
391 |
|
493
|
392 Matrix result; |
|
511
|
393 Octave_object retval (nargout, result); |
|
455
|
394 |
|
493
|
395 if (count == 0) |
|
511
|
396 return retval; |
|
|
397 |
|
|
398 ColumnVector i_idx (count); |
|
|
399 ColumnVector j_idx (count); |
|
|
400 ColumnVector v (count); |
|
493
|
401 |
|
511
|
402 count = 0; |
|
|
403 for (j = 0; j < m_nc; j++) |
|
|
404 for (i = 0; i < m_nr; i++) |
|
|
405 { |
|
|
406 double d = m.elem (i, j); |
|
|
407 if (d != 0.0) |
|
493
|
408 { |
|
511
|
409 i_idx (count) = i + 1; |
|
|
410 j_idx (count) = j + 1; |
|
|
411 v (count) = d; |
|
493
|
412 count++; |
|
|
413 } |
|
511
|
414 } |
|
|
415 |
|
513
|
416 tree_constant tmp (v, 1); |
|
|
417 return find_to_fortran_idx (i_idx, j_idx, tmp, m_nr, m_nc, nargout); |
|
455
|
418 } |
|
|
419 |
|
511
|
420 static Octave_object |
|
|
421 find_nonzero_elem_idx (const ComplexMatrix& m, int nargout) |
|
455
|
422 { |
|
511
|
423 int count = 0; |
|
|
424 int m_nr = m.rows (); |
|
|
425 int m_nc = m.columns (); |
|
|
426 |
|
|
427 int i, j; |
|
|
428 for (j = 0; j < m_nc; j++) |
|
|
429 for (i = 0; i < m_nr; i++) |
|
|
430 if (m.elem (i, j) != 0.0) |
|
|
431 count++; |
|
|
432 |
|
|
433 Matrix result; |
|
|
434 Octave_object retval (nargout, result); |
|
|
435 |
|
|
436 if (count == 0) |
|
|
437 return retval; |
|
|
438 |
|
|
439 ColumnVector i_idx (count); |
|
|
440 ColumnVector j_idx (count); |
|
|
441 ComplexColumnVector v (count); |
|
|
442 |
|
|
443 count = 0; |
|
|
444 for (j = 0; j < m_nc; j++) |
|
|
445 for (i = 0; i < m_nr; i++) |
|
|
446 { |
|
|
447 Complex c = m.elem (i, j); |
|
|
448 if (c != 0.0) |
|
|
449 { |
|
|
450 i_idx (count) = i; |
|
|
451 j_idx (count) = j; |
|
|
452 v (count) = c; |
|
|
453 count++; |
|
|
454 } |
|
|
455 } |
|
|
456 |
|
513
|
457 tree_constant tmp (v, 1); |
|
|
458 return find_to_fortran_idx (i_idx, j_idx, tmp, m_nr, m_nc, nargout); |
|
511
|
459 } |
|
|
460 |
|
|
461 Octave_object |
|
|
462 find_nonzero_elem_idx (const tree_constant& a, int nargout) |
|
|
463 { |
|
|
464 Matrix result; |
|
|
465 |
|
|
466 nargout = (nargout == 0) ? 1 : nargout; |
|
|
467 Octave_object retval (nargout, result); |
|
455
|
468 |
|
493
|
469 tree_constant tmp = a.make_numeric (); |
|
455
|
470 |
|
493
|
471 switch (tmp.const_type ()) |
|
455
|
472 { |
|
493
|
473 case tree_constant_rep::matrix_constant: |
|
|
474 { |
|
|
475 Matrix m = tmp.matrix_value (); |
|
511
|
476 return find_nonzero_elem_idx (m, nargout); |
|
493
|
477 } |
|
|
478 break; |
|
|
479 case tree_constant_rep::scalar_constant: |
|
|
480 { |
|
|
481 double d = tmp.double_value (); |
|
|
482 if (d != 0.0) |
|
511
|
483 { |
|
|
484 retval(0) = 1.0; |
|
|
485 if (nargout > 1) |
|
|
486 retval(1) = 1.0; |
|
|
487 if (nargout > 2) |
|
|
488 retval(2) = d; |
|
|
489 } |
|
493
|
490 } |
|
|
491 break; |
|
|
492 case tree_constant_rep::complex_matrix_constant: |
|
|
493 { |
|
|
494 ComplexMatrix m = tmp.complex_matrix_value (); |
|
511
|
495 return find_nonzero_elem_idx (m, nargout); |
|
493
|
496 } |
|
|
497 break; |
|
|
498 case tree_constant_rep::complex_scalar_constant: |
|
|
499 { |
|
|
500 Complex c = tmp.complex_value (); |
|
|
501 if (c != 0.0) |
|
511
|
502 { |
|
|
503 retval(0) = 1.0; |
|
|
504 if (nargout > 1) |
|
|
505 retval(1) = 1.0; |
|
|
506 if (nargout > 2) |
|
|
507 retval(2) = c; |
|
|
508 } |
|
493
|
509 } |
|
|
510 break; |
|
|
511 default: |
|
|
512 break; |
|
455
|
513 } |
|
493
|
514 return retval; |
|
455
|
515 } |
|
|
516 |
|
493
|
517 // XXX FIXME XXX -- the next two functions (and expm) should really be just |
|
|
518 // one... |
|
|
519 |
|
500
|
520 Octave_object |
|
493
|
521 matrix_log (const tree_constant& a) |
|
455
|
522 { |
|
500
|
523 Octave_object retval (1); |
|
455
|
524 |
|
493
|
525 tree_constant tmp = a.make_numeric ();; |
|
|
526 |
|
|
527 if (tmp.rows () == 0 || tmp.columns () == 0) |
|
455
|
528 { |
|
493
|
529 int flag = user_pref.propagate_empty_matrices; |
|
|
530 if (flag != 0) |
|
455
|
531 { |
|
493
|
532 if (flag < 0) |
|
|
533 gripe_empty_arg ("logm", 0); |
|
|
534 Matrix m; |
|
500
|
535 retval(0) = tree_constant (m); |
|
493
|
536 return retval; |
|
455
|
537 } |
|
|
538 else |
|
493
|
539 gripe_empty_arg ("logm", 1); |
|
455
|
540 } |
|
|
541 |
|
493
|
542 switch (tmp.const_type ()) |
|
|
543 { |
|
|
544 case tree_constant_rep::matrix_constant: |
|
|
545 { |
|
|
546 Matrix m = tmp.matrix_value (); |
|
455
|
547 |
|
493
|
548 int nr = m.rows (); |
|
|
549 int nc = m.columns (); |
|
455
|
550 |
|
493
|
551 if (nr == 0 || nc == 0 || nr != nc) |
|
|
552 gripe_square_matrix_required ("logm"); |
|
|
553 else |
|
|
554 { |
|
|
555 EIG m_eig (m); |
|
|
556 ComplexColumnVector lambda (m_eig.eigenvalues ()); |
|
|
557 ComplexMatrix Q (m_eig.eigenvectors ()); |
|
455
|
558 |
|
493
|
559 for (int i = 0; i < nr; i++) |
|
|
560 { |
|
|
561 Complex elt = lambda.elem (i); |
|
|
562 if (imag (elt) == 0.0 && real (elt) > 0.0) |
|
|
563 lambda.elem (i) = log (real (elt)); |
|
|
564 else |
|
|
565 lambda.elem (i) = log (elt); |
|
|
566 } |
|
455
|
567 |
|
493
|
568 ComplexDiagMatrix D (lambda); |
|
|
569 ComplexMatrix result = Q * D * Q.inverse (); |
|
455
|
570 |
|
500
|
571 retval(0) = tree_constant (result); |
|
493
|
572 } |
|
455
|
573 } |
|
|
574 break; |
|
493
|
575 case tree_constant_rep::complex_matrix_constant: |
|
455
|
576 { |
|
493
|
577 ComplexMatrix m = tmp.complex_matrix_value (); |
|
455
|
578 |
|
493
|
579 int nr = m.rows (); |
|
|
580 int nc = m.columns (); |
|
|
581 |
|
|
582 if (nr == 0 || nc == 0 || nr != nc) |
|
|
583 gripe_square_matrix_required ("logm"); |
|
455
|
584 else |
|
|
585 { |
|
493
|
586 EIG m_eig (m); |
|
|
587 ComplexColumnVector lambda (m_eig.eigenvalues ()); |
|
|
588 ComplexMatrix Q (m_eig.eigenvectors ()); |
|
|
589 |
|
|
590 for (int i = 0; i < nr; i++) |
|
|
591 { |
|
|
592 Complex elt = lambda.elem (i); |
|
|
593 if (imag (elt) == 0.0 && real (elt) > 0.0) |
|
|
594 lambda.elem (i) = log (real (elt)); |
|
|
595 else |
|
|
596 lambda.elem (i) = log (elt); |
|
|
597 } |
|
|
598 |
|
|
599 ComplexDiagMatrix D (lambda); |
|
|
600 ComplexMatrix result = Q * D * Q.inverse (); |
|
|
601 |
|
500
|
602 retval(0) = tree_constant (result); |
|
455
|
603 } |
|
|
604 } |
|
|
605 break; |
|
493
|
606 case tree_constant_rep::scalar_constant: |
|
455
|
607 { |
|
493
|
608 double d = tmp.double_value (); |
|
|
609 if (d > 0.0) |
|
500
|
610 retval(0) = tree_constant (log (d)); |
|
455
|
611 else |
|
|
612 { |
|
493
|
613 Complex dtmp (d); |
|
500
|
614 retval(0) = tree_constant (log (dtmp)); |
|
455
|
615 } |
|
|
616 } |
|
|
617 break; |
|
493
|
618 case tree_constant_rep::complex_scalar_constant: |
|
455
|
619 { |
|
493
|
620 Complex c = tmp.complex_value (); |
|
500
|
621 retval(0) = tree_constant (log (c)); |
|
455
|
622 } |
|
|
623 break; |
|
|
624 default: |
|
|
625 break; |
|
|
626 } |
|
493
|
627 return retval; |
|
455
|
628 } |
|
|
629 |
|
500
|
630 Octave_object |
|
493
|
631 matrix_sqrt (const tree_constant& a) |
|
455
|
632 { |
|
500
|
633 Octave_object retval (1); |
|
455
|
634 |
|
493
|
635 tree_constant tmp = a.make_numeric ();; |
|
|
636 |
|
|
637 if (tmp.rows () == 0 || tmp.columns () == 0) |
|
455
|
638 { |
|
493
|
639 int flag = user_pref.propagate_empty_matrices; |
|
|
640 if (flag != 0) |
|
|
641 { |
|
|
642 if (flag < 0) |
|
|
643 gripe_empty_arg ("sqrtm", 0); |
|
|
644 Matrix m; |
|
500
|
645 retval(0) = tree_constant (m); |
|
493
|
646 return retval; |
|
|
647 } |
|
|
648 else |
|
|
649 gripe_empty_arg ("sqrtm", 1); |
|
|
650 } |
|
455
|
651 |
|
493
|
652 switch (tmp.const_type ()) |
|
|
653 { |
|
|
654 case tree_constant_rep::matrix_constant: |
|
455
|
655 { |
|
493
|
656 Matrix m = tmp.matrix_value (); |
|
455
|
657 |
|
493
|
658 int nr = m.rows (); |
|
|
659 int nc = m.columns (); |
|
455
|
660 |
|
493
|
661 if (nr == 0 || nc == 0 || nr != nc) |
|
|
662 gripe_square_matrix_required ("sqrtm"); |
|
455
|
663 else |
|
|
664 { |
|
493
|
665 EIG m_eig (m); |
|
|
666 ComplexColumnVector lambda (m_eig.eigenvalues ()); |
|
|
667 ComplexMatrix Q (m_eig.eigenvectors ()); |
|
455
|
668 |
|
493
|
669 for (int i = 0; i < nr; i++) |
|
|
670 { |
|
|
671 Complex elt = lambda.elem (i); |
|
|
672 if (imag (elt) == 0.0 && real (elt) > 0.0) |
|
|
673 lambda.elem (i) = sqrt (real (elt)); |
|
|
674 else |
|
|
675 lambda.elem (i) = sqrt (elt); |
|
|
676 } |
|
|
677 |
|
|
678 ComplexDiagMatrix D (lambda); |
|
|
679 ComplexMatrix result = Q * D * Q.inverse (); |
|
|
680 |
|
500
|
681 retval(0) = tree_constant (result); |
|
455
|
682 } |
|
|
683 } |
|
|
684 break; |
|
493
|
685 case tree_constant_rep::complex_matrix_constant: |
|
455
|
686 { |
|
493
|
687 ComplexMatrix m = tmp.complex_matrix_value (); |
|
455
|
688 |
|
493
|
689 int nr = m.rows (); |
|
|
690 int nc = m.columns (); |
|
|
691 |
|
|
692 if (nr == 0 || nc == 0 || nr != nc) |
|
|
693 gripe_square_matrix_required ("sqrtm"); |
|
|
694 else |
|
455
|
695 { |
|
493
|
696 EIG m_eig (m); |
|
|
697 ComplexColumnVector lambda (m_eig.eigenvalues ()); |
|
|
698 ComplexMatrix Q (m_eig.eigenvectors ()); |
|
|
699 |
|
|
700 for (int i = 0; i < nr; i++) |
|
|
701 { |
|
|
702 Complex elt = lambda.elem (i); |
|
|
703 if (imag (elt) == 0.0 && real (elt) > 0.0) |
|
|
704 lambda.elem (i) = sqrt (real (elt)); |
|
|
705 else |
|
|
706 lambda.elem (i) = sqrt (elt); |
|
|
707 } |
|
|
708 |
|
|
709 ComplexDiagMatrix D (lambda); |
|
|
710 ComplexMatrix result = Q * D * Q.inverse (); |
|
|
711 |
|
500
|
712 retval(0) = tree_constant (result); |
|
455
|
713 } |
|
493
|
714 } |
|
|
715 break; |
|
|
716 case tree_constant_rep::scalar_constant: |
|
|
717 { |
|
|
718 double d = tmp.double_value (); |
|
|
719 if (d > 0.0) |
|
500
|
720 retval(0) = tree_constant (sqrt (d)); |
|
455
|
721 else |
|
|
722 { |
|
493
|
723 Complex dtmp (d); |
|
500
|
724 retval(0) = tree_constant (sqrt (dtmp)); |
|
455
|
725 } |
|
493
|
726 } |
|
|
727 break; |
|
|
728 case tree_constant_rep::complex_scalar_constant: |
|
|
729 { |
|
|
730 Complex c = tmp.complex_value (); |
|
500
|
731 retval(0) = tree_constant (log (c)); |
|
493
|
732 } |
|
|
733 break; |
|
|
734 default: |
|
|
735 break; |
|
|
736 } |
|
|
737 return retval; |
|
|
738 } |
|
455
|
739 |
|
500
|
740 Octave_object |
|
506
|
741 column_max (const Octave_object& args, int nargout) |
|
493
|
742 { |
|
500
|
743 Octave_object retval; |
|
493
|
744 |
|
|
745 tree_constant arg1; |
|
|
746 tree_constant arg2; |
|
|
747 tree_constant_rep::constant_type arg1_type = |
|
|
748 tree_constant_rep::unknown_constant; |
|
|
749 tree_constant_rep::constant_type arg2_type = |
|
|
750 tree_constant_rep::unknown_constant; |
|
|
751 |
|
506
|
752 int nargin = args.length (); |
|
|
753 |
|
493
|
754 switch (nargin) |
|
|
755 { |
|
|
756 case 3: |
|
500
|
757 arg2 = args(2).make_numeric (); |
|
493
|
758 arg2_type = arg2.const_type (); |
|
|
759 // Fall through... |
|
|
760 case 2: |
|
500
|
761 arg1 = args(1).make_numeric (); |
|
493
|
762 arg1_type = arg1.const_type (); |
|
455
|
763 break; |
|
|
764 default: |
|
|
765 panic_impossible (); |
|
|
766 break; |
|
|
767 } |
|
|
768 |
|
500
|
769 if (nargin == 2 && (nargout == 1 || nargout == 0)) |
|
455
|
770 { |
|
500
|
771 retval.resize (1); |
|
493
|
772 switch (arg1_type) |
|
|
773 { |
|
|
774 case tree_constant_rep::scalar_constant: |
|
500
|
775 retval(0) = tree_constant (arg1.double_value ()); |
|
493
|
776 break; |
|
|
777 case tree_constant_rep::complex_scalar_constant: |
|
500
|
778 retval(0) = tree_constant (arg1.complex_value ()); |
|
493
|
779 break; |
|
|
780 case tree_constant_rep::matrix_constant: |
|
|
781 { |
|
|
782 Matrix m = arg1.matrix_value (); |
|
|
783 if (m.rows () == 1) |
|
500
|
784 retval(0) = tree_constant (m.row_max ()); |
|
493
|
785 else |
|
500
|
786 retval(0) = tree_constant (m.column_max (), 0); |
|
493
|
787 } |
|
|
788 break; |
|
|
789 case tree_constant_rep::complex_matrix_constant: |
|
|
790 { |
|
|
791 ComplexMatrix m = arg1.complex_matrix_value (); |
|
|
792 if (m.rows () == 1) |
|
500
|
793 retval(0) = tree_constant (m.row_max ()); |
|
493
|
794 else |
|
500
|
795 retval(0) = tree_constant (m.column_max (), 0); |
|
493
|
796 } |
|
|
797 break; |
|
|
798 default: |
|
|
799 panic_impossible (); |
|
|
800 break; |
|
|
801 } |
|
|
802 } |
|
|
803 else if (nargin == 2 && nargout == 2) |
|
|
804 { |
|
500
|
805 retval.resize (2); |
|
493
|
806 switch (arg1_type) |
|
|
807 { |
|
|
808 case tree_constant_rep::scalar_constant: |
|
455
|
809 { |
|
500
|
810 retval(0) = tree_constant (arg1.double_value ()); |
|
|
811 retval(1) = tree_constant (1); |
|
455
|
812 } |
|
493
|
813 break; |
|
|
814 case tree_constant_rep::complex_scalar_constant: |
|
455
|
815 { |
|
500
|
816 retval(0) = tree_constant (arg1.complex_value ()); |
|
|
817 retval(1) = tree_constant (1); |
|
455
|
818 } |
|
493
|
819 break; |
|
|
820 case tree_constant_rep::matrix_constant: |
|
455
|
821 { |
|
493
|
822 Matrix m = arg1.matrix_value (); |
|
|
823 if (m.rows () == 1) |
|
455
|
824 { |
|
500
|
825 retval(0) = tree_constant (m.row_max ()); |
|
|
826 retval(1) = tree_constant (m.row_max_loc ()); |
|
493
|
827 } |
|
|
828 else |
|
|
829 { |
|
500
|
830 retval(0) = tree_constant (m.column_max (), 0); |
|
|
831 retval(1) = tree_constant (m.column_max_loc (), 0); |
|
455
|
832 } |
|
|
833 } |
|
493
|
834 break; |
|
|
835 case tree_constant_rep::complex_matrix_constant: |
|
455
|
836 { |
|
493
|
837 ComplexMatrix m = arg1.complex_matrix_value (); |
|
|
838 if (m.rows () == 1) |
|
455
|
839 { |
|
500
|
840 retval(0) = tree_constant (m.row_max ()); |
|
|
841 retval(1) = tree_constant (m.row_max_loc ()); |
|
455
|
842 } |
|
|
843 else |
|
|
844 { |
|
500
|
845 retval(0) = tree_constant (m.column_max (), 0); |
|
|
846 retval(1) = tree_constant (m.column_max_loc (), 0); |
|
455
|
847 } |
|
|
848 } |
|
493
|
849 break; |
|
|
850 default: |
|
|
851 panic_impossible (); |
|
|
852 break; |
|
|
853 } |
|
|
854 } |
|
|
855 else if (nargin == 3) |
|
|
856 { |
|
|
857 if (arg1.rows () == arg2.rows () |
|
|
858 && arg1.columns () == arg2.columns ()) |
|
|
859 { |
|
500
|
860 retval.resize (1); |
|
493
|
861 switch (arg1_type) |
|
|
862 { |
|
|
863 case tree_constant_rep::scalar_constant: |
|
455
|
864 { |
|
493
|
865 double result; |
|
|
866 double a_elem = arg1.double_value (); |
|
|
867 double b_elem = arg2.double_value (); |
|
|
868 result = MAX (a_elem, b_elem); |
|
500
|
869 retval(0) = tree_constant (result); |
|
455
|
870 } |
|
493
|
871 break; |
|
|
872 case tree_constant_rep::complex_scalar_constant: |
|
455
|
873 { |
|
493
|
874 Complex result; |
|
|
875 Complex a_elem = arg1.complex_value (); |
|
|
876 Complex b_elem = arg2.complex_value (); |
|
|
877 if (abs (a_elem) > abs (b_elem)) |
|
|
878 result = a_elem; |
|
|
879 else |
|
|
880 result = b_elem; |
|
500
|
881 retval(0) = tree_constant (result); |
|
493
|
882 } |
|
|
883 break; |
|
|
884 case tree_constant_rep::matrix_constant: |
|
|
885 { |
|
|
886 Matrix result; |
|
|
887 result = max (arg1.matrix_value (), arg2.matrix_value ()); |
|
500
|
888 retval(0) = tree_constant (result); |
|
455
|
889 } |
|
493
|
890 break; |
|
|
891 case tree_constant_rep::complex_matrix_constant: |
|
455
|
892 { |
|
493
|
893 ComplexMatrix result; |
|
|
894 result = max (arg1.complex_matrix_value (), |
|
|
895 arg2.complex_matrix_value ()); |
|
500
|
896 retval(0) = tree_constant (result); |
|
455
|
897 } |
|
493
|
898 break; |
|
|
899 default: |
|
|
900 panic_impossible (); |
|
|
901 break; |
|
|
902 } |
|
|
903 } |
|
|
904 else |
|
|
905 error ("max: nonconformant matrices"); |
|
|
906 } |
|
|
907 else |
|
|
908 panic_impossible (); |
|
|
909 |
|
|
910 return retval; |
|
|
911 } |
|
455
|
912 |
|
500
|
913 Octave_object |
|
506
|
914 column_min (const Octave_object& args, int nargout) |
|
493
|
915 { |
|
500
|
916 Octave_object retval; |
|
493
|
917 |
|
|
918 tree_constant arg1; |
|
|
919 tree_constant arg2; |
|
|
920 tree_constant_rep::constant_type arg1_type = |
|
|
921 tree_constant_rep::unknown_constant; |
|
|
922 tree_constant_rep::constant_type arg2_type = |
|
|
923 tree_constant_rep::unknown_constant; |
|
|
924 |
|
506
|
925 int nargin = args.length (); |
|
|
926 |
|
493
|
927 switch (nargin) |
|
|
928 { |
|
|
929 case 3: |
|
500
|
930 arg2 = args(2).make_numeric (); |
|
493
|
931 arg2_type = arg2.const_type (); |
|
|
932 // Fall through... |
|
|
933 case 2: |
|
500
|
934 arg1 = args(1).make_numeric (); |
|
493
|
935 arg1_type = arg1.const_type (); |
|
455
|
936 break; |
|
|
937 default: |
|
|
938 panic_impossible (); |
|
|
939 break; |
|
|
940 } |
|
|
941 |
|
500
|
942 if (nargin == 2 && (nargout == 1 || nargout == 0)) |
|
455
|
943 { |
|
500
|
944 retval.resize (1); |
|
493
|
945 switch (arg1_type) |
|
455
|
946 { |
|
493
|
947 case tree_constant_rep::scalar_constant: |
|
500
|
948 retval(0) = tree_constant (arg1.double_value ()); |
|
493
|
949 break; |
|
|
950 case tree_constant_rep::complex_scalar_constant: |
|
500
|
951 retval(0) = tree_constant (arg1.complex_value ()); |
|
493
|
952 break; |
|
|
953 case tree_constant_rep::matrix_constant: |
|
|
954 { |
|
|
955 Matrix m = arg1.matrix_value (); |
|
|
956 if (m.rows () == 1) |
|
500
|
957 retval(0) = tree_constant (m.row_min ()); |
|
493
|
958 else |
|
500
|
959 retval(0) = tree_constant (m.column_min (), 0); |
|
493
|
960 } |
|
|
961 break; |
|
|
962 case tree_constant_rep::complex_matrix_constant: |
|
|
963 { |
|
|
964 ComplexMatrix m = arg1.complex_matrix_value (); |
|
|
965 if (m.rows () == 1) |
|
500
|
966 retval(0) = tree_constant (m.row_min ()); |
|
493
|
967 else |
|
500
|
968 retval(0) = tree_constant (m.column_min (), 0); |
|
493
|
969 } |
|
|
970 break; |
|
|
971 default: |
|
|
972 panic_impossible (); |
|
|
973 break; |
|
455
|
974 } |
|
|
975 } |
|
493
|
976 else if (nargin == 2 && nargout == 2) |
|
455
|
977 { |
|
500
|
978 retval.resize (2); |
|
493
|
979 switch (arg1_type) |
|
|
980 { |
|
|
981 case tree_constant_rep::scalar_constant: |
|
|
982 { |
|
500
|
983 retval(0) = tree_constant (arg1.double_value ()); |
|
|
984 retval(1) = tree_constant (1); |
|
493
|
985 } |
|
|
986 break; |
|
|
987 case tree_constant_rep::complex_scalar_constant: |
|
|
988 { |
|
500
|
989 retval(0) = tree_constant (arg1.complex_value ()); |
|
|
990 retval(1) = tree_constant (1); |
|
493
|
991 } |
|
|
992 break; |
|
|
993 case tree_constant_rep::matrix_constant: |
|
|
994 { |
|
|
995 Matrix m = arg1.matrix_value (); |
|
|
996 if (m.rows () == 1) |
|
|
997 { |
|
500
|
998 retval(0) = tree_constant (m.row_min ()); |
|
|
999 retval(1) = tree_constant (m.row_min_loc ()); |
|
493
|
1000 } |
|
|
1001 else |
|
|
1002 { |
|
500
|
1003 retval(0) = tree_constant (m.column_min (), 0); |
|
|
1004 retval(1) = tree_constant (m.column_min_loc (), 0); |
|
493
|
1005 } |
|
|
1006 } |
|
|
1007 break; |
|
|
1008 case tree_constant_rep::complex_matrix_constant: |
|
|
1009 { |
|
|
1010 ComplexMatrix m = arg1.complex_matrix_value (); |
|
|
1011 if (m.rows () == 1) |
|
|
1012 { |
|
500
|
1013 retval(0) = tree_constant (m.row_min ()); |
|
|
1014 retval(1) = tree_constant (m.row_min_loc ()); |
|
493
|
1015 } |
|
|
1016 else |
|
|
1017 { |
|
500
|
1018 retval(0) = tree_constant (m.column_min (), 0); |
|
|
1019 retval(1) = tree_constant (m.column_min_loc (), 0); |
|
493
|
1020 } |
|
|
1021 } |
|
|
1022 break; |
|
|
1023 default: |
|
|
1024 panic_impossible (); |
|
|
1025 break; |
|
|
1026 } |
|
455
|
1027 } |
|
493
|
1028 else if (nargin == 3) |
|
455
|
1029 { |
|
493
|
1030 if (arg1.rows () == arg2.rows () |
|
|
1031 && arg1.columns () == arg2.columns ()) |
|
455
|
1032 { |
|
500
|
1033 retval.resize (1); |
|
493
|
1034 switch (arg1_type) |
|
|
1035 { |
|
|
1036 case tree_constant_rep::scalar_constant: |
|
|
1037 { |
|
|
1038 double result; |
|
|
1039 double a_elem = arg1.double_value (); |
|
|
1040 double b_elem = arg2.double_value (); |
|
|
1041 result = MIN (a_elem, b_elem); |
|
500
|
1042 retval(0) = tree_constant (result); |
|
493
|
1043 } |
|
|
1044 break; |
|
|
1045 case tree_constant_rep::complex_scalar_constant: |
|
|
1046 { |
|
|
1047 Complex result; |
|
|
1048 Complex a_elem = arg1.complex_value (); |
|
|
1049 Complex b_elem = arg2.complex_value (); |
|
|
1050 if (abs (a_elem) < abs (b_elem)) |
|
|
1051 result = a_elem; |
|
|
1052 else |
|
|
1053 result = b_elem; |
|
500
|
1054 retval(0) = tree_constant (result); |
|
493
|
1055 } |
|
|
1056 break; |
|
|
1057 case tree_constant_rep::matrix_constant: |
|
|
1058 { |
|
|
1059 Matrix result; |
|
|
1060 result = min (arg1.matrix_value (), arg2.matrix_value ()); |
|
500
|
1061 retval(0) = tree_constant (result); |
|
493
|
1062 } |
|
|
1063 break; |
|
|
1064 case tree_constant_rep::complex_matrix_constant: |
|
|
1065 { |
|
|
1066 ComplexMatrix result; |
|
|
1067 result = min (arg1.complex_matrix_value (), |
|
|
1068 arg2.complex_matrix_value ()); |
|
500
|
1069 retval(0) = tree_constant (result); |
|
493
|
1070 } |
|
|
1071 break; |
|
|
1072 default: |
|
|
1073 panic_impossible (); |
|
|
1074 break; |
|
455
|
1075 } |
|
|
1076 } |
|
493
|
1077 else |
|
|
1078 error ("min: nonconformant matrices"); |
|
455
|
1079 } |
|
|
1080 else |
|
|
1081 panic_impossible (); |
|
|
1082 |
|
|
1083 return retval; |
|
|
1084 } |
|
|
1085 |
|
493
|
1086 static void |
|
|
1087 mx_sort (Matrix& m, Matrix& idx, int return_idx) |
|
455
|
1088 { |
|
493
|
1089 int nr = m.rows (); |
|
|
1090 int nc = m.columns (); |
|
|
1091 idx.resize (nr, nc); |
|
|
1092 int i, j; |
|
|
1093 |
|
|
1094 if (return_idx) |
|
|
1095 { |
|
|
1096 for (j = 0; j < nc; j++) |
|
|
1097 for (i = 0; i < nr; i++) |
|
|
1098 idx.elem (i, j) = i+1; |
|
|
1099 } |
|
|
1100 |
|
|
1101 for (j = 0; j < nc; j++) |
|
|
1102 { |
|
|
1103 for (int gap = nr/2; gap > 0; gap /= 2) |
|
|
1104 for (i = gap; i < nr; i++) |
|
|
1105 for (int k = i - gap; |
|
|
1106 k >= 0 && m.elem (k, j) > m.elem (k+gap, j); |
|
|
1107 k -= gap) |
|
|
1108 { |
|
|
1109 double tmp = m.elem (k, j); |
|
|
1110 m.elem (k, j) = m.elem (k+gap, j); |
|
|
1111 m.elem (k+gap, j) = tmp; |
|
|
1112 |
|
|
1113 if (return_idx) |
|
|
1114 { |
|
|
1115 double tmp = idx.elem (k, j); |
|
|
1116 idx.elem (k, j) = idx.elem (k+gap, j); |
|
|
1117 idx.elem (k+gap, j) = tmp; |
|
|
1118 } |
|
|
1119 } |
|
|
1120 } |
|
|
1121 } |
|
455
|
1122 |
|
493
|
1123 static void |
|
|
1124 mx_sort (RowVector& v, RowVector& idx, int return_idx) |
|
|
1125 { |
|
|
1126 int n = v.capacity (); |
|
|
1127 idx.resize (n); |
|
|
1128 int i; |
|
|
1129 |
|
|
1130 if (return_idx) |
|
|
1131 for (i = 0; i < n; i++) |
|
|
1132 idx.elem (i) = i+1; |
|
455
|
1133 |
|
493
|
1134 for (int gap = n/2; gap > 0; gap /= 2) |
|
|
1135 for (i = gap; i < n; i++) |
|
|
1136 for (int k = i - gap; |
|
|
1137 k >= 0 && v.elem (k) > v.elem (k+gap); |
|
|
1138 k -= gap) |
|
|
1139 { |
|
|
1140 double tmp = v.elem (k); |
|
|
1141 v.elem (k) = v.elem (k+gap); |
|
|
1142 v.elem (k+gap) = tmp; |
|
455
|
1143 |
|
493
|
1144 if (return_idx) |
|
|
1145 { |
|
|
1146 double tmp = idx.elem (k); |
|
|
1147 idx.elem (k) = idx.elem (k+gap); |
|
|
1148 idx.elem (k+gap) = tmp; |
|
|
1149 } |
|
|
1150 } |
|
|
1151 } |
|
|
1152 |
|
|
1153 static void |
|
|
1154 mx_sort (ComplexMatrix& cm, Matrix& idx, int return_idx) |
|
|
1155 { |
|
|
1156 int nr = cm.rows (); |
|
|
1157 int nc = cm.columns (); |
|
|
1158 idx.resize (nr, nc); |
|
|
1159 int i, j; |
|
455
|
1160 |
|
493
|
1161 if (return_idx) |
|
|
1162 { |
|
|
1163 for (j = 0; j < nc; j++) |
|
|
1164 for (i = 0; i < nr; i++) |
|
|
1165 idx.elem (i, j) = i+1; |
|
|
1166 } |
|
|
1167 |
|
|
1168 for (j = 0; j < nc; j++) |
|
455
|
1169 { |
|
493
|
1170 for (int gap = nr/2; gap > 0; gap /= 2) |
|
|
1171 for (i = gap; i < nr; i++) |
|
|
1172 for (int k = i - gap; |
|
|
1173 k >= 0 && abs (cm.elem (k, j)) > abs (cm.elem (k+gap, j)); |
|
|
1174 k -= gap) |
|
|
1175 { |
|
|
1176 Complex ctmp = cm.elem (k, j); |
|
|
1177 cm.elem (k, j) = cm.elem (k+gap, j); |
|
|
1178 cm.elem (k+gap, j) = ctmp; |
|
|
1179 |
|
|
1180 if (return_idx) |
|
|
1181 { |
|
|
1182 double tmp = idx.elem (k, j); |
|
|
1183 idx.elem (k, j) = idx.elem (k+gap, j); |
|
|
1184 idx.elem (k+gap, j) = tmp; |
|
|
1185 } |
|
|
1186 } |
|
|
1187 } |
|
|
1188 } |
|
|
1189 |
|
|
1190 static void |
|
|
1191 mx_sort (ComplexRowVector& cv, RowVector& idx, int return_idx) |
|
|
1192 { |
|
|
1193 int n = cv.capacity (); |
|
|
1194 idx.resize (n); |
|
|
1195 int i; |
|
|
1196 |
|
|
1197 if (return_idx) |
|
|
1198 for (i = 0; i < n; i++) |
|
|
1199 idx.elem (i) = i+1; |
|
|
1200 |
|
|
1201 for (int gap = n/2; gap > 0; gap /= 2) |
|
|
1202 for (i = gap; i < n; i++) |
|
|
1203 for (int k = i - gap; |
|
|
1204 k >= 0 && abs (cv.elem (k)) > abs (cv.elem (k+gap)); |
|
|
1205 k -= gap) |
|
|
1206 { |
|
|
1207 Complex tmp = cv.elem (k); |
|
|
1208 cv.elem (k) = cv.elem (k+gap); |
|
|
1209 cv.elem (k+gap) = tmp; |
|
|
1210 |
|
|
1211 if (return_idx) |
|
|
1212 { |
|
|
1213 double tmp = idx.elem (k); |
|
|
1214 idx.elem (k) = idx.elem (k+gap); |
|
|
1215 idx.elem (k+gap) = tmp; |
|
|
1216 } |
|
|
1217 } |
|
|
1218 } |
|
|
1219 |
|
500
|
1220 Octave_object |
|
506
|
1221 sort (const Octave_object& args, int nargout) |
|
493
|
1222 { |
|
|
1223 // Assumes that we have been given the correct number of arguments. |
|
|
1224 |
|
500
|
1225 Octave_object retval; |
|
493
|
1226 |
|
|
1227 int return_idx = nargout > 1; |
|
|
1228 if (return_idx) |
|
500
|
1229 retval.resize (2); |
|
493
|
1230 else |
|
500
|
1231 retval.resize (1); |
|
493
|
1232 |
|
500
|
1233 switch (args(1).const_type ()) |
|
493
|
1234 { |
|
|
1235 case tree_constant_rep::scalar_constant: |
|
455
|
1236 { |
|
500
|
1237 retval(0) = tree_constant (args(1).double_value ()); |
|
493
|
1238 if (return_idx) |
|
500
|
1239 retval(1) = tree_constant (1.0); |
|
455
|
1240 } |
|
|
1241 break; |
|
493
|
1242 case tree_constant_rep::complex_scalar_constant: |
|
|
1243 { |
|
500
|
1244 retval(0) = tree_constant (args(1).complex_value ()); |
|
493
|
1245 if (return_idx) |
|
500
|
1246 retval(1) = tree_constant (1.0); |
|
493
|
1247 } |
|
|
1248 break; |
|
|
1249 case tree_constant_rep::string_constant: |
|
|
1250 case tree_constant_rep::range_constant: |
|
|
1251 case tree_constant_rep::matrix_constant: |
|
455
|
1252 { |
|
500
|
1253 Matrix m = args(1).to_matrix (); |
|
493
|
1254 if (m.rows () == 1) |
|
455
|
1255 { |
|
493
|
1256 int nc = m.columns (); |
|
|
1257 RowVector v (nc); |
|
|
1258 for (int i = 0; i < nc; i++) |
|
|
1259 v.elem (i) = m.elem (0, i); |
|
|
1260 RowVector idx; |
|
|
1261 mx_sort (v, idx, return_idx); |
|
|
1262 |
|
500
|
1263 retval(0) = tree_constant (v, 0); |
|
493
|
1264 if (return_idx) |
|
500
|
1265 retval(1) = tree_constant (idx, 0); |
|
455
|
1266 } |
|
|
1267 else |
|
|
1268 { |
|
493
|
1269 // Sorts m in place, optionally computes index Matrix. |
|
|
1270 Matrix idx; |
|
|
1271 mx_sort (m, idx, return_idx); |
|
|
1272 |
|
500
|
1273 retval(0) = tree_constant (m); |
|
493
|
1274 if (return_idx) |
|
500
|
1275 retval(1) = tree_constant (idx); |
|
455
|
1276 } |
|
|
1277 } |
|
|
1278 break; |
|
493
|
1279 case tree_constant_rep::complex_matrix_constant: |
|
|
1280 { |
|
500
|
1281 ComplexMatrix cm = args(1).complex_matrix_value (); |
|
493
|
1282 if (cm.rows () == 1) |
|
|
1283 { |
|
|
1284 int nc = cm.columns (); |
|
|
1285 ComplexRowVector cv (nc); |
|
|
1286 for (int i = 0; i < nc; i++) |
|
|
1287 cv.elem (i) = cm.elem (0, i); |
|
|
1288 RowVector idx; |
|
|
1289 mx_sort (cv, idx, return_idx); |
|
|
1290 |
|
500
|
1291 retval(0) = tree_constant (cv, 0); |
|
493
|
1292 if (return_idx) |
|
500
|
1293 retval(1) = tree_constant (idx, 0); |
|
493
|
1294 } |
|
|
1295 else |
|
|
1296 { |
|
|
1297 // Sorts cm in place, optionally computes index Matrix. |
|
|
1298 Matrix idx; |
|
|
1299 mx_sort (cm, idx, return_idx); |
|
|
1300 |
|
500
|
1301 retval(0) = tree_constant (cm); |
|
493
|
1302 if (return_idx) |
|
500
|
1303 retval(1) = tree_constant (idx); |
|
493
|
1304 } |
|
|
1305 } |
|
455
|
1306 break; |
|
|
1307 default: |
|
|
1308 panic_impossible (); |
|
|
1309 break; |
|
|
1310 } |
|
|
1311 |
|
|
1312 return retval; |
|
|
1313 } |
|
|
1314 |
|
500
|
1315 Octave_object |
|
506
|
1316 feval (const Octave_object& args, int nargout) |
|
455
|
1317 { |
|
493
|
1318 // Assumes that we have been given the correct number of arguments. |
|
455
|
1319 |
|
500
|
1320 Octave_object retval; |
|
455
|
1321 |
|
500
|
1322 tree_fvc *fcn = is_valid_function (args(1), "feval", 1); |
|
493
|
1323 if (fcn != (tree_fvc *) NULL) |
|
455
|
1324 { |
|
506
|
1325 int nargin = args.length () - 1; |
|
|
1326 Octave_object tmp_args (nargin); |
|
500
|
1327 for (int i = 0; i < nargin; i++) |
|
|
1328 tmp_args(i) = args(i+1); |
|
506
|
1329 retval = fcn->eval (0, nargout, tmp_args); |
|
455
|
1330 } |
|
|
1331 |
|
|
1332 return retval; |
|
|
1333 } |
|
|
1334 |
|
|
1335 tree_constant |
|
493
|
1336 eval_string (const char *string, int print, int ans_assign, |
|
|
1337 int& parse_status) |
|
455
|
1338 { |
|
493
|
1339 begin_unwind_frame ("eval_string"); |
|
455
|
1340 |
|
493
|
1341 unwind_protect_int (get_input_from_eval_string); |
|
|
1342 unwind_protect_ptr (global_command); |
|
|
1343 unwind_protect_ptr (current_eval_string); |
|
|
1344 |
|
|
1345 get_input_from_eval_string = 1; |
|
|
1346 current_eval_string = string; |
|
|
1347 |
|
|
1348 YY_BUFFER_STATE old_buf = current_buffer (); |
|
|
1349 YY_BUFFER_STATE new_buf = create_buffer ((FILE *) NULL); |
|
455
|
1350 |
|
493
|
1351 add_unwind_protect (restore_input_buffer, (void *) old_buf); |
|
|
1352 add_unwind_protect (delete_input_buffer, (void *) new_buf); |
|
|
1353 |
|
|
1354 switch_to_buffer (new_buf); |
|
|
1355 |
|
|
1356 unwind_protect_ptr (curr_sym_tab); |
|
|
1357 |
|
|
1358 reset_parser (); |
|
455
|
1359 |
|
493
|
1360 parse_status = yyparse (); |
|
455
|
1361 |
|
493
|
1362 // Important to reset the idea of where input is coming from before |
|
|
1363 // trying to eval the command we just parsed -- it might contain the |
|
|
1364 // name of an function file that still needs to be parsed! |
|
|
1365 |
|
|
1366 tree *command = global_command; |
|
|
1367 |
|
|
1368 run_unwind_frame ("eval_string"); |
|
|
1369 |
|
455
|
1370 tree_constant retval; |
|
|
1371 |
|
493
|
1372 if (parse_status == 0 && command != NULL_TREE) |
|
455
|
1373 { |
|
493
|
1374 retval = command->eval (print); |
|
|
1375 delete command; |
|
455
|
1376 } |
|
|
1377 |
|
|
1378 return retval; |
|
|
1379 } |
|
|
1380 |
|
|
1381 tree_constant |
|
493
|
1382 eval_string (const tree_constant& arg, int& parse_status) |
|
455
|
1383 { |
|
493
|
1384 if (! arg.is_string_type ()) |
|
455
|
1385 { |
|
493
|
1386 error ("eval: expecting string argument"); |
|
|
1387 return -1; |
|
455
|
1388 } |
|
|
1389 |
|
493
|
1390 char *string = arg.string_value (); |
|
|
1391 |
|
|
1392 // Yes Virginia, we always print here... |
|
|
1393 |
|
|
1394 return eval_string (string, 1, 1, parse_status); |
|
|
1395 } |
|
|
1396 |
|
|
1397 static int |
|
|
1398 match_sans_spaces (const char *standard, const char *test) |
|
|
1399 { |
|
|
1400 const char *tp = test; |
|
|
1401 while (*tp == ' ' || *tp == '\t') |
|
|
1402 tp++; |
|
|
1403 |
|
|
1404 const char *ep = test + strlen (test) - 1; |
|
|
1405 while (*ep == ' ' || *ep == '\t') |
|
|
1406 ep--; |
|
|
1407 |
|
|
1408 int len = ep - tp + 1; |
|
|
1409 |
|
|
1410 return (strncmp (standard, tp, len) == 0); |
|
455
|
1411 } |
|
|
1412 |
|
|
1413 tree_constant |
|
506
|
1414 get_user_input (const Octave_object& args, int nargout, int debug = 0) |
|
455
|
1415 { |
|
|
1416 tree_constant retval; |
|
|
1417 |
|
506
|
1418 int nargin = args.length (); |
|
|
1419 |
|
493
|
1420 int read_as_string = 0; |
|
506
|
1421 |
|
493
|
1422 if (nargin == 3) |
|
455
|
1423 { |
|
500
|
1424 if (args(2).is_string_type () |
|
|
1425 && strcmp ("s", args(2).string_value ()) == 0) |
|
493
|
1426 read_as_string++; |
|
|
1427 else |
|
455
|
1428 { |
|
493
|
1429 error ("input: unrecognized second argument"); |
|
|
1430 return retval; |
|
455
|
1431 } |
|
|
1432 } |
|
|
1433 |
|
493
|
1434 char *prompt = "debug> "; |
|
|
1435 if (nargin > 1) |
|
|
1436 { |
|
500
|
1437 if (args(1).is_string_type ()) |
|
|
1438 prompt = args(1).string_value (); |
|
493
|
1439 else |
|
455
|
1440 { |
|
493
|
1441 error ("input: unrecognized argument"); |
|
|
1442 return retval; |
|
455
|
1443 } |
|
|
1444 } |
|
|
1445 |
|
493
|
1446 again: |
|
455
|
1447 |
|
493
|
1448 flush_output_to_pager (); |
|
|
1449 |
|
|
1450 char *input_buf = gnu_readline (prompt); |
|
455
|
1451 |
|
493
|
1452 if (input_buf != (char *) NULL) |
|
|
1453 { |
|
|
1454 if (input_buf) |
|
|
1455 maybe_save_history (input_buf); |
|
455
|
1456 |
|
493
|
1457 int len = strlen (input_buf); |
|
455
|
1458 |
|
493
|
1459 if (len < 1) |
|
|
1460 { |
|
|
1461 if (debug) |
|
|
1462 goto again; |
|
|
1463 else |
|
|
1464 return retval; |
|
|
1465 } |
|
455
|
1466 |
|
493
|
1467 if (match_sans_spaces ("exit", input_buf) |
|
|
1468 || match_sans_spaces ("quit", input_buf) |
|
|
1469 || match_sans_spaces ("return", input_buf)) |
|
|
1470 return tree_constant (); |
|
|
1471 else if (read_as_string) |
|
|
1472 retval = tree_constant (input_buf); |
|
|
1473 else |
|
455
|
1474 { |
|
493
|
1475 int parse_status; |
|
|
1476 retval = eval_string (input_buf, 0, 0, parse_status); |
|
|
1477 if (debug && retval.is_defined ()) |
|
|
1478 retval.eval (1); |
|
455
|
1479 } |
|
|
1480 } |
|
493
|
1481 else |
|
|
1482 error ("input: reading user-input failed!"); |
|
455
|
1483 |
|
493
|
1484 if (debug) |
|
|
1485 goto again; |
|
455
|
1486 |
|
|
1487 return retval; |
|
|
1488 } |
|
|
1489 |
|
96
|
1490 /* |
|
1
|
1491 ;;; Local Variables: *** |
|
|
1492 ;;; mode: C++ *** |
|
|
1493 ;;; page-delimiter: "^/\\*" *** |
|
|
1494 ;;; End: *** |
|
|
1495 */ |
