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