1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 2006-2021 The Octave Project Developers
4 //
5 // See the file COPYRIGHT.md in the top-level directory of this
6 // distribution or <https://octave.org/copyright/>.
7 //
8 // This file is part of Octave.
9 //
10 // Octave is free software: you can redistribute it and/or modify it
11 // under the terms of the GNU General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // Octave is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU General Public License for more details.
19 //
20 // You should have received a copy of the GNU General Public License
21 // along with Octave; see the file COPYING. If not, see
22 // <https://www.gnu.org/licenses/>.
23 //
24 ////////////////////////////////////////////////////////////////////////
25
26 #if defined (HAVE_CONFIG_H)
27 # include "config.h"
28 #endif
29
30 #include <cinttypes>
31 #include <vector>
32
33 #include "MatrixType.h"
34 #include "dMatrix.h"
35 #include "fMatrix.h"
36 #include "CMatrix.h"
37 #include "fCMatrix.h"
38 #include "dSparse.h"
39 #include "CSparse.h"
40 #include "oct-spparms.h"
41 #include "oct-locbuf.h"
42
43 static void
warn_cached(void)44 warn_cached (void)
45 {
46 (*current_liboctave_warning_with_id_handler)
47 ("Octave:matrix-type-info", "using cached matrix type");
48 }
49
50 static void
warn_invalid(void)51 warn_invalid (void)
52 {
53 (*current_liboctave_warning_with_id_handler)
54 ("Octave:matrix-type-info", "invalid matrix type");
55 }
56
57 static void
warn_calculating_sparse_type(void)58 warn_calculating_sparse_type (void)
59 {
60 (*current_liboctave_warning_with_id_handler)
61 ("Octave:matrix-type-info", "calculating sparse matrix type");
62 }
63
64 // FIXME: There is a large code duplication here
65
MatrixType(void)66 MatrixType::MatrixType (void)
67 : typ (MatrixType::Unknown),
68 sp_bandden (octave_sparse_params::get_bandden ()),
69 bandden (0), upper_band (0),
70 lower_band (0), dense (false), full (false), nperm (0), perm (nullptr) { }
71
MatrixType(const MatrixType & a)72 MatrixType::MatrixType (const MatrixType& a)
73 : typ (a.typ), sp_bandden (a.sp_bandden), bandden (a.bandden),
74 upper_band (a.upper_band), lower_band (a.lower_band),
75 dense (a.dense), full (a.full), nperm (a.nperm), perm (nullptr)
76 {
77 if (nperm != 0)
78 {
79 perm = new octave_idx_type [nperm];
80 for (octave_idx_type i = 0; i < nperm; i++)
81 perm[i] = a.perm[i];
82 }
83 }
84
85 template <typename T>
86 MatrixType::matrix_type
matrix_real_probe(const MArray<T> & a)87 matrix_real_probe (const MArray<T>& a)
88 {
89 MatrixType::matrix_type typ;
90 octave_idx_type nrows = a.rows ();
91 octave_idx_type ncols = a.cols ();
92
93 const T zero = 0;
94
95 if (ncols == nrows)
96 {
97 bool upper = true;
98 bool lower = true;
99 bool hermitian = true;
100
101 // do the checks for lower/upper/hermitian all in one pass.
102 OCTAVE_LOCAL_BUFFER (T, diag, ncols);
103
104 for (octave_idx_type j = 0;
105 j < ncols && upper; j++)
106 {
107 T d = a.elem (j,j);
108 upper = upper && (d != zero);
109 lower = lower && (d != zero);
110 hermitian = hermitian && (d > zero);
111 diag[j] = d;
112 }
113
114 for (octave_idx_type j = 0;
115 j < ncols && (upper || lower || hermitian); j++)
116 {
117 for (octave_idx_type i = 0; i < j; i++)
118 {
119 double aij = a.elem (i,j);
120 double aji = a.elem (j,i);
121 lower = lower && (aij == zero);
122 upper = upper && (aji == zero);
123 hermitian = hermitian && (aij == aji
124 && aij*aij < diag[i]*diag[j]);
125 }
126 }
127
128 if (upper)
129 typ = MatrixType::Upper;
130 else if (lower)
131 typ = MatrixType::Lower;
132 else if (hermitian)
133 typ = MatrixType::Hermitian;
134 else
135 typ = MatrixType::Full;
136 }
137 else
138 typ = MatrixType::Rectangular;
139
140 return typ;
141 }
142
143 template <typename T>
144 MatrixType::matrix_type
matrix_complex_probe(const MArray<std::complex<T>> & a)145 matrix_complex_probe (const MArray<std::complex<T>>& a)
146 {
147 MatrixType::matrix_type typ = MatrixType::Unknown;
148 octave_idx_type nrows = a.rows ();
149 octave_idx_type ncols = a.cols ();
150
151 const T zero = 0;
152 // get the real type
153
154 if (ncols == nrows)
155 {
156 bool upper = true;
157 bool lower = true;
158 bool hermitian = true;
159
160 // do the checks for lower/upper/hermitian all in one pass.
161 OCTAVE_LOCAL_BUFFER (T, diag, ncols);
162
163 for (octave_idx_type j = 0;
164 j < ncols && upper; j++)
165 {
166 std::complex<T> d = a.elem (j,j);
167 upper = upper && (d != zero);
168 lower = lower && (d != zero);
169 hermitian = hermitian && (d.real () > zero && d.imag () == zero);
170 diag[j] = d.real ();
171 }
172
173 for (octave_idx_type j = 0;
174 j < ncols && (upper || lower || hermitian); j++)
175 {
176 for (octave_idx_type i = 0; i < j; i++)
177 {
178 std::complex<T> aij = a.elem (i,j);
179 std::complex<T> aji = a.elem (j,i);
180 lower = lower && (aij == zero);
181 upper = upper && (aji == zero);
182 hermitian = hermitian && (aij == octave::math::conj (aji)
183 && std::norm (aij) < diag[i]*diag[j]);
184 }
185 }
186
187 if (upper)
188 typ = MatrixType::Upper;
189 else if (lower)
190 typ = MatrixType::Lower;
191 else if (hermitian)
192 typ = MatrixType::Hermitian;
193 else
194 typ = MatrixType::Full;
195 }
196 else
197 typ = MatrixType::Rectangular;
198
199 return typ;
200 }
201
MatrixType(const Matrix & a)202 MatrixType::MatrixType (const Matrix& a)
203 : typ (MatrixType::Unknown),
204 sp_bandden (0), bandden (0), upper_band (0), lower_band (0),
205 dense (false), full (true), nperm (0), perm (nullptr)
206 {
207 typ = matrix_real_probe (a);
208 }
209
MatrixType(const ComplexMatrix & a)210 MatrixType::MatrixType (const ComplexMatrix& a)
211 : typ (MatrixType::Unknown),
212 sp_bandden (0), bandden (0), upper_band (0), lower_band (0),
213 dense (false), full (true), nperm (0), perm (nullptr)
214 {
215 typ = matrix_complex_probe (a);
216 }
217
MatrixType(const FloatMatrix & a)218 MatrixType::MatrixType (const FloatMatrix& a)
219 : typ (MatrixType::Unknown),
220 sp_bandden (0), bandden (0), upper_band (0), lower_band (0),
221 dense (false), full (true), nperm (0), perm (nullptr)
222 {
223 typ = matrix_real_probe (a);
224 }
225
MatrixType(const FloatComplexMatrix & a)226 MatrixType::MatrixType (const FloatComplexMatrix& a)
227 : typ (MatrixType::Unknown),
228 sp_bandden (0), bandden (0), upper_band (0), lower_band (0),
229 dense (false), full (true), nperm (0), perm (nullptr)
230 {
231 typ = matrix_complex_probe (a);
232 }
233
234
235 template <typename T>
MatrixType(const MSparse<T> & a)236 MatrixType::MatrixType (const MSparse<T>& a)
237 : typ (MatrixType::Unknown),
238 sp_bandden (0), bandden (0), upper_band (0), lower_band (0),
239 dense (false), full (false), nperm (0), perm (nullptr)
240 {
241 octave_idx_type nrows = a.rows ();
242 octave_idx_type ncols = a.cols ();
243 octave_idx_type nm = (ncols < nrows ? ncols : nrows);
244 octave_idx_type nnz = a.nnz ();
245
246 if (octave_sparse_params::get_key ("spumoni") != 0.)
247 warn_calculating_sparse_type ();
248
249 sp_bandden = octave_sparse_params::get_bandden ();
250 bool maybe_hermitian = false;
251 typ = MatrixType::Full;
252
253 if (nnz == nm)
254 {
255 matrix_type tmp_typ = MatrixType::Diagonal;
256 octave_idx_type i;
257 // Maybe the matrix is diagonal
258 for (i = 0; i < nm; i++)
259 {
260 if (a.cidx (i+1) != a.cidx (i) + 1)
261 {
262 tmp_typ = MatrixType::Full;
263 break;
264 }
265 if (a.ridx (i) != i)
266 {
267 tmp_typ = MatrixType::Permuted_Diagonal;
268 break;
269 }
270 }
271
272 if (tmp_typ == MatrixType::Permuted_Diagonal)
273 {
274 std::vector<bool> found (nrows);
275
276 for (octave_idx_type j = 0; j < i; j++)
277 found[j] = true;
278 for (octave_idx_type j = i; j < nrows; j++)
279 found[j] = false;
280
281 for (octave_idx_type j = i; j < nm; j++)
282 {
283 if ((a.cidx (j+1) > a.cidx (j) + 1)
284 || ((a.cidx (j+1) == a.cidx (j) + 1) && found[a.ridx (j)]))
285 {
286 tmp_typ = MatrixType::Full;
287 break;
288 }
289 found[a.ridx (j)] = true;
290 }
291 }
292 typ = tmp_typ;
293 }
294
295 if (typ == MatrixType::Full)
296 {
297 // Search for banded, upper and lower triangular matrices
298 bool singular = false;
299 upper_band = 0;
300 lower_band = 0;
301 for (octave_idx_type j = 0; j < ncols; j++)
302 {
303 bool zero_on_diagonal = false;
304 if (j < nrows)
305 {
306 zero_on_diagonal = true;
307 for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++)
308 if (a.ridx (i) == j)
309 {
310 zero_on_diagonal = false;
311 break;
312 }
313 }
314
315 if (zero_on_diagonal)
316 {
317 singular = true;
318 break;
319 }
320
321 if (a.cidx (j+1) != a.cidx (j))
322 {
323 octave_idx_type ru = a.ridx (a.cidx (j));
324 octave_idx_type rl = a.ridx (a.cidx (j+1)-1);
325
326 if (j - ru > upper_band)
327 upper_band = j - ru;
328
329 if (rl - j > lower_band)
330 lower_band = rl - j;
331 }
332 }
333
334 if (! singular)
335 {
336 bandden = double (nnz) /
337 (double (ncols) * (double (lower_band) +
338 double (upper_band)) -
339 0.5 * double (upper_band + 1) * double (upper_band) -
340 0.5 * double (lower_band + 1) * double (lower_band));
341
342 if (nrows == ncols && sp_bandden != 1. && bandden > sp_bandden)
343 {
344 if (upper_band == 1 && lower_band == 1)
345 typ = MatrixType::Tridiagonal;
346 else
347 typ = MatrixType::Banded;
348
349 octave_idx_type nnz_in_band
350 = ((upper_band + lower_band + 1) * nrows
351 - (1 + upper_band) * upper_band / 2
352 - (1 + lower_band) * lower_band / 2);
353
354 if (nnz_in_band == nnz)
355 dense = true;
356 else
357 dense = false;
358 }
359
360 // If a matrix is Banded but also Upper/Lower, set to the latter.
361 if (upper_band == 0)
362 typ = MatrixType::Lower;
363 else if (lower_band == 0)
364 typ = MatrixType::Upper;
365
366 if (upper_band == lower_band && nrows == ncols)
367 maybe_hermitian = true;
368 }
369
370 if (typ == MatrixType::Full)
371 {
372 // Search for a permuted triangular matrix, and test if
373 // permutation is singular
374
375 // FIXME: Perhaps this should be based on a dmperm algorithm?
376 bool found = false;
377
378 nperm = ncols;
379 perm = new octave_idx_type [ncols];
380
381 for (octave_idx_type i = 0; i < ncols; i++)
382 perm[i] = -1;
383
384 for (octave_idx_type i = 0; i < nm; i++)
385 {
386 found = false;
387
388 for (octave_idx_type j = 0; j < ncols; j++)
389 {
390 if ((a.cidx (j+1) - a.cidx (j)) > 0
391 && (a.ridx (a.cidx (j+1)-1) == i))
392 {
393 perm[i] = j;
394 found = true;
395 break;
396 }
397 }
398
399 if (! found)
400 break;
401 }
402
403 if (found)
404 {
405 typ = MatrixType::Permuted_Upper;
406 if (ncols > nrows)
407 {
408 octave_idx_type k = nrows;
409 for (octave_idx_type i = 0; i < ncols; i++)
410 if (perm[i] == -1)
411 perm[i] = k++;
412 }
413 }
414 else if (a.cidx (nm) == a.cidx (ncols))
415 {
416 nperm = nrows;
417 delete [] perm;
418 perm = new octave_idx_type [nrows];
419 OCTAVE_LOCAL_BUFFER (octave_idx_type, tmp, nrows);
420
421 for (octave_idx_type i = 0; i < nrows; i++)
422 {
423 perm[i] = -1;
424 tmp[i] = -1;
425 }
426
427 for (octave_idx_type j = 0; j < ncols; j++)
428 for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++)
429 perm[a.ridx (i)] = j;
430
431 found = true;
432 for (octave_idx_type i = 0; i < nm; i++)
433 if (perm[i] == -1)
434 {
435 found = false;
436 break;
437 }
438 else
439 {
440 tmp[perm[i]] = 1;
441 }
442
443 if (found)
444 {
445 octave_idx_type k = ncols;
446 for (octave_idx_type i = 0; i < nrows; i++)
447 {
448 if (tmp[i] == -1)
449 {
450 if (k < nrows)
451 {
452 perm[k++] = i;
453 }
454 else
455 {
456 found = false;
457 break;
458 }
459 }
460 }
461 }
462
463 if (found)
464 typ = MatrixType::Permuted_Lower;
465 else
466 {
467 delete [] perm;
468 nperm = 0;
469 }
470 }
471 else
472 {
473 delete [] perm;
474 nperm = 0;
475 }
476 }
477
478 // FIXME: Disable lower under-determined and upper over-determined
479 // problems as being detected, and force to treat as singular
480 // as this seems to cause issues.
481 if (((typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
482 && nrows > ncols)
483 || ((typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
484 && nrows < ncols))
485 {
486 if (typ == MatrixType::Permuted_Upper
487 || typ == MatrixType::Permuted_Lower)
488 delete [] perm;
489 nperm = 0;
490 typ = MatrixType::Rectangular;
491 }
492
493 if (typ == MatrixType::Full && ncols != nrows)
494 typ = MatrixType::Rectangular;
495
496 if (maybe_hermitian && (typ == MatrixType::Full
497 || typ == MatrixType::Tridiagonal
498 || typ == MatrixType::Banded))
499 {
500 bool is_herm = true;
501
502 // first, check whether the diagonal is positive & extract it
503 ColumnVector diag (ncols);
504
505 for (octave_idx_type j = 0; is_herm && j < ncols; j++)
506 {
507 is_herm = false;
508 for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++)
509 {
510 if (a.ridx (i) == j)
511 {
512 T d = a.data (i);
513 is_herm = (std::real (d) > 0.0
514 && std::imag (d) == 0.0);
515 diag(j) = std::real (d);
516 break;
517 }
518 }
519 }
520
521 // next, check symmetry and 2x2 positiveness
522
523 for (octave_idx_type j = 0; is_herm && j < ncols; j++)
524 for (octave_idx_type i = a.cidx (j); is_herm && i < a.cidx (j+1); i++)
525 {
526 octave_idx_type k = a.ridx (i);
527 is_herm = k == j;
528 if (is_herm)
529 continue;
530
531 T d = a.data (i);
532 if (std::norm (d) < diag(j)*diag(k))
533 {
534 d = octave::math::conj (d);
535 for (octave_idx_type l = a.cidx (k); l < a.cidx (k+1); l++)
536 {
537 if (a.ridx (l) == j)
538 {
539 is_herm = a.data (l) == d;
540 break;
541 }
542 }
543 }
544 }
545
546 if (is_herm)
547 {
548 if (typ == MatrixType::Full)
549 typ = MatrixType::Hermitian;
550 else if (typ == MatrixType::Banded)
551 typ = MatrixType::Banded_Hermitian;
552 else
553 typ = MatrixType::Tridiagonal_Hermitian;
554 }
555 }
556 }
557 }
558
559
MatrixType(const matrix_type t,bool _full)560 MatrixType::MatrixType (const matrix_type t, bool _full)
561 : typ (MatrixType::Unknown),
562 sp_bandden (octave_sparse_params::get_bandden ()),
563 bandden (0), upper_band (0), lower_band (0),
564 dense (false), full (_full), nperm (0), perm (nullptr)
565 {
566 if (t == MatrixType::Unknown || t == MatrixType::Full
567 || t == MatrixType::Diagonal || t == MatrixType::Permuted_Diagonal
568 || t == MatrixType::Upper || t == MatrixType::Lower
569 || t == MatrixType::Tridiagonal || t == MatrixType::Tridiagonal_Hermitian
570 || t == MatrixType::Rectangular)
571 typ = t;
572 else
573 warn_invalid ();
574 }
575
MatrixType(const matrix_type t,const octave_idx_type np,const octave_idx_type * p,bool _full)576 MatrixType::MatrixType (const matrix_type t, const octave_idx_type np,
577 const octave_idx_type *p, bool _full)
578 : typ (MatrixType::Unknown),
579 sp_bandden (octave_sparse_params::get_bandden ()),
580 bandden (0), upper_band (0), lower_band (0),
581 dense (false), full (_full), nperm (0), perm (nullptr)
582 {
583 if ((t == MatrixType::Permuted_Upper || t == MatrixType::Permuted_Lower)
584 && np > 0 && p != nullptr)
585 {
586 typ = t;
587 nperm = np;
588 perm = new octave_idx_type [nperm];
589 for (octave_idx_type i = 0; i < nperm; i++)
590 perm[i] = p[i];
591 }
592 else
593 warn_invalid ();
594 }
595
MatrixType(const matrix_type t,const octave_idx_type ku,const octave_idx_type kl,bool _full)596 MatrixType::MatrixType (const matrix_type t, const octave_idx_type ku,
597 const octave_idx_type kl, bool _full)
598 : typ (MatrixType::Unknown),
599 sp_bandden (octave_sparse_params::get_bandden ()),
600 bandden (0), upper_band (0), lower_band (0),
601 dense (false), full (_full), nperm (0), perm (nullptr)
602 {
603 if (t == MatrixType::Banded || t == MatrixType::Banded_Hermitian)
604 {
605 typ = t;
606 upper_band = ku;
607 lower_band = kl;
608 }
609 else
610 warn_invalid ();
611 }
612
~MatrixType(void)613 MatrixType::~MatrixType (void)
614 {
615 if (nperm != 0)
616 {
617 delete [] perm;
618 }
619 }
620
621 MatrixType&
operator =(const MatrixType & a)622 MatrixType::operator = (const MatrixType& a)
623 {
624 if (this != &a)
625 {
626 typ = a.typ;
627 sp_bandden = a.sp_bandden;
628 bandden = a.bandden;
629 upper_band = a.upper_band;
630 lower_band = a.lower_band;
631 dense = a.dense;
632 full = a.full;
633
634 if (nperm)
635 {
636 delete[] perm;
637 }
638
639 if (a.nperm != 0)
640 {
641 perm = new octave_idx_type [a.nperm];
642 for (octave_idx_type i = 0; i < a.nperm; i++)
643 perm[i] = a.perm[i];
644 }
645
646 nperm = a.nperm;
647 }
648
649 return *this;
650 }
651
652 int
type(bool quiet)653 MatrixType::type (bool quiet)
654 {
655 if (typ != MatrixType::Unknown
656 && (full || sp_bandden == octave_sparse_params::get_bandden ()))
657 {
658 if (! quiet && octave_sparse_params::get_key ("spumoni") != 0.)
659 warn_cached ();
660
661 return typ;
662 }
663
664 if (typ != MatrixType::Unknown
665 && octave_sparse_params::get_key ("spumoni") != 0.)
666 (*current_liboctave_warning_with_id_handler)
667 ("Octave:matrix-type-info", "invalidating matrix type");
668
669 typ = MatrixType::Unknown;
670
671 return typ;
672 }
673
674 int
type(const SparseMatrix & a)675 MatrixType::type (const SparseMatrix& a)
676 {
677 if (typ != MatrixType::Unknown
678 && (full || sp_bandden == octave_sparse_params::get_bandden ()))
679 {
680 if (octave_sparse_params::get_key ("spumoni") != 0.)
681 warn_cached ();
682
683 return typ;
684 }
685
686 MatrixType tmp_typ (a);
687 typ = tmp_typ.typ;
688 sp_bandden = tmp_typ.sp_bandden;
689 bandden = tmp_typ.bandden;
690 upper_band = tmp_typ.upper_band;
691 lower_band = tmp_typ.lower_band;
692 dense = tmp_typ.dense;
693 full = tmp_typ.full;
694 nperm = tmp_typ.nperm;
695
696 if (nperm != 0)
697 {
698 perm = new octave_idx_type [nperm];
699 for (octave_idx_type i = 0; i < nperm; i++)
700 perm[i] = tmp_typ.perm[i];
701 }
702
703 return typ;
704 }
705
706 int
type(const SparseComplexMatrix & a)707 MatrixType::type (const SparseComplexMatrix& a)
708 {
709 if (typ != MatrixType::Unknown
710 && (full || sp_bandden == octave_sparse_params::get_bandden ()))
711 {
712 if (octave_sparse_params::get_key ("spumoni") != 0.)
713 warn_cached ();
714
715 return typ;
716 }
717
718 MatrixType tmp_typ (a);
719 typ = tmp_typ.typ;
720 sp_bandden = tmp_typ.sp_bandden;
721 bandden = tmp_typ.bandden;
722 upper_band = tmp_typ.upper_band;
723 lower_band = tmp_typ.lower_band;
724 dense = tmp_typ.dense;
725 full = tmp_typ.full;
726 nperm = tmp_typ.nperm;
727
728 if (nperm != 0)
729 {
730 perm = new octave_idx_type [nperm];
731 for (octave_idx_type i = 0; i < nperm; i++)
732 perm[i] = tmp_typ.perm[i];
733 }
734
735 return typ;
736 }
737
738 int
type(const Matrix & a)739 MatrixType::type (const Matrix& a)
740 {
741 if (typ != MatrixType::Unknown)
742 {
743 if (octave_sparse_params::get_key ("spumoni") != 0.)
744 warn_cached ();
745
746 return typ;
747 }
748
749 MatrixType tmp_typ (a);
750 typ = tmp_typ.typ;
751 full = tmp_typ.full;
752 nperm = tmp_typ.nperm;
753
754 if (nperm != 0)
755 {
756 perm = new octave_idx_type [nperm];
757 for (octave_idx_type i = 0; i < nperm; i++)
758 perm[i] = tmp_typ.perm[i];
759 }
760
761 return typ;
762 }
763
764 int
type(const ComplexMatrix & a)765 MatrixType::type (const ComplexMatrix& a)
766 {
767 if (typ != MatrixType::Unknown)
768 {
769 if (octave_sparse_params::get_key ("spumoni") != 0.)
770 warn_cached ();
771
772 return typ;
773 }
774
775 MatrixType tmp_typ (a);
776 typ = tmp_typ.typ;
777 full = tmp_typ.full;
778 nperm = tmp_typ.nperm;
779
780 if (nperm != 0)
781 {
782 perm = new octave_idx_type [nperm];
783 for (octave_idx_type i = 0; i < nperm; i++)
784 perm[i] = tmp_typ.perm[i];
785 }
786
787 return typ;
788 }
789
790 int
type(const FloatMatrix & a)791 MatrixType::type (const FloatMatrix& a)
792 {
793 if (typ != MatrixType::Unknown)
794 {
795 if (octave_sparse_params::get_key ("spumoni") != 0.)
796 warn_cached ();
797
798 return typ;
799 }
800
801 MatrixType tmp_typ (a);
802 typ = tmp_typ.typ;
803 full = tmp_typ.full;
804 nperm = tmp_typ.nperm;
805
806 if (nperm != 0)
807 {
808 perm = new octave_idx_type [nperm];
809 for (octave_idx_type i = 0; i < nperm; i++)
810 perm[i] = tmp_typ.perm[i];
811 }
812
813 return typ;
814 }
815
816 int
type(const FloatComplexMatrix & a)817 MatrixType::type (const FloatComplexMatrix& a)
818 {
819 if (typ != MatrixType::Unknown)
820 {
821 if (octave_sparse_params::get_key ("spumoni") != 0.)
822 warn_cached ();
823
824 return typ;
825 }
826
827 MatrixType tmp_typ (a);
828 typ = tmp_typ.typ;
829 full = tmp_typ.full;
830 nperm = tmp_typ.nperm;
831
832 if (nperm != 0)
833 {
834 perm = new octave_idx_type [nperm];
835 for (octave_idx_type i = 0; i < nperm; i++)
836 perm[i] = tmp_typ.perm[i];
837 }
838
839 return typ;
840 }
841
842 void
info() const843 MatrixType::info () const
844 {
845 if (octave_sparse_params::get_key ("spumoni") != 0.)
846 {
847 if (typ == MatrixType::Unknown)
848 (*current_liboctave_warning_with_id_handler)
849 ("Octave:matrix-type-info", "unknown matrix type");
850 else if (typ == MatrixType::Diagonal)
851 (*current_liboctave_warning_with_id_handler)
852 ("Octave:matrix-type-info", "diagonal sparse matrix");
853 else if (typ == MatrixType::Permuted_Diagonal)
854 (*current_liboctave_warning_with_id_handler)
855 ("Octave:matrix-type-info", "permuted diagonal sparse matrix");
856 else if (typ == MatrixType::Upper)
857 (*current_liboctave_warning_with_id_handler)
858 ("Octave:matrix-type-info", "upper triangular matrix");
859 else if (typ == MatrixType::Lower)
860 (*current_liboctave_warning_with_id_handler)
861 ("Octave:matrix-type-info", "lower triangular matrix");
862 else if (typ == MatrixType::Permuted_Upper)
863 (*current_liboctave_warning_with_id_handler)
864 ("Octave:matrix-type-info", "permuted upper triangular matrix");
865 else if (typ == MatrixType::Permuted_Lower)
866 (*current_liboctave_warning_with_id_handler)
867 ("Octave:matrix-type-info", "permuted lower triangular Matrix");
868 else if (typ == MatrixType::Banded)
869 (*current_liboctave_warning_with_id_handler)
870 ("Octave:matrix-type-info",
871 "banded sparse matrix %" OCTAVE_IDX_TYPE_FORMAT "-1-"
872 "%" OCTAVE_IDX_TYPE_FORMAT " (density %f)",
873 lower_band, upper_band, bandden);
874 else if (typ == MatrixType::Banded_Hermitian)
875 (*current_liboctave_warning_with_id_handler)
876 ("Octave:matrix-type-info",
877 "banded hermitian/symmetric sparse matrix %" OCTAVE_IDX_TYPE_FORMAT
878 "-1-%" OCTAVE_IDX_TYPE_FORMAT " (density %f)",
879 lower_band, upper_band, bandden);
880 else if (typ == MatrixType::Hermitian)
881 (*current_liboctave_warning_with_id_handler)
882 ("Octave:matrix-type-info", "hermitian/symmetric matrix");
883 else if (typ == MatrixType::Tridiagonal)
884 (*current_liboctave_warning_with_id_handler)
885 ("Octave:matrix-type-info", "tridiagonal sparse matrix");
886 else if (typ == MatrixType::Tridiagonal_Hermitian)
887 (*current_liboctave_warning_with_id_handler)
888 ("Octave:matrix-type-info",
889 "hermitian/symmetric tridiagonal sparse matrix");
890 else if (typ == MatrixType::Rectangular)
891 (*current_liboctave_warning_with_id_handler)
892 ("Octave:matrix-type-info", "rectangular/singular matrix");
893 else if (typ == MatrixType::Full)
894 (*current_liboctave_warning_with_id_handler)
895 ("Octave:matrix-type-info", "full matrix");
896 }
897 }
898
899 void
mark_as_symmetric(void)900 MatrixType::mark_as_symmetric (void)
901 {
902 if (typ == MatrixType::Tridiagonal
903 || typ == MatrixType::Tridiagonal_Hermitian)
904 typ = MatrixType::Tridiagonal_Hermitian;
905 else if (typ == MatrixType::Banded || typ == MatrixType::Banded_Hermitian)
906 typ = MatrixType::Banded_Hermitian;
907 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian
908 || typ == MatrixType::Unknown)
909 typ = MatrixType::Hermitian;
910 else
911 (*current_liboctave_error_handler)
912 ("Can not mark current matrix type as symmetric");
913 }
914
915 void
mark_as_unsymmetric(void)916 MatrixType::mark_as_unsymmetric (void)
917 {
918 if (typ == MatrixType::Tridiagonal
919 || typ == MatrixType::Tridiagonal_Hermitian)
920 typ = MatrixType::Tridiagonal;
921 else if (typ == MatrixType::Banded || typ == MatrixType::Banded_Hermitian)
922 typ = MatrixType::Banded;
923 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian
924 || typ == MatrixType::Unknown)
925 typ = MatrixType::Full;
926 }
927
928 void
mark_as_permuted(const octave_idx_type np,const octave_idx_type * p)929 MatrixType::mark_as_permuted (const octave_idx_type np,
930 const octave_idx_type *p)
931 {
932 nperm = np;
933 perm = new octave_idx_type [nperm];
934 for (octave_idx_type i = 0; i < nperm; i++)
935 perm[i] = p[i];
936
937 if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal)
938 typ = MatrixType::Permuted_Diagonal;
939 else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
940 typ = MatrixType::Permuted_Upper;
941 else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
942 typ = MatrixType::Permuted_Lower;
943 else
944 (*current_liboctave_error_handler)
945 ("Can not mark current matrix type as symmetric");
946 }
947
948 void
mark_as_unpermuted(void)949 MatrixType::mark_as_unpermuted (void)
950 {
951 if (nperm)
952 {
953 nperm = 0;
954 delete [] perm;
955 }
956
957 if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal)
958 typ = MatrixType::Diagonal;
959 else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
960 typ = MatrixType::Upper;
961 else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
962 typ = MatrixType::Lower;
963 }
964
965 MatrixType
transpose(void) const966 MatrixType::transpose (void) const
967 {
968 MatrixType retval (*this);
969 if (typ == MatrixType::Upper)
970 retval.typ = MatrixType::Lower;
971 else if (typ == MatrixType::Permuted_Upper)
972 retval.typ = MatrixType::Permuted_Lower;
973 else if (typ == MatrixType::Lower)
974 retval.typ = MatrixType::Upper;
975 else if (typ == MatrixType::Permuted_Lower)
976 retval.typ = MatrixType::Permuted_Upper;
977 else if (typ == MatrixType::Banded)
978 {
979 retval.upper_band = lower_band;
980 retval.lower_band = upper_band;
981 }
982
983 return retval;
984 }
985
986 // Instantiate MatrixType template constructors that we need.
987
988 template MatrixType::MatrixType (const MSparse<double>&);
989 template MatrixType::MatrixType (const MSparse<Complex>&);
990