24 #if defined (HAVE_CONFIG_H) 51 (*current_liboctave_error_handler)
52 (
"svd: U not computed because type == svd::sigma_only");
62 (*current_liboctave_error_handler)
63 (
"svd: V not computed because type == svd::sigma_only");
70 #define GESVD_REAL_STEP(f, F) \ 71 F77_XFCN (f, F, (F77_CONST_CHAR_ARG2 (&jobu, 1), \ 72 F77_CONST_CHAR_ARG2 (&jobv, 1), \ 73 m, n, tmp_data, m1, s_vec, u, m1, vt, \ 74 nrow_vt1, work.data (), lwork, info \ 75 F77_CHAR_ARG_LEN (1) \ 76 F77_CHAR_ARG_LEN (1))) 78 #define GESVD_COMPLEX_STEP(f, F, CMPLX_ARG) \ 79 F77_XFCN (f, F, (F77_CONST_CHAR_ARG2 (&jobu, 1), \ 80 F77_CONST_CHAR_ARG2 (&jobv, 1), \ 81 m, n, CMPLX_ARG (tmp_data), \ 82 m1, s_vec, CMPLX_ARG (u), m1, \ 83 CMPLX_ARG (vt), nrow_vt1, \ 84 CMPLX_ARG (work.data ()), \ 85 lwork, rwork.data (), info \ 86 F77_CHAR_ARG_LEN (1) \ 87 F77_CHAR_ARG_LEN (1))) 93 double *tmp_data,
F77_INT m1,
double *s_vec,
94 double *
u,
double *vt,
F77_INT nrow_vt1,
95 std::vector<double>& work,
F77_INT& lwork,
100 lwork =
static_cast<F77_INT> (work[0]);
101 work.reserve (lwork);
110 float *tmp_data,
F77_INT m1,
float *s_vec,
111 float *
u,
float *vt,
F77_INT nrow_vt1,
112 std::vector<float>& work,
F77_INT& lwork,
117 lwork =
static_cast<F77_INT> (work[0]);
118 work.reserve (lwork);
129 std::vector<Complex>& work,
F77_INT& lwork,
132 std::vector<double> rwork (5 *
std::max (m, n));
136 lwork =
static_cast<F77_INT> (work[0].real ());
137 work.reserve (lwork);
149 std::vector<FloatComplex>& work,
152 std::vector<float> rwork (5 *
std::max (m, n));
156 lwork =
static_cast<F77_INT> (work[0].real ());
157 work.reserve (lwork);
162 #undef GESVD_REAL_STEP 163 #undef GESVD_COMPLEX_STEP 167 #define GESDD_REAL_STEP(f, F) \ 168 F77_XFCN (f, F, (F77_CONST_CHAR_ARG2 (&jobz, 1), \ 169 m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, \ 170 work.data (), lwork, iwork, info \ 171 F77_CHAR_ARG_LEN (1))) 173 #define GESDD_COMPLEX_STEP(f, F, CMPLX_ARG) \ 174 F77_XFCN (f, F, (F77_CONST_CHAR_ARG2 (&jobz, 1), m, n, \ 175 CMPLX_ARG (tmp_data), m1, \ 176 s_vec, CMPLX_ARG (u), m1, \ 177 CMPLX_ARG (vt), nrow_vt1, \ 178 CMPLX_ARG (work.data ()), lwork, \ 179 rwork.data (), iwork, info \ 180 F77_CHAR_ARG_LEN (1))) 186 F77_INT m1,
double *s_vec,
double *
u,
double *vt,
187 F77_INT nrow_vt1, std::vector<double>& work,
192 lwork =
static_cast<F77_INT> (work[0]);
193 work.reserve (lwork);
202 F77_INT m1,
float *s_vec,
float *
u,
float *vt,
203 F77_INT nrow_vt1, std::vector<float>& work,
208 lwork =
static_cast<F77_INT> (work[0]);
209 work.reserve (lwork);
220 std::vector<Complex>& work,
F77_INT& lwork,
231 lrwork = min_mn *
std::max (5*min_mn+5, 2*max_mn+2*min_mn+1);
233 std::vector<double> rwork (lrwork);
237 lwork =
static_cast<F77_INT> (work[0].real ());
238 work.reserve (lwork);
250 std::vector<FloatComplex>& work,
261 lrwork = min_mn *
std::max (5*min_mn+5, 2*max_mn+2*min_mn+1);
262 std::vector<float> rwork (lrwork);
266 lwork =
static_cast<F77_INT> (work[0].real ());
267 work.reserve (lwork);
272 #undef GESDD_REAL_STEP 273 #undef GESDD_COMPLEX_STEP 278 : m_type (
type), m_driver (driver), left_sm (), sigma (), right_sm ()
285 if (m == 0 || n == 0)
314 P *tmp_data = atmp.fortran_vec ();
316 F77_INT min_mn = (m < n ? m : n);
330 ncol_u = nrow_vt = nrow_s = ncol_s = min_mn;
343 ncol_u = nrow_vt = 1;
350 if (! (jobu ==
'N' || jobu ==
'O'))
355 sigma.resize (nrow_s, ncol_s);
358 if (! (jobv ==
'N' || jobv ==
'O'))
367 std::vector<P> work (1);
373 gesvd (jobu, jobv, m, n, tmp_data, m1, s_vec,
u, vt, nrow_vt1,
377 assert (jobu == jobv);
380 std::vector<F77_INT> iwork (8 *
std::min (m, n));
382 gesdd (jobz, m, n, tmp_data, m1, s_vec,
u, vt, nrow_vt1,
383 work, lwork, iwork.data (), info);
388 if (! (jobv ==
'N' || jobv ==
'O'))
OCTAVE_NORETURN liboctave_error_handler current_liboctave_error_handler
#define F77_DBLE_CMPLX_ARG(x)
#define GESDD_COMPLEX_STEP(f, F, CMPLX_ARG)
#define GESVD_COMPLEX_STEP(f, F, CMPLX_ARG)
T right_singular_matrix(void) const
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
T::real_diag_matrix_type DM_T
#define GESDD_REAL_STEP(f, F)
void gesdd(char &jobz, octave_f77_int_type m, octave_f77_int_type n, P *tmp_data, octave_f77_int_type m1, DM_P *s_vec, P *u, P *vt, octave_f77_int_type nrow_vt1, std::vector< P > &work, octave_f77_int_type &lwork, octave_f77_int_type *iwork, octave_f77_int_type &info)
void gesvd(char &jobu, char &jobv, octave_f77_int_type m, octave_f77_int_type n, P *tmp_data, octave_f77_int_type m1, DM_P *s_vec, P *u, P *vt, octave_f77_int_type nrow_vt1, std::vector< P > &work, octave_f77_int_type &lwork, octave_f77_int_type &info)
charNDArray max(char d, const charNDArray &m)
octave_f77_int_type F77_INT
#define GESVD_REAL_STEP(f, F)
std::complex< float > FloatComplex
std::complex< double > Complex
charNDArray min(char d, const charNDArray &m)
T left_singular_matrix(void) const