# HG changeset patch # User jwe # Date 1045613630 0 # Node ID d3278845f764a90bc1f19a5ace521ec67ac2937e # Parent 9f86c2055b5826123da7db835d6d8c858b2417b7 [project @ 2003-02-19 00:12:31 by jwe] diff --git a/src/ChangeLog b/src/ChangeLog --- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,3 +1,14 @@ +2003-02-18 John W. Eaton + + * 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 + + * DLD-FUNCTIONS/sqrtm.cc: New file. + 2003-02-18 David Bateman * DLD-FUNCTIONS/lu.cc (Flu): Allow non-square matrices. diff --git a/src/DLD-FUNCTIONS/log.cc b/src/DLD-FUNCTIONS/log.cc deleted file mode 100644 --- a/src/DLD-FUNCTIONS/log.cc +++ /dev/null @@ -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 -#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: *** -*/ diff --git a/src/DLD-FUNCTIONS/sqrtm.cc b/src/DLD-FUNCTIONS/sqrtm.cc new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/sqrtm.cc @@ -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 +#endif + +#include + +#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; +} diff --git a/src/Makefile.in b/src/Makefile.in --- a/src/Makefile.in +++ b/src/Makefile.in @@ -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))