# HG changeset patch # User Rik # Date 1238086055 25200 # Node ID cfad8f9a77fac4b44a183631f055b4d9fb647a33 # Parent 97aa01a85ea423314eb2f6ac219e7f388d97246d# Parent 2df28ad88b0e75468a8f5cd0de6f32c6ab520340 Documentation merge to main branch diff --git a/.hgtags b/.hgtags --- a/.hgtags +++ b/.hgtags @@ -41,3 +41,4 @@ 739141cde75a22b9ac9b5b327ea30c61fd357748 ss-3-1-52 bd1b1fe9c6e952202bf5878cd042a168f6df6ddd ss-3-1-53 3c3cbe8f18e03c73cc34494c8037d8a45a30bf7d ss-3-1-54 +5e276a0b999747b72af34e72eb63f64640da6399 ss-3-1-55 diff --git a/liboctave/Array.cc b/liboctave/Array.cc --- a/liboctave/Array.cc +++ b/liboctave/Array.cc @@ -2487,6 +2487,69 @@ return idx; } +template +Array +Array::find (octave_idx_type n, bool backward) const +{ + Array retval; + const T *src = data (); + octave_idx_type nel = nelem (); + const T zero = T (); + if (n < 0) + { + // We want all elements, which means we'll almost surely need + // to resize. So count first, then allocate array of exact size. + octave_idx_type cnt = 0; + for (octave_idx_type i = 0; i < nel; i++) + cnt += src[i] != zero; + + retval = Array (cnt); + octave_idx_type *dest = retval.fortran_vec (); + for (octave_idx_type i = 0; i < nel; i++) + if (src[i] != zero) *dest++ = i; + } + else + { + // We want a fixed max number of elements, usually small. So be + // optimistic, alloc the array in advance, and then resize if + // needed. + retval = Array (n); + if (backward) + { + // Do the search as a series of successive single-element searches. + octave_idx_type k = 0, l = nel - 1; + for (; k < n; k++) + { + for (;l >= 0 && src[l] == zero; l--) ; + if (l >= 0) + retval(k) = l--; + else + break; + } + if (k < n) + retval.resize (k); + } + else + { + // Do the search as a series of successive single-element searches. + octave_idx_type k = 0, l = 0; + for (; k < n; k++) + { + for (;l != nel && src[l] == zero; l++) ; + if (l != nel) + retval(k) = l++; + else + break; + } + if (k < n) + retval.resize (k); + } + } + + return retval; +} + + #define INSTANTIATE_ARRAY_SORT(T) template class octave_sort; #define NO_INSTANTIATE_ARRAY_SORT(T) \ @@ -2520,7 +2583,9 @@ template <> Array \ Array::lookup (const Array&, sortmode, bool, bool) const \ { return Array (); } \ - +template <> Array \ +Array::find (octave_idx_type, bool) const\ +{ return Array (); } \ template diff --git a/liboctave/Array.h b/liboctave/Array.h --- a/liboctave/Array.h +++ b/liboctave/Array.h @@ -255,7 +255,8 @@ size_t byte_size (void) const { return numel () * sizeof (T); } - dim_vector dims (void) const { return dimensions; } + // Return a const-reference so that dims ()(i) works efficiently. + const dim_vector& dims (void) const { return dimensions; } Array squeeze (void) const; @@ -428,6 +429,8 @@ bool is_empty (void) const { return numel () == 0; } + bool is_vector (void) const { return dimensions.is_vector (); } + Array transpose (void) const; Array hermitian (T (*fcn) (const T&) = 0) const; @@ -579,6 +582,10 @@ Array lookup (const Array& values, sortmode mode = UNSORTED, bool linf = false, bool rinf = false) const; + // Find indices of (at most n) nonzero elements. If n is specified, backward + // specifies search from backward. + Array find (octave_idx_type n = -1, bool backward = false) const; + Array diag (octave_idx_type k = 0) const; template diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog --- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,5 +1,23 @@ +2009-03-26 Jaroslav Hajek + + * dim-vector.h (dim_vector::numel): Add optional argument, simplify. + +2009-03-26 Jaroslav Hajek + + * Array.h (Array::dims): Return a const reference. + (Array::is_vector): New method. + +2009-03-26 Jaroslav Hajek + + * Array.cc (Array::find): New method. + * Array.h: Declare it. + 2009-03-25 John W. Eaton + * EIG.cc (EIG::init (const Matrix&, bool), + EIG::init (const Matrix&, const Matrix&, bool)): + Avoid volatile declaration for tmp variable. + * Makefile.in (MATRIX_INC): Add Sparse-diag-op-defs.h and Sparse-perm-op-defs.h to the list. diff --git a/liboctave/EIG.cc b/liboctave/EIG.cc --- a/liboctave/EIG.cc +++ b/liboctave/EIG.cc @@ -161,8 +161,8 @@ Array wi (n); double *pwi = wi.fortran_vec (); - volatile octave_idx_type nvr = calc_ev ? n : 0; - Matrix vr (nvr, nvr); + octave_idx_type tnvr = calc_ev ? n : 0; + Matrix vr (tnvr, tnvr); double *pvr = vr.fortran_vec (); octave_idx_type lwork = -1; @@ -204,6 +204,7 @@ } lambda.resize (n); + octave_idx_type nvr = calc_ev ? n : 0; v.resize (nvr, nvr); for (octave_idx_type j = 0; j < n; j++) @@ -507,8 +508,8 @@ Array beta (n); double *pbeta = beta.fortran_vec (); - volatile octave_idx_type nvr = calc_ev ? n : 0; - Matrix vr (nvr, nvr); + octave_idx_type tnvr = calc_ev ? n : 0; + Matrix vr (tnvr, tnvr); double *pvr = vr.fortran_vec (); octave_idx_type lwork = -1; @@ -554,6 +555,7 @@ } lambda.resize (n); + octave_idx_type nvr = calc_ev ? n : 0; v.resize (nvr, nvr); for (octave_idx_type j = 0; j < n; j++) diff --git a/liboctave/dim-vector.h b/liboctave/dim-vector.h --- a/liboctave/dim-vector.h +++ b/liboctave/dim-vector.h @@ -310,13 +310,13 @@ // vector would have, NOT the number of dimensions (elements in the // dimension vector). - octave_idx_type numel (void) const + octave_idx_type numel (int n = 0) const { int n_dims = length (); - octave_idx_type retval = n_dims > 0 ? elem (0) : 0; + octave_idx_type retval = 1; - for (int i = 1; i < n_dims; i++) + for (int i = n; i < n_dims; i++) retval *= elem (i); return retval; diff --git a/src/ChangeLog b/src/ChangeLog --- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,5 +1,21 @@ +2009-03-26 Jaroslav Hajek + + * DLD-FUNCTIONS/find.cc + (find_nonzero_elem_idx (const Array&, ...)): Simplify. + Instantiate for bool and octave_int types. + (find_nonzero_elem_idx (const Sparse&, ...)): + Instantiate for bool. + (Ffind): Handle bool and octave_int cases. + 2009-03-25 John W. Eaton + * version.h (OCTAVE_VERSION): Now 3.1.55+. + (OCTAVE_API_VERSION): Now api-v37+. + + * version.h (OCTAVE_VERSION): Now 3.1.55. + (OCTAVE_API_VERSION): Now api-v37. + (OCTAVE_RELEASE_DATE): Now 2009-03-25. + * DLD-FUNCTIONS/find.cc (find_nonzero_elem_idx): Also return [](0x0) if the array has 0 rows and it is not a column vector. diff --git a/src/DLD-FUNCTIONS/find.cc b/src/DLD-FUNCTIONS/find.cc --- a/src/DLD-FUNCTIONS/find.cc +++ b/src/DLD-FUNCTIONS/find.cc @@ -43,155 +43,75 @@ { octave_value_list retval ((nargout == 0 ? 1 : nargout), Matrix ()); - octave_idx_type count = 0; - - octave_idx_type nel = nda.nelem (); - - // Set the starting element to the correct value based on the - // direction to search. - octave_idx_type k = 0; - if (direction == -1) - k = nel - 1; - - // Search in the default range. - octave_idx_type start_el = -1; - octave_idx_type end_el = -1; - - // Search for the number of elements to return. - while (k < nel && k > -1 && n_to_find != count) - { - OCTAVE_QUIT; - - if (nda(k) != T ()) - { - end_el = k; - if (start_el == -1) - start_el = k; - count++; - } - k = k + direction; - } - - // Reverse the range if we're looking backward. - if (direction == -1) - { - octave_idx_type tmp_el = start_el; - start_el = end_el; - end_el = tmp_el; - } - // Fix an off by one error. - end_el++; - - // If the original argument was a row vector, force a row vector of - // the overall indices to be returned. But see below for scalar - // case... - - octave_idx_type result_nr = count; - octave_idx_type result_nc = 1; - - bool column_vector_arg = false; - bool scalar_arg = false; - - if (nda.ndims () == 2) - { - octave_idx_type nr = nda.rows (); - octave_idx_type nc = nda.columns (); + Array idx; + if (n_to_find >= 0) + idx = nda.find (n_to_find, direction == -1); + else + idx = nda.find (); - if (nr == 1) - { - result_nr = 1; - result_nc = count; - - scalar_arg = (nc == 1); - } - else if (nc == 1) - column_vector_arg = true; - } - - Matrix idx (result_nr, result_nc); - - Matrix i_idx (result_nr, result_nc); - Matrix j_idx (result_nr, result_nc); - - ArrayN val (dim_vector (result_nr, result_nc)); - - if (count > 0) - { - count = 0; - - octave_idx_type nr = nda.rows (); - - octave_idx_type i = 0; - - // Search for elements to return. Only search the region where - // there are elements to be found using the count that we want - // to find. + // Fixup idx dimensions, for Matlab compatibility. + // find(zeros(0,0)) -> zeros(0,0) + // find(zeros(1,0)) -> zeros(1,0) + // find(zeros(0,1)) -> zeros(0,1) + // find(zeros(0,X)) -> zeros(0,1) + // find(zeros(1,1)) -> zeros(0,0) !!!! WHY? + // find(zeros(0,1,0)) -> zeros(0,0) + // find(zeros(0,1,0,1)) -> zeros(0,0) etc + // FIXME: I don't believe this is right. Matlab seems to violate its own docs + // here, because a scalar *is* a row vector. - // For compatibility, all N-d arrays are handled as if they are - // 2-d, with the number of columns equal to "prod (dims (2:end))". - - for (k = start_el; k < end_el; k++) - { - OCTAVE_QUIT; - - if (nda(k) != T ()) - { - idx(count) = k + 1; - - octave_idx_type xr = k % nr; - i_idx(count) = xr + 1; - j_idx(count) = (k - xr) / nr + 1; - - val(count) = nda(k); - - count++; - } - - i++; - } - } - else if (scalar_arg || (nda.rows () == 0 && ! column_vector_arg)) - { - idx.resize (0, 0); - - i_idx.resize (0, 0); - j_idx.resize (0, 0); - - val.resize (dim_vector (0, 0)); - } + if ((nda.numel () == 1 && idx.is_empty ()) + || (nda.rows () == 0 && nda.dims ().numel (1) == 0)) + idx = idx.reshape (dim_vector (0, 0)); + else if (nda.rows () == 1 && nda.ndims () == 2) + idx = idx.reshape (dim_vector (1, idx.length ())); switch (nargout) { default: case 3: - retval(2) = val; + retval(2) = ArrayN (nda.index (idx_vector (idx))); // Fall through! case 2: - retval(1) = j_idx; - retval(0) = i_idx; - break; + { + Array jdx (idx.dims ()); + octave_idx_type n = idx.length (), nr = nda.rows (); + for (octave_idx_type i = 0; i < n; i++) + { + jdx.xelem (i) = idx.xelem (i) / nr; + idx.xelem (i) %= nr; + } + retval(1) = NDArray (jdx, true); + } + // Fall through! case 1: case 0: - retval(0) = idx; + retval(0) = NDArray (idx, true); break; } return retval; } -template octave_value_list find_nonzero_elem_idx (const Array&, int, - octave_idx_type, int); - -template octave_value_list find_nonzero_elem_idx (const Array&, int, - octave_idx_type, int); +#define INSTANTIATE_FIND_ARRAY(T) \ +template octave_value_list find_nonzero_elem_idx (const Array&, int, \ + octave_idx_type, int) -template octave_value_list find_nonzero_elem_idx (const Array&, int, - octave_idx_type, int); - -template octave_value_list find_nonzero_elem_idx (const Array&, - int, octave_idx_type, int); +INSTANTIATE_FIND_ARRAY(double); +INSTANTIATE_FIND_ARRAY(float); +INSTANTIATE_FIND_ARRAY(Complex); +INSTANTIATE_FIND_ARRAY(FloatComplex); +INSTANTIATE_FIND_ARRAY(bool); +INSTANTIATE_FIND_ARRAY(octave_int8); +INSTANTIATE_FIND_ARRAY(octave_int16); +INSTANTIATE_FIND_ARRAY(octave_int32); +INSTANTIATE_FIND_ARRAY(octave_int64); +INSTANTIATE_FIND_ARRAY(octave_uint8); +INSTANTIATE_FIND_ARRAY(octave_uint16); +INSTANTIATE_FIND_ARRAY(octave_uint32); +INSTANTIATE_FIND_ARRAY(octave_uint64); template octave_value_list @@ -342,6 +262,9 @@ template octave_value_list find_nonzero_elem_idx (const Sparse&, int, octave_idx_type, int); +template octave_value_list find_nonzero_elem_idx (const Sparse&, int, + octave_idx_type, int); + octave_value_list find_nonzero_elem_idx (const PermMatrix& v, int nargout, octave_idx_type n_to_find, int direction) @@ -561,7 +484,51 @@ octave_value arg = args(0); - if (arg.is_sparse_type ()) + if (arg.is_bool_type ()) + { + if (arg.is_sparse_type ()) + { + SparseBoolMatrix v = arg.sparse_bool_matrix_value (); + + if (! error_state) + retval = find_nonzero_elem_idx (v, nargout, + n_to_find, direction); + } + else + { + boolNDArray v = arg.bool_array_value (); + + if (! error_state) + retval = find_nonzero_elem_idx (v, nargout, + n_to_find, direction); + } + } + else if (arg.is_integer_type ()) + { +#define DO_INT_BRANCH(INTT) \ + else if (arg.is_ ## INTT ## _type ()) \ + { \ + INTT ## NDArray v = arg.INTT ## _array_value (); \ + \ + if (! error_state) \ + retval = find_nonzero_elem_idx (v, nargout, \ + n_to_find, direction);\ + } + + if (false) + ; + DO_INT_BRANCH (int8) + DO_INT_BRANCH (int16) + DO_INT_BRANCH (int32) + DO_INT_BRANCH (int64) + DO_INT_BRANCH (uint8) + DO_INT_BRANCH (uint16) + DO_INT_BRANCH (uint32) + DO_INT_BRANCH (uint64) + else + panic_impossible (); + } + else if (arg.is_sparse_type ()) { if (arg.is_real_type ()) { diff --git a/src/version.h b/src/version.h --- a/src/version.h +++ b/src/version.h @@ -25,11 +25,11 @@ #if !defined (octave_version_h) #define octave_version_h 1 -#define OCTAVE_VERSION "3.1.54+" +#define OCTAVE_VERSION "3.1.55" -#define OCTAVE_API_VERSION "api-v36+" +#define OCTAVE_API_VERSION "api-v37" -#define OCTAVE_RELEASE_DATE "2009-03-07" +#define OCTAVE_RELEASE_DATE "2009-03-25" #define OCTAVE_COPYRIGHT "Copyright (C) 2009 John W. Eaton and others."