File: | liboctave/numeric/fCmplxHESS.cc |
Location: | line 126, 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 "fCmplxHESS.h" | |||
28 | #include "f77-fcn.h" | |||
29 | #include "lo-error.h" | |||
30 | ||||
31 | extern "C" | |||
32 | { | |||
33 | F77_RET_Tint | |||
34 | F77_FUNC (cgebal, CGEBAL)cgebal_ (F77_CONST_CHAR_ARG_DECLconst char *, | |||
35 | const octave_idx_type&, FloatComplex*, | |||
36 | const octave_idx_type&, octave_idx_type&, | |||
37 | octave_idx_type&, float*, octave_idx_type& | |||
38 | F77_CHAR_ARG_LEN_DECL, long); | |||
39 | ||||
40 | F77_RET_Tint | |||
41 | F77_FUNC (cgehrd, CGEHRD)cgehrd_ (const octave_idx_type&, const octave_idx_type&, | |||
42 | const octave_idx_type&, FloatComplex*, | |||
43 | const octave_idx_type&, FloatComplex*, | |||
44 | FloatComplex*, const octave_idx_type&, | |||
45 | octave_idx_type&); | |||
46 | ||||
47 | F77_RET_Tint | |||
48 | F77_FUNC (cunghr, CUNGHR)cunghr_ (const octave_idx_type&, const octave_idx_type&, | |||
49 | const octave_idx_type&, FloatComplex*, | |||
50 | const octave_idx_type&, FloatComplex*, | |||
51 | FloatComplex*, const octave_idx_type&, | |||
52 | octave_idx_type&); | |||
53 | ||||
54 | F77_RET_Tint | |||
55 | F77_FUNC (cgebak, CGEBAK)cgebak_ (F77_CONST_CHAR_ARG_DECLconst char *, | |||
56 | F77_CONST_CHAR_ARG_DECLconst char *, | |||
57 | const octave_idx_type&, const octave_idx_type&, | |||
58 | const octave_idx_type&, float*, | |||
59 | const octave_idx_type&, FloatComplex*, | |||
60 | const octave_idx_type&, octave_idx_type& | |||
61 | F77_CHAR_ARG_LEN_DECL, long | |||
62 | F77_CHAR_ARG_LEN_DECL, long); | |||
63 | } | |||
64 | ||||
65 | octave_idx_type | |||
66 | FloatComplexHESS::init (const FloatComplexMatrix& a) | |||
67 | { | |||
68 | octave_idx_type a_nr = a.rows (); | |||
69 | octave_idx_type a_nc = a.cols (); | |||
70 | ||||
71 | if (a_nr != a_nc) | |||
| ||||
72 | { | |||
73 | (*current_liboctave_error_handler) | |||
74 | ("FloatComplexHESS requires square matrix"); | |||
75 | return -1; | |||
76 | } | |||
77 | ||||
78 | char job = 'N'; | |||
79 | char side = 'R'; | |||
80 | ||||
81 | octave_idx_type n = a_nc; | |||
82 | octave_idx_type lwork = 32 * n; | |||
83 | octave_idx_type info; | |||
84 | octave_idx_type ilo; | |||
85 | octave_idx_type ihi; | |||
86 | ||||
87 | hess_mat = a; | |||
88 | FloatComplex *h = hess_mat.fortran_vec (); | |||
89 | ||||
90 | Array<float> scale (dim_vector (n, 1)); | |||
91 | float *pscale = scale.fortran_vec (); | |||
92 | ||||
93 | F77_XFCN (cgebal, CGEBAL, (F77_CONST_CHAR_ARG2 (&job, 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" , "cgebal_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cgebal_ (&job, n, h, n, ilo, ihi, pscale, info , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
94 | n, h, n, ilo, ihi, pscale, 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" , "cgebal_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cgebal_ (&job, n, h, n, ilo, ihi, pscale, info , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
95 | 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" , "cgebal_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cgebal_ (&job, n, h, n, ilo, ihi, pscale, info , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0); | |||
96 | ||||
97 | Array<FloatComplex> tau (dim_vector (n-1, 1)); | |||
98 | FloatComplex *ptau = tau.fortran_vec (); | |||
99 | ||||
100 | Array<FloatComplex> work (dim_vector (lwork, 1)); | |||
101 | FloatComplex *pwork = work.fortran_vec (); | |||
102 | ||||
103 | F77_XFCN (cgehrd, CGEHRD, (n, ilo, ihi, h, n, ptau, pwork, 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" , "cgehrd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cgehrd_ (n, ilo, ihi, h, n, ptau, pwork, lwork, info); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0); | |||
104 | ||||
105 | unitary_hess_mat = hess_mat; | |||
106 | FloatComplex *z = unitary_hess_mat.fortran_vec (); | |||
107 | ||||
108 | F77_XFCN (cunghr, CUNGHR, (n, ilo, ihi, z, n, ptau, pwork,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" , "cunghr_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cunghr_ (n, ilo, ihi, z, n, ptau, pwork, lwork, info); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
109 | 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" , "cunghr_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cunghr_ (n, ilo, ihi, z, n, ptau, pwork, lwork, info); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0); | |||
110 | ||||
111 | F77_XFCN (cgebak, CGEBAK, (F77_CONST_CHAR_ARG2 (&job, 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" , "cgebak_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cgebak_ (&job, &side, n, ilo, ihi, pscale, n, z, n , info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
112 | F77_CONST_CHAR_ARG2 (&side, 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" , "cgebak_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cgebak_ (&job, &side, n, ilo, ihi, pscale, n, z, n , info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
113 | n, ilo, ihi, pscale, n, z, n, 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" , "cgebak_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cgebak_ (&job, &side, n, ilo, ihi, pscale, n, z, n , info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
114 | 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" , "cgebak_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cgebak_ (&job, &side, n, ilo, ihi, pscale, n, z, n , info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
115 | 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" , "cgebak_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cgebak_ (&job, &side, n, ilo, ihi, pscale, n, z, n , info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0); | |||
116 | ||||
117 | // If someone thinks of a more graceful way of | |||
118 | // doing this (or faster for that matter :-)), | |||
119 | // please let me know! | |||
120 | ||||
121 | if (n > 2) | |||
122 | for (octave_idx_type j = 0; j < a_nc; j++) | |||
123 | for (octave_idx_type i = j+2; i < a_nr; i++) | |||
124 | hess_mat.elem (i, j) = 0; | |||
125 | ||||
126 | return info; | |||
| ||||
127 | } |