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