File: | liboctave/numeric/fCmplxSCHUR.cc |
Location: | line 147, 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 "fCmplxSCHUR.h" | |||
28 | #include "f77-fcn.h" | |||
29 | #include "lo-error.h" | |||
30 | #include "oct-locbuf.h" | |||
31 | ||||
32 | extern "C" | |||
33 | { | |||
34 | F77_RET_Tint | |||
35 | F77_FUNC (cgeesx, CGEESX)cgeesx_ (F77_CONST_CHAR_ARG_DECLconst char *, | |||
36 | F77_CONST_CHAR_ARG_DECLconst char *, | |||
37 | FloatComplexSCHUR::select_function, | |||
38 | F77_CONST_CHAR_ARG_DECLconst char *, | |||
39 | const octave_idx_type&, FloatComplex*, | |||
40 | const octave_idx_type&, octave_idx_type&, | |||
41 | FloatComplex*, FloatComplex*, | |||
42 | const octave_idx_type&, float&, float&, | |||
43 | FloatComplex*, const octave_idx_type&, | |||
44 | float*, octave_idx_type*, octave_idx_type& | |||
45 | F77_CHAR_ARG_LEN_DECL, long | |||
46 | F77_CHAR_ARG_LEN_DECL, long | |||
47 | F77_CHAR_ARG_LEN_DECL, long); | |||
48 | F77_RET_Tint | |||
49 | F77_FUNC (crsf2csf, CRSF2CSF)crsf2csf_ (const octave_idx_type&, FloatComplex *, | |||
50 | FloatComplex *, float *, float *); | |||
51 | } | |||
52 | ||||
53 | static octave_idx_type | |||
54 | select_ana (const FloatComplex& a) | |||
55 | { | |||
56 | return a.real () < 0.0; | |||
57 | } | |||
58 | ||||
59 | static octave_idx_type | |||
60 | select_dig (const FloatComplex& a) | |||
61 | { | |||
62 | return (abs (a) < 1.0); | |||
63 | } | |||
64 | ||||
65 | octave_idx_type | |||
66 | FloatComplexSCHUR::init (const FloatComplexMatrix& a, const std::string& ord, | |||
67 | bool calc_unitary) | |||
68 | { | |||
69 | octave_idx_type a_nr = a.rows (); | |||
70 | octave_idx_type a_nc = a.cols (); | |||
71 | ||||
72 | if (a_nr != a_nc) | |||
| ||||
73 | { | |||
74 | (*current_liboctave_error_handler) | |||
75 | ("FloatComplexSCHUR requires square matrix"); | |||
76 | return -1; | |||
77 | } | |||
78 | else if (a_nr == 0) | |||
79 | { | |||
80 | schur_mat.clear (); | |||
81 | unitary_mat.clear (); | |||
82 | return 0; | |||
83 | } | |||
84 | ||||
85 | // Workspace requirements may need to be fixed if any of the | |||
86 | // following change. | |||
87 | ||||
88 | char jobvs; | |||
89 | char sense = 'N'; | |||
90 | char sort = 'N'; | |||
91 | ||||
92 | if (calc_unitary) | |||
93 | jobvs = 'V'; | |||
94 | else | |||
95 | jobvs = 'N'; | |||
96 | ||||
97 | char ord_char = ord.empty () ? 'U' : ord[0]; | |||
98 | ||||
99 | if (ord_char == 'A' || ord_char == 'D' || ord_char == 'a' || ord_char == 'd') | |||
100 | sort = 'S'; | |||
101 | ||||
102 | if (ord_char == 'A' || ord_char == 'a') | |||
103 | selector = select_ana; | |||
104 | else if (ord_char == 'D' || ord_char == 'd') | |||
105 | selector = select_dig; | |||
106 | else | |||
107 | selector = 0; | |||
108 | ||||
109 | octave_idx_type n = a_nc; | |||
110 | octave_idx_type lwork = 8 * n; | |||
111 | octave_idx_type info; | |||
112 | octave_idx_type sdim; | |||
113 | float rconde; | |||
114 | float rcondv; | |||
115 | ||||
116 | schur_mat = a; | |||
117 | if (calc_unitary) | |||
118 | unitary_mat.clear (n, n); | |||
119 | ||||
120 | FloatComplex *s = schur_mat.fortran_vec (); | |||
121 | FloatComplex *q = unitary_mat.fortran_vec (); | |||
122 | ||||
123 | Array<float> rwork (dim_vector (n, 1)); | |||
124 | float *prwork = rwork.fortran_vec (); | |||
125 | ||||
126 | Array<FloatComplex> w (dim_vector (n, 1)); | |||
127 | FloatComplex *pw = w.fortran_vec (); | |||
128 | ||||
129 | Array<FloatComplex> work (dim_vector (lwork, 1)); | |||
130 | FloatComplex *pwork = work.fortran_vec (); | |||
131 | ||||
132 | // BWORK is not referenced for non-ordered Schur. | |||
133 | octave_idx_type ntmp = (ord_char == 'N' || ord_char == 'n') ? 0 : n; | |||
134 | Array<octave_idx_type> bwork (dim_vector (ntmp, 1)); | |||
135 | octave_idx_type *pbwork = bwork.fortran_vec (); | |||
136 | ||||
137 | F77_XFCN (cgeesx, CGEESX, (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" , "cgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cgeesx_ (&jobvs, &sort, selector, &sense, n, s , n, sdim, pw, q, n, rconde, rcondv, pwork, lwork, prwork, pbwork , info , 1 , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
138 | 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" , "cgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cgeesx_ (&jobvs, &sort, selector, &sense, n, s , n, sdim, pw, q, n, rconde, rcondv, pwork, lwork, prwork, pbwork , info , 1 , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
139 | 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" , "cgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cgeesx_ (&jobvs, &sort, selector, &sense, n, s , n, sdim, pw, q, n, rconde, rcondv, pwork, lwork, prwork, pbwork , info , 1 , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
140 | 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" , "cgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cgeesx_ (&jobvs, &sort, selector, &sense, n, s , n, sdim, pw, q, n, rconde, rcondv, pwork, lwork, prwork, pbwork , info , 1 , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
141 | n, s, n, sdim, pw, 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" , "cgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cgeesx_ (&jobvs, &sort, selector, &sense, n, s , n, sdim, pw, q, n, rconde, rcondv, pwork, lwork, prwork, pbwork , info , 1 , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
142 | pwork, lwork, prwork, 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" , "cgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cgeesx_ (&jobvs, &sort, selector, &sense, n, s , n, sdim, pw, q, n, rconde, rcondv, pwork, lwork, prwork, pbwork , info , 1 , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
143 | 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" , "cgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cgeesx_ (&jobvs, &sort, selector, &sense, n, s , n, sdim, pw, q, n, rconde, rcondv, pwork, lwork, prwork, pbwork , info , 1 , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
144 | 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" , "cgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cgeesx_ (&jobvs, &sort, selector, &sense, n, s , n, sdim, pw, q, n, rconde, rcondv, pwork, lwork, prwork, 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" , "cgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cgeesx_ (&jobvs, &sort, selector, &sense, n, s , n, sdim, pw, q, n, rconde, rcondv, pwork, lwork, prwork, pbwork , info , 1 , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0); | |||
146 | ||||
147 | return info; | |||
| ||||
148 | } | |||
149 | ||||
150 | FloatComplexSCHUR::FloatComplexSCHUR (const FloatComplexMatrix& s, | |||
151 | const FloatComplexMatrix& u) | |||
152 | : schur_mat (s), unitary_mat (u), selector (0) | |||
153 | { | |||
154 | octave_idx_type n = s.rows (); | |||
155 | if (s.columns () != n || u.rows () != n || u.columns () != n) | |||
156 | (*current_liboctave_error_handler) | |||
157 | ("schur: inconsistent matrix dimensions"); | |||
158 | } | |||
159 | ||||
160 | FloatComplexSCHUR::FloatComplexSCHUR (const FloatSCHUR& s) | |||
161 | : schur_mat (s.schur_matrix ()), unitary_mat (s.unitary_matrix ()), | |||
162 | selector (0) | |||
163 | { | |||
164 | octave_idx_type n = schur_mat.rows (); | |||
165 | if (n > 0) | |||
166 | { | |||
167 | OCTAVE_LOCAL_BUFFER (float, c, n-1)octave_local_buffer<float> _buffer_c (n-1); float *c = _buffer_c; | |||
168 | OCTAVE_LOCAL_BUFFER (float, sx, n-1)octave_local_buffer<float> _buffer_sx (n-1); float *sx = _buffer_sx; | |||
169 | ||||
170 | F77_XFCN (crsf2csf, CRSF2CSF, (n, schur_mat.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" , "crsf2csf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; crsf2csf_ (n, schur_mat.fortran_vec (), unitary_mat.fortran_vec (), c, sx); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
171 | unitary_mat.fortran_vec (), c, sx))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "crsf2csf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; crsf2csf_ (n, schur_mat.fortran_vec (), unitary_mat.fortran_vec (), c, sx); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0); | |||
172 | } | |||
173 | } |