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