Mercurial > octave-nkf

comparison liboctave/SparseType.cc @ 5164:57077d0ddc8e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2005-02-25 19:55:24 by jwe]
author jwe
date Fri, 25 Feb 2005 19:55:28 +0000
parents
children 23b37da9fd5b
comparison
equal deleted inserted replaced
5163:9f3299378193 5164:57077d0ddc8e
1 /*
2
3 Copyright (C) 2004 David Bateman
4 Copyright (C) 1998-2004 Andy Adler
5
6 Octave is free software; you can redistribute it and/or modify it
7 under the terms of the GNU General Public License as published by the
8 Free Software Foundation; either version 2, or (at your option) any
9 later version.
10
11 Octave is distributed in the hope that it will be useful, but WITHOUT
12 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
14 for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with this program; see the file COPYING. If not, write to the Free
18 Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
19
20 */
21
22 #ifdef HAVE_CONFIG_H
23 #include <config.h>
24 #endif
25
26 #include "SparseType.h"
27 #include "dSparse.h"
28 #include "CSparse.h"
29 #include "oct-spparms.h"
30
31 // XXX FIXME XXX There is a large code duplication here
32
33 SparseType::SparseType (const SparseType &a) : typ (a.typ),
34 sp_bandden (a.sp_bandden), bandden (a.bandden),
35 upper_band (a.upper_band), lower_band (a.lower_band),
36 dense (a.dense), nperm (a.nperm)
37 {
38 if (nperm != 0)
39 {
40 row_perm = new int [nperm];
41 col_perm = new int [nperm];
42 for (int i = 0; i < nperm; i++)
43 {
44 row_perm[i] = a.row_perm[i];
45 col_perm[i] = a.col_perm[i];
46 }
47 }
48 }
49
50 SparseType::SparseType (const SparseMatrix &a)
51 {
52 int nrows = a.rows ();
53 int ncols = a.cols ();
54 int nnz = a.nnz ();
55
56 nperm = 0;
57
58 if (nrows != ncols)
59 typ = SparseType::Rectangular;
60 else
61 {
62 sp_bandden = Voctave_sparse_controls.get_key ("bandden");
63 bool maybe_hermitian = false;
64 typ = SparseType::Full;
65
66 if (nnz == ncols)
67 {
68 matrix_type tmp_typ = SparseType::Diagonal;
69 int i;
70 // Maybe the matrix is diagonal
71 for (i = 0; i < ncols; i++)
72 {
73 if (a.cidx(i+1) != a.cidx(i) + 1)
74 {
75 tmp_typ = Full;
76 break;
77 }
78 if (a.ridx(i) != i)
79 {
80 tmp_typ = SparseType::Permuted_Diagonal;
81 break;
82 }
83 }
84
85 if (tmp_typ == SparseType::Permuted_Diagonal)
86 {
87 bool found [ncols];
88
89 for (int j = 0; j < i; j++)
90 found [j] = true;
91 for (int j = i; j < ncols; j++)
92 found [j] = false;
93
94 for (int j = i; j < ncols; j++)
95 {
96 if ((a.cidx(j+1) != a.cidx(j) + 1) || found [a.ridx(j)])
97 {
98 tmp_typ = Full;
99 break;
100 }
101 found [a.ridx(j)] = true;
102 }
103 }
104 typ = tmp_typ;
105 }
106
107 if (typ == Full)
108 {
109 // Search for banded, upper and lower triangular matrices
110 bool singular = false;
111 upper_band = 0;
112 lower_band = 0;
113 for (int j = 0; j < ncols; j++)
114 {
115 bool zero_on_diagonal = true;
116 for (int i = a.cidx(j); i < a.cidx(j+1); i++)
117 if (a.ridx(i) == j)
118 {
119 zero_on_diagonal = false;
120 break;
121 }
122
123 if (zero_on_diagonal)
124 {
125 singular = true;
126 break;
127 }
128
129 if (a.cidx(j+1) - a.cidx(j) > 0)
130 {
131 int ru = a.ridx(a.cidx(j));
132 int rl = a.ridx(a.cidx(j+1)-1);
133
134 if (j - ru > upper_band)
135 upper_band = j - ru;
136
137 if (rl - j > lower_band)
138 lower_band = rl - j;
139 }
140 }
141
142 if (!singular)
143 {
144 bandden = double (nnz) /
145 (double (ncols) * (double (lower_band) +
146 double (upper_band)) -
147 0.5 * double (upper_band + 1) * double (upper_band) -
148 0.5 * double (lower_band + 1) * double (lower_band));
149
150 if (sp_bandden != 1. && bandden > sp_bandden)
151 {
152 if (upper_band == 1 && lower_band == 1)
153 typ = SparseType::Tridiagonal;
154 else
155 typ = SparseType::Banded;
156
157 int nnz_in_band = (upper_band + lower_band + 1) * nrows -
158 (1 + upper_band) * upper_band / 2 -
159 (1 + lower_band) * lower_band / 2;
160 if (nnz_in_band == nnz)
161 dense = true;
162 else
163 dense = false;
164 }
165 else if (upper_band == 0)
166 typ = SparseType::Lower;
167 else if (lower_band == 0)
168 typ = SparseType::Upper;
169
170 if (upper_band == lower_band)
171 maybe_hermitian = true;
172 }
173
174 if (typ == Full)
175 {
176 // Search for a permuted triangular matrix, and test if
177 // permutation is singular
178
179 // XXX FIXME XXX Write this test based on dmperm
180
181 }
182 }
183
184 if (maybe_hermitian && (typ == Full || typ == Tridiagonal ||
185 typ == Banded))
186 {
187 // Check for symmetry, with positive real diagonal, which
188 // has a very good chance of being symmetric positive
189 // definite..
190 bool is_herm = true;
191
192 for (int j = 0; j < ncols; j++)
193 {
194 bool diag_positive = false;
195
196 for (int i = a.cidx(j); i < a.cidx(j+1); i++)
197 {
198 int ri = a.ridx(i);
199
200 if (ri == j)
201 {
202 if (a.data(i) == std::abs(a.data(i)))
203 diag_positive = true;
204 else
205 break;
206 }
207 else
208 {
209 bool found = false;
210
211 for (int k = a.cidx(ri); k < a.cidx(ri+1); k++)
212 {
213 if (a.ridx(k) == j)
214 {
215 if (a.data(i) == conj (a.data(k)))
216 found = true;
217 break;
218 }
219 }
220
221 if (! found)
222 {
223 is_herm = false;
224 break;
225 }
226 }
227 }
228
229 if (! diag_positive || ! is_herm)
230 {
231 is_herm = false;
232 break;
233 }
234 }
235
236 if (is_herm)
237 {
238 if (typ == Full)
239 typ = Hermitian;
240 else if (typ == Banded)
241 typ = Banded_Hermitian;
242 else
243 typ = Tridiagonal_Hermitian;
244 }
245 }
246 }
247 }
248
249 SparseType::SparseType (const SparseComplexMatrix &a)
250 {
251 int nrows = a.rows ();
252 int ncols = a.cols ();
253 int nnz = a.nnz ();
254
255 nperm = 0;
256
257 if (nrows != ncols)
258 typ = SparseType::Rectangular;
259 else
260 {
261 sp_bandden = Voctave_sparse_controls.get_key ("bandden");
262 bool maybe_hermitian = false;
263 typ = SparseType::Full;
264
265 if (nnz == ncols)
266 {
267 matrix_type tmp_typ = SparseType::Diagonal;
268 int i;
269 // Maybe the matrix is diagonal
270 for (i = 0; i < ncols; i++)
271 {
272 if (a.cidx(i+1) != a.cidx(i) + 1)
273 {
274 tmp_typ = Full;
275 break;
276 }
277 if (a.ridx(i) != i)
278 {
279 tmp_typ = SparseType::Permuted_Diagonal;
280 break;
281 }
282 }
283
284 if (tmp_typ == SparseType::Permuted_Diagonal)
285 {
286 bool found [ncols];
287
288 for (int j = 0; j < i; j++)
289 found [j] = true;
290 for (int j = i; j < ncols; j++)
291 found [j] = false;
292
293 for (int j = i; j < ncols; j++)
294 {
295 if ((a.cidx(j+1) != a.cidx(j) + 1) || found [a.ridx(j)])
296 {
297 tmp_typ = Full;
298 break;
299 }
300 found [a.ridx(j)] = true;
301 }
302 }
303 typ = tmp_typ;
304 }
305
306 if (typ == Full)
307 {
308 // Search for banded, upper and lower triangular matrices
309 bool singular = false;
310 upper_band = 0;
311 lower_band = 0;
312 for (int j = 0; j < ncols; j++)
313 {
314 bool zero_on_diagonal = true;
315 for (int i = a.cidx(j); i < a.cidx(j+1); i++)
316 if (a.ridx(i) == j)
317 {
318 zero_on_diagonal = false;
319 break;
320 }
321
322 if (zero_on_diagonal)
323 {
324 singular = true;
325 break;
326 }
327
328 if (a.cidx(j+1) - a.cidx(j) > 0)
329 {
330 int ru = a.ridx(a.cidx(j));
331 int rl = a.ridx(a.cidx(j+1)-1);
332
333 if (j - ru > upper_band)
334 upper_band = j - ru;
335
336 if (rl - j > lower_band)
337 lower_band = rl - j;
338 }
339 }
340
341 if (!singular)
342 {
343 bandden = double (nnz) /
344 (double (ncols) * (double (lower_band) +
345 double (upper_band)) -
346 0.5 * double (upper_band + 1) * double (upper_band) -
347 0.5 * double (lower_band + 1) * double (lower_band));
348
349 if (sp_bandden != 1. && bandden > sp_bandden)
350 {
351 if (upper_band == 1 && lower_band == 1)
352 typ = SparseType::Tridiagonal;
353 else
354 typ = SparseType::Banded;
355
356 int nnz_in_band = (upper_band + lower_band + 1) * nrows -
357 (1 + upper_band) * upper_band / 2 -
358 (1 + lower_band) * lower_band / 2;
359 if (nnz_in_band == nnz)
360 dense = true;
361 else
362 dense = false;
363 }
364 else if (upper_band == 0)
365 typ = SparseType::Lower;
366 else if (lower_band == 0)
367 typ = SparseType::Upper;
368
369 if (upper_band == lower_band)
370 maybe_hermitian = true;
371 }
372
373 if (typ == Full)
374 {
375 // Search for a permuted triangular matrix, and test if
376 // permutation is singular
377
378 // XXX FIXME XXX Write this test based on dmperm
379
380 }
381 }
382
383 if (maybe_hermitian && (typ == Full || typ == Tridiagonal ||
384 typ == Banded))
385 {
386 // Check for symmetry, with positive real diagonal, which
387 // has a very good chance of being symmetric positive
388 // definite..
389 bool is_herm = true;
390
391 for (int j = 0; j < ncols; j++)
392 {
393 bool diag_positive = false;
394
395 for (int i = a.cidx(j); i < a.cidx(j+1); i++)
396 {
397 int ri = a.ridx(i);
398
399 if (ri == j)
400 {
401 if (a.data(i) == std::abs(a.data(i)))
402 diag_positive = true;
403 else
404 break;
405 }
406 else
407 {
408 bool found = false;
409
410 for (int k = a.cidx(ri); k < a.cidx(ri+1); k++)
411 {
412 if (a.ridx(k) == j)
413 {
414 if (a.data(i) == a.data(k))
415 found = true;
416 break;
417 }
418 }
419
420 if (! found)
421 {
422 is_herm = false;
423 break;
424 }
425 }
426 }
427
428 if (! diag_positive || ! is_herm)
429 {
430 is_herm = false;
431 break;
432 }
433 }
434
435 if (is_herm)
436 {
437 if (typ == Full)
438 typ = Hermitian;
439 else if (typ == Banded)
440 typ = Banded_Hermitian;
441 else
442 typ = Tridiagonal_Hermitian;
443 }
444 }
445 }
446 }
447
448 SparseType::~SparseType (void)
449 {
450 if (nperm != 0)
451 {
452 delete [] row_perm;
453 delete [] col_perm;
454 }
455 }
456
457 SparseType&
458 SparseType::operator = (const SparseType& a)
459 {
460 if (this != &a)
461 {
462 typ = a.typ;
463 sp_bandden = a.sp_bandden;
464 bandden = a.bandden;
465 upper_band = a.upper_band;
466 lower_band = a.lower_band;
467 dense = a.dense;
468 nperm = a.nperm;
469
470 if (nperm != 0)
471 {
472 row_perm = new int [nperm];
473 col_perm = new int [nperm];
474 for (int i = 0; i < nperm; i++)
475 {
476 row_perm[i] = a.row_perm[i];
477 col_perm[i] = a.col_perm[i];
478 }
479 }
480
481 }
482 return *this;
483 }
484
485 int
486 SparseType::type (const SparseMatrix &a)
487 {
488 if (typ != SparseType::Unknown &&
489 sp_bandden == Voctave_sparse_controls.get_key ("bandden"))
490 {
491 if (Voctave_sparse_controls.get_key ("spumoni") != 0.)
492 (*current_liboctave_warning_handler)
493 ("Using Cached Sparse Matrix Type");
494
495 return typ;
496 }
497
498 if (Voctave_sparse_controls.get_key ("spumoni") != 0.)
499 (*current_liboctave_warning_handler)
500 ("Calculating Sparse Matrix Type");
501
502
503 SparseType tmp_typ (a);
504 typ = tmp_typ.typ;
505 sp_bandden = tmp_typ.sp_bandden;
506 bandden = tmp_typ.bandden;
507 upper_band = tmp_typ.upper_band;
508 lower_band = tmp_typ.lower_band;
509 dense = tmp_typ.dense;
510 nperm = tmp_typ.nperm;
511
512 if (nperm != 0)
513 {
514 row_perm = new int [nperm];
515 col_perm = new int [nperm];
516 for (int i = 0; i < nperm; i++)
517 {
518 row_perm[i] = tmp_typ.row_perm[i];
519 col_perm[i] = tmp_typ.col_perm[i];
520 }
521 }
522
523 return typ;
524 }
525
526 int
527 SparseType::type (const SparseComplexMatrix &a)
528 {
529 if (typ != SparseType::Unknown &&
530 sp_bandden == Voctave_sparse_controls.get_key ("bandden"))
531 {
532 if (Voctave_sparse_controls.get_key ("spumoni") != 0.)
533 (*current_liboctave_warning_handler)
534 ("Using Cached Sparse Matrix Type");
535
536 return typ;
537 }
538
539 if (Voctave_sparse_controls.get_key ("spumoni") != 0.)
540 (*current_liboctave_warning_handler)
541 ("Calculating Sparse Matrix Type");
542
543
544 SparseType tmp_typ (a);
545 typ = tmp_typ.typ;
546 sp_bandden = tmp_typ.sp_bandden;
547 bandden = tmp_typ.bandden;
548 upper_band = tmp_typ.upper_band;
549 lower_band = tmp_typ.lower_band;
550 dense = tmp_typ.dense;
551 nperm = tmp_typ.nperm;
552
553 if (nperm != 0)
554 {
555 row_perm = new int [nperm];
556 col_perm = new int [nperm];
557 for (int i = 0; i < nperm; i++)
558 {
559 row_perm[i] = tmp_typ.row_perm[i];
560 col_perm[i] = tmp_typ.col_perm[i];
561 }
562 }
563
564 return typ;
565 }
566
567 void
568 SparseType::info (void) const
569 {
570 if (Voctave_sparse_controls.get_key ("spumoni") != 0.)
571 {
572 if (typ == SparseType::Unknown)
573 (*current_liboctave_warning_handler)
574 ("Unknown Sparse Matrix Type");
575 else if (typ == SparseType::Diagonal)
576 (*current_liboctave_warning_handler)
577 ("Diagonal Sparse Matrix");
578 else if (typ == SparseType::Permuted_Diagonal)
579 (*current_liboctave_warning_handler)
580 ("Permuted Diagonal Sparse Matrix");
581 else if (typ == SparseType::Upper)
582 (*current_liboctave_warning_handler)
583 ("Upper Triangular Sparse Matrix");
584 else if (typ == SparseType::Lower)
585 (*current_liboctave_warning_handler)
586 ("Lower Triangular Sparse Matrix");
587 else if (typ == SparseType::Permuted_Upper)
588 (*current_liboctave_warning_handler)
589 ("Permuted Upper Triangular Sparse Matrix");
590 else if (typ == SparseType::Permuted_Lower)
591 (*current_liboctave_warning_handler)
592 ("Permuted Lower Triangular Sparse Matrix");
593 else if (typ == SparseType::Banded)
594 (*current_liboctave_warning_handler)
595 ("Banded Sparse Matrix %g-1-%g (Density %g)", lower_band,
596 upper_band, bandden);
597 else if (typ == SparseType::Banded_Hermitian)
598 (*current_liboctave_warning_handler)
599 ("Banded Hermitian/Symmetric Sparse Matrix %g-1-%g (Density %g)",
600 lower_band, upper_band, bandden);
601 else if (typ == SparseType::Hermitian)
602 (*current_liboctave_warning_handler)
603 ("Hermitian/Symmetric Sparse Matrix");
604 else if (typ == SparseType::Tridiagonal)
605 (*current_liboctave_warning_handler)
606 ("Tridiagonal Sparse Matrix");
607 else if (typ == SparseType::Tridiagonal_Hermitian)
608 (*current_liboctave_warning_handler)
609 ("Hermitian/Symmetric Tridiagonal Sparse Matrix");
610 else if (typ == SparseType::Rectangular)
611 (*current_liboctave_warning_handler)
612 ("Rectangular Sparse Matrix");
613 else if (typ == SparseType::Full)
614 (*current_liboctave_warning_handler)
615 ("Full Sparse Matrix");
616 }
617 }
618
619 void
620 SparseType::mark_as_symmetric (void)
621 {
622 if (typ == SparseType::Tridiagonal ||
623 typ == SparseType::Tridiagonal_Hermitian)
624 typ = SparseType::Tridiagonal_Hermitian;
625 else if (typ == SparseType::Banded ||
626 typ == SparseType::Banded_Hermitian)
627 typ = SparseType::Banded_Hermitian;
628 else if (typ == SparseType::Full || typ == SparseType::Hermitian ||
629 typ == SparseType::Unknown)
630 typ = SparseType::Hermitian;
631 else
632 (*current_liboctave_error_handler)
633 ("Can not mark current matrix type as symmetric");
634 }
635
636 void
637 SparseType::mark_as_unsymmetric (void)
638 {
639 if (typ == SparseType::Tridiagonal ||
640 typ == SparseType::Tridiagonal_Hermitian)
641 typ = SparseType::Tridiagonal;
642 else if (typ == SparseType::Banded ||
643 typ == SparseType::Banded_Hermitian)
644 typ = SparseType::Banded;
645 else if (typ == SparseType::Full || typ == SparseType::Hermitian ||
646 typ == SparseType::Unknown)
647 typ = SparseType::Full;
648 }
649
650 void
651 SparseType::mark_as_permuted (const int np, const int *pr, const int *pc)
652 {
653 nperm = np;
654 row_perm = new int [nperm];
655 col_perm = new int [nperm];
656 for (int i = 0; i < nperm; i++)
657 {
658 row_perm[i] = pr[i];
659 col_perm[i] = pc[i];
660 }
661
662 if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal)
663 typ = SparseType::Permuted_Diagonal;
664 else if (typ == SparseType::Upper || typ == SparseType::Permuted_Upper)
665 typ = SparseType::Permuted_Upper;
666 else if (typ == SparseType::Lower || typ == SparseType::Permuted_Lower)
667 typ = SparseType::Permuted_Lower;
668 else
669 (*current_liboctave_error_handler)
670 ("Can not mark current matrix type as symmetric");
671 }
672
673 void
674 SparseType::mark_as_unpermuted (void)
675 {
676 if (nperm)
677 {
678 nperm = 0;
679 delete [] row_perm;
680 delete [] col_perm;
681 }
682
683 if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal)
684 typ = SparseType::Diagonal;
685 else if (typ == SparseType::Upper || typ == SparseType::Permuted_Upper)
686 typ = SparseType::Upper;
687 else if (typ == SparseType::Lower || typ == SparseType::Permuted_Lower)
688 typ = SparseType::Lower;
689 }
690
691 /*
692 ;;; Local Variables: ***
693 ;;; mode: C++ ***
694 ;;; End: ***
695 */
696