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