1 | |
2 | |
3 | |
4 | |
5 | |
6 | |
7 | |
8 | |
9 | |
10 | |
11 | |
12 | |
13 | |
14 | |
15 | |
16 | |
17 | |
18 | |
19 | |
20 | |
21 | |
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 | |
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 | |
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 | |
117 | template<class T> |
118 | MatrixType::matrix_type |
119 | matrix_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 | |
127 | |
128 | if (ncols == nrows) |
| |
129 | { |
130 | bool upper = true; |
131 | bool lower = true; |
132 | bool hermitian = true; |
133 | |
134 | |
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) |
| |
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) |
| 13 | | Assuming 'ncols' is not equal to 'nrows' | |
|
| |
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 | |
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); |
| 3 | | Calling 'matrix_complex_probe' | |
|
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 | |
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 | |
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 | |
344 | |
345 | |
346 | |
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 | |
450 | |
451 | |
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 | |
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 | |
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 | |
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 | |
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 | |
663 | |
664 | |
665 | |
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 | |
769 | |
770 | |
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 | |
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 | |
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); |
| 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 | |
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 | } |