File: | liboctave/numeric/CmplxSCHUR.cc |
Location: | line 148, 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 "CmplxSCHUR.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 (zgeesx, ZGEESX)zgeesx_ (F77_CONST_CHAR_ARG_DECLconst char *, | |||
36 | F77_CONST_CHAR_ARG_DECLconst char *, | |||
37 | ComplexSCHUR::select_function, | |||
38 | F77_CONST_CHAR_ARG_DECLconst char *, | |||
39 | const octave_idx_type&, Complex*, | |||
40 | const octave_idx_type&, octave_idx_type&, | |||
41 | Complex*, Complex*, const octave_idx_type&, | |||
42 | double&, double&, Complex*, | |||
43 | const octave_idx_type&, double*, | |||
44 | 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 | ||||
49 | F77_RET_Tint | |||
50 | F77_FUNC (zrsf2csf, ZRSF2CSF)zrsf2csf_ (const octave_idx_type&, Complex *, | |||
51 | Complex *, double *, double *); | |||
52 | } | |||
53 | ||||
54 | static octave_idx_type | |||
55 | select_ana (const Complex& a) | |||
56 | { | |||
57 | return a.real () < 0.0; | |||
58 | } | |||
59 | ||||
60 | static octave_idx_type | |||
61 | select_dig (const Complex& a) | |||
62 | { | |||
63 | return (abs (a) < 1.0); | |||
64 | } | |||
65 | ||||
66 | octave_idx_type | |||
67 | ComplexSCHUR::init (const ComplexMatrix& a, const std::string& ord, | |||
68 | bool calc_unitary) | |||
69 | { | |||
70 | octave_idx_type a_nr = a.rows (); | |||
71 | octave_idx_type a_nc = a.cols (); | |||
72 | ||||
73 | if (a_nr != a_nc) | |||
| ||||
74 | { | |||
75 | (*current_liboctave_error_handler) | |||
76 | ("ComplexSCHUR requires square matrix"); | |||
77 | return -1; | |||
78 | } | |||
79 | else if (a_nr == 0) | |||
80 | { | |||
81 | schur_mat.clear (); | |||
82 | unitary_mat.clear (); | |||
83 | return 0; | |||
84 | } | |||
85 | ||||
86 | // Workspace requirements may need to be fixed if any of the | |||
87 | // following change. | |||
88 | ||||
89 | char jobvs; | |||
90 | char sense = 'N'; | |||
91 | char sort = 'N'; | |||
92 | ||||
93 | if (calc_unitary) | |||
94 | jobvs = 'V'; | |||
95 | else | |||
96 | jobvs = 'N'; | |||
97 | ||||
98 | char ord_char = ord.empty () ? 'U' : ord[0]; | |||
99 | ||||
100 | if (ord_char == 'A' || ord_char == 'D' || ord_char == 'a' || ord_char == 'd') | |||
101 | sort = 'S'; | |||
102 | ||||
103 | if (ord_char == 'A' || ord_char == 'a') | |||
104 | selector = select_ana; | |||
105 | else if (ord_char == 'D' || ord_char == 'd') | |||
106 | selector = select_dig; | |||
107 | else | |||
108 | selector = 0; | |||
109 | ||||
110 | octave_idx_type n = a_nc; | |||
111 | octave_idx_type lwork = 8 * n; | |||
112 | octave_idx_type info; | |||
113 | octave_idx_type sdim; | |||
114 | double rconde; | |||
115 | double rcondv; | |||
116 | ||||
117 | schur_mat = a; | |||
118 | if (calc_unitary) | |||
119 | unitary_mat.clear (n, n); | |||
120 | ||||
121 | Complex *s = schur_mat.fortran_vec (); | |||
122 | Complex *q = unitary_mat.fortran_vec (); | |||
123 | ||||
124 | Array<double> rwork (dim_vector (n, 1)); | |||
125 | double *prwork = rwork.fortran_vec (); | |||
126 | ||||
127 | Array<Complex> w (dim_vector (n, 1)); | |||
128 | Complex *pw = w.fortran_vec (); | |||
129 | ||||
130 | Array<Complex> work (dim_vector (lwork, 1)); | |||
131 | Complex *pwork = work.fortran_vec (); | |||
132 | ||||
133 | // BWORK is not referenced for non-ordered Schur. | |||
134 | octave_idx_type ntmp = (ord_char == 'N' || ord_char == 'n') ? 0 : n; | |||
135 | Array<octave_idx_type> bwork (dim_vector (ntmp, 1)); | |||
136 | octave_idx_type *pbwork = bwork.fortran_vec (); | |||
137 | ||||
138 | F77_XFCN (zgeesx, ZGEESX, (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" , "zgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; zgeesx_ (&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 | 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" , "zgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; zgeesx_ (&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 | 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" , "zgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; zgeesx_ (&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 | 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" , "zgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; zgeesx_ (&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 | 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" , "zgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; zgeesx_ (&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 | 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" , "zgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; zgeesx_ (&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" , "zgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; zgeesx_ (&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" , "zgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; zgeesx_ (&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 | 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" , "zgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; zgeesx_ (&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); | |||
147 | ||||
148 | return info; | |||
| ||||
149 | } | |||
150 | ||||
151 | ComplexSCHUR::ComplexSCHUR (const ComplexMatrix& s, const ComplexMatrix& 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 | ComplexSCHUR::ComplexSCHUR (const SCHUR& 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 (double, c, n-1)octave_local_buffer<double> _buffer_c (n-1); double *c = _buffer_c; | |||
168 | OCTAVE_LOCAL_BUFFER (double, sx, n-1)octave_local_buffer<double> _buffer_sx (n-1); double *sx = _buffer_sx; | |||
169 | ||||
170 | F77_XFCN (zrsf2csf, ZRSF2CSF, (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" , "zrsf2csf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; zrsf2csf_ (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" , "zrsf2csf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; zrsf2csf_ (n, schur_mat.fortran_vec (), unitary_mat.fortran_vec (), c, sx); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0); | |||
172 | } | |||
173 | } |