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