Mercurial > hg > octave-nkf
annotate liboctave/CMatrix.cc @ 10182:0522a65bcd56
assume unistd.h and sys/types.h exist
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Thu, 21 Jan 2010 15:41:19 -0500 |
parents | 4c0cdbe0acca |
children | f7ba6cfe7fb7 |
rev | line source |
---|---|
1993 | 1 // Matrix manipulations. |
458 | 2 /* |
3 | |
7017 | 4 Copyright (C) 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002, |
5 2003, 2004, 2005, 2006, 2007 John W. Eaton | |
8920 | 6 Copyright (C) 2008, 2009 Jaroslav Hajek |
9601
a9b37bae1802
add a couple of missing copyright statements
Jaroslav Hajek <highegg@gmail.com>
parents:
9578
diff
changeset
|
7 Copyright (C) 2009 VZLU Prague, a.s. |
458 | 8 |
9 This file is part of Octave. | |
10 | |
11 Octave is free software; you can redistribute it and/or modify it | |
12 under the terms of the GNU General Public License as published by the | |
7016 | 13 Free Software Foundation; either version 3 of the License, or (at your |
14 option) any later version. | |
458 | 15 |
16 Octave is distributed in the hope that it will be useful, but WITHOUT | |
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
18 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
19 for more details. | |
20 | |
21 You should have received a copy of the GNU General Public License | |
7016 | 22 along with Octave; see the file COPYING. If not, see |
23 <http://www.gnu.org/licenses/>. | |
458 | 24 |
25 */ | |
26 | |
27 #ifdef HAVE_CONFIG_H | |
1192 | 28 #include <config.h> |
458 | 29 #endif |
30 | |
1367 | 31 #include <cfloat> |
32 | |
3503 | 33 #include <iostream> |
6209 | 34 #include <vector> |
1367 | 35 |
5775 | 36 // FIXME |
2443 | 37 #include <sys/types.h> |
458 | 38 |
4669 | 39 #include "Array-util.h" |
2828 | 40 #include "CMatrix.h" |
1819 | 41 #include "CmplxAEPBAL.h" |
9523
0ce82753dd72
more configure changes for libraries
John W. Eaton <jwe@octave.org>
parents:
9469
diff
changeset
|
42 #include "CmplxCHOL.h" |
1819 | 43 #include "CmplxSCHUR.h" |
740 | 44 #include "CmplxSVD.h" |
9523
0ce82753dd72
more configure changes for libraries
John W. Eaton <jwe@octave.org>
parents:
9469
diff
changeset
|
45 #include "DET.h" |
1847 | 46 #include "f77-fcn.h" |
7503
8c32f95c2639
convert mapper functions to new format
David Bateman <dbateman@free.fr>
parents:
7488
diff
changeset
|
47 #include "functor.h" |
458 | 48 #include "lo-error.h" |
2354 | 49 #include "lo-ieee.h" |
50 #include "lo-mappers.h" | |
1968 | 51 #include "lo-utils.h" |
1367 | 52 #include "mx-base.h" |
2828 | 53 #include "mx-cm-dm.h" |
9523
0ce82753dd72
more configure changes for libraries
John W. Eaton <jwe@octave.org>
parents:
9469
diff
changeset
|
54 #include "mx-cm-s.h" |
3176 | 55 #include "mx-dm-cm.h" |
1367 | 56 #include "mx-inlines.cc" |
8774
b756ce0002db
split implementation and interface in mx-op-defs and MArray-defs
Jaroslav Hajek <highegg@gmail.com>
parents:
8743
diff
changeset
|
57 #include "mx-op-defs.h" |
1650 | 58 #include "oct-cmplx.h" |
9523
0ce82753dd72
more configure changes for libraries
John W. Eaton <jwe@octave.org>
parents:
9469
diff
changeset
|
59 #include "oct-fftw.h" |
0ce82753dd72
more configure changes for libraries
John W. Eaton <jwe@octave.org>
parents:
9469
diff
changeset
|
60 #include "oct-locbuf.h" |
8336
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
61 #include "oct-norm.h" |
458 | 62 |
63 // Fortran functions we call. | |
64 | |
65 extern "C" | |
66 { | |
7477 | 67 F77_RET_T |
68 F77_FUNC (xilaenv, XILAENV) (const octave_idx_type&, F77_CONST_CHAR_ARG_DECL, | |
69 F77_CONST_CHAR_ARG_DECL, | |
70 const octave_idx_type&, const octave_idx_type&, | |
7478 | 71 const octave_idx_type&, const octave_idx_type&, |
72 octave_idx_type& | |
73 F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL); | |
7476 | 74 |
4552 | 75 F77_RET_T |
76 F77_FUNC (zgebal, ZGEBAL) (F77_CONST_CHAR_ARG_DECL, | |
5275 | 77 const octave_idx_type&, Complex*, const octave_idx_type&, octave_idx_type&, |
78 octave_idx_type&, double*, octave_idx_type& | |
4552 | 79 F77_CHAR_ARG_LEN_DECL); |
80 | |
81 F77_RET_T | |
82 F77_FUNC (dgebak, DGEBAK) (F77_CONST_CHAR_ARG_DECL, | |
83 F77_CONST_CHAR_ARG_DECL, | |
5275 | 84 const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, double*, |
85 const octave_idx_type&, double*, const octave_idx_type&, octave_idx_type& | |
4552 | 86 F77_CHAR_ARG_LEN_DECL |
87 F77_CHAR_ARG_LEN_DECL); | |
88 | |
89 F77_RET_T | |
90 F77_FUNC (zgemm, ZGEMM) (F77_CONST_CHAR_ARG_DECL, | |
91 F77_CONST_CHAR_ARG_DECL, | |
5275 | 92 const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, |
93 const Complex&, const Complex*, const octave_idx_type&, | |
94 const Complex*, const octave_idx_type&, const Complex&, | |
95 Complex*, const octave_idx_type& | |
4552 | 96 F77_CHAR_ARG_LEN_DECL |
97 F77_CHAR_ARG_LEN_DECL); | |
98 | |
99 F77_RET_T | |
5983 | 100 F77_FUNC (zgemv, ZGEMV) (F77_CONST_CHAR_ARG_DECL, |
101 const octave_idx_type&, const octave_idx_type&, const Complex&, | |
102 const Complex*, const octave_idx_type&, const Complex*, | |
103 const octave_idx_type&, const Complex&, Complex*, const octave_idx_type& | |
104 F77_CHAR_ARG_LEN_DECL); | |
105 | |
106 F77_RET_T | |
107 F77_FUNC (xzdotu, XZDOTU) (const octave_idx_type&, const Complex*, const octave_idx_type&, | |
108 const Complex*, const octave_idx_type&, Complex&); | |
109 | |
110 F77_RET_T | |
7800
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
111 F77_FUNC (xzdotc, XZDOTC) (const octave_idx_type&, const Complex*, const octave_idx_type&, |
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
112 const Complex*, const octave_idx_type&, Complex&); |
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
113 |
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
114 F77_RET_T |
7801
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
115 F77_FUNC (zsyrk, ZSYRK) (F77_CONST_CHAR_ARG_DECL, |
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
116 F77_CONST_CHAR_ARG_DECL, |
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
117 const octave_idx_type&, const octave_idx_type&, |
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
118 const Complex&, const Complex*, const octave_idx_type&, |
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
119 const Complex&, Complex*, const octave_idx_type& |
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
120 F77_CHAR_ARG_LEN_DECL |
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
121 F77_CHAR_ARG_LEN_DECL); |
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
122 |
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
123 F77_RET_T |
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
124 F77_FUNC (zherk, ZHERK) (F77_CONST_CHAR_ARG_DECL, |
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
125 F77_CONST_CHAR_ARG_DECL, |
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
126 const octave_idx_type&, const octave_idx_type&, |
10105
018d92ad0730
fix prototypes for zherk and cherk
John W. Eaton <jwe@octave.org>
parents:
9825
diff
changeset
|
127 const double&, const Complex*, const octave_idx_type&, |
018d92ad0730
fix prototypes for zherk and cherk
John W. Eaton <jwe@octave.org>
parents:
9825
diff
changeset
|
128 const double&, Complex*, const octave_idx_type& |
7801
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
129 F77_CHAR_ARG_LEN_DECL |
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
130 F77_CHAR_ARG_LEN_DECL); |
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
131 |
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
132 F77_RET_T |
5275 | 133 F77_FUNC (zgetrf, ZGETRF) (const octave_idx_type&, const octave_idx_type&, Complex*, const octave_idx_type&, |
134 octave_idx_type*, octave_idx_type&); | |
4552 | 135 |
136 F77_RET_T | |
137 F77_FUNC (zgetrs, ZGETRS) (F77_CONST_CHAR_ARG_DECL, | |
5275 | 138 const octave_idx_type&, const octave_idx_type&, Complex*, const octave_idx_type&, |
139 const octave_idx_type*, Complex*, const octave_idx_type&, octave_idx_type& | |
4552 | 140 F77_CHAR_ARG_LEN_DECL); |
141 | |
142 F77_RET_T | |
5275 | 143 F77_FUNC (zgetri, ZGETRI) (const octave_idx_type&, Complex*, const octave_idx_type&, const octave_idx_type*, |
144 Complex*, const octave_idx_type&, octave_idx_type&); | |
4552 | 145 |
146 F77_RET_T | |
147 F77_FUNC (zgecon, ZGECON) (F77_CONST_CHAR_ARG_DECL, | |
5275 | 148 const octave_idx_type&, Complex*, |
149 const octave_idx_type&, const double&, double&, | |
150 Complex*, double*, octave_idx_type& | |
4552 | 151 F77_CHAR_ARG_LEN_DECL); |
152 | |
153 F77_RET_T | |
7072 | 154 F77_FUNC (zgelsy, ZGELSY) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, |
155 Complex*, const octave_idx_type&, Complex*, | |
156 const octave_idx_type&, octave_idx_type*, double&, octave_idx_type&, | |
157 Complex*, const octave_idx_type&, double*, octave_idx_type&); | |
158 | |
159 F77_RET_T | |
160 F77_FUNC (zgelsd, ZGELSD) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, | |
5275 | 161 Complex*, const octave_idx_type&, Complex*, |
7071 | 162 const octave_idx_type&, double*, double&, octave_idx_type&, |
7072 | 163 Complex*, const octave_idx_type&, double*, |
164 octave_idx_type*, octave_idx_type&); | |
458 | 165 |
5785 | 166 F77_RET_T |
167 F77_FUNC (zpotrf, ZPOTRF) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, | |
168 Complex*, const octave_idx_type&, | |
169 octave_idx_type& F77_CHAR_ARG_LEN_DECL); | |
170 | |
171 F77_RET_T | |
172 F77_FUNC (zpocon, ZPOCON) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, | |
173 Complex*, const octave_idx_type&, const double&, | |
174 double&, Complex*, double*, | |
175 octave_idx_type& F77_CHAR_ARG_LEN_DECL); | |
176 | |
177 F77_RET_T | |
178 F77_FUNC (zpotrs, ZPOTRS) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, | |
179 const octave_idx_type&, const Complex*, | |
180 const octave_idx_type&, Complex*, | |
181 const octave_idx_type&, octave_idx_type& | |
182 F77_CHAR_ARG_LEN_DECL); | |
183 | |
184 F77_RET_T | |
6207 | 185 F77_FUNC (ztrtri, ZTRTRI) (F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, |
186 const octave_idx_type&, const Complex*, | |
187 const octave_idx_type&, octave_idx_type& | |
188 F77_CHAR_ARG_LEN_DECL | |
189 F77_CHAR_ARG_LEN_DECL); | |
190 | |
191 F77_RET_T | |
5785 | 192 F77_FUNC (ztrcon, ZTRCON) (F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, |
193 F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, | |
194 const Complex*, const octave_idx_type&, double&, | |
195 Complex*, double*, octave_idx_type& | |
196 F77_CHAR_ARG_LEN_DECL | |
197 F77_CHAR_ARG_LEN_DECL | |
198 F77_CHAR_ARG_LEN_DECL); | |
199 | |
200 F77_RET_T | |
201 F77_FUNC (ztrtrs, ZTRTRS) (F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, | |
202 F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, | |
203 const octave_idx_type&, const Complex*, | |
204 const octave_idx_type&, Complex*, | |
205 const octave_idx_type&, octave_idx_type& | |
206 F77_CHAR_ARG_LEN_DECL | |
207 F77_CHAR_ARG_LEN_DECL | |
208 F77_CHAR_ARG_LEN_DECL); | |
209 | |
4552 | 210 F77_RET_T |
211 F77_FUNC (zlartg, ZLARTG) (const Complex&, const Complex&, | |
212 double&, Complex&, Complex&); | |
213 | |
214 F77_RET_T | |
215 F77_FUNC (ztrsyl, ZTRSYL) (F77_CONST_CHAR_ARG_DECL, | |
216 F77_CONST_CHAR_ARG_DECL, | |
5275 | 217 const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, |
218 const Complex*, const octave_idx_type&, | |
219 const Complex*, const octave_idx_type&, | |
220 const Complex*, const octave_idx_type&, double&, octave_idx_type& | |
4552 | 221 F77_CHAR_ARG_LEN_DECL |
222 F77_CHAR_ARG_LEN_DECL); | |
223 | |
224 F77_RET_T | |
225 F77_FUNC (xzlange, XZLANGE) (F77_CONST_CHAR_ARG_DECL, | |
5275 | 226 const octave_idx_type&, const octave_idx_type&, const Complex*, |
227 const octave_idx_type&, double*, double& | |
4552 | 228 F77_CHAR_ARG_LEN_DECL); |
458 | 229 } |
230 | |
2354 | 231 static const Complex Complex_NaN_result (octave_NaN, octave_NaN); |
232 | |
1360 | 233 // Complex Matrix class |
458 | 234 |
235 ComplexMatrix::ComplexMatrix (const Matrix& a) | |
1214 | 236 : MArray2<Complex> (a.rows (), a.cols ()) |
458 | 237 { |
5275 | 238 for (octave_idx_type j = 0; j < cols (); j++) |
239 for (octave_idx_type i = 0; i < rows (); i++) | |
458 | 240 elem (i, j) = a.elem (i, j); |
241 } | |
242 | |
2349 | 243 ComplexMatrix::ComplexMatrix (const RowVector& rv) |
244 : MArray2<Complex> (1, rv.length (), 0.0) | |
245 { | |
5275 | 246 for (octave_idx_type i = 0; i < rv.length (); i++) |
2349 | 247 elem (0, i) = rv.elem (i); |
248 } | |
249 | |
250 ComplexMatrix::ComplexMatrix (const ColumnVector& cv) | |
251 : MArray2<Complex> (cv.length (), 1, 0.0) | |
252 { | |
5275 | 253 for (octave_idx_type i = 0; i < cv.length (); i++) |
2349 | 254 elem (i, 0) = cv.elem (i); |
255 } | |
256 | |
458 | 257 ComplexMatrix::ComplexMatrix (const DiagMatrix& a) |
1214 | 258 : MArray2<Complex> (a.rows (), a.cols (), 0.0) |
458 | 259 { |
5275 | 260 for (octave_idx_type i = 0; i < a.length (); i++) |
458 | 261 elem (i, i) = a.elem (i, i); |
262 } | |
263 | |
2349 | 264 ComplexMatrix::ComplexMatrix (const ComplexRowVector& rv) |
8614
5114ea5a41b5
use shallow copying in Matrix/RowVector/ColumnVector conversions
Jaroslav Hajek <highegg@gmail.com>
parents:
8392
diff
changeset
|
265 : MArray2<Complex> (Array2<Complex> (rv, 1, rv.length ())) |
2349 | 266 { |
267 } | |
268 | |
269 ComplexMatrix::ComplexMatrix (const ComplexColumnVector& cv) | |
8614
5114ea5a41b5
use shallow copying in Matrix/RowVector/ColumnVector conversions
Jaroslav Hajek <highegg@gmail.com>
parents:
8392
diff
changeset
|
270 : MArray2<Complex> (Array2<Complex> (cv, cv.length (), 1)) |
2349 | 271 { |
272 } | |
273 | |
458 | 274 ComplexMatrix::ComplexMatrix (const ComplexDiagMatrix& a) |
1214 | 275 : MArray2<Complex> (a.rows (), a.cols (), 0.0) |
458 | 276 { |
5275 | 277 for (octave_idx_type i = 0; i < a.length (); i++) |
458 | 278 elem (i, i) = a.elem (i, i); |
279 } | |
280 | |
5775 | 281 // FIXME -- could we use a templated mixed-type copy function |
1574 | 282 // here? |
283 | |
2828 | 284 ComplexMatrix::ComplexMatrix (const boolMatrix& a) |
3180 | 285 : MArray2<Complex> (a.rows (), a.cols (), 0.0) |
2828 | 286 { |
5275 | 287 for (octave_idx_type i = 0; i < a.rows (); i++) |
288 for (octave_idx_type j = 0; j < a.cols (); j++) | |
2828 | 289 elem (i, j) = a.elem (i, j); |
290 } | |
291 | |
1574 | 292 ComplexMatrix::ComplexMatrix (const charMatrix& a) |
3180 | 293 : MArray2<Complex> (a.rows (), a.cols (), 0.0) |
1574 | 294 { |
5275 | 295 for (octave_idx_type i = 0; i < a.rows (); i++) |
296 for (octave_idx_type j = 0; j < a.cols (); j++) | |
8956
d91fa4b20bbb
ensure nonnegative char -> real conversion
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
297 elem (i, j) = static_cast<unsigned char> (a.elem (i, j)); |
1574 | 298 } |
299 | |
9663
7e5b4de5fbfe
improve mixed real x complex ops
Jaroslav Hajek <highegg@gmail.com>
parents:
9661
diff
changeset
|
300 ComplexMatrix::ComplexMatrix (const Matrix& re, const Matrix& im) |
7e5b4de5fbfe
improve mixed real x complex ops
Jaroslav Hajek <highegg@gmail.com>
parents:
9661
diff
changeset
|
301 : MArray2<Complex> (re.rows (), re.cols ()) |
7e5b4de5fbfe
improve mixed real x complex ops
Jaroslav Hajek <highegg@gmail.com>
parents:
9661
diff
changeset
|
302 { |
7e5b4de5fbfe
improve mixed real x complex ops
Jaroslav Hajek <highegg@gmail.com>
parents:
9661
diff
changeset
|
303 if (im.rows () != rows () || im.cols () != cols ()) |
7e5b4de5fbfe
improve mixed real x complex ops
Jaroslav Hajek <highegg@gmail.com>
parents:
9661
diff
changeset
|
304 (*current_liboctave_error_handler) ("complex: internal error"); |
7e5b4de5fbfe
improve mixed real x complex ops
Jaroslav Hajek <highegg@gmail.com>
parents:
9661
diff
changeset
|
305 |
7e5b4de5fbfe
improve mixed real x complex ops
Jaroslav Hajek <highegg@gmail.com>
parents:
9661
diff
changeset
|
306 octave_idx_type nel = numel (); |
7e5b4de5fbfe
improve mixed real x complex ops
Jaroslav Hajek <highegg@gmail.com>
parents:
9661
diff
changeset
|
307 for (octave_idx_type i = 0; i < nel; i++) |
7e5b4de5fbfe
improve mixed real x complex ops
Jaroslav Hajek <highegg@gmail.com>
parents:
9661
diff
changeset
|
308 xelem (i) = Complex (re(i), im(i)); |
7e5b4de5fbfe
improve mixed real x complex ops
Jaroslav Hajek <highegg@gmail.com>
parents:
9661
diff
changeset
|
309 } |
7e5b4de5fbfe
improve mixed real x complex ops
Jaroslav Hajek <highegg@gmail.com>
parents:
9661
diff
changeset
|
310 |
2384 | 311 bool |
458 | 312 ComplexMatrix::operator == (const ComplexMatrix& a) const |
313 { | |
314 if (rows () != a.rows () || cols () != a.cols ()) | |
2384 | 315 return false; |
458 | 316 |
9550
3d6a9aea2aea
refactor binary & bool ops in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9528
diff
changeset
|
317 return mx_inline_equal (length (), data (), a.data ()); |
458 | 318 } |
319 | |
2384 | 320 bool |
458 | 321 ComplexMatrix::operator != (const ComplexMatrix& a) const |
322 { | |
323 return !(*this == a); | |
324 } | |
325 | |
2815 | 326 bool |
327 ComplexMatrix::is_hermitian (void) const | |
328 { | |
5275 | 329 octave_idx_type nr = rows (); |
330 octave_idx_type nc = cols (); | |
2815 | 331 |
332 if (is_square () && nr > 0) | |
333 { | |
5275 | 334 for (octave_idx_type i = 0; i < nr; i++) |
335 for (octave_idx_type j = i; j < nc; j++) | |
2815 | 336 if (elem (i, j) != conj (elem (j, i))) |
337 return false; | |
338 | |
339 return true; | |
340 } | |
341 | |
342 return false; | |
343 } | |
344 | |
458 | 345 // destructive insert/delete/reorder operations |
346 | |
347 ComplexMatrix& | |
5275 | 348 ComplexMatrix::insert (const Matrix& a, octave_idx_type r, octave_idx_type c) |
458 | 349 { |
5275 | 350 octave_idx_type a_nr = a.rows (); |
351 octave_idx_type a_nc = a.cols (); | |
1699 | 352 |
353 if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ()) | |
458 | 354 { |
355 (*current_liboctave_error_handler) ("range error for insert"); | |
356 return *this; | |
357 } | |
358 | |
4316 | 359 if (a_nr >0 && a_nc > 0) |
360 { | |
361 make_unique (); | |
362 | |
5275 | 363 for (octave_idx_type j = 0; j < a_nc; j++) |
364 for (octave_idx_type i = 0; i < a_nr; i++) | |
4316 | 365 xelem (r+i, c+j) = a.elem (i, j); |
366 } | |
458 | 367 |
368 return *this; | |
369 } | |
370 | |
371 ComplexMatrix& | |
5275 | 372 ComplexMatrix::insert (const RowVector& a, octave_idx_type r, octave_idx_type c) |
458 | 373 { |
5275 | 374 octave_idx_type a_len = a.length (); |
4316 | 375 |
1699 | 376 if (r < 0 || r >= rows () || c < 0 || c + a_len > cols ()) |
458 | 377 { |
378 (*current_liboctave_error_handler) ("range error for insert"); | |
379 return *this; | |
380 } | |
381 | |
4316 | 382 if (a_len > 0) |
383 { | |
384 make_unique (); | |
385 | |
5275 | 386 for (octave_idx_type i = 0; i < a_len; i++) |
4316 | 387 xelem (r, c+i) = a.elem (i); |
388 } | |
458 | 389 |
390 return *this; | |
391 } | |
392 | |
393 ComplexMatrix& | |
5275 | 394 ComplexMatrix::insert (const ColumnVector& a, octave_idx_type r, octave_idx_type c) |
458 | 395 { |
5275 | 396 octave_idx_type a_len = a.length (); |
4316 | 397 |
1699 | 398 if (r < 0 || r + a_len > rows () || c < 0 || c >= cols ()) |
458 | 399 { |
400 (*current_liboctave_error_handler) ("range error for insert"); | |
401 return *this; | |
402 } | |
403 | |
4316 | 404 if (a_len > 0) |
405 { | |
406 make_unique (); | |
407 | |
5275 | 408 for (octave_idx_type i = 0; i < a_len; i++) |
4316 | 409 xelem (r+i, c) = a.elem (i); |
410 } | |
458 | 411 |
412 return *this; | |
413 } | |
414 | |
415 ComplexMatrix& | |
5275 | 416 ComplexMatrix::insert (const DiagMatrix& a, octave_idx_type r, octave_idx_type c) |
458 | 417 { |
5275 | 418 octave_idx_type a_nr = a.rows (); |
419 octave_idx_type a_nc = a.cols (); | |
1699 | 420 |
421 if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ()) | |
458 | 422 { |
423 (*current_liboctave_error_handler) ("range error for insert"); | |
424 return *this; | |
425 } | |
426 | |
1699 | 427 fill (0.0, r, c, r + a_nr - 1, c + a_nc - 1); |
428 | |
5275 | 429 octave_idx_type a_len = a.length (); |
4316 | 430 |
431 if (a_len > 0) | |
432 { | |
433 make_unique (); | |
434 | |
5275 | 435 for (octave_idx_type i = 0; i < a_len; i++) |
4316 | 436 xelem (r+i, c+i) = a.elem (i, i); |
437 } | |
458 | 438 |
439 return *this; | |
440 } | |
441 | |
442 ComplexMatrix& | |
5275 | 443 ComplexMatrix::insert (const ComplexMatrix& a, octave_idx_type r, octave_idx_type c) |
458 | 444 { |
1561 | 445 Array2<Complex>::insert (a, r, c); |
458 | 446 return *this; |
447 } | |
448 | |
449 ComplexMatrix& | |
5275 | 450 ComplexMatrix::insert (const ComplexRowVector& a, octave_idx_type r, octave_idx_type c) |
458 | 451 { |
5275 | 452 octave_idx_type a_len = a.length (); |
1699 | 453 if (r < 0 || r >= rows () || c < 0 || c + a_len > cols ()) |
458 | 454 { |
455 (*current_liboctave_error_handler) ("range error for insert"); | |
456 return *this; | |
457 } | |
458 | |
5275 | 459 for (octave_idx_type i = 0; i < a_len; i++) |
458 | 460 elem (r, c+i) = a.elem (i); |
461 | |
462 return *this; | |
463 } | |
464 | |
465 ComplexMatrix& | |
5275 | 466 ComplexMatrix::insert (const ComplexColumnVector& a, octave_idx_type r, octave_idx_type c) |
458 | 467 { |
5275 | 468 octave_idx_type a_len = a.length (); |
4316 | 469 |
1699 | 470 if (r < 0 || r + a_len > rows () || c < 0 || c >= cols ()) |
458 | 471 { |
472 (*current_liboctave_error_handler) ("range error for insert"); | |
473 return *this; | |
474 } | |
475 | |
4316 | 476 if (a_len > 0) |
477 { | |
478 make_unique (); | |
479 | |
5275 | 480 for (octave_idx_type i = 0; i < a_len; i++) |
4316 | 481 xelem (r+i, c) = a.elem (i); |
482 } | |
458 | 483 |
484 return *this; | |
485 } | |
486 | |
487 ComplexMatrix& | |
5275 | 488 ComplexMatrix::insert (const ComplexDiagMatrix& a, octave_idx_type r, octave_idx_type c) |
458 | 489 { |
5275 | 490 octave_idx_type a_nr = a.rows (); |
491 octave_idx_type a_nc = a.cols (); | |
1699 | 492 |
493 if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ()) | |
458 | 494 { |
495 (*current_liboctave_error_handler) ("range error for insert"); | |
496 return *this; | |
497 } | |
498 | |
1699 | 499 fill (0.0, r, c, r + a_nr - 1, c + a_nc - 1); |
500 | |
5275 | 501 octave_idx_type a_len = a.length (); |
4316 | 502 |
503 if (a_len > 0) | |
504 { | |
505 make_unique (); | |
506 | |
5275 | 507 for (octave_idx_type i = 0; i < a_len; i++) |
4316 | 508 xelem (r+i, c+i) = a.elem (i, i); |
509 } | |
458 | 510 |
511 return *this; | |
512 } | |
513 | |
514 ComplexMatrix& | |
515 ComplexMatrix::fill (double val) | |
516 { | |
5275 | 517 octave_idx_type nr = rows (); |
518 octave_idx_type nc = cols (); | |
4316 | 519 |
458 | 520 if (nr > 0 && nc > 0) |
4316 | 521 { |
522 make_unique (); | |
523 | |
5275 | 524 for (octave_idx_type j = 0; j < nc; j++) |
525 for (octave_idx_type i = 0; i < nr; i++) | |
4316 | 526 xelem (i, j) = val; |
527 } | |
458 | 528 |
529 return *this; | |
530 } | |
531 | |
532 ComplexMatrix& | |
533 ComplexMatrix::fill (const Complex& val) | |
534 { | |
5275 | 535 octave_idx_type nr = rows (); |
536 octave_idx_type nc = cols (); | |
4316 | 537 |
458 | 538 if (nr > 0 && nc > 0) |
4316 | 539 { |
540 make_unique (); | |
541 | |
5275 | 542 for (octave_idx_type j = 0; j < nc; j++) |
543 for (octave_idx_type i = 0; i < nr; i++) | |
4316 | 544 xelem (i, j) = val; |
545 } | |
458 | 546 |
547 return *this; | |
548 } | |
549 | |
550 ComplexMatrix& | |
5275 | 551 ComplexMatrix::fill (double val, octave_idx_type r1, octave_idx_type c1, octave_idx_type r2, octave_idx_type c2) |
458 | 552 { |
5275 | 553 octave_idx_type nr = rows (); |
554 octave_idx_type nc = cols (); | |
4316 | 555 |
458 | 556 if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0 |
557 || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc) | |
558 { | |
559 (*current_liboctave_error_handler) ("range error for fill"); | |
560 return *this; | |
561 } | |
562 | |
5275 | 563 if (r1 > r2) { octave_idx_type tmp = r1; r1 = r2; r2 = tmp; } |
564 if (c1 > c2) { octave_idx_type tmp = c1; c1 = c2; c2 = tmp; } | |
458 | 565 |
4316 | 566 if (r2 >= r1 && c2 >= c1) |
567 { | |
568 make_unique (); | |
569 | |
5275 | 570 for (octave_idx_type j = c1; j <= c2; j++) |
571 for (octave_idx_type i = r1; i <= r2; i++) | |
4316 | 572 xelem (i, j) = val; |
573 } | |
458 | 574 |
575 return *this; | |
576 } | |
577 | |
578 ComplexMatrix& | |
5275 | 579 ComplexMatrix::fill (const Complex& val, octave_idx_type r1, octave_idx_type c1, octave_idx_type r2, octave_idx_type c2) |
458 | 580 { |
5275 | 581 octave_idx_type nr = rows (); |
582 octave_idx_type nc = cols (); | |
4316 | 583 |
458 | 584 if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0 |
585 || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc) | |
586 { | |
587 (*current_liboctave_error_handler) ("range error for fill"); | |
588 return *this; | |
589 } | |
590 | |
5275 | 591 if (r1 > r2) { octave_idx_type tmp = r1; r1 = r2; r2 = tmp; } |
592 if (c1 > c2) { octave_idx_type tmp = c1; c1 = c2; c2 = tmp; } | |
458 | 593 |
4316 | 594 if (r2 >= r1 && c2 >=c1) |
595 { | |
596 make_unique (); | |
597 | |
5275 | 598 for (octave_idx_type j = c1; j <= c2; j++) |
599 for (octave_idx_type i = r1; i <= r2; i++) | |
4316 | 600 xelem (i, j) = val; |
601 } | |
458 | 602 |
603 return *this; | |
604 } | |
605 | |
606 ComplexMatrix | |
607 ComplexMatrix::append (const Matrix& a) const | |
608 { | |
5275 | 609 octave_idx_type nr = rows (); |
610 octave_idx_type nc = cols (); | |
458 | 611 if (nr != a.rows ()) |
612 { | |
613 (*current_liboctave_error_handler) ("row dimension mismatch for append"); | |
614 return *this; | |
615 } | |
616 | |
5275 | 617 octave_idx_type nc_insert = nc; |
458 | 618 ComplexMatrix retval (nr, nc + a.cols ()); |
619 retval.insert (*this, 0, 0); | |
620 retval.insert (a, 0, nc_insert); | |
621 return retval; | |
622 } | |
623 | |
624 ComplexMatrix | |
625 ComplexMatrix::append (const RowVector& a) const | |
626 { | |
5275 | 627 octave_idx_type nr = rows (); |
628 octave_idx_type nc = cols (); | |
458 | 629 if (nr != 1) |
630 { | |
631 (*current_liboctave_error_handler) ("row dimension mismatch for append"); | |
632 return *this; | |
633 } | |
634 | |
5275 | 635 octave_idx_type nc_insert = nc; |
458 | 636 ComplexMatrix retval (nr, nc + a.length ()); |
637 retval.insert (*this, 0, 0); | |
638 retval.insert (a, 0, nc_insert); | |
639 return retval; | |
640 } | |
641 | |
642 ComplexMatrix | |
643 ComplexMatrix::append (const ColumnVector& a) const | |
644 { | |
5275 | 645 octave_idx_type nr = rows (); |
646 octave_idx_type nc = cols (); | |
458 | 647 if (nr != a.length ()) |
648 { | |
649 (*current_liboctave_error_handler) ("row dimension mismatch for append"); | |
650 return *this; | |
651 } | |
652 | |
5275 | 653 octave_idx_type nc_insert = nc; |
458 | 654 ComplexMatrix retval (nr, nc + 1); |
655 retval.insert (*this, 0, 0); | |
656 retval.insert (a, 0, nc_insert); | |
657 return retval; | |
658 } | |
659 | |
660 ComplexMatrix | |
661 ComplexMatrix::append (const DiagMatrix& a) const | |
662 { | |
5275 | 663 octave_idx_type nr = rows (); |
664 octave_idx_type nc = cols (); | |
458 | 665 if (nr != a.rows ()) |
666 { | |
667 (*current_liboctave_error_handler) ("row dimension mismatch for append"); | |
668 return *this; | |
669 } | |
670 | |
5275 | 671 octave_idx_type nc_insert = nc; |
458 | 672 ComplexMatrix retval (nr, nc + a.cols ()); |
673 retval.insert (*this, 0, 0); | |
674 retval.insert (a, 0, nc_insert); | |
675 return retval; | |
676 } | |
677 | |
678 ComplexMatrix | |
679 ComplexMatrix::append (const ComplexMatrix& a) const | |
680 { | |
5275 | 681 octave_idx_type nr = rows (); |
682 octave_idx_type nc = cols (); | |
458 | 683 if (nr != a.rows ()) |
684 { | |
685 (*current_liboctave_error_handler) ("row dimension mismatch for append"); | |
686 return *this; | |
687 } | |
688 | |
5275 | 689 octave_idx_type nc_insert = nc; |
458 | 690 ComplexMatrix retval (nr, nc + a.cols ()); |
691 retval.insert (*this, 0, 0); | |
692 retval.insert (a, 0, nc_insert); | |
693 return retval; | |
694 } | |
695 | |
696 ComplexMatrix | |
697 ComplexMatrix::append (const ComplexRowVector& a) const | |
698 { | |
5275 | 699 octave_idx_type nr = rows (); |
700 octave_idx_type nc = cols (); | |
458 | 701 if (nr != 1) |
702 { | |
703 (*current_liboctave_error_handler) ("row dimension mismatch for append"); | |
704 return *this; | |
705 } | |
706 | |
5275 | 707 octave_idx_type nc_insert = nc; |
458 | 708 ComplexMatrix retval (nr, nc + a.length ()); |
709 retval.insert (*this, 0, 0); | |
710 retval.insert (a, 0, nc_insert); | |
711 return retval; | |
712 } | |
713 | |
714 ComplexMatrix | |
715 ComplexMatrix::append (const ComplexColumnVector& a) const | |
716 { | |
5275 | 717 octave_idx_type nr = rows (); |
718 octave_idx_type nc = cols (); | |
458 | 719 if (nr != a.length ()) |
720 { | |
721 (*current_liboctave_error_handler) ("row dimension mismatch for append"); | |
722 return *this; | |
723 } | |
724 | |
5275 | 725 octave_idx_type nc_insert = nc; |
458 | 726 ComplexMatrix retval (nr, nc + 1); |
727 retval.insert (*this, 0, 0); | |
728 retval.insert (a, 0, nc_insert); | |
729 return retval; | |
730 } | |
731 | |
732 ComplexMatrix | |
733 ComplexMatrix::append (const ComplexDiagMatrix& a) const | |
734 { | |
5275 | 735 octave_idx_type nr = rows (); |
736 octave_idx_type nc = cols (); | |
458 | 737 if (nr != a.rows ()) |
738 { | |
739 (*current_liboctave_error_handler) ("row dimension mismatch for append"); | |
740 return *this; | |
741 } | |
742 | |
5275 | 743 octave_idx_type nc_insert = nc; |
458 | 744 ComplexMatrix retval (nr, nc + a.cols ()); |
745 retval.insert (*this, 0, 0); | |
746 retval.insert (a, 0, nc_insert); | |
747 return retval; | |
748 } | |
749 | |
750 ComplexMatrix | |
751 ComplexMatrix::stack (const Matrix& a) const | |
752 { | |
5275 | 753 octave_idx_type nr = rows (); |
754 octave_idx_type nc = cols (); | |
458 | 755 if (nc != a.cols ()) |
756 { | |
757 (*current_liboctave_error_handler) | |
758 ("column dimension mismatch for stack"); | |
759 return *this; | |
760 } | |
761 | |
5275 | 762 octave_idx_type nr_insert = nr; |
458 | 763 ComplexMatrix retval (nr + a.rows (), nc); |
764 retval.insert (*this, 0, 0); | |
765 retval.insert (a, nr_insert, 0); | |
766 return retval; | |
767 } | |
768 | |
769 ComplexMatrix | |
770 ComplexMatrix::stack (const RowVector& a) const | |
771 { | |
5275 | 772 octave_idx_type nr = rows (); |
773 octave_idx_type nc = cols (); | |
458 | 774 if (nc != a.length ()) |
775 { | |
776 (*current_liboctave_error_handler) | |
777 ("column dimension mismatch for stack"); | |
778 return *this; | |
779 } | |
780 | |
5275 | 781 octave_idx_type nr_insert = nr; |
458 | 782 ComplexMatrix retval (nr + 1, nc); |
783 retval.insert (*this, 0, 0); | |
784 retval.insert (a, nr_insert, 0); | |
785 return retval; | |
786 } | |
787 | |
788 ComplexMatrix | |
789 ComplexMatrix::stack (const ColumnVector& a) const | |
790 { | |
5275 | 791 octave_idx_type nr = rows (); |
792 octave_idx_type nc = cols (); | |
458 | 793 if (nc != 1) |
794 { | |
795 (*current_liboctave_error_handler) | |
796 ("column dimension mismatch for stack"); | |
797 return *this; | |
798 } | |
799 | |
5275 | 800 octave_idx_type nr_insert = nr; |
458 | 801 ComplexMatrix retval (nr + a.length (), nc); |
802 retval.insert (*this, 0, 0); | |
803 retval.insert (a, nr_insert, 0); | |
804 return retval; | |
805 } | |
806 | |
807 ComplexMatrix | |
808 ComplexMatrix::stack (const DiagMatrix& a) const | |
809 { | |
5275 | 810 octave_idx_type nr = rows (); |
811 octave_idx_type nc = cols (); | |
458 | 812 if (nc != a.cols ()) |
813 { | |
814 (*current_liboctave_error_handler) | |
815 ("column dimension mismatch for stack"); | |
816 return *this; | |
817 } | |
818 | |
5275 | 819 octave_idx_type nr_insert = nr; |
458 | 820 ComplexMatrix retval (nr + a.rows (), nc); |
821 retval.insert (*this, 0, 0); | |
822 retval.insert (a, nr_insert, 0); | |
823 return retval; | |
824 } | |
825 | |
826 ComplexMatrix | |
827 ComplexMatrix::stack (const ComplexMatrix& a) const | |
828 { | |
5275 | 829 octave_idx_type nr = rows (); |
830 octave_idx_type nc = cols (); | |
458 | 831 if (nc != a.cols ()) |
832 { | |
833 (*current_liboctave_error_handler) | |
834 ("column dimension mismatch for stack"); | |
835 return *this; | |
836 } | |
837 | |
5275 | 838 octave_idx_type nr_insert = nr; |
458 | 839 ComplexMatrix retval (nr + a.rows (), nc); |
840 retval.insert (*this, 0, 0); | |
841 retval.insert (a, nr_insert, 0); | |
842 return retval; | |
843 } | |
844 | |
845 ComplexMatrix | |
846 ComplexMatrix::stack (const ComplexRowVector& a) const | |
847 { | |
5275 | 848 octave_idx_type nr = rows (); |
849 octave_idx_type nc = cols (); | |
458 | 850 if (nc != a.length ()) |
851 { | |
852 (*current_liboctave_error_handler) | |
853 ("column dimension mismatch for stack"); | |
854 return *this; | |
855 } | |
856 | |
5275 | 857 octave_idx_type nr_insert = nr; |
458 | 858 ComplexMatrix retval (nr + 1, nc); |
859 retval.insert (*this, 0, 0); | |
860 retval.insert (a, nr_insert, 0); | |
861 return retval; | |
862 } | |
863 | |
864 ComplexMatrix | |
865 ComplexMatrix::stack (const ComplexColumnVector& a) const | |
866 { | |
5275 | 867 octave_idx_type nr = rows (); |
868 octave_idx_type nc = cols (); | |
458 | 869 if (nc != 1) |
870 { | |
871 (*current_liboctave_error_handler) | |
872 ("column dimension mismatch for stack"); | |
873 return *this; | |
874 } | |
875 | |
5275 | 876 octave_idx_type nr_insert = nr; |
458 | 877 ComplexMatrix retval (nr + a.length (), nc); |
878 retval.insert (*this, 0, 0); | |
879 retval.insert (a, nr_insert, 0); | |
880 return retval; | |
881 } | |
882 | |
883 ComplexMatrix | |
884 ComplexMatrix::stack (const ComplexDiagMatrix& a) const | |
885 { | |
5275 | 886 octave_idx_type nr = rows (); |
887 octave_idx_type nc = cols (); | |
458 | 888 if (nc != a.cols ()) |
889 { | |
890 (*current_liboctave_error_handler) | |
891 ("column dimension mismatch for stack"); | |
892 return *this; | |
893 } | |
894 | |
5275 | 895 octave_idx_type nr_insert = nr; |
458 | 896 ComplexMatrix retval (nr + a.rows (), nc); |
897 retval.insert (*this, 0, 0); | |
898 retval.insert (a, nr_insert, 0); | |
899 return retval; | |
900 } | |
901 | |
902 ComplexMatrix | |
903 conj (const ComplexMatrix& a) | |
904 { | |
8650
a1ae2aae903e
abs,real,imag,conj: use code from mx-inlines rather than the generic map
Jaroslav Hajek <highegg@gmail.com>
parents:
8614
diff
changeset
|
905 return ComplexMatrix (mx_inline_conj_dup (a.data (), a.length ()), |
a1ae2aae903e
abs,real,imag,conj: use code from mx-inlines rather than the generic map
Jaroslav Hajek <highegg@gmail.com>
parents:
8614
diff
changeset
|
906 a.rows (), a.cols ()); |
458 | 907 } |
908 | |
909 // resize is the destructive equivalent for this one | |
910 | |
911 ComplexMatrix | |
5275 | 912 ComplexMatrix::extract (octave_idx_type r1, octave_idx_type c1, octave_idx_type r2, octave_idx_type c2) const |
458 | 913 { |
5275 | 914 if (r1 > r2) { octave_idx_type tmp = r1; r1 = r2; r2 = tmp; } |
915 if (c1 > c2) { octave_idx_type tmp = c1; c1 = c2; c2 = tmp; } | |
916 | |
917 octave_idx_type new_r = r2 - r1 + 1; | |
918 octave_idx_type new_c = c2 - c1 + 1; | |
458 | 919 |
920 ComplexMatrix result (new_r, new_c); | |
921 | |
5275 | 922 for (octave_idx_type j = 0; j < new_c; j++) |
923 for (octave_idx_type i = 0; i < new_r; i++) | |
4316 | 924 result.xelem (i, j) = elem (r1+i, c1+j); |
925 | |
926 return result; | |
927 } | |
928 | |
929 ComplexMatrix | |
5275 | 930 ComplexMatrix::extract_n (octave_idx_type r1, octave_idx_type c1, octave_idx_type nr, octave_idx_type nc) const |
4316 | 931 { |
932 ComplexMatrix result (nr, nc); | |
933 | |
5275 | 934 for (octave_idx_type j = 0; j < nc; j++) |
935 for (octave_idx_type i = 0; i < nr; i++) | |
4316 | 936 result.xelem (i, j) = elem (r1+i, c1+j); |
458 | 937 |
938 return result; | |
939 } | |
940 | |
941 // extract row or column i. | |
942 | |
943 ComplexRowVector | |
5275 | 944 ComplexMatrix::row (octave_idx_type i) const |
458 | 945 { |
8614
5114ea5a41b5
use shallow copying in Matrix/RowVector/ColumnVector conversions
Jaroslav Hajek <highegg@gmail.com>
parents:
8392
diff
changeset
|
946 return MArray<Complex> (index (idx_vector (i), idx_vector::colon)); |
458 | 947 } |
948 | |
949 ComplexColumnVector | |
5275 | 950 ComplexMatrix::column (octave_idx_type i) const |
458 | 951 { |
8614
5114ea5a41b5
use shallow copying in Matrix/RowVector/ColumnVector conversions
Jaroslav Hajek <highegg@gmail.com>
parents:
8392
diff
changeset
|
952 return MArray<Complex> (index (idx_vector::colon, idx_vector (i))); |
458 | 953 } |
954 | |
955 ComplexMatrix | |
956 ComplexMatrix::inverse (void) const | |
957 { | |
5275 | 958 octave_idx_type info; |
7788 | 959 double rcon; |
6207 | 960 MatrixType mattype (*this); |
7788 | 961 return inverse (mattype, info, rcon, 0, 0); |
6207 | 962 } |
963 | |
964 ComplexMatrix | |
6479 | 965 ComplexMatrix::inverse (octave_idx_type& info) const |
966 { | |
7788 | 967 double rcon; |
6479 | 968 MatrixType mattype (*this); |
7788 | 969 return inverse (mattype, info, rcon, 0, 0); |
6479 | 970 } |
971 | |
972 ComplexMatrix | |
7788 | 973 ComplexMatrix::inverse (octave_idx_type& info, double& rcon, int force, |
6479 | 974 int calc_cond) const |
975 { | |
976 MatrixType mattype (*this); | |
7788 | 977 return inverse (mattype, info, rcon, force, calc_cond); |
6479 | 978 } |
979 | |
980 ComplexMatrix | |
6207 | 981 ComplexMatrix::inverse (MatrixType &mattype) const |
982 { | |
983 octave_idx_type info; | |
7788 | 984 double rcon; |
985 return inverse (mattype, info, rcon, 0, 0); | |
6207 | 986 } |
987 | |
988 ComplexMatrix | |
989 ComplexMatrix::inverse (MatrixType &mattype, octave_idx_type& info) const | |
990 { | |
7788 | 991 double rcon; |
992 return inverse (mattype, info, rcon, 0, 0); | |
458 | 993 } |
994 | |
995 ComplexMatrix | |
6207 | 996 ComplexMatrix::tinverse (MatrixType &mattype, octave_idx_type& info, |
7788 | 997 double& rcon, int force, int calc_cond) const |
458 | 998 { |
6207 | 999 ComplexMatrix retval; |
1000 | |
1001 octave_idx_type nr = rows (); | |
1002 octave_idx_type nc = cols (); | |
1003 | |
1004 if (nr != nc || nr == 0 || nc == 0) | |
1005 (*current_liboctave_error_handler) ("inverse requires square matrix"); | |
1006 else | |
1007 { | |
1008 int typ = mattype.type (); | |
1009 char uplo = (typ == MatrixType::Lower ? 'L' : 'U'); | |
1010 char udiag = 'N'; | |
1011 retval = *this; | |
1012 Complex *tmp_data = retval.fortran_vec (); | |
1013 | |
1014 F77_XFCN (ztrtri, ZTRTRI, (F77_CONST_CHAR_ARG2 (&uplo, 1), | |
1015 F77_CONST_CHAR_ARG2 (&udiag, 1), | |
1016 nr, tmp_data, nr, info | |
1017 F77_CHAR_ARG_LEN (1) | |
1018 F77_CHAR_ARG_LEN (1))); | |
1019 | |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1020 // Throw-away extra info LAPACK gives so as to not change output. |
7788 | 1021 rcon = 0.0; |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1022 if (info != 0) |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1023 info = -1; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1024 else if (calc_cond) |
6207 | 1025 { |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1026 octave_idx_type ztrcon_info = 0; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1027 char job = '1'; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1028 |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1029 OCTAVE_LOCAL_BUFFER (Complex, cwork, 2*nr); |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1030 OCTAVE_LOCAL_BUFFER (double, rwork, nr); |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1031 |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1032 F77_XFCN (ztrcon, ZTRCON, (F77_CONST_CHAR_ARG2 (&job, 1), |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1033 F77_CONST_CHAR_ARG2 (&uplo, 1), |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1034 F77_CONST_CHAR_ARG2 (&udiag, 1), |
7788 | 1035 nr, tmp_data, nr, rcon, |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1036 cwork, rwork, ztrcon_info |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1037 F77_CHAR_ARG_LEN (1) |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1038 F77_CHAR_ARG_LEN (1) |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1039 F77_CHAR_ARG_LEN (1))); |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1040 |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1041 if (ztrcon_info != 0) |
6207 | 1042 info = -1; |
1043 } | |
1044 | |
1045 if (info == -1 && ! force) | |
1046 retval = *this; // Restore matrix contents. | |
1047 } | |
1048 | |
1049 return retval; | |
458 | 1050 } |
1051 | |
1052 ComplexMatrix | |
6207 | 1053 ComplexMatrix::finverse (MatrixType &mattype, octave_idx_type& info, |
7788 | 1054 double& rcon, int force, int calc_cond) const |
458 | 1055 { |
1948 | 1056 ComplexMatrix retval; |
1057 | |
5275 | 1058 octave_idx_type nr = rows (); |
1059 octave_idx_type nc = cols (); | |
1948 | 1060 |
458 | 1061 if (nr != nc) |
1948 | 1062 (*current_liboctave_error_handler) ("inverse requires square matrix"); |
458 | 1063 else |
1064 { | |
5275 | 1065 Array<octave_idx_type> ipvt (nr); |
1066 octave_idx_type *pipvt = ipvt.fortran_vec (); | |
1948 | 1067 |
1068 retval = *this; | |
1069 Complex *tmp_data = retval.fortran_vec (); | |
1070 | |
4329 | 1071 Array<Complex> z(1); |
5275 | 1072 octave_idx_type lwork = -1; |
4330 | 1073 |
1074 // Query the optimum work array size. | |
4329 | 1075 |
1076 F77_XFCN (zgetri, ZGETRI, (nc, tmp_data, nr, pipvt, | |
1077 z.fortran_vec (), lwork, info)); | |
1078 | |
5315 | 1079 lwork = static_cast<octave_idx_type> (std::real(z(0))); |
4329 | 1080 lwork = (lwork < 2 *nc ? 2*nc : lwork); |
1081 z.resize (lwork); | |
1082 Complex *pz = z.fortran_vec (); | |
1083 | |
1084 info = 0; | |
1085 | |
4330 | 1086 // Calculate the norm of the matrix, for later use. |
4329 | 1087 double anorm; |
1088 if (calc_cond) | |
5275 | 1089 anorm = retval.abs().sum().row(static_cast<octave_idx_type>(0)).max(); |
4329 | 1090 |
1091 F77_XFCN (zgetrf, ZGETRF, (nc, nc, tmp_data, nr, pipvt, info)); | |
1948 | 1092 |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1093 // Throw-away extra info LAPACK gives so as to not change output. |
7788 | 1094 rcon = 0.0; |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1095 if (info != 0) |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1096 info = -1; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1097 else if (calc_cond) |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1098 { |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1099 // Now calculate the condition number for non-singular matrix. |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1100 octave_idx_type zgecon_info = 0; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1101 char job = '1'; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1102 Array<double> rz (2 * nc); |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1103 double *prz = rz.fortran_vec (); |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1104 F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1), |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1105 nc, tmp_data, nr, anorm, |
7788 | 1106 rcon, pz, prz, zgecon_info |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1107 F77_CHAR_ARG_LEN (1))); |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1108 |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1109 if (zgecon_info != 0) |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1110 info = -1; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1111 } |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1112 |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1113 if (info == -1 && ! force) |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1114 retval = *this; // Restore contents. |
1948 | 1115 else |
1116 { | |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1117 octave_idx_type zgetri_info = 0; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1118 |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1119 F77_XFCN (zgetri, ZGETRI, (nc, tmp_data, nr, pipvt, |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1120 pz, lwork, zgetri_info)); |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1121 |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
1122 if (zgetri_info != 0) |
1948 | 1123 info = -1; |
1124 } | |
6207 | 1125 |
1126 if (info != 0) | |
1127 mattype.mark_as_rectangular(); | |
458 | 1128 } |
4329 | 1129 |
1948 | 1130 return retval; |
458 | 1131 } |
1132 | |
1133 ComplexMatrix | |
6207 | 1134 ComplexMatrix::inverse (MatrixType &mattype, octave_idx_type& info, |
7788 | 1135 double& rcon, int force, int calc_cond) const |
6207 | 1136 { |
1137 int typ = mattype.type (false); | |
1138 ComplexMatrix ret; | |
1139 | |
1140 if (typ == MatrixType::Unknown) | |
1141 typ = mattype.type (*this); | |
1142 | |
1143 if (typ == MatrixType::Upper || typ == MatrixType::Lower) | |
7788 | 1144 ret = tinverse (mattype, info, rcon, force, calc_cond); |
6840 | 1145 else |
6207 | 1146 { |
1147 if (mattype.is_hermitian ()) | |
1148 { | |
6486 | 1149 ComplexCHOL chol (*this, info, calc_cond); |
6207 | 1150 if (info == 0) |
6486 | 1151 { |
1152 if (calc_cond) | |
7788 | 1153 rcon = chol.rcond(); |
6486 | 1154 else |
7788 | 1155 rcon = 1.0; |
6486 | 1156 ret = chol.inverse (); |
1157 } | |
6207 | 1158 else |
1159 mattype.mark_as_unsymmetric (); | |
1160 } | |
1161 | |
1162 if (!mattype.is_hermitian ()) | |
7788 | 1163 ret = finverse(mattype, info, rcon, force, calc_cond); |
1164 | |
1165 if ((mattype.is_hermitian () || calc_cond) && rcon == 0.) | |
6840 | 1166 ret = ComplexMatrix (rows (), columns (), Complex (octave_Inf, 0.)); |
6207 | 1167 } |
1168 | |
1169 return ret; | |
1170 } | |
1171 | |
1172 ComplexMatrix | |
4384 | 1173 ComplexMatrix::pseudo_inverse (double tol) const |
740 | 1174 { |
1549 | 1175 ComplexMatrix retval; |
1176 | |
3480 | 1177 ComplexSVD result (*this, SVD::economy); |
740 | 1178 |
1179 DiagMatrix S = result.singular_values (); | |
1180 ComplexMatrix U = result.left_singular_matrix (); | |
1181 ComplexMatrix V = result.right_singular_matrix (); | |
1182 | |
1183 ColumnVector sigma = S.diag (); | |
1184 | |
5275 | 1185 octave_idx_type r = sigma.length () - 1; |
1186 octave_idx_type nr = rows (); | |
1187 octave_idx_type nc = cols (); | |
740 | 1188 |
1189 if (tol <= 0.0) | |
1190 { | |
1191 if (nr > nc) | |
1192 tol = nr * sigma.elem (0) * DBL_EPSILON; | |
1193 else | |
1194 tol = nc * sigma.elem (0) * DBL_EPSILON; | |
1195 } | |
1196 | |
1197 while (r >= 0 && sigma.elem (r) < tol) | |
1198 r--; | |
1199 | |
1200 if (r < 0) | |
1549 | 1201 retval = ComplexMatrix (nc, nr, 0.0); |
740 | 1202 else |
1203 { | |
1204 ComplexMatrix Ur = U.extract (0, 0, nr-1, r); | |
1205 DiagMatrix D = DiagMatrix (sigma.extract (0, r)) . inverse (); | |
1206 ComplexMatrix Vr = V.extract (0, 0, nc-1, r); | |
1549 | 1207 retval = Vr * D * Ur.hermitian (); |
740 | 1208 } |
1549 | 1209 |
1210 return retval; | |
740 | 1211 } |
1212 | |
9523
0ce82753dd72
more configure changes for libraries
John W. Eaton <jwe@octave.org>
parents:
9469
diff
changeset
|
1213 #if defined (HAVE_FFTW) |
3827 | 1214 |
1215 ComplexMatrix | |
1216 ComplexMatrix::fourier (void) const | |
1217 { | |
1218 size_t nr = rows (); | |
1219 size_t nc = cols (); | |
1220 | |
1221 ComplexMatrix retval (nr, nc); | |
1222 | |
1223 size_t npts, nsamples; | |
1224 | |
1225 if (nr == 1 || nc == 1) | |
1226 { | |
1227 npts = nr > nc ? nr : nc; | |
1228 nsamples = 1; | |
1229 } | |
1230 else | |
1231 { | |
1232 npts = nr; | |
1233 nsamples = nc; | |
1234 } | |
1235 | |
1236 const Complex *in (data ()); | |
1237 Complex *out (retval.fortran_vec ()); | |
1238 | |
4773 | 1239 octave_fftw::fft (in, out, npts, nsamples); |
3827 | 1240 |
1241 return retval; | |
1242 } | |
1243 | |
1244 ComplexMatrix | |
1245 ComplexMatrix::ifourier (void) const | |
1246 { | |
1247 size_t nr = rows (); | |
1248 size_t nc = cols (); | |
1249 | |
1250 ComplexMatrix retval (nr, nc); | |
1251 | |
1252 size_t npts, nsamples; | |
1253 | |
1254 if (nr == 1 || nc == 1) | |
1255 { | |
1256 npts = nr > nc ? nr : nc; | |
1257 nsamples = 1; | |
1258 } | |
1259 else | |
1260 { | |
1261 npts = nr; | |
1262 nsamples = nc; | |
1263 } | |
1264 | |
1265 const Complex *in (data ()); | |
1266 Complex *out (retval.fortran_vec ()); | |
1267 | |
4773 | 1268 octave_fftw::ifft (in, out, npts, nsamples); |
3827 | 1269 |
1270 return retval; | |
1271 } | |
1272 | |
1273 ComplexMatrix | |
1274 ComplexMatrix::fourier2d (void) const | |
1275 { | |
4773 | 1276 dim_vector dv(rows (), cols ()); |
1277 | |
1278 ComplexMatrix retval (rows (), cols ()); | |
1279 const Complex *in (data ()); | |
1280 Complex *out (retval.fortran_vec ()); | |
1281 | |
1282 octave_fftw::fftNd (in, out, 2, dv); | |
3827 | 1283 |
1284 return retval; | |
1285 } | |
1286 | |
1287 ComplexMatrix | |
1288 ComplexMatrix::ifourier2d (void) const | |
1289 { | |
4773 | 1290 dim_vector dv(rows (), cols ()); |
1291 | |
1292 ComplexMatrix retval (rows (), cols ()); | |
1293 const Complex *in (data ()); | |
1294 Complex *out (retval.fortran_vec ()); | |
1295 | |
1296 octave_fftw::ifftNd (in, out, 2, dv); | |
3827 | 1297 |
1298 return retval; | |
1299 } | |
1300 | |
1301 #else | |
1302 | |
9523
0ce82753dd72
more configure changes for libraries
John W. Eaton <jwe@octave.org>
parents:
9469
diff
changeset
|
1303 extern "C" |
0ce82753dd72
more configure changes for libraries
John W. Eaton <jwe@octave.org>
parents:
9469
diff
changeset
|
1304 { |
0ce82753dd72
more configure changes for libraries
John W. Eaton <jwe@octave.org>
parents:
9469
diff
changeset
|
1305 // Note that the original complex fft routines were not written for |
0ce82753dd72
more configure changes for libraries
John W. Eaton <jwe@octave.org>
parents:
9469
diff
changeset
|
1306 // double complex arguments. They have been modified by adding an |
0ce82753dd72
more configure changes for libraries
John W. Eaton <jwe@octave.org>
parents:
9469
diff
changeset
|
1307 // implicit double precision (a-h,o-z) statement at the beginning of |
0ce82753dd72
more configure changes for libraries
John W. Eaton <jwe@octave.org>
parents:
9469
diff
changeset
|
1308 // each subroutine. |
0ce82753dd72
more configure changes for libraries
John W. Eaton <jwe@octave.org>
parents:
9469
diff
changeset
|
1309 |
0ce82753dd72
more configure changes for libraries
John W. Eaton <jwe@octave.org>
parents:
9469
diff
changeset
|
1310 F77_RET_T |
0ce82753dd72
more configure changes for libraries
John W. Eaton <jwe@octave.org>
parents:
9469
diff
changeset
|
1311 F77_FUNC (zffti, ZFFTI) (const octave_idx_type&, Complex*); |
0ce82753dd72
more configure changes for libraries
John W. Eaton <jwe@octave.org>
parents:
9469
diff
changeset
|
1312 |
0ce82753dd72
more configure changes for libraries
John W. Eaton <jwe@octave.org>
parents:
9469
diff
changeset
|
1313 F77_RET_T |
0ce82753dd72
more configure changes for libraries
John W. Eaton <jwe@octave.org>
parents:
9469
diff
changeset
|
1314 F77_FUNC (zfftf, ZFFTF) (const octave_idx_type&, Complex*, Complex*); |
0ce82753dd72
more configure changes for libraries
John W. Eaton <jwe@octave.org>
parents:
9469
diff
changeset
|
1315 |
0ce82753dd72
more configure changes for libraries
John W. Eaton <jwe@octave.org>
parents:
9469
diff
changeset
|
1316 F77_RET_T |
0ce82753dd72
more configure changes for libraries
John W. Eaton <jwe@octave.org>
parents:
9469
diff
changeset
|
1317 F77_FUNC (zfftb, ZFFTB) (const octave_idx_type&, Complex*, Complex*); |
0ce82753dd72
more configure changes for libraries
John W. Eaton <jwe@octave.org>
parents:
9469
diff
changeset
|
1318 } |
0ce82753dd72
more configure changes for libraries
John W. Eaton <jwe@octave.org>
parents:
9469
diff
changeset
|
1319 |
740 | 1320 ComplexMatrix |
458 | 1321 ComplexMatrix::fourier (void) const |
1322 { | |
1948 | 1323 ComplexMatrix retval; |
1324 | |
5275 | 1325 octave_idx_type nr = rows (); |
1326 octave_idx_type nc = cols (); | |
1327 | |
1328 octave_idx_type npts, nsamples; | |
1948 | 1329 |
458 | 1330 if (nr == 1 || nc == 1) |
1331 { | |
1332 npts = nr > nc ? nr : nc; | |
1333 nsamples = 1; | |
1334 } | |
1335 else | |
1336 { | |
1337 npts = nr; | |
1338 nsamples = nc; | |
1339 } | |
1340 | |
5275 | 1341 octave_idx_type nn = 4*npts+15; |
1948 | 1342 |
1343 Array<Complex> wsave (nn); | |
1344 Complex *pwsave = wsave.fortran_vec (); | |
1345 | |
1346 retval = *this; | |
1347 Complex *tmp_data = retval.fortran_vec (); | |
1348 | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7788
diff
changeset
|
1349 F77_FUNC (zffti, ZFFTI) (npts, pwsave); |
458 | 1350 |
5275 | 1351 for (octave_idx_type j = 0; j < nsamples; j++) |
4153 | 1352 { |
10142
829e69ec3110
make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents:
10124
diff
changeset
|
1353 octave_quit (); |
4153 | 1354 |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7788
diff
changeset
|
1355 F77_FUNC (zfftf, ZFFTF) (npts, &tmp_data[npts*j], pwsave); |
4153 | 1356 } |
1948 | 1357 |
1358 return retval; | |
458 | 1359 } |
1360 | |
1361 ComplexMatrix | |
1362 ComplexMatrix::ifourier (void) const | |
1363 { | |
1948 | 1364 ComplexMatrix retval; |
1365 | |
5275 | 1366 octave_idx_type nr = rows (); |
1367 octave_idx_type nc = cols (); | |
1368 | |
1369 octave_idx_type npts, nsamples; | |
1948 | 1370 |
458 | 1371 if (nr == 1 || nc == 1) |
1372 { | |
1373 npts = nr > nc ? nr : nc; | |
1374 nsamples = 1; | |
1375 } | |
1376 else | |
1377 { | |
1378 npts = nr; | |
1379 nsamples = nc; | |
1380 } | |
1381 | |
5275 | 1382 octave_idx_type nn = 4*npts+15; |
1948 | 1383 |
1384 Array<Complex> wsave (nn); | |
1385 Complex *pwsave = wsave.fortran_vec (); | |
1386 | |
1387 retval = *this; | |
1388 Complex *tmp_data = retval.fortran_vec (); | |
1389 | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7788
diff
changeset
|
1390 F77_FUNC (zffti, ZFFTI) (npts, pwsave); |
458 | 1391 |
5275 | 1392 for (octave_idx_type j = 0; j < nsamples; j++) |
4153 | 1393 { |
10142
829e69ec3110
make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents:
10124
diff
changeset
|
1394 octave_quit (); |
4153 | 1395 |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7788
diff
changeset
|
1396 F77_FUNC (zfftb, ZFFTB) (npts, &tmp_data[npts*j], pwsave); |
4153 | 1397 } |
458 | 1398 |
5275 | 1399 for (octave_idx_type j = 0; j < npts*nsamples; j++) |
3572 | 1400 tmp_data[j] = tmp_data[j] / static_cast<double> (npts); |
458 | 1401 |
1948 | 1402 return retval; |
458 | 1403 } |
1404 | |
677 | 1405 ComplexMatrix |
1406 ComplexMatrix::fourier2d (void) const | |
1407 { | |
1948 | 1408 ComplexMatrix retval; |
1409 | |
5275 | 1410 octave_idx_type nr = rows (); |
1411 octave_idx_type nc = cols (); | |
1412 | |
1413 octave_idx_type npts, nsamples; | |
1948 | 1414 |
677 | 1415 if (nr == 1 || nc == 1) |
1416 { | |
1417 npts = nr > nc ? nr : nc; | |
1418 nsamples = 1; | |
1419 } | |
1420 else | |
1421 { | |
1422 npts = nr; | |
1423 nsamples = nc; | |
1424 } | |
1425 | |
5275 | 1426 octave_idx_type nn = 4*npts+15; |
1948 | 1427 |
1428 Array<Complex> wsave (nn); | |
1429 Complex *pwsave = wsave.fortran_vec (); | |
1430 | |
1431 retval = *this; | |
1432 Complex *tmp_data = retval.fortran_vec (); | |
1433 | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7788
diff
changeset
|
1434 F77_FUNC (zffti, ZFFTI) (npts, pwsave); |
677 | 1435 |
5275 | 1436 for (octave_idx_type j = 0; j < nsamples; j++) |
4153 | 1437 { |
10142
829e69ec3110
make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents:
10124
diff
changeset
|
1438 octave_quit (); |
4153 | 1439 |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7788
diff
changeset
|
1440 F77_FUNC (zfftf, ZFFTF) (npts, &tmp_data[npts*j], pwsave); |
4153 | 1441 } |
677 | 1442 |
1443 npts = nc; | |
1444 nsamples = nr; | |
1445 nn = 4*npts+15; | |
1948 | 1446 |
1447 wsave.resize (nn); | |
1448 pwsave = wsave.fortran_vec (); | |
1449 | |
4773 | 1450 Array<Complex> tmp (npts); |
1451 Complex *prow = tmp.fortran_vec (); | |
1948 | 1452 |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7788
diff
changeset
|
1453 F77_FUNC (zffti, ZFFTI) (npts, pwsave); |
677 | 1454 |
5275 | 1455 for (octave_idx_type j = 0; j < nsamples; j++) |
677 | 1456 { |
10142
829e69ec3110
make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents:
10124
diff
changeset
|
1457 octave_quit (); |
4153 | 1458 |
5275 | 1459 for (octave_idx_type i = 0; i < npts; i++) |
1948 | 1460 prow[i] = tmp_data[i*nr + j]; |
1461 | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7788
diff
changeset
|
1462 F77_FUNC (zfftf, ZFFTF) (npts, prow, pwsave); |
677 | 1463 |
5275 | 1464 for (octave_idx_type i = 0; i < npts; i++) |
1948 | 1465 tmp_data[i*nr + j] = prow[i]; |
677 | 1466 } |
1467 | |
1948 | 1468 return retval; |
677 | 1469 } |
1470 | |
1471 ComplexMatrix | |
1472 ComplexMatrix::ifourier2d (void) const | |
1473 { | |
1948 | 1474 ComplexMatrix retval; |
1475 | |
5275 | 1476 octave_idx_type nr = rows (); |
1477 octave_idx_type nc = cols (); | |
1478 | |
1479 octave_idx_type npts, nsamples; | |
1948 | 1480 |
677 | 1481 if (nr == 1 || nc == 1) |
1482 { | |
1483 npts = nr > nc ? nr : nc; | |
1484 nsamples = 1; | |
1485 } | |
1486 else | |
1487 { | |
1488 npts = nr; | |
1489 nsamples = nc; | |
1490 } | |
1491 | |
5275 | 1492 octave_idx_type nn = 4*npts+15; |
1948 | 1493 |
1494 Array<Complex> wsave (nn); | |
1495 Complex *pwsave = wsave.fortran_vec (); | |
1496 | |
1497 retval = *this; | |
1498 Complex *tmp_data = retval.fortran_vec (); | |
1499 | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7788
diff
changeset
|
1500 F77_FUNC (zffti, ZFFTI) (npts, pwsave); |
677 | 1501 |
5275 | 1502 for (octave_idx_type j = 0; j < nsamples; j++) |
4153 | 1503 { |
10142
829e69ec3110
make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents:
10124
diff
changeset
|
1504 octave_quit (); |
4153 | 1505 |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7788
diff
changeset
|
1506 F77_FUNC (zfftb, ZFFTB) (npts, &tmp_data[npts*j], pwsave); |
4153 | 1507 } |
677 | 1508 |
5275 | 1509 for (octave_idx_type j = 0; j < npts*nsamples; j++) |
3572 | 1510 tmp_data[j] = tmp_data[j] / static_cast<double> (npts); |
677 | 1511 |
1512 npts = nc; | |
1513 nsamples = nr; | |
1514 nn = 4*npts+15; | |
1948 | 1515 |
1516 wsave.resize (nn); | |
1517 pwsave = wsave.fortran_vec (); | |
1518 | |
4773 | 1519 Array<Complex> tmp (npts); |
1520 Complex *prow = tmp.fortran_vec (); | |
1948 | 1521 |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7788
diff
changeset
|
1522 F77_FUNC (zffti, ZFFTI) (npts, pwsave); |
677 | 1523 |
5275 | 1524 for (octave_idx_type j = 0; j < nsamples; j++) |
677 | 1525 { |
10142
829e69ec3110
make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents:
10124
diff
changeset
|
1526 octave_quit (); |
4153 | 1527 |
5275 | 1528 for (octave_idx_type i = 0; i < npts; i++) |
1948 | 1529 prow[i] = tmp_data[i*nr + j]; |
1530 | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7788
diff
changeset
|
1531 F77_FUNC (zfftb, ZFFTB) (npts, prow, pwsave); |
677 | 1532 |
5275 | 1533 for (octave_idx_type i = 0; i < npts; i++) |
3572 | 1534 tmp_data[i*nr + j] = prow[i] / static_cast<double> (npts); |
677 | 1535 } |
1536 | |
1948 | 1537 return retval; |
677 | 1538 } |
1539 | |
3827 | 1540 #endif |
1541 | |
458 | 1542 ComplexDET |
1543 ComplexMatrix::determinant (void) const | |
1544 { | |
5275 | 1545 octave_idx_type info; |
7788 | 1546 double rcon; |
1547 return determinant (info, rcon, 0); | |
458 | 1548 } |
1549 | |
1550 ComplexDET | |
5275 | 1551 ComplexMatrix::determinant (octave_idx_type& info) const |
458 | 1552 { |
7788 | 1553 double rcon; |
1554 return determinant (info, rcon, 0); | |
458 | 1555 } |
1556 | |
1557 ComplexDET | |
7788 | 1558 ComplexMatrix::determinant (octave_idx_type& info, double& rcon, int calc_cond) const |
458 | 1559 { |
8336
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1560 MatrixType mattype (*this); |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1561 return determinant (mattype, info, rcon, calc_cond); |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1562 } |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1563 |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1564 ComplexDET |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1565 ComplexMatrix::determinant (MatrixType& mattype, |
8806 | 1566 octave_idx_type& info, double& rcon, |
1567 int calc_cond) const | |
8336
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1568 { |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1569 ComplexDET retval (1.0); |
458 | 1570 |
5275 | 1571 octave_idx_type nr = rows (); |
1572 octave_idx_type nc = cols (); | |
458 | 1573 |
8336
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1574 if (nr != nc) |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1575 (*current_liboctave_error_handler) ("matrix must be square"); |
458 | 1576 else |
1577 { | |
8806 | 1578 volatile int typ = mattype.type (); |
8336
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1579 |
8337
e02242c54c49
reuse matrix type detected in det
Jaroslav Hajek <highegg@gmail.com>
parents:
8336
diff
changeset
|
1580 if (typ == MatrixType::Unknown) |
e02242c54c49
reuse matrix type detected in det
Jaroslav Hajek <highegg@gmail.com>
parents:
8336
diff
changeset
|
1581 typ = mattype.type (*this); |
e02242c54c49
reuse matrix type detected in det
Jaroslav Hajek <highegg@gmail.com>
parents:
8336
diff
changeset
|
1582 |
8336
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1583 if (typ == MatrixType::Lower || typ == MatrixType::Upper) |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1584 { |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1585 for (octave_idx_type i = 0; i < nc; i++) |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1586 retval *= elem (i,i); |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1587 } |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1588 else if (typ == MatrixType::Hermitian) |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1589 { |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1590 ComplexMatrix atmp = *this; |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1591 Complex *tmp_data = atmp.fortran_vec (); |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1592 |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1593 info = 0; |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1594 double anorm = 0; |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1595 if (calc_cond) anorm = xnorm (*this, 1); |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1596 |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1597 |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1598 char job = 'L'; |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1599 F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr, |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1600 tmp_data, nr, info |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1601 F77_CHAR_ARG_LEN (1))); |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1602 |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1603 if (info != 0) |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1604 { |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1605 rcon = 0.0; |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1606 mattype.mark_as_unsymmetric (); |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1607 typ = MatrixType::Full; |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1608 } |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1609 else |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1610 { |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1611 Array<Complex> z (2 * nc); |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1612 Complex *pz = z.fortran_vec (); |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1613 Array<double> rz (nc); |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1614 double *prz = rz.fortran_vec (); |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1615 |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1616 F77_XFCN (zpocon, ZPOCON, (F77_CONST_CHAR_ARG2 (&job, 1), |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1617 nr, tmp_data, nr, anorm, |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1618 rcon, pz, prz, info |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1619 F77_CHAR_ARG_LEN (1))); |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1620 |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1621 if (info != 0) |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1622 rcon = 0.0; |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1623 |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1624 for (octave_idx_type i = 0; i < nc; i++) |
8337
e02242c54c49
reuse matrix type detected in det
Jaroslav Hajek <highegg@gmail.com>
parents:
8336
diff
changeset
|
1625 retval *= atmp (i,i); |
8336
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1626 |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1627 retval = retval.square (); |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1628 } |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1629 } |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1630 else if (typ != MatrixType::Full) |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1631 (*current_liboctave_error_handler) ("det: invalid dense matrix type"); |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1632 |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1633 if (typ == MatrixType::Full) |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1634 { |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1635 Array<octave_idx_type> ipvt (nr); |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1636 octave_idx_type *pipvt = ipvt.fortran_vec (); |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1637 |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1638 ComplexMatrix atmp = *this; |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1639 Complex *tmp_data = atmp.fortran_vec (); |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1640 |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1641 info = 0; |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1642 |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1643 // Calculate the norm of the matrix, for later use. |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1644 double anorm = 0; |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1645 if (calc_cond) anorm = xnorm (*this, 1); |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1646 |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1647 F77_XFCN (zgetrf, ZGETRF, (nr, nr, tmp_data, nr, pipvt, info)); |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1648 |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1649 // Throw-away extra info LAPACK gives so as to not change output. |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1650 rcon = 0.0; |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1651 if (info != 0) |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1652 { |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1653 info = -1; |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1654 retval = ComplexDET (); |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1655 } |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1656 else |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1657 { |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1658 if (calc_cond) |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1659 { |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1660 // Now calc the condition number for non-singular matrix. |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1661 char job = '1'; |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1662 Array<Complex> z (2 * nc); |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1663 Complex *pz = z.fortran_vec (); |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1664 Array<double> rz (2 * nc); |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1665 double *prz = rz.fortran_vec (); |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1666 |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1667 F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1), |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1668 nc, tmp_data, nr, anorm, |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1669 rcon, pz, prz, info |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1670 F77_CHAR_ARG_LEN (1))); |
8335 | 1671 } |
8336
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1672 |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1673 if (info != 0) |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1674 { |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1675 info = -1; |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1676 retval = ComplexDET (); |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1677 } |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1678 else |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1679 { |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1680 for (octave_idx_type i = 0; i < nc; i++) |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1681 { |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1682 Complex c = atmp(i,i); |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1683 retval *= (ipvt(i) != (i+1)) ? -c : c; |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1684 } |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1685 } |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1686 } |
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1687 } |
458 | 1688 } |
8336
9813c07ca946
make det take advantage of matrix type
Jaroslav Hajek <highegg@gmail.com>
parents:
8335
diff
changeset
|
1689 |
458 | 1690 return retval; |
1691 } | |
1692 | |
7788 | 1693 double |
1694 ComplexMatrix::rcond (void) const | |
1695 { | |
1696 MatrixType mattype (*this); | |
1697 return rcond (mattype); | |
1698 } | |
1699 | |
1700 double | |
1701 ComplexMatrix::rcond (MatrixType &mattype) const | |
1702 { | |
1703 double rcon; | |
1704 octave_idx_type nr = rows (); | |
1705 octave_idx_type nc = cols (); | |
1706 | |
1707 if (nr != nc) | |
1708 (*current_liboctave_error_handler) ("matrix must be square"); | |
1709 else if (nr == 0 || nc == 0) | |
1710 rcon = octave_Inf; | |
1711 else | |
1712 { | |
1713 int typ = mattype.type (); | |
1714 | |
1715 if (typ == MatrixType::Unknown) | |
1716 typ = mattype.type (*this); | |
1717 | |
1718 // Only calculate the condition number for LU/Cholesky | |
1719 if (typ == MatrixType::Upper) | |
1720 { | |
1721 const Complex *tmp_data = fortran_vec (); | |
1722 octave_idx_type info = 0; | |
1723 char norm = '1'; | |
1724 char uplo = 'U'; | |
1725 char dia = 'N'; | |
1726 | |
1727 Array<Complex> z (2 * nc); | |
1728 Complex *pz = z.fortran_vec (); | |
1729 Array<double> rz (nc); | |
1730 double *prz = rz.fortran_vec (); | |
1731 | |
1732 F77_XFCN (ztrcon, ZTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), | |
1733 F77_CONST_CHAR_ARG2 (&uplo, 1), | |
1734 F77_CONST_CHAR_ARG2 (&dia, 1), | |
1735 nr, tmp_data, nr, rcon, | |
1736 pz, prz, info | |
1737 F77_CHAR_ARG_LEN (1) | |
1738 F77_CHAR_ARG_LEN (1) | |
1739 F77_CHAR_ARG_LEN (1))); | |
1740 | |
1741 if (info != 0) | |
1742 rcon = 0; | |
1743 } | |
1744 else if (typ == MatrixType::Permuted_Upper) | |
1745 (*current_liboctave_error_handler) | |
1746 ("permuted triangular matrix not implemented"); | |
1747 else if (typ == MatrixType::Lower) | |
1748 { | |
1749 const Complex *tmp_data = fortran_vec (); | |
1750 octave_idx_type info = 0; | |
1751 char norm = '1'; | |
1752 char uplo = 'L'; | |
1753 char dia = 'N'; | |
1754 | |
1755 Array<Complex> z (2 * nc); | |
1756 Complex *pz = z.fortran_vec (); | |
1757 Array<double> rz (nc); | |
1758 double *prz = rz.fortran_vec (); | |
1759 | |
1760 F77_XFCN (ztrcon, ZTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), | |
1761 F77_CONST_CHAR_ARG2 (&uplo, 1), | |
1762 F77_CONST_CHAR_ARG2 (&dia, 1), | |
1763 nr, tmp_data, nr, rcon, | |
1764 pz, prz, info | |
1765 F77_CHAR_ARG_LEN (1) | |
1766 F77_CHAR_ARG_LEN (1) | |
1767 F77_CHAR_ARG_LEN (1))); | |
1768 | |
1769 if (info != 0) | |
1770 rcon = 0.0; | |
1771 } | |
1772 else if (typ == MatrixType::Permuted_Lower) | |
1773 (*current_liboctave_error_handler) | |
1774 ("permuted triangular matrix not implemented"); | |
1775 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian) | |
1776 { | |
1777 double anorm = -1.0; | |
1778 ComplexMatrix atmp = *this; | |
1779 Complex *tmp_data = atmp.fortran_vec (); | |
1780 | |
1781 if (typ == MatrixType::Hermitian) | |
1782 { | |
1783 octave_idx_type info = 0; | |
1784 char job = 'L'; | |
1785 anorm = atmp.abs().sum(). | |
1786 row(static_cast<octave_idx_type>(0)).max(); | |
1787 | |
1788 F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr, | |
1789 tmp_data, nr, info | |
1790 F77_CHAR_ARG_LEN (1))); | |
1791 | |
1792 if (info != 0) | |
1793 { | |
1794 rcon = 0.0; | |
1795 | |
1796 mattype.mark_as_unsymmetric (); | |
1797 typ = MatrixType::Full; | |
1798 } | |
1799 else | |
1800 { | |
1801 Array<Complex> z (2 * nc); | |
1802 Complex *pz = z.fortran_vec (); | |
1803 Array<double> rz (nc); | |
1804 double *prz = rz.fortran_vec (); | |
1805 | |
1806 F77_XFCN (zpocon, ZPOCON, (F77_CONST_CHAR_ARG2 (&job, 1), | |
1807 nr, tmp_data, nr, anorm, | |
1808 rcon, pz, prz, info | |
1809 F77_CHAR_ARG_LEN (1))); | |
1810 | |
1811 if (info != 0) | |
1812 rcon = 0.0; | |
1813 } | |
1814 } | |
1815 | |
1816 | |
1817 if (typ == MatrixType::Full) | |
1818 { | |
1819 octave_idx_type info = 0; | |
1820 | |
1821 Array<octave_idx_type> ipvt (nr); | |
1822 octave_idx_type *pipvt = ipvt.fortran_vec (); | |
1823 | |
1824 if(anorm < 0.) | |
1825 anorm = atmp.abs().sum(). | |
1826 row(static_cast<octave_idx_type>(0)).max(); | |
1827 | |
1828 Array<Complex> z (2 * nc); | |
1829 Complex *pz = z.fortran_vec (); | |
1830 Array<double> rz (2 * nc); | |
1831 double *prz = rz.fortran_vec (); | |
1832 | |
1833 F77_XFCN (zgetrf, ZGETRF, (nr, nr, tmp_data, nr, pipvt, info)); | |
1834 | |
1835 if (info != 0) | |
1836 { | |
1837 rcon = 0.0; | |
1838 mattype.mark_as_rectangular (); | |
1839 } | |
1840 else | |
1841 { | |
1842 char job = '1'; | |
1843 F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1), | |
1844 nc, tmp_data, nr, anorm, | |
1845 rcon, pz, prz, info | |
1846 F77_CHAR_ARG_LEN (1))); | |
1847 | |
1848 if (info != 0) | |
1849 rcon = 0.0; | |
1850 } | |
1851 } | |
1852 } | |
1853 else | |
1854 rcon = 0.0; | |
1855 } | |
1856 | |
1857 return rcon; | |
1858 } | |
1859 | |
458 | 1860 ComplexMatrix |
5785 | 1861 ComplexMatrix::utsolve (MatrixType &mattype, const ComplexMatrix& b, |
7788 | 1862 octave_idx_type& info, double& rcon, |
5785 | 1863 solve_singularity_handler sing_handler, |
9661
afcf852256d2
optimize / and '\ for triangular matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
9658
diff
changeset
|
1864 bool calc_cond, blas_trans_type transt) const |
5785 | 1865 { |
1866 ComplexMatrix retval; | |
1867 | |
1868 octave_idx_type nr = rows (); | |
1869 octave_idx_type nc = cols (); | |
1870 | |
6924 | 1871 if (nr != b.rows ()) |
5785 | 1872 (*current_liboctave_error_handler) |
1873 ("matrix dimension mismatch solution of linear equations"); | |
6924 | 1874 else if (nr == 0 || nc == 0 || b.cols () == 0) |
1875 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0)); | |
5785 | 1876 else |
1877 { | |
1878 volatile int typ = mattype.type (); | |
1879 | |
1880 if (typ == MatrixType::Permuted_Upper || | |
1881 typ == MatrixType::Upper) | |
1882 { | |
1883 octave_idx_type b_nc = b.cols (); | |
7788 | 1884 rcon = 1.; |
5785 | 1885 info = 0; |
1886 | |
1887 if (typ == MatrixType::Permuted_Upper) | |
1888 { | |
1889 (*current_liboctave_error_handler) | |
6390 | 1890 ("permuted triangular matrix not implemented"); |
5785 | 1891 } |
1892 else | |
1893 { | |
1894 const Complex *tmp_data = fortran_vec (); | |
1895 | |
1896 if (calc_cond) | |
1897 { | |
1898 char norm = '1'; | |
1899 char uplo = 'U'; | |
1900 char dia = 'N'; | |
1901 | |
1902 Array<Complex> z (2 * nc); | |
1903 Complex *pz = z.fortran_vec (); | |
1904 Array<double> rz (nc); | |
1905 double *prz = rz.fortran_vec (); | |
1906 | |
1907 F77_XFCN (ztrcon, ZTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), | |
1908 F77_CONST_CHAR_ARG2 (&uplo, 1), | |
1909 F77_CONST_CHAR_ARG2 (&dia, 1), | |
7788 | 1910 nr, tmp_data, nr, rcon, |
5785 | 1911 pz, prz, info |
1912 F77_CHAR_ARG_LEN (1) | |
1913 F77_CHAR_ARG_LEN (1) | |
1914 F77_CHAR_ARG_LEN (1))); | |
1915 | |
1916 if (info != 0) | |
1917 info = -2; | |
1918 | |
7788 | 1919 volatile double rcond_plus_one = rcon + 1.0; |
1920 | |
1921 if (rcond_plus_one == 1.0 || xisnan (rcon)) | |
5785 | 1922 { |
1923 info = -2; | |
1924 | |
1925 if (sing_handler) | |
7788 | 1926 sing_handler (rcon); |
5785 | 1927 else |
1928 (*current_liboctave_error_handler) | |
1929 ("matrix singular to machine precision, rcond = %g", | |
7788 | 1930 rcon); |
5785 | 1931 } |
1932 } | |
1933 | |
1934 if (info == 0) | |
1935 { | |
1936 retval = b; | |
1937 Complex *result = retval.fortran_vec (); | |
1938 | |
1939 char uplo = 'U'; | |
9661
afcf852256d2
optimize / and '\ for triangular matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
9658
diff
changeset
|
1940 char trans = get_blas_char (transt); |
5785 | 1941 char dia = 'N'; |
1942 | |
1943 F77_XFCN (ztrtrs, ZTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1), | |
1944 F77_CONST_CHAR_ARG2 (&trans, 1), | |
1945 F77_CONST_CHAR_ARG2 (&dia, 1), | |
1946 nr, b_nc, tmp_data, nr, | |
1947 result, nr, info | |
1948 F77_CHAR_ARG_LEN (1) | |
1949 F77_CHAR_ARG_LEN (1) | |
1950 F77_CHAR_ARG_LEN (1))); | |
1951 } | |
1952 } | |
1953 } | |
1954 else | |
1955 (*current_liboctave_error_handler) ("incorrect matrix type"); | |
1956 } | |
1957 | |
1958 return retval; | |
1959 } | |
1960 | |
1961 ComplexMatrix | |
1962 ComplexMatrix::ltsolve (MatrixType &mattype, const ComplexMatrix& b, | |
7788 | 1963 octave_idx_type& info, double& rcon, |
5785 | 1964 solve_singularity_handler sing_handler, |
9661
afcf852256d2
optimize / and '\ for triangular matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
9658
diff
changeset
|
1965 bool calc_cond, blas_trans_type transt) const |
5785 | 1966 { |
1967 ComplexMatrix retval; | |
1968 | |
1969 octave_idx_type nr = rows (); | |
1970 octave_idx_type nc = cols (); | |
1971 | |
6924 | 1972 if (nr != b.rows ()) |
5785 | 1973 (*current_liboctave_error_handler) |
1974 ("matrix dimension mismatch solution of linear equations"); | |
6924 | 1975 else if (nr == 0 || nc == 0 || b.cols () == 0) |
1976 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0)); | |
5785 | 1977 else |
1978 { | |
1979 volatile int typ = mattype.type (); | |
1980 | |
1981 if (typ == MatrixType::Permuted_Lower || | |
1982 typ == MatrixType::Lower) | |
1983 { | |
1984 octave_idx_type b_nc = b.cols (); | |
7788 | 1985 rcon = 1.; |
5785 | 1986 info = 0; |
1987 | |
1988 if (typ == MatrixType::Permuted_Lower) | |
1989 { | |
1990 (*current_liboctave_error_handler) | |
6390 | 1991 ("permuted triangular matrix not implemented"); |
5785 | 1992 } |
1993 else | |
1994 { | |
1995 const Complex *tmp_data = fortran_vec (); | |
1996 | |
1997 if (calc_cond) | |
1998 { | |
1999 char norm = '1'; | |
2000 char uplo = 'L'; | |
2001 char dia = 'N'; | |
2002 | |
2003 Array<Complex> z (2 * nc); | |
2004 Complex *pz = z.fortran_vec (); | |
2005 Array<double> rz (nc); | |
2006 double *prz = rz.fortran_vec (); | |
2007 | |
2008 F77_XFCN (ztrcon, ZTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), | |
2009 F77_CONST_CHAR_ARG2 (&uplo, 1), | |
2010 F77_CONST_CHAR_ARG2 (&dia, 1), | |
7788 | 2011 nr, tmp_data, nr, rcon, |
5785 | 2012 pz, prz, info |
2013 F77_CHAR_ARG_LEN (1) | |
2014 F77_CHAR_ARG_LEN (1) | |
2015 F77_CHAR_ARG_LEN (1))); | |
2016 | |
2017 if (info != 0) | |
2018 info = -2; | |
2019 | |
7788 | 2020 volatile double rcond_plus_one = rcon + 1.0; |
2021 | |
2022 if (rcond_plus_one == 1.0 || xisnan (rcon)) | |
5785 | 2023 { |
2024 info = -2; | |
2025 | |
2026 if (sing_handler) | |
7788 | 2027 sing_handler (rcon); |
5785 | 2028 else |
2029 (*current_liboctave_error_handler) | |
2030 ("matrix singular to machine precision, rcond = %g", | |
7788 | 2031 rcon); |
5785 | 2032 } |
2033 } | |
2034 | |
2035 if (info == 0) | |
2036 { | |
2037 retval = b; | |
2038 Complex *result = retval.fortran_vec (); | |
2039 | |
2040 char uplo = 'L'; | |
9661
afcf852256d2
optimize / and '\ for triangular matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
9658
diff
changeset
|
2041 char trans = get_blas_char (transt); |
5785 | 2042 char dia = 'N'; |
2043 | |
2044 F77_XFCN (ztrtrs, ZTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1), | |
2045 F77_CONST_CHAR_ARG2 (&trans, 1), | |
2046 F77_CONST_CHAR_ARG2 (&dia, 1), | |
2047 nr, b_nc, tmp_data, nr, | |
2048 result, nr, info | |
2049 F77_CHAR_ARG_LEN (1) | |
2050 F77_CHAR_ARG_LEN (1) | |
2051 F77_CHAR_ARG_LEN (1))); | |
2052 } | |
2053 } | |
2054 } | |
2055 else | |
2056 (*current_liboctave_error_handler) ("incorrect matrix type"); | |
2057 } | |
2058 | |
2059 return retval; | |
2060 } | |
2061 | |
2062 ComplexMatrix | |
2063 ComplexMatrix::fsolve (MatrixType &mattype, const ComplexMatrix& b, | |
7788 | 2064 octave_idx_type& info, double& rcon, |
5785 | 2065 solve_singularity_handler sing_handler, |
2066 bool calc_cond) const | |
2067 { | |
2068 ComplexMatrix retval; | |
2069 | |
2070 octave_idx_type nr = rows (); | |
2071 octave_idx_type nc = cols (); | |
2072 | |
6924 | 2073 |
2074 if (nr != nc || nr != b.rows ()) | |
5785 | 2075 (*current_liboctave_error_handler) |
6924 | 2076 ("matrix dimension mismatch solution of linear equations"); |
2077 else if (nr == 0 || b.cols () == 0) | |
2078 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0)); | |
5785 | 2079 else |
2080 { | |
2081 volatile int typ = mattype.type (); | |
2082 | |
2083 // Calculate the norm of the matrix, for later use. | |
2084 double anorm = -1.; | |
2085 | |
2086 if (typ == MatrixType::Hermitian) | |
2087 { | |
2088 info = 0; | |
2089 char job = 'L'; | |
2090 ComplexMatrix atmp = *this; | |
2091 Complex *tmp_data = atmp.fortran_vec (); | |
2092 anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max(); | |
2093 | |
2094 F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr, | |
2095 tmp_data, nr, info | |
2096 F77_CHAR_ARG_LEN (1))); | |
2097 | |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2098 // Throw-away extra info LAPACK gives so as to not change output. |
7788 | 2099 rcon = 0.0; |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2100 if (info != 0) |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2101 { |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2102 info = -2; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2103 |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2104 mattype.mark_as_unsymmetric (); |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2105 typ = MatrixType::Full; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2106 } |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2107 else |
5785 | 2108 { |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2109 if (calc_cond) |
5785 | 2110 { |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2111 Array<Complex> z (2 * nc); |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2112 Complex *pz = z.fortran_vec (); |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2113 Array<double> rz (nc); |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2114 double *prz = rz.fortran_vec (); |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2115 |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2116 F77_XFCN (zpocon, ZPOCON, (F77_CONST_CHAR_ARG2 (&job, 1), |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2117 nr, tmp_data, nr, anorm, |
7788 | 2118 rcon, pz, prz, info |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2119 F77_CHAR_ARG_LEN (1))); |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2120 |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2121 if (info != 0) |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2122 info = -2; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2123 |
7788 | 2124 volatile double rcond_plus_one = rcon + 1.0; |
2125 | |
2126 if (rcond_plus_one == 1.0 || xisnan (rcon)) | |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2127 { |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2128 info = -2; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2129 |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2130 if (sing_handler) |
7788 | 2131 sing_handler (rcon); |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2132 else |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2133 (*current_liboctave_error_handler) |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2134 ("matrix singular to machine precision, rcond = %g", |
7788 | 2135 rcon); |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2136 } |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2137 } |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2138 |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2139 if (info == 0) |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2140 { |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2141 retval = b; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2142 Complex *result = retval.fortran_vec (); |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2143 |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2144 octave_idx_type b_nc = b.cols (); |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2145 |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2146 F77_XFCN (zpotrs, ZPOTRS, (F77_CONST_CHAR_ARG2 (&job, 1), |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2147 nr, b_nc, tmp_data, nr, |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2148 result, b.rows(), info |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2149 F77_CHAR_ARG_LEN (1))); |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2150 } |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2151 else |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2152 { |
5785 | 2153 mattype.mark_as_unsymmetric (); |
2154 typ = MatrixType::Full; | |
2155 } | |
2156 } | |
2157 } | |
2158 | |
2159 if (typ == MatrixType::Full) | |
2160 { | |
2161 info = 0; | |
2162 | |
2163 Array<octave_idx_type> ipvt (nr); | |
2164 octave_idx_type *pipvt = ipvt.fortran_vec (); | |
2165 | |
2166 ComplexMatrix atmp = *this; | |
2167 Complex *tmp_data = atmp.fortran_vec (); | |
2168 | |
2169 Array<Complex> z (2 * nc); | |
2170 Complex *pz = z.fortran_vec (); | |
2171 Array<double> rz (2 * nc); | |
2172 double *prz = rz.fortran_vec (); | |
2173 | |
2174 // Calculate the norm of the matrix, for later use. | |
2175 if (anorm < 0.) | |
2176 anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max(); | |
2177 | |
2178 F77_XFCN (zgetrf, ZGETRF, (nr, nr, tmp_data, nr, pipvt, info)); | |
2179 | |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2180 // Throw-away extra info LAPACK gives so as to not change output. |
7788 | 2181 rcon = 0.0; |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2182 if (info != 0) |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2183 { |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2184 info = -2; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2185 |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2186 if (sing_handler) |
7788 | 2187 sing_handler (rcon); |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2188 else |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2189 (*current_liboctave_error_handler) |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2190 ("matrix singular to machine precision"); |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2191 |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2192 mattype.mark_as_rectangular (); |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2193 } |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2194 else |
5785 | 2195 { |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2196 if (calc_cond) |
5785 | 2197 { |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2198 // Now calculate the condition number for |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2199 // non-singular matrix. |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2200 char job = '1'; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2201 F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1), |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2202 nc, tmp_data, nr, anorm, |
7788 | 2203 rcon, pz, prz, info |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2204 F77_CHAR_ARG_LEN (1))); |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2205 |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2206 if (info != 0) |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2207 info = -2; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2208 |
7788 | 2209 volatile double rcond_plus_one = rcon + 1.0; |
2210 | |
2211 if (rcond_plus_one == 1.0 || xisnan (rcon)) | |
5785 | 2212 { |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2213 info = -2; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2214 |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2215 if (sing_handler) |
7788 | 2216 sing_handler (rcon); |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2217 else |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2218 (*current_liboctave_error_handler) |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2219 ("matrix singular to machine precision, rcond = %g", |
7788 | 2220 rcon); |
5785 | 2221 } |
2222 } | |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2223 |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2224 if (info == 0) |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2225 { |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2226 retval = b; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2227 Complex *result = retval.fortran_vec (); |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2228 |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2229 octave_idx_type b_nc = b.cols (); |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2230 |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2231 char job = 'N'; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2232 F77_XFCN (zgetrs, ZGETRS, (F77_CONST_CHAR_ARG2 (&job, 1), |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2233 nr, b_nc, tmp_data, nr, |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2234 pipvt, result, b.rows(), info |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2235 F77_CHAR_ARG_LEN (1))); |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2236 } |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2237 else |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2238 mattype.mark_as_rectangular (); |
5785 | 2239 } |
2240 } | |
2241 } | |
2242 | |
2243 return retval; | |
2244 } | |
2245 | |
2246 ComplexMatrix | |
2247 ComplexMatrix::solve (MatrixType &typ, const Matrix& b) const | |
2248 { | |
2249 octave_idx_type info; | |
7788 | 2250 double rcon; |
2251 return solve (typ, b, info, rcon, 0); | |
5785 | 2252 } |
2253 | |
2254 ComplexMatrix | |
2255 ComplexMatrix::solve (MatrixType &typ, const Matrix& b, | |
2256 octave_idx_type& info) const | |
2257 { | |
7788 | 2258 double rcon; |
2259 return solve (typ, b, info, rcon, 0); | |
5785 | 2260 } |
2261 | |
2262 ComplexMatrix | |
2263 ComplexMatrix::solve (MatrixType &typ, const Matrix& b, octave_idx_type& info, | |
7788 | 2264 double& rcon) const |
5785 | 2265 { |
7788 | 2266 return solve (typ, b, info, rcon, 0); |
5785 | 2267 } |
2268 | |
2269 ComplexMatrix | |
2270 ComplexMatrix::solve (MatrixType &typ, const Matrix& b, octave_idx_type& info, | |
7788 | 2271 double& rcon, solve_singularity_handler sing_handler, |
9661
afcf852256d2
optimize / and '\ for triangular matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
9658
diff
changeset
|
2272 bool singular_fallback, blas_trans_type transt) const |
5785 | 2273 { |
2274 ComplexMatrix tmp (b); | |
9661
afcf852256d2
optimize / and '\ for triangular matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
9658
diff
changeset
|
2275 return solve (typ, tmp, info, rcon, sing_handler, singular_fallback, transt); |
5785 | 2276 } |
2277 | |
2278 ComplexMatrix | |
2279 ComplexMatrix::solve (MatrixType &typ, const ComplexMatrix& b) const | |
2280 { | |
2281 octave_idx_type info; | |
7788 | 2282 double rcon; |
2283 return solve (typ, b, info, rcon, 0); | |
5785 | 2284 } |
2285 | |
2286 ComplexMatrix | |
2287 ComplexMatrix::solve (MatrixType &typ, const ComplexMatrix& b, | |
2288 octave_idx_type& info) const | |
2289 { | |
7788 | 2290 double rcon; |
2291 return solve (typ, b, info, rcon, 0); | |
5785 | 2292 } |
2293 | |
2294 ComplexMatrix | |
2295 ComplexMatrix::solve (MatrixType &typ, const ComplexMatrix& b, | |
7788 | 2296 octave_idx_type& info, double& rcon) const |
5785 | 2297 { |
7788 | 2298 return solve (typ, b, info, rcon, 0); |
5785 | 2299 } |
2300 | |
2301 ComplexMatrix | |
2302 ComplexMatrix::solve (MatrixType &mattype, const ComplexMatrix& b, | |
7788 | 2303 octave_idx_type& info, double& rcon, |
5785 | 2304 solve_singularity_handler sing_handler, |
9661
afcf852256d2
optimize / and '\ for triangular matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
9658
diff
changeset
|
2305 bool singular_fallback, blas_trans_type transt) const |
5785 | 2306 { |
2307 ComplexMatrix retval; | |
2308 int typ = mattype.type (); | |
2309 | |
2310 if (typ == MatrixType::Unknown) | |
2311 typ = mattype.type (*this); | |
2312 | |
2313 // Only calculate the condition number for LU/Cholesky | |
2314 if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper) | |
9661
afcf852256d2
optimize / and '\ for triangular matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
9658
diff
changeset
|
2315 retval = utsolve (mattype, b, info, rcon, sing_handler, false, transt); |
5785 | 2316 else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower) |
9661
afcf852256d2
optimize / and '\ for triangular matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
9658
diff
changeset
|
2317 retval = ltsolve (mattype, b, info, rcon, sing_handler, false, transt); |
afcf852256d2
optimize / and '\ for triangular matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
9658
diff
changeset
|
2318 else if (transt == blas_trans) |
afcf852256d2
optimize / and '\ for triangular matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
9658
diff
changeset
|
2319 return transpose ().solve (mattype, b, info, rcon, sing_handler, singular_fallback); |
afcf852256d2
optimize / and '\ for triangular matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
9658
diff
changeset
|
2320 else if (transt == blas_conj_trans) |
afcf852256d2
optimize / and '\ for triangular matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
9658
diff
changeset
|
2321 retval = hermitian ().solve (mattype, b, info, rcon, sing_handler, singular_fallback); |
5785 | 2322 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian) |
7788 | 2323 retval = fsolve (mattype, b, info, rcon, sing_handler, true); |
5785 | 2324 else if (typ != MatrixType::Rectangular) |
2325 { | |
2326 (*current_liboctave_error_handler) ("unknown matrix type"); | |
2327 return ComplexMatrix (); | |
2328 } | |
2329 | |
2330 // Rectangular or one of the above solvers flags a singular matrix | |
2331 if (singular_fallback && mattype.type () == MatrixType::Rectangular) | |
2332 { | |
2333 octave_idx_type rank; | |
7788 | 2334 retval = lssolve (b, info, rank, rcon); |
5785 | 2335 } |
2336 | |
2337 return retval; | |
2338 } | |
2339 | |
2340 ComplexColumnVector | |
2341 ComplexMatrix::solve (MatrixType &typ, const ColumnVector& b) const | |
2342 { | |
2343 octave_idx_type info; | |
7788 | 2344 double rcon; |
2345 return solve (typ, ComplexColumnVector (b), info, rcon, 0); | |
5785 | 2346 } |
2347 | |
2348 ComplexColumnVector | |
2349 ComplexMatrix::solve (MatrixType &typ, const ColumnVector& b, | |
2350 octave_idx_type& info) const | |
2351 { | |
7788 | 2352 double rcon; |
2353 return solve (typ, ComplexColumnVector (b), info, rcon, 0); | |
5785 | 2354 } |
2355 | |
2356 ComplexColumnVector | |
2357 ComplexMatrix::solve (MatrixType &typ, const ColumnVector& b, | |
7788 | 2358 octave_idx_type& info, double& rcon) const |
5785 | 2359 { |
7788 | 2360 return solve (typ, ComplexColumnVector (b), info, rcon, 0); |
5785 | 2361 } |
2362 | |
2363 ComplexColumnVector | |
2364 ComplexMatrix::solve (MatrixType &typ, const ColumnVector& b, | |
7788 | 2365 octave_idx_type& info, double& rcon, |
9661
afcf852256d2
optimize / and '\ for triangular matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
9658
diff
changeset
|
2366 solve_singularity_handler sing_handler, blas_trans_type transt) const |
5785 | 2367 { |
9661
afcf852256d2
optimize / and '\ for triangular matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
9658
diff
changeset
|
2368 return solve (typ, ComplexColumnVector (b), info, rcon, sing_handler, transt); |
5785 | 2369 } |
2370 | |
2371 ComplexColumnVector | |
2372 ComplexMatrix::solve (MatrixType &typ, const ComplexColumnVector& b) const | |
2373 { | |
2374 octave_idx_type info; | |
7788 | 2375 double rcon; |
2376 return solve (typ, b, info, rcon, 0); | |
5785 | 2377 } |
2378 | |
2379 ComplexColumnVector | |
2380 ComplexMatrix::solve (MatrixType &typ, const ComplexColumnVector& b, | |
2381 octave_idx_type& info) const | |
2382 { | |
7788 | 2383 double rcon; |
2384 return solve (typ, b, info, rcon, 0); | |
5785 | 2385 } |
2386 | |
2387 ComplexColumnVector | |
2388 ComplexMatrix::solve (MatrixType &typ, const ComplexColumnVector& b, | |
7788 | 2389 octave_idx_type& info, double& rcon) const |
5785 | 2390 { |
7788 | 2391 return solve (typ, b, info, rcon, 0); |
5785 | 2392 } |
2393 | |
2394 ComplexColumnVector | |
2395 ComplexMatrix::solve (MatrixType &typ, const ComplexColumnVector& b, | |
7788 | 2396 octave_idx_type& info, double& rcon, |
9661
afcf852256d2
optimize / and '\ for triangular matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
9658
diff
changeset
|
2397 solve_singularity_handler sing_handler, blas_trans_type transt) const |
5785 | 2398 { |
2399 | |
2400 ComplexMatrix tmp (b); | |
9661
afcf852256d2
optimize / and '\ for triangular matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
9658
diff
changeset
|
2401 return solve (typ, tmp, info, rcon, sing_handler, transt).column(static_cast<octave_idx_type> (0)); |
5785 | 2402 } |
2403 | |
2404 ComplexMatrix | |
458 | 2405 ComplexMatrix::solve (const Matrix& b) const |
2406 { | |
5275 | 2407 octave_idx_type info; |
7788 | 2408 double rcon; |
2409 return solve (b, info, rcon, 0); | |
458 | 2410 } |
2411 | |
2412 ComplexMatrix | |
5275 | 2413 ComplexMatrix::solve (const Matrix& b, octave_idx_type& info) const |
458 | 2414 { |
7788 | 2415 double rcon; |
2416 return solve (b, info, rcon, 0); | |
458 | 2417 } |
2418 | |
2419 ComplexMatrix | |
7788 | 2420 ComplexMatrix::solve (const Matrix& b, octave_idx_type& info, double& rcon) const |
458 | 2421 { |
7788 | 2422 return solve (b, info, rcon, 0); |
3480 | 2423 } |
2424 | |
2425 ComplexMatrix | |
7788 | 2426 ComplexMatrix::solve (const Matrix& b, octave_idx_type& info, double& rcon, |
9661
afcf852256d2
optimize / and '\ for triangular matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
9658
diff
changeset
|
2427 solve_singularity_handler sing_handler, blas_trans_type transt) const |
3480 | 2428 { |
458 | 2429 ComplexMatrix tmp (b); |
9661
afcf852256d2
optimize / and '\ for triangular matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
9658
diff
changeset
|
2430 return solve (tmp, info, rcon, sing_handler, transt); |
458 | 2431 } |
2432 | |
2433 ComplexMatrix | |
2434 ComplexMatrix::solve (const ComplexMatrix& b) const | |
2435 { | |
5275 | 2436 octave_idx_type info; |
7788 | 2437 double rcon; |
2438 return solve (b, info, rcon, 0); | |
458 | 2439 } |
2440 | |
2441 ComplexMatrix | |
5275 | 2442 ComplexMatrix::solve (const ComplexMatrix& b, octave_idx_type& info) const |
458 | 2443 { |
7788 | 2444 double rcon; |
2445 return solve (b, info, rcon, 0); | |
458 | 2446 } |
3480 | 2447 |
458 | 2448 ComplexMatrix |
7788 | 2449 ComplexMatrix::solve (const ComplexMatrix& b, octave_idx_type& info, double& rcon) const |
458 | 2450 { |
7788 | 2451 return solve (b, info, rcon, 0); |
3480 | 2452 } |
2453 | |
2454 ComplexMatrix | |
7788 | 2455 ComplexMatrix::solve (const ComplexMatrix& b, octave_idx_type& info, double& rcon, |
9661
afcf852256d2
optimize / and '\ for triangular matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
9658
diff
changeset
|
2456 solve_singularity_handler sing_handler, blas_trans_type transt) const |
3480 | 2457 { |
5785 | 2458 MatrixType mattype (*this); |
9661
afcf852256d2
optimize / and '\ for triangular matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
9658
diff
changeset
|
2459 return solve (mattype, b, info, rcon, sing_handler, true, transt); |
458 | 2460 } |
2461 | |
2462 ComplexColumnVector | |
3585 | 2463 ComplexMatrix::solve (const ColumnVector& b) const |
2464 { | |
5275 | 2465 octave_idx_type info; |
7788 | 2466 double rcon; |
2467 return solve (ComplexColumnVector (b), info, rcon, 0); | |
3585 | 2468 } |
2469 | |
2470 ComplexColumnVector | |
5275 | 2471 ComplexMatrix::solve (const ColumnVector& b, octave_idx_type& info) const |
3585 | 2472 { |
7788 | 2473 double rcon; |
2474 return solve (ComplexColumnVector (b), info, rcon, 0); | |
3585 | 2475 } |
2476 | |
2477 ComplexColumnVector | |
5785 | 2478 ComplexMatrix::solve (const ColumnVector& b, octave_idx_type& info, |
7788 | 2479 double& rcon) const |
3585 | 2480 { |
7788 | 2481 return solve (ComplexColumnVector (b), info, rcon, 0); |
3585 | 2482 } |
2483 | |
2484 ComplexColumnVector | |
5785 | 2485 ComplexMatrix::solve (const ColumnVector& b, octave_idx_type& info, |
7788 | 2486 double& rcon, |
9661
afcf852256d2
optimize / and '\ for triangular matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
9658
diff
changeset
|
2487 solve_singularity_handler sing_handler, blas_trans_type transt) const |
3585 | 2488 { |
9661
afcf852256d2
optimize / and '\ for triangular matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
9658
diff
changeset
|
2489 return solve (ComplexColumnVector (b), info, rcon, sing_handler, transt); |
3585 | 2490 } |
2491 | |
2492 ComplexColumnVector | |
458 | 2493 ComplexMatrix::solve (const ComplexColumnVector& b) const |
2494 { | |
5275 | 2495 octave_idx_type info; |
7788 | 2496 double rcon; |
2497 return solve (b, info, rcon, 0); | |
458 | 2498 } |
2499 | |
2500 ComplexColumnVector | |
5275 | 2501 ComplexMatrix::solve (const ComplexColumnVector& b, octave_idx_type& info) const |
458 | 2502 { |
7788 | 2503 double rcon; |
2504 return solve (b, info, rcon, 0); | |
458 | 2505 } |
2506 | |
2507 ComplexColumnVector | |
5275 | 2508 ComplexMatrix::solve (const ComplexColumnVector& b, octave_idx_type& info, |
7788 | 2509 double& rcon) const |
458 | 2510 { |
7788 | 2511 return solve (b, info, rcon, 0); |
3480 | 2512 } |
2513 | |
2514 ComplexColumnVector | |
5275 | 2515 ComplexMatrix::solve (const ComplexColumnVector& b, octave_idx_type& info, |
7788 | 2516 double& rcon, |
9661
afcf852256d2
optimize / and '\ for triangular matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
9658
diff
changeset
|
2517 solve_singularity_handler sing_handler, blas_trans_type transt) const |
3480 | 2518 { |
5785 | 2519 MatrixType mattype (*this); |
9661
afcf852256d2
optimize / and '\ for triangular matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
9658
diff
changeset
|
2520 return solve (mattype, b, info, rcon, sing_handler, transt); |
458 | 2521 } |
2522 | |
2523 ComplexMatrix | |
3585 | 2524 ComplexMatrix::lssolve (const Matrix& b) const |
2525 { | |
5275 | 2526 octave_idx_type info; |
2527 octave_idx_type rank; | |
7788 | 2528 double rcon; |
2529 return lssolve (ComplexMatrix (b), info, rank, rcon); | |
3585 | 2530 } |
2531 | |
2532 ComplexMatrix | |
5275 | 2533 ComplexMatrix::lssolve (const Matrix& b, octave_idx_type& info) const |
3585 | 2534 { |
5275 | 2535 octave_idx_type rank; |
7788 | 2536 double rcon; |
2537 return lssolve (ComplexMatrix (b), info, rank, rcon); | |
3585 | 2538 } |
2539 | |
2540 ComplexMatrix | |
7076 | 2541 ComplexMatrix::lssolve (const Matrix& b, octave_idx_type& info, |
2542 octave_idx_type& rank) const | |
3585 | 2543 { |
7788 | 2544 double rcon; |
2545 return lssolve (ComplexMatrix (b), info, rank, rcon); | |
7076 | 2546 } |
2547 | |
2548 ComplexMatrix | |
2549 ComplexMatrix::lssolve (const Matrix& b, octave_idx_type& info, | |
7788 | 2550 octave_idx_type& rank, double& rcon) const |
7076 | 2551 { |
7788 | 2552 return lssolve (ComplexMatrix (b), info, rank, rcon); |
3585 | 2553 } |
2554 | |
2555 ComplexMatrix | |
458 | 2556 ComplexMatrix::lssolve (const ComplexMatrix& b) const |
2557 { | |
5275 | 2558 octave_idx_type info; |
2559 octave_idx_type rank; | |
7788 | 2560 double rcon; |
2561 return lssolve (b, info, rank, rcon); | |
458 | 2562 } |
2563 | |
2564 ComplexMatrix | |
5275 | 2565 ComplexMatrix::lssolve (const ComplexMatrix& b, octave_idx_type& info) const |
458 | 2566 { |
5275 | 2567 octave_idx_type rank; |
7788 | 2568 double rcon; |
2569 return lssolve (b, info, rank, rcon); | |
458 | 2570 } |
2571 | |
2572 ComplexMatrix | |
7076 | 2573 ComplexMatrix::lssolve (const ComplexMatrix& b, octave_idx_type& info, |
2574 octave_idx_type& rank) const | |
2575 { | |
7788 | 2576 double rcon; |
2577 return lssolve (b, info, rank, rcon); | |
7076 | 2578 } |
2579 | |
2580 ComplexMatrix | |
2581 ComplexMatrix::lssolve (const ComplexMatrix& b, octave_idx_type& info, | |
7788 | 2582 octave_idx_type& rank, double& rcon) const |
458 | 2583 { |
1948 | 2584 ComplexMatrix retval; |
2585 | |
5275 | 2586 octave_idx_type nrhs = b.cols (); |
2587 | |
2588 octave_idx_type m = rows (); | |
2589 octave_idx_type n = cols (); | |
458 | 2590 |
6924 | 2591 if (m != b.rows ()) |
1948 | 2592 (*current_liboctave_error_handler) |
2593 ("matrix dimension mismatch solution of linear equations"); | |
6924 | 2594 else if (m== 0 || n == 0 || b.cols () == 0) |
2595 retval = ComplexMatrix (n, b.cols (), Complex (0.0, 0.0)); | |
1948 | 2596 else |
458 | 2597 { |
7072 | 2598 volatile octave_idx_type minmn = (m < n ? m : n); |
2599 octave_idx_type maxmn = m > n ? m : n; | |
7788 | 2600 rcon = -1.0; |
7072 | 2601 |
2602 if (m != n) | |
2603 { | |
2604 retval = ComplexMatrix (maxmn, nrhs); | |
2605 | |
2606 for (octave_idx_type j = 0; j < nrhs; j++) | |
2607 for (octave_idx_type i = 0; i < m; i++) | |
2608 retval.elem (i, j) = b.elem (i, j); | |
2609 } | |
2610 else | |
2611 retval = b; | |
2612 | |
1948 | 2613 ComplexMatrix atmp = *this; |
2614 Complex *tmp_data = atmp.fortran_vec (); | |
2615 | |
7072 | 2616 Complex *pretval = retval.fortran_vec (); |
2617 Array<double> s (minmn); | |
7071 | 2618 double *ps = s.fortran_vec (); |
2563 | 2619 |
7072 | 2620 // Ask ZGELSD what the dimension of WORK should be. |
5275 | 2621 octave_idx_type lwork = -1; |
3752 | 2622 |
2623 Array<Complex> work (1); | |
7079 | 2624 |
7477 | 2625 octave_idx_type smlsiz; |
2626 F77_FUNC (xilaenv, XILAENV) (9, F77_CONST_CHAR_ARG2 ("ZGELSD", 6), | |
2627 F77_CONST_CHAR_ARG2 (" ", 1), | |
7478 | 2628 0, 0, 0, 0, smlsiz |
7477 | 2629 F77_CHAR_ARG_LEN (6) |
7478 | 2630 F77_CHAR_ARG_LEN (1)); |
7079 | 2631 |
7486
6a6d2abe51ff
more xGELSD workspace fixes
John W. Eaton <jwe@octave.org>
parents:
7482
diff
changeset
|
2632 octave_idx_type mnthr; |
6a6d2abe51ff
more xGELSD workspace fixes
John W. Eaton <jwe@octave.org>
parents:
7482
diff
changeset
|
2633 F77_FUNC (xilaenv, XILAENV) (6, F77_CONST_CHAR_ARG2 ("ZGELSD", 6), |
6a6d2abe51ff
more xGELSD workspace fixes
John W. Eaton <jwe@octave.org>
parents:
7482
diff
changeset
|
2634 F77_CONST_CHAR_ARG2 (" ", 1), |
6a6d2abe51ff
more xGELSD workspace fixes
John W. Eaton <jwe@octave.org>
parents:
7482
diff
changeset
|
2635 m, n, nrhs, -1, mnthr |
6a6d2abe51ff
more xGELSD workspace fixes
John W. Eaton <jwe@octave.org>
parents:
7482
diff
changeset
|
2636 F77_CHAR_ARG_LEN (6) |
6a6d2abe51ff
more xGELSD workspace fixes
John W. Eaton <jwe@octave.org>
parents:
7482
diff
changeset
|
2637 F77_CHAR_ARG_LEN (1)); |
6a6d2abe51ff
more xGELSD workspace fixes
John W. Eaton <jwe@octave.org>
parents:
7482
diff
changeset
|
2638 |
7079 | 2639 // We compute the size of rwork and iwork because ZGELSD in |
2640 // older versions of LAPACK does not return them on a query | |
2641 // call. | |
7124 | 2642 double dminmn = static_cast<double> (minmn); |
2643 double dsmlsizp1 = static_cast<double> (smlsiz+1); | |
7079 | 2644 #if defined (HAVE_LOG2) |
7544
f9983d2761df
more xGELSD workspace fixes
Jaroslav Hajek <highegg@gmail.com>
parents:
7532
diff
changeset
|
2645 double tmp = log2 (dminmn / dsmlsizp1); |
7079 | 2646 #else |
7544
f9983d2761df
more xGELSD workspace fixes
Jaroslav Hajek <highegg@gmail.com>
parents:
7532
diff
changeset
|
2647 double tmp = log (dminmn / dsmlsizp1) / log (2.0); |
7079 | 2648 #endif |
7544
f9983d2761df
more xGELSD workspace fixes
Jaroslav Hajek <highegg@gmail.com>
parents:
7532
diff
changeset
|
2649 octave_idx_type nlvl = static_cast<octave_idx_type> (tmp) + 1; |
7079 | 2650 if (nlvl < 0) |
2651 nlvl = 0; | |
2652 | |
2653 octave_idx_type lrwork = minmn*(10 + 2*smlsiz + 8*nlvl) | |
2654 + 3*smlsiz*nrhs + (smlsiz+1)*(smlsiz+1); | |
2655 if (lrwork < 1) | |
2656 lrwork = 1; | |
2657 Array<double> rwork (lrwork); | |
2658 double *prwork = rwork.fortran_vec (); | |
2659 | |
2660 octave_idx_type liwork = 3 * minmn * nlvl + 11 * minmn; | |
2661 if (liwork < 1) | |
2662 liwork = 1; | |
2663 Array<octave_idx_type> iwork (liwork); | |
2664 octave_idx_type* piwork = iwork.fortran_vec (); | |
7072 | 2665 |
2666 F77_XFCN (zgelsd, ZGELSD, (m, n, nrhs, tmp_data, m, pretval, maxmn, | |
7788 | 2667 ps, rcon, rank, work.fortran_vec (), |
7079 | 2668 lwork, prwork, piwork, info)); |
1948 | 2669 |
7476 | 2670 // The workspace query is broken in at least LAPACK 3.0.0 |
7488
6470f946a425
another small xGELSD workspace fix
John W. Eaton <jwe@octave.org>
parents:
7486
diff
changeset
|
2671 // through 3.1.1 when n >= mnthr. The obtuse formula below |
7486
6a6d2abe51ff
more xGELSD workspace fixes
John W. Eaton <jwe@octave.org>
parents:
7482
diff
changeset
|
2672 // should provide sufficient workspace for ZGELSD to operate |
7476 | 2673 // efficiently. |
7488
6470f946a425
another small xGELSD workspace fix
John W. Eaton <jwe@octave.org>
parents:
7486
diff
changeset
|
2674 if (n >= mnthr) |
7476 | 2675 { |
2676 octave_idx_type addend = m; | |
2677 | |
2678 if (2*m-4 > addend) | |
2679 addend = 2*m-4; | |
2680 | |
2681 if (nrhs > addend) | |
2682 addend = nrhs; | |
2683 | |
2684 if (n-3*m > addend) | |
2685 addend = n-3*m; | |
2686 | |
2687 const octave_idx_type lworkaround = 4*m + m*m + addend; | |
2688 | |
2689 if (std::real (work(0)) < lworkaround) | |
2690 work(0) = lworkaround; | |
2691 } | |
7532
493bb0de3199
avoid another xGELSD workspace query bug
John W. Eaton <jwe@octave.org>
parents:
7503
diff
changeset
|
2692 else if (m >= n) |
493bb0de3199
avoid another xGELSD workspace query bug
John W. Eaton <jwe@octave.org>
parents:
7503
diff
changeset
|
2693 { |
493bb0de3199
avoid another xGELSD workspace query bug
John W. Eaton <jwe@octave.org>
parents:
7503
diff
changeset
|
2694 octave_idx_type lworkaround = 2*m + m*nrhs; |
493bb0de3199
avoid another xGELSD workspace query bug
John W. Eaton <jwe@octave.org>
parents:
7503
diff
changeset
|
2695 |
493bb0de3199
avoid another xGELSD workspace query bug
John W. Eaton <jwe@octave.org>
parents:
7503
diff
changeset
|
2696 if (std::real (work(0)) < lworkaround) |
493bb0de3199
avoid another xGELSD workspace query bug
John W. Eaton <jwe@octave.org>
parents:
7503
diff
changeset
|
2697 work(0) = lworkaround; |
493bb0de3199
avoid another xGELSD workspace query bug
John W. Eaton <jwe@octave.org>
parents:
7503
diff
changeset
|
2698 } |
7476 | 2699 |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2700 lwork = static_cast<octave_idx_type> (std::real (work(0))); |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2701 work.resize (lwork); |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2702 |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2703 F77_XFCN (zgelsd, ZGELSD, (m, n, nrhs, tmp_data, m, pretval, |
7788 | 2704 maxmn, ps, rcon, rank, |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2705 work.fortran_vec (), lwork, |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2706 prwork, piwork, info)); |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2707 |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2708 if (rank < minmn) |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2709 (*current_liboctave_warning_handler) |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2710 ("zgelsd: rank deficient %dx%d matrix, rank = %d, tol = %e", |
7788 | 2711 m, n, rank, rcon); |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2712 |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2713 if (s.elem (0) == 0.0) |
7788 | 2714 rcon = 0.0; |
1948 | 2715 else |
7788 | 2716 rcon = s.elem (minmn - 1) / s.elem (0); |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2717 |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2718 retval.resize (n, nrhs); |
458 | 2719 } |
2720 | |
2721 return retval; | |
2722 } | |
2723 | |
2724 ComplexColumnVector | |
3585 | 2725 ComplexMatrix::lssolve (const ColumnVector& b) const |
2726 { | |
5275 | 2727 octave_idx_type info; |
2728 octave_idx_type rank; | |
7788 | 2729 double rcon; |
2730 return lssolve (ComplexColumnVector (b), info, rank, rcon); | |
3585 | 2731 } |
2732 | |
2733 ComplexColumnVector | |
5275 | 2734 ComplexMatrix::lssolve (const ColumnVector& b, octave_idx_type& info) const |
3585 | 2735 { |
5275 | 2736 octave_idx_type rank; |
7788 | 2737 double rcon; |
2738 return lssolve (ComplexColumnVector (b), info, rank, rcon); | |
3585 | 2739 } |
2740 | |
2741 ComplexColumnVector | |
7076 | 2742 ComplexMatrix::lssolve (const ColumnVector& b, octave_idx_type& info, |
2743 octave_idx_type& rank) const | |
3585 | 2744 { |
7788 | 2745 double rcon; |
2746 return lssolve (ComplexColumnVector (b), info, rank, rcon); | |
7076 | 2747 } |
2748 | |
2749 ComplexColumnVector | |
2750 ComplexMatrix::lssolve (const ColumnVector& b, octave_idx_type& info, | |
7788 | 2751 octave_idx_type& rank, double& rcon) const |
7076 | 2752 { |
7788 | 2753 return lssolve (ComplexColumnVector (b), info, rank, rcon); |
3585 | 2754 } |
2755 | |
2756 ComplexColumnVector | |
458 | 2757 ComplexMatrix::lssolve (const ComplexColumnVector& b) const |
2758 { | |
5275 | 2759 octave_idx_type info; |
2760 octave_idx_type rank; | |
7788 | 2761 double rcon; |
2762 return lssolve (b, info, rank, rcon); | |
458 | 2763 } |
2764 | |
2765 ComplexColumnVector | |
5275 | 2766 ComplexMatrix::lssolve (const ComplexColumnVector& b, octave_idx_type& info) const |
458 | 2767 { |
5275 | 2768 octave_idx_type rank; |
7788 | 2769 double rcon; |
2770 return lssolve (b, info, rank, rcon); | |
458 | 2771 } |
2772 | |
2773 ComplexColumnVector | |
5275 | 2774 ComplexMatrix::lssolve (const ComplexColumnVector& b, octave_idx_type& info, |
2775 octave_idx_type& rank) const | |
458 | 2776 { |
7788 | 2777 double rcon; |
2778 return lssolve (b, info, rank, rcon); | |
7076 | 2779 |
2780 } | |
2781 | |
2782 ComplexColumnVector | |
2783 ComplexMatrix::lssolve (const ComplexColumnVector& b, octave_idx_type& info, | |
7788 | 2784 octave_idx_type& rank, double& rcon) const |
7076 | 2785 { |
1948 | 2786 ComplexColumnVector retval; |
2787 | |
5275 | 2788 octave_idx_type nrhs = 1; |
2789 | |
2790 octave_idx_type m = rows (); | |
2791 octave_idx_type n = cols (); | |
458 | 2792 |
6924 | 2793 if (m != b.length ()) |
1948 | 2794 (*current_liboctave_error_handler) |
6924 | 2795 ("matrix dimension mismatch solution of linear equations"); |
2796 else if (m == 0 || n == 0 || b.cols () == 0) | |
2797 retval = ComplexColumnVector (n, Complex (0.0, 0.0)); | |
1948 | 2798 else |
458 | 2799 { |
7072 | 2800 volatile octave_idx_type minmn = (m < n ? m : n); |
2801 octave_idx_type maxmn = m > n ? m : n; | |
7788 | 2802 rcon = -1.0; |
7072 | 2803 |
2804 if (m != n) | |
2805 { | |
2806 retval = ComplexColumnVector (maxmn); | |
2807 | |
2808 for (octave_idx_type i = 0; i < m; i++) | |
2809 retval.elem (i) = b.elem (i); | |
2810 } | |
2811 else | |
2812 retval = b; | |
2813 | |
1948 | 2814 ComplexMatrix atmp = *this; |
2815 Complex *tmp_data = atmp.fortran_vec (); | |
2816 | |
7072 | 2817 Complex *pretval = retval.fortran_vec (); |
2818 Array<double> s (minmn); | |
7071 | 2819 double *ps = s.fortran_vec (); |
1948 | 2820 |
7072 | 2821 // Ask ZGELSD what the dimension of WORK should be. |
5275 | 2822 octave_idx_type lwork = -1; |
3752 | 2823 |
2824 Array<Complex> work (1); | |
7079 | 2825 |
7544
f9983d2761df
more xGELSD workspace fixes
Jaroslav Hajek <highegg@gmail.com>
parents:
7532
diff
changeset
|
2826 octave_idx_type smlsiz; |
f9983d2761df
more xGELSD workspace fixes
Jaroslav Hajek <highegg@gmail.com>
parents:
7532
diff
changeset
|
2827 F77_FUNC (xilaenv, XILAENV) (9, F77_CONST_CHAR_ARG2 ("ZGELSD", 6), |
f9983d2761df
more xGELSD workspace fixes
Jaroslav Hajek <highegg@gmail.com>
parents:
7532
diff
changeset
|
2828 F77_CONST_CHAR_ARG2 (" ", 1), |
f9983d2761df
more xGELSD workspace fixes
Jaroslav Hajek <highegg@gmail.com>
parents:
7532
diff
changeset
|
2829 0, 0, 0, 0, smlsiz |
f9983d2761df
more xGELSD workspace fixes
Jaroslav Hajek <highegg@gmail.com>
parents:
7532
diff
changeset
|
2830 F77_CHAR_ARG_LEN (6) |
f9983d2761df
more xGELSD workspace fixes
Jaroslav Hajek <highegg@gmail.com>
parents:
7532
diff
changeset
|
2831 F77_CHAR_ARG_LEN (1)); |
7079 | 2832 |
2833 // We compute the size of rwork and iwork because ZGELSD in | |
2834 // older versions of LAPACK does not return them on a query | |
2835 // call. | |
7124 | 2836 double dminmn = static_cast<double> (minmn); |
2837 double dsmlsizp1 = static_cast<double> (smlsiz+1); | |
7079 | 2838 #if defined (HAVE_LOG2) |
7544
f9983d2761df
more xGELSD workspace fixes
Jaroslav Hajek <highegg@gmail.com>
parents:
7532
diff
changeset
|
2839 double tmp = log2 (dminmn / dsmlsizp1); |
7079 | 2840 #else |
7544
f9983d2761df
more xGELSD workspace fixes
Jaroslav Hajek <highegg@gmail.com>
parents:
7532
diff
changeset
|
2841 double tmp = log (dminmn / dsmlsizp1) / log (2.0); |
7079 | 2842 #endif |
7544
f9983d2761df
more xGELSD workspace fixes
Jaroslav Hajek <highegg@gmail.com>
parents:
7532
diff
changeset
|
2843 octave_idx_type nlvl = static_cast<octave_idx_type> (tmp) + 1; |
7079 | 2844 if (nlvl < 0) |
2845 nlvl = 0; | |
2846 | |
2847 octave_idx_type lrwork = minmn*(10 + 2*smlsiz + 8*nlvl) | |
2848 + 3*smlsiz*nrhs + (smlsiz+1)*(smlsiz+1); | |
2849 if (lrwork < 1) | |
2850 lrwork = 1; | |
2851 Array<double> rwork (lrwork); | |
2852 double *prwork = rwork.fortran_vec (); | |
2853 | |
2854 octave_idx_type liwork = 3 * minmn * nlvl + 11 * minmn; | |
2855 if (liwork < 1) | |
2856 liwork = 1; | |
2857 Array<octave_idx_type> iwork (liwork); | |
2858 octave_idx_type* piwork = iwork.fortran_vec (); | |
7072 | 2859 |
2860 F77_XFCN (zgelsd, ZGELSD, (m, n, nrhs, tmp_data, m, pretval, maxmn, | |
7788 | 2861 ps, rcon, rank, work.fortran_vec (), |
7079 | 2862 lwork, prwork, piwork, info)); |
1948 | 2863 |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2864 lwork = static_cast<octave_idx_type> (std::real (work(0))); |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2865 work.resize (lwork); |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2866 rwork.resize (static_cast<octave_idx_type> (rwork(0))); |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2867 iwork.resize (iwork(0)); |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2868 |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2869 F77_XFCN (zgelsd, ZGELSD, (m, n, nrhs, tmp_data, m, pretval, |
7788 | 2870 maxmn, ps, rcon, rank, |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2871 work.fortran_vec (), lwork, |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2872 prwork, piwork, info)); |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2873 |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2874 if (rank < minmn) |
1948 | 2875 { |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2876 if (rank < minmn) |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2877 (*current_liboctave_warning_handler) |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2878 ("zgelsd: rank deficient %dx%d matrix, rank = %d, tol = %e", |
7788 | 2879 m, n, rank, rcon); |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2880 |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2881 if (s.elem (0) == 0.0) |
7788 | 2882 rcon = 0.0; |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2883 else |
7788 | 2884 rcon = s.elem (minmn - 1) / s.elem (0); |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2885 |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
2886 retval.resize (n, nrhs); |
1948 | 2887 } |
458 | 2888 } |
2889 | |
2890 return retval; | |
2891 } | |
2892 | |
1205 | 2893 // column vector by row vector -> matrix operations |
2894 | |
2895 ComplexMatrix | |
2896 operator * (const ColumnVector& v, const ComplexRowVector& a) | |
2897 { | |
2898 ComplexColumnVector tmp (v); | |
2899 return tmp * a; | |
2900 } | |
2901 | |
2902 ComplexMatrix | |
2903 operator * (const ComplexColumnVector& a, const RowVector& b) | |
2904 { | |
2905 ComplexRowVector tmp (b); | |
2906 return a * tmp; | |
2907 } | |
2908 | |
2909 ComplexMatrix | |
2910 operator * (const ComplexColumnVector& v, const ComplexRowVector& a) | |
2911 { | |
1948 | 2912 ComplexMatrix retval; |
2913 | |
5275 | 2914 octave_idx_type len = v.length (); |
3233 | 2915 |
2916 if (len != 0) | |
1205 | 2917 { |
5275 | 2918 octave_idx_type a_len = a.length (); |
3233 | 2919 |
9359
be6867ba8104
avoid useless zero initialization when doing matrix multiply
Jaroslav Hajek <highegg@gmail.com>
parents:
9227
diff
changeset
|
2920 retval = ComplexMatrix (len, a_len); |
3233 | 2921 Complex *c = retval.fortran_vec (); |
2922 | |
4552 | 2923 F77_XFCN (zgemm, ZGEMM, (F77_CONST_CHAR_ARG2 ("N", 1), |
2924 F77_CONST_CHAR_ARG2 ("N", 1), | |
2925 len, a_len, 1, 1.0, v.data (), len, | |
2926 a.data (), 1, 0.0, c, len | |
2927 F77_CHAR_ARG_LEN (1) | |
2928 F77_CHAR_ARG_LEN (1))); | |
1205 | 2929 } |
2930 | |
1948 | 2931 return retval; |
1205 | 2932 } |
2933 | |
458 | 2934 // matrix by diagonal matrix -> matrix operations |
2935 | |
2936 ComplexMatrix& | |
2937 ComplexMatrix::operator += (const DiagMatrix& a) | |
2938 { | |
5275 | 2939 octave_idx_type nr = rows (); |
2940 octave_idx_type nc = cols (); | |
2941 | |
2942 octave_idx_type a_nr = rows (); | |
2943 octave_idx_type a_nc = cols (); | |
2384 | 2944 |
2945 if (nr != a_nr || nc != a_nc) | |
458 | 2946 { |
2384 | 2947 gripe_nonconformant ("operator +=", nr, nc, a_nr, a_nc); |
889 | 2948 return *this; |
458 | 2949 } |
2950 | |
5275 | 2951 for (octave_idx_type i = 0; i < a.length (); i++) |
458 | 2952 elem (i, i) += a.elem (i, i); |
2953 | |
2954 return *this; | |
2955 } | |
2956 | |
2957 ComplexMatrix& | |
2958 ComplexMatrix::operator -= (const DiagMatrix& a) | |
2959 { | |
5275 | 2960 octave_idx_type nr = rows (); |
2961 octave_idx_type nc = cols (); | |
2962 | |
2963 octave_idx_type a_nr = rows (); | |
2964 octave_idx_type a_nc = cols (); | |
2384 | 2965 |
2966 if (nr != a_nr || nc != a_nc) | |
458 | 2967 { |
2384 | 2968 gripe_nonconformant ("operator -=", nr, nc, a_nr, a_nc); |
889 | 2969 return *this; |
458 | 2970 } |
2971 | |
5275 | 2972 for (octave_idx_type i = 0; i < a.length (); i++) |
458 | 2973 elem (i, i) -= a.elem (i, i); |
2974 | |
2975 return *this; | |
2976 } | |
2977 | |
2978 ComplexMatrix& | |
2979 ComplexMatrix::operator += (const ComplexDiagMatrix& a) | |
2980 { | |
5275 | 2981 octave_idx_type nr = rows (); |
2982 octave_idx_type nc = cols (); | |
2983 | |
2984 octave_idx_type a_nr = rows (); | |
2985 octave_idx_type a_nc = cols (); | |
2384 | 2986 |
2987 if (nr != a_nr || nc != a_nc) | |
458 | 2988 { |
2384 | 2989 gripe_nonconformant ("operator +=", nr, nc, a_nr, a_nc); |
889 | 2990 return *this; |
458 | 2991 } |
2992 | |
5275 | 2993 for (octave_idx_type i = 0; i < a.length (); i++) |
458 | 2994 elem (i, i) += a.elem (i, i); |
2995 | |
2996 return *this; | |
2997 } | |
2998 | |
2999 ComplexMatrix& | |
3000 ComplexMatrix::operator -= (const ComplexDiagMatrix& a) | |
3001 { | |
5275 | 3002 octave_idx_type nr = rows (); |
3003 octave_idx_type nc = cols (); | |
3004 | |
3005 octave_idx_type a_nr = rows (); | |
3006 octave_idx_type a_nc = cols (); | |
2384 | 3007 |
3008 if (nr != a_nr || nc != a_nc) | |
458 | 3009 { |
2384 | 3010 gripe_nonconformant ("operator -=", nr, nc, a_nr, a_nc); |
889 | 3011 return *this; |
458 | 3012 } |
3013 | |
5275 | 3014 for (octave_idx_type i = 0; i < a.length (); i++) |
458 | 3015 elem (i, i) -= a.elem (i, i); |
3016 | |
3017 return *this; | |
3018 } | |
3019 | |
3020 // matrix by matrix -> matrix operations | |
3021 | |
3022 ComplexMatrix& | |
3023 ComplexMatrix::operator += (const Matrix& a) | |
3024 { | |
5275 | 3025 octave_idx_type nr = rows (); |
3026 octave_idx_type nc = cols (); | |
3027 | |
3028 octave_idx_type a_nr = a.rows (); | |
3029 octave_idx_type a_nc = a.cols (); | |
2384 | 3030 |
3031 if (nr != a_nr || nc != a_nc) | |
458 | 3032 { |
2384 | 3033 gripe_nonconformant ("operator +=", nr, nc, a_nr, a_nc); |
458 | 3034 return *this; |
3035 } | |
3036 | |
3037 if (nr == 0 || nc == 0) | |
3038 return *this; | |
3039 | |
3040 Complex *d = fortran_vec (); // Ensures only one reference to my privates! | |
3041 | |
9550
3d6a9aea2aea
refactor binary & bool ops in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9528
diff
changeset
|
3042 mx_inline_add2 (length (), d, a.data ()); |
458 | 3043 return *this; |
3044 } | |
3045 | |
3046 ComplexMatrix& | |
3047 ComplexMatrix::operator -= (const Matrix& a) | |
3048 { | |
5275 | 3049 octave_idx_type nr = rows (); |
3050 octave_idx_type nc = cols (); | |
3051 | |
3052 octave_idx_type a_nr = a.rows (); | |
3053 octave_idx_type a_nc = a.cols (); | |
2384 | 3054 |
3055 if (nr != a_nr || nc != a_nc) | |
458 | 3056 { |
2384 | 3057 gripe_nonconformant ("operator -=", nr, nc, a_nr, a_nc); |
458 | 3058 return *this; |
3059 } | |
3060 | |
3061 if (nr == 0 || nc == 0) | |
3062 return *this; | |
3063 | |
3064 Complex *d = fortran_vec (); // Ensures only one reference to my privates! | |
3065 | |
9550
3d6a9aea2aea
refactor binary & bool ops in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9528
diff
changeset
|
3066 mx_inline_sub2 (length (), d, a.data ()); |
458 | 3067 return *this; |
3068 } | |
3069 | |
3070 // unary operations | |
3071 | |
2964 | 3072 boolMatrix |
458 | 3073 ComplexMatrix::operator ! (void) const |
3074 { | |
9553
0c72d9284087
further bool ops tweaks
Jaroslav Hajek <highegg@gmail.com>
parents:
9551
diff
changeset
|
3075 return do_mx_unary_op<boolMatrix, ComplexMatrix> (*this, mx_inline_not); |
458 | 3076 } |
3077 | |
3078 // other operations | |
3079 | |
2384 | 3080 bool |
7922
935be827eaf8
error for NaN values in & and | expressions
John W. Eaton <jwe@octave.org>
parents:
7803
diff
changeset
|
3081 ComplexMatrix::any_element_is_nan (void) const |
935be827eaf8
error for NaN values in & and | expressions
John W. Eaton <jwe@octave.org>
parents:
7803
diff
changeset
|
3082 { |
935be827eaf8
error for NaN values in & and | expressions
John W. Eaton <jwe@octave.org>
parents:
7803
diff
changeset
|
3083 octave_idx_type nr = rows (); |
935be827eaf8
error for NaN values in & and | expressions
John W. Eaton <jwe@octave.org>
parents:
7803
diff
changeset
|
3084 octave_idx_type nc = cols (); |
935be827eaf8
error for NaN values in & and | expressions
John W. Eaton <jwe@octave.org>
parents:
7803
diff
changeset
|
3085 |
935be827eaf8
error for NaN values in & and | expressions
John W. Eaton <jwe@octave.org>
parents:
7803
diff
changeset
|
3086 for (octave_idx_type j = 0; j < nc; j++) |
935be827eaf8
error for NaN values in & and | expressions
John W. Eaton <jwe@octave.org>
parents:
7803
diff
changeset
|
3087 for (octave_idx_type i = 0; i < nr; i++) |
935be827eaf8
error for NaN values in & and | expressions
John W. Eaton <jwe@octave.org>
parents:
7803
diff
changeset
|
3088 { |
935be827eaf8
error for NaN values in & and | expressions
John W. Eaton <jwe@octave.org>
parents:
7803
diff
changeset
|
3089 Complex val = elem (i, j); |
935be827eaf8
error for NaN values in & and | expressions
John W. Eaton <jwe@octave.org>
parents:
7803
diff
changeset
|
3090 if (xisnan (val)) |
935be827eaf8
error for NaN values in & and | expressions
John W. Eaton <jwe@octave.org>
parents:
7803
diff
changeset
|
3091 return true; |
935be827eaf8
error for NaN values in & and | expressions
John W. Eaton <jwe@octave.org>
parents:
7803
diff
changeset
|
3092 } |
935be827eaf8
error for NaN values in & and | expressions
John W. Eaton <jwe@octave.org>
parents:
7803
diff
changeset
|
3093 |
935be827eaf8
error for NaN values in & and | expressions
John W. Eaton <jwe@octave.org>
parents:
7803
diff
changeset
|
3094 return false; |
935be827eaf8
error for NaN values in & and | expressions
John W. Eaton <jwe@octave.org>
parents:
7803
diff
changeset
|
3095 } |
935be827eaf8
error for NaN values in & and | expressions
John W. Eaton <jwe@octave.org>
parents:
7803
diff
changeset
|
3096 |
935be827eaf8
error for NaN values in & and | expressions
John W. Eaton <jwe@octave.org>
parents:
7803
diff
changeset
|
3097 bool |
2384 | 3098 ComplexMatrix::any_element_is_inf_or_nan (void) const |
3099 { | |
5275 | 3100 octave_idx_type nr = rows (); |
3101 octave_idx_type nc = cols (); | |
3102 | |
3103 for (octave_idx_type j = 0; j < nc; j++) | |
3104 for (octave_idx_type i = 0; i < nr; i++) | |
2384 | 3105 { |
3106 Complex val = elem (i, j); | |
3107 if (xisinf (val) || xisnan (val)) | |
3108 return true; | |
3109 } | |
3110 | |
3111 return false; | |
3112 } | |
3113 | |
2408 | 3114 // Return true if no elements have imaginary components. |
3115 | |
3116 bool | |
3117 ComplexMatrix::all_elements_are_real (void) const | |
3118 { | |
9825
7483fe200fab
narrow complex values with negative zero imaginary parts
Jaroslav Hajek <highegg@gmail.com>
parents:
9812
diff
changeset
|
3119 return mx_inline_all_real (numel (), data ()); |
2408 | 3120 } |
3121 | |
1968 | 3122 // Return nonzero if any element of CM has a non-integer real or |
3123 // imaginary part. Also extract the largest and smallest (real or | |
3124 // imaginary) values and return them in MAX_VAL and MIN_VAL. | |
3125 | |
2384 | 3126 bool |
1968 | 3127 ComplexMatrix::all_integers (double& max_val, double& min_val) const |
3128 { | |
5275 | 3129 octave_idx_type nr = rows (); |
3130 octave_idx_type nc = cols (); | |
1968 | 3131 |
3132 if (nr > 0 && nc > 0) | |
3133 { | |
3134 Complex val = elem (0, 0); | |
3135 | |
5315 | 3136 double r_val = std::real (val); |
3137 double i_val = std::imag (val); | |
1968 | 3138 |
3139 max_val = r_val; | |
3140 min_val = r_val; | |
3141 | |
3142 if (i_val > max_val) | |
3143 max_val = i_val; | |
3144 | |
3145 if (i_val < max_val) | |
3146 min_val = i_val; | |
3147 } | |
3148 else | |
2384 | 3149 return false; |
1968 | 3150 |
5275 | 3151 for (octave_idx_type j = 0; j < nc; j++) |
3152 for (octave_idx_type i = 0; i < nr; i++) | |
1968 | 3153 { |
3154 Complex val = elem (i, j); | |
3155 | |
5315 | 3156 double r_val = std::real (val); |
3157 double i_val = std::imag (val); | |
1968 | 3158 |
3159 if (r_val > max_val) | |
3160 max_val = r_val; | |
3161 | |
3162 if (i_val > max_val) | |
3163 max_val = i_val; | |
3164 | |
3165 if (r_val < min_val) | |
3166 min_val = r_val; | |
3167 | |
3168 if (i_val < min_val) | |
3169 min_val = i_val; | |
3170 | |
3171 if (D_NINT (r_val) != r_val || D_NINT (i_val) != i_val) | |
2384 | 3172 return false; |
1968 | 3173 } |
2384 | 3174 |
3175 return true; | |
1968 | 3176 } |
3177 | |
2384 | 3178 bool |
1968 | 3179 ComplexMatrix::too_large_for_float (void) const |
3180 { | |
5275 | 3181 octave_idx_type nr = rows (); |
3182 octave_idx_type nc = cols (); | |
3183 | |
3184 for (octave_idx_type j = 0; j < nc; j++) | |
3185 for (octave_idx_type i = 0; i < nr; i++) | |
1968 | 3186 { |
3187 Complex val = elem (i, j); | |
3188 | |
5315 | 3189 double r_val = std::real (val); |
3190 double i_val = std::imag (val); | |
1968 | 3191 |
5389 | 3192 if ((! (xisnan (r_val) || xisinf (r_val)) |
5387 | 3193 && fabs (r_val) > FLT_MAX) |
5389 | 3194 || (! (xisnan (i_val) || xisinf (i_val)) |
5387 | 3195 && fabs (i_val) > FLT_MAX)) |
2384 | 3196 return true; |
1968 | 3197 } |
3198 | |
2384 | 3199 return false; |
1968 | 3200 } |
3201 | |
5775 | 3202 // FIXME Do these really belong here? Maybe they should be |
4015 | 3203 // in a base class? |
3204 | |
2832 | 3205 boolMatrix |
4015 | 3206 ComplexMatrix::all (int dim) const |
458 | 3207 { |
9227
8145f2255276
use explicit template qualifs to please Intel C++ and MSVC++
Jaroslav Hajek <highegg@gmail.com>
parents:
8999
diff
changeset
|
3208 return do_mx_red_op<boolMatrix, Complex> (*this, dim, mx_inline_all); |
458 | 3209 } |
3210 | |
2832 | 3211 boolMatrix |
4015 | 3212 ComplexMatrix::any (int dim) const |
458 | 3213 { |
9227
8145f2255276
use explicit template qualifs to please Intel C++ and MSVC++
Jaroslav Hajek <highegg@gmail.com>
parents:
8999
diff
changeset
|
3214 return do_mx_red_op<boolMatrix, Complex> (*this, dim, mx_inline_any); |
458 | 3215 } |
3216 | |
3217 ComplexMatrix | |
3723 | 3218 ComplexMatrix::cumprod (int dim) const |
458 | 3219 { |
9227
8145f2255276
use explicit template qualifs to please Intel C++ and MSVC++
Jaroslav Hajek <highegg@gmail.com>
parents:
8999
diff
changeset
|
3220 return do_mx_cum_op<ComplexMatrix, Complex> (*this, dim, mx_inline_cumprod); |
458 | 3221 } |
3222 | |
3223 ComplexMatrix | |
3723 | 3224 ComplexMatrix::cumsum (int dim) const |
458 | 3225 { |
9227
8145f2255276
use explicit template qualifs to please Intel C++ and MSVC++
Jaroslav Hajek <highegg@gmail.com>
parents:
8999
diff
changeset
|
3226 return do_mx_cum_op<ComplexMatrix, Complex> (*this, dim, mx_inline_cumsum); |
458 | 3227 } |
3228 | |
3229 ComplexMatrix | |
3723 | 3230 ComplexMatrix::prod (int dim) const |
458 | 3231 { |
9227
8145f2255276
use explicit template qualifs to please Intel C++ and MSVC++
Jaroslav Hajek <highegg@gmail.com>
parents:
8999
diff
changeset
|
3232 return do_mx_red_op<ComplexMatrix, Complex> (*this, dim, mx_inline_prod); |
458 | 3233 } |
3234 | |
3235 ComplexMatrix | |
3723 | 3236 ComplexMatrix::sum (int dim) const |
458 | 3237 { |
9227
8145f2255276
use explicit template qualifs to please Intel C++ and MSVC++
Jaroslav Hajek <highegg@gmail.com>
parents:
8999
diff
changeset
|
3238 return do_mx_red_op<ComplexMatrix, Complex> (*this, dim, mx_inline_sum); |
458 | 3239 } |
3240 | |
3241 ComplexMatrix | |
3723 | 3242 ComplexMatrix::sumsq (int dim) const |
458 | 3243 { |
9227
8145f2255276
use explicit template qualifs to please Intel C++ and MSVC++
Jaroslav Hajek <highegg@gmail.com>
parents:
8999
diff
changeset
|
3244 return do_mx_red_op<Matrix, Complex> (*this, dim, mx_inline_sumsq); |
458 | 3245 } |
3246 | |
4329 | 3247 Matrix ComplexMatrix::abs (void) const |
3248 { | |
8650
a1ae2aae903e
abs,real,imag,conj: use code from mx-inlines rather than the generic map
Jaroslav Hajek <highegg@gmail.com>
parents:
8614
diff
changeset
|
3249 return Matrix (mx_inline_cabs_dup (data (), length ()), |
a1ae2aae903e
abs,real,imag,conj: use code from mx-inlines rather than the generic map
Jaroslav Hajek <highegg@gmail.com>
parents:
8614
diff
changeset
|
3250 rows (), cols ()); |
4329 | 3251 } |
3252 | |
7620
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7544
diff
changeset
|
3253 ComplexMatrix |
5275 | 3254 ComplexMatrix::diag (octave_idx_type k) const |
458 | 3255 { |
7620
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7544
diff
changeset
|
3256 return MArray2<Complex>::diag (k); |
458 | 3257 } |
3258 | |
2354 | 3259 bool |
5275 | 3260 ComplexMatrix::row_is_real_only (octave_idx_type i) const |
2354 | 3261 { |
3262 bool retval = true; | |
3263 | |
5275 | 3264 octave_idx_type nc = columns (); |
3265 | |
3266 for (octave_idx_type j = 0; j < nc; j++) | |
2354 | 3267 { |
5315 | 3268 if (std::imag (elem (i, j)) != 0.0) |
2354 | 3269 { |
3270 retval = false; | |
3271 break; | |
3272 } | |
3273 } | |
3274 | |
3275 return retval; | |
3276 } | |
3277 | |
3278 bool | |
5275 | 3279 ComplexMatrix::column_is_real_only (octave_idx_type j) const |
2354 | 3280 { |
3281 bool retval = true; | |
3282 | |
5275 | 3283 octave_idx_type nr = rows (); |
3284 | |
3285 for (octave_idx_type i = 0; i < nr; i++) | |
2354 | 3286 { |
5315 | 3287 if (std::imag (elem (i, j)) != 0.0) |
2354 | 3288 { |
3289 retval = false; | |
3290 break; | |
3291 } | |
3292 } | |
3293 | |
3294 return retval; | |
3295 } | |
891 | 3296 |
458 | 3297 ComplexColumnVector |
3298 ComplexMatrix::row_min (void) const | |
3299 { | |
5275 | 3300 Array<octave_idx_type> dummy_idx; |
4587 | 3301 return row_min (dummy_idx); |
458 | 3302 } |
3303 | |
3304 ComplexColumnVector | |
5275 | 3305 ComplexMatrix::row_min (Array<octave_idx_type>& idx_arg) const |
458 | 3306 { |
3307 ComplexColumnVector result; | |
3308 | |
5275 | 3309 octave_idx_type nr = rows (); |
3310 octave_idx_type nc = cols (); | |
458 | 3311 |
3312 if (nr > 0 && nc > 0) | |
3313 { | |
3314 result.resize (nr); | |
4587 | 3315 idx_arg.resize (nr); |
458 | 3316 |
5275 | 3317 for (octave_idx_type i = 0; i < nr; i++) |
458 | 3318 { |
2354 | 3319 bool real_only = row_is_real_only (i); |
3320 | |
5275 | 3321 octave_idx_type idx_j; |
4469 | 3322 |
3323 Complex tmp_min; | |
3324 | |
3325 double abs_min = octave_NaN; | |
3326 | |
3327 for (idx_j = 0; idx_j < nc; idx_j++) | |
3328 { | |
3329 tmp_min = elem (i, idx_j); | |
3330 | |
5389 | 3331 if (! xisnan (tmp_min)) |
4469 | 3332 { |
5315 | 3333 abs_min = real_only ? std::real (tmp_min) : std::abs (tmp_min); |
4469 | 3334 break; |
3335 } | |
3336 } | |
3337 | |
5275 | 3338 for (octave_idx_type j = idx_j+1; j < nc; j++) |
4469 | 3339 { |
3340 Complex tmp = elem (i, j); | |
3341 | |
5389 | 3342 if (xisnan (tmp)) |
4469 | 3343 continue; |
3344 | |
5315 | 3345 double abs_tmp = real_only ? std::real (tmp) : std::abs (tmp); |
4469 | 3346 |
3347 if (abs_tmp < abs_min) | |
3348 { | |
3349 idx_j = j; | |
3350 tmp_min = tmp; | |
3351 abs_min = abs_tmp; | |
3352 } | |
3353 } | |
3354 | |
5389 | 3355 if (xisnan (tmp_min)) |
4469 | 3356 { |
3357 result.elem (i) = Complex_NaN_result; | |
4587 | 3358 idx_arg.elem (i) = 0; |
4469 | 3359 } |
891 | 3360 else |
3361 { | |
4469 | 3362 result.elem (i) = tmp_min; |
4587 | 3363 idx_arg.elem (i) = idx_j; |
891 | 3364 } |
458 | 3365 } |
3366 } | |
3367 | |
3368 return result; | |
3369 } | |
3370 | |
3371 ComplexColumnVector | |
3372 ComplexMatrix::row_max (void) const | |
3373 { | |
5275 | 3374 Array<octave_idx_type> dummy_idx; |
4587 | 3375 return row_max (dummy_idx); |
458 | 3376 } |
3377 | |
3378 ComplexColumnVector | |
5275 | 3379 ComplexMatrix::row_max (Array<octave_idx_type>& idx_arg) const |
458 | 3380 { |
3381 ComplexColumnVector result; | |
3382 | |
5275 | 3383 octave_idx_type nr = rows (); |
3384 octave_idx_type nc = cols (); | |
458 | 3385 |
3386 if (nr > 0 && nc > 0) | |
3387 { | |
3388 result.resize (nr); | |
4587 | 3389 idx_arg.resize (nr); |
458 | 3390 |
5275 | 3391 for (octave_idx_type i = 0; i < nr; i++) |
458 | 3392 { |
2354 | 3393 bool real_only = row_is_real_only (i); |
3394 | |
5275 | 3395 octave_idx_type idx_j; |
4469 | 3396 |
3397 Complex tmp_max; | |
3398 | |
3399 double abs_max = octave_NaN; | |
3400 | |
3401 for (idx_j = 0; idx_j < nc; idx_j++) | |
3402 { | |
3403 tmp_max = elem (i, idx_j); | |
3404 | |
5389 | 3405 if (! xisnan (tmp_max)) |
4469 | 3406 { |
5315 | 3407 abs_max = real_only ? std::real (tmp_max) : std::abs (tmp_max); |
4469 | 3408 break; |
3409 } | |
3410 } | |
3411 | |
5275 | 3412 for (octave_idx_type j = idx_j+1; j < nc; j++) |
4469 | 3413 { |
3414 Complex tmp = elem (i, j); | |
3415 | |
5389 | 3416 if (xisnan (tmp)) |
4469 | 3417 continue; |
3418 | |
5315 | 3419 double abs_tmp = real_only ? std::real (tmp) : std::abs (tmp); |
4469 | 3420 |
3421 if (abs_tmp > abs_max) | |
3422 { | |
3423 idx_j = j; | |
3424 tmp_max = tmp; | |
3425 abs_max = abs_tmp; | |
3426 } | |
3427 } | |
3428 | |
5389 | 3429 if (xisnan (tmp_max)) |
4469 | 3430 { |
3431 result.elem (i) = Complex_NaN_result; | |
4587 | 3432 idx_arg.elem (i) = 0; |
4469 | 3433 } |
891 | 3434 else |
3435 { | |
4469 | 3436 result.elem (i) = tmp_max; |
4587 | 3437 idx_arg.elem (i) = idx_j; |
891 | 3438 } |
458 | 3439 } |
3440 } | |
3441 | |
3442 return result; | |
3443 } | |
3444 | |
3445 ComplexRowVector | |
3446 ComplexMatrix::column_min (void) const | |
3447 { | |
5275 | 3448 Array<octave_idx_type> dummy_idx; |
4587 | 3449 return column_min (dummy_idx); |
458 | 3450 } |
3451 | |
3452 ComplexRowVector | |
5275 | 3453 ComplexMatrix::column_min (Array<octave_idx_type>& idx_arg) const |
458 | 3454 { |
3455 ComplexRowVector result; | |
3456 | |
5275 | 3457 octave_idx_type nr = rows (); |
3458 octave_idx_type nc = cols (); | |
458 | 3459 |
3460 if (nr > 0 && nc > 0) | |
3461 { | |
3462 result.resize (nc); | |
4587 | 3463 idx_arg.resize (nc); |
458 | 3464 |
5275 | 3465 for (octave_idx_type j = 0; j < nc; j++) |
458 | 3466 { |
2354 | 3467 bool real_only = column_is_real_only (j); |
3468 | |
5275 | 3469 octave_idx_type idx_i; |
4469 | 3470 |
3471 Complex tmp_min; | |
3472 | |
3473 double abs_min = octave_NaN; | |
3474 | |
3475 for (idx_i = 0; idx_i < nr; idx_i++) | |
3476 { | |
3477 tmp_min = elem (idx_i, j); | |
3478 | |
5389 | 3479 if (! xisnan (tmp_min)) |
4469 | 3480 { |
5315 | 3481 abs_min = real_only ? std::real (tmp_min) : std::abs (tmp_min); |
4469 | 3482 break; |
3483 } | |
3484 } | |
3485 | |
5275 | 3486 for (octave_idx_type i = idx_i+1; i < nr; i++) |
4469 | 3487 { |
3488 Complex tmp = elem (i, j); | |
3489 | |
5389 | 3490 if (xisnan (tmp)) |
4469 | 3491 continue; |
3492 | |
5315 | 3493 double abs_tmp = real_only ? std::real (tmp) : std::abs (tmp); |
4469 | 3494 |
3495 if (abs_tmp < abs_min) | |
3496 { | |
3497 idx_i = i; | |
3498 tmp_min = tmp; | |
3499 abs_min = abs_tmp; | |
3500 } | |
3501 } | |
3502 | |
5389 | 3503 if (xisnan (tmp_min)) |
4469 | 3504 { |
3505 result.elem (j) = Complex_NaN_result; | |
4587 | 3506 idx_arg.elem (j) = 0; |
4469 | 3507 } |
891 | 3508 else |
3509 { | |
4469 | 3510 result.elem (j) = tmp_min; |
4587 | 3511 idx_arg.elem (j) = idx_i; |
891 | 3512 } |
458 | 3513 } |
3514 } | |
3515 | |
3516 return result; | |
3517 } | |
3518 | |
3519 ComplexRowVector | |
3520 ComplexMatrix::column_max (void) const | |
3521 { | |
5275 | 3522 Array<octave_idx_type> dummy_idx; |
4587 | 3523 return column_max (dummy_idx); |
458 | 3524 } |
3525 | |
3526 ComplexRowVector | |
5275 | 3527 ComplexMatrix::column_max (Array<octave_idx_type>& idx_arg) const |
458 | 3528 { |
3529 ComplexRowVector result; | |
3530 | |
5275 | 3531 octave_idx_type nr = rows (); |
3532 octave_idx_type nc = cols (); | |
458 | 3533 |
3534 if (nr > 0 && nc > 0) | |
3535 { | |
3536 result.resize (nc); | |
4587 | 3537 idx_arg.resize (nc); |
458 | 3538 |
5275 | 3539 for (octave_idx_type j = 0; j < nc; j++) |
458 | 3540 { |
2354 | 3541 bool real_only = column_is_real_only (j); |
3542 | |
5275 | 3543 octave_idx_type idx_i; |
4469 | 3544 |
3545 Complex tmp_max; | |
3546 | |
3547 double abs_max = octave_NaN; | |
3548 | |
3549 for (idx_i = 0; idx_i < nr; idx_i++) | |
3550 { | |
3551 tmp_max = elem (idx_i, j); | |
3552 | |
5389 | 3553 if (! xisnan (tmp_max)) |
4469 | 3554 { |
5315 | 3555 abs_max = real_only ? std::real (tmp_max) : std::abs (tmp_max); |
4469 | 3556 break; |
3557 } | |
3558 } | |
3559 | |
5275 | 3560 for (octave_idx_type i = idx_i+1; i < nr; i++) |
4469 | 3561 { |
3562 Complex tmp = elem (i, j); | |
3563 | |
5389 | 3564 if (xisnan (tmp)) |
4469 | 3565 continue; |
3566 | |
5315 | 3567 double abs_tmp = real_only ? std::real (tmp) : std::abs (tmp); |
4469 | 3568 |
3569 if (abs_tmp > abs_max) | |
3570 { | |
3571 idx_i = i; | |
3572 tmp_max = tmp; | |
3573 abs_max = abs_tmp; | |
3574 } | |
3575 } | |
3576 | |
5389 | 3577 if (xisnan (tmp_max)) |
4469 | 3578 { |
3579 result.elem (j) = Complex_NaN_result; | |
4587 | 3580 idx_arg.elem (j) = 0; |
4469 | 3581 } |
891 | 3582 else |
3583 { | |
4469 | 3584 result.elem (j) = tmp_max; |
4587 | 3585 idx_arg.elem (j) = idx_i; |
891 | 3586 } |
458 | 3587 } |
3588 } | |
3589 | |
3590 return result; | |
3591 } | |
3592 | |
3593 // i/o | |
3594 | |
3504 | 3595 std::ostream& |
3596 operator << (std::ostream& os, const ComplexMatrix& a) | |
458 | 3597 { |
5275 | 3598 for (octave_idx_type i = 0; i < a.rows (); i++) |
458 | 3599 { |
5275 | 3600 for (octave_idx_type j = 0; j < a.cols (); j++) |
4130 | 3601 { |
3602 os << " "; | |
3603 octave_write_complex (os, a.elem (i, j)); | |
3604 } | |
458 | 3605 os << "\n"; |
3606 } | |
3607 return os; | |
3608 } | |
3609 | |
3504 | 3610 std::istream& |
3611 operator >> (std::istream& is, ComplexMatrix& a) | |
458 | 3612 { |
5275 | 3613 octave_idx_type nr = a.rows (); |
3614 octave_idx_type nc = a.cols (); | |
458 | 3615 |
8999
dc07bc4157b8
allow empty matrices in stream input operators
Jaroslav Hajek <highegg@gmail.com>
parents:
8956
diff
changeset
|
3616 if (nr > 0 && nc > 0) |
458 | 3617 { |
3618 Complex tmp; | |
5275 | 3619 for (octave_idx_type i = 0; i < nr; i++) |
3620 for (octave_idx_type j = 0; j < nc; j++) | |
458 | 3621 { |
9469
c6edba80dfae
sanity checks for loading sparse matrices
John W. Eaton <jwe@octave.org>
parents:
9359
diff
changeset
|
3622 tmp = octave_read_value<Complex> (is); |
458 | 3623 if (is) |
3624 a.elem (i, j) = tmp; | |
3625 else | |
2993 | 3626 goto done; |
458 | 3627 } |
3628 } | |
3629 | |
2993 | 3630 done: |
3631 | |
458 | 3632 return is; |
3633 } | |
3634 | |
1819 | 3635 ComplexMatrix |
3636 Givens (const Complex& x, const Complex& y) | |
3637 { | |
3638 double cc; | |
3639 Complex cs, temp_r; | |
3640 | |
3887 | 3641 F77_FUNC (zlartg, ZLARTG) (x, y, cc, cs, temp_r); |
1819 | 3642 |
3643 ComplexMatrix g (2, 2); | |
3644 | |
3645 g.elem (0, 0) = cc; | |
3646 g.elem (1, 1) = cc; | |
3647 g.elem (0, 1) = cs; | |
3648 g.elem (1, 0) = -conj (cs); | |
3649 | |
3650 return g; | |
3651 } | |
3652 | |
3653 ComplexMatrix | |
3654 Sylvester (const ComplexMatrix& a, const ComplexMatrix& b, | |
3655 const ComplexMatrix& c) | |
3656 { | |
3657 ComplexMatrix retval; | |
3658 | |
5775 | 3659 // FIXME -- need to check that a, b, and c are all the same |
1819 | 3660 // size. |
3661 | |
3662 // Compute Schur decompositions | |
3663 | |
3664 ComplexSCHUR as (a, "U"); | |
3665 ComplexSCHUR bs (b, "U"); | |
3666 | |
3667 // Transform c to new coordinates. | |
3668 | |
3669 ComplexMatrix ua = as.unitary_matrix (); | |
3670 ComplexMatrix sch_a = as.schur_matrix (); | |
3671 | |
3672 ComplexMatrix ub = bs.unitary_matrix (); | |
3673 ComplexMatrix sch_b = bs.schur_matrix (); | |
3674 | |
3675 ComplexMatrix cx = ua.hermitian () * c * ub; | |
3676 | |
3677 // Solve the sylvester equation, back-transform, and return the | |
3678 // solution. | |
3679 | |
5275 | 3680 octave_idx_type a_nr = a.rows (); |
3681 octave_idx_type b_nr = b.rows (); | |
1819 | 3682 |
3683 double scale; | |
5275 | 3684 octave_idx_type info; |
1950 | 3685 |
3686 Complex *pa = sch_a.fortran_vec (); | |
3687 Complex *pb = sch_b.fortran_vec (); | |
3688 Complex *px = cx.fortran_vec (); | |
1819 | 3689 |
4552 | 3690 F77_XFCN (ztrsyl, ZTRSYL, (F77_CONST_CHAR_ARG2 ("N", 1), |
3691 F77_CONST_CHAR_ARG2 ("N", 1), | |
3692 1, a_nr, b_nr, pa, a_nr, pb, | |
3693 b_nr, px, a_nr, scale, info | |
3694 F77_CHAR_ARG_LEN (1) | |
3695 F77_CHAR_ARG_LEN (1))); | |
1950 | 3696 |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
3697 // FIXME -- check info? |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
3698 |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7478
diff
changeset
|
3699 retval = -ua * cx * ub.hermitian (); |
1819 | 3700 |
3701 return retval; | |
3702 } | |
3703 | |
2828 | 3704 ComplexMatrix |
3705 operator * (const ComplexMatrix& m, const Matrix& a) | |
3706 { | |
9663
7e5b4de5fbfe
improve mixed real x complex ops
Jaroslav Hajek <highegg@gmail.com>
parents:
9661
diff
changeset
|
3707 if (m.columns () > std::min (m.rows (), a.columns ()) / 10) |
7e5b4de5fbfe
improve mixed real x complex ops
Jaroslav Hajek <highegg@gmail.com>
parents:
9661
diff
changeset
|
3708 return ComplexMatrix (real (m) * a, imag (m) * a); |
7e5b4de5fbfe
improve mixed real x complex ops
Jaroslav Hajek <highegg@gmail.com>
parents:
9661
diff
changeset
|
3709 else |
7e5b4de5fbfe
improve mixed real x complex ops
Jaroslav Hajek <highegg@gmail.com>
parents:
9661
diff
changeset
|
3710 return m * ComplexMatrix (a); |
2828 | 3711 } |
3712 | |
3713 ComplexMatrix | |
3714 operator * (const Matrix& m, const ComplexMatrix& a) | |
3715 { | |
9663
7e5b4de5fbfe
improve mixed real x complex ops
Jaroslav Hajek <highegg@gmail.com>
parents:
9661
diff
changeset
|
3716 if (a.rows () > std::min (m.rows (), a.columns ()) / 10) |
7e5b4de5fbfe
improve mixed real x complex ops
Jaroslav Hajek <highegg@gmail.com>
parents:
9661
diff
changeset
|
3717 return ComplexMatrix (m * real (a), m * imag (a)); |
7e5b4de5fbfe
improve mixed real x complex ops
Jaroslav Hajek <highegg@gmail.com>
parents:
9661
diff
changeset
|
3718 else |
7e5b4de5fbfe
improve mixed real x complex ops
Jaroslav Hajek <highegg@gmail.com>
parents:
9661
diff
changeset
|
3719 return m * ComplexMatrix (a); |
2828 | 3720 } |
3721 | |
6162 | 3722 /* Simple Dot Product, Matrix-Vector and Matrix-Matrix Unit tests |
3723 %!assert([1+i 2+i 3+i] * [ 4+i ; 5+i ; 6+i], 29+21i, 1e-14) | |
3724 %!assert([1+i 2+i ; 3+i 4+i ] * [5+i ; 6+i], [15 + 14i ; 37 + 18i], 1e-14) | |
3725 %!assert([1+i 2+i ; 3+i 4+i ] * [5+i 6+i ; 7+i 8+i], [17 + 15i 20 + 17i; 41 + 19i 48 + 21i], 1e-14) | |
9526
f3ebc728ffd7
fix typos in complex xgemm
Jaroslav Hajek <highegg@gmail.com>
parents:
9523
diff
changeset
|
3726 %!assert([1 i]*[i 0]', -i); |
6162 | 3727 */ |
3728 | |
3729 /* Test some simple identities | |
3730 %!shared M, cv, rv | |
3731 %! M = randn(10,10)+i*rand(10,10); | |
3732 %! cv = randn(10,1)+i*rand(10,1); | |
3733 %! rv = randn(1,10)+i*rand(1,10); | |
3734 %!assert([M*cv,M*cv],M*[cv,cv],1e-14) | |
9528
ec066ba012c8
more fixes & tests for matrix multiply
Jaroslav Hajek <highegg@gmail.com>
parents:
9526
diff
changeset
|
3735 %!assert([M.'*cv,M.'*cv],M.'*[cv,cv],1e-14) |
ec066ba012c8
more fixes & tests for matrix multiply
Jaroslav Hajek <highegg@gmail.com>
parents:
9526
diff
changeset
|
3736 %!assert([M'*cv,M'*cv],M'*[cv,cv],1e-14) |
6162 | 3737 %!assert([rv*M;rv*M],[rv;rv]*M,1e-14) |
9528
ec066ba012c8
more fixes & tests for matrix multiply
Jaroslav Hajek <highegg@gmail.com>
parents:
9526
diff
changeset
|
3738 %!assert([rv*M.';rv*M.'],[rv;rv]*M.',1e-14) |
ec066ba012c8
more fixes & tests for matrix multiply
Jaroslav Hajek <highegg@gmail.com>
parents:
9526
diff
changeset
|
3739 %!assert([rv*M';rv*M'],[rv;rv]*M',1e-14) |
6162 | 3740 %!assert(2*rv*cv,[rv,rv]*[cv;cv],1e-14) |
3741 */ | |
3742 | |
7800
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3743 static const char * |
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3744 get_blas_trans_arg (bool trans, bool conj) |
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3745 { |
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3746 static char blas_notrans = 'N', blas_trans = 'T', blas_conj_trans = 'C'; |
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3747 return trans ? (conj ? &blas_conj_trans : &blas_trans) : &blas_notrans; |
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3748 } |
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3749 |
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3750 // the general GEMM operation |
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3751 |
2828 | 3752 ComplexMatrix |
9665
1dba57e9d08d
use blas_trans_type for xgemm
Jaroslav Hajek <highegg@gmail.com>
parents:
9663
diff
changeset
|
3753 xgemm (const ComplexMatrix& a, const ComplexMatrix& b, |
1dba57e9d08d
use blas_trans_type for xgemm
Jaroslav Hajek <highegg@gmail.com>
parents:
9663
diff
changeset
|
3754 blas_trans_type transa, blas_trans_type transb) |
2828 | 3755 { |
3756 ComplexMatrix retval; | |
3757 | |
9665
1dba57e9d08d
use blas_trans_type for xgemm
Jaroslav Hajek <highegg@gmail.com>
parents:
9663
diff
changeset
|
3758 bool tra = transa != blas_no_trans, trb = transb != blas_no_trans; |
1dba57e9d08d
use blas_trans_type for xgemm
Jaroslav Hajek <highegg@gmail.com>
parents:
9663
diff
changeset
|
3759 bool cja = transa == blas_conj_trans, cjb = transb == blas_conj_trans; |
1dba57e9d08d
use blas_trans_type for xgemm
Jaroslav Hajek <highegg@gmail.com>
parents:
9663
diff
changeset
|
3760 |
1dba57e9d08d
use blas_trans_type for xgemm
Jaroslav Hajek <highegg@gmail.com>
parents:
9663
diff
changeset
|
3761 octave_idx_type a_nr = tra ? a.cols () : a.rows (); |
1dba57e9d08d
use blas_trans_type for xgemm
Jaroslav Hajek <highegg@gmail.com>
parents:
9663
diff
changeset
|
3762 octave_idx_type a_nc = tra ? a.rows () : a.cols (); |
1dba57e9d08d
use blas_trans_type for xgemm
Jaroslav Hajek <highegg@gmail.com>
parents:
9663
diff
changeset
|
3763 |
1dba57e9d08d
use blas_trans_type for xgemm
Jaroslav Hajek <highegg@gmail.com>
parents:
9663
diff
changeset
|
3764 octave_idx_type b_nr = trb ? b.cols () : b.rows (); |
1dba57e9d08d
use blas_trans_type for xgemm
Jaroslav Hajek <highegg@gmail.com>
parents:
9663
diff
changeset
|
3765 octave_idx_type b_nc = trb ? b.rows () : b.cols (); |
7800
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3766 |
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3767 if (a_nc != b_nr) |
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3768 gripe_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc); |
2828 | 3769 else |
3770 { | |
7800
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3771 if (a_nr == 0 || a_nc == 0 || b_nc == 0) |
9359
be6867ba8104
avoid useless zero initialization when doing matrix multiply
Jaroslav Hajek <highegg@gmail.com>
parents:
9227
diff
changeset
|
3772 retval = ComplexMatrix (a_nr, b_nc, 0.0); |
9665
1dba57e9d08d
use blas_trans_type for xgemm
Jaroslav Hajek <highegg@gmail.com>
parents:
9663
diff
changeset
|
3773 else if (a.data () == b.data () && a_nr == b_nc && tra != trb) |
7801
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
3774 { |
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
3775 octave_idx_type lda = a.rows (); |
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
3776 |
9359
be6867ba8104
avoid useless zero initialization when doing matrix multiply
Jaroslav Hajek <highegg@gmail.com>
parents:
9227
diff
changeset
|
3777 retval = ComplexMatrix (a_nr, b_nc); |
7801
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
3778 Complex *c = retval.fortran_vec (); |
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
3779 |
9665
1dba57e9d08d
use blas_trans_type for xgemm
Jaroslav Hajek <highegg@gmail.com>
parents:
9663
diff
changeset
|
3780 const char *ctra = get_blas_trans_arg (tra, cja); |
1dba57e9d08d
use blas_trans_type for xgemm
Jaroslav Hajek <highegg@gmail.com>
parents:
9663
diff
changeset
|
3781 if (cja || cjb) |
7801
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
3782 { |
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
3783 F77_XFCN (zherk, ZHERK, (F77_CONST_CHAR_ARG2 ("U", 1), |
9665
1dba57e9d08d
use blas_trans_type for xgemm
Jaroslav Hajek <highegg@gmail.com>
parents:
9663
diff
changeset
|
3784 F77_CONST_CHAR_ARG2 (ctra, 1), |
7801
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
3785 a_nr, a_nc, 1.0, |
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
3786 a.data (), lda, 0.0, c, a_nr |
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
3787 F77_CHAR_ARG_LEN (1) |
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
3788 F77_CHAR_ARG_LEN (1))); |
10124
e5c6600e3400
use octave_idx_type where needed
Jaroslav Hajek <highegg@gmail.com>
parents:
10105
diff
changeset
|
3789 for (octave_idx_type j = 0; j < a_nr; j++) |
e5c6600e3400
use octave_idx_type where needed
Jaroslav Hajek <highegg@gmail.com>
parents:
10105
diff
changeset
|
3790 for (octave_idx_type i = 0; i < j; i++) |
7801
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
3791 retval.xelem (j,i) = std::conj (retval.xelem (i,j)); |
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
3792 } |
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
3793 else |
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
3794 { |
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
3795 F77_XFCN (zsyrk, ZSYRK, (F77_CONST_CHAR_ARG2 ("U", 1), |
9665
1dba57e9d08d
use blas_trans_type for xgemm
Jaroslav Hajek <highegg@gmail.com>
parents:
9663
diff
changeset
|
3796 F77_CONST_CHAR_ARG2 (ctra, 1), |
7801
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
3797 a_nr, a_nc, 1.0, |
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
3798 a.data (), lda, 0.0, c, a_nr |
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
3799 F77_CHAR_ARG_LEN (1) |
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
3800 F77_CHAR_ARG_LEN (1))); |
10124
e5c6600e3400
use octave_idx_type where needed
Jaroslav Hajek <highegg@gmail.com>
parents:
10105
diff
changeset
|
3801 for (octave_idx_type j = 0; j < a_nr; j++) |
e5c6600e3400
use octave_idx_type where needed
Jaroslav Hajek <highegg@gmail.com>
parents:
10105
diff
changeset
|
3802 for (octave_idx_type i = 0; i < j; i++) |
7801
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
3803 retval.xelem (j,i) = retval.xelem (i,j); |
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
3804 |
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
3805 } |
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
3806 |
776791438957
map symmetric cases to xHERK, xSYRK
Jaroslav Hajek <highegg@gmail.com>
parents:
7800
diff
changeset
|
3807 } |
2828 | 3808 else |
3809 { | |
7800
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3810 octave_idx_type lda = a.rows (), tda = a.cols (); |
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3811 octave_idx_type ldb = b.rows (), tdb = b.cols (); |
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3812 |
9359
be6867ba8104
avoid useless zero initialization when doing matrix multiply
Jaroslav Hajek <highegg@gmail.com>
parents:
9227
diff
changeset
|
3813 retval = ComplexMatrix (a_nr, b_nc); |
2828 | 3814 Complex *c = retval.fortran_vec (); |
3815 | |
7800
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3816 if (b_nc == 1 && a_nr == 1) |
5983 | 3817 { |
9665
1dba57e9d08d
use blas_trans_type for xgemm
Jaroslav Hajek <highegg@gmail.com>
parents:
9663
diff
changeset
|
3818 if (cja == cjb) |
7800
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3819 { |
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3820 F77_FUNC (xzdotu, XZDOTU) (a_nc, a.data (), 1, b.data (), 1, *c); |
9665
1dba57e9d08d
use blas_trans_type for xgemm
Jaroslav Hajek <highegg@gmail.com>
parents:
9663
diff
changeset
|
3821 if (cja) *c = std::conj (*c); |
7800
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3822 } |
9665
1dba57e9d08d
use blas_trans_type for xgemm
Jaroslav Hajek <highegg@gmail.com>
parents:
9663
diff
changeset
|
3823 else if (cja) |
7800
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3824 F77_FUNC (xzdotc, XZDOTC) (a_nc, a.data (), 1, b.data (), 1, *c); |
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3825 else |
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3826 F77_FUNC (xzdotc, XZDOTC) (a_nc, b.data (), 1, a.data (), 1, *c); |
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3827 } |
9665
1dba57e9d08d
use blas_trans_type for xgemm
Jaroslav Hajek <highegg@gmail.com>
parents:
9663
diff
changeset
|
3828 else if (b_nc == 1 && ! cjb) |
7800
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3829 { |
9665
1dba57e9d08d
use blas_trans_type for xgemm
Jaroslav Hajek <highegg@gmail.com>
parents:
9663
diff
changeset
|
3830 const char *ctra = get_blas_trans_arg (tra, cja); |
1dba57e9d08d
use blas_trans_type for xgemm
Jaroslav Hajek <highegg@gmail.com>
parents:
9663
diff
changeset
|
3831 F77_XFCN (zgemv, ZGEMV, (F77_CONST_CHAR_ARG2 (ctra, 1), |
7800
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3832 lda, tda, 1.0, a.data (), lda, |
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3833 b.data (), 1, 0.0, c, 1 |
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3834 F77_CHAR_ARG_LEN (1))); |
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3835 } |
9665
1dba57e9d08d
use blas_trans_type for xgemm
Jaroslav Hajek <highegg@gmail.com>
parents:
9663
diff
changeset
|
3836 else if (a_nr == 1 && ! cja && ! cjb) |
7800
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3837 { |
9665
1dba57e9d08d
use blas_trans_type for xgemm
Jaroslav Hajek <highegg@gmail.com>
parents:
9663
diff
changeset
|
3838 const char *crevtrb = get_blas_trans_arg (! trb, cjb); |
1dba57e9d08d
use blas_trans_type for xgemm
Jaroslav Hajek <highegg@gmail.com>
parents:
9663
diff
changeset
|
3839 F77_XFCN (zgemv, ZGEMV, (F77_CONST_CHAR_ARG2 (crevtrb, 1), |
7800
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3840 ldb, tdb, 1.0, b.data (), ldb, |
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3841 a.data (), 1, 0.0, c, 1 |
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3842 F77_CHAR_ARG_LEN (1))); |
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3843 } |
5983 | 3844 else |
6390 | 3845 { |
9665
1dba57e9d08d
use blas_trans_type for xgemm
Jaroslav Hajek <highegg@gmail.com>
parents:
9663
diff
changeset
|
3846 const char *ctra = get_blas_trans_arg (tra, cja); |
1dba57e9d08d
use blas_trans_type for xgemm
Jaroslav Hajek <highegg@gmail.com>
parents:
9663
diff
changeset
|
3847 const char *ctrb = get_blas_trans_arg (trb, cjb); |
1dba57e9d08d
use blas_trans_type for xgemm
Jaroslav Hajek <highegg@gmail.com>
parents:
9663
diff
changeset
|
3848 F77_XFCN (zgemm, ZGEMM, (F77_CONST_CHAR_ARG2 (ctra, 1), |
1dba57e9d08d
use blas_trans_type for xgemm
Jaroslav Hajek <highegg@gmail.com>
parents:
9663
diff
changeset
|
3849 F77_CONST_CHAR_ARG2 (ctrb, 1), |
7800
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3850 a_nr, b_nc, a_nc, 1.0, a.data (), |
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3851 lda, b.data (), ldb, 0.0, c, a_nr |
6390 | 3852 F77_CHAR_ARG_LEN (1) |
3853 F77_CHAR_ARG_LEN (1))); | |
3854 } | |
2828 | 3855 } |
3856 } | |
3857 | |
3858 return retval; | |
3859 } | |
3860 | |
7800
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3861 ComplexMatrix |
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3862 operator * (const ComplexMatrix& a, const ComplexMatrix& b) |
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3863 { |
9665
1dba57e9d08d
use blas_trans_type for xgemm
Jaroslav Hajek <highegg@gmail.com>
parents:
9663
diff
changeset
|
3864 return xgemm (a, b); |
7800
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3865 } |
5861b95e9879
support for compound operators, implement trans_mul, mul_trans, herm_mul and mul_herm
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
3866 |
5775 | 3867 // FIXME -- it would be nice to share code among the min/max |
4309 | 3868 // functions below. |
3869 | |
3870 #define EMPTY_RETURN_CHECK(T) \ | |
3871 if (nr == 0 || nc == 0) \ | |
3872 return T (nr, nc); | |
3873 | |
3874 ComplexMatrix | |
3875 min (const Complex& c, const ComplexMatrix& m) | |
3876 { | |
5275 | 3877 octave_idx_type nr = m.rows (); |
3878 octave_idx_type nc = m.columns (); | |
4309 | 3879 |
3880 EMPTY_RETURN_CHECK (ComplexMatrix); | |
3881 | |
3882 ComplexMatrix result (nr, nc); | |
3883 | |
5275 | 3884 for (octave_idx_type j = 0; j < nc; j++) |
3885 for (octave_idx_type i = 0; i < nr; i++) | |
4309 | 3886 { |
10142
829e69ec3110
make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents:
10124
diff
changeset
|
3887 octave_quit (); |
4309 | 3888 result (i, j) = xmin (c, m (i, j)); |
3889 } | |
3890 | |
3891 return result; | |
3892 } | |
3893 | |
3894 ComplexMatrix | |
3895 min (const ComplexMatrix& m, const Complex& c) | |
3896 { | |
5275 | 3897 octave_idx_type nr = m.rows (); |
3898 octave_idx_type nc = m.columns (); | |
4309 | 3899 |
3900 EMPTY_RETURN_CHECK (ComplexMatrix); | |
3901 | |
3902 ComplexMatrix result (nr, nc); | |
3903 | |
5275 | 3904 for (octave_idx_type j = 0; j < nc; j++) |
3905 for (octave_idx_type i = 0; i < nr; i++) | |
4309 | 3906 { |
10142
829e69ec3110
make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents:
10124
diff
changeset
|
3907 octave_quit (); |
4309 | 3908 result (i, j) = xmin (m (i, j), c); |
3909 } | |
3910 | |
3911 return result; | |
3912 } | |
3913 | |
3914 ComplexMatrix | |
3915 min (const ComplexMatrix& a, const ComplexMatrix& b) | |
3916 { | |
5275 | 3917 octave_idx_type nr = a.rows (); |
3918 octave_idx_type nc = a.columns (); | |
4309 | 3919 |
3920 if (nr != b.rows () || nc != b.columns ()) | |
3921 { | |
3922 (*current_liboctave_error_handler) | |
3923 ("two-arg min expecting args of same size"); | |
3924 return ComplexMatrix (); | |
3925 } | |
3926 | |
3927 EMPTY_RETURN_CHECK (ComplexMatrix); | |
3928 | |
3929 ComplexMatrix result (nr, nc); | |
3930 | |
5275 | 3931 for (octave_idx_type j = 0; j < nc; j++) |
4309 | 3932 { |
3933 int columns_are_real_only = 1; | |
5275 | 3934 for (octave_idx_type i = 0; i < nr; i++) |
4309 | 3935 { |
10142
829e69ec3110
make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents:
10124
diff
changeset
|
3936 octave_quit (); |
5315 | 3937 if (std::imag (a (i, j)) != 0.0 || std::imag (b (i, j)) != 0.0) |
4309 | 3938 { |
3939 columns_are_real_only = 0; | |
3940 break; | |
3941 } | |
3942 } | |
3943 | |
3944 if (columns_are_real_only) | |
3945 { | |
5275 | 3946 for (octave_idx_type i = 0; i < nr; i++) |
5315 | 3947 result (i, j) = xmin (std::real (a (i, j)), std::real (b (i, j))); |
4309 | 3948 } |
3949 else | |
3950 { | |
5275 | 3951 for (octave_idx_type i = 0; i < nr; i++) |
4309 | 3952 { |
10142
829e69ec3110
make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents:
10124
diff
changeset
|
3953 octave_quit (); |
4309 | 3954 result (i, j) = xmin (a (i, j), b (i, j)); |
3955 } | |
3956 } | |
3957 } | |
3958 | |
3959 return result; | |
3960 } | |
3961 | |
3962 ComplexMatrix | |
3963 max (const Complex& c, const ComplexMatrix& m) | |
3964 { | |
5275 | 3965 octave_idx_type nr = m.rows (); |
3966 octave_idx_type nc = m.columns (); | |
4309 | 3967 |
3968 EMPTY_RETURN_CHECK (ComplexMatrix); | |
3969 | |
3970 ComplexMatrix result (nr, nc); | |
3971 | |
5275 | 3972 for (octave_idx_type j = 0; j < nc; j++) |
3973 for (octave_idx_type i = 0; i < nr; i++) | |
4309 | 3974 { |
10142
829e69ec3110
make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents:
10124
diff
changeset
|
3975 octave_quit (); |
4309 | 3976 result (i, j) = xmax (c, m (i, j)); |
3977 } | |
3978 | |
3979 return result; | |
3980 } | |
3981 | |
3982 ComplexMatrix | |
3983 max (const ComplexMatrix& m, const Complex& c) | |
3984 { | |
5275 | 3985 octave_idx_type nr = m.rows (); |
3986 octave_idx_type nc = m.columns (); | |
4309 | 3987 |
3988 EMPTY_RETURN_CHECK (ComplexMatrix); | |
3989 | |
3990 ComplexMatrix result (nr, nc); | |
3991 | |
5275 | 3992 for (octave_idx_type j = 0; j < nc; j++) |
3993 for (octave_idx_type i = 0; i < nr; i++) | |
4309 | 3994 { |
10142
829e69ec3110
make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents:
10124
diff
changeset
|
3995 octave_quit (); |
4309 | 3996 result (i, j) = xmax (m (i, j), c); |
3997 } | |
3998 | |
3999 return result; | |
4000 } | |
4001 | |
4002 ComplexMatrix | |
4003 max (const ComplexMatrix& a, const ComplexMatrix& b) | |
4004 { | |
5275 | 4005 octave_idx_type nr = a.rows (); |
4006 octave_idx_type nc = a.columns (); | |
4309 | 4007 |
4008 if (nr != b.rows () || nc != b.columns ()) | |
4009 { | |
4010 (*current_liboctave_error_handler) | |
4011 ("two-arg max expecting args of same size"); | |
4012 return ComplexMatrix (); | |
4013 } | |
4014 | |
4015 EMPTY_RETURN_CHECK (ComplexMatrix); | |
4016 | |
4017 ComplexMatrix result (nr, nc); | |
4018 | |
5275 | 4019 for (octave_idx_type j = 0; j < nc; j++) |
4309 | 4020 { |
4021 int columns_are_real_only = 1; | |
5275 | 4022 for (octave_idx_type i = 0; i < nr; i++) |
4309 | 4023 { |
10142
829e69ec3110
make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents:
10124
diff
changeset
|
4024 octave_quit (); |
5315 | 4025 if (std::imag (a (i, j)) != 0.0 || std::imag (b (i, j)) != 0.0) |
4309 | 4026 { |
4027 columns_are_real_only = 0; | |
4028 break; | |
4029 } | |
4030 } | |
4031 | |
4032 if (columns_are_real_only) | |
4033 { | |
5275 | 4034 for (octave_idx_type i = 0; i < nr; i++) |
4309 | 4035 { |
10142
829e69ec3110
make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents:
10124
diff
changeset
|
4036 octave_quit (); |
5315 | 4037 result (i, j) = xmax (std::real (a (i, j)), std::real (b (i, j))); |
4309 | 4038 } |
4039 } | |
4040 else | |
4041 { | |
5275 | 4042 for (octave_idx_type i = 0; i < nr; i++) |
4309 | 4043 { |
10142
829e69ec3110
make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents:
10124
diff
changeset
|
4044 octave_quit (); |
4309 | 4045 result (i, j) = xmax (a (i, j), b (i, j)); |
4046 } | |
4047 } | |
4048 } | |
4049 | |
4050 return result; | |
4051 } | |
4052 | |
9653
e087d7c77ff9
improve linspace in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9601
diff
changeset
|
4053 ComplexMatrix linspace (const ComplexColumnVector& x1, |
e087d7c77ff9
improve linspace in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9601
diff
changeset
|
4054 const ComplexColumnVector& x2, |
e087d7c77ff9
improve linspace in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9601
diff
changeset
|
4055 octave_idx_type n) |
e087d7c77ff9
improve linspace in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9601
diff
changeset
|
4056 |
e087d7c77ff9
improve linspace in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9601
diff
changeset
|
4057 { |
e087d7c77ff9
improve linspace in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9601
diff
changeset
|
4058 if (n < 1) n = 1; |
e087d7c77ff9
improve linspace in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9601
diff
changeset
|
4059 |
e087d7c77ff9
improve linspace in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9601
diff
changeset
|
4060 octave_idx_type m = x1.length (); |
e087d7c77ff9
improve linspace in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9601
diff
changeset
|
4061 |
e087d7c77ff9
improve linspace in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9601
diff
changeset
|
4062 if (x2.length () != m) |
e087d7c77ff9
improve linspace in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9601
diff
changeset
|
4063 (*current_liboctave_error_handler) ("linspace: vectors must be of equal length"); |
e087d7c77ff9
improve linspace in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9601
diff
changeset
|
4064 |
e087d7c77ff9
improve linspace in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9601
diff
changeset
|
4065 NoAlias<ComplexMatrix> retval; |
e087d7c77ff9
improve linspace in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9601
diff
changeset
|
4066 |
e087d7c77ff9
improve linspace in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9601
diff
changeset
|
4067 retval.clear (m, n); |
e087d7c77ff9
improve linspace in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9601
diff
changeset
|
4068 for (octave_idx_type i = 0; i < m; i++) |
e087d7c77ff9
improve linspace in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9601
diff
changeset
|
4069 retval(i, 0) = x1(i); |
e087d7c77ff9
improve linspace in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9601
diff
changeset
|
4070 |
e087d7c77ff9
improve linspace in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9601
diff
changeset
|
4071 // The last column is not needed while using delta. |
9658
3429c956de6f
extend linspace & fix up liboctave rewrite
Jaroslav Hajek <highegg@gmail.com>
parents:
9653
diff
changeset
|
4072 Complex *delta = &retval(0, n-1); |
9653
e087d7c77ff9
improve linspace in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9601
diff
changeset
|
4073 for (octave_idx_type i = 0; i < m; i++) |
e087d7c77ff9
improve linspace in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9601
diff
changeset
|
4074 delta[i] = (x2(i) - x1(i)) / (n - 1.0); |
e087d7c77ff9
improve linspace in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9601
diff
changeset
|
4075 |
e087d7c77ff9
improve linspace in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9601
diff
changeset
|
4076 for (octave_idx_type j = 1; j < n-1; j++) |
e087d7c77ff9
improve linspace in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9601
diff
changeset
|
4077 for (octave_idx_type i = 0; i < m; i++) |
9658
3429c956de6f
extend linspace & fix up liboctave rewrite
Jaroslav Hajek <highegg@gmail.com>
parents:
9653
diff
changeset
|
4078 retval(i, j) = x1(i) + static_cast<double> (j)*delta[i]; |
9653
e087d7c77ff9
improve linspace in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9601
diff
changeset
|
4079 |
e087d7c77ff9
improve linspace in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9601
diff
changeset
|
4080 for (octave_idx_type i = 0; i < m; i++) |
e087d7c77ff9
improve linspace in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9601
diff
changeset
|
4081 retval(i, n-1) = x2(i); |
e087d7c77ff9
improve linspace in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9601
diff
changeset
|
4082 |
e087d7c77ff9
improve linspace in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9601
diff
changeset
|
4083 return retval; |
e087d7c77ff9
improve linspace in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9601
diff
changeset
|
4084 } |
e087d7c77ff9
improve linspace in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9601
diff
changeset
|
4085 |
9578
7dafdb8b062f
refactor comparison ops implementations
Jaroslav Hajek <highegg@gmail.com>
parents:
9553
diff
changeset
|
4086 MS_CMP_OPS (ComplexMatrix, Complex) |
9550
3d6a9aea2aea
refactor binary & bool ops in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9528
diff
changeset
|
4087 MS_BOOL_OPS (ComplexMatrix, Complex) |
2870 | 4088 |
9578
7dafdb8b062f
refactor comparison ops implementations
Jaroslav Hajek <highegg@gmail.com>
parents:
9553
diff
changeset
|
4089 SM_CMP_OPS (Complex, ComplexMatrix) |
9550
3d6a9aea2aea
refactor binary & bool ops in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9528
diff
changeset
|
4090 SM_BOOL_OPS (Complex, ComplexMatrix) |
2870 | 4091 |
9578
7dafdb8b062f
refactor comparison ops implementations
Jaroslav Hajek <highegg@gmail.com>
parents:
9553
diff
changeset
|
4092 MM_CMP_OPS (ComplexMatrix, ComplexMatrix) |
9550
3d6a9aea2aea
refactor binary & bool ops in liboctave
Jaroslav Hajek <highegg@gmail.com>
parents:
9528
diff
changeset
|
4093 MM_BOOL_OPS (ComplexMatrix, ComplexMatrix) |