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