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