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