Mercurial > octave-nkf
comparison src/ov-flt-cx-mat.cc @ 7789:82be108cc558
First attempt at single precision tyeps
* * *
corrections to qrupdate single precision routines
* * *
prefer demotion to single over promotion to double
* * *
Add single precision support to log2 function
* * *
Trivial PROJECT file update
* * *
Cache optimized hermitian/transpose methods
* * *
Add tests for tranpose/hermitian and ChangeLog entry for new transpose code
author | David Bateman <dbateman@free.fr> |
---|---|
date | Sun, 27 Apr 2008 22:34:17 +0200 |
parents | |
children | a41df65f3f00 |
comparison
equal
deleted
inserted
replaced
7788:45f5faba05a2 | 7789:82be108cc558 |
---|---|
1 /* | |
2 | |
3 Copyright (C) 1996, 1997, 1998, 2000, 2001, 2002, 2003, 2004, 2005, | |
4 2006, 2007 John W. Eaton | |
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 3 of the License, or (at your | |
11 option) any 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, see | |
20 <http://www.gnu.org/licenses/>. | |
21 | |
22 */ | |
23 | |
24 #ifdef HAVE_CONFIG_H | |
25 #include <config.h> | |
26 #endif | |
27 | |
28 #include <iostream> | |
29 #include <vector> | |
30 | |
31 #include "data-conv.h" | |
32 #include "lo-ieee.h" | |
33 #include "lo-specfun.h" | |
34 #include "lo-mappers.h" | |
35 #include "mx-base.h" | |
36 #include "mach-info.h" | |
37 | |
38 #include "gripes.h" | |
39 #include "oct-obj.h" | |
40 #include "oct-stream.h" | |
41 #include "ops.h" | |
42 #include "ov-base.h" | |
43 #include "ov-base-mat.h" | |
44 #include "ov-base-mat.cc" | |
45 #include "ov-complex.h" | |
46 #include "ov-flt-complex.h" | |
47 #include "ov-cx-mat.h" | |
48 #include "ov-flt-cx-mat.h" | |
49 #include "ov-re-mat.h" | |
50 #include "ov-flt-re-mat.h" | |
51 #include "ov-scalar.h" | |
52 #include "ov-float.h" | |
53 #include "pr-output.h" | |
54 #include "ops.h" | |
55 | |
56 #include "byte-swap.h" | |
57 #include "ls-oct-ascii.h" | |
58 #include "ls-hdf5.h" | |
59 #include "ls-utils.h" | |
60 | |
61 template class octave_base_matrix<FloatComplexNDArray>; | |
62 | |
63 DEFINE_OCTAVE_ALLOCATOR (octave_float_complex_matrix); | |
64 | |
65 DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA (octave_float_complex_matrix, | |
66 "float complex matrix", "single"); | |
67 | |
68 octave_base_value * | |
69 octave_float_complex_matrix::try_narrowing_conversion (void) | |
70 { | |
71 octave_base_value *retval = 0; | |
72 | |
73 if (matrix.ndims () == 2) | |
74 { | |
75 FloatComplexMatrix cm = matrix.matrix_value (); | |
76 | |
77 octave_idx_type nr = cm.rows (); | |
78 octave_idx_type nc = cm.cols (); | |
79 | |
80 if (nr == 1 && nc == 1) | |
81 { | |
82 FloatComplex c = matrix (0, 0); | |
83 | |
84 float im = std::imag (c); | |
85 | |
86 if (im == 0.0 && ! lo_ieee_signbit (im)) | |
87 retval = new octave_float_scalar (std::real (c)); | |
88 else | |
89 retval = new octave_float_complex (c); | |
90 } | |
91 else if (nr == 0 || nc == 0) | |
92 retval = new octave_float_matrix (FloatMatrix (nr, nc)); | |
93 else if (cm.all_elements_are_real ()) | |
94 retval = new octave_float_matrix (::real (cm)); | |
95 } | |
96 else if (matrix.all_elements_are_real ()) | |
97 retval = new octave_float_matrix (::real (matrix)); | |
98 | |
99 return retval; | |
100 } | |
101 | |
102 void | |
103 octave_float_complex_matrix::assign (const octave_value_list& idx, | |
104 const FloatComplexNDArray& rhs) | |
105 { | |
106 octave_base_matrix<FloatComplexNDArray>::assign (idx, rhs); | |
107 } | |
108 | |
109 void | |
110 octave_float_complex_matrix::assign (const octave_value_list& idx, | |
111 const FloatNDArray& rhs) | |
112 { | |
113 octave_idx_type len = idx.length (); | |
114 | |
115 for (octave_idx_type i = 0; i < len; i++) | |
116 matrix.set_index (idx(i).index_vector ()); | |
117 | |
118 ::assign (matrix, rhs); | |
119 } | |
120 | |
121 bool | |
122 octave_float_complex_matrix::valid_as_scalar_index (void) const | |
123 { | |
124 // FIXME | |
125 return false; | |
126 } | |
127 | |
128 double | |
129 octave_float_complex_matrix::double_value (bool force_conversion) const | |
130 { | |
131 double retval = lo_ieee_nan_value (); | |
132 | |
133 if (! force_conversion) | |
134 gripe_implicit_conversion ("Octave:imag-to-real", | |
135 "complex matrix", "real scalar"); | |
136 | |
137 if (rows () > 0 && columns () > 0) | |
138 { | |
139 gripe_implicit_conversion ("Octave:array-as-scalar", | |
140 "complex matrix", "real scalar"); | |
141 | |
142 retval = std::real (matrix (0, 0)); | |
143 } | |
144 else | |
145 gripe_invalid_conversion ("complex matrix", "real scalar"); | |
146 | |
147 return retval; | |
148 } | |
149 | |
150 float | |
151 octave_float_complex_matrix::float_value (bool force_conversion) const | |
152 { | |
153 float retval = lo_ieee_float_nan_value (); | |
154 | |
155 if (! force_conversion) | |
156 gripe_implicit_conversion ("Octave:imag-to-real", | |
157 "complex matrix", "real scalar"); | |
158 | |
159 if (rows () > 0 && columns () > 0) | |
160 { | |
161 gripe_implicit_conversion ("Octave:array-as-scalar", | |
162 "complex matrix", "real scalar"); | |
163 | |
164 retval = std::real (matrix (0, 0)); | |
165 } | |
166 else | |
167 gripe_invalid_conversion ("complex matrix", "real scalar"); | |
168 | |
169 return retval; | |
170 } | |
171 | |
172 Matrix | |
173 octave_float_complex_matrix::matrix_value (bool force_conversion) const | |
174 { | |
175 Matrix retval; | |
176 | |
177 if (! force_conversion) | |
178 gripe_implicit_conversion ("Octave:imag-to-real", | |
179 "complex matrix", "real matrix"); | |
180 | |
181 retval = ::real (matrix.matrix_value ()); | |
182 | |
183 return retval; | |
184 } | |
185 | |
186 FloatMatrix | |
187 octave_float_complex_matrix::float_matrix_value (bool force_conversion) const | |
188 { | |
189 FloatMatrix retval; | |
190 | |
191 if (! force_conversion) | |
192 gripe_implicit_conversion ("Octave:imag-to-real", | |
193 "complex matrix", "real matrix"); | |
194 | |
195 retval = ::real (matrix.matrix_value ()); | |
196 | |
197 return retval; | |
198 } | |
199 | |
200 Complex | |
201 octave_float_complex_matrix::complex_value (bool) const | |
202 { | |
203 double tmp = lo_ieee_nan_value (); | |
204 | |
205 Complex retval (tmp, tmp); | |
206 | |
207 if (rows () > 0 && columns () > 0) | |
208 { | |
209 gripe_implicit_conversion ("Octave:array-as-scalar", | |
210 "complex matrix", "complex scalar"); | |
211 | |
212 retval = matrix (0, 0); | |
213 } | |
214 else | |
215 gripe_invalid_conversion ("complex matrix", "complex scalar"); | |
216 | |
217 return retval; | |
218 } | |
219 | |
220 FloatComplex | |
221 octave_float_complex_matrix::float_complex_value (bool) const | |
222 { | |
223 float tmp = lo_ieee_float_nan_value (); | |
224 | |
225 FloatComplex retval (tmp, tmp); | |
226 | |
227 if (rows () > 0 && columns () > 0) | |
228 { | |
229 gripe_implicit_conversion ("Octave:array-as-scalar", | |
230 "complex matrix", "complex scalar"); | |
231 | |
232 retval = matrix (0, 0); | |
233 } | |
234 else | |
235 gripe_invalid_conversion ("complex matrix", "complex scalar"); | |
236 | |
237 return retval; | |
238 } | |
239 | |
240 ComplexMatrix | |
241 octave_float_complex_matrix::complex_matrix_value (bool) const | |
242 { | |
243 return matrix.matrix_value (); | |
244 } | |
245 | |
246 FloatComplexMatrix | |
247 octave_float_complex_matrix::float_complex_matrix_value (bool) const | |
248 { | |
249 return FloatComplexMatrix (matrix.matrix_value ()); | |
250 } | |
251 | |
252 charNDArray | |
253 octave_float_complex_matrix::char_array_value (bool frc_str_conv) const | |
254 { | |
255 charNDArray retval; | |
256 | |
257 if (! frc_str_conv) | |
258 gripe_implicit_conversion ("Octave:num-to-str", | |
259 "complex matrix", "string"); | |
260 else | |
261 { | |
262 retval = charNDArray (dims ()); | |
263 octave_idx_type nel = numel (); | |
264 | |
265 for (octave_idx_type i = 0; i < nel; i++) | |
266 retval.elem (i) = static_cast<char>(std::real (matrix.elem (i))); | |
267 } | |
268 | |
269 return retval; | |
270 } | |
271 | |
272 FloatComplexNDArray | |
273 octave_float_complex_matrix::float_complex_array_value (bool) const | |
274 { | |
275 return FloatComplexNDArray (matrix); | |
276 } | |
277 | |
278 SparseMatrix | |
279 octave_float_complex_matrix::sparse_matrix_value (bool force_conversion) const | |
280 { | |
281 SparseMatrix retval; | |
282 | |
283 if (! force_conversion) | |
284 gripe_implicit_conversion ("Octave:imag-to-real", | |
285 "complex matrix", "real matrix"); | |
286 | |
287 retval = SparseMatrix (::real (matrix.matrix_value ())); | |
288 | |
289 return retval; | |
290 } | |
291 | |
292 SparseComplexMatrix | |
293 octave_float_complex_matrix::sparse_complex_matrix_value (bool) const | |
294 { | |
295 return SparseComplexMatrix (matrix.matrix_value ()); | |
296 } | |
297 | |
298 bool | |
299 octave_float_complex_matrix::save_ascii (std::ostream& os) | |
300 { | |
301 dim_vector d = dims (); | |
302 if (d.length () > 2) | |
303 { | |
304 FloatComplexNDArray tmp = complex_array_value (); | |
305 | |
306 os << "# ndims: " << d.length () << "\n"; | |
307 | |
308 for (int i = 0; i < d.length (); i++) | |
309 os << " " << d (i); | |
310 | |
311 os << "\n" << tmp; | |
312 } | |
313 else | |
314 { | |
315 // Keep this case, rather than use generic code above for backward | |
316 // compatiability. Makes load_ascii much more complex!! | |
317 os << "# rows: " << rows () << "\n" | |
318 << "# columns: " << columns () << "\n"; | |
319 | |
320 os << complex_matrix_value (); | |
321 } | |
322 | |
323 return true; | |
324 } | |
325 | |
326 bool | |
327 octave_float_complex_matrix::load_ascii (std::istream& is) | |
328 { | |
329 bool success = true; | |
330 | |
331 string_vector keywords(2); | |
332 | |
333 keywords[0] = "ndims"; | |
334 keywords[1] = "rows"; | |
335 | |
336 std::string kw; | |
337 octave_idx_type val = 0; | |
338 | |
339 if (extract_keyword (is, keywords, kw, val, true)) | |
340 { | |
341 if (kw == "ndims") | |
342 { | |
343 int mdims = static_cast<int> (val); | |
344 | |
345 if (mdims >= 0) | |
346 { | |
347 dim_vector dv; | |
348 dv.resize (mdims); | |
349 | |
350 for (int i = 0; i < mdims; i++) | |
351 is >> dv(i); | |
352 | |
353 if (is) | |
354 { | |
355 FloatComplexNDArray tmp(dv); | |
356 | |
357 if (tmp.is_empty ()) | |
358 matrix = tmp; | |
359 else | |
360 { | |
361 is >> tmp; | |
362 | |
363 if (is) | |
364 matrix = tmp; | |
365 else | |
366 { | |
367 error ("load: failed to load matrix constant"); | |
368 success = false; | |
369 } | |
370 } | |
371 } | |
372 else | |
373 { | |
374 error ("load: failed to read dimensions"); | |
375 success = false; | |
376 } | |
377 } | |
378 else | |
379 { | |
380 error ("load: failed to extract number of dimensions"); | |
381 success = false; | |
382 } | |
383 } | |
384 else if (kw == "rows") | |
385 { | |
386 octave_idx_type nr = val; | |
387 octave_idx_type nc = 0; | |
388 | |
389 if (nr >= 0 && extract_keyword (is, "columns", nc) && nc >= 0) | |
390 { | |
391 if (nr > 0 && nc > 0) | |
392 { | |
393 FloatComplexMatrix tmp (nr, nc); | |
394 is >> tmp; | |
395 if (is) | |
396 matrix = tmp; | |
397 else | |
398 { | |
399 error ("load: failed to load matrix constant"); | |
400 success = false; | |
401 } | |
402 } | |
403 else if (nr == 0 || nc == 0) | |
404 matrix = FloatComplexMatrix (nr, nc); | |
405 else | |
406 panic_impossible (); | |
407 } | |
408 else | |
409 { | |
410 error ("load: failed to extract number of rows and columns"); | |
411 success = false; | |
412 } | |
413 } | |
414 else | |
415 panic_impossible (); | |
416 } | |
417 else | |
418 { | |
419 error ("load: failed to extract number of rows and columns"); | |
420 success = false; | |
421 } | |
422 | |
423 return success; | |
424 } | |
425 | |
426 bool | |
427 octave_float_complex_matrix::save_binary (std::ostream& os, bool&) | |
428 { | |
429 dim_vector d = dims (); | |
430 if (d.length() < 1) | |
431 return false; | |
432 | |
433 // Use negative value for ndims to differentiate with old format!! | |
434 int32_t tmp = - d.length(); | |
435 os.write (reinterpret_cast<char *> (&tmp), 4); | |
436 for (int i = 0; i < d.length (); i++) | |
437 { | |
438 tmp = d(i); | |
439 os.write (reinterpret_cast<char *> (&tmp), 4); | |
440 } | |
441 | |
442 FloatComplexNDArray m = complex_array_value (); | |
443 save_type st = LS_FLOAT; | |
444 if (d.numel () > 4096) // FIXME -- make this configurable. | |
445 { | |
446 float max_val, min_val; | |
447 if (m.all_integers (max_val, min_val)) | |
448 st = get_save_type (max_val, min_val); | |
449 } | |
450 | |
451 const FloatComplex *mtmp = m.data (); | |
452 write_floats (os, reinterpret_cast<const float *> (mtmp), st, 2 * d.numel ()); | |
453 | |
454 return true; | |
455 } | |
456 | |
457 bool | |
458 octave_float_complex_matrix::load_binary (std::istream& is, bool swap, | |
459 oct_mach_info::float_format fmt) | |
460 { | |
461 char tmp; | |
462 int32_t mdims; | |
463 if (! is.read (reinterpret_cast<char *> (&mdims), 4)) | |
464 return false; | |
465 if (swap) | |
466 swap_bytes<4> (&mdims); | |
467 if (mdims < 0) | |
468 { | |
469 mdims = - mdims; | |
470 int32_t di; | |
471 dim_vector dv; | |
472 dv.resize (mdims); | |
473 | |
474 for (int i = 0; i < mdims; i++) | |
475 { | |
476 if (! is.read (reinterpret_cast<char *> (&di), 4)) | |
477 return false; | |
478 if (swap) | |
479 swap_bytes<4> (&di); | |
480 dv(i) = di; | |
481 } | |
482 | |
483 // Convert an array with a single dimension to be a row vector. | |
484 // Octave should never write files like this, other software | |
485 // might. | |
486 | |
487 if (mdims == 1) | |
488 { | |
489 mdims = 2; | |
490 dv.resize (mdims); | |
491 dv(1) = dv(0); | |
492 dv(0) = 1; | |
493 } | |
494 | |
495 if (! is.read (reinterpret_cast<char *> (&tmp), 1)) | |
496 return false; | |
497 | |
498 FloatComplexNDArray m(dv); | |
499 FloatComplex *im = m.fortran_vec (); | |
500 read_floats (is, reinterpret_cast<float *> (im), | |
501 static_cast<save_type> (tmp), 2 * dv.numel (), swap, fmt); | |
502 if (error_state || ! is) | |
503 return false; | |
504 matrix = m; | |
505 } | |
506 else | |
507 { | |
508 int32_t nr, nc; | |
509 nr = mdims; | |
510 if (! is.read (reinterpret_cast<char *> (&nc), 4)) | |
511 return false; | |
512 if (swap) | |
513 swap_bytes<4> (&nc); | |
514 if (! is.read (reinterpret_cast<char *> (&tmp), 1)) | |
515 return false; | |
516 FloatComplexMatrix m (nr, nc); | |
517 FloatComplex *im = m.fortran_vec (); | |
518 octave_idx_type len = nr * nc; | |
519 read_floats (is, reinterpret_cast<float *> (im), | |
520 static_cast<save_type> (tmp), 2*len, swap, fmt); | |
521 if (error_state || ! is) | |
522 return false; | |
523 matrix = m; | |
524 } | |
525 return true; | |
526 } | |
527 | |
528 #if defined (HAVE_HDF5) | |
529 | |
530 bool | |
531 octave_float_complex_matrix::save_hdf5 (hid_t loc_id, const char *name, bool) | |
532 { | |
533 dim_vector dv = dims (); | |
534 int empty = save_hdf5_empty (loc_id, name, dv); | |
535 if (empty) | |
536 return (empty > 0); | |
537 | |
538 int rank = dv.length (); | |
539 hid_t space_hid = -1, data_hid = -1, type_hid = -1; | |
540 bool retval = true; | |
541 FloatComplexNDArray m = complex_array_value (); | |
542 | |
543 OCTAVE_LOCAL_BUFFER (hsize_t, hdims, rank); | |
544 | |
545 // Octave uses column-major, while HDF5 uses row-major ordering | |
546 for (int i = 0; i < rank; i++) | |
547 hdims[i] = dv (rank-i-1); | |
548 | |
549 space_hid = H5Screate_simple (rank, hdims, 0); | |
550 if (space_hid < 0) return false; | |
551 | |
552 hid_t save_type_hid = H5T_NATIVE_FLOAT; | |
553 | |
554 #if HAVE_HDF5_INT2FLOAT_CONVERSIONS | |
555 // hdf5 currently doesn't support float/integer conversions | |
556 else | |
557 { | |
558 float max_val, min_val; | |
559 | |
560 if (m.all_integers (max_val, min_val)) | |
561 save_type_hid | |
562 = save_type_to_hdf5 (get_save_type (max_val, min_val)); | |
563 } | |
564 #endif /* HAVE_HDF5_INT2FLOAT_CONVERSIONS */ | |
565 | |
566 type_hid = hdf5_make_complex_type (save_type_hid); | |
567 if (type_hid < 0) | |
568 { | |
569 H5Sclose (space_hid); | |
570 return false; | |
571 } | |
572 | |
573 data_hid = H5Dcreate (loc_id, name, type_hid, space_hid, H5P_DEFAULT); | |
574 if (data_hid < 0) | |
575 { | |
576 H5Sclose (space_hid); | |
577 H5Tclose (type_hid); | |
578 return false; | |
579 } | |
580 | |
581 hid_t complex_type_hid = hdf5_make_complex_type (H5T_NATIVE_FLOAT); | |
582 if (complex_type_hid < 0) retval = false; | |
583 | |
584 if (retval) | |
585 { | |
586 FloatComplex *mtmp = m.fortran_vec (); | |
587 if (H5Dwrite (data_hid, complex_type_hid, H5S_ALL, H5S_ALL, H5P_DEFAULT, | |
588 mtmp) < 0) | |
589 { | |
590 H5Tclose (complex_type_hid); | |
591 retval = false; | |
592 } | |
593 } | |
594 | |
595 H5Tclose (complex_type_hid); | |
596 H5Dclose (data_hid); | |
597 H5Tclose (type_hid); | |
598 H5Sclose (space_hid); | |
599 | |
600 return retval; | |
601 } | |
602 | |
603 bool | |
604 octave_float_complex_matrix::load_hdf5 (hid_t loc_id, const char *name, | |
605 bool /* have_h5giterate_bug */) | |
606 { | |
607 bool retval = false; | |
608 | |
609 dim_vector dv; | |
610 int empty = load_hdf5_empty (loc_id, name, dv); | |
611 if (empty > 0) | |
612 matrix.resize(dv); | |
613 if (empty) | |
614 return (empty > 0); | |
615 | |
616 hid_t data_hid = H5Dopen (loc_id, name); | |
617 hid_t type_hid = H5Dget_type (data_hid); | |
618 | |
619 hid_t complex_type = hdf5_make_complex_type (H5T_NATIVE_FLOAT); | |
620 | |
621 if (! hdf5_types_compatible (type_hid, complex_type)) | |
622 { | |
623 H5Tclose (complex_type); | |
624 H5Dclose (data_hid); | |
625 return false; | |
626 } | |
627 | |
628 hid_t space_id = H5Dget_space (data_hid); | |
629 | |
630 hsize_t rank = H5Sget_simple_extent_ndims (space_id); | |
631 | |
632 if (rank < 1) | |
633 { | |
634 H5Tclose (complex_type); | |
635 H5Sclose (space_id); | |
636 H5Dclose (data_hid); | |
637 return false; | |
638 } | |
639 | |
640 OCTAVE_LOCAL_BUFFER (hsize_t, hdims, rank); | |
641 OCTAVE_LOCAL_BUFFER (hsize_t, maxdims, rank); | |
642 | |
643 H5Sget_simple_extent_dims (space_id, hdims, maxdims); | |
644 | |
645 // Octave uses column-major, while HDF5 uses row-major ordering | |
646 if (rank == 1) | |
647 { | |
648 dv.resize (2); | |
649 dv(0) = 1; | |
650 dv(1) = hdims[0]; | |
651 } | |
652 else | |
653 { | |
654 dv.resize (rank); | |
655 for (hsize_t i = 0, j = rank - 1; i < rank; i++, j--) | |
656 dv(j) = hdims[i]; | |
657 } | |
658 | |
659 FloatComplexNDArray m (dv); | |
660 FloatComplex *reim = m.fortran_vec (); | |
661 if (H5Dread (data_hid, complex_type, H5S_ALL, H5S_ALL, H5P_DEFAULT, | |
662 reim) >= 0) | |
663 { | |
664 retval = true; | |
665 matrix = m; | |
666 } | |
667 | |
668 H5Tclose (complex_type); | |
669 H5Sclose (space_id); | |
670 H5Dclose (data_hid); | |
671 | |
672 return retval; | |
673 } | |
674 | |
675 #endif | |
676 | |
677 void | |
678 octave_float_complex_matrix::print_raw (std::ostream& os, | |
679 bool pr_as_read_syntax) const | |
680 { | |
681 octave_print_internal (os, matrix, pr_as_read_syntax, | |
682 current_print_indent_level ()); | |
683 } | |
684 | |
685 mxArray * | |
686 octave_float_complex_matrix::as_mxArray (void) const | |
687 { | |
688 mxArray *retval = new mxArray (mxSINGLE_CLASS, dims (), mxCOMPLEX); | |
689 | |
690 float *pr = static_cast<float *> (retval->get_data ()); | |
691 float *pi = static_cast<float *> (retval->get_imag_data ()); | |
692 | |
693 mwSize nel = numel (); | |
694 | |
695 const FloatComplex *p = matrix.data (); | |
696 | |
697 for (mwIndex i = 0; i < nel; i++) | |
698 { | |
699 pr[i] = std::real (p[i]); | |
700 pi[i] = std::imag (p[i]); | |
701 } | |
702 | |
703 return retval; | |
704 } | |
705 | |
706 static float | |
707 xabs (const FloatComplex& x) | |
708 { | |
709 return (xisinf (x.real ()) || xisinf (x.imag ())) ? octave_Inf : abs (x); | |
710 } | |
711 | |
712 static float | |
713 ximag (const FloatComplex& x) | |
714 { | |
715 return x.imag (); | |
716 } | |
717 | |
718 static float | |
719 xreal (const FloatComplex& x) | |
720 { | |
721 return x.real (); | |
722 } | |
723 | |
724 static bool | |
725 any_element_less_than (const FloatNDArray& a, float val) | |
726 { | |
727 octave_idx_type len = a.length (); | |
728 const float *m = a.fortran_vec (); | |
729 | |
730 for (octave_idx_type i = 0; i < len; i++) | |
731 { | |
732 OCTAVE_QUIT; | |
733 | |
734 if (m[i] < val) | |
735 return true; | |
736 } | |
737 | |
738 return false; | |
739 } | |
740 | |
741 static bool | |
742 any_element_greater_than (const FloatNDArray& a, float val) | |
743 { | |
744 octave_idx_type len = a.length (); | |
745 const float *m = a.fortran_vec (); | |
746 | |
747 for (octave_idx_type i = 0; i < len; i++) | |
748 { | |
749 OCTAVE_QUIT; | |
750 | |
751 if (m[i] > val) | |
752 return true; | |
753 } | |
754 | |
755 return false; | |
756 } | |
757 | |
758 #define ARRAY_MAPPER(MAP, AMAP, FCN) \ | |
759 octave_value \ | |
760 octave_float_complex_matrix::MAP (void) const \ | |
761 { \ | |
762 static AMAP cmap = FCN; \ | |
763 return matrix.map (cmap); \ | |
764 } | |
765 | |
766 #define DARRAY_MAPPER(MAP, AMAP, FCN) \ | |
767 octave_value \ | |
768 octave_float_complex_matrix::MAP (void) const \ | |
769 { \ | |
770 static FloatComplexNDArray::dmapper dmap = ximag; \ | |
771 NDArray m = matrix.map (dmap); \ | |
772 if (m.all_elements_are_zero ()) \ | |
773 { \ | |
774 dmap = xreal; \ | |
775 m = matrix.map (dmap); \ | |
776 static AMAP cmap = FCN; \ | |
777 return m.map (cmap); \ | |
778 } \ | |
779 else \ | |
780 { \ | |
781 error ("%s: not defined for complex arguments", #MAP); \ | |
782 return octave_value (); \ | |
783 } \ | |
784 } | |
785 | |
786 #define CD_ARRAY_MAPPER(MAP, RFCN, CFCN, L1, L2) \ | |
787 octave_value \ | |
788 octave_float_complex_matrix::MAP (void) const \ | |
789 { \ | |
790 static FloatComplexNDArray::dmapper idmap = ximag; \ | |
791 NDArray m = matrix.map (idmap); \ | |
792 if (m.all_elements_are_zero ()) \ | |
793 { \ | |
794 static FloatComplexNDArray::dmapper rdmap = xreal; \ | |
795 m = matrix.map (rdmap); \ | |
796 static NDArray::dmapper dmap = RFCN; \ | |
797 static NDArray::cmapper cmap = CFCN; \ | |
798 return (any_element_less_than (m, L1) \ | |
799 ? octave_value (m.map (cmap)) \ | |
800 : (any_element_greater_than (m, L2) \ | |
801 ? octave_value (m.map (cmap)) \ | |
802 : octave_value (m.map (dmap)))); \ | |
803 } \ | |
804 else \ | |
805 { \ | |
806 /*error ("%s: not defined for complex arguments", #MAP); */ \ | |
807 return octave_value (m); \ | |
808 } \ | |
809 } | |
810 | |
811 DARRAY_MAPPER (erf, NDArray::dmapper, ::erf) | |
812 DARRAY_MAPPER (erfc, NDArray::dmapper, ::erfc) | |
813 DARRAY_MAPPER (gamma, NDArray::dmapper, xgamma) | |
814 CD_ARRAY_MAPPER (lgamma, xlgamma, xlgamma, 0.0, octave_Inf) | |
815 | |
816 ARRAY_MAPPER (abs, FloatComplexNDArray::dmapper, xabs) | |
817 ARRAY_MAPPER (acos, FloatComplexNDArray::cmapper, ::acos) | |
818 ARRAY_MAPPER (acosh, FloatComplexNDArray::cmapper, ::acosh) | |
819 ARRAY_MAPPER (angle, FloatComplexNDArray::dmapper, std::arg) | |
820 ARRAY_MAPPER (arg, FloatComplexNDArray::dmapper, std::arg) | |
821 ARRAY_MAPPER (asin, FloatComplexNDArray::cmapper, ::asin) | |
822 ARRAY_MAPPER (asinh, FloatComplexNDArray::cmapper, ::asinh) | |
823 ARRAY_MAPPER (atan, FloatComplexNDArray::cmapper, ::atan) | |
824 ARRAY_MAPPER (atanh, FloatComplexNDArray::cmapper, ::atanh) | |
825 ARRAY_MAPPER (ceil, FloatComplexNDArray::cmapper, ::ceil) | |
826 ARRAY_MAPPER (conj, FloatComplexNDArray::cmapper, std::conj) | |
827 ARRAY_MAPPER (cos, FloatComplexNDArray::cmapper, std::cos) | |
828 ARRAY_MAPPER (cosh, FloatComplexNDArray::cmapper, std::cosh) | |
829 ARRAY_MAPPER (exp, FloatComplexNDArray::cmapper, std::exp) | |
830 ARRAY_MAPPER (expm1, FloatComplexNDArray::cmapper, ::expm1f) | |
831 ARRAY_MAPPER (fix, FloatComplexNDArray::cmapper, ::fix) | |
832 ARRAY_MAPPER (floor, FloatComplexNDArray::cmapper, ::floor) | |
833 ARRAY_MAPPER (imag, FloatComplexNDArray::dmapper, ximag) | |
834 ARRAY_MAPPER (log, FloatComplexNDArray::cmapper, std::log) | |
835 ARRAY_MAPPER (log2, FloatComplexNDArray::cmapper, xlog2) | |
836 ARRAY_MAPPER (log10, FloatComplexNDArray::cmapper, std::log10) | |
837 ARRAY_MAPPER (log1p, FloatComplexNDArray::cmapper, ::log1pf) | |
838 ARRAY_MAPPER (real, FloatComplexNDArray::dmapper, xreal) | |
839 ARRAY_MAPPER (round, FloatComplexNDArray::cmapper, xround) | |
840 ARRAY_MAPPER (roundb, FloatComplexNDArray::cmapper, xroundb) | |
841 ARRAY_MAPPER (signum, FloatComplexNDArray::cmapper, ::signum) | |
842 ARRAY_MAPPER (sin, FloatComplexNDArray::cmapper, std::sin) | |
843 ARRAY_MAPPER (sinh, FloatComplexNDArray::cmapper, std::sinh) | |
844 ARRAY_MAPPER (sqrt, FloatComplexNDArray::cmapper, std::sqrt) | |
845 ARRAY_MAPPER (tan, FloatComplexNDArray::cmapper, std::tan) | |
846 ARRAY_MAPPER (tanh, FloatComplexNDArray::cmapper, std::tanh) | |
847 ARRAY_MAPPER (finite, FloatComplexNDArray::bmapper, xfinite) | |
848 ARRAY_MAPPER (isinf, FloatComplexNDArray::bmapper, xisinf) | |
849 ARRAY_MAPPER (isna, FloatComplexNDArray::bmapper, octave_is_NA) | |
850 ARRAY_MAPPER (isnan, FloatComplexNDArray::bmapper, xisnan) | |
851 | |
852 /* | |
853 ;;; Local Variables: *** | |
854 ;;; mode: C++ *** | |
855 ;;; End: *** | |
856 */ |