Mercurial > octave

comparison liboctave/numeric/gsvd.cc @ 22236:065a44375723

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
gsvd: reduce code duplication with templates. * CmplxGSVD.cc, CmplxGSVD.h, dbleGSVD.cc, dbleGSVD.h: Remove files for no longer existing classes. Replaced by gsvd template class. This classes never existed in an Octave release, this was freshly imported from Octave Forge so backwards compatibility is not an issue. * liboctave/numeric/gsvd.h, liboctave/numeric/gsvd.cc: New files for gsvd class template generated from CmplxGSVD.cc, CmplxGSVD.h, dbleGSVD.cc, and dbleGSVD.h and converted to template. Removed unused << operator, unused constructor with &info, and commented code. Only instantiated for Matrix and ComplexMatrix, providing interface to DGGSVD and ZGGSVD. * liboctave/numeric/module.mk: Update. * mx-defs.h, mx-ext.h: Use new classes.
author Barbara Locsi <locsi.barbara@gmail.com>
date Tue, 09 Aug 2016 18:02:11 +0200
parents
children da201af35c97
comparison
equal deleted inserted replaced
22235:63b41167ef1e 22236:065a44375723
1 // Copyright (C) 1996, 1997 John W. Eaton
2 // Copyright (C) 2006 Pascal Dupuis <Pascal.Dupuis@uclouvain.be>
3 // Copyright (C) 2016 Barbara Lócsi
4 //
5 // This program is free software; you can redistribute it and/or modify it under
6 // the terms of the GNU General Public License as published by the Free Software
7 // Foundation; either version 3 of the License, or (at your option) any later
8 // version.
9 //
10 // This program is distributed in the hope that it will be useful, but WITHOUT
11 // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 // FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 // details.
14 //
15 // You should have received a copy of the GNU General Public License along with
16 // this program; if not, see <http://www.gnu.org/licenses/>.
17
18 #ifdef HAVE_CONFIG_H
19 # include <config.h>
20 #endif
21
22 #include "gsvd.h"
23 #include "f77-fcn.h"
24 #include "lo-error.h"
25 #include "CMatrix.h"
26 #include "dDiagMatrix.h"
27 #include "dMatrix.h"
28
29 extern "C"
30 {
31 F77_RET_T
32 F77_FUNC (dggsvd, DGGSVD)
33 (
34 F77_CONST_CHAR_ARG_DECL, // JOBU (input) CHARACTER*1
35 F77_CONST_CHAR_ARG_DECL, // JOBV (input) CHARACTER*1
36 F77_CONST_CHAR_ARG_DECL, // JOBQ (input) CHARACTER*1
37 const F77_INT&, // M (input) INTEGER
38 const F77_INT&, // N (input) INTEGER
39 const F77_INT&, // P (input) INTEGER
40 F77_INT &, // K (output) INTEGER
41 F77_INT &, // L (output) INTEGER
42 F77_DBLE*, // A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
43 const F77_INT&, // LDA (input) INTEGER
44 F77_DBLE*, // B (input/output) DOUBLE PRECISION array, dimension (LDB,N)
45 const F77_INT&, // LDB (input) INTEGER
46 F77_DBLE*, // ALPHA (output) DOUBLE PRECISION array, dimension (N)
47 F77_DBLE*, // BETA (output) DOUBLE PRECISION array, dimension (N)
48 F77_DBLE*, // U (output) DOUBLE PRECISION array, dimension (LDU,M)
49 const F77_INT&, // LDU (input) INTEGER
50 F77_DBLE*, // V (output) DOUBLE PRECISION array, dimension (LDV,P)
51 const F77_INT&, // LDV (input) INTEGER
52 F77_DBLE*, // Q (output) DOUBLE PRECISION array, dimension (LDQ,N)
53 const F77_INT&, // LDQ (input) INTEGER
54 F77_DBLE*, // WORK (workspace) DOUBLE PRECISION array
55 int*, // IWORK (workspace/output) INTEGER array, dimension (N)
56 F77_INT& // INFO (output)INTEGER
57 F77_CHAR_ARG_LEN_DECL
58 F77_CHAR_ARG_LEN_DECL
59 F77_CHAR_ARG_LEN_DECL
60 );
61
62 F77_RET_T
63 F77_FUNC (zggsvd, ZGGSVD)
64 (
65 F77_CONST_CHAR_ARG_DECL, // JOBU (input) CHARACTER*1
66 F77_CONST_CHAR_ARG_DECL, // JOBV (input) CHARACTER*1
67 F77_CONST_CHAR_ARG_DECL, // JOBQ (input) CHARACTER*1
68 const F77_INT&, // M (input) INTEGER
69 const F77_INT&, // N (input) INTEGER
70 const F77_INT&, // P (input) INTEGER
71 F77_INT &, // K (output) INTEGER
72 F77_INT &, // L (output) INTEGER
73 F77_DBLE_CMPLX*, // A (input/output) COMPLEX*16 array, dimension (LDA,N)
74 const F77_INT&, // LDA (input) INTEGER
75 F77_DBLE_CMPLX*, // B (input/output) COMPLEX*16 array, dimension (LDB,N)
76 const F77_INT&, // LDB (input) INTEGER
77 F77_DBLE*, // ALPHA (output) DOUBLE PRECISION array, dimension (N)
78 F77_DBLE*, // BETA (output) DOUBLE PRECISION array, dimension (N)
79 F77_DBLE_CMPLX*, // U (output) COMPLEX*16 array, dimension (LDU,M)
80 const F77_INT&, // LDU (input) INTEGER
81 F77_DBLE_CMPLX*, // V (output) COMPLEX*16 array, dimension (LDV,P)
82 const F77_INT&, // LDV (input) INTEGER
83 F77_DBLE_CMPLX*, // Q (output) COMPLEX*16 array, dimension (LDQ,N)
84 const F77_INT&, // LDQ (input) INTEGER
85 F77_DBLE_CMPLX*, // WORK (workspace) COMPLEX*16 array
86 F77_DBLE*, // RWORK (workspace) DOUBLE PRECISION array
87 int*, // IWORK (workspace/output) INTEGER array, dimension (N)
88 F77_INT& // INFO (output)INTEGER
89 F77_CHAR_ARG_LEN_DECL
90 F77_CHAR_ARG_LEN_DECL
91 F77_CHAR_ARG_LEN_DECL
92 );
93 }
94
95 template <>
96 void
97 gsvd<Matrix>::ggsvd (char& jobu, char& jobv, char& jobq, octave_idx_type m,
98 octave_idx_type n, octave_idx_type p, octave_idx_type& k,
99 octave_idx_type& l, double *tmp_dataA, octave_idx_type m1,
100 double *tmp_dataB, octave_idx_type p1, Matrix& alpha,
101 Matrix& beta, double *u, octave_idx_type nrow_u, double *v,
102 octave_idx_type nrow_v, double *q, octave_idx_type nrow_q,
103 Matrix& work, octave_idx_type* iwork, octave_idx_type& info)
104 {
105 F77_XFCN (dggsvd, DGGSVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
106 F77_CONST_CHAR_ARG2 (&jobv, 1),
107 F77_CONST_CHAR_ARG2 (&jobq, 1),
108 m, n, p, k, l, tmp_dataA, m1,
109 tmp_dataB, p1, alpha.fortran_vec (),
110 beta.fortran_vec (), u, nrow_u,
111 v, nrow_v, q, nrow_q, work.fortran_vec (),
112 iwork, info
113 F77_CHAR_ARG_LEN (1)
114 F77_CHAR_ARG_LEN (1)
115 F77_CHAR_ARG_LEN (1)));
116 }
117
118
119 template <>
120 void
121 gsvd<ComplexMatrix>::ggsvd (char& jobu, char& jobv, char& jobq,
122 octave_idx_type m, octave_idx_type n,
123 octave_idx_type p, octave_idx_type& k,
124 octave_idx_type& l, Complex *tmp_dataA,
125 octave_idx_type m1, Complex *tmp_dataB,
126 octave_idx_type p1, Matrix& alpha,
127 Matrix& beta, Complex *u, octave_idx_type nrow_u,
128 Complex *v, octave_idx_type nrow_v, Complex *q,
129 octave_idx_type nrow_q, ComplexMatrix& work,
130 octave_idx_type* iwork, octave_idx_type& info)
131 {
132 Matrix rwork(2*n, 1);
133 F77_XFCN (zggsvd, ZGGSVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
134 F77_CONST_CHAR_ARG2 (&jobv, 1),
135 F77_CONST_CHAR_ARG2 (&jobq, 1),
136 m, n, p, k, l, F77_DBLE_CMPLX_ARG (tmp_dataA), m1,
137 F77_DBLE_CMPLX_ARG (tmp_dataB), p1,
138 alpha.fortran_vec (), beta.fortran_vec (),
139 F77_DBLE_CMPLX_ARG (u), nrow_u,
140 F77_DBLE_CMPLX_ARG (v), nrow_v,
141 F77_DBLE_CMPLX_ARG (q), nrow_q,
142 F77_DBLE_CMPLX_ARG (work.fortran_vec ()),
143 rwork.fortran_vec (), iwork, info
144 F77_CHAR_ARG_LEN (1)
145 F77_CHAR_ARG_LEN (1)
146 F77_CHAR_ARG_LEN (1)));
147 }
148
149 template <typename T>
150 T
151 gsvd<T>::left_singular_matrix_A (void) const
152 {
153 if (type == gsvd::Type::sigma_only)
154 {
155 (*current_liboctave_error_handler)
156 ("gsvd: U not computed because type == gsvd::sigma_only");
157 return T ();
158 }
159 else
160 return left_smA;
161 }
162
163 template <typename T>
164 T
165 gsvd<T>::left_singular_matrix_B (void) const
166 {
167 if (type == gsvd::Type::sigma_only)
168 {
169 (*current_liboctave_error_handler)
170 ("gsvd: V not computed because type == gsvd::sigma_only");
171 return T ();
172 }
173 else
174 return left_smB;
175 }
176
177 template <typename T>
178 T
179 gsvd<T>::right_singular_matrix (void) const
180 {
181 if (type == gsvd::Type::sigma_only)
182 {
183 (*current_liboctave_error_handler)
184 ("gsvd: X not computed because type == gsvd::sigma_only");
185 return T ();
186 }
187 else
188 return right_sm;
189 }
190
191 template <typename T>
192 T
193 gsvd<T>::R_matrix (void) const
194 {
195 if (type != gsvd::Type::std)
196 {
197 (*current_liboctave_error_handler)
198 ("gsvd: R not computed because type != gsvd::std");
199 return T ();
200 }
201 else
202 return R;
203 }
204
205 template <typename T>
206 gsvd<T>::gsvd (const T& a, const T& b, gsvd::Type gsvd_type)
207 {
208 octave_idx_type info;
209
210 octave_idx_type m = a.rows ();
211 octave_idx_type n = a.cols ();
212 octave_idx_type p = b.rows ();
213
214 T atmp = a;
215 P *tmp_dataA = atmp.fortran_vec ();
216
217 T btmp = b;
218 P *tmp_dataB = btmp.fortran_vec ();
219
220 char jobu = 'U';
221 char jobv = 'V';
222 char jobq = 'Q';
223
224 octave_idx_type nrow_u = m;
225 octave_idx_type nrow_v = p;
226 octave_idx_type nrow_q = n;
227
228 octave_idx_type k, l;
229
230 switch (gsvd_type)
231 {
232 case gsvd<T>::Type::sigma_only:
233
234 // Note: for this case, both jobu and jobv should be 'N', but
235 // there seems to be a bug in dgesvd from Lapack V2.0. To
236 // demonstrate the bug, set both jobu and jobv to 'N' and find
237 // the singular values of [eye(3), eye(3)]. The result is
238 // [-sqrt(2), -sqrt(2), -sqrt(2)].
239 //
240 // For Lapack 3.0, this problem seems to be fixed.
241
242 jobu = 'N';
243 jobv = 'N';
244 jobq = 'N';
245 nrow_u = nrow_v = nrow_q = 1;
246 break;
247
248 default:
249 break;
250 }
251
252 type = gsvd_type;
253
254 if (! (jobu == 'N' || jobu == 'O'))
255 left_smA.resize (nrow_u, m);
256
257 P *u = left_smA.fortran_vec ();
258
259 if (! (jobv == 'N' || jobv == 'O'))
260 left_smB.resize (nrow_v, p);
261
262 P *v = left_smB.fortran_vec ();
263
264 if (! (jobq == 'N' || jobq == 'O'))
265 right_sm.resize (nrow_q, n);
266
267 P *q = right_sm.fortran_vec ();
268
269 octave_idx_type lwork = 3*n;
270 lwork = lwork > m ? lwork : m;
271 lwork = (lwork > p ? lwork : p) + n;
272
273 T work (lwork, 1);
274 Matrix alpha (n, 1);
275 Matrix beta (n, 1);
276
277 std::vector<octave_idx_type> iwork (n);
278
279 gsvd<T>::ggsvd (jobu, jobv, jobq, m, n, p, k, l,
280 tmp_dataA, m, tmp_dataB, p, alpha, beta, u,
281 nrow_u, v, nrow_v, q, nrow_q, work, iwork.data (), info);
282
283 if (f77_exception_encountered)
284 (*current_liboctave_error_handler) ("unrecoverable error in *ggsvd");
285
286 if (info < 0)
287 (*current_liboctave_error_handler) ("*ggsvd.f: argument %d illegal", -info);
288 else
289 {
290 if (info > 0)
291 (*current_liboctave_error_handler)
292 ("*ggsvd.f: Jacobi-type procedure failed to converge.");
293 else
294 {
295 octave_idx_type i, j;
296
297 if (gsvd<T>::Type::std == gsvd_type)
298 {
299 R.resize(k+l, k+l);
300 int astart = n-k-l;
301 if (m - k - l >= 0)
302 {
303 astart = n-k-l;
304 // R is stored in A(1:K+L,N-K-L+1:N)
305 for (i = 0; i < k+l; i++)
306 for (j = 0; j < k+l; j++)
307 R.xelem (i, j) = atmp.xelem (i, astart + j);
308 }
309 else
310 {
311 // (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N),
312 // ( 0 R22 R23 )
313
314 for (i = 0; i < m; i++)
315 for (j = 0; j < k+l; j++)
316 R.xelem (i, j) = atmp.xelem (i, astart + j);
317 // and R33 is stored in B(M-K+1:L,N+M-K-L+1:N)
318 for (i = k+l-1; i >=m; i--) {
319 for (j = 0; j < m; j++)
320 R.xelem(i, j) = 0.0;
321 for (j = m; j < k+l; j++)
322 R.xelem (i, j) = btmp.xelem (i - k, astart + j);
323 }
324 }
325 }
326
327 if (m-k-l >= 0)
328 {
329 // Fills in C and S
330 sigmaA.resize (l, l);
331 sigmaB.resize (l, l);
332 for (i = 0; i < l; i++)
333 {
334 sigmaA.dgxelem(i) = alpha.elem(k+i);
335 sigmaB.dgxelem(i) = beta.elem(k+i);
336 }
337 }
338 else
339 {
340 // Fills in C and S
341 sigmaA.resize (m-k, m-k);
342 sigmaB.resize (m-k, m-k);
343 for (i = 0; i < m-k; i++)
344 {
345 sigmaA.dgxelem(i) = alpha.elem(k+i);
346 sigmaB.dgxelem(i) = beta.elem(k+i);
347 }
348 }
349 }
350 }
351 }
352
353 // Instantiations we need.
354
355 template class gsvd<Matrix>;
356
357 template class gsvd<ComplexMatrix>;