File: | liboctave/numeric/floatSCHUR.cc |
Location: | line 149, column 3 |
Description: | Undefined or garbage value returned to caller |
1 | /* | |||
2 | ||||
3 | Copyright (C) 1994-2013 John W. Eaton | |||
4 | ||||
5 | This file is part of Octave. | |||
6 | ||||
7 | Octave is free software; you can redistribute it and/or modify it | |||
8 | under the terms of the GNU General Public License as published by the | |||
9 | Free Software Foundation; either version 3 of the License, or (at your | |||
10 | option) any later version. | |||
11 | ||||
12 | Octave is distributed in the hope that it will be useful, but WITHOUT | |||
13 | ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |||
14 | FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |||
15 | for more details. | |||
16 | ||||
17 | You should have received a copy of the GNU General Public License | |||
18 | along with Octave; see the file COPYING. If not, see | |||
19 | <http://www.gnu.org/licenses/>. | |||
20 | ||||
21 | */ | |||
22 | ||||
23 | #ifdef HAVE_CONFIG_H1 | |||
24 | #include <config.h> | |||
25 | #endif | |||
26 | ||||
27 | #include <iostream> | |||
28 | ||||
29 | #include "floatSCHUR.h" | |||
30 | #include "f77-fcn.h" | |||
31 | #include "lo-error.h" | |||
32 | ||||
33 | extern "C" | |||
34 | { | |||
35 | F77_RET_Tint | |||
36 | F77_FUNC (sgeesx, SGEESX)sgeesx_ (F77_CONST_CHAR_ARG_DECLconst char *, | |||
37 | F77_CONST_CHAR_ARG_DECLconst char *, | |||
38 | FloatSCHUR::select_function, | |||
39 | F77_CONST_CHAR_ARG_DECLconst char *, | |||
40 | const octave_idx_type&, float*, | |||
41 | const octave_idx_type&, octave_idx_type&, | |||
42 | float*, float*, float*, const octave_idx_type&, | |||
43 | float&, float&, float*, const octave_idx_type&, | |||
44 | octave_idx_type*, const octave_idx_type&, | |||
45 | octave_idx_type*, octave_idx_type& | |||
46 | F77_CHAR_ARG_LEN_DECL, long | |||
47 | F77_CHAR_ARG_LEN_DECL, long | |||
48 | F77_CHAR_ARG_LEN_DECL, long); | |||
49 | } | |||
50 | ||||
51 | static octave_idx_type | |||
52 | select_ana (const float& a, const float&) | |||
53 | { | |||
54 | return (a < 0.0); | |||
55 | } | |||
56 | ||||
57 | static octave_idx_type | |||
58 | select_dig (const float& a, const float& b) | |||
59 | { | |||
60 | return (hypot (a, b) < 1.0); | |||
61 | } | |||
62 | ||||
63 | octave_idx_type | |||
64 | FloatSCHUR::init (const FloatMatrix& a, const std::string& ord, | |||
65 | bool calc_unitary) | |||
66 | { | |||
67 | octave_idx_type a_nr = a.rows (); | |||
68 | octave_idx_type a_nc = a.cols (); | |||
69 | ||||
70 | if (a_nr != a_nc) | |||
| ||||
71 | { | |||
72 | (*current_liboctave_error_handler) ("FloatSCHUR requires square matrix"); | |||
73 | return -1; | |||
74 | } | |||
75 | else if (a_nr == 0) | |||
76 | { | |||
77 | schur_mat.clear (); | |||
78 | unitary_mat.clear (); | |||
79 | return 0; | |||
80 | } | |||
81 | ||||
82 | // Workspace requirements may need to be fixed if any of the | |||
83 | // following change. | |||
84 | ||||
85 | char jobvs; | |||
86 | char sense = 'N'; | |||
87 | char sort = 'N'; | |||
88 | ||||
89 | if (calc_unitary) | |||
90 | jobvs = 'V'; | |||
91 | else | |||
92 | jobvs = 'N'; | |||
93 | ||||
94 | char ord_char = ord.empty () ? 'U' : ord[0]; | |||
95 | ||||
96 | if (ord_char == 'A' || ord_char == 'D' || ord_char == 'a' || ord_char == 'd') | |||
97 | sort = 'S'; | |||
98 | ||||
99 | if (ord_char == 'A' || ord_char == 'a') | |||
100 | selector = select_ana; | |||
101 | else if (ord_char == 'D' || ord_char == 'd') | |||
102 | selector = select_dig; | |||
103 | else | |||
104 | selector = 0; | |||
105 | ||||
106 | octave_idx_type n = a_nc; | |||
107 | octave_idx_type lwork = 8 * n; | |||
108 | octave_idx_type liwork = 1; | |||
109 | octave_idx_type info; | |||
110 | octave_idx_type sdim; | |||
111 | float rconde; | |||
112 | float rcondv; | |||
113 | ||||
114 | schur_mat = a; | |||
115 | ||||
116 | if (calc_unitary) | |||
117 | unitary_mat.clear (n, n); | |||
118 | ||||
119 | float *s = schur_mat.fortran_vec (); | |||
120 | float *q = unitary_mat.fortran_vec (); | |||
121 | ||||
122 | Array<float> wr (dim_vector (n, 1)); | |||
123 | float *pwr = wr.fortran_vec (); | |||
124 | ||||
125 | Array<float> wi (dim_vector (n, 1)); | |||
126 | float *pwi = wi.fortran_vec (); | |||
127 | ||||
128 | Array<float> work (dim_vector (lwork, 1)); | |||
129 | float *pwork = work.fortran_vec (); | |||
130 | ||||
131 | // BWORK is not referenced for the non-ordered Schur routine. | |||
132 | octave_idx_type ntmp = (ord_char == 'N' || ord_char == 'n') ? 0 : n; | |||
133 | Array<octave_idx_type> bwork (dim_vector (ntmp, 1)); | |||
134 | octave_idx_type *pbwork = bwork.fortran_vec (); | |||
135 | ||||
136 | Array<octave_idx_type> iwork (dim_vector (liwork, 1)); | |||
137 | octave_idx_type *piwork = iwork.fortran_vec (); | |||
138 | ||||
139 | F77_XFCN (sgeesx, SGEESX, (F77_CONST_CHAR_ARG2 (&jobvs, 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" , "sgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgeesx_ (&jobvs, &sort, selector, &sense, n, s , n, sdim, pwr, pwi, q, n, rconde, rcondv, pwork, lwork, piwork , liwork, pbwork, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
140 | F77_CONST_CHAR_ARG2 (&sort, 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" , "sgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgeesx_ (&jobvs, &sort, selector, &sense, n, s , n, sdim, pwr, pwi, q, n, rconde, rcondv, pwork, lwork, piwork , liwork, pbwork, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
141 | selector,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "sgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgeesx_ (&jobvs, &sort, selector, &sense, n, s , n, sdim, pwr, pwi, q, n, rconde, rcondv, pwork, lwork, piwork , liwork, pbwork, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
142 | F77_CONST_CHAR_ARG2 (&sense, 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" , "sgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgeesx_ (&jobvs, &sort, selector, &sense, n, s , n, sdim, pwr, pwi, q, n, rconde, rcondv, pwork, lwork, piwork , liwork, pbwork, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
143 | n, s, n, sdim, pwr, pwi, q, n, rconde, rcondv,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "sgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgeesx_ (&jobvs, &sort, selector, &sense, n, s , n, sdim, pwr, pwi, q, n, rconde, rcondv, pwork, lwork, piwork , liwork, pbwork, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
144 | pwork, lwork, piwork, liwork, pbwork, 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" , "sgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgeesx_ (&jobvs, &sort, selector, &sense, n, s , n, sdim, pwr, pwi, q, n, rconde, rcondv, pwork, lwork, piwork , liwork, pbwork, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
145 | 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" , "sgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgeesx_ (&jobvs, &sort, selector, &sense, n, s , n, sdim, pwr, pwi, q, n, rconde, rcondv, pwork, lwork, piwork , liwork, pbwork, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
146 | 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" , "sgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgeesx_ (&jobvs, &sort, selector, &sense, n, s , n, sdim, pwr, pwi, q, n, rconde, rcondv, pwork, lwork, piwork , liwork, pbwork, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
147 | 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" , "sgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgeesx_ (&jobvs, &sort, selector, &sense, n, s , n, sdim, pwr, pwi, q, n, rconde, rcondv, pwork, lwork, piwork , liwork, pbwork, info , 1 , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0); | |||
148 | ||||
149 | return info; | |||
| ||||
150 | } | |||
151 | ||||
152 | FloatSCHUR::FloatSCHUR (const FloatMatrix& s, const FloatMatrix& u) | |||
153 | : schur_mat (s), unitary_mat (u), selector (0) | |||
154 | { | |||
155 | octave_idx_type n = s.rows (); | |||
156 | if (s.columns () != n || u.rows () != n || u.columns () != n) | |||
157 | (*current_liboctave_error_handler) | |||
158 | ("schur: inconsistent matrix dimensions"); | |||
159 | } | |||
160 | ||||
161 | std::ostream& | |||
162 | operator << (std::ostream& os, const FloatSCHUR& a) | |||
163 | { | |||
164 | os << a.schur_matrix () << "\n"; | |||
165 | os << a.unitary_matrix () << "\n"; | |||
166 | ||||
167 | return os; | |||
168 | } |