home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / src / f-expm.cc < prev    next >
C/C++ Source or Header  |  1996-09-28  |  7KB  |  315 lines
  1. // f-expm.cc                                           -*- C++ -*-
  2. /*
  3.  
  4. Copyright (C) 1993, 1994, 1995 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 HAVE_CONFIG_H
  27. #include "config.h"
  28. #endif
  29.  
  30. #include <math.h>
  31.  
  32. #include "dMatrix.h"
  33. #include "CMatrix.h"
  34. #include "CColVector.h"
  35. #include "dbleAEPBAL.h"
  36. #include "CmplxAEPBAL.h"
  37. #include "f77-uscore.h"
  38.  
  39. #include "tree-const.h"
  40. #include "user-prefs.h"
  41. #include "gripes.h"
  42. #include "error.h"
  43. #include "utils.h"
  44. #include "help.h"
  45. #include "defun-dld.h"
  46.  
  47. extern "C"
  48. {
  49.   double F77_FCN (dlange) (const char*, const int*, const int*,
  50.                const double*, const int*, double*);
  51.  
  52.   double F77_FCN (zlange) (const char*, const int*, const int*,
  53.                const Complex*, const int*, double*);
  54. }
  55.  
  56. DEFUN_DLD_BUILTIN ("expm", Fexpm, Sexpm, 2, 1,
  57.   "expm (X): matrix exponential, e^A")
  58. {
  59.   Octave_object retval;
  60.  
  61.   int nargin = args.length ();
  62.  
  63.   if (nargin != 1)
  64.     {
  65.       print_usage ("expm");
  66.       return retval;
  67.     }
  68.  
  69.   tree_constant arg = args(0);
  70.  
  71. // Constants for matrix exponential calculation.
  72.  
  73.   static double padec [] =
  74.     {
  75.       5.0000000000000000e-1,
  76.       1.1666666666666667e-1,
  77.       1.6666666666666667e-2,
  78.       1.6025641025641026e-3,
  79.       1.0683760683760684e-4,
  80.       4.8562548562548563e-6,
  81.       1.3875013875013875e-7,
  82.       1.9270852604185938e-9,
  83.     };
  84.  
  85.   int nr = arg.rows ();
  86.   int nc = arg.columns ();
  87.  
  88.   int arg_is_empty = empty_arg ("expm", nr, nc);
  89.  
  90.   if (arg_is_empty < 0)
  91.     return retval;
  92.   if (arg_is_empty > 0)
  93.     return Matrix ();
  94.  
  95.   if (nr != nc)
  96.     {
  97.       gripe_square_matrix_required ("expm");
  98.       return retval;
  99.     }
  100.  
  101.   int i, j;
  102.  
  103.   char* balance_job = "B";    // variables for balancing
  104.  
  105.   int sqpow;        // power for scaling and squaring
  106.   double inf_norm;    // norm of preconditioned matrix
  107.   int minus_one_j;    // used in computing pade approx
  108.  
  109.   if (arg.is_real_type ())
  110.     {
  111.  
  112. // Compute the exponential.
  113.  
  114.       Matrix m = arg.matrix_value ();
  115.  
  116.       if (error_state)
  117.     return retval;
  118.  
  119.       double trshift = 0;        // trace shift value
  120.  
  121. // Preconditioning step 1: trace normalization.
  122.  
  123.       for (i = 0; i < nc; i++)
  124.     trshift += m.elem (i, i);
  125.       trshift /= nc;
  126.       for (i = 0; i < nc; i++)
  127.     m.elem (i, i) -= trshift;
  128.  
  129. // Preconditioning step 2: balancing.
  130.  
  131.       AEPBALANCE mbal (m, balance_job);
  132.       m = mbal.balanced_matrix ();
  133.       Matrix d = mbal.balancing_matrix ();
  134.  
  135. // Preconditioning step 3: scaling.
  136.  
  137.       ColumnVector work(nc);
  138.       inf_norm = F77_FCN (dlange) ("I", &nc, &nc,
  139.                    m.fortran_vec (), &nc,
  140.                    work.fortran_vec ());
  141.  
  142.       sqpow = (int) (1.0 + log (inf_norm) / log (2.0));
  143.  
  144. // Check whether we need to square at all.
  145.  
  146.       if (sqpow < 0)
  147.     sqpow = 0;
  148.       else
  149.     {
  150.       for (inf_norm = 1.0, i = 0; i < sqpow; i++)
  151.         inf_norm *= 2.0;
  152.  
  153.       m = m / inf_norm;
  154.     }
  155.  
  156. // npp, dpp: pade' approx polynomial matrices.
  157.  
  158.       Matrix npp (nc, nc, 0.0);
  159.       Matrix dpp = npp;
  160.  
  161. // now powers a^8 ... a^1.
  162.  
  163.       minus_one_j = -1;
  164.       for (j = 7; j >= 0; j--)
  165.     {
  166.       npp = m * npp + m * padec[j];
  167.       dpp = m * dpp + m * (minus_one_j * padec[j]);
  168.       minus_one_j *= -1;
  169.     }
  170. // Zero power.
  171.  
  172.       dpp = -dpp;
  173.       for(j = 0; j < nc; j++)
  174.     {
  175.       npp.elem (j, j) += 1.0;
  176.       dpp.elem (j, j) += 1.0;
  177.     }
  178.  
  179. // Compute pade approximation = inverse (dpp) * npp.
  180.  
  181.       Matrix result = dpp.solve (npp);
  182.  
  183. // Reverse preconditioning step 3: repeated squaring.
  184.  
  185.       while (sqpow)
  186.     {
  187.       result = result * result;
  188.       sqpow--;
  189.     }
  190.  
  191. // Reverse preconditioning step 2: inverse balancing.
  192.  
  193.       result = result.transpose();
  194.       d = d.transpose ();
  195.       result = result * d;
  196.       result = d.solve (result);
  197.       result = result.transpose ();
  198.  
  199. // Reverse preconditioning step 1: fix trace normalization.
  200.  
  201.       result = result * exp (trshift);
  202.  
  203.       retval = result;
  204.     }
  205.   else if (arg.is_complex_type ())
  206.     {
  207.       ComplexMatrix m = arg.complex_matrix_value ();
  208.  
  209.       if (error_state)
  210.     return retval;
  211.  
  212.       Complex trshift = 0.0;        // trace shift value
  213.  
  214. // Preconditioning step 1: trace normalization.
  215.  
  216.       for (i = 0; i < nc; i++)
  217.     trshift += m.elem (i, i);
  218.       trshift /= nc;
  219.       for (i = 0; i < nc; i++)
  220.     m.elem (i, i) -= trshift;
  221.  
  222. // Preconditioning step 2: eigenvalue balancing.
  223.  
  224.       ComplexAEPBALANCE mbal (m, balance_job);
  225.       m = mbal.balanced_matrix ();
  226.       ComplexMatrix d = mbal.balancing_matrix ();
  227.  
  228. // Preconditioning step 3: scaling.
  229.  
  230.       ColumnVector work (nc);
  231.       inf_norm = F77_FCN (zlange) ("I", &nc, &nc, m.
  232.                    fortran_vec (), &nc,
  233.                    work.fortran_vec ());
  234.  
  235.       sqpow = (int) (1.0 + log (inf_norm) / log (2.0));
  236.  
  237. // Check whether we need to square at all.
  238.  
  239.       if (sqpow < 0)
  240.     sqpow = 0;
  241.       else
  242.     {
  243.       for (inf_norm = 1.0, i = 0; i < sqpow; i++)
  244.         inf_norm *= 2.0;
  245.  
  246.       m = m / inf_norm;
  247.     }
  248.  
  249. // npp, dpp: pade' approx polynomial matrices.
  250.  
  251.       ComplexMatrix npp (nc, nc, 0.0);
  252.       ComplexMatrix dpp = npp;
  253.  
  254. // Now powers a^8 ... a^1.
  255.  
  256.       minus_one_j = -1;
  257.       for (j = 7; j >= 0; j--)
  258.     {
  259.       npp = m * npp + m * padec[j];
  260.       dpp = m * dpp + m * (minus_one_j * padec[j]);
  261.       minus_one_j *= -1;
  262.     }
  263.  
  264. // Zero power.
  265.  
  266.       dpp = -dpp;
  267.       for (j = 0; j < nc; j++)
  268.     {
  269.       npp.elem (j, j) += 1.0;
  270.       dpp.elem (j, j) += 1.0;
  271.     }
  272.  
  273. // Compute pade approximation = inverse (dpp) * npp.
  274.  
  275.       ComplexMatrix result = dpp.solve (npp);
  276.     
  277. // Reverse preconditioning step 3: repeated squaring.
  278.  
  279.       while (sqpow)
  280.     {
  281.       result = result * result;
  282.       sqpow--;
  283.     }
  284.  
  285. // reverse preconditioning step 2: inverse balancing XXX FIXME XXX:
  286. // should probably do this with lapack calls instead of a complete
  287. // matrix inversion.
  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 = result;
  300.     }
  301.   else
  302.     {
  303.       gripe_wrong_type_arg ("expm", arg);
  304.     }
  305.  
  306.   return retval;
  307. }
  308.  
  309. /*
  310. ;;; Local Variables: ***
  311. ;;; mode: C++ ***
  312. ;;; page-delimiter: "^/\\*" ***
  313. ;;; End: ***
  314. */
  315.