Bug Summary

File:liboctave/array/MatrixType.cc
Location:line 173, column 3
Description:Undefined or garbage value returned to caller

Annotated Source Code

1/*
2
3Copyright (C) 2006-2013 David Bateman
4Copyright (C) 2006 Andy Adler
5Copyright (C) 2009 VZLU Prague
6
7This file is part of Octave.
8
9Octave is free software; you can redistribute it and/or modify it
10under the terms of the GNU General Public License as published by the
11Free Software Foundation; either version 3 of the License, or (at your
12option) any later version.
13
14Octave is distributed in the hope that it will be useful, but WITHOUT
15ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
16FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
17for more details.
18
19You should have received a copy of the GNU General Public License
20along with Octave; see the file COPYING. If not, see
21<http://www.gnu.org/licenses/>.
22
23*/
24
25#ifdef HAVE_CONFIG_H1
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
41MatrixType::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
47MatrixType::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
60template<class T>
61MatrixType::matrix_type
62matrix_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)octave_local_buffer<T> _buffer_diag (ncols); T *diag = _buffer_diag;
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
117template<class T>
118MatrixType::matrix_type
119matrix_complex_probe (const MArray<std::complex<T> >& a)
120{
121 MatrixType::matrix_type typ;
4
'typ' declared without an initial value
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)
5
Taking true branch
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)octave_local_buffer<T> _buffer_diag (ncols); T *diag = _buffer_diag;
136
137 for (octave_idx_type j = 0;
7
Loop condition is true. Entering loop body
138 j < ncols && upper; j++)
6
Assuming 'j' is < 'ncols'
8
Assuming 'j' is >= 'ncols'
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;
9
Loop condition is false. Execution continues on line 161
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)
10
Taking false branch
162 typ = MatrixType::Upper;
163 else if (lower)
11
Taking false branch
164 typ = MatrixType::Lower;
165 else if (hermitian)
12
Taking false branch
166 typ = MatrixType::Hermitian;
167 else if (ncols == nrows)
13
Assuming 'ncols' is not equal to 'nrows'
14
Taking false branch
168 typ = MatrixType::Full;
169 }
170 else
171 typ = MatrixType::Rectangular;
172
173 return typ;
15
Undefined or garbage value returned to caller
174}
175
176MatrixType::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
184MatrixType::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
193MatrixType::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
201MatrixType::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);
3
Calling 'matrix_complex_probe'
207}
208
209MatrixType::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)octave_local_buffer<octave_idx_type> _buffer_tmp (nrows
); octave_idx_type *tmp = _buffer_tmp
;
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
528MatrixType::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)octave_local_buffer<octave_idx_type> _buffer_tmp (nrows
); octave_idx_type *tmp = _buffer_tmp
;
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}
847MatrixType::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
863MatrixType::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
883MatrixType::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
900MatrixType::~MatrixType (void)
901{
902 if (nperm != 0)
903 {
904 delete [] perm;
905 }
906}
907
908MatrixType&
909MatrixType::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
939int
940MatrixType::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
963int
964MatrixType::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
996int
997MatrixType::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
1029int
1030MatrixType::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
1056int
1057MatrixType::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
1083int
1084MatrixType::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
1110int
1111MatrixType::type (const FloatComplexMatrix &a)
1112{
1113 if (typ != MatrixType::Unknown)
1
Taking false branch
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);
2
Calling constructor for 'MatrixType'
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
1137void
1138MatrixType::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
1189void
1190MatrixType::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
1206void
1207MatrixType::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
1220void
1221MatrixType::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
1240void
1241MatrixType::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
1257MatrixType
1258MatrixType::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}