Mercurial > hg > octave-lyh

changeset 6508:184ab67c3bc1

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-04-07 00:43:09 by jwe]
author jwe
date Sat, 07 Apr 2007 00:43:10 +0000
parents 679c2f987943
children 84f2d0253aea
files libcruft/ChangeLog libcruft/blas-xtra/Makefile.in libcruft/blas-xtra/xdnrm2.f libcruft/blas-xtra/xdznrm2.f liboctave/ChangeLog liboctave/MArray-C.cc liboctave/MArray-d.cc liboctave/MArray-defs.h liboctave/MArray.cc liboctave/MArray.h scripts/ChangeLog scripts/linear-algebra/norm.m src/ChangeLog src/data.cc
diffstat 14 files changed, 266 insertions(+), 30 deletions(-) [+]
line wrap: on
line diff
--- a/libcruft/ChangeLog
+++ b/libcruft/ChangeLog
@@ -1,5 +1,8 @@
 2007-04-06  John W. Eaton  <jwe@octave.org>
 
+	* blas-xtra/xdnrm2.f, blas-xtra/xdznrm2.f: New functions.
+	* blas-xtra/Makefile.in (FSRC): Add them to the list.
+
 	* ranlib/phrtsd.f (phrtsd): Ensure that the types of the arguments
 	passed to mod are the same even on 64-bit systems.
 
--- a/libcruft/blas-xtra/Makefile.in
+++ b/libcruft/blas-xtra/Makefile.in
@@ -14,7 +14,7 @@
 
 EXTERNAL_DISTFILES = $(DISTFILES)
 
-FSRC = xddot.f xerbla.f xzdotu.f
+FSRC = xddot.f xdnrm2.f xdznrm2.f xerbla.f xzdotu.f
 
 include $(TOPDIR)/Makeconf
 
new file mode 100644
--- /dev/null
+++ b/libcruft/blas-xtra/xdnrm2.f
@@ -0,0 +1,6 @@
+      subroutine xdnrm2 (n, x, incx, retval)
+      double precision dnrm2, x(*), retval
+      integer n, incx
+      retval = dnrm2 (n, x, incx);
+      return
+      end
new file mode 100644
--- /dev/null
+++ b/libcruft/blas-xtra/xdznrm2.f
@@ -0,0 +1,7 @@
+      subroutine xdznrm2 (n, x, incx, retval)
+      double precision dznrm2, retval
+      double complex x(*)
+      integer n, incx
+      retval = dznrm2 (n, x, incx);
+      return
+      end
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,11 @@
+2007-04-06  John W. Eaton  <jwe@octave.org>
+
+	* MArray-defs.h (MARRAY_NORM_BODY): New macro.
+	* MArray.h (MArray<T>::norm): New function.
+	* MArray.cc: Provide decl.
+	* MArray-d.cc (MArray<double>::norm): Define double specialization.
+	* MArray-C.cc (MArray<Complex>::norm): Define Complex specialization.
+
 2007-04-04  John W. Eaton  <jwe@octave.org>
 
 	* Range.cc (Range::nelem_internal): Likewise.
--- a/liboctave/MArray-C.cc
+++ b/liboctave/MArray-C.cc
@@ -28,10 +28,25 @@
 // 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"
 
+template <>
+double
+MArray<Complex>::norm (double p) const
+{
+  MARRAY_NORM_BODY (Complex, xdznrm2, XDZNRM2);
+}
+
 template class OCTAVE_API MArray<Complex>;
 
 INSTANTIATE_MARRAY_FRIENDS (Complex)
--- a/liboctave/MArray-d.cc
+++ b/liboctave/MArray-d.cc
@@ -27,9 +27,25 @@
 
 // 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"
 
+template <>
+double
+MArray<double>::norm (double p) const
+{
+  MARRAY_NORM_BODY (double, xdnrm2, XDNRM2);
+}
+
 template class MArray<double>;
 
 INSTANTIATE_MARRAY_FRIENDS (double)
--- a/liboctave/MArray-defs.h
+++ b/liboctave/MArray-defs.h
@@ -321,6 +321,78 @@
   MDIAGARRAY2_DADA_BINOP_FWD_DEFS \
     (R, T, dynamic_cast<const B<T>&>, R, dynamic_cast<const B<T>&>, R)
 
+#define MARRAY_NORM_BODY(TYPE, blas_norm, BLAS_NORM)	\
+ \
+  double retval = octave_NaN; \
+ \
+  octave_idx_type len = length (); \
+ \
+  if (len > 0) \
+    { \
+      const TYPE *d = data (); \
+ \
+      if (p == -1) \
+	{ \
+	  /* Frobenius norm.  */ \
+	  retval = 0; \
+ \
+	  for (octave_idx_type i = 0; i < len; i++) \
+	    { \
+	      double d_abs = std::abs (d[i]); \
+	      retval += d_abs * d_abs; \
+	    } \
+ \
+	  retval = ::sqrt (retval); \
+	} \
+      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) \
+		{ \
+		  double d_abs = std::abs (d[i++]); \
+ \
+		  if (d_abs > retval) \
+		    retval = d_abs; \
+		} \
+	    } \
+	  else \
+	    { \
+	      while (i < len) \
+		{ \
+		  double 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++) \
+	    { \
+	      double 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.cc
+++ b/liboctave/MArray.cc
@@ -33,6 +33,16 @@
 
 // One dimensional array with math ops.
 
+template <class T>
+double
+MArray<T>::norm (double) const
+{
+  (*current_liboctave_error_handler)
+    ("norm: only implemented for double and complex values");
+
+  return 0;
+}
+
 // Element by element MArray by scalar ops.
 
 template <class T>
--- a/liboctave/MArray.h
+++ b/liboctave/MArray.h
@@ -80,6 +80,7 @@
       return retval;
     }
 
+  double norm (double p) const;
 
   // Currently, the OPS functions don't need to be friends, but that
   // may change.
--- a/scripts/ChangeLog
+++ b/scripts/ChangeLog
@@ -1,3 +1,11 @@
+2007-04-06  John W. Eaton  <jwe@octave.org>
+
+	* linear-algebra/norm.m: Use new __vnorm__ function for vector args.
+
+2007-04-06  Daniel J Sebald  <daniel.sebald@ieee.org>
+
+	* plot/stem.m: Use plot instead of a series of calls to line.
+
 2007-04-05  John W. Eaton  <jwe@octave.org>
 
 	* sparse/nonzeros.m, sparse/normest.m, sparse/spconvert.m,
--- a/scripts/linear-algebra/norm.m
+++ b/scripts/linear-algebra/norm.m
@@ -69,14 +69,17 @@
   endif
 
   if (ndims (x) > 2)
-    error ("norm: Only valid on 2-D objects")
+    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 (rows (x) == 1 || columns (x) == 1)
-
-    if (nargin == 2)
+  if (is_vector (x))
+    if (isinteger (x) || issparse (x))
       if (ischar (p))
         if (strcmp (p, "fro"))
 	  retval = sqrt (sum (abs (x) .^ 2));
@@ -94,36 +97,48 @@
           retval = sum (abs (x) .^ p) ^ (1/p);
         endif
       endif
-    elseif (nargin == 1)
-      retval = sqrt (sum (abs (x) .^ 2));
+    else
+      retval = __vnorm__ (x, p);
     endif
-
   else
-
-    if (nargin == 2)
-      if (ischar (p))
-        if (strcmp (p, "fro"))
-	  retval = sqrt (sum (sum (abs (x) .^ 2)));
-        elseif (strcmp (p, "inf"))
-          retval = max (sum (abs (x')));
-        else
-          error ("norm: unrecognized norm");
-        endif
+    if (ischar (p))
+      if (strcmp (p, "fro"))
+	retval = sqrt (sum (sum (abs (x) .^ 2)));
+      elseif (strcmp (p, "inf"))
+        retval = max (sum (abs (x')));
       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')));
-        endif
+        error ("norm: unrecognized vector norm");
       endif
-    elseif (nargin == 1)
-      s = svd (x);
-      retval = s (1);
+    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')));
+      else
+	error ("norm: unrecognized matrix norm");
+      endif
     endif
-
   endif
 
 endfunction
+
+%!shared x
+%! x = [1, -3, 4, 5, -7];
+%!assert(norm(x,1), 20);
+%!assert(norm(x,2), 10);
+%!assert(norm(x,3), 8.24257059961711, -4*eps);
+%!assert(norm(x,Inf), 7);
+%!assert(norm(x,-Inf), 1);
+%!assert(norm(x,"inf"), 7);
+%!assert(norm(x,"fro"), 10);
+%!assert(norm(x), 10);
+%!assert(norm([1e200, 1]), 1e200);
+%!shared m
+%! m = magic (4);
+%!assert(norm(m,1), 34);
+%!assert(norm(m,2), 34);
+%!assert(norm(m,Inf), 34);
+%!assert(norm(m,"inf"), 34);
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,5 +1,7 @@
 2007-04-06  John W. Eaton  <jwe@octave.org>
 
+	* data.cc (F__vnorm__): New function.
+
 	* pt-fcn-handle.cc (tree_anon_fcn_handle::param_list, 
 	tree_anon_fcn_handle::cmd_list,	tree_anon_fcn_handle::ret_list,
 	tree_anon_fcn_handle::sym_tab): Delete.  Remove all uses.
--- a/src/data.cc
+++ b/src/data.cc
@@ -2583,6 +2583,79 @@
   return retval;
 }
 
+DEFUN (__vnorm__, args, ,
+  "-*- texinfo -*-\n\
+@deftypefn {Built-in Function} {} __vnorm__ (@var{x}, @var{p})\n\
+Compute various norms of the vector @var{x}.\n\
+@end deftypefn")
+{
+  octave_value retval;
+
+  int nargin = args.length ();
+
+  if (nargin == 1 || nargin == 2)
+    {
+      double p_val;
+
+      octave_value p_arg;
+
+      if (nargin == 1)
+	p_arg = 2;
+      else
+	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 ());
+	      return retval;
+	    }
+	}
+      else
+	{
+	  p_val = args(1).double_value ();
+
+	  if (error_state)
+	    {
+	      error ("norm: unrecognized norm value");
+	      return retval;
+	    }
+	}
+
+      octave_value x_arg = args(0);
+
+      if (x_arg.is_real_type ())
+	{
+	  ColumnVector x (x_arg.vector_value ());
+
+	  if (! error_state)
+	    retval = x.norm (p_val);
+	  else
+	    error ("norm: expecting real vector");
+	}
+      else
+	{
+	  ComplexColumnVector x (x_arg.complex_vector_value ());
+
+	  if (! error_state)
+	    retval = x.norm (p_val);
+	  else
+	    error ("norm: expecting complex vector");
+	}
+    }
+  else
+    print_usage ();
+
+  return retval;
+}
+
 /*
 ;;; Local Variables: ***
 ;;; mode: C++ ***