# HG changeset patch
# User Christof Kaufmann <christofkaufmann.dev@gmail.com>
# Date 1643126838 -3600
#      Tue Jan 25 17:07:18 2022 +0100
# Node ID b7c7d2ef8c75b80f661e38ceb4d52027b5e02eab
# Parent  20b05c5c41a1599e1fc58f756b5ad14b07966f4d
Fix interp1 and interp2 for complex matrix (bug #61863, bug #61903)

diff --git a/libinterp/corefcn/__pchip_deriv__.cc b/libinterp/corefcn/__pchip_deriv__.cc
--- a/libinterp/corefcn/__pchip_deriv__.cc
+++ b/libinterp/corefcn/__pchip_deriv__.cc
@@ -57,79 +57,167 @@ Undocumented internal function.
       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 @@ Undocumented internal function.
 }
 
 /*
-## 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
diff --git a/scripts/general/interp1.m b/scripts/general/interp1.m
--- a/scripts/general/interp1.m
+++ b/scripts/general/interp1.m
@@ -123,7 +123,7 @@ function yi = interp1 (x, y, varargin)
   endif
 
   method = "linear";
-  extrap = NA;
+  extrap = [];
   xi = [];
   ispp = false;
   have_xi = false;
@@ -167,6 +167,14 @@ function yi = interp1 (x, y, varargin)
     endif
   endif
 
+  if (isempty (extrap))
+    if (iscomplex (y))
+      extrap = complex (NA, NA);
+    else
+      extrap = NA;
+    endif
+  endif
+
   ## reshape matrices for convenience
   x = x(:);
   nx = rows (x);
@@ -503,6 +511,8 @@ endfunction
 %!assert (isempty (interp1 (xp,yp,[],style)))
 %!assert (interp1 (xp,[yp',yp'],xi(:),style),...
 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
+%!assert <*61903> (interp1 (xp, yp + 1i * yp.^2, xi, style),...
+%!                 interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style))
 %!assert (interp1 (xp,yp,xi,style),...
 %!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
 %!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
@@ -520,6 +530,8 @@ endfunction
 %!assert (isempty (interp1 (xp,yp,[],style)))
 %!assert (interp1 (xp,[yp',yp'],xi(:),style),...
 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
+%!assert <*61903> (interp1 (xp, yp + 1i * yp.^2, xi, style),...
+%!                 interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style))
 %!assert (interp1 (xp,yp,xi,style),...
 %!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
 %!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
@@ -540,6 +552,8 @@ endfunction
 %!assert (isempty (interp1 (xp,yp,[],style)))
 %!assert (interp1 (xp,[yp',yp'],xi(:),style),...
 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
+%!assert <*61903> (interp1 (xp, yp + 1i * yp.^2, xi, style),...
+%!                 interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style))
 ## This test is expected to fail, so commented out.
 ## "previous" and "next" options are not symmetric w.r.t to flipping xp,yp
 #%!assert (interp1 (xp,yp,xi,style),...
@@ -559,6 +573,8 @@ endfunction
 %!assert (isempty (interp1 (xp,yp,[],style)))
 %!assert (interp1 (xp,[yp',yp'],xi(:),style),...
 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
+%!assert <*61903> (interp1 (xp, yp + 1i * yp.^2, xi, style),...
+%!                 interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style))
 #%!assert (interp1 (xp,yp,xi,style),...
 #%!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
 %!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
@@ -579,6 +595,8 @@ endfunction
 %!assert (isempty (interp1 (xp,yp,[],style)))
 %!assert (interp1 (xp,[yp',yp'],xi(:),style),...
 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
+%!assert <*61903> (interp1 (xp, yp + 1i * yp.^2, xi, style),...
+%!                 interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style))
 #%!assert (interp1 (xp,yp,xi,style),...
 #%!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
 %!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
@@ -596,6 +614,8 @@ endfunction
 %!assert (isempty (interp1 (xp,yp,[],style)))
 %!assert (interp1 (xp,[yp',yp'],xi(:),style),...
 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
+%!assert <*61903> (interp1 (xp, yp + 1i * yp.^2, xi, style),...
+%!                 interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style))
 #%!assert (interp1 (xp,yp,xi,style),...
 #%!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
 %!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
@@ -616,6 +636,8 @@ endfunction
 %!assert (isempty (interp1 (xp,yp,[],style)))
 %!assert (interp1 (xp,[yp',yp'],xi(:),style),...
 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
+%!assert <*61903> (interp1 (xp, yp + 1i * yp.^2, xi, style),...
+%!                 interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style))
 %!assert (interp1 (xp,yp,xi,style),...
 %!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
 %!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
@@ -633,6 +655,8 @@ endfunction
 %!assert (isempty (interp1 (xp,yp,[],style)))
 %!assert (interp1 (xp,[yp',yp'],xi(:),style),...
 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
+%!assert <*61903> (interp1 (xp, yp + 1i * yp.^2, xi, style),...
+%!                 interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style))
 %!assert (interp1 (xp,yp,xi,style),...
 %!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
 %!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
@@ -655,6 +679,8 @@ endfunction
 %!assert (isempty (interp1 (xp,yp,[],style)))
 %!assert (interp1 (xp,[yp',yp'],xi(:),style),...
 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
+%!assert <*61903> (interp1 (xp, yp + 1i * yp.^2, xi, style),...
+%!                 interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style))
 %!assert (interp1 (xp,yp,xi,style),...
 %!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
 %!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
@@ -672,6 +698,8 @@ endfunction
 %!assert (isempty (interp1 (xp,yp,[],style)))
 %!assert (interp1 (xp,[yp',yp'],xi(:),style),...
 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
+%!assert <*61903> (interp1 (xp, yp + 1i * yp.^2, xi, style),...
+%!                 interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style))
 %!assert (interp1 (xp,yp,xi,style),...
 %!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
 %!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
@@ -692,6 +720,8 @@ endfunction
 %!assert (isempty (interp1 (xp,yp,[],style)))
 %!assert (interp1 (xp,[yp',yp'],xi(:),style),...
 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
+%!assert <*61903> (interp1 (xp, yp + 1i * yp.^2, xi, style),...
+%!                 interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style))
 %!assert (interp1 (xp,yp,xi,style),...
 %!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
 %!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
@@ -709,6 +739,8 @@ endfunction
 %!assert (isempty (interp1 (xp,yp,[],style)))
 %!assert (interp1 (xp,[yp',yp'],xi(:),style),...
 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
+%!assert <*61903> (interp1 (xp, yp + 1i * yp.^2, xi, style),...
+%!                 interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style))
 %!assert (interp1 (xp,yp,xi,style),...
 %!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
 %!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
@@ -729,6 +761,8 @@ endfunction
 %!assert (isempty (interp1 (xp,yp,[],style)))
 %!assert (interp1 (xp,[yp',yp'],xi(:),style),...
 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
+%!assert <*61903> (interp1 (xp, yp + 1i * yp.^2, xi, style),...
+%!                 interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style))
 %!assert (interp1 (xp,yp,xi,style),...
 %!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
 %!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
@@ -746,6 +780,8 @@ endfunction
 %!assert (isempty (interp1 (xp,yp,[],style)))
 %!assert (interp1 (xp,[yp',yp'],xi(:),style),...
 %!        [interp1(xp,yp,xi(:),style),interp1(xp,yp,xi(:),style)])
+%!assert <*61903> (interp1 (xp, yp + 1i * yp.^2, xi, style),...
+%!                 interp1 (xp,yp,xi,style) + 1i * interp1 (xp,yp.^2,xi,style))
 %!assert (interp1 (xp,yp,xi,style),...
 %!        interp1 (fliplr (xp),fliplr (yp),xi,style),100*eps)
 %!assert (ppval (interp1 (xp,yp,style,"pp"),xi),
diff --git a/scripts/general/interp2.m b/scripts/general/interp2.m
--- a/scripts/general/interp2.m
+++ b/scripts/general/interp2.m
@@ -257,18 +257,6 @@ function ZI = interp2 (varargin)
       ## 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 @@ endfunction
 
 ## 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);
