Mercurial > octave-nkf
annotate 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 |
rev | line source |
---|---|
2928 | 1 /* |
2 | |
7017 | 3 Copyright (C) 1996, 1997, 1999, 2000, 2002, 2005, 2006, 2007 |
4 John W. Eaton | |
2928 | 5 |
6 This file is part of Octave. | |
7 | |
8 Octave is free software; you can redistribute it and/or modify it | |
9 under the terms of the GNU General Public License as published by the | |
7016 | 10 Free Software Foundation; either version 3 of the License, or (at your |
11 option) any later version. | |
2928 | 12 |
13 Octave is distributed in the hope that it will be useful, but WITHOUT | |
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
16 for more details. | |
17 | |
18 You should have received a copy of the GNU General Public License | |
7016 | 19 along with Octave; see the file COPYING. If not, see |
20 <http://www.gnu.org/licenses/>. | |
2928 | 21 |
22 */ | |
23 | |
3911 | 24 // Author: A. S. Hodel <scotte@eng.auburn.edu> |
2928 | 25 |
26 #ifdef HAVE_CONFIG_H | |
27 #include <config.h> | |
28 #endif | |
29 | |
30 #include "defun-dld.h" | |
31 #include "error.h" | |
32 #include "gripes.h" | |
33 #include "oct-obj.h" | |
34 #include "utils.h" | |
35 | |
36 DEFUN_DLD (expm, args, , | |
3548 | 37 "-*- texinfo -*-\n\ |
3372 | 38 @deftypefn {Loadable Function} {} expm (@var{a})\n\ |
39 Return the exponential of a matrix, defined as the infinite Taylor\n\ | |
40 series\n\ | |
41 @iftex\n\ | |
42 @tex\n\ | |
43 $$\n\ | |
44 \\exp (A) = I + A + {A^2 \\over 2!} + {A^3 \\over 3!} + \\cdots\n\ | |
45 $$\n\ | |
46 @end tex\n\ | |
47 @end iftex\n\ | |
48 @ifinfo\n\ | |
49 \n\ | |
50 @example\n\ | |
51 expm(a) = I + a + a^2/2! + a^3/3! + ...\n\ | |
52 @end example\n\ | |
53 \n\ | |
54 @end ifinfo\n\ | |
55 The Taylor series is @emph{not} the way to compute the matrix\n\ | |
56 exponential; see Moler and Van Loan, @cite{Nineteen Dubious Ways to\n\ | |
57 Compute the Exponential of a Matrix}, SIAM Review, 1978. This routine\n\ | |
58 uses Ward's diagonal\n\ | |
59 @iftex\n\ | |
60 @tex\n\ | |
61 Pad\\'e\n\ | |
62 @end tex\n\ | |
63 @end iftex\n\ | |
64 @ifinfo\n\ | |
65 Pade'\n\ | |
66 @end ifinfo\n\ | |
67 approximation method with three step preconditioning (SIAM Journal on\n\ | |
68 Numerical Analysis, 1977). Diagonal\n\ | |
69 @iftex\n\ | |
70 @tex\n\ | |
71 Pad\\'e\n\ | |
72 @end tex\n\ | |
73 @end iftex\n\ | |
74 @ifinfo\n\ | |
75 Pade'\n\ | |
76 @end ifinfo\n\ | |
77 approximations are rational polynomials of matrices\n\ | |
78 @iftex\n\ | |
79 @tex\n\ | |
80 $D_q(a)^{-1}N_q(a)$\n\ | |
81 @end tex\n\ | |
82 @end iftex\n\ | |
83 @ifinfo\n\ | |
84 \n\ | |
85 @example\n\ | |
86 -1\n\ | |
87 D (a) N (a)\n\ | |
88 @end example\n\ | |
89 \n\ | |
90 @end ifinfo\n\ | |
91 whose Taylor series matches the first\n\ | |
92 @iftex\n\ | |
93 @tex\n\ | |
94 $2 q + 1 $\n\ | |
95 @end tex\n\ | |
96 @end iftex\n\ | |
97 @ifinfo\n\ | |
98 @code{2q+1}\n\ | |
99 @end ifinfo\n\ | |
100 terms of the Taylor series above; direct evaluation of the Taylor series\n\ | |
101 (with the same preconditioning steps) may be desirable in lieu of the\n\ | |
102 @iftex\n\ | |
103 @tex\n\ | |
104 Pad\\'e\n\ | |
105 @end tex\n\ | |
106 @end iftex\n\ | |
107 @ifinfo\n\ | |
108 Pade'\n\ | |
109 @end ifinfo\n\ | |
110 approximation when\n\ | |
111 @iftex\n\ | |
112 @tex\n\ | |
113 $D_q(a)$\n\ | |
114 @end tex\n\ | |
115 @end iftex\n\ | |
116 @ifinfo\n\ | |
117 @code{Dq(a)}\n\ | |
118 @end ifinfo\n\ | |
119 is ill-conditioned.\n\ | |
120 @end deftypefn") | |
2928 | 121 { |
4233 | 122 octave_value retval; |
2928 | 123 |
124 int nargin = args.length (); | |
125 | |
126 if (nargin != 1) | |
127 { | |
5823 | 128 print_usage (); |
2928 | 129 return retval; |
130 } | |
131 | |
132 octave_value arg = args(0); | |
133 | |
5275 | 134 octave_idx_type nr = arg.rows (); |
135 octave_idx_type nc = arg.columns (); | |
2928 | 136 |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
137 bool isfloat = arg.is_single_type (); |
2928 | 138 int arg_is_empty = empty_arg ("expm", nr, nc); |
139 | |
140 if (arg_is_empty < 0) | |
141 return retval; | |
7910
b2f212b51488
DLD-FUNCTIONS/expm.cc (expm): Avoid GCC warning
John W. Eaton <jwe@octave.org>
parents:
7814
diff
changeset
|
142 |
2928 | 143 if (arg_is_empty > 0) |
7910
b2f212b51488
DLD-FUNCTIONS/expm.cc (expm): Avoid GCC warning
John W. Eaton <jwe@octave.org>
parents:
7814
diff
changeset
|
144 return isfloat ? octave_value (FloatMatrix ()) : octave_value (Matrix ()); |
2928 | 145 |
146 if (nr != nc) | |
147 { | |
148 gripe_square_matrix_required ("expm"); | |
149 return retval; | |
150 } | |
151 | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
152 if (isfloat) |
2928 | 153 { |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
154 if (arg.is_real_type ()) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
155 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
156 FloatMatrix m = arg.float_matrix_value (); |
2928 | 157 |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
158 if (error_state) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
159 return retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
160 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
161 retval = m.expm (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
162 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
163 else if (arg.is_complex_type ()) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
164 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
165 FloatComplexMatrix m = arg.float_complex_matrix_value (); |
2928 | 166 |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
167 if (error_state) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
168 return retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
169 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
170 retval = m.expm (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
171 } |
2928 | 172 } |
173 else | |
174 { | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
175 if (arg.is_real_type ()) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
176 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
177 Matrix m = arg.matrix_value (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
178 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
179 if (error_state) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
180 return retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
181 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
182 retval = m.expm (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
183 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
184 else if (arg.is_complex_type ()) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
185 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
186 ComplexMatrix m = arg.complex_matrix_value (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
187 |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
188 if (error_state) |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
189 return retval; |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
190 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
191 retval = m.expm (); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
192 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
193 else |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
194 { |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
195 gripe_wrong_type_arg ("expm", arg); |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
196 } |
2928 | 197 } |
198 | |
199 return retval; | |
200 } | |
201 | |
202 /* | |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
203 |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
204 %!assert(expm ([-49, 24; -64, 31]), [-0.735758758144742, 0.551819099658089; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
205 %! -1.471517599088239, 1.103638240715556], 128*eps); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
206 |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
207 %!assert(expm ([1, 1; 0, 1]), [2.718281828459045, 2.718281828459045; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
208 %! 0.000000000000000, 2.718281828459045],4 * eps); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
209 |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
210 %!test |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
211 %! arg = diag ([6, 6, 6], 1); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
212 %! result = [1, 6, 18, 36; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
213 %! 0, 1, 6, 18; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
214 %! 0, 0, 1, 6; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
215 %! 0, 0, 0, 1]; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
216 %! assert(expm (arg), result); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
217 |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
218 %!assert(expm (single([-49, 24; -64, 31])), single([-0.735758758144742, ... |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
219 %! 0.551819099658089; -1.471517599088239, 1.103638240715556]), ... |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
220 %! 512*eps('single')); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
221 |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
222 %!assert(expm (single([1, 1; 0, 1])), single([2.718281828459045, ... |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
223 %! 2.718281828459045; 0.000000000000000, 2.718281828459045]), ... |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
224 %! 4 * eps('single')); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
225 |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
226 %!test |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
227 %! arg = single(diag ([6, 6, 6], 1)); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
228 %! result = single([1, 6, 18, 36; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
229 %! 0, 1, 6, 18; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
230 %! 0, 0, 1, 6; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
231 %! 0, 0, 0, 1]); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
232 %! assert(expm (arg), result); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
233 |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
234 %!error <Invalid call to expm.*> expm(); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
235 %!error <Invalid call to expm.*> expm(1,2); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
236 |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
237 */ |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
238 |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
239 /* |
2928 | 240 ;;; Local Variables: *** |
241 ;;; mode: C++ *** | |
242 ;;; End: *** | |
243 */ |