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