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