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