# HG changeset patch # User Jaroslav Hajek # Date 1218275908 -7200 # Node ID 17c72510835710d2dd8589528a48a3bcd407ebe5 # Parent f0fbf47c914c61d63b45bda33339ffdf7b3c7b0b improve norm computation capabilities diff --git a/liboctave/CSparse.cc b/liboctave/CSparse.cc --- a/liboctave/CSparse.cc +++ b/liboctave/CSparse.cc @@ -518,6 +518,37 @@ } return result; +} + +ComplexRowVector +SparseComplexMatrix::row (octave_idx_type i) const +{ + octave_idx_type nc = columns (); + ComplexRowVector retval (nc, 0); + + for (octave_idx_type j = 0; j < nc; j++) + for (octave_idx_type k = cidx (j); k < cidx (j+1); k++) + { + if (ridx (k) == i) + { + retval(j) = data (k); + break; + } + } + + return retval; +} + +ComplexColumnVector +SparseComplexMatrix::column (octave_idx_type i) const +{ + octave_idx_type nr = rows (); + ComplexColumnVector retval (nr); + + for (octave_idx_type k = cidx (i); k < cidx (i+1); k++) + retval(ridx (k)) = data (k); + + return retval; } // destructive insert/delete/reorder operations diff --git a/liboctave/CSparse.h b/liboctave/CSparse.h --- a/liboctave/CSparse.h +++ b/liboctave/CSparse.h @@ -126,6 +126,12 @@ { return MSparse::transpose (); } friend SparseComplexMatrix conj (const SparseComplexMatrix& a); + + // extract row or column i. + + ComplexRowVector row (octave_idx_type i) const; + + ComplexColumnVector column (octave_idx_type i) const; private: SparseComplexMatrix dinverse (MatrixType &mattyp, octave_idx_type& info, diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog --- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,20 @@ +2008-08-10 Jaroslav Hajek + + * oct-norm.h: New header file. + * oct-norm.cc: New source. + * CSparse.cc (SparseComplexMatrix::row, SparseComplexMatrix::column): + New member functions. + * CSparse.h (SparseComplexMatrix): Declare them. + * dSparse.cc (SparseMatrix::row, SparseMatrix::column): + New member functions. + * dSparse.h (SparseMatrix): Declare them. + * MArray-C.cc (MArray::norm), + MArray-d.cc (MArray::norm), + MArray-fC.cc (MArray::norm), + MArray-f.cc (MArray::norm): Wrap a call to xnorm. + + * MArray-defs.h (MARRAY_NORM_BODY): Remove. + 2008-08-11 Jaroslav Hajek * Array.cc (no_op_fcn): New static function. diff --git a/liboctave/MArray-C.cc b/liboctave/MArray-C.cc --- a/liboctave/MArray-C.cc +++ b/liboctave/MArray-C.cc @@ -28,23 +28,17 @@ // Instantiate MArrays of Complex values. #include "oct-cmplx.h" -#include "f77-fcn.h" - -extern "C" -{ - F77_RET_T - F77_FUNC (xdznrm2, XDZNRM2) (const octave_idx_type&, const Complex*, - const octave_idx_type&, double&); -} #include "MArray.h" #include "MArray.cc" +#include "CColVector.h" +#include "oct-norm.h" template <> OCTAVE_API double MArray::norm (double p) const { - MARRAY_NORM_BODY (Complex, double, xdznrm2, XDZNRM2, octave_NaN); + return xnorm (ComplexColumnVector (*this), p); } template class OCTAVE_API MArray; diff --git a/liboctave/MArray-d.cc b/liboctave/MArray-d.cc --- a/liboctave/MArray-d.cc +++ b/liboctave/MArray-d.cc @@ -26,23 +26,16 @@ // Instantiate MArrays of double values. -#include "f77-fcn.h" - -extern "C" -{ - F77_RET_T - F77_FUNC (xdnrm2, XDNRM2) (const octave_idx_type&, const double*, - const octave_idx_type&, double&); -} - #include "MArray.h" #include "MArray.cc" +#include "dColVector.h" +#include "oct-norm.h" template <> OCTAVE_API double MArray::norm (double p) const { - MARRAY_NORM_BODY (double, double, xdnrm2, XDNRM2, octave_NaN); + return xnorm (ColumnVector (*this), p); } template class OCTAVE_API MArray; diff --git a/liboctave/MArray-defs.h b/liboctave/MArray-defs.h --- a/liboctave/MArray-defs.h +++ b/liboctave/MArray-defs.h @@ -343,90 +343,6 @@ MDIAGARRAY2_DADA_BINOP_FWD_DEFS \ (R, T, dynamic_cast&>, R, dynamic_cast&>, R) -#define MARRAY_NORM_BODY(TYPE, RTYPE, blas_norm, BLAS_NORM, NAN_VALUE) \ - \ - RTYPE retval = NAN_VALUE; \ - \ - octave_idx_type len = length (); \ - \ - if (len > 0) \ - { \ - const TYPE *d = data (); \ - \ - if (p == -1) \ - { \ - /* Frobenius norm. */ \ - retval = 0; \ - \ - /* precondition */ \ - RTYPE inf_norm = 0.; \ - for (octave_idx_type i = 0; i < len; i++) \ - { \ - RTYPE d_abs = std::abs (d[i]); \ - if (d_abs > inf_norm) \ - inf_norm = d_abs; \ - } \ - inf_norm = (inf_norm == octave_Inf || inf_norm == 0. ? 1.0 : \ - inf_norm); \ - RTYPE scale = 1. / inf_norm; \ -\ - for (octave_idx_type i = 0; i < len; i++) \ - { \ - RTYPE d_abs = std::abs (d[i]) * scale; \ - retval += d_abs * d_abs; \ - } \ - \ - retval = ::sqrt (retval) * inf_norm; \ - } \ - else if (p == 2) \ - F77_FCN (blas_norm, BLAS_NORM) (len, d, 1, retval); \ - else if (xisinf (p)) \ - { \ - octave_idx_type i = 0; \ - \ - while (i < len && xisnan (d[i])) \ - i++; \ - \ - if (i < len) \ - retval = std::abs (d[i]); \ - \ - if (p > 0) \ - { \ - while (i < len) \ - { \ - RTYPE d_abs = std::abs (d[i++]); \ - \ - if (d_abs > retval) \ - retval = d_abs; \ - } \ - } \ - else \ - { \ - while (i < len) \ - { \ - RTYPE d_abs = std::abs (d[i++]); \ - \ - if (d_abs < retval) \ - retval = d_abs; \ - } \ - } \ - } \ - else \ - { \ - retval = 0; \ - \ - for (octave_idx_type i = 0; i < len; i++) \ - { \ - RTYPE d_abs = std::abs (d[i]); \ - retval += pow (d_abs, p); \ - } \ - \ - retval = pow (retval, 1/p); \ - } \ - } \ - \ - return retval - // Now we have all the definitions we need. #endif diff --git a/liboctave/MArray-f.cc b/liboctave/MArray-f.cc --- a/liboctave/MArray-f.cc +++ b/liboctave/MArray-f.cc @@ -26,23 +26,16 @@ // Instantiate MArrays of float values. -#include "f77-fcn.h" - -extern "C" -{ - F77_RET_T - F77_FUNC (xsnrm2, XSNRM2) (const octave_idx_type&, const float*, - const octave_idx_type&, float&); -} - #include "MArray.h" #include "MArray.cc" +#include "fColVector.h" +#include "oct-norm.h" template <> OCTAVE_API float MArray::norm (float p) const { - MARRAY_NORM_BODY (float, float, xsnrm2, XSNRM2, octave_Float_NaN); + return xnorm (FloatColumnVector (*this)); } template class OCTAVE_API MArray; diff --git a/liboctave/MArray-fC.cc b/liboctave/MArray-fC.cc --- a/liboctave/MArray-fC.cc +++ b/liboctave/MArray-fC.cc @@ -28,23 +28,17 @@ // Instantiate MArrays of FloatComplex values. #include "oct-cmplx.h" -#include "f77-fcn.h" - -extern "C" -{ - F77_RET_T - F77_FUNC (xscnrm2, XSCNRM2) (const octave_idx_type&, const FloatComplex*, - const octave_idx_type&, float&); -} #include "MArray.h" #include "MArray.cc" +#include "fCColVector.h" +#include "oct-norm.h" template <> OCTAVE_API float MArray::norm (float p) const { - MARRAY_NORM_BODY (FloatComplex, float, xscnrm2, XSCNRM2, octave_Float_NaN); + return xnorm (FloatComplexColumnVector (*this)); } template class OCTAVE_API MArray; diff --git a/liboctave/Makefile.in b/liboctave/Makefile.in --- a/liboctave/Makefile.in +++ b/liboctave/Makefile.in @@ -88,9 +88,9 @@ lo-ieee.h lo-mappers.h lo-math.h lo-specfun.h \ lo-sysdep.h lo-utils.h mach-info.h md5.h oct-alloc.h oct-cmplx.h \ oct-env.h oct-fftw.h oct-getopt.h oct-group.h oct-inttypes.h \ - oct-lookup.h oct-md5.h oct-mutex.h oct-passwd.h oct-rand.h oct-rl-edit.h \ - oct-rl-hist.h oct-shlib.h oct-sort.h oct-spparms.h oct-syscalls.h \ - oct-sparse.h oct-time.h oct-uname.h \ + oct-lookup.h oct-md5.h oct-mutex.h oct-norm.h oct-passwd.h \ + oct-rand.h oct-rl-edit.h oct-rl-hist.h oct-shlib.h oct-sort.h \ + oct-spparms.h oct-syscalls.h oct-sparse.h oct-time.h oct-uname.h \ pathlen.h pathsearch.h prog-args.h \ randgamma.h randmtzig.h randpoisson.h regex-match.h \ sparse-sort.h statdefs.h str-vec.h \ @@ -148,7 +148,8 @@ file-ops.cc file-stat.cc glob-match.cc idx-vector.cc \ lo-ieee.cc lo-mappers.cc lo-specfun.cc lo-sysdep.cc \ lo-utils.cc mach-info.cc oct-alloc.cc oct-env.cc \ - oct-fftw.cc oct-group.cc oct-mutex.cc oct-md5.cc oct-passwd.cc oct-rand.cc \ + oct-fftw.cc oct-group.cc oct-mutex.cc oct-md5.cc \ + oct-norm.cc oct-passwd.cc oct-rand.cc \ oct-shlib.cc oct-spparms.cc oct-syscalls.cc oct-time.cc oct-uname.cc \ prog-args.cc regex-match.cc \ sparse-sort.cc sparse-util.cc str-vec.cc \ diff --git a/liboctave/dSparse.cc b/liboctave/dSparse.cc --- a/liboctave/dSparse.cc +++ b/liboctave/dSparse.cc @@ -513,6 +513,37 @@ } return result; +} + +RowVector +SparseMatrix::row (octave_idx_type i) const +{ + octave_idx_type nc = columns (); + RowVector retval (nc, 0); + + for (octave_idx_type j = 0; j < nc; j++) + for (octave_idx_type k = cidx (j); k < cidx (j+1); k++) + { + if (ridx (k) == i) + { + retval(j) = data (k); + break; + } + } + + return retval; +} + +ColumnVector +SparseMatrix::column (octave_idx_type i) const +{ + octave_idx_type nr = rows (); + ColumnVector retval (nr); + + for (octave_idx_type k = cidx (i); k < cidx (i+1); k++) + retval(ridx (k)) = data (k); + + return retval; } SparseMatrix diff --git a/liboctave/dSparse.h b/liboctave/dSparse.h --- a/liboctave/dSparse.h +++ b/liboctave/dSparse.h @@ -120,6 +120,12 @@ return MSparse::transpose (); } SparseMatrix hermitian (void) const { return transpose (); } + + // extract row or column i. + + RowVector row (octave_idx_type i) const; + + ColumnVector column (octave_idx_type i) const; private: SparseMatrix dinverse (MatrixType &mattyp, octave_idx_type& info, diff --git a/liboctave/oct-norm.cc b/liboctave/oct-norm.cc new file mode 100644 --- /dev/null +++ b/liboctave/oct-norm.cc @@ -0,0 +1,510 @@ +/* + +Copyright (C) 2008 VZLU Prague, a.s. + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, see +. + +*/ + +// author: Jaroslav Hajek + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "oct-cmplx.h" +#include "lo-error.h" +#include "lo-ieee.h" +#include "Array.h" +#include "Array.cc" +#include "Array-util.h" +#include "CMatrix.h" +#include "dMatrix.h" +#include "fCMatrix.h" +#include "fMatrix.h" +#include "CColVector.h" +#include "dColVector.h" +#include "CRowVector.h" +#include "dRowVector.h" +#include "fCColVector.h" +#include "fColVector.h" +#include "fCRowVector.h" +#include "fRowVector.h" +#include "CSparse.h" +#include "dSparse.h" +#include "dbleSVD.h" +#include "CmplxSVD.h" +#include "floatSVD.h" +#include "fCmplxSVD.h" + +// Theory: norm accumulator is an object that has an accum method able to handle +// both real and complex element, and a cast operator returning the intermediate +// norm. + +// norm accumulator for the p-norm +template +class norm_accumulator_p +{ + R p,scl,sum; +public: + norm_accumulator_p () {} // we need this one for Array + norm_accumulator_p (R pp) : p(pp), scl(0), sum(1) {} + + template + void accum (U val) + { + R t = std::abs (val); + if (scl == t) // we need this to handle Infs properly + sum += 1; + else if (scl < t) + { + sum *= std::pow (scl/t, p); + sum += 1; + scl = t; + } + else if (t != 0) + sum += std::pow (t/scl, p); + } + operator R () { return scl * std::pow (sum, 1/p); } +}; + +// norm accumulator for the minus p-pseudonorm +template +class norm_accumulator_mp +{ + R p,scl,sum; +public: + norm_accumulator_mp () {} // we need this one for Array + norm_accumulator_mp (R pp) : p(pp), scl(0), sum(1) {} + + template + void accum (U val) + { + R t = 1 / std::abs (val); + if (scl == t) + sum += 1; + else if (scl < t) + { + sum *= std::pow (scl/t, p); + sum += 1; + scl = t; + } + else if (t != 0) + sum += std::pow (t/scl, p); + } + operator R () { return scl * std::pow (sum, -1/p); } +}; + +// norm accumulator for the 2-norm (euclidean) +template +class norm_accumulator_2 +{ + R scl,sum; + static R pow2 (R x) { return x*x; } +public: + norm_accumulator_2 () : scl(0), sum(1) {} + + void accum (R val) + { + R t = std::abs (val); + if (scl == t) + sum += 1; + else if (scl < t) + { + sum *= pow2 (scl/t); + sum += 1; + scl = t; + } + else if (t != 0) + sum += pow2 (t/scl); + } + + void accum (std::complex val) + { + accum (val.real ()); + accum (val.imag ()); + } + + operator R () { return scl * std::sqrt (sum); } +}; + +// norm accumulator for the 1-norm (city metric) +template +class norm_accumulator_1 +{ + R sum; +public: + norm_accumulator_1 () : sum (0) {} + template + void accum (U val) + { + sum += std::abs (val); + } + operator R () { return sum; } +}; + +// norm accumulator for the inf-norm (max metric) +template +class norm_accumulator_inf +{ + R max; +public: + norm_accumulator_inf () : max (0) {} + template + void accum (U val) + { + max = std::max (max, std::abs (val)); + } + operator R () { return max; } +}; + +// norm accumulator for the -inf pseudonorm (min abs value) +template +class norm_accumulator_minf +{ + R min; +public: + norm_accumulator_minf () : min (octave_Inf) {} + template + void accum (U val) + { + min = std::min (min, std::abs (val)); + } + operator R () { return min; } +}; + +// norm accumulator for the 0-pseudonorm (hamming distance) +template +class norm_accumulator_0 +{ + unsigned int num; +public: + norm_accumulator_0 () : num (0) {} + template + void accum (U val) + { + if (val != static_cast (0)) ++num; + } + operator R () { return num; } +}; + + +// OK, we're armed :) Now let's go for the fun + +template +inline void vector_norm (const Array& v, R& res, ACC acc) +{ + for (octave_idx_type i = 0; i < v.numel (); i++) + acc.accum (v(i)); + + res = acc; +} + +// dense versions +template +void column_norms (const MArray2& m, MArray& res, ACC acc) +{ + res.resize (m.columns ()); + for (octave_idx_type j = 0; j < m.columns (); j++) + { + ACC accj = acc; + for (octave_idx_type i = 0; i < m.rows (); i++) + accj.accum (m(i, j)); + + res.xelem (j) = accj; + } +} + +template +void row_norms (const MArray2& m, MArray& res, ACC acc) +{ + res.resize (m.rows ()); + Array acci (m.rows (), acc); + for (octave_idx_type j = 0; j < m.columns (); j++) + { + for (octave_idx_type i = 0; i < m.rows (); i++) + acci.xelem (i).accum (m(i, j)); + } + + for (octave_idx_type i = 0; i < m.rows (); i++) + res.xelem (i) = acci(i); +} + +// sparse versions +template +void column_norms (const MSparse& m, MArray& res, ACC acc) +{ + res.resize (m.columns ()); + for (octave_idx_type j = 0; j < m.columns (); j++) + { + ACC accj = acc; + for (octave_idx_type k = m.cidx (j); k < m.cidx (j+1); k++) + accj.accum (m.data (k)); + + res.xelem (j) = accj; + } +} + +template +void row_norms (const MSparse& m, MArray& res, ACC acc) +{ + res.resize (m.rows ()); + Array acci (m.rows (), acc); + for (octave_idx_type j = 0; j < m.columns (); j++) + { + for (octave_idx_type k = m.cidx (j); k < m.cidx (j+1); k++) + acci(m.ridx (k)).accum (m.data (k)); + } + + for (octave_idx_type i = 0; i < m.rows (); i++) + res.xelem (i) = acci(i); +} + +// now the dispatchers +#define DEFINE_DISPATCHER(FUNC_NAME, ARG_TYPE, RES_TYPE) \ +template \ +RES_TYPE FUNC_NAME (const ARG_TYPE& v, R p) \ +{ \ + RES_TYPE res; \ + if (p == 2) \ + FUNC_NAME (v, res, norm_accumulator_2 ()); \ + else if (p == 1) \ + FUNC_NAME (v, res, norm_accumulator_1 ()); \ + else if (lo_ieee_isinf (p)) \ + { \ + if (p > 0) \ + FUNC_NAME (v, res, norm_accumulator_inf ()); \ + else \ + FUNC_NAME (v, res, norm_accumulator_minf ()); \ + } \ + else if (p == 0) \ + FUNC_NAME (v, res, norm_accumulator_0 ()); \ + else if (p > 0) \ + FUNC_NAME (v, res, norm_accumulator_p (p)); \ + else \ + FUNC_NAME (v, res, norm_accumulator_mp (p)); \ + return res; \ +} + +DEFINE_DISPATCHER (vector_norm, MArray, R) +DEFINE_DISPATCHER (vector_norm, MArray2, R) +DEFINE_DISPATCHER (column_norms, MArray2, MArray) +DEFINE_DISPATCHER (row_norms, MArray2, MArray) +DEFINE_DISPATCHER (column_norms, MSparse, MArray) +DEFINE_DISPATCHER (row_norms, MSparse, MArray) + +// the subproblem in Higham's method +// TODO: works for complex matrices? +template +void higham_subp (const ColVectorT& y, const ColVectorT& col, + octave_idx_type nsamp, R p, R& lambda, R& mu) +{ + R nrm = 0; + for (octave_idx_type i = 0; i < nsamp; i++) + { + R fi = i*M_PI/nsamp, lambda1 = cos (fi), mu1 = sin (fi); + R lmnr = std::pow (std::pow (std::abs (lambda1), p) + + std::pow (std::abs (mu1), p), 1/p); + lambda1 /= lmnr; mu1 /= lmnr; + R nrm1 = vector_norm (lambda1 * y + mu1 * col, p); + if (nrm1 > nrm) + { + lambda = lambda1; + mu = mu1; + nrm = nrm1; + } + } +} + +// the p-dual element (should work for both real and complex) +template +inline T elem_dual_p (T x, R p) +{ + return conj (signum (x)) * std::pow (std::abs (x), p-1); +} + +// the VectorT is used for vectors, but actually it has to be +// a Matrix type to allow all the operations. For instance SparseMatrix +// does not support multiplication with column/row vectors. +// the dual vector +template +VectorT dual_p (const VectorT& x, R p, R q) +{ + VectorT res (x.dims ()); + for (octave_idx_type i = 0; i < x.numel (); i++) + res.xelem (i) = elem_dual_p (x(i), p); + return res / vector_norm (res, q); +} + +// Higham's hybrid method +template +R higham (const MatrixT& m, R p, R tol, int maxiter, + VectorT& x) +{ + x.resize (m.columns (), 1); + // the OSE part + VectorT y(m.rows (), 1, 0), z(m.rows (), 1); + R lambda = 0, mu = 0; + for (octave_idx_type k = 0; k < m.columns (); k++) + { + VectorT col (m.column (k)); + if (k > 0) + higham_subp (y, col, 4*k, p, lambda, mu); + for (octave_idx_type i = 0; i < k; i++) + x(i) *= lambda; + x(k) = mu; + y = lambda * y + mu * col; + } + + // the PM part + x = x / vector_norm (x, p); + R q = p/(p-1); + + R gamma = 0, gamma1; + int iter = 0; + while (iter < maxiter) + { + y = m*x; + gamma1 = gamma; + gamma = vector_norm (y, p); + z = dual_p (y, p, q); + z = z.hermitian (); + z = z * m; + + if (iter > 0 && (vector_norm (z, q) <= std::abs ((z*x)(0)) + || (gamma - gamma1) <= tol*gamma)) + break; + + z = z.hermitian (); + x = dual_p (z, q, p); + iter ++; + } + + return gamma; +} + +// TODO: is there a better way? +inline float get_eps (float) { return FLT_EPSILON; } +inline double get_eps (double) { return DBL_EPSILON; } +// derive column vector and SVD types + +static const char *p_less1_gripe = "xnorm: p must be at least 1"; + +// version with SVD for dense matrices +template +R matrix_norm (const MatrixT& m, R p, VectorT, SVDT) +{ + R res = 0; + if (p == 2) + { + SVDT svd (m, SVD::sigma_only); + res = svd.singular_values () (0,0); + } + else if (p == 1) + res = xcolnorms (m, 1).max (); + else if (lo_ieee_isinf (p)) + res = xrownorms (m, 1).max (); + else if (p > 1) + { + VectorT x; + res = higham (m, p, std::sqrt (get_eps (p)), 100, x); + } + else + (*current_liboctave_error_handler) (p_less1_gripe); + + return res; +} + +// SVD-free version for sparse matrices +template +R matrix_norm (const MatrixT& m, R p, VectorT) +{ + R res = 0; + if (p == 1) + res = xcolnorms (m, 1).max (); + else if (lo_ieee_isinf (p)) + res = xrownorms (m, 1).max (); + else if (p > 1) + { + VectorT x; + res = higham (m, p, std::sqrt (get_eps (p)), 100, x); + } + else + (*current_liboctave_error_handler) (p_less1_gripe); + + return res; +} + +// and finally, here's what we've promised in the header file + +#define DEFINE_XNORM_FUNCS(PREFIX, RTYPE) \ + RTYPE xnorm (const PREFIX##ColumnVector& x, RTYPE p) \ + { return vector_norm (x, p); } \ + RTYPE xnorm (const PREFIX##RowVector& x, RTYPE p) \ + { return vector_norm (x, p); } \ + RTYPE xnorm (const PREFIX##Matrix& x, RTYPE p) \ + { return matrix_norm (x, p, PREFIX##Matrix (), PREFIX##SVD ()); } \ + RTYPE xfrobnorm (const PREFIX##Matrix& x) \ + { return vector_norm (x, static_cast (2)); } + +DEFINE_XNORM_FUNCS(, double) +DEFINE_XNORM_FUNCS(Complex, double) +DEFINE_XNORM_FUNCS(Float, float) +DEFINE_XNORM_FUNCS(FloatComplex, float) + +// this is needed to avoid copying the sparse matrix for xfrobnorm +template +inline void array_norm_2 (const T* v, octave_idx_type n, R& res) +{ + norm_accumulator_2 acc; + for (octave_idx_type i = 0; i < n; i++) + acc.accum (v[i]); + + res = acc; +} + +#define DEFINE_XNORM_SPARSE_FUNCS(PREFIX, RTYPE) \ + RTYPE xnorm (const Sparse##PREFIX##Matrix& x, RTYPE p) \ + { return matrix_norm (x, p, PREFIX##Matrix ()); } \ + RTYPE xfrobnorm (const Sparse##PREFIX##Matrix& x) \ + { \ + RTYPE res; \ + array_norm_2 (x.data (), x.nnz (), res); \ + return res; \ + } + +DEFINE_XNORM_SPARSE_FUNCS(, double) +DEFINE_XNORM_SPARSE_FUNCS(Complex, double) + +#define DEFINE_COLROW_NORM_FUNCS(PREFIX, RPREFIX, RTYPE) \ + extern RPREFIX##RowVector xcolnorms (const PREFIX##Matrix& m, RTYPE p) \ + { return column_norms (m, p); } \ + extern RPREFIX##ColumnVector xrownorms (const PREFIX##Matrix& m, RTYPE p) \ + { return row_norms (m, p); } \ + +DEFINE_COLROW_NORM_FUNCS(, , double) +DEFINE_COLROW_NORM_FUNCS(Complex, , double) +DEFINE_COLROW_NORM_FUNCS(Float, Float, float) +DEFINE_COLROW_NORM_FUNCS(FloatComplex, Float, float) + +DEFINE_COLROW_NORM_FUNCS(Sparse, , double) +DEFINE_COLROW_NORM_FUNCS(SparseComplex, , double) + diff --git a/liboctave/oct-norm.h b/liboctave/oct-norm.h new file mode 100644 --- /dev/null +++ b/liboctave/oct-norm.h @@ -0,0 +1,60 @@ +/* + +Copyright (C) 2008 VZLU Prague, a.s. + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, see +. + +*/ + +// author: Jaroslav Hajek + +#if !defined (octave_xnorm_h) +#define octave_xnorm_h 1 + +#include "oct-cmplx.h" + +#define DECLARE_XNORM_FUNCS(PREFIX, RTYPE) \ + class PREFIX##Matrix; \ + class PREFIX##ColumnVector; \ + class PREFIX##RowVector; \ + \ + extern RTYPE xnorm (const PREFIX##ColumnVector&, RTYPE p = 2); \ + extern RTYPE xnorm (const PREFIX##RowVector&, RTYPE p = 2); \ + extern RTYPE xnorm (const PREFIX##Matrix&, RTYPE p = 2); \ + extern RTYPE xfrobnorm (const PREFIX##Matrix&); + +DECLARE_XNORM_FUNCS(, double) +DECLARE_XNORM_FUNCS(Complex, double) +DECLARE_XNORM_FUNCS(Float, float) +DECLARE_XNORM_FUNCS(FloatComplex, float) + +DECLARE_XNORM_FUNCS(Sparse, double) +DECLARE_XNORM_FUNCS(SparseComplex, double) + +#define DECLARE_COLROW_NORM_FUNCS(PREFIX, RPREFIX, RTYPE) \ + extern RPREFIX##RowVector xcolnorms (const PREFIX##Matrix&, RTYPE p = 2); \ + extern RPREFIX##ColumnVector xrownorms (const PREFIX##Matrix&, RTYPE p = 2); \ + +DECLARE_COLROW_NORM_FUNCS(, , double) +DECLARE_COLROW_NORM_FUNCS(Complex, , double) +DECLARE_COLROW_NORM_FUNCS(Float, Float, float) +DECLARE_COLROW_NORM_FUNCS(FloatComplex, Float, float) + +DECLARE_COLROW_NORM_FUNCS(Sparse, , double) +DECLARE_COLROW_NORM_FUNCS(SparseComplex, , double) + +#endif diff --git a/src/ChangeLog b/src/ChangeLog --- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,3 +1,12 @@ +2008-08-10 Jaroslav Hajek + + * xnorm.cc: New source. + * xnorm.h: New header file. + * Makefile.in: Include xnorm.cc in the build process. + * data.cc (Fnorm): Call xnorm, xcolnorms, xrownorms or xfrobnorm + to do the actual work. + + 2008-08-07 John W. Eaton * ov.cc (octave_value::idx_type_value): Don't include default diff --git a/src/Makefile.in b/src/Makefile.in --- a/src/Makefile.in +++ b/src/Makefile.in @@ -128,7 +128,7 @@ pager.h parse.h pr-output.h procstream.h sighandlers.h \ siglist.h sparse-xdiv.h sparse-xpow.h symtab.h sysdep.h \ token.h toplev.h unwind-prot.h utils.h variables.h \ - version.h xdiv.h xpow.h \ + version.h xdiv.h xnorm.h xpow.h \ $(OV_INCLUDES) \ $(PT_INCLUDES) \ $(OV_SPARSE_INCLUDES) @@ -207,7 +207,7 @@ parse.y pr-output.cc procstream.cc sighandlers.cc \ siglist.c sparse-xdiv.cc sparse-xpow.cc strfns.cc \ syscalls.cc symtab.cc sysdep.cc token.cc toplev.cc \ - unwind-prot.cc utils.cc variables.cc xdiv.cc xpow.cc \ + unwind-prot.cc utils.cc variables.cc xdiv.cc xnorm.cc xpow.cc \ $(OV_SRC) \ $(PT_SRC) diff --git a/src/data.cc b/src/data.cc --- a/src/data.cc +++ b/src/data.cc @@ -62,6 +62,7 @@ #include "utils.h" #include "variables.h" #include "pager.h" +#include "xnorm.h" #if ! defined (HAVE_HYPOTF) && defined (HAVE__HYPOTF) #define hypotf _hypotf @@ -4464,6 +4465,10 @@ Frobenius norm of @var{a}, @code{sqrt (sum (diag (@var{a}' * @var{a})))}.\n\ @end table\n\ \n\ address@hidden otherwise\n\ address@hidden general p-norm, i.e. maximum @code{norm (A*x, p)} such that \n\ address@hidden (x, p) == 1}\n\ +\n\ If @var{a} is a vector or a scalar:\n\ \n\ @table @asis\n\ @@ -4482,14 +4487,11 @@ @seealso{cond, svd}\n\ @end deftypefn") { - // Currently only handles vector norms for full double/complex - // vectors internally. Other cases are handled by __norm__.m. - octave_value_list retval; int nargin = args.length (); - if (nargin == 1 || nargin == 2) + if (nargin >= 1 && nargin <= 3) { octave_value x_arg = args(0); @@ -4502,103 +4504,54 @@ } else if (x_arg.ndims () == 2) { - if ((x_arg.rows () == 1 || x_arg.columns () == 1) - && ! (x_arg.is_sparse_type () || x_arg.is_integer_type ())) - { - double p_val = 2; - - if (nargin == 2) - { - octave_value p_arg = args(1); - - if (p_arg.is_string ()) - { - std::string p = args(1).string_value (); - - if (p == "inf") - p_val = octave_Inf; - else if (p == "fro") - p_val = -1; - else - error ("norm: unrecognized norm `%s'", p.c_str ()); - } - else - { - p_val = p_arg.double_value (); - - if (error_state) - error ("norm: unrecognized norm value"); - } - } - - if (x_arg.is_single_type ()) - { - if (! error_state) - { - if (x_arg.is_real_type ()) - { - MArray x (x_arg.float_array_value ()); - - if (! error_state) - retval(0) = x.norm (static_cast(p_val)); - else - error ("norm: expecting real vector"); - } - else - { - MArray x (x_arg.float_complex_array_value ()); - - if (! error_state) - retval(0) = x.norm (static_cast(p_val)); - else - error ("norm: expecting complex vector"); - } - } - } - else - { - if (! error_state) - { - if (x_arg.is_real_type ()) - { - MArray x (x_arg.array_value ()); - - if (! error_state) - retval(0) = x.norm (p_val); - else - error ("norm: expecting real vector"); - } - else - { - MArray x (x_arg.complex_array_value ()); - - if (! error_state) - retval(0) = x.norm (p_val); - else - error ("norm: expecting complex vector"); - } - } - } - } - else - retval = feval ("__norm__", args); + enum { sfmatrix, sfcols, sfrows, sffrob, sfinf } strflag = sfmatrix; + if (nargin > 1 && args(nargin-1).is_string ()) + { + std::string str = args(nargin-1).string_value (); + if (str == "cols" || str == "columns") + strflag = sfcols; + else if (str == "rows") + strflag = sfrows; + else if (str == "fro") + strflag = sffrob; + else if (str == "inf") + strflag = sfinf; + else + error ("norm: unrecognized option: %s", str.c_str ()); + // we've handled the last parameter, so act as if it was removed + nargin --; + } + else if (nargin > 1 && ! args(1).is_scalar_type ()) + gripe_wrong_type_arg ("norm", args(1), true); + + if (! error_state) + { + octave_value p_arg = (nargin > 1) ? args(1) : octave_value (2); + switch (strflag) + { + case sfmatrix: + retval(0) = xnorm (x_arg, p_arg); + break; + case sfcols: + retval(0) = xcolnorms (x_arg, p_arg); + break; + case sfrows: + retval(0) = xrownorms (x_arg, p_arg); + break; + case sffrob: + retval(0) = xfrobnorm (x_arg); + break; + case sfinf: + retval(0) = xnorm (x_arg, octave_Inf); + break; + } + } } else error ("norm: only valid for 2-D objects"); } else print_usage (); - - // Should not return a sparse type - if (retval(0).is_sparse_type ()) - { - if (retval(0).type_name () == "sparse matrix") - retval(0) = retval(0).matrix_value (); - else if (retval(0).type_name () == "sparse complex matrix") - retval(0) = retval(0).complex_matrix_value (); - else if (retval(0).type_name () == "sparse bool matrix") - retval(0) = retval(0).bool_matrix_value (); - } return retval; } diff --git a/src/xnorm.cc b/src/xnorm.cc new file mode 100644 --- /dev/null +++ b/src/xnorm.cc @@ -0,0 +1,212 @@ +/* + +Copyright (C) 2008 VZLU Prague, a.s. + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, see +. + +*/ + +// author: Jaroslav Hajek + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "oct-norm.h" + +#include "error.h" +#include "xnorm.h" +#include "ov.h" +#include "gripes.h" + +octave_value xnorm (const octave_value& x, const octave_value& p) +{ + octave_value retval; + + bool isvector = (x.columns () == 1 || x.rows () == 1); + bool iscomplex = x.is_complex_type (); + bool issparse = x.is_sparse_type (); + bool isfloat = x.is_single_type (); + + if (isfloat || x.is_double_type ()) + { + if (isvector) + { + if (isfloat & iscomplex) + retval = xnorm (x.float_complex_column_vector_value (), + p.float_value ()); + else if (isfloat) + retval = xnorm (x.float_column_vector_value (), + p.float_value ()); + else if (iscomplex) + retval = xnorm (x.complex_column_vector_value (), + p.double_value ()); + else + retval = xnorm (x.column_vector_value (), + p.double_value ()); + } + else if (issparse) + { + if (iscomplex) + retval = xnorm (x.sparse_complex_matrix_value (), + p.double_value ()); + else + retval = xnorm (x.sparse_matrix_value (), + p.double_value ()); + } + else + { + if (isfloat & iscomplex) + retval = xnorm (x.float_complex_matrix_value (), + p.float_value ()); + else if (isfloat) + retval = xnorm (x.float_matrix_value (), + p.float_value ()); + else if (iscomplex) + retval = xnorm (x.complex_matrix_value (), + p.double_value ()); + else + retval = xnorm (x.matrix_value (), + p.double_value ()); + } + } + else + gripe_wrong_type_arg ("xnorm", x, true); + + return retval; +} + +octave_value xcolnorms (const octave_value& x, const octave_value& p) +{ + octave_value retval; + + bool iscomplex = x.is_complex_type (); + bool issparse = x.is_sparse_type (); + bool isfloat = x.is_single_type (); + + if (isfloat || x.is_double_type ()) + { + if (issparse) + { + if (iscomplex) + retval = xcolnorms (x.sparse_complex_matrix_value (), + p.double_value ()); + else + retval = xcolnorms (x.sparse_matrix_value (), + p.double_value ()); + } + else + { + if (isfloat & iscomplex) + retval = xcolnorms (x.float_complex_matrix_value (), + p.float_value ()); + else if (isfloat) + retval = xcolnorms (x.float_matrix_value (), + p.float_value ()); + else if (iscomplex) + retval = xcolnorms (x.complex_matrix_value (), + p.double_value ()); + else + retval = xcolnorms (x.matrix_value (), + p.double_value ()); + } + } + else + gripe_wrong_type_arg ("xcolnorms", x, true); + + return retval; +} + +octave_value xrownorms (const octave_value& x, const octave_value& p) +{ + octave_value retval; + + bool iscomplex = x.is_complex_type (); + bool issparse = x.is_sparse_type (); + bool isfloat = x.is_single_type (); + + if (isfloat || x.is_double_type ()) + { + if (issparse) + { + if (iscomplex) + retval = xrownorms (x.sparse_complex_matrix_value (), + p.double_value ()); + else + retval = xrownorms (x.sparse_matrix_value (), + p.double_value ()); + } + else + { + if (isfloat & iscomplex) + retval = xrownorms (x.float_complex_matrix_value (), + p.float_value ()); + else if (isfloat) + retval = xrownorms (x.float_matrix_value (), + p.float_value ()); + else if (iscomplex) + retval = xrownorms (x.complex_matrix_value (), + p.double_value ()); + else + retval = xrownorms (x.matrix_value (), + p.double_value ()); + } + } + else + gripe_wrong_type_arg ("xrownorms", x, true); + + return retval; +} + +octave_value xfrobnorm (const octave_value& x) +{ + octave_value retval; + + bool iscomplex = x.is_complex_type (); + bool issparse = x.is_sparse_type (); + bool isfloat = x.is_single_type (); + + if (isfloat || x.is_double_type ()) + { + if (issparse) + { + if (iscomplex) + retval = xfrobnorm (x.sparse_complex_matrix_value ()); + else + retval = xfrobnorm (x.sparse_matrix_value ()); + } + else + { + if (isfloat & iscomplex) + retval = xfrobnorm (x.float_complex_matrix_value ()); + else if (isfloat) + retval = xfrobnorm (x.float_matrix_value ()); + else if (iscomplex) + retval = xfrobnorm (x.complex_matrix_value ()); + else + retval = xfrobnorm (x.matrix_value ()); + } + } + else + gripe_wrong_type_arg ("xfrobnorm", x, true); + + return retval; +} diff --git a/src/xnorm.h b/src/xnorm.h new file mode 100644 --- /dev/null +++ b/src/xnorm.h @@ -0,0 +1,35 @@ +/* + +Copyright (C) 2008 VZLU Prague, a.s. + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, see +. + +*/ + +// author: Jaroslav Hajek + +#if !defined (octave_xnorm_h) +#define octave_xnorm_h 1 + +class octave_value; + +extern octave_value xnorm (const octave_value& x, const octave_value& p); +extern octave_value xcolnorms (const octave_value& x, const octave_value& p); +extern octave_value xrownorms (const octave_value& x, const octave_value& p); +extern octave_value xfrobnorm (const octave_value& x); + +#endif