Bug Summary

File:liboctave/numeric/dbleSVD.cc
Location:line 200, column 3
Description:Undefined or garbage value returned to caller

Annotated Source Code

1/*
2
3Copyright (C) 1994-2013 John W. Eaton
4
5This file is part of Octave.
6
7Octave is free software; you can redistribute it and/or modify it
8under the terms of the GNU General Public License as published by the
9Free Software Foundation; either version 3 of the License, or (at your
10option) any later version.
11
12Octave is distributed in the hope that it will be useful, but WITHOUT
13ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15for more details.
16
17You should have received a copy of the GNU General Public License
18along 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
33extern "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
57Matrix
58SVD::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
70Matrix
71SVD::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
83octave_idx_type
84SVD::init (const Matrix& a, SVD::type svd_type, SVD::driver svd_driver)
85{
86 octave_idx_type info;
1
'info' declared without an initial value
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;
2
'?' condition is false
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)
3
Control jumps to the 'default' case at line 125
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;
4
Execution continues on line 129
127 }
128
129 type_computed = svd_type;
130
131 if (! (jobu == 'N' || jobu == 'O'))
5
Taking true branch
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'))
6
Taking true branch
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)
7
Assuming 'svd_driver' is equal to GESVD
8
Taking true branch
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'))
9
Taking true branch
198 right_sm = right_sm.transpose ();
199
200 return info;
10
Undefined or garbage value returned to caller
201}
202
203std::ostream&
204operator << (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}