diff --git a/libinterp/corefcn/sparse-xpow.cc b/libinterp/corefcn/sparse-xpow.cc --- a/libinterp/corefcn/sparse-xpow.cc +++ b/libinterp/corefcn/sparse-xpow.cc @@ -403,11 +403,116 @@ elem_xpow (const SparseMatrix& a, double } // -*- 4 -*- + +static inline double +real_pow (double a, double b) +{ + return std::pow (a, b); +} + +static inline Complex +real_to_complex_pow (double a, double b) +{ + return std::pow (Complex (a), Complex (b)); +} + +static inline Complex +complex_real_pow (const Complex& a, double b) +{ + if (xisint (b)) + return std::pow (Complex (a), static_cast (b)); + else + return std::pow (Complex (a), b); +} + +static inline Complex +complex_pow (const Complex& a, const Complex& b) +{ + return std::pow (Complex (a), Complex (b)); +} + +template +RT +do_elem_xpow (const T_LHS& a, const T_RHS& b, POW_FCN pow_fcn) +{ + // Expect the worst, almost all elements are one. This will happen + // if sparsity is high for both matrices and there is little + // correspondence in the locations of the non-zero elements. + + octave_idx_type nr = a.rows (); + octave_idx_type nc = a.cols (); + + RT r (nr, nc, 1.0); + + for (octave_idx_type j = 0; j < nc; j++) + { + octave_idx_type i1 = a.cidx(j); + octave_idx_type e1 = a.cidx(j+1); + octave_idx_type i2 = b.cidx(j); + octave_idx_type e2 = b.cidx(j+1); + + // If i1 == e1, then the LHS has no non-zero elements in this column. + // If i2 == e2, then the RHS has no non-zero elements in this column. + // If neither matrix has elements in this column, then the + // result for all elements is the initial value (1.0) and there + // is no need to execute the following loop. + + while (i1 < e1 || i2 < e2) + { + if (i1 == e1 || (i2 < e2 && a.ridx (i1) > b.ridx (i2))) + { + // LHS has no non-zero elements in this column but RHS + // has at least one OR both columns have elements but + // the row index of the LHS element is greater than the + // row index of the RHS element, so we are looking at + // 0^Y. If Y is not 0.0, then the result is 0.0. + // Otherwise, keep the initial result value (1.0). + + if (b.data (i2) != 0.0) + r.xelem (b.ridx (i2), j) = 0.0; + + i2++; + } + else if (i2 == e2 || b.ridx (i2) > a.ridx (i1)) + { + // RHS has no non-zero elements in this column but LHS + // has at least one OR both columns have elements but + // the row index of the RHS element is greater than the + // row index of the LHS element, so we are looking at + // X^0. Regardless of the value of X, the result is 1.0, + // so keep the initial result value (1.0). + + i1++; + } + else + { + // We have elements for both the RHS and LHS at the same + // row index so we are looking at X^Y. If the value of + // Y is 0.0, then for any value of X, keep the initial + // result (1.0). If the values of X and Y are both + // non-zero, call pow_fcn to compute the result. + // Otherwise, we have 0^Y with a non-zero value of Y so + // the result is 0.0. + + if (b.data (i2) != 0.0) + r.xelem (b.ridx (i2), j) + = (a.data (i1) != 0.0 + ? pow_fcn (a.data (i1), b.data (i2)) : 0.0); + + i1++; + i2++; + } + } + } + + r.maybe_compress (true); + + return r; +} + octave_value elem_xpow (const SparseMatrix& a, const SparseMatrix& b) { - octave_value retval; - octave_idx_type nr = a.rows (); octave_idx_type nc = a.cols (); @@ -437,45 +542,10 @@ elem_xpow (const SparseMatrix& a, const done: - // This is a dumb operator for sparse matrices anyway, and there is - // no sensible way to handle the 0.^0 versus the 0.^x cases. Therefore - // allocate a full matrix filled for the 0.^0 case and shrink it later - // as needed. - if (convert_to_complex) - { - SparseComplexMatrix complex_result (nr, nc, Complex (1.0, 0.0)); - - for (octave_idx_type j = 0; j < nc; j++) - { - for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++) - { - octave_quit (); - complex_result.xelem (a.ridx (i), j) - = std::pow (Complex (a.data (i)), Complex (b(a.ridx (i), j))); - } - } - complex_result.maybe_compress (true); - retval = complex_result; - } + return do_elem_xpow (a, b, real_to_complex_pow); else - { - SparseMatrix result (nr, nc, 1.0); - - for (octave_idx_type j = 0; j < nc; j++) - { - for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++) - { - octave_quit (); - result.xelem (a.ridx (i), j) = std::pow (a.data (i), - b(a.ridx (i), j)); - } - } - result.maybe_compress (true); - retval = result; - } - - return retval; + return do_elem_xpow (a, b, real_pow); } // -*- 5 -*- @@ -522,19 +592,7 @@ elem_xpow (const SparseMatrix& a, const if (nr != b_nr || nc != b_nc) octave::err_nonconformant ("operator .^", nr, nc, b_nr, b_nc); - SparseComplexMatrix result (nr, nc, Complex (1.0, 0.0)); - for (octave_idx_type j = 0; j < nc; j++) - { - for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++) - { - octave_quit (); - result.xelem (a.ridx(i), j) = std::pow (a.data (i), b(a.ridx (i), j)); - } - } - - result.maybe_compress (true); - - return result; + return do_elem_xpow (a, b, complex_pow); } // -*- 7 -*- @@ -662,25 +720,7 @@ elem_xpow (const SparseComplexMatrix& a, if (nr != b_nr || nc != b_nc) octave::err_nonconformant ("operator .^", nr, nc, b_nr, b_nc); - SparseComplexMatrix result (nr, nc, Complex (1.0, 0.0)); - for (octave_idx_type j = 0; j < nc; j++) - { - for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++) - { - octave_quit (); - double btmp = b(a.ridx (i), j); - - if (xisint (btmp)) - result.xelem (a.ridx (i), j) = std::pow (a.data (i), - static_cast (btmp)); - else - result.xelem (a.ridx (i), j) = std::pow (a.data (i), btmp); - } - } - - result.maybe_compress (true); - - return result; + return do_elem_xpow (a, b, complex_real_pow); } // -*- 11 -*- @@ -729,17 +769,5 @@ elem_xpow (const SparseComplexMatrix& a, if (nr != b_nr || nc != b_nc) octave::err_nonconformant ("operator .^", nr, nc, b_nr, b_nc); - SparseComplexMatrix result (nr, nc, Complex (1.0, 0.0)); - for (octave_idx_type j = 0; j < nc; j++) - { - for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++) - { - octave_quit (); - result.xelem (a.ridx (i), j) = std::pow (a.data (i), - b(a.ridx (i), j)); - } - } - result.maybe_compress (true); - - return result; + return do_elem_xpow (a, b, complex_pow); }