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