Mercurial > octave-nkf

changeset 4331:d3278845f764

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2003-02-19 00:12:31 by jwe]
author jwe
date Wed, 19 Feb 2003 00:13:50 +0000
parents 9f86c2055b58
children e41906608e0f
files src/ChangeLog src/DLD-FUNCTIONS/log.cc src/DLD-FUNCTIONS/sqrtm.cc src/Makefile.in
diffstat 4 files changed, 285 insertions(+), 284 deletions(-) [+]
line wrap: on
line diff
--- a/src/ChangeLog	Tue Feb 18 22:31:55 2003 +0000
+++ b/src/ChangeLog	Wed Feb 19 00:13:50 2003 +0000
@@ -1,3 +1,14 @@
+2003-02-18  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* Makefile.in (DLD_XSRC): Delete log.cc from the list.
+	Add sqrtm.cc to the list.
+
+	* DLD-FUNCTIONS/log.cc: Delete.
+
+2003-02-18  Paul Kienzle <pkienzle@users.sf.net>
+
+	* DLD-FUNCTIONS/sqrtm.cc: New file.
+
 2003-02-18  David Bateman <dbateman@free.fr>
 
 	* DLD-FUNCTIONS/lu.cc (Flu): Allow non-square matrices.
--- a/src/DLD-FUNCTIONS/log.cc	Tue Feb 18 22:31:55 2003 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,282 +0,0 @@
-/*
-
-Copyright (C) 1996, 1997 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 2, 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, write to the Free
-Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
-
-*/
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include "EIG.h"
-#include "mx-cm-cdm.h"
-#include "quit.h"
-
-#include "defun-dld.h"
-#include "error.h"
-#include "gripes.h"
-#include "oct-obj.h"
-#include "utils.h"
-
-// XXX FIXME XXX -- the next two functions should really be just
-// one...
-
-DEFUN_DLD (logm, args, ,
-  "-*- texinfo -*-\n\
-@deftypefn {Loadable Function} {} logm (@var{a})\n\
-Compute the matrix logarithm of the square matrix @var{a}.  Note that\n\
-this is currently implemented in terms of an eigenvalue expansion and\n\
-needs to be improved to be more robust.\n\
-@end deftypefn")
-{
-  octave_value_list retval;
-
-  int nargin = args.length ();
-
-  if (nargin != 1)
-    {
-      print_usage ("logm");
-      return retval;
-    }
-
-  octave_value arg = args(0);
-
-  int arg_is_empty = empty_arg ("logm", arg.rows (), arg.columns ());
-
-  if (arg_is_empty < 0)
-    return retval;
-  else if (arg_is_empty > 0)
-    return octave_value (Matrix ());
-
-  if (arg.is_real_scalar ())
-    {
-      double d = arg.double_value ();
-      if (d > 0.0)
-	retval(0) = log (d);
-      else
-	{
-	  Complex dtmp (d);
-	  retval(0) = log (dtmp);
-	}
-    }
-  else if (arg.is_complex_scalar ())
-    {
-      Complex c = arg.complex_value ();
-      retval(0) = log (c);
-    }
-  else if (arg.is_real_type ())
-    {
-      Matrix m = arg.matrix_value ();
-
-      if (! error_state)
-	{
-	  int nr = m.rows ();
-	  int nc = m.columns ();
-
-	  if (nr == 0 || nc == 0 || nr != nc)
-	    gripe_square_matrix_required ("logm");
-	  else
-	    {
-	      EIG m_eig (m);
-	      ComplexColumnVector lambda (m_eig.eigenvalues ());
-	      ComplexMatrix Q (m_eig.eigenvectors ());
-
-	      for (int i = 0; i < nr; i++)
-		{
-		  OCTAVE_QUIT;
-		  Complex elt = lambda (i);
-		  if (imag (elt) == 0.0 && real (elt) > 0.0)
-		    lambda (i) = log (real (elt));
-		  else
-		    lambda (i) = log (elt);
-		}
-
-	      ComplexDiagMatrix D (lambda);
-	      ComplexMatrix result = Q * D * Q.inverse ();
-
-	      retval(0) = result;
-	    }
-	}
-    }
-  else if (arg.is_complex_type ())
-    {
-      ComplexMatrix m = arg.complex_matrix_value ();
-
-      if (! error_state)
-	{
-	  int nr = m.rows ();
-	  int nc = m.columns ();
-
-	  if (nr == 0 || nc == 0 || nr != nc)
-	    gripe_square_matrix_required ("logm");
-	  else
-	    {
-	      EIG m_eig (m);
-	      ComplexColumnVector lambda (m_eig.eigenvalues ());
-	      ComplexMatrix Q (m_eig.eigenvectors ());
-
-	      for (int i = 0; i < nr; i++)
-		{
-		  OCTAVE_QUIT;
-		  Complex elt = lambda (i);
-		  if (imag (elt) == 0.0 && real (elt) > 0.0)
-		    lambda (i) = log (real (elt));
-		  else
-		    lambda (i) = log (elt);
-		}
-
-	      ComplexDiagMatrix D (lambda);
-	      ComplexMatrix result = Q * D * Q.inverse ();
-
-	      retval(0) = result;
-	    }
-	}
-    }
-  else
-    {
-      gripe_wrong_type_arg ("logm", arg);
-    }
-
-  return retval;
-}
-
-DEFUN_DLD (sqrtm, args, ,
- "-*- texinfo -*-\n\
-@deftypefn {Loadable Function} {} sqrtm (@var{a})\n\
-Compute the matrix square root of the square matrix @var{a}.  Note that\n\
-this is currently implemented in terms of an eigenvalue expansion and\n\
-needs to be improved to be more robust.\n\
-@end deftypefn")
-{
-  octave_value_list retval;
-
-  int nargin = args.length ();
-
-  if (nargin != 1)
-    {
-      print_usage ("sqrtm");
-      return retval;
-    }
-
-  octave_value arg = args(0);
-
-  int arg_is_empty = empty_arg ("sqrtm", arg.rows (), arg.columns ());
-
-  if (arg_is_empty < 0)
-    return retval;
-  else if (arg_is_empty > 0)
-    return octave_value (Matrix ());
-
-  if (arg.is_real_scalar ())
-    {
-      double d = arg.double_value ();
-      if (d > 0.0)
-	retval(0) = sqrt (d);
-      else
-	{
-	  Complex dtmp (d);
-	  retval(0) = sqrt (dtmp);
-	}
-    }
-  else if (arg.is_complex_scalar ())
-    {
-      Complex c = arg.complex_value ();
-      retval(0) = sqrt (c);
-    }
-  else if (arg.is_real_type ())
-    {
-      Matrix m = arg.matrix_value ();
-
-      if (! error_state)
-	{
-	  int nr = m.rows ();
-	  int nc = m.columns ();
-
-	  if (nr == 0 || nc == 0 || nr != nc)
-	    gripe_square_matrix_required ("sqrtm");
-	  else
-	    {
-	      EIG m_eig (m);
-	      ComplexColumnVector lambda (m_eig.eigenvalues ());
-	      ComplexMatrix Q (m_eig.eigenvectors ());
-
-	      for (int i = 0; i < nr; i++)
-		{
-		  OCTAVE_QUIT;
-		  Complex elt = lambda (i);
-		  if (imag (elt) == 0.0 && real (elt) > 0.0)
-		    lambda (i) = sqrt (real (elt));
-		  else
-		    lambda (i) = sqrt (elt);
-		}
-
-	      ComplexDiagMatrix D (lambda);
-	      ComplexMatrix result = Q * D * Q.inverse ();
-
-	      retval(0) = result;
-	    }
-	}
-    }
-  else if (arg.is_complex_type ())
-    {
-      ComplexMatrix m = arg.complex_matrix_value ();
-
-      if (! error_state)
-	{
-	  int nr = m.rows ();
-	  int nc = m.columns ();
-
-	  if (nr == 0 || nc == 0 || nr != nc)
-	    gripe_square_matrix_required ("sqrtm");
-	  else
-	    {
-	      EIG m_eig (m);
-	      ComplexColumnVector lambda (m_eig.eigenvalues ());
-	      ComplexMatrix Q (m_eig.eigenvectors ());
-
-	      for (int i = 0; i < nr; i++)
-		{
-		  OCTAVE_QUIT;
-		  Complex elt = lambda (i);
-		  if (imag (elt) == 0.0 && real (elt) > 0.0)
-		    lambda (i) = sqrt (real (elt));
-		  else
-		    lambda (i) = sqrt (elt);
-		}
-
-	      ComplexDiagMatrix D (lambda);
-	      ComplexMatrix result = Q * D * Q.inverse ();
-
-	      retval(0) = result;
-	    }
-	}
-    }
-  else
-    {
-      gripe_wrong_type_arg ("sqrtm", arg);
-    }
-
-  return retval;
-}
-
-/*
-;;; Local Variables: ***
-;;; mode: C++ ***
-;;; End: ***
-*/
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/DLD-FUNCTIONS/sqrtm.cc	Wed Feb 19 00:13:50 2003 +0000
@@ -0,0 +1,272 @@
+/*
+
+Copyright (C) 2001 Ross Lippert and Paul Kienzle
+
+This program 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 2 of the License, or
+(at your option) any later version.
+
+This program 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 this program; if not, write to the Free Software
+Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <float.h>
+
+#include "CmplxSCHUR.h"
+#include "lo-ieee.h"
+#include "lo-mappers.h"
+
+#include "defun-dld.h"
+#include "error.h"
+#include "gripes.h"
+#include "utils.h"
+
+static inline double
+getmin (double x, double y)
+{
+  return x < y ? x : y;
+}
+
+static inline double
+getmax (double x, double y)
+{
+  return x > y ? x : y;
+}
+
+static double
+frobnorm (const ComplexMatrix& A)
+{
+  double sum = 0;
+
+  for (int i = 0; i < A.rows (); i++)
+    for (int j = 0; j < A.columns (); j++)
+      sum += real (A(i,j) * conj (A(i,j)));
+
+  return sqrt (sum);
+}
+
+static double
+frobnorm (const Matrix& A)
+{
+  double sum = 0;
+  for (int i = 0; i < A.rows (); i++)
+    for (int j = 0; j < A.columns (); j++)
+      sum += A(i,j) * A(i,j);
+
+  return sqrt (sum);
+}
+
+
+static ComplexMatrix
+sqrtm_from_schur (const ComplexMatrix& U, const ComplexMatrix& T)
+{
+  const int n = U.rows ();
+
+  ComplexMatrix R (n, n, 0.0);
+
+  for (int j = 0; j < n; j++)
+    R(j,j) = sqrt (T(j,j));
+
+  const double fudge = sqrt (DBL_MIN);
+
+  for (int p = 0; p < n-1; p++)
+    {
+      for (int i = 0; i < n-(p+1); i++)
+	{
+	  const int j = i + p + 1;
+
+	  Complex s = T(i,j);
+
+	  for (int k = i+1; k < j; k++)
+	    s -= R(i,k) * R(k,j);
+
+	  // dividing
+	  //     R(i,j) = s/(R(i,i)+R(j,j));
+	  // screwing around to not / 0
+
+	  const Complex d = R(i,i) + R(j,j) + fudge;
+	  const Complex conjd = conj (d);
+
+	  R(i,j) =  (s*conjd)/(d*conjd);
+	}
+    }
+
+  return U * R * U.hermitian ();
+}
+
+DEFUN_DLD (sqrtm, args, nargout,
+ "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {[@var{result}, @var{error_estimate}] =} sqrtm (@var{a})\n\
+Compute the matrix square root of the square matrix @var{a}.\n\
+\n\
+Ref: Nicholas J. Higham. A new sqrtm for MATLAB. Numerical Analysis\n\
+Report No. 336, Manchester Centre for Computational Mathematics,\n\
+Manchester, England, January 1999.\n\
+@end deftypefn\n\
+@seealso{expm, logm, and funm}")
+{
+  octave_value_list retval;
+
+  int nargin = args.length ();
+
+  if (nargin != 1)
+    {
+      print_usage ("sqrtm");
+      return retval;
+    }
+
+  octave_value arg = args(0);
+
+  int n = arg.rows ();
+  int nc = arg.columns ();
+
+  int arg_is_empty = empty_arg ("sqrtm", n, nc);
+
+  if (arg_is_empty < 0)
+    return retval;
+  else if (arg_is_empty > 0)
+    return octave_value (Matrix ());
+
+  if (n != nc)
+    {
+      gripe_square_matrix_required ("sqrtm");
+      return retval;
+    }
+
+  retval(1) = lo_ieee_inf_value ();
+  retval(0) = lo_ieee_nan_value ();
+
+  if (arg.is_real_scalar ())
+    {
+      double d = arg.double_value ();
+      if (d > 0.0)
+	{
+	  retval(0) = sqrt (d);
+	  retval(1) = 0.0;
+	}
+      else
+	{
+	  retval(0) = Complex (0.0, sqrt (d));
+	  retval(1) = 0.0;
+	}
+    }
+  else if (arg.is_complex_scalar ())
+    {
+      Complex c = arg.complex_value ();
+      retval(0) = sqrt (c);
+      retval(1) = 0.0;
+    }
+  else if (arg.is_matrix_type ())
+    {
+      double err, minT;
+
+      if (arg.is_real_matrix ())
+	{
+	  Matrix A = arg.matrix_value();
+
+	  if (error_state)
+	    return retval;
+
+	  // XXX FIXME XXX -- eventually, ComplexSCHUR will accept a
+	  // real matrix arg.
+
+	  ComplexMatrix Ac (A);
+
+	  const ComplexSCHUR schur (Ac, std::string ());
+
+	  if (error_state)
+	    return retval;
+
+	  const ComplexMatrix U (schur.unitary_matrix ());
+	  const ComplexMatrix T (schur.schur_matrix ());
+	  const ComplexMatrix X (sqrtm_from_schur (U, T));
+
+	  // Check for minimal imaginary part
+	  double normX = 0.0;
+	  double imagX = 0.0;
+	  for (int i = 0; i < n; i++)
+	    for (int j = 0; j < n; j++)
+	      {
+		imagX = getmax (imagX, imag (X(i,j)));
+		normX = getmax (normX, abs (X(i,j)));
+	      }
+
+	  if (imagX < normX * 100 * DBL_EPSILON)
+	    retval(0) = real (X);
+	  else
+	    retval(0) = X;
+
+	  // Compute error
+	  // XXX FIXME XXX can we estimate the error without doing the
+	  // matrix multiply?
+
+	  err = frobnorm (X*X - ComplexMatrix (A)) / frobnorm (A);
+
+	  if (xisnan (err))
+	    err = lo_ieee_inf_value ();
+
+	  // Find min diagonal
+	  minT = lo_ieee_inf_value ();
+	  for (int i=0; i < n; i++)
+	    minT = getmin(minT, abs(T(i,i)));
+	}
+      else
+	{
+	  ComplexMatrix A = arg.complex_matrix_value ();
+
+	  if (error_state)
+	    return retval;
+
+	  const ComplexSCHUR schur (A, std::string ());
+
+	  if (error_state)
+	    return retval;
+
+	  const ComplexMatrix U (schur.unitary_matrix ());
+	  const ComplexMatrix T (schur.schur_matrix ());
+	  const ComplexMatrix X (sqrtm_from_schur (U, T));
+
+	  retval(0) = X;
+
+	  err = frobnorm (X*X - A) / frobnorm (A);
+
+	  if (xisnan (err))
+	    err = lo_ieee_inf_value ();
+
+	  minT = lo_ieee_inf_value ();
+	  for (int i = 0; i < n; i++)
+	    minT = getmin (minT, abs (T(i,i)));
+	}
+
+      retval(1) = err;
+
+      if (nargout < 2)
+	{
+	  if (err > 100*(minT+DBL_EPSILON)*n)
+	    {
+	      if (minT == 0.0)
+		error ("sqrtm: A is singular, sqrt may not exist");
+	      else if (minT <= sqrt (DBL_MIN))
+		error ("sqrtm: A is nearly singular, failed to find sqrt");
+	      else
+		error ("sqrtm: failed to find sqrt");
+	    }
+	}
+    }
+  else
+    gripe_wrong_type_arg ("sqrtm", arg);
+
+  return retval;
+}
--- a/src/Makefile.in	Tue Feb 18 22:31:55 2003 +0000
+++ b/src/Makefile.in	Wed Feb 19 00:13:50 2003 +0000
@@ -47,9 +47,9 @@
 	daspk.cc dasrt.cc dassl.cc det.cc eig.cc expm.cc fft.cc fft2.cc \
 	filter.cc find.cc fsolve.cc gammainc.cc getgrent.cc \
 	getpwent.cc getrusage.cc givens.cc hess.cc ifft.cc ifft2.cc \
-	inv.cc kron.cc log.cc lpsolve.cc lsode.cc lu.cc minmax.cc \
+	inv.cc kron.cc lpsolve.cc lsode.cc lu.cc minmax.cc \
 	odessa.cc pinv.cc qr.cc quad.cc qz.cc rand.cc schur.cc \
-	sort.cc svd.cc syl.cc time.cc
+	sort.cc sqrtm.cc svd.cc syl.cc time.cc
 
 DLD_SRC := $(addprefix DLD-FUNCTIONS/, $(DLD_XSRC))