new file mode 100644
--- /dev/null
+++ b/src/sparse-xdiv.cc
@@ -0,0 +1,640 @@
+/*
+
+Copyright (C) 2004 David Bateman
+Copyright (C) 1998-2004 Andy Adler
+
+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 this program; see the file COPYING.  If not, write to the Free
+Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <cassert>
+
+#include "Array-util.h"
+#include "oct-cmplx.h"
+#include "quit.h"
+#include "error.h"
+
+#include "dSparse.h"
+#include "CSparse.h"
+#include "oct-spparms.h"
+#include "sparse-xdiv.h"
+
+static inline bool
+result_ok (int info)
+{
+#ifdef HAVE_LSSOLVE
+  return (info != -2 && info != -1);
+#else
+  // If the matrix is singular, who cares as we don't have QR based solver yet
+  return (info != -1);
+#endif
+}
+
+static void
+solve_singularity_warning (double rcond)
+{
+  warning ("matrix singular to machine precision, rcond = %g", rcond);
+  warning ("attempting to find minimum norm solution");
+}
+
+template <class T1, class T2>
+bool
+mx_leftdiv_conform (const T1& a, const T2& b)
+{
+  int a_nr = a.rows ();
+  int b_nr = b.rows ();
+
+  if (a_nr != b_nr)
+    {
+      int a_nc = a.cols ();
+      int b_nc = b.cols ();
+
+      gripe_nonconformant ("operator \\", a_nr, a_nc, b_nr, b_nc);
+      return false;
+    }
+
+  return true;
+}
+
+#define INSTANTIATE_MX_LEFTDIV_CONFORM(T1, T2) \
+  template bool mx_leftdiv_conform (const T1&, const T2&)
+
+INSTANTIATE_MX_LEFTDIV_CONFORM (SparseMatrix, SparseMatrix);
+INSTANTIATE_MX_LEFTDIV_CONFORM (SparseMatrix, SparseComplexMatrix);
+INSTANTIATE_MX_LEFTDIV_CONFORM (SparseComplexMatrix, SparseMatrix);
+INSTANTIATE_MX_LEFTDIV_CONFORM (SparseComplexMatrix, SparseComplexMatrix);
+INSTANTIATE_MX_LEFTDIV_CONFORM (SparseMatrix, Matrix);
+INSTANTIATE_MX_LEFTDIV_CONFORM (SparseMatrix, ComplexMatrix);
+INSTANTIATE_MX_LEFTDIV_CONFORM (SparseComplexMatrix, Matrix);
+INSTANTIATE_MX_LEFTDIV_CONFORM (SparseComplexMatrix, ComplexMatrix);
+
+template <class T1, class T2>
+bool
+mx_div_conform (const T1& a, const T2& b)
+{
+  int a_nc = a.cols ();
+  int b_nc = b.cols ();
+
+  if (a_nc != b_nc)
+    {
+      int a_nr = a.rows ();
+      int b_nr = b.rows ();
+
+      gripe_nonconformant ("operator /", a_nr, a_nc, b_nr, b_nc);
+      return false;
+    }
+
+  return true;
+}
+
+#define INSTANTIATE_MX_DIV_CONFORM(T1, T2) \
+  template bool mx_div_conform (const T1&, const T2&)
+
+INSTANTIATE_MX_DIV_CONFORM (SparseMatrix, SparseMatrix);
+INSTANTIATE_MX_DIV_CONFORM (SparseMatrix, SparseComplexMatrix);
+INSTANTIATE_MX_DIV_CONFORM (SparseComplexMatrix, SparseMatrix);
+INSTANTIATE_MX_DIV_CONFORM (SparseComplexMatrix, SparseComplexMatrix);
+INSTANTIATE_MX_DIV_CONFORM (Matrix, SparseMatrix);
+INSTANTIATE_MX_DIV_CONFORM (Matrix, SparseComplexMatrix);
+INSTANTIATE_MX_DIV_CONFORM (ComplexMatrix, SparseMatrix);
+INSTANTIATE_MX_DIV_CONFORM (ComplexMatrix, SparseComplexMatrix);
+
+// Right division functions.
+//
+//              op2 / op1:   m   cm   sm  scm
+//                   +--   +---+----+----+----+
+//   sparse matrix         | 1 |  3 |  5 |  7 |
+//                         +---+----+----+----+
+//   sparse complex_matrix | 2 |  4 |  6 |  8 |
+//                         +---+----+----+----+
+
+// -*- 1 -*-
+Matrix
+xdiv (const Matrix& a, const SparseMatrix& b)
+{
+  if (! mx_div_conform (a, b))
+    return Matrix ();
+
+  Matrix atmp = a.transpose ();
+  SparseMatrix btmp = b.transpose ();
+
+  int info;
+  if (btmp.rows () == btmp.columns ())
+    {
+      double rcond = 0.0;
+
+      Matrix result = btmp.solve (atmp, info, rcond, 
+				  solve_singularity_warning);
+
+      if (result_ok (info))
+	return Matrix (result.transpose ());
+    }
+
+  int rank;
+  Matrix result = btmp.lssolve (atmp, info, rank);
+
+  return result.transpose ();
+}
+
+// -*- 2 -*-
+ComplexMatrix
+xdiv (const Matrix& a, const SparseComplexMatrix& b)
+{
+  if (! mx_div_conform (a, b))
+    return ComplexMatrix ();
+
+  Matrix atmp = a.transpose ();
+  SparseComplexMatrix btmp = b.hermitian ();
+
+  int info;
+  if (btmp.rows () == btmp.columns ())
+    {
+      double rcond = 0.0;
+
+      ComplexMatrix result
+	= btmp.solve (atmp, info, rcond, solve_singularity_warning);
+
+      if (result_ok (info))
+	return result.hermitian ();
+    }
+
+  int rank;
+  ComplexMatrix result = btmp.lssolve (atmp, info, rank);
+
+  return result.hermitian ();
+}
+
+// -*- 3 -*-
+ComplexMatrix
+xdiv (const ComplexMatrix& a, const SparseMatrix& b)
+{
+  if (! mx_div_conform (a, b))
+    return ComplexMatrix ();
+
+  ComplexMatrix atmp = a.hermitian ();
+  SparseMatrix btmp = b.transpose ();
+
+  int info;
+  if (btmp.rows () == btmp.columns ())
+    {
+      double rcond = 0.0;
+
+      ComplexMatrix result
+	= btmp.solve (atmp, info, rcond, solve_singularity_warning);
+
+      if (result_ok (info))
+	return result.hermitian ();
+    }
+
+  int rank;
+  ComplexMatrix result = btmp.lssolve (atmp, info, rank);
+
+  return result.hermitian ();
+}
+
+// -*- 4 -*-
+ComplexMatrix
+xdiv (const ComplexMatrix& a, const SparseComplexMatrix& b)
+{
+  if (! mx_div_conform (a, b))
+    return ComplexMatrix ();
+
+  ComplexMatrix atmp = a.hermitian ();
+  SparseComplexMatrix btmp = b.hermitian ();
+
+  int info;
+  if (btmp.rows () == btmp.columns ())
+    {
+      double rcond = 0.0;
+
+      ComplexMatrix result
+	= btmp.solve (atmp, info, rcond, solve_singularity_warning);
+
+      if (result_ok (info))
+	return result.hermitian ();
+    }
+
+  int rank;
+  ComplexMatrix result = btmp.lssolve (atmp, info, rank);
+
+  return result.hermitian ();
+}
+
+// -*- 5 -*-
+SparseMatrix
+xdiv (const SparseMatrix& a, const SparseMatrix& b)
+{
+  if (! mx_div_conform (a, b))
+    return SparseMatrix ();
+
+  SparseMatrix atmp = a.transpose ();
+  SparseMatrix btmp = b.transpose ();
+
+  int info;
+  if (btmp.rows () == btmp.columns ())
+    {
+      double rcond = 0.0;
+
+      SparseMatrix result = btmp.solve (atmp, info, rcond, 
+					solve_singularity_warning);
+
+      if (result_ok (info))
+	return SparseMatrix (result.transpose ());
+    }
+
+  int rank;
+  SparseMatrix result = btmp.lssolve (atmp, info, rank);
+
+  return result.transpose ();
+}
+
+// -*- 6 -*-
+SparseComplexMatrix
+xdiv (const SparseMatrix& a, const SparseComplexMatrix& b)
+{
+  if (! mx_div_conform (a, b))
+    return SparseComplexMatrix ();
+
+  SparseMatrix atmp = a.transpose ();
+  SparseComplexMatrix btmp = b.hermitian ();
+
+  int info;
+  if (btmp.rows () == btmp.columns ())
+    {
+      double rcond = 0.0;
+
+      SparseComplexMatrix result
+	= btmp.solve (atmp, info, rcond, solve_singularity_warning);
+
+      if (result_ok (info))
+	return result.hermitian ();
+    }
+
+  int rank;
+  SparseComplexMatrix result = btmp.lssolve (atmp, info, rank);
+
+  return result.hermitian ();
+}
+
+// -*- 7 -*-
+SparseComplexMatrix
+xdiv (const SparseComplexMatrix& a, const SparseMatrix& b)
+{
+  if (! mx_div_conform (a, b))
+    return SparseComplexMatrix ();
+
+  SparseComplexMatrix atmp = a.hermitian ();
+  SparseMatrix btmp = b.transpose ();
+
+  int info;
+  if (btmp.rows () == btmp.columns ())
+    {
+      double rcond = 0.0;
+
+      SparseComplexMatrix result
+	= btmp.solve (atmp, info, rcond, solve_singularity_warning);
+
+      if (result_ok (info))
+	return result.hermitian ();
+    }
+
+  int rank;
+  SparseComplexMatrix result = btmp.lssolve (atmp, info, rank);
+
+  return result.hermitian ();
+}
+
+// -*- 8 -*-
+SparseComplexMatrix
+xdiv (const SparseComplexMatrix& a, const SparseComplexMatrix& b)
+{
+  if (! mx_div_conform (a, b))
+    return SparseComplexMatrix ();
+
+  SparseComplexMatrix atmp = a.hermitian ();
+  SparseComplexMatrix btmp = b.hermitian ();
+
+  int info;
+  if (btmp.rows () == btmp.columns ())
+    {
+      double rcond = 0.0;
+
+      SparseComplexMatrix result
+	= btmp.solve (atmp, info, rcond, solve_singularity_warning);
+
+      if (result_ok (info))
+	return result.hermitian ();
+    }
+
+  int rank;
+  SparseComplexMatrix result = btmp.lssolve (atmp, info, rank);
+
+  return result.hermitian ();
+}
+
+// Funny element by element division operations.
+//
+//       op2 \ op1:   s   cs
+//            +--   +---+----+
+//   matrix         | 1 |  3 |
+//                  +---+----+
+//   complex_matrix | 2 |  4 |
+//                  +---+----+
+
+Matrix
+x_el_div (double a, const SparseMatrix& b)
+{
+  int nr = b.rows ();
+  int nc = b.cols ();
+
+  Matrix result;
+  if (a == 0.)
+    result = Matrix (nr, nc, octave_NaN);
+  else if (a > 0.)
+    result = Matrix (nr, nc, octave_Inf);
+  else
+    result = Matrix (nr, nc, -octave_Inf);
+
+
+  for (int j = 0; j < nc; j++)
+    for (int i = b.cidx(j); i < b.cidx(j+1); i++)
+      {
+	OCTAVE_QUIT;
+	result.elem (b.ridx(i), j) = a / b.data (i);
+      }
+
+  return result;
+}
+
+ComplexMatrix
+x_el_div (double a, const SparseComplexMatrix& b)
+{
+  int nr = b.rows ();
+  int nc = b.cols ();
+
+  ComplexMatrix  result (nr, nc, Complex(octave_NaN, octave_NaN));
+
+  for (int j = 0; j < nc; j++)
+    for (int i = b.cidx(j); i < b.cidx(j+1); i++)
+      {
+	OCTAVE_QUIT;
+	result.elem (b.ridx(i), j) = a / b.data (i);
+      }
+
+  return result;
+}
+
+ComplexMatrix
+x_el_div (const Complex a, const SparseMatrix& b)
+{
+  int nr = b.rows ();
+  int nc = b.cols ();
+
+  ComplexMatrix result (nr, nc, (a / 0.0));
+
+  for (int j = 0; j < nc; j++)
+    for (int i = b.cidx(j); i < b.cidx(j+1); i++)
+      {
+	OCTAVE_QUIT;
+	result.elem (b.ridx(i), j) = a / b.data (i);
+      }
+
+  return result;
+}
+
+ComplexMatrix
+x_el_div (const Complex a, const SparseComplexMatrix& b)
+{
+  int nr = b.rows ();
+  int nc = b.cols ();
+
+  ComplexMatrix result (nr, nc, (a / 0.0));
+
+  for (int j = 0; j < nc; j++)
+    for (int i = b.cidx(j); i < b.cidx(j+1); i++)
+      {
+	OCTAVE_QUIT;
+	result.elem (b.ridx(i), j) = a / b.data (i);
+      }
+
+  return result;
+}
+
+// Left division functions.
+//
+//              op2 \ op1:   m   cm
+//                   +--   +---+----+
+//   matrix                | 1 |  5 |
+//                         +---+----+
+//   complex_matrix        | 2 |  6 |
+//                         +---+----+
+//   sparse matrix         | 3 |  7 |
+//                         +---+----+
+//   sparse complex_matrix | 4 |  8 |
+//                         +---+----+
+
+// -*- 1 -*-
+Matrix
+xleftdiv (const SparseMatrix& a, const Matrix& b)
+{
+  if (! mx_leftdiv_conform (a, b))
+    return Matrix ();
+
+  int info;
+  if (a.rows () == a.columns ())
+    {
+      double rcond = 0.0;
+
+      Matrix result
+	= a.solve (b, info, rcond, solve_singularity_warning);
+
+      if (result_ok (info))
+	return result;
+    }
+
+  int rank;
+  return a.lssolve (b, info, rank);
+}
+
+// -*- 2 -*-
+ComplexMatrix
+xleftdiv (const SparseMatrix& a, const ComplexMatrix& b)
+{
+  if (! mx_leftdiv_conform (a, b))
+    return ComplexMatrix ();
+
+  int info;
+  if (a.rows () == a.columns ())
+    {
+      double rcond = 0.0;
+
+      ComplexMatrix result
+	= a.solve (b, info, rcond, solve_singularity_warning);
+
+      if (result_ok (info))
+	return result;
+    }
+
+  int rank;
+  return a.lssolve (b, info, rank);
+}
+
+// -*- 3 -*-
+SparseMatrix
+xleftdiv (const SparseMatrix& a, const SparseMatrix& b)
+{
+  if (! mx_leftdiv_conform (a, b))
+    return SparseMatrix ();
+
+  int info;
+  if (a.rows () == a.columns ())
+    {
+      double rcond = 0.0;
+
+      SparseMatrix result
+	= a.solve (b, info, rcond, solve_singularity_warning);
+
+      if (result_ok (info))
+	return result;
+    }
+
+  int rank;
+  return a.lssolve (b, info, rank);
+}
+
+// -*- 4 -*-
+SparseComplexMatrix
+xleftdiv (const SparseMatrix& a, const SparseComplexMatrix& b)
+{
+  if (! mx_leftdiv_conform (a, b))
+    return SparseComplexMatrix ();
+
+  int info;
+  if (a.rows () == a.columns ())
+    {
+      double rcond = 0.0;
+
+      SparseComplexMatrix result
+	= a.solve (b, info, rcond, solve_singularity_warning);
+
+      if (result_ok (info))
+	return result;
+    }
+
+  int rank;
+  return a.lssolve (b, info, rank);
+}
+
+// -*- 5 -*-
+ComplexMatrix
+xleftdiv (const SparseComplexMatrix& a, const Matrix& b)
+{
+  if (! mx_leftdiv_conform (a, b))
+    return ComplexMatrix ();
+
+  int info;
+  if (a.rows () == a.columns ())
+    {
+      double rcond = 0.0;
+
+      ComplexMatrix result
+	= a.solve (b, info, rcond, solve_singularity_warning);
+
+      if (result_ok (info))
+	return result;
+    }
+
+  int rank;
+  return a.lssolve (b, info, rank);
+}
+
+// -*- 6 -*-
+ComplexMatrix
+xleftdiv (const SparseComplexMatrix& a, const ComplexMatrix& b)
+{
+  if (! mx_leftdiv_conform (a, b))
+    return ComplexMatrix ();
+
+  int info;
+  if (a.rows () == a.columns ())
+    {
+      double rcond = 0.0;
+
+      ComplexMatrix result
+	= a.solve (b, info, rcond, solve_singularity_warning);
+
+      if (result_ok (info))
+	return result;
+    }
+
+  int rank;
+  return a.lssolve (b, info, rank);
+}
+
+// -*- 7 -*-
+SparseComplexMatrix
+xleftdiv (const SparseComplexMatrix& a, const SparseMatrix& b)
+{
+  if (! mx_leftdiv_conform (a, b))
+    return SparseComplexMatrix ();
+
+  int info;
+  if (a.rows () == a.columns ())
+    {
+      double rcond = 0.0;
+
+      SparseComplexMatrix result
+	= a.solve (b, info, rcond, solve_singularity_warning);
+
+      if (result_ok (info))
+	return result;
+    }
+
+  int rank;
+  return a.lssolve (b, info, rank);
+}
+
+// -*- 8 -*-
+SparseComplexMatrix
+xleftdiv (const SparseComplexMatrix& a, const SparseComplexMatrix& b)
+{
+  if (! mx_leftdiv_conform (a, b))
+    return SparseComplexMatrix ();
+
+  int info;
+  if (a.rows () == a.columns ())
+    {
+      double rcond = 0.0;
+
+      SparseComplexMatrix result
+	= a.solve (b, info, rcond, solve_singularity_warning);
+
+      if (result_ok (info))
+	return result;
+    }
+
+  int rank;
+  return a.lssolve (b, info, rank);
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/