# 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