# HG changeset patch # User dbateman # Date 1138708667 0 # Node ID 2c66c36d2698676ee17d0d09fa75f748f095c36c # Parent dfa2da0563e5e93dd6539ccdda8e56331c3c5af3 [project @ 2006-01-31 11:57:47 by dbateman] diff --git a/configure.in b/configure.in --- 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= diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog --- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,18 @@ +2006-01-21 David Bateman + + * sparse-sort.cc (bool octave_sparse_sidxl_comp): 64-bit fix. + (bool octave_idx_vector_comp): New function. + (template class octave_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&, Sparse&)): Treat cases of + unordered LHS indexes in assignment using new octave_idx_vector_sort + class. + (int assign(Sparse&, Sparse&)): ditto. + 2006-01-30 John W. Eaton * so-array.h (streamoff_array::nnz): New funtion. diff --git a/liboctave/Sparse.cc b/liboctave/Sparse.cc --- 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 + sort (octave_idx_vector_comp); + + sort.sort (sidx, n); + + intNDArray 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 + sort (octave_idx_vector_comp); + + sort.sort (sidx, n); + + intNDArray 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 + sort (octave_idx_vector_comp); + + sort.sort (sidx, m); + + intNDArray 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 + sort (octave_idx_vector_comp); + + sort.sort (sidx, len); + + intNDArray 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 ()) diff --git a/liboctave/sparse-sort.cc b/liboctave/sparse-sort.cc --- 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; +// 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; + // Instantiate the sorting class of octave_idx_type, need in MUL macro template class octave_sort; diff --git a/liboctave/sparse-sort.h b/liboctave/sparse-sort.h --- 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 /* diff --git a/src/ChangeLog b/src/ChangeLog --- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,3 +1,7 @@ +2006-01-31 Kim Hansen + + * __glpk__.cc (F_glpk__): Fix for sparse matrices. + 2006-01-30 John W. Eaton * gripes.cc (gripe_wrong_type_arg (const char*, const char*, bool)): 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 @@ -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); } } diff --git a/test/ChangeLog b/test/ChangeLog --- a/test/ChangeLog +++ b/test/ChangeLog @@ -1,3 +1,8 @@ +2006-01-21 David Bateman + + * build_sparsetest.sh: Add new un-ordered indexing, assignment and + deletion tests. + 2006-01-13 Bill Denney * test_system.m: Use filesep instead of "/" where needed. diff --git a/test/build_sparse_tests.sh b/test/build_sparse_tests.sh --- 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))