Mercurial > octave

comparison liboctave/CDiagMatrix.cc @ 1205:8302fab9fe24

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1995-04-04 02:05:01 by jwe]
author jwe
date Tue, 04 Apr 1995 02:05:01 +0000
parents b6360f2d4fa6
children 0bf4d2b7def4
comparison
equal deleted inserted replaced
1204:68d147abe7ca 1205:8302fab9fe24
64 { 64 {
65 for (int i = 0; i < length (); i++) 65 for (int i = 0; i < length (); i++)
66 elem (i, i) = a.elem (i, i); 66 elem (i, i) = a.elem (i, i);
67 } 67 }
68 68
69 #if 0
70 ComplexDiagMatrix&
71 ComplexDiagMatrix::resize (int r, int c)
72 {
73 if (r < 0 || c < 0)
74 {
75 (*current_liboctave_error_handler)
76 ("can't resize to negative dimensions");
77 return *this;
78 }
79
80 int new_len = r < c ? r : c;
81 Complex *new_data = 0;
82 if (new_len > 0)
83 {
84 new_data = new Complex [new_len];
85
86 int min_len = new_len < len ? new_len : len;
87
88 for (int i = 0; i < min_len; i++)
89 new_data[i] = data[i];
90 }
91
92 delete [] data;
93 nr = r;
94 nc = c;
95 len = new_len;
96 data = new_data;
97
98 return *this;
99 }
100
101 ComplexDiagMatrix&
102 ComplexDiagMatrix::resize (int r, int c, double val)
103 {
104 if (r < 0 || c < 0)
105 {
106 (*current_liboctave_error_handler)
107 ("can't resize to negative dimensions");
108 return *this;
109 }
110
111 int new_len = r < c ? r : c;
112 Complex *new_data = 0;
113 if (new_len > 0)
114 {
115 new_data = new Complex [new_len];
116
117 int min_len = new_len < len ? new_len : len;
118
119 for (int i = 0; i < min_len; i++)
120 new_data[i] = data[i];
121
122 for (i = min_len; i < new_len; i++)
123 new_data[i] = val;
124 }
125
126 delete [] data;
127 nr = r;
128 nc = c;
129 len = new_len;
130 data = new_data;
131
132 return *this;
133 }
134
135 ComplexDiagMatrix&
136 ComplexDiagMatrix::resize (int r, int c, const Complex& val)
137 {
138 if (r < 0 || c < 0)
139 {
140 (*current_liboctave_error_handler)
141 ("can't resize to negative dimensions");
142 return *this;
143 }
144
145 int new_len = r < c ? r : c;
146 Complex *new_data = 0;
147 if (new_len > 0)
148 {
149 new_data = new Complex [new_len];
150
151 int min_len = new_len < len ? new_len : len;
152
153 for (int i = 0; i < min_len; i++)
154 new_data[i] = data[i];
155
156 for (i = min_len; i < new_len; i++)
157 new_data[i] = val;
158 }
159
160 delete [] data;
161 nr = r;
162 nc = c;
163 len = new_len;
164 data = new_data;
165
166 return *this;
167 }
168 #endif
169
170 int 69 int
171 ComplexDiagMatrix::operator == (const ComplexDiagMatrix& a) const 70 ComplexDiagMatrix::operator == (const ComplexDiagMatrix& a) const
172 { 71 {
173 if (rows () != a.rows () || cols () != a.cols ()) 72 if (rows () != a.rows () || cols () != a.cols ())
174 return 0; 73 return 0;
180 ComplexDiagMatrix::operator != (const ComplexDiagMatrix& a) const 79 ComplexDiagMatrix::operator != (const ComplexDiagMatrix& a) const
181 { 80 {
182 return !(*this == a); 81 return !(*this == a);
183 } 82 }
184 83
185 ComplexDiagMatrix
186 ComplexDiagMatrix::hermitian (void) const
187 {
188 return ComplexDiagMatrix (conj_dup (data (), length ()), cols (), rows ());
189 }
190
191 ComplexDiagMatrix& 84 ComplexDiagMatrix&
192 ComplexDiagMatrix::fill (double val) 85 ComplexDiagMatrix::fill (double val)
193 { 86 {
194 for (int i = 0; i < length (); i++) 87 for (int i = 0; i < length (); i++)
195 elem (i, i) = val; 88 elem (i, i) = val;
361 254
362 return *this; 255 return *this;
363 } 256 }
364 257
365 ComplexDiagMatrix 258 ComplexDiagMatrix
259 ComplexDiagMatrix::hermitian (void) const
260 {
261 return ComplexDiagMatrix (conj_dup (data (), length ()), cols (), rows ());
262 }
263
264 ComplexDiagMatrix
366 ComplexDiagMatrix::transpose (void) const 265 ComplexDiagMatrix::transpose (void) const
367 { 266 {
368 return ComplexDiagMatrix (dup (data (), length ()), cols (), rows ()); 267 return ComplexDiagMatrix (dup (data (), length ()), cols (), rows ());
369 }
370
371 DiagMatrix
372 real (const ComplexDiagMatrix& a)
373 {
374 DiagMatrix retval;
375 int a_len = a.length ();
376 if (a_len > 0)
377 retval = DiagMatrix (real_dup (a.data (), a_len), a.rows (),
378 a.cols ());
379 return retval;
380 }
381
382 DiagMatrix
383 imag (const ComplexDiagMatrix& a)
384 {
385 DiagMatrix retval;
386 int a_len = a.length ();
387 if (a_len > 0)
388 retval = DiagMatrix (imag_dup (a.data (), a_len), a.rows (),
389 a.cols ());
390 return retval;
391 } 268 }
392 269
393 ComplexDiagMatrix 270 ComplexDiagMatrix
394 conj (const ComplexDiagMatrix& a) 271 conj (const ComplexDiagMatrix& a)
395 { 272 {
620 497
621 subtract2 (d, a.data (), length ()); 498 subtract2 (d, a.data (), length ());
622 return *this; 499 return *this;
623 } 500 }
624 501
625 // diagonal matrix by scalar -> matrix operations
626
627 ComplexMatrix
628 operator + (const ComplexDiagMatrix& a, double s)
629 {
630 ComplexMatrix tmp (a.rows (), a.cols (), s);
631 return a + tmp;
632 }
633
634 ComplexMatrix
635 operator - (const ComplexDiagMatrix& a, double s)
636 {
637 ComplexMatrix tmp (a.rows (), a.cols (), -s);
638 return a + tmp;
639 }
640
641 ComplexMatrix
642 operator + (const ComplexDiagMatrix& a, const Complex& s)
643 {
644 ComplexMatrix tmp (a.rows (), a.cols (), s);
645 return a + tmp;
646 }
647
648 ComplexMatrix
649 operator - (const ComplexDiagMatrix& a, const Complex& s)
650 {
651 ComplexMatrix tmp (a.rows (), a.cols (), -s);
652 return a + tmp;
653 }
654
655 // diagonal matrix by scalar -> diagonal matrix operations 502 // diagonal matrix by scalar -> diagonal matrix operations
656 503
657 ComplexDiagMatrix 504 ComplexDiagMatrix
658 operator * (const ComplexDiagMatrix& a, double s) 505 operator * (const ComplexDiagMatrix& a, double s)
659 { 506 {
666 { 513 {
667 return ComplexDiagMatrix (divide (a.data (), a.length (), s), 514 return ComplexDiagMatrix (divide (a.data (), a.length (), s),
668 a.rows (), a.cols ()); 515 a.rows (), a.cols ());
669 } 516 }
670 517
671 // scalar by diagonal matrix -> matrix operations 518 ComplexDiagMatrix
672 519 operator * (const DiagMatrix& a, const Complex& s)
673 ComplexMatrix
674 operator + (double s, const ComplexDiagMatrix& a)
675 {
676 ComplexMatrix tmp (a.rows (), a.cols (), s);
677 return tmp + a;
678 }
679
680 ComplexMatrix
681 operator - (double s, const ComplexDiagMatrix& a)
682 {
683 ComplexMatrix tmp (a.rows (), a.cols (), s);
684 return tmp - a;
685 }
686
687 ComplexMatrix
688 operator + (const Complex& s, const ComplexDiagMatrix& a)
689 {
690 ComplexMatrix tmp (a.rows (), a.cols (), s);
691 return tmp + a;
692 }
693
694 ComplexMatrix
695 operator - (const Complex& s, const ComplexDiagMatrix& a)
696 {
697 ComplexMatrix tmp (a.rows (), a.cols (), s);
698 return tmp - a;
699 }
700
701 // scalar by diagonal matrix -> diagonal matrix operations
702
703 ComplexDiagMatrix
704 operator * (double s, const ComplexDiagMatrix& a)
705 { 520 {
706 return ComplexDiagMatrix (multiply (a.data (), a.length (), s), 521 return ComplexDiagMatrix (multiply (a.data (), a.length (), s),
707 a.rows (), a.cols ()); 522 a.rows (), a.cols ());
708 } 523 }
709 524
710 // diagonal matrix by column vector -> column vector operations 525 ComplexDiagMatrix
711 526 operator / (const DiagMatrix& a, const Complex& s)
712 ComplexColumnVector 527 {
713 operator * (const ComplexDiagMatrix& m, const ColumnVector& a) 528 return ComplexDiagMatrix (divide (a.data (), a.length (), s),
714 { 529 a.rows (), a.cols ());
715 int nr = m.rows (); 530 }
716 int nc = m.cols (); 531
717 int a_len = a.length (); 532 // scalar by diagonal matrix -> diagonal matrix operations
718 if (nc != a_len) 533
719 { 534 ComplexDiagMatrix
720 (*current_liboctave_error_handler) 535 operator * (double s, const ComplexDiagMatrix& a)
721 ("nonconformant matrix muliplication attempted"); 536 {
722 return ComplexColumnVector (); 537 return ComplexDiagMatrix (multiply (a.data (), a.length (), s),
723 } 538 a.rows (), a.cols ());
724 539 }
725 if (nc == 0 || nr == 0) 540
726 return ComplexColumnVector (0); 541 ComplexDiagMatrix
727 542 operator * (const Complex& s, const DiagMatrix& a)
728 ComplexColumnVector result (nr); 543 {
729 544 return ComplexDiagMatrix (multiply (a.data (), a.length (), s),
730 for (int i = 0; i < a_len; i++) 545 a.rows (), a.cols ());
731 result.elem (i) = a.elem (i) * m.elem (i, i);
732
733 for (i = a_len; i < nr; i++)
734 result.elem (i) = 0.0;
735
736 return result;
737 }
738
739 ComplexColumnVector
740 operator * (const ComplexDiagMatrix& m, const ComplexColumnVector& a)
741 {
742 int nr = m.rows ();
743 int nc = m.cols ();
744 int a_len = a.length ();
745 if (nc != a_len)
746 {
747 (*current_liboctave_error_handler)
748 ("nonconformant matrix muliplication attempted");
749 return ComplexColumnVector ();
750 }
751
752 if (nc == 0 || nr == 0)
753 return ComplexColumnVector (0);
754
755 ComplexColumnVector result (nr);
756
757 for (int i = 0; i < a_len; i++)
758 result.elem (i) = a.elem (i) * m.elem (i, i);
759
760 for (i = a_len; i < nr; i++)
761 result.elem (i) = 0.0;
762
763 return result;
764 } 546 }
765 547
766 // diagonal matrix by diagonal matrix -> diagonal matrix operations 548 // diagonal matrix by diagonal matrix -> diagonal matrix operations
767 549
768 ComplexDiagMatrix 550 ComplexDiagMatrix
879 661
880 return c; 662 return c;
881 } 663 }
882 664
883 ComplexDiagMatrix 665 ComplexDiagMatrix
884 product (const ComplexDiagMatrix& m, const DiagMatrix& a) 666 operator + (const DiagMatrix& m, const ComplexDiagMatrix& a)
885 { 667 {
886 int nr = m.rows (); 668 int nr = m.rows ();
887 int nc = m.cols (); 669 int nc = m.cols ();
888 if (nr != a.rows () || nc != a.cols ()) 670 if (nr != a.rows () || nc != a.cols ())
889 { 671 {
890 (*current_liboctave_error_handler) 672 (*current_liboctave_error_handler)
673 ("nonconformant matrix addition attempted");
674 return ComplexDiagMatrix ();
675 }
676
677 if (nc == 0 || nr == 0)
678 return ComplexDiagMatrix (nr, nc);
679
680 return ComplexDiagMatrix (add (m.data (), a.data (), m.length ()), nr, nc);
681 }
682
683 ComplexDiagMatrix
684 operator - (const DiagMatrix& m, const ComplexDiagMatrix& a)
685 {
686 int nr = m.rows ();
687 int nc = m.cols ();
688 if (nr != a.rows () || nc != a.cols ())
689 {
690 (*current_liboctave_error_handler)
691 ("nonconformant matrix subtraction attempted");
692 return ComplexDiagMatrix ();
693 }
694
695 if (nc == 0 || nr == 0)
696 return ComplexDiagMatrix (nr, nc);
697
698 return ComplexDiagMatrix (subtract (m.data (), a.data (), m.length ()),
699 nr, nc);
700 }
701
702 ComplexDiagMatrix
703 operator * (const DiagMatrix& a, const ComplexDiagMatrix& b)
704 {
705 int nr_a = a.rows ();
706 int nc_a = a.cols ();
707 int nr_b = b.rows ();
708 int nc_b = b.cols ();
709 if (nc_a != nr_b)
710 {
711 (*current_liboctave_error_handler)
712 ("nonconformant matrix multiplication attempted");
713 return ComplexDiagMatrix ();
714 }
715
716 if (nr_a == 0 || nc_a == 0 || nc_b == 0)
717 return ComplexDiagMatrix (nr_a, nc_a, 0.0);
718
719 ComplexDiagMatrix c (nr_a, nc_b);
720
721 int len = nr_a < nc_b ? nr_a : nc_b;
722
723 for (int i = 0; i < len; i++)
724 {
725 double a_element = a.elem (i, i);
726 Complex b_element = b.elem (i, i);
727
728 if (a_element == 0.0 || b_element == 0.0)
729 c.elem (i, i) = 0.0;
730 else if (a_element == 1.0)
731 c.elem (i, i) = b_element;
732 else if (b_element == 1.0)
733 c.elem (i, i) = a_element;
734 else
735 c.elem (i, i) = a_element * b_element;
736 }
737
738 return c;
739 }
740
741 ComplexDiagMatrix
742 product (const ComplexDiagMatrix& m, const DiagMatrix& a)
743 {
744 int nr = m.rows ();
745 int nc = m.cols ();
746 if (nr != a.rows () || nc != a.cols ())
747 {
748 (*current_liboctave_error_handler)
891 ("nonconformant matrix product attempted"); 749 ("nonconformant matrix product attempted");
892 return ComplexDiagMatrix (); 750 return ComplexDiagMatrix ();
893 } 751 }
894 752
895 if (nr == 0 || nc == 0) 753 if (nr == 0 || nc == 0)
897 755
898 return ComplexDiagMatrix (multiply (m.data (), a.data (), m.length ()), 756 return ComplexDiagMatrix (multiply (m.data (), a.data (), m.length ()),
899 nr, nc); 757 nr, nc);
900 } 758 }
901 759
902 // diagonal matrix by matrix -> matrix operations 760 ComplexDiagMatrix
903 761 product (const DiagMatrix& m, const ComplexDiagMatrix& a)
904 ComplexMatrix
905 operator + (const ComplexDiagMatrix& m, const Matrix& a)
906 { 762 {
907 int nr = m.rows (); 763 int nr = m.rows ();
908 int nc = m.cols (); 764 int nc = m.cols ();
909 if (nr != a.rows () || nc != a.cols ()) 765 if (nr != a.rows () || nc != a.cols ())
910 { 766 {
911 (*current_liboctave_error_handler) 767 (*current_liboctave_error_handler)
912 ("nonconformant matrix addition attempted"); 768 ("nonconformant matrix product attempted");
913 return ComplexMatrix (); 769 return ComplexDiagMatrix ();
914 } 770 }
915 771
916 if (nr == 0 || nc == 0) 772 if (nc == 0 || nr == 0)
917 return ComplexMatrix (nr, nc); 773 return ComplexDiagMatrix (nr, nc);
918 774
919 ComplexMatrix result (a); 775 return ComplexDiagMatrix (multiply (m.data (), a.data (), m.length ()),
920 for (int i = 0; i < m.length (); i++) 776 nr, nc);
921 result.elem (i, i) += m.elem (i, i);
922
923 return result;
924 }
925
926 ComplexMatrix
927 operator - (const ComplexDiagMatrix& m, const Matrix& a)
928 {
929 int nr = m.rows ();
930 int nc = m.cols ();
931 if (nr != a.rows () || nc != a.cols ())
932 {
933 (*current_liboctave_error_handler)
934 ("nonconformant matrix subtraction attempted");
935 return ComplexMatrix ();
936 }
937
938 if (nr == 0 || nc == 0)
939 return ComplexMatrix (nr, nc);
940
941 ComplexMatrix result (-a);
942 for (int i = 0; i < m.length (); i++)
943 result.elem (i, i) += m.elem (i, i);
944
945 return result;
946 }
947
948 ComplexMatrix
949 operator * (const ComplexDiagMatrix& m, const Matrix& a)
950 {
951 int nr = m.rows ();
952 int nc = m.cols ();
953 int a_nr = a.rows ();
954 int a_nc = a.cols ();
955 if (nc != a_nr)
956 {
957 (*current_liboctave_error_handler)
958 ("nonconformant matrix multiplication attempted");
959 return ComplexMatrix ();
960 }
961
962 if (nr == 0 || nc == 0 || a_nc == 0)
963 return ComplexMatrix (nr, a_nc, 0.0);
964
965 ComplexMatrix c (nr, a_nc);
966
967 for (int i = 0; i < m.length (); i++)
968 {
969 if (m.elem (i, i) == 1.0)
970 {
971 for (int j = 0; j < a_nc; j++)
972 c.elem (i, j) = a.elem (i, j);
973 }
974 else if (m.elem (i, i) == 0.0)
975 {
976 for (int j = 0; j < a_nc; j++)
977 c.elem (i, j) = 0.0;
978 }
979 else
980 {
981 for (int j = 0; j < a_nc; j++)
982 c.elem (i, j) = m.elem (i, i) * a.elem (i, j);
983 }
984 }
985
986 if (nr > nc)
987 {
988 for (int j = 0; j < a_nc; j++)
989 for (int i = a_nr; i < nr; i++)
990 c.elem (i, j) = 0.0;
991 }
992
993 return c;
994 }
995
996 ComplexMatrix
997 operator + (const ComplexDiagMatrix& m, const ComplexMatrix& a)
998 {
999 int nr = m.rows ();
1000 int nc = m.cols ();
1001 if (nr != a.rows () || nc != a.cols ())
1002 {
1003 (*current_liboctave_error_handler)
1004 ("nonconformant matrix addition attempted");
1005 return ComplexMatrix ();
1006 }
1007
1008 if (nr == 0 || nc == 0)
1009 return ComplexMatrix (nr, nc);
1010
1011 ComplexMatrix result (a);
1012 for (int i = 0; i < m.length (); i++)
1013 result.elem (i, i) += m.elem (i, i);
1014
1015 return result;
1016 }
1017
1018 ComplexMatrix
1019 operator - (const ComplexDiagMatrix& m, const ComplexMatrix& a)
1020 {
1021 int nr = m.rows ();
1022 int nc = m.cols ();
1023 if (nr != a.rows () || nc != a.cols ())
1024 {
1025 (*current_liboctave_error_handler)
1026 ("nonconformant matrix subtraction attempted");
1027 return ComplexMatrix ();
1028 }
1029
1030 if (nr == 0 || nc == 0)
1031 return ComplexMatrix (nr, nc);
1032
1033 ComplexMatrix result (-a);
1034 for (int i = 0; i < m.length (); i++)
1035 result.elem (i, i) += m.elem (i, i);
1036
1037 return result;
1038 }
1039
1040 ComplexMatrix
1041 operator * (const ComplexDiagMatrix& m, const ComplexMatrix& a)
1042 {
1043 int nr = m.rows ();
1044 int nc = m.cols ();
1045 int a_nr = a.rows ();
1046 int a_nc = a.cols ();
1047 if (nc != a_nr)
1048 {
1049 (*current_liboctave_error_handler)
1050 ("nonconformant matrix multiplication attempted");
1051 return ComplexMatrix ();
1052 }
1053
1054 if (nr == 0 || nc == 0 || a_nc == 0)
1055 return ComplexMatrix (nr, a_nc, 0.0);
1056
1057 ComplexMatrix c (nr, a_nc);
1058
1059 for (int i = 0; i < m.length (); i++)
1060 {
1061 if (m.elem (i, i) == 1.0)
1062 {
1063 for (int j = 0; j < a_nc; j++)
1064 c.elem (i, j) = a.elem (i, j);
1065 }
1066 else if (m.elem (i, i) == 0.0)
1067 {
1068 for (int j = 0; j < a_nc; j++)
1069 c.elem (i, j) = 0.0;
1070 }
1071 else
1072 {
1073 for (int j = 0; j < a_nc; j++)
1074 c.elem (i, j) = m.elem (i, i) * a.elem (i, j);
1075 }
1076 }
1077
1078 if (nr > nc)
1079 {
1080 for (int j = 0; j < a_nc; j++)
1081 for (int i = a_nr; i < nr; i++)
1082 c.elem (i, j) = 0.0;
1083 }
1084
1085 return c;
1086 } 777 }
1087 778
1088 // other operations 779 // other operations
1089 780
1090 ComplexColumnVector 781 ComplexColumnVector