File: | liboctave/array/dMatrix.cc |
Location: | line 1542, column 3 |
Description: | Undefined or garbage value returned to caller |
1 | // Matrix manipulations. | |||
2 | /* | |||
3 | ||||
4 | Copyright (C) 1994-2013 John W. Eaton | |||
5 | Copyright (C) 2008 Jaroslav Hajek | |||
6 | Copyright (C) 2009 VZLU Prague, a.s. | |||
7 | ||||
8 | This file is part of Octave. | |||
9 | ||||
10 | Octave is free software; you can redistribute it and/or modify it | |||
11 | under the terms of the GNU General Public License as published by the | |||
12 | Free Software Foundation; either version 3 of the License, or (at your | |||
13 | option) any later version. | |||
14 | ||||
15 | Octave is distributed in the hope that it will be useful, but WITHOUT | |||
16 | ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |||
17 | FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |||
18 | for more details. | |||
19 | ||||
20 | You should have received a copy of the GNU General Public License | |||
21 | along with Octave; see the file COPYING. If not, see | |||
22 | <http://www.gnu.org/licenses/>. | |||
23 | ||||
24 | */ | |||
25 | ||||
26 | #ifdef HAVE_CONFIG_H1 | |||
27 | #include <config.h> | |||
28 | #endif | |||
29 | ||||
30 | #include <cfloat> | |||
31 | ||||
32 | #include <iostream> | |||
33 | #include <vector> | |||
34 | ||||
35 | #include "Array-util.h" | |||
36 | #include "byte-swap.h" | |||
37 | #include "dMatrix.h" | |||
38 | #include "dbleAEPBAL.h" | |||
39 | #include "DET.h" | |||
40 | #include "dbleSCHUR.h" | |||
41 | #include "dbleSVD.h" | |||
42 | #include "dbleCHOL.h" | |||
43 | #include "f77-fcn.h" | |||
44 | #include "functor.h" | |||
45 | #include "lo-error.h" | |||
46 | #include "oct-locbuf.h" | |||
47 | #include "lo-ieee.h" | |||
48 | #include "lo-mappers.h" | |||
49 | #include "lo-utils.h" | |||
50 | #include "mx-base.h" | |||
51 | #include "mx-m-dm.h" | |||
52 | #include "mx-dm-m.h" | |||
53 | #include "mx-inlines.cc" | |||
54 | #include "mx-op-defs.h" | |||
55 | #include "oct-cmplx.h" | |||
56 | #include "oct-fftw.h" | |||
57 | #include "oct-norm.h" | |||
58 | #include "quit.h" | |||
59 | ||||
60 | // Fortran functions we call. | |||
61 | ||||
62 | extern "C" | |||
63 | { | |||
64 | F77_RET_Tint | |||
65 | F77_FUNC (xilaenv, XILAENV)xilaenv_ (const octave_idx_type&, | |||
66 | F77_CONST_CHAR_ARG_DECLconst char *, | |||
67 | F77_CONST_CHAR_ARG_DECLconst char *, | |||
68 | const octave_idx_type&, const octave_idx_type&, | |||
69 | const octave_idx_type&, const octave_idx_type&, | |||
70 | octave_idx_type& | |||
71 | F77_CHAR_ARG_LEN_DECL, long | |||
72 | F77_CHAR_ARG_LEN_DECL, long); | |||
73 | ||||
74 | F77_RET_Tint | |||
75 | F77_FUNC (dgebal, DGEBAL)dgebal_ (F77_CONST_CHAR_ARG_DECLconst char *, | |||
76 | const octave_idx_type&, double*, | |||
77 | const octave_idx_type&, octave_idx_type&, | |||
78 | octave_idx_type&, double*, octave_idx_type& | |||
79 | F77_CHAR_ARG_LEN_DECL, long); | |||
80 | ||||
81 | F77_RET_Tint | |||
82 | F77_FUNC (dgebak, DGEBAK)dgebak_ (F77_CONST_CHAR_ARG_DECLconst char *, | |||
83 | F77_CONST_CHAR_ARG_DECLconst char *, | |||
84 | const octave_idx_type&, const octave_idx_type&, | |||
85 | const octave_idx_type&, double*, | |||
86 | const octave_idx_type&, double*, | |||
87 | const octave_idx_type&, octave_idx_type& | |||
88 | F77_CHAR_ARG_LEN_DECL, long | |||
89 | F77_CHAR_ARG_LEN_DECL, long); | |||
90 | ||||
91 | ||||
92 | F77_RET_Tint | |||
93 | F77_FUNC (dgemm, DGEMM)dgemm_ (F77_CONST_CHAR_ARG_DECLconst char *, | |||
94 | F77_CONST_CHAR_ARG_DECLconst char *, | |||
95 | const octave_idx_type&, const octave_idx_type&, | |||
96 | const octave_idx_type&, const double&, | |||
97 | const double*, const octave_idx_type&, | |||
98 | const double*, const octave_idx_type&, | |||
99 | const double&, double*, const octave_idx_type& | |||
100 | F77_CHAR_ARG_LEN_DECL, long | |||
101 | F77_CHAR_ARG_LEN_DECL, long); | |||
102 | ||||
103 | F77_RET_Tint | |||
104 | F77_FUNC (dgemv, DGEMV)dgemv_ (F77_CONST_CHAR_ARG_DECLconst char *, | |||
105 | const octave_idx_type&, const octave_idx_type&, | |||
106 | const double&, const double*, | |||
107 | const octave_idx_type&, const double*, | |||
108 | const octave_idx_type&, const double&, double*, | |||
109 | const octave_idx_type& | |||
110 | F77_CHAR_ARG_LEN_DECL, long); | |||
111 | ||||
112 | F77_RET_Tint | |||
113 | F77_FUNC (xddot, XDDOT)xddot_ (const octave_idx_type&, const double*, | |||
114 | const octave_idx_type&, const double*, | |||
115 | const octave_idx_type&, double&); | |||
116 | ||||
117 | F77_RET_Tint | |||
118 | F77_FUNC (dsyrk, DSYRK)dsyrk_ (F77_CONST_CHAR_ARG_DECLconst char *, | |||
119 | F77_CONST_CHAR_ARG_DECLconst char *, | |||
120 | const octave_idx_type&, const octave_idx_type&, | |||
121 | const double&, const double*, const octave_idx_type&, | |||
122 | const double&, double*, const octave_idx_type& | |||
123 | F77_CHAR_ARG_LEN_DECL, long | |||
124 | F77_CHAR_ARG_LEN_DECL, long); | |||
125 | ||||
126 | F77_RET_Tint | |||
127 | F77_FUNC (dgetrf, DGETRF)dgetrf_ (const octave_idx_type&, const octave_idx_type&, | |||
128 | double*, const octave_idx_type&, | |||
129 | octave_idx_type*, octave_idx_type&); | |||
130 | ||||
131 | F77_RET_Tint | |||
132 | F77_FUNC (dgetrs, DGETRS)dgetrs_ (F77_CONST_CHAR_ARG_DECLconst char *, | |||
133 | const octave_idx_type&, const octave_idx_type&, | |||
134 | const double*, const octave_idx_type&, | |||
135 | const octave_idx_type*, double*, | |||
136 | const octave_idx_type&, octave_idx_type& | |||
137 | F77_CHAR_ARG_LEN_DECL, long); | |||
138 | ||||
139 | F77_RET_Tint | |||
140 | F77_FUNC (dgetri, DGETRI)dgetri_ (const octave_idx_type&, double*, | |||
141 | const octave_idx_type&, const octave_idx_type*, | |||
142 | double*, const octave_idx_type&, | |||
143 | octave_idx_type&); | |||
144 | ||||
145 | F77_RET_Tint | |||
146 | F77_FUNC (dgecon, DGECON)dgecon_ (F77_CONST_CHAR_ARG_DECLconst char *, | |||
147 | const octave_idx_type&, double*, | |||
148 | const octave_idx_type&, const double&, double&, | |||
149 | double*, octave_idx_type*, octave_idx_type& | |||
150 | F77_CHAR_ARG_LEN_DECL, long); | |||
151 | ||||
152 | F77_RET_Tint | |||
153 | F77_FUNC (dgelsy, DGELSY)dgelsy_ (const octave_idx_type&, const octave_idx_type&, | |||
154 | const octave_idx_type&, double*, | |||
155 | const octave_idx_type&, double*, | |||
156 | const octave_idx_type&, octave_idx_type*, | |||
157 | double&, octave_idx_type&, double*, | |||
158 | const octave_idx_type&, octave_idx_type&); | |||
159 | ||||
160 | F77_RET_Tint | |||
161 | F77_FUNC (dgelsd, DGELSD)dgelsd_ (const octave_idx_type&, const octave_idx_type&, | |||
162 | const octave_idx_type&, double*, | |||
163 | const octave_idx_type&, double*, | |||
164 | const octave_idx_type&, double*, double&, | |||
165 | octave_idx_type&, double*, | |||
166 | const octave_idx_type&, octave_idx_type*, | |||
167 | octave_idx_type&); | |||
168 | ||||
169 | F77_RET_Tint | |||
170 | F77_FUNC (dpotrf, DPOTRF)dpotrf_ (F77_CONST_CHAR_ARG_DECLconst char *, | |||
171 | const octave_idx_type&, double *, | |||
172 | const octave_idx_type&, octave_idx_type& | |||
173 | F77_CHAR_ARG_LEN_DECL, long); | |||
174 | ||||
175 | F77_RET_Tint | |||
176 | F77_FUNC (dpocon, DPOCON)dpocon_ (F77_CONST_CHAR_ARG_DECLconst char *, | |||
177 | const octave_idx_type&, double*, | |||
178 | const octave_idx_type&, const double&, | |||
179 | double&, double*, octave_idx_type*, | |||
180 | octave_idx_type& | |||
181 | F77_CHAR_ARG_LEN_DECL, long); | |||
182 | F77_RET_Tint | |||
183 | F77_FUNC (dpotrs, DPOTRS)dpotrs_ (F77_CONST_CHAR_ARG_DECLconst char *, | |||
184 | const octave_idx_type&, const octave_idx_type&, | |||
185 | const double*, const octave_idx_type&, double*, | |||
186 | const octave_idx_type&, octave_idx_type& | |||
187 | F77_CHAR_ARG_LEN_DECL, long); | |||
188 | ||||
189 | F77_RET_Tint | |||
190 | F77_FUNC (dtrtri, DTRTRI)dtrtri_ (F77_CONST_CHAR_ARG_DECLconst char *, | |||
191 | F77_CONST_CHAR_ARG_DECLconst char *, | |||
192 | const octave_idx_type&, const double*, | |||
193 | const octave_idx_type&, octave_idx_type& | |||
194 | F77_CHAR_ARG_LEN_DECL, long | |||
195 | F77_CHAR_ARG_LEN_DECL, long); | |||
196 | F77_RET_Tint | |||
197 | F77_FUNC (dtrcon, DTRCON)dtrcon_ (F77_CONST_CHAR_ARG_DECLconst char *, | |||
198 | F77_CONST_CHAR_ARG_DECLconst char *, | |||
199 | F77_CONST_CHAR_ARG_DECLconst char *, | |||
200 | const octave_idx_type&, const double*, | |||
201 | const octave_idx_type&, double&, | |||
202 | double*, octave_idx_type*, octave_idx_type& | |||
203 | F77_CHAR_ARG_LEN_DECL, long | |||
204 | F77_CHAR_ARG_LEN_DECL, long | |||
205 | F77_CHAR_ARG_LEN_DECL, long); | |||
206 | F77_RET_Tint | |||
207 | F77_FUNC (dtrtrs, DTRTRS)dtrtrs_ (F77_CONST_CHAR_ARG_DECLconst char *, | |||
208 | F77_CONST_CHAR_ARG_DECLconst char *, | |||
209 | F77_CONST_CHAR_ARG_DECLconst char *, | |||
210 | const octave_idx_type&, const octave_idx_type&, | |||
211 | const double*, const octave_idx_type&, double*, | |||
212 | const octave_idx_type&, octave_idx_type& | |||
213 | F77_CHAR_ARG_LEN_DECL, long | |||
214 | F77_CHAR_ARG_LEN_DECL, long | |||
215 | F77_CHAR_ARG_LEN_DECL, long); | |||
216 | ||||
217 | F77_RET_Tint | |||
218 | F77_FUNC (dlartg, DLARTG)dlartg_ (const double&, const double&, double&, | |||
219 | double&, double&); | |||
220 | ||||
221 | F77_RET_Tint | |||
222 | F77_FUNC (dtrsyl, DTRSYL)dtrsyl_ (F77_CONST_CHAR_ARG_DECLconst char *, | |||
223 | F77_CONST_CHAR_ARG_DECLconst char *, | |||
224 | const octave_idx_type&, const octave_idx_type&, | |||
225 | const octave_idx_type&, const double*, | |||
226 | const octave_idx_type&, const double*, | |||
227 | const octave_idx_type&, const double*, | |||
228 | const octave_idx_type&, double&, octave_idx_type& | |||
229 | F77_CHAR_ARG_LEN_DECL, long | |||
230 | F77_CHAR_ARG_LEN_DECL, long); | |||
231 | ||||
232 | F77_RET_Tint | |||
233 | F77_FUNC (xdlange, XDLANGE)xdlange_ (F77_CONST_CHAR_ARG_DECLconst char *, | |||
234 | const octave_idx_type&, const octave_idx_type&, | |||
235 | const double*, const octave_idx_type&, | |||
236 | double*, double& | |||
237 | F77_CHAR_ARG_LEN_DECL, long); | |||
238 | } | |||
239 | ||||
240 | // Matrix class. | |||
241 | ||||
242 | Matrix::Matrix (const RowVector& rv) | |||
243 | : MArray<double> (rv) | |||
244 | { | |||
245 | } | |||
246 | ||||
247 | Matrix::Matrix (const ColumnVector& cv) | |||
248 | : MArray<double> (cv) | |||
249 | { | |||
250 | } | |||
251 | ||||
252 | Matrix::Matrix (const DiagMatrix& a) | |||
253 | : MArray<double> (a.dims (), 0.0) | |||
254 | { | |||
255 | for (octave_idx_type i = 0; i < a.length (); i++) | |||
256 | elem (i, i) = a.elem (i, i); | |||
257 | } | |||
258 | ||||
259 | Matrix::Matrix (const PermMatrix& a) | |||
260 | : MArray<double> (a.dims (), 0.0) | |||
261 | { | |||
262 | const Array<octave_idx_type> ia (a.pvec ()); | |||
263 | octave_idx_type len = a.rows (); | |||
264 | if (a.is_col_perm ()) | |||
265 | for (octave_idx_type i = 0; i < len; i++) | |||
266 | elem (ia(i), i) = 1.0; | |||
267 | else | |||
268 | for (octave_idx_type i = 0; i < len; i++) | |||
269 | elem (i, ia(i)) = 1.0; | |||
270 | } | |||
271 | ||||
272 | // FIXME: could we use a templated mixed-type copy function here? | |||
273 | ||||
274 | Matrix::Matrix (const boolMatrix& a) | |||
275 | : MArray<double> (a) | |||
276 | { | |||
277 | } | |||
278 | ||||
279 | Matrix::Matrix (const charMatrix& a) | |||
280 | : MArray<double> (a.dims ()) | |||
281 | { | |||
282 | for (octave_idx_type i = 0; i < a.rows (); i++) | |||
283 | for (octave_idx_type j = 0; j < a.cols (); j++) | |||
284 | elem (i, j) = static_cast<unsigned char> (a.elem (i, j)); | |||
285 | } | |||
286 | ||||
287 | bool | |||
288 | Matrix::operator == (const Matrix& a) const | |||
289 | { | |||
290 | if (rows () != a.rows () || cols () != a.cols ()) | |||
291 | return false; | |||
292 | ||||
293 | return mx_inline_equal (length (), data (), a.data ()); | |||
294 | } | |||
295 | ||||
296 | bool | |||
297 | Matrix::operator != (const Matrix& a) const | |||
298 | { | |||
299 | return !(*this == a); | |||
300 | } | |||
301 | ||||
302 | bool | |||
303 | Matrix::is_symmetric (void) const | |||
304 | { | |||
305 | if (is_square () && rows () > 0) | |||
306 | { | |||
307 | for (octave_idx_type i = 0; i < rows (); i++) | |||
308 | for (octave_idx_type j = i+1; j < cols (); j++) | |||
309 | if (elem (i, j) != elem (j, i)) | |||
310 | return false; | |||
311 | ||||
312 | return true; | |||
313 | } | |||
314 | ||||
315 | return false; | |||
316 | } | |||
317 | ||||
318 | Matrix& | |||
319 | Matrix::insert (const Matrix& a, octave_idx_type r, octave_idx_type c) | |||
320 | { | |||
321 | Array<double>::insert (a, r, c); | |||
322 | return *this; | |||
323 | } | |||
324 | ||||
325 | Matrix& | |||
326 | Matrix::insert (const RowVector& a, octave_idx_type r, octave_idx_type c) | |||
327 | { | |||
328 | octave_idx_type a_len = a.length (); | |||
329 | ||||
330 | if (r < 0 || r >= rows () || c < 0 || c + a_len > cols ()) | |||
331 | { | |||
332 | (*current_liboctave_error_handler) ("range error for insert"); | |||
333 | return *this; | |||
334 | } | |||
335 | ||||
336 | if (a_len > 0) | |||
337 | { | |||
338 | make_unique (); | |||
339 | ||||
340 | for (octave_idx_type i = 0; i < a_len; i++) | |||
341 | xelem (r, c+i) = a.elem (i); | |||
342 | } | |||
343 | ||||
344 | return *this; | |||
345 | } | |||
346 | ||||
347 | Matrix& | |||
348 | Matrix::insert (const ColumnVector& a, octave_idx_type r, octave_idx_type c) | |||
349 | { | |||
350 | octave_idx_type a_len = a.length (); | |||
351 | ||||
352 | if (r < 0 || r + a_len > rows () || c < 0 || c >= cols ()) | |||
353 | { | |||
354 | (*current_liboctave_error_handler) ("range error for insert"); | |||
355 | return *this; | |||
356 | } | |||
357 | ||||
358 | if (a_len > 0) | |||
359 | { | |||
360 | make_unique (); | |||
361 | ||||
362 | for (octave_idx_type i = 0; i < a_len; i++) | |||
363 | xelem (r+i, c) = a.elem (i); | |||
364 | } | |||
365 | ||||
366 | return *this; | |||
367 | } | |||
368 | ||||
369 | Matrix& | |||
370 | Matrix::insert (const DiagMatrix& a, octave_idx_type r, octave_idx_type c) | |||
371 | { | |||
372 | octave_idx_type a_nr = a.rows (); | |||
373 | octave_idx_type a_nc = a.cols (); | |||
374 | ||||
375 | if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ()) | |||
376 | { | |||
377 | (*current_liboctave_error_handler) ("range error for insert"); | |||
378 | return *this; | |||
379 | } | |||
380 | ||||
381 | fill (0.0, r, c, r + a_nr - 1, c + a_nc - 1); | |||
382 | ||||
383 | octave_idx_type a_len = a.length (); | |||
384 | ||||
385 | if (a_len > 0) | |||
386 | { | |||
387 | make_unique (); | |||
388 | ||||
389 | for (octave_idx_type i = 0; i < a_len; i++) | |||
390 | xelem (r+i, c+i) = a.elem (i, i); | |||
391 | } | |||
392 | ||||
393 | return *this; | |||
394 | } | |||
395 | ||||
396 | Matrix& | |||
397 | Matrix::fill (double val) | |||
398 | { | |||
399 | octave_idx_type nr = rows (); | |||
400 | octave_idx_type nc = cols (); | |||
401 | ||||
402 | if (nr > 0 && nc > 0) | |||
403 | { | |||
404 | make_unique (); | |||
405 | ||||
406 | for (octave_idx_type j = 0; j < nc; j++) | |||
407 | for (octave_idx_type i = 0; i < nr; i++) | |||
408 | xelem (i, j) = val; | |||
409 | } | |||
410 | ||||
411 | return *this; | |||
412 | } | |||
413 | ||||
414 | Matrix& | |||
415 | Matrix::fill (double val, octave_idx_type r1, octave_idx_type c1, | |||
416 | octave_idx_type r2, octave_idx_type c2) | |||
417 | { | |||
418 | octave_idx_type nr = rows (); | |||
419 | octave_idx_type nc = cols (); | |||
420 | ||||
421 | if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0 | |||
422 | || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc) | |||
423 | { | |||
424 | (*current_liboctave_error_handler) ("range error for fill"); | |||
425 | return *this; | |||
426 | } | |||
427 | ||||
428 | if (r1 > r2) { std::swap (r1, r2); } | |||
429 | if (c1 > c2) { std::swap (c1, c2); } | |||
430 | ||||
431 | if (r2 >= r1 && c2 >= c1) | |||
432 | { | |||
433 | make_unique (); | |||
434 | ||||
435 | for (octave_idx_type j = c1; j <= c2; j++) | |||
436 | for (octave_idx_type i = r1; i <= r2; i++) | |||
437 | xelem (i, j) = val; | |||
438 | } | |||
439 | ||||
440 | return *this; | |||
441 | } | |||
442 | ||||
443 | Matrix | |||
444 | Matrix::append (const Matrix& a) const | |||
445 | { | |||
446 | octave_idx_type nr = rows (); | |||
447 | octave_idx_type nc = cols (); | |||
448 | if (nr != a.rows ()) | |||
449 | { | |||
450 | (*current_liboctave_error_handler) ("row dimension mismatch for append"); | |||
451 | return Matrix (); | |||
452 | } | |||
453 | ||||
454 | octave_idx_type nc_insert = nc; | |||
455 | Matrix retval (nr, nc + a.cols ()); | |||
456 | retval.insert (*this, 0, 0); | |||
457 | retval.insert (a, 0, nc_insert); | |||
458 | return retval; | |||
459 | } | |||
460 | ||||
461 | Matrix | |||
462 | Matrix::append (const RowVector& a) const | |||
463 | { | |||
464 | octave_idx_type nr = rows (); | |||
465 | octave_idx_type nc = cols (); | |||
466 | if (nr != 1) | |||
467 | { | |||
468 | (*current_liboctave_error_handler) ("row dimension mismatch for append"); | |||
469 | return Matrix (); | |||
470 | } | |||
471 | ||||
472 | octave_idx_type nc_insert = nc; | |||
473 | Matrix retval (nr, nc + a.length ()); | |||
474 | retval.insert (*this, 0, 0); | |||
475 | retval.insert (a, 0, nc_insert); | |||
476 | return retval; | |||
477 | } | |||
478 | ||||
479 | Matrix | |||
480 | Matrix::append (const ColumnVector& a) const | |||
481 | { | |||
482 | octave_idx_type nr = rows (); | |||
483 | octave_idx_type nc = cols (); | |||
484 | if (nr != a.length ()) | |||
485 | { | |||
486 | (*current_liboctave_error_handler) ("row dimension mismatch for append"); | |||
487 | return Matrix (); | |||
488 | } | |||
489 | ||||
490 | octave_idx_type nc_insert = nc; | |||
491 | Matrix retval (nr, nc + 1); | |||
492 | retval.insert (*this, 0, 0); | |||
493 | retval.insert (a, 0, nc_insert); | |||
494 | return retval; | |||
495 | } | |||
496 | ||||
497 | Matrix | |||
498 | Matrix::append (const DiagMatrix& a) const | |||
499 | { | |||
500 | octave_idx_type nr = rows (); | |||
501 | octave_idx_type nc = cols (); | |||
502 | if (nr != a.rows ()) | |||
503 | { | |||
504 | (*current_liboctave_error_handler) ("row dimension mismatch for append"); | |||
505 | return *this; | |||
506 | } | |||
507 | ||||
508 | octave_idx_type nc_insert = nc; | |||
509 | Matrix retval (nr, nc + a.cols ()); | |||
510 | retval.insert (*this, 0, 0); | |||
511 | retval.insert (a, 0, nc_insert); | |||
512 | return retval; | |||
513 | } | |||
514 | ||||
515 | Matrix | |||
516 | Matrix::stack (const Matrix& a) const | |||
517 | { | |||
518 | octave_idx_type nr = rows (); | |||
519 | octave_idx_type nc = cols (); | |||
520 | if (nc != a.cols ()) | |||
521 | { | |||
522 | (*current_liboctave_error_handler) | |||
523 | ("column dimension mismatch for stack"); | |||
524 | return Matrix (); | |||
525 | } | |||
526 | ||||
527 | octave_idx_type nr_insert = nr; | |||
528 | Matrix retval (nr + a.rows (), nc); | |||
529 | retval.insert (*this, 0, 0); | |||
530 | retval.insert (a, nr_insert, 0); | |||
531 | return retval; | |||
532 | } | |||
533 | ||||
534 | Matrix | |||
535 | Matrix::stack (const RowVector& a) const | |||
536 | { | |||
537 | octave_idx_type nr = rows (); | |||
538 | octave_idx_type nc = cols (); | |||
539 | if (nc != a.length ()) | |||
540 | { | |||
541 | (*current_liboctave_error_handler) | |||
542 | ("column dimension mismatch for stack"); | |||
543 | return Matrix (); | |||
544 | } | |||
545 | ||||
546 | octave_idx_type nr_insert = nr; | |||
547 | Matrix retval (nr + 1, nc); | |||
548 | retval.insert (*this, 0, 0); | |||
549 | retval.insert (a, nr_insert, 0); | |||
550 | return retval; | |||
551 | } | |||
552 | ||||
553 | Matrix | |||
554 | Matrix::stack (const ColumnVector& a) const | |||
555 | { | |||
556 | octave_idx_type nr = rows (); | |||
557 | octave_idx_type nc = cols (); | |||
558 | if (nc != 1) | |||
559 | { | |||
560 | (*current_liboctave_error_handler) | |||
561 | ("column dimension mismatch for stack"); | |||
562 | return Matrix (); | |||
563 | } | |||
564 | ||||
565 | octave_idx_type nr_insert = nr; | |||
566 | Matrix retval (nr + a.length (), nc); | |||
567 | retval.insert (*this, 0, 0); | |||
568 | retval.insert (a, nr_insert, 0); | |||
569 | return retval; | |||
570 | } | |||
571 | ||||
572 | Matrix | |||
573 | Matrix::stack (const DiagMatrix& a) const | |||
574 | { | |||
575 | octave_idx_type nr = rows (); | |||
576 | octave_idx_type nc = cols (); | |||
577 | if (nc != a.cols ()) | |||
578 | { | |||
579 | (*current_liboctave_error_handler) | |||
580 | ("column dimension mismatch for stack"); | |||
581 | return Matrix (); | |||
582 | } | |||
583 | ||||
584 | octave_idx_type nr_insert = nr; | |||
585 | Matrix retval (nr + a.rows (), nc); | |||
586 | retval.insert (*this, 0, 0); | |||
587 | retval.insert (a, nr_insert, 0); | |||
588 | return retval; | |||
589 | } | |||
590 | ||||
591 | Matrix | |||
592 | real (const ComplexMatrix& a) | |||
593 | { | |||
594 | return do_mx_unary_op<double, Complex> (a, mx_inline_real); | |||
595 | } | |||
596 | ||||
597 | Matrix | |||
598 | imag (const ComplexMatrix& a) | |||
599 | { | |||
600 | return do_mx_unary_op<double, Complex> (a, mx_inline_imag); | |||
601 | } | |||
602 | ||||
603 | Matrix | |||
604 | Matrix::extract (octave_idx_type r1, octave_idx_type c1, | |||
605 | octave_idx_type r2, octave_idx_type c2) const | |||
606 | { | |||
607 | if (r1 > r2) { std::swap (r1, r2); } | |||
608 | if (c1 > c2) { std::swap (c1, c2); } | |||
609 | ||||
610 | return index (idx_vector (r1, r2+1), idx_vector (c1, c2+1)); | |||
611 | } | |||
612 | ||||
613 | Matrix | |||
614 | Matrix::extract_n (octave_idx_type r1, octave_idx_type c1, octave_idx_type nr, | |||
615 | octave_idx_type nc) const | |||
616 | { | |||
617 | return index (idx_vector (r1, r1 + nr), idx_vector (c1, c1 + nc)); | |||
618 | } | |||
619 | ||||
620 | // extract row or column i. | |||
621 | ||||
622 | RowVector | |||
623 | Matrix::row (octave_idx_type i) const | |||
624 | { | |||
625 | return index (idx_vector (i), idx_vector::colon); | |||
626 | } | |||
627 | ||||
628 | ColumnVector | |||
629 | Matrix::column (octave_idx_type i) const | |||
630 | { | |||
631 | return index (idx_vector::colon, idx_vector (i)); | |||
632 | } | |||
633 | ||||
634 | Matrix | |||
635 | Matrix::inverse (void) const | |||
636 | { | |||
637 | octave_idx_type info; | |||
638 | double rcon; | |||
639 | MatrixType mattype (*this); | |||
640 | return inverse (mattype, info, rcon, 0, 0); | |||
641 | } | |||
642 | ||||
643 | Matrix | |||
644 | Matrix::inverse (octave_idx_type& info) const | |||
645 | { | |||
646 | double rcon; | |||
647 | MatrixType mattype (*this); | |||
648 | return inverse (mattype, info, rcon, 0, 0); | |||
649 | } | |||
650 | ||||
651 | Matrix | |||
652 | Matrix::inverse (octave_idx_type& info, double& rcon, int force, | |||
653 | int calc_cond) const | |||
654 | { | |||
655 | MatrixType mattype (*this); | |||
656 | return inverse (mattype, info, rcon, force, calc_cond); | |||
657 | } | |||
658 | ||||
659 | Matrix | |||
660 | Matrix::inverse (MatrixType& mattype) const | |||
661 | { | |||
662 | octave_idx_type info; | |||
663 | double rcon; | |||
664 | return inverse (mattype, info, rcon, 0, 0); | |||
665 | } | |||
666 | ||||
667 | Matrix | |||
668 | Matrix::inverse (MatrixType &mattype, octave_idx_type& info) const | |||
669 | { | |||
670 | double rcon; | |||
671 | return inverse (mattype, info, rcon, 0, 0); | |||
672 | } | |||
673 | ||||
674 | Matrix | |||
675 | Matrix::tinverse (MatrixType &mattype, octave_idx_type& info, double& rcon, | |||
676 | int force, int calc_cond) const | |||
677 | { | |||
678 | Matrix retval; | |||
679 | ||||
680 | octave_idx_type nr = rows (); | |||
681 | octave_idx_type nc = cols (); | |||
682 | ||||
683 | if (nr != nc || nr == 0 || nc == 0) | |||
684 | (*current_liboctave_error_handler) ("inverse requires square matrix"); | |||
685 | else | |||
686 | { | |||
687 | int typ = mattype.type (); | |||
688 | char uplo = (typ == MatrixType::Lower ? 'L' : 'U'); | |||
689 | char udiag = 'N'; | |||
690 | retval = *this; | |||
691 | double *tmp_data = retval.fortran_vec (); | |||
692 | ||||
693 | F77_XFCN (dtrtri, DTRTRI, (F77_CONST_CHAR_ARG2 (&uplo, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrtri_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrtri_ (&uplo, &udiag, nr, tmp_data, nr, info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
694 | F77_CONST_CHAR_ARG2 (&udiag, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrtri_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrtri_ (&uplo, &udiag, nr, tmp_data, nr, info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
695 | nr, tmp_data, nr, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrtri_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrtri_ (&uplo, &udiag, nr, tmp_data, nr, info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
696 | F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrtri_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrtri_ (&uplo, &udiag, nr, tmp_data, nr, info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
697 | F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrtri_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrtri_ (&uplo, &udiag, nr, tmp_data, nr, info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0); | |||
698 | ||||
699 | // Throw-away extra info LAPACK gives so as to not change output. | |||
700 | rcon = 0.0; | |||
701 | if (info != 0) | |||
702 | info = -1; | |||
703 | else if (calc_cond) | |||
704 | { | |||
705 | octave_idx_type dtrcon_info = 0; | |||
706 | char job = '1'; | |||
707 | ||||
708 | OCTAVE_LOCAL_BUFFER (double, work, 3 * nr)octave_local_buffer<double> _buffer_work (3 * nr); double *work = _buffer_work; | |||
709 | OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, nr)octave_local_buffer<octave_idx_type> _buffer_iwork (nr) ; octave_idx_type *iwork = _buffer_iwork; | |||
710 | ||||
711 | F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&job, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&job, &uplo, &udiag, nr, tmp_data, nr , rcon, work, iwork, dtrcon_info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
712 | F77_CONST_CHAR_ARG2 (&uplo, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&job, &uplo, &udiag, nr, tmp_data, nr , rcon, work, iwork, dtrcon_info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
713 | F77_CONST_CHAR_ARG2 (&udiag, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&job, &uplo, &udiag, nr, tmp_data, nr , rcon, work, iwork, dtrcon_info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
714 | nr, tmp_data, nr, rcon,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&job, &uplo, &udiag, nr, tmp_data, nr , rcon, work, iwork, dtrcon_info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
715 | work, iwork, dtrcon_infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&job, &uplo, &udiag, nr, tmp_data, nr , rcon, work, iwork, dtrcon_info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
716 | F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&job, &uplo, &udiag, nr, tmp_data, nr , rcon, work, iwork, dtrcon_info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
717 | F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&job, &uplo, &udiag, nr, tmp_data, nr , rcon, work, iwork, dtrcon_info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
718 | F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&job, &uplo, &udiag, nr, tmp_data, nr , rcon, work, iwork, dtrcon_info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0); | |||
719 | ||||
720 | if (dtrcon_info != 0) | |||
721 | info = -1; | |||
722 | } | |||
723 | ||||
724 | if (info == -1 && ! force) | |||
725 | retval = *this; // Restore matrix contents. | |||
726 | } | |||
727 | ||||
728 | return retval; | |||
729 | } | |||
730 | ||||
731 | ||||
732 | Matrix | |||
733 | Matrix::finverse (MatrixType &mattype, octave_idx_type& info, double& rcon, | |||
734 | int force, int calc_cond) const | |||
735 | { | |||
736 | Matrix retval; | |||
737 | ||||
738 | octave_idx_type nr = rows (); | |||
739 | octave_idx_type nc = cols (); | |||
740 | ||||
741 | if (nr != nc || nr == 0 || nc == 0) | |||
742 | (*current_liboctave_error_handler) ("inverse requires square matrix"); | |||
743 | else | |||
744 | { | |||
745 | Array<octave_idx_type> ipvt (dim_vector (nr, 1)); | |||
746 | octave_idx_type *pipvt = ipvt.fortran_vec (); | |||
747 | ||||
748 | retval = *this; | |||
749 | double *tmp_data = retval.fortran_vec (); | |||
750 | ||||
751 | Array<double> z (dim_vector (1, 1)); | |||
752 | octave_idx_type lwork = -1; | |||
753 | ||||
754 | // Query the optimum work array size. | |||
755 | F77_XFCN (dgetri, DGETRI, (nc, tmp_data, nr, pipvt,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgetri_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgetri_ (nc, tmp_data, nr, pipvt, z.fortran_vec (), lwork , info); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
756 | z.fortran_vec (), lwork, info))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgetri_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgetri_ (nc, tmp_data, nr, pipvt, z.fortran_vec (), lwork , info); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0); | |||
757 | ||||
758 | lwork = static_cast<octave_idx_type> (z(0)); | |||
759 | lwork = (lwork < 2 *nc ? 2*nc : lwork); | |||
760 | z.resize (dim_vector (lwork, 1)); | |||
761 | double *pz = z.fortran_vec (); | |||
762 | ||||
763 | info = 0; | |||
764 | ||||
765 | // Calculate the norm of the matrix, for later use. | |||
766 | double anorm = 0; | |||
767 | if (calc_cond) | |||
768 | anorm = retval.abs ().sum ().row (static_cast<octave_idx_type>(0)) | |||
769 | .max (); | |||
770 | ||||
771 | F77_XFCN (dgetrf, DGETRF, (nc, nc, tmp_data, nr, pipvt, info))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgetrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgetrf_ (nc, nc, tmp_data, nr, pipvt, info); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0); | |||
772 | ||||
773 | // Throw-away extra info LAPACK gives so as to not change output. | |||
774 | rcon = 0.0; | |||
775 | if (info != 0) | |||
776 | info = -1; | |||
777 | else if (calc_cond) | |||
778 | { | |||
779 | octave_idx_type dgecon_info = 0; | |||
780 | ||||
781 | // Now calculate the condition number for non-singular matrix. | |||
782 | char job = '1'; | |||
783 | Array<octave_idx_type> iz (dim_vector (nc, 1)); | |||
784 | octave_idx_type *piz = iz.fortran_vec (); | |||
785 | F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgecon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgecon_ (&job, nc, tmp_data, nr, anorm, rcon, pz, piz , dgecon_info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
786 | nc, tmp_data, nr, anorm,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgecon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgecon_ (&job, nc, tmp_data, nr, anorm, rcon, pz, piz , dgecon_info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
787 | rcon, pz, piz, dgecon_infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgecon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgecon_ (&job, nc, tmp_data, nr, anorm, rcon, pz, piz , dgecon_info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
788 | F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgecon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgecon_ (&job, nc, tmp_data, nr, anorm, rcon, pz, piz , dgecon_info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0); | |||
789 | ||||
790 | if (dgecon_info != 0) | |||
791 | info = -1; | |||
792 | } | |||
793 | ||||
794 | if (info == -1 && ! force) | |||
795 | retval = *this; // Restore matrix contents. | |||
796 | else | |||
797 | { | |||
798 | octave_idx_type dgetri_info = 0; | |||
799 | ||||
800 | F77_XFCN (dgetri, DGETRI, (nc, tmp_data, nr, pipvt,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgetri_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgetri_ (nc, tmp_data, nr, pipvt, pz, lwork, dgetri_info) ; octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
801 | pz, lwork, dgetri_info))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgetri_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgetri_ (nc, tmp_data, nr, pipvt, pz, lwork, dgetri_info) ; octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0); | |||
802 | ||||
803 | if (dgetri_info != 0) | |||
804 | info = -1; | |||
805 | } | |||
806 | ||||
807 | if (info != 0) | |||
808 | mattype.mark_as_rectangular (); | |||
809 | } | |||
810 | ||||
811 | return retval; | |||
812 | } | |||
813 | ||||
814 | Matrix | |||
815 | Matrix::inverse (MatrixType &mattype, octave_idx_type& info, double& rcon, | |||
816 | int force, int calc_cond) const | |||
817 | { | |||
818 | int typ = mattype.type (false); | |||
819 | Matrix ret; | |||
820 | ||||
821 | if (typ == MatrixType::Unknown) | |||
822 | typ = mattype.type (*this); | |||
823 | ||||
824 | if (typ == MatrixType::Upper || typ == MatrixType::Lower) | |||
825 | ret = tinverse (mattype, info, rcon, force, calc_cond); | |||
826 | else | |||
827 | { | |||
828 | if (mattype.is_hermitian ()) | |||
829 | { | |||
830 | CHOL chol (*this, info, calc_cond); | |||
831 | if (info == 0) | |||
832 | { | |||
833 | if (calc_cond) | |||
834 | rcon = chol.rcond (); | |||
835 | else | |||
836 | rcon = 1.0; | |||
837 | ret = chol.inverse (); | |||
838 | } | |||
839 | else | |||
840 | mattype.mark_as_unsymmetric (); | |||
841 | } | |||
842 | ||||
843 | if (!mattype.is_hermitian ()) | |||
844 | ret = finverse (mattype, info, rcon, force, calc_cond); | |||
845 | ||||
846 | if ((mattype.is_hermitian () || calc_cond) && rcon == 0.) | |||
847 | ret = Matrix (rows (), columns (), octave_Inf); | |||
848 | } | |||
849 | ||||
850 | return ret; | |||
851 | } | |||
852 | ||||
853 | Matrix | |||
854 | Matrix::pseudo_inverse (double tol) const | |||
855 | { | |||
856 | SVD result (*this, SVD::economy); | |||
857 | ||||
858 | DiagMatrix S = result.singular_values (); | |||
859 | Matrix U = result.left_singular_matrix (); | |||
860 | Matrix V = result.right_singular_matrix (); | |||
861 | ||||
862 | ColumnVector sigma = S.extract_diag (); | |||
863 | ||||
864 | octave_idx_type r = sigma.length () - 1; | |||
865 | octave_idx_type nr = rows (); | |||
866 | octave_idx_type nc = cols (); | |||
867 | ||||
868 | if (tol <= 0.0) | |||
869 | { | |||
870 | if (nr > nc) | |||
871 | tol = nr * sigma.elem (0) * std::numeric_limits<double>::epsilon (); | |||
872 | else | |||
873 | tol = nc * sigma.elem (0) * std::numeric_limits<double>::epsilon (); | |||
874 | } | |||
875 | ||||
876 | while (r >= 0 && sigma.elem (r) < tol) | |||
877 | r--; | |||
878 | ||||
879 | if (r < 0) | |||
880 | return Matrix (nc, nr, 0.0); | |||
881 | else | |||
882 | { | |||
883 | Matrix Ur = U.extract (0, 0, nr-1, r); | |||
884 | DiagMatrix D = DiagMatrix (sigma.extract (0, r)) . inverse (); | |||
885 | Matrix Vr = V.extract (0, 0, nc-1, r); | |||
886 | return Vr * D * Ur.transpose (); | |||
887 | } | |||
888 | } | |||
889 | ||||
890 | #if defined (HAVE_FFTW1) | |||
891 | ||||
892 | ComplexMatrix | |||
893 | Matrix::fourier (void) const | |||
894 | { | |||
895 | size_t nr = rows (); | |||
896 | size_t nc = cols (); | |||
897 | ||||
898 | ComplexMatrix retval (nr, nc); | |||
899 | ||||
900 | size_t npts, nsamples; | |||
901 | ||||
902 | if (nr == 1 || nc == 1) | |||
903 | { | |||
904 | npts = nr > nc ? nr : nc; | |||
905 | nsamples = 1; | |||
906 | } | |||
907 | else | |||
908 | { | |||
909 | npts = nr; | |||
910 | nsamples = nc; | |||
911 | } | |||
912 | ||||
913 | const double *in (fortran_vec ()); | |||
914 | Complex *out (retval.fortran_vec ()); | |||
915 | ||||
916 | octave_fftw::fft (in, out, npts, nsamples); | |||
917 | ||||
918 | return retval; | |||
919 | } | |||
920 | ||||
921 | ComplexMatrix | |||
922 | Matrix::ifourier (void) const | |||
923 | { | |||
924 | size_t nr = rows (); | |||
925 | size_t nc = cols (); | |||
926 | ||||
927 | ComplexMatrix retval (nr, nc); | |||
928 | ||||
929 | size_t npts, nsamples; | |||
930 | ||||
931 | if (nr == 1 || nc == 1) | |||
932 | { | |||
933 | npts = nr > nc ? nr : nc; | |||
934 | nsamples = 1; | |||
935 | } | |||
936 | else | |||
937 | { | |||
938 | npts = nr; | |||
939 | nsamples = nc; | |||
940 | } | |||
941 | ||||
942 | ComplexMatrix tmp (*this); | |||
943 | Complex *in (tmp.fortran_vec ()); | |||
944 | Complex *out (retval.fortran_vec ()); | |||
945 | ||||
946 | octave_fftw::ifft (in, out, npts, nsamples); | |||
947 | ||||
948 | return retval; | |||
949 | } | |||
950 | ||||
951 | ComplexMatrix | |||
952 | Matrix::fourier2d (void) const | |||
953 | { | |||
954 | dim_vector dv(rows (), cols ()); | |||
955 | ||||
956 | const double *in = fortran_vec (); | |||
957 | ComplexMatrix retval (rows (), cols ()); | |||
958 | octave_fftw::fftNd (in, retval.fortran_vec (), 2, dv); | |||
959 | ||||
960 | return retval; | |||
961 | } | |||
962 | ||||
963 | ComplexMatrix | |||
964 | Matrix::ifourier2d (void) const | |||
965 | { | |||
966 | dim_vector dv(rows (), cols ()); | |||
967 | ||||
968 | ComplexMatrix retval (*this); | |||
969 | Complex *out (retval.fortran_vec ()); | |||
970 | ||||
971 | octave_fftw::ifftNd (out, out, 2, dv); | |||
972 | ||||
973 | return retval; | |||
974 | } | |||
975 | ||||
976 | #else | |||
977 | ||||
978 | extern "C" | |||
979 | { | |||
980 | // Note that the original complex fft routines were not written for | |||
981 | // double complex arguments. They have been modified by adding an | |||
982 | // implicit double precision (a-h,o-z) statement at the beginning of | |||
983 | // each subroutine. | |||
984 | ||||
985 | F77_RET_Tint | |||
986 | F77_FUNC (zffti, ZFFTI)zffti_ (const octave_idx_type&, Complex*); | |||
987 | ||||
988 | F77_RET_Tint | |||
989 | F77_FUNC (zfftf, ZFFTF)zfftf_ (const octave_idx_type&, Complex*, Complex*); | |||
990 | ||||
991 | F77_RET_Tint | |||
992 | F77_FUNC (zfftb, ZFFTB)zfftb_ (const octave_idx_type&, Complex*, Complex*); | |||
993 | } | |||
994 | ||||
995 | ComplexMatrix | |||
996 | Matrix::fourier (void) const | |||
997 | { | |||
998 | ComplexMatrix retval; | |||
999 | ||||
1000 | octave_idx_type nr = rows (); | |||
1001 | octave_idx_type nc = cols (); | |||
1002 | ||||
1003 | octave_idx_type npts, nsamples; | |||
1004 | ||||
1005 | if (nr == 1 || nc == 1) | |||
1006 | { | |||
1007 | npts = nr > nc ? nr : nc; | |||
1008 | nsamples = 1; | |||
1009 | } | |||
1010 | else | |||
1011 | { | |||
1012 | npts = nr; | |||
1013 | nsamples = nc; | |||
1014 | } | |||
1015 | ||||
1016 | octave_idx_type nn = 4*npts+15; | |||
1017 | ||||
1018 | Array<Complex> wsave (dim_vector (nn, 1)); | |||
1019 | Complex *pwsave = wsave.fortran_vec (); | |||
1020 | ||||
1021 | retval = ComplexMatrix (*this); | |||
1022 | Complex *tmp_data = retval.fortran_vec (); | |||
1023 | ||||
1024 | F77_FUNC (zffti, ZFFTI)zffti_ (npts, pwsave); | |||
1025 | ||||
1026 | for (octave_idx_type j = 0; j < nsamples; j++) | |||
1027 | { | |||
1028 | octave_quit (); | |||
1029 | ||||
1030 | F77_FUNC (zfftf, ZFFTF)zfftf_ (npts, &tmp_data[npts*j], pwsave); | |||
1031 | } | |||
1032 | ||||
1033 | return retval; | |||
1034 | } | |||
1035 | ||||
1036 | ComplexMatrix | |||
1037 | Matrix::ifourier (void) const | |||
1038 | { | |||
1039 | ComplexMatrix retval; | |||
1040 | ||||
1041 | octave_idx_type nr = rows (); | |||
1042 | octave_idx_type nc = cols (); | |||
1043 | ||||
1044 | octave_idx_type npts, nsamples; | |||
1045 | ||||
1046 | if (nr == 1 || nc == 1) | |||
1047 | { | |||
1048 | npts = nr > nc ? nr : nc; | |||
1049 | nsamples = 1; | |||
1050 | } | |||
1051 | else | |||
1052 | { | |||
1053 | npts = nr; | |||
1054 | nsamples = nc; | |||
1055 | } | |||
1056 | ||||
1057 | octave_idx_type nn = 4*npts+15; | |||
1058 | ||||
1059 | Array<Complex> wsave (dim_vector (nn, 1)); | |||
1060 | Complex *pwsave = wsave.fortran_vec (); | |||
1061 | ||||
1062 | retval = ComplexMatrix (*this); | |||
1063 | Complex *tmp_data = retval.fortran_vec (); | |||
1064 | ||||
1065 | F77_FUNC (zffti, ZFFTI)zffti_ (npts, pwsave); | |||
1066 | ||||
1067 | for (octave_idx_type j = 0; j < nsamples; j++) | |||
1068 | { | |||
1069 | octave_quit (); | |||
1070 | ||||
1071 | F77_FUNC (zfftb, ZFFTB)zfftb_ (npts, &tmp_data[npts*j], pwsave); | |||
1072 | } | |||
1073 | ||||
1074 | for (octave_idx_type j = 0; j < npts*nsamples; j++) | |||
1075 | tmp_data[j] = tmp_data[j] / static_cast<double> (npts); | |||
1076 | ||||
1077 | return retval; | |||
1078 | } | |||
1079 | ||||
1080 | ComplexMatrix | |||
1081 | Matrix::fourier2d (void) const | |||
1082 | { | |||
1083 | ComplexMatrix retval; | |||
1084 | ||||
1085 | octave_idx_type nr = rows (); | |||
1086 | octave_idx_type nc = cols (); | |||
1087 | ||||
1088 | octave_idx_type npts, nsamples; | |||
1089 | ||||
1090 | if (nr == 1 || nc == 1) | |||
1091 | { | |||
1092 | npts = nr > nc ? nr : nc; | |||
1093 | nsamples = 1; | |||
1094 | } | |||
1095 | else | |||
1096 | { | |||
1097 | npts = nr; | |||
1098 | nsamples = nc; | |||
1099 | } | |||
1100 | ||||
1101 | octave_idx_type nn = 4*npts+15; | |||
1102 | ||||
1103 | Array<Complex> wsave (dim_vector (nn, 1)); | |||
1104 | Complex *pwsave = wsave.fortran_vec (); | |||
1105 | ||||
1106 | retval = ComplexMatrix (*this); | |||
1107 | Complex *tmp_data = retval.fortran_vec (); | |||
1108 | ||||
1109 | F77_FUNC (zffti, ZFFTI)zffti_ (npts, pwsave); | |||
1110 | ||||
1111 | for (octave_idx_type j = 0; j < nsamples; j++) | |||
1112 | { | |||
1113 | octave_quit (); | |||
1114 | ||||
1115 | F77_FUNC (zfftf, ZFFTF)zfftf_ (npts, &tmp_data[npts*j], pwsave); | |||
1116 | } | |||
1117 | ||||
1118 | npts = nc; | |||
1119 | nsamples = nr; | |||
1120 | nn = 4*npts+15; | |||
1121 | ||||
1122 | wsave.resize (dim_vector (nn, 1)); | |||
1123 | pwsave = wsave.fortran_vec (); | |||
1124 | ||||
1125 | Array<Complex> tmp (dim_vector (npts, 1)); | |||
1126 | Complex *prow = tmp.fortran_vec (); | |||
1127 | ||||
1128 | F77_FUNC (zffti, ZFFTI)zffti_ (npts, pwsave); | |||
1129 | ||||
1130 | for (octave_idx_type j = 0; j < nsamples; j++) | |||
1131 | { | |||
1132 | octave_quit (); | |||
1133 | ||||
1134 | for (octave_idx_type i = 0; i < npts; i++) | |||
1135 | prow[i] = tmp_data[i*nr + j]; | |||
1136 | ||||
1137 | F77_FUNC (zfftf, ZFFTF)zfftf_ (npts, prow, pwsave); | |||
1138 | ||||
1139 | for (octave_idx_type i = 0; i < npts; i++) | |||
1140 | tmp_data[i*nr + j] = prow[i]; | |||
1141 | } | |||
1142 | ||||
1143 | return retval; | |||
1144 | } | |||
1145 | ||||
1146 | ComplexMatrix | |||
1147 | Matrix::ifourier2d (void) const | |||
1148 | { | |||
1149 | ComplexMatrix retval; | |||
1150 | ||||
1151 | octave_idx_type nr = rows (); | |||
1152 | octave_idx_type nc = cols (); | |||
1153 | ||||
1154 | octave_idx_type npts, nsamples; | |||
1155 | ||||
1156 | if (nr == 1 || nc == 1) | |||
1157 | { | |||
1158 | npts = nr > nc ? nr : nc; | |||
1159 | nsamples = 1; | |||
1160 | } | |||
1161 | else | |||
1162 | { | |||
1163 | npts = nr; | |||
1164 | nsamples = nc; | |||
1165 | } | |||
1166 | ||||
1167 | octave_idx_type nn = 4*npts+15; | |||
1168 | ||||
1169 | Array<Complex> wsave (dim_vector (nn, 1)); | |||
1170 | Complex *pwsave = wsave.fortran_vec (); | |||
1171 | ||||
1172 | retval = ComplexMatrix (*this); | |||
1173 | Complex *tmp_data = retval.fortran_vec (); | |||
1174 | ||||
1175 | F77_FUNC (zffti, ZFFTI)zffti_ (npts, pwsave); | |||
1176 | ||||
1177 | for (octave_idx_type j = 0; j < nsamples; j++) | |||
1178 | { | |||
1179 | octave_quit (); | |||
1180 | ||||
1181 | F77_FUNC (zfftb, ZFFTB)zfftb_ (npts, &tmp_data[npts*j], pwsave); | |||
1182 | } | |||
1183 | ||||
1184 | for (octave_idx_type j = 0; j < npts*nsamples; j++) | |||
1185 | tmp_data[j] = tmp_data[j] / static_cast<double> (npts); | |||
1186 | ||||
1187 | npts = nc; | |||
1188 | nsamples = nr; | |||
1189 | nn = 4*npts+15; | |||
1190 | ||||
1191 | wsave.resize (dim_vector (nn, 1)); | |||
1192 | pwsave = wsave.fortran_vec (); | |||
1193 | ||||
1194 | Array<Complex> tmp (dim_vector (npts, 1)); | |||
1195 | Complex *prow = tmp.fortran_vec (); | |||
1196 | ||||
1197 | F77_FUNC (zffti, ZFFTI)zffti_ (npts, pwsave); | |||
1198 | ||||
1199 | for (octave_idx_type j = 0; j < nsamples; j++) | |||
1200 | { | |||
1201 | octave_quit (); | |||
1202 | ||||
1203 | for (octave_idx_type i = 0; i < npts; i++) | |||
1204 | prow[i] = tmp_data[i*nr + j]; | |||
1205 | ||||
1206 | F77_FUNC (zfftb, ZFFTB)zfftb_ (npts, prow, pwsave); | |||
1207 | ||||
1208 | for (octave_idx_type i = 0; i < npts; i++) | |||
1209 | tmp_data[i*nr + j] = prow[i] / static_cast<double> (npts); | |||
1210 | } | |||
1211 | ||||
1212 | return retval; | |||
1213 | } | |||
1214 | ||||
1215 | #endif | |||
1216 | ||||
1217 | DET | |||
1218 | Matrix::determinant (void) const | |||
1219 | { | |||
1220 | octave_idx_type info; | |||
1221 | double rcon; | |||
1222 | return determinant (info, rcon, 0); | |||
1223 | } | |||
1224 | ||||
1225 | DET | |||
1226 | Matrix::determinant (octave_idx_type& info) const | |||
1227 | { | |||
1228 | double rcon; | |||
1229 | return determinant (info, rcon, 0); | |||
1230 | } | |||
1231 | ||||
1232 | DET | |||
1233 | Matrix::determinant (octave_idx_type& info, double& rcon, int calc_cond) const | |||
1234 | { | |||
1235 | MatrixType mattype (*this); | |||
1236 | return determinant (mattype, info, rcon, calc_cond); | |||
1237 | } | |||
1238 | ||||
1239 | DET | |||
1240 | Matrix::determinant (MatrixType& mattype, | |||
1241 | octave_idx_type& info, double& rcon, int calc_cond) const | |||
1242 | { | |||
1243 | DET retval (1.0); | |||
1244 | ||||
1245 | info = 0; | |||
1246 | rcon = 0.0; | |||
1247 | ||||
1248 | octave_idx_type nr = rows (); | |||
1249 | octave_idx_type nc = cols (); | |||
1250 | ||||
1251 | if (nr != nc) | |||
1252 | (*current_liboctave_error_handler) ("matrix must be square"); | |||
1253 | else | |||
1254 | { | |||
1255 | volatile int typ = mattype.type (); | |||
1256 | ||||
1257 | // Even though the matrix is marked as singular (Rectangular), we may | |||
1258 | // still get a useful number from the LU factorization, because it always | |||
1259 | // completes. | |||
1260 | ||||
1261 | if (typ == MatrixType::Unknown) | |||
1262 | typ = mattype.type (*this); | |||
1263 | else if (typ == MatrixType::Rectangular) | |||
1264 | typ = MatrixType::Full; | |||
1265 | ||||
1266 | if (typ == MatrixType::Lower || typ == MatrixType::Upper) | |||
1267 | { | |||
1268 | for (octave_idx_type i = 0; i < nc; i++) | |||
1269 | retval *= elem (i,i); | |||
1270 | } | |||
1271 | else if (typ == MatrixType::Hermitian) | |||
1272 | { | |||
1273 | Matrix atmp = *this; | |||
1274 | double *tmp_data = atmp.fortran_vec (); | |||
1275 | ||||
1276 | double anorm = 0; | |||
1277 | if (calc_cond) anorm = xnorm (*this, 1); | |||
1278 | ||||
1279 | ||||
1280 | char job = 'L'; | |||
1281 | F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dpotrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dpotrf_ (&job, nr, tmp_data, nr, info , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1282 | tmp_data, nr, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dpotrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dpotrf_ (&job, nr, tmp_data, nr, info , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1283 | F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dpotrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dpotrf_ (&job, nr, tmp_data, nr, info , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0); | |||
1284 | ||||
1285 | if (info != 0) | |||
1286 | { | |||
1287 | rcon = 0.0; | |||
1288 | mattype.mark_as_unsymmetric (); | |||
1289 | typ = MatrixType::Full; | |||
1290 | } | |||
1291 | else | |||
1292 | { | |||
1293 | Array<double> z (dim_vector (3 * nc, 1)); | |||
1294 | double *pz = z.fortran_vec (); | |||
1295 | Array<octave_idx_type> iz (dim_vector (nc, 1)); | |||
1296 | octave_idx_type *piz = iz.fortran_vec (); | |||
1297 | ||||
1298 | F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dpocon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dpocon_ (&job, nr, tmp_data, nr, anorm, rcon, pz, piz , info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
1299 | nr, tmp_data, nr, anorm,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dpocon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dpocon_ (&job, nr, tmp_data, nr, anorm, rcon, pz, piz , info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
1300 | rcon, pz, piz, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dpocon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dpocon_ (&job, nr, tmp_data, nr, anorm, rcon, pz, piz , info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
1301 | F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dpocon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dpocon_ (&job, nr, tmp_data, nr, anorm, rcon, pz, piz , info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0); | |||
1302 | ||||
1303 | if (info != 0) | |||
1304 | rcon = 0.0; | |||
1305 | ||||
1306 | for (octave_idx_type i = 0; i < nc; i++) | |||
1307 | retval *= atmp (i,i); | |||
1308 | ||||
1309 | retval = retval.square (); | |||
1310 | } | |||
1311 | } | |||
1312 | else if (typ != MatrixType::Full) | |||
1313 | (*current_liboctave_error_handler) ("det: invalid dense matrix type"); | |||
1314 | ||||
1315 | if (typ == MatrixType::Full) | |||
1316 | { | |||
1317 | Array<octave_idx_type> ipvt (dim_vector (nr, 1)); | |||
1318 | octave_idx_type *pipvt = ipvt.fortran_vec (); | |||
1319 | ||||
1320 | Matrix atmp = *this; | |||
1321 | double *tmp_data = atmp.fortran_vec (); | |||
1322 | ||||
1323 | info = 0; | |||
1324 | ||||
1325 | // Calculate the norm of the matrix, for later use. | |||
1326 | double anorm = 0; | |||
1327 | if (calc_cond) anorm = xnorm (*this, 1); | |||
1328 | ||||
1329 | F77_XFCN (dgetrf, DGETRF, (nr, nr, tmp_data, nr, pipvt, info))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgetrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgetrf_ (nr, nr, tmp_data, nr, pipvt, info); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0); | |||
1330 | ||||
1331 | // Throw-away extra info LAPACK gives so as to not change output. | |||
1332 | rcon = 0.0; | |||
1333 | if (info != 0) | |||
1334 | { | |||
1335 | info = -1; | |||
1336 | retval = DET (); | |||
1337 | } | |||
1338 | else | |||
1339 | { | |||
1340 | if (calc_cond) | |||
1341 | { | |||
1342 | // Now calc the condition number for non-singular matrix. | |||
1343 | char job = '1'; | |||
1344 | Array<double> z (dim_vector (4 * nc, 1)); | |||
1345 | double *pz = z.fortran_vec (); | |||
1346 | Array<octave_idx_type> iz (dim_vector (nc, 1)); | |||
1347 | octave_idx_type *piz = iz.fortran_vec (); | |||
1348 | ||||
1349 | F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgecon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgecon_ (&job, nc, tmp_data, nr, anorm, rcon, pz, piz , info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
1350 | nc, tmp_data, nr, anorm,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgecon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgecon_ (&job, nc, tmp_data, nr, anorm, rcon, pz, piz , info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
1351 | rcon, pz, piz, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgecon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgecon_ (&job, nc, tmp_data, nr, anorm, rcon, pz, piz , info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
1352 | F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgecon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgecon_ (&job, nc, tmp_data, nr, anorm, rcon, pz, piz , info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0); | |||
1353 | } | |||
1354 | ||||
1355 | if (info != 0) | |||
1356 | { | |||
1357 | info = -1; | |||
1358 | retval = DET (); | |||
1359 | } | |||
1360 | else | |||
1361 | { | |||
1362 | for (octave_idx_type i = 0; i < nc; i++) | |||
1363 | { | |||
1364 | double c = atmp(i,i); | |||
1365 | retval *= (ipvt(i) != (i+1)) ? -c : c; | |||
1366 | } | |||
1367 | } | |||
1368 | } | |||
1369 | } | |||
1370 | } | |||
1371 | ||||
1372 | return retval; | |||
1373 | } | |||
1374 | ||||
1375 | double | |||
1376 | Matrix::rcond (void) const | |||
1377 | { | |||
1378 | MatrixType mattype (*this); | |||
1379 | return rcond (mattype); | |||
1380 | } | |||
1381 | ||||
1382 | double | |||
1383 | Matrix::rcond (MatrixType &mattype) const | |||
1384 | { | |||
1385 | double rcon; | |||
| ||||
1386 | octave_idx_type nr = rows (); | |||
1387 | octave_idx_type nc = cols (); | |||
1388 | ||||
1389 | if (nr != nc) | |||
1390 | (*current_liboctave_error_handler) ("matrix must be square"); | |||
1391 | else if (nr == 0 || nc == 0) | |||
1392 | rcon = octave_Inf; | |||
1393 | else | |||
1394 | { | |||
1395 | volatile int typ = mattype.type (); | |||
1396 | ||||
1397 | if (typ == MatrixType::Unknown) | |||
1398 | typ = mattype.type (*this); | |||
1399 | ||||
1400 | // Only calculate the condition number for LU/Cholesky | |||
1401 | if (typ == MatrixType::Upper) | |||
1402 | { | |||
1403 | const double *tmp_data = fortran_vec (); | |||
1404 | octave_idx_type info = 0; | |||
1405 | char norm = '1'; | |||
1406 | char uplo = 'U'; | |||
1407 | char dia = 'N'; | |||
1408 | ||||
1409 | Array<double> z (dim_vector (3 * nc, 1)); | |||
1410 | double *pz = z.fortran_vec (); | |||
1411 | Array<octave_idx_type> iz (dim_vector (nc, 1)); | |||
1412 | octave_idx_type *piz = iz.fortran_vec (); | |||
1413 | ||||
1414 | F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr , rcon, pz, piz, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1415 | F77_CONST_CHAR_ARG2 (&uplo, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr , rcon, pz, piz, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1416 | F77_CONST_CHAR_ARG2 (&dia, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr , rcon, pz, piz, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1417 | nr, tmp_data, nr, rcon,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr , rcon, pz, piz, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1418 | pz, piz, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr , rcon, pz, piz, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1419 | F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr , rcon, pz, piz, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1420 | F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr , rcon, pz, piz, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1421 | F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr , rcon, pz, piz, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0); | |||
1422 | ||||
1423 | if (info != 0) | |||
1424 | rcon = 0.0; | |||
1425 | } | |||
1426 | else if (typ == MatrixType::Permuted_Upper) | |||
1427 | (*current_liboctave_error_handler) | |||
1428 | ("permuted triangular matrix not implemented"); | |||
1429 | else if (typ == MatrixType::Lower) | |||
1430 | { | |||
1431 | const double *tmp_data = fortran_vec (); | |||
1432 | octave_idx_type info = 0; | |||
1433 | char norm = '1'; | |||
1434 | char uplo = 'L'; | |||
1435 | char dia = 'N'; | |||
1436 | ||||
1437 | Array<double> z (dim_vector (3 * nc, 1)); | |||
1438 | double *pz = z.fortran_vec (); | |||
1439 | Array<octave_idx_type> iz (dim_vector (nc, 1)); | |||
1440 | octave_idx_type *piz = iz.fortran_vec (); | |||
1441 | ||||
1442 | F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr , rcon, pz, piz, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1443 | F77_CONST_CHAR_ARG2 (&uplo, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr , rcon, pz, piz, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1444 | F77_CONST_CHAR_ARG2 (&dia, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr , rcon, pz, piz, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1445 | nr, tmp_data, nr, rcon,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr , rcon, pz, piz, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1446 | pz, piz, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr , rcon, pz, piz, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1447 | F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr , rcon, pz, piz, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1448 | F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr , rcon, pz, piz, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1449 | F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr , rcon, pz, piz, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0); | |||
1450 | ||||
1451 | if (info != 0) | |||
1452 | rcon = 0.0; | |||
1453 | } | |||
1454 | else if (typ == MatrixType::Permuted_Lower) | |||
1455 | (*current_liboctave_error_handler) | |||
1456 | ("permuted triangular matrix not implemented"); | |||
1457 | else if (typ == MatrixType::Full || typ == MatrixType::Hermitian) | |||
1458 | { | |||
1459 | double anorm = -1.0; | |||
1460 | ||||
1461 | if (typ == MatrixType::Hermitian) | |||
1462 | { | |||
1463 | octave_idx_type info = 0; | |||
1464 | char job = 'L'; | |||
1465 | ||||
1466 | Matrix atmp = *this; | |||
1467 | double *tmp_data = atmp.fortran_vec (); | |||
1468 | ||||
1469 | anorm = atmp.abs().sum(). | |||
1470 | row(static_cast<octave_idx_type>(0)).max(); | |||
1471 | ||||
1472 | F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dpotrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dpotrf_ (&job, nr, tmp_data, nr, info , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1473 | tmp_data, nr, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dpotrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dpotrf_ (&job, nr, tmp_data, nr, info , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1474 | F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dpotrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dpotrf_ (&job, nr, tmp_data, nr, info , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0); | |||
1475 | ||||
1476 | if (info != 0) | |||
1477 | { | |||
1478 | rcon = 0.0; | |||
1479 | mattype.mark_as_unsymmetric (); | |||
1480 | typ = MatrixType::Full; | |||
1481 | } | |||
1482 | else | |||
1483 | { | |||
1484 | Array<double> z (dim_vector (3 * nc, 1)); | |||
1485 | double *pz = z.fortran_vec (); | |||
1486 | Array<octave_idx_type> iz (dim_vector (nc, 1)); | |||
1487 | octave_idx_type *piz = iz.fortran_vec (); | |||
1488 | ||||
1489 | F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dpocon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dpocon_ (&job, nr, tmp_data, nr, anorm, rcon, pz, piz , info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
1490 | nr, tmp_data, nr, anorm,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dpocon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dpocon_ (&job, nr, tmp_data, nr, anorm, rcon, pz, piz , info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
1491 | rcon, pz, piz, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dpocon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dpocon_ (&job, nr, tmp_data, nr, anorm, rcon, pz, piz , info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
1492 | F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dpocon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dpocon_ (&job, nr, tmp_data, nr, anorm, rcon, pz, piz , info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0); | |||
1493 | ||||
1494 | if (info != 0) | |||
1495 | rcon = 0.0; | |||
1496 | } | |||
1497 | } | |||
1498 | ||||
1499 | if (typ == MatrixType::Full) | |||
1500 | { | |||
1501 | octave_idx_type info = 0; | |||
1502 | ||||
1503 | Matrix atmp = *this; | |||
1504 | double *tmp_data = atmp.fortran_vec (); | |||
1505 | ||||
1506 | Array<octave_idx_type> ipvt (dim_vector (nr, 1)); | |||
1507 | octave_idx_type *pipvt = ipvt.fortran_vec (); | |||
1508 | ||||
1509 | if (anorm < 0.) | |||
1510 | anorm = atmp.abs ().sum (). | |||
1511 | row(static_cast<octave_idx_type>(0)).max (); | |||
1512 | ||||
1513 | Array<double> z (dim_vector (4 * nc, 1)); | |||
1514 | double *pz = z.fortran_vec (); | |||
1515 | Array<octave_idx_type> iz (dim_vector (nc, 1)); | |||
1516 | octave_idx_type *piz = iz.fortran_vec (); | |||
1517 | ||||
1518 | F77_XFCN (dgetrf, DGETRF, (nr, nr, tmp_data, nr, pipvt, info))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgetrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgetrf_ (nr, nr, tmp_data, nr, pipvt, info); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0); | |||
1519 | ||||
1520 | if (info != 0) | |||
1521 | { | |||
1522 | rcon = 0.0; | |||
1523 | mattype.mark_as_rectangular (); | |||
1524 | } | |||
1525 | else | |||
1526 | { | |||
1527 | char job = '1'; | |||
1528 | F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgecon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgecon_ (&job, nc, tmp_data, nr, anorm, rcon, pz, piz , info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
1529 | nc, tmp_data, nr, anorm,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgecon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgecon_ (&job, nc, tmp_data, nr, anorm, rcon, pz, piz , info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
1530 | rcon, pz, piz, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgecon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgecon_ (&job, nc, tmp_data, nr, anorm, rcon, pz, piz , info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
1531 | F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgecon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgecon_ (&job, nc, tmp_data, nr, anorm, rcon, pz, piz , info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0); | |||
1532 | ||||
1533 | if (info != 0) | |||
1534 | rcon = 0.0; | |||
1535 | } | |||
1536 | } | |||
1537 | } | |||
1538 | else | |||
1539 | rcon = 0.0; | |||
1540 | } | |||
1541 | ||||
1542 | return rcon; | |||
| ||||
1543 | } | |||
1544 | ||||
1545 | Matrix | |||
1546 | Matrix::utsolve (MatrixType &mattype, const Matrix& b, octave_idx_type& info, | |||
1547 | double& rcon, solve_singularity_handler sing_handler, | |||
1548 | bool calc_cond, blas_trans_type transt) const | |||
1549 | { | |||
1550 | Matrix retval; | |||
1551 | ||||
1552 | octave_idx_type nr = rows (); | |||
1553 | octave_idx_type nc = cols (); | |||
1554 | ||||
1555 | if (nr != b.rows ()) | |||
1556 | (*current_liboctave_error_handler) | |||
1557 | ("matrix dimension mismatch solution of linear equations"); | |||
1558 | else if (nr == 0 || nc == 0 || b.cols () == 0) | |||
1559 | retval = Matrix (nc, b.cols (), 0.0); | |||
1560 | else | |||
1561 | { | |||
1562 | volatile int typ = mattype.type (); | |||
1563 | ||||
1564 | if (typ == MatrixType::Permuted_Upper || | |||
1565 | typ == MatrixType::Upper) | |||
1566 | { | |||
1567 | octave_idx_type b_nc = b.cols (); | |||
1568 | rcon = 1.; | |||
1569 | info = 0; | |||
1570 | ||||
1571 | if (typ == MatrixType::Permuted_Upper) | |||
1572 | { | |||
1573 | (*current_liboctave_error_handler) | |||
1574 | ("permuted triangular matrix not implemented"); | |||
1575 | } | |||
1576 | else | |||
1577 | { | |||
1578 | const double *tmp_data = fortran_vec (); | |||
1579 | ||||
1580 | if (calc_cond) | |||
1581 | { | |||
1582 | char norm = '1'; | |||
1583 | char uplo = 'U'; | |||
1584 | char dia = 'N'; | |||
1585 | ||||
1586 | Array<double> z (dim_vector (3 * nc, 1)); | |||
1587 | double *pz = z.fortran_vec (); | |||
1588 | Array<octave_idx_type> iz (dim_vector (nc, 1)); | |||
1589 | octave_idx_type *piz = iz.fortran_vec (); | |||
1590 | ||||
1591 | F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr , rcon, pz, piz, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1592 | F77_CONST_CHAR_ARG2 (&uplo, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr , rcon, pz, piz, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1593 | F77_CONST_CHAR_ARG2 (&dia, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr , rcon, pz, piz, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1594 | nr, tmp_data, nr, rcon,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr , rcon, pz, piz, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1595 | pz, piz, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr , rcon, pz, piz, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1596 | F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr , rcon, pz, piz, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1597 | F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr , rcon, pz, piz, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1598 | F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr , rcon, pz, piz, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0); | |||
1599 | ||||
1600 | if (info != 0) | |||
1601 | info = -2; | |||
1602 | ||||
1603 | volatile double rcond_plus_one = rcon + 1.0; | |||
1604 | ||||
1605 | if (rcond_plus_one == 1.0 || xisnan (rcon)) | |||
1606 | { | |||
1607 | info = -2; | |||
1608 | ||||
1609 | if (sing_handler) | |||
1610 | sing_handler (rcon); | |||
1611 | else | |||
1612 | (*current_liboctave_error_handler) | |||
1613 | ("matrix singular to machine precision, rcond = %g", | |||
1614 | rcon); | |||
1615 | } | |||
1616 | } | |||
1617 | ||||
1618 | if (info == 0) | |||
1619 | { | |||
1620 | retval = b; | |||
1621 | double *result = retval.fortran_vec (); | |||
1622 | ||||
1623 | char uplo = 'U'; | |||
1624 | char trans = get_blas_char (transt); | |||
1625 | char dia = 'N'; | |||
1626 | ||||
1627 | F77_XFCN (dtrtrs, DTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrtrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrtrs_ (&uplo, &trans, &dia, nr, b_nc, tmp_data , nr, result, nr, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1628 | F77_CONST_CHAR_ARG2 (&trans, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrtrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrtrs_ (&uplo, &trans, &dia, nr, b_nc, tmp_data , nr, result, nr, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1629 | F77_CONST_CHAR_ARG2 (&dia, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrtrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrtrs_ (&uplo, &trans, &dia, nr, b_nc, tmp_data , nr, result, nr, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1630 | nr, b_nc, tmp_data, nr,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrtrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrtrs_ (&uplo, &trans, &dia, nr, b_nc, tmp_data , nr, result, nr, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1631 | result, nr, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrtrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrtrs_ (&uplo, &trans, &dia, nr, b_nc, tmp_data , nr, result, nr, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1632 | F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrtrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrtrs_ (&uplo, &trans, &dia, nr, b_nc, tmp_data , nr, result, nr, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1633 | F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrtrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrtrs_ (&uplo, &trans, &dia, nr, b_nc, tmp_data , nr, result, nr, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1634 | F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrtrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrtrs_ (&uplo, &trans, &dia, nr, b_nc, tmp_data , nr, result, nr, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0); | |||
1635 | } | |||
1636 | } | |||
1637 | } | |||
1638 | else | |||
1639 | (*current_liboctave_error_handler) ("incorrect matrix type"); | |||
1640 | } | |||
1641 | ||||
1642 | return retval; | |||
1643 | } | |||
1644 | ||||
1645 | Matrix | |||
1646 | Matrix::ltsolve (MatrixType &mattype, const Matrix& b, octave_idx_type& info, | |||
1647 | double& rcon, solve_singularity_handler sing_handler, | |||
1648 | bool calc_cond, blas_trans_type transt) const | |||
1649 | { | |||
1650 | Matrix retval; | |||
1651 | ||||
1652 | octave_idx_type nr = rows (); | |||
1653 | octave_idx_type nc = cols (); | |||
1654 | ||||
1655 | if (nr != b.rows ()) | |||
1656 | (*current_liboctave_error_handler) | |||
1657 | ("matrix dimension mismatch solution of linear equations"); | |||
1658 | else if (nr == 0 || nc == 0 || b.cols () == 0) | |||
1659 | retval = Matrix (nc, b.cols (), 0.0); | |||
1660 | else | |||
1661 | { | |||
1662 | volatile int typ = mattype.type (); | |||
1663 | ||||
1664 | if (typ == MatrixType::Permuted_Lower || | |||
1665 | typ == MatrixType::Lower) | |||
1666 | { | |||
1667 | octave_idx_type b_nc = b.cols (); | |||
1668 | rcon = 1.; | |||
1669 | info = 0; | |||
1670 | ||||
1671 | if (typ == MatrixType::Permuted_Lower) | |||
1672 | { | |||
1673 | (*current_liboctave_error_handler) | |||
1674 | ("permuted triangular matrix not implemented"); | |||
1675 | } | |||
1676 | else | |||
1677 | { | |||
1678 | const double *tmp_data = fortran_vec (); | |||
1679 | ||||
1680 | if (calc_cond) | |||
1681 | { | |||
1682 | char norm = '1'; | |||
1683 | char uplo = 'L'; | |||
1684 | char dia = 'N'; | |||
1685 | ||||
1686 | Array<double> z (dim_vector (3 * nc, 1)); | |||
1687 | double *pz = z.fortran_vec (); | |||
1688 | Array<octave_idx_type> iz (dim_vector (nc, 1)); | |||
1689 | octave_idx_type *piz = iz.fortran_vec (); | |||
1690 | ||||
1691 | F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr , rcon, pz, piz, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1692 | F77_CONST_CHAR_ARG2 (&uplo, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr , rcon, pz, piz, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1693 | F77_CONST_CHAR_ARG2 (&dia, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr , rcon, pz, piz, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1694 | nr, tmp_data, nr, rcon,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr , rcon, pz, piz, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1695 | pz, piz, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr , rcon, pz, piz, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1696 | F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr , rcon, pz, piz, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1697 | F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr , rcon, pz, piz, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1698 | F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrcon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrcon_ (&norm, &uplo, &dia, nr, tmp_data, nr , rcon, pz, piz, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0); | |||
1699 | ||||
1700 | if (info != 0) | |||
1701 | info = -2; | |||
1702 | ||||
1703 | volatile double rcond_plus_one = rcon + 1.0; | |||
1704 | ||||
1705 | if (rcond_plus_one == 1.0 || xisnan (rcon)) | |||
1706 | { | |||
1707 | info = -2; | |||
1708 | ||||
1709 | if (sing_handler) | |||
1710 | sing_handler (rcon); | |||
1711 | else | |||
1712 | (*current_liboctave_error_handler) | |||
1713 | ("matrix singular to machine precision, rcond = %g", | |||
1714 | rcon); | |||
1715 | } | |||
1716 | } | |||
1717 | ||||
1718 | if (info == 0) | |||
1719 | { | |||
1720 | retval = b; | |||
1721 | double *result = retval.fortran_vec (); | |||
1722 | ||||
1723 | char uplo = 'L'; | |||
1724 | char trans = get_blas_char (transt); | |||
1725 | char dia = 'N'; | |||
1726 | ||||
1727 | F77_XFCN (dtrtrs, DTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrtrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrtrs_ (&uplo, &trans, &dia, nr, b_nc, tmp_data , nr, result, nr, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1728 | F77_CONST_CHAR_ARG2 (&trans, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrtrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrtrs_ (&uplo, &trans, &dia, nr, b_nc, tmp_data , nr, result, nr, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1729 | F77_CONST_CHAR_ARG2 (&dia, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrtrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrtrs_ (&uplo, &trans, &dia, nr, b_nc, tmp_data , nr, result, nr, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1730 | nr, b_nc, tmp_data, nr,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrtrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrtrs_ (&uplo, &trans, &dia, nr, b_nc, tmp_data , nr, result, nr, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1731 | result, nr, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrtrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrtrs_ (&uplo, &trans, &dia, nr, b_nc, tmp_data , nr, result, nr, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1732 | F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrtrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrtrs_ (&uplo, &trans, &dia, nr, b_nc, tmp_data , nr, result, nr, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1733 | F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrtrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrtrs_ (&uplo, &trans, &dia, nr, b_nc, tmp_data , nr, result, nr, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1734 | F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrtrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrtrs_ (&uplo, &trans, &dia, nr, b_nc, tmp_data , nr, result, nr, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0); | |||
1735 | } | |||
1736 | } | |||
1737 | } | |||
1738 | else | |||
1739 | (*current_liboctave_error_handler) ("incorrect matrix type"); | |||
1740 | } | |||
1741 | ||||
1742 | return retval; | |||
1743 | } | |||
1744 | ||||
1745 | Matrix | |||
1746 | Matrix::fsolve (MatrixType &mattype, const Matrix& b, octave_idx_type& info, | |||
1747 | double& rcon, solve_singularity_handler sing_handler, | |||
1748 | bool calc_cond) const | |||
1749 | { | |||
1750 | Matrix retval; | |||
1751 | ||||
1752 | octave_idx_type nr = rows (); | |||
1753 | octave_idx_type nc = cols (); | |||
1754 | ||||
1755 | if (nr != nc || nr != b.rows ()) | |||
1756 | (*current_liboctave_error_handler) | |||
1757 | ("matrix dimension mismatch solution of linear equations"); | |||
1758 | else if (nr == 0 || b.cols () == 0) | |||
1759 | retval = Matrix (nc, b.cols (), 0.0); | |||
1760 | else | |||
1761 | { | |||
1762 | volatile int typ = mattype.type (); | |||
1763 | ||||
1764 | // Calculate the norm of the matrix, for later use. | |||
1765 | double anorm = -1.; | |||
1766 | ||||
1767 | if (typ == MatrixType::Hermitian) | |||
1768 | { | |||
1769 | info = 0; | |||
1770 | char job = 'L'; | |||
1771 | ||||
1772 | Matrix atmp = *this; | |||
1773 | double *tmp_data = atmp.fortran_vec (); | |||
1774 | ||||
1775 | anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max(); | |||
1776 | ||||
1777 | F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dpotrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dpotrf_ (&job, nr, tmp_data, nr, info , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1778 | tmp_data, nr, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dpotrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dpotrf_ (&job, nr, tmp_data, nr, info , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
1779 | F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dpotrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dpotrf_ (&job, nr, tmp_data, nr, info , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0); | |||
1780 | ||||
1781 | // Throw-away extra info LAPACK gives so as to not change output. | |||
1782 | rcon = 0.0; | |||
1783 | if (info != 0) | |||
1784 | { | |||
1785 | info = -2; | |||
1786 | ||||
1787 | mattype.mark_as_unsymmetric (); | |||
1788 | typ = MatrixType::Full; | |||
1789 | } | |||
1790 | else | |||
1791 | { | |||
1792 | if (calc_cond) | |||
1793 | { | |||
1794 | Array<double> z (dim_vector (3 * nc, 1)); | |||
1795 | double *pz = z.fortran_vec (); | |||
1796 | Array<octave_idx_type> iz (dim_vector (nc, 1)); | |||
1797 | octave_idx_type *piz = iz.fortran_vec (); | |||
1798 | ||||
1799 | F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dpocon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dpocon_ (&job, nr, tmp_data, nr, anorm, rcon, pz, piz , info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
1800 | nr, tmp_data, nr, anorm,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dpocon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dpocon_ (&job, nr, tmp_data, nr, anorm, rcon, pz, piz , info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
1801 | rcon, pz, piz, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dpocon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dpocon_ (&job, nr, tmp_data, nr, anorm, rcon, pz, piz , info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
1802 | F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dpocon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dpocon_ (&job, nr, tmp_data, nr, anorm, rcon, pz, piz , info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0); | |||
1803 | ||||
1804 | if (info != 0) | |||
1805 | info = -2; | |||
1806 | ||||
1807 | volatile double rcond_plus_one = rcon + 1.0; | |||
1808 | ||||
1809 | if (rcond_plus_one == 1.0 || xisnan (rcon)) | |||
1810 | { | |||
1811 | info = -2; | |||
1812 | ||||
1813 | if (sing_handler) | |||
1814 | sing_handler (rcon); | |||
1815 | else | |||
1816 | (*current_liboctave_error_handler) | |||
1817 | ("matrix singular to machine precision, rcond = %g", | |||
1818 | rcon); | |||
1819 | } | |||
1820 | } | |||
1821 | ||||
1822 | if (info == 0) | |||
1823 | { | |||
1824 | retval = b; | |||
1825 | double *result = retval.fortran_vec (); | |||
1826 | ||||
1827 | octave_idx_type b_nc = b.cols (); | |||
1828 | ||||
1829 | F77_XFCN (dpotrs, DPOTRS, (F77_CONST_CHAR_ARG2 (&job, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dpotrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dpotrs_ (&job, nr, b_nc, tmp_data, nr, result, b.rows (), info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
1830 | nr, b_nc, tmp_data, nr,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dpotrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dpotrs_ (&job, nr, b_nc, tmp_data, nr, result, b.rows (), info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
1831 | result, b.rows (), infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dpotrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dpotrs_ (&job, nr, b_nc, tmp_data, nr, result, b.rows (), info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
1832 | F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dpotrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dpotrs_ (&job, nr, b_nc, tmp_data, nr, result, b.rows (), info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0); | |||
1833 | } | |||
1834 | else | |||
1835 | { | |||
1836 | mattype.mark_as_unsymmetric (); | |||
1837 | typ = MatrixType::Full; | |||
1838 | } | |||
1839 | } | |||
1840 | } | |||
1841 | ||||
1842 | if (typ == MatrixType::Full) | |||
1843 | { | |||
1844 | info = 0; | |||
1845 | ||||
1846 | Array<octave_idx_type> ipvt (dim_vector (nr, 1)); | |||
1847 | octave_idx_type *pipvt = ipvt.fortran_vec (); | |||
1848 | ||||
1849 | Matrix atmp = *this; | |||
1850 | double *tmp_data = atmp.fortran_vec (); | |||
1851 | ||||
1852 | if (anorm < 0.) | |||
1853 | anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max(); | |||
1854 | ||||
1855 | Array<double> z (dim_vector (4 * nc, 1)); | |||
1856 | double *pz = z.fortran_vec (); | |||
1857 | Array<octave_idx_type> iz (dim_vector (nc, 1)); | |||
1858 | octave_idx_type *piz = iz.fortran_vec (); | |||
1859 | ||||
1860 | F77_XFCN (dgetrf, DGETRF, (nr, nr, tmp_data, nr, pipvt, info))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgetrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgetrf_ (nr, nr, tmp_data, nr, pipvt, info); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0); | |||
1861 | ||||
1862 | // Throw-away extra info LAPACK gives so as to not change output. | |||
1863 | rcon = 0.0; | |||
1864 | if (info != 0) | |||
1865 | { | |||
1866 | info = -2; | |||
1867 | ||||
1868 | if (sing_handler) | |||
1869 | sing_handler (rcon); | |||
1870 | else | |||
1871 | (*current_liboctave_error_handler) | |||
1872 | ("matrix singular to machine precision"); | |||
1873 | ||||
1874 | mattype.mark_as_rectangular (); | |||
1875 | } | |||
1876 | else | |||
1877 | { | |||
1878 | if (calc_cond) | |||
1879 | { | |||
1880 | // Now calculate the condition number for | |||
1881 | // non-singular matrix. | |||
1882 | char job = '1'; | |||
1883 | F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgecon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgecon_ (&job, nc, tmp_data, nr, anorm, rcon, pz, piz , info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
1884 | nc, tmp_data, nr, anorm,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgecon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgecon_ (&job, nc, tmp_data, nr, anorm, rcon, pz, piz , info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
1885 | rcon, pz, piz, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgecon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgecon_ (&job, nc, tmp_data, nr, anorm, rcon, pz, piz , info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
1886 | F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgecon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgecon_ (&job, nc, tmp_data, nr, anorm, rcon, pz, piz , info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0); | |||
1887 | ||||
1888 | if (info != 0) | |||
1889 | info = -2; | |||
1890 | ||||
1891 | volatile double rcond_plus_one = rcon + 1.0; | |||
1892 | ||||
1893 | if (rcond_plus_one == 1.0 || xisnan (rcon)) | |||
1894 | { | |||
1895 | info = -2; | |||
1896 | ||||
1897 | if (sing_handler) | |||
1898 | sing_handler (rcon); | |||
1899 | else | |||
1900 | (*current_liboctave_error_handler) | |||
1901 | ("matrix singular to machine precision, rcond = %g", | |||
1902 | rcon); | |||
1903 | } | |||
1904 | } | |||
1905 | ||||
1906 | if (info == 0) | |||
1907 | { | |||
1908 | retval = b; | |||
1909 | double *result = retval.fortran_vec (); | |||
1910 | ||||
1911 | octave_idx_type b_nc = b.cols (); | |||
1912 | ||||
1913 | char job = 'N'; | |||
1914 | F77_XFCN (dgetrs, DGETRS, (F77_CONST_CHAR_ARG2 (&job, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgetrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgetrs_ (&job, nr, b_nc, tmp_data, nr, pipvt, result, b.rows (), info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
1915 | nr, b_nc, tmp_data, nr,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgetrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgetrs_ (&job, nr, b_nc, tmp_data, nr, pipvt, result, b.rows (), info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
1916 | pipvt, result, b.rows (), infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgetrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgetrs_ (&job, nr, b_nc, tmp_data, nr, pipvt, result, b.rows (), info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
1917 | F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgetrs_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgetrs_ (&job, nr, b_nc, tmp_data, nr, pipvt, result, b.rows (), info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0); | |||
1918 | } | |||
1919 | else | |||
1920 | mattype.mark_as_rectangular (); | |||
1921 | } | |||
1922 | } | |||
1923 | else if (typ != MatrixType::Hermitian) | |||
1924 | (*current_liboctave_error_handler) ("incorrect matrix type"); | |||
1925 | } | |||
1926 | ||||
1927 | return retval; | |||
1928 | } | |||
1929 | ||||
1930 | Matrix | |||
1931 | Matrix::solve (MatrixType &typ, const Matrix& b) const | |||
1932 | { | |||
1933 | octave_idx_type info; | |||
1934 | double rcon; | |||
1935 | return solve (typ, b, info, rcon, 0); | |||
1936 | } | |||
1937 | ||||
1938 | Matrix | |||
1939 | Matrix::solve (MatrixType &typ, const Matrix& b, octave_idx_type& info) const | |||
1940 | { | |||
1941 | double rcon; | |||
1942 | return solve (typ, b, info, rcon, 0); | |||
1943 | } | |||
1944 | ||||
1945 | Matrix | |||
1946 | Matrix::solve (MatrixType &typ, const Matrix& b, octave_idx_type& info, | |||
1947 | double& rcon) const | |||
1948 | { | |||
1949 | return solve (typ, b, info, rcon, 0); | |||
1950 | } | |||
1951 | ||||
1952 | Matrix | |||
1953 | Matrix::solve (MatrixType &mattype, const Matrix& b, octave_idx_type& info, | |||
1954 | double& rcon, solve_singularity_handler sing_handler, | |||
1955 | bool singular_fallback, blas_trans_type transt) const | |||
1956 | { | |||
1957 | Matrix retval; | |||
1958 | int typ = mattype.type (); | |||
1959 | ||||
1960 | if (typ == MatrixType::Unknown) | |||
1961 | typ = mattype.type (*this); | |||
1962 | ||||
1963 | // Only calculate the condition number for LU/Cholesky | |||
1964 | if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper) | |||
1965 | retval = utsolve (mattype, b, info, rcon, sing_handler, false, transt); | |||
1966 | else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower) | |||
1967 | retval = ltsolve (mattype, b, info, rcon, sing_handler, false, transt); | |||
1968 | else if (transt == blas_trans || transt == blas_conj_trans) | |||
1969 | return transpose ().solve (mattype, b, info, rcon, sing_handler, | |||
1970 | singular_fallback); | |||
1971 | else if (typ == MatrixType::Full || typ == MatrixType::Hermitian) | |||
1972 | retval = fsolve (mattype, b, info, rcon, sing_handler, true); | |||
1973 | else if (typ != MatrixType::Rectangular) | |||
1974 | { | |||
1975 | (*current_liboctave_error_handler) ("unknown matrix type"); | |||
1976 | return Matrix (); | |||
1977 | } | |||
1978 | ||||
1979 | // Rectangular or one of the above solvers flags a singular matrix | |||
1980 | if (singular_fallback && mattype.type () == MatrixType::Rectangular) | |||
1981 | { | |||
1982 | octave_idx_type rank; | |||
1983 | retval = lssolve (b, info, rank, rcon); | |||
1984 | } | |||
1985 | ||||
1986 | return retval; | |||
1987 | } | |||
1988 | ||||
1989 | ComplexMatrix | |||
1990 | Matrix::solve (MatrixType &typ, const ComplexMatrix& b) const | |||
1991 | { | |||
1992 | octave_idx_type info; | |||
1993 | double rcon; | |||
1994 | return solve (typ, b, info, rcon, 0); | |||
1995 | } | |||
1996 | ||||
1997 | ComplexMatrix | |||
1998 | Matrix::solve (MatrixType &typ, const ComplexMatrix& b, | |||
1999 | octave_idx_type& info) const | |||
2000 | { | |||
2001 | double rcon; | |||
2002 | return solve (typ, b, info, rcon, 0); | |||
2003 | } | |||
2004 | ||||
2005 | ComplexMatrix | |||
2006 | Matrix::solve (MatrixType &typ, const ComplexMatrix& b, octave_idx_type& info, | |||
2007 | double& rcon) const | |||
2008 | { | |||
2009 | return solve (typ, b, info, rcon, 0); | |||
2010 | } | |||
2011 | ||||
2012 | static Matrix | |||
2013 | stack_complex_matrix (const ComplexMatrix& cm) | |||
2014 | { | |||
2015 | octave_idx_type m = cm.rows (), n = cm.cols (), nel = m*n; | |||
2016 | Matrix retval (m, 2*n); | |||
2017 | const Complex *cmd = cm.data (); | |||
2018 | double *rd = retval.fortran_vec (); | |||
2019 | for (octave_idx_type i = 0; i < nel; i++) | |||
2020 | { | |||
2021 | rd[i] = std::real (cmd[i]); | |||
2022 | rd[nel+i] = std::imag (cmd[i]); | |||
2023 | } | |||
2024 | return retval; | |||
2025 | } | |||
2026 | ||||
2027 | static ComplexMatrix | |||
2028 | unstack_complex_matrix (const Matrix& sm) | |||
2029 | { | |||
2030 | octave_idx_type m = sm.rows (), n = sm.cols () / 2, nel = m*n; | |||
2031 | ComplexMatrix retval (m, n); | |||
2032 | const double *smd = sm.data (); | |||
2033 | Complex *rd = retval.fortran_vec (); | |||
2034 | for (octave_idx_type i = 0; i < nel; i++) | |||
2035 | rd[i] = Complex (smd[i], smd[nel+i]); | |||
2036 | return retval; | |||
2037 | } | |||
2038 | ||||
2039 | ComplexMatrix | |||
2040 | Matrix::solve (MatrixType &typ, const ComplexMatrix& b, octave_idx_type& info, | |||
2041 | double& rcon, solve_singularity_handler sing_handler, | |||
2042 | bool singular_fallback, blas_trans_type transt) const | |||
2043 | { | |||
2044 | Matrix tmp = stack_complex_matrix (b); | |||
2045 | tmp = solve (typ, tmp, info, rcon, sing_handler, singular_fallback, transt); | |||
2046 | return unstack_complex_matrix (tmp); | |||
2047 | } | |||
2048 | ||||
2049 | ColumnVector | |||
2050 | Matrix::solve (MatrixType &typ, const ColumnVector& b) const | |||
2051 | { | |||
2052 | octave_idx_type info; double rcon; | |||
2053 | return solve (typ, b, info, rcon); | |||
2054 | } | |||
2055 | ||||
2056 | ColumnVector | |||
2057 | Matrix::solve (MatrixType &typ, const ColumnVector& b, | |||
2058 | octave_idx_type& info) const | |||
2059 | { | |||
2060 | double rcon; | |||
2061 | return solve (typ, b, info, rcon); | |||
2062 | } | |||
2063 | ||||
2064 | ColumnVector | |||
2065 | Matrix::solve (MatrixType &typ, const ColumnVector& b, octave_idx_type& info, | |||
2066 | double& rcon) const | |||
2067 | { | |||
2068 | return solve (typ, b, info, rcon, 0); | |||
2069 | } | |||
2070 | ||||
2071 | ColumnVector | |||
2072 | Matrix::solve (MatrixType &typ, const ColumnVector& b, octave_idx_type& info, | |||
2073 | double& rcon, solve_singularity_handler sing_handler, | |||
2074 | blas_trans_type transt) const | |||
2075 | { | |||
2076 | Matrix tmp (b); | |||
2077 | tmp = solve (typ, tmp, info, rcon, sing_handler, true, transt); | |||
2078 | return tmp.column (static_cast<octave_idx_type> (0)); | |||
2079 | } | |||
2080 | ||||
2081 | ComplexColumnVector | |||
2082 | Matrix::solve (MatrixType &typ, const ComplexColumnVector& b) const | |||
2083 | { | |||
2084 | ComplexMatrix tmp (*this); | |||
2085 | return tmp.solve (typ, b); | |||
2086 | } | |||
2087 | ||||
2088 | ComplexColumnVector | |||
2089 | Matrix::solve (MatrixType &typ, const ComplexColumnVector& b, | |||
2090 | octave_idx_type& info) const | |||
2091 | { | |||
2092 | ComplexMatrix tmp (*this); | |||
2093 | return tmp.solve (typ, b, info); | |||
2094 | } | |||
2095 | ||||
2096 | ComplexColumnVector | |||
2097 | Matrix::solve (MatrixType &typ, const ComplexColumnVector& b, | |||
2098 | octave_idx_type& info, double& rcon) const | |||
2099 | { | |||
2100 | ComplexMatrix tmp (*this); | |||
2101 | return tmp.solve (typ, b, info, rcon); | |||
2102 | } | |||
2103 | ||||
2104 | ComplexColumnVector | |||
2105 | Matrix::solve (MatrixType &typ, const ComplexColumnVector& b, | |||
2106 | octave_idx_type& info, double& rcon, | |||
2107 | solve_singularity_handler sing_handler, | |||
2108 | blas_trans_type transt) const | |||
2109 | { | |||
2110 | ComplexMatrix tmp (*this); | |||
2111 | return tmp.solve (typ, b, info, rcon, sing_handler, transt); | |||
2112 | } | |||
2113 | ||||
2114 | Matrix | |||
2115 | Matrix::solve (const Matrix& b) const | |||
2116 | { | |||
2117 | octave_idx_type info; | |||
2118 | double rcon; | |||
2119 | return solve (b, info, rcon, 0); | |||
2120 | } | |||
2121 | ||||
2122 | Matrix | |||
2123 | Matrix::solve (const Matrix& b, octave_idx_type& info) const | |||
2124 | { | |||
2125 | double rcon; | |||
2126 | return solve (b, info, rcon, 0); | |||
2127 | } | |||
2128 | ||||
2129 | Matrix | |||
2130 | Matrix::solve (const Matrix& b, octave_idx_type& info, double& rcon) const | |||
2131 | { | |||
2132 | return solve (b, info, rcon, 0); | |||
2133 | } | |||
2134 | ||||
2135 | Matrix | |||
2136 | Matrix::solve (const Matrix& b, octave_idx_type& info, | |||
2137 | double& rcon, solve_singularity_handler sing_handler, | |||
2138 | blas_trans_type transt) const | |||
2139 | { | |||
2140 | MatrixType mattype (*this); | |||
2141 | return solve (mattype, b, info, rcon, sing_handler, true, transt); | |||
2142 | } | |||
2143 | ||||
2144 | ComplexMatrix | |||
2145 | Matrix::solve (const ComplexMatrix& b) const | |||
2146 | { | |||
2147 | ComplexMatrix tmp (*this); | |||
2148 | return tmp.solve (b); | |||
2149 | } | |||
2150 | ||||
2151 | ComplexMatrix | |||
2152 | Matrix::solve (const ComplexMatrix& b, octave_idx_type& info) const | |||
2153 | { | |||
2154 | ComplexMatrix tmp (*this); | |||
2155 | return tmp.solve (b, info); | |||
2156 | } | |||
2157 | ||||
2158 | ComplexMatrix | |||
2159 | Matrix::solve (const ComplexMatrix& b, octave_idx_type& info, | |||
2160 | double& rcon) const | |||
2161 | { | |||
2162 | ComplexMatrix tmp (*this); | |||
2163 | return tmp.solve (b, info, rcon); | |||
2164 | } | |||
2165 | ||||
2166 | ComplexMatrix | |||
2167 | Matrix::solve (const ComplexMatrix& b, octave_idx_type& info, double& rcon, | |||
2168 | solve_singularity_handler sing_handler, | |||
2169 | blas_trans_type transt) const | |||
2170 | { | |||
2171 | ComplexMatrix tmp (*this); | |||
2172 | return tmp.solve (b, info, rcon, sing_handler, transt); | |||
2173 | } | |||
2174 | ||||
2175 | ColumnVector | |||
2176 | Matrix::solve (const ColumnVector& b) const | |||
2177 | { | |||
2178 | octave_idx_type info; double rcon; | |||
2179 | return solve (b, info, rcon); | |||
2180 | } | |||
2181 | ||||
2182 | ColumnVector | |||
2183 | Matrix::solve (const ColumnVector& b, octave_idx_type& info) const | |||
2184 | { | |||
2185 | double rcon; | |||
2186 | return solve (b, info, rcon); | |||
2187 | } | |||
2188 | ||||
2189 | ColumnVector | |||
2190 | Matrix::solve (const ColumnVector& b, octave_idx_type& info, double& rcon) const | |||
2191 | { | |||
2192 | return solve (b, info, rcon, 0); | |||
2193 | } | |||
2194 | ||||
2195 | ColumnVector | |||
2196 | Matrix::solve (const ColumnVector& b, octave_idx_type& info, double& rcon, | |||
2197 | solve_singularity_handler sing_handler, | |||
2198 | blas_trans_type transt) const | |||
2199 | { | |||
2200 | MatrixType mattype (*this); | |||
2201 | return solve (mattype, b, info, rcon, sing_handler, transt); | |||
2202 | } | |||
2203 | ||||
2204 | ComplexColumnVector | |||
2205 | Matrix::solve (const ComplexColumnVector& b) const | |||
2206 | { | |||
2207 | ComplexMatrix tmp (*this); | |||
2208 | return tmp.solve (b); | |||
2209 | } | |||
2210 | ||||
2211 | ComplexColumnVector | |||
2212 | Matrix::solve (const ComplexColumnVector& b, octave_idx_type& info) const | |||
2213 | { | |||
2214 | ComplexMatrix tmp (*this); | |||
2215 | return tmp.solve (b, info); | |||
2216 | } | |||
2217 | ||||
2218 | ComplexColumnVector | |||
2219 | Matrix::solve (const ComplexColumnVector& b, octave_idx_type& info, | |||
2220 | double& rcon) const | |||
2221 | { | |||
2222 | ComplexMatrix tmp (*this); | |||
2223 | return tmp.solve (b, info, rcon); | |||
2224 | } | |||
2225 | ||||
2226 | ComplexColumnVector | |||
2227 | Matrix::solve (const ComplexColumnVector& b, octave_idx_type& info, | |||
2228 | double& rcon, | |||
2229 | solve_singularity_handler sing_handler, | |||
2230 | blas_trans_type transt) const | |||
2231 | { | |||
2232 | ComplexMatrix tmp (*this); | |||
2233 | return tmp.solve (b, info, rcon, sing_handler, transt); | |||
2234 | } | |||
2235 | ||||
2236 | Matrix | |||
2237 | Matrix::lssolve (const Matrix& b) const | |||
2238 | { | |||
2239 | octave_idx_type info; | |||
2240 | octave_idx_type rank; | |||
2241 | double rcon; | |||
2242 | return lssolve (b, info, rank, rcon); | |||
2243 | } | |||
2244 | ||||
2245 | Matrix | |||
2246 | Matrix::lssolve (const Matrix& b, octave_idx_type& info) const | |||
2247 | { | |||
2248 | octave_idx_type rank; | |||
2249 | double rcon; | |||
2250 | return lssolve (b, info, rank, rcon); | |||
2251 | } | |||
2252 | ||||
2253 | Matrix | |||
2254 | Matrix::lssolve (const Matrix& b, octave_idx_type& info, | |||
2255 | octave_idx_type& rank) const | |||
2256 | { | |||
2257 | double rcon; | |||
2258 | return lssolve (b, info, rank, rcon); | |||
2259 | } | |||
2260 | ||||
2261 | Matrix | |||
2262 | Matrix::lssolve (const Matrix& b, octave_idx_type& info, | |||
2263 | octave_idx_type& rank, double &rcon) const | |||
2264 | { | |||
2265 | Matrix retval; | |||
2266 | ||||
2267 | octave_idx_type nrhs = b.cols (); | |||
2268 | ||||
2269 | octave_idx_type m = rows (); | |||
2270 | octave_idx_type n = cols (); | |||
2271 | ||||
2272 | if (m != b.rows ()) | |||
2273 | (*current_liboctave_error_handler) | |||
2274 | ("matrix dimension mismatch solution of linear equations"); | |||
2275 | else if (m == 0 || n == 0 || b.cols () == 0) | |||
2276 | retval = Matrix (n, b.cols (), 0.0); | |||
2277 | else | |||
2278 | { | |||
2279 | volatile octave_idx_type minmn = (m < n ? m : n); | |||
2280 | octave_idx_type maxmn = m > n ? m : n; | |||
2281 | rcon = -1.0; | |||
2282 | if (m != n) | |||
2283 | { | |||
2284 | retval = Matrix (maxmn, nrhs, 0.0); | |||
2285 | ||||
2286 | for (octave_idx_type j = 0; j < nrhs; j++) | |||
2287 | for (octave_idx_type i = 0; i < m; i++) | |||
2288 | retval.elem (i, j) = b.elem (i, j); | |||
2289 | } | |||
2290 | else | |||
2291 | retval = b; | |||
2292 | ||||
2293 | Matrix atmp = *this; | |||
2294 | double *tmp_data = atmp.fortran_vec (); | |||
2295 | ||||
2296 | double *pretval = retval.fortran_vec (); | |||
2297 | Array<double> s (dim_vector (minmn, 1)); | |||
2298 | double *ps = s.fortran_vec (); | |||
2299 | ||||
2300 | // Ask DGELSD what the dimension of WORK should be. | |||
2301 | octave_idx_type lwork = -1; | |||
2302 | ||||
2303 | Array<double> work (dim_vector (1, 1)); | |||
2304 | ||||
2305 | octave_idx_type smlsiz; | |||
2306 | F77_FUNC (xilaenv, XILAENV)xilaenv_ (9, F77_CONST_CHAR_ARG2 ("DGELSD", 6)"DGELSD", | |||
2307 | F77_CONST_CHAR_ARG2 (" ", 1)" ", | |||
2308 | 0, 0, 0, 0, smlsiz | |||
2309 | F77_CHAR_ARG_LEN (6), 6 | |||
2310 | F77_CHAR_ARG_LEN (1), 1); | |||
2311 | ||||
2312 | octave_idx_type mnthr; | |||
2313 | F77_FUNC (xilaenv, XILAENV)xilaenv_ (6, F77_CONST_CHAR_ARG2 ("DGELSD", 6)"DGELSD", | |||
2314 | F77_CONST_CHAR_ARG2 (" ", 1)" ", | |||
2315 | m, n, nrhs, -1, mnthr | |||
2316 | F77_CHAR_ARG_LEN (6), 6 | |||
2317 | F77_CHAR_ARG_LEN (1), 1); | |||
2318 | ||||
2319 | // We compute the size of iwork because DGELSD in older versions | |||
2320 | // of LAPACK does not return it on a query call. | |||
2321 | double dminmn = static_cast<double> (minmn); | |||
2322 | double dsmlsizp1 = static_cast<double> (smlsiz+1); | |||
2323 | #if defined (HAVE_LOG21) | |||
2324 | double tmp = log2 (dminmn / dsmlsizp1); | |||
2325 | #else | |||
2326 | double tmp = log (dminmn / dsmlsizp1) / log (2.0); | |||
2327 | #endif | |||
2328 | octave_idx_type nlvl = static_cast<octave_idx_type> (tmp) + 1; | |||
2329 | if (nlvl < 0) | |||
2330 | nlvl = 0; | |||
2331 | ||||
2332 | octave_idx_type liwork = 3 * minmn * nlvl + 11 * minmn; | |||
2333 | if (liwork < 1) | |||
2334 | liwork = 1; | |||
2335 | Array<octave_idx_type> iwork (dim_vector (liwork, 1)); | |||
2336 | octave_idx_type* piwork = iwork.fortran_vec (); | |||
2337 | ||||
2338 | F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval, maxmn,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgelsd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgelsd_ (m, n, nrhs, tmp_data, m, pretval, maxmn, ps, rcon , rank, work.fortran_vec (), lwork, piwork, info); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
2339 | ps, rcon, rank, work.fortran_vec (),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgelsd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgelsd_ (m, n, nrhs, tmp_data, m, pretval, maxmn, ps, rcon , rank, work.fortran_vec (), lwork, piwork, info); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
2340 | lwork, piwork, info))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgelsd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgelsd_ (m, n, nrhs, tmp_data, m, pretval, maxmn, ps, rcon , rank, work.fortran_vec (), lwork, piwork, info); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0); | |||
2341 | ||||
2342 | // The workspace query is broken in at least LAPACK 3.0.0 | |||
2343 | // through 3.1.1 when n >= mnthr. The obtuse formula below | |||
2344 | // should provide sufficient workspace for DGELSD to operate | |||
2345 | // efficiently. | |||
2346 | if (n > m && n >= mnthr) | |||
2347 | { | |||
2348 | const octave_idx_type wlalsd | |||
2349 | = 9*m + 2*m*smlsiz + 8*m*nlvl + m*nrhs + (smlsiz+1)*(smlsiz+1); | |||
2350 | ||||
2351 | octave_idx_type addend = m; | |||
2352 | ||||
2353 | if (2*m-4 > addend) | |||
2354 | addend = 2*m-4; | |||
2355 | ||||
2356 | if (nrhs > addend) | |||
2357 | addend = nrhs; | |||
2358 | ||||
2359 | if (n-3*m > addend) | |||
2360 | addend = n-3*m; | |||
2361 | ||||
2362 | if (wlalsd > addend) | |||
2363 | addend = wlalsd; | |||
2364 | ||||
2365 | const octave_idx_type lworkaround = 4*m + m*m + addend; | |||
2366 | ||||
2367 | if (work(0) < lworkaround) | |||
2368 | work(0) = lworkaround; | |||
2369 | } | |||
2370 | else if (m >= n) | |||
2371 | { | |||
2372 | octave_idx_type lworkaround | |||
2373 | = 12*n + 2*n*smlsiz + 8*n*nlvl + n*nrhs + (smlsiz+1)*(smlsiz+1); | |||
2374 | ||||
2375 | if (work(0) < lworkaround) | |||
2376 | work(0) = lworkaround; | |||
2377 | } | |||
2378 | ||||
2379 | lwork = static_cast<octave_idx_type> (work(0)); | |||
2380 | work.resize (dim_vector (lwork, 1)); | |||
2381 | ||||
2382 | F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgelsd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgelsd_ (m, n, nrhs, tmp_data, m, pretval, maxmn, ps, rcon , rank, work.fortran_vec (), lwork, piwork, info); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
2383 | maxmn, ps, rcon, rank,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgelsd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgelsd_ (m, n, nrhs, tmp_data, m, pretval, maxmn, ps, rcon , rank, work.fortran_vec (), lwork, piwork, info); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
2384 | work.fortran_vec (), lwork,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgelsd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgelsd_ (m, n, nrhs, tmp_data, m, pretval, maxmn, ps, rcon , rank, work.fortran_vec (), lwork, piwork, info); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
2385 | piwork, info))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgelsd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgelsd_ (m, n, nrhs, tmp_data, m, pretval, maxmn, ps, rcon , rank, work.fortran_vec (), lwork, piwork, info); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0); | |||
2386 | ||||
2387 | if (s.elem (0) == 0.0) | |||
2388 | rcon = 0.0; | |||
2389 | else | |||
2390 | rcon = s.elem (minmn - 1) / s.elem (0); | |||
2391 | ||||
2392 | retval.resize (n, nrhs); | |||
2393 | } | |||
2394 | ||||
2395 | return retval; | |||
2396 | } | |||
2397 | ||||
2398 | ComplexMatrix | |||
2399 | Matrix::lssolve (const ComplexMatrix& b) const | |||
2400 | { | |||
2401 | ComplexMatrix tmp (*this); | |||
2402 | octave_idx_type info; | |||
2403 | octave_idx_type rank; | |||
2404 | double rcon; | |||
2405 | return tmp.lssolve (b, info, rank, rcon); | |||
2406 | } | |||
2407 | ||||
2408 | ComplexMatrix | |||
2409 | Matrix::lssolve (const ComplexMatrix& b, octave_idx_type& info) const | |||
2410 | { | |||
2411 | ComplexMatrix tmp (*this); | |||
2412 | octave_idx_type rank; | |||
2413 | double rcon; | |||
2414 | return tmp.lssolve (b, info, rank, rcon); | |||
2415 | } | |||
2416 | ||||
2417 | ComplexMatrix | |||
2418 | Matrix::lssolve (const ComplexMatrix& b, octave_idx_type& info, | |||
2419 | octave_idx_type& rank) const | |||
2420 | { | |||
2421 | ComplexMatrix tmp (*this); | |||
2422 | double rcon; | |||
2423 | return tmp.lssolve (b, info, rank, rcon); | |||
2424 | } | |||
2425 | ||||
2426 | ComplexMatrix | |||
2427 | Matrix::lssolve (const ComplexMatrix& b, octave_idx_type& info, | |||
2428 | octave_idx_type& rank, double& rcon) const | |||
2429 | { | |||
2430 | ComplexMatrix tmp (*this); | |||
2431 | return tmp.lssolve (b, info, rank, rcon); | |||
2432 | } | |||
2433 | ||||
2434 | ColumnVector | |||
2435 | Matrix::lssolve (const ColumnVector& b) const | |||
2436 | { | |||
2437 | octave_idx_type info; | |||
2438 | octave_idx_type rank; | |||
2439 | double rcon; | |||
2440 | return lssolve (b, info, rank, rcon); | |||
2441 | } | |||
2442 | ||||
2443 | ColumnVector | |||
2444 | Matrix::lssolve (const ColumnVector& b, octave_idx_type& info) const | |||
2445 | { | |||
2446 | octave_idx_type rank; | |||
2447 | double rcon; | |||
2448 | return lssolve (b, info, rank, rcon); | |||
2449 | } | |||
2450 | ||||
2451 | ColumnVector | |||
2452 | Matrix::lssolve (const ColumnVector& b, octave_idx_type& info, | |||
2453 | octave_idx_type& rank) const | |||
2454 | { | |||
2455 | double rcon; | |||
2456 | return lssolve (b, info, rank, rcon); | |||
2457 | } | |||
2458 | ||||
2459 | ColumnVector | |||
2460 | Matrix::lssolve (const ColumnVector& b, octave_idx_type& info, | |||
2461 | octave_idx_type& rank, double &rcon) const | |||
2462 | { | |||
2463 | ColumnVector retval; | |||
2464 | ||||
2465 | octave_idx_type nrhs = 1; | |||
2466 | ||||
2467 | octave_idx_type m = rows (); | |||
2468 | octave_idx_type n = cols (); | |||
2469 | ||||
2470 | if (m != b.length ()) | |||
2471 | (*current_liboctave_error_handler) | |||
2472 | ("matrix dimension mismatch solution of linear equations"); | |||
2473 | else if (m == 0 || n == 0) | |||
2474 | retval = ColumnVector (n, 0.0); | |||
2475 | else | |||
2476 | { | |||
2477 | volatile octave_idx_type minmn = (m < n ? m : n); | |||
2478 | octave_idx_type maxmn = m > n ? m : n; | |||
2479 | rcon = -1.0; | |||
2480 | ||||
2481 | if (m != n) | |||
2482 | { | |||
2483 | retval = ColumnVector (maxmn, 0.0); | |||
2484 | ||||
2485 | for (octave_idx_type i = 0; i < m; i++) | |||
2486 | retval.elem (i) = b.elem (i); | |||
2487 | } | |||
2488 | else | |||
2489 | retval = b; | |||
2490 | ||||
2491 | Matrix atmp = *this; | |||
2492 | double *tmp_data = atmp.fortran_vec (); | |||
2493 | ||||
2494 | double *pretval = retval.fortran_vec (); | |||
2495 | Array<double> s (dim_vector (minmn, 1)); | |||
2496 | double *ps = s.fortran_vec (); | |||
2497 | ||||
2498 | // Ask DGELSD what the dimension of WORK should be. | |||
2499 | octave_idx_type lwork = -1; | |||
2500 | ||||
2501 | Array<double> work (dim_vector (1, 1)); | |||
2502 | ||||
2503 | octave_idx_type smlsiz; | |||
2504 | F77_FUNC (xilaenv, XILAENV)xilaenv_ (9, F77_CONST_CHAR_ARG2 ("DGELSD", 6)"DGELSD", | |||
2505 | F77_CONST_CHAR_ARG2 (" ", 1)" ", | |||
2506 | 0, 0, 0, 0, smlsiz | |||
2507 | F77_CHAR_ARG_LEN (6), 6 | |||
2508 | F77_CHAR_ARG_LEN (1), 1); | |||
2509 | ||||
2510 | // We compute the size of iwork because DGELSD in older versions | |||
2511 | // of LAPACK does not return it on a query call. | |||
2512 | double dminmn = static_cast<double> (minmn); | |||
2513 | double dsmlsizp1 = static_cast<double> (smlsiz+1); | |||
2514 | #if defined (HAVE_LOG21) | |||
2515 | double tmp = log2 (dminmn / dsmlsizp1); | |||
2516 | #else | |||
2517 | double tmp = log (dminmn / dsmlsizp1) / log (2.0); | |||
2518 | #endif | |||
2519 | octave_idx_type nlvl = static_cast<octave_idx_type> (tmp) + 1; | |||
2520 | if (nlvl < 0) | |||
2521 | nlvl = 0; | |||
2522 | ||||
2523 | octave_idx_type liwork = 3 * minmn * nlvl + 11 * minmn; | |||
2524 | if (liwork < 1) | |||
2525 | liwork = 1; | |||
2526 | Array<octave_idx_type> iwork (dim_vector (liwork, 1)); | |||
2527 | octave_idx_type* piwork = iwork.fortran_vec (); | |||
2528 | ||||
2529 | F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval, maxmn,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgelsd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgelsd_ (m, n, nrhs, tmp_data, m, pretval, maxmn, ps, rcon , rank, work.fortran_vec (), lwork, piwork, info); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
2530 | ps, rcon, rank, work.fortran_vec (),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgelsd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgelsd_ (m, n, nrhs, tmp_data, m, pretval, maxmn, ps, rcon , rank, work.fortran_vec (), lwork, piwork, info); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
2531 | lwork, piwork, info))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgelsd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgelsd_ (m, n, nrhs, tmp_data, m, pretval, maxmn, ps, rcon , rank, work.fortran_vec (), lwork, piwork, info); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0); | |||
2532 | ||||
2533 | lwork = static_cast<octave_idx_type> (work(0)); | |||
2534 | work.resize (dim_vector (lwork, 1)); | |||
2535 | ||||
2536 | F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgelsd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgelsd_ (m, n, nrhs, tmp_data, m, pretval, maxmn, ps, rcon , rank, work.fortran_vec (), lwork, piwork, info); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
2537 | maxmn, ps, rcon, rank,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgelsd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgelsd_ (m, n, nrhs, tmp_data, m, pretval, maxmn, ps, rcon , rank, work.fortran_vec (), lwork, piwork, info); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
2538 | work.fortran_vec (), lwork,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgelsd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgelsd_ (m, n, nrhs, tmp_data, m, pretval, maxmn, ps, rcon , rank, work.fortran_vec (), lwork, piwork, info); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
2539 | piwork, info))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgelsd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgelsd_ (m, n, nrhs, tmp_data, m, pretval, maxmn, ps, rcon , rank, work.fortran_vec (), lwork, piwork, info); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0); | |||
2540 | ||||
2541 | if (rank < minmn) | |||
2542 | { | |||
2543 | if (s.elem (0) == 0.0) | |||
2544 | rcon = 0.0; | |||
2545 | else | |||
2546 | rcon = s.elem (minmn - 1) / s.elem (0); | |||
2547 | } | |||
2548 | ||||
2549 | retval.resize (n, nrhs); | |||
2550 | } | |||
2551 | ||||
2552 | return retval; | |||
2553 | } | |||
2554 | ||||
2555 | ComplexColumnVector | |||
2556 | Matrix::lssolve (const ComplexColumnVector& b) const | |||
2557 | { | |||
2558 | ComplexMatrix tmp (*this); | |||
2559 | octave_idx_type info; | |||
2560 | octave_idx_type rank; | |||
2561 | double rcon; | |||
2562 | return tmp.lssolve (b, info, rank, rcon); | |||
2563 | } | |||
2564 | ||||
2565 | ComplexColumnVector | |||
2566 | Matrix::lssolve (const ComplexColumnVector& b, octave_idx_type& info) const | |||
2567 | { | |||
2568 | ComplexMatrix tmp (*this); | |||
2569 | octave_idx_type rank; | |||
2570 | double rcon; | |||
2571 | return tmp.lssolve (b, info, rank, rcon); | |||
2572 | } | |||
2573 | ||||
2574 | ComplexColumnVector | |||
2575 | Matrix::lssolve (const ComplexColumnVector& b, octave_idx_type& info, | |||
2576 | octave_idx_type& rank) const | |||
2577 | { | |||
2578 | ComplexMatrix tmp (*this); | |||
2579 | double rcon; | |||
2580 | return tmp.lssolve (b, info, rank, rcon); | |||
2581 | } | |||
2582 | ||||
2583 | ComplexColumnVector | |||
2584 | Matrix::lssolve (const ComplexColumnVector& b, octave_idx_type& info, | |||
2585 | octave_idx_type& rank, double &rcon) const | |||
2586 | { | |||
2587 | ComplexMatrix tmp (*this); | |||
2588 | return tmp.lssolve (b, info, rank, rcon); | |||
2589 | } | |||
2590 | ||||
2591 | Matrix& | |||
2592 | Matrix::operator += (const DiagMatrix& a) | |||
2593 | { | |||
2594 | octave_idx_type nr = rows (); | |||
2595 | octave_idx_type nc = cols (); | |||
2596 | ||||
2597 | octave_idx_type a_nr = a.rows (); | |||
2598 | octave_idx_type a_nc = a.cols (); | |||
2599 | ||||
2600 | if (nr != a_nr || nc != a_nc) | |||
2601 | { | |||
2602 | gripe_nonconformant ("operator +=", nr, nc, a_nr, a_nc); | |||
2603 | return *this; | |||
2604 | } | |||
2605 | ||||
2606 | for (octave_idx_type i = 0; i < a.length (); i++) | |||
2607 | elem (i, i) += a.elem (i, i); | |||
2608 | ||||
2609 | return *this; | |||
2610 | } | |||
2611 | ||||
2612 | Matrix& | |||
2613 | Matrix::operator -= (const DiagMatrix& a) | |||
2614 | { | |||
2615 | octave_idx_type nr = rows (); | |||
2616 | octave_idx_type nc = cols (); | |||
2617 | ||||
2618 | octave_idx_type a_nr = a.rows (); | |||
2619 | octave_idx_type a_nc = a.cols (); | |||
2620 | ||||
2621 | if (nr != a_nr || nc != a_nc) | |||
2622 | { | |||
2623 | gripe_nonconformant ("operator -=", nr, nc, a_nr, a_nc); | |||
2624 | return *this; | |||
2625 | } | |||
2626 | ||||
2627 | for (octave_idx_type i = 0; i < a.length (); i++) | |||
2628 | elem (i, i) -= a.elem (i, i); | |||
2629 | ||||
2630 | return *this; | |||
2631 | } | |||
2632 | ||||
2633 | // unary operations | |||
2634 | ||||
2635 | boolMatrix | |||
2636 | Matrix::operator ! (void) const | |||
2637 | { | |||
2638 | if (any_element_is_nan ()) | |||
2639 | gripe_nan_to_logical_conversion (); | |||
2640 | ||||
2641 | return do_mx_unary_op<bool, double> (*this, mx_inline_not); | |||
2642 | } | |||
2643 | ||||
2644 | // column vector by row vector -> matrix operations | |||
2645 | ||||
2646 | Matrix | |||
2647 | operator * (const ColumnVector& v, const RowVector& a) | |||
2648 | { | |||
2649 | Matrix retval; | |||
2650 | ||||
2651 | octave_idx_type len = v.length (); | |||
2652 | ||||
2653 | if (len != 0) | |||
2654 | { | |||
2655 | octave_idx_type a_len = a.length (); | |||
2656 | ||||
2657 | retval = Matrix (len, a_len); | |||
2658 | double *c = retval.fortran_vec (); | |||
2659 | ||||
2660 | F77_XFCN (dgemm, DGEMM, (F77_CONST_CHAR_ARG2 ("N", 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgemm_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgemm_ ("N", "N", len, a_len, 1, 1.0, v.data (), len, a.data (), 1, 0.0, c, len , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0 ) | |||
2661 | F77_CONST_CHAR_ARG2 ("N", 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgemm_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgemm_ ("N", "N", len, a_len, 1, 1.0, v.data (), len, a.data (), 1, 0.0, c, len , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0 ) | |||
2662 | len, a_len, 1, 1.0, v.data (), len,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgemm_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgemm_ ("N", "N", len, a_len, 1, 1.0, v.data (), len, a.data (), 1, 0.0, c, len , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0 ) | |||
2663 | a.data (), 1, 0.0, c, lendo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgemm_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgemm_ ("N", "N", len, a_len, 1, 1.0, v.data (), len, a.data (), 1, 0.0, c, len , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0 ) | |||
2664 | F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgemm_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgemm_ ("N", "N", len, a_len, 1, 1.0, v.data (), len, a.data (), 1, 0.0, c, len , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0 ) | |||
2665 | F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgemm_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgemm_ ("N", "N", len, a_len, 1, 1.0, v.data (), len, a.data (), 1, 0.0, c, len , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0 ); | |||
2666 | } | |||
2667 | ||||
2668 | return retval; | |||
2669 | } | |||
2670 | ||||
2671 | // other operations. | |||
2672 | ||||
2673 | bool | |||
2674 | Matrix::any_element_is_negative (bool neg_zero) const | |||
2675 | { | |||
2676 | return (neg_zero ? test_all (xnegative_sign) | |||
2677 | : do_mx_check<double> (*this, mx_inline_any_negative)); | |||
2678 | } | |||
2679 | ||||
2680 | bool | |||
2681 | Matrix::any_element_is_positive (bool neg_zero) const | |||
2682 | { | |||
2683 | return (neg_zero ? test_all (xpositive_sign) | |||
2684 | : do_mx_check<double> (*this, mx_inline_any_positive)); | |||
2685 | } | |||
2686 | ||||
2687 | bool | |||
2688 | Matrix::any_element_is_nan (void) const | |||
2689 | { | |||
2690 | return do_mx_check<double> (*this, mx_inline_any_nan); | |||
2691 | } | |||
2692 | ||||
2693 | bool | |||
2694 | Matrix::any_element_is_inf_or_nan (void) const | |||
2695 | { | |||
2696 | return ! do_mx_check<double> (*this, mx_inline_all_finite); | |||
2697 | } | |||
2698 | ||||
2699 | bool | |||
2700 | Matrix::any_element_not_one_or_zero (void) const | |||
2701 | { | |||
2702 | return ! test_all (xis_one_or_zero); | |||
2703 | } | |||
2704 | ||||
2705 | bool | |||
2706 | Matrix::all_elements_are_int_or_inf_or_nan (void) const | |||
2707 | { | |||
2708 | return test_all (xis_int_or_inf_or_nan); | |||
2709 | } | |||
2710 | ||||
2711 | // Return nonzero if any element of M is not an integer. Also extract | |||
2712 | // the largest and smallest values and return them in MAX_VAL and MIN_VAL. | |||
2713 | ||||
2714 | bool | |||
2715 | Matrix::all_integers (double& max_val, double& min_val) const | |||
2716 | { | |||
2717 | octave_idx_type nel = nelem (); | |||
2718 | ||||
2719 | if (nel > 0) | |||
2720 | { | |||
2721 | max_val = elem (0); | |||
2722 | min_val = elem (0); | |||
2723 | } | |||
2724 | else | |||
2725 | return false; | |||
2726 | ||||
2727 | for (octave_idx_type i = 0; i < nel; i++) | |||
2728 | { | |||
2729 | double val = elem (i); | |||
2730 | ||||
2731 | if (val > max_val) | |||
2732 | max_val = val; | |||
2733 | ||||
2734 | if (val < min_val) | |||
2735 | min_val = val; | |||
2736 | ||||
2737 | if (! xisinteger (val)) | |||
2738 | return false; | |||
2739 | } | |||
2740 | ||||
2741 | return true; | |||
2742 | } | |||
2743 | ||||
2744 | bool | |||
2745 | Matrix::too_large_for_float (void) const | |||
2746 | { | |||
2747 | return test_any (xtoo_large_for_float); | |||
2748 | } | |||
2749 | ||||
2750 | // FIXME: Do these really belong here? Maybe they should be in a base class? | |||
2751 | ||||
2752 | boolMatrix | |||
2753 | Matrix::all (int dim) const | |||
2754 | { | |||
2755 | return do_mx_red_op<bool, double> (*this, dim, mx_inline_all); | |||
2756 | } | |||
2757 | ||||
2758 | boolMatrix | |||
2759 | Matrix::any (int dim) const | |||
2760 | { | |||
2761 | return do_mx_red_op<bool, double> (*this, dim, mx_inline_any); | |||
2762 | } | |||
2763 | ||||
2764 | Matrix | |||
2765 | Matrix::cumprod (int dim) const | |||
2766 | { | |||
2767 | return do_mx_cum_op<double, double> (*this, dim, mx_inline_cumprod); | |||
2768 | } | |||
2769 | ||||
2770 | Matrix | |||
2771 | Matrix::cumsum (int dim) const | |||
2772 | { | |||
2773 | return do_mx_cum_op<double, double> (*this, dim, mx_inline_cumsum); | |||
2774 | } | |||
2775 | ||||
2776 | Matrix | |||
2777 | Matrix::prod (int dim) const | |||
2778 | { | |||
2779 | return do_mx_red_op<double, double> (*this, dim, mx_inline_prod); | |||
2780 | } | |||
2781 | ||||
2782 | Matrix | |||
2783 | Matrix::sum (int dim) const | |||
2784 | { | |||
2785 | return do_mx_red_op<double, double> (*this, dim, mx_inline_sum); | |||
2786 | } | |||
2787 | ||||
2788 | Matrix | |||
2789 | Matrix::sumsq (int dim) const | |||
2790 | { | |||
2791 | return do_mx_red_op<double, double> (*this, dim, mx_inline_sumsq); | |||
2792 | } | |||
2793 | ||||
2794 | Matrix | |||
2795 | Matrix::abs (void) const | |||
2796 | { | |||
2797 | return do_mx_unary_map<double, double, std::abs> (*this); | |||
2798 | } | |||
2799 | ||||
2800 | Matrix | |||
2801 | Matrix::diag (octave_idx_type k) const | |||
2802 | { | |||
2803 | return MArray<double>::diag (k); | |||
2804 | } | |||
2805 | ||||
2806 | DiagMatrix | |||
2807 | Matrix::diag (octave_idx_type m, octave_idx_type n) const | |||
2808 | { | |||
2809 | DiagMatrix retval; | |||
2810 | ||||
2811 | octave_idx_type nr = rows (); | |||
2812 | octave_idx_type nc = cols (); | |||
2813 | ||||
2814 | if (nr == 1 || nc == 1) | |||
2815 | retval = DiagMatrix (*this, m, n); | |||
2816 | else | |||
2817 | (*current_liboctave_error_handler) | |||
2818 | ("diag: expecting vector argument"); | |||
2819 | ||||
2820 | return retval; | |||
2821 | } | |||
2822 | ||||
2823 | ColumnVector | |||
2824 | Matrix::row_min (void) const | |||
2825 | { | |||
2826 | Array<octave_idx_type> dummy_idx; | |||
2827 | return row_min (dummy_idx); | |||
2828 | } | |||
2829 | ||||
2830 | ColumnVector | |||
2831 | Matrix::row_min (Array<octave_idx_type>& idx_arg) const | |||
2832 | { | |||
2833 | ColumnVector result; | |||
2834 | ||||
2835 | octave_idx_type nr = rows (); | |||
2836 | octave_idx_type nc = cols (); | |||
2837 | ||||
2838 | if (nr > 0 && nc > 0) | |||
2839 | { | |||
2840 | result.resize (nr); | |||
2841 | idx_arg.resize (dim_vector (nr, 1)); | |||
2842 | ||||
2843 | for (octave_idx_type i = 0; i < nr; i++) | |||
2844 | { | |||
2845 | octave_idx_type idx_j; | |||
2846 | ||||
2847 | double tmp_min = octave_NaN; | |||
2848 | ||||
2849 | for (idx_j = 0; idx_j < nc; idx_j++) | |||
2850 | { | |||
2851 | tmp_min = elem (i, idx_j); | |||
2852 | ||||
2853 | if (! xisnan (tmp_min)) | |||
2854 | break; | |||
2855 | } | |||
2856 | ||||
2857 | for (octave_idx_type j = idx_j+1; j < nc; j++) | |||
2858 | { | |||
2859 | double tmp = elem (i, j); | |||
2860 | ||||
2861 | if (xisnan (tmp)) | |||
2862 | continue; | |||
2863 | else if (tmp < tmp_min) | |||
2864 | { | |||
2865 | idx_j = j; | |||
2866 | tmp_min = tmp; | |||
2867 | } | |||
2868 | } | |||
2869 | ||||
2870 | result.elem (i) = tmp_min; | |||
2871 | idx_arg.elem (i) = xisnan (tmp_min) ? 0 : idx_j; | |||
2872 | } | |||
2873 | } | |||
2874 | ||||
2875 | return result; | |||
2876 | } | |||
2877 | ||||
2878 | ColumnVector | |||
2879 | Matrix::row_max (void) const | |||
2880 | { | |||
2881 | Array<octave_idx_type> dummy_idx; | |||
2882 | return row_max (dummy_idx); | |||
2883 | } | |||
2884 | ||||
2885 | ColumnVector | |||
2886 | Matrix::row_max (Array<octave_idx_type>& idx_arg) const | |||
2887 | { | |||
2888 | ColumnVector result; | |||
2889 | ||||
2890 | octave_idx_type nr = rows (); | |||
2891 | octave_idx_type nc = cols (); | |||
2892 | ||||
2893 | if (nr > 0 && nc > 0) | |||
2894 | { | |||
2895 | result.resize (nr); | |||
2896 | idx_arg.resize (dim_vector (nr, 1)); | |||
2897 | ||||
2898 | for (octave_idx_type i = 0; i < nr; i++) | |||
2899 | { | |||
2900 | octave_idx_type idx_j; | |||
2901 | ||||
2902 | double tmp_max = octave_NaN; | |||
2903 | ||||
2904 | for (idx_j = 0; idx_j < nc; idx_j++) | |||
2905 | { | |||
2906 | tmp_max = elem (i, idx_j); | |||
2907 | ||||
2908 | if (! xisnan (tmp_max)) | |||
2909 | break; | |||
2910 | } | |||
2911 | ||||
2912 | for (octave_idx_type j = idx_j+1; j < nc; j++) | |||
2913 | { | |||
2914 | double tmp = elem (i, j); | |||
2915 | ||||
2916 | if (xisnan (tmp)) | |||
2917 | continue; | |||
2918 | else if (tmp > tmp_max) | |||
2919 | { | |||
2920 | idx_j = j; | |||
2921 | tmp_max = tmp; | |||
2922 | } | |||
2923 | } | |||
2924 | ||||
2925 | result.elem (i) = tmp_max; | |||
2926 | idx_arg.elem (i) = xisnan (tmp_max) ? 0 : idx_j; | |||
2927 | } | |||
2928 | } | |||
2929 | ||||
2930 | return result; | |||
2931 | } | |||
2932 | ||||
2933 | RowVector | |||
2934 | Matrix::column_min (void) const | |||
2935 | { | |||
2936 | Array<octave_idx_type> dummy_idx; | |||
2937 | return column_min (dummy_idx); | |||
2938 | } | |||
2939 | ||||
2940 | RowVector | |||
2941 | Matrix::column_min (Array<octave_idx_type>& idx_arg) const | |||
2942 | { | |||
2943 | RowVector result; | |||
2944 | ||||
2945 | octave_idx_type nr = rows (); | |||
2946 | octave_idx_type nc = cols (); | |||
2947 | ||||
2948 | if (nr > 0 && nc > 0) | |||
2949 | { | |||
2950 | result.resize (nc); | |||
2951 | idx_arg.resize (dim_vector (1, nc)); | |||
2952 | ||||
2953 | for (octave_idx_type j = 0; j < nc; j++) | |||
2954 | { | |||
2955 | octave_idx_type idx_i; | |||
2956 | ||||
2957 | double tmp_min = octave_NaN; | |||
2958 | ||||
2959 | for (idx_i = 0; idx_i < nr; idx_i++) | |||
2960 | { | |||
2961 | tmp_min = elem (idx_i, j); | |||
2962 | ||||
2963 | if (! xisnan (tmp_min)) | |||
2964 | break; | |||
2965 | } | |||
2966 | ||||
2967 | for (octave_idx_type i = idx_i+1; i < nr; i++) | |||
2968 | { | |||
2969 | double tmp = elem (i, j); | |||
2970 | ||||
2971 | if (xisnan (tmp)) | |||
2972 | continue; | |||
2973 | else if (tmp < tmp_min) | |||
2974 | { | |||
2975 | idx_i = i; | |||
2976 | tmp_min = tmp; | |||
2977 | } | |||
2978 | } | |||
2979 | ||||
2980 | result.elem (j) = tmp_min; | |||
2981 | idx_arg.elem (j) = xisnan (tmp_min) ? 0 : idx_i; | |||
2982 | } | |||
2983 | } | |||
2984 | ||||
2985 | return result; | |||
2986 | } | |||
2987 | ||||
2988 | RowVector | |||
2989 | Matrix::column_max (void) const | |||
2990 | { | |||
2991 | Array<octave_idx_type> dummy_idx; | |||
2992 | return column_max (dummy_idx); | |||
2993 | } | |||
2994 | ||||
2995 | RowVector | |||
2996 | Matrix::column_max (Array<octave_idx_type>& idx_arg) const | |||
2997 | { | |||
2998 | RowVector result; | |||
2999 | ||||
3000 | octave_idx_type nr = rows (); | |||
3001 | octave_idx_type nc = cols (); | |||
3002 | ||||
3003 | if (nr > 0 && nc > 0) | |||
3004 | { | |||
3005 | result.resize (nc); | |||
3006 | idx_arg.resize (dim_vector (1, nc)); | |||
3007 | ||||
3008 | for (octave_idx_type j = 0; j < nc; j++) | |||
3009 | { | |||
3010 | octave_idx_type idx_i; | |||
3011 | ||||
3012 | double tmp_max = octave_NaN; | |||
3013 | ||||
3014 | for (idx_i = 0; idx_i < nr; idx_i++) | |||
3015 | { | |||
3016 | tmp_max = elem (idx_i, j); | |||
3017 | ||||
3018 | if (! xisnan (tmp_max)) | |||
3019 | break; | |||
3020 | } | |||
3021 | ||||
3022 | for (octave_idx_type i = idx_i+1; i < nr; i++) | |||
3023 | { | |||
3024 | double tmp = elem (i, j); | |||
3025 | ||||
3026 | if (xisnan (tmp)) | |||
3027 | continue; | |||
3028 | else if (tmp > tmp_max) | |||
3029 | { | |||
3030 | idx_i = i; | |||
3031 | tmp_max = tmp; | |||
3032 | } | |||
3033 | } | |||
3034 | ||||
3035 | result.elem (j) = tmp_max; | |||
3036 | idx_arg.elem (j) = xisnan (tmp_max) ? 0 : idx_i; | |||
3037 | } | |||
3038 | } | |||
3039 | ||||
3040 | return result; | |||
3041 | } | |||
3042 | ||||
3043 | std::ostream& | |||
3044 | operator << (std::ostream& os, const Matrix& a) | |||
3045 | { | |||
3046 | for (octave_idx_type i = 0; i < a.rows (); i++) | |||
3047 | { | |||
3048 | for (octave_idx_type j = 0; j < a.cols (); j++) | |||
3049 | { | |||
3050 | os << " "; | |||
3051 | octave_write_double (os, a.elem (i, j)); | |||
3052 | } | |||
3053 | os << "\n"; | |||
3054 | } | |||
3055 | return os; | |||
3056 | } | |||
3057 | ||||
3058 | std::istream& | |||
3059 | operator >> (std::istream& is, Matrix& a) | |||
3060 | { | |||
3061 | octave_idx_type nr = a.rows (); | |||
3062 | octave_idx_type nc = a.cols (); | |||
3063 | ||||
3064 | if (nr > 0 && nc > 0) | |||
3065 | { | |||
3066 | double tmp; | |||
3067 | for (octave_idx_type i = 0; i < nr; i++) | |||
3068 | for (octave_idx_type j = 0; j < nc; j++) | |||
3069 | { | |||
3070 | tmp = octave_read_value<double> (is); | |||
3071 | if (is) | |||
3072 | a.elem (i, j) = tmp; | |||
3073 | else | |||
3074 | goto done; | |||
3075 | } | |||
3076 | } | |||
3077 | ||||
3078 | done: | |||
3079 | ||||
3080 | return is; | |||
3081 | } | |||
3082 | ||||
3083 | Matrix | |||
3084 | Givens (double x, double y) | |||
3085 | { | |||
3086 | double cc, s, temp_r; | |||
3087 | ||||
3088 | F77_FUNC (dlartg, DLARTG)dlartg_ (x, y, cc, s, temp_r); | |||
3089 | ||||
3090 | Matrix g (2, 2); | |||
3091 | ||||
3092 | g.elem (0, 0) = cc; | |||
3093 | g.elem (1, 1) = cc; | |||
3094 | g.elem (0, 1) = s; | |||
3095 | g.elem (1, 0) = -s; | |||
3096 | ||||
3097 | return g; | |||
3098 | } | |||
3099 | ||||
3100 | Matrix | |||
3101 | Sylvester (const Matrix& a, const Matrix& b, const Matrix& c) | |||
3102 | { | |||
3103 | Matrix retval; | |||
3104 | ||||
3105 | // FIXME: need to check that a, b, and c are all the same size. | |||
3106 | ||||
3107 | // Compute Schur decompositions. | |||
3108 | ||||
3109 | SCHUR as (a, "U"); | |||
3110 | SCHUR bs (b, "U"); | |||
3111 | ||||
3112 | // Transform c to new coordinates. | |||
3113 | ||||
3114 | Matrix ua = as.unitary_matrix (); | |||
3115 | Matrix sch_a = as.schur_matrix (); | |||
3116 | ||||
3117 | Matrix ub = bs.unitary_matrix (); | |||
3118 | Matrix sch_b = bs.schur_matrix (); | |||
3119 | ||||
3120 | Matrix cx = ua.transpose () * c * ub; | |||
3121 | ||||
3122 | // Solve the sylvester equation, back-transform, and return the solution. | |||
3123 | ||||
3124 | octave_idx_type a_nr = a.rows (); | |||
3125 | octave_idx_type b_nr = b.rows (); | |||
3126 | ||||
3127 | double scale; | |||
3128 | octave_idx_type info; | |||
3129 | ||||
3130 | double *pa = sch_a.fortran_vec (); | |||
3131 | double *pb = sch_b.fortran_vec (); | |||
3132 | double *px = cx.fortran_vec (); | |||
3133 | ||||
3134 | F77_XFCN (dtrsyl, DTRSYL, (F77_CONST_CHAR_ARG2 ("N", 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrsyl_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrsyl_ ("N", "N", 1, a_nr, b_nr, pa, a_nr, pb, b_nr, px, a_nr, scale, info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
3135 | F77_CONST_CHAR_ARG2 ("N", 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrsyl_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrsyl_ ("N", "N", 1, a_nr, b_nr, pa, a_nr, pb, b_nr, px, a_nr, scale, info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
3136 | 1, a_nr, b_nr, pa, a_nr, pb,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrsyl_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrsyl_ ("N", "N", 1, a_nr, b_nr, pa, a_nr, pb, b_nr, px, a_nr, scale, info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
3137 | b_nr, px, a_nr, scale, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrsyl_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrsyl_ ("N", "N", 1, a_nr, b_nr, pa, a_nr, pb, b_nr, px, a_nr, scale, info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
3138 | F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrsyl_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrsyl_ ("N", "N", 1, a_nr, b_nr, pa, a_nr, pb, b_nr, px, a_nr, scale, info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
3139 | F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dtrsyl_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dtrsyl_ ("N", "N", 1, a_nr, b_nr, pa, a_nr, pb, b_nr, px, a_nr, scale, info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0); | |||
3140 | ||||
3141 | ||||
3142 | // FIXME: check info? | |||
3143 | ||||
3144 | retval = -ua*cx*ub.transpose (); | |||
3145 | ||||
3146 | return retval; | |||
3147 | } | |||
3148 | ||||
3149 | // matrix by matrix -> matrix operations | |||
3150 | ||||
3151 | /* | |||
3152 | ||||
3153 | ## Simple Dot Product, Matrix-Vector and Matrix-Matrix Unit tests | |||
3154 | %!assert ([1 2 3] * [ 4 ; 5 ; 6], 32, 1e-14) | |||
3155 | %!assert ([1 2 ; 3 4 ] * [5 ; 6], [17 ; 39 ], 1e-14) | |||
3156 | %!assert ([1 2 ; 3 4 ] * [5 6 ; 7 8], [19 22; 43 50], 1e-14) | |||
3157 | ||||
3158 | ## Test some simple identities | |||
3159 | %!shared M, cv, rv, Mt, rvt | |||
3160 | %! M = randn (10,10) + 100*eye (10,10); | |||
3161 | %! Mt = M'; | |||
3162 | %! cv = randn (10,1); | |||
3163 | %! rv = randn (1,10); | |||
3164 | %! rvt = rv'; | |||
3165 | %!assert ([M*cv,M*cv], M*[cv,cv], 1e-13) | |||
3166 | %!assert ([M'*cv,M'*cv], M'*[cv,cv], 1e-13) | |||
3167 | %!assert ([rv*M;rv*M], [rv;rv]*M, 1e-13) | |||
3168 | %!assert ([rv*M';rv*M'], [rv;rv]*M', 1e-13) | |||
3169 | %!assert (2*rv*cv, [rv,rv]*[cv;cv], 1e-13) | |||
3170 | %!assert (M'\cv, Mt\cv, 1e-14) | |||
3171 | %!assert (M'\rv', Mt\rvt, 1e-14) | |||
3172 | ||||
3173 | */ | |||
3174 | ||||
3175 | static inline char | |||
3176 | get_blas_trans_arg (bool trans) | |||
3177 | { | |||
3178 | return trans ? 'T' : 'N'; | |||
3179 | } | |||
3180 | ||||
3181 | // the general GEMM operation | |||
3182 | ||||
3183 | Matrix | |||
3184 | xgemm (const Matrix& a, const Matrix& b, | |||
3185 | blas_trans_type transa, blas_trans_type transb) | |||
3186 | { | |||
3187 | Matrix retval; | |||
3188 | ||||
3189 | bool tra = transa != blas_no_trans, trb = transb != blas_no_trans; | |||
3190 | ||||
3191 | octave_idx_type a_nr = tra ? a.cols () : a.rows (); | |||
3192 | octave_idx_type a_nc = tra ? a.rows () : a.cols (); | |||
3193 | ||||
3194 | octave_idx_type b_nr = trb ? b.cols () : b.rows (); | |||
3195 | octave_idx_type b_nc = trb ? b.rows () : b.cols (); | |||
3196 | ||||
3197 | if (a_nc != b_nr) | |||
3198 | gripe_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc); | |||
3199 | else | |||
3200 | { | |||
3201 | if (a_nr == 0 || a_nc == 0 || b_nc == 0) | |||
3202 | retval = Matrix (a_nr, b_nc, 0.0); | |||
3203 | else if (a.data () == b.data () && a_nr == b_nc && tra != trb) | |||
3204 | { | |||
3205 | octave_idx_type lda = a.rows (); | |||
3206 | ||||
3207 | retval = Matrix (a_nr, b_nc); | |||
3208 | double *c = retval.fortran_vec (); | |||
3209 | ||||
3210 | const char ctra = get_blas_trans_arg (tra); | |||
3211 | F77_XFCN (dsyrk, DSYRK, (F77_CONST_CHAR_ARG2 ("U", 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dsyrk_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dsyrk_ ("U", &ctra, a_nr, a_nc, 1.0, a.data (), lda, 0.0 , c, a_nr , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
3212 | F77_CONST_CHAR_ARG2 (&ctra, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dsyrk_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dsyrk_ ("U", &ctra, a_nr, a_nc, 1.0, a.data (), lda, 0.0 , c, a_nr , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
3213 | a_nr, a_nc, 1.0,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dsyrk_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dsyrk_ ("U", &ctra, a_nr, a_nc, 1.0, a.data (), lda, 0.0 , c, a_nr , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
3214 | a.data (), lda, 0.0, c, a_nrdo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dsyrk_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dsyrk_ ("U", &ctra, a_nr, a_nc, 1.0, a.data (), lda, 0.0 , c, a_nr , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
3215 | F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dsyrk_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dsyrk_ ("U", &ctra, a_nr, a_nc, 1.0, a.data (), lda, 0.0 , c, a_nr , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
3216 | F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dsyrk_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dsyrk_ ("U", &ctra, a_nr, a_nc, 1.0, a.data (), lda, 0.0 , c, a_nr , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0); | |||
3217 | for (int j = 0; j < a_nr; j++) | |||
3218 | for (int i = 0; i < j; i++) | |||
3219 | retval.xelem (j,i) = retval.xelem (i,j); | |||
3220 | ||||
3221 | } | |||
3222 | else | |||
3223 | { | |||
3224 | octave_idx_type lda = a.rows (), tda = a.cols (); | |||
3225 | octave_idx_type ldb = b.rows (), tdb = b.cols (); | |||
3226 | ||||
3227 | retval = Matrix (a_nr, b_nc); | |||
3228 | double *c = retval.fortran_vec (); | |||
3229 | ||||
3230 | if (b_nc == 1) | |||
3231 | { | |||
3232 | if (a_nr == 1) | |||
3233 | F77_FUNC (xddot, XDDOT)xddot_ (a_nc, a.data (), 1, b.data (), 1, *c); | |||
3234 | else | |||
3235 | { | |||
3236 | const char ctra = get_blas_trans_arg (tra); | |||
3237 | F77_XFCN (dgemv, DGEMV, (F77_CONST_CHAR_ARG2 (&ctra, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgemv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgemv_ (&ctra, lda, tda, 1.0, a.data (), lda, b.data ( ), 1, 0.0, c, 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
3238 | lda, tda, 1.0, a.data (), lda,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgemv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgemv_ (&ctra, lda, tda, 1.0, a.data (), lda, b.data ( ), 1, 0.0, c, 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
3239 | b.data (), 1, 0.0, c, 1do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgemv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgemv_ (&ctra, lda, tda, 1.0, a.data (), lda, b.data ( ), 1, 0.0, c, 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
3240 | F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgemv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgemv_ (&ctra, lda, tda, 1.0, a.data (), lda, b.data ( ), 1, 0.0, c, 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0); | |||
3241 | } | |||
3242 | } | |||
3243 | else if (a_nr == 1) | |||
3244 | { | |||
3245 | const char crevtrb = get_blas_trans_arg (! trb); | |||
3246 | F77_XFCN (dgemv, DGEMV, (F77_CONST_CHAR_ARG2 (&crevtrb, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgemv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgemv_ (&crevtrb, ldb, tdb, 1.0, b.data (), ldb, a.data (), 1, 0.0, c, 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
3247 | ldb, tdb, 1.0, b.data (), ldb,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgemv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgemv_ (&crevtrb, ldb, tdb, 1.0, b.data (), ldb, a.data (), 1, 0.0, c, 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
3248 | a.data (), 1, 0.0, c, 1do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgemv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgemv_ (&crevtrb, ldb, tdb, 1.0, b.data (), ldb, a.data (), 1, 0.0, c, 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
3249 | F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgemv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgemv_ (&crevtrb, ldb, tdb, 1.0, b.data (), ldb, a.data (), 1, 0.0, c, 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0); | |||
3250 | } | |||
3251 | else | |||
3252 | { | |||
3253 | const char ctra = get_blas_trans_arg (tra); | |||
3254 | const char ctrb = get_blas_trans_arg (trb); | |||
3255 | F77_XFCN (dgemm, DGEMM, (F77_CONST_CHAR_ARG2 (&ctra, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgemm_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgemm_ (&ctra, &ctrb, a_nr, b_nc, a_nc, 1.0, a.data (), lda, b.data (), ldb, 0.0, c, a_nr , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
3256 | F77_CONST_CHAR_ARG2 (&ctrb, 1),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgemm_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgemm_ (&ctra, &ctrb, a_nr, b_nc, a_nc, 1.0, a.data (), lda, b.data (), ldb, 0.0, c, a_nr , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
3257 | a_nr, b_nc, a_nc, 1.0, a.data (),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgemm_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgemm_ (&ctra, &ctrb, a_nr, b_nc, a_nc, 1.0, a.data (), lda, b.data (), ldb, 0.0, c, a_nr , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
3258 | lda, b.data (), ldb, 0.0, c, a_nrdo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgemm_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgemm_ (&ctra, &ctrb, a_nr, b_nc, a_nc, 1.0, a.data (), lda, b.data (), ldb, 0.0, c, a_nr , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
3259 | F77_CHAR_ARG_LEN (1)do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgemm_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgemm_ (&ctra, &ctrb, a_nr, b_nc, a_nc, 1.0, a.data (), lda, b.data (), ldb, 0.0, c, a_nr , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
3260 | F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dgemm_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgemm_ (&ctra, &ctrb, a_nr, b_nc, a_nc, 1.0, a.data (), lda, b.data (), ldb, 0.0, c, a_nr , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0); | |||
3261 | } | |||
3262 | } | |||
3263 | } | |||
3264 | ||||
3265 | return retval; | |||
3266 | } | |||
3267 | ||||
3268 | Matrix | |||
3269 | operator * (const Matrix& a, const Matrix& b) | |||
3270 | { | |||
3271 | return xgemm (a, b); | |||
3272 | } | |||
3273 | ||||
3274 | // FIXME: it would be nice to share code among the min/max functions below. | |||
3275 | ||||
3276 | #define EMPTY_RETURN_CHECK(T)if (nr == 0 || nc == 0) return T (nr, nc); \ | |||
3277 | if (nr == 0 || nc == 0) \ | |||
3278 | return T (nr, nc); | |||
3279 | ||||
3280 | Matrix | |||
3281 | min (double d, const Matrix& m) | |||
3282 | { | |||
3283 | octave_idx_type nr = m.rows (); | |||
3284 | octave_idx_type nc = m.columns (); | |||
3285 | ||||
3286 | EMPTY_RETURN_CHECK (Matrix)if (nr == 0 || nc == 0) return Matrix (nr, nc);; | |||
3287 | ||||
3288 | Matrix result (nr, nc); | |||
3289 | ||||
3290 | for (octave_idx_type j = 0; j < nc; j++) | |||
3291 | for (octave_idx_type i = 0; i < nr; i++) | |||
3292 | { | |||
3293 | octave_quit (); | |||
3294 | result (i, j) = xmin (d, m (i, j)); | |||
3295 | } | |||
3296 | ||||
3297 | return result; | |||
3298 | } | |||
3299 | ||||
3300 | Matrix | |||
3301 | min (const Matrix& m, double d) | |||
3302 | { | |||
3303 | octave_idx_type nr = m.rows (); | |||
3304 | octave_idx_type nc = m.columns (); | |||
3305 | ||||
3306 | EMPTY_RETURN_CHECK (Matrix)if (nr == 0 || nc == 0) return Matrix (nr, nc);; | |||
3307 | ||||
3308 | Matrix result (nr, nc); | |||
3309 | ||||
3310 | for (octave_idx_type j = 0; j < nc; j++) | |||
3311 | for (octave_idx_type i = 0; i < nr; i++) | |||
3312 | { | |||
3313 | octave_quit (); | |||
3314 | result (i, j) = xmin (m (i, j), d); | |||
3315 | } | |||
3316 | ||||
3317 | return result; | |||
3318 | } | |||
3319 | ||||
3320 | Matrix | |||
3321 | min (const Matrix& a, const Matrix& b) | |||
3322 | { | |||
3323 | octave_idx_type nr = a.rows (); | |||
3324 | octave_idx_type nc = a.columns (); | |||
3325 | ||||
3326 | if (nr != b.rows () || nc != b.columns ()) | |||
3327 | { | |||
3328 | (*current_liboctave_error_handler) | |||
3329 | ("two-arg min expecting args of same size"); | |||
3330 | return Matrix (); | |||
3331 | } | |||
3332 | ||||
3333 | EMPTY_RETURN_CHECK (Matrix)if (nr == 0 || nc == 0) return Matrix (nr, nc);; | |||
3334 | ||||
3335 | Matrix result (nr, nc); | |||
3336 | ||||
3337 | for (octave_idx_type j = 0; j < nc; j++) | |||
3338 | for (octave_idx_type i = 0; i < nr; i++) | |||
3339 | { | |||
3340 | octave_quit (); | |||
3341 | result (i, j) = xmin (a (i, j), b (i, j)); | |||
3342 | } | |||
3343 | ||||
3344 | return result; | |||
3345 | } | |||
3346 | ||||
3347 | Matrix | |||
3348 | max (double d, const Matrix& m) | |||
3349 | { | |||
3350 | octave_idx_type nr = m.rows (); | |||
3351 | octave_idx_type nc = m.columns (); | |||
3352 | ||||
3353 | EMPTY_RETURN_CHECK (Matrix)if (nr == 0 || nc == 0) return Matrix (nr, nc);; | |||
3354 | ||||
3355 | Matrix result (nr, nc); | |||
3356 | ||||
3357 | for (octave_idx_type j = 0; j < nc; j++) | |||
3358 | for (octave_idx_type i = 0; i < nr; i++) | |||
3359 | { | |||
3360 | octave_quit (); | |||
3361 | result (i, j) = xmax (d, m (i, j)); | |||
3362 | } | |||
3363 | ||||
3364 | return result; | |||
3365 | } | |||
3366 | ||||
3367 | Matrix | |||
3368 | max (const Matrix& m, double d) | |||
3369 | { | |||
3370 | octave_idx_type nr = m.rows (); | |||
3371 | octave_idx_type nc = m.columns (); | |||
3372 | ||||
3373 | EMPTY_RETURN_CHECK (Matrix)if (nr == 0 || nc == 0) return Matrix (nr, nc);; | |||
3374 | ||||
3375 | Matrix result (nr, nc); | |||
3376 | ||||
3377 | for (octave_idx_type j = 0; j < nc; j++) | |||
3378 | for (octave_idx_type i = 0; i < nr; i++) | |||
3379 | { | |||
3380 | octave_quit (); | |||
3381 | result (i, j) = xmax (m (i, j), d); | |||
3382 | } | |||
3383 | ||||
3384 | return result; | |||
3385 | } | |||
3386 | ||||
3387 | Matrix | |||
3388 | max (const Matrix& a, const Matrix& b) | |||
3389 | { | |||
3390 | octave_idx_type nr = a.rows (); | |||
3391 | octave_idx_type nc = a.columns (); | |||
3392 | ||||
3393 | if (nr != b.rows () || nc != b.columns ()) | |||
3394 | { | |||
3395 | (*current_liboctave_error_handler) | |||
3396 | ("two-arg max expecting args of same size"); | |||
3397 | return Matrix (); | |||
3398 | } | |||
3399 | ||||
3400 | EMPTY_RETURN_CHECK (Matrix)if (nr == 0 || nc == 0) return Matrix (nr, nc);; | |||
3401 | ||||
3402 | Matrix result (nr, nc); | |||
3403 | ||||
3404 | for (octave_idx_type j = 0; j < nc; j++) | |||
3405 | for (octave_idx_type i = 0; i < nr; i++) | |||
3406 | { | |||
3407 | octave_quit (); | |||
3408 | result (i, j) = xmax (a (i, j), b (i, j)); | |||
3409 | } | |||
3410 | ||||
3411 | return result; | |||
3412 | } | |||
3413 | ||||
3414 | Matrix linspace (const ColumnVector& x1, | |||
3415 | const ColumnVector& x2, | |||
3416 | octave_idx_type n) | |||
3417 | ||||
3418 | { | |||
3419 | if (n < 1) n = 1; | |||
3420 | ||||
3421 | octave_idx_type m = x1.length (); | |||
3422 | ||||
3423 | if (x2.length () != m) | |||
3424 | (*current_liboctave_error_handler) | |||
3425 | ("linspace: vectors must be of equal length"); | |||
3426 | ||||
3427 | NoAlias<Matrix> retval; | |||
3428 | ||||
3429 | retval.clear (m, n); | |||
3430 | for (octave_idx_type i = 0; i < m; i++) | |||
3431 | retval(i, 0) = x1(i); | |||
3432 | ||||
3433 | // The last column is not needed while using delta. | |||
3434 | double *delta = &retval(0, n-1); | |||
3435 | for (octave_idx_type i = 0; i < m; i++) | |||
3436 | delta[i] = (x2(i) - x1(i)) / (n - 1); | |||
3437 | ||||
3438 | for (octave_idx_type j = 1; j < n-1; j++) | |||
3439 | for (octave_idx_type i = 0; i < m; i++) | |||
3440 | retval(i, j) = x1(i) + j*delta[i]; | |||
3441 | ||||
3442 | for (octave_idx_type i = 0; i < m; i++) | |||
3443 | retval(i, n-1) = x2(i); | |||
3444 | ||||
3445 | return retval; | |||
3446 | } | |||
3447 | ||||
3448 | MS_CMP_OPS (Matrix, double)boolMatrix mx_el_lt (const Matrix& m, const double& s ) { return do_ms_binary_op<bool, Matrix::element_type, double > (m, s, mx_inline_lt); } boolMatrix mx_el_le (const Matrix & m, const double& s) { return do_ms_binary_op<bool , Matrix::element_type, double> (m, s, mx_inline_le); } boolMatrix mx_el_ge (const Matrix& m, const double& s) { return do_ms_binary_op<bool, Matrix::element_type, double> (m , s, mx_inline_ge); } boolMatrix mx_el_gt (const Matrix& m , const double& s) { return do_ms_binary_op<bool, Matrix ::element_type, double> (m, s, mx_inline_gt); } boolMatrix mx_el_eq (const Matrix& m, const double& s) { return do_ms_binary_op<bool, Matrix::element_type, double> (m , s, mx_inline_eq); } boolMatrix mx_el_ne (const Matrix& m , const double& s) { return do_ms_binary_op<bool, Matrix ::element_type, double> (m, s, mx_inline_ne); } | |||
3449 | MS_BOOL_OPS (Matrix, double)boolMatrix mx_el_and (const Matrix& m, const double& s ) { if (do_mx_check (m, mx_inline_any_nan<Matrix::element_type >)) gripe_nan_to_logical_conversion (); if (xisnan (s)) gripe_nan_to_logical_conversion (); return do_ms_binary_op<bool, Matrix::element_type, double > (m, s, mx_inline_and); } boolMatrix mx_el_or (const Matrix & m, const double& s) { if (do_mx_check (m, mx_inline_any_nan <Matrix::element_type>)) gripe_nan_to_logical_conversion (); if (xisnan (s)) gripe_nan_to_logical_conversion (); return do_ms_binary_op<bool, Matrix::element_type, double> (m , s, mx_inline_or); } | |||
3450 | ||||
3451 | SM_CMP_OPS (double, Matrix)boolMatrix mx_el_lt (const double& s, const Matrix& m ) { return do_sm_binary_op<bool, double, Matrix::element_type > (s, m, mx_inline_lt); } boolMatrix mx_el_le (const double & s, const Matrix& m) { return do_sm_binary_op<bool , double, Matrix::element_type> (s, m, mx_inline_le); } boolMatrix mx_el_ge (const double& s, const Matrix& m) { return do_sm_binary_op<bool, double, Matrix::element_type> (s , m, mx_inline_ge); } boolMatrix mx_el_gt (const double& s , const Matrix& m) { return do_sm_binary_op<bool, double , Matrix::element_type> (s, m, mx_inline_gt); } boolMatrix mx_el_eq (const double& s, const Matrix& m) { return do_sm_binary_op<bool, double, Matrix::element_type> (s , m, mx_inline_eq); } boolMatrix mx_el_ne (const double& s , const Matrix& m) { return do_sm_binary_op<bool, double , Matrix::element_type> (s, m, mx_inline_ne); } | |||
3452 | SM_BOOL_OPS (double, Matrix)boolMatrix mx_el_and (const double& s, const Matrix& m ) { if (xisnan (s)) gripe_nan_to_logical_conversion (); if (do_mx_check (m, mx_inline_any_nan<Matrix::element_type>)) gripe_nan_to_logical_conversion (); return do_sm_binary_op<bool, double, Matrix::element_type > (s, m, mx_inline_and); } boolMatrix mx_el_or (const double & s, const Matrix& m) { if (xisnan (s)) gripe_nan_to_logical_conversion (); if (do_mx_check (m, mx_inline_any_nan<Matrix::element_type >)) gripe_nan_to_logical_conversion (); return do_sm_binary_op <bool, double, Matrix::element_type> (s, m, mx_inline_or ); } | |||
3453 | ||||
3454 | MM_CMP_OPS (Matrix, Matrix)boolMatrix mx_el_lt (const Matrix& m1, const Matrix& m2 ) { return do_mm_binary_op<bool, Matrix::element_type, Matrix ::element_type> (m1, m2, mx_inline_lt, mx_inline_lt, mx_inline_lt , "mx_el_lt"); } boolMatrix mx_el_le (const Matrix& m1, const Matrix& m2) { return do_mm_binary_op<bool, Matrix::element_type , Matrix::element_type> (m1, m2, mx_inline_le, mx_inline_le , mx_inline_le, "mx_el_le"); } boolMatrix mx_el_ge (const Matrix & m1, const Matrix& m2) { return do_mm_binary_op<bool , Matrix::element_type, Matrix::element_type> (m1, m2, mx_inline_ge , mx_inline_ge, mx_inline_ge, "mx_el_ge"); } boolMatrix mx_el_gt (const Matrix& m1, const Matrix& m2) { return do_mm_binary_op <bool, Matrix::element_type, Matrix::element_type> (m1, m2, mx_inline_gt, mx_inline_gt, mx_inline_gt, "mx_el_gt"); } boolMatrix mx_el_eq (const Matrix& m1, const Matrix& m2) { return do_mm_binary_op<bool, Matrix::element_type, Matrix ::element_type> (m1, m2, mx_inline_eq, mx_inline_eq, mx_inline_eq , "mx_el_eq"); } boolMatrix mx_el_ne (const Matrix& m1, const Matrix& m2) { return do_mm_binary_op<bool, Matrix::element_type , Matrix::element_type> (m1, m2, mx_inline_ne, mx_inline_ne , mx_inline_ne, "mx_el_ne"); } | |||
3455 | MM_BOOL_OPS (Matrix, Matrix)boolMatrix mx_el_and (const Matrix& m1, const Matrix& m2) { if (do_mx_check (m1, mx_inline_any_nan<Matrix::element_type >)) gripe_nan_to_logical_conversion (); if (do_mx_check (m2 , mx_inline_any_nan<Matrix::element_type>)) gripe_nan_to_logical_conversion (); return do_mm_binary_op<bool, Matrix::element_type, Matrix ::element_type> (m1, m2, mx_inline_and, mx_inline_and, mx_inline_and , "mx_el_and"); } boolMatrix mx_el_or (const Matrix& m1, const Matrix& m2) { if (do_mx_check (m1, mx_inline_any_nan< Matrix::element_type>)) gripe_nan_to_logical_conversion () ; if (do_mx_check (m2, mx_inline_any_nan<Matrix::element_type >)) gripe_nan_to_logical_conversion (); return do_mm_binary_op <bool, Matrix::element_type, Matrix::element_type> (m1, m2, mx_inline_or, mx_inline_or, mx_inline_or, "mx_el_or"); } |