48 #define xSTRINGIZE(x) #x 49 #define STRINGIZE(x) xSTRINGIZE(x) 149 template<
class T1,
class T2>
182 template<
class T1,
class T2>
234 double *tmp_dataA,
F77_INT m1,
double *tmp_dataB,
246 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
247 F77_CONST_CHAR_ARG2 (&jobv, 1),
248 F77_CONST_CHAR_ARG2 (&jobq, 1),
249 m, n,
p,
k, l, tmp_dataA, m1, tmp_dataB, p1,
251 u, nrow_u, v, nrow_v, q, nrow_q,
255 F77_CHAR_ARG_LEN (1));
260 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
261 F77_CONST_CHAR_ARG2 (&jobv, 1),
262 F77_CONST_CHAR_ARG2 (&jobq, 1),
263 m, n,
p,
k, l, tmp_dataA, m1, tmp_dataB, p1,
265 u, nrow_u, v, nrow_v, q, nrow_q,
269 F77_CHAR_ARG_LEN (1));
277 float *tmp_dataA,
F77_INT m1,
float *tmp_dataB,
280 float *v,
F77_INT nrow_v,
float *q,
290 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
291 F77_CONST_CHAR_ARG2 (&jobv, 1),
292 F77_CONST_CHAR_ARG2 (&jobq, 1),
293 m, n,
p,
k, l, tmp_dataA, m1, tmp_dataB, p1,
295 u, nrow_u, v, nrow_v, q, nrow_q,
299 F77_CHAR_ARG_LEN (1));
304 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
305 F77_CONST_CHAR_ARG2 (&jobv, 1),
306 F77_CONST_CHAR_ARG2 (&jobq, 1),
307 m, n,
p,
k, l, tmp_dataA, m1, tmp_dataB, p1,
309 u, nrow_u, v, nrow_v, q, nrow_q,
313 F77_CHAR_ARG_LEN (1));
335 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
336 F77_CONST_CHAR_ARG2 (&jobv, 1),
337 F77_CONST_CHAR_ARG2 (&jobq, 1),
349 F77_CHAR_ARG_LEN (1));
354 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
355 F77_CONST_CHAR_ARG2 (&jobv, 1),
356 F77_CONST_CHAR_ARG2 (&jobq, 1),
368 F77_CHAR_ARG_LEN (1));
393 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
394 F77_CONST_CHAR_ARG2 (&jobv, 1),
395 F77_CONST_CHAR_ARG2 (&jobq, 1),
407 F77_CHAR_ARG_LEN (1));
412 f_ptr (F77_CONST_CHAR_ARG2 (&jobu, 1),
413 F77_CONST_CHAR_ARG2 (&jobv, 1),
414 F77_CONST_CHAR_ARG2 (&jobq, 1),
426 F77_CHAR_ARG_LEN (1));
430 template <
typename T>
436 (*current_liboctave_error_handler)
437 (
"gsvd: U not computed because type == gsvd::sigma_only");
444 template <
typename T>
450 (*current_liboctave_error_handler)
451 (
"gsvd: V not computed because type == gsvd::sigma_only");
458 template <
typename T>
464 (*current_liboctave_error_handler)
465 (
"gsvd: X not computed because type == gsvd::sigma_only");
472 template <
typename T>
478 (*current_liboctave_error_handler)
479 (
"gsvd: R not computed because type != gsvd::std");
486 template <
typename T>
496 P *tmp_dataA = atmp.fortran_vec ();
499 P *tmp_dataB = btmp.fortran_vec ();
526 nrow_u = nrow_v = nrow_q = 1;
535 if (! (jobu ==
'N' || jobu ==
'O'))
536 left_smA.resize (nrow_u, m);
538 P *
u = left_smA.fortran_vec ();
540 if (! (jobv ==
'N' || jobv ==
'O'))
541 left_smB.resize (nrow_v,
p);
543 P *v = left_smB.fortran_vec ();
545 if (! (jobq ==
'N' || jobq ==
'O'))
546 right_sm.resize (nrow_q, n);
548 P *q = right_sm.fortran_vec ();
553 std::vector<F77_INT> iwork (n);
565 tmp_dataA, m, tmp_dataB,
p, alpha, beta,
u,
566 nrow_u, v, nrow_v, q, nrow_q, work_tmp, lwork,
567 iwork.data (), info);
574 lwork = (lwork > m ? lwork : m);
575 lwork = (lwork >
p ? lwork :
p) + n;
582 tmp_dataA, m, tmp_dataB,
p, alpha, beta,
u,
583 nrow_u, v, nrow_v, q, nrow_q, work, lwork, iwork.data (),
587 (*current_liboctave_error_handler) (
"*ggsvd.f: argument %d illegal",
592 (*current_liboctave_error_handler)
593 (
"*ggsvd.f: Jacobi-type procedure failed to converge.");
606 for (
i = 0;
i <
k+l;
i++)
607 for (j = 0; j <
k+l; j++)
608 R.xelem (
i, j) = atmp.xelem (
i, astart + j);
615 for (
i = 0;
i < m;
i++)
616 for (j = 0; j <
k+l; j++)
617 R.xelem (
i, j) = atmp.xelem (
i, astart + j);
619 for (
i =
k+l-1;
i >=m;
i--)
621 for (j = 0; j < m; j++)
623 for (j = m; j <
k+l; j++)
624 R.xelem (
i, j) = btmp.xelem (
i -
k, astart + j);
632 sigmaA.resize (l, l);
633 sigmaB.resize (l, l);
634 for (
i = 0;
i < l;
i++)
636 sigmaA.dgxelem(
i) = alpha.elem(
k+
i);
637 sigmaB.dgxelem(
i) = beta.elem(
k+
i);
643 sigmaA.resize (m-
k, m-
k);
644 sigmaB.resize (m-
k, m-
k);
645 for (
i = 0;
i < m-
k;
i++)
647 sigmaA.dgxelem(
i) = alpha.elem(
k+
i);
648 sigmaB.dgxelem(
i) = beta.elem(
k+
i);
static bool gsvd_initialized
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)
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
#define F77_DBLE_CMPLX_ARG(x)
real_ggsvd3_ptr< F77_REAL >::type sggsvd3_type
comp_ggsvd3_ptr< F77_CMPLX, F77_REAL >::type cggsvd3_type
const T * fortran_vec(void) const
T::real_matrix_type real_matrix
comp_ggsvd_ptr< F77_CMPLX, F77_REAL >::type cggsvd_type
void initialize_gsvd(void)
F77_RET_T F77_CONST_CHAR_ARG_DECL
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)
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)
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
comp_ggsvd3_ptr< F77_DBLE_CMPLX, F77_DBLE >::type zggsvd3_type
static octave::math::gsvd< T >::Type gsvd_type(int nargout)
T left_singular_matrix_B(void) const
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)
F77_RET_T F77_FUNC(xerbla, XERBLA)
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)
octave_f77_int_type F77_INT
T left_singular_matrix_A(void) const
real_ggsvd3_ptr< F77_DBLE >::type dggsvd3_type
T right_singular_matrix(void) const
comp_ggsvd_ptr< F77_DBLE_CMPLX, F77_DBLE >::type zggsvd_type
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)
real_ggsvd_ptr< F77_DBLE >::type dggsvd_type
std::complex< float > FloatComplex
std::complex< double > Complex
real_ggsvd_ptr< F77_REAL >::type sggsvd_type
void * search(const std::string &nm, name_mangler mangler=nullptr) const
static std::map< std::string, void * > gsvd_fcn