File: | liboctave/numeric/dbleSCHUR.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 <iostream> | |||
28 | ||||
29 | #include "dbleSCHUR.h" | |||
30 | #include "f77-fcn.h" | |||
31 | #include "lo-error.h" | |||
32 | ||||
33 | extern "C" | |||
34 | { | |||
35 | F77_RET_Tint | |||
36 | F77_FUNC (dgeesx, DGEESX)dgeesx_ (F77_CONST_CHAR_ARG_DECLconst char *, | |||
37 | F77_CONST_CHAR_ARG_DECLconst char *, | |||
38 | SCHUR::select_function, | |||
39 | F77_CONST_CHAR_ARG_DECLconst char *, | |||
40 | const octave_idx_type&, double*, | |||
41 | const octave_idx_type&, octave_idx_type&, | |||
42 | double*, double*, double*, const octave_idx_type&, | |||
43 | double&, double&, double*, 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 double& a, const double&) | |||
53 | { | |||
54 | return (a < 0.0); | |||
55 | } | |||
56 | ||||
57 | static octave_idx_type | |||
58 | select_dig (const double& a, const double& b) | |||
59 | { | |||
60 | return (hypot (a, b) < 1.0); | |||
61 | } | |||
62 | ||||
63 | octave_idx_type | |||
64 | SCHUR::init (const Matrix& a, const std::string& ord, bool calc_unitary) | |||
65 | { | |||
66 | octave_idx_type a_nr = a.rows (); | |||
67 | octave_idx_type a_nc = a.cols (); | |||
68 | ||||
69 | if (a_nr != a_nc) | |||
| ||||
70 | { | |||
71 | (*current_liboctave_error_handler) ("SCHUR requires square matrix"); | |||
72 | return -1; | |||
73 | } | |||
74 | else if (a_nr == 0) | |||
75 | { | |||
76 | schur_mat.clear (); | |||
77 | unitary_mat.clear (); | |||
78 | return 0; | |||
79 | } | |||
80 | ||||
81 | // Workspace requirements may need to be fixed if any of the | |||
82 | // following change. | |||
83 | ||||
84 | char jobvs; | |||
85 | char sense = 'N'; | |||
86 | char sort = 'N'; | |||
87 | ||||
88 | if (calc_unitary) | |||
89 | jobvs = 'V'; | |||
90 | else | |||
91 | jobvs = 'N'; | |||
92 | ||||
93 | char ord_char = ord.empty () ? 'U' : ord[0]; | |||
94 | ||||
95 | if (ord_char == 'A' || ord_char == 'D' || ord_char == 'a' || ord_char == 'd') | |||
96 | sort = 'S'; | |||
97 | ||||
98 | if (ord_char == 'A' || ord_char == 'a') | |||
99 | selector = select_ana; | |||
100 | else if (ord_char == 'D' || ord_char == 'd') | |||
101 | selector = select_dig; | |||
102 | else | |||
103 | selector = 0; | |||
104 | ||||
105 | octave_idx_type n = a_nc; | |||
106 | octave_idx_type lwork = 8 * n; | |||
107 | octave_idx_type liwork = 1; | |||
108 | octave_idx_type info; | |||
109 | octave_idx_type sdim; | |||
110 | double rconde; | |||
111 | double rcondv; | |||
112 | ||||
113 | schur_mat = a; | |||
114 | ||||
115 | if (calc_unitary) | |||
116 | unitary_mat.clear (n, n); | |||
117 | ||||
118 | double *s = schur_mat.fortran_vec (); | |||
119 | double *q = unitary_mat.fortran_vec (); | |||
120 | ||||
121 | Array<double> wr (dim_vector (n, 1)); | |||
122 | double *pwr = wr.fortran_vec (); | |||
123 | ||||
124 | Array<double> wi (dim_vector (n, 1)); | |||
125 | double *pwi = wi.fortran_vec (); | |||
126 | ||||
127 | Array<double> work (dim_vector (lwork, 1)); | |||
128 | double *pwork = work.fortran_vec (); | |||
129 | ||||
130 | // BWORK is not referenced for the non-ordered Schur routine. | |||
131 | octave_idx_type ntmp = (ord_char == 'N' || ord_char == 'n') ? 0 : n; | |||
132 | Array<octave_idx_type> bwork (dim_vector (ntmp, 1)); | |||
133 | octave_idx_type *pbwork = bwork.fortran_vec (); | |||
134 | ||||
135 | Array<octave_idx_type> iwork (dim_vector (liwork, 1)); | |||
136 | octave_idx_type *piwork = iwork.fortran_vec (); | |||
137 | ||||
138 | F77_XFCN (dgeesx, DGEESX, (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" , "dgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgeesx_ (&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) | |||
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" , "dgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgeesx_ (&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 | 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" , "dgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgeesx_ (&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 | 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" , "dgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgeesx_ (&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 | 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" , "dgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgeesx_ (&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 | 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" , "dgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgeesx_ (&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 | 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" , "dgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgeesx_ (&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" , "dgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgeesx_ (&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" , "dgeesx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgeesx_ (&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 | ||||
148 | return info; | |||
| ||||
149 | } | |||
150 | ||||
151 | std::ostream& | |||
152 | operator << (std::ostream& os, const SCHUR& a) | |||
153 | { | |||
154 | os << a.schur_matrix () << "\n"; | |||
155 | os << a.unitary_matrix () << "\n"; | |||
156 | ||||
157 | return os; | |||
158 | } | |||
159 | ||||
160 | SCHUR::SCHUR (const Matrix& s, const Matrix& u) | |||
161 | : schur_mat (s), unitary_mat (u), selector (0) | |||
162 | { | |||
163 | octave_idx_type n = s.rows (); | |||
164 | if (s.columns () != n || u.rows () != n || u.columns () != n) | |||
165 | (*current_liboctave_error_handler) | |||
166 | ("schur: inconsistent matrix dimensions"); | |||
167 | } | |||
168 |