GNU Octave  4.0.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
__pchip_deriv__.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2002-2015 Kai Habel
4 Copyright (C) 2008-2009 Jaroslav Hajek
5 
6 This file is part of Octave.
7 
8 Octave is free software; you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by the
10 Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12 
13 Octave is distributed in the hope that it will be useful, but WITHOUT
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 for more details.
17 
18 You should have received a copy of the GNU General Public License
19 along with Octave; see the file COPYING. If not, see
20 <http://www.gnu.org/licenses/>.
21 
22 */
23 
24 #ifdef HAVE_CONFIG_H
25 #include <config.h>
26 #endif
27 
28 #include "defun.h"
29 #include "error.h"
30 #include "gripes.h"
31 #include "oct-obj.h"
32 #include "utils.h"
33 #include "f77-fcn.h"
34 
35 extern "C"
36 {
37  F77_RET_T
38  F77_FUNC (dpchim, DPCHIM) (const octave_idx_type& n, const double *x,
39  const double *f, double *d,
40  const octave_idx_type &incfd,
42 
43  F77_RET_T
44  F77_FUNC (pchim, PCHIM) (const octave_idx_type& n, const float *x,
45  const float *f, float *d,
46  const octave_idx_type& incfd,
48 }
49 
50 // Wrapper for SLATEC/PCHIP function DPCHIM to calculate the derivates
51 // for piecewise polynomials.
52 
53 DEFUN (__pchip_deriv__, args, ,
54  "-*- texinfo -*-\n\
55 @deftypefn {Built-in Function} {} __pchip_deriv__ (@var{x}, @var{y}, @var{dim})\n\
56 Undocumented internal function.\n\
57 @end deftypefn")
58 {
59  octave_value retval;
60  const int nargin = args.length ();
61 
62  bool rows = (nargin == 3 && args(2).uint_value () == 2);
63 
64  if (nargin >= 2)
65  {
66  if (args(0).is_single_type () || args(1).is_single_type ())
67  {
68  FloatColumnVector xvec (args(0).float_vector_value ());
69  FloatMatrix ymat (args(1).float_matrix_value ());
70 
71  octave_idx_type nx = xvec.length ();
72 
73  if (nx < 2)
74  {
75  error ("__pchip_deriv__: X must be at least of length 2");
76  return retval;
77  }
78 
79  octave_idx_type nyr = ymat.rows ();
80  octave_idx_type nyc = ymat.columns ();
81 
82  if (nx != (rows ? nyc : nyr))
83  {
84  error ("__pchip_deriv__: X and Y dimension mismatch");
85  return retval;
86  }
87 
88  FloatMatrix dmat (nyr, nyc);
89 
91  const octave_idx_type incfd = rows ? nyr : 1;
92  volatile const octave_idx_type inc = rows ? 1 : nyr;
93  volatile octave_idx_type k = 0;
94 
95  for (volatile octave_idx_type i = (rows ? nyr : nyc); i > 0; i--)
96  {
97  F77_XFCN (pchim, PCHIM, (nx, xvec.data (),
98  ymat.data () + k * inc,
99  dmat.fortran_vec () + k * inc,
100  incfd, &ierr));
101 
102  k++;
103 
104  if (ierr < 0)
105  {
106  error ("__pchip_deriv__: PCHIM failed with ierr = %i", ierr);
107  return retval;
108  }
109  }
110 
111  retval = dmat;
112  }
113  else
114  {
115  ColumnVector xvec (args(0).vector_value ());
116  Matrix ymat (args(1).matrix_value ());
117 
118  octave_idx_type nx = xvec.length ();
119 
120  if (nx < 2)
121  {
122  error ("__pchip_deriv__: X must be at least of length 2");
123  return retval;
124  }
125 
126  octave_idx_type nyr = ymat.rows ();
127  octave_idx_type nyc = ymat.columns ();
128 
129  if (nx != (rows ? nyc : nyr))
130  {
131  error ("__pchip_deriv__: X and Y dimension mismatch");
132  return retval;
133  }
134 
135  Matrix dmat (nyr, nyc);
136 
138  const octave_idx_type incfd = rows ? nyr : 1;
139  volatile const octave_idx_type inc = rows ? 1 : nyr;
140  volatile octave_idx_type k = 0;
141 
142  for (volatile octave_idx_type i = (rows ? nyr : nyc); i > 0; i--)
143  {
144  F77_XFCN (dpchim, DPCHIM, (nx, xvec.data (),
145  ymat.data () + k * inc,
146  dmat.fortran_vec () + k * inc,
147  incfd, &ierr));
148  k++;
149 
150  if (ierr < 0)
151  {
152  error ("__pchip_deriv__: DPCHIM failed with ierr = %i", ierr);
153  return retval;
154  }
155  }
156 
157  retval = dmat;
158  }
159  }
160 
161  return retval;
162 }
163 
164 /*
165 ## No test needed for internal helper function.
166 %!assert (1)
167 */
subroutine dpchim(N, X, F, D, INCFD, IERR)
Definition: dpchim.f:2
F77_RET_T F77_FUNC(dpchim, DPCHIM)(const octave_idx_type &n
F77_RET_T const double const double double const octave_idx_type octave_idx_type * ierr
#define DEFUN(name, args_name, nargout_name, doc)
Definition: defun.h:44
void error(const char *fmt,...)
Definition: error.cc:476
#define F77_XFCN(f, F, args)
Definition: f77-fcn.h:51
F77_RET_T const double const double double const octave_idx_type & incfd
octave_idx_type rows(void) const
Definition: Array.h:313
F77_RET_T const double const double double * d
F77_RET_T const double const double * f
const T * data(void) const
Definition: Array.h:479
octave_idx_type length(void) const
Definition: ov.cc:1525
#define F77_RET_T
Definition: f77-fcn.h:264
Definition: dMatrix.h:35
subroutine pchim(N, X, F, D, INCFD, IERR)
Definition: pchim.f:2
octave_idx_type length(void) const
Number of elements in the array.
Definition: Array.h:267
const T * fortran_vec(void) const
Definition: Array.h:481
octave_idx_type columns(void) const
Definition: Array.h:322
F77_RET_T const double * x