GNU Octave 7.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
gsvd.cc
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////
2//
3// Copyright (C) 1997-2022 The Octave Project Developers
4//
5// See the file COPYRIGHT.md in the top-level directory of this
6// distribution or <https://octave.org/copyright/>.
7//
8// This file is part of Octave.
9//
10// Octave is free software: you can redistribute it and/or modify it
11// under the terms of the GNU General Public License as published by
12// the Free Software Foundation, either version 3 of the License, or
13// (at your option) any later version.
14//
15// Octave is distributed in the hope that it will be useful, but
16// WITHOUT ANY WARRANTY; without even the implied warranty of
17// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18// GNU General Public License for more details.
19//
20// You should have received a copy of the GNU General Public License
21// along with Octave; see the file COPYING. If not, see
22// <https://www.gnu.org/licenses/>.
23//
24////////////////////////////////////////////////////////////////////////
25
26#ifdef HAVE_CONFIG_H
27# include <config.h>
28#endif
29
30#include <algorithm>
31#include <unordered_map>
32
33#include "CMatrix.h"
34#include "dDiagMatrix.h"
35#include "dMatrix.h"
36#include "fCMatrix.h"
37#include "fDiagMatrix.h"
38#include "fMatrix.h"
39#include "gsvd.h"
40#include "lo-error.h"
41#include "lo-lapack-proto.h"
42#include "oct-locbuf.h"
43#include "oct-shlib.h"
44
45namespace octave
46{
47 static std::unordered_map<std::string, void *> gsvd_fcn;
48
49 static bool have_DGGSVD3 = false;
50 static bool gsvd_initialized = false;
51
52 /* Hack to stringize results of F77_FUNC macro. */
53 #define xSTRINGIZE(x) #x
54 #define STRINGIZE(x) xSTRINGIZE(x)
55
56 static void initialize_gsvd (void)
57 {
58 if (gsvd_initialized)
59 return;
60
61 dynamic_library libs ("");
62 if (! libs)
63 (*current_liboctave_error_handler)
64 ("gsvd: unable to query LAPACK library");
65
66 have_DGGSVD3 = (libs.search (STRINGIZE (F77_FUNC (dggsvd3, DGGSVD3)))
67 != nullptr);
68
69 if (have_DGGSVD3)
70 {
71 gsvd_fcn["dg"] = libs.search (STRINGIZE (F77_FUNC (dggsvd3, DGGSVD3)));
72 gsvd_fcn["sg"] = libs.search (STRINGIZE (F77_FUNC (sggsvd3, SGGSVD3)));
73 gsvd_fcn["zg"] = libs.search (STRINGIZE (F77_FUNC (zggsvd3, ZGGSVD3)));
74 gsvd_fcn["cg"] = libs.search (STRINGIZE (F77_FUNC (cggsvd3, CGGSVD3)));
75 }
76 else
77 {
78 gsvd_fcn["dg"] = libs.search (STRINGIZE (F77_FUNC (dggsvd, DGGSVD)));
79 gsvd_fcn["sg"] = libs.search (STRINGIZE (F77_FUNC (sggsvd, SGGSVD)));
80 gsvd_fcn["zg"] = libs.search (STRINGIZE (F77_FUNC (zggsvd, ZGGSVD)));
81 gsvd_fcn["cg"] = libs.search (STRINGIZE (F77_FUNC (cggsvd, CGGSVD)));
82 }
83
84 gsvd_initialized = true;
85 }
86
87 /* Clean up macro namespace as soon as we are done using them */
88 #undef xSTRINGIZE
89 #undef STRINGIZE
90
91 template<class T1>
92 struct real_ggsvd_ptr
93 {
94 typedef F77_RET_T (*type)
95 (F77_CONST_CHAR_ARG_DECL, // JOBU
96 F77_CONST_CHAR_ARG_DECL, // JOBV
97 F77_CONST_CHAR_ARG_DECL, // JOBQ
98 const F77_INT&, // M
99 const F77_INT&, // N
100 const F77_INT&, // P
101 F77_INT&, // K
102 F77_INT&, // L
103 T1 *, // A(LDA,N)
104 const F77_INT&, // LDA
105 T1 *, // B(LDB,N)
106 const F77_INT&, // LDB
107 T1 *, // ALPHA(N)
108 T1 *, // BETA(N)
109 T1 *, // U(LDU,M)
110 const F77_INT&, // LDU
111 T1 *, // V(LDV,P)
112 const F77_INT&, // LDV
113 T1 *, // Q(LDQ,N)
114 const F77_INT&, // LDQ
115 T1 *, // WORK
116 F77_INT *, // IWORK(N)
117 F77_INT& // INFO
118 F77_CHAR_ARG_LEN_DECL
119 F77_CHAR_ARG_LEN_DECL
120 F77_CHAR_ARG_LEN_DECL);
121 };
122
123 template<class T1>
124 struct real_ggsvd3_ptr
125 {
126 typedef F77_RET_T (*type)
127 (F77_CONST_CHAR_ARG_DECL, // JOBU
128 F77_CONST_CHAR_ARG_DECL, // JOBV
129 F77_CONST_CHAR_ARG_DECL, // JOBQ
130 const F77_INT&, // M
131 const F77_INT&, // N
132 const F77_INT&, // P
133 F77_INT&, // K
134 F77_INT&, // L
135 T1 *, // A(LDA,N)
136 const F77_INT&, // LDA
137 T1 *, // B(LDB,N)
138 const F77_INT&, // LDB
139 T1 *, // ALPHA(N)
140 T1 *, // BETA(N)
141 T1 *, // U(LDU,M)
142 const F77_INT&, // LDU
143 T1 *, // V(LDV,P)
144 const F77_INT&, // LDV
145 T1 *, // Q(LDQ,N)
146 const F77_INT&, // LDQ
147 T1 *, // WORK
148 const F77_INT&, // LWORK
149 F77_INT *, // IWORK(N)
150 F77_INT& // INFO
151 F77_CHAR_ARG_LEN_DECL
152 F77_CHAR_ARG_LEN_DECL
153 F77_CHAR_ARG_LEN_DECL);
154 };
155
156 template<class T1, class T2>
157 struct comp_ggsvd_ptr
158 {
159 typedef F77_RET_T (*type)
160 (F77_CONST_CHAR_ARG_DECL, // JOBU
161 F77_CONST_CHAR_ARG_DECL, // JOBV
162 F77_CONST_CHAR_ARG_DECL, // JOBQ
163 const F77_INT&, // M
164 const F77_INT&, // N
165 const F77_INT&, // P
166 F77_INT&, // K
167 F77_INT&, // L
168 T1 *, // A(LDA,N)
169 const F77_INT&, // LDA
170 T1 *, // B(LDB,N)
171 const F77_INT&, // LDB
172 T2 *, // ALPHA(N)
173 T2 *, // BETA(N)
174 T1 *, // U(LDU,M)
175 const F77_INT&, // LDU
176 T1 *, // V(LDV,P)
177 const F77_INT&, // LDV
178 T1 *, // Q(LDQ,N)
179 const F77_INT&, // LDQ
180 T1 *, // WORK
181 T2 *, // RWORK
182 F77_INT *, // IWORK(N)
183 F77_INT& // INFO
184 F77_CHAR_ARG_LEN_DECL
185 F77_CHAR_ARG_LEN_DECL
186 F77_CHAR_ARG_LEN_DECL);
187 };
188
189 template<class T1, class T2>
190 struct comp_ggsvd3_ptr
191 {
192 typedef F77_RET_T (*type)
193 (F77_CONST_CHAR_ARG_DECL, // JOBU
194 F77_CONST_CHAR_ARG_DECL, // JOBV
195 F77_CONST_CHAR_ARG_DECL, // JOBQ
196 const F77_INT&, // M
197 const F77_INT&, // N
198 const F77_INT&, // P
199 F77_INT&, // K
200 F77_INT&, // L
201 T1 *, // A(LDA,N)
202 const F77_INT&, // LDA
203 T1 *, // B(LDB,N)
204 const F77_INT&, // LDB
205 T2 *, // ALPHA(N)
206 T2 *, // BETA(N)
207 T1 *, // U(LDU,M)
208 const F77_INT&, // LDU
209 T1 *, // V(LDV,P)
210 const F77_INT&, // LDV
211 T1 *, // Q(LDQ,N)
212 const F77_INT&, // LDQ
213 T1 *, // WORK
214 const F77_INT&, // LWORK
215 T2 *, // RWORK
216 F77_INT *, // IWORK(N)
217 F77_INT& // INFO
218 F77_CHAR_ARG_LEN_DECL
219 F77_CHAR_ARG_LEN_DECL
220 F77_CHAR_ARG_LEN_DECL);
221 };
222
223 // template specializations
224 typedef real_ggsvd_ptr<F77_DBLE>::type dggsvd_type;
225 typedef real_ggsvd3_ptr<F77_DBLE>::type dggsvd3_type;
226 typedef real_ggsvd_ptr<F77_REAL>::type sggsvd_type;
227 typedef real_ggsvd3_ptr<F77_REAL>::type sggsvd3_type;
228 typedef comp_ggsvd_ptr<F77_DBLE_CMPLX, F77_DBLE>::type zggsvd_type;
229 typedef comp_ggsvd3_ptr<F77_DBLE_CMPLX, F77_DBLE>::type zggsvd3_type;
230 typedef comp_ggsvd_ptr<F77_CMPLX, F77_REAL>::type cggsvd_type;
231 typedef comp_ggsvd3_ptr<F77_CMPLX, F77_REAL>::type cggsvd3_type;
232
233 namespace math
234 {
235 template <>
236 void
237 gsvd<Matrix>::ggsvd (char& jobu, char& jobv, char& jobq, F77_INT m,
238 F77_INT n, F77_INT p, F77_INT& k, F77_INT& l,
239 double *tmp_dataA, F77_INT m1, double *tmp_dataB,
240 F77_INT p1, Matrix& alpha, Matrix& beta, double *u,
241 F77_INT nrow_u, double *v, F77_INT nrow_v, double *q,
242 F77_INT nrow_q, double *work, F77_INT lwork,
243 F77_INT *iwork, F77_INT& info)
244 {
245 if (! gsvd_initialized)
246 initialize_gsvd ();
247
248 if (have_DGGSVD3)
249 {
250 dggsvd3_type f_ptr = reinterpret_cast<dggsvd3_type> (gsvd_fcn["dg"]);
251 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
252 F77_CONST_CHAR_ARG2 (&jobv, 1),
253 F77_CONST_CHAR_ARG2 (&jobq, 1),
254 m, n, p, k, l, tmp_dataA, m1, tmp_dataB, p1,
255 alpha.fortran_vec (), beta.fortran_vec (),
256 u, nrow_u, v, nrow_v, q, nrow_q,
257 work, lwork, iwork, info
258 F77_CHAR_ARG_LEN (1)
259 F77_CHAR_ARG_LEN (1)
260 F77_CHAR_ARG_LEN (1));
261 }
262 else
263 {
264 dggsvd_type f_ptr = reinterpret_cast<dggsvd_type> (gsvd_fcn["dg"]);
265 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
266 F77_CONST_CHAR_ARG2 (&jobv, 1),
267 F77_CONST_CHAR_ARG2 (&jobq, 1),
268 m, n, p, k, l, tmp_dataA, m1, tmp_dataB, p1,
269 alpha.fortran_vec (), beta.fortran_vec (),
270 u, nrow_u, v, nrow_v, q, nrow_q,
271 work, iwork, info
272 F77_CHAR_ARG_LEN (1)
273 F77_CHAR_ARG_LEN (1)
274 F77_CHAR_ARG_LEN (1));
275 }
276 }
277
278 template <>
279 void
280 gsvd<FloatMatrix>::ggsvd (char& jobu, char& jobv, char& jobq, F77_INT m,
281 F77_INT n, F77_INT p, F77_INT& k, F77_INT& l,
282 float *tmp_dataA, F77_INT m1, float *tmp_dataB,
283 F77_INT p1, FloatMatrix& alpha,
284 FloatMatrix& beta, float *u, F77_INT nrow_u,
285 float *v, F77_INT nrow_v, float *q,
286 F77_INT nrow_q, float *work,
287 F77_INT lwork, F77_INT *iwork, F77_INT& info)
288 {
289 if (! gsvd_initialized)
290 initialize_gsvd ();
291
292 if (have_DGGSVD3)
293 {
294 sggsvd3_type f_ptr = reinterpret_cast<sggsvd3_type> (gsvd_fcn["sg"]);
295 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
296 F77_CONST_CHAR_ARG2 (&jobv, 1),
297 F77_CONST_CHAR_ARG2 (&jobq, 1),
298 m, n, p, k, l, tmp_dataA, m1, tmp_dataB, p1,
299 alpha.fortran_vec (), beta.fortran_vec (),
300 u, nrow_u, v, nrow_v, q, nrow_q,
301 work, lwork, iwork, info
302 F77_CHAR_ARG_LEN (1)
303 F77_CHAR_ARG_LEN (1)
304 F77_CHAR_ARG_LEN (1));
305 }
306 else
307 {
308 sggsvd_type f_ptr = reinterpret_cast<sggsvd_type> (gsvd_fcn["sg"]);
309 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
310 F77_CONST_CHAR_ARG2 (&jobv, 1),
311 F77_CONST_CHAR_ARG2 (&jobq, 1),
312 m, n, p, k, l, tmp_dataA, m1, tmp_dataB, p1,
313 alpha.fortran_vec (), beta.fortran_vec (),
314 u, nrow_u, v, nrow_v, q, nrow_q,
315 work, iwork, info
316 F77_CHAR_ARG_LEN (1)
317 F77_CHAR_ARG_LEN (1)
318 F77_CHAR_ARG_LEN (1));
319 }
320 }
321
322 template <>
323 void
324 gsvd<ComplexMatrix>::ggsvd (char& jobu, char& jobv, char& jobq,
325 F77_INT m, F77_INT n, F77_INT p, F77_INT& k,
326 F77_INT& l, Complex *tmp_dataA, F77_INT m1,
327 Complex *tmp_dataB, F77_INT p1, Matrix& alpha,
328 Matrix& beta, Complex *u, F77_INT nrow_u,
329 Complex *v, F77_INT nrow_v, Complex *q,
330 F77_INT nrow_q, Complex *work,
331 F77_INT lwork, F77_INT *iwork, F77_INT& info)
332 {
333 if (! gsvd_initialized)
334 initialize_gsvd ();
335
336 OCTAVE_LOCAL_BUFFER(double, rwork, 2*n);
337
338 if (have_DGGSVD3)
339 {
340 zggsvd3_type f_ptr = reinterpret_cast<zggsvd3_type> (gsvd_fcn["zg"]);
341 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
342 F77_CONST_CHAR_ARG2 (&jobv, 1),
343 F77_CONST_CHAR_ARG2 (&jobq, 1),
344 m, n, p, k, l,
345 F77_DBLE_CMPLX_ARG (tmp_dataA), m1,
346 F77_DBLE_CMPLX_ARG (tmp_dataB), p1,
347 alpha.fortran_vec (), beta.fortran_vec (),
348 F77_DBLE_CMPLX_ARG (u), nrow_u,
349 F77_DBLE_CMPLX_ARG (v), nrow_v,
350 F77_DBLE_CMPLX_ARG (q), nrow_q,
351 F77_DBLE_CMPLX_ARG (work),
352 lwork, rwork, iwork, info
353 F77_CHAR_ARG_LEN (1)
354 F77_CHAR_ARG_LEN (1)
355 F77_CHAR_ARG_LEN (1));
356 }
357 else
358 {
359 zggsvd_type f_ptr = reinterpret_cast<zggsvd_type> (gsvd_fcn["zg"]);
360 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
361 F77_CONST_CHAR_ARG2 (&jobv, 1),
362 F77_CONST_CHAR_ARG2 (&jobq, 1),
363 m, n, p, k, l,
364 F77_DBLE_CMPLX_ARG (tmp_dataA), m1,
365 F77_DBLE_CMPLX_ARG (tmp_dataB), p1,
366 alpha.fortran_vec (), beta.fortran_vec (),
367 F77_DBLE_CMPLX_ARG (u), nrow_u,
368 F77_DBLE_CMPLX_ARG (v), nrow_v,
369 F77_DBLE_CMPLX_ARG (q), nrow_q,
370 F77_DBLE_CMPLX_ARG (work),
371 rwork, iwork, info
372 F77_CHAR_ARG_LEN (1)
373 F77_CHAR_ARG_LEN (1)
374 F77_CHAR_ARG_LEN (1));
375 }
376 }
377
378 template <>
379 void
380 gsvd<FloatComplexMatrix>::ggsvd (char& jobu, char& jobv, char& jobq,
381 F77_INT m, F77_INT n, F77_INT p,
382 F77_INT& k, F77_INT& l,
383 FloatComplex *tmp_dataA, F77_INT m1,
384 FloatComplex *tmp_dataB, F77_INT p1,
385 FloatMatrix& alpha, FloatMatrix& beta,
386 FloatComplex *u, F77_INT nrow_u,
387 FloatComplex *v, F77_INT nrow_v,
388 FloatComplex *q, F77_INT nrow_q,
389 FloatComplex *work, F77_INT lwork,
390 F77_INT *iwork, F77_INT& info)
391 {
392 if (! gsvd_initialized)
393 initialize_gsvd ();
394
395 OCTAVE_LOCAL_BUFFER(float, rwork, 2*n);
396
397 if (have_DGGSVD3)
398 {
399 cggsvd3_type f_ptr = reinterpret_cast<cggsvd3_type> (gsvd_fcn["cg"]);
400 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
401 F77_CONST_CHAR_ARG2 (&jobv, 1),
402 F77_CONST_CHAR_ARG2 (&jobq, 1),
403 m, n, p, k, l,
404 F77_CMPLX_ARG (tmp_dataA), m1,
405 F77_CMPLX_ARG (tmp_dataB), p1,
406 alpha.fortran_vec (), beta.fortran_vec (),
407 F77_CMPLX_ARG (u), nrow_u,
408 F77_CMPLX_ARG (v), nrow_v,
409 F77_CMPLX_ARG (q), nrow_q,
410 F77_CMPLX_ARG (work), lwork,
411 rwork, iwork, info
412 F77_CHAR_ARG_LEN (1)
413 F77_CHAR_ARG_LEN (1)
414 F77_CHAR_ARG_LEN (1));
415 }
416 else
417 {
418 cggsvd_type f_ptr = reinterpret_cast<cggsvd_type> (gsvd_fcn["cg"]);
419 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
420 F77_CONST_CHAR_ARG2 (&jobv, 1),
421 F77_CONST_CHAR_ARG2 (&jobq, 1),
422 m, n, p, k, l,
423 F77_CMPLX_ARG (tmp_dataA), m1,
424 F77_CMPLX_ARG (tmp_dataB), p1,
425 alpha.fortran_vec (), beta.fortran_vec (),
426 F77_CMPLX_ARG (u), nrow_u,
427 F77_CMPLX_ARG (v), nrow_v,
428 F77_CMPLX_ARG (q), nrow_q,
429 F77_CMPLX_ARG (work),
430 rwork, iwork, info
431 F77_CHAR_ARG_LEN (1)
432 F77_CHAR_ARG_LEN (1)
433 F77_CHAR_ARG_LEN (1));
434 }
435 }
436
437 template <typename T>
438 T
439 gsvd<T>::left_singular_matrix_A (void) const
440 {
441 if (m_type == gsvd::Type::sigma_only)
442 (*current_liboctave_error_handler)
443 ("gsvd: U not computed because type == gsvd::sigma_only");
444
445 return m_left_smA;
446 }
447
448 template <typename T>
449 T
450 gsvd<T>::left_singular_matrix_B (void) const
451 {
452 if (m_type == gsvd::Type::sigma_only)
453 (*current_liboctave_error_handler)
454 ("gsvd: V not computed because type == gsvd::sigma_only");
455
456 return m_left_smB;
457 }
458
459 template <typename T>
460 T
461 gsvd<T>::right_singular_matrix (void) const
462 {
463 if (m_type == gsvd::Type::sigma_only)
464 (*current_liboctave_error_handler)
465 ("gsvd: X not computed because type == gsvd::sigma_only");
466
467 return m_right_sm;
468 }
469
470 template <typename T>
471 gsvd<T>::gsvd (const T& a, const T& b, gsvd::Type gsvd_type)
472 {
473 if (a.isempty () || b.isempty ())
474 (*current_liboctave_error_handler)
475 ("gsvd: A and B cannot be empty matrices");
476
477 F77_INT info;
478
479 F77_INT m = to_f77_int (a.rows ());
480 F77_INT n = to_f77_int (a.cols ());
481 F77_INT p = to_f77_int (b.rows ());
482
483 T atmp = a;
484 P *tmp_dataA = atmp.fortran_vec ();
485
486 T btmp = b;
487 P *tmp_dataB = btmp.fortran_vec ();
488
489 char jobu = 'U';
490 char jobv = 'V';
491 char jobq = 'Q';
492
493 F77_INT nrow_u = m;
494 F77_INT nrow_v = p;
495 F77_INT nrow_q = n;
496
497 F77_INT k, l;
498
499 switch (gsvd_type)
500 {
501 case gsvd<T>::Type::sigma_only:
502
503 // FIXME: In LAPACK 3.0, problem below seems to be fixed,
504 // so we now set jobX = 'N'.
505 //
506 // For calculating sigma_only, both jobu and jobv should be 'N', but
507 // there seems to be a bug in dgesvd from LAPACK V2.0. To
508 // demonstrate the bug, set both jobu and jobv to 'N' and find the
509 // singular values of [eye(3), eye(3)]. The result is
510 // [-sqrt(2), -sqrt(2), -sqrt(2)].
511
512 jobu = jobv = jobq = 'N';
513 nrow_u = nrow_v = nrow_q = 1;
514 break;
515
516 default:
517 break;
518 }
519
520 m_type = gsvd_type;
521
522 if (jobu != 'N')
523 m_left_smA.resize (nrow_u, m);
524
525 P *u = m_left_smA.fortran_vec ();
526
527 if (jobv != 'N')
528 m_left_smB.resize (nrow_v, p);
529
530 P *v = m_left_smB.fortran_vec ();
531
532 if (jobq != 'N')
533 m_right_sm.resize (nrow_q, n);
534
535 P *q = m_right_sm.fortran_vec ();
536
537 real_matrix alpha (n, 1);
538 real_matrix beta (n, 1);
539
540 OCTAVE_LOCAL_BUFFER(F77_INT, iwork, n);
541
542 if (! gsvd_initialized)
543 initialize_gsvd ();
544
545 F77_INT lwork;
546 if (have_DGGSVD3)
547 {
548 lwork = -1;
549 P work_tmp;
550
551 gsvd<T>::ggsvd (jobu, jobv, jobq, m, n, p, k, l,
552 tmp_dataA, m, tmp_dataB, p,
553 alpha, beta, u, nrow_u, v, nrow_v, q, nrow_q,
554 &work_tmp, lwork, iwork, info);
555
556 lwork = static_cast<F77_INT> (std::abs (work_tmp));
557 }
558 else
559 {
560 lwork = std::max ({3*n, m, p}) + n;
561 }
562 info = 0;
563
564 OCTAVE_LOCAL_BUFFER(P, work, lwork);
565
566 gsvd<T>::ggsvd (jobu, jobv, jobq, m, n, p, k, l,
567 tmp_dataA, m, tmp_dataB, p,
568 alpha, beta, u, nrow_u, v, nrow_v, q, nrow_q,
569 work, lwork, iwork, info);
570
571 if (info < 0)
572 (*current_liboctave_error_handler)
573 ("*ggsvd.f: argument %d illegal", -info);
574
575 if (info > 0)
576 (*current_liboctave_error_handler)
577 ("*ggsvd.f: Jacobi-type procedure failed to converge");
578
579 F77_INT i, j;
580
581 if (gsvd_type != gsvd<T>::Type::sigma_only)
582 {
583 // Size according to LAPACK is k+l,k+l, but this needs
584 // to be nxn for Matlab compatibility.
585 T R (n, n, 0.0);
586 int astart = n-k-l;
587 if (m - k - l >= 0)
588 {
589 // R is stored in A(1:K+L,N-K-L+1:N)
590 for (i = 0; i < k+l; i++)
591 for (j = 0; j < k+l; j++)
592 R.xelem (i, j) = atmp.xelem (i, astart + j);
593 }
594 else
595 {
596 // (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N)
597 // ( 0 R22 R23 )
598 for (i = 0; i < m; i++)
599 for (j = 0; j < k+l; j++)
600 R.xelem (i, j) = atmp.xelem (i, astart + j);
601 // and R33 is stored in B(M-K+1:L,N+M-K-L+1:N)
602 for (i = m; i < k + l; i++)
603 for (j = n - l - k + m; j < n; j++)
604 R.xelem (i, j) = btmp.xelem (i - k, j);
605 }
606
607 // Output X = Q*R'
608 // FIXME: Is there a way to call BLAS multiply directly
609 // with flags so that R is transposed?
610 m_right_sm = m_right_sm * R.hermitian ();
611 }
612
613 // Fill in C and S
614 F77_INT rank;
615 bool fill_ptn;
616 if (m-k-l >= 0)
617 {
618 rank = l;
619 fill_ptn = true;
620 }
621 else
622 {
623 rank = m-k;
624 fill_ptn = false;
625 }
626
627 if (gsvd_type == gsvd<T>::Type::sigma_only)
628 {
629 // Return column vector with results
630 m_sigmaA.resize (k+l, 1);
631 m_sigmaB.resize (k+l, 1);
632
633 if (fill_ptn)
634 {
635 for (i = 0; i < k; i++)
636 {
637 m_sigmaA.xelem (i) = 1.0;
638 m_sigmaB.xelem (i) = 0.0;
639 }
640 for (i = k, j = k+l-1; i < k+l; i++, j--)
641 {
642 m_sigmaA.xelem (i) = alpha.xelem (i);
643 m_sigmaB.xelem (i) = beta.xelem (i);
644 }
645 }
646 else
647 {
648 for (i = 0; i < k; i++)
649 {
650 m_sigmaA.xelem (i) = 1.0;
651 m_sigmaB.xelem (i) = 0.0;
652 }
653 for (i = k; i < m; i++)
654 {
655 m_sigmaA.xelem (i) = alpha.xelem (i);
656 m_sigmaB.xelem (i) = beta.xelem (i);
657 }
658 for (i = m; i < k+l; i++)
659 {
660 m_sigmaA.xelem (i) = 0.0;
661 m_sigmaB.xelem (i) = 1.0;
662 }
663 }
664 }
665 else // returning all matrices
666 {
667 // Number of columns according to LAPACK is k+l, but this needs
668 // to be n for Matlab compatibility.
669 m_sigmaA.resize (m, n);
670 m_sigmaB.resize (p, n);
671
672 for (i = 0; i < k; i++)
673 m_sigmaA.xelem (i, i) = 1.0;
674
675 for (i = 0; i < rank; i++)
676 {
677 m_sigmaA.xelem (k+i, k+i) = alpha.xelem (k+i);
678 m_sigmaB.xelem (i, k+i) = beta.xelem (k+i);
679 }
680
681 if (! fill_ptn)
682 {
683 for (i = m; i < n; i++)
684 m_sigmaB.xelem (i-k, i) = 1.0;
685 }
686
687 }
688 }
689
690 // Instantiations needed in octave::math namespace.
691 template class gsvd<Matrix>;
692 template class gsvd<FloatMatrix>;
693 template class gsvd<ComplexMatrix>;
694 template class gsvd<FloatComplexMatrix>;
695 }
696}
max
charNDArray max(char d, const charNDArray &m)
Definition: chNDArray.cc:230
Array::fortran_vec
OCTARRAY_API T * fortran_vec(void)
Size of the specified dimension.
Definition: Array.cc:1744
Matrix
Definition: dMatrix.h:42
octave::dynamic_library
Definition: oct-shlib.h:43
octave::dynamic_library::search
void * search(const std::string &nm, const name_mangler &mangler=name_mangler()) const
Definition: oct-shlib.h:177
octave::math::gsvd::real_matrix
T::real_matrix_type real_matrix
Definition: gsvd.h:91
octave::math::gsvd::right_singular_matrix
T right_singular_matrix(void) const
Definition: gsvd.cc:461
octave::math::gsvd::P
T::value_type P
Definition: gsvd.h:90
octave::math::gsvd::left_singular_matrix_A
T left_singular_matrix_A(void) const
Definition: gsvd.cc:439
octave::math::gsvd::left_singular_matrix_B
T left_singular_matrix_B(void) const
Definition: gsvd.cc:450
octave::math::gsvd::ggsvd
void ggsvd(char &jobu, char &jobv, char &jobq, octave_f77_int_type m, octave_f77_int_type n, octave_f77_int_type p, octave_f77_int_type &k, octave_f77_int_type &l, P *tmp_dataA, octave_f77_int_type m1, P *tmp_dataB, octave_f77_int_type p1, real_matrix &alpha, real_matrix &beta, P *u, octave_f77_int_type nrow_u, P *v, octave_f77_int_type nrow_v, P *q, octave_f77_int_type nrow_q, P *work, octave_f77_int_type lwork, octave_f77_int_type *iwork, octave_f77_int_type &info)
octave::math::gsvd::gsvd
gsvd(void)
Definition: gsvd.h:49
F77_DBLE_CMPLX_ARG
#define F77_DBLE_CMPLX_ARG(x)
Definition: f77-fcn.h:316
F77_CMPLX_ARG
#define F77_CMPLX_ARG(x)
Definition: f77-fcn.h:310
F77_INT
octave_f77_int_type F77_INT
Definition: f77-fcn.h:306
gsvd_type
static OCTAVE_NAMESPACE_BEGIN math::gsvd< T >::Type gsvd_type(int nargout, int nargin)
Definition: gsvd.cc:47
STRINGIZE
#define STRINGIZE(x)
Definition: gsvd.cc:54
F77_CONST_CHAR_ARG_DECL
F77_RET_T F77_CONST_CHAR_ARG_DECL
Definition: lo-arpack-proto.h:42
F77_CHAR_ARG_LEN_DECL
F77_RET_T const F77_INT F77_INT const F77_DBLE F77_DBLE const F77_INT F77_DBLE const F77_INT F77_INT F77_INT F77_DBLE F77_DBLE const F77_INT F77_INT &F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL
Definition: lo-arpack-proto.h:45
F77_RET_T
F77_RET_T(F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, const F77_INT &, const F77_INT &, const F77_INT &, F77_INT &, F77_INT &, F77_DBLE *, const F77_INT &, F77_DBLE *, const F77_INT &, F77_DBLE *, F77_DBLE *, F77_DBLE *, const F77_INT &, F77_DBLE *, const F77_INT &, F77_DBLE *, const F77_INT &, F77_DBLE *, F77_INT *, F77_INT &F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL)
octave::initialize_gsvd
static void initialize_gsvd(void)
Definition: gsvd.cc:56
octave::gsvd_initialized
static bool gsvd_initialized
Definition: gsvd.cc:50
octave::cggsvd3_type
comp_ggsvd3_ptr< F77_CMPLX, F77_REAL >::type cggsvd3_type
Definition: gsvd.cc:231
octave::zggsvd_type
comp_ggsvd_ptr< F77_DBLE_CMPLX, F77_DBLE >::type zggsvd_type
Definition: gsvd.cc:228
octave::F77_FUNC
F77_RET_T F77_FUNC(dconv2o, DCONV2O)(const F77_INT &
octave::have_DGGSVD3
static bool have_DGGSVD3
Definition: gsvd.cc:49
octave::zggsvd3_type
comp_ggsvd3_ptr< F77_DBLE_CMPLX, F77_DBLE >::type zggsvd3_type
Definition: gsvd.cc:229
octave::gsvd_fcn
static std::unordered_map< std::string, void * > gsvd_fcn
Definition: gsvd.cc:47
octave::sggsvd_type
real_ggsvd_ptr< F77_REAL >::type sggsvd_type
Definition: gsvd.cc:226
octave::dggsvd3_type
real_ggsvd3_ptr< F77_DBLE >::type dggsvd3_type
Definition: gsvd.cc:225
octave::dggsvd_type
real_ggsvd_ptr< F77_DBLE >::type dggsvd_type
Definition: gsvd.cc:224
octave::cggsvd_type
comp_ggsvd_ptr< F77_CMPLX, F77_REAL >::type cggsvd_type
Definition: gsvd.cc:230
octave::sggsvd3_type
real_ggsvd3_ptr< F77_REAL >::type sggsvd3_type
Definition: gsvd.cc:227
Complex
std::complex< double > Complex
Definition: oct-cmplx.h:33
FloatComplex
std::complex< float > FloatComplex
Definition: oct-cmplx.h:34
OCTAVE_LOCAL_BUFFER
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:44
abs
static T abs(T x)
Definition: pr-output.cc:1678
octave::comp_ggsvd3_ptr
Definition: gsvd.cc:191
octave::comp_ggsvd3_ptr::type
F77_RET_T(* type)(F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, const F77_INT &, const F77_INT &, const F77_INT &, F77_INT &, F77_INT &, T1 *, const F77_INT &, T1 *, const F77_INT &, T2 *, T2 *, T1 *, const F77_INT &, T1 *, const F77_INT &, T1 *, const F77_INT &, T1 *, const F77_INT &, T2 *, F77_INT *, F77_INT &F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL)
Definition: gsvd.cc:193
octave::comp_ggsvd_ptr
Definition: gsvd.cc:158
octave::comp_ggsvd_ptr::type
F77_RET_T(* type)(F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, const F77_INT &, const F77_INT &, const F77_INT &, F77_INT &, F77_INT &, T1 *, const F77_INT &, T1 *, const F77_INT &, T2 *, T2 *, T1 *, const F77_INT &, T1 *, const F77_INT &, T1 *, const F77_INT &, T1 *, T2 *, F77_INT *, F77_INT &F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL)
Definition: gsvd.cc:160
octave::real_ggsvd3_ptr
Definition: gsvd.cc:125
octave::real_ggsvd3_ptr::type
F77_RET_T(* type)(F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, const F77_INT &, const F77_INT &, const F77_INT &, F77_INT &, F77_INT &, T1 *, const F77_INT &, T1 *, const F77_INT &, T1 *, T1 *, T1 *, const F77_INT &, T1 *, const F77_INT &, T1 *, const F77_INT &, T1 *, const F77_INT &, F77_INT *, F77_INT &F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL)
Definition: gsvd.cc:127
octave::real_ggsvd_ptr
Definition: gsvd.cc:93
octave::real_ggsvd_ptr::type
F77_RET_T(* type)(F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, const F77_INT &, const F77_INT &, const F77_INT &, F77_INT &, F77_INT &, T1 *, const F77_INT &, T1 *, const F77_INT &, T1 *, T1 *, T1 *, const F77_INT &, T1 *, const F77_INT &, T1 *, const F77_INT &, T1 *, F77_INT *, F77_INT &F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL)
Definition: gsvd.cc:95