Mercurial > hg > octave-max

comparison src/data.cc @ 7789:82be108cc558

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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 fbe27e477578
children df9519e9990c
comparison
equal deleted inserted replaced
7788:45f5faba05a2 7789:82be108cc558
49 #include "error.h" 49 #include "error.h"
50 #include "gripes.h" 50 #include "gripes.h"
51 #include "oct-map.h" 51 #include "oct-map.h"
52 #include "oct-obj.h" 52 #include "oct-obj.h"
53 #include "ov.h" 53 #include "ov.h"
54 #include "ov-float.h"
54 #include "ov-complex.h" 55 #include "ov-complex.h"
56 #include "ov-flt-complex.h"
55 #include "ov-cx-mat.h" 57 #include "ov-cx-mat.h"
58 #include "ov-flt-cx-mat.h"
59 #include "ov-cx-sparse.h"
56 #include "parse.h" 60 #include "parse.h"
57 #include "pt-mat.h" 61 #include "pt-mat.h"
58 #include "utils.h" 62 #include "utils.h"
59 #include "variables.h" 63 #include "variables.h"
60 #include "pager.h" 64 #include "pager.h"
127 } 131 }
128 132
129 // These mapping functions may also be useful in other places, eh? 133 // These mapping functions may also be useful in other places, eh?
130 134
131 typedef double (*d_dd_fcn) (double, double); 135 typedef double (*d_dd_fcn) (double, double);
136 typedef float (*f_ff_fcn) (float, float);
132 137
133 static NDArray 138 static NDArray
134 map_d_m (d_dd_fcn f, double x, const NDArray& y) 139 map_d_m (d_dd_fcn f, double x, const NDArray& y)
135 { 140 {
136 NDArray retval (y.dims ()); 141 NDArray retval (y.dims ());
147 } 152 }
148 153
149 return retval; 154 return retval;
150 } 155 }
151 156
157 static FloatNDArray
158 map_f_fm (f_ff_fcn f, float x, const FloatNDArray& y)
159 {
160 FloatNDArray retval (y.dims ());
161 float *r_data = retval.fortran_vec ();
162
163 const float *y_data = y.data ();
164
165 octave_idx_type nel = y.numel ();
166
167 for (octave_idx_type i = 0; i < nel; i++)
168 {
169 OCTAVE_QUIT;
170 r_data[i] = f (x, y_data[i]);
171 }
172
173 return retval;
174 }
175
152 static NDArray 176 static NDArray
153 map_m_d (d_dd_fcn f, const NDArray& x, double y) 177 map_m_d (d_dd_fcn f, const NDArray& x, double y)
154 { 178 {
155 NDArray retval (x.dims ()); 179 NDArray retval (x.dims ());
156 double *r_data = retval.fortran_vec (); 180 double *r_data = retval.fortran_vec ();
166 } 190 }
167 191
168 return retval; 192 return retval;
169 } 193 }
170 194
195 static FloatNDArray
196 map_fm_f (f_ff_fcn f, const FloatNDArray& x, float y)
197 {
198 FloatNDArray retval (x.dims ());
199 float *r_data = retval.fortran_vec ();
200
201 const float *x_data = x.data ();
202
203 octave_idx_type nel = x.numel ();
204
205 for (octave_idx_type i = 0; i < nel; i++)
206 {
207 OCTAVE_QUIT;
208 r_data[i] = f (x_data[i], y);
209 }
210
211 return retval;
212 }
213
171 static NDArray 214 static NDArray
172 map_m_m (d_dd_fcn f, const NDArray& x, const NDArray& y) 215 map_m_m (d_dd_fcn f, const NDArray& x, const NDArray& y)
173 { 216 {
174 assert (x.dims () == y.dims ()); 217 assert (x.dims () == y.dims ());
175 218
176 NDArray retval (x.dims ()); 219 NDArray retval (x.dims ());
177 double *r_data = retval.fortran_vec (); 220 double *r_data = retval.fortran_vec ();
178 221
179 const double *x_data = x.data (); 222 const double *x_data = x.data ();
180 const double *y_data = y.data (); 223 const double *y_data = y.data ();
224
225 octave_idx_type nel = x.numel ();
226
227 for (octave_idx_type i = 0; i < nel; i++)
228 {
229 OCTAVE_QUIT;
230 r_data[i] = f (x_data[i], y_data[i]);
231 }
232
233 return retval;
234 }
235
236 static FloatNDArray
237 map_fm_fm (f_ff_fcn f, const FloatNDArray& x, const FloatNDArray& y)
238 {
239 assert (x.dims () == y.dims ());
240
241 FloatNDArray retval (x.dims ());
242 float *r_data = retval.fortran_vec ();
243
244 const float *x_data = x.data ();
245 const float *y_data = y.data ();
181 246
182 octave_idx_type nel = x.numel (); 247 octave_idx_type nel = x.numel ();
183 248
184 for (octave_idx_type i = 0; i < nel; i++) 249 for (octave_idx_type i = 0; i < nel; i++)
185 { 250 {
436 dim_vector x_dims = arg_x.dims (); 501 dim_vector x_dims = arg_x.dims ();
437 502
438 bool y_is_scalar = y_dims.all_ones (); 503 bool y_is_scalar = y_dims.all_ones ();
439 bool x_is_scalar = x_dims.all_ones (); 504 bool x_is_scalar = x_dims.all_ones ();
440 505
506 bool is_float = arg_y.is_single_type () || arg_x.is_single_type ();
507
441 if (y_is_scalar && x_is_scalar) 508 if (y_is_scalar && x_is_scalar)
442 { 509 {
443 double y = arg_y.double_value (); 510 if (is_float)
444 511 {
445 if (! error_state) 512 float y = arg_y.float_value ();
446 {
447 double x = arg_x.double_value ();
448 513
449 if (! error_state) 514 if (! error_state)
450 retval = atan2 (y, x); 515 {
516 float x = arg_x.float_value ();
517
518 if (! error_state)
519 retval = atan2f (y, x);
520 }
521 }
522 else
523 {
524 double y = arg_y.double_value ();
525
526 if (! error_state)
527 {
528 double x = arg_x.double_value ();
529
530 if (! error_state)
531 retval = atan2 (y, x);
532 }
451 } 533 }
452 } 534 }
453 else if (y_is_scalar) 535 else if (y_is_scalar)
454 { 536 {
455 double y = arg_y.double_value (); 537 if (is_float)
456 538 {
457 if (! error_state) 539 float y = arg_y.float_value ();
458 {
459 // Even if x is sparse return a full matrix here
460 NDArray x = arg_x.array_value ();
461 540
462 if (! error_state) 541 if (! error_state)
463 retval = map_d_m (atan2, y, x); 542 {
543 // Even if x is sparse return a full matrix here
544 FloatNDArray x = arg_x.float_array_value ();
545
546 if (! error_state)
547 retval = map_f_fm (atan2f, y, x);
548 }
549 }
550 else
551 {
552 double y = arg_y.double_value ();
553
554 if (! error_state)
555 {
556 // Even if x is sparse return a full matrix here
557 NDArray x = arg_x.array_value ();
558
559 if (! error_state)
560 retval = map_d_m (atan2, y, x);
561 }
464 } 562 }
465 } 563 }
466 else if (x_is_scalar) 564 else if (x_is_scalar)
467 { 565 {
468 if (arg_y.is_sparse_type ()) 566 if (arg_y.is_sparse_type ())
475 573
476 if (! error_state) 574 if (! error_state)
477 retval = map_s_d (atan2, y, x); 575 retval = map_s_d (atan2, y, x);
478 } 576 }
479 } 577 }
578 else if (is_float)
579 {
580 FloatNDArray y = arg_y.float_array_value ();
581
582 if (! error_state)
583 {
584 float x = arg_x.float_value ();
585
586 if (! error_state)
587 retval = map_fm_f (atan2f, y, x);
588 }
589 }
480 else 590 else
481 { 591 {
482 NDArray y = arg_y.array_value (); 592 NDArray y = arg_y.array_value ();
483 593
484 if (! error_state) 594 if (! error_state)
501 { 611 {
502 SparseMatrix x = arg_x.sparse_matrix_value (); 612 SparseMatrix x = arg_x.sparse_matrix_value ();
503 613
504 if (! error_state) 614 if (! error_state)
505 retval = map_s_s (atan2, y, x); 615 retval = map_s_s (atan2, y, x);
616 }
617 }
618 else if (is_float)
619 {
620 FloatNDArray y = arg_y.array_value ();
621
622 if (! error_state)
623 {
624 FloatNDArray x = arg_x.array_value ();
625
626 if (! error_state)
627 retval = map_fm_fm (atan2f, y, x);
506 } 628 }
507 } 629 }
508 else 630 else
509 { 631 {
510 NDArray y = arg_y.array_value (); 632 NDArray y = arg_y.array_value ();
562 dim_vector y_dims = arg_y.dims (); 684 dim_vector y_dims = arg_y.dims ();
563 685
564 bool x_is_scalar = x_dims.all_ones (); 686 bool x_is_scalar = x_dims.all_ones ();
565 bool y_is_scalar = y_dims.all_ones (); 687 bool y_is_scalar = y_dims.all_ones ();
566 688
689 bool is_float = arg_y.is_single_type () || arg_x.is_single_type ();
690
567 if (y_is_scalar && x_is_scalar) 691 if (y_is_scalar && x_is_scalar)
568 { 692 {
569 double x; 693 if (is_float)
570 if (arg_x.is_complex_type ()) 694 {
571 x = abs (arg_x.complex_value ()); 695 float x;
696 if (arg_x.is_complex_type ())
697 x = abs (arg_x.float_complex_value ());
698 else
699 x = arg_x.float_value ();
700
701 if (! error_state)
702 {
703 float y;
704 if (arg_y.is_complex_type ())
705 y = abs (arg_y.float_complex_value ());
706 else
707 y = arg_y.float_value ();
708
709 if (! error_state)
710 retval = hypotf (x, y);
711 }
712 }
572 else 713 else
573 x = arg_x.double_value (); 714 {
574 715 double x;
575 if (! error_state) 716 if (arg_x.is_complex_type ())
576 { 717 x = abs (arg_x.complex_value ());
577 double y;
578 if (arg_y.is_complex_type ())
579 y = abs (arg_y.complex_value ());
580 else 718 else
581 y = arg_y.double_value (); 719 x = arg_x.double_value ();
582 720
583 if (! error_state) 721 if (! error_state)
584 retval = hypot (x, y); 722 {
723 double y;
724 if (arg_y.is_complex_type ())
725 y = abs (arg_y.complex_value ());
726 else
727 y = arg_y.double_value ();
728
729 if (! error_state)
730 retval = hypot (x, y);
731 }
585 } 732 }
586 } 733 }
587 else if (y_is_scalar) 734 else if (y_is_scalar)
588 { 735 {
589 NDArray x; 736 if (is_float)
590 if (arg_x.is_complex_type ()) 737 {
591 x = arg_x.complex_array_value ().abs (); 738 FloatNDArray x;
739 if (arg_x.is_complex_type ())
740 x = arg_x.float_complex_array_value ().abs ();
741 else
742 x = arg_x.float_array_value ();
743
744 if (! error_state)
745 {
746 float y;
747 if (arg_y.is_complex_type ())
748 y = abs (arg_y.float_complex_value ());
749 else
750 y = arg_y.float_value ();
751
752 if (! error_state)
753 retval = map_fm_f (hypotf, x, y);
754 }
755 }
592 else 756 else
593 x = arg_x.array_value (); 757 {
594 758 NDArray x;
595 if (! error_state) 759 if (arg_x.is_complex_type ())
596 { 760 x = arg_x.complex_array_value ().abs ();
597 double y;
598 if (arg_y.is_complex_type ())
599 y = abs (arg_y.complex_value ());
600 else 761 else
601 y = arg_y.double_value (); 762 x = arg_x.array_value ();
602 763
603 if (! error_state) 764 if (! error_state)
604 retval = map_m_d (hypot, x, y); 765 {
766 double y;
767 if (arg_y.is_complex_type ())
768 y = abs (arg_y.complex_value ());
769 else
770 y = arg_y.double_value ();
771
772 if (! error_state)
773 retval = map_m_d (hypot, x, y);
774 }
605 } 775 }
606 } 776 }
607 else if (x_is_scalar) 777 else if (x_is_scalar)
608 { 778 {
609 double x; 779 if (is_float)
610 if (arg_x.is_complex_type ()) 780 {
611 x = abs (arg_x.complex_value ()); 781 float x;
782 if (arg_x.is_complex_type ())
783 x = abs (arg_x.float_complex_value ());
784 else
785 x = arg_x.float_value ();
786
787 if (! error_state)
788 {
789 FloatNDArray y;
790 if (arg_y.is_complex_type ())
791 y = arg_y.float_complex_array_value ().abs ();
792 else
793 y = arg_y.float_array_value ();
794
795 if (! error_state)
796 retval = map_f_fm (hypotf, x, y);
797 }
798 }
612 else 799 else
613 x = arg_x.double_value (); 800 {
614 801 double x;
615 if (! error_state) 802 if (arg_x.is_complex_type ())
616 { 803 x = abs (arg_x.complex_value ());
617 NDArray y;
618 if (arg_y.is_complex_type ())
619 y = arg_y.complex_array_value ().abs ();
620 else 804 else
621 y = arg_y.array_value (); 805 x = arg_x.double_value ();
622 806
623 if (! error_state) 807 if (! error_state)
624 retval = map_d_m (hypot, x, y); 808 {
809 NDArray y;
810 if (arg_y.is_complex_type ())
811 y = arg_y.complex_array_value ().abs ();
812 else
813 y = arg_y.array_value ();
814
815 if (! error_state)
816 retval = map_d_m (hypot, x, y);
817 }
625 } 818 }
626 } 819 }
627 else if (y_dims == x_dims) 820 else if (y_dims == x_dims)
628 { 821 {
629 if (arg_x.is_sparse_type () && arg_y.is_sparse_type ()) 822 if (arg_x.is_sparse_type () && arg_y.is_sparse_type ())
644 837
645 if (! error_state) 838 if (! error_state)
646 retval = map_s_s (hypot, x, y); 839 retval = map_s_s (hypot, x, y);
647 } 840 }
648 } 841 }
842 else if (is_float)
843 {
844 FloatNDArray x;
845 if (arg_x.is_complex_type ())
846 x = arg_x.float_complex_array_value ().abs ();
847 else
848 x = arg_x.float_array_value ();
849
850 if (! error_state)
851 {
852 FloatNDArray y;
853 if (arg_y.is_complex_type ())
854 y = arg_y.float_complex_array_value ().abs ();
855 else
856 y = arg_y.float_array_value ();
857
858 if (! error_state)
859 retval = map_fm_fm (hypotf, x, y);
860 }
861 }
649 else 862 else
650 { 863 {
651 NDArray x; 864 NDArray x;
652 if (arg_x.is_complex_type ()) 865 if (arg_x.is_complex_type ())
653 x = arg_x.complex_array_value ().abs (); 866 x = arg_x.complex_array_value ().abs ();
682 %!assert (size (hypot (rand (2, 3, 4), zeros (2, 3, 4))), [2, 3, 4]) 895 %!assert (size (hypot (rand (2, 3, 4), zeros (2, 3, 4))), [2, 3, 4])
683 %!assert (size (hypot (rand (2, 3, 4), 1)), [2, 3, 4]) 896 %!assert (size (hypot (rand (2, 3, 4), 1)), [2, 3, 4])
684 %!assert (size (hypot (1, rand (2, 3, 4))), [2, 3, 4]) 897 %!assert (size (hypot (1, rand (2, 3, 4))), [2, 3, 4])
685 %!assert (size (hypot (1, 2)), [1, 1]) 898 %!assert (size (hypot (1, 2)), [1, 1])
686 %!assert (hypot (1:10, 1:10), sqrt(2) * [1:10], 16*eps) 899 %!assert (hypot (1:10, 1:10), sqrt(2) * [1:10], 16*eps)
900 %!assert (hypot (single(1:10), single(1:10)), single(sqrt(2) * [1:10]));
687 */ 901 */
688 902
689 template<typename T, typename ET> 903 template<typename T, typename ET>
690 void 904 void
691 map_2_xlog2 (const Array<T>& x, Array<T>& f, Array<ET>& e) 905 map_2_xlog2 (const Array<T>& x, Array<T>& f, Array<ET>& e)
715 929
716 if (args.length () == 1) 930 if (args.length () == 1)
717 { 931 {
718 if (nargout < 2) 932 if (nargout < 2)
719 retval(0) = args(0).log2 (); 933 retval(0) = args(0).log2 ();
934 else if (args(0).is_single_type ())
935 {
936 if (args(0).is_real_type ())
937 {
938 FloatNDArray f;
939 FloatNDArray x = args(0).float_array_value ();
940 // FIXME -- should E be an int value?
941 FloatMatrix e;
942 map_2_xlog2 (x, f, e);
943 retval (1) = e;
944 retval (0) = f;
945 }
946 else if (args(0).is_complex_type ())
947 {
948 FloatComplexNDArray f;
949 FloatComplexNDArray x = args(0).float_complex_array_value ();
950 // FIXME -- should E be an int value?
951 FloatNDArray e;
952 map_2_xlog2 (x, f, e);
953 retval (1) = e;
954 retval (0) = f;
955 }
956 }
720 else if (args(0).is_real_type ()) 957 else if (args(0).is_real_type ())
721 { 958 {
722 NDArray f; 959 NDArray f;
723 NDArray x = args(0).array_value (); 960 NDArray x = args(0).array_value ();
724 // FIXME -- should E be an int value? 961 // FIXME -- should E be an int value?
785 dim_vector x_dims = arg_x.dims (); 1022 dim_vector x_dims = arg_x.dims ();
786 1023
787 bool y_is_scalar = y_dims.all_ones (); 1024 bool y_is_scalar = y_dims.all_ones ();
788 bool x_is_scalar = x_dims.all_ones (); 1025 bool x_is_scalar = x_dims.all_ones ();
789 1026
1027 bool is_float = arg_y.is_single_type () || arg_x.is_single_type ();
1028
790 if (y_is_scalar && x_is_scalar) 1029 if (y_is_scalar && x_is_scalar)
791 { 1030 {
792 double y = arg_y.double_value (); 1031 if (is_float)
793 1032 {
794 if (! error_state) 1033 float y = arg_y.float_value ();
795 {
796 double x = arg_x.double_value ();
797 1034
798 if (! error_state) 1035 if (! error_state)
799 retval = fmod (x, y); 1036 {
1037 float x = arg_x.float_value ();
1038
1039 if (! error_state)
1040 retval = fmod (x, y);
1041 }
1042 }
1043 else
1044 {
1045 double y = arg_y.double_value ();
1046
1047 if (! error_state)
1048 {
1049 double x = arg_x.double_value ();
1050
1051 if (! error_state)
1052 retval = fmod (x, y);
1053 }
800 } 1054 }
801 } 1055 }
802 else if (y_is_scalar) 1056 else if (y_is_scalar)
803 { 1057 {
804 double y = arg_y.double_value (); 1058 if (is_float)
805 1059 {
806 if (! error_state) 1060 float y = arg_y.float_value ();
807 { 1061
808 if (arg_x.is_sparse_type ()) 1062 if (! error_state)
809 { 1063 {
810 SparseMatrix x = arg_x.sparse_matrix_value (); 1064 FloatNDArray x = arg_x.float_array_value ();
811 1065
812 if (! error_state) 1066 if (! error_state)
813 retval = map_s_d (fmod, x, y); 1067 retval = map_fm_f (fmodf, x, y);
814 } 1068 }
815 else 1069 }
816 { 1070 else
817 NDArray x = arg_x.array_value (); 1071 {
818 1072 double y = arg_y.double_value ();
819 if (! error_state) 1073
820 retval = map_m_d (fmod, x, y); 1074 if (! error_state)
1075 {
1076 if (arg_x.is_sparse_type ())
1077 {
1078 SparseMatrix x = arg_x.sparse_matrix_value ();
1079
1080 if (! error_state)
1081 retval = map_s_d (fmod, x, y);
1082 }
1083 else
1084 {
1085 NDArray x = arg_x.array_value ();
1086
1087 if (! error_state)
1088 retval = map_m_d (fmod, x, y);
1089 }
821 } 1090 }
822 } 1091 }
823 } 1092 }
824 else if (x_is_scalar) 1093 else if (x_is_scalar)
825 { 1094 {
833 1102
834 if (! error_state) 1103 if (! error_state)
835 retval = map_d_s (fmod, x, y); 1104 retval = map_d_s (fmod, x, y);
836 } 1105 }
837 } 1106 }
1107 else if (is_float)
1108 {
1109 FloatNDArray y = arg_y.float_array_value ();
1110
1111 if (! error_state)
1112 {
1113 float x = arg_x.float_value ();
1114
1115 if (! error_state)
1116 retval = map_f_fm (fmodf, x, y);
1117 }
1118 }
838 else 1119 else
839 { 1120 {
840 NDArray y = arg_y.array_value (); 1121 NDArray y = arg_y.array_value ();
841 1122
842 if (! error_state) 1123 if (! error_state)
858 { 1139 {
859 SparseMatrix x = arg_x.sparse_matrix_value (); 1140 SparseMatrix x = arg_x.sparse_matrix_value ();
860 1141
861 if (! error_state) 1142 if (! error_state)
862 retval = map_s_s (fmod, x, y); 1143 retval = map_s_s (fmod, x, y);
1144 }
1145 }
1146 else if (is_float)
1147 {
1148 FloatNDArray y = arg_y.float_array_value ();
1149
1150 if (! error_state)
1151 {
1152 FloatNDArray x = arg_x.float_array_value ();
1153
1154 if (! error_state)
1155 retval = map_fm_fm (fmodf, x, y);
863 } 1156 }
864 } 1157 }
865 else 1158 else
866 { 1159 {
867 NDArray y = arg_y.array_value (); 1160 NDArray y = arg_y.array_value ();
889 %!assert (size (fmod (rand (2, 3, 4), zeros (2, 3, 4))), [2, 3, 4]) 1182 %!assert (size (fmod (rand (2, 3, 4), zeros (2, 3, 4))), [2, 3, 4])
890 %!assert (size (fmod (rand (2, 3, 4), 1)), [2, 3, 4]) 1183 %!assert (size (fmod (rand (2, 3, 4), 1)), [2, 3, 4])
891 %!assert (size (fmod (1, rand (2, 3, 4))), [2, 3, 4]) 1184 %!assert (size (fmod (1, rand (2, 3, 4))), [2, 3, 4])
892 %!assert (size (fmod (1, 2)), [1, 1]) 1185 %!assert (size (fmod (1, 2)), [1, 1])
893 */ 1186 */
1187
1188 // FIXME Need to convert the reduction functions of this file for single precision
894 1189
895 #define NATIVE_REDUCTION_1(FCN, TYPE, DIM) \ 1190 #define NATIVE_REDUCTION_1(FCN, TYPE, DIM) \
896 (arg.is_ ## TYPE ## _type ()) \ 1191 (arg.is_ ## TYPE ## _type ()) \
897 { \ 1192 { \
898 TYPE ## NDArray tmp = arg. TYPE ##_array_value (); \ 1193 TYPE ## NDArray tmp = arg. TYPE ##_array_value (); \
965 else if NATIVE_REDUCTION_1 (FCN, bool, dim) \ 1260 else if NATIVE_REDUCTION_1 (FCN, bool, dim) \
966 else if (arg.is_char_matrix ()) \ 1261 else if (arg.is_char_matrix ()) \
967 { \ 1262 { \
968 error (#FCN, ": invalid char type"); \ 1263 error (#FCN, ": invalid char type"); \
969 } \ 1264 } \
1265 else if (arg.is_single_type ()) \
1266 { \
1267 if (arg.is_complex_type ()) \
1268 { \
1269 FloatComplexNDArray tmp = \
1270 arg.float_complex_array_value (); \
1271 \
1272 if (! error_state) \
1273 retval = tmp.FCN (dim); \
1274 } \
1275 else if (arg.is_real_type ()) \
1276 { \
1277 FloatNDArray tmp = arg.float_array_value (); \
1278 \
1279 if (! error_state) \
1280 retval = tmp.FCN (dim); \
1281 } \
1282 } \
970 else if (arg.is_complex_type ()) \ 1283 else if (arg.is_complex_type ()) \
971 { \ 1284 { \
972 ComplexNDArray tmp = arg.complex_array_value (); \ 1285 ComplexNDArray tmp = arg.complex_array_value (); \
973 \ 1286 \
974 if (! error_state) \ 1287 if (! error_state) \
985 { \ 1298 { \
986 gripe_wrong_type_arg (#FCN, arg); \ 1299 gripe_wrong_type_arg (#FCN, arg); \
987 return retval; \ 1300 return retval; \
988 } \ 1301 } \
989 } \ 1302 } \
1303 else if (arg.is_single_type ()) \
1304 { \
1305 if (arg.is_real_type ()) \
1306 { \
1307 FloatNDArray tmp = arg.float_array_value (); \
1308 \
1309 if (! error_state) \
1310 retval = tmp.FCN (dim); \
1311 } \
1312 else if (arg.is_complex_type ()) \
1313 { \
1314 FloatComplexNDArray tmp = \
1315 arg.float_complex_array_value (); \
1316 \
1317 if (! error_state) \
1318 retval = tmp.FCN (dim); \
1319 } \
1320 } \
990 else if (arg.is_real_type ()) \ 1321 else if (arg.is_real_type ()) \
991 { \ 1322 { \
992 NDArray tmp = arg.array_value (); \ 1323 NDArray tmp = arg.array_value (); \
993 \ 1324 \
994 if (! error_state) \ 1325 if (! error_state) \
1041 SparseMatrix tmp = arg.sparse_matrix_value (); \ 1372 SparseMatrix tmp = arg.sparse_matrix_value (); \
1042 \ 1373 \
1043 if (! error_state) \ 1374 if (! error_state) \
1044 retval = tmp.FCN (dim); \ 1375 retval = tmp.FCN (dim); \
1045 } \ 1376 } \
1377 else if (arg.is_single_type ()) \
1378 { \
1379 FloatNDArray tmp = arg.float_array_value (); \
1380 \
1381 if (! error_state) \
1382 retval = tmp.FCN (dim); \
1383 } \
1046 else \ 1384 else \
1047 { \ 1385 { \
1048 NDArray tmp = arg.array_value (); \ 1386 NDArray tmp = arg.array_value (); \
1049 \ 1387 \
1050 if (! error_state) \ 1388 if (! error_state) \
1054 else if (arg.is_complex_type ()) \ 1392 else if (arg.is_complex_type ()) \
1055 { \ 1393 { \
1056 if (arg.is_sparse_type ()) \ 1394 if (arg.is_sparse_type ()) \
1057 { \ 1395 { \
1058 SparseComplexMatrix tmp = arg.sparse_complex_matrix_value (); \ 1396 SparseComplexMatrix tmp = arg.sparse_complex_matrix_value (); \
1397 \
1398 if (! error_state) \
1399 retval = tmp.FCN (dim); \
1400 } \
1401 else if (arg.is_single_type ()) \
1402 { \
1403 FloatComplexNDArray tmp = arg.float_complex_array_value (); \
1059 \ 1404 \
1060 if (! error_state) \ 1405 if (! error_state) \
1061 retval = tmp.FCN (dim); \ 1406 retval = tmp.FCN (dim); \
1062 } \ 1407 } \
1063 else \ 1408 else \
1848 2193
1849 if (arg.is_complex_type ()) 2194 if (arg.is_complex_type ())
1850 retval = arg; 2195 retval = arg;
1851 else 2196 else
1852 { 2197 {
1853 if (arg.numel () == 1) 2198 if (arg.is_sparse_type ())
1854 { 2199 {
1855 Complex val = arg.complex_value (); 2200 SparseComplexMatrix val = arg.sparse_complex_matrix_value ();
1856 2201
1857 if (! error_state) 2202 if (! error_state)
1858 retval = octave_value (new octave_complex (val)); 2203 retval = octave_value (new octave_sparse_complex_matrix (val));
2204 }
2205 else if (arg.is_single_type ())
2206 {
2207 if (arg.numel () == 1)
2208 {
2209 FloatComplex val = arg.float_complex_value ();
2210
2211 if (! error_state)
2212 retval = octave_value (new octave_float_complex (val));
2213 }
2214 else
2215 {
2216 FloatComplexNDArray val = arg.float_complex_array_value ();
2217
2218 if (! error_state)
2219 retval = octave_value (new octave_float_complex_matrix (val));
2220 }
1859 } 2221 }
1860 else 2222 else
1861 { 2223 {
1862 ComplexNDArray val = arg.complex_array_value (); 2224 if (arg.numel () == 1)
1863 2225 {
1864 if (! error_state) 2226 Complex val = arg.complex_value ();
1865 retval = octave_value (new octave_complex_matrix (val)); 2227
2228 if (! error_state)
2229 retval = octave_value (new octave_complex (val));
2230 }
2231 else
2232 {
2233 ComplexNDArray val = arg.complex_array_value ();
2234
2235 if (! error_state)
2236 retval = octave_value (new octave_complex_matrix (val));
2237 }
1866 } 2238 }
1867 2239
1868 if (error_state) 2240 if (error_state)
1869 error ("complex: invalid conversion"); 2241 error ("complex: invalid conversion");
1870 } 2242 }
1872 else if (nargin == 2) 2244 else if (nargin == 2)
1873 { 2245 {
1874 octave_value re = args(0); 2246 octave_value re = args(0);
1875 octave_value im = args(1); 2247 octave_value im = args(1);
1876 2248
1877 if (re.numel () == 1) 2249 if (re.is_sparse_type () && im.is_sparse_type ())
2250 {
2251 const SparseMatrix re_val = re.sparse_matrix_value ();
2252 const SparseMatrix im_val = im.sparse_matrix_value ();
2253
2254 if (!error_state)
2255 {
2256 if (re.numel () == 1)
2257 {
2258 SparseComplexMatrix result;
2259 if (re_val.nnz () == 0)
2260 result = Complex(0, 1) * SparseComplexMatrix (im_val);
2261 else
2262 {
2263 result = SparseComplexMatrix (im_val.dims (), re_val (0));
2264 octave_idx_type nr = im_val.rows ();
2265 octave_idx_type nc = im_val.cols ();
2266
2267 for (octave_idx_type j = 0; j < nc; j++)
2268 {
2269 octave_idx_type off = j * nr;
2270 for (octave_idx_type i = im_val.cidx(j);
2271 i < im_val.cidx(j + 1); i++)
2272 result.data (im_val.ridx(i) + off) =
2273 result.data (im_val.ridx(i) + off) +
2274 Complex (0, im_val.data (i));
2275 }
2276 }
2277 retval = octave_value (new octave_sparse_complex_matrix (result));
2278 }
2279 else if (im.numel () == 1)
2280 {
2281 SparseComplexMatrix result;
2282 if (im_val.nnz () == 0)
2283 result = SparseComplexMatrix (re_val);
2284 else
2285 {
2286 result = SparseComplexMatrix (re_val.rows(), re_val.cols(), Complex(0, im_val (0)));
2287 octave_idx_type nr = re_val.rows ();
2288 octave_idx_type nc = re_val.cols ();
2289
2290 for (octave_idx_type j = 0; j < nc; j++)
2291 {
2292 octave_idx_type off = j * nr;
2293 for (octave_idx_type i = re_val.cidx(j);
2294 i < re_val.cidx(j + 1); i++)
2295 result.data (re_val.ridx(i) + off) =
2296 result.data (re_val.ridx(i) + off) +
2297 re_val.data (i);
2298 }
2299 }
2300 retval = octave_value (new octave_sparse_complex_matrix (result));
2301 }
2302 else
2303 {
2304 if (re_val.dims () == im_val.dims ())
2305 {
2306 SparseComplexMatrix result = SparseComplexMatrix(re_val)
2307 + Complex(0, 1) * SparseComplexMatrix (im_val);
2308 retval = octave_value (new octave_sparse_complex_matrix (result));
2309 }
2310 else
2311 error ("complex: dimension mismatch");
2312 }
2313 }
2314 }
2315 else if (re.is_single_type () || im.is_single_type ())
2316 {
2317 if (re.numel () == 1)
2318 {
2319 float re_val = re.float_value ();
2320
2321 if (im.numel () == 1)
2322 {
2323 float im_val = im.double_value ();
2324
2325 if (! error_state)
2326 retval = octave_value (new octave_float_complex (FloatComplex (re_val, im_val)));
2327 }
2328 else
2329 {
2330 const FloatNDArray im_val = im.float_array_value ();
2331
2332 if (! error_state)
2333 {
2334 FloatComplexNDArray result (im_val.dims (), FloatComplex ());
2335
2336 for (octave_idx_type i = 0; i < im_val.numel (); i++)
2337 result.xelem (i) = FloatComplex (re_val, im_val(i));
2338
2339 retval = octave_value (new octave_float_complex_matrix (result));
2340 }
2341 }
2342 }
2343 else
2344 {
2345 const FloatNDArray re_val = re.float_array_value ();
2346
2347 if (im.numel () == 1)
2348 {
2349 float im_val = im.float_value ();
2350
2351 if (! error_state)
2352 {
2353 FloatComplexNDArray result (re_val.dims (), FloatComplex ());
2354
2355 for (octave_idx_type i = 0; i < re_val.numel (); i++)
2356 result.xelem (i) = FloatComplex (re_val(i), im_val);
2357
2358 retval = octave_value (new octave_float_complex_matrix (result));
2359 }
2360 }
2361 else
2362 {
2363 const FloatNDArray im_val = im.float_array_value ();
2364
2365 if (! error_state)
2366 {
2367 if (re_val.dims () == im_val.dims ())
2368 {
2369 FloatComplexNDArray result (re_val.dims (), FloatComplex ());
2370
2371 for (octave_idx_type i = 0; i < re_val.numel (); i++)
2372 result.xelem (i) = FloatComplex (re_val(i), im_val(i));
2373
2374 retval = octave_value (new octave_float_complex_matrix (result));
2375 }
2376 else
2377 error ("complex: dimension mismatch");
2378 }
2379 }
2380 }
2381 }
2382 else if (re.numel () == 1)
1878 { 2383 {
1879 double re_val = re.double_value (); 2384 double re_val = re.double_value ();
1880 2385
1881 if (im.numel () == 1) 2386 if (im.numel () == 1)
1882 { 2387 {
2132 2637
2133 case oct_data_conv::dt_uint64: 2638 case oct_data_conv::dt_uint64:
2134 retval = uint64NDArray (dims, val); 2639 retval = uint64NDArray (dims, val);
2135 break; 2640 break;
2136 2641
2137 case oct_data_conv::dt_single: // FIXME 2642 case oct_data_conv::dt_single:
2643 retval = FloatNDArray (dims, val);
2644 break;
2645
2138 case oct_data_conv::dt_double: 2646 case oct_data_conv::dt_double:
2139 retval = NDArray (dims, val); 2647 retval = NDArray (dims, val);
2140 break; 2648 break;
2141 2649
2142 case oct_data_conv::dt_logical: 2650 case oct_data_conv::dt_logical:
2213 2721
2214 if (! error_state) 2722 if (! error_state)
2215 { 2723 {
2216 switch (dt) 2724 switch (dt)
2217 { 2725 {
2218 case oct_data_conv::dt_single: // FIXME 2726 case oct_data_conv::dt_single:
2727 retval = FloatNDArray (dims, val);
2728 break;
2729
2219 case oct_data_conv::dt_double: 2730 case oct_data_conv::dt_double:
2220 retval = NDArray (dims, val); 2731 retval = NDArray (dims, val);
2221 break; 2732 break;
2222 2733
2223 default: 2734 default:
2291 2802
2292 if (! error_state) 2803 if (! error_state)
2293 { 2804 {
2294 switch (dt) 2805 switch (dt)
2295 { 2806 {
2296 case oct_data_conv::dt_single: // FIXME 2807 case oct_data_conv::dt_single:
2808 retval = FloatComplexNDArray (dims, static_cast <FloatComplex> (val));
2809 break;
2810
2297 case oct_data_conv::dt_double: 2811 case oct_data_conv::dt_double:
2298 retval = ComplexNDArray (dims, val); 2812 retval = ComplexNDArray (dims, val);
2299 break; 2813 break;
2300 2814
2301 default: 2815 default:
2690 INSTANTIATE_EYE (uint16NDArray); 3204 INSTANTIATE_EYE (uint16NDArray);
2691 INSTANTIATE_EYE (int32NDArray); 3205 INSTANTIATE_EYE (int32NDArray);
2692 INSTANTIATE_EYE (uint32NDArray); 3206 INSTANTIATE_EYE (uint32NDArray);
2693 INSTANTIATE_EYE (int64NDArray); 3207 INSTANTIATE_EYE (int64NDArray);
2694 INSTANTIATE_EYE (uint64NDArray); 3208 INSTANTIATE_EYE (uint64NDArray);
3209 INSTANTIATE_EYE (FloatNDArray);
2695 INSTANTIATE_EYE (NDArray); 3210 INSTANTIATE_EYE (NDArray);
2696 INSTANTIATE_EYE (boolNDArray); 3211 INSTANTIATE_EYE (boolNDArray);
2697 3212
2698 static octave_value 3213 static octave_value
2699 identity_matrix (int nr, int nc, oct_data_conv::data_type dt) 3214 identity_matrix (int nr, int nc, oct_data_conv::data_type dt)
2738 3253
2739 case oct_data_conv::dt_uint64: 3254 case oct_data_conv::dt_uint64:
2740 retval = identity_matrix<uint64NDArray> (nr, nc); 3255 retval = identity_matrix<uint64NDArray> (nr, nc);
2741 break; 3256 break;
2742 3257
2743 case oct_data_conv::dt_single: // FIXME 3258 case oct_data_conv::dt_single:
3259 retval = identity_matrix<FloatNDArray> (nr, nc);
3260 break;
3261
2744 case oct_data_conv::dt_double: 3262 case oct_data_conv::dt_double:
2745 retval = identity_matrix<NDArray> (nr, nc); 3263 retval = identity_matrix<NDArray> (nr, nc);
2746 break; 3264 break;
2747 3265
2748 case oct_data_conv::dt_logical: 3266 case oct_data_conv::dt_logical:
2892 if (! error_state) 3410 if (! error_state)
2893 { 3411 {
2894 octave_value arg_1 = args(0); 3412 octave_value arg_1 = args(0);
2895 octave_value arg_2 = args(1); 3413 octave_value arg_2 = args(1);
2896 3414
2897 if (arg_1.is_complex_type () || arg_2.is_complex_type ()) 3415 if (arg_1.is_single_type () || arg_2.is_single_type ())
2898 { 3416 {
2899 Complex x1 = arg_1.complex_value (); 3417 if (arg_1.is_complex_type () || arg_2.is_complex_type ())
2900 Complex x2 = arg_2.complex_value (); 3418 {
2901 3419 FloatComplex x1 = arg_1.float_complex_value ();
2902 if (! error_state) 3420 FloatComplex x2 = arg_2.float_complex_value ();
2903 {
2904 ComplexRowVector rv = linspace (x1, x2, npoints);
2905 3421
2906 if (! error_state) 3422 if (! error_state)
2907 retval = rv; 3423 {
3424 FloatComplexRowVector rv = linspace (x1, x2, npoints);
3425
3426 if (! error_state)
3427 retval = rv;
3428 }
3429 }
3430 else
3431 {
3432 float x1 = arg_1.float_value ();
3433 float x2 = arg_2.float_value ();
3434
3435 if (! error_state)
3436 {
3437 FloatRowVector rv = linspace (x1, x2, npoints);
3438
3439 if (! error_state)
3440 retval = rv;
3441 }
2908 } 3442 }
2909 } 3443 }
2910 else 3444 else
2911 { 3445 {
2912 double x1 = arg_1.double_value (); 3446 if (arg_1.is_complex_type () || arg_2.is_complex_type ())
2913 double x2 = arg_2.double_value (); 3447 {
2914 3448 Complex x1 = arg_1.complex_value ();
2915 if (! error_state) 3449 Complex x2 = arg_2.complex_value ();
2916 {
2917 RowVector rv = linspace (x1, x2, npoints);
2918 3450
2919 if (! error_state) 3451 if (! error_state)
2920 retval = rv; 3452 {
3453 ComplexRowVector rv = linspace (x1, x2, npoints);
3454
3455 if (! error_state)
3456 retval = rv;
3457 }
3458 }
3459 else
3460 {
3461 double x1 = arg_1.double_value ();
3462 double x2 = arg_2.double_value ();
3463
3464 if (! error_state)
3465 {
3466 RowVector rv = linspace (x1, x2, npoints);
3467
3468 if (! error_state)
3469 retval = rv;
3470 }
2921 } 3471 }
2922 } 3472 }
2923 } 3473 }
2924 else 3474 else
2925 error ("linspace: expecting third argument to be an integer"); 3475 error ("linspace: expecting third argument to be an integer");