File: | liboctave/numeric/dbleSVD.cc |
Location: | line 200, 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 "dbleSVD.h" | |||
30 | #include "f77-fcn.h" | |||
31 | #include "oct-locbuf.h" | |||
32 | ||||
33 | extern "C" | |||
34 | { | |||
35 | F77_RET_Tint | |||
36 | F77_FUNC (dgesvd, DGESVD)dgesvd_ (F77_CONST_CHAR_ARG_DECLconst char *, | |||
37 | F77_CONST_CHAR_ARG_DECLconst char *, | |||
38 | const octave_idx_type&, const octave_idx_type&, | |||
39 | double*, const octave_idx_type&, double*, | |||
40 | double*, const octave_idx_type&, double*, | |||
41 | const octave_idx_type&, double*, | |||
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 (dgesdd, DGESDD)dgesdd_ (F77_CONST_CHAR_ARG_DECLconst char *, | |||
48 | const octave_idx_type&, const octave_idx_type&, | |||
49 | double*, const octave_idx_type&, double*, | |||
50 | double*, const octave_idx_type&, double*, | |||
51 | const octave_idx_type&, double*, | |||
52 | const octave_idx_type&, octave_idx_type *, | |||
53 | octave_idx_type& | |||
54 | F77_CHAR_ARG_LEN_DECL, long); | |||
55 | } | |||
56 | ||||
57 | Matrix | |||
58 | SVD::left_singular_matrix (void) const | |||
59 | { | |||
60 | if (type_computed == SVD::sigma_only) | |||
61 | { | |||
62 | (*current_liboctave_error_handler) | |||
63 | ("SVD: U not computed because type == SVD::sigma_only"); | |||
64 | return Matrix (); | |||
65 | } | |||
66 | else | |||
67 | return left_sm; | |||
68 | } | |||
69 | ||||
70 | Matrix | |||
71 | SVD::right_singular_matrix (void) const | |||
72 | { | |||
73 | if (type_computed == SVD::sigma_only) | |||
74 | { | |||
75 | (*current_liboctave_error_handler) | |||
76 | ("SVD: V not computed because type == SVD::sigma_only"); | |||
77 | return Matrix (); | |||
78 | } | |||
79 | else | |||
80 | return right_sm; | |||
81 | } | |||
82 | ||||
83 | octave_idx_type | |||
84 | SVD::init (const Matrix& a, SVD::type svd_type, 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 | Matrix atmp = a; | |||
92 | double *tmp_data = atmp.fortran_vec (); | |||
93 | ||||
94 | octave_idx_type min_mn = m < n ? m : n; | |||
95 | ||||
96 | char jobu = 'A'; | |||
97 | char jobv = 'A'; | |||
98 | ||||
99 | octave_idx_type ncol_u = m; | |||
100 | octave_idx_type nrow_vt = n; | |||
101 | octave_idx_type nrow_s = m; | |||
102 | octave_idx_type ncol_s = n; | |||
103 | ||||
104 | switch (svd_type) | |||
105 | { | |||
106 | case SVD::economy: | |||
107 | jobu = jobv = 'S'; | |||
108 | ncol_u = nrow_vt = nrow_s = ncol_s = min_mn; | |||
109 | break; | |||
110 | ||||
111 | case SVD::sigma_only: | |||
112 | ||||
113 | // Note: for this case, both jobu and jobv should be 'N', but | |||
114 | // there seems to be a bug in dgesvd from Lapack V2.0. To | |||
115 | // demonstrate the bug, set both jobu and jobv to 'N' and find | |||
116 | // the singular values of [eye(3), eye(3)]. The result is | |||
117 | // [-sqrt(2), -sqrt(2), -sqrt(2)]. | |||
118 | // | |||
119 | // For Lapack 3.0, this problem seems to be fixed. | |||
120 | ||||
121 | jobu = jobv = 'N'; | |||
122 | ncol_u = nrow_vt = 1; | |||
123 | break; | |||
124 | ||||
125 | default: | |||
126 | break; | |||
127 | } | |||
128 | ||||
129 | type_computed = svd_type; | |||
130 | ||||
131 | if (! (jobu == 'N' || jobu == 'O')) | |||
132 | left_sm.resize (m, ncol_u); | |||
133 | ||||
134 | double *u = left_sm.fortran_vec (); | |||
135 | ||||
136 | sigma.resize (nrow_s, ncol_s); | |||
137 | double *s_vec = sigma.fortran_vec (); | |||
138 | ||||
139 | if (! (jobv == 'N' || jobv == 'O')) | |||
140 | right_sm.resize (nrow_vt, n); | |||
141 | ||||
142 | double *vt = right_sm.fortran_vec (); | |||
143 | ||||
144 | // Query DGESVD for the correct dimension of WORK. | |||
145 | ||||
146 | octave_idx_type lwork = -1; | |||
147 | ||||
148 | Array<double> work (dim_vector (1, 1)); | |||
149 | ||||
150 | octave_idx_type one = 1; | |||
151 | octave_idx_type m1 = std::max (m, one); | |||
152 | octave_idx_type nrow_vt1 = std::max (nrow_vt, one); | |||
153 | ||||
154 | if (svd_driver == SVD::GESVD) | |||
155 | { | |||
156 | F77_XFCN (dgesvd, DGESVD, (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" , "dgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgesvd_ (&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) | |||
157 | 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" , "dgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgesvd_ (&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 | 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" , "dgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgesvd_ (&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 | 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" , "dgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgesvd_ (&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 | 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" , "dgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgesvd_ (&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" , "dgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgesvd_ (&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 | ||||
163 | lwork = static_cast<octave_idx_type> (work(0)); | |||
164 | work.resize (dim_vector (lwork, 1)); | |||
165 | ||||
166 | F77_XFCN (dgesvd, DGESVD, (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" , "dgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgesvd_ (&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) | |||
167 | 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" , "dgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgesvd_ (&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 | 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" , "dgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgesvd_ (&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 | 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" , "dgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgesvd_ (&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 | 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" , "dgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgesvd_ (&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" , "dgesvd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgesvd_ (&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 | ||||
173 | } | |||
174 | else if (svd_driver == SVD::GESDD) | |||
175 | { | |||
176 | assert (jobu == jobv)((jobu == jobv) ? static_cast<void> (0) : __assert_fail ("jobu == jobv", "numeric/dbleSVD.cc", 176, __PRETTY_FUNCTION__ )); | |||
177 | char jobz = jobu; | |||
178 | 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; | |||
179 | ||||
180 | F77_XFCN (dgesdd, DGESDD, (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" , "dgesdd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgesdd_ (&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) | |||
181 | 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" , "dgesdd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgesdd_ (&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 | 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" , "dgesdd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgesdd_ (&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 | 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" , "dgesdd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgesdd_ (&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 | ||||
185 | lwork = static_cast<octave_idx_type> (work(0)); | |||
186 | work.resize (dim_vector (lwork, 1)); | |||
187 | ||||
188 | F77_XFCN (dgesdd, DGESDD, (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" , "dgesdd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgesdd_ (&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) | |||
189 | 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" , "dgesdd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgesdd_ (&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 | 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" , "dgesdd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgesdd_ (&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 | 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" , "dgesdd_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dgesdd_ (&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 | ||||
193 | } | |||
194 | else | |||
195 | assert (0)((0) ? static_cast<void> (0) : __assert_fail ("0", "numeric/dbleSVD.cc" , 195, __PRETTY_FUNCTION__)); // impossible | |||
196 | ||||
197 | if (! (jobv == 'N' || jobv == 'O')) | |||
198 | right_sm = right_sm.transpose (); | |||
199 | ||||
200 | return info; | |||
| ||||
201 | } | |||
202 | ||||
203 | std::ostream& | |||
204 | operator << (std::ostream& os, const SVD& a) | |||
205 | { | |||
206 | os << a.left_singular_matrix () << "\n"; | |||
207 | os << a.singular_values () << "\n"; | |||
208 | os << a.right_singular_matrix () << "\n"; | |||
209 | ||||
210 | return os; | |||
211 | } |