# HG changeset patch # User jwe # Date 745102080 0 # Node ID 6028dcac27efe55282cec599a67a9eff335fd3fd # Parent 445ea777560a6704fc835710c1d0d19cc1687c6e [project @ 1993-08-11 20:48:00 by jwe] Initial revision diff -r 445ea777560a -r 6028dcac27ef src/expm.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/expm.cc Wed Aug 11 20:48:00 1993 +0000 @@ -0,0 +1,321 @@ +// tc-expm.cc -*- C++ -*- +/* + +Copyright (C) 1993 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, 675 Mass Ave, Cambridge, MA 02139, USA. + +*/ + +// Written by A. S. Hodel + +#ifdef __GNUG__ +#pragma implementation +#endif + +#include + +#include "Matrix.h" + +#include "tree-const.h" +#include "user-prefs.h" +#include "gripes.h" +#include "error.h" +#include "f-expm.h" + +#ifdef WITH_DLD +tree_constant * +builtin_matrix_exp_2 (tree_constant *args, int nargin, int nargout) +{ + tree_constant *retval = new tree_constant [2]; + retval[0] = matrix_exp (args[1]); + return retval; +} +#endif + +extern "C" +{ + double F77_FCN (dlange) (const char*, const int*, const int*, + const double*, const int*, double*); + + double F77_FCN (zlange) (const char*, const int*, const int*, + const Complex*, const int*, double*); +} + +// XXX FIXME XXX +extern int empty_arg (tree_constant&); + +tree_constant +matrix_exp (const tree_constant& a) +{ + tree_constant retval; + tree_constant tmp = a.make_numeric (); + +// Constants for matrix exponential calculation. + + static double padec[] = + { + 5.0000000000000000e-1, + 1.1666666666666667e-1, + 1.6666666666666667e-2, + 1.6025641025641026e-3, + 1.0683760683760684e-4, + 4.8562548562548563e-6, + 1.3875013875013875e-7, + 1.9270852604185938e-9, + }; + + if (empty_arg (tmp)) + { + int flag = user_pref.propagate_empty_matrices; + if (flag != 0) + { + if (flag < 0) + gripe_empty_arg ("expm", 0); + Matrix m; + retval = tree_constant (m); + } + else gripe_empty_arg ("expm", 1); + } + else if (tmp.rows () != tmp.columns ()) + gripe_square_matrix_required ("expm"); + else + { + int i, j; + int n_cols = tmp.columns (); + + char* balance_job = "B"; // variables for balancing + + int sqpow; // power for scaling and squaring + double inf_norm; // norm of preconditioned matrix + int minus_one_j; // used in computing pade approx + + switch (tmp.const_type ()) + { + case tree_constant_rep::complex_matrix_constant: + { + ComplexMatrix m = tmp.complex_matrix_value (); + Complex trshift = 0.0; // trace shift value + +// Preconditioning step 1: trace normalization. + + for (i = 0; i < n_cols; i++) + trshift += m.elem (i, i); + trshift /= n_cols; + for (i = 0; i < n_cols; i++) + m.elem (i, i) -= trshift; + +// Preconditioning step 2: eigenvalue balancing. + + ComplexAEPBALANCE mbal (m, balance_job); + m = mbal.balanced_matrix (); + ComplexMatrix d = mbal.balancing_matrix (); + +// Preconditioning step 3: scaling. + + ColumnVector work (n_cols); + inf_norm = F77_FCN (zlange) ("I", &n_cols, &n_cols, m. + fortran_vec (), &n_cols, + work.fortran_vec ()); + + sqpow = 1.0 + log (inf_norm) / log (2.0); + +// Check whether we need to square at all. + + if (sqpow <= 0.0) + sqpow = 0.0; + else + { + for (inf_norm = 1.0, i = 0; i < sqpow; i++) + inf_norm *= 2.0; + + m = m / inf_norm; + } + +// npp, dpp: pade' approx polynomial matrices. + + ComplexMatrix npp (n_cols, n_cols, 0.0); + ComplexMatrix dpp = npp; + +// Now powers a^8 ... a^1. + + minus_one_j = -1; + for (j = 7; j >= 0; j--) + { + npp = m * npp + m * padec[j]; + dpp = m * dpp + m * (minus_one_j * padec[j]); + minus_one_j *= -1; + } + +// Zero power. + + dpp = -dpp; + for (j = 0; j < n_cols; j++) + { + npp.elem (j, j) += 1.0; + dpp.elem (j, j) += 1.0; + } + +// Compute pade approximation = inverse (dpp) * npp. + + ComplexMatrix result = dpp.solve (npp); + +// Reverse preconditioning step 3: repeated squaring. + + while (sqpow) + { + result = result * result; + sqpow--; + } + +// reverse preconditioning step 2: inverse balancing XXX FIXME XXX: +// should probably do this with lapack calls instead of a complete +// matrix inversion. + + result = result.transpose (); + d = d.transpose (); + result = result * d; + result = d.solve (result); + result = result.transpose (); + +// Reverse preconditioning step 1: fix trace normalization. + + result = result * exp (trshift); + + retval = tree_constant (result); + } + break; + case tree_constant_rep::complex_scalar_constant: + { + Complex c = tmp.complex_value (); + retval = tree_constant (exp (c)); + } + break; + case tree_constant_rep::matrix_constant: + { + +// Compute the exponential. + + Matrix m = tmp.matrix_value (); + + double trshift = 0; // trace shift value + +// Preconditioning step 1: trace normalization. + + for (i = 0; i < n_cols; i++) + trshift += m.elem (i, i); + trshift /= n_cols; + for (i = 0; i < n_cols; i++) + m.elem (i, i) -= trshift; + +// Preconditioning step 2: balancing. + + AEPBALANCE mbal (m, balance_job); + m = mbal.balanced_matrix (); + Matrix d = mbal.balancing_matrix (); + +// Preconditioning step 3: scaling. + + ColumnVector work(n_cols); + inf_norm = F77_FCN (dlange) ("I", &n_cols, &n_cols, + m.fortran_vec (), &n_cols, + work.fortran_vec ()); + + sqpow = 1.0 + log (inf_norm) / log (2.0); + +// Check whether we need to square at all. + + if (sqpow <= 0.0) + sqpow = 0.0; + else + { + for (inf_norm = 1.0, i = 0; i < sqpow; i++) + inf_norm *= 2.0; + + m = m / inf_norm; + } + +// npp, dpp: pade' approx polynomial matrices. + + Matrix npp (n_cols, n_cols, 0.0); + Matrix dpp = npp; + +// now powers a^8 ... a^1. + + minus_one_j = -1; + for (j = 7; j >= 0; j--) + { + npp = m * npp + m * padec[j]; + dpp = m * dpp + m * (minus_one_j * padec[j]); + minus_one_j *= -1; + } +// Zero power. + + dpp = -dpp; + for(j = 0; j < n_cols; j++) + { + npp.elem (j, j) += 1.0; + dpp.elem (j, j) += 1.0; + } + +// Compute pade approximation = inverse (dpp) * npp. + + Matrix result = dpp.solve (npp); + +// Reverse preconditioning step 3: repeated squaring. + + while(sqpow) + { + result = result * result; + sqpow--; + } + +// Reverse preconditioning step 2: inverse balancing. + + result = result.transpose(); + d = d.transpose (); + result = result * d; + result = d.solve (result); + result = result.transpose (); + +// Reverse preconditioning step 1: fix trace normalization. + + result = result * exp (trshift); + + retval = tree_constant (result); + } + break; + case tree_constant_rep::scalar_constant: + { + double d = tmp.double_value (); + retval = tree_constant (exp (d)); + } + break; + default: + panic_impossible(); + break; + } + } + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/