Mercurial > octave

comparison liboctave/MatrixType.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 a1dbe9d80eee
children 42c42c640108
comparison
equal deleted inserted replaced
7788:45f5faba05a2 7789:82be108cc558
112 else 112 else
113 typ = MatrixType::Rectangular; 113 typ = MatrixType::Rectangular;
114 } 114 }
115 115
116 MatrixType::MatrixType (const ComplexMatrix &a) 116 MatrixType::MatrixType (const ComplexMatrix &a)
117 : typ (MatrixType::Unknown),
118 sp_bandden (0), bandden (0), upper_band (0), lower_band (0),
119 dense (false), full (true), nperm (0), perm (0)
120 {
121 octave_idx_type nrows = a.rows ();
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 }
177
178
179 MatrixType::MatrixType (const FloatMatrix &a)
180 : typ (MatrixType::Unknown),
181 sp_bandden (0), bandden (0), upper_band (0), lower_band (0),
182 dense (false), full (true), nperm (0), perm (0)
183 {
184 octave_idx_type nrows = a.rows ();
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 }
236
237 MatrixType::MatrixType (const FloatComplexMatrix &a)
117 : typ (MatrixType::Unknown), 238 : typ (MatrixType::Unknown),
118 sp_bandden (0), bandden (0), upper_band (0), lower_band (0), 239 sp_bandden (0), bandden (0), upper_band (0), lower_band (0),
119 dense (false), full (true), nperm (0), perm (0) 240 dense (false), full (true), nperm (0), perm (0)
120 { 241 {
121 octave_idx_type nrows = a.rows (); 242 octave_idx_type nrows = a.rows ();
998 perm[i] = tmp_typ.perm[i]; 1119 perm[i] = tmp_typ.perm[i];
999 } 1120 }
1000 1121
1001 return typ; 1122 return typ;
1002 } 1123 }
1124
1003 int 1125 int
1004 MatrixType::type (const Matrix &a) 1126 MatrixType::type (const Matrix &a)
1005 { 1127 {
1006 if (typ != MatrixType::Unknown) 1128 if (typ != MatrixType::Unknown)
1007 { 1129 {
1027 return typ; 1149 return typ;
1028 } 1150 }
1029 1151
1030 int 1152 int
1031 MatrixType::type (const ComplexMatrix &a) 1153 MatrixType::type (const ComplexMatrix &a)
1154 {
1155 if (typ != MatrixType::Unknown)
1156 {
1157 if (octave_sparse_params::get_key ("spumoni") != 0.)
1158 (*current_liboctave_warning_handler)
1159 ("Using Cached Matrix Type");
1160
1161 return typ;
1162 }
1163
1164 MatrixType tmp_typ (a);
1165 typ = tmp_typ.typ;
1166 full = tmp_typ.full;
1167 nperm = tmp_typ.nperm;
1168
1169 if (nperm != 0)
1170 {
1171 perm = new octave_idx_type [nperm];
1172 for (octave_idx_type i = 0; i < nperm; i++)
1173 perm[i] = tmp_typ.perm[i];
1174 }
1175
1176 return typ;
1177 }
1178
1179 int
1180 MatrixType::type (const FloatMatrix &a)
1181 {
1182 if (typ != MatrixType::Unknown)
1183 {
1184 if (octave_sparse_params::get_key ("spumoni") != 0.)
1185 (*current_liboctave_warning_handler)
1186 ("Using Cached Matrix Type");
1187
1188 return typ;
1189 }
1190
1191 MatrixType tmp_typ (a);
1192 typ = tmp_typ.typ;
1193 full = tmp_typ.full;
1194 nperm = tmp_typ.nperm;
1195
1196 if (nperm != 0)
1197 {
1198 perm = new octave_idx_type [nperm];
1199 for (octave_idx_type i = 0; i < nperm; i++)
1200 perm[i] = tmp_typ.perm[i];
1201 }
1202
1203 return typ;
1204 }
1205
1206 int
1207 MatrixType::type (const FloatComplexMatrix &a)
1032 { 1208 {
1033 if (typ != MatrixType::Unknown) 1209 if (typ != MatrixType::Unknown)
1034 { 1210 {
1035 if (octave_sparse_params::get_key ("spumoni") != 0.) 1211 if (octave_sparse_params::get_key ("spumoni") != 0.)
1036 (*current_liboctave_warning_handler) 1212 (*current_liboctave_warning_handler)