Mercurial > octave-nkf
view src/DLD-FUNCTIONS/expm.cc @ 7910:b2f212b51488
DLD-FUNCTIONS/expm.cc (expm): Avoid GCC warning
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Wed, 09 Jul 2008 12:26:37 -0400 |
parents | 87865ed7405f |
children |
line wrap: on
line source
/* Copyright (C) 1996, 1997, 1999, 2000, 2002, 2005, 2006, 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/>. */ // Author: A. S. Hodel <scotte@eng.auburn.edu> #ifdef HAVE_CONFIG_H #include <config.h> #endif #include "defun-dld.h" #include "error.h" #include "gripes.h" #include "oct-obj.h" #include "utils.h" DEFUN_DLD (expm, args, , "-*- texinfo -*-\n\ @deftypefn {Loadable Function} {} expm (@var{a})\n\ Return the exponential of a matrix, defined as the infinite Taylor\n\ series\n\ @iftex\n\ @tex\n\ $$\n\ \\exp (A) = I + A + {A^2 \\over 2!} + {A^3 \\over 3!} + \\cdots\n\ $$\n\ @end tex\n\ @end iftex\n\ @ifinfo\n\ \n\ @example\n\ expm(a) = I + a + a^2/2! + a^3/3! + ...\n\ @end example\n\ \n\ @end ifinfo\n\ The Taylor series is @emph{not} the way to compute the matrix\n\ exponential; see Moler and Van Loan, @cite{Nineteen Dubious Ways to\n\ Compute the Exponential of a Matrix}, SIAM Review, 1978. This routine\n\ uses Ward's diagonal\n\ @iftex\n\ @tex\n\ Pad\\'e\n\ @end tex\n\ @end iftex\n\ @ifinfo\n\ Pade'\n\ @end ifinfo\n\ approximation method with three step preconditioning (SIAM Journal on\n\ Numerical Analysis, 1977). Diagonal\n\ @iftex\n\ @tex\n\ Pad\\'e\n\ @end tex\n\ @end iftex\n\ @ifinfo\n\ Pade'\n\ @end ifinfo\n\ approximations are rational polynomials of matrices\n\ @iftex\n\ @tex\n\ $D_q(a)^{-1}N_q(a)$\n\ @end tex\n\ @end iftex\n\ @ifinfo\n\ \n\ @example\n\ -1\n\ D (a) N (a)\n\ @end example\n\ \n\ @end ifinfo\n\ whose Taylor series matches the first\n\ @iftex\n\ @tex\n\ $2 q + 1 $\n\ @end tex\n\ @end iftex\n\ @ifinfo\n\ @code{2q+1}\n\ @end ifinfo\n\ terms of the Taylor series above; direct evaluation of the Taylor series\n\ (with the same preconditioning steps) may be desirable in lieu of the\n\ @iftex\n\ @tex\n\ Pad\\'e\n\ @end tex\n\ @end iftex\n\ @ifinfo\n\ Pade'\n\ @end ifinfo\n\ approximation when\n\ @iftex\n\ @tex\n\ $D_q(a)$\n\ @end tex\n\ @end iftex\n\ @ifinfo\n\ @code{Dq(a)}\n\ @end ifinfo\n\ is ill-conditioned.\n\ @end deftypefn") { octave_value retval; int nargin = args.length (); if (nargin != 1) { print_usage (); return retval; } octave_value arg = args(0); octave_idx_type nr = arg.rows (); octave_idx_type nc = arg.columns (); bool isfloat = arg.is_single_type (); int arg_is_empty = empty_arg ("expm", nr, nc); if (arg_is_empty < 0) return retval; if (arg_is_empty > 0) return isfloat ? octave_value (FloatMatrix ()) : octave_value (Matrix ()); if (nr != nc) { gripe_square_matrix_required ("expm"); return retval; } if (isfloat) { if (arg.is_real_type ()) { FloatMatrix m = arg.float_matrix_value (); if (error_state) return retval; else retval = m.expm (); } else if (arg.is_complex_type ()) { FloatComplexMatrix m = arg.float_complex_matrix_value (); if (error_state) return retval; else retval = m.expm (); } } else { if (arg.is_real_type ()) { Matrix m = arg.matrix_value (); if (error_state) return retval; else retval = m.expm (); } else if (arg.is_complex_type ()) { ComplexMatrix m = arg.complex_matrix_value (); if (error_state) return retval; else retval = m.expm (); } else { gripe_wrong_type_arg ("expm", arg); } } return retval; } /* %!assert(expm ([-49, 24; -64, 31]), [-0.735758758144742, 0.551819099658089; %! -1.471517599088239, 1.103638240715556], 128*eps); %!assert(expm ([1, 1; 0, 1]), [2.718281828459045, 2.718281828459045; %! 0.000000000000000, 2.718281828459045],4 * eps); %!test %! arg = diag ([6, 6, 6], 1); %! result = [1, 6, 18, 36; %! 0, 1, 6, 18; %! 0, 0, 1, 6; %! 0, 0, 0, 1]; %! assert(expm (arg), result); %!assert(expm (single([-49, 24; -64, 31])), single([-0.735758758144742, ... %! 0.551819099658089; -1.471517599088239, 1.103638240715556]), ... %! 512*eps('single')); %!assert(expm (single([1, 1; 0, 1])), single([2.718281828459045, ... %! 2.718281828459045; 0.000000000000000, 2.718281828459045]), ... %! 4 * eps('single')); %!test %! arg = single(diag ([6, 6, 6], 1)); %! result = single([1, 6, 18, 36; %! 0, 1, 6, 18; %! 0, 0, 1, 6; %! 0, 0, 0, 1]); %! assert(expm (arg), result); %!error <Invalid call to expm.*> expm(); %!error <Invalid call to expm.*> expm(1,2); */ /* ;;; Local Variables: *** ;;; mode: C++ *** ;;; End: *** */