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