--- a/liboctave/sparse-dmsolve.cc	Thu Feb 11 12:16:43 2010 -0500
+++ b/liboctave/sparse-dmsolve.cc	Thu Feb 11 12:23:32 2010 -0500
@@ -37,10 +37,10 @@
 template <class T>
 static MSparse<T>
 dmsolve_extract (const MSparse<T> &A, const octave_idx_type *Pinv, 
-		const octave_idx_type *Q, octave_idx_type rst, 
-		octave_idx_type rend, octave_idx_type cst, 
-		octave_idx_type cend, octave_idx_type maxnz = -1,
-		bool lazy = false)
+                const octave_idx_type *Q, octave_idx_type rst, 
+                octave_idx_type rend, octave_idx_type cst, 
+                octave_idx_type cend, octave_idx_type maxnz = -1,
+                bool lazy = false)
 {
   octave_idx_type nz = (rend - rst) * (cend - cst);
   maxnz = (maxnz < 0 ? A.nnz () : maxnz);
@@ -53,20 +53,20 @@
     {
       nz = 0;
       for (octave_idx_type j = cst ; j < cend ; j++)
-	{
-	  octave_idx_type qq = (Q ? Q [j] : j);
-	  B.xcidx (j - cst) = nz;
-	  for (octave_idx_type p = A.cidx(qq) ; p < A.cidx (qq+1) ; p++)
-	    {
-	      octave_quit ();
-	      octave_idx_type r = (Pinv ? Pinv [A.ridx (p)] : A.ridx (p));
-	      if (r >= rst && r < rend)
-		{
-		  B.xdata (nz) = A.data (p);
-		  B.xridx (nz++) =  r - rst ;
-		}
-	    }
-	}
+        {
+          octave_idx_type qq = (Q ? Q [j] : j);
+          B.xcidx (j - cst) = nz;
+          for (octave_idx_type p = A.cidx(qq) ; p < A.cidx (qq+1) ; p++)
+            {
+              octave_quit ();
+              octave_idx_type r = (Pinv ? Pinv [A.ridx (p)] : A.ridx (p));
+              if (r >= rst && r < rend)
+                {
+                  B.xdata (nz) = A.data (p);
+                  B.xridx (nz++) =  r - rst ;
+                }
+            }
+        }
       B.xcidx (cend - cst) = nz ;
     }
   else
@@ -76,23 +76,23 @@
       octave_idx_type *ri = B.xridx();
       nz = 0;
       for (octave_idx_type j = cst ; j < cend ; j++)
-	{
-	  octave_idx_type qq = (Q ? Q [j] : j);
-	  B.xcidx (j - cst) = nz;
-	  for (octave_idx_type p = A.cidx(qq) ; p < A.cidx (qq+1) ; p++)
-	    {
-	      octave_quit ();
-	      octave_idx_type r = (Pinv ? Pinv [A.ridx (p)] : A.ridx (p));
-	      if (r >= rst && r < rend)
-		{
-		  X [r-rst] = A.data (p);
-		  B.xridx (nz++) =  r - rst ;
-		}
-	    }
-	  sort.sort (ri + B.xcidx (j - cst), nz - B.xcidx (j - cst));
-	  for (octave_idx_type p = B.cidx (j - cst); p < nz; p++)
-	    B.xdata (p) = X [B.xridx (p)]; 
-	}
+        {
+          octave_idx_type qq = (Q ? Q [j] : j);
+          B.xcidx (j - cst) = nz;
+          for (octave_idx_type p = A.cidx(qq) ; p < A.cidx (qq+1) ; p++)
+            {
+              octave_quit ();
+              octave_idx_type r = (Pinv ? Pinv [A.ridx (p)] : A.ridx (p));
+              if (r >= rst && r < rend)
+                {
+                  X [r-rst] = A.data (p);
+                  B.xridx (nz++) =  r - rst ;
+                }
+            }
+          sort.sort (ri + B.xcidx (j - cst), nz - B.xcidx (j - cst));
+          for (octave_idx_type p = B.cidx (j - cst); p < nz; p++)
+            B.xdata (p) = X [B.xridx (p)]; 
+        }
       B.xcidx (cend - cst) = nz ;
     }
 
@@ -102,25 +102,25 @@
 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
 static MSparse<double>
 dmsolve_extract (const MSparse<double> &A, const octave_idx_type *Pinv, 
-		const octave_idx_type *Q, octave_idx_type rst, 
-		octave_idx_type rend, octave_idx_type cst, 
-		octave_idx_type cend, octave_idx_type maxnz,
-		bool lazy);
+                const octave_idx_type *Q, octave_idx_type rst, 
+                octave_idx_type rend, octave_idx_type cst, 
+                octave_idx_type cend, octave_idx_type maxnz,
+                bool lazy);
 
 static MSparse<Complex>
 dmsolve_extract (const MSparse<Complex> &A, const octave_idx_type *Pinv, 
-		const octave_idx_type *Q, octave_idx_type rst, 
-		octave_idx_type rend, octave_idx_type cst, 
-		octave_idx_type cend, octave_idx_type maxnz,
-		bool lazy);
+                const octave_idx_type *Q, octave_idx_type rst, 
+                octave_idx_type rend, octave_idx_type cst, 
+                octave_idx_type cend, octave_idx_type maxnz,
+                bool lazy);
 #endif
 
 template <class T>
 static MArray2<T>
 dmsolve_extract (const MArray2<T> &m, const octave_idx_type *, 
-		 const octave_idx_type *, octave_idx_type r1, 
-		 octave_idx_type r2, octave_idx_type c1, 
-		 octave_idx_type c2)
+                 const octave_idx_type *, octave_idx_type r1, 
+                 octave_idx_type r2, octave_idx_type c1, 
+                 octave_idx_type c2)
 {
   r2 -= 1;
   c2 -= 1;
@@ -142,21 +142,21 @@
 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
 static MArray2<double>
 dmsolve_extract (const MArray2<double> &m, const octave_idx_type *, 
-		 const octave_idx_type *, octave_idx_type r1, 
-		 octave_idx_type r2, octave_idx_type c1, 
-		 octave_idx_type c2)
+                 const octave_idx_type *, octave_idx_type r1, 
+                 octave_idx_type r2, octave_idx_type c1, 
+                 octave_idx_type c2)
 
 static MArray2<Complex>
 dmsolve_extract (const MArray2<Complex> &m, const octave_idx_type *, 
-		 const octave_idx_type *, octave_idx_type r1, 
-		 octave_idx_type r2, octave_idx_type c1, 
-		 octave_idx_type c2)
+                 const octave_idx_type *, octave_idx_type r1, 
+                 octave_idx_type r2, octave_idx_type c1, 
+                 octave_idx_type c2)
 #endif
 
 template <class T>
 static void
 dmsolve_insert (MArray2<T> &a, const MArray2<T> &b, const octave_idx_type *Q,
-	       octave_idx_type r, octave_idx_type c)
+               octave_idx_type r, octave_idx_type c)
 {
   T *ax = a.fortran_vec();
   const T *bx = b.fortran_vec();
@@ -168,27 +168,27 @@
       octave_idx_type aoff = (c + j) * anr;
       octave_idx_type boff = j * nr;
       for (octave_idx_type i = 0; i < nr; i++)
-	{
-	  octave_quit ();
-	  ax [Q [r + i] + aoff] = bx [i + boff];
-	}
+        {
+          octave_quit ();
+          ax [Q [r + i] + aoff] = bx [i + boff];
+        }
     }
 }
 
 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
 static void
 dmsolve_insert (MArray2<double> &a, const MArray2<double> &b, 
-	       const octave_idx_type *Q, octave_idx_type r, octave_idx_type c);
+               const octave_idx_type *Q, octave_idx_type r, octave_idx_type c);
 
 static void
 dmsolve_insert (MArray2<Complex> &a, const MArray2<Complex> &b,
-	       const octave_idx_type *Q, octave_idx_type r, octave_idx_type c);
+               const octave_idx_type *Q, octave_idx_type r, octave_idx_type c);
 #endif
 
 template <class T>
 static void
 dmsolve_insert (MSparse<T> &a, const MSparse<T> &b, const octave_idx_type *Q,
-	       octave_idx_type r, octave_idx_type c)
+               octave_idx_type r, octave_idx_type c)
 {
   octave_idx_type b_rows = b.rows ();
   octave_idx_type b_cols = b.cols ();
@@ -208,7 +208,7 @@
   for (octave_idx_type i = c; i < c + b_cols; i++)
     for (octave_idx_type j = a.xcidx(i); j < a.xcidx(i+1); j++)
       if (Qinv [a.xridx(j)] < r || Qinv [a.xridx(j)] >= r + b_rows)
-	nel++;
+        nel++;
 
   OCTAVE_LOCAL_BUFFER (T, X, nr);
   octave_sort<octave_idx_type> sort;
@@ -231,33 +231,33 @@
       octave_quit ();
 
       for (octave_idx_type j = tmp.xcidx(i); j < tmp.xcidx(i+1); j++)
-	if (Qinv [tmp.xridx(j)] < r || 	Qinv [tmp.xridx(j)] >= r + b_rows)
-	  {
-	    X [tmp.xridx(j)] = tmp.xdata(j);
-	    a.xridx(ii++) = tmp.xridx(j);
-	  }
+        if (Qinv [tmp.xridx(j)] < r ||  Qinv [tmp.xridx(j)] >= r + b_rows)
+          {
+            X [tmp.xridx(j)] = tmp.xdata(j);
+            a.xridx(ii++) = tmp.xridx(j);
+          }
 
       octave_quit ();
 
       for (octave_idx_type j = b.cidx(i-c); j < b.cidx(i-c+1); j++)
-	{
-	  X [Q [r + b.ridx(j)]] = b.data(j);
-	  a.xridx(ii++) = Q [r + b.ridx(j)];
-	}
+        {
+          X [Q [r + b.ridx(j)]] = b.data(j);
+          a.xridx(ii++) = Q [r + b.ridx(j)];
+        }
 
       sort.sort (ri + a.xcidx (i), ii - a.xcidx (i));
       for (octave_idx_type p = a.xcidx (i); p < ii; p++)
-	a.xdata (p) = X [a.xridx (p)]; 
+        a.xdata (p) = X [a.xridx (p)]; 
       a.xcidx(i+1) = ii;
     }
 
   for (octave_idx_type i = c + b_cols; i < nc; i++)
     {
       for (octave_idx_type j = tmp.xcidx(i); j < tmp.cidx(i+1); j++)
-	{
-	  a.xdata(ii) = tmp.xdata(j);
-	  a.xridx(ii++) = tmp.xridx(j);
-	}
+        {
+          a.xdata(ii) = tmp.xdata(j);
+          a.xridx(ii++) = tmp.xridx(j);
+        }
       a.xcidx(i+1) = ii;
     }
 }
@@ -265,11 +265,11 @@
 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
 static void
 dmsolve_insert (MSparse<double> &a, const SparseMatrix &b, 
-	       const octave_idx_type *Q, octave_idx_type r, octave_idx_type c);
+               const octave_idx_type *Q, octave_idx_type r, octave_idx_type c);
 
 static void
 dmsolve_insert (MSparse<Complex> &a, const MSparse<Complex> &b,
-	       const octave_idx_type *Q, octave_idx_type r, octave_idx_type c);
+               const octave_idx_type *Q, octave_idx_type r, octave_idx_type c);
 #endif
 
 template <class T, class RT>
@@ -285,25 +285,25 @@
     {
       octave_idx_type off = j * b_nr;
       for (octave_idx_type i = 0; i < b_nr; i++)
-	{
-	  octave_quit ();
-	  Btx [p [i] + off] = Bx [ i + off];
-	}
+        {
+          octave_quit ();
+          Btx [p [i] + off] = Bx [ i + off];
+        }
     }
 }
 
 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
 static void
 dmsolve_permute (MArray2<double> &a, const MArray2<double>& b,
-		 const octave_idx_type *p);
+                 const octave_idx_type *p);
 
 static void
 dmsolve_permute (MArray2<Complex> &a, const MArray2<double>& b,
-		 const octave_idx_type *p);
+                 const octave_idx_type *p);
 
 static void
 dmsolve_permute (MArray2<Complex> &a, const MArray2<Complex>& b,
-		 const octave_idx_type *p);
+                 const octave_idx_type *p);
 #endif
 
 template <class T, class RT>
@@ -322,18 +322,18 @@
   for (octave_idx_type j = 0; j < b_nc; j++)
     {
       for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++)
-	{
-	  octave_quit ();
-	  octave_idx_type r = p [b.ridx (i)];
-	  X [r] = b.data (i);
-	  a.xridx(nz++) = p [b.ridx (i)];
-	}
+        {
+          octave_quit ();
+          octave_idx_type r = p [b.ridx (i)];
+          X [r] = b.data (i);
+          a.xridx(nz++) = p [b.ridx (i)];
+        }
       sort.sort (ri + a.xcidx (j), nz - a.xcidx (j));
       for (octave_idx_type i = a.cidx (j); i < nz; i++)
-	{
-	  octave_quit ();
-	  a.xdata (i) = X [a.xridx (i)]; 
-	}
+        {
+          octave_quit ();
+          a.xdata (i) = X [a.xridx (i)]; 
+        }
       a.xcidx(j+1) = nz;
     }
 }
@@ -341,15 +341,15 @@
 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
 static void
 dmsolve_permute (MSparse<double> &a, const MSparse<double>& b, 
-		 const octave_idx_type *p);
+                 const octave_idx_type *p);
 
 static void
 dmsolve_permute (MSparse<Complex> &a, const MSparse<double>& b,
-		 const octave_idx_type *p);
+                 const octave_idx_type *p);
 
 static void
 dmsolve_permute (MSparse<Complex> &a, const MSparse<Complex>& b,
-		 const octave_idx_type *p);
+                 const octave_idx_type *p);
 #endif
 
 static void
@@ -400,7 +400,7 @@
 #endif
       OCTAVE_LOCAL_BUFFER (octave_idx_type, pinv, nr);
       for (octave_idx_type i = 0; i < nr; i++)
-	pinv [p [i]] = i;
+        pinv [p [i]] = i;
       RT btmp;
       dmsolve_permute (btmp, b, pinv);
       info = 0;
@@ -408,66 +408,66 @@
 
       // Leading over-determined block
       if (dm->rr [2] < nr && dm->cc [3] < nc)
-	{
-	  ST m = dmsolve_extract (a, pinv, q, dm->rr [2], nr, dm->cc [3], nc, 
-				  nnz_remaining, true);
-	  nnz_remaining -= m.nnz();
-	  RT mtmp = 
-	    qrsolve (m, dmsolve_extract (btmp, 0, 0, dm->rr[2], b_nr, 0,
-					 b_nc), info);
-	  dmsolve_insert (retval, mtmp, q, dm->cc [3], 0);
-	  if (dm->rr [2] > 0 && !info)
-	    {
-	      m = dmsolve_extract (a, pinv, q, 0, dm->rr [2], 
-				   dm->cc [3], nc, nnz_remaining, true);
-	      nnz_remaining -= m.nnz();
-	      RT ctmp = dmsolve_extract (btmp, 0, 0, 0, 
-					 dm->rr[2], 0, b_nc);
-	      btmp.insert (ctmp - m * mtmp, 0, 0);
-	    }
-	}
+        {
+          ST m = dmsolve_extract (a, pinv, q, dm->rr [2], nr, dm->cc [3], nc, 
+                                  nnz_remaining, true);
+          nnz_remaining -= m.nnz();
+          RT mtmp = 
+            qrsolve (m, dmsolve_extract (btmp, 0, 0, dm->rr[2], b_nr, 0,
+                                         b_nc), info);
+          dmsolve_insert (retval, mtmp, q, dm->cc [3], 0);
+          if (dm->rr [2] > 0 && !info)
+            {
+              m = dmsolve_extract (a, pinv, q, 0, dm->rr [2], 
+                                   dm->cc [3], nc, nnz_remaining, true);
+              nnz_remaining -= m.nnz();
+              RT ctmp = dmsolve_extract (btmp, 0, 0, 0, 
+                                         dm->rr[2], 0, b_nc);
+              btmp.insert (ctmp - m * mtmp, 0, 0);
+            }
+        }
       
       // Structurally non-singular blocks
       // FIXME Should use fine Dulmange-Mendelsohn decomposition here.
       if (dm->rr [1] < dm->rr [2] && dm->cc [2] < dm->cc [3] && !info)
-	{
-	  ST m = dmsolve_extract (a, pinv, q, dm->rr [1], dm->rr [2], 
-				  dm->cc [2], dm->cc [3], nnz_remaining, false);
-	  nnz_remaining -= m.nnz();
-	  RT btmp2 = dmsolve_extract (btmp, 0, 0, dm->rr [1], dm->rr [2], 
-				      0, b_nc);
-	  double rcond = 0.0;
-	  MatrixType mtyp (MatrixType::Full);
-	  RT mtmp = m.solve (mtyp, btmp2, info, rcond, 
-			     solve_singularity_warning, false);	
-	  if (info != 0)
-	    {
-	      info = 0;
-	      mtmp = qrsolve (m, btmp2, info);
-	    }
+        {
+          ST m = dmsolve_extract (a, pinv, q, dm->rr [1], dm->rr [2], 
+                                  dm->cc [2], dm->cc [3], nnz_remaining, false);
+          nnz_remaining -= m.nnz();
+          RT btmp2 = dmsolve_extract (btmp, 0, 0, dm->rr [1], dm->rr [2], 
+                                      0, b_nc);
+          double rcond = 0.0;
+          MatrixType mtyp (MatrixType::Full);
+          RT mtmp = m.solve (mtyp, btmp2, info, rcond, 
+                             solve_singularity_warning, false); 
+          if (info != 0)
+            {
+              info = 0;
+              mtmp = qrsolve (m, btmp2, info);
+            }
 
-	  dmsolve_insert (retval, mtmp, q, dm->cc [2], 0);
-	  if (dm->rr [1] > 0 && !info)
-	    {
-	      m = dmsolve_extract (a, pinv, q, 0, dm->rr [1], dm->cc [2],
-				   dm->cc [3], nnz_remaining, true);
-	      nnz_remaining -= m.nnz();
-	      RT ctmp = dmsolve_extract (btmp, 0, 0, 0,
-					 dm->rr[1], 0, b_nc);
-	      btmp.insert (ctmp - m * mtmp, 0, 0);
-	    }
-	}
+          dmsolve_insert (retval, mtmp, q, dm->cc [2], 0);
+          if (dm->rr [1] > 0 && !info)
+            {
+              m = dmsolve_extract (a, pinv, q, 0, dm->rr [1], dm->cc [2],
+                                   dm->cc [3], nnz_remaining, true);
+              nnz_remaining -= m.nnz();
+              RT ctmp = dmsolve_extract (btmp, 0, 0, 0,
+                                         dm->rr[1], 0, b_nc);
+              btmp.insert (ctmp - m * mtmp, 0, 0);
+            }
+        }
 
       // Trailing under-determined block
       if (dm->rr [1] > 0 && dm->cc [2] > 0 && !info)
-	{
-	  ST m = dmsolve_extract (a, pinv, q, 0, dm->rr [1], 0, 
-				  dm->cc [2], nnz_remaining, true);
-	  RT mtmp = 
-	    qrsolve (m, dmsolve_extract(btmp, 0, 0, 0, dm->rr [1] , 0, 
-					b_nc), info);
-	  dmsolve_insert (retval, mtmp, q, 0, 0);
-	}
+        {
+          ST m = dmsolve_extract (a, pinv, q, 0, dm->rr [1], 0, 
+                                  dm->cc [2], nnz_remaining, true);
+          RT mtmp = 
+            qrsolve (m, dmsolve_extract(btmp, 0, 0, 0, dm->rr [1] , 0, 
+                                        b_nc), info);
+          dmsolve_insert (retval, mtmp, q, 0, 0);
+        }
 
       CXSPARSE_DNAME (_dfree) (dm);
     }
@@ -480,33 +480,33 @@
 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
 extern Matrix
 dmsolve (const SparseMatrix &a, const Matrix &b, 
-	 octave_idx_type &info);
+         octave_idx_type &info);
 
 extern ComplexMatrix
 dmsolve (const SparseMatrix &a, const ComplexMatrix &b, 
-	 octave_idx_type &info);
+         octave_idx_type &info);
 
 extern ComplexMatrix
 dmsolve (const SparseComplexMatrix &a, const Matrix &b, 
-	 octave_idx_type &info);
+         octave_idx_type &info);
 
 extern ComplexMatrix
 dmsolve (const SparseComplexMatrix &a, const ComplexMatrix &b, 
-	 octave_idx_type &info);
+         octave_idx_type &info);
 
 extern SparseMatrix
 dmsolve (const SparseMatrix &a, const SparseMatrix &b, 
-	 octave_idx_type &info);
+         octave_idx_type &info);
 
 extern SparseComplexMatrix
 dmsolve (const SparseMatrix &a, const SparseComplexMatrix &b, 
-	 octave_idx_type &info);
+         octave_idx_type &info);
 
 extern SparseComplexMatrix
 dmsolve (const SparseComplexMatrix &a, const SparseMatrix &b, 
-	 octave_idx_type &info);
+         octave_idx_type &info);
 
 extern SparseComplexMatrix
 dmsolve (const SparseComplexMatrix &a, const SparseComplexMatrix &b, 
-	 octave_idx_type &info);
+         octave_idx_type &info);
 #endif