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