94 operator R () {
return scl *
std::pow (sum, 1/p); }
122 operator R () {
return scl *
std::pow (sum, -1/p); }
155 operator R () {
return scl * std::sqrt (sum); }
170 operator R () {
return sum; }
185 operator R () {
return max; }
200 operator R () {
return min; }
213 if (val != static_cast<U> (0)) ++
num;
215 operator R () {
return num; }
221 template <
class T,
class R,
class ACC>
231 template <
class T,
class R,
class ACC>
239 accj.accum (m(i, j));
241 res.
xelem (j) = accj;
245 template <
class T,
class R,
class ACC>
249 std::vector<ACC> acci (m.
rows (), acc);
253 acci[i].accum (m(i, j));
257 res.
xelem (i) = acci[i];
261 template <
class T,
class R,
class ACC>
269 accj.accum (m.
data (k));
271 res.
xelem (j) = accj;
275 template <
class T,
class R,
class ACC>
279 std::vector<ACC> acci (m.
rows (), acc);
283 acci[m.
ridx (k)].accum (m.
data (k));
287 res.
xelem (i) = acci[i];
291 #define DEFINE_DISPATCHER(FUNC_NAME, ARG_TYPE, RES_TYPE) \
292 template <class T, class R> \
293 RES_TYPE FUNC_NAME (const ARG_TYPE& v, R p) \
297 FUNC_NAME (v, res, norm_accumulator_2<R> ()); \
299 FUNC_NAME (v, res, norm_accumulator_1<R> ()); \
300 else if (lo_ieee_isinf (p)) \
303 FUNC_NAME (v, res, norm_accumulator_inf<R> ()); \
305 FUNC_NAME (v, res, norm_accumulator_minf<R> ()); \
308 FUNC_NAME (v, res, norm_accumulator_0<R> ()); \
310 FUNC_NAME (v, res, norm_accumulator_p<R> (p)); \
312 FUNC_NAME (v, res, norm_accumulator_mp<R> (p)); \
325 template <class ColVectorT, class R>
334 R
fi = i *
static_cast<R
> (M_PI) / nsamp;
335 R lambda1 = cos (fi);
339 lambda1 /= lmnr; mu1 /= lmnr;
353 template <
class ColVectorT,
class R>
357 std::complex<R>& lambda, std::complex<R>& mu)
359 typedef std::complex<R> CR;
362 CR lamcu = lambda /
std::abs (lambda);
367 R
fi = i *
static_cast<R
> (M_PI) / nsamp;
368 R lambda1 = cos (fi);
372 lambda1 /= lmnr; mu1 /= lmnr;
373 R nrm1 =
vector_norm (lambda1 * lamcu * y + mu1 * col, p);
376 lambda = lambda1 * lamcu;
386 R
fi = i *
static_cast<R
> (M_PI) / nsamp;
387 lamcu = CR (cos (fi), sin (fi));
388 R nrm1 =
vector_norm (lama * lamcu * y + mu * col, p);
391 lambda = lama * lamcu;
398 template <
class T,
class R>
408 template <
class VectorT,
class R>
411 VectorT res (x.dims ());
418 template <
class MatrixT,
class VectorT,
class R>
419 R
higham (
const MatrixT& m, R p, R tol,
int maxiter,
422 x.resize (m.columns (), 1);
424 VectorT y(m.rows (), 1, 0), z(m.rows (), 1);
425 typedef typename VectorT::element_type RR;
431 VectorT col (m.column (k));
437 y = lambda * y + mu * col;
446 while (iter < maxiter)
457 || (gamma - gamma1) <= tol*gamma))
478 template <
class MatrixT,
class VectorT,
class SVDT,
class R>
485 res = svd.singular_values () (0,0);
494 const R sqrteps = std::sqrt (std::numeric_limits<R>::epsilon ());
495 res =
higham (m, p, sqrteps, max_norm_iter, x);
504 template <
class MatrixT,
class VectorT,
class R>
515 const R sqrteps = std::sqrt (std::numeric_limits<R>::epsilon ());
516 res =
higham (m, p, sqrteps, max_norm_iter, x);
526 #define DEFINE_XNORM_FUNCS(PREFIX, RTYPE) \
527 OCTAVE_API RTYPE xnorm (const PREFIX##ColumnVector& x, RTYPE p) \
528 { return vector_norm (x, p); } \
529 OCTAVE_API RTYPE xnorm (const PREFIX##RowVector& x, RTYPE p) \
530 { return vector_norm (x, p); } \
531 OCTAVE_API RTYPE xnorm (const PREFIX##Matrix& x, RTYPE p) \
532 { return matrix_norm (x, p, PREFIX##Matrix (), PREFIX##SVD ()); } \
533 OCTAVE_API RTYPE xfrobnorm (const PREFIX##Matrix& x) \
534 { return vector_norm (x, static_cast<RTYPE> (2)); }
542 template <class T, class R>
552 #define DEFINE_XNORM_SPARSE_FUNCS(PREFIX, RTYPE) \
553 OCTAVE_API RTYPE xnorm (const Sparse##PREFIX##Matrix& x, RTYPE p) \
554 { return matrix_norm (x, p, PREFIX##Matrix ()); } \
555 OCTAVE_API RTYPE xfrobnorm (const Sparse##PREFIX##Matrix& x) \
558 array_norm_2 (x.data (), x.nnz (), res); \
565 #define DEFINE_COLROW_NORM_FUNCS(PREFIX, RPREFIX, RTYPE) \
566 extern OCTAVE_API RPREFIX##RowVector xcolnorms (const PREFIX##Matrix& m, RTYPE p) \
567 { return column_norms (m, p); } \
568 extern OCTAVE_API RPREFIX##ColumnVector xrownorms (const PREFIX##Matrix& m, RTYPE p) \
569 { return row_norms (m, p); } \
OCTAVE_API ColumnVector xrownorms(const Matrix &m, double p)
#define DEFINE_XNORM_FUNCS(PREFIX, RTYPE)
octave_idx_type rows(void) const
#define DEFINE_DISPATCHER(FUNC_NAME, ARG_TYPE, RES_TYPE)
octave_idx_type numel(void) const
Number of elements in the array.
VectorT dual_p(const VectorT &x, R p, R q)
octave_idx_type * cidx(void)
octave_idx_type columns(void) const
R matrix_norm(const MatrixT &m, R p, VectorT, SVDT)
octave_idx_type rows(void) const
liboctave_error_handler current_liboctave_error_handler
static void higham_subp(const ColVectorT &y, const ColVectorT &col, octave_idx_type nsamp, R p, R &lambda, R &mu)
#define DEFINE_COLROW_NORM_FUNCS(PREFIX, RPREFIX, RTYPE)
OCTAVE_API RowVector xcolnorms(const Matrix &m, double p)
octave_int< T > pow(const octave_int< T > &a, const octave_int< T > &b)
static const char * p_less1_gripe
void accum(std::complex< R > val)
T & xelem(octave_idx_type n)
Handles the reference counting for all the derived classes.
charNDArray max(char d, const charNDArray &m)
octave_idx_type * ridx(void)
#define DEFINE_XNORM_SPARSE_FUNCS(PREFIX, RTYPE)
norm_accumulator_mp(R pp)
void array_norm_2(const T *v, octave_idx_type n, R &res)
void row_norms(const MArray< T > &m, MArray< R > &res, ACC acc)
void vector_norm(const Array< T > &v, R &res, ACC acc)
void column_norms(const MArray< T > &m, MArray< R > &res, ACC acc)
std::complex< float > FloatComplex
R higham(const MatrixT &m, R p, R tol, int maxiter, VectorT &x)
std::complex< double > Complex
octave_idx_type columns(void) const
F77_RET_T const double * x
charNDArray min(char d, const charNDArray &m)