File: | liboctave/numeric/CmplxSVD.cc |
Location: | line 214, 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 "CmplxSVD.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 (zgesvd, ZGESVD)zgesvd_ (F77_CONST_CHAR_ARG_DECLconst char *, | |||
36 | F77_CONST_CHAR_ARG_DECLconst char *, | |||
37 | const octave_idx_type&, const octave_idx_type&, | |||
38 | Complex*, const octave_idx_type&, | |||
39 | double*, Complex*, const octave_idx_type&, | |||
40 | Complex*, const octave_idx_type&, Complex*, | |||
41 | const octave_idx_type&, double*, octave_idx_type& | |||
42 | F77_CHAR_ARG_LEN_DECL, long | |||
43 | F77_CHAR_ARG_LEN_DECL, long); | |||
44 | ||||
45 | F77_RET_Tint | |||
46 | F77_FUNC (zgesdd, ZGESDD)zgesdd_ (F77_CONST_CHAR_ARG_DECLconst char *, | |||
47 | const octave_idx_type&, const octave_idx_type&, | |||
48 | Complex*, const octave_idx_type&, | |||
49 | double*, Complex*, const octave_idx_type&, | |||
50 | Complex*, const octave_idx_type&, Complex*, | |||
51 | const octave_idx_type&, double*, | |||
52 | octave_idx_type *, octave_idx_type& | |||
53 | F77_CHAR_ARG_LEN_DECL, long); | |||
54 | } | |||
55 | ||||
56 | ComplexMatrix | |||
57 | ComplexSVD::left_singular_matrix (void) const | |||
58 | { | |||
59 | if (type_computed == SVD::sigma_only) | |||
60 | { | |||
61 | (*current_liboctave_error_handler) | |||
62 | ("ComplexSVD: U not computed because type == SVD::sigma_only"); | |||
63 | return ComplexMatrix (); | |||
64 | } | |||
65 | else | |||
66 | return left_sm; | |||
67 | } | |||
68 | ||||
69 | ComplexMatrix | |||
70 | ComplexSVD::right_singular_matrix (void) const | |||
71 | { | |||
72 | if (type_computed == SVD::sigma_only) | |||
73 | { | |||
74 | (*current_liboctave_error_handler) | |||
75 | ("ComplexSVD: V not computed because type == SVD::sigma_only"); | |||
76 | return ComplexMatrix (); | |||
77 | } | |||
78 | else | |||
79 | return right_sm; | |||
80 | } | |||
81 | ||||
82 | octave_idx_type | |||
83 | ComplexSVD::init (const ComplexMatrix& a, SVD::type svd_type, | |||
84 | SVD::driver svd_driver) | |||
85 | { | |||
86 | octave_idx_type info; | |||
| ||||
87 | ||||
88 | octave_idx_type m = a.rows (); | |||
89 | octave_idx_type n = a.cols (); | |||
90 | ||||
91 | ComplexMatrix atmp = a; | |||
92 | Complex *tmp_data = atmp.fortran_vec (); | |||
93 | ||||
94 | octave_idx_type min_mn = m < n ? m : n; | |||
95 | octave_idx_type max_mn = m > n ? m : n; | |||
96 | ||||
97 | char jobu = 'A'; | |||
98 | char jobv = 'A'; | |||
99 | ||||
100 | octave_idx_type ncol_u = m; | |||
101 | octave_idx_type nrow_vt = n; | |||
102 | octave_idx_type nrow_s = m; | |||
103 | octave_idx_type ncol_s = n; | |||
104 | ||||
105 | switch (svd_type) | |||
106 | { | |||
107 | case SVD::economy: | |||
108 | jobu = jobv = 'S'; | |||
109 | ncol_u = nrow_vt = nrow_s = ncol_s = min_mn; | |||
110 | break; | |||
111 | ||||
112 | case SVD::sigma_only: | |||
113 | ||||
114 | // Note: for this case, both jobu and jobv should be 'N', but | |||
115 | // there seems to be a bug in dgesvd from Lapack V2.0. To | |||
116 | // demonstrate the bug, set both jobu and jobv to 'N' and find | |||
117 | // the singular values of [eye(3), eye(3)]. The result is | |||
118 | // [-sqrt(2), -sqrt(2), -sqrt(2)]. | |||
119 | // | |||
120 | // For Lapack 3.0, this problem seems to be fixed. | |||
121 | ||||
122 | jobu = jobv = 'N'; | |||
123 | ncol_u = nrow_vt = 1; | |||
124 | break; | |||
125 | ||||
126 | default: | |||
127 | break; | |||
128 | } | |||
129 | ||||
130 | type_computed = svd_type; | |||
131 | ||||
132 | if (! (jobu == 'N' || jobu == 'O')) | |||
133 | left_sm.resize (m, ncol_u); | |||
134 | ||||
135 | Complex *u = left_sm.fortran_vec (); | |||
136 | ||||
137 | sigma.resize (nrow_s, ncol_s); | |||
138 | double *s_vec = sigma.fortran_vec (); | |||
139 | ||||
140 | if (! (jobv == 'N' || jobv == 'O')) | |||
141 | right_sm.resize (nrow_vt, n); | |||
142 | ||||
143 | Complex *vt = right_sm.fortran_vec (); | |||
144 | ||||
145 | // Query ZGESVD for the correct dimension of WORK. | |||
146 | ||||
147 | octave_idx_type lwork = -1; | |||
148 | ||||
149 | Array<Complex> work (dim_vector (1, 1)); | |||
150 | ||||
151 | octave_idx_type one = 1; | |||
152 | octave_idx_type m1 = std::max (m, one); | |||
153 | octave_idx_type nrow_vt1 = std::max (nrow_vt, one); | |||
154 | ||||
155 | if (svd_driver == SVD::GESVD) | |||
156 | { | |||
157 | octave_idx_type lrwork = 5*max_mn; | |||
158 | Array<double> rwork (dim_vector (lrwork, 1)); | |||
159 | ||||
160 | F77_XFCN (zgesvd, ZGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 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" , "zgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; zgesvd_ (&jobu, &jobv, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, rwork.fortran_vec (), info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
161 | F77_CONST_CHAR_ARG2 (&jobv, 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" , "zgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; zgesvd_ (&jobu, &jobv, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, rwork.fortran_vec (), info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
162 | m, n, tmp_data, m1, s_vec, u, m1, vt,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" , "zgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; zgesvd_ (&jobu, &jobv, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, rwork.fortran_vec (), info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
163 | nrow_vt1, work.fortran_vec (), lwork,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" , "zgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; zgesvd_ (&jobu, &jobv, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, rwork.fortran_vec (), info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
164 | rwork.fortran_vec (), 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" , "zgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; zgesvd_ (&jobu, &jobv, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, rwork.fortran_vec (), info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
165 | 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" , "zgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; zgesvd_ (&jobu, &jobv, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, rwork.fortran_vec (), info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
166 | 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" , "zgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; zgesvd_ (&jobu, &jobv, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, rwork.fortran_vec (), info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0); | |||
167 | ||||
168 | lwork = static_cast<octave_idx_type> (work(0).real ()); | |||
169 | work.resize (dim_vector (lwork, 1)); | |||
170 | ||||
171 | F77_XFCN (zgesvd, ZGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 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" , "zgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; zgesvd_ (&jobu, &jobv, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, rwork.fortran_vec (), info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
172 | F77_CONST_CHAR_ARG2 (&jobv, 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" , "zgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; zgesvd_ (&jobu, &jobv, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, rwork.fortran_vec (), info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
173 | m, n, tmp_data, m1, s_vec, u, m1, vt,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" , "zgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; zgesvd_ (&jobu, &jobv, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, rwork.fortran_vec (), info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
174 | nrow_vt1, work.fortran_vec (), lwork,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" , "zgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; zgesvd_ (&jobu, &jobv, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, rwork.fortran_vec (), info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
175 | rwork.fortran_vec (), 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" , "zgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; zgesvd_ (&jobu, &jobv, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, rwork.fortran_vec (), info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
176 | 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" , "zgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; zgesvd_ (&jobu, &jobv, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, rwork.fortran_vec (), info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
177 | 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" , "zgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; zgesvd_ (&jobu, &jobv, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, rwork.fortran_vec (), info , 1 , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0); | |||
178 | } | |||
179 | else if (svd_driver == SVD::GESDD) | |||
180 | { | |||
181 | assert (jobu == jobv)((jobu == jobv) ? static_cast<void> (0) : __assert_fail ("jobu == jobv", "numeric/CmplxSVD.cc", 181, __PRETTY_FUNCTION__ )); | |||
182 | char jobz = jobu; | |||
183 | ||||
184 | octave_idx_type lrwork; | |||
185 | if (jobz == 'N') | |||
186 | lrwork = 7*min_mn; | |||
187 | else | |||
188 | lrwork = 5*min_mn*min_mn + 5*min_mn; | |||
189 | Array<double> rwork (dim_vector (lrwork, 1)); | |||
190 | ||||
191 | OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, 8*min_mn)octave_local_buffer<octave_idx_type> _buffer_iwork (8*min_mn ); octave_idx_type *iwork = _buffer_iwork; | |||
192 | ||||
193 | F77_XFCN (zgesdd, ZGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 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" , "zgesdd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; zgesdd_ (&jobz, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, rwork.fortran_vec (), iwork , info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
194 | m, n, tmp_data, m1, s_vec, u, m1, vt,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" , "zgesdd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; zgesdd_ (&jobz, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, rwork.fortran_vec (), iwork , info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
195 | nrow_vt1, work.fortran_vec (), lwork,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" , "zgesdd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; zgesdd_ (&jobz, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, rwork.fortran_vec (), iwork , info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
196 | rwork.fortran_vec (), iwork, 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" , "zgesdd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; zgesdd_ (&jobz, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, rwork.fortran_vec (), iwork , info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
197 | 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" , "zgesdd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; zgesdd_ (&jobz, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, rwork.fortran_vec (), iwork , info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0); | |||
198 | ||||
199 | lwork = static_cast<octave_idx_type> (work(0).real ()); | |||
200 | work.resize (dim_vector (lwork, 1)); | |||
201 | ||||
202 | F77_XFCN (zgesdd, ZGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 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" , "zgesdd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; zgesdd_ (&jobz, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, rwork.fortran_vec (), iwork , info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
203 | m, n, tmp_data, m1, s_vec, u, m1, vt,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" , "zgesdd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; zgesdd_ (&jobz, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, rwork.fortran_vec (), iwork , info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
204 | nrow_vt1, work.fortran_vec (), lwork,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" , "zgesdd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; zgesdd_ (&jobz, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, rwork.fortran_vec (), iwork , info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
205 | rwork.fortran_vec (), iwork, 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" , "zgesdd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; zgesdd_ (&jobz, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, rwork.fortran_vec (), iwork , info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
206 | 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" , "zgesdd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; zgesdd_ (&jobz, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, rwork.fortran_vec (), iwork , info , 1); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0); | |||
207 | } | |||
208 | else | |||
209 | assert (0)((0) ? static_cast<void> (0) : __assert_fail ("0", "numeric/CmplxSVD.cc" , 209, __PRETTY_FUNCTION__)); // impossible | |||
210 | ||||
211 | if (! (jobv == 'N' || jobv == 'O')) | |||
212 | right_sm = right_sm.hermitian (); | |||
213 | ||||
214 | return info; | |||
| ||||
215 | } |