Bug Summary

File:liboctave/array/fMatrix.cc
Location:line 1550, column 3
Description:Undefined or garbage value returned to caller

Annotated Source Code

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