GNU Octave  6.2.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-2021 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 <vector>
31 
32 #include "CMatrix.h"
33 #include "dDiagMatrix.h"
34 #include "dMatrix.h"
35 #include "fCMatrix.h"
36 #include "fDiagMatrix.h"
37 #include "fMatrix.h"
38 #include "gsvd.h"
39 #include "lo-error.h"
40 #include "lo-lapack-proto.h"
41 #include "oct-shlib.h"
42 
43 namespace octave
44 {
45  static std::map<std::string, void *> gsvd_fcn;
46 
47  static bool have_DGGSVD3 = false;
48  static bool gsvd_initialized = false;
49 
50  /* Hack to stringize macro results. */
51 #define xSTRINGIZE(x) #x
52 #define STRINGIZE(x) xSTRINGIZE(x)
53 
54  static void initialize_gsvd (void)
55  {
56  if (gsvd_initialized)
57  return;
58 
59  dynamic_library libs ("");
60  if (! libs)
61  {
62  // FIXME: Should we throw an error if we cannot check the libraries?
63  have_DGGSVD3 = false;
64  return;
65  }
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  gsvd_initialized = true;
85  }
86 
87  template<class T1>
88  struct real_ggsvd_ptr
89  {
90  typedef F77_RET_T (*type)
91  (F77_CONST_CHAR_ARG_DECL, // JOBU
92  F77_CONST_CHAR_ARG_DECL, // JOBV
93  F77_CONST_CHAR_ARG_DECL, // JOBQ
94  const F77_INT&, // M
95  const F77_INT&, // N
96  const F77_INT&, // P
97  F77_INT &, // K
98  F77_INT &, // L
99  T1*, // A(LDA,N)
100  const F77_INT&, // LDA
101  T1*, // B(LDB,N)
102  const F77_INT&, // LDB
103  T1*, // ALPHA(N)
104  T1*, // BETA(N)
105  T1*, // U(LDU,M)
106  const F77_INT&, // LDU
107  T1*, // V(LDV,P)
108  const F77_INT&, // LDV
109  T1*, // Q(LDQ,N)
110  const F77_INT&, // LDQ
111  T1*, // WORK
112  F77_INT*, // IWORK(N)
113  F77_INT& // INFO
114  F77_CHAR_ARG_LEN_DECL
115  F77_CHAR_ARG_LEN_DECL
116  F77_CHAR_ARG_LEN_DECL);
117  };
118 
119  template<class T1>
120  struct real_ggsvd3_ptr
121  {
122  typedef F77_RET_T (*type)
123  (F77_CONST_CHAR_ARG_DECL, // JOBU
124  F77_CONST_CHAR_ARG_DECL, // JOBV
125  F77_CONST_CHAR_ARG_DECL, // JOBQ
126  const F77_INT&, // M
127  const F77_INT&, // N
128  const F77_INT&, // P
129  F77_INT &, // K
130  F77_INT &, // L
131  T1*, // A(LDA,N)
132  const F77_INT&, // LDA
133  T1*, // B(LDB,N)
134  const F77_INT&, // LDB
135  T1*, // ALPHA(N)
136  T1*, // BETA(N)
137  T1*, // U(LDU,M)
138  const F77_INT&, // LDU
139  T1*, // V(LDV,P)
140  const F77_INT&, // LDV
141  T1*, // Q(LDQ,N)
142  const F77_INT&, // LDQ
143  T1*, // WORK
144  const F77_INT&, // LWORK
145  F77_INT*, // IWORK(N)
146  F77_INT& // INFO
147  F77_CHAR_ARG_LEN_DECL
148  F77_CHAR_ARG_LEN_DECL
149  F77_CHAR_ARG_LEN_DECL);
150  };
151 
152  template<class T1, class T2>
153  struct comp_ggsvd_ptr
154  {
155  typedef F77_RET_T (*type)
156  (F77_CONST_CHAR_ARG_DECL, // JOBU
157  F77_CONST_CHAR_ARG_DECL, // JOBV
158  F77_CONST_CHAR_ARG_DECL, // JOBQ
159  const F77_INT&, // M
160  const F77_INT&, // N
161  const F77_INT&, // P
162  F77_INT &, // K
163  F77_INT &, // L
164  T1*, // A(LDA,N)
165  const F77_INT&, // LDA
166  T1*, // B(LDB,N)
167  const F77_INT&, // LDB
168  T2*, // ALPHA(N)
169  T2*, // BETA(N)
170  T1*, // U(LDU,M)
171  const F77_INT&, // LDU
172  T1*, // V(LDV,P)
173  const F77_INT&, // LDV
174  T1*, // Q(LDQ,N)
175  const F77_INT&, // LDQ
176  T1*, // WORK
177  T2*, // RWORK
178  F77_INT*, // IWORK(N)
179  F77_INT& // INFO
180  F77_CHAR_ARG_LEN_DECL
181  F77_CHAR_ARG_LEN_DECL
182  F77_CHAR_ARG_LEN_DECL);
183  };
184 
185  template<class T1, class T2>
186  struct comp_ggsvd3_ptr
187  {
188  typedef F77_RET_T (*type)
189  (F77_CONST_CHAR_ARG_DECL, // JOBU
190  F77_CONST_CHAR_ARG_DECL, // JOBV
191  F77_CONST_CHAR_ARG_DECL, // JOBQ
192  const F77_INT&, // M
193  const F77_INT&, // N
194  const F77_INT&, // P
195  F77_INT &, // K
196  F77_INT &, // L
197  T1*, // A(LDA,N)
198  const F77_INT&, // LDA
199  T1*, // B(LDB,N)
200  const F77_INT&, // LDB
201  T2*, // ALPHA(N)
202  T2*, // BETA(N)
203  T1*, // U(LDU,M)
204  const F77_INT&, // LDU
205  T1*, // V(LDV,P)
206  const F77_INT&, // LDV
207  T1*, // Q(LDQ,N)
208  const F77_INT&, // LDQ
209  T1*, // WORK
210  const F77_INT&, // LWORK
211  T2*, // RWORK
212  F77_INT*, // IWORK(N)
213  F77_INT& // INFO
214  F77_CHAR_ARG_LEN_DECL
215  F77_CHAR_ARG_LEN_DECL
216  F77_CHAR_ARG_LEN_DECL);
217  };
218 
219  // template specializations
220  typedef real_ggsvd3_ptr<F77_DBLE>::type dggsvd3_type;
221  typedef real_ggsvd_ptr<F77_DBLE>::type dggsvd_type;
222  typedef real_ggsvd3_ptr<F77_REAL>::type sggsvd3_type;
223  typedef real_ggsvd_ptr<F77_REAL>::type sggsvd_type;
224  typedef comp_ggsvd3_ptr<F77_DBLE_CMPLX, F77_DBLE>::type zggsvd3_type;
225  typedef comp_ggsvd_ptr<F77_DBLE_CMPLX, F77_DBLE>::type zggsvd_type;
226  typedef comp_ggsvd3_ptr<F77_CMPLX, F77_REAL>::type cggsvd3_type;
227  typedef comp_ggsvd_ptr<F77_CMPLX, F77_REAL>::type cggsvd_type;
228 
229  namespace math
230  {
231  template <>
232  void
233  gsvd<Matrix>::ggsvd (char& jobu, char& jobv, char& jobq, F77_INT m,
234  F77_INT n, F77_INT p, F77_INT& k, F77_INT& l,
235  double *tmp_dataA, F77_INT m1, double *tmp_dataB,
236  F77_INT p1, Matrix& alpha, Matrix& beta, double *u,
237  F77_INT nrow_u, double *v, F77_INT nrow_v, double *q,
238  F77_INT nrow_q, Matrix& work, F77_INT lwork,
239  F77_INT *iwork, F77_INT& info)
240  {
241  if (! gsvd_initialized)
242  initialize_gsvd ();
243 
244  if (have_DGGSVD3)
245  {
246  dggsvd3_type f_ptr = reinterpret_cast<dggsvd3_type> (gsvd_fcn["dg"]);
247  f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
248  F77_CONST_CHAR_ARG2 (&jobv, 1),
249  F77_CONST_CHAR_ARG2 (&jobq, 1),
250  m, n, p, k, l, tmp_dataA, m1, tmp_dataB, p1,
251  alpha.fortran_vec (), beta.fortran_vec (),
252  u, nrow_u, v, nrow_v, q, nrow_q,
253  work.fortran_vec (), lwork, iwork, info
254  F77_CHAR_ARG_LEN (1)
255  F77_CHAR_ARG_LEN (1)
256  F77_CHAR_ARG_LEN (1));
257  }
258  else
259  {
260  dggsvd_type f_ptr = reinterpret_cast<dggsvd_type> (gsvd_fcn["dg"]);
261  f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
262  F77_CONST_CHAR_ARG2 (&jobv, 1),
263  F77_CONST_CHAR_ARG2 (&jobq, 1),
264  m, n, p, k, l, tmp_dataA, m1, tmp_dataB, p1,
265  alpha.fortran_vec (), beta.fortran_vec (),
266  u, nrow_u, v, nrow_v, q, nrow_q,
267  work.fortran_vec (), iwork, info
268  F77_CHAR_ARG_LEN (1)
269  F77_CHAR_ARG_LEN (1)
270  F77_CHAR_ARG_LEN (1));
271  }
272  }
273 
274  template <>
275  void
276  gsvd<FloatMatrix>::ggsvd (char& jobu, char& jobv, char& jobq, F77_INT m,
277  F77_INT n, F77_INT p, F77_INT& k, F77_INT& l,
278  float *tmp_dataA, F77_INT m1, float *tmp_dataB,
279  F77_INT p1, FloatMatrix& alpha,
280  FloatMatrix& beta, float *u, F77_INT nrow_u,
281  float *v, F77_INT nrow_v, float *q,
282  F77_INT nrow_q, FloatMatrix& work,
283  F77_INT lwork, F77_INT *iwork, F77_INT& info)
284  {
285  if (! gsvd_initialized)
286  initialize_gsvd ();
287 
288  if (have_DGGSVD3)
289  {
290  sggsvd3_type f_ptr = reinterpret_cast<sggsvd3_type> (gsvd_fcn["sg"]);
291  f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
292  F77_CONST_CHAR_ARG2 (&jobv, 1),
293  F77_CONST_CHAR_ARG2 (&jobq, 1),
294  m, n, p, k, l, tmp_dataA, m1, tmp_dataB, p1,
295  alpha.fortran_vec (), beta.fortran_vec (),
296  u, nrow_u, v, nrow_v, q, nrow_q,
297  work.fortran_vec (), lwork, iwork, info
298  F77_CHAR_ARG_LEN (1)
299  F77_CHAR_ARG_LEN (1)
300  F77_CHAR_ARG_LEN (1));
301  }
302  else
303  {
304  sggsvd_type f_ptr = reinterpret_cast<sggsvd_type> (gsvd_fcn["sg"]);
305  f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
306  F77_CONST_CHAR_ARG2 (&jobv, 1),
307  F77_CONST_CHAR_ARG2 (&jobq, 1),
308  m, n, p, k, l, tmp_dataA, m1, tmp_dataB, p1,
309  alpha.fortran_vec (), beta.fortran_vec (),
310  u, nrow_u, v, nrow_v, q, nrow_q,
311  work.fortran_vec (), iwork, info
312  F77_CHAR_ARG_LEN (1)
313  F77_CHAR_ARG_LEN (1)
314  F77_CHAR_ARG_LEN (1));
315  }
316  }
317 
318  template <>
319  void
320  gsvd<ComplexMatrix>::ggsvd (char& jobu, char& jobv, char& jobq,
321  F77_INT m, F77_INT n, F77_INT p, F77_INT& k,
322  F77_INT& l, Complex *tmp_dataA, F77_INT m1,
323  Complex *tmp_dataB, F77_INT p1, Matrix& alpha,
324  Matrix& beta, Complex *u, F77_INT nrow_u,
325  Complex *v, F77_INT nrow_v, Complex *q,
326  F77_INT nrow_q, ComplexMatrix& work,
327  F77_INT lwork, F77_INT *iwork, F77_INT& info)
328  {
329  if (! gsvd_initialized)
330  initialize_gsvd ();
331 
332  Matrix rwork(2*n, 1);
333  if (have_DGGSVD3)
334  {
335  zggsvd3_type f_ptr = reinterpret_cast<zggsvd3_type> (gsvd_fcn["zg"]);
336  f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
337  F77_CONST_CHAR_ARG2 (&jobv, 1),
338  F77_CONST_CHAR_ARG2 (&jobq, 1),
339  m, n, p, k, l,
340  F77_DBLE_CMPLX_ARG (tmp_dataA), m1,
341  F77_DBLE_CMPLX_ARG (tmp_dataB), p1,
342  alpha.fortran_vec (), beta.fortran_vec (),
343  F77_DBLE_CMPLX_ARG (u), nrow_u,
344  F77_DBLE_CMPLX_ARG (v), nrow_v,
345  F77_DBLE_CMPLX_ARG (q), nrow_q,
346  F77_DBLE_CMPLX_ARG (work.fortran_vec ()),
347  lwork, rwork.fortran_vec (), iwork, info
348  F77_CHAR_ARG_LEN (1)
349  F77_CHAR_ARG_LEN (1)
350  F77_CHAR_ARG_LEN (1));
351  }
352  else
353  {
354  zggsvd_type f_ptr = reinterpret_cast<zggsvd_type> (gsvd_fcn["zg"]);
355  f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
356  F77_CONST_CHAR_ARG2 (&jobv, 1),
357  F77_CONST_CHAR_ARG2 (&jobq, 1),
358  m, n, p, k, l,
359  F77_DBLE_CMPLX_ARG (tmp_dataA), m1,
360  F77_DBLE_CMPLX_ARG (tmp_dataB), p1,
361  alpha.fortran_vec (), beta.fortran_vec (),
362  F77_DBLE_CMPLX_ARG (u), nrow_u,
363  F77_DBLE_CMPLX_ARG (v), nrow_v,
364  F77_DBLE_CMPLX_ARG (q), nrow_q,
365  F77_DBLE_CMPLX_ARG (work.fortran_vec ()),
366  rwork.fortran_vec (), iwork, info
367  F77_CHAR_ARG_LEN (1)
368  F77_CHAR_ARG_LEN (1)
369  F77_CHAR_ARG_LEN (1));
370  }
371  }
372 
373  template <>
374  void
375  gsvd<FloatComplexMatrix>::ggsvd (char& jobu, char& jobv, char& jobq,
376  F77_INT m, F77_INT n, F77_INT p,
377  F77_INT& k, F77_INT& l,
378  FloatComplex *tmp_dataA, F77_INT m1,
379  FloatComplex *tmp_dataB, F77_INT p1,
380  FloatMatrix& alpha, FloatMatrix& beta,
381  FloatComplex *u, F77_INT nrow_u,
382  FloatComplex *v, F77_INT nrow_v,
383  FloatComplex *q, F77_INT nrow_q,
384  FloatComplexMatrix& work, F77_INT lwork,
385  F77_INT *iwork, F77_INT& info)
386  {
387  if (! gsvd_initialized)
388  initialize_gsvd ();
389 
390  FloatMatrix rwork(2*n, 1);
391  if (have_DGGSVD3)
392  {
393  cggsvd3_type f_ptr = reinterpret_cast<cggsvd3_type> (gsvd_fcn["cg"]);
394  f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
395  F77_CONST_CHAR_ARG2 (&jobv, 1),
396  F77_CONST_CHAR_ARG2 (&jobq, 1),
397  m, n, p, k, l,
398  F77_CMPLX_ARG (tmp_dataA), m1,
399  F77_CMPLX_ARG (tmp_dataB), p1,
400  alpha.fortran_vec (), beta.fortran_vec (),
401  F77_CMPLX_ARG (u), nrow_u,
402  F77_CMPLX_ARG (v), nrow_v,
403  F77_CMPLX_ARG (q), nrow_q,
404  F77_CMPLX_ARG (work.fortran_vec ()), lwork,
405  rwork.fortran_vec (), iwork, info
406  F77_CHAR_ARG_LEN (1)
407  F77_CHAR_ARG_LEN (1)
408  F77_CHAR_ARG_LEN (1));
409  }
410  else
411  {
412  cggsvd_type f_ptr = reinterpret_cast<cggsvd_type> (gsvd_fcn["cg"]);
413  f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
414  F77_CONST_CHAR_ARG2 (&jobv, 1),
415  F77_CONST_CHAR_ARG2 (&jobq, 1),
416  m, n, p, k, l,
417  F77_CMPLX_ARG (tmp_dataA), m1,
418  F77_CMPLX_ARG (tmp_dataB), p1,
419  alpha.fortran_vec (), beta.fortran_vec (),
420  F77_CMPLX_ARG (u), nrow_u,
421  F77_CMPLX_ARG (v), nrow_v,
422  F77_CMPLX_ARG (q), nrow_q,
423  F77_CMPLX_ARG (work.fortran_vec ()),
424  rwork.fortran_vec (), iwork, info
425  F77_CHAR_ARG_LEN (1)
426  F77_CHAR_ARG_LEN (1)
427  F77_CHAR_ARG_LEN (1));
428  }
429  }
430 
431  template <typename T>
432  T
433  gsvd<T>::left_singular_matrix_A (void) const
434  {
435  if (type == gsvd::Type::sigma_only)
436  {
437  (*current_liboctave_error_handler)
438  ("gsvd: U not computed because type == gsvd::sigma_only");
439  return T ();
440  }
441  else
442  return left_smA;
443  }
444 
445  template <typename T>
446  T
447  gsvd<T>::left_singular_matrix_B (void) const
448  {
449  if (type == gsvd::Type::sigma_only)
450  {
451  (*current_liboctave_error_handler)
452  ("gsvd: V not computed because type == gsvd::sigma_only");
453  return T ();
454  }
455  else
456  return left_smB;
457  }
458 
459  template <typename T>
460  T
461  gsvd<T>::right_singular_matrix (void) const
462  {
463  if (type == gsvd::Type::sigma_only)
464  {
465  (*current_liboctave_error_handler)
466  ("gsvd: X not computed because type == gsvd::sigma_only");
467  return T ();
468  }
469  else
470  return right_sm;
471  }
472 
473  template <typename T>
474  T
475  gsvd<T>::R_matrix (void) const
476  {
477  if (type != gsvd::Type::std)
478  {
479  (*current_liboctave_error_handler)
480  ("gsvd: R not computed because type != gsvd::std");
481  return T ();
482  }
483  else
484  return R;
485  }
486 
487  template <typename T>
488  gsvd<T>::gsvd (const T& a, const T& b, gsvd::Type gsvd_type)
489  {
490  F77_INT info;
491 
492  F77_INT m = to_f77_int (a.rows ());
493  F77_INT n = to_f77_int (a.cols ());
494  F77_INT p = to_f77_int (b.rows ());
495 
496  T atmp = a;
497  P *tmp_dataA = atmp.fortran_vec ();
498 
499  T btmp = b;
500  P *tmp_dataB = btmp.fortran_vec ();
501 
502  char jobu = 'U';
503  char jobv = 'V';
504  char jobq = 'Q';
505 
506  F77_INT nrow_u = m;
507  F77_INT nrow_v = p;
508  F77_INT nrow_q = n;
509 
510  F77_INT k, l;
511 
512  switch (gsvd_type)
513  {
514  case gsvd<T>::Type::sigma_only:
515 
516  // Note: for this case, both jobu and jobv should be 'N', but
517  // there seems to be a bug in dgesvd from Lapack V2.0. To
518  // demonstrate the bug, set both jobu and jobv to 'N' and find
519  // the singular values of [eye(3), eye(3)]. The result is
520  // [-sqrt(2), -sqrt(2), -sqrt(2)].
521  //
522  // For Lapack 3.0, this problem seems to be fixed.
523 
524  jobu = 'N';
525  jobv = 'N';
526  jobq = 'N';
527  nrow_u = nrow_v = nrow_q = 1;
528  break;
529 
530  default:
531  break;
532  }
533 
534  type = gsvd_type;
535 
536  if (! (jobu == 'N' || jobu == 'O'))
537  left_smA.resize (nrow_u, m);
538 
539  P *u = left_smA.fortran_vec ();
540 
541  if (! (jobv == 'N' || jobv == 'O'))
542  left_smB.resize (nrow_v, p);
543 
544  P *v = left_smB.fortran_vec ();
545 
546  if (! (jobq == 'N' || jobq == 'O'))
547  right_sm.resize (nrow_q, n);
548 
549  P *q = right_sm.fortran_vec ();
550 
551  real_matrix alpha (n, 1);
552  real_matrix beta (n, 1);
553 
554  std::vector<F77_INT> iwork (n);
555 
556  if (! gsvd_initialized)
557  initialize_gsvd ();
558 
559  F77_INT lwork;
560  if (have_DGGSVD3)
561  {
562  lwork = -1;
563  T work_tmp (1, 1);
564 
565  gsvd<T>::ggsvd (jobu, jobv, jobq, m, n, p, k, l,
566  tmp_dataA, m, tmp_dataB, p, alpha, beta, u,
567  nrow_u, v, nrow_v, q, nrow_q, work_tmp, lwork,
568  iwork.data (), info);
569 
570  lwork = static_cast<F77_INT> (std::abs (work_tmp(0, 0)));
571  }
572  else
573  {
574  lwork = 3*n;
575  lwork = (lwork > m ? lwork : m);
576  lwork = (lwork > p ? lwork : p) + n;
577  }
578  info = 0;
579 
580  T work (lwork, 1);
581 
582  gsvd<T>::ggsvd (jobu, jobv, jobq, m, n, p, k, l,
583  tmp_dataA, m, tmp_dataB, p, alpha, beta, u,
584  nrow_u, v, nrow_v, q, nrow_q, work, lwork, iwork.data (),
585  info);
586 
587  if (info < 0)
588  (*current_liboctave_error_handler) ("*ggsvd.f: argument %d illegal",
589  -info);
590  else
591  {
592  if (info > 0)
593  (*current_liboctave_error_handler)
594  ("*ggsvd.f: Jacobi-type procedure failed to converge.");
595  else
596  {
597  F77_INT i, j;
598 
599  if (gsvd<T>::Type::std == gsvd_type)
600  {
601  R.resize(k+l, k+l);
602  int astart = n-k-l;
603  if (m - k - l >= 0)
604  {
605  astart = n-k-l;
606  // R is stored in A(1:K+L,N-K-L+1:N)
607  for (i = 0; i < k+l; i++)
608  for (j = 0; j < k+l; j++)
609  R.xelem (i, j) = atmp.xelem (i, astart + j);
610  }
611  else
612  {
613  // (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N),
614  // ( 0 R22 R23 )
615 
616  for (i = 0; i < m; i++)
617  for (j = 0; j < k+l; j++)
618  R.xelem (i, j) = atmp.xelem (i, astart + j);
619  // and R33 is stored in B(M-K+1:L,N+M-K-L+1:N)
620  for (i = k+l-1; i >=m; i--)
621  {
622  for (j = 0; j < m; j++)
623  R.xelem(i, j) = 0.0;
624  for (j = m; j < k+l; j++)
625  R.xelem (i, j) = btmp.xelem (i - k, astart + j);
626  }
627  }
628  }
629 
630  if (m-k-l >= 0)
631  {
632  // Fills in C and S
633  sigmaA.resize (l, l);
634  sigmaB.resize (l, l);
635  for (i = 0; i < l; i++)
636  {
637  sigmaA.dgxelem(i) = alpha.elem(k+i);
638  sigmaB.dgxelem(i) = beta.elem(k+i);
639  }
640  }
641  else
642  {
643  // Fills in C and S
644  sigmaA.resize (m-k, m-k);
645  sigmaB.resize (m-k, m-k);
646  for (i = 0; i < m-k; i++)
647  {
648  sigmaA.dgxelem(i) = alpha.elem(k+i);
649  sigmaB.dgxelem(i) = beta.elem(k+i);
650  }
651  }
652  }
653  }
654  }
655 
656  // Instantiations we need.
657  template class gsvd<Matrix>;
658  template class gsvd<FloatMatrix>;
659  template class gsvd<ComplexMatrix>;
660  template class gsvd<FloatComplexMatrix>;
661  }
662 }
Array::fortran_vec
const T * fortran_vec(void) const
Size of the specified dimension.
Definition: Array.h:583
Matrix
Definition: dMatrix.h:42
octave::dynamic_library
Definition: oct-shlib.h:44
octave::dynamic_library::search
void * search(const std::string &nm, const name_mangler &mangler=name_mangler()) const
Definition: oct-shlib.h:175
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, T &work, octave_f77_int_type lwork, octave_f77_int_type *iwork, octave_f77_int_type &info)
octave::math::gsvd::real_matrix
T::real_matrix_type real_matrix
Definition: gsvd.h:92
octave::math::gsvd::right_singular_matrix
T right_singular_matrix(void) const
Definition: gsvd.cc:461
octave::math::gsvd::R_matrix
T R_matrix(void) const
Definition: gsvd.cc:475
octave::math::gsvd::P
T::value_type P
Definition: gsvd.h:91
octave::math::gsvd::left_singular_matrix_A
T left_singular_matrix_A(void) const
Definition: gsvd.cc:433
octave::math::gsvd::left_singular_matrix_B
T left_singular_matrix_B(void) const
Definition: gsvd.cc:447
octave::math::gsvd::gsvd
gsvd(void)
Definition: gsvd.h:48
F77_DBLE_CMPLX_ARG
#define F77_DBLE_CMPLX_ARG(x)
Definition: f77-fcn.h:315
F77_CMPLX_ARG
#define F77_CMPLX_ARG(x)
Definition: f77-fcn.h:309
F77_INT
octave_f77_int_type F77_INT
Definition: f77-fcn.h:305
gsvd_type
static octave::math::gsvd< T >::Type gsvd_type(int nargout)
Definition: gsvd.cc:46
STRINGIZE
#define STRINGIZE(x)
Definition: gsvd.cc:52
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
octave::jit_convention::type
type
Definition: jit-typeinfo.h:124
octave::initialize_gsvd
static void initialize_gsvd(void)
Definition: gsvd.cc:54
octave::gsvd_initialized
static bool gsvd_initialized
Definition: gsvd.cc:48
octave::have_DGGSVD3
static bool have_DGGSVD3
Definition: gsvd.cc:47
octave::gsvd_fcn
static std::map< std::string, void * > gsvd_fcn
Definition: gsvd.cc:45
octave::zggsvd_type
comp_ggsvd_ptr< F77_DBLE_CMPLX, F77_DBLE >::type zggsvd_type
Definition: gsvd.cc:225
octave::sggsvd_type
real_ggsvd_ptr< F77_REAL >::type sggsvd_type
Definition: gsvd.cc:223
octave::zggsvd3_type
comp_ggsvd3_ptr< F77_DBLE_CMPLX, F77_DBLE >::type zggsvd3_type
Definition: gsvd.cc:224
octave::dggsvd3_type
real_ggsvd3_ptr< F77_DBLE >::type dggsvd3_type
Definition: gsvd.cc:220
octave::cggsvd_type
comp_ggsvd_ptr< F77_CMPLX, F77_REAL >::type cggsvd_type
Definition: gsvd.cc:227
octave::dggsvd_type
real_ggsvd_ptr< F77_DBLE >::type dggsvd_type
Definition: gsvd.cc:221
octave::cggsvd3_type
comp_ggsvd3_ptr< F77_CMPLX, F77_REAL >::type cggsvd3_type
Definition: gsvd.cc:226
octave::sggsvd3_type
real_ggsvd3_ptr< F77_REAL >::type sggsvd3_type
Definition: gsvd.cc:222
Complex
std::complex< double > Complex
Definition: oct-cmplx.h:33
FloatComplex
std::complex< float > FloatComplex
Definition: oct-cmplx.h:34
abs
static T abs(T x)
Definition: pr-output.cc:1678
octave::comp_ggsvd3_ptr
Definition: gsvd.cc:187
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:189
octave::comp_ggsvd_ptr
Definition: gsvd.cc:154
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:156
octave::real_ggsvd3_ptr
Definition: gsvd.cc:121
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:123
octave::real_ggsvd_ptr
Definition: gsvd.cc:89
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:91
F77_FUNC
F77_RET_T F77_FUNC(xerbla, XERBLA)(F77_CONST_CHAR_ARG_DEF(s_arg