File: | liboctave/numeric/floatQRP.cc |
Location: | line 80, column 7 |
Description: | Assigned value is garbage or undefined |
1 | /* | |||
2 | ||||
3 | Copyright (C) 1994-2013 John W. Eaton | |||
4 | Copyright (C) 2009 VZLU Prague | |||
5 | ||||
6 | This file is part of Octave. | |||
7 | ||||
8 | Octave is free software; you can redistribute it and/or modify it | |||
9 | under the terms of the GNU General Public License as published by the | |||
10 | Free Software Foundation; either version 3 of the License, or (at your | |||
11 | option) any later version. | |||
12 | ||||
13 | Octave is distributed in the hope that it will be useful, but WITHOUT | |||
14 | ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |||
15 | FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |||
16 | for more details. | |||
17 | ||||
18 | You should have received a copy of the GNU General Public License | |||
19 | along with Octave; see the file COPYING. If not, see | |||
20 | <http://www.gnu.org/licenses/>. | |||
21 | ||||
22 | */ | |||
23 | ||||
24 | #ifdef HAVE_CONFIG_H1 | |||
25 | #include <config.h> | |||
26 | #endif | |||
27 | ||||
28 | #include <cassert> | |||
29 | ||||
30 | #include "floatQRP.h" | |||
31 | #include "f77-fcn.h" | |||
32 | #include "lo-error.h" | |||
33 | #include "oct-locbuf.h" | |||
34 | ||||
35 | extern "C" | |||
36 | { | |||
37 | F77_RET_Tint | |||
38 | F77_FUNC (sgeqp3, SGEQP3)sgeqp3_ (const octave_idx_type&, const octave_idx_type&, | |||
39 | float*, const octave_idx_type&, octave_idx_type*, | |||
40 | float*, float*, const octave_idx_type&, | |||
41 | octave_idx_type&); | |||
42 | } | |||
43 | ||||
44 | // It would be best to share some of this code with QR class... | |||
45 | ||||
46 | FloatQRP::FloatQRP (const FloatMatrix& a, qr_type_t qr_type) | |||
47 | : FloatQR (), p () | |||
48 | { | |||
49 | init (a, qr_type); | |||
50 | } | |||
51 | ||||
52 | void | |||
53 | FloatQRP::init (const FloatMatrix& a, qr_type_t qr_type) | |||
54 | { | |||
55 | assert (qr_type != qr_type_raw)((qr_type != qr_type_raw) ? static_cast<void> (0) : __assert_fail ("qr_type != qr_type_raw", "numeric/floatQRP.cc", 55, __PRETTY_FUNCTION__ )); | |||
56 | ||||
57 | octave_idx_type m = a.rows (); | |||
58 | octave_idx_type n = a.cols (); | |||
59 | ||||
60 | octave_idx_type min_mn = m < n ? m : n; | |||
| ||||
61 | OCTAVE_LOCAL_BUFFER (float, tau, min_mn)octave_local_buffer<float> _buffer_tau (min_mn); float * tau = _buffer_tau; | |||
62 | ||||
63 | octave_idx_type info = 0; | |||
64 | ||||
65 | FloatMatrix afact = a; | |||
66 | if (m > n && qr_type == qr_type_std) | |||
67 | afact.resize (m, m); | |||
68 | ||||
69 | MArray<octave_idx_type> jpvt (dim_vector (n, 1), 0); | |||
70 | ||||
71 | if (m > 0) | |||
72 | { | |||
73 | // workspace query. | |||
74 | float rlwork; | |||
75 | F77_XFCN (sgeqp3, SGEQP3, (m, n, afact.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" , "sgeqp3_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgeqp3_ (m, n, afact.fortran_vec (), m, jpvt.fortran_vec ( ), tau, &rlwork, -1, info); octave_interrupt_immediately-- ; octave_restore_current_context (saved_context); } } while ( 0) | |||
76 | m, jpvt.fortran_vec (), tau,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "sgeqp3_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgeqp3_ (m, n, afact.fortran_vec (), m, jpvt.fortran_vec ( ), tau, &rlwork, -1, info); octave_interrupt_immediately-- ; octave_restore_current_context (saved_context); } } while ( 0) | |||
77 | &rlwork, -1, 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" , "sgeqp3_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgeqp3_ (m, n, afact.fortran_vec (), m, jpvt.fortran_vec ( ), tau, &rlwork, -1, info); octave_interrupt_immediately-- ; octave_restore_current_context (saved_context); } } while ( 0); | |||
78 | ||||
79 | // allocate buffer and do the job. | |||
80 | octave_idx_type lwork = rlwork; | |||
| ||||
81 | lwork = std::max (lwork, static_cast<octave_idx_type> (1)); | |||
82 | OCTAVE_LOCAL_BUFFER (float, work, lwork)octave_local_buffer<float> _buffer_work (lwork); float * work = _buffer_work; | |||
83 | F77_XFCN (sgeqp3, SGEQP3, (m, n, afact.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" , "sgeqp3_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgeqp3_ (m, n, afact.fortran_vec (), m, jpvt.fortran_vec ( ), tau, work, lwork, info); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
84 | m, jpvt.fortran_vec (), tau,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "sgeqp3_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgeqp3_ (m, n, afact.fortran_vec (), m, jpvt.fortran_vec ( ), tau, work, lwork, info); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
85 | work, 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" , "sgeqp3_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgeqp3_ (m, n, afact.fortran_vec (), m, jpvt.fortran_vec ( ), tau, work, lwork, info); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0); | |||
86 | } | |||
87 | else | |||
88 | for (octave_idx_type i = 0; i < n; i++) jpvt(i) = i+1; | |||
89 | ||||
90 | // Form Permutation matrix (if economy is requested, return the | |||
91 | // indices only!) | |||
92 | ||||
93 | jpvt -= static_cast<octave_idx_type> (1); | |||
94 | p = PermMatrix (jpvt, true); | |||
95 | ||||
96 | ||||
97 | form (n, afact, tau, qr_type); | |||
98 | } | |||
99 | ||||
100 | FloatRowVector | |||
101 | FloatQRP::Pvec (void) const | |||
102 | { | |||
103 | Array<float> pa (p.pvec ()); | |||
104 | FloatRowVector pv (MArray<float> (pa) + 1.0f); | |||
105 | return pv; | |||
106 | } |