Re: FYI: sparse indexed assignment rewritten

[David Bateman]

Ben Abbott wrote:
> On Friday, April 16, 2010, at 11:18AM, "David Bateman" <address@hidden> wrote:
>
>   
> > Ok, then consider the attached changeset which has the following properties
> >
> >     
>
> Wrong changeset I think.
>
> The one you attached is for the background color bug seen when displaying 
> images, correct?
>
> Ben
>
>
>   
Damn.. Changeset attached

D.

# HG changeset patch
# User David Bateman <address@hidden>
# Date 1271427285 -7200
# Node ID b4d2080b6df714365c77f734481ce2a663f74ca3
# Parent  660c244d3206ebbf9261416b202bd2bf15642502
Replace nzmax by nnz as needed

diff --git a/ChangeLog b/ChangeLog
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,7 @@
+2010-04-16  David Bateman  <address@hidden>
+
+       * PROJECTS: Update for new sparse functionality.
+
 2010-04-14  Shai Ayal  <address@hidden>
 
        * NEWS: Update.
diff --git a/PROJECTS b/PROJECTS
--- a/PROJECTS
+++ b/PROJECTS
@@ -61,41 +61,21 @@
   * Improve QR factorization functions, using idea based on CSPARSE
     cs_dmsol.m 
 
+  * Improve QR fqctorization by replace CXSPARSE code with SPQR code, and make
+    the linear solve return 2-norm solutions for ill-conditioned matrices
+    based on this new code
+
   * Implement fourth argument to the sprand and sprandn, and addition
     arguments to sprandsym that the leading brand implements.
 
   * Sparse logical indexing in idx_vector class so that something like
-    "a=sprandn(1e6,1e6,1e-6); a(a<1) = 0" won't cause a memory overflow.
-
-  * Make spalloc (r, c, n) actually create an empty sparse with n
-    non-zero elements?  This allows something like
-
-    sm = spalloc (r, c, n)
-    for j=1:c
-      for i=1:r
-        tmp = foo (i,j);
-        if (tmp != 0.)
-          sm (i,j) = tmp;
-        endif
-      endfor
-    endfor
-
-    actually make sense.  Otherwise the above will cause massive amounts
-    of memory reallocation.
-
-    The fact is that this doesn't make sense in any case as the assign
-    function makes another copy of the sparse matrix.  So although spalloc
-    might easily be made to have the correct behavior, the first assign
-    will cause the matrix to be resized!  There seems to be no simple
-    way to treat this but a complete rewrite of the sparse assignment
-    functions...
+    'a=sprandn(1e6,1e6,1e-6); a(a<1) = 0' won't cause a memory overflow.
 
   * Other missing Functions
       - symmmd      Superseded by symamd
       - colmmd      Superseded by colamd
       - cholinc
       - bicg        Can this be taken from octave-forge?
-      - cgs
       - gmres
       - lsqr
       - minres
diff --git a/liboctave/CSparse.h b/liboctave/CSparse.h
--- a/liboctave/CSparse.h
+++ b/liboctave/CSparse.h
@@ -80,8 +80,9 @@
 
   SparseComplexMatrix (const Array<Complex>& a, const idx_vector& r, 
                        const idx_vector& c, octave_idx_type nr = -1, 
-                       octave_idx_type nc = -1, bool sum_terms = true)
-    : MSparse<Complex> (a, r, c, nr, nc, sum_terms) { }
+                       octave_idx_type nc = -1, bool sum_terms = true,
+                       octave_idx_type nzm = -1)
+    : MSparse<Complex> (a, r, c, nr, nc, sum_terms, nzm) { }
 
   explicit SparseComplexMatrix (const SparseMatrix& a);
 
diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,17 @@
+2010-04-16  David Bateman  <address@hidden>
+
+       * Sparse.cc (template <class T> Sparse<T>::Sparse (const Array<T>&,
+       const idx_vector&, const idx_vector&, octave_idx_type,
+        octave_idx_type, bool, octave_idx_type)): Add argument defining the
+       minimum storage to allocate for the sparse matrix.
+       * Sparse.h (template <class T> Sparse (const Array<T>&,
+       const idx_vector&, const idx_vector&, octave_idx_type,
+        octave_idx_type, bool, octave_idx_type)): ditto.
+       * MSparse.h : ditto
+       * CSparse.h : ditto
+       * dSparse.h : ditto
+       * boolSparse.h : ditto
+
 2010-04-14  Jaroslav Hajek  <address@hidden>
 
        * Sparse.cc: Update failing tests.
diff --git a/liboctave/MSparse.h b/liboctave/MSparse.h
--- a/liboctave/MSparse.h
+++ b/liboctave/MSparse.h
@@ -59,8 +59,9 @@
   MSparse (const Sparse<U>& a) : Sparse<T> (a) { }
 
   MSparse (const Array<T>& a, const idx_vector& r, const idx_vector& c,
-           octave_idx_type nr = -1, octave_idx_type nc = -1, bool sum_terms = 
true)
-    : Sparse<T> (a, r, c, nr, nc, sum_terms) { }
+           octave_idx_type nr = -1, octave_idx_type nc = -1, 
+           bool sum_terms = true, octave_idx_type nzm = -1)
+    : Sparse<T> (a, r, c, nr, nc, sum_terms, nzm) { }
 
   explicit MSparse (octave_idx_type r, octave_idx_type c, T val) : Sparse<T> 
(r, c, val) { }
 
diff --git a/liboctave/Sparse.cc b/liboctave/Sparse.cc
--- a/liboctave/Sparse.cc
+++ b/liboctave/Sparse.cc
@@ -257,7 +257,8 @@
 template <class T>
 Sparse<T>::Sparse (const Array<T>& a, const idx_vector& r, 
                    const idx_vector& c, octave_idx_type nr,
-                   octave_idx_type nc, bool sum_terms)
+                   octave_idx_type nc, bool sum_terms,
+                   octave_idx_type nzm)
   : rep (nil_rep ()), dimensions ()
 {
   if (nr < 0)
@@ -296,7 +297,7 @@
     {
       if (n == 1 && a(0) != T ())
         {
-          change_capacity (1);
+          change_capacity (nzm > 1 ? nzm : 1);
           xridx(0) = r(0);
           xdata(0) = a(0);
           for (octave_idx_type j = 0; j < nc; j++)
@@ -324,7 +325,7 @@
           for (octave_idx_type i = 1; i < n; i++)
             new_nz += rd[i-1] != rd[i];
           // Allocate result.
-          change_capacity (new_nz);
+          change_capacity (nzm > new_nz ? nzm : new_nz);
           xcidx (1) = new_nz;
           octave_idx_type *rri = ridx ();
           T *rrd = data ();
@@ -407,7 +408,7 @@
               xcidx(j+1) = xcidx(j) + nzj;
             }
 
-          change_capacity (xcidx (nc));
+          change_capacity (nzm > xcidx (nc) ? nzm : xcidx (nc));
           octave_idx_type *rri = ridx ();
           T *rrd = data ();
 
@@ -463,7 +464,7 @@
       for (octave_idx_type i = 1; i < n; i++)
         new_nz += rd[i-1] != rd[i];
       // Allocate result.
-      change_capacity (new_nz);
+      change_capacity (nzm > new_nz ? nzm : new_nz);
       xcidx(1) = new_nz;
       octave_idx_type *rri = ridx ();
       T *rrd = data ();
@@ -550,7 +551,7 @@
           xcidx(j+1) = xcidx(j) + nzj;
         }
 
-      change_capacity (xcidx (nc));
+      change_capacity (nzm > xcidx (nc) ? nzm : xcidx (nc));
       octave_idx_type *rri = ridx ();
       T *rrd = data ();
 
diff --git a/liboctave/Sparse.h b/liboctave/Sparse.h
--- a/liboctave/Sparse.h
+++ b/liboctave/Sparse.h
@@ -221,7 +221,8 @@
   Sparse (const Sparse<T>& a, const dim_vector& dv);
 
   Sparse (const Array<T>& a, const idx_vector& r, const idx_vector& c,
-          octave_idx_type nr = -1, octave_idx_type nc = -1, bool sum_terms = 
true);
+          octave_idx_type nr = -1, octave_idx_type nc = -1,
+          bool sum_terms = true, octave_idx_type nzm = -1);
 
   // Sparsify a normal matrix
   Sparse (const Array<T>& a);
diff --git a/liboctave/boolSparse.h b/liboctave/boolSparse.h
--- a/liboctave/boolSparse.h
+++ b/liboctave/boolSparse.h
@@ -57,8 +57,9 @@
 
   SparseBoolMatrix (const Array<bool>& a, const idx_vector& r, 
                     const idx_vector& c, octave_idx_type nr = -1, 
-                    octave_idx_type nc = -1, bool sum_terms = true)
-    : Sparse<bool> (a, r, c, nr, nc, sum_terms) { }
+                    octave_idx_type nc = -1, bool sum_terms = true,
+                    octave_idx_type nzm = -1)
+    : Sparse<bool> (a, r, c, nr, nc, sum_terms, nzm) { }
 
   SparseBoolMatrix (octave_idx_type r, octave_idx_type c, octave_idx_type 
num_nz) : Sparse<bool> (r, c, num_nz) { }
 
diff --git a/liboctave/dSparse.h b/liboctave/dSparse.h
--- a/liboctave/dSparse.h
+++ b/liboctave/dSparse.h
@@ -74,8 +74,9 @@
 
   SparseMatrix (const Array<double>& a, const idx_vector& r, 
                 const idx_vector& c, octave_idx_type nr = -1, 
-                octave_idx_type nc = -1, bool sum_terms = true)
-    : MSparse<double> (a, r, c, nr, nc, sum_terms) { }
+                octave_idx_type nc = -1, bool sum_terms = true,
+                octave_idx_type nzm = -1)
+    : MSparse<double> (a, r, c, nr, nc, sum_terms, nzm) { }
 
   explicit SparseMatrix (const DiagMatrix& a);
 
diff --git a/src/DLD-FUNCTIONS/__glpk__.cc b/src/DLD-FUNCTIONS/__glpk__.cc
--- a/src/DLD-FUNCTIONS/__glpk__.cc
+++ b/src/DLD-FUNCTIONS/__glpk__.cc
@@ -519,7 +519,7 @@
 
       mrowsA = A.rows ();
       octave_idx_type Anc = A.cols ();
-      octave_idx_type Anz = A.nzmax ();
+      octave_idx_type Anz = A.nnz ();
       rn.resize (Anz+1, 1);
       cn.resize (Anz+1, 1);
       a.resize (Anz+1, 0.0);
diff --git a/src/DLD-FUNCTIONS/amd.cc b/src/DLD-FUNCTIONS/amd.cc
--- a/src/DLD-FUNCTIONS/amd.cc
+++ b/src/DLD-FUNCTIONS/amd.cc
@@ -91,7 +91,7 @@
     print_usage ();
   else
     {
-      octave_idx_type n_row, n_col, nnz;
+      octave_idx_type n_row, n_col;
       const octave_idx_type *ridx, *cidx;
       SparseMatrix sm;
       SparseComplexMatrix scm;
@@ -103,7 +103,6 @@
               scm = args(0).sparse_complex_matrix_value ();
               n_row = scm.rows ();
               n_col = scm.cols ();
-              nnz = scm.nzmax ();
               ridx = scm.xridx ();
               cidx = scm.xcidx ();
             }
@@ -112,7 +111,6 @@
               sm = args(0).sparse_matrix_value ();
               n_row = sm.rows ();
               n_col = sm.cols ();
-              nnz = sm.nzmax ();
               ridx = sm.xridx ();
               cidx = sm.xcidx ();
             }
@@ -126,7 +124,6 @@
           
           n_row = sm.rows ();
           n_col = sm.cols ();
-          nnz = sm.nzmax ();
           ridx = sm.xridx ();
           cidx = sm.xcidx ();
         }
diff --git a/src/DLD-FUNCTIONS/ccolamd.cc b/src/DLD-FUNCTIONS/ccolamd.cc
--- a/src/DLD-FUNCTIONS/ccolamd.cc
+++ b/src/DLD-FUNCTIONS/ccolamd.cc
@@ -223,7 +223,7 @@
               scm = args(0). sparse_complex_matrix_value ();
               n_row = scm.rows ();
               n_col = scm.cols ();
-              nnz = scm.nzmax ();
+              nnz = scm.nnz ();
               ridx = scm.xridx ();
               cidx = scm.xcidx ();
             }
@@ -233,7 +233,7 @@
 
               n_row = sm.rows ();
               n_col = sm.cols ();
-              nnz = sm.nzmax ();
+              nnz = sm.nnz ();
               ridx = sm.xridx ();
               cidx = sm.xcidx ();
             }
@@ -247,7 +247,7 @@
 
           n_row = sm.rows ();
           n_col = sm.cols ();
-          nnz = sm.nzmax ();
+          nnz = sm.nnz ();
           ridx = sm.xridx ();
           cidx = sm.xcidx ();
         }
@@ -461,7 +461,7 @@
               scm = args(0).sparse_complex_matrix_value ();
               n_row = scm.rows ();
               n_col = scm.cols ();
-              nnz = scm.nzmax ();
+              nnz = scm.nnz ();
               ridx = scm.xridx ();
               cidx = scm.xcidx ();
             }
@@ -470,7 +470,7 @@
               sm = args(0).sparse_matrix_value ();
               n_row = sm.rows ();
               n_col = sm.cols ();
-              nnz = sm.nzmax ();
+              nnz = sm.nnz ();
               ridx = sm.xridx ();
               cidx = sm.xcidx ();
             }
@@ -484,7 +484,7 @@
           
           n_row = sm.rows ();
           n_col = sm.cols ();
-          nnz = sm.nzmax ();
+          nnz = sm.nnz ();
           ridx = sm.xridx ();
           cidx = sm.xcidx ();
         }
diff --git a/src/DLD-FUNCTIONS/colamd.cc b/src/DLD-FUNCTIONS/colamd.cc
--- a/src/DLD-FUNCTIONS/colamd.cc
+++ b/src/DLD-FUNCTIONS/colamd.cc
@@ -347,7 +347,7 @@
               scm = args(0). sparse_complex_matrix_value ();
               n_row = scm.rows ();
               n_col = scm.cols ();
-              nnz = scm.nzmax ();
+              nnz = scm.nnz ();
               ridx = scm.xridx ();
               cidx = scm.xcidx ();
             }
@@ -357,7 +357,7 @@
 
               n_row = sm.rows ();
               n_col = sm.cols ();
-              nnz = sm.nzmax ();
+              nnz = sm.nnz ();
               ridx = sm.xridx ();
               cidx = sm.xcidx ();
             }
@@ -371,7 +371,7 @@
 
           n_row = sm.rows ();
           n_col = sm.cols ();
-          nnz = sm.nzmax ();
+          nnz = sm.nnz ();
           ridx = sm.xridx ();
           cidx = sm.xcidx ();
         }
@@ -556,7 +556,7 @@
               scm = args(0).sparse_complex_matrix_value ();
               n_row = scm.rows ();
               n_col = scm.cols ();
-              nnz = scm.nzmax ();
+              nnz = scm.nnz ();
               ridx = scm.xridx ();
               cidx = scm.xcidx ();
             }
@@ -565,7 +565,7 @@
               sm = args(0).sparse_matrix_value ();
               n_row = sm.rows ();
               n_col = sm.cols ();
-              nnz = sm.nzmax ();
+              nnz = sm.nnz ();
               ridx = sm.xridx ();
               cidx = sm.xcidx ();
             }
@@ -579,7 +579,7 @@
           
           n_row = sm.rows ();
           n_col = sm.cols ();
-          nnz = sm.nzmax ();
+          nnz = sm.nnz ();
           ridx = sm.xridx ();
           cidx = sm.xcidx ();
         }
@@ -681,7 +681,7 @@
               scm = args(0).sparse_complex_matrix_value ();
               n_row = scm.rows ();
               n_col = scm.cols ();
-              nnz = scm.nzmax ();
+              nnz = scm.nnz ();
               ridx = scm.xridx ();
               cidx = scm.xcidx ();
             }
@@ -690,7 +690,7 @@
               sm = args(0).sparse_matrix_value ();
               n_row = sm.rows ();
               n_col = sm.cols ();
-              nnz = sm.nzmax ();
+              nnz = sm.nnz ();
               ridx = sm.xridx ();
               cidx = sm.xcidx ();
             }
diff --git a/src/DLD-FUNCTIONS/kron.cc b/src/DLD-FUNCTIONS/kron.cc
--- a/src/DLD-FUNCTIONS/kron.cc
+++ b/src/DLD-FUNCTIONS/kron.cc
@@ -100,7 +100,7 @@
 {
   octave_idx_type idx = 0;
   MSparse<T> C (A.rows () * B.rows (), A.columns () * B.columns (), 
-                A.nzmax () * B.nzmax ());
+                A.nnz () * B.nnz ());
 
   C.cidx (0) = 0;
 
diff --git a/src/DLD-FUNCTIONS/sparse.cc b/src/DLD-FUNCTIONS/sparse.cc
--- a/src/DLD-FUNCTIONS/sparse.cc
+++ b/src/DLD-FUNCTIONS/sparse.cc
@@ -156,7 +156,13 @@
 
        if (! error_state)
          {
-           octave_idx_type m = -1, n = -1;
+           octave_idx_type m = -1, n = -1, nzmax = -1;
+           if (nargin == 6)
+             {
+               nzmax = args(5).idx_type_value ();
+               nargin --;
+             }
+
            if (nargin == 5)
              {
                if (args(3).is_scalar_type () && args(4).is_scalar_type ())
@@ -181,13 +187,13 @@
 
                if (args(2).is_bool_type ())
                  retval = SparseBoolMatrix (args(2).bool_array_value (), i, j,
-                                            m, n, summation);
+                                            m, n, summation, nzmax);
                else if (args(2).is_complex_type ())
-                 retval = SparseComplexMatrix (args(2).complex_array_value (), 
i, j,
-                                               m, n, summation);
+                 retval = SparseComplexMatrix (args(2).complex_array_value (),
+                                               i, j, m, n, summation, nzmax);
                else if (args(2).is_numeric_type ())
                  retval = SparseMatrix (args(2).array_value (), i, j,
-                                        m, n, summation);
+                                        m, n, summation, nzmax);
                else
                  gripe_wrong_type_arg ("sparse", args(2));
              }
diff --git a/src/ov-base-sparse.cc b/src/ov-base-sparse.cc
--- a/src/ov-base-sparse.cc
+++ b/src/ov-base-sparse.cc
@@ -361,7 +361,7 @@
   // Ensure that additional memory is deallocated
   matrix.maybe_compress ();
 
-  os << "# nnz: "      << nzmax () << "\n";
+  os << "# nnz: "      << nnz () << "\n";
   os << "# rows: "     << dv (0) << "\n";
   os << "# columns: "  << dv (1) << "\n";
 
diff --git a/src/ov-bool-sparse.cc b/src/ov-bool-sparse.cc
--- a/src/ov-bool-sparse.cc
+++ b/src/ov-bool-sparse.cc
@@ -217,7 +217,7 @@
 
   int nr = d(0);
   int nc = d(1);
-  int nz = nzmax ();
+  int nz = nnz ();
 
   int32_t itmp;
   // Use negative value for ndims to be consistent with other formats
@@ -427,7 +427,7 @@
       return false;
     }
   
-  tmp = m.nzmax ();
+  tmp = m.nnz ();
   retval = H5Dwrite (data_hid, H5T_NATIVE_IDX, H5S_ALL, H5S_ALL,
                      H5P_DEFAULT, &tmp) >= 0;
   H5Dclose (data_hid);
@@ -478,7 +478,7 @@
 
   H5Sclose (space_hid);
 
-  hdims[0] = m.nzmax ();
+  hdims[0] = m.nnz ();
   hdims[1] = 1;
 
   space_hid = H5Screate_simple (2, hdims, 0);
@@ -528,8 +528,8 @@
       return false;
     }
 
-  OCTAVE_LOCAL_BUFFER (hbool_t, htmp, m.nzmax ());  
-  for (int i = 0; i < m.nzmax (); i++)
+  OCTAVE_LOCAL_BUFFER (hbool_t, htmp, m.nnz ());  
+  for (int i = 0; i < m.nnz (); i++)
     htmp[i] = m.xdata(i);
 
   retval = H5Dwrite (data_hid, H5T_NATIVE_HBOOL, H5S_ALL, H5S_ALL,
diff --git a/src/ov-cx-sparse.cc b/src/ov-cx-sparse.cc
--- a/src/ov-cx-sparse.cc
+++ b/src/ov-cx-sparse.cc
@@ -236,7 +236,7 @@
 
   int nr = d(0);
   int nc = d(1);
-  int nz = nzmax ();
+  int nz = nnz ();
 
   int32_t itmp;
   // Use negative value for ndims to be consistent with other formats
@@ -263,7 +263,7 @@
       else
         st = LS_FLOAT;
     }
-  else if (matrix.nzmax () > 8192) // FIXME -- make this configurable.
+  else if (matrix.nnz () > 8192) // FIXME -- make this configurable.
     {
       double max_val, min_val;
       if (matrix.all_integers (max_val, min_val))
@@ -463,7 +463,7 @@
       return false;
     }
   
-  tmp = m.nzmax ();
+  tmp = m.nnz ();
   retval = H5Dwrite (data_hid, H5T_NATIVE_IDX, H5S_ALL, H5S_ALL,
                      H5P_DEFAULT, &tmp) >= 0;
   H5Dclose (data_hid);
@@ -514,7 +514,7 @@
 
   H5Sclose (space_hid);
 
-  hdims[0] = m.nzmax ();
+  hdims[0] = m.nnz ();
   hdims[1] = 1;
 
   space_hid = H5Screate_simple (2, hdims, 0);
diff --git a/src/ov-re-sparse.cc b/src/ov-re-sparse.cc
--- a/src/ov-re-sparse.cc
+++ b/src/ov-re-sparse.cc
@@ -269,7 +269,7 @@
 
   int nr = d(0);
   int nc = d(1);
-  int nz = nzmax ();
+  int nz = nnz ();
 
   int32_t itmp;
   // Use negative value for ndims to be consistent with other formats
@@ -296,7 +296,7 @@
       else
         st = LS_FLOAT;
     }
-  else if (matrix.nzmax () > 8192) // FIXME -- make this configurable.
+  else if (matrix.nnz () > 8192) // FIXME -- make this configurable.
     {
       double max_val, min_val;
       if (matrix.all_integers (max_val, min_val))
@@ -493,7 +493,7 @@
       return false;
     }
   
-  tmp = m.nzmax ();
+  tmp = m.nnz ();
   retval = H5Dwrite (data_hid, H5T_NATIVE_IDX, H5S_ALL, H5S_ALL, H5P_DEFAULT,
                      &tmp) >= 0;
   H5Dclose (data_hid);
@@ -544,7 +544,7 @@
 
   H5Sclose (space_hid);
 
-  hdims[0] = m.nzmax ();
+  hdims[0] = m.nnz ();
   hdims[1] = 1;
 
   space_hid = H5Screate_simple (2, hdims, 0);
diff --git a/src/sparse-xpow.cc b/src/sparse-xpow.cc
--- a/src/sparse-xpow.cc
+++ b/src/sparse-xpow.cc
@@ -307,7 +307,7 @@
 
   octave_value retval;
 
-  octave_idx_type nz = a.nzmax ();
+  octave_idx_type nz = a.nnz ();
 
   if (b <= 0.0)
     {
@@ -474,7 +474,7 @@
     retval = octave_value (NDArray (a.dims (), 1));
   else
     {
-      octave_idx_type nz = a.nzmax ();
+      octave_idx_type nz = a.nnz ();
       SparseComplexMatrix result (a);
       
       for (octave_idx_type i = 0; i < nz; i++)
@@ -602,7 +602,7 @@
     }
   else
     {
-      octave_idx_type nz = a.nzmax ();
+      octave_idx_type nz = a.nnz ();
 
       SparseComplexMatrix result (a);
   
@@ -681,7 +681,7 @@
   else
     {
 
-      octave_idx_type nz = a.nzmax ();
+      octave_idx_type nz = a.nnz ();
 
       SparseComplexMatrix result (a);