File: | liboctave/numeric/floatSVD.cc |
Location: | line 201, 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 "floatSVD.h" | |||
30 | #include "f77-fcn.h" | |||
31 | #include "oct-locbuf.h" | |||
32 | ||||
33 | extern "C" | |||
34 | { | |||
35 | F77_RET_Tint | |||
36 | F77_FUNC (sgesvd, SGESVD)sgesvd_ (F77_CONST_CHAR_ARG_DECLconst char *, | |||
37 | F77_CONST_CHAR_ARG_DECLconst char *, | |||
38 | const octave_idx_type&, const octave_idx_type&, | |||
39 | float*, const octave_idx_type&, float*, | |||
40 | float*, const octave_idx_type&, float*, | |||
41 | const octave_idx_type&, float*, | |||
42 | const octave_idx_type&, 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 (sgesdd, SGESDD)sgesdd_ (F77_CONST_CHAR_ARG_DECLconst char *, | |||
48 | const octave_idx_type&, const octave_idx_type&, | |||
49 | float*, const octave_idx_type&, float*, | |||
50 | float*, const octave_idx_type&, float*, | |||
51 | const octave_idx_type&, float*, | |||
52 | const octave_idx_type&, octave_idx_type *, | |||
53 | octave_idx_type& | |||
54 | F77_CHAR_ARG_LEN_DECL, long); | |||
55 | } | |||
56 | ||||
57 | FloatMatrix | |||
58 | FloatSVD::left_singular_matrix (void) const | |||
59 | { | |||
60 | if (type_computed == SVD::sigma_only) | |||
61 | { | |||
62 | (*current_liboctave_error_handler) | |||
63 | ("FloatSVD: U not computed because type == SVD::sigma_only"); | |||
64 | return FloatMatrix (); | |||
65 | } | |||
66 | else | |||
67 | return left_sm; | |||
68 | } | |||
69 | ||||
70 | FloatMatrix | |||
71 | FloatSVD::right_singular_matrix (void) const | |||
72 | { | |||
73 | if (type_computed == SVD::sigma_only) | |||
74 | { | |||
75 | (*current_liboctave_error_handler) | |||
76 | ("FloatSVD: V not computed because type == SVD::sigma_only"); | |||
77 | return FloatMatrix (); | |||
78 | } | |||
79 | else | |||
80 | return right_sm; | |||
81 | } | |||
82 | ||||
83 | octave_idx_type | |||
84 | FloatSVD::init (const FloatMatrix& 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 | FloatMatrix atmp = a; | |||
93 | float *tmp_data = atmp.fortran_vec (); | |||
94 | ||||
95 | octave_idx_type min_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 | float *u = left_sm.fortran_vec (); | |||
136 | ||||
137 | sigma.resize (nrow_s, ncol_s); | |||
138 | float *s_vec = sigma.fortran_vec (); | |||
139 | ||||
140 | if (! (jobv == 'N' || jobv == 'O')) | |||
141 | right_sm.resize (nrow_vt, n); | |||
142 | ||||
143 | float *vt = right_sm.fortran_vec (); | |||
144 | ||||
145 | // Query SGESVD for the correct dimension of WORK. | |||
146 | ||||
147 | octave_idx_type lwork = -1; | |||
148 | ||||
149 | Array<float> 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 | F77_XFCN (sgesvd, SGESVD, (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" , "sgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgesvd_ (&jobu, &jobv, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, info , 1 , 1 ); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
158 | 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" , "sgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgesvd_ (&jobu, &jobv, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, info , 1 , 1 ); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
159 | 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" , "sgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgesvd_ (&jobu, &jobv, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, info , 1 , 1 ); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
160 | nrow_vt1, work.fortran_vec (), lwork, 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" , "sgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgesvd_ (&jobu, &jobv, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, info , 1 , 1 ); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
161 | 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" , "sgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgesvd_ (&jobu, &jobv, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, info , 1 , 1 ); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
162 | 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" , "sgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgesvd_ (&jobu, &jobv, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, info , 1 , 1 ); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0); | |||
163 | ||||
164 | lwork = static_cast<octave_idx_type> (work(0)); | |||
165 | work.resize (dim_vector (lwork, 1)); | |||
166 | ||||
167 | F77_XFCN (sgesvd, SGESVD, (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" , "sgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgesvd_ (&jobu, &jobv, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, info , 1 , 1 ); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
168 | 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" , "sgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgesvd_ (&jobu, &jobv, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, info , 1 , 1 ); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
169 | 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" , "sgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgesvd_ (&jobu, &jobv, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, info , 1 , 1 ); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
170 | nrow_vt1, work.fortran_vec (), lwork, 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" , "sgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgesvd_ (&jobu, &jobv, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, info , 1 , 1 ); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
171 | 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" , "sgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgesvd_ (&jobu, &jobv, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, info , 1 , 1 ); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0) | |||
172 | 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" , "sgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgesvd_ (&jobu, &jobv, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, info , 1 , 1 ); octave_interrupt_immediately--; octave_restore_current_context (saved_context); } } while (0); | |||
173 | ||||
174 | } | |||
175 | else if (svd_driver == SVD::GESDD) | |||
176 | { | |||
177 | assert (jobu == jobv)((jobu == jobv) ? static_cast<void> (0) : __assert_fail ("jobu == jobv", "numeric/floatSVD.cc", 177, __PRETTY_FUNCTION__ )); | |||
178 | char jobz = jobu; | |||
179 | 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; | |||
180 | ||||
181 | F77_XFCN (sgesdd, SGESDD, (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" , "sgesdd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgesdd_ (&jobz, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, iwork, info , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
182 | m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1,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" , "sgesdd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgesdd_ (&jobz, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, iwork, info , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
183 | work.fortran_vec (), lwork, 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" , "sgesdd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgesdd_ (&jobz, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, iwork, info , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
184 | 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" , "sgesdd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgesdd_ (&jobz, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, iwork, info , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0); | |||
185 | ||||
186 | lwork = static_cast<octave_idx_type> (work(0)); | |||
187 | work.resize (dim_vector (lwork, 1)); | |||
188 | ||||
189 | F77_XFCN (sgesdd, SGESDD, (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" , "sgesdd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgesdd_ (&jobz, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, iwork, info , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
190 | m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1,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" , "sgesdd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgesdd_ (&jobz, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, iwork, info , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
191 | work.fortran_vec (), lwork, 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" , "sgesdd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgesdd_ (&jobz, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, iwork, info , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
192 | 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" , "sgesdd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgesdd_ (&jobz, m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, work.fortran_vec (), lwork, iwork, info , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0); | |||
193 | ||||
194 | } | |||
195 | else | |||
196 | assert (0)((0) ? static_cast<void> (0) : __assert_fail ("0", "numeric/floatSVD.cc" , 196, __PRETTY_FUNCTION__)); // impossible | |||
197 | ||||
198 | if (! (jobv == 'N' || jobv == 'O')) | |||
199 | right_sm = right_sm.transpose (); | |||
200 | ||||
201 | return info; | |||
| ||||
202 | } | |||
203 | ||||
204 | std::ostream& | |||
205 | operator << (std::ostream& os, const FloatSVD& a) | |||
206 | { | |||
207 | os << a.left_singular_matrix () << "\n"; | |||
208 | os << a.singular_values () << "\n"; | |||
209 | os << a.right_singular_matrix () << "\n"; | |||
210 | ||||
211 | return os; | |||
212 | } |