Mercurial > octave

comparison liboctave/numeric/CmplxGSVD.cc @ 22235:63b41167ef1e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
gsvd: new function imported from Octave-Forge linear-algebra package. * libinterp/corefcn/gsvd.cc: New function to the interpreter. Imported from the linear-algebra package. * CmplxGSVD.cc, CmplxGSVD.h, dbleGSVD.cc, dbleGSVD.h: new classes imported from the linear-algebra package to compute gsvd of Matrix and ComplexMatrix. * liboctave/operators/mx-defs.h, liboctave/operators/mx-ext.h: use new classes. * libinterp/corefcn/module.mk, liboctave/numeric/module.mk: Add to the * scripts/help/__unimplemented__.m: Remove "gsvd" from list. * doc/interpreter/linalg.txi: Add to manual. build system. * NEWS: Add function to list of new functions for 4.2.
author Barbara Locsi <locsi.barbara@gmail.com>
date Thu, 04 Aug 2016 07:50:31 +0200
parents
children
comparison
equal deleted inserted replaced
22234:66dd260512a4 22235:63b41167ef1e
1 // Copyright (C) 1996, 1997 John W. Eaton
2 // Copyright (C) 2006 Pascal Dupuis <Pascal.Dupuis@uclouvain.be>
3 //
4 // This program is free software; you can redistribute it and/or modify it under
5 // the terms of the GNU General Public License as published by the Free Software
6 // Foundation; either version 3 of the License, or (at your option) any later
7 // version.
8 //
9 // This program is distributed in the hope that it will be useful, but WITHOUT
10 // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
11 // FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
12 // details.
13 //
14 // You should have received a copy of the GNU General Public License along with
15 // this program; if not, see <http://www.gnu.org/licenses/>.
16
17 #ifdef HAVE_CONFIG_H
18 # include <config.h>
19 #endif
20
21 #include "CmplxGSVD.h"
22 #include "f77-fcn.h"
23 #include "lo-error.h"
24
25
26 /*
27 uncomment those lines to monitor k and l
28 #include "oct-obj.h"
29 #include "pager.h"
30 */
31
32 extern "C"
33 {
34 F77_RET_T
35 F77_FUNC (zggsvd, ZGGSVD)
36 (
37 F77_CONST_CHAR_ARG_DECL, // JOBU (input) CHARACTER*1
38 F77_CONST_CHAR_ARG_DECL, // JOBV (input) CHARACTER*1
39 F77_CONST_CHAR_ARG_DECL, // JOBQ (input) CHARACTER*1
40 const octave_idx_type&, // M (input) INTEGER
41 const octave_idx_type&, // N (input) INTEGER
42 const octave_idx_type&, // P (input) INTEGER
43 octave_idx_type &, // K (output) INTEGER
44 octave_idx_type &, // L (output) INTEGER
45 Complex*, // A (input/output) COMPLEX*16 array, dimension (LDA,N)
46 const octave_idx_type&, // LDA (input) INTEGER
47 Complex*, // B (input/output) COMPLEX*16 array, dimension (LDB,N)
48 const octave_idx_type&, // LDB (input) INTEGER
49 double*, // ALPHA (output) DOUBLE PRECISION array, dimension (N)
50 double*, // BETA (output) DOUBLE PRECISION array, dimension (N)
51 Complex*, // U (output) COMPLEX*16 array, dimension (LDU,M)
52 const octave_idx_type&, // LDU (input) INTEGER
53 Complex*, // V (output) COMPLEX*16 array, dimension (LDV,P)
54 const octave_idx_type&, // LDV (input) INTEGER
55 Complex*, // Q (output) COMPLEX*16 array, dimension (LDQ,N)
56 const octave_idx_type&, // LDQ (input) INTEGER
57 Complex*, // WORK (workspace) COMPLEX*16 array
58 double*, // RWORK (workspace) DOUBLE PRECISION array
59 int*, // IWORK (workspace/output) INTEGER array, dimension (N)
60 octave_idx_type& // INFO (output)INTEGER
61 F77_CHAR_ARG_LEN_DECL
62 F77_CHAR_ARG_LEN_DECL
63 F77_CHAR_ARG_LEN_DECL
64 );
65 }
66
67 ComplexMatrix
68 ComplexGSVD::left_singular_matrix_A (void) const
69 {
70 if (type_computed == GSVD::sigma_only)
71 {
72 (*current_liboctave_error_handler)
73 ("CmplxGSVD: U not computed because type == GSVD::sigma_only");
74 return ComplexMatrix ();
75 }
76 else
77 return left_smA;
78 }
79
80 ComplexMatrix
81 ComplexGSVD::left_singular_matrix_B (void) const
82 {
83 if (type_computed == GSVD::sigma_only)
84 {
85 (*current_liboctave_error_handler)
86 ("CmplxGSVD: V not computed because type == GSVD::sigma_only");
87 return ComplexMatrix ();
88 }
89 else
90 return left_smB;
91 }
92
93 ComplexMatrix
94 ComplexGSVD::right_singular_matrix (void) const
95 {
96 if (type_computed == GSVD::sigma_only)
97 {
98 (*current_liboctave_error_handler)
99 ("CmplxGSVD: X not computed because type == GSVD::sigma_only");
100 return ComplexMatrix ();
101 }
102 else
103 return right_sm;
104 }
105
106 ComplexMatrix
107 ComplexGSVD::R_matrix (void) const
108 {
109 if (type_computed != GSVD::std)
110 {
111 (*current_liboctave_error_handler)
112 ("CmplxGSVD: R not computed because type != GSVD::std");
113 return ComplexMatrix ();
114 }
115 else
116 return R;
117 }
118
119 octave_idx_type
120 ComplexGSVD::init (const ComplexMatrix& a, const ComplexMatrix& b,
121 GSVD::type gsvd_type)
122 {
123 octave_idx_type info;
124
125 octave_idx_type m = a.rows ();
126 octave_idx_type n = a.cols ();
127 octave_idx_type p = b.rows ();
128
129 ComplexMatrix atmp = a;
130 Complex *tmp_dataA = atmp.fortran_vec ();
131
132 ComplexMatrix btmp = b;
133 Complex *tmp_dataB = btmp.fortran_vec ();
134
135 // octave_idx_type min_mn = m < n ? m : n;
136
137 char jobu = 'U';
138 char jobv = 'V';
139 char jobq = 'Q';
140
141 octave_idx_type nrow_u = m;
142 octave_idx_type nrow_v = p;
143 octave_idx_type nrow_q = n;
144
145 octave_idx_type k, l;
146
147 switch (gsvd_type)
148 {
149
150 case GSVD::sigma_only:
151
152 // Note: for this case, both jobu and jobv should be 'N', but
153 // there seems to be a bug in dgesvd from Lapack V2.0. To
154 // demonstrate the bug, set both jobu and jobv to 'N' and find
155 // the singular values of [eye(3), eye(3)]. The result is
156 // [-sqrt(2), -sqrt(2), -sqrt(2)].
157 //
158 // For Lapack 3.0, this problem seems to be fixed.
159
160 jobu = 'N';
161 jobv = 'N';
162 jobq = 'N';
163 nrow_u = nrow_v = nrow_q = 1;
164 break;
165
166 default:
167 break;
168 }
169
170 type_computed = gsvd_type;
171
172 if (! (jobu == 'N' || jobu == 'O')) {
173 left_smA.resize (nrow_u, m);
174 }
175
176 Complex *u = left_smA.fortran_vec ();
177
178 if (! (jobv == 'N' || jobv == 'O')) {
179 left_smB.resize (nrow_v, p);
180 }
181
182 Complex *v = left_smB.fortran_vec ();
183
184 if (! (jobq == 'N' || jobq == 'O')) {
185 right_sm.resize (nrow_q, n);
186 }
187 Complex *q = right_sm.fortran_vec ();
188
189 octave_idx_type lwork = 3*n;
190 lwork = lwork > m ? lwork : m;
191 lwork = (lwork > p ? lwork : p) + n;
192
193 Array<Complex> work (dim_vector (lwork, 1));
194 Array<double> alpha (dim_vector (n, 1));
195 Array<double> beta (dim_vector (n, 1));
196 Array<double> rwork(dim_vector (2*n, 1));
197 Array<int> iwork (dim_vector (n, 1));
198
199 F77_XFCN (zggsvd, ZGGSVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
200 F77_CONST_CHAR_ARG2 (&jobv, 1),
201 F77_CONST_CHAR_ARG2 (&jobq, 1),
202 m, n, p, k, l, tmp_dataA, m,
203 tmp_dataB, p, alpha.fortran_vec (),
204 beta.fortran_vec (), u, nrow_u,
205 v, nrow_v, q, nrow_q, work.fortran_vec (),
206 rwork.fortran_vec (), iwork.fortran_vec (),
207 info
208 F77_CHAR_ARG_LEN (1)
209 F77_CHAR_ARG_LEN (1)
210 F77_CHAR_ARG_LEN (1)));
211
212 if (f77_exception_encountered)
213 (*current_liboctave_error_handler) ("unrecoverable error in zggsvd");
214
215 if (info < 0) {
216 (*current_liboctave_error_handler) ("zggsvd.f: argument %d illegal", -info);
217 } else {
218 if (info > 0) {
219 (*current_liboctave_error_handler) ("zggsvd.f: Jacobi-type procedure failed to converge.");
220 } else {
221 octave_idx_type i, j;
222
223 if (GSVD::std == gsvd_type) {
224 R.resize(k+l, k+l);
225 int astart = n-k-l;
226 if (m - k - l >= 0) {
227 int astart = n-k-l;
228 /*
229 * R is stored in A(1:K+L,N-K-L+1:N)
230 */
231 for (i = 0; i < k+l; i++)
232 for (j = 0; j < k+l; j++)
233 R.xelem(i, j) = atmp.xelem(i, astart + j);
234 } else {
235 /*
236 * (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N),
237 * ( 0 R22 R23 )
238 */
239
240 for (i = 0; i < m; i++)
241 for (j = 0; j < k+l; j++)
242 R.xelem(i, j) = atmp.xelem(i, astart + j);
243 /*
244 * and R33 is stored in B(M-K+1:L,N+M-K-L+1:N)
245 */
246 for (i = k+l-1; i >=m; i--) {
247 for (j = 0; j < m; j++)
248 R.xelem(i, j) = 0.0;
249 for (j = m; j < k+l; j++)
250 R.xelem(i, j) = btmp.xelem(i - k, astart + j);
251 }
252 }
253 }
254 /*
255 uncomment this to monitor k and l
256 octave_value tmp;
257 octave_stdout << "CmplxGSVD k: ";
258 tmp = k;
259 tmp.print(octave_stdout);
260 octave_stdout << "\n";
261 octave_stdout << "CmplxGSVD l: ";
262 tmp = l;
263 tmp.print(octave_stdout);
264 octave_stdout << "\n";
265 */
266
267 if (m-k-l >= 0) {
268 // Fills in C and S
269 sigmaA.resize (l, l);
270 sigmaB.resize (l, l);
271 for (i = 0; i < l; i++) {
272 sigmaA.dgxelem(i) = alpha.elem(k+i);
273 sigmaB.dgxelem(i) = beta.elem(k+i);
274 }
275 } else {
276 // Fills in C and S
277 sigmaA.resize (m-k, m-k);
278 sigmaB.resize (m-k, m-k);
279 for (i = 0; i < m-k; i++) {
280 sigmaA.dgxelem(i) = alpha.elem(k+i);
281 sigmaB.dgxelem(i) = beta.elem(k+i);
282 }
283 }
284 }
285 }
286 return info;
287 }