Mercurial > hg > octave-nkf
annotate liboctave/numeric/CmplxSVD.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 | |
| 27 #include "CmplxSVD.h" | |
| 1847 | 28 #include "f77-fcn.h" |
| 1543 | 29 #include "lo-error.h" |
| 10601 | 30 #include "oct-locbuf.h" |
| 457 | 31 |
| 32 extern "C" | |
| 33 { | |
| 4552 | 34 F77_RET_T |
| 35 F77_FUNC (zgesvd, ZGESVD) (F77_CONST_CHAR_ARG_DECL, | |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10258
diff
changeset
|
36 F77_CONST_CHAR_ARG_DECL, |
| 11495 | 37 const octave_idx_type&, const octave_idx_type&, |
| 38 Complex*, const octave_idx_type&, | |
| 39 double*, Complex*, const octave_idx_type&, | |
| 40 Complex*, const octave_idx_type&, Complex*, | |
| 41 const octave_idx_type&, double*, octave_idx_type& | |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10258
diff
changeset
|
42 F77_CHAR_ARG_LEN_DECL |
|
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10258
diff
changeset
|
43 F77_CHAR_ARG_LEN_DECL); |
| 10601 | 44 |
| 45 F77_RET_T | |
| 46 F77_FUNC (zgesdd, ZGESDD) (F77_CONST_CHAR_ARG_DECL, | |
| 11495 | 47 const octave_idx_type&, const octave_idx_type&, |
| 48 Complex*, const octave_idx_type&, | |
| 49 double*, Complex*, const octave_idx_type&, | |
| 50 Complex*, const octave_idx_type&, Complex*, | |
| 51 const octave_idx_type&, double*, | |
| 52 octave_idx_type *, octave_idx_type& | |
| 10601 | 53 F77_CHAR_ARG_LEN_DECL); |
| 457 | 54 } |
| 55 | |
| 1543 | 56 ComplexMatrix |
| 57 ComplexSVD::left_singular_matrix (void) const | |
| 58 { | |
| 1544 | 59 if (type_computed == SVD::sigma_only) |
| 1543 | 60 { |
| 61 (*current_liboctave_error_handler) | |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10258
diff
changeset
|
62 ("ComplexSVD: U not computed because type == SVD::sigma_only"); |
| 1543 | 63 return ComplexMatrix (); |
| 64 } | |
| 65 else | |
| 66 return left_sm; | |
| 67 } | |
| 68 | |
| 69 ComplexMatrix | |
| 70 ComplexSVD::right_singular_matrix (void) const | |
| 71 { | |
| 1544 | 72 if (type_computed == SVD::sigma_only) |
| 1543 | 73 { |
| 74 (*current_liboctave_error_handler) | |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10258
diff
changeset
|
75 ("ComplexSVD: V not computed because type == SVD::sigma_only"); |
| 1543 | 76 return ComplexMatrix (); |
| 77 } | |
| 78 else | |
| 79 return right_sm; | |
| 80 } | |
| 81 | |
| 5275 | 82 octave_idx_type |
|
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
83 ComplexSVD::init (const ComplexMatrix& a, SVD::type svd_type, |
|
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
84 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 |
| 1946 | 91 ComplexMatrix atmp = a; |
| 92 Complex *tmp_data = atmp.fortran_vec (); | |
| 457 | 93 |
| 5275 | 94 octave_idx_type min_mn = m < n ? m : n; |
| 95 octave_idx_type max_mn = m > n ? m : n; | |
| 457 | 96 |
| 1930 | 97 char jobu = 'A'; |
| 98 char jobv = 'A'; | |
| 537 | 99 |
| 5275 | 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; | |
| 537 | 104 |
| 1543 | 105 switch (svd_type) |
| 537 | 106 { |
| 1543 | 107 case SVD::economy: |
| 1930 | 108 jobu = jobv = 'S'; |
| 537 | 109 ncol_u = nrow_vt = nrow_s = ncol_s = min_mn; |
| 1543 | 110 break; |
| 111 | |
| 112 case SVD::sigma_only: | |
| 2621 | 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)]. | |
| 3335 | 119 // |
| 120 // For Lapack 3.0, this problem seems to be fixed. | |
| 2621 | 121 |
|
15887
8ced82e96b48
Fix segfaults with gesdd driver for svd (bug #37998).
Rik <rik@octave.org>
parents:
14138
diff
changeset
|
122 jobu = jobv = 'N'; |
| 1545 | 123 ncol_u = nrow_vt = 1; |
| 1543 | 124 break; |
| 125 | |
| 126 default: | |
| 127 break; | |
| 537 | 128 } |
| 129 | |
| 1544 | 130 type_computed = svd_type; |
| 131 | |
| 2621 | 132 if (! (jobu == 'N' || jobu == 'O')) |
| 1930 | 133 left_sm.resize (m, ncol_u); |
| 134 | |
| 135 Complex *u = left_sm.fortran_vec (); | |
| 136 | |
| 137 sigma.resize (nrow_s, ncol_s); | |
| 138 double *s_vec = sigma.fortran_vec (); | |
| 139 | |
| 2621 | 140 if (! (jobv == 'N' || jobv == 'O')) |
| 1930 | 141 right_sm.resize (nrow_vt, n); |
| 142 | |
| 143 Complex *vt = right_sm.fortran_vec (); | |
| 457 | 144 |
|
15887
8ced82e96b48
Fix segfaults with gesdd driver for svd (bug #37998).
Rik <rik@octave.org>
parents:
14138
diff
changeset
|
145 // Query ZGESVD for the correct dimension of WORK. |
| 3336 | 146 |
| 5275 | 147 octave_idx_type lwork = -1; |
| 457 | 148 |
|
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
149 Array<Complex> work (dim_vector (1, 1)); |
| 3336 | 150 |
| 10258 | 151 octave_idx_type one = 1; |
|
15887
8ced82e96b48
Fix segfaults with gesdd driver for svd (bug #37998).
Rik <rik@octave.org>
parents:
14138
diff
changeset
|
152 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
|
153 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
|
154 |
| 10601 | 155 if (svd_driver == SVD::GESVD) |
| 156 { | |
|
15887
8ced82e96b48
Fix segfaults with gesdd driver for svd (bug #37998).
Rik <rik@octave.org>
parents:
14138
diff
changeset
|
157 octave_idx_type lrwork = 5*max_mn; |
|
8ced82e96b48
Fix segfaults with gesdd driver for svd (bug #37998).
Rik <rik@octave.org>
parents:
14138
diff
changeset
|
158 Array<double> rwork (dim_vector (lrwork, 1)); |
|
8ced82e96b48
Fix segfaults with gesdd driver for svd (bug #37998).
Rik <rik@octave.org>
parents:
14138
diff
changeset
|
159 |
| 10601 | 160 F77_XFCN (zgesvd, ZGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), |
| 161 F77_CONST_CHAR_ARG2 (&jobv, 1), | |
| 162 m, n, tmp_data, m1, s_vec, u, m1, vt, | |
| 163 nrow_vt1, work.fortran_vec (), lwork, | |
| 164 rwork.fortran_vec (), info | |
| 165 F77_CHAR_ARG_LEN (1) | |
| 166 F77_CHAR_ARG_LEN (1))); | |
| 167 | |
| 168 lwork = static_cast<octave_idx_type> (work(0).real ()); | |
|
11574
a83bad07f7e3
attempt better backward compatibility for Array resize functions
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
169 work.resize (dim_vector (lwork, 1)); |
| 1543 | 170 |
| 10601 | 171 F77_XFCN (zgesvd, ZGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), |
| 172 F77_CONST_CHAR_ARG2 (&jobv, 1), | |
| 173 m, n, tmp_data, m1, s_vec, u, m1, vt, | |
| 174 nrow_vt1, work.fortran_vec (), lwork, | |
| 175 rwork.fortran_vec (), info | |
| 176 F77_CHAR_ARG_LEN (1) | |
| 177 F77_CHAR_ARG_LEN (1))); | |
| 178 } | |
| 179 else if (svd_driver == SVD::GESDD) | |
| 180 { | |
| 181 assert (jobu == jobv); | |
| 182 char jobz = jobu; | |
|
15887
8ced82e96b48
Fix segfaults with gesdd driver for svd (bug #37998).
Rik <rik@octave.org>
parents:
14138
diff
changeset
|
183 |
|
8ced82e96b48
Fix segfaults with gesdd driver for svd (bug #37998).
Rik <rik@octave.org>
parents:
14138
diff
changeset
|
184 octave_idx_type lrwork; |
|
8ced82e96b48
Fix segfaults with gesdd driver for svd (bug #37998).
Rik <rik@octave.org>
parents:
14138
diff
changeset
|
185 if (jobz == 'N') |
|
8ced82e96b48
Fix segfaults with gesdd driver for svd (bug #37998).
Rik <rik@octave.org>
parents:
14138
diff
changeset
|
186 lrwork = 7*min_mn; |
|
8ced82e96b48
Fix segfaults with gesdd driver for svd (bug #37998).
Rik <rik@octave.org>
parents:
14138
diff
changeset
|
187 else |
|
8ced82e96b48
Fix segfaults with gesdd driver for svd (bug #37998).
Rik <rik@octave.org>
parents:
14138
diff
changeset
|
188 lrwork = 5*min_mn*min_mn + 5*min_mn; |
|
8ced82e96b48
Fix segfaults with gesdd driver for svd (bug #37998).
Rik <rik@octave.org>
parents:
14138
diff
changeset
|
189 Array<double> rwork (dim_vector (lrwork, 1)); |
|
8ced82e96b48
Fix segfaults with gesdd driver for svd (bug #37998).
Rik <rik@octave.org>
parents:
14138
diff
changeset
|
190 |
| 10601 | 191 OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, 8*min_mn); |
| 3336 | 192 |
| 10601 | 193 F77_XFCN (zgesdd, ZGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1), |
| 194 m, n, tmp_data, m1, s_vec, u, m1, vt, | |
| 195 nrow_vt1, work.fortran_vec (), lwork, | |
| 196 rwork.fortran_vec (), iwork, info | |
| 197 F77_CHAR_ARG_LEN (1))); | |
| 198 | |
| 199 lwork = static_cast<octave_idx_type> (work(0).real ()); | |
|
11574
a83bad07f7e3
attempt better backward compatibility for Array resize functions
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
200 work.resize (dim_vector (lwork, 1)); |
| 10601 | 201 |
| 202 F77_XFCN (zgesdd, ZGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1), | |
| 203 m, n, tmp_data, m1, s_vec, u, m1, vt, | |
| 204 nrow_vt1, work.fortran_vec (), lwork, | |
| 205 rwork.fortran_vec (), iwork, info | |
| 206 F77_CHAR_ARG_LEN (1))); | |
| 207 } | |
| 208 else | |
| 209 assert (0); // impossible | |
| 3336 | 210 |
|
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
211 if (! (jobv == 'N' || jobv == 'O')) |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
212 right_sm = right_sm.hermitian (); |
| 457 | 213 |
| 214 return info; | |
| 215 } |