Mercurial > hg > octave-nkf

comparison liboctave/MatrixType.cc @ 7798:42c42c640108

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
improved matrix_type check
author Jaroslav Hajek <highegg@gmail.com>
date Wed, 21 May 2008 17:14:08 +0200
parents 82be108cc558
children 25bc2d31e1bf
comparison
equal deleted inserted replaced
7797:f42c6f8d6d8e 7798:42c42c640108
53 for (octave_idx_type i = 0; i < nperm; i++) 53 for (octave_idx_type i = 0; i < nperm; i++)
54 perm[i] = a.perm[i]; 54 perm[i] = a.perm[i];
55 } 55 }
56 } 56 }
57 57
58 MatrixType::MatrixType (const Matrix &a) 58 template<class T>
59 : typ (MatrixType::Unknown), 59 MatrixType::matrix_type
60 sp_bandden (0), bandden (0), upper_band (0), lower_band (0), 60 matrix_real_probe (const MArray2<T>& a)
61 dense (false), full (true), nperm (0), perm (0) 61 {
62 { 62 MatrixType::matrix_type typ;
63 octave_idx_type nrows = a.rows (); 63 octave_idx_type nrows = a.rows ();
64 octave_idx_type ncols = a.cols (); 64 octave_idx_type ncols = a.cols ();
65
66 const T zero = 0;
65 67
66 if (ncols == nrows) 68 if (ncols == nrows)
67 { 69 {
68 bool upper = true; 70 bool upper = true;
69 bool lower = true; 71 bool lower = true;
70 bool hermitian = true; 72 bool hermitian = true;
71 73
72 for (octave_idx_type j = 0; j < ncols; j++) 74 // do the checks for lower/upper/hermitian all in one pass.
73 { 75 ColumnVector diag(ncols);
74 if (j < nrows) 76
75 { 77 for (octave_idx_type j = 0;
76 if (a.elem (j,j) == 0.) 78 j < ncols && upper; j++)
77 { 79 {
78 upper = false; 80 T d = a.elem (j,j);
79 lower = false; 81 upper = upper && (d != zero);
80 hermitian = false; 82 lower = lower && (d != zero);
81 break; 83 hermitian = hermitian && (d > zero);
82 } 84 diag(j) = d;
83 if (a.elem (j,j) < 0.) 85 }
84 hermitian = false; 86
85 } 87 for (octave_idx_type j = 0;
86 for (octave_idx_type i = 0; i < j; i++) 88 j < ncols && (upper || lower || hermitian); j++)
87 if (lower && a.elem (i,j) != 0.) 89 {
88 { 90 for (octave_idx_type i = 0; i < j; i++)
89 lower = false; 91 {
90 break; 92 double aij = a.elem (i,j), aji = a.elem (j,i);
91 } 93 lower = lower && (aij == zero);
92 for (octave_idx_type i = j+1; i < nrows; i++) 94 upper = upper && (aji == zero);
93 { 95 hermitian = hermitian && (aij == aji
94 if (hermitian && a.elem (i, j) != a.elem (j, i)) 96 && aij*aij < diag(i)*diag(j));
95 hermitian = false; 97 }
96 if (upper && a.elem (i,j) != 0) 98 }
97 upper = false; 99
98 } 100 if (upper)
99 if (!upper && !lower && !hermitian) 101 typ = MatrixType::Upper;
100 break; 102 else if (lower)
101 } 103 typ = MatrixType::Lower;
104 else if (hermitian)
105 typ = MatrixType::Hermitian;
106 else
107 typ = MatrixType::Full;
108 }
109 else
110 typ = MatrixType::Rectangular;
111
112 return typ;
113 }
114
115 template<class T>
116 MatrixType::matrix_type
117 matrix_complex_probe (const MArray2<T>& a)
118 {
119 MatrixType::matrix_type typ;
120 octave_idx_type nrows = a.rows ();
121 octave_idx_type ncols = a.cols ();
122
123 const typename T::value_type zero = 0;
124
125 if (ncols == nrows)
126 {
127 bool upper = true;
128 bool lower = true;
129 bool hermitian = true;
130
131 // do the checks for lower/upper/hermitian all in one pass.
132 ColumnVector diag(ncols);
133
134 for (octave_idx_type j = 0;
135 j < ncols && upper; j++)
136 {
137 T d = a.elem (j,j);
138 upper = upper && (d != zero);
139 lower = lower && (d != zero);
140 hermitian = hermitian && (d.real() > zero && d.imag() == zero);
141 diag (j) = d.real();
142 }
143
144 for (octave_idx_type j = 0;
145 j < ncols && (upper || lower || hermitian); j++)
146 {
147 for (octave_idx_type i = 0; i < j; i++)
148 {
149 T aij = a.elem (i,j), aji = a.elem (j,i);
150 lower = lower && (aij == zero);
151 upper = upper && (aji == zero);
152 hermitian = hermitian && (aij == std::conj (aji)
153 && std::norm (aij) < diag(i)*diag(j));
154 }
155 }
156
102 157
103 if (upper) 158 if (upper)
104 typ = MatrixType::Upper; 159 typ = MatrixType::Upper;
105 else if (lower) 160 else if (lower)
106 typ = MatrixType::Lower; 161 typ = MatrixType::Lower;
109 else if (ncols == nrows) 164 else if (ncols == nrows)
110 typ = MatrixType::Full; 165 typ = MatrixType::Full;
111 } 166 }
112 else 167 else
113 typ = MatrixType::Rectangular; 168 typ = MatrixType::Rectangular;
169
170 return typ;
171 }
172
173 MatrixType::MatrixType (const Matrix &a)
174 : typ (MatrixType::Unknown),
175 sp_bandden (0), bandden (0), upper_band (0), lower_band (0),
176 dense (false), full (true), nperm (0), perm (0)
177 {
178 typ = matrix_real_probe (a);
114 } 179 }
115 180
116 MatrixType::MatrixType (const ComplexMatrix &a) 181 MatrixType::MatrixType (const ComplexMatrix &a)
117 : typ (MatrixType::Unknown), 182 : typ (MatrixType::Unknown),
118 sp_bandden (0), bandden (0), upper_band (0), lower_band (0), 183 sp_bandden (0), bandden (0), upper_band (0), lower_band (0),
119 dense (false), full (true), nperm (0), perm (0) 184 dense (false), full (true), nperm (0), perm (0)
120 { 185 {
121 octave_idx_type nrows = a.rows (); 186 typ = matrix_complex_probe (a);
122 octave_idx_type ncols = a.cols ();
123
124 if (ncols == nrows)
125 {
126 bool upper = true;
127 bool lower = true;
128 bool hermitian = true;
129
130 for (octave_idx_type j = 0; j < ncols; j++)
131 {
132 if (j < ncols)
133 {
134 if (imag(a.elem (j,j)) == 0. &&
135 real(a.elem (j,j)) == 0.)
136 {
137 upper = false;
138 lower = false;
139 hermitian = false;
140 break;
141 }
142
143 if (imag(a.elem (j,j)) != 0. ||
144 real(a.elem (j,j)) < 0.)
145 hermitian = false;
146 }
147 for (octave_idx_type i = 0; i < j; i++)
148 if (lower && (real(a.elem (i,j)) != 0 || imag(a.elem (i,j)) != 0))
149 {
150 lower = false;
151 break;
152 }
153 for (octave_idx_type i = j+1; i < nrows; i++)
154 {
155 if (hermitian && a.elem (i, j) != conj(a.elem (j, i)))
156 hermitian = false;
157 if (upper && (real(a.elem (i,j)) != 0 ||
158 imag(a.elem (i,j)) != 0))
159 upper = false;
160 }
161 if (!upper && !lower && !hermitian)
162 break;
163 }
164
165 if (upper)
166 typ = MatrixType::Upper;
167 else if (lower)
168 typ = MatrixType::Lower;
169 else if (hermitian)
170 typ = MatrixType::Hermitian;
171 else if (ncols == nrows)
172 typ = MatrixType::Full;
173 }
174 else
175 typ = MatrixType::Rectangular;
176 } 187 }
177 188
178 189
179 MatrixType::MatrixType (const FloatMatrix &a) 190 MatrixType::MatrixType (const FloatMatrix &a)
180 : typ (MatrixType::Unknown), 191 : typ (MatrixType::Unknown),
181 sp_bandden (0), bandden (0), upper_band (0), lower_band (0), 192 sp_bandden (0), bandden (0), upper_band (0), lower_band (0),
182 dense (false), full (true), nperm (0), perm (0) 193 dense (false), full (true), nperm (0), perm (0)
183 { 194 {
184 octave_idx_type nrows = a.rows (); 195 typ = matrix_real_probe (a);
185 octave_idx_type ncols = a.cols ();
186
187 if (ncols == nrows)
188 {
189 bool upper = true;
190 bool lower = true;
191 bool hermitian = true;
192
193 for (octave_idx_type j = 0; j < ncols; j++)
194 {
195 if (j < nrows)
196 {
197 if (a.elem (j,j) == 0.)
198 {
199 upper = false;
200 lower = false;
201 hermitian = false;
202 break;
203 }
204 if (a.elem (j,j) < 0.)
205 hermitian = false;
206 }
207 for (octave_idx_type i = 0; i < j; i++)
208 if (lower && a.elem (i,j) != 0.)
209 {
210 lower = false;
211 break;
212 }
213 for (octave_idx_type i = j+1; i < nrows; i++)
214 {
215 if (hermitian && a.elem (i, j) != a.elem (j, i))
216 hermitian = false;
217 if (upper && a.elem (i,j) != 0)
218 upper = false;
219 }
220 if (!upper && !lower && !hermitian)
221 break;
222 }
223
224 if (upper)
225 typ = MatrixType::Upper;
226 else if (lower)
227 typ = MatrixType::Lower;
228 else if (hermitian)
229 typ = MatrixType::Hermitian;
230 else if (ncols == nrows)
231 typ = MatrixType::Full;
232 }
233 else
234 typ = MatrixType::Rectangular;
235 } 196 }
236 197
237 MatrixType::MatrixType (const FloatComplexMatrix &a) 198 MatrixType::MatrixType (const FloatComplexMatrix &a)
238 : typ (MatrixType::Unknown), 199 : typ (MatrixType::Unknown),
239 sp_bandden (0), bandden (0), upper_band (0), lower_band (0), 200 sp_bandden (0), bandden (0), upper_band (0), lower_band (0),
240 dense (false), full (true), nperm (0), perm (0) 201 dense (false), full (true), nperm (0), perm (0)
241 { 202 {
242 octave_idx_type nrows = a.rows (); 203 typ = matrix_complex_probe (a);
243 octave_idx_type ncols = a.cols ();
244
245 if (ncols == nrows)
246 {
247 bool upper = true;
248 bool lower = true;
249 bool hermitian = true;
250
251 for (octave_idx_type j = 0; j < ncols; j++)
252 {
253 if (j < ncols)
254 {
255 if (imag(a.elem (j,j)) == 0. &&
256 real(a.elem (j,j)) == 0.)
257 {
258 upper = false;
259 lower = false;
260 hermitian = false;
261 break;
262 }
263
264 if (imag(a.elem (j,j)) != 0. ||
265 real(a.elem (j,j)) < 0.)
266 hermitian = false;
267 }
268 for (octave_idx_type i = 0; i < j; i++)
269 if (lower && (real(a.elem (i,j)) != 0 || imag(a.elem (i,j)) != 0))
270 {
271 lower = false;
272 break;
273 }
274 for (octave_idx_type i = j+1; i < nrows; i++)
275 {
276 if (hermitian && a.elem (i, j) != conj(a.elem (j, i)))
277 hermitian = false;
278 if (upper && (real(a.elem (i,j)) != 0 ||
279 imag(a.elem (i,j)) != 0))
280 upper = false;
281 }
282 if (!upper && !lower && !hermitian)
283 break;
284 }
285
286 if (upper)
287 typ = MatrixType::Upper;
288 else if (lower)
289 typ = MatrixType::Lower;
290 else if (hermitian)
291 typ = MatrixType::Hermitian;
292 else if (ncols == nrows)
293 typ = MatrixType::Full;
294 }
295 else
296 typ = MatrixType::Rectangular;
297 } 204 }
298 205
299 MatrixType::MatrixType (const SparseMatrix &a) 206 MatrixType::MatrixType (const SparseMatrix &a)
300 : typ (MatrixType::Unknown), 207 : typ (MatrixType::Unknown),
301 sp_bandden (0), bandden (0), upper_band (0), lower_band (0), 208 sp_bandden (0), bandden (0), upper_band (0), lower_band (0),
558 465
559 if (maybe_hermitian && (typ == MatrixType::Full || 466 if (maybe_hermitian && (typ == MatrixType::Full ||
560 typ == MatrixType::Tridiagonal || 467 typ == MatrixType::Tridiagonal ||
561 typ == MatrixType::Banded)) 468 typ == MatrixType::Banded))
562 { 469 {
563 // Check for symmetry, with positive real diagonal, which
564 // has a very good chance of being symmetric positive
565 // definite..
566 bool is_herm = true; 470 bool is_herm = true;
567 471
568 for (octave_idx_type j = 0; j < ncols; j++) 472 // first, check whether the diagonal is positive & extract it
569 { 473 ColumnVector diag (ncols);
570 bool diag_positive = false; 474
571 475 for (octave_idx_type j = 0; is_herm && j < ncols; j++)
572 for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++) 476 {
573 { 477 is_herm = false;
574 octave_idx_type ri = a.ridx(i); 478 for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++)
575 479 {
576 if (ri == j) 480 if (a.ridx(i) == j)
577 { 481 {
578 if (a.data(i) == std::abs(a.data(i))) 482 double d = a.data(i);
579 diag_positive = true; 483 is_herm = d > 0.;
580 else 484 diag(j) = d;
581 break; 485 break;
582 } 486 }
583 else 487 }
584 { 488 }
585 bool found = false; 489
586 490
587 for (octave_idx_type k = a.cidx(ri); k < a.cidx(ri+1); k++) 491 // next, check symmetry and 2x2 positiveness
588 { 492
589 if (a.ridx(k) == j) 493 for (octave_idx_type j = 0; is_herm && j < ncols; j++)
590 { 494 for (octave_idx_type i = a.cidx(j); is_herm && i < a.cidx(j+1); i++)
591 if (a.data(i) == a.data(k)) 495 {
592 found = true; 496 octave_idx_type k = a.ridx(i);
593 break; 497 is_herm = k == j;
594 } 498 if (is_herm)
595 } 499 continue;
596 500 double d = a.data(i);
597 if (! found) 501 if (d*d < diag(j)*diag(k))
598 { 502 {
599 is_herm = false; 503 for (octave_idx_type l = a.cidx(k); l < a.cidx(k+1); l++)
600 break; 504 {
601 } 505 if (a.ridx(l) == j)
602 } 506 {
603 } 507 is_herm = a.data(l) == d;
604 508 break;
605 if (! diag_positive || ! is_herm) 509 }
606 { 510 }
607 is_herm = false; 511 }
608 break; 512 }
609 }
610 }
611 513
612 if (is_herm) 514 if (is_herm)
613 { 515 {
614 if (typ == MatrixType::Full) 516 if (typ == MatrixType::Full)
615 typ = MatrixType::Hermitian; 517 typ = MatrixType::Hermitian;
884 786
885 if (maybe_hermitian && (typ == MatrixType::Full || 787 if (maybe_hermitian && (typ == MatrixType::Full ||
886 typ == MatrixType::Tridiagonal || 788 typ == MatrixType::Tridiagonal ||
887 typ == MatrixType::Banded)) 789 typ == MatrixType::Banded))
888 { 790 {
889 // Check for symmetry, with positive real diagonal, which
890 // has a very good chance of being symmetric positive
891 // definite..
892 bool is_herm = true; 791 bool is_herm = true;
893 792
894 for (octave_idx_type j = 0; j < ncols; j++) 793 // first, check whether the diagonal is positive & extract it
895 { 794 ColumnVector diag (ncols);
896 bool diag_positive = false; 795
897 796 for (octave_idx_type j = 0; is_herm && j < ncols; j++)
898 for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++) 797 {
899 { 798 is_herm = false;
900 octave_idx_type ri = a.ridx(i); 799 for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++)
901 800 {
902 if (ri == j) 801 if (a.ridx(i) == j)
903 { 802 {
904 if (a.data(i) == std::abs(a.data(i))) 803 Complex d = a.data(i);
905 diag_positive = true; 804 is_herm = d.real() > 0. && d.imag() == 0.;
906 else 805 diag(j) = d.real();
907 break; 806 break;
908 } 807 }
909 else 808 }
910 { 809 }
911 bool found = false; 810
912 811 // next, check symmetry and 2x2 positiveness
913 for (octave_idx_type k = a.cidx(ri); k < a.cidx(ri+1); k++) 812
914 { 813 for (octave_idx_type j = 0; is_herm && j < ncols; j++)
915 if (a.ridx(k) == j) 814 for (octave_idx_type i = a.cidx(j); is_herm && i < a.cidx(j+1); i++)
916 { 815 {
917 if (a.data(i) == conj(a.data(k))) 816 octave_idx_type k = a.ridx(i);
918 found = true; 817 is_herm = k == j;
919 break; 818 if (is_herm)
920 } 819 continue;
921 } 820 Complex d = a.data(i);
922 821 if (std::norm (d) < diag(j)*diag(k))
923 if (! found) 822 {
924 { 823 d = std::conj (d);
925 is_herm = false; 824 for (octave_idx_type l = a.cidx(k); l < a.cidx(k+1); l++)
926 break; 825 {
927 } 826 if (a.ridx(l) == j)
928 } 827 {
929 } 828 is_herm = a.data(l) == d;
930 829 break;
931 if (! diag_positive || ! is_herm) 830 }
932 { 831 }
933 is_herm = false; 832 }
934 break; 833 }
935 } 834
936 }
937 835
938 if (is_herm) 836 if (is_herm)
939 { 837 {
940 if (typ == MatrixType::Full) 838 if (typ == MatrixType::Full)
941 typ = MatrixType::Hermitian; 839 typ = MatrixType::Hermitian;
951 : typ (MatrixType::Unknown), 849 : typ (MatrixType::Unknown),
952 sp_bandden (octave_sparse_params::get_bandden()), 850 sp_bandden (octave_sparse_params::get_bandden()),
953 bandden (0), upper_band (0), lower_band (0), 851 bandden (0), upper_band (0), lower_band (0),
954 dense (false), full (_full), nperm (0), perm (0) 852 dense (false), full (_full), nperm (0), perm (0)
955 { 853 {
956 if (t == MatrixType::Full || t == MatrixType::Diagonal || 854 if (t == MatrixType::Unknown || t == MatrixType::Full
957 t == MatrixType::Permuted_Diagonal || t == MatrixType::Upper || 855 || t == MatrixType::Diagonal || t == MatrixType::Permuted_Diagonal
958 t == MatrixType::Lower || t == MatrixType::Tridiagonal || 856 || t == MatrixType::Upper || t == MatrixType::Lower
959 t == MatrixType::Tridiagonal_Hermitian || t == MatrixType::Rectangular) 857 || t == MatrixType::Tridiagonal || t == MatrixType::Tridiagonal_Hermitian
858 || t == MatrixType::Rectangular)
960 typ = t; 859 typ = t;
961 else 860 else
962 (*current_liboctave_warning_handler) ("Invalid matrix type"); 861 (*current_liboctave_warning_handler) ("Invalid matrix type");
963 } 862 }
964 863