Mercurial > octave-nkf
changeset 8303:b11c31849b44
improve norm computation capabilities
| author | Jaroslav Hajek <highegg@gmail.com> |
|---|---|
| date | Fri, 31 Oct 2008 08:05:32 +0100 |
| parents | f2e050b62199 |
| children | eeaee297c0da |
| files | liboctave/CSparse.cc liboctave/CSparse.h liboctave/ChangeLog liboctave/MArray-C.cc liboctave/MArray-d.cc liboctave/MArray-defs.h liboctave/MArray-f.cc liboctave/MArray-fC.cc liboctave/Makefile.in liboctave/dSparse.cc liboctave/dSparse.h liboctave/oct-norm.cc liboctave/oct-norm.h scripts/ChangeLog scripts/linear-algebra/Makefile.in scripts/linear-algebra/__norm__.m src/ChangeLog src/Makefile.in src/data.cc src/xnorm.cc src/xnorm.h |
| diffstat | 21 files changed, 1054 insertions(+), 369 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/CSparse.cc Sun Nov 02 06:44:12 2008 +0100 +++ b/liboctave/CSparse.cc Fri Oct 31 08:05:32 2008 +0100 @@ -520,6 +520,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 SparseComplexMatrix&
--- a/liboctave/CSparse.h Sun Nov 02 06:44:12 2008 +0100 +++ b/liboctave/CSparse.h Fri Oct 31 08:05:32 2008 +0100 @@ -127,6 +127,12 @@ 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, double& rcond, const bool force = false,
--- a/liboctave/ChangeLog Sun Nov 02 06:44:12 2008 +0100 +++ b/liboctave/ChangeLog Fri Oct 31 08:05:32 2008 +0100 @@ -1,3 +1,20 @@ +2008-10-31 Jaroslav Hajek <highegg@gmail.com> + + * 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<Complex>::norm), + MArray-d.cc (MArray<double>::norm), + MArray-fC.cc (MArray<FloatComplex>::norm), + MArray-f.cc (MArray<float>::norm): Wrap a call to xnorm. + + * MArray-defs.h (MARRAY_NORM_BODY): Remove. + 2008-11-02 Jaroslav Hajek <highegg@gmail.com> * idx-vector.cc (idx_vector::is_complement): Set resulting extent
--- a/liboctave/MArray-C.cc Sun Nov 02 06:44:12 2008 +0100 +++ b/liboctave/MArray-C.cc Fri Oct 31 08:05:32 2008 +0100 @@ -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<Complex>::norm (double p) const { - MARRAY_NORM_BODY (Complex, double, xdznrm2, XDZNRM2, octave_NaN); + return xnorm (ComplexColumnVector (*this), p); } template class OCTAVE_API MArray<Complex>;
--- a/liboctave/MArray-d.cc Sun Nov 02 06:44:12 2008 +0100 +++ b/liboctave/MArray-d.cc Fri Oct 31 08:05:32 2008 +0100 @@ -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<double>::norm (double p) const { - MARRAY_NORM_BODY (double, double, xdnrm2, XDNRM2, octave_NaN); + return xnorm (ColumnVector (*this), p); } template class OCTAVE_API MArray<double>;
--- a/liboctave/MArray-defs.h Sun Nov 02 06:44:12 2008 +0100 +++ b/liboctave/MArray-defs.h Fri Oct 31 08:05:32 2008 +0100 @@ -343,90 +343,6 @@ MDIAGARRAY2_DADA_BINOP_FWD_DEFS \ (R, T, dynamic_cast<const B<T>&>, R, dynamic_cast<const B<T>&>, 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
--- a/liboctave/MArray-f.cc Sun Nov 02 06:44:12 2008 +0100 +++ b/liboctave/MArray-f.cc Fri Oct 31 08:05:32 2008 +0100 @@ -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<float>::norm (float p) const { - MARRAY_NORM_BODY (float, float, xsnrm2, XSNRM2, octave_Float_NaN); + return xnorm (FloatColumnVector (*this)); } template class OCTAVE_API MArray<float>;
--- a/liboctave/MArray-fC.cc Sun Nov 02 06:44:12 2008 +0100 +++ b/liboctave/MArray-fC.cc Fri Oct 31 08:05:32 2008 +0100 @@ -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<FloatComplex>::norm (float p) const { - MARRAY_NORM_BODY (FloatComplex, float, xscnrm2, XSCNRM2, octave_Float_NaN); + return xnorm (FloatComplexColumnVector (*this)); } template class OCTAVE_API MArray<FloatComplex>;
--- a/liboctave/Makefile.in Sun Nov 02 06:44:12 2008 +0100 +++ b/liboctave/Makefile.in Fri Oct 31 08:05:32 2008 +0100 @@ -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 \
--- a/liboctave/dSparse.cc Sun Nov 02 06:44:12 2008 +0100 +++ b/liboctave/dSparse.cc Fri Oct 31 08:05:32 2008 +0100 @@ -515,6 +515,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 SparseMatrix::concat (const SparseMatrix& rb, const Array<octave_idx_type>& ra_idx) {
--- a/liboctave/dSparse.h Sun Nov 02 06:44:12 2008 +0100 +++ b/liboctave/dSparse.h Fri Oct 31 08:05:32 2008 +0100 @@ -121,6 +121,12 @@ } 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, double& rcond, const bool force = false,
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/oct-norm.cc Fri Oct 31 08:05:32 2008 +0100 @@ -0,0 +1,562 @@ +/* + +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 +<http://www.gnu.org/licenses/>. + +*/ + +// author: Jaroslav Hajek <highegg@gmail.com> + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <cassert> +#include <cfloat> +#include <cmath> + +#include <iostream> + +#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 R> +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<class U> + 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 R> +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<class U> + 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 R> +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<R> val) + { + accum (val.real ()); + accum (val.imag ()); + } + + operator R () { return scl * std::sqrt (sum); } +}; + +// norm accumulator for the 1-norm (city metric) +template <class R> +class norm_accumulator_1 +{ + R sum; +public: + norm_accumulator_1 () : sum (0) {} + template<class U> + void accum (U val) + { + sum += std::abs (val); + } + operator R () { return sum; } +}; + +// norm accumulator for the inf-norm (max metric) +template <class R> +class norm_accumulator_inf +{ + R max; +public: + norm_accumulator_inf () : max (0) {} + template<class U> + 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 R> +class norm_accumulator_minf +{ + R min; +public: + norm_accumulator_minf () : min (octave_Inf) {} + template<class U> + 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 R> +class norm_accumulator_0 +{ + unsigned int num; +public: + norm_accumulator_0 () : num (0) {} + template<class U> + void accum (U val) + { + if (val != static_cast<U> (0)) ++num; + } + operator R () { return num; } +}; + + +// OK, we're armed :) Now let's go for the fun + +template <class T, class R, class ACC> +inline void vector_norm (const Array<T>& v, R& res, ACC acc) +{ + for (octave_idx_type i = 0; i < v.numel (); i++) + acc.accum (v(i)); + + res = acc; +} + +// dense versions +template <class T, class R, class ACC> +void column_norms (const MArray2<T>& m, MArray<R>& 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 <class T, class R, class ACC> +void row_norms (const MArray2<T>& m, MArray<R>& res, ACC acc) +{ + res.resize (m.rows ()); + Array<ACC> 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 <class T, class R, class ACC> +void column_norms (const MSparse<T>& m, MArray<R>& 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 <class T, class R, class ACC> +void row_norms (const MSparse<T>& m, MArray<R>& res, ACC acc) +{ + res.resize (m.rows ()); + Array<ACC> 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 <class T, class R> \ +RES_TYPE FUNC_NAME (const ARG_TYPE& v, R p) \ +{ \ + RES_TYPE res; \ + if (p == 2) \ + FUNC_NAME (v, res, norm_accumulator_2<R> ()); \ + else if (p == 1) \ + FUNC_NAME (v, res, norm_accumulator_1<R> ()); \ + else if (lo_ieee_isinf (p)) \ + { \ + if (p > 0) \ + FUNC_NAME (v, res, norm_accumulator_inf<R> ()); \ + else \ + FUNC_NAME (v, res, norm_accumulator_minf<R> ()); \ + } \ + else if (p == 0) \ + FUNC_NAME (v, res, norm_accumulator_0<R> ()); \ + else if (p > 0) \ + FUNC_NAME (v, res, norm_accumulator_p<R> (p)); \ + else \ + FUNC_NAME (v, res, norm_accumulator_mp<R> (p)); \ + return res; \ +} + +DEFINE_DISPATCHER (vector_norm, MArray<T>, R) +DEFINE_DISPATCHER (vector_norm, MArray2<T>, R) +DEFINE_DISPATCHER (column_norms, MArray2<T>, MArray<R>) +DEFINE_DISPATCHER (row_norms, MArray2<T>, MArray<R>) +DEFINE_DISPATCHER (column_norms, MSparse<T>, MArray<R>) +DEFINE_DISPATCHER (row_norms, MSparse<T>, MArray<R>) + +// The approximate subproblem in Higham's method. Find lambda and mu such that +// norm ([lambda, mu], p) == 1 and norm (y*lambda + col*mu, p) is maximized. +// Real version. As in Higham's paper. +template <class ColVectorT, class R> +static 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; + } + } +} + +// Complex version. Higham's paper does not deal with complex case, so we use a simple +// extension. First, guess the magnitudes as in real version, then try to rotate lambda +// to improve further. +template <class ColVectorT, class R> +static void +higham_subp (const ColVectorT& y, const ColVectorT& col, + octave_idx_type nsamp, R p, + std::complex<R>& lambda, std::complex<R>& mu) +{ + typedef std::complex<R> CR; + R nrm = 0; + lambda = 1.0; + CR lamcu = lambda / std::abs (lambda); + // Probe magnitudes + 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 * lamcu * y + mu1 * col, p); + if (nrm1 > nrm) + { + lambda = lambda1 * lamcu; + mu = mu1; + nrm = nrm1; + } + } + R lama = std::abs (lambda); + // Probe orientation + for (octave_idx_type i = 0; i < nsamp; i++) + { + R fi = i*M_PI/nsamp; + lamcu = CR (cos (fi), sin (fi)); + R nrm1 = vector_norm (lama * lamcu * y + mu * col, p); + if (nrm1 > nrm) + { + lambda = lama * lamcu; + nrm = nrm1; + } + } +} + +// the p-dual element (should work for both real and complex) +template <class T, class R> +inline T elem_dual_p (T x, R p) +{ + return 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 <class VectorT, class R> +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 <class MatrixT, class VectorT, class R> +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); + typedef typename VectorT::element_type RR; + RR 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) <= gamma + || (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"; + +// Static constant to control the maximum number of iterations. 100 seems to be a good value. +// Eventually, we can provide a means to change this constant from Octave. +static int max_norm_iter = 100; + +// version with SVD for dense matrices +template <class MatrixT, class VectorT, class SVDT, class R> +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)), max_norm_iter, x); + } + else + (*current_liboctave_error_handler) (p_less1_gripe); + + return res; +} + +// SVD-free version for sparse matrices +template <class MatrixT, class VectorT, class R> +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)), max_norm_iter, 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<RTYPE> (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 <class T, class R> +inline void array_norm_2 (const T* v, octave_idx_type n, R& res) +{ + norm_accumulator_2<R> 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) +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/oct-norm.h Fri Oct 31 08:05:32 2008 +0100 @@ -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 +<http://www.gnu.org/licenses/>. + +*/ + +// author: Jaroslav Hajek <highegg@gmail.com> + +#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
--- a/scripts/ChangeLog Sun Nov 02 06:44:12 2008 +0100 +++ b/scripts/ChangeLog Fri Oct 31 08:05:32 2008 +0100 @@ -1,3 +1,7 @@ +2008-10-31 Jaroslav Hajek <highegg@gmail.com> + + * linear-algebra/__norm__.m: Remove. + 2008-10-31 David Bateman <dbateman@free.fr> * plot/__contour__.m: Exclude infinite values when calculating contour
--- a/scripts/linear-algebra/Makefile.in Sun Nov 02 06:44:12 2008 +0100 +++ b/scripts/linear-algebra/Makefile.in Fri Oct 31 08:05:32 2008 +0100 @@ -33,7 +33,7 @@ INSTALL_PROGRAM = @INSTALL_PROGRAM@ INSTALL_DATA = @INSTALL_DATA@ -SOURCES = __norm__.m commutation_matrix.m cond.m condest.m cross.m \ +SOURCES = commutation_matrix.m cond.m condest.m cross.m \ dmult.m dot.m duplication_matrix.m housh.m krylov.m krylovb.m logm.m \ null.m onenormest.m orth.m planerot.m qzhess.m rank.m rref.m subspace.m \ trace.m vec.m vech.m
--- a/scripts/linear-algebra/__norm__.m Sun Nov 02 06:44:12 2008 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,142 +0,0 @@ -## Copyright (C) 1996, 1997, 2007 John W. Eaton -## -## 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 -## <http://www.gnu.org/licenses/>. - -## Undocumented internal function. - -## Author: jwe - -function retval = __norm__ (x, p) - - if (nargin < 1 || nargin > 2) - print_usage (); - endif - - if (isempty (x)) - retval = []; - return; - endif - - if (ndims (x) > 2) - error ("norm: only valid on 2-D objects") - endif - - if (nargin == 1) - p = 2; - endif - - ## Do we have a vector or matrix as the first argument? - if (ndims(x) == 2 && (rows (x) == 1 || columns (x) == 1)) - if (ischar (p)) - if (strcmp (p, "fro")) - inf_norm = norm (x, "inf"); - if (inf_norm && finite (inf_norm)) - retval = inf_norm .* sqrt (sum (abs (x ./ inf_norm) .^ 2)); - else - retval = inf_norm; - endif - elseif (strcmp (p, "inf")) - retval = max (abs (x)); - else - error ("norm: unrecognized norm"); - endif - else - if (p == Inf) - retval = max (abs (x)); - elseif (p == -Inf) - retval = min (abs (x)); - elseif (p == 1) - retval = sum (abs (x)); - elseif (p == 2) - if (iscomplex (x)) - x = abs (x); - endif - retval = sqrt (sumsq (x)); - else - retval = sum (abs (x) .^ p) ^ (1/p); - endif - endif - else - if (ischar (p)) - if (strcmp (p, "fro")) - inf_norm = norm (x, "inf"); - if (inf_norm && finite (inf_norm)) - retval = inf_norm .* sqrt (sum (sum (abs (x ./ inf_norm) .^ 2))); - else - retval = inf_norm; - endif - elseif (strcmp (p, "inf")) - retval = max (sum (abs (x), 2)); - else - error ("norm: unrecognized vector norm"); - endif - else - if (p == 1) - retval = max (sum (abs (x))); - elseif (p == 2) - s = svd (x); - retval = s (1); - elseif (p == Inf) - retval = max (sum (abs (x), 2)); - else - error ("norm: unrecognized matrix norm"); - endif - endif - endif - -endfunction - -%!test -%! assert (__norm__ (magic (3)), 15, -2*eps); -%! assert (__norm__ (magic (3) * i), 15, -2*eps); - -%!test -%! assert (__norm__ (zeros (5), "fro"), 0); -%! assert (__norm__ (ones (5), "fro"), 5); -%! assert (__norm__ (zeros (5,1), "fro"), 0); -%! assert (__norm__ (2*ones (5,3), "fro"), sqrt (60)); - -%!test -%! assert (__norm__ (zeros (5), "inf"), 0); -%! assert (__norm__ (ones (5), "inf"), 5); -%! assert (__norm__ (2*ones (5,1), "inf"), 2); -%! assert (__norm__ (2*ones (5,3), "inf"), 6); - -%!test -%! assert (__norm__ (zeros (5), 1), 0); -%! assert (__norm__ (ones (5), 1), 5); -%! assert (__norm__ (2*ones (1,5), 1), 10); -%! assert (__norm__ (2*ones (3,5), 1), 6); - - -%!test -%! assert (__norm__ (1e304 * ones (5, 3), "fro"), 1e304 * sqrt (15)); -%! assert (__norm__ (1e-320 * ones (5, 3), "fro"), 1e-320 * sqrt (15)); -%! assert (x = __norm__ ([1, 2; 3, Inf], "fro"), Inf); -%! assert (x = __norm__ ([1, 2, 3, Inf], "fro"), Inf); -%! assert (x = __norm__ ([1, -Inf; 3, 4], "fro"), Inf); -%! assert (x = __norm__ ([1, 2; 3, NaN], "fro"), NaN); - -%!test -%! assert (__norm__ (1e304 * ones (5, 3), "inf"), 3e304); -%! assert (__norm__ (1e-320 * ones (5, 3), "inf"), 3e-320); -%! assert (x = __norm__ ([1, 2; 3, Inf], "inf"), Inf); -%! assert (x = __norm__ ([1, 2, 3, Inf], "inf"), Inf); -%! assert (x = __norm__ ([1, -Inf; 3, 4], "inf"), Inf); -%! assert (x = __norm__ ([1, 2; 3, NaN], "inf"), 3); - -
--- a/src/ChangeLog Sun Nov 02 06:44:12 2008 +0100 +++ b/src/ChangeLog Fri Oct 31 08:05:32 2008 +0100 @@ -1,3 +1,12 @@ +2008-10-31 Jaroslav Hajek <highegg@gmail.com> + + * 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-10-30 John W. Eaton <jwe@octave.org> * oct-map.cc (Octave_map::index): Copy key_list.
--- a/src/Makefile.in Sun Nov 02 06:44:12 2008 +0100 +++ b/src/Makefile.in Fri Oct 31 08:05:32 2008 +0100 @@ -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)
--- a/src/data.cc Sun Nov 02 06:44:12 2008 +0100 +++ b/src/data.cc Fri Oct 31 08:05:32 2008 +0100 @@ -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 @@ -4574,11 +4575,11 @@ DEFUN (norm, args, , "-*- texinfo -*-\n\ -@deftypefn {Function File} {} norm (@var{a}, @var{p})\n\ +@deftypefn {Function File} {} norm (@var{a}, @var{p}, @var{opt})\n\ Compute the p-norm of the matrix @var{a}. If the second argument is\n\ missing, @code{p = 2} is assumed.\n\ \n\ -If @var{a} is a matrix:\n\ +If @var{a} is a matrix (or sparse matrix):\n\ \n\ @table @asis\n\ @item @var{p} = @code{1}\n\ @@ -4594,6 +4595,10 @@ @item @var{p} = @code{\"fro\"}\n\ @cindex Frobenius norm\n\ Frobenius norm of @var{a}, @code{sqrt (sum (diag (@var{a}' * @var{a})))}.\n\ +\n\ +@item other @var{p}, @code{@var{p} > 1}\n\ +@cindex general p-norm \n\ +maximum @code{norm (A*x, p)} such that @code{norm (x, p) == 1}\n\ @end table\n\ \n\ If @var{a} is a vector or a scalar:\n\ @@ -4608,20 +4613,27 @@ @item @var{p} = @code{\"fro\"}\n\ Frobenius norm of @var{a}, @code{sqrt (sumsq (abs (a)))}.\n\ \n\ -@item other\n\ +@item @var{p} = 0\n\ +Hamming norm - the number of nonzero elements.\n\ +\n\ +@item other @var{p}, @code{@var{p} > 1}\n\ p-norm of @var{a}, @code{(sum (abs (@var{a}) .^ @var{p})) ^ (1/@var{p})}.\n\ +\n\ +@item other @var{p} @code{@var{p} < 1}\n\ +the p-pseudonorm defined as above.\n\ @end table\n\ +\n\ +If @code{\"rows\"} is given as @var{opt}, the norms of all rows of the matrix @var{a} are\n\ +returned as a column vector. Similarly, if @code{\"columns\"} or @code{\"cols\"} is passed\n\ +column norms are computed.\n\ @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); @@ -4634,86 +4646,48 @@ } 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<float> x (x_arg.float_array_value ()); - - if (! error_state) - retval(0) = x.norm (static_cast<float>(p_val)); - else - error ("norm: expecting real vector"); - } - else - { - MArray<FloatComplex> x (x_arg.float_complex_array_value ()); - - if (! error_state) - retval(0) = x.norm (static_cast<float>(p_val)); - else - error ("norm: expecting complex vector"); - } - } - } - else - { - if (! error_state) - { - if (x_arg.is_real_type ()) - { - MArray<double> x (x_arg.array_value ()); - - if (! error_state) - retval(0) = x.norm (p_val); - else - error ("norm: expecting real vector"); - } - else - { - MArray<Complex> 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"); @@ -4721,17 +4695,6 @@ 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; }
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/xnorm.cc Fri Oct 31 08:05:32 2008 +0100 @@ -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 +<http://www.gnu.org/licenses/>. + +*/ + +// author: Jaroslav Hajek <highegg@gmail.com> + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <cassert> +#include <cfloat> +#include <cmath> + +#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; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/xnorm.h Fri Oct 31 08:05:32 2008 +0100 @@ -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 +<http://www.gnu.org/licenses/>. + +*/ + +// author: Jaroslav Hajek <highegg@gmail.com> + +#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