GNU Octave  8.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-2023 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 
45 OCTAVE_BEGIN_NAMESPACE(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 OCTAVE_BEGIN_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 %" OCTAVE_F77_INT_TYPE_FORMAT " illegal",
574  -info);
575 
576  if (info > 0)
577  (*current_liboctave_error_handler)
578  ("*ggsvd.f: Jacobi-type procedure failed to converge");
579 
580  F77_INT i, j;
581 
582  if (gsvd_type != gsvd<T>::Type::sigma_only)
583  {
584  // Size according to LAPACK is k+l,k+l, but this needs
585  // to be nxn for Matlab compatibility.
586  T R (n, n, 0.0);
587  int astart = n-k-l;
588  if (m - k - l >= 0)
589  {
590  // R is stored in A(1:K+L,N-K-L+1:N)
591  for (i = 0; i < k+l; i++)
592  for (j = 0; j < k+l; j++)
593  R.xelem (i, j) = atmp.xelem (i, astart + j);
594  }
595  else
596  {
597  // (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N)
598  // ( 0 R22 R23 )
599  for (i = 0; i < m; i++)
600  for (j = 0; j < k+l; j++)
601  R.xelem (i, j) = atmp.xelem (i, astart + j);
602  // and R33 is stored in B(M-K+1:L,N+M-K-L+1:N)
603  for (i = m; i < k + l; i++)
604  for (j = n - l - k + m; j < n; j++)
605  R.xelem (i, j) = btmp.xelem (i - k, j);
606  }
607 
608  // Output X = Q*R'
609  // FIXME: Is there a way to call BLAS multiply directly
610  // with flags so that R is transposed?
611  m_right_sm = m_right_sm * R.hermitian ();
612  }
613 
614  // Fill in C and S
615  F77_INT rank;
616  bool fill_ptn;
617  if (m-k-l >= 0)
618  {
619  rank = l;
620  fill_ptn = true;
621  }
622  else
623  {
624  rank = m-k;
625  fill_ptn = false;
626  }
627 
628  if (gsvd_type == gsvd<T>::Type::sigma_only)
629  {
630  // Return column vector with results
631  m_sigmaA.resize (k+l, 1);
632  m_sigmaB.resize (k+l, 1);
633 
634  if (fill_ptn)
635  {
636  for (i = 0; i < k; i++)
637  {
638  m_sigmaA.xelem (i) = 1.0;
639  m_sigmaB.xelem (i) = 0.0;
640  }
641  for (i = k, j = k+l-1; i < k+l; i++, j--)
642  {
643  m_sigmaA.xelem (i) = alpha.xelem (i);
644  m_sigmaB.xelem (i) = beta.xelem (i);
645  }
646  }
647  else
648  {
649  for (i = 0; i < k; i++)
650  {
651  m_sigmaA.xelem (i) = 1.0;
652  m_sigmaB.xelem (i) = 0.0;
653  }
654  for (i = k; i < m; i++)
655  {
656  m_sigmaA.xelem (i) = alpha.xelem (i);
657  m_sigmaB.xelem (i) = beta.xelem (i);
658  }
659  for (i = m; i < k+l; i++)
660  {
661  m_sigmaA.xelem (i) = 0.0;
662  m_sigmaB.xelem (i) = 1.0;
663  }
664  }
665  }
666  else // returning all matrices
667  {
668  // Number of columns according to LAPACK is k+l, but this needs
669  // to be n for Matlab compatibility.
670  m_sigmaA.resize (m, n);
671  m_sigmaB.resize (p, n);
672 
673  for (i = 0; i < k; i++)
674  m_sigmaA.xelem (i, i) = 1.0;
675 
676  for (i = 0; i < rank; i++)
677  {
678  m_sigmaA.xelem (k+i, k+i) = alpha.xelem (k+i);
679  m_sigmaB.xelem (i, k+i) = beta.xelem (k+i);
680  }
681 
682  if (! fill_ptn)
683  {
684  for (i = m; i < n; i++)
685  m_sigmaB.xelem (i-k, i) = 1.0;
686  }
687 
688  }
689 }
690 
691 // Instantiations needed in octave::math namespace.
692 template class gsvd<Matrix>;
693 template class gsvd<FloatMatrix>;
694 template class gsvd<ComplexMatrix>;
695 template class gsvd<FloatComplexMatrix>;
696 
697 OCTAVE_END_NAMESPACE(math)
698 OCTAVE_END_NAMESPACE(octave)
OCTAVE_END_NAMESPACE
OCTAVE_END_NAMESPACE(octave)
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-base.cc:1766
Matrix
Definition: dMatrix.h:42
dynamic_library
Definition: oct-shlib.h:43
dynamic_library::search
void * search(const std::string &nm, const name_mangler &mangler=name_mangler()) const
Definition: oct-shlib.h:177
gsvd
Definition: gsvd.h:39
gsvd::P
T::value_type P
Definition: gsvd.h:90
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)
gsvd::right_singular_matrix
T right_singular_matrix(void) const
Definition: gsvd.cc:461
gsvd::real_matrix
T::real_matrix_type real_matrix
Definition: gsvd.h:91
gsvd::left_singular_matrix_A
T left_singular_matrix_A(void) const
Definition: gsvd.cc:439
gsvd::left_singular_matrix_B
T left_singular_matrix_B(void) const
Definition: gsvd.cc:450
gsvd::gsvd
gsvd(void)
Definition: gsvd.h:49
gsvd::Type
Type
Definition: gsvd.h:43
OCTAVE_BEGIN_NAMESPACE
OCTAVE_BEGIN_NAMESPACE(octave) static octave_value daspk_fcn
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 math::gsvd< T >::Type gsvd_type(int nargout, int nargin)
Definition: gsvd.cc:47
gsvd_fcn
static std::unordered_map< std::string, void * > gsvd_fcn
Definition: gsvd.cc:47
zggsvd3_type
comp_ggsvd3_ptr< F77_DBLE_CMPLX, F77_DBLE >::type zggsvd3_type
Definition: gsvd.cc:229
cggsvd_type
comp_ggsvd_ptr< F77_CMPLX, F77_REAL >::type cggsvd_type
Definition: gsvd.cc:230
dggsvd_type
real_ggsvd_ptr< F77_DBLE >::type dggsvd_type
Definition: gsvd.cc:224
cggsvd3_type
comp_ggsvd3_ptr< F77_CMPLX, F77_REAL >::type cggsvd3_type
Definition: gsvd.cc:231
STRINGIZE
#define STRINGIZE(x)
Definition: gsvd.cc:54
sggsvd3_type
real_ggsvd3_ptr< F77_REAL >::type sggsvd3_type
Definition: gsvd.cc:227
have_DGGSVD3
static bool have_DGGSVD3
Definition: gsvd.cc:49
initialize_gsvd
static void initialize_gsvd(void)
Definition: gsvd.cc:56
dggsvd3_type
real_ggsvd3_ptr< F77_DBLE >::type dggsvd3_type
Definition: gsvd.cc:225
zggsvd_type
comp_ggsvd_ptr< F77_DBLE_CMPLX, F77_DBLE >::type zggsvd_type
Definition: gsvd.cc:228
sggsvd_type
real_ggsvd_ptr< F77_REAL >::type sggsvd_type
Definition: gsvd.cc:226
gsvd_initialized
static bool gsvd_initialized
Definition: gsvd.cc:50
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)
m
T octave_idx_type m
Definition: mx-inlines.cc:773
n
octave_idx_type n
Definition: mx-inlines.cc:753
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
comp_ggsvd3_ptr
Definition: gsvd.cc:191
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
comp_ggsvd_ptr
Definition: gsvd.cc:158
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
real_ggsvd3_ptr
Definition: gsvd.cc:125
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
real_ggsvd_ptr
Definition: gsvd.cc:93
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
F77_FUNC
F77_RET_T F77_FUNC(xerbla, XERBLA)(F77_CONST_CHAR_ARG_DEF(s_arg