Mercurial > octave
changeset 30689:4b367bf5eb16
Fix __pchip_deriv__ for complex matrix (bug #61863, bug #61903)
* libinterp/corefcn/__pchip_deriv__.cc (F__pchip_deriv__): Generalize
implementation for complex valued input.
* scripts/general/interp2.m, scripts/polynomial/pchip.m: Remove work-around.
| author | Christof Kaufmann <christofkaufmann.dev@gmail.com> |
|---|---|
| date | Tue, 25 Jan 2022 17:07:18 +0100 |
| parents | d1fe2cb16d95 |
| children | e2c8c852399e |
| files | libinterp/corefcn/__pchip_deriv__.cc scripts/general/interp2.m scripts/polynomial/pchip.m |
| diffstat | 3 files changed, 166 insertions(+), 66 deletions(-) [+] |
line wrap: on
line diff
--- a/libinterp/corefcn/__pchip_deriv__.cc Wed Jan 26 14:42:09 2022 +0100 +++ b/libinterp/corefcn/__pchip_deriv__.cc Tue Jan 25 17:07:18 2022 +0100 @@ -57,79 +57,167 @@ if (args(0).is_single_type () || args(1).is_single_type ()) { FloatColumnVector xvec (args(0).float_vector_value ()); - FloatMatrix ymat (args(1).float_matrix_value ()); - F77_INT nx = to_f77_int (xvec.numel ()); - if (nx < 2) error ("__pchip_deriv__: X must be at least of length 2"); - octave_idx_type nyr = ymat.rows (); - octave_idx_type nyc = ymat.columns (); + if (args(1).iscomplex ()) + { + FloatComplexMatrix ymat (args(1).float_complex_matrix_value ()); + + octave_idx_type nyr = ymat.rows (); + octave_idx_type nyc = ymat.columns (); - if (nx != (rows ? nyc : nyr)) - error ("__pchip_deriv__: X and Y dimension mismatch"); + if (nx != (rows ? nyc : nyr)) + error ("__pchip_deriv__: X and Y dimension mismatch"); + + FloatComplexMatrix dmat (nyr, nyc); + + F77_INT ierr; + const F77_INT incfd = (rows ? to_f77_int (2*nyr) : 2); + volatile const octave_idx_type inc = (rows ? 2 : 2*nyr); + volatile octave_idx_type k = 0; - FloatMatrix dmat (nyr, nyc); + for (volatile octave_idx_type i = (rows ? nyr : nyc); i > 0; i--) + { + F77_XFCN (pchim, PCHIM, (nx, xvec.data (), + reinterpret_cast<float const*>(ymat.data ()) + k * inc, + reinterpret_cast<float*>(dmat.fortran_vec ()) + k * inc, + incfd, ierr)); - F77_INT ierr; - const F77_INT incfd = (rows ? to_f77_int (nyr) : 1); - volatile const octave_idx_type inc = (rows ? 1 : nyr); - volatile octave_idx_type k = 0; + if (ierr < 0) + error ("__pchip_deriv__: PCHIM failed for real part with ierr = %" + OCTAVE_F77_INT_TYPE_FORMAT, ierr); + + F77_XFCN (pchim, PCHIM, (nx, xvec.data (), + reinterpret_cast<float const*>(ymat.data ()) + 1 + k * inc, + reinterpret_cast<float*>(dmat.fortran_vec ()) + 1 + k * inc, + incfd, ierr)); + + if (ierr < 0) + error ("__pchip_deriv__: PCHIM failed for imaginary part with ierr = %" + OCTAVE_F77_INT_TYPE_FORMAT, ierr); - for (volatile octave_idx_type i = (rows ? nyr : nyc); i > 0; i--) + k++; + } + + retval = dmat; + } + else { - F77_XFCN (pchim, PCHIM, (nx, xvec.data (), - ymat.data () + k * inc, - dmat.fortran_vec () + k * inc, - incfd, ierr)); + FloatMatrix ymat (args(1).float_matrix_value ()); + + octave_idx_type nyr = ymat.rows (); + octave_idx_type nyc = ymat.columns (); + + if (nx != (rows ? nyc : nyr)) + error ("__pchip_deriv__: X and Y dimension mismatch"); + + FloatMatrix dmat (nyr, nyc); - k++; + F77_INT ierr; + const F77_INT incfd = (rows ? to_f77_int (nyr) : 1); + volatile const octave_idx_type inc = (rows ? 1 : nyr); + volatile octave_idx_type k = 0; - if (ierr < 0) - error ("__pchip_deriv__: PCHIM failed with ierr = %" - OCTAVE_F77_INT_TYPE_FORMAT, ierr); + for (volatile octave_idx_type i = (rows ? nyr : nyc); i > 0; i--) + { + F77_XFCN (pchim, PCHIM, (nx, xvec.data (), + ymat.data () + k * inc, + dmat.fortran_vec () + k * inc, + incfd, ierr)); + + k++; + + if (ierr < 0) + error ("__pchip_deriv__: PCHIM failed with ierr = %" + OCTAVE_F77_INT_TYPE_FORMAT, ierr); + } + + retval = dmat; } - - retval = dmat; } else { ColumnVector xvec (args(0).vector_value ()); - Matrix ymat (args(1).matrix_value ()); - F77_INT nx = to_f77_int (xvec.numel ()); - if (nx < 2) error ("__pchip_deriv__: X must be at least of length 2"); - octave_idx_type nyr = ymat.rows (); - octave_idx_type nyc = ymat.columns (); + if (args(1).iscomplex ()) + { + ComplexMatrix ymat (args(1).complex_matrix_value ()); + + octave_idx_type nyr = ymat.rows (); + octave_idx_type nyc = ymat.columns (); - if (nx != (rows ? nyc : nyr)) - error ("__pchip_deriv__: X and Y dimension mismatch"); + if (nx != (rows ? nyc : nyr)) + error ("__pchip_deriv__: X and Y dimension mismatch"); + + ComplexMatrix dmat (nyr, nyc); + + F77_INT ierr; + const F77_INT incfd = (rows ? to_f77_int (2*nyr) : 2); + volatile const octave_idx_type inc = (rows ? 2 : 2*nyr); + volatile octave_idx_type k = 0; - Matrix dmat (nyr, nyc); + for (volatile octave_idx_type i = (rows ? nyr : nyc); i > 0; i--) + { + F77_XFCN (dpchim, DPCHIM, (nx, xvec.data (), + reinterpret_cast<double const*>(ymat.data ()) + k * inc, + reinterpret_cast<double*>(dmat.fortran_vec ()) + k * inc, + incfd, ierr)); - F77_INT ierr; - const F77_INT incfd = (rows ? to_f77_int (nyr) : 1); - volatile const octave_idx_type inc = (rows ? 1 : nyr); - volatile octave_idx_type k = 0; + if (ierr < 0) + error ("__pchip_deriv__: DPCHIM failed for real part with ierr = %" + OCTAVE_F77_INT_TYPE_FORMAT, ierr); + + F77_XFCN (dpchim, DPCHIM, (nx, xvec.data (), + reinterpret_cast<double const*>(ymat.data ()) + 1 + k * inc, + reinterpret_cast<double*>(dmat.fortran_vec ()) + 1 + k * inc, + incfd, ierr)); + + if (ierr < 0) + error ("__pchip_deriv__: DPCHIM failed for imaginary part with ierr = %" + OCTAVE_F77_INT_TYPE_FORMAT, ierr); - for (volatile octave_idx_type i = (rows ? nyr : nyc); i > 0; i--) + k++; + } + + retval = dmat; + } + else { - F77_XFCN (dpchim, DPCHIM, (nx, xvec.data (), - ymat.data () + k * inc, - dmat.fortran_vec () + k * inc, - incfd, ierr)); - k++; + Matrix ymat (args(1).matrix_value ()); + + octave_idx_type nyr = ymat.rows (); + octave_idx_type nyc = ymat.columns (); + + if (nx != (rows ? nyc : nyr)) + error ("__pchip_deriv__: X and Y dimension mismatch"); + + Matrix dmat (nyr, nyc); + + F77_INT ierr; + const F77_INT incfd = (rows ? to_f77_int (nyr) : 1); + volatile const octave_idx_type inc = (rows ? 1 : nyr); + volatile octave_idx_type k = 0; - if (ierr < 0) - error ("__pchip_deriv__: DPCHIM failed with ierr = %" - OCTAVE_F77_INT_TYPE_FORMAT, ierr); + for (volatile octave_idx_type i = (rows ? nyr : nyc); i > 0; i--) + { + F77_XFCN (dpchim, DPCHIM, (nx, xvec.data (), + ymat.data () + k * inc, + dmat.fortran_vec () + k * inc, + incfd, ierr)); + k++; + + if (ierr < 0) + error ("__pchip_deriv__: DPCHIM failed with ierr = %" + OCTAVE_F77_INT_TYPE_FORMAT, ierr); + } + + retval = dmat; } - - retval = dmat; } } @@ -137,8 +225,37 @@ } /* -## No test needed for internal helper function. -%!assert (1) +%!shared x, y +%! x = 0:3; +%! y = x.^2 + 1i * x.^3; + +%!test +%! d_complex = __pchip_deriv__ (x, y, 2); +%! d_real = __pchip_deriv__ (x, real (y), 2); +%! d_imag = __pchip_deriv__ (x, imag (y), 2); +%! assert (real (d_complex), d_real); +%! assert (imag (d_complex), d_imag); + +%!test +%! d_complex = __pchip_deriv__ (x.', y.'); +%! d_real = __pchip_deriv__ (x.', real (y.')); +%! d_imag = __pchip_deriv__ (x.', imag (y.')); +%! assert (real (d_complex), d_real); +%! assert (imag (d_complex), d_imag); + +%!test +%! d_complex = __pchip_deriv__ (single (x), single (y), 2); +%! d_real = __pchip_deriv__ (single (x), real (single (y)), 2); +%! d_imag = __pchip_deriv__ (single (x), imag (single (y)), 2); +%! assert (real (d_complex), d_real); +%! assert (imag (d_complex), d_imag); + +%!test +%! d_complex = __pchip_deriv__ (single (x'), single (y')); +%! d_real = __pchip_deriv__ (single (x'), real (single (y'))); +%! d_imag = __pchip_deriv__ (single (x'), imag (single (y'))); +%! assert (real (d_complex), d_real); +%! assert (imag (d_complex), d_imag); */ OCTAVE_NAMESPACE_END
--- a/scripts/general/interp2.m Wed Jan 26 14:42:09 2022 +0100 +++ b/scripts/general/interp2.m Tue Jan 25 17:07:18 2022 +0100 @@ -257,18 +257,6 @@ ## Compute mixed derivatives row-wise and column-wise. Use the average. DXY = (__pchip_deriv__ (X, DY, 2) + __pchip_deriv__ (Y, DX, 1))/2; - if (iscomplex (Z)) - ## __pchip_deriv__ works only on real part. Do it again for imag part. - ## FIXME: Adapt __pchip_deriv__ to correctly handle complex input. - - ## first order derivatives - DX += 1i * __pchip_deriv__ (X, imag (Z), 2); - DY += 1i * __pchip_deriv__ (Y, imag (Z), 1); - ## Compute mixed derivatives row-wise and column-wise. Use the average. - DXY += 1i * (__pchip_deriv__ (X, imag (DY), 2) - + __pchip_deriv__ (Y, imag (DX), 1))/2; - endif - ## do the bicubic interpolation hx = diff (X); hx = hx(xidx); hy = diff (Y); hy = hy(yidx); @@ -498,7 +486,7 @@ ## Test that interpolating a complex matrix is equivalent to interpolating its ## real and imaginary parts separately. -%!test <61863> +%!test <*61863> %! xi = [2.5, 3.5]; %! yi = [0.5, 1.5]'; %! orig = rand (2, 3) + 1i * rand (2, 3);
--- a/scripts/polynomial/pchip.m Wed Jan 26 14:42:09 2022 +0100 +++ b/scripts/polynomial/pchip.m Tue Jan 25 17:07:18 2022 +0100 @@ -111,11 +111,6 @@ ## Compute derivatives. d = __pchip_deriv__ (x, y, 2); - if (iscomplex (y)) - ## __pchip_deriv__ ignores imaginary part. Do it again for imag part. - ## FIXME: Adapt __pchip_deriv__ to correctly handle complex input. - d += 1i * __pchip_deriv__ (x, imag (y), 2); - endif d1 = d(:, 1:n-1); d2 = d(:, 2:n);