Mercurial > octave

comparison src/expm.cc @ 50:6028dcac27ef

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1993-08-11 20:48:00 by jwe] Initial revision
author jwe
date Wed, 11 Aug 1993 20:48:00 +0000
parents
children 9eda1009cf95
comparison
equal deleted inserted replaced
49:445ea777560a 50:6028dcac27ef
1 // tc-expm.cc -*- C++ -*-
2 /*
3
4 Copyright (C) 1993 John W. Eaton
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
10 Free Software Foundation; either version 2, or (at your option) any
11 later version.
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
19 along with Octave; see the file COPYING. If not, write to the Free
20 Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
21
22 */
23
24 // Written by A. S. Hodel <scotte@eng.auburn.edu>
25
26 #ifdef __GNUG__
27 #pragma implementation
28 #endif
29
30 #include <math.h>
31
32 #include "Matrix.h"
33
34 #include "tree-const.h"
35 #include "user-prefs.h"
36 #include "gripes.h"
37 #include "error.h"
38 #include "f-expm.h"
39
40 #ifdef WITH_DLD
41 tree_constant *
42 builtin_matrix_exp_2 (tree_constant *args, int nargin, int nargout)
43 {
44 tree_constant *retval = new tree_constant [2];
45 retval[0] = matrix_exp (args[1]);
46 return retval;
47 }
48 #endif
49
50 extern "C"
51 {
52 double F77_FCN (dlange) (const char*, const int*, const int*,
53 const double*, const int*, double*);
54
55 double F77_FCN (zlange) (const char*, const int*, const int*,
56 const Complex*, const int*, double*);
57 }
58
59 // XXX FIXME XXX
60 extern int empty_arg (tree_constant&);
61
62 tree_constant
63 matrix_exp (const tree_constant& a)
64 {
65 tree_constant retval;
66 tree_constant tmp = a.make_numeric ();
67
68 // Constants for matrix exponential calculation.
69
70 static double padec[] =
71 {
72 5.0000000000000000e-1,
73 1.1666666666666667e-1,
74 1.6666666666666667e-2,
75 1.6025641025641026e-3,
76 1.0683760683760684e-4,
77 4.8562548562548563e-6,
78 1.3875013875013875e-7,
79 1.9270852604185938e-9,
80 };
81
82 if (empty_arg (tmp))
83 {
84 int flag = user_pref.propagate_empty_matrices;
85 if (flag != 0)
86 {
87 if (flag < 0)
88 gripe_empty_arg ("expm", 0);
89 Matrix m;
90 retval = tree_constant (m);
91 }
92 else gripe_empty_arg ("expm", 1);
93 }
94 else if (tmp.rows () != tmp.columns ())
95 gripe_square_matrix_required ("expm");
96 else
97 {
98 int i, j;
99 int n_cols = tmp.columns ();
100
101 char* balance_job = "B"; // variables for balancing
102
103 int sqpow; // power for scaling and squaring
104 double inf_norm; // norm of preconditioned matrix
105 int minus_one_j; // used in computing pade approx
106
107 switch (tmp.const_type ())
108 {
109 case tree_constant_rep::complex_matrix_constant:
110 {
111 ComplexMatrix m = tmp.complex_matrix_value ();
112 Complex trshift = 0.0; // trace shift value
113
114 // Preconditioning step 1: trace normalization.
115
116 for (i = 0; i < n_cols; i++)
117 trshift += m.elem (i, i);
118 trshift /= n_cols;
119 for (i = 0; i < n_cols; i++)
120 m.elem (i, i) -= trshift;
121
122 // Preconditioning step 2: eigenvalue balancing.
123
124 ComplexAEPBALANCE mbal (m, balance_job);
125 m = mbal.balanced_matrix ();
126 ComplexMatrix d = mbal.balancing_matrix ();
127
128 // Preconditioning step 3: scaling.
129
130 ColumnVector work (n_cols);
131 inf_norm = F77_FCN (zlange) ("I", &n_cols, &n_cols, m.
132 fortran_vec (), &n_cols,
133 work.fortran_vec ());
134
135 sqpow = 1.0 + log (inf_norm) / log (2.0);
136
137 // Check whether we need to square at all.
138
139 if (sqpow <= 0.0)
140 sqpow = 0.0;
141 else
142 {
143 for (inf_norm = 1.0, i = 0; i < sqpow; i++)
144 inf_norm *= 2.0;
145
146 m = m / inf_norm;
147 }
148
149 // npp, dpp: pade' approx polynomial matrices.
150
151 ComplexMatrix npp (n_cols, n_cols, 0.0);
152 ComplexMatrix dpp = npp;
153
154 // Now powers a^8 ... a^1.
155
156 minus_one_j = -1;
157 for (j = 7; j >= 0; j--)
158 {
159 npp = m * npp + m * padec[j];
160 dpp = m * dpp + m * (minus_one_j * padec[j]);
161 minus_one_j *= -1;
162 }
163
164 // Zero power.
165
166 dpp = -dpp;
167 for (j = 0; j < n_cols; j++)
168 {
169 npp.elem (j, j) += 1.0;
170 dpp.elem (j, j) += 1.0;
171 }
172
173 // Compute pade approximation = inverse (dpp) * npp.
174
175 ComplexMatrix result = dpp.solve (npp);
176
177 // Reverse preconditioning step 3: repeated squaring.
178
179 while (sqpow)
180 {
181 result = result * result;
182 sqpow--;
183 }
184
185 // reverse preconditioning step 2: inverse balancing XXX FIXME XXX:
186 // should probably do this with lapack calls instead of a complete
187 // matrix inversion.
188
189 result = result.transpose ();
190 d = d.transpose ();
191 result = result * d;
192 result = d.solve (result);
193 result = result.transpose ();
194
195 // Reverse preconditioning step 1: fix trace normalization.
196
197 result = result * exp (trshift);
198
199 retval = tree_constant (result);
200 }
201 break;
202 case tree_constant_rep::complex_scalar_constant:
203 {
204 Complex c = tmp.complex_value ();
205 retval = tree_constant (exp (c));
206 }
207 break;
208 case tree_constant_rep::matrix_constant:
209 {
210
211 // Compute the exponential.
212
213 Matrix m = tmp.matrix_value ();
214
215 double trshift = 0; // trace shift value
216
217 // Preconditioning step 1: trace normalization.
218
219 for (i = 0; i < n_cols; i++)
220 trshift += m.elem (i, i);
221 trshift /= n_cols;
222 for (i = 0; i < n_cols; i++)
223 m.elem (i, i) -= trshift;
224
225 // Preconditioning step 2: balancing.
226
227 AEPBALANCE mbal (m, balance_job);
228 m = mbal.balanced_matrix ();
229 Matrix d = mbal.balancing_matrix ();
230
231 // Preconditioning step 3: scaling.
232
233 ColumnVector work(n_cols);
234 inf_norm = F77_FCN (dlange) ("I", &n_cols, &n_cols,
235 m.fortran_vec (), &n_cols,
236 work.fortran_vec ());
237
238 sqpow = 1.0 + log (inf_norm) / log (2.0);
239
240 // Check whether we need to square at all.
241
242 if (sqpow <= 0.0)
243 sqpow = 0.0;
244 else
245 {
246 for (inf_norm = 1.0, i = 0; i < sqpow; i++)
247 inf_norm *= 2.0;
248
249 m = m / inf_norm;
250 }
251
252 // npp, dpp: pade' approx polynomial matrices.
253
254 Matrix npp (n_cols, n_cols, 0.0);
255 Matrix dpp = npp;
256
257 // now powers a^8 ... a^1.
258
259 minus_one_j = -1;
260 for (j = 7; j >= 0; j--)
261 {
262 npp = m * npp + m * padec[j];
263 dpp = m * dpp + m * (minus_one_j * padec[j]);
264 minus_one_j *= -1;
265 }
266 // Zero power.
267
268 dpp = -dpp;
269 for(j = 0; j < n_cols; j++)
270 {
271 npp.elem (j, j) += 1.0;
272 dpp.elem (j, j) += 1.0;
273 }
274
275 // Compute pade approximation = inverse (dpp) * npp.
276
277 Matrix result = dpp.solve (npp);
278
279 // Reverse preconditioning step 3: repeated squaring.
280
281 while(sqpow)
282 {
283 result = result * result;
284 sqpow--;
285 }
286
287 // Reverse preconditioning step 2: inverse balancing.
288
289 result = result.transpose();
290 d = d.transpose ();
291 result = result * d;
292 result = d.solve (result);
293 result = result.transpose ();
294
295 // Reverse preconditioning step 1: fix trace normalization.
296
297 result = result * exp (trshift);
298
299 retval = tree_constant (result);
300 }
301 break;
302 case tree_constant_rep::scalar_constant:
303 {
304 double d = tmp.double_value ();
305 retval = tree_constant (exp (d));
306 }
307 break;
308 default:
309 panic_impossible();
310 break;
311 }
312 }
313 return retval;
314 }
315
316 /*
317 ;;; Local Variables: ***
318 ;;; mode: C++ ***
319 ;;; page-delimiter: "^/\\*" ***
320 ;;; End: ***
321 */