Mercurial > hg > octave-nkf
annotate liboctave/numeric/dbleSVD.cc @ 20685:7fa1970a655d
pkg.m: drop check of nargout value, the interpreter already does that.
* scripts/pkg/pkg.m: the interpreter already checks if there was any variable
that got no value assigned, there's no need to make the code more
complicated to cover that. Also, there's no point in calling describe()
with different nargout since it doesn't check nargout.
| author | Carnë Draug <carandraug@octave.org> |
|---|---|
| date | Thu, 03 Sep 2015 16:21:08 +0100 |
| parents | 4197fc428c7d |
| children |
| rev | line source |
|---|---|
| 457 | 1 /* |
| 2 | |
|
19898
4197fc428c7d
maint: Update copyright notices for 2015.
John W. Eaton <jwe@octave.org>
parents:
17769
diff
changeset
|
3 Copyright (C) 1994-2015 John W. Eaton |
| 457 | 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 | |
| 7016 | 9 Free Software Foundation; either version 3 of the License, or (at your |
| 10 option) any later version. | |
| 457 | 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 | |
| 7016 | 18 along with Octave; see the file COPYING. If not, see |
| 19 <http://www.gnu.org/licenses/>. | |
| 457 | 20 |
| 21 */ | |
| 22 | |
| 23 #ifdef HAVE_CONFIG_H | |
| 1192 | 24 #include <config.h> |
| 457 | 25 #endif |
| 26 | |
| 3503 | 27 #include <iostream> |
| 1631 | 28 |
| 457 | 29 #include "dbleSVD.h" |
| 1847 | 30 #include "f77-fcn.h" |
| 10601 | 31 #include "oct-locbuf.h" |
| 457 | 32 |
| 33 extern "C" | |
| 34 { | |
| 4552 | 35 F77_RET_T |
| 36 F77_FUNC (dgesvd, DGESVD) (F77_CONST_CHAR_ARG_DECL, | |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10258
diff
changeset
|
37 F77_CONST_CHAR_ARG_DECL, |
| 11495 | 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& | |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10258
diff
changeset
|
43 F77_CHAR_ARG_LEN_DECL |
|
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10258
diff
changeset
|
44 F77_CHAR_ARG_LEN_DECL); |
| 10601 | 45 |
| 46 F77_RET_T | |
| 47 F77_FUNC (dgesdd, DGESDD) (F77_CONST_CHAR_ARG_DECL, | |
| 11495 | 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& | |
| 10601 | 54 F77_CHAR_ARG_LEN_DECL); |
| 457 | 55 } |
| 56 | |
| 1543 | 57 Matrix |
| 1545 | 58 SVD::left_singular_matrix (void) const |
| 1543 | 59 { |
| 1544 | 60 if (type_computed == SVD::sigma_only) |
| 1543 | 61 { |
| 62 (*current_liboctave_error_handler) | |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10258
diff
changeset
|
63 ("SVD: U not computed because type == SVD::sigma_only"); |
| 1543 | 64 return Matrix (); |
| 65 } | |
| 66 else | |
| 67 return left_sm; | |
| 68 } | |
| 69 | |
| 70 Matrix | |
| 1545 | 71 SVD::right_singular_matrix (void) const |
| 1543 | 72 { |
| 1544 | 73 if (type_computed == SVD::sigma_only) |
| 1543 | 74 { |
| 75 (*current_liboctave_error_handler) | |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10258
diff
changeset
|
76 ("SVD: V not computed because type == SVD::sigma_only"); |
| 1543 | 77 return Matrix (); |
| 78 } | |
| 79 else | |
| 80 return right_sm; | |
| 81 } | |
| 82 | |
| 5275 | 83 octave_idx_type |
| 10601 | 84 SVD::init (const Matrix& a, SVD::type svd_type, SVD::driver svd_driver) |
| 457 | 85 { |
| 5275 | 86 octave_idx_type info; |
| 457 | 87 |
| 5275 | 88 octave_idx_type m = a.rows (); |
| 89 octave_idx_type n = a.cols (); | |
| 457 | 90 |
| 1930 | 91 Matrix atmp = a; |
| 92 double *tmp_data = atmp.fortran_vec (); | |
| 457 | 93 |
| 5275 | 94 octave_idx_type min_mn = m < n ? m : n; |
| 457 | 95 |
| 1930 | 96 char jobu = 'A'; |
| 97 char jobv = 'A'; | |
| 537 | 98 |
| 5275 | 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; | |
| 537 | 103 |
| 1543 | 104 switch (svd_type) |
| 537 | 105 { |
| 1543 | 106 case SVD::economy: |
| 1930 | 107 jobu = jobv = 'S'; |
| 537 | 108 ncol_u = nrow_vt = nrow_s = ncol_s = min_mn; |
| 1543 | 109 break; |
| 110 | |
| 111 case SVD::sigma_only: | |
| 2621 | 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)]. | |
| 3335 | 118 // |
| 119 // For Lapack 3.0, this problem seems to be fixed. | |
| 2621 | 120 |
|
15887
8ced82e96b48
Fix segfaults with gesdd driver for svd (bug #37998).
Rik <rik@octave.org>
parents:
14138
diff
changeset
|
121 jobu = jobv = 'N'; |
| 1545 | 122 ncol_u = nrow_vt = 1; |
| 1543 | 123 break; |
| 124 | |
| 125 default: | |
| 126 break; | |
| 537 | 127 } |
| 128 | |
| 1544 | 129 type_computed = svd_type; |
| 130 | |
| 2621 | 131 if (! (jobu == 'N' || jobu == 'O')) |
| 1931 | 132 left_sm.resize (m, ncol_u); |
| 1930 | 133 |
| 134 double *u = left_sm.fortran_vec (); | |
| 135 | |
| 136 sigma.resize (nrow_s, ncol_s); | |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14138
diff
changeset
|
137 double *s_vec = sigma.fortran_vec (); |
| 1930 | 138 |
| 2621 | 139 if (! (jobv == 'N' || jobv == 'O')) |
| 1930 | 140 right_sm.resize (nrow_vt, n); |
| 141 | |
| 142 double *vt = right_sm.fortran_vec (); | |
| 457 | 143 |
|
15887
8ced82e96b48
Fix segfaults with gesdd driver for svd (bug #37998).
Rik <rik@octave.org>
parents:
14138
diff
changeset
|
144 // Query DGESVD for the correct dimension of WORK. |
| 1930 | 145 |
| 5275 | 146 octave_idx_type lwork = -1; |
| 3336 | 147 |
|
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
148 Array<double> work (dim_vector (1, 1)); |
| 457 | 149 |
| 10258 | 150 octave_idx_type one = 1; |
|
15887
8ced82e96b48
Fix segfaults with gesdd driver for svd (bug #37998).
Rik <rik@octave.org>
parents:
14138
diff
changeset
|
151 octave_idx_type m1 = std::max (m, one); |
|
8ced82e96b48
Fix segfaults with gesdd driver for svd (bug #37998).
Rik <rik@octave.org>
parents:
14138
diff
changeset
|
152 octave_idx_type nrow_vt1 = std::max (nrow_vt, one); |
|
10185
455759a5fcbe
fix norm and svd on empty matrices
Jaroslav Hajek <highegg@gmail.com>
parents:
10158
diff
changeset
|
153 |
| 10601 | 154 if (svd_driver == SVD::GESVD) |
| 155 { | |
| 156 F77_XFCN (dgesvd, DGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), | |
| 157 F77_CONST_CHAR_ARG2 (&jobv, 1), | |
| 158 m, n, tmp_data, m1, s_vec, u, m1, vt, | |
| 159 nrow_vt1, work.fortran_vec (), lwork, info | |
| 160 F77_CHAR_ARG_LEN (1) | |
| 161 F77_CHAR_ARG_LEN (1))); | |
| 162 | |
| 163 lwork = static_cast<octave_idx_type> (work(0)); | |
|
11574
a83bad07f7e3
attempt better backward compatibility for Array resize functions
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
164 work.resize (dim_vector (lwork, 1)); |
| 10601 | 165 |
| 166 F77_XFCN (dgesvd, DGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), | |
| 167 F77_CONST_CHAR_ARG2 (&jobv, 1), | |
| 168 m, n, tmp_data, m1, s_vec, u, m1, vt, | |
| 169 nrow_vt1, work.fortran_vec (), lwork, info | |
| 170 F77_CHAR_ARG_LEN (1) | |
| 171 F77_CHAR_ARG_LEN (1))); | |
| 1543 | 172 |
| 10601 | 173 } |
| 174 else if (svd_driver == SVD::GESDD) | |
| 175 { | |
| 176 assert (jobu == jobv); | |
| 177 char jobz = jobu; | |
| 178 OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, 8*min_mn); | |
| 179 | |
| 180 F77_XFCN (dgesdd, DGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1), | |
|
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
181 m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, |
|
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
182 work.fortran_vec (), lwork, iwork, info |
| 10601 | 183 F77_CHAR_ARG_LEN (1))); |
| 3336 | 184 |
| 10601 | 185 lwork = static_cast<octave_idx_type> (work(0)); |
|
11574
a83bad07f7e3
attempt better backward compatibility for Array resize functions
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
186 work.resize (dim_vector (lwork, 1)); |
| 10601 | 187 |
| 188 F77_XFCN (dgesdd, DGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1), | |
|
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
189 m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, |
|
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
190 work.fortran_vec (), lwork, iwork, info |
| 10601 | 191 F77_CHAR_ARG_LEN (1))); |
| 192 | |
| 193 } | |
| 194 else | |
| 195 assert (0); // impossible | |
| 3336 | 196 |
|
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
197 if (! (jobv == 'N' || jobv == 'O')) |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
198 right_sm = right_sm.transpose (); |
| 457 | 199 |
| 200 return info; | |
| 201 } | |
| 202 | |
| 3504 | 203 std::ostream& |
| 204 operator << (std::ostream& os, const SVD& a) | |
| 457 | 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 } |