|
5785
|
1 /* |
|
|
2 |
|
|
3 Copyright (C) 2006 David Bateman |
|
|
4 Copyright (C) 2006 Andy Adler |
|
|
5 |
|
|
6 Octave is free software; you can redistribute it and/or modify it |
|
|
7 under the terms of the GNU General Public License as published by the |
|
|
8 Free Software Foundation; either version 2, or (at your option) any |
|
|
9 later version. |
|
|
10 |
|
|
11 Octave is distributed in the hope that it will be useful, but WITHOUT |
|
|
12 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
|
13 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
|
14 for more details. |
|
|
15 |
|
|
16 You should have received a copy of the GNU General Public License |
|
|
17 along with this program; see the file COPYING. If not, write to the |
|
|
18 Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, |
|
|
19 Boston, MA 02110-1301, USA. |
|
|
20 |
|
|
21 */ |
|
|
22 |
|
|
23 #ifdef HAVE_CONFIG_H |
|
|
24 #include <config.h> |
|
|
25 #endif |
|
|
26 |
|
|
27 #include <vector> |
|
|
28 |
|
|
29 #include "MatrixType.h" |
|
|
30 #include "dMatrix.h" |
|
|
31 #include "CMatrix.h" |
|
|
32 #include "dSparse.h" |
|
|
33 #include "CSparse.h" |
|
|
34 #include "oct-spparms.h" |
|
|
35 |
|
|
36 // FIXME There is a large code duplication here |
|
|
37 |
|
5892
|
38 MatrixType::MatrixType (void) |
|
6452
|
39 : typ (MatrixType::Unknown), sp_bandden (0), bandden (0), upper_band (0), |
|
|
40 lower_band (0), dense (false), full (false), nperm (0), perm (0) { } |
|
5785
|
41 |
|
5892
|
42 MatrixType::MatrixType (const MatrixType &a) |
|
|
43 : typ (a.typ), sp_bandden (a.sp_bandden), bandden (a.bandden), |
|
5785
|
44 upper_band (a.upper_band), lower_band (a.lower_band), |
|
|
45 dense (a.dense), full (a.full), nperm (a.nperm) |
|
|
46 { |
|
|
47 if (nperm != 0) |
|
|
48 { |
|
|
49 perm = new octave_idx_type [nperm]; |
|
|
50 for (octave_idx_type i = 0; i < nperm; i++) |
|
|
51 perm[i] = a.perm[i]; |
|
|
52 } |
|
|
53 } |
|
|
54 |
|
|
55 MatrixType::MatrixType (const Matrix &a) |
|
5892
|
56 : typ (MatrixType::Unknown), |
|
|
57 sp_bandden (0), bandden (0), upper_band (0), lower_band (0), |
|
|
58 dense (false), full (true), nperm (0), perm (0) |
|
5785
|
59 { |
|
|
60 octave_idx_type nrows = a.rows (); |
|
|
61 octave_idx_type ncols = a.cols (); |
|
5997
|
62 |
|
5785
|
63 if (ncols == nrows) |
|
|
64 { |
|
|
65 bool upper = true; |
|
|
66 bool lower = true; |
|
|
67 bool hermitian = true; |
|
|
68 |
|
|
69 for (octave_idx_type j = 0; j < ncols; j++) |
|
|
70 { |
|
|
71 if (j < nrows) |
|
|
72 { |
|
|
73 if (a.elem (j,j) == 0.) |
|
|
74 { |
|
|
75 upper = false; |
|
|
76 lower = false; |
|
|
77 hermitian = false; |
|
|
78 break; |
|
|
79 } |
|
|
80 if (a.elem (j,j) < 0.) |
|
|
81 hermitian = false; |
|
|
82 } |
|
|
83 for (octave_idx_type i = 0; i < j; i++) |
|
|
84 if (lower && a.elem (i,j) != 0.) |
|
|
85 { |
|
|
86 lower = false; |
|
|
87 break; |
|
|
88 } |
|
|
89 for (octave_idx_type i = j+1; i < nrows; i++) |
|
|
90 { |
|
|
91 if (hermitian && a.elem (i, j) != a.elem (j, i)) |
|
|
92 hermitian = false; |
|
|
93 if (upper && a.elem (i,j) != 0) |
|
|
94 upper = false; |
|
|
95 } |
|
|
96 if (!upper && !lower && !hermitian) |
|
|
97 break; |
|
|
98 } |
|
|
99 |
|
|
100 if (upper) |
|
|
101 typ = MatrixType::Upper; |
|
|
102 else if (lower) |
|
|
103 typ = MatrixType::Lower; |
|
|
104 else if (hermitian) |
|
|
105 typ = MatrixType::Hermitian; |
|
|
106 else if (ncols == nrows) |
|
|
107 typ = MatrixType::Full; |
|
|
108 } |
|
|
109 else |
|
|
110 typ = MatrixType::Rectangular; |
|
|
111 } |
|
|
112 |
|
|
113 MatrixType::MatrixType (const ComplexMatrix &a) |
|
5892
|
114 : typ (MatrixType::Unknown), |
|
|
115 sp_bandden (0), bandden (0), upper_band (0), lower_band (0), |
|
|
116 dense (false), full (true), nperm (0), perm (0) |
|
5785
|
117 { |
|
|
118 octave_idx_type nrows = a.rows (); |
|
|
119 octave_idx_type ncols = a.cols (); |
|
|
120 |
|
|
121 if (ncols == nrows) |
|
|
122 { |
|
|
123 bool upper = true; |
|
|
124 bool lower = true; |
|
|
125 bool hermitian = true; |
|
|
126 |
|
|
127 for (octave_idx_type j = 0; j < ncols; j++) |
|
|
128 { |
|
|
129 if (j < ncols) |
|
|
130 { |
|
|
131 if (imag(a.elem (j,j)) == 0. && |
|
|
132 real(a.elem (j,j)) == 0.) |
|
|
133 { |
|
|
134 upper = false; |
|
|
135 lower = false; |
|
|
136 hermitian = false; |
|
|
137 break; |
|
|
138 } |
|
|
139 |
|
|
140 if (imag(a.elem (j,j)) != 0. || |
|
|
141 real(a.elem (j,j)) < 0.) |
|
|
142 hermitian = false; |
|
|
143 } |
|
|
144 for (octave_idx_type i = 0; i < j; i++) |
|
|
145 if (lower && (real(a.elem (i,j)) != 0 || imag(a.elem (i,j)) != 0)) |
|
|
146 { |
|
|
147 lower = false; |
|
|
148 break; |
|
|
149 } |
|
|
150 for (octave_idx_type i = j+1; i < nrows; i++) |
|
|
151 { |
|
|
152 if (hermitian && a.elem (i, j) != conj(a.elem (j, i))) |
|
|
153 hermitian = false; |
|
|
154 if (upper && (real(a.elem (i,j)) != 0 || |
|
|
155 imag(a.elem (i,j)) != 0)) |
|
|
156 upper = false; |
|
|
157 } |
|
|
158 if (!upper && !lower && !hermitian) |
|
|
159 break; |
|
|
160 } |
|
|
161 |
|
|
162 if (upper) |
|
|
163 typ = MatrixType::Upper; |
|
|
164 else if (lower) |
|
|
165 typ = MatrixType::Lower; |
|
|
166 else if (hermitian) |
|
|
167 typ = MatrixType::Hermitian; |
|
|
168 else if (ncols == nrows) |
|
|
169 typ = MatrixType::Full; |
|
|
170 } |
|
|
171 else |
|
|
172 typ = MatrixType::Rectangular; |
|
|
173 } |
|
|
174 |
|
|
175 MatrixType::MatrixType (const SparseMatrix &a) |
|
5892
|
176 : typ (MatrixType::Unknown), |
|
|
177 sp_bandden (0), bandden (0), upper_band (0), lower_band (0), |
|
|
178 dense (false), full (false), nperm (0), perm (0) |
|
5785
|
179 { |
|
|
180 octave_idx_type nrows = a.rows (); |
|
|
181 octave_idx_type ncols = a.cols (); |
|
|
182 octave_idx_type nm = (ncols < nrows ? ncols : nrows); |
|
|
183 octave_idx_type nnz = a.nzmax (); |
|
|
184 |
|
5893
|
185 if (octave_sparse_params::get_key ("spumoni") != 0.) |
|
5785
|
186 (*current_liboctave_warning_handler) |
|
|
187 ("Calculating Sparse Matrix Type"); |
|
|
188 |
|
5893
|
189 sp_bandden = octave_sparse_params::get_key ("bandden"); |
|
5785
|
190 bool maybe_hermitian = false; |
|
|
191 typ = MatrixType::Full; |
|
|
192 |
|
|
193 if (nnz == nm) |
|
|
194 { |
|
|
195 matrix_type tmp_typ = MatrixType::Diagonal; |
|
|
196 octave_idx_type i; |
|
|
197 // Maybe the matrix is diagonal |
|
|
198 for (i = 0; i < nm; i++) |
|
|
199 { |
|
|
200 if (a.cidx(i+1) != a.cidx(i) + 1) |
|
|
201 { |
|
|
202 tmp_typ = MatrixType::Full; |
|
|
203 break; |
|
|
204 } |
|
|
205 if (a.ridx(i) != i) |
|
|
206 { |
|
|
207 tmp_typ = MatrixType::Permuted_Diagonal; |
|
|
208 break; |
|
|
209 } |
|
|
210 } |
|
|
211 |
|
|
212 if (tmp_typ == MatrixType::Permuted_Diagonal) |
|
|
213 { |
|
|
214 std::vector<bool> found (nrows); |
|
|
215 |
|
|
216 for (octave_idx_type j = 0; j < i; j++) |
|
|
217 found [j] = true; |
|
|
218 for (octave_idx_type j = i; j < nrows; j++) |
|
|
219 found [j] = false; |
|
|
220 |
|
|
221 for (octave_idx_type j = i; j < nm; j++) |
|
|
222 { |
|
|
223 if ((a.cidx(j+1) > a.cidx(j) + 1) || |
|
|
224 ((a.cidx(j+1) == a.cidx(j) + 1) && found [a.ridx(j)])) |
|
|
225 { |
|
|
226 tmp_typ = MatrixType::Full; |
|
|
227 break; |
|
|
228 } |
|
|
229 found [a.ridx(j)] = true; |
|
|
230 } |
|
|
231 } |
|
|
232 typ = tmp_typ; |
|
|
233 } |
|
|
234 |
|
|
235 if (typ == MatrixType::Full) |
|
|
236 { |
|
|
237 // Search for banded, upper and lower triangular matrices |
|
|
238 bool singular = false; |
|
|
239 upper_band = 0; |
|
|
240 lower_band = 0; |
|
|
241 for (octave_idx_type j = 0; j < ncols; j++) |
|
|
242 { |
|
|
243 bool zero_on_diagonal = false; |
|
|
244 if (j < nrows) |
|
|
245 { |
|
|
246 zero_on_diagonal = true; |
|
|
247 for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++) |
|
|
248 if (a.ridx(i) == j) |
|
|
249 { |
|
|
250 zero_on_diagonal = false; |
|
|
251 break; |
|
|
252 } |
|
|
253 } |
|
|
254 |
|
|
255 if (zero_on_diagonal) |
|
|
256 { |
|
|
257 singular = true; |
|
|
258 break; |
|
|
259 } |
|
|
260 |
|
|
261 if (a.cidx(j+1) != a.cidx(j)) |
|
|
262 { |
|
|
263 octave_idx_type ru = a.ridx(a.cidx(j)); |
|
|
264 octave_idx_type rl = a.ridx(a.cidx(j+1)-1); |
|
|
265 |
|
|
266 if (j - ru > upper_band) |
|
|
267 upper_band = j - ru; |
|
|
268 |
|
|
269 if (rl - j > lower_band) |
|
|
270 lower_band = rl - j; |
|
|
271 } |
|
|
272 } |
|
|
273 |
|
|
274 if (!singular) |
|
|
275 { |
|
|
276 bandden = double (nnz) / |
|
|
277 (double (ncols) * (double (lower_band) + |
|
|
278 double (upper_band)) - |
|
|
279 0.5 * double (upper_band + 1) * double (upper_band) - |
|
|
280 0.5 * double (lower_band + 1) * double (lower_band)); |
|
|
281 |
|
|
282 if (nrows == ncols && sp_bandden != 1. && bandden > sp_bandden) |
|
|
283 { |
|
|
284 if (upper_band == 1 && lower_band == 1) |
|
|
285 typ = MatrixType::Tridiagonal; |
|
|
286 else |
|
|
287 typ = MatrixType::Banded; |
|
|
288 |
|
|
289 octave_idx_type nnz_in_band = |
|
|
290 (upper_band + lower_band + 1) * nrows - |
|
|
291 (1 + upper_band) * upper_band / 2 - |
|
|
292 (1 + lower_band) * lower_band / 2; |
|
|
293 if (nnz_in_band == nnz) |
|
|
294 dense = true; |
|
|
295 else |
|
|
296 dense = false; |
|
|
297 } |
|
|
298 else if (upper_band == 0) |
|
|
299 typ = MatrixType::Lower; |
|
|
300 else if (lower_band == 0) |
|
|
301 typ = MatrixType::Upper; |
|
|
302 |
|
|
303 if (upper_band == lower_band && nrows == ncols) |
|
|
304 maybe_hermitian = true; |
|
|
305 } |
|
|
306 |
|
|
307 if (typ == MatrixType::Full) |
|
|
308 { |
|
|
309 // Search for a permuted triangular matrix, and test if |
|
|
310 // permutation is singular |
|
|
311 |
|
|
312 // FIXME |
|
|
313 // Perhaps this should be based on a dmperm algorithm |
|
|
314 bool found = false; |
|
|
315 |
|
|
316 nperm = ncols; |
|
|
317 perm = new octave_idx_type [ncols]; |
|
|
318 |
|
|
319 for (octave_idx_type i = 0; i < ncols; i++) |
|
|
320 perm [i] = -1; |
|
|
321 |
|
|
322 for (octave_idx_type i = 0; i < nm; i++) |
|
|
323 { |
|
|
324 found = false; |
|
|
325 |
|
|
326 for (octave_idx_type j = 0; j < ncols; j++) |
|
|
327 { |
|
|
328 if ((a.cidx(j+1) - a.cidx(j)) > 0 && |
|
|
329 (a.ridx(a.cidx(j+1)-1) == i)) |
|
|
330 { |
|
|
331 perm [i] = j; |
|
|
332 found = true; |
|
|
333 break; |
|
|
334 } |
|
|
335 } |
|
|
336 |
|
|
337 if (!found) |
|
|
338 break; |
|
|
339 } |
|
|
340 |
|
|
341 if (found) |
|
|
342 { |
|
|
343 typ = MatrixType::Permuted_Upper; |
|
|
344 if (ncols > nrows) |
|
|
345 { |
|
|
346 octave_idx_type k = nrows; |
|
|
347 for (octave_idx_type i = 0; i < ncols; i++) |
|
|
348 if (perm [i] == -1) |
|
|
349 perm[i] = k++; |
|
|
350 } |
|
|
351 } |
|
|
352 else if (a.cidx(nm) == a.cidx(ncols)) |
|
|
353 { |
|
|
354 nperm = nrows; |
|
|
355 delete [] perm; |
|
|
356 perm = new octave_idx_type [nrows]; |
|
|
357 OCTAVE_LOCAL_BUFFER (octave_idx_type, tmp, nrows); |
|
|
358 |
|
|
359 for (octave_idx_type i = 0; i < nrows; i++) |
|
|
360 { |
|
|
361 perm [i] = -1; |
|
|
362 tmp [i] = -1; |
|
|
363 } |
|
|
364 |
|
|
365 for (octave_idx_type j = 0; j < ncols; j++) |
|
|
366 for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++) |
|
|
367 perm [a.ridx(i)] = j; |
|
|
368 |
|
|
369 found = true; |
|
|
370 for (octave_idx_type i = 0; i < nm; i++) |
|
|
371 if (perm[i] == -1) |
|
|
372 { |
|
|
373 found = false; |
|
|
374 break; |
|
|
375 } |
|
|
376 else |
|
|
377 { |
|
|
378 tmp[perm[i]] = 1; |
|
|
379 } |
|
|
380 |
|
|
381 if (found) |
|
|
382 { |
|
|
383 octave_idx_type k = ncols; |
|
|
384 for (octave_idx_type i = 0; i < nrows; i++) |
|
|
385 { |
|
|
386 if (tmp[i] == -1) |
|
|
387 { |
|
|
388 if (k < nrows) |
|
|
389 { |
|
|
390 perm[k++] = i; |
|
|
391 } |
|
|
392 else |
|
|
393 { |
|
|
394 found = false; |
|
|
395 break; |
|
|
396 } |
|
|
397 } |
|
|
398 } |
|
|
399 } |
|
|
400 |
|
|
401 if (found) |
|
|
402 typ = MatrixType::Permuted_Lower; |
|
|
403 else |
|
|
404 { |
|
|
405 delete [] perm; |
|
|
406 nperm = 0; |
|
|
407 } |
|
|
408 } |
|
|
409 else |
|
|
410 { |
|
|
411 delete [] perm; |
|
|
412 nperm = 0; |
|
|
413 } |
|
|
414 } |
|
|
415 |
|
|
416 // FIXME |
|
|
417 // Disable lower under-determined and upper over-determined problems |
|
|
418 // as being detected, and force to treat as singular. As this seems |
|
|
419 // to cause issues |
|
|
420 if (((typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower) |
|
|
421 && nrows > ncols) || |
|
|
422 ((typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper) |
|
|
423 && nrows < ncols)) |
|
|
424 { |
|
|
425 typ = MatrixType::Rectangular; |
|
|
426 if (typ == MatrixType::Permuted_Upper || |
|
|
427 typ == MatrixType::Permuted_Lower) |
|
|
428 delete [] perm; |
|
|
429 nperm = 0; |
|
|
430 } |
|
|
431 |
|
|
432 if (typ == MatrixType::Full && ncols != nrows) |
|
|
433 typ = MatrixType::Rectangular; |
|
|
434 |
|
|
435 if (maybe_hermitian && (typ == MatrixType::Full || |
|
|
436 typ == MatrixType::Tridiagonal || |
|
|
437 typ == MatrixType::Banded)) |
|
|
438 { |
|
|
439 // Check for symmetry, with positive real diagonal, which |
|
|
440 // has a very good chance of being symmetric positive |
|
|
441 // definite.. |
|
|
442 bool is_herm = true; |
|
|
443 |
|
|
444 for (octave_idx_type j = 0; j < ncols; j++) |
|
|
445 { |
|
|
446 bool diag_positive = false; |
|
|
447 |
|
|
448 for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++) |
|
|
449 { |
|
|
450 octave_idx_type ri = a.ridx(i); |
|
|
451 |
|
|
452 if (ri == j) |
|
|
453 { |
|
|
454 if (a.data(i) == std::abs(a.data(i))) |
|
|
455 diag_positive = true; |
|
|
456 else |
|
|
457 break; |
|
|
458 } |
|
|
459 else |
|
|
460 { |
|
|
461 bool found = false; |
|
|
462 |
|
|
463 for (octave_idx_type k = a.cidx(ri); k < a.cidx(ri+1); k++) |
|
|
464 { |
|
|
465 if (a.ridx(k) == j) |
|
|
466 { |
|
|
467 if (a.data(i) == a.data(k)) |
|
|
468 found = true; |
|
|
469 break; |
|
|
470 } |
|
|
471 } |
|
|
472 |
|
|
473 if (! found) |
|
|
474 { |
|
|
475 is_herm = false; |
|
|
476 break; |
|
|
477 } |
|
|
478 } |
|
|
479 } |
|
|
480 |
|
|
481 if (! diag_positive || ! is_herm) |
|
|
482 { |
|
|
483 is_herm = false; |
|
|
484 break; |
|
|
485 } |
|
|
486 } |
|
|
487 |
|
|
488 if (is_herm) |
|
|
489 { |
|
|
490 if (typ == MatrixType::Full) |
|
|
491 typ = MatrixType::Hermitian; |
|
|
492 else if (typ == MatrixType::Banded) |
|
|
493 typ = MatrixType::Banded_Hermitian; |
|
|
494 else |
|
|
495 typ = MatrixType::Tridiagonal_Hermitian; |
|
|
496 } |
|
|
497 } |
|
|
498 } |
|
|
499 } |
|
|
500 |
|
|
501 MatrixType::MatrixType (const SparseComplexMatrix &a) |
|
5892
|
502 : typ (MatrixType::Unknown), |
|
|
503 sp_bandden (0), bandden (0), upper_band (0), lower_band (0), |
|
|
504 dense (false), full (false), nperm (0), perm (0) |
|
5785
|
505 { |
|
|
506 octave_idx_type nrows = a.rows (); |
|
|
507 octave_idx_type ncols = a.cols (); |
|
|
508 octave_idx_type nm = (ncols < nrows ? ncols : nrows); |
|
|
509 octave_idx_type nnz = a.nzmax (); |
|
|
510 |
|
5996
|
511 if (octave_sparse_params::get_key ("spumoni") != 0.) |
|
5785
|
512 (*current_liboctave_warning_handler) |
|
|
513 ("Calculating Sparse Matrix Type"); |
|
|
514 |
|
5893
|
515 sp_bandden = octave_sparse_params::get_key ("bandden"); |
|
5785
|
516 bool maybe_hermitian = false; |
|
|
517 typ = MatrixType::Full; |
|
|
518 |
|
|
519 if (nnz == nm) |
|
|
520 { |
|
|
521 matrix_type tmp_typ = MatrixType::Diagonal; |
|
|
522 octave_idx_type i; |
|
|
523 // Maybe the matrix is diagonal |
|
|
524 for (i = 0; i < nm; i++) |
|
|
525 { |
|
|
526 if (a.cidx(i+1) != a.cidx(i) + 1) |
|
|
527 { |
|
|
528 tmp_typ = MatrixType::Full; |
|
|
529 break; |
|
|
530 } |
|
|
531 if (a.ridx(i) != i) |
|
|
532 { |
|
|
533 tmp_typ = MatrixType::Permuted_Diagonal; |
|
|
534 break; |
|
|
535 } |
|
|
536 } |
|
|
537 |
|
|
538 if (tmp_typ == MatrixType::Permuted_Diagonal) |
|
|
539 { |
|
|
540 std::vector<bool> found (nrows); |
|
|
541 |
|
|
542 for (octave_idx_type j = 0; j < i; j++) |
|
|
543 found [j] = true; |
|
|
544 for (octave_idx_type j = i; j < nrows; j++) |
|
|
545 found [j] = false; |
|
|
546 |
|
|
547 for (octave_idx_type j = i; j < nm; j++) |
|
|
548 { |
|
|
549 if ((a.cidx(j+1) > a.cidx(j) + 1) || |
|
|
550 ((a.cidx(j+1) == a.cidx(j) + 1) && found [a.ridx(j)])) |
|
|
551 { |
|
|
552 tmp_typ = MatrixType::Full; |
|
|
553 break; |
|
|
554 } |
|
|
555 found [a.ridx(j)] = true; |
|
|
556 } |
|
|
557 } |
|
|
558 typ = tmp_typ; |
|
|
559 } |
|
|
560 |
|
|
561 if (typ == MatrixType::Full) |
|
|
562 { |
|
|
563 // Search for banded, upper and lower triangular matrices |
|
|
564 bool singular = false; |
|
|
565 upper_band = 0; |
|
|
566 lower_band = 0; |
|
|
567 for (octave_idx_type j = 0; j < ncols; j++) |
|
|
568 { |
|
|
569 bool zero_on_diagonal = false; |
|
|
570 if (j < nrows) |
|
|
571 { |
|
|
572 zero_on_diagonal = true; |
|
|
573 for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++) |
|
|
574 if (a.ridx(i) == j) |
|
|
575 { |
|
|
576 zero_on_diagonal = false; |
|
|
577 break; |
|
|
578 } |
|
|
579 } |
|
|
580 |
|
|
581 if (zero_on_diagonal) |
|
|
582 { |
|
|
583 singular = true; |
|
|
584 break; |
|
|
585 } |
|
|
586 |
|
|
587 if (a.cidx(j+1) != a.cidx(j)) |
|
|
588 { |
|
|
589 octave_idx_type ru = a.ridx(a.cidx(j)); |
|
|
590 octave_idx_type rl = a.ridx(a.cidx(j+1)-1); |
|
|
591 |
|
|
592 if (j - ru > upper_band) |
|
|
593 upper_band = j - ru; |
|
|
594 |
|
|
595 if (rl - j > lower_band) |
|
|
596 lower_band = rl - j; |
|
|
597 } |
|
|
598 } |
|
|
599 |
|
|
600 if (!singular) |
|
|
601 { |
|
|
602 bandden = double (nnz) / |
|
|
603 (double (ncols) * (double (lower_band) + |
|
|
604 double (upper_band)) - |
|
|
605 0.5 * double (upper_band + 1) * double (upper_band) - |
|
|
606 0.5 * double (lower_band + 1) * double (lower_band)); |
|
|
607 |
|
|
608 if (nrows == ncols && sp_bandden != 1. && bandden > sp_bandden) |
|
|
609 { |
|
|
610 if (upper_band == 1 && lower_band == 1) |
|
|
611 typ = MatrixType::Tridiagonal; |
|
|
612 else |
|
|
613 typ = MatrixType::Banded; |
|
|
614 |
|
|
615 octave_idx_type nnz_in_band = |
|
|
616 (upper_band + lower_band + 1) * nrows - |
|
|
617 (1 + upper_band) * upper_band / 2 - |
|
|
618 (1 + lower_band) * lower_band / 2; |
|
|
619 if (nnz_in_band == nnz) |
|
|
620 dense = true; |
|
|
621 else |
|
|
622 dense = false; |
|
|
623 } |
|
|
624 else if (upper_band == 0) |
|
|
625 typ = MatrixType::Lower; |
|
|
626 else if (lower_band == 0) |
|
|
627 typ = MatrixType::Upper; |
|
|
628 |
|
|
629 if (upper_band == lower_band && nrows == ncols) |
|
|
630 maybe_hermitian = true; |
|
|
631 } |
|
|
632 |
|
|
633 if (typ == MatrixType::Full) |
|
|
634 { |
|
|
635 // Search for a permuted triangular matrix, and test if |
|
|
636 // permutation is singular |
|
|
637 |
|
|
638 // FIXME |
|
|
639 // Perhaps this should be based on a dmperm algorithm |
|
|
640 bool found = false; |
|
|
641 |
|
|
642 nperm = ncols; |
|
|
643 perm = new octave_idx_type [ncols]; |
|
|
644 |
|
|
645 for (octave_idx_type i = 0; i < ncols; i++) |
|
|
646 perm [i] = -1; |
|
|
647 |
|
|
648 for (octave_idx_type i = 0; i < nm; i++) |
|
|
649 { |
|
|
650 found = false; |
|
|
651 |
|
|
652 for (octave_idx_type j = 0; j < ncols; j++) |
|
|
653 { |
|
|
654 if ((a.cidx(j+1) - a.cidx(j)) > 0 && |
|
|
655 (a.ridx(a.cidx(j+1)-1) == i)) |
|
|
656 { |
|
|
657 perm [i] = j; |
|
|
658 found = true; |
|
|
659 break; |
|
|
660 } |
|
|
661 } |
|
|
662 |
|
|
663 if (!found) |
|
|
664 break; |
|
|
665 } |
|
|
666 |
|
|
667 if (found) |
|
|
668 { |
|
|
669 typ = MatrixType::Permuted_Upper; |
|
|
670 if (ncols > nrows) |
|
|
671 { |
|
|
672 octave_idx_type k = nrows; |
|
|
673 for (octave_idx_type i = 0; i < ncols; i++) |
|
|
674 if (perm [i] == -1) |
|
|
675 perm[i] = k++; |
|
|
676 } |
|
|
677 } |
|
|
678 else if (a.cidx(nm) == a.cidx(ncols)) |
|
|
679 { |
|
|
680 nperm = nrows; |
|
|
681 delete [] perm; |
|
|
682 perm = new octave_idx_type [nrows]; |
|
|
683 OCTAVE_LOCAL_BUFFER (octave_idx_type, tmp, nrows); |
|
|
684 |
|
|
685 for (octave_idx_type i = 0; i < nrows; i++) |
|
|
686 { |
|
|
687 perm [i] = -1; |
|
|
688 tmp [i] = -1; |
|
|
689 } |
|
|
690 |
|
|
691 for (octave_idx_type j = 0; j < ncols; j++) |
|
|
692 for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++) |
|
|
693 perm [a.ridx(i)] = j; |
|
|
694 |
|
|
695 found = true; |
|
|
696 for (octave_idx_type i = 0; i < nm; i++) |
|
|
697 if (perm[i] == -1) |
|
|
698 { |
|
|
699 found = false; |
|
|
700 break; |
|
|
701 } |
|
|
702 else |
|
|
703 { |
|
|
704 tmp[perm[i]] = 1; |
|
|
705 } |
|
|
706 |
|
|
707 if (found) |
|
|
708 { |
|
|
709 octave_idx_type k = ncols; |
|
|
710 for (octave_idx_type i = 0; i < nrows; i++) |
|
|
711 { |
|
|
712 if (tmp[i] == -1) |
|
|
713 { |
|
|
714 if (k < nrows) |
|
|
715 { |
|
|
716 perm[k++] = i; |
|
|
717 } |
|
|
718 else |
|
|
719 { |
|
|
720 found = false; |
|
|
721 break; |
|
|
722 } |
|
|
723 } |
|
|
724 } |
|
|
725 } |
|
|
726 |
|
|
727 if (found) |
|
|
728 typ = MatrixType::Permuted_Lower; |
|
|
729 else |
|
|
730 { |
|
|
731 delete [] perm; |
|
|
732 nperm = 0; |
|
|
733 } |
|
|
734 } |
|
|
735 else |
|
|
736 { |
|
|
737 delete [] perm; |
|
|
738 nperm = 0; |
|
|
739 } |
|
|
740 } |
|
|
741 |
|
|
742 // FIXME |
|
|
743 // Disable lower under-determined and upper over-determined problems |
|
|
744 // as being detected, and force to treat as singular. As this seems |
|
|
745 // to cause issues |
|
|
746 if (((typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower) |
|
|
747 && nrows > ncols) || |
|
|
748 ((typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper) |
|
|
749 && nrows < ncols)) |
|
|
750 { |
|
|
751 typ = MatrixType::Rectangular; |
|
|
752 if (typ == MatrixType::Permuted_Upper || |
|
|
753 typ == MatrixType::Permuted_Lower) |
|
|
754 delete [] perm; |
|
|
755 nperm = 0; |
|
|
756 } |
|
|
757 |
|
|
758 if (typ == MatrixType::Full && ncols != nrows) |
|
|
759 typ = MatrixType::Rectangular; |
|
|
760 |
|
|
761 if (maybe_hermitian && (typ == MatrixType::Full || |
|
|
762 typ == MatrixType::Tridiagonal || |
|
|
763 typ == MatrixType::Banded)) |
|
|
764 { |
|
|
765 // Check for symmetry, with positive real diagonal, which |
|
|
766 // has a very good chance of being symmetric positive |
|
|
767 // definite.. |
|
|
768 bool is_herm = true; |
|
|
769 |
|
|
770 for (octave_idx_type j = 0; j < ncols; j++) |
|
|
771 { |
|
|
772 bool diag_positive = false; |
|
|
773 |
|
|
774 for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++) |
|
|
775 { |
|
|
776 octave_idx_type ri = a.ridx(i); |
|
|
777 |
|
|
778 if (ri == j) |
|
|
779 { |
|
|
780 if (a.data(i) == std::abs(a.data(i))) |
|
|
781 diag_positive = true; |
|
|
782 else |
|
|
783 break; |
|
|
784 } |
|
|
785 else |
|
|
786 { |
|
|
787 bool found = false; |
|
|
788 |
|
|
789 for (octave_idx_type k = a.cidx(ri); k < a.cidx(ri+1); k++) |
|
|
790 { |
|
|
791 if (a.ridx(k) == j) |
|
|
792 { |
|
|
793 if (a.data(i) == conj(a.data(k))) |
|
|
794 found = true; |
|
|
795 break; |
|
|
796 } |
|
|
797 } |
|
|
798 |
|
|
799 if (! found) |
|
|
800 { |
|
|
801 is_herm = false; |
|
|
802 break; |
|
|
803 } |
|
|
804 } |
|
|
805 } |
|
|
806 |
|
|
807 if (! diag_positive || ! is_herm) |
|
|
808 { |
|
|
809 is_herm = false; |
|
|
810 break; |
|
|
811 } |
|
|
812 } |
|
|
813 |
|
|
814 if (is_herm) |
|
|
815 { |
|
|
816 if (typ == MatrixType::Full) |
|
|
817 typ = MatrixType::Hermitian; |
|
|
818 else if (typ == MatrixType::Banded) |
|
|
819 typ = MatrixType::Banded_Hermitian; |
|
|
820 else |
|
|
821 typ = MatrixType::Tridiagonal_Hermitian; |
|
|
822 } |
|
|
823 } |
|
|
824 } |
|
|
825 } |
|
5892
|
826 MatrixType::MatrixType (const matrix_type t, bool _full) |
|
|
827 : typ (MatrixType::Unknown), |
|
5893
|
828 sp_bandden (octave_sparse_params::get_key ("bandden")), |
|
5892
|
829 bandden (0), upper_band (0), lower_band (0), |
|
|
830 dense (false), full (_full), nperm (0), perm (0) |
|
5785
|
831 { |
|
|
832 if (t == MatrixType::Full || t == MatrixType::Diagonal || |
|
|
833 t == MatrixType::Permuted_Diagonal || t == MatrixType::Upper || |
|
|
834 t == MatrixType::Lower || t == MatrixType::Tridiagonal || |
|
|
835 t == MatrixType::Tridiagonal_Hermitian || t == MatrixType::Rectangular) |
|
|
836 typ = t; |
|
|
837 else |
|
|
838 (*current_liboctave_warning_handler) ("Invalid matrix type"); |
|
|
839 } |
|
|
840 |
|
|
841 MatrixType::MatrixType (const matrix_type t, const octave_idx_type np, |
|
5892
|
842 const octave_idx_type *p, bool _full) |
|
|
843 : typ (MatrixType::Unknown), |
|
5893
|
844 sp_bandden (octave_sparse_params::get_key ("bandden")), |
|
5892
|
845 bandden (0), upper_band (0), lower_band (0), |
|
|
846 dense (false), full (_full), nperm (0), perm (0) |
|
5785
|
847 { |
|
6027
|
848 if ((t == MatrixType::Permuted_Upper || t == MatrixType::Permuted_Lower) && |
|
|
849 np > 0 && p != 0) |
|
5785
|
850 { |
|
|
851 typ = t; |
|
|
852 nperm = np; |
|
|
853 perm = new octave_idx_type [nperm]; |
|
|
854 for (octave_idx_type i = 0; i < nperm; i++) |
|
|
855 perm[i] = p[i]; |
|
|
856 } |
|
|
857 else |
|
|
858 (*current_liboctave_warning_handler) ("Invalid matrix type"); |
|
|
859 } |
|
|
860 |
|
|
861 MatrixType::MatrixType (const matrix_type t, const octave_idx_type ku, |
|
5892
|
862 const octave_idx_type kl, bool _full) |
|
|
863 : typ (MatrixType::Unknown), |
|
5893
|
864 sp_bandden (octave_sparse_params::get_key ("bandden")), |
|
5892
|
865 bandden (0), upper_band (0), lower_band (0), |
|
|
866 dense (false), full (_full), nperm (0), perm (0) |
|
5785
|
867 { |
|
|
868 if (t == MatrixType::Banded || t == MatrixType::Banded_Hermitian) |
|
|
869 { |
|
|
870 typ = t; |
|
|
871 upper_band = ku; |
|
|
872 lower_band = kl; |
|
|
873 } |
|
|
874 else |
|
|
875 (*current_liboctave_warning_handler) ("Invalid sparse matrix type"); |
|
|
876 } |
|
|
877 |
|
|
878 MatrixType::~MatrixType (void) |
|
|
879 { |
|
|
880 if (nperm != 0) |
|
|
881 { |
|
|
882 delete [] perm; |
|
|
883 } |
|
|
884 } |
|
|
885 |
|
|
886 MatrixType& |
|
|
887 MatrixType::operator = (const MatrixType& a) |
|
|
888 { |
|
|
889 if (this != &a) |
|
|
890 { |
|
|
891 typ = a.typ; |
|
|
892 sp_bandden = a.sp_bandden; |
|
|
893 bandden = a.bandden; |
|
|
894 upper_band = a.upper_band; |
|
|
895 lower_band = a.lower_band; |
|
|
896 dense = a.dense; |
|
|
897 full = a.full; |
|
|
898 nperm = a.nperm; |
|
|
899 |
|
|
900 if (nperm != 0) |
|
|
901 { |
|
|
902 perm = new octave_idx_type [nperm]; |
|
|
903 for (octave_idx_type i = 0; i < nperm; i++) |
|
|
904 perm[i] = a.perm[i]; |
|
|
905 } |
|
|
906 } |
|
|
907 |
|
|
908 return *this; |
|
|
909 } |
|
|
910 |
|
|
911 int |
|
|
912 MatrixType::type (bool quiet) |
|
|
913 { |
|
|
914 if (typ != MatrixType::Unknown && (full || |
|
5893
|
915 sp_bandden == octave_sparse_params::get_key ("bandden"))) |
|
5785
|
916 { |
|
|
917 if (!quiet && |
|
5893
|
918 octave_sparse_params::get_key ("spumoni") != 0.) |
|
5785
|
919 (*current_liboctave_warning_handler) |
|
|
920 ("Using Cached Matrix Type"); |
|
|
921 |
|
|
922 return typ; |
|
|
923 } |
|
|
924 |
|
|
925 if (typ != MatrixType::Unknown && |
|
5893
|
926 octave_sparse_params::get_key ("spumoni") != 0.) |
|
5785
|
927 (*current_liboctave_warning_handler) |
|
|
928 ("Invalidating Matrix Type"); |
|
|
929 |
|
|
930 typ = MatrixType::Unknown; |
|
|
931 |
|
|
932 return typ; |
|
|
933 } |
|
|
934 |
|
|
935 int |
|
|
936 MatrixType::type (const SparseMatrix &a) |
|
|
937 { |
|
|
938 if (typ != MatrixType::Unknown && (full || |
|
5893
|
939 sp_bandden == octave_sparse_params::get_key ("bandden"))) |
|
5785
|
940 { |
|
5893
|
941 if (octave_sparse_params::get_key ("spumoni") != 0.) |
|
5785
|
942 (*current_liboctave_warning_handler) |
|
|
943 ("Using Cached Matrix Type"); |
|
|
944 |
|
|
945 return typ; |
|
|
946 } |
|
|
947 |
|
|
948 MatrixType tmp_typ (a); |
|
|
949 typ = tmp_typ.typ; |
|
|
950 sp_bandden = tmp_typ.sp_bandden; |
|
|
951 bandden = tmp_typ.bandden; |
|
|
952 upper_band = tmp_typ.upper_band; |
|
|
953 lower_band = tmp_typ.lower_band; |
|
|
954 dense = tmp_typ.dense; |
|
|
955 full = tmp_typ.full; |
|
|
956 nperm = tmp_typ.nperm; |
|
|
957 |
|
|
958 if (nperm != 0) |
|
|
959 { |
|
|
960 perm = new octave_idx_type [nperm]; |
|
|
961 for (octave_idx_type i = 0; i < nperm; i++) |
|
|
962 perm[i] = tmp_typ.perm[i]; |
|
|
963 } |
|
|
964 |
|
|
965 return typ; |
|
|
966 } |
|
|
967 |
|
|
968 int |
|
|
969 MatrixType::type (const SparseComplexMatrix &a) |
|
|
970 { |
|
|
971 if (typ != MatrixType::Unknown && (full || |
|
5893
|
972 sp_bandden == octave_sparse_params::get_key ("bandden"))) |
|
5785
|
973 { |
|
5893
|
974 if (octave_sparse_params::get_key ("spumoni") != 0.) |
|
5785
|
975 (*current_liboctave_warning_handler) |
|
|
976 ("Using Cached Matrix Type"); |
|
|
977 |
|
|
978 return typ; |
|
|
979 } |
|
|
980 |
|
|
981 MatrixType tmp_typ (a); |
|
|
982 typ = tmp_typ.typ; |
|
|
983 sp_bandden = tmp_typ.sp_bandden; |
|
|
984 bandden = tmp_typ.bandden; |
|
|
985 upper_band = tmp_typ.upper_band; |
|
|
986 lower_band = tmp_typ.lower_band; |
|
|
987 dense = tmp_typ.dense; |
|
|
988 full = tmp_typ.full; |
|
|
989 nperm = tmp_typ.nperm; |
|
|
990 |
|
|
991 if (nperm != 0) |
|
|
992 { |
|
|
993 perm = new octave_idx_type [nperm]; |
|
|
994 for (octave_idx_type i = 0; i < nperm; i++) |
|
|
995 perm[i] = tmp_typ.perm[i]; |
|
|
996 } |
|
|
997 |
|
|
998 return typ; |
|
|
999 } |
|
|
1000 int |
|
|
1001 MatrixType::type (const Matrix &a) |
|
|
1002 { |
|
|
1003 if (typ != MatrixType::Unknown) |
|
|
1004 { |
|
5893
|
1005 if (octave_sparse_params::get_key ("spumoni") != 0.) |
|
5785
|
1006 (*current_liboctave_warning_handler) |
|
|
1007 ("Using Cached Matrix Type"); |
|
|
1008 |
|
|
1009 return typ; |
|
|
1010 } |
|
|
1011 |
|
|
1012 MatrixType tmp_typ (a); |
|
|
1013 typ = tmp_typ.typ; |
|
|
1014 full = tmp_typ.full; |
|
|
1015 nperm = tmp_typ.nperm; |
|
|
1016 |
|
|
1017 if (nperm != 0) |
|
|
1018 { |
|
|
1019 perm = new octave_idx_type [nperm]; |
|
|
1020 for (octave_idx_type i = 0; i < nperm; i++) |
|
|
1021 perm[i] = tmp_typ.perm[i]; |
|
|
1022 } |
|
|
1023 |
|
|
1024 return typ; |
|
|
1025 } |
|
|
1026 |
|
|
1027 int |
|
|
1028 MatrixType::type (const ComplexMatrix &a) |
|
|
1029 { |
|
|
1030 if (typ != MatrixType::Unknown) |
|
|
1031 { |
|
5893
|
1032 if (octave_sparse_params::get_key ("spumoni") != 0.) |
|
5785
|
1033 (*current_liboctave_warning_handler) |
|
|
1034 ("Using Cached Matrix Type"); |
|
|
1035 |
|
|
1036 return typ; |
|
|
1037 } |
|
|
1038 |
|
|
1039 MatrixType tmp_typ (a); |
|
|
1040 typ = tmp_typ.typ; |
|
|
1041 full = tmp_typ.full; |
|
|
1042 nperm = tmp_typ.nperm; |
|
|
1043 |
|
|
1044 if (nperm != 0) |
|
|
1045 { |
|
|
1046 perm = new octave_idx_type [nperm]; |
|
|
1047 for (octave_idx_type i = 0; i < nperm; i++) |
|
|
1048 perm[i] = tmp_typ.perm[i]; |
|
|
1049 } |
|
|
1050 |
|
|
1051 return typ; |
|
|
1052 } |
|
|
1053 |
|
|
1054 void |
|
|
1055 MatrixType::info () const |
|
|
1056 { |
|
5893
|
1057 if (octave_sparse_params::get_key ("spumoni") != 0.) |
|
5785
|
1058 { |
|
|
1059 if (typ == MatrixType::Unknown) |
|
|
1060 (*current_liboctave_warning_handler) |
|
|
1061 ("Unknown Matrix Type"); |
|
|
1062 else if (typ == MatrixType::Diagonal) |
|
|
1063 (*current_liboctave_warning_handler) |
|
|
1064 ("Diagonal Sparse Matrix"); |
|
|
1065 else if (typ == MatrixType::Permuted_Diagonal) |
|
|
1066 (*current_liboctave_warning_handler) |
|
|
1067 ("Permuted Diagonal Sparse Matrix"); |
|
|
1068 else if (typ == MatrixType::Upper) |
|
|
1069 (*current_liboctave_warning_handler) |
|
|
1070 ("Upper Triangular Matrix"); |
|
|
1071 else if (typ == MatrixType::Lower) |
|
|
1072 (*current_liboctave_warning_handler) |
|
|
1073 ("Lower Triangular Matrix"); |
|
|
1074 else if (typ == MatrixType::Permuted_Upper) |
|
|
1075 (*current_liboctave_warning_handler) |
|
|
1076 ("Permuted Upper Triangular Matrix"); |
|
|
1077 else if (typ == MatrixType::Permuted_Lower) |
|
|
1078 (*current_liboctave_warning_handler) |
|
|
1079 ("Permuted Lower Triangular Matrix"); |
|
|
1080 else if (typ == MatrixType::Banded) |
|
|
1081 (*current_liboctave_warning_handler) |
|
|
1082 ("Banded Sparse Matrix %d-1-%d (Density %f)", lower_band, |
|
|
1083 upper_band, bandden); |
|
|
1084 else if (typ == MatrixType::Banded_Hermitian) |
|
|
1085 (*current_liboctave_warning_handler) |
|
|
1086 ("Banded Hermitian/Symmetric Sparse Matrix %d-1-%d (Density %f)", |
|
|
1087 lower_band, upper_band, bandden); |
|
|
1088 else if (typ == MatrixType::Hermitian) |
|
|
1089 (*current_liboctave_warning_handler) |
|
|
1090 ("Hermitian/Symmetric Matrix"); |
|
|
1091 else if (typ == MatrixType::Tridiagonal) |
|
|
1092 (*current_liboctave_warning_handler) |
|
|
1093 ("Tridiagonal Sparse Matrix"); |
|
|
1094 else if (typ == MatrixType::Tridiagonal_Hermitian) |
|
|
1095 (*current_liboctave_warning_handler) |
|
|
1096 ("Hermitian/Symmetric Tridiagonal Sparse Matrix"); |
|
|
1097 else if (typ == MatrixType::Rectangular) |
|
|
1098 (*current_liboctave_warning_handler) |
|
|
1099 ("Rectangular/Singular Matrix"); |
|
|
1100 else if (typ == MatrixType::Full) |
|
|
1101 (*current_liboctave_warning_handler) |
|
|
1102 ("Full Matrix"); |
|
|
1103 } |
|
|
1104 } |
|
|
1105 |
|
|
1106 void |
|
|
1107 MatrixType::mark_as_symmetric (void) |
|
|
1108 { |
|
|
1109 if (typ == MatrixType::Tridiagonal || |
|
|
1110 typ == MatrixType::Tridiagonal_Hermitian) |
|
|
1111 typ = MatrixType::Tridiagonal_Hermitian; |
|
|
1112 else if (typ == MatrixType::Banded || |
|
|
1113 typ == MatrixType::Banded_Hermitian) |
|
|
1114 typ = MatrixType::Banded_Hermitian; |
|
|
1115 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian || |
|
|
1116 typ == MatrixType::Unknown) |
|
|
1117 typ = MatrixType::Hermitian; |
|
|
1118 else |
|
|
1119 (*current_liboctave_error_handler) |
|
|
1120 ("Can not mark current matrix type as symmetric"); |
|
|
1121 } |
|
|
1122 |
|
|
1123 void |
|
|
1124 MatrixType::mark_as_unsymmetric (void) |
|
|
1125 { |
|
|
1126 if (typ == MatrixType::Tridiagonal || |
|
|
1127 typ == MatrixType::Tridiagonal_Hermitian) |
|
|
1128 typ = MatrixType::Tridiagonal; |
|
|
1129 else if (typ == MatrixType::Banded || |
|
|
1130 typ == MatrixType::Banded_Hermitian) |
|
|
1131 typ = MatrixType::Banded; |
|
|
1132 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian || |
|
|
1133 typ == MatrixType::Unknown) |
|
|
1134 typ = MatrixType::Full; |
|
|
1135 } |
|
|
1136 |
|
|
1137 void |
|
|
1138 MatrixType::mark_as_permuted (const octave_idx_type np, const octave_idx_type *p) |
|
|
1139 { |
|
|
1140 nperm = np; |
|
|
1141 perm = new octave_idx_type [nperm]; |
|
|
1142 for (octave_idx_type i = 0; i < nperm; i++) |
|
|
1143 perm[i] = p[i]; |
|
|
1144 |
|
|
1145 if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal) |
|
|
1146 typ = MatrixType::Permuted_Diagonal; |
|
|
1147 else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper) |
|
|
1148 typ = MatrixType::Permuted_Upper; |
|
|
1149 else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower) |
|
|
1150 typ = MatrixType::Permuted_Lower; |
|
|
1151 else |
|
|
1152 (*current_liboctave_error_handler) |
|
|
1153 ("Can not mark current matrix type as symmetric"); |
|
|
1154 } |
|
|
1155 |
|
|
1156 void |
|
|
1157 MatrixType::mark_as_unpermuted (void) |
|
|
1158 { |
|
|
1159 if (nperm) |
|
|
1160 { |
|
|
1161 nperm = 0; |
|
|
1162 delete [] perm; |
|
|
1163 } |
|
|
1164 |
|
|
1165 if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal) |
|
|
1166 typ = MatrixType::Diagonal; |
|
|
1167 else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper) |
|
|
1168 typ = MatrixType::Upper; |
|
|
1169 else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower) |
|
|
1170 typ = MatrixType::Lower; |
|
|
1171 } |
|
|
1172 |
|
|
1173 MatrixType |
|
|
1174 MatrixType::transpose (void) const |
|
|
1175 { |
|
|
1176 MatrixType retval (*this); |
|
|
1177 if (typ == MatrixType::Upper) |
|
|
1178 retval.typ = MatrixType::Lower; |
|
|
1179 else if (typ == MatrixType::Permuted_Upper) |
|
|
1180 retval.typ = MatrixType::Permuted_Lower; |
|
|
1181 else if (typ == MatrixType::Lower) |
|
|
1182 retval.typ = MatrixType::Upper; |
|
|
1183 else if (typ == MatrixType::Permuted_Lower) |
|
|
1184 retval.typ = MatrixType::Permuted_Upper; |
|
|
1185 else if (typ == MatrixType::Banded) |
|
|
1186 { |
|
|
1187 retval.upper_band = lower_band; |
|
|
1188 retval.lower_band = upper_band; |
|
|
1189 } |
|
|
1190 |
|
|
1191 return retval; |
|
|
1192 } |
|
|
1193 |
|
|
1194 /* |
|
|
1195 ;;; Local Variables: *** |
|
|
1196 ;;; mode: C++ *** |
|
|
1197 ;;; End: *** |
|
|
1198 */ |
|
|
1199 |
