Mercurial > octave-nkf

changeset 8392:c187f0e3a7ee

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
use m-file implementation for expm
author Jaroslav Hajek <highegg@gmail.com>
date Wed, 10 Dec 2008 11:55:23 +0100
parents 343f0fbca6eb
children c8a785d0e867
files liboctave/CMatrix.cc liboctave/CMatrix.h liboctave/ChangeLog liboctave/dMatrix.cc liboctave/dMatrix.h liboctave/fCMatrix.cc liboctave/fCMatrix.h liboctave/fMatrix.cc liboctave/fMatrix.h scripts/linear-algebra/Makefile.in src/ChangeLog src/DLD-FUNCTIONS/expm.cc src/Makefile.in
diffstat 13 files changed, 14 insertions(+), 1002 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/CMatrix.cc	Wed Dec 10 11:55:21 2008 +0100
+++ b/liboctave/CMatrix.cc	Wed Dec 10 11:55:23 2008 +0100
@@ -2926,194 +2926,6 @@
      rcon);
 }
 
-ComplexMatrix
-ComplexMatrix::expm (void) const
-{
-  ComplexMatrix retval;
-
-  ComplexMatrix m = *this;
-
-  octave_idx_type nc = columns ();
-
-  // Preconditioning step 1: trace normalization to reduce dynamic
-  // range of poles, but avoid making stable eigenvalues unstable.
-
-  // trace shift value
-  Complex trshift = 0.0;
-
-  for (octave_idx_type i = 0; i < nc; i++)
-    trshift += m.elem (i, i);
-
-  trshift /= nc;
-
-  if (trshift.real () < 0.0)
-    {
-      trshift = trshift.imag ();
-      if (trshift.real () > 709.0)
-	trshift = 709.0;
-    }
-
-  for (octave_idx_type i = 0; i < nc; i++)
-    m.elem (i, i) -= trshift;
-
-  // Preconditioning step 2: eigenvalue balancing.
-  // code follows development in AEPBAL
-
-  Complex *mp = m.fortran_vec ();
-
-  octave_idx_type info, ilo, ihi,ilos,ihis;
-  Array<double> dpermute (nc);
-  Array<double> dscale (nc);
-
-  // FIXME -- should pass job as a parameter in expm
-
-  // Permute first
-  char job = 'P';
-  F77_XFCN (zgebal, ZGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
-			     nc, mp, nc, ilo, ihi,
-			     dpermute.fortran_vec (), info
-			     F77_CHAR_ARG_LEN (1)));
-
-  // then scale
-  job = 'S';
-  F77_XFCN (zgebal, ZGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
-			     nc, mp, nc, ilos, ihis,
-			     dscale.fortran_vec (), info
-			     F77_CHAR_ARG_LEN (1)));
-
-  // Preconditioning step 3: scaling.
-
-  ColumnVector work (nc);
-  double inf_norm;
-
-  F77_XFCN (xzlange, XZLANGE, (F77_CONST_CHAR_ARG2 ("I", 1),
-			       nc, nc, m.fortran_vec (), nc,
-			       work.fortran_vec (), inf_norm
-			       F77_CHAR_ARG_LEN (1)));
-
-  int sqpow = (inf_norm > 0.0
-	       ? static_cast<int> (1.0 + log (inf_norm) / log (2.0)) : 0);
-
-  // Check whether we need to square at all.
-
-  if (sqpow < 0)
-    sqpow = 0;
-
-  if (sqpow > 0)
-    {
-      if (sqpow > 1023)
-	sqpow = 1023;
-
-      double scale_factor = 1.0;
-      for (octave_idx_type i = 0; i < sqpow; i++)
-	scale_factor *= 2.0;
-
-      m = m / scale_factor;
-    }
-
-  // npp, dpp: pade' approx polynomial matrices.
-
-  ComplexMatrix npp (nc, nc, 0.0);
-  Complex *pnpp = npp.fortran_vec ();
-  ComplexMatrix dpp = npp;
-  Complex *pdpp = dpp.fortran_vec ();
-
-  // Now powers a^8 ... a^1.
-
-  int minus_one_j = -1;
-  for (octave_idx_type j = 7; j >= 0; j--)
-    {
-      for (octave_idx_type i = 0; i < nc; i++)
-	{
-	  octave_idx_type k = i * nc + i;
-	  pnpp[k] += padec[j];
-	  pdpp[k] += minus_one_j * padec[j];
-	}      
-
-      npp = m * npp;
-      pnpp = npp.fortran_vec ();
-
-      dpp = m * dpp;
-      pdpp = dpp.fortran_vec ();
-
-      minus_one_j *= -1;
-    }
-
-  // Zero power.
-
-  dpp = -dpp;
-  for (octave_idx_type j = 0; j < nc; j++)
-    {
-      npp.elem (j, j) += 1.0;
-      dpp.elem (j, j) += 1.0;
-    }
-
-  // Compute pade approximation = inverse (dpp) * npp.
-
-  double rcon;
-  retval = dpp.solve (npp, info, rcon, solve_singularity_warning);
-
-  if (info < 0)
-    return retval;
-
-  // Reverse preconditioning step 3: repeated squaring.
-
-  while (sqpow)
-    {
-      retval = retval * retval;
-      sqpow--;
-    }
-
-  // Reverse preconditioning step 2: inverse balancing.
-  // Done in two steps: inverse scaling, then inverse permutation
-
-  // inverse scaling (diagonal transformation)
-  for (octave_idx_type i = 0; i < nc; i++)
-    for (octave_idx_type j = 0; j < nc; j++)
-      retval(i,j) *= dscale(i) / dscale(j);
-
-  OCTAVE_QUIT;
-
-  // construct balancing permutation vector
-  Array<octave_idx_type> iperm (nc);
-  for (octave_idx_type i = 0; i < nc; i++)
-    iperm(i) = i;  // identity permutation
-
-  // trailing permutations must be done in reverse order
-  for (octave_idx_type i = nc - 1; i >= ihi; i--)
-    {
-      octave_idx_type swapidx = static_cast<octave_idx_type> (dpermute(i)) - 1;
-      octave_idx_type tmp = iperm(i);
-      iperm(i) = iperm(swapidx);
-      iperm(swapidx) = tmp;
-    }
-
-  // leading permutations in forward order
-  for (octave_idx_type i = 0; i < (ilo-1); i++)
-    {
-      octave_idx_type swapidx = static_cast<octave_idx_type> (dpermute(i)) - 1;
-      octave_idx_type tmp = iperm(i);
-      iperm(i) = iperm (swapidx);
-      iperm(swapidx) = tmp;
-    }
-
-  // construct inverse balancing permutation vector
-  Array<octave_idx_type> invpvec (nc);
-  for (octave_idx_type i = 0; i < nc; i++)
-    invpvec(iperm(i)) = i;     // Thanks to R. A. Lippert for this method
-
-  OCTAVE_QUIT;
-
-  ComplexMatrix tmpMat = retval;
-  for (octave_idx_type i = 0; i < nc; i++)
-    for (octave_idx_type j = 0; j < nc; j++)
-      retval(i,j) = tmpMat(invpvec(i),invpvec(j));
-
-  // Reverse preconditioning step 1: fix trace normalization. 
-
-  return exp (trshift) * retval;
-}
-
 // column vector by row vector -> matrix operations
 
 ComplexMatrix
--- a/liboctave/CMatrix.h	Wed Dec 10 11:55:21 2008 +0100
+++ b/liboctave/CMatrix.h	Wed Dec 10 11:55:23 2008 +0100
@@ -301,8 +301,6 @@
 			       octave_idx_type& info,
 			       octave_idx_type& rank, double& rcon) const;
 
-  ComplexMatrix expm (void) const;
-
   // matrix by diagonal matrix -> matrix operations
 
   ComplexMatrix& operator += (const DiagMatrix& a);
--- a/liboctave/ChangeLog	Wed Dec 10 11:55:21 2008 +0100
+++ b/liboctave/ChangeLog	Wed Dec 10 11:55:23 2008 +0100
@@ -1,3 +1,10 @@
+2008-12-10  Jaroslav Hajek  <highegg@gmail.com>
+
+	* dMatrix.h, dMatrix.cc (Matrix::expm): Remove.
+	* fMatrix.h, fMatrix.cc (FloatMatrix::expm): Remove.
+	* CMatrix.h, CMatrix.cc (ComplexMatrix::expm): Remove.
+	* fCMatrix.h, fCMatrix.cc (FloatComplexMatrix::expm): Remove.
+
 2008-12-09  Jaroslav Hajek  <highegg@gmail.com>
 
 	* base-aepbal.h: New source.
--- a/liboctave/dMatrix.cc	Wed Dec 10 11:55:21 2008 +0100
+++ b/liboctave/dMatrix.cc	Wed Dec 10 11:55:23 2008 +0100
@@ -2575,192 +2575,6 @@
      rcon);
 }
 
-Matrix
-Matrix::expm (void) const
-{
-  Matrix retval;
-
-  Matrix m = *this;
-
-  if (numel () == 1)
-    return Matrix (1, 1, exp (m(0)));
-
-  octave_idx_type nc = columns ();
-
-  // Preconditioning step 1: trace normalization to reduce dynamic
-  // range of poles, but avoid making stable eigenvalues unstable.
-
-  // trace shift value
-  volatile double trshift = 0.0;
-
-  for (octave_idx_type i = 0; i < nc; i++)
-    trshift += m.elem (i, i);
-
-  trshift /= nc;
-
-  if (trshift > 0.0)
-    {
-      for (octave_idx_type i = 0; i < nc; i++)
-	m.elem (i, i) -= trshift;
-    }
-
-  // Preconditioning step 2: balancing; code follows development
-  // in AEPBAL
-
-  double *p_m = m.fortran_vec ();
-
-  octave_idx_type info, ilo, ihi, ilos, ihis;
-  Array<double> dpermute (nc);
-  Array<double> dscale (nc);
-
-  // permutation first
-  char job = 'P';
-  F77_XFCN (dgebal, DGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
-			     nc, p_m, nc, ilo, ihi,
-			     dpermute.fortran_vec (), info
-			     F77_CHAR_ARG_LEN (1)));
-
-  // then scaling
-  job = 'S';
-  F77_XFCN (dgebal, DGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
-			     nc, p_m, nc, ilos, ihis,
-			     dscale.fortran_vec (), info
-			     F77_CHAR_ARG_LEN (1)));
-
-  // Preconditioning step 3: scaling.
-  
-  ColumnVector work(nc);
-  double inf_norm;
-  
-  F77_XFCN (xdlange, XDLANGE, (F77_CONST_CHAR_ARG2 ("I", 1),
-			       nc, nc, m.fortran_vec (), nc,
-			       work.fortran_vec (), inf_norm
-			       F77_CHAR_ARG_LEN (1)));
-  
-  octave_idx_type sqpow = static_cast<octave_idx_type> (inf_norm > 0.0
-		     ? (1.0 + log (inf_norm) / log (2.0))
-		     : 0.0);
-  
-  // Check whether we need to square at all.
-  
-  if (sqpow < 0)
-    sqpow = 0;
-  
-  if (sqpow > 0)
-    {
-      if (sqpow > 1023)
-	sqpow = 1023;
-
-      double scale_factor = 1.0;
-      for (octave_idx_type i = 0; i < sqpow; i++)
-	scale_factor *= 2.0;
-
-      m = m / scale_factor;
-    }
-  
-  // npp, dpp: pade' approx polynomial matrices.
-  
-  Matrix npp (nc, nc, 0.0);
-  double *pnpp = npp.fortran_vec ();
-  Matrix dpp = npp;
-  double *pdpp = dpp.fortran_vec ();
-  
-  // Now powers a^8 ... a^1.
-  
-  octave_idx_type minus_one_j = -1;
-  for (octave_idx_type j = 7; j >= 0; j--)
-    {
-      for (octave_idx_type i = 0; i < nc; i++)
-	{
-	  octave_idx_type k = i * nc + i;
-	  pnpp[k] += padec[j];
-	  pdpp[k] += minus_one_j * padec[j];
-	}      
-
-      npp = m * npp;
-      pnpp = npp.fortran_vec ();
-
-      dpp = m * dpp;
-      pdpp = dpp.fortran_vec ();
-
-      minus_one_j *= -1;
-    }
-  
-  // Zero power.
-  
-  dpp = -dpp;
-  for (octave_idx_type j = 0; j < nc; j++)
-    {
-      npp.elem (j, j) += 1.0;
-      dpp.elem (j, j) += 1.0;
-    }
-  
-  // Compute pade approximation = inverse (dpp) * npp.
-
-  double rcon;
-  retval = dpp.solve (npp, info, rcon, solve_singularity_warning);
-
-  if (info < 0)
-    return retval;
-
-  // Reverse preconditioning step 3: repeated squaring.
-  
-  while (sqpow)
-    {
-      retval = retval * retval;
-      sqpow--;
-    }
-  
-  // Reverse preconditioning step 2: inverse balancing.
-  // apply inverse scaling to computed exponential
-  for (octave_idx_type i = 0; i < nc; i++)
-    for (octave_idx_type j = 0; j < nc; j++)
-      retval(i,j) *= dscale(i) / dscale(j);
-
-  OCTAVE_QUIT;
-
-  // construct balancing permutation vector
-  Array<octave_idx_type> iperm (nc);
-  for (octave_idx_type i = 0; i < nc; i++)
-    iperm(i) = i;  // identity permutation
-
-  // trailing permutations must be done in reverse order
-  for (octave_idx_type i = nc - 1; i >= ihi; i--)
-    {
-      octave_idx_type swapidx = static_cast<octave_idx_type> (dpermute(i)) - 1;
-      octave_idx_type tmp = iperm(i);
-      iperm(i) = iperm(swapidx);
-      iperm(swapidx) = tmp;
-    }
-
-  // leading permutations in forward order
-  for (octave_idx_type i = 0; i < (ilo-1); i++)
-    {
-      octave_idx_type swapidx = static_cast<octave_idx_type> (dpermute(i)) - 1;
-      octave_idx_type tmp = iperm(i);
-      iperm(i) = iperm (swapidx);
-      iperm(swapidx) = tmp;
-    }
-
-  // construct inverse balancing permutation vector
-  Array<octave_idx_type> invpvec (nc);
-  for (octave_idx_type i = 0; i < nc; i++)
-    invpvec(iperm(i)) = i;     // Thanks to R. A. Lippert for this method
-
-  OCTAVE_QUIT;
-
-  Matrix tmpMat = retval;
-  for (octave_idx_type i = 0; i < nc; i++)
-    for (octave_idx_type j = 0; j < nc; j++)
-      retval(i,j) = tmpMat(invpvec(i),invpvec(j));
-
-  // Reverse preconditioning step 1: fix trace normalization. 
-  if (trshift > 0.0)
-    retval = exp (trshift) * retval;
-
-  return retval;
-}
-
 Matrix&
 Matrix::operator += (const DiagMatrix& a)
 {
--- a/liboctave/dMatrix.h	Wed Dec 10 11:55:21 2008 +0100
+++ b/liboctave/dMatrix.h	Wed Dec 10 11:55:23 2008 +0100
@@ -267,8 +267,6 @@
 			       octave_idx_type& info,
 			       octave_idx_type& rank, double& rcon) const;
 
-  Matrix expm (void) const;
-
   Matrix& operator += (const DiagMatrix& a);
   Matrix& operator -= (const DiagMatrix& a);
 
--- a/liboctave/fCMatrix.cc	Wed Dec 10 11:55:21 2008 +0100
+++ b/liboctave/fCMatrix.cc	Wed Dec 10 11:55:23 2008 +0100
@@ -2920,194 +2920,6 @@
      rcon);
 }
 
-FloatComplexMatrix
-FloatComplexMatrix::expm (void) const
-{
-  FloatComplexMatrix retval;
-
-  FloatComplexMatrix m = *this;
-
-  octave_idx_type nc = columns ();
-
-  // Preconditioning step 1: trace normalization to reduce dynamic
-  // range of poles, but avoid making stable eigenvalues unstable.
-
-  // trace shift value
-  FloatComplex trshift = 0.0;
-
-  for (octave_idx_type i = 0; i < nc; i++)
-    trshift += m.elem (i, i);
-
-  trshift /= nc;
-
-  if (trshift.real () < 0.0)
-    {
-      trshift = trshift.imag ();
-      if (trshift.real () > 709.0)
-	trshift = 709.0;
-    }
-
-  for (octave_idx_type i = 0; i < nc; i++)
-    m.elem (i, i) -= trshift;
-
-  // Preconditioning step 2: eigenvalue balancing.
-  // code follows development in AEPBAL
-
-  FloatComplex *mp = m.fortran_vec ();
-
-  octave_idx_type info, ilo, ihi,ilos,ihis;
-  Array<float> dpermute (nc);
-  Array<float> dscale (nc);
-
-  // FIXME -- should pass job as a parameter in expm
-
-  // Permute first
-  char job = 'P';
-  F77_XFCN (cgebal, CGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
-			     nc, mp, nc, ilo, ihi,
-			     dpermute.fortran_vec (), info
-			     F77_CHAR_ARG_LEN (1)));
-
-  // then scale
-  job = 'S';
-  F77_XFCN (cgebal, CGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
-			     nc, mp, nc, ilos, ihis,
-			     dscale.fortran_vec (), info
-			     F77_CHAR_ARG_LEN (1)));
-
-  // Preconditioning step 3: scaling.
-
-  FloatColumnVector work (nc);
-  float inf_norm;
-
-  F77_XFCN (xclange, XCLANGE, (F77_CONST_CHAR_ARG2 ("I", 1),
-			       nc, nc, m.fortran_vec (), nc,
-			       work.fortran_vec (), inf_norm
-			       F77_CHAR_ARG_LEN (1)));
-
-  int sqpow = (inf_norm > 0.0
-	       ? static_cast<int> (1.0 + log (inf_norm) / log (2.0)) : 0);
-
-  // Check whether we need to square at all.
-
-  if (sqpow < 0)
-    sqpow = 0;
-
-  if (sqpow > 0)
-    {
-      if (sqpow > 1023)
-	sqpow = 1023;
-
-      float scale_factor = 1.0;
-      for (octave_idx_type i = 0; i < sqpow; i++)
-	scale_factor *= 2.0;
-
-      m = m / scale_factor;
-    }
-
-  // npp, dpp: pade' approx polynomial matrices.
-
-  FloatComplexMatrix npp (nc, nc, 0.0);
-  FloatComplex *pnpp = npp.fortran_vec ();
-  FloatComplexMatrix dpp = npp;
-  FloatComplex *pdpp = dpp.fortran_vec ();
-
-  // Now powers a^8 ... a^1.
-
-  int minus_one_j = -1;
-  for (octave_idx_type j = 7; j >= 0; j--)
-    {
-      for (octave_idx_type i = 0; i < nc; i++)
-	{
-	  octave_idx_type k = i * nc + i;
-	  pnpp[k] += padec[j];
-	  pdpp[k] += minus_one_j * padec[j];
-	}      
-
-      npp = m * npp;
-      pnpp = npp.fortran_vec ();
-
-      dpp = m * dpp;
-      pdpp = dpp.fortran_vec ();
-
-      minus_one_j *= -1;
-    }
-
-  // Zero power.
-
-  dpp = -dpp;
-  for (octave_idx_type j = 0; j < nc; j++)
-    {
-      npp.elem (j, j) += 1.0;
-      dpp.elem (j, j) += 1.0;
-    }
-
-  // Compute pade approximation = inverse (dpp) * npp.
-
-  float rcon;
-  retval = dpp.solve (npp, info, rcon, solve_singularity_warning);
-
-  if (info < 0)
-    return retval;
-
-  // Reverse preconditioning step 3: repeated squaring.
-
-  while (sqpow)
-    {
-      retval = retval * retval;
-      sqpow--;
-    }
-
-  // Reverse preconditioning step 2: inverse balancing.
-  // Done in two steps: inverse scaling, then inverse permutation
-
-  // inverse scaling (diagonal transformation)
-  for (octave_idx_type i = 0; i < nc; i++)
-    for (octave_idx_type j = 0; j < nc; j++)
-      retval(i,j) *= dscale(i) / dscale(j);
-
-  OCTAVE_QUIT;
-
-  // construct balancing permutation vector
-  Array<octave_idx_type> iperm (nc);
-  for (octave_idx_type i = 0; i < nc; i++)
-    iperm(i) = i;  // identity permutation
-
-  // trailing permutations must be done in reverse order
-  for (octave_idx_type i = nc - 1; i >= ihi; i--)
-    {
-      octave_idx_type swapidx = static_cast<octave_idx_type> (dpermute(i)) - 1;
-      octave_idx_type tmp = iperm(i);
-      iperm(i) = iperm(swapidx);
-      iperm(swapidx) = tmp;
-    }
-
-  // leading permutations in forward order
-  for (octave_idx_type i = 0; i < (ilo-1); i++)
-    {
-      octave_idx_type swapidx = static_cast<octave_idx_type> (dpermute(i)) - 1;
-      octave_idx_type tmp = iperm(i);
-      iperm(i) = iperm (swapidx);
-      iperm(swapidx) = tmp;
-    }
-
-  // construct inverse balancing permutation vector
-  Array<octave_idx_type> invpvec (nc);
-  for (octave_idx_type i = 0; i < nc; i++)
-    invpvec(iperm(i)) = i;     // Thanks to R. A. Lippert for this method
-
-  OCTAVE_QUIT;
-
-  FloatComplexMatrix tmpMat = retval;
-  for (octave_idx_type i = 0; i < nc; i++)
-    for (octave_idx_type j = 0; j < nc; j++)
-      retval(i,j) = tmpMat(invpvec(i),invpvec(j));
-
-  // Reverse preconditioning step 1: fix trace normalization. 
-
-  return exp (trshift) * retval;
-}
-
 // column vector by row vector -> matrix operations
 
 FloatComplexMatrix
--- a/liboctave/fCMatrix.h	Wed Dec 10 11:55:21 2008 +0100
+++ b/liboctave/fCMatrix.h	Wed Dec 10 11:55:23 2008 +0100
@@ -301,8 +301,6 @@
 			       octave_idx_type& info,
 			       octave_idx_type& rank, float& rcon) const;
 
-  FloatComplexMatrix expm (void) const;
-
   // matrix by diagonal matrix -> matrix operations
 
   FloatComplexMatrix& operator += (const FloatDiagMatrix& a);
--- a/liboctave/fMatrix.cc	Wed Dec 10 11:55:21 2008 +0100
+++ b/liboctave/fMatrix.cc	Wed Dec 10 11:55:23 2008 +0100
@@ -2574,193 +2574,6 @@
      rcon);
 }
 
-FloatMatrix
-FloatMatrix::expm (void) const
-{
-  FloatMatrix retval;
-
-  FloatMatrix m = *this;
-
-  if (numel () == 1)
-    return FloatMatrix (1, 1, exp (m(0)));
-
-  octave_idx_type nc = columns ();
-
-  // Preconditioning step 1: trace normalization to reduce dynamic
-  // range of poles, but avoid making stable eigenvalues unstable.
-
-  // trace shift value
-  volatile float trshift = 0.0;
-
-  for (octave_idx_type i = 0; i < nc; i++)
-    trshift += m.elem (i, i);
-
-  trshift /= nc;
-
-  if (trshift > 0.0)
-    {
-      for (octave_idx_type i = 0; i < nc; i++)
-	m.elem (i, i) -= trshift;
-    }
-
-  // Preconditioning step 2: balancing; code follows development
-  // in AEPBAL
-
-  float *p_m = m.fortran_vec ();
-
-  octave_idx_type info, ilo, ihi, ilos, ihis;
-  Array<float> dpermute (nc);
-  Array<float> dscale (nc);
-
-  // permutation first
-  char job = 'P';
-  F77_XFCN (sgebal, SGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
-			     nc, p_m, nc, ilo, ihi,
-			     dpermute.fortran_vec (), info
-			     F77_CHAR_ARG_LEN (1)));
-
-  // then scaling
-  job = 'S';
-  F77_XFCN (sgebal, SGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
-			     nc, p_m, nc, ilos, ihis,
-			     dscale.fortran_vec (), info
-			     F77_CHAR_ARG_LEN (1)));
-
-  // Preconditioning step 3: scaling.
-  
-  FloatColumnVector work(nc);
-  float inf_norm;
-  
-  F77_XFCN (xslange, XSLANGE, (F77_CONST_CHAR_ARG2 ("I", 1),
-			       nc, nc, m.fortran_vec (), nc,
-			       work.fortran_vec (), inf_norm
-			       F77_CHAR_ARG_LEN (1)));
-  
-  octave_idx_type sqpow = static_cast<octave_idx_type> (inf_norm > 0.0
-		     ? (1.0 + log (inf_norm) / log (2.0))
-		     : 0.0);
-  
-  // Check whether we need to square at all.
-  
-  if (sqpow < 0)
-    sqpow = 0;
-  
-  if (sqpow > 0)
-    {
-      if (sqpow > 1023)
-	sqpow = 1023;
-
-      float scale_factor = 1.0;
-      for (octave_idx_type i = 0; i < sqpow; i++)
-	scale_factor *= 2.0;
-
-      m = m / scale_factor;
-    }
-  
-  // npp, dpp: pade' approx polynomial matrices.
-  
-  FloatMatrix npp (nc, nc, 0.0);
-  float *pnpp = npp.fortran_vec ();
-  FloatMatrix dpp = npp;
-  float *pdpp = dpp.fortran_vec ();
-  
-  // Now powers a^8 ... a^1.
-  
-  octave_idx_type minus_one_j = -1;
-  for (octave_idx_type j = 7; j >= 0; j--)
-    {
-      for (octave_idx_type i = 0; i < nc; i++)
-	{
-	  octave_idx_type k = i * nc + i;
-	  pnpp[k] += padec[j];
-	  pdpp[k] += minus_one_j * padec[j];
-	}      
-
-      npp = m * npp;
-      pnpp = npp.fortran_vec ();
-
-      dpp = m * dpp;
-      pdpp = dpp.fortran_vec ();
-
-      minus_one_j *= -1;
-    }
-  
-  // Zero power.
-  
-  dpp = -dpp;
-  for (octave_idx_type j = 0; j < nc; j++)
-    {
-      npp.elem (j, j) += 1.0;
-      dpp.elem (j, j) += 1.0;
-    }
-  
-  // Compute pade approximation = inverse (dpp) * npp.
-
-  float rcon;
-  retval = dpp.solve (npp, info, rcon, solve_singularity_warning);
-
-  if (info < 0)
-    return retval;
-
-  // Reverse preconditioning step 3: repeated squaring.
-  
-  while (sqpow)
-    {
-      retval = retval * retval;
-      sqpow--;
-    }
-  
-  // Reverse preconditioning step 2: inverse balancing.
-  // apply inverse scaling to computed exponential
-  for (octave_idx_type i = 0; i < nc; i++)
-    for (octave_idx_type j = 0; j < nc; j++)
-      retval(i,j) *= dscale(i) / dscale(j);
-
-  OCTAVE_QUIT;
-
-  // construct balancing permutation vector
-  Array<octave_idx_type> iperm (nc);
-  for (octave_idx_type i = 0; i < nc; i++)
-    iperm(i) = i;  // identity permutation
-
-  // trailing permutations must be done in reverse order
-  for (octave_idx_type i = nc - 1; i >= ihi; i--)
-    {
-      octave_idx_type swapidx = static_cast<octave_idx_type> (dpermute(i)) - 1;
-      octave_idx_type tmp = iperm(i);
-      iperm(i) = iperm(swapidx);
-      iperm(swapidx) = tmp;
-    }
-
-  // leading permutations in forward order
-  for (octave_idx_type i = 0; i < (ilo-1); i++)
-    {
-      octave_idx_type swapidx = static_cast<octave_idx_type> (dpermute(i)) - 1;
-      octave_idx_type tmp = iperm(i);
-      iperm(i) = iperm (swapidx);
-      iperm(swapidx) = tmp;
-    }
-
-  // construct inverse balancing permutation vector
-  Array<octave_idx_type> invpvec (nc);
-  for (octave_idx_type i = 0; i < nc; i++)
-    invpvec(iperm(i)) = i;     // Thanks to R. A. Lippert for this method
-
-  OCTAVE_QUIT;
-
-  FloatMatrix tmpMat = retval;
-  for (octave_idx_type i = 0; i < nc; i++)
-    for (octave_idx_type j = 0; j < nc; j++)
-      retval(i,j) = tmpMat(invpvec(i),invpvec(j));
-
-  // Reverse preconditioning step 1: fix trace normalization. 
-  
-  if (trshift > 0.0)
-    retval = expf (trshift) * retval;
-
-  return retval;
-}
-
 FloatMatrix&
 FloatMatrix::operator += (const FloatDiagMatrix& a)
 {
--- a/liboctave/fMatrix.h	Wed Dec 10 11:55:21 2008 +0100
+++ b/liboctave/fMatrix.h	Wed Dec 10 11:55:23 2008 +0100
@@ -268,8 +268,6 @@
 			       octave_idx_type& info,
 			       octave_idx_type& rank, float& rcon) const;
 
-  FloatMatrix expm (void) const;
-
   FloatMatrix& operator += (const FloatDiagMatrix& a);
   FloatMatrix& operator -= (const FloatDiagMatrix& a);
 
--- a/scripts/linear-algebra/Makefile.in	Wed Dec 10 11:55:21 2008 +0100
+++ b/scripts/linear-algebra/Makefile.in	Wed Dec 10 11:55:23 2008 +0100
@@ -34,7 +34,7 @@
 INSTALL_DATA = @INSTALL_DATA@
 
 SOURCES = commutation_matrix.m cond.m condest.m cross.m \
-  dmult.m dot.m duplication_matrix.m housh.m krylov.m krylovb.m logm.m \
+  dmult.m dot.m duplication_matrix.m expm.m housh.m krylov.m krylovb.m logm.m \
   null.m onenormest.m orth.m planerot.m qzhess.m rank.m rref.m subspace.m \
   trace.m vec.m vech.m
 
--- a/src/ChangeLog	Wed Dec 10 11:55:21 2008 +0100
+++ b/src/ChangeLog	Wed Dec 10 11:55:23 2008 +0100
@@ -1,3 +1,8 @@
+2008-12-10  Jaroslav Hajek  <highegg@gmail.com>
+
+	* DLD-FUNCTIONS/expm.cc: Remove.
+	* Makefile.in: Update.
+
 2008-12-10  Jaroslav Hajek  <highegg@gmail.com>
 
 	* ov-intx.h (OCTAVE_VALUE_INT_SCALAR_T::empty_clone): Construct an
--- a/src/DLD-FUNCTIONS/expm.cc	Wed Dec 10 11:55:21 2008 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,243 +0,0 @@
-/*
-
-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: ***
-*/
--- a/src/Makefile.in	Wed Dec 10 11:55:21 2008 +0100
+++ b/src/Makefile.in	Wed Dec 10 11:55:23 2008 +0100
@@ -65,7 +65,7 @@
 DLD_XSRC := amd.cc balance.cc besselj.cc betainc.cc bsxfun.cc cellfun.cc \
 	chol.cc ccolamd.cc colamd.cc colloc.cc conv2.cc convhulln.cc daspk.cc \
 	dasrt.cc dassl.cc det.cc dispatch.cc dlmread.cc dmperm.cc eig.cc \
-	expm.cc fft.cc fft2.cc fftn.cc fftw.cc filter.cc find.cc \
+	fft.cc fft2.cc fftn.cc fftw.cc filter.cc find.cc \
 	fltk_backend.cc \
 	gammainc.cc gcd.cc getgrent.cc getpwent.cc getrusage.cc \
 	givens.cc hess.cc hex2num.cc inv.cc kron.cc lookup.cc lsode.cc \