--- a/liboctave/fMatrix.cc
+++ b/liboctave/fMatrix.cc
@@ -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)
 {