Mercurial > octave

annotate src/DLD-FUNCTIONS/__pchip_deriv__.cc @ 8712:010e15c7be9a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
support pchip method in interp2
author Jaroslav Hajek <highegg@gmail.com>
date Fri, 14 Mar 2008 12:54:56 +0100
parents 82be108cc558
children 02d4c764de61
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
1 /*
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
2
7017
a1dbe9d80eee [project @ 2007-10-12 21:27:11 by jwe]
jwe
parents: 7016
diff changeset
3 Copyright (C) 2002, 2006, 2007 Kai Habel
8712
010e15c7be9a support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents: 7789
diff changeset
4 Copyright (C) 2008, 2009 Jaroslav Hajek
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
5
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
6 This file is part of Octave.
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
7
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
8 Octave is free software; you can redistribute it and/or modify it
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
9 under the terms of the GNU General Public License as published by the
7016
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6945
diff changeset
10 Free Software Foundation; either version 3 of the License, or (at your
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6945
diff changeset
11 option) any later version.
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
12
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
13 Octave is distributed in the hope that it will be useful, but WITHOUT
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
16 for more details.
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
17
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
18 You should have received a copy of the GNU General Public License
7016
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6945
diff changeset
19 along with Octave; see the file COPYING. If not, see
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6945
diff changeset
20 <http://www.gnu.org/licenses/>.
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
21
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
22 */
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
23
8712
010e15c7be9a support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents: 7789
diff changeset
24 // Jaroslav Hajek, Feb 2008: handle row-wise derivatives,
010e15c7be9a support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents: 7789
diff changeset
25 // use const pointers to avoid unnecessary copying.
010e15c7be9a support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents: 7789
diff changeset
26
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
27 #ifdef HAVE_CONFIG_H
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
28 #include <config.h>
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
29 #endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
30
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
31 #include "defun-dld.h"
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
32 #include "error.h"
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
33 #include "gripes.h"
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
34 #include "oct-obj.h"
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
35 #include "utils.h"
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
36 #include "f77-fcn.h"
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
37
5838
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe]
jwe
parents: 5837
diff changeset
38 extern "C"
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe]
jwe
parents: 5837
diff changeset
39 {
6242
64bad7c6a607 [project @ 2007-01-16 07:03:51 by jwe]
jwe
parents: 5838
diff changeset
40 F77_RET_T
8712
010e15c7be9a support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents: 7789
diff changeset
41 F77_FUNC (dpchim, DPCHIM) (const octave_idx_type& n, const double *x, const double *f,
6242
64bad7c6a607 [project @ 2007-01-16 07:03:51 by jwe]
jwe
parents: 5838
diff changeset
42 double *d, const octave_idx_type &incfd,
64bad7c6a607 [project @ 2007-01-16 07:03:51 by jwe]
jwe
parents: 5838
diff changeset
43 octave_idx_type *ierr);
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
44
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
45 F77_RET_T
8712
010e15c7be9a support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents: 7789
diff changeset
46 F77_FUNC (pchim, PCHIM) (const octave_idx_type& n, const float *x, const float *f,
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
47 float *d, const octave_idx_type &incfd,
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
48 octave_idx_type *ierr);
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
49 }
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
50
6945
6bbf56a9718a [project @ 2007-10-02 20:47:22 by jwe]
jwe
parents: 6549
diff changeset
51 // Wrapper for SLATEC/PCHIP function DPCHIM to calculate the derivates
6bbf56a9718a [project @ 2007-10-02 20:47:22 by jwe]
jwe
parents: 6549
diff changeset
52 // for piecewise polynomials.
6bbf56a9718a [project @ 2007-10-02 20:47:22 by jwe]
jwe
parents: 6549
diff changeset
53
6549
5a5a09d7deb8 [project @ 2007-04-20 06:55:29 by jwe]
jwe
parents: 6242
diff changeset
54 DEFUN_DLD (__pchip_deriv__, args, ,
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
55 "-*- texinfo -*-\n\
8712
010e15c7be9a support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents: 7789
diff changeset
56 @deftypefn {Loadable Function} {} __pchip_deriv__ (@var{x}, @var{y}, @var{dim})\n\
6945
6bbf56a9718a [project @ 2007-10-02 20:47:22 by jwe]
jwe
parents: 6549
diff changeset
57 Undocumented internal function.\n\
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
58 @end deftypefn")
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
59 {
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
60 octave_value retval;
5838
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe]
jwe
parents: 5837
diff changeset
61 const int nargin = args.length ();
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
62
8712
010e15c7be9a support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents: 7789
diff changeset
63 bool rows = (nargin == 3 && args (2).uint_value() == 2);
010e15c7be9a support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents: 7789
diff changeset
64
010e15c7be9a support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents: 7789
diff changeset
65 if (nargin >= 2)
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
66 {
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
67 if (args(0).is_single_type () || args(1).is_single_type ())
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
68 {
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
69 FloatColumnVector xvec (args(0).float_vector_value ());
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
70 FloatMatrix ymat (args(1).float_matrix_value ());
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
71
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
72 octave_idx_type nx = xvec.length ();
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
73 octave_idx_type nyr = ymat.rows ();
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
74 octave_idx_type nyc = ymat.columns ();
5838
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe]
jwe
parents: 5837
diff changeset
75
8712
010e15c7be9a support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents: 7789
diff changeset
76 if (nx != (rows ? nyc : nyr))
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
77 {
8712
010e15c7be9a support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents: 7789
diff changeset
78 error ("__pchip_deriv__: dimension mismatch");
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
79 return retval;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
80 }
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
81
8712
010e15c7be9a support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents: 7789
diff changeset
82 const float *yvec = ymat.data ();
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
83 FloatMatrix dmat (nyr, nyc);
8712
010e15c7be9a support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents: 7789
diff changeset
84 float *dvec = dmat.fortran_vec ();
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
85
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
86 octave_idx_type ierr;
8712
010e15c7be9a support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents: 7789
diff changeset
87 const octave_idx_type incfd = rows ? nyr : 1;
010e15c7be9a support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents: 7789
diff changeset
88 const octave_idx_type inc = rows ? 1 : nyc;
010e15c7be9a support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents: 7789
diff changeset
89
010e15c7be9a support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents: 7789
diff changeset
90 for (octave_idx_type i = (rows ? nyr : nyc); i > 0; i--)
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
91 {
8712
010e15c7be9a support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents: 7789
diff changeset
92 F77_FUNC (pchim, PCHIM) (nx, xvec.data (),
010e15c7be9a support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents: 7789
diff changeset
93 yvec, dvec, incfd, &ierr);
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
94
8712
010e15c7be9a support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents: 7789
diff changeset
95 yvec += inc;
010e15c7be9a support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents: 7789
diff changeset
96 dvec += inc;
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
97
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
98 if (ierr < 0)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
99 {
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
100 error ("PCHIM error: %i\n", ierr);
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
101 return retval;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
102 }
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
103 }
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
104
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
105 retval = dmat;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
106 }
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
107 else
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
108 {
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
109 ColumnVector xvec (args(0).vector_value ());
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
110 Matrix ymat (args(1).matrix_value ());
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
111
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
112 octave_idx_type nx = xvec.length ();
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
113 octave_idx_type nyr = ymat.rows ();
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
114 octave_idx_type nyc = ymat.columns ();
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
115
8712
010e15c7be9a support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents: 7789
diff changeset
116 if (nx != (rows ? nyc : nyr))
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
117 {
8712
010e15c7be9a support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents: 7789
diff changeset
118 error ("__pchip_deriv__: dimension mismatch");
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
119 return retval;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
120 }
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
121
8712
010e15c7be9a support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents: 7789
diff changeset
122 const double *yvec = ymat.data ();
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
123 Matrix dmat (nyr, nyc);
8712
010e15c7be9a support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents: 7789
diff changeset
124 double *dvec = dmat.fortran_vec ();
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
125
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
126 octave_idx_type ierr;
8712
010e15c7be9a support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents: 7789
diff changeset
127 const octave_idx_type incfd = rows ? nyr : 1;
010e15c7be9a support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents: 7789
diff changeset
128 const octave_idx_type inc = rows ? 1 : nyc;
010e15c7be9a support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents: 7789
diff changeset
129
010e15c7be9a support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents: 7789
diff changeset
130 for (octave_idx_type i = (rows ? nyr : nyc); i > 0; i--)
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
131 {
8712
010e15c7be9a support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents: 7789
diff changeset
132 F77_FUNC (dpchim, DPCHIM) (nx, xvec.data (),
010e15c7be9a support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents: 7789
diff changeset
133 yvec, dvec, incfd, &ierr);
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
134
8712
010e15c7be9a support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents: 7789
diff changeset
135 yvec += inc;
010e15c7be9a support pchip method in interp2
Jaroslav Hajek <highegg@gmail.com>
parents: 7789
diff changeset
136 dvec += inc;
5838
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe]
jwe
parents: 5837
diff changeset
137
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
138 if (ierr < 0)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
139 {
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
140 error ("DPCHIM error: %i\n", ierr);
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
141 return retval;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
142 }
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
143 }
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
144
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
145 retval = dmat;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
146 }
5837
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
147 }
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
148
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
149 return retval;
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe]
jwe
parents:
diff changeset
150 }