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