Mercurial > hg > octave-max

changeset 5603:2c66c36d2698

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2006-01-31 11:57:47 by dbateman]
author dbateman
date Tue, 31 Jan 2006 11:57:47 +0000
parents dfa2da0563e5
children 2857357f9d3c
files configure.in liboctave/ChangeLog liboctave/Sparse.cc liboctave/sparse-sort.cc liboctave/sparse-sort.h src/ChangeLog src/DLD-FUNCTIONS/__glpk__.cc test/ChangeLog test/build_sparse_tests.sh
diffstat 9 files changed, 273 insertions(+), 40 deletions(-) [+]
line wrap: on
line diff
--- a/configure.in
+++ b/configure.in
@@ -29,7 +29,7 @@
 EXTERN_CXXFLAGS="$CXXFLAGS"
 
 AC_INIT
-AC_REVISION($Revision: 1.497 $)
+AC_REVISION($Revision: 1.498 $)
 AC_PREREQ(2.57)
 AC_CONFIG_SRCDIR([src/octave.cc])
 AC_CONFIG_HEADER(config.h)
@@ -841,7 +841,7 @@
 fi 
 
 if test -n "$warn_ccolamd"; then
-  AC_MSG_WARN($warn_colamd)
+  AC_MSG_WARN($warn_ccolamd)
 fi
 
 CHOLMOD_LIBS=
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,18 @@
+2006-01-21  David Bateman  <dbateman@free.fr>
+
+        * sparse-sort.cc (bool octave_sparse_sidxl_comp): 64-bit fix.
+        (bool octave_idx_vector_comp): New function.
+        (template class octave_sort<octave_idx_vector_sort *>): Instantiate
+        indexed idx_vector sorting function.
+        * sparse-sort.h (class octave_sparse_sort_idxl): 64-bit fix.
+        (class octave_idx_vector_sort): New class for indexed idx_vector
+        sorting.
+        (bool octave_idx_vector_comp): Declaration.
+        * Sparse.cc (int assign1(Sparse<LT>&, Sparse<RT>&)): Treat cases of
+        unordered LHS indexes in assignment using new octave_idx_vector_sort
+        class.
+        (int assign(Sparse<LT>&, Sparse<RT>&)): ditto.
+
 2006-01-30  John W. Eaton  <jwe@octave.org>
 
 	* so-array.h (streamoff_array::nnz): New funtion.
--- a/liboctave/Sparse.cc
+++ b/liboctave/Sparse.cc
@@ -1938,7 +1938,8 @@
   octave_idx_type nc = lhs.cols ();
   octave_idx_type nz = lhs.nnz ();
 
-  octave_idx_type n = lhs_idx.freeze (lhs_len, "vector", true, liboctave_wrore_flag);
+  octave_idx_type n = lhs_idx.freeze (lhs_len, "vector", true, 
+				      liboctave_wrore_flag);
 
   if (n != 0)
     {
@@ -1953,6 +1954,43 @@
 	{
 	  octave_idx_type new_nnz = lhs.nnz ();
 
+	  OCTAVE_LOCAL_BUFFER (octave_idx_type, rhs_idx, n);
+	  if (! lhs_idx.is_colon ())
+	    {
+	      // Ok here we have to be careful with the indexing,
+	      // to treat cases like "a([3,2,1]) = b", and still 
+	      // handle the need for strict sorting of the sparse 
+	      // elements.
+	      OCTAVE_LOCAL_BUFFER (octave_idx_vector_sort *, sidx, n);
+	      OCTAVE_LOCAL_BUFFER (octave_idx_vector_sort, sidxX, n);
+
+	      for (octave_idx_type i = 0; i < n; i++)
+		{
+		  sidx[i] = &sidxX[i];
+		  sidx[i]->i = lhs_idx.elem(i);
+		  sidx[i]->idx = i;
+		}
+			  
+	      OCTAVE_QUIT;
+	      octave_sort<octave_idx_vector_sort *> 
+		sort (octave_idx_vector_comp);
+
+	      sort.sort (sidx, n);
+
+	      intNDArray<octave_idx_type> new_idx (dim_vector (n,1));
+
+	      for (octave_idx_type i = 0; i < n; i++)
+		{
+		  new_idx.xelem(i) = sidx[i]->i + 1;
+		  rhs_idx[i] = sidx[i]->idx;
+		}
+
+	      lhs_idx = idx_vector (new_idx);
+	    }
+	  else
+	    for (octave_idx_type i = 0; i < n; i++)
+	      rhs_idx[i] = i;
+
 	  // First count the number of non-zero elements
 	  for (octave_idx_type i = 0; i < n; i++)
 	    {
@@ -1961,7 +1999,7 @@
 	      octave_idx_type ii = lhs_idx.elem (i);
 	      if (ii < lhs_len && c_lhs.elem(ii) != LT ())
 		new_nnz--;
-	      if (rhs.elem(i) != RT ())
+	      if (rhs.elem(rhs_idx[i]) != RT ())
 		new_nnz++;
 	    }
 
@@ -1992,7 +2030,7 @@
 		    }
 		  else
 		    {
-		      RT rtmp = rhs.elem (j);
+		      RT rtmp = rhs.elem (rhs_idx[j]);
 		      if (rtmp != RT ())
 			{
 			  tmp.xdata (kk) = rtmp;
@@ -2031,6 +2069,7 @@
 		      while (ic <= ii)
 			tmp.xcidx (ic++) = kk;
 		      tmp.xdata (kk) = c_lhs.data (i);
+		      tmp.xridx (kk++) = 0;
 		      i++;
 		      while (ii < nc && c_lhs.cidx(ii+1) <= i)
 			ii++;
@@ -2040,9 +2079,12 @@
 		      while (ic <= jj)
 			tmp.xcidx (ic++) = kk;
 
-		      RT rtmp = rhs.elem (j);
+		      RT rtmp = rhs.elem (rhs_idx[j]);
 		      if (rtmp != RT ())
-			tmp.xdata (kk) = rtmp;
+			{
+			  tmp.xdata (kk) = rtmp;
+			  tmp.xridx (kk++) = 0;
+			}
 		      if (ii == jj)
 			{
 			  i++;
@@ -2053,7 +2095,6 @@
 		      if (j < n)
 			jj = lhs_idx.elem(j);
 		    }
-		  tmp.xridx (kk++) = 0;
 		}
 
 	      for (octave_idx_type iidx = ic; iidx < max_idx+1; iidx++)
@@ -2067,6 +2108,7 @@
 	  octave_idx_type new_nnz = lhs.nnz ();
 	  RT scalar = rhs.elem (0);
 	  bool scalar_non_zero = (scalar != RT ());
+	  lhs_idx.sort (true);
 
 	  // First count the number of non-zero elements
 	  if (scalar != RT ())
@@ -2260,12 +2302,10 @@
 
   if (n_idx == 2)
     {
-      octave_idx_type n = idx_i.freeze (lhs_nr, "row", true, liboctave_wrore_flag);
-      idx_i.sort (true);
-
-      octave_idx_type m = idx_j.freeze (lhs_nc, "column", true, liboctave_wrore_flag);
-      idx_j.sort (true);
-
+      octave_idx_type n = idx_i.freeze (lhs_nr, "row", true, 
+					liboctave_wrore_flag);
+      octave_idx_type m = idx_j.freeze (lhs_nc, "column", true, 
+					liboctave_wrore_flag);
 
       int idx_i_is_colon = idx_i.is_colon ();
       int idx_j_is_colon = idx_j.is_colon ();
@@ -2291,14 +2331,17 @@
 
 		  if (n > 0 && m > 0)
 		    {
+		      idx_i.sort (true);
+		      idx_j.sort (true);
+
 		      octave_idx_type max_row_idx = idx_i_is_colon ? rhs_nr : 
 			idx_i.max () + 1;
 		      octave_idx_type max_col_idx = idx_j_is_colon ? rhs_nc : 
 			idx_j.max () + 1;
-		      octave_idx_type new_nr = max_row_idx > lhs_nr ? max_row_idx : 
-			lhs_nr;
-		      octave_idx_type new_nc = max_col_idx > lhs_nc ? max_col_idx : 
-			lhs_nc;
+		      octave_idx_type new_nr = max_row_idx > lhs_nr ? 
+			max_row_idx : lhs_nr;
+		      octave_idx_type new_nc = max_col_idx > lhs_nc ? 
+			max_col_idx : lhs_nc;
 		      RT scalar = rhs.elem (0, 0);
 
 		      // Count the number of non-zero terms
@@ -2399,10 +2442,88 @@
 			idx_i.max () + 1;
 		      octave_idx_type max_col_idx = idx_j_is_colon ? rhs_nc : 
 			idx_j.max () + 1;
-		      octave_idx_type new_nr = max_row_idx > lhs_nr ? max_row_idx : 
-			lhs_nr;
-		      octave_idx_type new_nc = max_col_idx > lhs_nc ? max_col_idx : 
-			lhs_nc;
+		      octave_idx_type new_nr = max_row_idx > lhs_nr ?
+			max_row_idx : lhs_nr;
+		      octave_idx_type new_nc = max_col_idx > lhs_nc ? 
+			max_col_idx : lhs_nc;
+
+		      OCTAVE_LOCAL_BUFFER (octave_idx_type, rhs_idx_i, n);
+		      if (! idx_i.is_colon ())
+			{
+			  // Ok here we have to be careful with the indexing,
+			  // to treat cases like "a([3,2,1],:) = b", and still 
+			  // handle the need for strict sorting of the sparse 
+			  // elements.
+			  OCTAVE_LOCAL_BUFFER (octave_idx_vector_sort *,
+					       sidx, n);
+			  OCTAVE_LOCAL_BUFFER (octave_idx_vector_sort,
+					       sidxX, n);
+
+			  for (octave_idx_type i = 0; i < n; i++)
+			    {
+			      sidx[i] = &sidxX[i];
+			      sidx[i]->i = idx_i.elem(i);
+			      sidx[i]->idx = i;
+			    }
+			  
+			  OCTAVE_QUIT;
+			  octave_sort<octave_idx_vector_sort *> 
+			    sort (octave_idx_vector_comp);
+
+			  sort.sort (sidx, n);
+
+			  intNDArray<octave_idx_type> new_idx (dim_vector (n,1));
+
+			  for (octave_idx_type i = 0; i < n; i++)
+			    {
+			      new_idx.xelem(i) = sidx[i]->i + 1;
+			      rhs_idx_i[i] = sidx[i]->idx;
+			    }
+
+			  idx_i = idx_vector (new_idx);
+			}
+		      else
+			for (octave_idx_type i = 0; i < n; i++)
+			  rhs_idx_i[i] = i;
+
+		      OCTAVE_LOCAL_BUFFER (octave_idx_type, rhs_idx_j, m);
+		      if (! idx_j.is_colon ())
+			{
+			  // Ok here we have to be careful with the indexing,
+			  // to treat cases like "a([3,2,1],:) = b", and still 
+			  // handle the need for strict sorting of the sparse 
+			  // elements.
+			  OCTAVE_LOCAL_BUFFER (octave_idx_vector_sort *,
+					       sidx, m);
+			  OCTAVE_LOCAL_BUFFER (octave_idx_vector_sort,
+					       sidxX, m);
+
+			  for (octave_idx_type i = 0; i < m; i++)
+			    {
+			      sidx[i] = &sidxX[i];
+			      sidx[i]->i = idx_j.elem(i);
+			      sidx[i]->idx = i;
+			    }
+			  
+			  OCTAVE_QUIT;
+			  octave_sort<octave_idx_vector_sort *> 
+			    sort (octave_idx_vector_comp);
+
+			  sort.sort (sidx, m);
+
+			  intNDArray<octave_idx_type> new_idx (dim_vector (m,1));
+
+			  for (octave_idx_type i = 0; i < m; i++)
+			    {
+			      new_idx.xelem(i) = sidx[i]->i + 1;
+			      rhs_idx_j[i] = sidx[i]->idx;
+			    }
+
+			  idx_j = idx_vector (new_idx);
+			}
+		      else
+			for (octave_idx_type i = 0; i < m; i++)
+			  rhs_idx_j[i] = i;
 
 		      // Count the number of non-zero terms
 		      octave_idx_type new_nnz = lhs.nnz ();
@@ -2430,7 +2551,7 @@
 				    }
 				}
 			      
-			      if (rhs.elem(i,j) != RT ())
+			      if (rhs.elem(rhs_idx_i[i],rhs_idx_j[j]) != RT ())
 				new_nnz++;
 			    }
 			}
@@ -2453,7 +2574,8 @@
 
 				  if (iii < n && ii == i)
 				    {
-				      RT rtmp = rhs.elem (iii, jji);
+				      RT rtmp = rhs.elem (rhs_idx_i[iii], 
+							  rhs_idx_j[jji]);
 				      if (rtmp != RT ())
 					{
 					  stmp.data(kk) = rtmp;
@@ -2529,8 +2651,8 @@
 	{
 	  octave_idx_type lhs_len = lhs.length ();
 
-	  octave_idx_type n = idx_i.freeze (lhs_len, 0, true, liboctave_wrore_flag);
-	  idx_i.sort (true);
+	  octave_idx_type n = idx_i.freeze (lhs_len, 0, true, 
+					    liboctave_wrore_flag);
 
 	  if (idx_i)
 	    {
@@ -2570,7 +2692,6 @@
       else if (lhs_nr == 1)
 	{
 	  idx_i.freeze (lhs_nc, "vector", true, liboctave_wrore_flag);
-	  idx_i.sort (true);
 
 	  if (idx_i)
 	    {
@@ -2584,7 +2705,6 @@
       else if (lhs_nc == 1)
 	{
 	  idx_i.freeze (lhs_nr, "vector", true, liboctave_wrore_flag);
-	  idx_i.sort (true);
 
 	  if (idx_i)
 	    {
@@ -2608,7 +2728,6 @@
 	  octave_idx_type lhs_len = lhs.length ();
 
 	  octave_idx_type len = idx_i.freeze (lhs_nr * lhs_nc, "matrix");
-	  idx_i.sort (true);
 
 	  if (idx_i)
 	    {
@@ -2628,6 +2747,45 @@
 	      else if (len == rhs_nr * rhs_nc)
 		{
 		  octave_idx_type new_nnz = lhs_nz;
+		  OCTAVE_LOCAL_BUFFER (octave_idx_type, rhs_idx, len);
+		  
+		  if (! idx_i.is_colon ())
+		    {
+		      // Ok here we have to be careful with the indexing, to
+		      // treat cases like "a([3,2,1]) = b", and still handle
+		      // the need for strict sorting of the sparse elements.
+
+		      OCTAVE_LOCAL_BUFFER (octave_idx_vector_sort *, sidx, 
+					   len);
+		      OCTAVE_LOCAL_BUFFER (octave_idx_vector_sort, sidxX, 
+					   len);
+
+		      for (octave_idx_type i = 0; i < len; i++)
+			{
+			  sidx[i] = &sidxX[i];
+			  sidx[i]->i = idx_i.elem(i);
+			  sidx[i]->idx = i;
+			}
+
+		      OCTAVE_QUIT;
+		      octave_sort<octave_idx_vector_sort *> 
+			sort (octave_idx_vector_comp);
+
+		      sort.sort (sidx, len);
+
+		      intNDArray<octave_idx_type> new_idx (dim_vector (len,1));
+
+		      for (octave_idx_type i = 0; i < len; i++)
+			{
+			  new_idx.xelem(i) = sidx[i]->i + 1;
+			  rhs_idx[i] = sidx[i]->idx;
+			}
+
+		      idx_i = idx_vector (new_idx);
+		    }
+		  else
+		    for (octave_idx_type i = 0; i < len; i++)
+		      rhs_idx[i] = i;
 
 		  // First count the number of non-zero elements
 		  for (octave_idx_type i = 0; i < len; i++)
@@ -2637,7 +2795,7 @@
 		      octave_idx_type ii = idx_i.elem (i);
 		      if (ii < lhs_len && c_lhs.elem(ii) != LT ())
 			new_nnz--;
-		      if (rhs.elem(i) != RT ())
+		      if (rhs.elem(rhs_idx[i]) != RT ())
 			new_nnz++;
 		    }
 
@@ -2679,7 +2837,7 @@
 			{
 			  while (kc <= jc)
 			    stmp.xcidx (kc++) = kk;
-			  RT rtmp = rhs.elem (j);
+			  RT rtmp = rhs.elem (rhs_idx[j]);
 			  if (rtmp != RT ())
 			    {
 			      stmp.xdata (kk) = rtmp;
@@ -2704,8 +2862,7 @@
 		    }
 
 		  for (octave_idx_type iidx = kc; iidx < lhs_nc+1; iidx++)
-		    stmp.xcidx(iidx) = kk;
-		  
+		    stmp.xcidx(iidx) = kk; 
 
 		  lhs = stmp;
 		}
@@ -2713,6 +2870,7 @@
 		{
 		  RT scalar = rhs.elem (0, 0);
 		  octave_idx_type new_nnz = lhs_nz;
+		  idx_i.sort (true);
 
 		  // First count the number of non-zero elements
 		  if (scalar != RT ())
--- a/liboctave/sparse-sort.cc
+++ b/liboctave/sparse-sort.cc
@@ -39,7 +39,7 @@
 octave_sparse_sidxl_comp (octave_sparse_sort_idxl* i, 
 			  octave_sparse_sort_idxl* j)
 {
-  int tmp = i->c - j->c;
+  octave_idx_type tmp = i->c - j->c;
   if (tmp < 0)
     return true;
   else if (tmp > 0)
@@ -50,6 +50,18 @@
 // Instantiate the sparse sorting class
 template class octave_sort<octave_sparse_sort_idxl *>;
 
+// Need to know the original order of the sorted indexes in
+// sparse assignments, and this class does that
+bool
+octave_idx_vector_comp (octave_idx_vector_sort* i,
+			octave_idx_vector_sort* j)
+{
+  return (i->i < j->i);
+}
+
+// Instantiate the sparse index sorting class
+template class octave_sort<octave_idx_vector_sort *>;
+
 // Instantiate the sorting class of octave_idx_type, need in MUL macro
 template class octave_sort<octave_idx_type>;
 
--- a/liboctave/sparse-sort.h
+++ b/liboctave/sparse-sort.h
@@ -28,15 +28,26 @@
 class
 octave_sparse_sort_idxl
 {
- public:
-  unsigned int r;
-  unsigned int c;
-  unsigned int idx; 
+public:
+  octave_idx_type r;
+  octave_idx_type c;
+  octave_idx_type idx; 
 };
 
 bool octave_sparse_sidxl_comp (octave_sparse_sort_idxl* i,
 			       octave_sparse_sort_idxl* j);
 
+class
+octave_idx_vector_sort
+{
+public:
+  octave_idx_type i;
+  octave_idx_type idx;
+};
+
+bool octave_idx_vector_comp (octave_idx_vector_sort* i,
+			     octave_idx_vector_sort* j);
+
 #endif
 
 /*
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,3 +1,7 @@
+2006-01-31  Kim Hansen  <kim@i9.dk>
+
+	* __glpk__.cc (F_glpk__): Fix for sparse matrices.
+
 2006-01-30  John W. Eaton  <jwe@octave.org>
 
 	* gripes.cc (gripe_wrong_type_arg (const char*, const char*, bool)):
--- a/src/DLD-FUNCTIONS/__glpk__.cc
+++ b/src/DLD-FUNCTIONS/__glpk__.cc
@@ -482,7 +482,7 @@
     }
   else
     {
-      SparseMatrix A (args(1).matrix_value ()); // get the sparse matrix
+      SparseMatrix A = args(1).sparse_matrix_value (); // get the sparse matrix
 
       if (error_state)
 	{
@@ -509,7 +509,7 @@
 	    nz++;
 	    rn(nz) = A.ridx(i) + 1;
 	    cn(nz) = j + 1;
-	    a(nz) = A(i,j);
+	    a(nz) = A.data(i);
 	  }
     }
 
--- a/test/ChangeLog
+++ b/test/ChangeLog
@@ -1,3 +1,8 @@
+2006-01-21  David Bateman  <dbateman@free.fr>
+
+        * build_sparsetest.sh: Add new un-ordered indexing, assignment and
+        deletion tests.
+
 2006-01-13  Bill Denney  <bill@givebillmoney.com>
 
 	* test_system.m: Use filesep instead of "/" where needed.
--- a/test/build_sparse_tests.sh
+++ b/test/build_sparse_tests.sh
@@ -772,6 +772,7 @@
 %!assert(sparse(as(idx),true),sparse(af(idx),true));
 %!assert(as(idx),sparse(af(idx),true));
 %!assert(as(idx'),sparse(af(idx'),true));
+%!assert(as(flipud(idx(:))),sparse(af(flipud(idx(:))),true))
 %!assert(as([idx,idx]),sparse(af([idx,idx]),true));
 %!error(as(reshape([idx;idx],[1,length(idx),2])));
 
@@ -780,6 +781,33 @@
 %!assert(as(ridx,:), sparse(af(ridx,:),true))
 %!assert(as(:,cidx), sparse(af(:,cidx),true))
 %!assert(as(:,:), sparse(af(:,:),true))
+%!assert(as((size(as,1):-1:1),:),sparse(af((size(af,1):-1:1),:),true))
+%!assert(as(:,(size(as,2):-1:1)),sparse(af(:,(size(af,2):-1:1)),true))
+
+%% Assignment test
+%!test
+%! ts=as;ts(:,:)=ts(fliplr(1:size(as,1)),:);
+%! tf=af;tf(:,:)=tf(fliplr(1:size(af,1)),:);
+%! assert(ts,sparse(tf,true));
+%!test
+%! ts=as;ts(fliplr(1:size(as,1)),:)=ts;
+%! tf=af;tf(fliplr(1:size(af,1)),:)=tf;
+%! assert(ts,sparse(tf,true));
+%!test
+%! ts=as;ts(:,fliplr(1:size(as,2)))=ts;
+%! tf=af;tf(:,fliplr(1:size(af,2)))=tf;
+%! assert(ts,sparse(tf,true));
+%!test
+%! ts(fliplr(1:size(as,1)))=as(:,1);tf(fliplr(1:size(af,1)))=af(:,1);
+%! assert(ts,sparse(tf,true));
+
+%% Deletion tests
+%!test
+%! ts=as;ts(1,:)=[];tf=af;tf(1,:)=[];
+%! assert(ts,sparse(tf,true));
+%!test
+%! ts=as;ts(:,1)=[];tf=af;tf(:,1)=[];
+%! assert(ts,sparse(tf,true));
 
 %% Test 'end' keyword
 %!assert(as(end),af(end))