# HG changeset patch # User David Bateman # Date 1206054504 -3600 # Node ID 36594d5bbe132bcb25f14fcd87cc524894e810fd # Parent 56012914972a9aef16f517407167c677bdd5f207 Move diag function into the octave_value class diff --git a/liboctave/Array.cc b/liboctave/Array.cc --- a/liboctave/Array.cc +++ b/liboctave/Array.cc @@ -2624,6 +2624,93 @@ #endif +template +Array +Array::diag (octave_idx_type k) const +{ + dim_vector dv = dims (); + octave_idx_type nd = dv.length (); + Array d; + + if (nd > 2) + (*current_liboctave_error_handler) ("Matrix must be 2-dimensional"); + else + { + octave_idx_type nnr = dv (0); + octave_idx_type nnc = dv (1); + + if (nnr == 0 || nnc == 0) + ; // do nothing + else if (nnr != 1 && nnc != 1) + { + if (k > 0) + nnc -= k; + else if (k < 0) + nnr += k; + + if (nnr > 0 && nnc > 0) + { + octave_idx_type ndiag = (nnr < nnc) ? nnr : nnc; + + d.resize (dim_vector (ndiag, 1)); + + if (k > 0) + { + for (octave_idx_type i = 0; i < ndiag; i++) + d.xelem (i) = elem (i, i+k); + } + else if (k < 0) + { + for (octave_idx_type i = 0; i < ndiag; i++) + d.xelem (i) = elem (i-k, i); + } + else + { + for (octave_idx_type i = 0; i < ndiag; i++) + d.xelem (i) = elem (i, i); + } + } + else + (*current_liboctave_error_handler) + ("diag: requested diagonal out of range"); + } + else if (nnr != 0 && nnc != 0) + { + octave_idx_type roff = 0; + octave_idx_type coff = 0; + if (k > 0) + { + roff = 0; + coff = k; + } + else if (k < 0) + { + roff = -k; + coff = 0; + } + + if (nnr == 1) + { + octave_idx_type n = nnc + std::abs (k); + d = Array (dim_vector (n, n), resize_fill_value (T ())); + + for (octave_idx_type i = 0; i < nnc; i++) + d.xelem (i+roff, i+coff) = elem (0, i); + } + else + { + octave_idx_type n = nnr + std::abs (k); + d = Array (dim_vector (n, n), resize_fill_value (T ())); + + for (octave_idx_type i = 0; i < nnr; i++) + d.xelem (i+roff, i+coff) = elem (i, 0); + } + } + } + + return d; +} + // FIXME -- this is a mess. template diff --git a/liboctave/Array.h b/liboctave/Array.h --- a/liboctave/Array.h +++ b/liboctave/Array.h @@ -550,6 +550,8 @@ Array sort (Array &sidx, octave_idx_type dim = 0, sortmode mode = ASCENDING) const; + Array diag (octave_idx_type k = 0) const; + template Array map (F fcn) const diff --git a/liboctave/Array2.h b/liboctave/Array2.h --- a/liboctave/Array2.h +++ b/liboctave/Array2.h @@ -136,6 +136,11 @@ return Array2 (tmp, tmp.rows (), tmp.columns ()); } + Array2 diag (octave_idx_type k) const + { + return Array::diag (k); + } + template Array2 map (F fcn) const { diff --git a/liboctave/ArrayN.h b/liboctave/ArrayN.h --- a/liboctave/ArrayN.h +++ b/liboctave/ArrayN.h @@ -149,6 +149,11 @@ return ArrayN (tmp, tmp.dims ()); } + ArrayN diag (octave_idx_type k) const + { + return Array::diag (k); + } + template ArrayN map (F fcn) const { diff --git a/liboctave/CDiagMatrix.cc b/liboctave/CDiagMatrix.cc --- a/liboctave/CDiagMatrix.cc +++ b/liboctave/CDiagMatrix.cc @@ -501,14 +501,6 @@ // other operations ComplexColumnVector -ComplexDiagMatrix::diag (void) const -{ - return diag (0); -} - -// Could be optimized... - -ComplexColumnVector ComplexDiagMatrix::diag (octave_idx_type k) const { octave_idx_type nnr = rows (); diff --git a/liboctave/CDiagMatrix.h b/liboctave/CDiagMatrix.h --- a/liboctave/CDiagMatrix.h +++ b/liboctave/CDiagMatrix.h @@ -114,8 +114,7 @@ // other operations - ComplexColumnVector diag (void) const; - ComplexColumnVector diag (octave_idx_type k) const; + ComplexColumnVector diag (octave_idx_type k = 0) const; // i/o diff --git a/liboctave/CMatrix.cc b/liboctave/CMatrix.cc --- a/liboctave/CMatrix.cc +++ b/liboctave/CMatrix.cc @@ -3316,51 +3316,10 @@ return retval; } -ComplexColumnVector -ComplexMatrix::diag (void) const -{ - return diag (0); -} - -ComplexColumnVector +ComplexMatrix ComplexMatrix::diag (octave_idx_type k) const { - octave_idx_type nnr = rows (); - octave_idx_type nnc = cols (); - if (k > 0) - nnc -= k; - else if (k < 0) - nnr += k; - - ComplexColumnVector d; - - if (nnr > 0 && nnc > 0) - { - octave_idx_type ndiag = (nnr < nnc) ? nnr : nnc; - - d.resize (ndiag); - - if (k > 0) - { - for (octave_idx_type i = 0; i < ndiag; i++) - d.elem (i) = elem (i, i+k); - } - else if (k < 0) - { - for (octave_idx_type i = 0; i < ndiag; i++) - d.elem (i) = elem (i-k, i); - } - else - { - for (octave_idx_type i = 0; i < ndiag; i++) - d.elem (i) = elem (i, i); - } - } - else - (*current_liboctave_error_handler) - ("diag: requested diagonal out of range"); - - return d; + return MArray2::diag (k); } bool diff --git a/liboctave/CMatrix.h b/liboctave/CMatrix.h --- a/liboctave/CMatrix.h +++ b/liboctave/CMatrix.h @@ -334,8 +334,7 @@ ComplexMatrix sumsq (int dim = -1) const; Matrix abs (void) const; - ComplexColumnVector diag (void) const; - ComplexColumnVector diag (octave_idx_type k) const; + ComplexMatrix diag (octave_idx_type k = 0) const; bool row_is_real_only (octave_idx_type) const; bool column_is_real_only (octave_idx_type) const; diff --git a/liboctave/CNDArray.cc b/liboctave/CNDArray.cc --- a/liboctave/CNDArray.cc +++ b/liboctave/CNDArray.cc @@ -991,6 +991,12 @@ return ::compute_index (ra_idx, dimensions); } +ComplexNDArray +ComplexNDArray::diag (octave_idx_type k) const +{ + return MArrayN::diag (k); +} + NDArray ComplexNDArray::map (dmapper fcn) const { diff --git a/liboctave/CNDArray.h b/liboctave/CNDArray.h --- a/liboctave/CNDArray.h +++ b/liboctave/CNDArray.h @@ -116,6 +116,8 @@ // bool all_elements_are_real (void) const; // bool all_integers (double& max_val, double& min_val) const; + ComplexNDArray diag (octave_idx_type k = 0) const; + typedef double (*dmapper) (const Complex&); typedef Complex (*cmapper) (const Complex&); typedef bool (*bmapper) (const Complex&); diff --git a/liboctave/CSparse.cc b/liboctave/CSparse.cc --- a/liboctave/CSparse.cc +++ b/liboctave/CSparse.cc @@ -7359,88 +7359,7 @@ SparseComplexMatrix SparseComplexMatrix::diag (octave_idx_type k) const { - octave_idx_type nnr = rows (); - octave_idx_type nnc = cols (); - - if (k > 0) - nnc -= k; - else if (k < 0) - nnr += k; - - SparseComplexMatrix d; - - if (nnr > 0 && nnc > 0) - { - octave_idx_type ndiag = (nnr < nnc) ? nnr : nnc; - - // Count the number of non-zero elements - octave_idx_type nel = 0; - if (k > 0) - { - for (octave_idx_type i = 0; i < ndiag; i++) - if (elem (i, i+k) != 0.) - nel++; - } - else if ( k < 0) - { - for (octave_idx_type i = 0; i < ndiag; i++) - if (elem (i-k, i) != 0.) - nel++; - } - else - { - for (octave_idx_type i = 0; i < ndiag; i++) - if (elem (i, i) != 0.) - nel++; - } - - d = SparseComplexMatrix (ndiag, 1, nel); - d.xcidx (0) = 0; - d.xcidx (1) = nel; - - octave_idx_type ii = 0; - if (k > 0) - { - for (octave_idx_type i = 0; i < ndiag; i++) - { - Complex tmp = elem (i, i+k); - if (tmp != 0.) - { - d.xdata (ii) = tmp; - d.xridx (ii++) = i; - } - } - } - else if ( k < 0) - { - for (octave_idx_type i = 0; i < ndiag; i++) - { - Complex tmp = elem (i-k, i); - if (tmp != 0.) - { - d.xdata (ii) = tmp; - d.xridx (ii++) = i; - } - } - } - else - { - for (octave_idx_type i = 0; i < ndiag; i++) - { - Complex tmp = elem (i, i); - if (tmp != 0.) - { - d.xdata (ii) = tmp; - d.xridx (ii++) = i; - } - } - } - } - else - (*current_liboctave_error_handler) - ("diag: requested diagonal out of range"); - - return d; + return MSparse::diag (k); } SparseMatrix diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog --- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -2,6 +2,77 @@ * oct-sparse.h: Add headers for amd.h. +2008-03-20 David Bateman + + * Array.cc (Array Array::diag (octave_idx_type) const): New + method for diag function. + * Array.h (Array diag (octave_idx_type) const): Declare it. + * Array2.h (Array2 diag (octave_idx_type) const): New method. + * MArray2.h (MArray2 diag (octave_idx_type) const): ditto. + * ArrayN.h (ArrayN diag (octave_idx_type) const): ditto. + * MArrayN.h (MArrayN diag (octave_idx_type) const): ditto. + + * Sparse.cc (Sparse Sparse::diag (octave_idx_type) const): + New method for the diag function. + * Sparse.h (Sparse diag (octave_idx_type) const): Declare it. + * MSparse.h (MSparse diag (octave_idx_type) const): New method. + + * Range.cc (Matrix Range::diag (octave_idx_type) const): + New method for the diag function. + * Range.h (Matrix diag (octave_idx_type) const): Declare it. + + * CDiagMatrix.cc (ComplexColumnVector ComplexDiagMatrix::diag + (void) const): delete. + * dDiagMatrix.cc (ColumnVector DiagMatrix::diag (void) const): delete. + * dDiagMatrix.h (ColumnVector diag (void) const): ditto. + * CMatrix.cc (ComplexColumnVector ComplexMatrix::diag (void) const): + delete. + * CMatrix.h (ComplexColumnVector diag (void) const): ditto. + * dMatrix.cc (ColumnVector Matrix::diag (void) const): ditto. + * dMatrix.h (ColumnVector diag (void) const): ditto. + * boolMatrix.cc (boolMatrix boolMatrix::diag (void) const): ditto. + * boolMatrix.h (boolMatrix diag (void) const): ditto. + * chMatrix.cc (charMatrix charMatrix::diag (void) const): ditto. + * chMatrix.h (charMatrix diag (void) const): ditto. + * intNDArray.cc (intNDArray intNDArray::diag (void) const): ditto. + * intNDArray.h (intNDArray diag (void) const): ditto. + + * CMatrix.cc (ComplexMatrix ComplexMatrix::diag (octave_idx_type) + const): Rewrite in terms of template classes function. + * CMatrix.h (ComplexMatrix diag (octave_idx_type)const ): Change + return type. + * dMatrix.cc (Matrix Matrix::diag (octave_idx_type) const): Rewrite in + terms of template classes function. + * dMatrix.h (Matrix diag (octave_idx_type) const): Change return type. + * boolMatrix.cc (boolMatrix boolMatrix::diag (octave_idx_type) const): + Rewrite in terms of template classes function. + * boolMatrix.h (boolMatrix diag (octave_idx_type) const): Change + return type. + * chMatrix.cc (charMatrix charMatrix::diag (octave_idx_type) + const): Rewrite in terms of template classes function. + + * dSparse.cc (SparseMatrix SparseMatrix::diag (octave_idx_type) const): + Rewrite in terms of template classes function. + * CSparse.cc (SparseComplexMatrix SparseComplexMatrix::diag + (octave_idx_type) const): ditto. + * boolSparse.cc (SparseBoolMatrix SparseBoolMatrix::diag + (octave_idx_type) const): ditto. + * intNDArray.cc (intNDArray intNDArray::diag + (octave_idx_type) const): ditto. + + * CNDArray.cc (ComplexNDArray ComplexNDArray::diag + (octave_idx_type) const): New method. + * CNDArray.h (ComplexNDArray diag (octave_idx_type) const): + Declare it. + * dNDArray.cc (NDArray NDArray::diag (octave_idx_type) const): New + method. + * dNDArray.h (NDArray diag (octave_idx_type) const): Declare it. + * chNDArray.cc (charNDArray charNDArray::diag + (octave_idx_type) const): New method. + * chNDArray.h (charNDArray diag (octave_idx_type) const): + Declare it. + + 2008-03-19 John W. Eaton * oct-env.cc (octave_env::do_base_pathname): Also handle rooted diff --git a/liboctave/MArray2.h b/liboctave/MArray2.h --- a/liboctave/MArray2.h +++ b/liboctave/MArray2.h @@ -81,6 +81,11 @@ MArray2 transpose (void) const { return Array2::transpose (); } + MArray2 diag (octave_idx_type k) const + { + return Array2::diag (k); + } + template MArray2 map (F fcn) const { diff --git a/liboctave/MArrayN.h b/liboctave/MArrayN.h --- a/liboctave/MArrayN.h +++ b/liboctave/MArrayN.h @@ -98,6 +98,11 @@ MArrayN squeeze (void) const { return ArrayN::squeeze (); } + MArrayN diag (octave_idx_type k) const + { + return ArrayN::diag (k); + } + template MArrayN map (F fcn) const { diff --git a/liboctave/MSparse.h b/liboctave/MSparse.h --- a/liboctave/MSparse.h +++ b/liboctave/MSparse.h @@ -112,7 +112,12 @@ { return Sparse::ipermute (vec); } - template + MSparse diag (octave_idx_type k = 0) const + { + return Sparse::diag (k); + } + + template MSparse map (F fcn) const { return Sparse::template map (fcn); diff --git a/liboctave/Range.cc b/liboctave/Range.cc --- a/liboctave/Range.cc +++ b/liboctave/Range.cc @@ -178,6 +178,40 @@ } +Matrix +Range::diag (octave_idx_type k) const +{ + octave_idx_type nnr = 1; + octave_idx_type nnc = nelem (); + Matrix d; + + if (nnr != 0) + { + octave_idx_type roff = 0; + octave_idx_type coff = 0; + if (k > 0) + { + roff = 0; + coff = k; + } + else if (k < 0) + { + roff = -k; + coff = 0; + } + + // Force cached matrix to be created + matrix_value (); + + octave_idx_type n = nnc + std::abs (k); + d = Matrix (n, n, Matrix::resize_fill_value ()); + for (octave_idx_type i = 0; i < nnc; i++) + d.xelem (i+roff, i+coff) = cache.xelem (i); + } + + return d; +} + Range Range::sort (octave_idx_type dim, sortmode mode) const { diff --git a/liboctave/Range.h b/liboctave/Range.h --- a/liboctave/Range.h +++ b/liboctave/Range.h @@ -65,6 +65,8 @@ void sort_internal (bool ascending = true); void sort_internal (Array& sidx, bool ascending = true); + Matrix diag (octave_idx_type k = 0) const; + Range sort (octave_idx_type dim = 0, sortmode mode = ASCENDING) const; Range sort (Array& sidx, octave_idx_type dim = 0, diff --git a/liboctave/Sparse.cc b/liboctave/Sparse.cc --- a/liboctave/Sparse.cc +++ b/liboctave/Sparse.cc @@ -2253,6 +2253,156 @@ return m; } +template +Sparse +Sparse::diag (octave_idx_type k) const +{ + octave_idx_type nnr = rows (); + octave_idx_type nnc = cols (); + Sparse d; + + if (nnr == 0 || nnc == 0) + ; // do nothing + else if (nnr != 1 && nnc != 1) + { + if (k > 0) + nnc -= k; + else if (k < 0) + nnr += k; + + if (nnr > 0 && nnc > 0) + { + octave_idx_type ndiag = (nnr < nnc) ? nnr : nnc; + + // Count the number of non-zero elements + octave_idx_type nel = 0; + if (k > 0) + { + for (octave_idx_type i = 0; i < ndiag; i++) + if (elem (i, i+k) != 0.) + nel++; + } + else if ( k < 0) + { + for (octave_idx_type i = 0; i < ndiag; i++) + if (elem (i-k, i) != 0.) + nel++; + } + else + { + for (octave_idx_type i = 0; i < ndiag; i++) + if (elem (i, i) != 0.) + nel++; + } + + d = Sparse (ndiag, 1, nel); + d.xcidx (0) = 0; + d.xcidx (1) = nel; + + octave_idx_type ii = 0; + if (k > 0) + { + for (octave_idx_type i = 0; i < ndiag; i++) + { + T tmp = elem (i, i+k); + if (tmp != 0.) + { + d.xdata (ii) = tmp; + d.xridx (ii++) = i; + } + } + } + else if ( k < 0) + { + for (octave_idx_type i = 0; i < ndiag; i++) + { + T tmp = elem (i-k, i); + if (tmp != 0.) + { + d.xdata (ii) = tmp; + d.xridx (ii++) = i; + } + } + } + else + { + for (octave_idx_type i = 0; i < ndiag; i++) + { + T tmp = elem (i, i); + if (tmp != 0.) + { + d.xdata (ii) = tmp; + d.xridx (ii++) = i; + } + } + } + } + else + (*current_liboctave_error_handler) + ("diag: requested diagonal out of range"); + } + else if (nnr != 0 && nnc != 0) + { + octave_idx_type roff = 0; + octave_idx_type coff = 0; + if (k > 0) + { + roff = 0; + coff = k; + } + else if (k < 0) + { + roff = -k; + coff = 0; + } + + if (nnr == 1) + { + octave_idx_type n = nnc + std::abs (k); + octave_idx_type nz = nzmax (); + d = Sparse (n, n, nz); + for (octave_idx_type i = 0; i < coff+1; i++) + d.xcidx (i) = 0; + for (octave_idx_type j = 0; j < nnc; j++) + { + for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) + { + d.xdata (i) = data (i); + d.xridx (i) = j + roff; + } + d.xcidx (j + coff + 1) = cidx(j+1); + } + for (octave_idx_type i = nnc + coff + 1; i < n + 1; i++) + d.xcidx (i) = nz; + } + else + { + octave_idx_type n = nnr + std::abs (k); + octave_idx_type nz = nzmax (); + octave_idx_type ii = 0; + octave_idx_type ir = ridx(0); + d = Sparse (n, n, nz); + for (octave_idx_type i = 0; i < coff+1; i++) + d.xcidx (i) = 0; + for (octave_idx_type i = 0; i < nnr; i++) + { + if (ir == i) + { + d.xdata (ii) = data (ii); + d.xridx (ii++) = ir + roff; + if (ii != nz) + ir = ridx (ii); + } + d.xcidx (i + coff + 1) = ii; + } + for (octave_idx_type i = nnr + coff + 1; i < n+1; i++) + d.xcidx (i) = nz; + } + } + + return d; +} + // FIXME // Unfortunately numel can overflow for very large but very sparse matrices. // For now just flag an error when this happens. diff --git a/liboctave/Sparse.h b/liboctave/Sparse.h --- a/liboctave/Sparse.h +++ b/liboctave/Sparse.h @@ -522,6 +522,8 @@ Sparse sort (Array &sidx, octave_idx_type dim = 0, sortmode mode = ASCENDING) const; + Sparse diag (octave_idx_type k = 0) const; + template Sparse map (F fcn) const diff --git a/liboctave/boolMatrix.cc b/liboctave/boolMatrix.cc --- a/liboctave/boolMatrix.cc +++ b/liboctave/boolMatrix.cc @@ -78,50 +78,9 @@ // other operations boolMatrix -boolMatrix::diag (void) const -{ - return diag (0); -} - -boolMatrix boolMatrix::diag (octave_idx_type k) const { - octave_idx_type nnr = rows (); - octave_idx_type nnc = cols (); - if (k > 0) - nnc -= k; - else if (k < 0) - nnr += k; - - boolMatrix d; - - if (nnr > 0 && nnc > 0) - { - octave_idx_type ndiag = (nnr < nnc) ? nnr : nnc; - - d.resize (ndiag, 1); - - if (k > 0) - { - for (octave_idx_type i = 0; i < ndiag; i++) - d.xelem (i) = elem (i, i+k); - } - else if (k < 0) - { - for (octave_idx_type i = 0; i < ndiag; i++) - d.xelem (i) = elem (i-k, i); - } - else - { - for (octave_idx_type i = 0; i < ndiag; i++) - d.xelem (i) = elem (i, i); - } - } - else - (*current_liboctave_error_handler) - ("diag: requested diagonal out of range"); - - return d; + return Array2::diag (k); } // FIXME Do these really belong here? Maybe they should be diff --git a/liboctave/boolMatrix.h b/liboctave/boolMatrix.h --- a/liboctave/boolMatrix.h +++ b/liboctave/boolMatrix.h @@ -64,8 +64,7 @@ // other operations - boolMatrix diag (void) const; - boolMatrix diag (octave_idx_type k) const; + boolMatrix diag (octave_idx_type k = 0) const; boolMatrix all (int dim = -1) const; boolMatrix any (int dim = -1) const; diff --git a/liboctave/boolNDArray.cc b/liboctave/boolNDArray.cc --- a/liboctave/boolNDArray.cc +++ b/liboctave/boolNDArray.cc @@ -129,6 +129,12 @@ return ::compute_index (ra_idx, dimensions); } +boolNDArray +boolNDArray::diag (octave_idx_type k) const +{ + return ArrayN::diag (k); +} + NDND_BOOL_OPS (boolNDArray, boolNDArray, false) NDND_CMP_OPS (boolNDArray, , boolNDArray, ) diff --git a/liboctave/boolNDArray.h b/liboctave/boolNDArray.h --- a/liboctave/boolNDArray.h +++ b/liboctave/boolNDArray.h @@ -109,6 +109,8 @@ return retval; } + boolNDArray diag (octave_idx_type k = 0) const; + private: boolNDArray (bool *d, dim_vector& dv) : ArrayN (d, dv) { } diff --git a/liboctave/boolSparse.cc b/liboctave/boolSparse.cc --- a/liboctave/boolSparse.cc +++ b/liboctave/boolSparse.cc @@ -143,88 +143,7 @@ SparseBoolMatrix SparseBoolMatrix::diag (octave_idx_type k) const { - octave_idx_type nnr = rows (); - octave_idx_type nnc = cols (); - - if (k > 0) - nnc -= k; - else if (k < 0) - nnr += k; - - SparseBoolMatrix d; - - if (nnr > 0 && nnc > 0) - { - octave_idx_type ndiag = (nnr < nnc) ? nnr : nnc; - - // Count the number of non-zero elements - octave_idx_type nel = 0; - if (k > 0) - { - for (octave_idx_type i = 0; i < ndiag; i++) - if (elem (i, i+k) != 0.) - nel++; - } - else if ( k < 0) - { - for (octave_idx_type i = 0; i < ndiag; i++) - if (elem (i-k, i) != 0.) - nel++; - } - else - { - for (octave_idx_type i = 0; i < ndiag; i++) - if (elem (i, i) != 0.) - nel++; - } - - d = SparseBoolMatrix (ndiag, 1, nel); - d.xcidx (0) = 0; - d.xcidx (1) = nel; - - octave_idx_type ii = 0; - if (k > 0) - { - for (octave_idx_type i = 0; i < ndiag; i++) - { - bool tmp = elem (i, i+k); - if (tmp != 0.) - { - d.xdata (ii) = tmp; - d.xridx (ii++) = i; - } - } - } - else if ( k < 0) - { - for (octave_idx_type i = 0; i < ndiag; i++) - { - bool tmp = elem (i-k, i); - if (tmp != 0.) - { - d.xdata (ii) = tmp; - d.xridx (ii++) = i; - } - } - } - else - { - for (octave_idx_type i = 0; i < ndiag; i++) - { - bool tmp = elem (i, i); - if (tmp != 0.) - { - d.xdata (ii) = tmp; - d.xridx (ii++) = i; - } - } - } - } - else - (*current_liboctave_error_handler) - ("diag: requested diagonal out of range"); - - return d; + return Sparse::diag (k); } boolMatrix diff --git a/liboctave/chMatrix.cc b/liboctave/chMatrix.cc --- a/liboctave/chMatrix.cc +++ b/liboctave/chMatrix.cc @@ -192,50 +192,9 @@ } charMatrix -charMatrix::diag (void) const -{ - return diag (0); -} - -charMatrix charMatrix::diag (octave_idx_type k) const { - octave_idx_type nnr = rows (); - octave_idx_type nnc = cols (); - if (k > 0) - nnc -= k; - else if (k < 0) - nnr += k; - - charMatrix d; - - if (nnr > 0 && nnc > 0) - { - octave_idx_type ndiag = (nnr < nnc) ? nnr : nnc; - - d.resize (ndiag, 1); - - if (k > 0) - { - for (octave_idx_type i = 0; i < ndiag; i++) - d.xelem (i) = elem (i, i+k); - } - else if (k < 0) - { - for (octave_idx_type i = 0; i < ndiag; i++) - d.xelem (i) = elem (i-k, i); - } - else - { - for (octave_idx_type i = 0; i < ndiag; i++) - d.xelem (i) = elem (i, i); - } - } - else - (*current_liboctave_error_handler) - ("diag: requested diagonal out of range"); - - return d; + return MArray2::diag (k); } // FIXME Do these really belong here? Maybe they should be diff --git a/liboctave/chMatrix.h b/liboctave/chMatrix.h --- a/liboctave/chMatrix.h +++ b/liboctave/chMatrix.h @@ -73,8 +73,7 @@ charMatrix extract (octave_idx_type r1, octave_idx_type c1, octave_idx_type r2, octave_idx_type c2) const; - charMatrix diag (void) const; - charMatrix diag (octave_idx_type k) const; + charMatrix diag (octave_idx_type k = 0) const; boolMatrix all (int dim = -1) const; boolMatrix any (int dim = -1) const; diff --git a/liboctave/chNDArray.cc b/liboctave/chNDArray.cc --- a/liboctave/chNDArray.cc +++ b/liboctave/chNDArray.cc @@ -145,6 +145,12 @@ return ::compute_index (ra_idx, dimensions); } +charNDArray +charNDArray::diag (octave_idx_type k) const +{ + return MArrayN::diag (k); +} + boolNDArray charNDArray::bmap (mapper fcn) const { diff --git a/liboctave/chNDArray.h b/liboctave/chNDArray.h --- a/liboctave/chNDArray.h +++ b/liboctave/chNDArray.h @@ -89,6 +89,8 @@ static char resize_fill_value (void) { return '\0'; } + charNDArray diag (octave_idx_type k = 0) const; + typedef int (*mapper) (int); boolNDArray bmap (mapper fcn) const; NDArray dmap (mapper fcn) const; diff --git a/liboctave/dDiagMatrix.cc b/liboctave/dDiagMatrix.cc --- a/liboctave/dDiagMatrix.cc +++ b/liboctave/dDiagMatrix.cc @@ -347,18 +347,13 @@ // other operations ColumnVector -DiagMatrix::diag (void) const -{ - return diag (0); -} - -// Could be optimized... - -ColumnVector DiagMatrix::diag (octave_idx_type k) const { octave_idx_type nnr = rows (); octave_idx_type nnc = cols (); + + if (nnr == 0 || nnc == 0) + if (k > 0) nnc -= k; else if (k < 0) diff --git a/liboctave/dDiagMatrix.h b/liboctave/dDiagMatrix.h --- a/liboctave/dDiagMatrix.h +++ b/liboctave/dDiagMatrix.h @@ -92,8 +92,7 @@ // other operations - ColumnVector diag (void) const; - ColumnVector diag (octave_idx_type k) const; + ColumnVector diag (octave_idx_type k = 0) const; // i/o diff --git a/liboctave/dMatrix.cc b/liboctave/dMatrix.cc --- a/liboctave/dMatrix.cc +++ b/liboctave/dMatrix.cc @@ -2843,51 +2843,10 @@ return retval; } -ColumnVector -Matrix::diag (void) const -{ - return diag (0); -} - -ColumnVector +Matrix Matrix::diag (octave_idx_type k) const { - octave_idx_type nnr = rows (); - octave_idx_type nnc = cols (); - if (k > 0) - nnc -= k; - else if (k < 0) - nnr += k; - - ColumnVector d; - - if (nnr > 0 && nnc > 0) - { - octave_idx_type ndiag = (nnr < nnc) ? nnr : nnc; - - d.resize (ndiag); - - if (k > 0) - { - for (octave_idx_type i = 0; i < ndiag; i++) - d.elem (i) = elem (i, i+k); - } - else if (k < 0) - { - for (octave_idx_type i = 0; i < ndiag; i++) - d.elem (i) = elem (i-k, i); - } - else - { - for (octave_idx_type i = 0; i < ndiag; i++) - d.elem (i) = elem (i, i); - } - } - else - (*current_liboctave_error_handler) - ("diag: requested diagonal out of range"); - - return d; + return MArray2::diag (k); } ColumnVector diff --git a/liboctave/dMatrix.h b/liboctave/dMatrix.h --- a/liboctave/dMatrix.h +++ b/liboctave/dMatrix.h @@ -291,8 +291,7 @@ Matrix sumsq (int dim = -1) const; Matrix abs (void) const; - ColumnVector diag (void) const; - ColumnVector diag (octave_idx_type k) const; + Matrix diag (octave_idx_type k = 0) const; ColumnVector row_min (void) const; ColumnVector row_max (void) const; diff --git a/liboctave/dNDArray.cc b/liboctave/dNDArray.cc --- a/liboctave/dNDArray.cc +++ b/liboctave/dNDArray.cc @@ -968,6 +968,12 @@ } NDArray +NDArray::diag (octave_idx_type k) const +{ + return MArrayN::diag (k); +} + +NDArray NDArray::map (dmapper fcn) const { return MArrayN::map (func_ptr (fcn)); diff --git a/liboctave/dNDArray.h b/liboctave/dNDArray.h --- a/liboctave/dNDArray.h +++ b/liboctave/dNDArray.h @@ -124,6 +124,8 @@ static double resize_fill_value (void) { return 0; } + NDArray diag (octave_idx_type k = 0) const; + typedef double (*dmapper) (double); typedef Complex (*cmapper) (const Complex&); typedef bool (*bmapper) (double); diff --git a/liboctave/dSparse.cc b/liboctave/dSparse.cc --- a/liboctave/dSparse.cc +++ b/liboctave/dSparse.cc @@ -7437,88 +7437,7 @@ SparseMatrix SparseMatrix::diag (octave_idx_type k) const { - octave_idx_type nnr = rows (); - octave_idx_type nnc = cols (); - - if (k > 0) - nnc -= k; - else if (k < 0) - nnr += k; - - SparseMatrix d; - - if (nnr > 0 && nnc > 0) - { - octave_idx_type ndiag = (nnr < nnc) ? nnr : nnc; - - // Count the number of non-zero elements - octave_idx_type nel = 0; - if (k > 0) - { - for (octave_idx_type i = 0; i < ndiag; i++) - if (elem (i, i+k) != 0.) - nel++; - } - else if ( k < 0) - { - for (octave_idx_type i = 0; i < ndiag; i++) - if (elem (i-k, i) != 0.) - nel++; - } - else - { - for (octave_idx_type i = 0; i < ndiag; i++) - if (elem (i, i) != 0.) - nel++; - } - - d = SparseMatrix (ndiag, 1, nel); - d.xcidx (0) = 0; - d.xcidx (1) = nel; - - octave_idx_type ii = 0; - if (k > 0) - { - for (octave_idx_type i = 0; i < ndiag; i++) - { - double tmp = elem (i, i+k); - if (tmp != 0.) - { - d.xdata (ii) = tmp; - d.xridx (ii++) = i; - } - } - } - else if ( k < 0) - { - for (octave_idx_type i = 0; i < ndiag; i++) - { - double tmp = elem (i-k, i); - if (tmp != 0.) - { - d.xdata (ii) = tmp; - d.xridx (ii++) = i; - } - } - } - else - { - for (octave_idx_type i = 0; i < ndiag; i++) - { - double tmp = elem (i, i); - if (tmp != 0.) - { - d.xdata (ii) = tmp; - d.xridx (ii++) = i; - } - } - } - } - else - (*current_liboctave_error_handler) - ("diag: requested diagonal out of range"); - - return d; + return MSparse::diag (k); } Matrix diff --git a/liboctave/intNDArray.cc b/liboctave/intNDArray.cc --- a/liboctave/intNDArray.cc +++ b/liboctave/intNDArray.cc @@ -60,66 +60,11 @@ return false; } - -template -intNDArray -intNDArray::diag (void) const -{ - return diag (0); -} - template intNDArray intNDArray::diag (octave_idx_type k) const { - dim_vector dv = this->dims (); - octave_idx_type nd = dv.length (); - - if (nd > 2) - { - (*current_liboctave_error_handler) ("Matrix must be 2-dimensional"); - return intNDArray(); - } - else - { - octave_idx_type nnr = dv (0); - octave_idx_type nnc = dv (1); - - if (k > 0) - nnc -= k; - else if (k < 0) - nnr += k; - - intNDArray d; - - if (nnr > 0 && nnc > 0) - { - octave_idx_type ndiag = (nnr < nnc) ? nnr : nnc; - - d.resize (dim_vector (ndiag, 1)); - - if (k > 0) - { - for (octave_idx_type i = 0; i < ndiag; i++) - d.xelem (i) = this->elem (i, i+k); - } - else if (k < 0) - { - for (octave_idx_type i = 0; i < ndiag; i++) - d.xelem (i) = this->elem (i-k, i); - } - else - { - for (octave_idx_type i = 0; i < ndiag; i++) - d.xelem (i) = this->elem (i, i); - } - } - else - (*current_liboctave_error_handler) - ("diag: requested diagonal out of range"); - - return d; - } + return MArrayN::diag (k); } // FIXME -- this is not quite the right thing. diff --git a/liboctave/intNDArray.h b/liboctave/intNDArray.h --- a/liboctave/intNDArray.h +++ b/liboctave/intNDArray.h @@ -65,8 +65,7 @@ bool any_element_not_one_or_zero (void) const; - intNDArray diag (void) const; - intNDArray diag (octave_idx_type k) const; + intNDArray diag (octave_idx_type k = 0) const; // FIXME -- this is not quite the right thing. diff --git a/src/Cell.cc b/src/Cell.cc --- a/src/Cell.cc +++ b/src/Cell.cc @@ -239,49 +239,9 @@ } Cell -Cell::diag (void) const -{ - return diag (0); -} - -Cell Cell::diag (octave_idx_type k) const { - octave_idx_type nnr = rows (); - octave_idx_type nnc = cols (); - if (k > 0) - nnc -= k; - else if (k < 0) - nnr += k; - - Cell d; - - if (nnr > 0 && nnc > 0) - { - octave_idx_type ndiag = (nnr < nnc) ? nnr : nnc; - - d.resize (dim_vector (ndiag, 1)); - - if (k > 0) - { - for (octave_idx_type i = 0; i < ndiag; i++) - d.elem (i) = elem (i, i+k); - } - else if (k < 0) - { - for (octave_idx_type i = 0; i < ndiag; i++) - d.elem (i) = elem (i-k, i); - } - else - { - for (octave_idx_type i = 0; i < ndiag; i++) - d.elem (i) = elem (i, i); - } - } - else - error ("diag: requested diagonal out of range"); - - return d; + return ArrayN::diag (k); } /* diff --git a/src/Cell.h b/src/Cell.h --- a/src/Cell.h +++ b/src/Cell.h @@ -115,8 +115,7 @@ static octave_value resize_fill_value (void) { return Matrix (); } - Cell diag (void) const; - Cell diag (octave_idx_type k) const; + Cell diag (octave_idx_type k = 0) const; Cell xisalnum (void) const { return map (&octave_value::xisalnum); } Cell xisalpha (void) const { return map (&octave_value::xisalpha); } diff --git a/src/ChangeLog b/src/ChangeLog --- a/src/ChangeLog +++ b/src/ChangeLog @@ -5,6 +5,29 @@ 2008-03-20 David Bateman + * Cell.cc (Cell Cell::diag (void) const): delete. + (Cell Cell::diag (octave__idx_type) const):Rewrite in terms of + template classes function. + * Cell.h (Cell diag (void) const): delete. + + * ov.h (octave_value diag (octave_idx_type) const): New method. + * ov-base.h (virtual octave_value diag (octave_idx_type) const): + New virtual method. + * ov-base.cc (octave_value octave_base_value::diag + (octave_idx_type) const): New default method. + + * ov-base-mat.h (octave_value diag (octave_idx_type) const): New + method. + * ov-base-sparse.h (octave_value diag (octave_idx_type) const): New + method. + * ov-base-scalar.h (octave_value diag (octave_idx_type) const): New + method. + * ov-range.h (octave_value diag (octave_idx_type) const): New + method. + + * data.cc (make_diag, make_spdiag): Delete. + (Fdiag): Rewrite in terms of octave_value diag function. + * data.cc (static octave_value make_diag (const Cell&, octave_idx_type)): New instantiation of template function. (static octave_value make_diag (const octave_value&, diff --git a/src/data.cc b/src/data.cc --- a/src/data.cc +++ b/src/data.cc @@ -883,305 +883,6 @@ DATA_REDUCTION (cumsum); } -template -static octave_value -make_diag (const T& v, octave_idx_type k) -{ - octave_value retval; - dim_vector dv = v.dims (); - octave_idx_type nd = dv.length (); - - if (nd > 2) - error ("diag: expecting 2-dimensional matrix"); - else - { - octave_idx_type nr = dv (0); - octave_idx_type nc = dv (1); - - if (nr == 0 || nc == 0) - retval = T (); - else if (nr != 1 && nc != 1) - retval = v.diag (k); - else - { - octave_idx_type roff = 0; - octave_idx_type coff = 0; - if (k > 0) - { - roff = 0; - coff = k; - } - else if (k < 0) - { - roff = -k; - coff = 0; - } - - if (nr == 1) - { - octave_idx_type n = nc + std::abs (k); - T m (dim_vector (n, n), T::resize_fill_value ()); - - for (octave_idx_type i = 0; i < nc; i++) - m (i+roff, i+coff) = v (0, i); - retval = m; - } - else - { - octave_idx_type n = nr + std::abs (k); - T m (dim_vector (n, n), T::resize_fill_value ()); - for (octave_idx_type i = 0; i < nr; i++) - m (i+roff, i+coff) = v (i, 0); - retval = m; - } - } - } - - return retval; -} - -#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL) -static octave_value -make_diag (const Matrix& v, octave_idx_type k); - -static octave_value -make_diag (const ComplexMatrix& v, octave_idx_type k); - -static octave_value -make_diag (const charMatrix& v, octave_idx_type k); - -static octave_value -make_diag (const boolMatrix& v, octave_idx_type k); - -static octave_value -make_diag (const int8NDArray& v, octave_idx_type k); - -static octave_value -make_diag (const int16NDArray& v, octave_idx_type k); - -static octave_value -make_diag (const int32NDArray& v, octave_idx_type k); - -static octave_value -make_diag (const int64NDArray& v, octave_idx_type k); - -static octave_value -make_diag (const uint8NDArray& v, octave_idx_type k); - -static octave_value -make_diag (const uint16NDArray& v, octave_idx_type k); - -static octave_value -make_diag (const uint32NDArray& v, octave_idx_type k); - -static octave_value -make_diag (const uint64NDArray& v, octave_idx_type k); - -static octave_value -make_diag (const Cell& v, octave_idx_type k); -#endif - -template -static octave_value -make_spdiag (const T& v, octave_idx_type k) -{ - octave_value retval; - dim_vector dv = v.dims (); - octave_idx_type nr = dv (0); - octave_idx_type nc = dv (1); - - if (nr == 0 || nc == 0) - retval = T (); - else if (nr != 1 && nc != 1) - retval = v.diag (k); - else - { - octave_idx_type roff = 0; - octave_idx_type coff = 0; - if (k > 0) - { - roff = 0; - coff = k; - } - else if (k < 0) - { - roff = -k; - coff = 0; - } - - if (nr == 1) - { - octave_idx_type n = nc + std::abs (k); - octave_idx_type nz = v.nzmax (); - T r (n, n, nz); - for (octave_idx_type i = 0; i < coff+1; i++) - r.xcidx (i) = 0; - for (octave_idx_type j = 0; j < nc; j++) - { - for (octave_idx_type i = v.cidx(j); i < v.cidx(j+1); i++) - { - r.xdata (i) = v.data (i); - r.xridx (i) = j + roff; - } - r.xcidx (j+coff+1) = v.cidx(j+1); - } - for (octave_idx_type i = nc+coff+1; i < n+1; i++) - r.xcidx (i) = nz; - retval = r; - } - else - { - octave_idx_type n = nr + std::abs (k); - octave_idx_type nz = v.nzmax (); - octave_idx_type ii = 0; - octave_idx_type ir = v.ridx(0); - T r (n, n, nz); - for (octave_idx_type i = 0; i < coff+1; i++) - r.xcidx (i) = 0; - for (octave_idx_type i = 0; i < nr; i++) - { - if (ir == i) - { - r.xdata (ii) = v.data (ii); - r.xridx (ii++) = ir + roff; - if (ii != nz) - ir = v.ridx (ii); - } - r.xcidx (i+coff+1) = ii; - } - for (octave_idx_type i = nr+coff+1; i < n+1; i++) - r.xcidx (i) = nz; - retval = r; - } - } - - return retval; -} - -#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL) -static octave_value -make_spdiag (const SparseMatrix& v, octave_idx_type k); - -static octave_value -make_spdiag (const SparseComplexMatrix& v, octave_idx_type k); - -static octave_value -make_spdiag (const SparseBoolMatrix& v, octave_idx_type k); -#endif - -static octave_value -make_diag (const octave_value& a, octave_idx_type k) -{ - octave_value retval; - std::string result_type = a.class_name (); - - if (result_type == "double") - { - if (a.is_sparse_type ()) - { - if (a.is_real_type ()) - { - SparseMatrix m = a.sparse_matrix_value (); - if (!error_state) - retval = make_spdiag (m, k); - } - else - { - SparseComplexMatrix m = a.sparse_complex_matrix_value (); - if (!error_state) - retval = make_spdiag (m, k); - } - } - else - { - if (a.is_real_type ()) - { - Matrix m = a.matrix_value (); - if (!error_state) - retval = make_diag (m, k); - } - else - { - ComplexMatrix m = a.complex_matrix_value (); - if (!error_state) - retval = make_diag (m, k); - } - } - } -#if 0 - else if (result_type == "single") - retval = make_diag (a.single_array_value (), k); -#endif - else if (result_type == "char") - { - charMatrix m = a.char_matrix_value (); - if (!error_state) - { - retval = make_diag (m, k); - if (a.is_sq_string ()) - retval = octave_value (retval.char_array_value (), true, '\''); - } - } - else if (result_type == "logical") - { - if (a.is_sparse_type ()) - { - SparseBoolMatrix m = a.sparse_bool_matrix_value (); - if (!error_state) - retval = make_spdiag (m, k); - } - else - { - boolMatrix m = a.bool_matrix_value (); - if (!error_state) - retval = make_diag (m, k); - } - } - else if (result_type == "int8") - retval = make_diag (a.int8_array_value (), k); - else if (result_type == "int16") - retval = make_diag (a.int16_array_value (), k); - else if (result_type == "int32") - retval = make_diag (a.int32_array_value (), k); - else if (result_type == "int64") - retval = make_diag (a.int64_array_value (), k); - else if (result_type == "uint8") - retval = make_diag (a.uint8_array_value (), k); - else if (result_type == "uint16") - retval = make_diag (a.uint16_array_value (), k); - else if (result_type == "uint32") - retval = make_diag (a.uint32_array_value (), k); - else if (result_type == "uint64") - retval = make_diag (a.uint64_array_value (), k); - else if (result_type == "cell") - retval = make_diag (a.cell_value (), k); - else - gripe_wrong_type_arg ("diag", a); - - return retval; -} - -static octave_value -make_diag (const octave_value& arg) -{ - return make_diag (arg, 0); -} - -static octave_value -make_diag (const octave_value& a, const octave_value& b) -{ - octave_value retval; - - octave_idx_type k = b.int_value (); - - if (error_state) - error ("diag: invalid second argument"); - else - retval = make_diag (a, k); - - return retval; -} - DEFUN (diag, args, , "-*- texinfo -*-\n\ @deftypefn {Built-in Function} {} diag (@var{v}, @var{k})\n\ @@ -1211,9 +912,16 @@ int nargin = args.length (); if (nargin == 1 && args(0).is_defined ()) - retval = make_diag (args(0)); + retval = args(0).diag(); else if (nargin == 2 && args(0).is_defined () && args(1).is_defined ()) - retval = make_diag (args(0), args(1)); + { + octave_idx_type k = args(1).int_value (); + + if (error_state) + error ("diag: invalid second argument"); + else + retval = args(0).diag(k); + } else print_usage (); diff --git a/src/ov-base-mat.h b/src/ov-base-mat.h --- a/src/ov-base-mat.h +++ b/src/ov-base-mat.h @@ -110,6 +110,9 @@ MatrixType matrix_type (const MatrixType& _typ) const { MatrixType ret = typ; typ = _typ; return ret; } + octave_value diag (octave_idx_type k = 0) const + { return octave_value (matrix.diag (k)); } + octave_value sort (octave_idx_type dim = 0, sortmode mode = ASCENDING) const { return octave_value (matrix.sort (dim, mode)); } octave_value sort (Array &sidx, octave_idx_type dim = 0, diff --git a/src/ov-base-scalar.h b/src/ov-base-scalar.h --- a/src/ov-base-scalar.h +++ b/src/ov-base-scalar.h @@ -93,6 +93,9 @@ octave_value any (int = 0) const { return (scalar != ST ()); } + octave_value diag (octave_idx_type k = 0) const + { return octave_value (matrix_value (). diag (k)); } + octave_value sort (octave_idx_type, sortmode) const { return octave_value (scalar); } octave_value sort (Array &sidx, octave_idx_type, diff --git a/src/ov-base-sparse.h b/src/ov-base-sparse.h --- a/src/ov-base-sparse.h +++ b/src/ov-base-sparse.h @@ -116,6 +116,9 @@ octave_value all (int dim = 0) const { return matrix.all (dim); } octave_value any (int dim = 0) const { return matrix.any (dim); } + octave_value diag (octave_idx_type k = 0) const + { return octave_value (matrix.diag (k)); } + octave_value sort (octave_idx_type dim = 0, sortmode mode = ASCENDING) const { return octave_value (matrix.sort (dim, mode)); } octave_value sort (Array &sidx, octave_idx_type dim = 0, diff --git a/src/ov-base.cc b/src/ov-base.cc --- a/src/ov-base.cc +++ b/src/ov-base.cc @@ -879,6 +879,14 @@ } octave_value +octave_base_value::diag (octave_idx_type) const +{ + gripe_wrong_type_arg ("octave_base_value::diag ()", type_name ()); + + return octave_value(); +} + +octave_value octave_base_value::sort (octave_idx_type, sortmode) const { gripe_wrong_type_arg ("octave_base_value::sort ()", type_name ()); diff --git a/src/ov-base.h b/src/ov-base.h --- a/src/ov-base.h +++ b/src/ov-base.h @@ -458,6 +458,8 @@ virtual mxArray *as_mxArray (void) const; + virtual octave_value diag (octave_idx_type k = 0) const; + virtual octave_value sort (octave_idx_type dim = 0, sortmode mode = ASCENDING) const; virtual octave_value sort (Array &sidx, diff --git a/src/ov-range.h b/src/ov-range.h --- a/src/ov-range.h +++ b/src/ov-range.h @@ -132,6 +132,9 @@ octave_value any (int dim = 0) const; + octave_value diag (octave_idx_type k = 0) const + { return octave_value (range.diag (k)); } + octave_value sort (octave_idx_type dim = 0, sortmode mode = ASCENDING) const { return range.sort (dim, mode); } diff --git a/src/ov.h b/src/ov.h --- a/src/ov.h +++ b/src/ov.h @@ -868,6 +868,9 @@ mxArray *as_mxArray (void) const { return rep->as_mxArray (); } + octave_value diag (octave_idx_type k = 0) const + { return rep->diag (k); } + octave_value sort (octave_idx_type dim = 0, sortmode mode = ASCENDING) const { return rep->sort (dim, mode); } octave_value sort (Array &sidx, octave_idx_type dim = 0,